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