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# 207 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
108
109# 232 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
110
111# 243 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
112
113# 245 "/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# 283 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
117
118# 293 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
119
120# 303 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
121
122# 312 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
123
124# 329 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
125
126# 339 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
127
128# 346 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
129
130# 352 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
131
132# 358 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
133
134# 364 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
135
136# 370 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
137
138# 376 "/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# 192 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
208
209# 213 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
210
211# 241 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
212
213# 256 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
214
215# 266 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
216
217# 275 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
218
219# 291 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
220
221# 301 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
222
223# 308 "/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# 21 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
228
229# 37 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
230
231# 50 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
232
233# 76 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
234
235# 91 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
236
237# 102 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
238
239# 115 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
240
241# 143 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
242
243# 154 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
244
245# 165 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
246
247# 176 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
248
249# 187 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
250
251# 198 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
252
253# 208 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
254
255# 214 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
256
257# 220 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
258
259# 226 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
260
261# 232 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
262
263# 234 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
264# 235 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
265! New line at end of file is required for FYPP
266# 2 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp" 2
267
268# 14 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
269
270! Caution:
271! This macro requires the use of a binding script to set CUDA_VISIBLE_DEVICES, such that we have one GPU device per MPI rank.
272! That's because for both cudaMemAdvise (preferred location) and cudaMemPrefetchAsync we use location = device_id = 0.
273! For an example see misc/nvidia_uvm/bind.sh.
274# 63 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
275
276# 81 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
277
278# 88 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
279
280# 111 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
281
282# 127 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
283
284# 153 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
285
286# 159 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
287
288# 167 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
289! New line at end of file is required for FYPP
290# 7 "/home/runner/work/MFC/MFC/src/common/m_model.fpp" 2
291
292!> @brief Binary STL file reader and processor for immersed boundary geometry
293module m_model
294
295 use m_helper
296 use m_mpi_proxy
297 use m_derived_types
298
299 use iso_c_binding, only: c_char, c_int32_t, c_int16_t, c_float
300
301 implicit none
302
303 private
304
305 public :: f_model_read, s_model_write, s_model_free, f_model_is_inside
306
307 ! Subroutines for STL immersed boundaries
308 public :: f_check_boundary, f_register_edge, f_check_interpolation_2d, &
309 f_check_interpolation_3d, f_interpolate_2d, f_interpolate_3d, &
310 f_interpolated_distance, f_normals, f_distance, f_distance_normals_3d, f_tri_area
311
312contains
313
314 !> This procedure reads a binary STL file.
315 !! @param filepath Path to the STL file.
316 !! @param model The binary of the STL file.
317 impure subroutine s_read_stl_binary(filepath, model)
318
319 character(LEN=*), intent(in) :: filepath
320 type(t_model), intent(out) :: model
321
322 integer :: i, iunit, iostat
323
324 character(kind=c_char, len=80) :: header
325 integer(kind=c_int32_t) :: ntriangles
326
327 real(kind=c_float) :: normal(3), v(3, 3)
328 integer(kind=c_int16_t) :: attribute
329
330 open (newunit=iunit, file=filepath, action='READ', &
331 form='UNFORMATTED', status='OLD', iostat=iostat, &
332 access='STREAM')
333
334 if (iostat /= 0) then
335 print *, "Error: could not open Binary STL file ", filepath
336
337 call s_mpi_abort()
338 end if
339
340 read (iunit, iostat=iostat) header, ntriangles
341
342 if (iostat /= 0) then
343 print *, "Error: could not read header from Binary STL file ", filepath
344
345 call s_mpi_abort()
346 end if
347
348 model%ntrs = ntriangles
349
350 allocate (model%trs(model%ntrs))
351
352 do i = 1, model%ntrs
353 read (iunit) normal(:), v(1, :), v(2, :), v(3, :), attribute
354
355 model%trs(i)%v = v
356 model%trs(i)%n = normal
357 end do
358
359 close (iunit)
360
361 end subroutine s_read_stl_binary
362
363 !> This procedure reads an ASCII STL file.
364 !! @param filepath Path to the STL file.
365 !! @param model the STL file.
366 impure subroutine s_read_stl_ascii(filepath, model)
367 character(LEN=*), intent(in) :: filepath
368 type(t_model), intent(out) :: model
369
370 integer :: i, j, iunit, iostat
371 character(80) :: line, buffered_line
372 logical :: is_buffered
373
374 is_buffered = .false.
375
376 open (newunit=iunit, file=filepath, action='READ', &
377 form='FORMATTED', status='OLD', iostat=iostat, &
378 access='STREAM')
379
380 if (iostat /= 0) then
381 print *, "Error: could not open ASCII STL file ", filepath
382 call s_mpi_abort()
383 end if
384
385 model%ntrs = 0
386 do
387 if (is_buffered) then
388 line = buffered_line
389 is_buffered = .false.
390 else
391 if (.not. f_read_line(iunit, line)) exit
392 end if
393
394 if (line(1:6) == "facet ") then
395 model%ntrs = model%ntrs + 1
396 end if
397 end do
398
399 allocate (model%trs(model%ntrs))
400
401 rewind(iunit)
402
403 i = 1
404 do
405 if (is_buffered) then
406 line = buffered_line
407 is_buffered = .false.
408 else
409 if (.not. f_read_line(iunit, line)) exit
410 end if
411
412 if (line(1:5) == "solid") cycle
413 if (line(1:8) == "endsolid") exit
414
415 if (line(1:12) /= "facet normal") then
416 print *, "Error: expected facet normal in STL file ", filepath
417 call s_mpi_abort()
418 end if
419
420 call s_skip_ignored_lines(iunit, buffered_line, is_buffered)
421 read (line(13:), *) model%trs(i)%n
422
423 call s_skip_ignored_lines(iunit, buffered_line, is_buffered)
424 if (is_buffered) then
425 line = buffered_line
426 is_buffered = .false.
427 else
428 read (iunit, '(A)') line
429 end if
430
431 do j = 1, 3
432 if (is_buffered) then
433 line = buffered_line
434 is_buffered = .false.
435 else
436 if (.not. f_read_line(iunit, line)) exit
437 end if
438
439 if (line(1:6) /= "vertex") then
440 print *, "Error: expected vertex in STL file ", filepath
441 call s_mpi_abort()
442 end if
443
444 call s_skip_ignored_lines(iunit, buffered_line, is_buffered)
445 read (line(7:), *) model%trs(i)%v(j, :)
446 end do
447
448 if (is_buffered) then
449 line = buffered_line
450 is_buffered = .false.
451 else
452 if (.not. f_read_line(iunit, line)) exit
453 end if
454
455 if (is_buffered) then
456 line = buffered_line
457 is_buffered = .false.
458 else
459 if (.not. f_read_line(iunit, line)) exit
460 end if
461
462 if (line(1:8) /= "endfacet") then
463 print *, "Error: expected endfacet in STL file ", filepath
464 call s_mpi_abort()
465 end if
466
467 i = i + 1
468 end do
469 end subroutine s_read_stl_ascii
470
471 !> This procedure reads an STL file.
472 !! @param filepath Path to the STL file.
473 !! @param model the STL file.
474 impure subroutine s_read_stl(filepath, model)
475
476 character(LEN=*), intent(in) :: filepath
477 type(t_model), intent(out) :: model
478
479 integer :: iunit, iostat
480
481 character(80) :: line
482
483 open (newunit=iunit, file=filepath, action='READ', &
484 form='FORMATTED', status='OLD', iostat=iostat, &
485 access='STREAM')
486
487 if (iostat /= 0) then
488 print *, "Error: could not open STL file ", filepath
489
490 call s_mpi_abort()
491 end if
492
493 read (iunit, '(A)') line
494
495 close (iunit)
496
497 if (line(1:5) == "solid") then
498 call s_read_stl_ascii(filepath, model)
499 else
500 call s_read_stl_binary(filepath, model)
501 end if
502
503 end subroutine s_read_stl
504
505 !> This procedure reads an OBJ file.
506 !! @param filepath Path to the odj file.
507 !! @param model The obj file.
508 impure subroutine s_read_obj(filepath, model)
509
510 character(LEN=*), intent(in) :: filepath
511 type(t_model), intent(out) :: model
512
513 integer :: i, j, k, l, iunit, iostat, nvertices
514
515 real(wp), dimension(1:3), allocatable :: vertices(:, :)
516
517 character(80) :: line
518
519 open (newunit=iunit, file=filepath, action='READ', &
520 form='FORMATTED', status='OLD', iostat=iostat, &
521 access='STREAM')
522
523 if (iostat /= 0) then
524 print *, "Error: could not open model file ", filepath
525
526 call s_mpi_abort()
527 end if
528
529 nvertices = 0
530 model%ntrs = 0
531 do
532 if (.not. f_read_line(iunit, line)) exit
533
534 select case (line(1:2))
535 case ("v ")
536 nvertices = nvertices + 1
537 case ("f ")
538 model%ntrs = model%ntrs + 1
539 end select
540 end do
541
542 rewind(iunit)
543
544 allocate (vertices(nvertices, 1:3))
545 allocate (model%trs(model%ntrs))
546
547 i = 1
548 j = 1
549
550 do
551 if (.not. f_read_line(iunit, line)) exit
552
553 select case (line(1:2))
554 case ("g ")
555 case ("vn")
556 case ("vt")
557 case ("l ")
558 case ("v ")
559 read (line(3:), *) vertices(i, :)
560 i = i + 1
561 case ("f ")
562 read (line(3:), *) k, l, j
563 model%trs(j)%v(1, :) = vertices(k, :)
564 model%trs(j)%v(2, :) = vertices(l, :)
565 model%trs(j)%v(3, :) = vertices(j, :)
566 j = j + 1
567 case default
568 print *, "Error: unknown line type in OBJ file ", filepath
569 print *, "Line: ", line
570
571 call s_mpi_abort()
572 end select
573 end do
574
575 deallocate (vertices)
576
577 close (iunit)
578
579 end subroutine s_read_obj
580
581 !> This procedure reads a mesh from a file.
582 !! @param filepath Path to the file to read.
583 !! @return The model read from the file.
584 impure function f_model_read(filepath) result(model)
585
586 character(LEN=*), intent(in) :: filepath
587
588 type(t_model) :: model
589
590 select case (filepath(len(trim(filepath)) - 3:len(trim(filepath))))
591 case (".stl")
592 call s_read_stl(filepath, model)
593 case (".obj")
594 call s_read_obj(filepath, model)
595 case default
596 print *, "Error: unknown model file format for file ", filepath
597
598 call s_mpi_abort()
599 end select
600
601 end function f_model_read
602
603 !> This procedure writes a binary STL file.
604 !! @param filepath Path to the STL file.
605 !! @param model STL to write
606 impure subroutine s_write_stl(filepath, model)
607
608 character(LEN=*), intent(in) :: filepath
609 type(t_model), intent(in) :: model
610
611 integer :: i, j, iunit, iostat
612
613 character(kind=c_char, len=80), parameter :: header = "Model file written by MFC."
614 integer(kind=c_int32_t) :: ntriangles
615 real(wp) :: normal(3), v(3)
616 integer(kind=c_int16_t) :: attribute
617
618 open (newunit=iunit, file=filepath, action='WRITE', &
619 form='UNFORMATTED', iostat=iostat, access='STREAM')
620
621 if (iostat /= 0) then
622 print *, "Error: could not open STL file ", filepath
623
624 call s_mpi_abort()
625 end if
626
627 ntriangles = model%ntrs
628 write (iunit, iostat=iostat) header, ntriangles
629
630 if (iostat /= 0) then
631 print *, "Error: could not write header to STL file ", filepath
632
633 call s_mpi_abort()
634 end if
635
636 do i = 1, model%ntrs
637 normal = model%trs(i)%n
638 write (iunit) normal
639
640 do j = 1, 3
641 v = model%trs(i)%v(j, :)
642 write (iunit) v(:)
643 end do
644
645 attribute = 0
646 write (iunit) attribute
647 end do
648
649 close (iunit)
650
651 end subroutine s_write_stl
652
653 !> This procedure writes an OBJ file.
654 !! @param filepath Path to the obj file.
655 !! @param model obj to write.
656 impure subroutine s_write_obj(filepath, model)
657
658 character(LEN=*), intent(in) :: filepath
659 type(t_model), intent(in) :: model
660
661 integer :: iunit, iostat
662
663 integer :: i, j
664
665 open (newunit=iunit, file=filepath, action='WRITE', &
666 form='FORMATTED', iostat=iostat, access='STREAM')
667
668 if (iostat /= 0) then
669 print *, "Error: could not open OBJ file ", filepath
670
671 call s_mpi_abort()
672 end if
673
674 write (iunit, '(A)') "# Model file written by MFC."
675
676 do i = 1, model%ntrs
677 do j = 1, 3
678 write (iunit, '(A, " ", (f30.20), " ", (f30.20), " ", (f30.20))') &
679 "v", model%trs(i)%v(j, 1), model%trs(i)%v(j, 2), model%trs(i)%v(j, 3)
680 end do
681
682 write (iunit, '(A, " ", I0, " ", I0, " ", I0)') &
683 "f", i*3 - 2, i*3 - 1, i*3
684 end do
685
686 close (iunit)
687
688 end subroutine s_write_obj
689
690 !> This procedure writes a binary STL file.
691 !! @param filepath Path to the file to write.
692 !! @param model Model to write.
693 impure subroutine s_model_write(filepath, model)
694
695 character(LEN=*), intent(in) :: filepath
696 type(t_model), intent(in) :: model
697
698 select case (filepath(len(trim(filepath)) - 3:len(trim(filepath))))
699 case (".stl")
700 call s_write_stl(filepath, model)
701 case (".obj")
702 call s_write_obj(filepath, model)
703 case default
704 print *, "Error: unknown model file format for file ", filepath
705
706 call s_mpi_abort()
707 end select
708
709 end subroutine s_model_write
710
711 !> This procedure frees the memory allocated for an STL mesh.
712 subroutine s_model_free(model)
713
714 type(t_model), intent(inout) :: model
715
716 deallocate (model%trs)
717
718 end subroutine s_model_free
719
720 impure function f_read_line(iunit, line) result(bIsLine)
721
722 integer, intent(in) :: iunit
723 character(80), intent(out) :: line
724
725 logical :: bisline
726 integer :: iostat
727
728 bisline = .true.
729
730 do
731 read (iunit, '(A)', iostat=iostat) line
732
733 if (iostat < 0) then
734 bisline = .false.
735 exit
736 end if
737
738 line = adjustl(trim(line))
739
740 if (len(trim(line)) == 0) cycle
741 if (line(1:5) == "solid") cycle
742 if (line(1:1) == "#") cycle
743
744 exit
745 end do
746
747 end function f_read_line
748
749 !> @brief Reads the next non-comment line from a model file, using a buffered look-ahead mechanism.
750 impure subroutine s_skip_ignored_lines(iunit, buffered_line, is_buffered)
751 integer, intent(in) :: iunit
752 character(80), intent(inout) :: buffered_line
753 logical, intent(inout) :: is_buffered
754
755 character(80) :: line
756
757 if (is_buffered) then
758 line = buffered_line
759 is_buffered = .false.
760 else
761 if (.not. f_read_line(iunit, line)) return
762 end if
763
764 buffered_line = line
765 is_buffered = .true.
766 end subroutine s_skip_ignored_lines
767
768 !> This procedure, recursively, finds whether a point is inside an octree.
769 !! @param model Model to search in.
770 !! @param point Point to test.
771 !! @param spacing Space around the point to search in (grid spacing).
772 !! @param spc Number of samples per cell.
773 !! @return True if the point is inside the octree, false otherwise.
774 impure function f_model_is_inside(model, point, spacing, spc) result(fraction)
775
776 ! $:GPU_ROUTINE(parallelism='[seq]')
777
778 type(t_model), intent(in) :: model
779 real(wp), dimension(1:3), intent(in) :: point
780 real(wp), dimension(1:3), intent(in) :: spacing
781 integer, intent(in) :: spc
782
783 real(wp) :: fraction
784
785 type(t_ray) :: ray
786 integer :: i, j, ninorout, nhits
787
788 real(wp), dimension(1:spc, 1:3) :: ray_origins, ray_dirs
789
790 ! TODO :: The random number generation prohibits GPU compute due to the subroutine not being able to be called in kernels
791 ! This should be swapped out with something that allows GPU compute. I recommend the fibonacci sphere:
792 ! do i = 1, spc
793 ! phi = acos(1.0 - 2.0*(i-1.0)/(spc-1.0))
794 ! theta = pi * (1.0 + sqrt(5.0)) * (i-1.0)
795 ! ray_dirs(i,:) = [cos(theta)*sin(phi), sin(theta)*sin(phi), cos(phi)]
796 ! ray_origins(i,:) = point
797 ! end do
798
799 do i = 1, spc
800 call random_number(ray_origins(i, :))
801 ray_origins(i, :) = point + (ray_origins(i, :) - 0.5_wp)*spacing(:)
802
803 call random_number(ray_dirs(i, :))
804 ray_dirs(i, :) = ray_dirs(i, :) - 0.5_wp
805 ray_dirs(i, :) = ray_dirs(i, :)/sqrt(sum(ray_dirs(i, :)*ray_dirs(i, :)))
806 end do
807
808 ninorout = 0
809 do i = 1, spc
810 ray%o = ray_origins(i, :)
811 ray%d = ray_dirs(i, :)
812
813 nhits = 0
814 do j = 1, model%ntrs
815 if (f_intersects_triangle(ray, model%trs(j))) then
816 nhits = nhits + 1
817 end if
818 end do
819
820 ninorout = ninorout + mod(nhits, 2)
821 end do
822
823 fraction = real(ninorout)/real(spc)
824
825 end function f_model_is_inside
826
827 ! From https://www.scratchapixel.com/lessons/3e-basic-rendering/ray-tracing-rendering-a-triangle/ray-triangle-intersection-geometric-solution.html
828 !> This procedure checks if a ray intersects a triangle.
829 !! @param ray Ray.
830 !! @param triangle Triangle.
831 !! @return True if the ray intersects the triangle, false otherwise.
832 elemental function f_intersects_triangle(ray, triangle) result(intersects)
833
834 type(t_ray), intent(in) :: ray
835 type(t_triangle), intent(in) :: triangle
836
837 logical :: intersects
838
839 real(wp) :: n(3), p(3), c(3), edge(3), vp(3)
840 real(wp) :: area2, d, t, ndotraydirection
841
842 intersects = .false.
843
844 n = triangle%n
845 area2 = sqrt(sum(n(:)*n(:)))
846
847 ndotraydirection = sum(n(:)*ray%d(:))
848
849 if (abs(ndotraydirection) < 0.0000001_wp) then
850 return
851 end if
852
853 d = -sum(n(:)*triangle%v(1, :))
854 t = -(sum(n(:)*ray%o(:)) + d)/ndotraydirection
855
856 if (t < 0) then
857 return
858 end if
859
860 p = ray%o + t*ray%d
861
862 edge = triangle%v(2, :) - triangle%v(1, :)
863 vp = p - triangle%v(1, :)
864 c = f_cross(edge, vp)
865 if (sum(n(:)*c(:)) < 0) then
866 return
867 end if
868
869 edge = triangle%v(3, :) - triangle%v(2, :)
870 vp = p - triangle%v(2, :)
871 c = f_cross(edge, vp)
872 if (sum(n(:)*c(:)) < 0) then
873 return
874 end if
875
876 edge = triangle%v(1, :) - triangle%v(3, :)
877 vp = p - triangle%v(3, :)
878 c = f_cross(edge, vp)
879 if (sum(n(:)*c(:)) < 0) then
880 return
881 end if
882
883 intersects = .true.
884
885 end function f_intersects_triangle
886
887 !> This procedure checks and labels edges shared by two or more triangles facets of the 2D STL model.
888 !! @param model Model to search in.
889 !! @param boundary_vertex_count Output total boundary vertex count
890 subroutine f_check_boundary(model, boundary_v, boundary_vertex_count, boundary_edge_count)
891
892 type(t_model), intent(in) :: model
893 real(wp), allocatable, intent(out), dimension(:, :, :) :: boundary_v !< Output boundary vertices/normals
894 integer, intent(out) :: boundary_vertex_count, boundary_edge_count !< Output boundary vertex/edge count
895
896 integer :: i, j !< Model index iterator
897 integer :: edge_count, edge_index, store_index !< Boundary edge index iterator
898 real(wp), dimension(1:2, 1:2) :: edge !< Edge end points buffer
899 real(wp), dimension(1:2) :: boundary_edge !< Boundary edge end points buffer
900 real(wp), dimension(1:(3*model%ntrs), 1:2, 1:2) :: temp_boundary_v !< Temporary boundary vertex buffer
901 integer, dimension(1:(3*model%ntrs)) :: edge_occurrence !< The manifoldness of the edges
902 real(wp) :: edgetan, initial, v_norm, xnormal, ynormal !< The manifoldness of the edges
903
904 ! Total number of edges in 2D STL
905 edge_count = 3*model%ntrs
906
907 ! Initialize edge_occurrence array to zero
908 edge_occurrence = 0
909 edge_index = 0
910
911 ! Collect all edges of all triangles and store them
912 do i = 1, model%ntrs
913 ! First edge (v1, v2)
914 edge(1, 1:2) = model%trs(i)%v(1, 1:2)
915 edge(2, 1:2) = model%trs(i)%v(2, 1:2)
916 call f_register_edge(temp_boundary_v, edge, edge_index, edge_count)
917
918 ! Second edge (v2, v3)
919 edge(1, 1:2) = model%trs(i)%v(2, 1:2)
920 edge(2, 1:2) = model%trs(i)%v(3, 1:2)
921 call f_register_edge(temp_boundary_v, edge, edge_index, edge_count)
922
923 ! Third edge (v3, v1)
924 edge(1, 1:2) = model%trs(i)%v(3, 1:2)
925 edge(2, 1:2) = model%trs(i)%v(1, 1:2)
926 call f_register_edge(temp_boundary_v, edge, edge_index, edge_count)
927 end do
928
929 ! Check all edges and count repeated edges
930 do i = 1, edge_count
931 do j = 1, edge_count
932 if (i /= j) then
933 if (((abs(temp_boundary_v(i, 1, 1) - temp_boundary_v(j, 1, 1)) < threshold_edge_zero) .and. &
934 (abs(temp_boundary_v(i, 1, 2) - temp_boundary_v(j, 1, 2)) < threshold_edge_zero) .and. &
935 (abs(temp_boundary_v(i, 2, 1) - temp_boundary_v(j, 2, 1)) < threshold_edge_zero) .and. &
936 (abs(temp_boundary_v(i, 2, 2) - temp_boundary_v(j, 2, 2)) < threshold_edge_zero)) .or. &
937 ((abs(temp_boundary_v(i, 1, 1) - temp_boundary_v(j, 2, 1)) < threshold_edge_zero) .and. &
938 (abs(temp_boundary_v(i, 1, 2) - temp_boundary_v(j, 2, 2)) < threshold_edge_zero) .and. &
939 (abs(temp_boundary_v(i, 2, 1) - temp_boundary_v(j, 1, 1)) < threshold_edge_zero) .and. &
940 (abs(temp_boundary_v(i, 2, 2) - temp_boundary_v(j, 1, 2)) < threshold_edge_zero))) then
941
942 edge_occurrence(i) = edge_occurrence(i) + 1
943 end if
944 end if
945 end do
946 end do
947
948 ! Count the number of boundary vertices/edges
949 boundary_vertex_count = 0
950 boundary_edge_count = 0
951
952 do i = 1, edge_count
953 if (edge_occurrence(i) == 0) then
954 boundary_vertex_count = boundary_vertex_count + 2
955 boundary_edge_count = boundary_edge_count + 1
956 end if
957 end do
958
959 ! Allocate the boundary_v array based on the number of boundary edges
960 allocate (boundary_v(boundary_edge_count, 1:3, 1:2))
961
962 ! Store boundary vertices
963 store_index = 0
964 do i = 1, edge_count
965 if (edge_occurrence(i) == 0) then
966 store_index = store_index + 1
967 boundary_v(store_index, 1, 1:2) = temp_boundary_v(i, 1, 1:2)
968 boundary_v(store_index, 2, 1:2) = temp_boundary_v(i, 2, 1:2)
969 end if
970 end do
971
972 ! Find/store the normal vector of the boundary edges
973 do i = 1, boundary_edge_count
974 boundary_edge(1) = boundary_v(i, 2, 1) - boundary_v(i, 1, 1)
975 boundary_edge(2) = boundary_v(i, 2, 2) - boundary_v(i, 1, 2)
976 edgetan = boundary_edge(1)/sign(max(sgm_eps, abs(boundary_edge(2))), boundary_edge(2))
977
978 if (abs(boundary_edge(2)) < threshold_vector_zero) then
979 if (edgetan > 0._wp) then
980 ynormal = -1
981 xnormal = 0._wp
982 else
983 ynormal = 1
984 xnormal = 0._wp
985 end if
986 else
987 initial = boundary_edge(2)
988 ynormal = -edgetan*initial
989 xnormal = initial
990 end if
991
992 v_norm = sqrt(xnormal**2 + ynormal**2)
993 boundary_v(i, 3, 1) = xnormal/v_norm
994 boundary_v(i, 3, 2) = ynormal/v_norm
995 end do
996
997 end subroutine f_check_boundary
998
999 !> This procedure appends the edge end vertices to a temporary buffer.
1000 subroutine f_register_edge(temp_boundary_v, edge, edge_index, edge_count)
1001
1002 integer, intent(inout) :: edge_index !< Edge index iterator
1003 integer, intent(inout) :: edge_count !< Total number of edges
1004 real(wp), intent(in), dimension(1:2, 1:2) :: edge !< Edges end points to be registered
1005 real(wp), dimension(1:edge_count, 1:2, 1:2), intent(inout) :: temp_boundary_v !< Temporary edge end vertex buffer
1006
1007 ! Increment edge index and store the edge
1008 edge_index = edge_index + 1
1009 temp_boundary_v(edge_index, 1, 1:2) = edge(1, 1:2)
1010 temp_boundary_v(edge_index, 2, 1:2) = edge(2, 1:2)
1011
1012 end subroutine f_register_edge
1013
1014 !> This procedure check if interpolates is needed for 2D models.
1015 !! @param boundary_v Temporary edge end vertex buffer
1016 !! @param spacing Dimensions of the current levelset cell
1017 subroutine f_check_interpolation_2d(boundary_v, boundary_edge_count, spacing, interpolate)
1018
1019 logical, intent(inout) :: interpolate !< Logical indicator of interpolation
1020 integer, intent(in) :: boundary_edge_count !< Number of boundary edges
1021 real(wp), intent(in), dimension(1:boundary_edge_count, 1:3, 1:2) :: boundary_v
1022 real(wp), dimension(1:3), intent(in) :: spacing
1023
1024 real(wp) :: l1, cell_width !< Length of each boundary edge and cell width
1025 integer :: j !< Boundary edge index iterator
1026
1027 cell_width = minval(spacing(1:2))
1028 interpolate = .false.
1029
1030 do j = 1, boundary_edge_count
1031 l1 = sqrt((boundary_v(j, 2, 1) - boundary_v(j, 1, 1))**2 + &
1032 (boundary_v(j, 2, 2) - boundary_v(j, 1, 2))**2)
1033
1034 if ((l1 > cell_width)) then
1035 interpolate = .true.
1036 else
1037 interpolate = .false.
1038 end if
1039 end do
1040
1041 end subroutine f_check_interpolation_2d
1042
1043 !> This procedure check if interpolates is needed for 3D models.
1044 !! @param model Model to search in.
1045 !! @param spacing Dimensions of the current levelset cell
1046 !! @param interpolate Logical output
1047 subroutine f_check_interpolation_3d(model, spacing, interpolate)
1048
1049 logical, intent(inout) :: interpolate
1050 type(t_model), intent(in) :: model
1051 real(wp), dimension(1:3), intent(in) :: spacing
1052 real(wp), dimension(1:3) :: edge_l
1053 real(wp) :: cell_width
1054 real(wp), dimension(1:3, 1:3) :: tri_v
1055 integer :: i, j !< Loop iterator
1056
1057 cell_width = minval(spacing)
1058 interpolate = .false.
1059
1060 do i = 1, model%ntrs
1061 do j = 1, 3
1062 tri_v(1, j) = model%trs(i)%v(1, j)
1063 tri_v(2, j) = model%trs(i)%v(2, j)
1064 tri_v(3, j) = model%trs(i)%v(3, j)
1065 end do
1066
1067 edge_l(1) = sqrt((tri_v(1, 2) - tri_v(1, 1))**2 + &
1068 (tri_v(2, 2) - tri_v(2, 1))**2 + &
1069 (tri_v(3, 2) - tri_v(3, 1))**2)
1070 edge_l(2) = sqrt((tri_v(1, 3) - tri_v(1, 2))**2 + &
1071 (tri_v(2, 3) - tri_v(2, 2))**2 + &
1072 (tri_v(3, 3) - tri_v(3, 2))**2)
1073 edge_l(3) = sqrt((tri_v(1, 1) - tri_v(1, 3))**2 + &
1074 (tri_v(2, 1) - tri_v(2, 3))**2 + &
1075 (tri_v(3, 1) - tri_v(3, 3))**2)
1076
1077 if ((edge_l(1) > cell_width) .or. &
1078 (edge_l(2) > cell_width) .or. &
1079 (edge_l(3) > cell_width)) then
1080 interpolate = .true.
1081 else
1082 interpolate = .false.
1083 end if
1084 end do
1085
1086 end subroutine f_check_interpolation_3d
1087
1088 !> This procedure interpolates 2D models.
1089 !! @param boundary_v Group of all the boundary vertices of the 2D model without interpolation
1090 !! @param boundary_edge_count Output total number of boundary edges
1091 !! @param spacing Dimensions of the current levelset cell
1092 !! @param interpolated_boundary_v Output all the boundary vertices of the interpolated 2D model
1093 !! @param total_vertices Total number of vertices after interpolation
1094 subroutine f_interpolate_2d(boundary_v, boundary_edge_count, spacing, interpolated_boundary_v, total_vertices)
1095
1096 real(wp), intent(in), dimension(:, :, :) :: boundary_v
1097 real(wp), dimension(1:3), intent(in) :: spacing
1098 real(wp), allocatable, intent(inout), dimension(:, :) :: interpolated_boundary_v
1099
1100 integer, intent(inout) :: total_vertices, boundary_edge_count
1101 integer :: num_segments
1102 integer :: i, j
1103
1104 real(wp) :: edge_length, cell_width
1105 real(wp), dimension(1:2) :: edge_x, edge_y, edge_del
1106
1107 ! Get the number of boundary edges
1108 cell_width = minval(spacing(1:2))
1109 num_segments = 0
1110
1111 ! First pass: Calculate the total number of vertices including interpolated ones
1112 total_vertices = 1
1113 do i = 1, boundary_edge_count
1114 ! Get the coordinates of the two ends of the current edge
1115 edge_x(1) = boundary_v(i, 1, 1)
1116 edge_y(1) = boundary_v(i, 1, 2)
1117 edge_x(2) = boundary_v(i, 2, 1)
1118 edge_y(2) = boundary_v(i, 2, 2)
1119
1120 ! Compute the length of the edge
1121 edge_length = sqrt((edge_x(2) - edge_x(1))**2 + &
1122 (edge_y(2) - edge_y(1))**2)
1123
1124 ! Determine the number of segments
1125 if (edge_length > cell_width) then
1126 num_segments = ifactor_2d*ceiling(edge_length/cell_width)
1127 else
1128 num_segments = 1
1129 end if
1130
1131 ! Each edge contributes num_segments vertices
1132 total_vertices = total_vertices + num_segments
1133 end do
1134
1135 ! Allocate memory for the new boundary vertices array
1136 allocate (interpolated_boundary_v(1:total_vertices, 1:3))
1137
1138 ! Fill the new boundary vertices array with original and interpolated vertices
1139 total_vertices = 1
1140 do i = 1, boundary_edge_count
1141 ! Get the coordinates of the two ends of the current edge
1142 edge_x(1) = boundary_v(i, 1, 1)
1143 edge_y(1) = boundary_v(i, 1, 2)
1144 edge_x(2) = boundary_v(i, 2, 1)
1145 edge_y(2) = boundary_v(i, 2, 2)
1146
1147 ! Compute the length of the edge
1148 edge_length = sqrt((edge_x(2) - edge_x(1))**2 + &
1149 (edge_y(2) - edge_y(1))**2)
1150
1151 ! Determine the number of segments and interpolation step
1152 if (edge_length > cell_width) then
1153 num_segments = ifactor_2d*ceiling(edge_length/cell_width)
1154 edge_del(1) = (edge_x(2) - edge_x(1))/num_segments
1155 edge_del(2) = (edge_y(2) - edge_y(1))/num_segments
1156 else
1157 num_segments = 1
1158 edge_del(1) = 0._wp
1159 edge_del(2) = 0._wp
1160 end if
1161
1162 interpolated_boundary_v(1, 1) = edge_x(1)
1163 interpolated_boundary_v(1, 2) = edge_y(1)
1164 interpolated_boundary_v(1, 3) = 0._wp
1165
1166 ! Add original and interpolated vertices to the output array
1167 do j = 1, num_segments - 1
1168 total_vertices = total_vertices + 1
1169 interpolated_boundary_v(total_vertices, 1) = edge_x(1) + j*edge_del(1)
1170 interpolated_boundary_v(total_vertices, 2) = edge_y(1) + j*edge_del(2)
1171 end do
1172
1173 ! Add the last vertex of the edge
1174 if (num_segments > 0) then
1175 total_vertices = total_vertices + 1
1176 interpolated_boundary_v(total_vertices, 1) = edge_x(2)
1177 interpolated_boundary_v(total_vertices, 2) = edge_y(2)
1178 end if
1179 end do
1180
1181 end subroutine f_interpolate_2d
1182
1183 !> This procedure interpolates 3D models.
1184 !! @param model Model to search in.
1185 !! @param spacing Dimensions of the current levelset cell
1186 !! @param interpolated_boundary_v Output all the boundary vertices of the interpolated 3D model
1187 !! @param total_vertices Total number of vertices after interpolation
1188 impure subroutine f_interpolate_3d(model, spacing, interpolated_boundary_v, total_vertices)
1189 real(wp), dimension(1:3), intent(in) :: spacing
1190 type(t_model), intent(in) :: model
1191 real(wp), allocatable, intent(inout), dimension(:, :) :: interpolated_boundary_v
1192 integer, intent(out) :: total_vertices
1193
1194 integer :: i, j, k, num_triangles, num_segments, num_inner_vertices
1195 real(wp), dimension(1:3, 1:3) :: tri
1196 real(wp), dimension(1:3) :: edge_del, cell_area
1197 real(wp), dimension(1:3) :: bary_coord !< Barycentric coordinates
1198 real(wp) :: edge_length, cell_width, cell_area_min, tri_area
1199
1200 ! Number of triangles in the model
1201 num_triangles = model%ntrs
1202 cell_width = minval(spacing)
1203
1204 ! Find the minimum surface area
1205 cell_area(1) = spacing(1)*spacing(2)
1206 cell_area(2) = spacing(1)*spacing(3)
1207 cell_area(3) = spacing(2)*spacing(3)
1208 cell_area_min = minval(cell_area)
1209 num_inner_vertices = 0
1210
1211 ! Calculate the total number of vertices including interpolated ones
1212 total_vertices = 0
1213 do i = 1, num_triangles
1214 do j = 1, 3
1215 ! Get the coordinates of the two vertices of the current edge
1216 tri(1, 1) = model%trs(i)%v(j, 1)
1217 tri(1, 2) = model%trs(i)%v(j, 2)
1218 tri(1, 3) = model%trs(i)%v(j, 3)
1219 ! Next vertex in the triangle (cyclic)
1220 tri(2, 1) = model%trs(i)%v(mod(j, 3) + 1, 1)
1221 tri(2, 2) = model%trs(i)%v(mod(j, 3) + 1, 2)
1222 tri(2, 3) = model%trs(i)%v(mod(j, 3) + 1, 3)
1223
1224 ! Compute the length of the edge
1225 edge_length = sqrt((tri(2, 1) - tri(1, 1))**2 + &
1226 (tri(2, 2) - tri(1, 2))**2 + &
1227 (tri(2, 3) - tri(1, 3))**2)
1228
1229 ! Determine the number of segments
1230 if (edge_length > cell_width) then
1231 num_segments = ifactor_3d*ceiling(edge_length/cell_width)
1232 else
1233 num_segments = 1
1234 end if
1235
1236 ! Each edge contributes num_segments vertices
1237 total_vertices = total_vertices + num_segments + 1
1238 end do
1239
1240 ! Add vertices inside the triangle
1241 do k = 1, 3
1242 tri(k, 1) = model%trs(i)%v(k, 1)
1243 tri(k, 2) = model%trs(i)%v(k, 2)
1244 tri(k, 3) = model%trs(i)%v(k, 3)
1245 end do
1246 call f_tri_area(tri, tri_area)
1247
1248 if (tri_area > threshold_bary*cell_area_min) then
1249 num_inner_vertices = ifactor_bary_3d*ceiling(tri_area/cell_area_min)
1250 total_vertices = total_vertices + num_inner_vertices
1251 end if
1252 end do
1253
1254 ! Allocate memory for the new boundary vertices array
1255 allocate (interpolated_boundary_v(1:total_vertices, 1:3))
1256
1257 ! Fill the new boundary vertices array with original and interpolated vertices
1258 total_vertices = 0
1259 do i = 1, num_triangles
1260 ! Loop through the 3 edges of each triangle
1261 do j = 1, 3
1262 ! Get the coordinates of the two vertices of the current edge
1263 tri(1, 1) = model%trs(i)%v(j, 1)
1264 tri(1, 2) = model%trs(i)%v(j, 2)
1265 tri(1, 3) = model%trs(i)%v(j, 3)
1266 ! Next vertex in the triangle (cyclic)
1267 tri(2, 1) = model%trs(i)%v(mod(j, 3) + 1, 1)
1268 tri(2, 2) = model%trs(i)%v(mod(j, 3) + 1, 2)
1269 tri(2, 3) = model%trs(i)%v(mod(j, 3) + 1, 3)
1270
1271 ! Compute the length of the edge
1272 edge_length = sqrt((tri(2, 1) - tri(1, 1))**2 + &
1273 (tri(2, 2) - tri(1, 2))**2 + &
1274 (tri(2, 3) - tri(1, 3))**2)
1275
1276 ! Determine the number of segments and interpolation step
1277 if (edge_length > cell_width) then
1278 num_segments = ifactor_3d*ceiling(edge_length/cell_width)
1279 edge_del(1) = (tri(2, 1) - tri(1, 1))/num_segments
1280 edge_del(2) = (tri(2, 2) - tri(1, 2))/num_segments
1281 edge_del(3) = (tri(2, 3) - tri(1, 3))/num_segments
1282 else
1283 num_segments = 1
1284 edge_del = 0._wp
1285 end if
1286
1287 ! Add original and interpolated vertices to the output array
1288 do k = 0, num_segments - 1
1289 total_vertices = total_vertices + 1
1290 interpolated_boundary_v(total_vertices, 1) = tri(1, 1) + k*edge_del(1)
1291 interpolated_boundary_v(total_vertices, 2) = tri(1, 2) + k*edge_del(2)
1292 interpolated_boundary_v(total_vertices, 3) = tri(1, 3) + k*edge_del(3)
1293 end do
1294
1295 ! Add the last vertex of the edge
1296 total_vertices = total_vertices + 1
1297 interpolated_boundary_v(total_vertices, 1) = tri(2, 1)
1298 interpolated_boundary_v(total_vertices, 2) = tri(2, 2)
1299 interpolated_boundary_v(total_vertices, 3) = tri(2, 3)
1300 end do
1301
1302 ! Interpolate verties that are not on edges
1303 do k = 1, 3
1304 tri(k, 1) = model%trs(i)%v(k, 1)
1305 tri(k, 2) = model%trs(i)%v(k, 2)
1306 tri(k, 3) = model%trs(i)%v(k, 3)
1307 end do
1308 call f_tri_area(tri, tri_area)
1309
1310 if (tri_area > threshold_bary*cell_area_min) then
1311 num_inner_vertices = ifactor_bary_3d*ceiling(tri_area/cell_area_min)
1312 !Use barycentric coordinates for randomly distributed points
1313 do k = 1, num_inner_vertices
1314 call random_number(bary_coord(1))
1315 call random_number(bary_coord(2))
1316
1317 if ((bary_coord(1) + bary_coord(2)) >= 1._wp) then
1318 bary_coord(1) = 1._wp - bary_coord(1)
1319 bary_coord(2) = 1._wp - bary_coord(2)
1320 end if
1321 bary_coord(3) = 1._wp - bary_coord(1) - bary_coord(2)
1322
1323 total_vertices = total_vertices + 1
1324 interpolated_boundary_v(total_vertices, 1) = dot_product(bary_coord, tri(1:3, 1))
1325 interpolated_boundary_v(total_vertices, 2) = dot_product(bary_coord, tri(1:3, 2))
1326 interpolated_boundary_v(total_vertices, 3) = dot_product(bary_coord, tri(1:3, 3))
1327 end do
1328 end if
1329 end do
1330
1331 end subroutine f_interpolate_3d
1332
1333 !> This procedure determines the levelset distance and normals of the 3D models without interpolation.
1334 !! @param model Model to search in.
1335 !! @param point The cell centers of the current level cell
1336 !! @param normals The output levelset normals
1337 !! @param distance The output levelset distance
1338 subroutine f_distance_normals_3d(model, point, normals, distance)
1339
1340
1341# 1056 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1342#if MFC_OpenACC
1343# 1056 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1344!$acc routine seq
1345# 1056 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1346#elif MFC_OpenMP
1347# 1056 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1348
1349# 1056 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1350
1351# 1056 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1352!$omp declare target device_type(any)
1353# 1056 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1354#endif
1355
1356 type(t_model), intent(IN) :: model
1357 real(wp), dimension(1:3), intent(in) :: point
1358 real(wp), dimension(1:3), intent(out) :: normals
1359 real(wp), intent(out) :: distance
1360
1361 real(wp), dimension(1:3, 1:3) :: tri
1362 real(wp) :: dist_min, dist_t_min
1363 real(wp) :: dist_min_normal, dist_buffer_normal
1364 real(wp), dimension(1:3) :: midp !< Centers of the triangle facets
1365 real(wp), dimension(1:3) :: dist_buffer !< Distance between the cell center and the vertices
1366 integer :: i, j, tri_idx !< Iterator
1367
1368 dist_min = 1.e12_wp
1369 dist_min_normal = 1.e12_wp
1370 distance = 0._wp
1371
1372 tri_idx = 0
1373 do i = 1, model%ntrs
1374 do j = 1, 3
1375 tri(j, 1) = model%trs(i)%v(j, 1)
1376 tri(j, 2) = model%trs(i)%v(j, 2)
1377 tri(j, 3) = model%trs(i)%v(j, 3)
1378 dist_buffer(j) = sqrt((point(1) - tri(j, 1))**2 + &
1379 (point(2) - tri(j, 2))**2 + &
1380 (point(3) - tri(j, 3))**2)
1381 end do
1382
1383 ! Get the surface center of each triangle facet
1384 do j = 1, 3
1385 midp(j) = (tri(1, j) + tri(2, j) + tri(3, j))/3
1386 end do
1387
1388 dist_t_min = minval(dist_buffer(1:3))
1389 dist_buffer_normal = sqrt((point(1) - midp(1))**2 + &
1390 (point(2) - midp(2))**2 + &
1391 (point(3) - midp(3))**2)
1392
1393 if (dist_t_min < dist_min) then
1394 dist_min = dist_t_min
1395 end if
1396
1397 if (dist_buffer_normal < dist_min_normal) then
1398 dist_min_normal = dist_buffer_normal
1399 tri_idx = i
1400 end if
1401 end do
1402
1403 normals(1) = model%trs(tri_idx)%n(1)
1404 normals(2) = model%trs(tri_idx)%n(2)
1405 normals(3) = model%trs(tri_idx)%n(3)
1406 distance = dist_min
1407
1408 end subroutine f_distance_normals_3d
1409
1410 !> This procedure determines the levelset distance of 2D models without interpolation.
1411 !! @param boundary_v Group of all the boundary vertices of the 2D model without interpolation
1412 !! @param boundary_edge_count Output the total number of boundary edges
1413 !! @param point The cell centers of the current levelset cell
1414 !! @return Distance which the levelset distance without interpolation
1415 function f_distance(boundary_v, boundary_edge_count, point) result(distance)
1416
1417
1418# 1119 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1419#if MFC_OpenACC
1420# 1119 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1421!$acc routine seq
1422# 1119 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1423#elif MFC_OpenMP
1424# 1119 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1425
1426# 1119 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1427
1428# 1119 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1429!$omp declare target device_type(any)
1430# 1119 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1431#endif
1432
1433 integer, intent(in) :: boundary_edge_count
1434 real(wp), intent(in), dimension(:, :, :) :: boundary_v
1435 ! real(wp), intent(in), dimension(1:boundary_edge_count, 1:3, 1:2) :: boundary_v
1436 real(wp), dimension(1:3), intent(in) :: point
1437
1438 integer :: i
1439 real(wp) :: dist_buffer1, dist_buffer2
1440 real(wp), dimension(1:boundary_edge_count) :: dist_buffer
1441 real(wp) :: distance
1442
1443 distance = 0._wp
1444 do i = 1, boundary_edge_count
1445 dist_buffer1 = sqrt((point(1) - boundary_v(i, 1, 1))**2 + &
1446 & (point(2) - boundary_v(i, 1, 2))**2)
1447
1448 dist_buffer2 = sqrt((point(1) - boundary_v(i, 2, 1))**2 + &
1449 & (point(2) - boundary_v(i, 2, 2))**2)
1450
1451 dist_buffer(i) = minval((/dist_buffer1, dist_buffer2/))
1452 end do
1453
1454 distance = minval(dist_buffer)
1455
1456 end function f_distance
1457
1458 !> This procedure determines the levelset normals of 2D models without interpolation.
1459 !! @param boundary_v Group of all the boundary vertices of the 2D model without interpolation
1460 !! @param boundary_edge_count Output the total number of boundary edges
1461 !! @param point The cell centers of the current levelset cell
1462 !! @param normals Output levelset normals without interpolation
1463 subroutine f_normals(boundary_v, boundary_edge_count, point, normals)
1464
1465
1466# 1153 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1467#if MFC_OpenACC
1468# 1153 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1469!$acc routine seq
1470# 1153 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1471#elif MFC_OpenMP
1472# 1153 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1473
1474# 1153 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1475
1476# 1153 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1477!$omp declare target device_type(any)
1478# 1153 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1479#endif
1480
1481 integer, intent(in) :: boundary_edge_count
1482 real(wp), intent(in), dimension(:, :, :) :: boundary_v
1483 real(wp), dimension(1:3), intent(in) :: point
1484 real(wp), dimension(1:3), intent(out) :: normals
1485
1486 integer :: i, idx_buffer
1487 real(wp) :: dist_min, dist_buffer
1488 real(wp) :: midp(1:3)
1489
1490 dist_buffer = 0._wp
1491 dist_min = initial_distance_buffer
1492 idx_buffer = 0
1493
1494 do i = 1, boundary_edge_count
1495 midp(1) = (boundary_v(i, 2, 1) + boundary_v(i, 1, 1))/2
1496 midp(2) = (boundary_v(i, 2, 2) + boundary_v(i, 1, 2))/2
1497 midp(3) = 0._wp
1498
1499 dist_buffer = sqrt((point(1) - midp(1))**2 + &
1500 & (point(2) - midp(2))**2)
1501
1502 if (dist_buffer < dist_min) then
1503 dist_min = dist_buffer
1504 idx_buffer = i
1505 end if
1506 end do
1507
1508 normals(1) = boundary_v(idx_buffer, 3, 1)
1509 normals(2) = boundary_v(idx_buffer, 3, 2)
1510 normals(3) = 0._wp
1511
1512 end subroutine f_normals
1513
1514 !> This procedure calculates the barycentric facet area
1515 subroutine f_tri_area(tri, tri_area)
1516
1517 real(wp), dimension(1:3, 1:3), intent(in) :: tri
1518 real(wp), intent(out) :: tri_area
1519 real(wp), dimension(1:3) :: ab, ac, cross
1520 integer :: i !< Loop iterator
1521
1522 do i = 1, 3
1523 ab(i) = tri(2, i) - tri(1, i)
1524 ac(i) = tri(3, i) - tri(1, i)
1525 end do
1526
1527 cross(1) = ab(2)*ac(3) - ab(3)*ac(2)
1528 cross(2) = ab(3)*ac(1) - ab(1)*ac(3)
1529 cross(3) = ab(1)*ac(2) - ab(2)*ac(1)
1530 tri_area = 0.5_wp*sqrt(cross(1)**2 + cross(2)**2 + cross(3)**2)
1531
1532 end subroutine f_tri_area
1533
1534 !> This procedure determines the levelset of interpolated 2D models.
1535 !! @param interpolated_boundary_v Group of all the boundary vertices of the interpolated 2D model
1536 !! @param total_vertices Total number of vertices after interpolation
1537 !! @param point The cell centers of the current levelset cell
1538 !! @return Distance which the levelset distance without interpolation
1539 function f_interpolated_distance(interpolated_boundary_v, total_vertices, point) result(distance)
1540
1541
1542# 1215 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1543#if MFC_OpenACC
1544# 1215 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1545!$acc routine seq
1546# 1215 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1547#elif MFC_OpenMP
1548# 1215 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1549
1550# 1215 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1551
1552# 1215 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1553!$omp declare target device_type(any)
1554# 1215 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1555#endif
1556
1557 integer, intent(in) :: total_vertices
1558 real(wp), intent(in), dimension(:, :) :: interpolated_boundary_v
1559 real(wp), dimension(1:3), intent(in) :: point
1560
1561 integer :: i !< Loop iterator
1562 real(wp) :: dist_buffer, min_dist
1563 real(wp) :: distance
1564
1565 distance = initial_distance_buffer
1566 dist_buffer = initial_distance_buffer
1567 min_dist = initial_distance_buffer
1568
1569 do i = 1, total_vertices
1570 dist_buffer = sqrt((point(1) - interpolated_boundary_v(i, 1))**2 + &
1571 (point(2) - interpolated_boundary_v(i, 2))**2 + &
1572 (point(3) - interpolated_boundary_v(i, 3))**2)
1573
1574 if (min_dist > dist_buffer) then
1575 min_dist = dist_buffer
1576 end if
1577 end do
1578
1579 distance = min_dist
1580
1581 end function f_interpolated_distance
1582
1583end module m_model
k
integer, intent(in) k
Definition m_phase_change.fpp.f90:1121
j
integer, intent(in) j
Definition m_phase_change.fpp.f90:1121
l
integer, intent(in) l
Definition m_phase_change.fpp.f90:1121
m_derived_types
Shared derived types for field data, patch geometry, bubble dynamics, and MPI I/O structures.
Definition m_derived_types.fpp.f90:292
m_helper
Utility routines for bubble model setup, coordinate transforms, array sampling, and special functions...
Definition m_helper.fpp.f90:303
m_helper::f_cross
pure real(wp) function, dimension(3), public f_cross(a, b)
This procedure computes the cross product of two vectors.
Definition m_helper.fpp.f90:803
m_model
Binary STL file reader and processor for immersed boundary geometry.
Definition m_model.fpp.f90:293
m_model::s_model_write
impure subroutine, public s_model_write(filepath, model)
This procedure writes a binary STL file.
Definition m_model.fpp.f90:694
m_model::s_write_stl
impure subroutine s_write_stl(filepath, model)
This procedure writes a binary STL file.
Definition m_model.fpp.f90:607
m_model::f_check_interpolation_2d
subroutine, public f_check_interpolation_2d(boundary_v, boundary_edge_count, spacing, interpolate)
This procedure check if interpolates is needed for 2D models.
Definition m_model.fpp.f90:1018
m_model::f_check_interpolation_3d
subroutine, public f_check_interpolation_3d(model, spacing, interpolate)
This procedure check if interpolates is needed for 3D models.
Definition m_model.fpp.f90:1048
m_model::f_distance_normals_3d
subroutine, public f_distance_normals_3d(model, point, normals, distance)
This procedure determines the levelset distance and normals of the 3D models without interpolation.
Definition m_model.fpp.f90:1339
m_model::f_read_line
impure logical function f_read_line(iunit, line)
Definition m_model.fpp.f90:721
m_model::f_interpolated_distance
real(wp) function, public f_interpolated_distance(interpolated_boundary_v, total_vertices, point)
This procedure determines the levelset of interpolated 2D models.
Definition m_model.fpp.f90:1540
m_model::s_read_obj
impure subroutine s_read_obj(filepath, model)
This procedure reads an OBJ file.
Definition m_model.fpp.f90:509
m_model::f_model_is_inside
impure real(wp) function, public f_model_is_inside(model, point, spacing, spc)
This procedure, recursively, finds whether a point is inside an octree.
Definition m_model.fpp.f90:775
m_model::s_skip_ignored_lines
impure subroutine s_skip_ignored_lines(iunit, buffered_line, is_buffered)
Reads the next non-comment line from a model file, using a buffered look-ahead mechanism.
Definition m_model.fpp.f90:751
m_model::f_check_boundary
subroutine, public f_check_boundary(model, boundary_v, boundary_vertex_count, boundary_edge_count)
This procedure checks and labels edges shared by two or more triangles facets of the 2D STL model.
Definition m_model.fpp.f90:891
m_model::s_model_free
subroutine, public s_model_free(model)
This procedure frees the memory allocated for an STL mesh.
Definition m_model.fpp.f90:713
m_model::f_interpolate_2d
subroutine, public f_interpolate_2d(boundary_v, boundary_edge_count, spacing, interpolated_boundary_v, total_vertices)
This procedure interpolates 2D models.
Definition m_model.fpp.f90:1095
m_model::s_read_stl
impure subroutine s_read_stl(filepath, model)
This procedure reads an STL file.
Definition m_model.fpp.f90:475
m_model::f_normals
subroutine, public f_normals(boundary_v, boundary_edge_count, point, normals)
This procedure determines the levelset normals of 2D models without interpolation.
Definition m_model.fpp.f90:1464
m_model::f_model_read
impure type(t_model) function, public f_model_read(filepath)
This procedure reads a mesh from a file.
Definition m_model.fpp.f90:585
m_model::s_write_obj
impure subroutine s_write_obj(filepath, model)
This procedure writes an OBJ file.
Definition m_model.fpp.f90:657
m_model::f_distance
real(wp) function, public f_distance(boundary_v, boundary_edge_count, point)
This procedure determines the levelset distance of 2D models without interpolation.
Definition m_model.fpp.f90:1416
m_model::f_register_edge
subroutine, public f_register_edge(temp_boundary_v, edge, edge_index, edge_count)
This procedure appends the edge end vertices to a temporary buffer.
Definition m_model.fpp.f90:1001
m_model::f_tri_area
subroutine, public f_tri_area(tri, tri_area)
This procedure calculates the barycentric facet area.
Definition m_model.fpp.f90:1516
m_model::f_interpolate_3d
impure subroutine, public f_interpolate_3d(model, spacing, interpolated_boundary_v, total_vertices)
This procedure interpolates 3D models.
Definition m_model.fpp.f90:1189
m_model::s_read_stl_ascii
impure subroutine s_read_stl_ascii(filepath, model)
This procedure reads an ASCII STL file.
Definition m_model.fpp.f90:367
m_model::s_read_stl_binary
impure subroutine s_read_stl_binary(filepath, model)
This procedure reads a binary STL file.
Definition m_model.fpp.f90:318
m_model::f_intersects_triangle
elemental logical function f_intersects_triangle(ray, triangle)
This procedure checks if a ray intersects a triangle.
Definition m_model.fpp.f90:833
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::t_model
Definition m_derived_types.fpp.f90:463
m_derived_types::t_ray
Definition m_derived_types.fpp.f90:453
m_derived_types::t_triangle
Definition m_derived_types.fpp.f90:448