m_model.fpp.f90

Go to the documentation of this file.

# 1 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
!>
!! @file
!! @author Henry Le Berre <hberre3@gatech.edu>
!! @brief  Contains module m_model
 
# 1 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp" 1
# 1 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp" 1
# 1 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp" 1
# 2 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 3 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 4 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 5 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 6 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 8 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 9 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 10 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 17 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 46 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 58 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 68 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 98 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 110 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 120 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
! New line at end of file is required for FYPP
# 2 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp" 2
# 1 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp" 1
# 1 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp" 1
# 2 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 3 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 4 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 5 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 6 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 8 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 9 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 10 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 17 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 46 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 58 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 68 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 98 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 110 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 120 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
! New line at end of file is required for FYPP
# 2 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp" 2
 
# 4 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
# 5 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
# 6 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
# 7 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
# 8 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 20 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 43 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 48 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 53 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 58 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 63 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 68 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 76 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 81 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 86 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 91 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 96 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 101 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 106 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 111 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 116 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 121 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 151 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 192 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 207 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 232 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 243 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 245 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
# 255 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 283 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 293 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 303 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 312 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 329 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 339 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 346 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 352 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 358 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 364 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 370 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 376 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
! New line at end of file is required for FYPP
# 3 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp" 2
# 1 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp" 1
# 1 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp" 1
# 2 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 3 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 4 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 5 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 6 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 8 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 9 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 10 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 17 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 46 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 58 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 68 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 98 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 110 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 120 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
! New line at end of file is required for FYPP
# 2 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp" 2
 
# 7 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 17 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 22 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 27 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 32 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 37 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 42 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 47 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 52 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 57 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 62 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 73 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 78 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 83 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 88 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 103 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 131 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 160 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 175 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 192 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 213 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 241 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 256 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 266 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 275 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 291 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 301 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 308 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
! New line at end of file is required for FYPP
# 4 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp" 2
 
# 21 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
# 37 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
# 50 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
# 104 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
# 119 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
# 130 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
# 143 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
# 171 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
# 182 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
# 193 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
# 204 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
# 214 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
# 225 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
# 236 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
# 246 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
# 252 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
# 258 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
# 264 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
# 270 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
# 272 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
# 273 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
! New line at end of file is required for FYPP
# 2 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp" 2
 
# 14 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
 
! Caution:
! This macro requires the use of a binding script to set CUDA_VISIBLE_DEVICES, such that we have one GPU device per MPI rank.
! That's because for both cudaMemAdvise (preferred location) and cudaMemPrefetchAsync we use location = device_id = 0.
! For an example see misc/nvidia_uvm/bind.sh.
# 63 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
 
# 81 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
 
# 88 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
 
# 111 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
 
# 127 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
 
# 153 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
 
# 159 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
 
# 167 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
! New line at end of file is required for FYPP
# 7 "/home/runner/work/MFC/MFC/src/common/m_model.fpp" 2
 
!> @brief Binary STL file reader and processor for immersed boundary geometry
module m_model
 
    use m_helper
    use m_mpi_proxy
    use m_derived_types
 
    use iso_c_binding, only: c_char, c_int32_t, c_int16_t, c_float
 
    implicit none
 
    private
 
    public :: f_model_read, s_model_write, s_model_free, f_model_is_inside, models, gpu_ntrs, &
              gpu_trs_v, gpu_trs_n, gpu_boundary_v, gpu_boundary_edge_count, &
              gpu_total_vertices, stl_bounding_boxes
 
    ! Subroutines for STL immersed boundaries
    public :: s_check_boundary, s_register_edge, f_model_is_inside_flat, &
              s_distance_normals_3d, s_distance_normals_2d, s_pack_model_for_gpu
 
#ifdef MFC_SIMULATION
    public :: s_instantiate_stl_models
#endif
 
    !! array of STL models that can be allocated and then used in IB marker and levelset compute
    type(t_model_array), allocatable, target :: models(:)
    !! GPU-friendly flat arrays for STL model data
    integer, allocatable :: gpu_ntrs(:)
    real(wp), allocatable, dimension(:, :, :, :) :: gpu_trs_v
    real(wp), allocatable, dimension(:, :, :) :: gpu_trs_n
    real(wp), allocatable, dimension(:, :, :, :) :: gpu_boundary_v
    integer, allocatable :: gpu_boundary_edge_count(:)
    integer, allocatable :: gpu_total_vertices(:)
    real(wp), allocatable :: stl_bounding_boxes(:, :, :)
 
# 43 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
#if defined(MFC_OpenACC)
# 43 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
!$acc declare create(gpu_ntrs, gpu_trs_v, gpu_trs_n, gpu_boundary_v, gpu_boundary_edge_count, gpu_total_vertices) 
# 43 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
#elif defined(MFC_OpenMP)
# 43 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
!$omp declare target (gpu_ntrs, gpu_trs_v, gpu_trs_n, gpu_boundary_v, gpu_boundary_edge_count, gpu_total_vertices) 
# 43 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
#endif
 
contains
 
    !> This procedure reads a binary STL file.
    !! @param filepath Path to the STL file.
    !! @param model The binary of the STL file.
    impure subroutine s_read_stl_binary(filepath, model)
 
        character(LEN=*), intent(in) :: filepath
        type(t_model), intent(out) :: model
 
        integer :: i, iunit, iostat
 
        character(kind=c_char, len=80) :: header
        integer(kind=c_int32_t) :: ntriangles
 
        real(kind=c_float) :: normal(3), v(3, 3), v_norm
        integer(kind=c_int16_t) :: attribute
 
        open (newunit=iunit, file=filepath, action='READ', &
              form='UNFORMATTED', status='OLD', iostat=iostat, &
              access='STREAM')
 
        if (iostat /= 0) then
            print *, "Error: could not open Binary STL file ", filepath
 
            call s_mpi_abort()
        end if
 
        read (iunit, iostat=iostat) header, ntriangles
 
        if (iostat /= 0) then
            print *, "Error: could not read header from Binary STL file ", filepath
 
            call s_mpi_abort()
        end if
 
        model%ntrs = ntriangles
 
        allocate (model%trs(model%ntrs))
 
        do i = 1, model%ntrs
            read (iunit) normal(:), v(1, :), v(2, :), v(3, :), attribute
 
            model%trs(i)%v = v
            model%trs(i)%n = normal
            v_norm = sqrt(normal(1)**2 + normal(2)**2 + normal(3)**2)
            if (v_norm > 0._wp) model%trs(i)%n = normal/v_norm
        end do
 
        close (iunit)
 
    end subroutine s_read_stl_binary
 
    !> This procedure reads an ASCII STL file.
    !! @param filepath Path to the STL file.
    !! @param model the STL file.
    impure subroutine s_read_stl_ascii(filepath, model)
        character(LEN=*), intent(in) :: filepath
        type(t_model), intent(out) :: model
 
        integer :: i, j, iunit, iostat
        character(80) :: line, buffered_line
        logical :: is_buffered
        real(wp) :: normal(3), v_norm
 
        is_buffered = .false.
 
        open (newunit=iunit, file=filepath, action='READ', &
              form='FORMATTED', status='OLD', iostat=iostat, &
              access='STREAM')
 
        if (iostat /= 0) then
            print *, "Error: could not open ASCII STL file ", filepath
            call s_mpi_abort()
        end if
 
        model%ntrs = 0
        do
            if (is_buffered) then
                line = buffered_line
                is_buffered = .false.
            else
                if (.not. f_read_line(iunit, line)) exit
            end if
 
            if (line(1:6) == "facet ") then
                model%ntrs = model%ntrs + 1
            end if
        end do
 
        allocate (model%trs(model%ntrs))
 
        rewind(iunit)
 
        i = 1
        do
            if (is_buffered) then
                line = buffered_line
                is_buffered = .false.
            else
                if (.not. f_read_line(iunit, line)) exit
            end if
 
            if (line(1:5) == "solid") cycle
            if (line(1:8) == "endsolid") exit
 
            if (line(1:12) /= "facet normal") then
                print *, "Error: expected facet normal in STL file ", filepath
                call s_mpi_abort()
            end if
 
            call s_skip_ignored_lines(iunit, buffered_line, is_buffered)
            read (line(13:), *) normal
            v_norm = sqrt(normal(1)**2 + normal(2)**2 + normal(3)**2)
            if (v_norm > 0._wp) model%trs(i)%n = normal/v_norm
 
            call s_skip_ignored_lines(iunit, buffered_line, is_buffered)
            if (is_buffered) then
                line = buffered_line
                is_buffered = .false.
            else
                read (iunit, '(A)') line
            end if
 
            do j = 1, 3
                if (is_buffered) then
                    line = buffered_line
                    is_buffered = .false.
                else
                    if (.not. f_read_line(iunit, line)) exit
                end if
 
                if (line(1:6) /= "vertex") then
                    print *, "Error: expected vertex in STL file ", filepath
                    call s_mpi_abort()
                end if
 
                call s_skip_ignored_lines(iunit, buffered_line, is_buffered)
                read (line(7:), *) model%trs(i)%v(j, :)
            end do
 
            if (is_buffered) then
                line = buffered_line
                is_buffered = .false.
            else
                if (.not. f_read_line(iunit, line)) exit
            end if
 
            if (is_buffered) then
                line = buffered_line
                is_buffered = .false.
            else
                if (.not. f_read_line(iunit, line)) exit
            end if
 
            if (line(1:8) /= "endfacet") then
                print *, "Error: expected endfacet in STL file ", filepath
                call s_mpi_abort()
            end if
 
            i = i + 1
        end do
    end subroutine s_read_stl_ascii
 
    !> This procedure reads an STL file.
    !! @param filepath Path to the STL file.
    !! @param model the STL file.
    impure subroutine s_read_stl(filepath, model)
 
        character(LEN=*), intent(in) :: filepath
        type(t_model), intent(out) :: model
 
        integer :: iunit, iostat
 
        character(80) :: line
 
        open (newunit=iunit, file=filepath, action='READ', &
              form='FORMATTED', status='OLD', iostat=iostat, &
              access='STREAM')
 
        if (iostat /= 0) then
            print *, "Error: could not open STL file ", filepath
 
            call s_mpi_abort()
        end if
 
        read (iunit, '(A)') line
 
        close (iunit)
 
        if (line(1:5) == "solid") then
            call s_read_stl_ascii(filepath, model)
        else
            call s_read_stl_binary(filepath, model)
        end if
 
    end subroutine s_read_stl
 
    !> This procedure reads an OBJ file.
    !! @param filepath Path to the odj file.
    !! @param model The obj file.
    impure subroutine s_read_obj(filepath, model)
 
        character(LEN=*), intent(in) :: filepath
        type(t_model), intent(out) :: model
 
        integer :: i, j, k, l, iunit, iostat, nvertices
 
        real(wp), dimension(1:3), allocatable :: vertices(:, :)
 
        character(80) :: line
 
        open (newunit=iunit, file=filepath, action='READ', &
              form='FORMATTED', status='OLD', iostat=iostat, &
              access='STREAM')
 
        if (iostat /= 0) then
            print *, "Error: could not open model file ", filepath
 
            call s_mpi_abort()
        end if
 
        nvertices = 0
        model%ntrs = 0
        do
            if (.not. f_read_line(iunit, line)) exit
 
            select case (line(1:2))
            case ("v ")
                nvertices = nvertices + 1
            case ("f ")
                model%ntrs = model%ntrs + 1
            end select
        end do
 
        rewind(iunit)
 
        allocate (vertices(nvertices, 1:3))
        allocate (model%trs(model%ntrs))
 
        i = 1
        j = 1
 
        do
            if (.not. f_read_line(iunit, line)) exit
 
            select case (line(1:2))
            case ("g ")
            case ("vn")
            case ("vt")
            case ("l ")
            case ("v ")
                read (line(3:), *) vertices(i, :)
                i = i + 1
            case ("f ")
                read (line(3:), *) k, l, j
                model%trs(j)%v(1, :) = vertices(k, :)
                model%trs(j)%v(2, :) = vertices(l, :)
                model%trs(j)%v(3, :) = vertices(j, :)
                j = j + 1
            case default
                print *, "Error: unknown line type in OBJ file ", filepath
                print *, "Line: ", line
 
                call s_mpi_abort()
            end select
        end do
 
        deallocate (vertices)
 
        close (iunit)
 
    end subroutine s_read_obj
 
    !> This procedure reads a mesh from a file.
    !! @param filepath Path to the file to read.
    !! @return The model read from the file.
    impure function f_model_read(filepath) result(model)
 
        character(LEN=*), intent(in) :: filepath
 
        type(t_model) :: model
 
        select case (filepath(len(trim(filepath)) - 3:len(trim(filepath))))
        case (".stl")
            call s_read_stl(filepath, model)
        case (".obj")
            call s_read_obj(filepath, model)
        case default
            print *, "Error: unknown model file format for file ", filepath
 
            call s_mpi_abort()
        end select
 
    end function f_model_read
 
    !> This procedure writes a binary STL file.
    !! @param filepath Path to the STL file.
    !! @param model STL to write
    impure subroutine s_write_stl(filepath, model)
 
        character(LEN=*), intent(in) :: filepath
        type(t_model), intent(in) :: model
 
        integer :: i, j, iunit, iostat
 
        character(kind=c_char, len=80), parameter :: header = "Model file written by MFC."
        integer(kind=c_int32_t) :: ntriangles
        real(wp) :: normal(3), v(3)
        integer(kind=c_int16_t) :: attribute
 
        open (newunit=iunit, file=filepath, action='WRITE', &
              form='UNFORMATTED', iostat=iostat, access='STREAM')
 
        if (iostat /= 0) then
            print *, "Error: could not open STL file ", filepath
 
            call s_mpi_abort()
        end if
 
        ntriangles = model%ntrs
        write (iunit, iostat=iostat) header, ntriangles
 
        if (iostat /= 0) then
            print *, "Error: could not write header to STL file ", filepath
 
            call s_mpi_abort()
        end if
 
        do i = 1, model%ntrs
            normal = model%trs(i)%n
            write (iunit) normal
 
            do j = 1, 3
                v = model%trs(i)%v(j, :)
                write (iunit) v(:)
            end do
 
            attribute = 0
            write (iunit) attribute
        end do
 
        close (iunit)
 
    end subroutine s_write_stl
 
    !> This procedure writes an OBJ file.
    !! @param filepath Path to the obj file.
    !! @param model obj to write.
    impure subroutine s_write_obj(filepath, model)
 
        character(LEN=*), intent(in) :: filepath
        type(t_model), intent(in) :: model
 
        integer :: iunit, iostat
 
        integer :: i, j
 
        open (newunit=iunit, file=filepath, action='WRITE', &
              form='FORMATTED', iostat=iostat, access='STREAM')
 
        if (iostat /= 0) then
            print *, "Error: could not open OBJ file ", filepath
 
            call s_mpi_abort()
        end if
 
        write (iunit, '(A)') "# Model file written by MFC."
 
        do i = 1, model%ntrs
            do j = 1, 3
                write (iunit, '(A, " ", (f30.20), " ", (f30.20), " ", (f30.20))') &
                    "v", model%trs(i)%v(j, 1), model%trs(i)%v(j, 2), model%trs(i)%v(j, 3)
            end do
 
            write (iunit, '(A, " ", I0, " ", I0, " ", I0)') &
                "f", i*3 - 2, i*3 - 1, i*3
        end do
 
        close (iunit)
 
    end subroutine s_write_obj
 
    !> This procedure writes a binary STL file.
    !! @param filepath  Path to the file to write.
    !! @param model Model to write.
    impure subroutine s_model_write(filepath, model)
 
        character(LEN=*), intent(in) :: filepath
        type(t_model), intent(in) :: model
 
        select case (filepath(len(trim(filepath)) - 3:len(trim(filepath))))
        case (".stl")
            call s_write_stl(filepath, model)
        case (".obj")
            call s_write_obj(filepath, model)
        case default
            print *, "Error: unknown model file format for file ", filepath
 
            call s_mpi_abort()
        end select
 
    end subroutine s_model_write
 
    !> This procedure frees the memory allocated for an STL mesh.
    subroutine s_model_free(model)
 
        type(t_model), intent(inout) :: model
 
        deallocate (model%trs)
 
    end subroutine s_model_free
 
    impure function f_read_line(iunit, line) result(bIsLine)
 
        integer, intent(in) :: iunit
        character(80), intent(out) :: line
 
        logical :: bisline
        integer :: iostat
 
        bisline = .true.
 
        do
            read (iunit, '(A)', iostat=iostat) line
 
            if (iostat < 0) then
                bisline = .false.
                exit
            end if
 
            line = adjustl(trim(line))
 
            if (len(trim(line)) == 0) cycle
            if (line(1:5) == "solid") cycle
            if (line(1:1) == "#") cycle
 
            exit
        end do
 
    end function f_read_line
 
    !> @brief Reads the next non-comment line from a model file, using a buffered look-ahead mechanism.
    impure subroutine s_skip_ignored_lines(iunit, buffered_line, is_buffered)
        integer, intent(in) :: iunit
        character(80), intent(inout) :: buffered_line
        logical, intent(inout) :: is_buffered
 
        character(80) :: line
 
        if (is_buffered) then
            line = buffered_line
            is_buffered = .false.
        else
            if (.not. f_read_line(iunit, line)) return
        end if
 
        buffered_line = line
        is_buffered = .true.
    end subroutine s_skip_ignored_lines
 
    !> This function is used to replace the fortran random number
    !! generator because the native generator is not compatible being called
    !! from GPU routines/functions
    function f_model_random_number(seed) result(rval)
 
        ! $:GPU_ROUTINE(parallelism='[seq]')
 
        integer, intent(inout) :: seed
        real(wp) :: rval
 
        seed = ieor(seed, ishft(seed, 13))
        seed = ieor(seed, ishft(seed, -17))
        seed = ieor(seed, ishft(seed, 5))
 
        rval = abs(real(seed, wp))/real(huge(seed), wp)
    end function f_model_random_number
 
    !> This procedure, recursively, finds whether a point is inside an octree.
    !! @param model    Model to search in.
    !! @param point    Point to test.
    !! @param spacing  Space around the point to search in (grid spacing).
    !! @param spc      Number of samples per cell.
    !! @return True if the point is inside the octree, false otherwise.
    impure function f_model_is_inside(model, point, spacing, spc) result(fraction)
 
        ! $:GPU_ROUTINE(parallelism='[seq]')
 
        type(t_model), intent(in) :: model
        real(wp), dimension(1:3), intent(in) :: point
        real(wp), dimension(1:3), intent(in) :: spacing
        integer, intent(in) :: spc
        real(wp) :: phi, theta
        integer :: rand_seed
 
        real(wp) :: fraction
 
        type(t_ray) :: ray
        integer :: i, j, k, ninorout, nhits
 
        real(wp), dimension(1:spc, 1:3) :: ray_origins, ray_dirs
 
        rand_seed = int(point(1)*73856093._wp) + &
                    int(point(2)*19349663._wp) + &
                    int(point(3)*83492791._wp)
        if (rand_seed == 0) rand_seed = 1
 
        ! generate our random collection or rays
        do i = 1, spc
            do k = 1, 3
                ! random jitter in the origin helps us estimate volume fraction instead of only at the cell center
                ray_origins(i, k) = point(k) + (f_model_random_number(rand_seed) - 0.5_wp)*spacing(k)
                ! cast sample rays in all directions
                ray_dirs(i, k) = f_model_random_number(rand_seed) - 0.5_wp
            end do
            ray_dirs(i, :) = ray_dirs(i, :)/sqrt(sum(ray_dirs(i, :)*ray_dirs(i, :)))
        end do
 
        ! ray trace
        ninorout = 0
        do i = 1, spc
            ray%o = ray_origins(i, :)
            ray%d = ray_dirs(i, :)
 
            nhits = 0
            do j = 1, model%ntrs
                ! count the number of triangles this ray intersects
                if (f_intersects_triangle(ray, model%trs(j)) == 1) then
                    nhits = nhits + 1
                end if
            end do
 
            ! if the ray hits an odd number of triangles on its way out, then
            ! it must be on the inside of the model
            ninorout = ninorout + mod(nhits, 2)
        end do
 
        fraction = real(ninorout)/real(spc)
 
    end function f_model_is_inside
 
    !> This procedure, given a cell center will determine if a point exists instide a surface
    !! @param ntrs     Number of triangles in the model
    !! @param pid      Patch ID od this model
    !! @param point    Point to test.
    !! @return fraction The perfentage of candidate rays cast indicate that we are inside the model
    function f_model_is_inside_flat(ntrs, pid, point) result(fraction)
 
 
# 593 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
#if MFC_OpenACC
# 593 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
!$acc routine seq 
# 593 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
#elif MFC_OpenMP
# 593 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
    
# 593 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
 
# 593 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
!$omp declare target device_type(any) 
# 593 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
#endif
 
        integer, intent(in) :: ntrs
        integer, intent(in) :: pid
        real(wp), dimension(1:3), intent(in) :: point
 
        real(wp) :: fraction
        type(t_ray) :: ray
        type(t_triangle) :: tri
        integer :: i, j, k, q, ninorout, nhits
 
        ! cast 26 rays from the point and count the number at leave the boundary
        ninorout = 0
        do i = -1, 1
            do j = -1, 1
                do k = -1, 1
                    if (i /= 0 .or. j /= 0 .or. k /= 0) then
                        ! We cannot get inersections if the ray is exactly in line with triangle plane
                        if (p == 0 .and. k == 0) cycle
 
                        ! generate the ray
                        ray%o = point
                        ray%d(:) = [real(i, wp), real(j, wp), real(k, wp)]
                        ray%d = ray%d/sqrt(real(abs(i) + abs(j) + abs(k), wp))
 
                        ! count the number of intersections
                        nhits = 0
                        do q = 1, ntrs
                            tri%v(:, :) = gpu_trs_v(:, :, q, pid)
                            tri%n(1:3) = gpu_trs_n(1:3, q, pid)
                            nhits = nhits + f_intersects_triangle(ray, tri)
                        end do
                        ! if the ray intersected an odd number of times, we must be inside
                        ninorout = ninorout + mod(nhits, 2)
                    end if
                end do
            end do
        end do
 
        if (p == 0) then
            ! in 2D, we skipped 8 rays
            fraction = real(ninorout)/18._wp
        else
            fraction = real(ninorout)/26._wp
        end if
 
    end function f_model_is_inside_flat
 
    ! From https://www.scratchapixel.com/lessons/3e-basic-rendering/ray-tracing-rendering-a-triangle/ray-triangle-intersection-geometric-solution.html
    !> This procedure checks if a ray intersects a triangle.
    !! @param ray      Ray.
    !! @param triangle Triangle.
    !! @return         True if the ray intersects the triangle, false otherwise.
    function f_intersects_triangle(ray, triangle) result(intersects)
 
 
# 648 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
#if MFC_OpenACC
# 648 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
!$acc routine seq 
# 648 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
#elif MFC_OpenMP
# 648 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
    
# 648 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
 
# 648 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
!$omp declare target device_type(any) 
# 648 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
#endif
 
        type(t_ray), intent(in) :: ray
        type(t_triangle), intent(in) :: triangle
 
        integer :: intersects
 
        real(wp) :: n(3), p(3), c(3), edge(3), vp(3)
        real(wp) :: d, t, ndotraydirection
 
        intersects = 0
 
        n(1:3) = triangle%n(1:3)
        ndotraydirection = dot_product(n(1:3), ray%d(1:3))
        if (abs(ndotraydirection) < 0.0000001_wp) then
            return
        end if
 
        d = -sum(n(:)*triangle%v(1, :))
        t = -(sum(n(:)*ray%o(:)) + d)/ndotraydirection
        if (t < 0) then
            return
        end if
 
        p = ray%o + t*ray%d
        edge = triangle%v(2, :) - triangle%v(1, :)
        vp = p - triangle%v(1, :)
        c = f_cross(edge, vp)
        if (sum(n(:)*c(:)) < 0) then
            return
        end if
 
        edge = triangle%v(3, :) - triangle%v(2, :)
        vp = p - triangle%v(2, :)
        c = f_cross(edge, vp)
        if (sum(n(:)*c(:)) < 0) then
            return
        end if
 
        edge = triangle%v(1, :) - triangle%v(3, :)
        vp = p - triangle%v(3, :)
        c = f_cross(edge, vp)
        if (sum(n(:)*c(:)) < 0) then
            return
        end if
 
        intersects = 1
 
    end function f_intersects_triangle
 
    !> This procedure checks and labels edges shared by two or more triangles facets of the 2D STL model.
    !! @param model                      Model to search in.
    !! @param boundary_vertex_count      Output total boundary vertex count
    subroutine s_check_boundary(model, boundary_v, boundary_vertex_count, boundary_edge_count)
 
        type(t_model), intent(in) :: model
        real(wp), allocatable, intent(out), dimension(:, :, :) :: boundary_v !< Output boundary vertices/normals
        integer, intent(out) :: boundary_vertex_count, boundary_edge_count !< Output boundary vertex/edge count
 
        integer :: i, j !< Model index iterator
        integer :: edge_count, edge_index, store_index !< Boundary edge index iterator
        real(wp), dimension(1:2, 1:2) :: edge !< Edge end points buffer
        real(wp), dimension(1:2) :: boundary_edge !< Boundary edge end points buffer
        real(wp), dimension(1:(3*model%ntrs), 1:2, 1:2) :: temp_boundary_v !< Temporary boundary vertex buffer
        integer, dimension(1:(3*model%ntrs)) :: edge_occurrence !< The manifoldness of the edges
        real(wp) :: edgetan, initial, v_norm, xnormal, ynormal !< The manifoldness of the edges
 
        ! Total number of edges in 2D STL
        edge_count = 3*model%ntrs
 
        ! Initialize edge_occurrence array to zero
        edge_occurrence = 0
        edge_index = 0
 
        ! Collect all edges of all triangles and store them
        do i = 1, model%ntrs
            ! First edge (v1, v2)
            edge(1, 1:2) = model%trs(i)%v(1, 1:2)
            edge(2, 1:2) = model%trs(i)%v(2, 1:2)
            call s_register_edge(temp_boundary_v, edge, edge_index, edge_count)
 
            ! Second edge (v2, v3)
            edge(1, 1:2) = model%trs(i)%v(2, 1:2)
            edge(2, 1:2) = model%trs(i)%v(3, 1:2)
            call s_register_edge(temp_boundary_v, edge, edge_index, edge_count)
 
            ! Third edge (v3, v1)
            edge(1, 1:2) = model%trs(i)%v(3, 1:2)
            edge(2, 1:2) = model%trs(i)%v(1, 1:2)
            call s_register_edge(temp_boundary_v, edge, edge_index, edge_count)
        end do
 
        ! Check all edges and count repeated edges
 
# 741 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
 
# 741 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
#if defined(MFC_OpenACC)
# 741 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
!$acc parallel loop collapse(2) gang vector default(present) private(i, j) copy(temp_boundary_v, edge_occurrence) 
# 741 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
#elif defined(MFC_OpenMP)
# 741 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
    
# 741 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
 
# 741 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
 
# 741 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
!$omp target teams loop defaultmap(firstprivate:scalar) bind(teams,parallel) collapse(2) defaultmap(tofrom:aggregate) defaultmap(tofrom:allocatable) defaultmap(tofrom:pointer) private(i, j) map(tofrom:temp_boundary_v, edge_occurrence) 
# 741 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
#endif
# 741 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
 
        do i = 1, edge_count
            do j = 1, edge_count
                if (i /= j) then
                    if (((abs(temp_boundary_v(i, 1, 1) - temp_boundary_v(j, 1, 1)) < threshold_edge_zero) .and. &
                         (abs(temp_boundary_v(i, 1, 2) - temp_boundary_v(j, 1, 2)) < threshold_edge_zero) .and. &
                         (abs(temp_boundary_v(i, 2, 1) - temp_boundary_v(j, 2, 1)) < threshold_edge_zero) .and. &
                         (abs(temp_boundary_v(i, 2, 2) - temp_boundary_v(j, 2, 2)) < threshold_edge_zero)) .or. &
                        ((abs(temp_boundary_v(i, 1, 1) - temp_boundary_v(j, 2, 1)) < threshold_edge_zero) .and. &
                         (abs(temp_boundary_v(i, 1, 2) - temp_boundary_v(j, 2, 2)) < threshold_edge_zero) .and. &
                         (abs(temp_boundary_v(i, 2, 1) - temp_boundary_v(j, 1, 1)) < threshold_edge_zero) .and. &
                         (abs(temp_boundary_v(i, 2, 2) - temp_boundary_v(j, 1, 2)) < threshold_edge_zero))) then
 
 
# 754 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
#if defined(MFC_OpenACC)
# 754 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
!$acc atomic update
# 754 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
#elif defined(MFC_OpenMP)
# 754 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
!$omp atomic update
# 754 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
#endif
                        edge_occurrence(i) = edge_occurrence(i) + 1
                    end if
                end if
            end do
        end do
 
# 760 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
 
# 760 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
#if defined(MFC_OpenACC)
# 760 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
!$acc end parallel loop
# 760 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
#elif defined(MFC_OpenMP)
# 760 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
    
# 760 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
 
# 760 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
!$omp end target teams loop
# 760 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
#endif
# 760 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
 
 
        ! Count the number of boundary vertices/edges
        boundary_vertex_count = 0
        boundary_edge_count = 0
 
        do i = 1, edge_count
            if (edge_occurrence(i) == 0) then
                boundary_vertex_count = boundary_vertex_count + 2
                boundary_edge_count = boundary_edge_count + 1
            end if
        end do
 
        ! Allocate the boundary_v array based on the number of boundary edges
        allocate (boundary_v(boundary_edge_count, 1:3, 1:2))
 
        ! Store boundary vertices
        store_index = 0
        do i = 1, edge_count
            if (edge_occurrence(i) == 0) then
                store_index = store_index + 1
                boundary_v(store_index, 1, 1:2) = temp_boundary_v(i, 1, 1:2)
                boundary_v(store_index, 2, 1:2) = temp_boundary_v(i, 2, 1:2)
            end if
        end do
 
        ! Find/store the normal vector of the boundary edges
        do i = 1, boundary_edge_count
            boundary_edge(1) = boundary_v(i, 2, 1) - boundary_v(i, 1, 1)
            boundary_edge(2) = boundary_v(i, 2, 2) - boundary_v(i, 1, 2)
            edgetan = boundary_edge(1)/sign(max(sgm_eps, abs(boundary_edge(2))), boundary_edge(2))
 
            if (abs(boundary_edge(2)) < threshold_vector_zero) then
                if (edgetan > 0._wp) then
                    ynormal = -1
                    xnormal = 0._wp
                else
                    ynormal = 1
                    xnormal = 0._wp
                end if
            else
                initial = boundary_edge(2)
                ynormal = -edgetan*initial
                xnormal = initial
            end if
 
            v_norm = sqrt(xnormal**2 + ynormal**2)
            boundary_v(i, 3, 1) = xnormal/v_norm
            boundary_v(i, 3, 2) = ynormal/v_norm
        end do
 
    end subroutine s_check_boundary
 
    !> This procedure appends the edge end vertices to a temporary buffer.
    subroutine s_register_edge(temp_boundary_v, edge, edge_index, edge_count)
 
        integer, intent(inout) :: edge_index !< Edge index iterator
        integer, intent(inout) :: edge_count !< Total number of edges
        real(wp), intent(in), dimension(1:2, 1:2) :: edge !< Edges end points to be registered
        real(wp), dimension(1:edge_count, 1:2, 1:2), intent(inout) :: temp_boundary_v !< Temporary edge end vertex buffer
 
        ! Increment edge index and store the edge
        edge_index = edge_index + 1
        temp_boundary_v(edge_index, 1, 1:2) = edge(1, 1:2)
        temp_boundary_v(edge_index, 2, 1:2) = edge(2, 1:2)
 
    end subroutine s_register_edge
 
    !> This procedure determines the levelset distance and normals of 3D models
    !! by computing the exact closest point via projection onto triangle surfaces.
    !! @param ntrs                  Number of triangles for this patch
    !! @param trs_v                 Flat GPU array of triangle vertices for all patches
    !! @param trs_n                 Flat GPU array of triangle normals for all patches
    !! @param pid                   Patch index into the arrays
    !! @param point                 The cell center of the current levelset cell
    !! @param normals               Output levelset normals
    !! @param distance              Output levelset distance
    subroutine s_distance_normals_3d(ntrs, pid, point, normals, distance)
 
# 838 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
#if MFC_OpenACC
# 838 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
!$acc routine seq 
# 838 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
#elif MFC_OpenMP
# 838 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
    
# 838 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
 
# 838 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
!$omp declare target device_type(any) 
# 838 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
#endif
 
        integer, intent(in) :: ntrs
        integer, intent(in) :: pid
        real(wp), dimension(1:3), intent(in) :: point
        real(wp), dimension(1:3), intent(out) :: normals
        real(wp), intent(out) :: distance
 
        integer :: i, j, l
        real(wp) :: dist_min, dist_proj, dist_v, dist_e, t
        real(wp) :: v1(1:3), v2(1:3), v3(1:3)
        real(wp) :: e0(1:3), e1(1:3), pv(1:3)
        real(wp) :: n(1:3), proj(1:3), norm_vec(1:3)
        real(wp) :: d, ndot, denom, norm_mag
        real(wp) :: u, v_bary, w
        real(wp) :: l00, l01, l11, l20, l21
        real(wp) :: edge(1:3), pe(1:3)
        real(wp) :: verts(1:3, 1:3)
 
        dist_min = initial_distance_buffer
        normals = 0._wp
 
        do i = 1, ntrs
            ! Triangle vertices
            v1(:) = gpu_trs_v(1, :, i, pid)
            v2(:) = gpu_trs_v(2, :, i, pid)
            v3(:) = gpu_trs_v(3, :, i, pid)
 
            ! Triangle normal
            n(:) = gpu_trs_n(:, i, pid)
 
            ! Project point onto triangle plane
            pv(:) = point(:) - v1(:)
            d = dot_product(pv, n)
            if (abs(d) >= dist_min) cycle ! minimum distance is not small enough, no need to check validity
            proj(:) = point(:) - d*n(:)
 
            ! Check if projection is inside triangle using barycentric coordinates
            e0(:) = v2(:) - v1(:)
            e1(:) = v3(:) - v1(:)
            pv(:) = proj(:) - v1(:)
 
            l00 = dot_product(e0, e0)
            l01 = dot_product(e0, e1)
            l11 = dot_product(e1, e1)
            l20 = dot_product(pv, e0)
            l21 = dot_product(pv, e1)
 
            denom = l00*l11 - l01*l01
 
            ! compute the barycentric coordinates of the projection in the triangle
            if (abs(denom) > 0._wp) then
                v_bary = (l11*l20 - l01*l21)/denom
                w = (l00*l21 - l01*l20)/denom
                u = 1._wp - v_bary - w
            else
                u = -1._wp
                v_bary = -1._wp
                w = -1._wp
            end if
 
            ! If projection is inside triangle
            if (u >= 0._wp .and. v_bary >= 0._wp .and. w >= 0._wp) then
                dist_proj = sqrt((point(1) - proj(1))**2 + &
                                 (point(2) - proj(2))**2 + &
                                 (point(3) - proj(3))**2)
 
                if (dist_proj < dist_min) then
                    dist_min = dist_proj
                    normals(:) = n(:)
                end if
            else
                ! Projection outside triangle: check edges and vertices
                verts(:, 1) = v1(:)
                verts(:, 2) = v2(:)
                verts(:, 3) = v3(:)
 
                ! Check three edges
                do j = 1, 3
                    edge(:) = verts(:, mod(j, 3) + 1) - verts(:, j)
                    pe(:) = point(:) - verts(:, j)
 
                    t = dot_product(pe, edge)/max(dot_product(edge, edge), 1.e-30_wp)
 
                    if (t >= 0._wp .and. t <= 1._wp) then
                        proj(:) = verts(:, j) + t*edge(:)
                        dist_e = sqrt((point(1) - proj(1))**2 + &
                                      (point(2) - proj(2))**2 + &
                                      (point(3) - proj(3))**2)
 
                        if (dist_e < dist_min) then
                            dist_min = dist_e
                            norm_vec(:) = proj(:) - point(:)
                            if (dist_e > 0._wp) norm_vec = norm_vec/dist_e
                            ! Snap to triangle normal if nearly parallel
                            if (f_approx_equal(dot_product(norm_vec, n), 1._wp)) then
                                normals(:) = n(:)
                            else
                                normals(:) = norm_vec(:)
                            end if
                        end if
                    else if (t < 0._wp) then
                        dist_v = sqrt((point(1) - verts(1, j))**2 + &
                                      (point(2) - verts(2, j))**2 + &
                                      (point(3) - verts(3, j))**2)
 
                        if (dist_v < dist_min) then
                            dist_min = dist_v
                            norm_vec(:) = verts(:, j) - point(:)
                            norm_mag = sqrt(dot_product(norm_vec, norm_vec))
                            if (norm_mag > 0._wp) norm_vec = norm_vec/norm_mag
                            normals(:) = norm_vec(:)
                        end if
                    else
                        dist_v = sqrt((point(1) - verts(1, mod(j, 3) + 1))**2 + &
                                      (point(2) - verts(2, mod(j, 3) + 1))**2 + &
                                      (point(3) - verts(3, mod(j, 3) + 1))**2)
 
                        if (dist_v < dist_min) then
                            dist_min = dist_v
                            norm_vec(:) = verts(:, mod(j, 3) + 1) - point(:)
                            norm_mag = sqrt(dot_product(norm_vec, norm_vec))
                            if (norm_mag > 0._wp) norm_vec = norm_vec/norm_mag
                            normals(:) = norm_vec(:)
                        end if
                    end if
                end do
            end if
        end do
 
        distance = dist_min
 
    end subroutine s_distance_normals_3d
 
    !> This procedure determines the levelset distance and normals of 2D models
    !! by computing the exact closest point via projection onto boundary edges.
    !! @param boundary_v            Flat GPU array of boundary vertices/normals for all patches
    !! @param pid                   Patch index into the boundary_v array
    !! @param boundary_edge_count   Total number of boundary edges for this patch
    !! @param point                 The cell center of the current levelset cell
    !! @param normals               Output levelset normals
    !! @param distance              Output levelset distance
    subroutine s_distance_normals_2d(pid, boundary_edge_count, point, normals, distance)
 
# 981 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
#if MFC_OpenACC
# 981 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
!$acc routine seq 
# 981 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
#elif MFC_OpenMP
# 981 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
    
# 981 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
 
# 981 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
!$omp declare target device_type(any) 
# 981 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
#endif
 
        integer, intent(in) :: pid
        integer, intent(in) :: boundary_edge_count
        real(wp), dimension(1:3), intent(in) :: point
        real(wp), dimension(1:3), intent(out) :: normals
        real(wp), intent(out) :: distance
 
        integer :: i
        real(wp) :: dist_min, dist, t, norm_mag
        real(wp) :: v1(1:2), v2(1:2), edge(1:2), pv(1:2)
        real(wp) :: edge_len_sq, proj(1:2), norm(1:2), c
 
        dist_min = initial_distance_buffer
        normals = 0._wp
        norm = 0._wp
 
        do i = 1, boundary_edge_count
            ! Edge endpoints
            v1(1) = gpu_boundary_v(i, 1, 1, pid)
            v1(2) = gpu_boundary_v(i, 1, 2, pid)
            v2(1) = gpu_boundary_v(i, 2, 1, pid)
            v2(2) = gpu_boundary_v(i, 2, 2, pid)
 
            ! Edge vector and point-to-v1 vector
            edge = v2 - v1
            pv(1) = point(1) - v1(1)
            pv(2) = point(2) - v1(2)
            edge_len_sq = dot_product(edge, edge)
 
            ! Parameter of projection onto the edge line
            if (edge_len_sq > 0._wp) then
                t = dot_product(pv, edge)/edge_len_sq
            else
                t = 0._wp
            end if
 
            ! Check if projection falls within the segment
            if (t >= 0._wp .and. t <= 1._wp) then
                proj = v1 + t*edge
                dist = sqrt((point(1) - proj(1))**2 + (point(2) - proj(2))**2)
                norm(1) = gpu_boundary_v(i, 3, 1, pid)
                norm(2) = gpu_boundary_v(i, 3, 2, pid)
            else if (t < 0._wp) then ! negative t means that v1 is the closest point on the edge
                dist = sqrt((point(1) - v1(1))**2 + (point(2) - v1(2))**2)
                norm(1) = v1(1) - point(1)
                norm(2) = v1(2) - point(2)
                norm = norm/dist
            else ! t > 1 means that v2 is the closest point on the line edge
                dist = sqrt((point(1) - v2(1))**2 + (point(2) - v2(2))**2)
                norm(1) = v2(1) - point(1)
                norm(2) = v2(2) - point(2)
                norm = norm/dist
            end if
 
            if (dist < dist_min) then
                dist_min = dist
                normals(1) = norm(1)
                normals(2) = norm(2)
            end if
        end do
 
        distance = dist_min
 
    end subroutine s_distance_normals_2d
 
#ifdef MFC_SIMULATION
 
    subroutine s_instantiate_stl_models()
 
        ! Variables for IBM+STL
        real(wp) :: normals(1:3) !< Boundary normal buffer
        integer :: boundary_vertex_count, boundary_edge_count, total_vertices !< Boundary vertex
        real(wp), allocatable, dimension(:, :, :) :: boundary_v !< Boundary vertex buffer
        real(wp) :: dx_local, dy_local, dz_local !< Levelset distance buffer
 
        integer :: i, j, k !< Generic loop iterators
        integer :: patch_id
 
        type(t_bbox) :: bbox, bbox_old
        type(t_model) :: model
        type(ic_model_parameters) :: params
 
        real(wp) :: eta
        real(wp), dimension(1:3) :: point, model_center
        real(wp) :: grid_mm(1:3, 1:2)
 
        real(wp), dimension(1:4, 1:4) :: transform, transform_n
 
        dx_local = minval(dx); dy_local = minval(dy)
        if (p /= 0) dz_local = minval(dz)
 
        allocate (stl_bounding_boxes(num_ibs, 1:3, 1:3))
 
        do patch_id = 1, num_ibs
            if (patch_ib(patch_id)%geometry == 5 .or. patch_ib(patch_id)%geometry == 12) then
                allocate (models(patch_id)%model)
                print *, " * Reading model: "//trim(patch_ib(patch_id)%model_filepath)
 
                model = f_model_read(patch_ib(patch_id)%model_filepath)
                params%scale(:) = patch_ib(patch_id)%model_scale(:)
                params%translate(:) = patch_ib(patch_id)%model_translate(:)
                params%rotate(:) = patch_ib(patch_id)%model_rotate(:)
                params%spc = patch_ib(patch_id)%model_spc
                params%threshold = patch_ib(patch_id)%model_threshold
 
                if (f_approx_equal(dot_product(params%scale, params%scale), 0._wp)) then
                    params%scale(:) = 1._wp
                end if
 
                if (proc_rank == 0) then
                    print *, " * Transforming model."
                end if
 
                ! Get the model center before transforming the model
                bbox_old = f_create_bbox(model)
                model_center(1:3) = (bbox_old%min(1:3) + bbox_old%max(1:3))/2._wp
 
                ! Compute the transform matrices for vertices and normals
                transform = f_create_transform_matrix(params, model_center)
                transform_n = f_create_transform_matrix(params)
 
                call s_transform_model(model, transform, transform_n)
 
                ! Recreate the bounding box after transformation
                bbox = f_create_bbox(model)
 
                ! Show the number of vertices in the original STL model
                if (proc_rank == 0) then
                    print *, ' * Number of input model vertices:', 3*model%ntrs
                end if
 
                ! Need the cells that form the boundary of the flat projection in 2D
                if (p == 0) call s_check_boundary(model, boundary_v, boundary_vertex_count, boundary_edge_count)
 
                ! Show the number of edges and boundary edges in 2D STL models
                if (proc_rank == 0 .and. p == 0) then
                    print *, ' * Number of 2D model boundary edges:', boundary_edge_count
                end if
 
                if (proc_rank == 0) then
                    write (*, "(A, 3(2X, F20.10))") "    > Model:  Min:", bbox%min(1:3)
                    write (*, "(A, 3(2X, F20.10))") "    >         Cen:", (bbox%min(1:3) + bbox%max(1:3))/2._wp
                    write (*, "(A, 3(2X, F20.10))") "    >         Max:", bbox%max(1:3)
 
                    grid_mm(1, :) = (/minval(x_cc(0:m)) - 0.5_wp*dx_local, maxval(x_cc(0:m)) + 0.5_wp*dx_local/)
                    grid_mm(2, :) = (/minval(y_cc(0:n)) - 0.5_wp*dy_local, maxval(y_cc(0:n)) + 0.5_wp*dy_local/)
 
                    if (p > 0) then
                        grid_mm(3, :) = (/minval(z_cc(0:p)) - 0.5_wp*dz_local, maxval(z_cc(0:p)) + 0.5_wp*dz_local/)
                    else
                        grid_mm(3, :) = (/0._wp, 0._wp/)
                    end if
 
                    write (*, "(A, 3(2X, F20.10))") "    > Domain: Min:", grid_mm(:, 1)
                    write (*, "(A, 3(2X, F20.10))") "    >         Cen:", (grid_mm(:, 1) + grid_mm(:, 2))/2._wp
                    write (*, "(A, 3(2X, F20.10))") "    >         Max:", grid_mm(:, 2)
                end if
 
                stl_bounding_boxes(patch_id, 1, 1:3) = [bbox%min(1), (bbox%min(1) + bbox%max(1))/2._wp, bbox%max(1)]
                stl_bounding_boxes(patch_id, 2, 1:3) = [bbox%min(2), (bbox%min(2) + bbox%max(2))/2._wp, bbox%max(2)]
                stl_bounding_boxes(patch_id, 3, 1:3) = [bbox%min(3), (bbox%min(3) + bbox%max(3))/2._wp, bbox%max(3)]
 
                models(patch_id)%model = model
                if (p == 0) then
                    models(patch_id)%boundary_v = boundary_v
                    models(patch_id)%boundary_edge_count = boundary_edge_count
                end if
            end if
        end do
 
        ! Pack and upload flat arrays for GPU (AFTER the loop)
        block
            integer :: pid, max_ntrs
            integer :: max_bv1, max_bv2, max_bv3, max_iv1, max_iv2
 
            max_ntrs = 0
            max_bv1 = 0; max_bv2 = 0; max_bv3 = 0
            max_iv1 = 0; max_iv2 = 0
 
            do pid = 1, num_ibs
                if (allocated(models(pid)%model)) then
                    call s_pack_model_for_gpu(models(pid))
                    max_ntrs = max(max_ntrs, models(pid)%ntrs)
                end if
                if (allocated(models(pid)%boundary_v)) then
                    max_bv1 = max(max_bv1, size(models(pid)%boundary_v, 1))
                    max_bv2 = max(max_bv2, size(models(pid)%boundary_v, 2))
                    max_bv3 = max(max_bv3, size(models(pid)%boundary_v, 3))
                end if
            end do
 
            if (max_ntrs > 0) then
#ifdef MFC_DEBUG
# 1174 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
    block
# 1174 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
        use iso_fortran_env, only: output_unit
# 1174 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
 
# 1174 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
        print *, 'm_model.fpp:1174: ', '@:ALLOCATE(gpu_ntrs(1:num_ibs))'
# 1174 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
 
# 1174 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
        call flush (output_unit)
# 1174 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
    end block
# 1174 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
#endif
# 1174 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
    allocate (gpu_ntrs(1:num_ibs))
# 1174 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
 
# 1174 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
 
# 1174 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
#if defined(MFC_OpenACC)
# 1174 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
!$acc enter data create(gpu_ntrs) 
# 1174 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
#elif defined(MFC_OpenMP)
# 1174 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
!$omp target enter data map(always,alloc:gpu_ntrs) 
# 1174 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
#endif
#ifdef MFC_DEBUG
# 1175 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
    block
# 1175 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
        use iso_fortran_env, only: output_unit
# 1175 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
 
# 1175 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
        print *, 'm_model.fpp:1175: ', '@:ALLOCATE(gpu_trs_v(1:3, 1:3, 1:max_ntrs, 1:num_ibs))'
# 1175 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
 
# 1175 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
        call flush (output_unit)
# 1175 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
    end block
# 1175 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
#endif
# 1175 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
    allocate (gpu_trs_v(1:3, 1:3, 1:max_ntrs, 1:num_ibs))
# 1175 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
 
# 1175 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
 
# 1175 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
#if defined(MFC_OpenACC)
# 1175 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
!$acc enter data create(gpu_trs_v) 
# 1175 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
#elif defined(MFC_OpenMP)
# 1175 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
!$omp target enter data map(always,alloc:gpu_trs_v) 
# 1175 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
#endif
#ifdef MFC_DEBUG
# 1176 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
    block
# 1176 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
        use iso_fortran_env, only: output_unit
# 1176 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
 
# 1176 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
        print *, 'm_model.fpp:1176: ', '@:ALLOCATE(gpu_trs_n(1:3, 1:max_ntrs, 1:num_ibs))'
# 1176 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
 
# 1176 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
        call flush (output_unit)
# 1176 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
    end block
# 1176 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
#endif
# 1176 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
    allocate (gpu_trs_n(1:3, 1:max_ntrs, 1:num_ibs))
# 1176 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
 
# 1176 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
 
# 1176 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
#if defined(MFC_OpenACC)
# 1176 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
!$acc enter data create(gpu_trs_n) 
# 1176 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
#elif defined(MFC_OpenMP)
# 1176 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
!$omp target enter data map(always,alloc:gpu_trs_n) 
# 1176 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
#endif
#ifdef MFC_DEBUG
# 1177 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
    block
# 1177 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
        use iso_fortran_env, only: output_unit
# 1177 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
 
# 1177 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
        print *, 'm_model.fpp:1177: ', '@:ALLOCATE(gpu_boundary_edge_count(1:num_ibs))'
# 1177 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
 
# 1177 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
        call flush (output_unit)
# 1177 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
    end block
# 1177 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
#endif
# 1177 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
    allocate (gpu_boundary_edge_count(1:num_ibs))
# 1177 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
 
# 1177 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
 
# 1177 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
#if defined(MFC_OpenACC)
# 1177 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
!$acc enter data create(gpu_boundary_edge_count) 
# 1177 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
#elif defined(MFC_OpenMP)
# 1177 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
!$omp target enter data map(always,alloc:gpu_boundary_edge_count) 
# 1177 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
#endif
#ifdef MFC_DEBUG
# 1178 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
    block
# 1178 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
        use iso_fortran_env, only: output_unit
# 1178 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
 
# 1178 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
        print *, 'm_model.fpp:1178: ', '@:ALLOCATE(gpu_total_vertices(1:num_ibs))'
# 1178 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
 
# 1178 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
        call flush (output_unit)
# 1178 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
    end block
# 1178 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
#endif
# 1178 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
    allocate (gpu_total_vertices(1:num_ibs))
# 1178 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
 
# 1178 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
 
# 1178 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
#if defined(MFC_OpenACC)
# 1178 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
!$acc enter data create(gpu_total_vertices) 
# 1178 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
#elif defined(MFC_OpenMP)
# 1178 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
!$omp target enter data map(always,alloc:gpu_total_vertices) 
# 1178 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
#endif
 
                gpu_ntrs = 0
                gpu_trs_v = 0._wp
                gpu_trs_n = 0._wp
                gpu_boundary_edge_count = 0
                gpu_total_vertices = 0
 
                if (max_bv1 > 0) then
#ifdef MFC_DEBUG
# 1187 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
    block
# 1187 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
        use iso_fortran_env, only: output_unit
# 1187 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
 
# 1187 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
        print *, 'm_model.fpp:1187: ', '@:ALLOCATE(gpu_boundary_v(1:max_bv1, 1:max_bv2, 1:max_bv3, 1:num_ibs))'
# 1187 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
 
# 1187 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
        call flush (output_unit)
# 1187 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
    end block
# 1187 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
#endif
# 1187 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
    allocate (gpu_boundary_v(1:max_bv1, 1:max_bv2, 1:max_bv3, 1:num_ibs))
# 1187 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
 
# 1187 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
 
# 1187 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
#if defined(MFC_OpenACC)
# 1187 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
!$acc enter data create(gpu_boundary_v) 
# 1187 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
#elif defined(MFC_OpenMP)
# 1187 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
!$omp target enter data map(always,alloc:gpu_boundary_v) 
# 1187 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
#endif
                    gpu_boundary_v = 0._wp
                end if
 
                do pid = 1, num_ibs
                    if (allocated(models(pid)%model)) then
                        gpu_ntrs(pid) = models(pid)%ntrs
                        gpu_trs_v(:, :, 1:models(pid)%ntrs, pid) = models(pid)%trs_v
                        gpu_trs_n(:, 1:models(pid)%ntrs, pid) = models(pid)%trs_n
                        gpu_boundary_edge_count(pid) = models(pid)%boundary_edge_count
                        gpu_total_vertices(pid) = models(pid)%total_vertices
                    end if
                    if (allocated(models(pid)%boundary_v) .and. p == 0) then
                        gpu_boundary_v(1:size(models(pid)%boundary_v, 1), &
                                       1:size(models(pid)%boundary_v, 2), &
                                       1:size(models(pid)%boundary_v, 3), pid) = models(pid)%boundary_v
                    end if
                end do
 
 
# 1206 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
#if defined(MFC_OpenACC)
# 1206 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
!$acc update device(gpu_ntrs, gpu_trs_v, gpu_trs_n, gpu_boundary_edge_count, gpu_total_vertices) 
# 1206 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
#elif defined(MFC_OpenMP)
# 1206 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
!$omp target update to(gpu_ntrs, gpu_trs_v, gpu_trs_n, gpu_boundary_edge_count, gpu_total_vertices) 
# 1206 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
#endif
                if (allocated(gpu_boundary_v)) then
 
# 1208 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
#if defined(MFC_OpenACC)
# 1208 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
!$acc update device(gpu_boundary_v) 
# 1208 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
#elif defined(MFC_OpenMP)
# 1208 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
!$omp target update to(gpu_boundary_v) 
# 1208 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
#endif
                end if
            end if
        end block
 
    end subroutine s_instantiate_stl_models
 
#endif
 
    subroutine s_pack_model_for_gpu(ma)
        type(t_model_array), intent(inout) :: ma
        integer :: i
 
        ma%ntrs = ma%model%ntrs
        allocate (ma%trs_v(1:3, 1:3, 1:ma%ntrs))
        allocate (ma%trs_n(1:3, 1:ma%ntrs))
 
        do i = 1, ma%ntrs
            ma%trs_v(:, :, i) = ma%model%trs(i)%v(:, :)
            ma%trs_n(:, i) = ma%model%trs(i)%n(:)
        end do
    end subroutine
 
end module m_model