MFC
Exascale flow solver
Loading...
Searching...
No Matches
m_model.fpp.f90
Go to the documentation of this file.
1# 1 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
2!>
3!! @file
4!! @author Henry Le Berre <hberre3@gatech.edu>
5!! @brief Contains module m_model
6
7# 1 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp" 1
8# 1 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp" 1
9# 1 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp" 1
10# 2 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
11# 3 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
12# 4 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
13# 5 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
14# 6 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
15
16# 8 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
17# 9 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
18# 10 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
19
20# 15 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
21# 16 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
22# 17 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
23
24# 24 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
25
26# 53 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
27
28# 65 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
29
30# 75 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
31
32# 105 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
33
34# 117 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
35
36# 127 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
37
38# 174 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
39! New line at end of file is required for FYPP
40# 2 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp" 2
41# 1 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp" 1
42# 1 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp" 1
43# 2 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
44# 3 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
45# 4 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
46# 5 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
47# 6 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
48
49# 8 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
50# 9 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
51# 10 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
52
53# 15 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
54# 16 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
55# 17 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
56
57# 24 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
58
59# 53 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
60
61# 65 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
62
63# 75 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
64
65# 105 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
66
67# 117 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
68
69# 127 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
70
71# 174 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
72! New line at end of file is required for FYPP
73# 2 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp" 2
74
75# 4 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
76# 5 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
77# 6 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
78# 7 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
79# 8 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
80
81# 20 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
82
83# 43 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
84
85# 48 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
86
87# 53 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
88
89# 58 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
90
91# 63 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
92
93# 68 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
94
95# 81 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
96
97# 86 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
98
99# 91 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
100
101# 96 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
102
103# 101 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
104
105# 106 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
106
107# 111 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
108
109# 116 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
110
111# 121 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
112
113# 126 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
114
115# 156 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
116
117# 197 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
118
119# 211 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
120
121# 236 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
122
123# 247 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
124
125# 249 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
126# 260 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
127
128# 310 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
129
130# 320 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
131
132# 330 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
133
134# 339 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
135
136# 356 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
137
138# 366 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
139
140# 373 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
141
142# 379 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
143
144# 385 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
145
146# 391 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
147
148# 397 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
149
150# 403 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
151! New line at end of file is required for FYPP
152# 3 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp" 2
153# 1 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp" 1
154# 1 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp" 1
155# 2 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
156# 3 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
157# 4 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
158# 5 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
159# 6 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
160
161# 8 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
162# 9 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
163# 10 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
164
165# 15 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
166# 16 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
167# 17 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
168
169# 24 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
170
171# 53 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
172
173# 65 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
174
175# 75 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
176
177# 105 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
178
179# 117 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
180
181# 127 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
182
183# 174 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
184! New line at end of file is required for FYPP
185# 2 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp" 2
186
187# 7 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
188
189# 17 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
190
191# 22 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
192
193# 27 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
194
195# 32 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
196
197# 37 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
198
199# 42 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
200
201# 47 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
202
203# 52 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
204
205# 57 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
206
207# 62 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
208
209# 73 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
210
211# 78 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
212
213# 83 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
214
215# 88 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
216
217# 103 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
218
219# 131 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
220
221# 160 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
222
223# 175 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
224
225# 193 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
226
227# 215 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
228
229# 244 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
230
231# 259 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
232
233# 269 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
234
235# 278 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
236
237# 294 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
238
239# 304 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
240
241# 311 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
242! New line at end of file is required for FYPP
243# 4 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp" 2
244
245! GPU parallel region (scalar reductions, maxval/minval)
246# 23 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
247
248! GPU parallel loop over threads (most common GPU macro)
249# 43 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
250
251! Required closing for GPU_PARALLEL_LOOP
252# 55 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
253
254! Mark routine for device compilation
255# 112 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
256
257! Declare device-resident data
258# 130 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
259
260! Inner loop within a GPU parallel region
261# 145 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
262
263! Scoped GPU data region
264# 164 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
265
266! Host code with device pointers (for MPI with GPU buffers)
267# 193 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
268
269! Allocate device memory (unscoped)
270# 207 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
271
272! Free device memory
273# 219 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
274
275! Atomic operation on device
276# 231 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
277
278! End atomic capture block
279# 242 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
280
281! Copy data between host and device
282# 254 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
283
284! Synchronization barrier
285# 266 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
286
287! Import GPU library module (openacc or omp_lib)
288# 275 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
289
290! Emit code only for AMD compiler
291# 282 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
292
293! Emit code for non-Cray compilers
294# 289 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
295
296! Emit code only for Cray compiler
297# 296 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
298
299! Emit code for non-NVIDIA compilers
300# 303 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
301
302# 305 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
303# 306 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
304! New line at end of file is required for FYPP
305# 2 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp" 2
306
307# 14 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
308
309! Caution: This macro requires the use of a binding script to set CUDA_VISIBLE_DEVICES, such that we have one GPU device per MPI
310! rank. That's because for both cudaMemAdvise (preferred location) and cudaMemPrefetchAsync we use location = device_id = 0. For an
311! example see misc/nvidia_uvm/bind.sh.
312# 52 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
313
314! Allocate and create GPU device memory
315# 72 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
316
317! Free GPU device memory and deallocate
318# 80 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
319
320! Cray-specific GPU pointer setup for vector fields
321# 104 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
322
323! Cray-specific GPU pointer setup for scalar fields
324# 120 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
325
326! Cray-specific GPU pointer setup for acoustic source spatials
327# 145 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
328
329# 151 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
330
331# 158 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
332! New line at end of file is required for FYPP
333# 7 "/home/runner/work/MFC/MFC/src/common/m_model.fpp" 2
334
335!> @brief Binary STL file reader and processor for immersed boundary geometry
336module m_model
337
338 use m_helper
339 use m_mpi_proxy
340 use m_derived_types
341 use iso_c_binding, only: c_char, c_int32_t, c_int16_t, c_float
342
343 implicit none
344
345 private
346
347 public :: f_model_read, s_model_write, s_model_free, models, gpu_ntrs, gpu_trs_v, gpu_trs_n, gpu_boundary_v, &
348 & gpu_boundary_edge_count, stl_bounding_boxes
349
350 ! Subroutines for STL immersed boundaries
351 public :: s_check_boundary, s_register_edge, f_model_is_inside, s_distance_normals_3d, s_distance_normals_2d, &
352 & s_pack_model_for_gpu
353
354 public :: s_instantiate_stl_models
355
356 type(t_model_array), allocatable, target :: models(:) !< STL/OBJ models for IB markers and levelset
357 integer, allocatable :: gpu_ntrs(:) !< GPU-friendly flat arrays for STL model data
358 real(wp), allocatable, dimension(:,:,:,:) :: gpu_trs_v
359 real(wp), allocatable, dimension(:,:,:) :: gpu_trs_n
360 real(wp), allocatable, dimension(:,:,:,:) :: gpu_boundary_v
361 integer, allocatable :: gpu_boundary_edge_count(:)
362 real(wp), allocatable :: stl_bounding_boxes(:,:,:)
363
364# 36 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
365#if defined(MFC_OpenACC)
366# 36 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
367!$acc declare create(gpu_ntrs, gpu_trs_v, gpu_trs_n, gpu_boundary_v, gpu_boundary_edge_count, stl_bounding_boxes)
368# 36 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
369#elif defined(MFC_OpenMP)
370# 36 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
371!$omp declare target (gpu_ntrs, gpu_trs_v, gpu_trs_n, gpu_boundary_v, gpu_boundary_edge_count, stl_bounding_boxes)
372# 36 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
373#endif
374
375contains
376
377 !> Read a binary STL file.
378 impure subroutine s_read_stl_binary(filepath, model)
379
380 character(LEN=*), intent(in) :: filepath
381 type(t_model), intent(out) :: model
382 integer :: i, iunit, iostat
383 character(kind=c_char, len=80) :: header
384 integer(kind=c_int32_t) :: ntriangles
385 real(kind=c_float) :: normal(3), v(3, 3), v_norm
386 integer(kind=c_int16_t) :: attribute
387
388 open (newunit=iunit, file=filepath, action='READ', form='UNFORMATTED', status='OLD', iostat=iostat, access='STREAM')
389
390 if (iostat /= 0) then
391 print *, "Error: could not open Binary STL file ", filepath
392
393 call s_mpi_abort()
394 end if
395
396 read (iunit, iostat=iostat) header, ntriangles
397
398 if (iostat /= 0) then
399 print *, "Error: could not read header from Binary STL file ", filepath
400
401 call s_mpi_abort()
402 end if
403
404 model%ntrs = ntriangles
405
406 allocate (model%trs(model%ntrs))
407
408 do i = 1, model%ntrs
409 read (iunit) normal(:), v(1,:), v(2,:), v(3,:), attribute
410
411 model%trs(i)%v = v
412 model%trs(i)%n = normal
413 v_norm = sqrt(normal(1)**2 + normal(2)**2 + normal(3)**2)
414 if (v_norm > 0._wp) model%trs(i)%n = normal/v_norm
415 end do
416
417 close (iunit)
418
419 end subroutine s_read_stl_binary
420
421 !> Read an ASCII STL file.
422 impure subroutine s_read_stl_ascii(filepath, model)
423
424 character(LEN=*), intent(in) :: filepath
425 type(t_model), intent(out) :: model
426 integer :: i, j, iunit, iostat
427 character(80) :: line, buffered_line
428 logical :: is_buffered
429 real(wp) :: normal(3), v_norm
430
431 is_buffered = .false.
432
433 open (newunit=iunit, file=filepath, action='READ', form='FORMATTED', status='OLD', iostat=iostat, access='STREAM')
434
435 if (iostat /= 0) then
436 print *, "Error: could not open ASCII STL file ", filepath
437 call s_mpi_abort()
438 end if
439
440 model%ntrs = 0
441 do
442 if (is_buffered) then
443 line = buffered_line
444 is_buffered = .false.
445 else
446 if (.not. f_read_line(iunit, line)) exit
447 end if
448
449 if (line(1:6) == "facet ") then
450 model%ntrs = model%ntrs + 1
451 end if
452 end do
453
454 allocate (model%trs(model%ntrs))
455
456 rewind(iunit)
457
458 i = 1
459 do
460 if (is_buffered) then
461 line = buffered_line
462 is_buffered = .false.
463 else
464 if (.not. f_read_line(iunit, line)) exit
465 end if
466
467 if (line(1:5) == "solid") cycle
468 if (line(1:8) == "endsolid") exit
469
470 if (line(1:12) /= "facet normal") then
471 print *, "Error: expected facet normal in STL file ", filepath
472 call s_mpi_abort()
473 end if
474
475 call s_skip_ignored_lines(iunit, buffered_line, is_buffered)
476 read (line(13:), *) normal
477 v_norm = sqrt(normal(1)**2 + normal(2)**2 + normal(3)**2)
478 if (v_norm > 0._wp) model%trs(i)%n = normal/v_norm
479
480 call s_skip_ignored_lines(iunit, buffered_line, is_buffered)
481 if (is_buffered) then
482 line = buffered_line
483 is_buffered = .false.
484 else
485 read (iunit, '(A)') line
486 end if
487
488 do j = 1, 3
489 if (is_buffered) then
490 line = buffered_line
491 is_buffered = .false.
492 else
493 if (.not. f_read_line(iunit, line)) exit
494 end if
495
496 if (line(1:6) /= "vertex") then
497 print *, "Error: expected vertex in STL file ", filepath
498 call s_mpi_abort()
499 end if
500
501 call s_skip_ignored_lines(iunit, buffered_line, is_buffered)
502 read (line(7:), *) model%trs(i)%v(j,:)
503 end do
504
505 if (is_buffered) then
506 line = buffered_line
507 is_buffered = .false.
508 else
509 if (.not. f_read_line(iunit, line)) exit
510 end if
511
512 if (is_buffered) then
513 line = buffered_line
514 is_buffered = .false.
515 else
516 if (.not. f_read_line(iunit, line)) exit
517 end if
518
519 if (line(1:8) /= "endfacet") then
520 print *, "Error: expected endfacet in STL file ", filepath
521 call s_mpi_abort()
522 end if
523
524 i = i + 1
525 end do
526
527 end subroutine s_read_stl_ascii
528
529 !> Read an STL file.
530 impure subroutine s_read_stl(filepath, model)
531
532 character(LEN=*), intent(in) :: filepath
533 type(t_model), intent(out) :: model
534 integer :: iunit, iostat
535 character(80) :: line
536
537 open (newunit=iunit, file=filepath, action='READ', form='FORMATTED', status='OLD', iostat=iostat, access='STREAM')
538
539 if (iostat /= 0) then
540 print *, "Error: could not open STL file ", filepath
541
542 call s_mpi_abort()
543 end if
544
545 read (iunit, '(A)') line
546
547 close (iunit)
548
549 if (line(1:5) == "solid") then
550 call s_read_stl_ascii(filepath, model)
551 else
552 call s_read_stl_binary(filepath, model)
553 end if
554
555 end subroutine s_read_stl
556
557 !> Read an OBJ file.
558 impure subroutine s_read_obj(filepath, model)
559
560 character(LEN=*), intent(in) :: filepath
561 type(t_model), intent(out) :: model
562 integer :: i, j, k, l, iv3, iunit, iostat, nvertices
563 real(wp), dimension(1:3), allocatable :: vertices(:,:)
564 character(80) :: line
565
566 open (newunit=iunit, file=filepath, action='READ', form='FORMATTED', status='OLD', iostat=iostat, access='STREAM')
567
568 if (iostat /= 0) then
569 print *, "Error: could not open model file ", filepath
570
571 call s_mpi_abort()
572 end if
573
574 nvertices = 0
575 model%ntrs = 0
576 do
577 if (.not. f_read_line(iunit, line)) exit
578
579 select case (line(1:2))
580 case ("v ")
581 nvertices = nvertices + 1
582 case ("f ")
583 model%ntrs = model%ntrs + 1
584 end select
585 end do
586
587 rewind(iunit)
588
589 allocate (vertices(nvertices,1:3))
590 allocate (model%trs(model%ntrs))
591
592 i = 1
593 j = 1
594
595 do
596 if (.not. f_read_line(iunit, line)) exit
597
598 select case (line(1:2))
599 case ("g ")
600 case ("vn")
601 case ("vt")
602 case ("l ")
603 case ("v ")
604 read (line(3:), *) vertices(i,:)
605 i = i + 1
606 case ("f ")
607 read (line(3:), *) k, l, iv3
608 model%trs(j)%v(1,:) = vertices(k,:)
609 model%trs(j)%v(2,:) = vertices(l,:)
610 model%trs(j)%v(3,:) = vertices(iv3,:)
611 j = j + 1
612 case default
613 print *, "Error: unknown line type in OBJ file ", filepath
614 print *, "Line: ", line
615
616 call s_mpi_abort()
617 end select
618 end do
619
620 deallocate (vertices)
621
622 close (iunit)
623
624 end subroutine s_read_obj
625
626 !> Read a mesh from a file.
627 impure function f_model_read(filepath) result(model)
628
629 character(LEN=*), intent(in) :: filepath
630 type(t_model) :: model
631
632 select case (filepath(len(trim(filepath)) - 3:len(trim(filepath))))
633 case (".stl")
634 call s_read_stl(filepath, model)
635 case (".obj")
636 call s_read_obj(filepath, model)
637 case default
638 print *, "Error: unknown model file format for file ", filepath
639
640 call s_mpi_abort()
641 end select
642
643 end function f_model_read
644
645 !> Write a binary STL file.
646 impure subroutine s_write_stl(filepath, model)
647
648 character(LEN=*), intent(in) :: filepath
649 type(t_model), intent(in) :: model
650 integer :: i, j, iunit, iostat
651 character(kind=c_char, len=80), parameter :: header = "Model file written by MFC."
652 integer(kind=c_int32_t) :: ntriangles
653 real(wp) :: normal(3), v(3)
654 integer(kind=c_int16_t) :: attribute
655
656 open (newunit=iunit, file=filepath, action='WRITE', form='UNFORMATTED', iostat=iostat, access='STREAM')
657
658 if (iostat /= 0) then
659 print *, "Error: could not open STL file ", filepath
660
661 call s_mpi_abort()
662 end if
663
664 ntriangles = model%ntrs
665 write (iunit, iostat=iostat) header, ntriangles
666
667 if (iostat /= 0) then
668 print *, "Error: could not write header to STL file ", filepath
669
670 call s_mpi_abort()
671 end if
672
673 do i = 1, model%ntrs
674 normal = model%trs(i)%n
675 write (iunit) normal
676
677 do j = 1, 3
678 v = model%trs(i)%v(j,:)
679 write (iunit) v(:)
680 end do
681
682 attribute = 0
683 write (iunit) attribute
684 end do
685
686 close (iunit)
687
688 end subroutine s_write_stl
689
690 !> Write an OBJ file.
691 impure subroutine s_write_obj(filepath, model)
692
693 character(LEN=*), intent(in) :: filepath
694 type(t_model), intent(in) :: model
695 integer :: iunit, iostat
696 integer :: i, j
697
698 open (newunit=iunit, file=filepath, action='WRITE', form='FORMATTED', iostat=iostat, access='STREAM')
699
700 if (iostat /= 0) then
701 print *, "Error: could not open OBJ file ", filepath
702
703 call s_mpi_abort()
704 end if
705
706 write (iunit, '(A)') "# Model file written by MFC."
707
708 do i = 1, model%ntrs
709 do j = 1, 3
710 write (iunit, '(A, " ", (f30.20), " ", (f30.20), " ", (f30.20))') "v", model%trs(i)%v(j, 1), model%trs(i)%v(j, &
711 & 2), model%trs(i)%v(j, 3)
712 end do
713
714 write (iunit, '(A, " ", I0, " ", I0, " ", I0)') "f", i*3 - 2, i*3 - 1, i*3
715 end do
716
717 close (iunit)
718
719 end subroutine s_write_obj
720
721 !> Write a mesh to a file.
722 impure subroutine s_model_write(filepath, model)
723
724 character(LEN=*), intent(in) :: filepath
725 type(t_model), intent(in) :: model
726
727 select case (filepath(len(trim(filepath)) - 3:len(trim(filepath))))
728 case (".stl")
729 call s_write_stl(filepath, model)
730 case (".obj")
731 call s_write_obj(filepath, model)
732 case default
733 print *, "Error: unknown model file format for file ", filepath
734
735 call s_mpi_abort()
736 end select
737
738 end subroutine s_model_write
739
740 !> Free the memory allocated for an STL mesh.
741 subroutine s_model_free(model)
742
743 type(t_model), intent(inout) :: model
744
745 deallocate (model%trs)
746
747 end subroutine s_model_free
748
749 !> Read the next non-blank, non-comment line from an STL or OBJ model file.
750 impure function f_read_line(iunit, line) result(bIsLine)
751
752 integer, intent(in) :: iunit
753 character(80), intent(out) :: line
754 logical :: bisline
755 integer :: iostat
756
757 bisline = .true.
758
759 do
760 read (iunit, '(A)', iostat=iostat) line
761
762 if (iostat < 0) then
763 bisline = .false.
764 exit
765 end if
766
767 line = adjustl(trim(line))
768
769 if (len(trim(line)) == 0) cycle
770 if (line(1:5) == "solid") cycle
771 if (line(1:1) == "#") cycle
772
773 exit
774 end do
775
776 end function f_read_line
777
778 !> Read the next non-comment line from a model file, using a buffered look-ahead mechanism.
779 impure subroutine s_skip_ignored_lines(iunit, buffered_line, is_buffered)
780
781 integer, intent(in) :: iunit
782 character(80), intent(inout) :: buffered_line
783 logical, intent(inout) :: is_buffered
784 character(80) :: line
785
786 if (is_buffered) then
787 line = buffered_line
788 is_buffered = .false.
789 else
790 if (.not. f_read_line(iunit, line)) return
791 end if
792
793 buffered_line = line
794 is_buffered = .true.
795
796 end subroutine s_skip_ignored_lines
797
798 !> Determine if a point is inside a surface using the generalized winding number (Jacobson et al., SIGGRAPH 2013). In 3D, sums
799 !! the solid angle subtended by each triangle (Van Oosterom-Strackee formula). In 2D (p==0), sums the signed angle subtended by
800 !! each boundary edge. Returns ~1.0 inside, ~0.0 outside. Unlike ray casting, this is robust to small triangles/edges and vertex
801 !! winding order.
802 !! @return fraction Winding number (~1.0 inside, ~0.0 outside).
803 function f_model_is_inside(ntrs, pid, point) result(fraction)
804
805
806# 468 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
807#if MFC_OpenACC
808# 468 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
809!$acc routine seq
810# 468 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
811#elif MFC_OpenMP
812# 468 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
813
814# 468 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
815
816# 468 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
817!$omp declare target device_type(any)
818# 468 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
819#endif
820
821 integer, intent(in) :: ntrs
822 integer, intent(in) :: pid
823 real(wp), dimension(1:3), intent(in) :: point
824 real(wp) :: fraction
825 real(wp) :: r1(3), r2(3), r3(3)
826 real(wp) :: r1_mag, r2_mag, r3_mag
827 real(wp) :: numerator, denominator
828 real(wp) :: d1(2), d2(2)
829 integer :: q
830
831 fraction = 0.0_wp
832
833 if (p == 0) then
834 ! 2D winding number: sum signed angles subtended by each boundary edge at the query point.
835 do q = 1, gpu_boundary_edge_count(pid)
836 d1(1) = gpu_boundary_v(q, 1, 1, pid) - point(1)
837 d1(2) = gpu_boundary_v(q, 1, 2, pid) - point(2)
838 d2(1) = gpu_boundary_v(q, 2, 1, pid) - point(1)
839 d2(2) = gpu_boundary_v(q, 2, 2, pid) - point(2)
840
841 ! Signed angle = atan2(d1 x d2, d1 . d2)
842 fraction = fraction + atan2(d1(1)*d2(2) - d1(2)*d2(1), d1(1)*d2(1) + d1(2)*d2(2))
843 end do
844
845 ! 2D winding number = total angle / (2*pi)
846 fraction = fraction/(2.0_wp*acos(-1.0_wp))
847 else
848 ! 3D winding number: sum solid angles via Van Oosterom-Strackee formula.
849 do q = 1, ntrs
850 r1 = gpu_trs_v(1,:,q, pid) - point
851 r2 = gpu_trs_v(2,:,q, pid) - point
852 r3 = gpu_trs_v(3,:,q, pid) - point
853
854 r1_mag = sqrt(dot_product(r1, r1))
855 r2_mag = sqrt(dot_product(r2, r2))
856 r3_mag = sqrt(dot_product(r3, r3))
857
858 ! Skip if query point is coincident with a vertex (magnitudes are zero/subnormal).
859 if (r1_mag*r2_mag*r3_mag < tiny(1.0_wp)) cycle
860
861 ! tan(Omega/2) = numerator / denominator numerator = scalar triple product r1 . (r2 x r3)
862 numerator = r1(1)*(r2(2)*r3(3) - r2(3)*r3(2)) + r1(2)*(r2(3)*r3(1) - r2(1)*r3(3)) + r1(3)*(r2(1)*r3(2) - r2(2) &
863 & *r3(1))
864
865 denominator = r1_mag*r2_mag*r3_mag + dot_product(r1, r2)*r3_mag + dot_product(r2, r3)*r1_mag + dot_product(r3, &
866 & r1)*r2_mag
867
868 fraction = fraction + atan2(numerator, denominator)
869 end do
870
871 ! Each atan2 returns Omega/2 per triangle; divide by 2*pi to get winding number = sum(Omega)/(4*pi).
872 fraction = fraction/(2.0_wp*acos(-1.0_wp))
873 end if
874
875 end function f_model_is_inside
876
877 !> Check and label edges shared by two or more triangle facets of the 2D STL model.
878 subroutine s_check_boundary(model, boundary_v, boundary_vertex_count, boundary_edge_count)
879
880 type(t_model), intent(in) :: model
881 real(wp), allocatable, intent(out), dimension(:,:,:) :: boundary_v !< Output boundary vertices/normals
882 integer, intent(out) :: boundary_vertex_count, boundary_edge_count !< Output boundary vertex/edge count
883 integer :: i, j !< Model index iterator
884 integer :: edge_count, edge_index, store_index !< Boundary edge index iterator
885 real(wp), dimension(1:2,1:2) :: edge !< Edge end points buffer
886 real(wp), dimension(1:2) :: boundary_edge !< Boundary edge end points buffer
887 real(wp), dimension(1:(3*model%ntrs),1:2,1:2) :: temp_boundary_v !< Temporary boundary vertex buffer
888 integer, dimension(1:(3*model%ntrs)) :: edge_occurrence !< The manifoldness of the edges
889 real(wp) :: edgetan, initial, v_norm, xnormal, ynormal !< The manifoldness of the edges
890 ! Total number of edges in 2D STL
891
892 edge_count = 3*model%ntrs
893
894 ! Initialize edge_occurrence array to zero
895 edge_occurrence = 0
896 edge_index = 0
897
898 ! Collect all edges of all triangles and store them
899 do i = 1, model%ntrs
900 ! First edge (v1, v2)
901 edge(1,1:2) = model%trs(i)%v(1,1:2)
902 edge(2,1:2) = model%trs(i)%v(2,1:2)
903 call s_register_edge(temp_boundary_v, edge, edge_index, edge_count)
904
905 ! Second edge (v2, v3)
906 edge(1,1:2) = model%trs(i)%v(2,1:2)
907 edge(2,1:2) = model%trs(i)%v(3,1:2)
908 call s_register_edge(temp_boundary_v, edge, edge_index, edge_count)
909
910 ! Third edge (v3, v1)
911 edge(1,1:2) = model%trs(i)%v(3,1:2)
912 edge(2,1:2) = model%trs(i)%v(1,1:2)
913 call s_register_edge(temp_boundary_v, edge, edge_index, edge_count)
914 end do
915
916 ! Check all edges and count repeated edges
917
918# 566 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
919
920# 566 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
921#if defined(MFC_OpenACC)
922# 566 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
923!$acc parallel loop collapse(2) gang vector default(present) private(i, j) copy(temp_boundary_v, edge_occurrence)
924# 566 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
925#elif defined(MFC_OpenMP)
926# 566 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
927
928# 566 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
929
930# 566 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
931
932# 566 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
933!$omp target teams loop defaultmap(firstprivate:scalar) bind(teams,parallel) collapse(2) defaultmap(tofrom:aggregate) defaultmap(tofrom:allocatable) defaultmap(tofrom:pointer) private(i, j) &
934# 566 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
935!$omp& map(tofrom:temp_boundary_v, edge_occurrence)
936# 566 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
937#endif
938 do i = 1, edge_count
939 do j = 1, edge_count
940 if (i /= j) then
941 if (((abs(temp_boundary_v(i, 1, 1) - temp_boundary_v(j, 1, &
942 & 1)) < threshold_edge_zero) .and. (abs(temp_boundary_v(i, 1, 2) - temp_boundary_v(j, 1, &
943 & 2)) < threshold_edge_zero) .and. (abs(temp_boundary_v(i, 2, 1) - temp_boundary_v(j, 2, &
944 & 1)) < threshold_edge_zero) .and. (abs(temp_boundary_v(i, 2, 2) - temp_boundary_v(j, 2, &
945 & 2)) < threshold_edge_zero)) .or. ((abs(temp_boundary_v(i, 1, 1) - temp_boundary_v(j, 2, &
946 & 1)) < threshold_edge_zero) .and. (abs(temp_boundary_v(i, 1, 2) - temp_boundary_v(j, 2, &
947 & 2)) < threshold_edge_zero) .and. (abs(temp_boundary_v(i, 2, 1) - temp_boundary_v(j, 1, &
948 & 1)) < threshold_edge_zero) .and. (abs(temp_boundary_v(i, 2, 2) - temp_boundary_v(j, 1, &
949 & 2)) < threshold_edge_zero))) then
950
951# 579 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
952#if defined(MFC_OpenACC)
953# 579 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
954!$acc atomic update
955# 579 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
956#elif defined(MFC_OpenMP)
957# 579 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
958!$omp atomic update
959# 579 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
960#endif
961 edge_occurrence(i) = edge_occurrence(i) + 1
962 end if
963 end if
964 end do
965 end do
966
967# 585 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
968#if defined(MFC_OpenACC)
969# 585 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
970!$acc end parallel loop
971# 585 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
972#elif defined(MFC_OpenMP)
973# 585 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
974
975# 585 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
976!$omp end target teams loop
977# 585 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
978#endif
979
980 ! Count the number of boundary vertices/edges
981 boundary_vertex_count = 0
982 boundary_edge_count = 0
983
984 do i = 1, edge_count
985 if (edge_occurrence(i) == 0) then
986 boundary_vertex_count = boundary_vertex_count + 2
987 boundary_edge_count = boundary_edge_count + 1
988 end if
989 end do
990
991 ! Allocate the boundary_v array based on the number of boundary edges
992 allocate (boundary_v(boundary_edge_count,1:3,1:2))
993
994 ! Store boundary vertices
995 store_index = 0
996 do i = 1, edge_count
997 if (edge_occurrence(i) == 0) then
998 store_index = store_index + 1
999 boundary_v(store_index, 1,1:2) = temp_boundary_v(i, 1,1:2)
1000 boundary_v(store_index, 2,1:2) = temp_boundary_v(i, 2,1:2)
1001 end if
1002 end do
1003
1004 ! Find/store the normal vector of the boundary edges
1005 do i = 1, boundary_edge_count
1006 boundary_edge(1) = boundary_v(i, 2, 1) - boundary_v(i, 1, 1)
1007 boundary_edge(2) = boundary_v(i, 2, 2) - boundary_v(i, 1, 2)
1008 edgetan = boundary_edge(1)/sign(max(sgm_eps, abs(boundary_edge(2))), boundary_edge(2))
1009
1010 if (abs(boundary_edge(2)) < threshold_vector_zero) then
1011 if (edgetan > 0._wp) then
1012 ynormal = -1
1013 xnormal = 0._wp
1014 else
1015 ynormal = 1
1016 xnormal = 0._wp
1017 end if
1018 else
1019 initial = boundary_edge(2)
1020 ynormal = -edgetan*initial
1021 xnormal = initial
1022 end if
1023
1024 v_norm = sqrt(xnormal**2 + ynormal**2)
1025 boundary_v(i, 3, 1) = xnormal/v_norm
1026 boundary_v(i, 3, 2) = ynormal/v_norm
1027 end do
1028
1029 end subroutine s_check_boundary
1030
1031 !> Append the edge end vertices to a temporary buffer.
1032 subroutine s_register_edge(temp_boundary_v, edge, edge_index, edge_count)
1033
1034 integer, intent(inout) :: edge_index !< Edge index iterator
1035 integer, intent(inout) :: edge_count !< Total number of edges
1036 real(wp), intent(in), dimension(1:2,1:2) :: edge !< Edges end points to be registered
1037 real(wp), dimension(1:edge_count,1:2,1:2), intent(inout) :: temp_boundary_v !< Temporary edge end vertex buffer
1038 ! Increment edge index and store the edge
1039
1040 edge_index = edge_index + 1
1041 temp_boundary_v(edge_index, 1,1:2) = edge(1,1:2)
1042 temp_boundary_v(edge_index, 2,1:2) = edge(2,1:2)
1043
1044 end subroutine s_register_edge
1045
1046 !> Determine the levelset distance and normals of 3D models by computing the exact closest point via projection onto triangle
1047 !! surfaces.
1048 subroutine s_distance_normals_3d(ntrs, pid, point, normals, distance)
1049
1050
1051# 657 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1052#if MFC_OpenACC
1053# 657 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1054!$acc routine seq
1055# 657 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1056#elif MFC_OpenMP
1057# 657 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1058
1059# 657 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1060
1061# 657 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1062!$omp declare target device_type(any)
1063# 657 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1064#endif
1065
1066 integer, intent(in) :: ntrs
1067 integer, intent(in) :: pid
1068 real(wp), dimension(1:3), intent(in) :: point
1069 real(wp), dimension(1:3), intent(out) :: normals
1070 real(wp), intent(out) :: distance
1071 integer :: i, j, l
1072 real(wp) :: dist_min, dist_proj, dist_v, dist_e, t
1073 real(wp) :: v1(1:3), v2(1:3), v3(1:3)
1074 real(wp) :: e0(1:3), e1(1:3), pv(1:3)
1075 real(wp) :: n(1:3), proj(1:3), norm_vec(1:3)
1076 real(wp) :: d, denom, norm_mag
1077 real(wp) :: u, v_bary, w
1078 real(wp) :: l00, l01, l11, l20, l21
1079 real(wp) :: edge(1:3), pe(1:3)
1080 real(wp) :: verts(1:3,1:3)
1081
1082 dist_min = initial_distance_buffer
1083 normals = 0._wp
1084
1085 do i = 1, ntrs
1086 ! Triangle vertices
1087 v1(:) = gpu_trs_v(1,:,i, pid)
1088 v2(:) = gpu_trs_v(2,:,i, pid)
1089 v3(:) = gpu_trs_v(3,:,i, pid)
1090
1091 ! Triangle normal
1092 n(:) = gpu_trs_n(:,i, pid)
1093
1094 ! Project point onto triangle plane
1095 pv(:) = point(:) - v1(:)
1096 d = dot_product(pv, n)
1097 if (abs(d) >= dist_min) cycle ! minimum distance is not small enough, no need to check validity
1098 proj(:) = point(:) - d*n(:)
1099
1100 ! Check if projection is inside triangle using barycentric coordinates
1101 e0(:) = v2(:) - v1(:)
1102 e1(:) = v3(:) - v1(:)
1103 pv(:) = proj(:) - v1(:)
1104
1105 l00 = dot_product(e0, e0)
1106 l01 = dot_product(e0, e1)
1107 l11 = dot_product(e1, e1)
1108 l20 = dot_product(pv, e0)
1109 l21 = dot_product(pv, e1)
1110
1111 denom = l00*l11 - l01*l01
1112
1113 ! compute the barycentric coordinates of the projection in the triangle
1114 if (abs(denom) > 0._wp) then
1115 v_bary = (l11*l20 - l01*l21)/denom
1116 w = (l00*l21 - l01*l20)/denom
1117 u = 1._wp - v_bary - w
1118 else
1119 u = -1._wp
1120 v_bary = -1._wp
1121 w = -1._wp
1122 end if
1123
1124 ! If projection is inside triangle
1125 if (u >= 0._wp .and. v_bary >= 0._wp .and. w >= 0._wp) then
1126 dist_proj = sqrt((point(1) - proj(1))**2 + (point(2) - proj(2))**2 + (point(3) - proj(3))**2)
1127
1128 if (dist_proj < dist_min) then
1129 dist_min = dist_proj
1130 normals(:) = n(:)
1131 end if
1132 else
1133 ! Projection outside triangle: check edges and vertices
1134 verts(:,1) = v1(:)
1135 verts(:,2) = v2(:)
1136 verts(:,3) = v3(:)
1137
1138 ! Check three edges
1139 do j = 1, 3
1140 edge(:) = verts(:,mod(j, 3) + 1) - verts(:,j)
1141 pe(:) = point(:) - verts(:,j)
1142
1143 t = dot_product(pe, edge)/max(dot_product(edge, edge), 1.e-30_wp)
1144
1145 if (t >= 0._wp .and. t <= 1._wp) then
1146 proj(:) = verts(:,j) + t*edge(:)
1147 dist_e = sqrt((point(1) - proj(1))**2 + (point(2) - proj(2))**2 + (point(3) - proj(3))**2)
1148
1149 if (dist_e < dist_min) then
1150 dist_min = dist_e
1151 norm_vec(:) = proj(:) - point(:)
1152 if (dist_e > 0._wp) norm_vec = norm_vec/dist_e
1153 ! Snap to triangle normal if nearly parallel
1154 if (f_approx_equal(dot_product(norm_vec, n), 1._wp)) then
1155 normals(:) = n(:)
1156 else
1157 normals(:) = norm_vec(:)
1158 end if
1159 end if
1160 else if (t < 0._wp) then
1161 dist_v = sqrt((point(1) - verts(1, j))**2 + (point(2) - verts(2, j))**2 + (point(3) - verts(3, j))**2)
1162
1163 if (dist_v < dist_min) then
1164 dist_min = dist_v
1165 norm_vec(:) = verts(:,j) - point(:)
1166 norm_mag = sqrt(dot_product(norm_vec, norm_vec))
1167 if (norm_mag > 0._wp) norm_vec = norm_vec/norm_mag
1168 normals(:) = norm_vec(:)
1169 end if
1170 else
1171 dist_v = sqrt((point(1) - verts(1, mod(j, 3) + 1))**2 + (point(2) - verts(2, mod(j, &
1172 & 3) + 1))**2 + (point(3) - verts(3, mod(j, 3) + 1))**2)
1173
1174 if (dist_v < dist_min) then
1175 dist_min = dist_v
1176 norm_vec(:) = verts(:,mod(j, 3) + 1) - point(:)
1177 norm_mag = sqrt(dot_product(norm_vec, norm_vec))
1178 if (norm_mag > 0._wp) norm_vec = norm_vec/norm_mag
1179 normals(:) = norm_vec(:)
1180 end if
1181 end if
1182 end do
1183 end if
1184 end do
1185
1186 distance = dist_min
1187
1188 end subroutine s_distance_normals_3d
1189
1190 !> Determine the levelset distance and normals of 2D models by computing the exact closest point via projection onto boundary
1191 !! edges.
1192 subroutine s_distance_normals_2d(pid, boundary_edge_count, point, normals, distance)
1193
1194
1195# 787 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1196#if MFC_OpenACC
1197# 787 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1198!$acc routine seq
1199# 787 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1200#elif MFC_OpenMP
1201# 787 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1202
1203# 787 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1204
1205# 787 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1206!$omp declare target device_type(any)
1207# 787 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1208#endif
1209
1210 integer, intent(in) :: pid
1211 integer, intent(in) :: boundary_edge_count
1212 real(wp), dimension(1:3), intent(in) :: point
1213 real(wp), dimension(1:3), intent(out) :: normals
1214 real(wp), intent(out) :: distance
1215 integer :: i
1216 real(wp) :: dist_min, dist, t
1217 real(wp) :: v1(1:2), v2(1:2), edge(1:2), pv(1:2)
1218 real(wp) :: edge_len_sq, proj(1:2), norm(1:2)
1219
1220 dist_min = initial_distance_buffer
1221 normals = 0._wp
1222 norm = 0._wp
1223
1224 do i = 1, boundary_edge_count
1225 ! Edge endpoints
1226 v1(1) = gpu_boundary_v(i, 1, 1, pid)
1227 v1(2) = gpu_boundary_v(i, 1, 2, pid)
1228 v2(1) = gpu_boundary_v(i, 2, 1, pid)
1229 v2(2) = gpu_boundary_v(i, 2, 2, pid)
1230
1231 ! Edge vector and point-to-v1 vector
1232 edge = v2 - v1
1233 pv(1) = point(1) - v1(1)
1234 pv(2) = point(2) - v1(2)
1235 edge_len_sq = dot_product(edge, edge)
1236
1237 ! Parameter of projection onto the edge line
1238 if (edge_len_sq > 0._wp) then
1239 t = dot_product(pv, edge)/edge_len_sq
1240 else
1241 t = 0._wp
1242 end if
1243
1244 ! Check if projection falls within the segment
1245 if (t >= 0._wp .and. t <= 1._wp) then
1246 proj = v1 + t*edge
1247 dist = sqrt((point(1) - proj(1))**2 + (point(2) - proj(2))**2)
1248 norm(1) = gpu_boundary_v(i, 3, 1, pid)
1249 norm(2) = gpu_boundary_v(i, 3, 2, pid)
1250 else if (t < 0._wp) then ! negative t means that v1 is the closest point on the edge
1251 dist = sqrt((point(1) - v1(1))**2 + (point(2) - v1(2))**2)
1252 norm(1) = v1(1) - point(1)
1253 norm(2) = v1(2) - point(2)
1254 norm = norm/dist
1255 else ! t > 1 means that v2 is the closest point on the line edge
1256 dist = sqrt((point(1) - v2(1))**2 + (point(2) - v2(2))**2)
1257 norm(1) = v2(1) - point(1)
1258 norm(2) = v2(2) - point(2)
1259 norm = norm/dist
1260 end if
1261
1262 if (dist < dist_min) then
1263 dist_min = dist
1264 normals(1) = norm(1)
1265 normals(2) = norm(2)
1266 end if
1267 end do
1268
1269 distance = dist_min
1270
1271 end subroutine s_distance_normals_2d
1272
1273 !> Load, transform, and register STL/OBJ immersed-boundary models onto the simulation grid.
1274 subroutine s_instantiate_stl_models()
1275
1276 ! Variables for IBM+STL
1277 integer :: boundary_vertex_count, boundary_edge_count !< Boundary vertex
1278 real(wp), allocatable, dimension(:,:,:) :: boundary_v !< Boundary vertex buffer
1279 integer :: i, j, k !< Generic loop iterators
1280 integer :: stl_id
1281 type(t_bbox) :: bbox, bbox_old
1282 type(t_model) :: model
1283 type(ic_model_parameters) :: params
1284 real(wp), dimension(1:3) :: point, model_center
1285 real(wp) :: grid_mm(1:3,1:2)
1286 real(wp), dimension(1:4,1:4) :: transform, transform_n
1287
1288 if (num_stl_models == 0) return
1289
1290#ifdef MFC_DEBUG
1291# 869 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1292 block
1293# 869 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1294 use iso_fortran_env, only: output_unit
1295# 869 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1296
1297# 869 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1298 print *, 'm_model.fpp:869: ', '@:ALLOCATE(stl_bounding_boxes(num_stl_models,1:3,1:3))'
1299# 869 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1300
1301# 869 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1302 call flush (output_unit)
1303# 869 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1304 end block
1305# 869 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1306#endif
1307# 869 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1308 allocate (stl_bounding_boxes(num_stl_models,1:3,1:3))
1309# 869 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1310
1311# 869 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1312
1313# 869 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1314#if defined(MFC_OpenACC)
1315# 869 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1316!$acc enter data create(stl_bounding_boxes)
1317# 869 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1318#elif defined(MFC_OpenMP)
1319# 869 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1320!$omp target enter data map(always,alloc:stl_bounding_boxes)
1321# 869 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1322#endif
1323#ifdef MFC_DEBUG
1324# 870 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1325 block
1326# 870 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1327 use iso_fortran_env, only: output_unit
1328# 870 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1329
1330# 870 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1331 print *, 'm_model.fpp:870: ', '@:ALLOCATE(models(num_stl_models))'
1332# 870 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1333
1334# 870 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1335 call flush (output_unit)
1336# 870 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1337 end block
1338# 870 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1339#endif
1340# 870 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1341 allocate (models(num_stl_models))
1342# 870 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1343
1344# 870 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1345
1346# 870 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1347#if defined(MFC_OpenACC)
1348# 870 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1349!$acc enter data create(models)
1350# 870 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1351#elif defined(MFC_OpenMP)
1352# 870 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1353!$omp target enter data map(always,alloc:models)
1354# 870 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1355#endif
1356
1357 do stl_id = 1, num_stl_models
1358 allocate (models(stl_id)%model)
1359 if (proc_rank == 0) print *, " * Reading model: " // trim(stl_models(stl_id)%model_filepath)
1360
1361 model = f_model_read(stl_models(stl_id)%model_filepath)
1362 params%scale(:) = stl_models(stl_id)%model_scale(:)
1363 params%translate(:) = stl_models(stl_id)%model_translate(:)
1364 params%rotate(:) = 0._wp
1365 params%spc = num_ray
1366 params%threshold = stl_models(stl_id)%model_threshold
1367
1368 if (f_approx_equal(dot_product(params%scale, params%scale), 0._wp)) then
1369 params%scale(:) = 1._wp
1370 end if
1371
1372 if (proc_rank == 0) print *, " * Transforming model."
1373
1374 ! Get the model center before transforming the model
1375 bbox_old = f_create_bbox(model)
1376 model_center(1:3) = (bbox_old%min(1:3) + bbox_old%max(1:3))/2._wp
1377
1378 ! Compute the transform matrices for vertices and normals
1379 transform = f_create_transform_matrix(params, model_center)
1380 transform_n = f_create_transform_matrix(params)
1381
1382 call s_transform_model(model, transform, transform_n)
1383
1384 ! Recreate the bounding box after transformation
1385 bbox = f_create_bbox(model)
1386
1387 ! Show the number of vertices in the original STL model
1388 if (proc_rank == 0) print *, ' * Number of input model vertices:', 3*model%ntrs
1389
1390 ! Need the cells that form the boundary of the flat projection in 2D
1391 if (p == 0) call s_check_boundary(model, boundary_v, boundary_vertex_count, boundary_edge_count)
1392
1393 ! Show the number of edges and boundary edges in 2D STL models
1394 if (proc_rank == 0 .and. p == 0) print *, ' * Number of 2D model boundary edges:', boundary_edge_count
1395
1396 if (proc_rank == 0) then
1397 write (*, "(A, 3(2X, F20.10))") " > Model: Min:", bbox%min(1:3)
1398 write (*, "(A, 3(2X, F20.10))") " > Cen:", (bbox%min(1:3) + bbox%max(1:3))/2._wp
1399 write (*, "(A, 3(2X, F20.10))") " > Max:", bbox%max(1:3)
1400
1401 grid_mm(1,:) = (/minval(x_cc(0:m)) - 0.5_wp*dx_min, maxval(x_cc(0:m)) + 0.5_wp*dx_min/)
1402 grid_mm(2,:) = (/minval(y_cc(0:n)) - 0.5_wp*dy_min, maxval(y_cc(0:n)) + 0.5_wp*dy_min/)
1403
1404 if (p > 0) then
1405 grid_mm(3,:) = (/minval(z_cc(0:p)) - 0.5_wp*dz_min, maxval(z_cc(0:p)) + 0.5_wp*dz_min/)
1406 else
1407 grid_mm(3,:) = (/0._wp, 0._wp/)
1408 end if
1409
1410 write (*, "(A, 3(2X, F20.10))") " > Domain: Min:", grid_mm(:,1)
1411 write (*, "(A, 3(2X, F20.10))") " > Cen:", (grid_mm(:,1) + grid_mm(:,2))/2._wp
1412 write (*, "(A, 3(2X, F20.10))") " > Max:", grid_mm(:,2)
1413 end if
1414 if (proc_rank == 0) print *, " * Transforming model."
1415
1416 stl_bounding_boxes(stl_id, 1,1:3) = [bbox%min(1), (bbox%min(1) + bbox%max(1))/2._wp, bbox%max(1)]
1417 stl_bounding_boxes(stl_id, 2,1:3) = [bbox%min(2), (bbox%min(2) + bbox%max(2))/2._wp, bbox%max(2)]
1418 stl_bounding_boxes(stl_id, 3,1:3) = [bbox%min(3), (bbox%min(3) + bbox%max(3))/2._wp, bbox%max(3)]
1419
1420 models(stl_id)%model = model
1421 if (p == 0) then
1422 models(stl_id)%boundary_v = boundary_v
1423 models(stl_id)%boundary_edge_count = boundary_edge_count
1424 end if
1425 end do
1426
1427
1428# 942 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1429#if defined(MFC_OpenACC)
1430# 942 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1431!$acc update device(stl_bounding_boxes)
1432# 942 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1433#elif defined(MFC_OpenMP)
1434# 942 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1435!$omp target update to(stl_bounding_boxes)
1436# 942 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1437#endif
1438
1439# 943 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1440#if defined(MFC_OpenACC)
1441# 943 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1442!$acc update device(stl_models(1:num_stl_models))
1443# 943 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1444#elif defined(MFC_OpenMP)
1445# 943 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1446!$omp target update to(stl_models(1:num_stl_models))
1447# 943 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1448#endif
1449
1450 ! Pack and upload flat arrays for GPU (AFTER the loop)
1451 block
1452 integer :: mid, max_ntrs
1453 integer :: max_bv1, max_bv2, max_bv3
1454
1455 max_ntrs = 0
1456 max_bv1 = 0; max_bv2 = 0; max_bv3 = 0
1457
1458 do mid = 1, num_stl_models
1459 if (allocated(models(mid)%model)) then
1460 call s_pack_model_for_gpu(models(mid))
1461 max_ntrs = max(max_ntrs, models(mid)%ntrs)
1462 end if
1463 if (allocated(models(mid)%boundary_v)) then
1464 max_bv1 = max(max_bv1, size(models(mid)%boundary_v, 1))
1465 max_bv2 = max(max_bv2, size(models(mid)%boundary_v, 2))
1466 max_bv3 = max(max_bv3, size(models(mid)%boundary_v, 3))
1467 end if
1468 end do
1469
1470 if (max_ntrs > 0) then
1471#ifdef MFC_DEBUG
1472# 966 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1473 block
1474# 966 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1475 use iso_fortran_env, only: output_unit
1476# 966 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1477
1478# 966 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1479 print *, 'm_model.fpp:966: ', '@:ALLOCATE(gpu_ntrs(1:num_stl_models))'
1480# 966 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1481
1482# 966 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1483 call flush (output_unit)
1484# 966 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1485 end block
1486# 966 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1487#endif
1488# 966 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1489 allocate (gpu_ntrs(1:num_stl_models))
1490# 966 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1491
1492# 966 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1493
1494# 966 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1495#if defined(MFC_OpenACC)
1496# 966 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1497!$acc enter data create(gpu_ntrs)
1498# 966 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1499#elif defined(MFC_OpenMP)
1500# 966 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1501!$omp target enter data map(always,alloc:gpu_ntrs)
1502# 966 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1503#endif
1504#ifdef MFC_DEBUG
1505# 967 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1506 block
1507# 967 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1508 use iso_fortran_env, only: output_unit
1509# 967 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1510
1511# 967 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1512 print *, 'm_model.fpp:967: ', '@:ALLOCATE(gpu_trs_v(1:3, 1:3, 1:max_ntrs, 1:num_stl_models))'
1513# 967 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1514
1515# 967 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1516 call flush (output_unit)
1517# 967 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1518 end block
1519# 967 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1520#endif
1521# 967 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1522 allocate (gpu_trs_v(1:3, 1:3, 1:max_ntrs, 1:num_stl_models))
1523# 967 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1524
1525# 967 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1526
1527# 967 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1528#if defined(MFC_OpenACC)
1529# 967 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1530!$acc enter data create(gpu_trs_v)
1531# 967 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1532#elif defined(MFC_OpenMP)
1533# 967 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1534!$omp target enter data map(always,alloc:gpu_trs_v)
1535# 967 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1536#endif
1537#ifdef MFC_DEBUG
1538# 968 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1539 block
1540# 968 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1541 use iso_fortran_env, only: output_unit
1542# 968 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1543
1544# 968 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1545 print *, 'm_model.fpp:968: ', '@:ALLOCATE(gpu_trs_n(1:3, 1:max_ntrs, 1:num_stl_models))'
1546# 968 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1547
1548# 968 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1549 call flush (output_unit)
1550# 968 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1551 end block
1552# 968 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1553#endif
1554# 968 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1555 allocate (gpu_trs_n(1:3, 1:max_ntrs, 1:num_stl_models))
1556# 968 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1557
1558# 968 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1559
1560# 968 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1561#if defined(MFC_OpenACC)
1562# 968 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1563!$acc enter data create(gpu_trs_n)
1564# 968 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1565#elif defined(MFC_OpenMP)
1566# 968 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1567!$omp target enter data map(always,alloc:gpu_trs_n)
1568# 968 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1569#endif
1570#ifdef MFC_DEBUG
1571# 969 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1572 block
1573# 969 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1574 use iso_fortran_env, only: output_unit
1575# 969 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1576
1577# 969 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1578 print *, 'm_model.fpp:969: ', '@:ALLOCATE(gpu_boundary_edge_count(1:num_stl_models))'
1579# 969 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1580
1581# 969 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1582 call flush (output_unit)
1583# 969 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1584 end block
1585# 969 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1586#endif
1587# 969 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1588 allocate (gpu_boundary_edge_count(1:num_stl_models))
1589# 969 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1590
1591# 969 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1592
1593# 969 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1594#if defined(MFC_OpenACC)
1595# 969 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1596!$acc enter data create(gpu_boundary_edge_count)
1597# 969 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1598#elif defined(MFC_OpenMP)
1599# 969 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1600!$omp target enter data map(always,alloc:gpu_boundary_edge_count)
1601# 969 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1602#endif
1603
1604 gpu_ntrs = 0
1605 gpu_trs_v = 0._wp
1606 gpu_trs_n = 0._wp
1607 gpu_boundary_edge_count = 0
1608
1609 if (max_bv1 > 0) then
1610#ifdef MFC_DEBUG
1611# 977 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1612 block
1613# 977 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1614 use iso_fortran_env, only: output_unit
1615# 977 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1616
1617# 977 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1618 print *, 'm_model.fpp:977: ', '@:ALLOCATE(gpu_boundary_v(1:max_bv1, 1:max_bv2, 1:max_bv3, 1:num_stl_models))'
1619# 977 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1620
1621# 977 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1622 call flush (output_unit)
1623# 977 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1624 end block
1625# 977 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1626#endif
1627# 977 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1628 allocate (gpu_boundary_v(1:max_bv1, 1:max_bv2, 1:max_bv3, 1:num_stl_models))
1629# 977 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1630
1631# 977 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1632
1633# 977 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1634#if defined(MFC_OpenACC)
1635# 977 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1636!$acc enter data create(gpu_boundary_v)
1637# 977 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1638#elif defined(MFC_OpenMP)
1639# 977 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1640!$omp target enter data map(always,alloc:gpu_boundary_v)
1641# 977 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1642#endif
1643 gpu_boundary_v = 0._wp
1644 end if
1645
1646 do mid = 1, num_stl_models
1647 if (allocated(models(mid)%model)) then
1648 gpu_ntrs(mid) = models(mid)%ntrs
1649 gpu_trs_v(:,:,1:models(mid)%ntrs,mid) = models(mid)%trs_v
1650 gpu_trs_n(:,1:models(mid)%ntrs,mid) = models(mid)%trs_n
1651 gpu_boundary_edge_count(mid) = models(mid)%boundary_edge_count
1652 end if
1653 if (allocated(models(mid)%boundary_v) .and. p == 0) then
1654 gpu_boundary_v(1:size(models(mid)%boundary_v, 1),1:size(models(mid)%boundary_v, 2), &
1655 & 1:size(models(mid)%boundary_v, 3),mid) = models(mid)%boundary_v
1656 end if
1657 end do
1658
1659
1660# 994 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1661#if defined(MFC_OpenACC)
1662# 994 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1663!$acc update device(gpu_ntrs, gpu_trs_v, gpu_trs_n, gpu_boundary_edge_count)
1664# 994 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1665#elif defined(MFC_OpenMP)
1666# 994 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1667!$omp target update to(gpu_ntrs, gpu_trs_v, gpu_trs_n, gpu_boundary_edge_count)
1668# 994 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1669#endif
1670 if (allocated(gpu_boundary_v)) then
1671
1672# 996 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1673#if defined(MFC_OpenACC)
1674# 996 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1675!$acc update device(gpu_boundary_v)
1676# 996 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1677#elif defined(MFC_OpenMP)
1678# 996 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1679!$omp target update to(gpu_boundary_v)
1680# 996 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1681#endif
1682 end if
1683 end if
1684 end block
1685
1686 end subroutine s_instantiate_stl_models
1687
1688 !> Pack triangle vertices and normals from a model into flat arrays for GPU transfer.
1689 subroutine s_pack_model_for_gpu(ma)
1690
1691 type(t_model_array), intent(inout) :: ma
1692 integer :: i
1693
1694 ma%ntrs = ma%model%ntrs
1695 allocate (ma%trs_v(1:3,1:3,1:ma%ntrs))
1696 allocate (ma%trs_n(1:3,1:ma%ntrs))
1697
1698 do i = 1, ma%ntrs
1699 ma%trs_v(:,:,i) = ma%model%trs(i)%v(:,:)
1700 ma%trs_n(:,i) = ma%model%trs(i)%n(:)
1701 end do
1702
1703 end subroutine s_pack_model_for_gpu
1704
1705end module m_model
k
integer, intent(in) k
Definition m_phase_change.fpp.f90:1007
j
integer, intent(in) j
Definition m_phase_change.fpp.f90:1007
l
integer, intent(in) l
Definition m_phase_change.fpp.f90:1007
m_derived_types
Shared derived types for field data, patch geometry, bubble dynamics, and MPI I/O structures.
Definition m_derived_types.fpp.f90:335
m_helper
Utility routines for bubble model setup, coordinate transforms, array sampling, and special functions...
Definition m_helper.fpp.f90:346
m_model
Binary STL file reader and processor for immersed boundary geometry.
Definition m_model.fpp.f90:336
m_model::s_model_write
impure subroutine, public s_model_write(filepath, model)
Write a mesh to a file.
Definition m_model.fpp.f90:723
m_model::s_write_stl
impure subroutine s_write_stl(filepath, model)
Write a binary STL file.
Definition m_model.fpp.f90:647
m_model::f_read_line
impure logical function f_read_line(iunit, line)
Read the next non-blank, non-comment line from an STL or OBJ model file.
Definition m_model.fpp.f90:751
m_model::s_read_obj
impure subroutine s_read_obj(filepath, model)
Read an OBJ file.
Definition m_model.fpp.f90:559
m_model::s_distance_normals_2d
subroutine, public s_distance_normals_2d(pid, boundary_edge_count, point, normals, distance)
Determine the levelset distance and normals of 2D models by computing the exact closest point via pro...
Definition m_model.fpp.f90:1193
m_model::s_skip_ignored_lines
impure subroutine s_skip_ignored_lines(iunit, buffered_line, is_buffered)
Read the next non-comment line from a model file, using a buffered look-ahead mechanism.
Definition m_model.fpp.f90:780
m_model::gpu_trs_v
real(wp), dimension(:,:,:,:), allocatable, public gpu_trs_v
Definition m_model.fpp.f90:358
m_model::s_model_free
subroutine, public s_model_free(model)
Free the memory allocated for an STL mesh.
Definition m_model.fpp.f90:742
m_model::gpu_ntrs
integer, dimension(:), allocatable, public gpu_ntrs
GPU-friendly flat arrays for STL model data.
Definition m_model.fpp.f90:357
m_model::gpu_boundary_v
real(wp), dimension(:,:,:,:), allocatable, public gpu_boundary_v
Definition m_model.fpp.f90:360
m_model::s_register_edge
subroutine, public s_register_edge(temp_boundary_v, edge, edge_index, edge_count)
Append the edge end vertices to a temporary buffer.
Definition m_model.fpp.f90:1033
m_model::s_read_stl
impure subroutine s_read_stl(filepath, model)
Read an STL file.
Definition m_model.fpp.f90:531
m_model::f_model_read
impure type(t_model) function, public f_model_read(filepath)
Read a mesh from a file.
Definition m_model.fpp.f90:628
m_model::gpu_trs_n
real(wp), dimension(:,:,:), allocatable, public gpu_trs_n
Definition m_model.fpp.f90:359
m_model::s_check_boundary
subroutine, public s_check_boundary(model, boundary_v, boundary_vertex_count, boundary_edge_count)
Check and label edges shared by two or more triangle facets of the 2D STL model.
Definition m_model.fpp.f90:879
m_model::s_distance_normals_3d
subroutine, public s_distance_normals_3d(ntrs, pid, point, normals, distance)
Determine the levelset distance and normals of 3D models by computing the exact closest point via pro...
Definition m_model.fpp.f90:1049
m_model::models
type(t_model_array), dimension(:), allocatable, target, public models
STL/OBJ models for IB markers and levelset.
Definition m_model.fpp.f90:356
m_model::s_write_obj
impure subroutine s_write_obj(filepath, model)
Write an OBJ file.
Definition m_model.fpp.f90:692
m_model::f_model_is_inside
real(wp) function, public f_model_is_inside(ntrs, pid, point)
Determine if a point is inside a surface using the generalized winding number (Jacobson et al....
Definition m_model.fpp.f90:804
m_model::stl_bounding_boxes
real(wp), dimension(:,:,:), allocatable, public stl_bounding_boxes
Definition m_model.fpp.f90:362
m_model::s_read_stl_ascii
impure subroutine s_read_stl_ascii(filepath, model)
Read an ASCII STL file.
Definition m_model.fpp.f90:423
m_model::s_read_stl_binary
impure subroutine s_read_stl_binary(filepath, model)
Read a binary STL file.
Definition m_model.fpp.f90:379
m_model::s_instantiate_stl_models
subroutine, public s_instantiate_stl_models()
Load, transform, and register STL/OBJ immersed-boundary models onto the simulation grid.
Definition m_model.fpp.f90:1275
m_model::s_pack_model_for_gpu
subroutine, public s_pack_model_for_gpu(ma)
Pack triangle vertices and normals from a model into flat arrays for GPU transfer.
Definition m_model.fpp.f90:1690
m_model::gpu_boundary_edge_count
integer, dimension(:), allocatable, public gpu_boundary_edge_count
Definition m_model.fpp.f90:361
m_mpi_proxy
Broadcasts user inputs and decomposes the domain across MPI ranks for pre-processing.
Definition m_mpi_proxy.fpp.f90:7
m_derived_types::ic_model_parameters
Defines parameters for a Model Patch.
Definition m_derived_types.fpp.f90:559
m_derived_types::t_bbox
Definition m_derived_types.fpp.f90:573
m_derived_types::t_model_array
Definition m_derived_types.fpp.f90:583
m_derived_types::t_model
Definition m_derived_types.fpp.f90:578