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