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