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