Removes the macros t_vec3 and the t_4x4 matrix one

[sbryngelson]

User description

Removes some macros and replaces them with normal text. Doesn't make the code longer or harder to read. removing macros seems like a good idea when appropriate.

---

PR Type

Enhancement

---

Description

• Remove unnecessary t_vec3 and t_mat4x4 macros

• Replace macro usage with explicit type declarations

• Improve code readability by using standard Fortran syntax

---

Changes diagram

flowchart LR
  A["macros.fpp"] -- "remove macros" --> B["Type definitions"]
  B -- "replace usage" --> C["m_derived_types.fpp"]
  B -- "replace usage" --> D["m_helper.fpp"]
  B -- "replace usage" --> E["m_model.fpp"]
  B -- "replace usage" --> F["m_patches.fpp"]

Loading

---

Changes walkthrough 📝

Relevant files

Enhancement

macros.fpp Remove vector and matrix macros                                                    src/common/include/macros.fpp Remove t_vec3 macro definition Remove t_mat4x4 macro definition +0/-3      m_derived_types.fpp Replace macros in derived types                                                    src/common/m_derived_types.fpp Replace t_vec3 with real(wp), dimension(1:3) in type definitions Update triangle, ray, and bbox type declarations Update model parameter type declarations +14/-14  m_helper.fpp Replace macros in helper functions                                              src/common/m_helper.fpp Replace t_vec3 and t_mat4x4 in function parameters Update transformation matrix function declarations Update vector and triangle transformation subroutines +6/-6      m_model.fpp Replace macros in model processing                                              src/pre_process/m_model.fpp Replace t_vec3 with explicit type in variable declarations Update function parameters and local variables Maintain same functionality with explicit types +19/-19  m_patches.fpp Replace macros in patch processing                                              src/pre_process/m_patches.fpp Replace t_vec3 and t_mat4x4 in local variables Update coordinate conversion function parameters Maintain geometric transformation functionality +4/-4

---

Need help?

Type /help how to ... in the comments thread for any questions about Qodo Merge usage.

Check out the documentation for more information.

[Copilot]

Pull Request Overview

This PR removes the t_vec3 and t_mat4x4 macros and replaces all their usages with explicit real(wp), dimension(...) declarations across the codebase.

• Replaces t_vec3 with real(wp), dimension(1:3) in various modules.
• Replaces t_mat4x4 with real(wp), dimension(1:4,1:4) in helper routines.
• Deletes the macro definitions from macros.fpp.

Reviewed Changes

Copilot reviewed 5 out of 5 changed files in this pull request and generated 1 comment.

Show a summary per file

File	Description
src/pre_process/m_patches.fpp	Swapped vector/matrix macros for explicit types
src/pre_process/m_model.fpp	Updated all t_vec3 usages to explicit arrays
src/common/m_helper.fpp	Updated transform routines to explicit matrices
src/common/m_derived_types.fpp	Replaced macro-based fields with explicit types
src/common/include/macros.fpp	Removed t_vec3/t_mat4x4 macro definitions

[Copilot]

[nitpick] The repeated explicit declarations for 3- and 4×4 arrays throughout the code increase duplication. Consider defining small derived types (e.g., type vec3 and type mat4x4) in a central module to improve readability and reduce boilerplate.

Suggested change

	real(wp), dimension(1:3), optional, intent(in) :: center
	real(wp), dimension(1:4,1:4) :: sc, rz, rx, ry, tr, t_back, t_to_origin, out_matrix
	type(vec3), optional, intent(in) :: center
	type(mat4x4) :: sc, rz, rx, ry, tr, t_back, t_to_origin, out_matrix

[qodo-code-review]

PR Reviewer Guide 🔍

Here are some key observations to aid the review process:

⏱️ Estimated effort to review: 2 🔵🔵⚪⚪⚪

🧪 No relevant tests

🔒 No security concerns identified

⚡ Recommended focus areas for review Code Duplication The same type declaration real(wp), dimension(1:3) is repeated multiple times across different type definitions (ic_model_parameters, t_triangle, t_ray, t_bbox, and patch types). This creates maintenance overhead and potential inconsistency if the type needs to be changed in the future. real(wp), dimension(1:3) :: translate !< !! Translation of the STL object. real(wp), dimension(1:3) :: scale !< !! Scale factor for the STL object. real(wp), dimension(1:3) :: rotate !< !! Angle to rotate the STL object along each cartesian coordinate axis, !! in radians. integer :: spc !< !! Number of samples per cell to use when discretizing the STL object. real(wp) :: threshold !< !! Threshold to turn on smoothen STL patch. end type ic_model_parameters type :: t_triangle real(wp), dimension(1:3, 1:3) :: v ! Vertices of the triangle real(wp), dimension(1:3) :: n ! Normal vector end type t_triangle type :: t_ray real(wp), dimension(1:3) :: o ! Origin real(wp), dimension(1:3) :: d ! Direction end type t_ray type :: t_bbox real(wp), dimension(1:3) :: min ! Minimum coordinates real(wp), dimension(1:3) :: max ! Maximum coordinates end type t_bbox type :: t_model integer :: ntrs ! Number of triangles type(t_triangle), allocatable :: trs(:) ! Triangles end type t_model !> Derived type adding initial condition (ic) patch parameters as attributes !! NOTE: The requirements for the specification of the above parameters !! are strongly dependent on both the choice of the multicomponent flow !! model as well as the choice of the patch geometry. type ic_patch_parameters integer :: geometry !< Type of geometry for the patch real(wp) :: x_centroid, y_centroid, z_centroid !< !! Location of the geometric center, i.e. the centroid, of the patch. It !! is specified through its x-, y- and z-coordinates, respectively. real(wp) :: length_x, length_y, length_z !< Dimensions of the patch. x,y,z Lengths. real(wp) :: radius !< Dimensions of the patch. radius. real(wp), dimension(3) :: radii !< !! Vector indicating the various radii for the elliptical and ellipsoidal !! patch geometries. It is specified through its x-, y-, and z-components !! respectively. real(wp) :: epsilon, beta !< !! The isentropic vortex parameters for the amplitude of the disturbance and !! domain of influence. real(wp), dimension(2:9) :: a !< !! The parameters needed for the spherical harmonic patch logical :: non_axis_sym real(wp), dimension(3) :: normal !< !! Normal vector indicating the orientation of the patch. It is specified !! through its x-, y- and z-components, respectively. logical, dimension(0:num_patches_max - 1) :: alter_patch !< !! List of permissions that indicate to the current patch which preceding !! patches it is allowed to overwrite when it is in process of being laid !! out in the domain logical :: smoothen !< !! Permission indicating to the current patch whether its boundaries will !! be smoothed out across a few cells or whether they are to remain sharp integer :: smooth_patch_id !< !! Identity (id) of the patch with which current patch is to get smoothed real(wp) :: smooth_coeff !< !! Smoothing coefficient (coeff) for the size of the stencil of !! cells across which boundaries of the current patch will be smeared out real(wp), dimension(num_fluids_max) :: alpha_rho real(wp) :: rho real(wp), dimension(3) :: vel real(wp) :: pres real(wp), dimension(num_fluids_max) :: alpha real(wp) :: gamma real(wp) :: pi_inf !< real(wp) :: cv !< real(wp) :: qv !< real(wp) :: qvp !< !! Primitive variables associated with the patch. In order, these include !! the partial densities, density, velocity, pressure, volume fractions, !! specific heat ratio function and the liquid stiffness function. real(wp) :: Bx, By, Bz !< !! Magnetic field components; B%x is not used for 1D real(wp), dimension(6) :: tau_e !< !! Elastic stresses added to primitive variables if hypoelasticity = True real(wp) :: R0 !< Bubble size real(wp) :: V0 !< Bubble velocity real(wp) :: p0 !< Bubble size real(wp) :: m0 !< Bubble velocity integer :: hcid !! id for hard coded initial condition real(wp) :: cf_val !! color function value real(wp) :: Y(1:num_species) !! STL or OBJ model input parameter character(LEN=pathlen_max) :: model_filepath !< !! Path the STL file relative to case_dir. real(wp), dimension(1:3) :: model_translate !< !! Translation of the STL object. real(wp), dimension(1:3) :: model_scale !< !! Scale factor for the STL object. real(wp), dimension(1:3) :: model_rotate !< !! Angle to rotate the STL object along each cartesian coordinate axis, !! in radians. integer :: model_spc !< !! Number of samples per cell to use when discretizing the STL object. real(wp) :: model_threshold !< !! Threshold to turn on smoothen STL patch. end type ic_patch_parameters type ib_patch_parameters integer :: geometry !< Type of geometry for the patch real(wp) :: x_centroid, y_centroid, z_centroid !< !! Location of the geometric center, i.e. the centroid, of the patch. It !! is specified through its x-, y- and z-coordinates, respectively. real(wp) :: c, p, t, m real(wp) :: length_x, length_y, length_z !< Dimensions of the patch. x,y,z Lengths. real(wp) :: radius !< Dimensions of the patch. radius. real(wp) :: theta logical :: slip !! STL or OBJ model input parameter character(LEN=pathlen_max) :: model_filepath !< !! Path the STL file relative to case_dir. real(wp), dimension(1:3) :: model_translate !< !! Translation of the STL object. real(wp), dimension(1:3) :: model_scale !< !! Scale factor for the STL object. real(wp), dimension(1:3) :: model_rotate !< !! Angle to rotate the STL object along each cartesian coordinate axis, !! in radians. Missing Tests The PR removes macros that were used in geometric transformation and model processing functions, but no tests are provided to verify that the explicit type declarations maintain the same functionality and precision as the original macro-based implementation. !! @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. pure 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 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 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) 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) :: fraction type(t_ray) :: ray integer :: i, j, nInOrOut, nHits real(wp), dimension(1:spc, 1:3) :: ray_origins, ray_dirs do i = 1, spc call random_number(ray_origins(i, :)) ray_origins(i, :) = point + (ray_origins(i, :) - 0.5_wp)*spacing(:) call random_number(ray_dirs(i, :)) ray_dirs(i, :) = ray_dirs(i, :) - 0.5_wp ray_dirs(i, :) = ray_dirs(i, :)/sqrt(sum(ray_dirs(i, :)*ray_dirs(i, :))) end do nInOrOut = 0 do i = 1, spc ray%o = ray_origins(i, :) ray%d = ray_dirs(i, :) nHits = 0 do j = 1, model%ntrs if (f_intersects_triangle(ray, model%trs(j))) then nHits = nHits + 1 end if end do nInOrOut = nInOrOut + mod(nHits, 2) end do fraction = real(nInOrOut)/real(spc) end function f_model_is_inside ! 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. pure elemental function f_intersects_triangle(ray, triangle) result(intersects) type(t_ray), intent(in) :: ray type(t_triangle), intent(in) :: triangle logical :: intersects real(wp) :: N(3), P(3), C(3), edge(3), vp(3) real(wp) :: area2, d, t, NdotRayDirection intersects = .false. N = triangle%n area2 = sqrt(sum(N(:)*N(:))) NdotRayDirection = sum(N(:)*ray%d(:)) 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 = .true. 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_v Output boundary vertices/normals. !! @param boundary_vertex_count Output total boundary vertex count !! @param boundary_edge_count Output total boundary edge counts pure subroutine f_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 f_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 f_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 f_register_edge(temp_boundary_v, edge, edge_index, edge_count) end do ! Check all edges and count repeated edges 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 edge_occurrence(i) = edge_occurrence(i) + 1 end if end if end do end do ! 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)/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 f_check_boundary !> This procedure appends the edge end vertices to a temporary buffer. !! @param temp_boundary_v Temporary edge end vertex buffer !! @param edge Edges end points to be registered !! @param edge_index Edge index iterator !! @param edge_count Total number of edges pure subroutine f_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 f_register_edge !> This procedure check if interpolates is needed for 2D models. !! @param boundary_v Temporary edge end vertex buffer !! @param boundary_edge_count Output total number of boundary edges !! @param spacing Dimensions of the current levelset cell !! @param interpolate Logical output pure subroutine f_check_interpolation_2D(boundary_v, boundary_edge_count, spacing, interpolate) logical, intent(inout) :: interpolate !< Logical indicator of interpolation integer, intent(in) :: boundary_edge_count !< Number of boundary edges real(wp), intent(in), dimension(1:boundary_edge_count, 1:3, 1:2) :: boundary_v real(wp), dimension(1:3), intent(in) :: spacing real(wp) :: l1, cell_width !< Length of each boundary edge and cell width integer :: j !< Boundary edge index iterator cell_width = minval(spacing(1:2)) interpolate = .false. do j = 1, boundary_edge_count l1 = sqrt((boundary_v(j, 2, 1) - boundary_v(j, 1, 1))**2 + & (boundary_v(j, 2, 2) - boundary_v(j, 1, 2))**2) if ((l1 > cell_width)) then interpolate = .true. else interpolate = .false. end if end do end subroutine f_check_interpolation_2D !> This procedure check if interpolates is needed for 3D models. !! @param model Model to search in. !! @param spacing Dimensions of the current levelset cell !! @param interpolate Logical output pure subroutine f_check_interpolation_3D(model, spacing, interpolate) logical, intent(inout) :: interpolate type(t_model), intent(in) :: model real(wp), dimension(1:3), intent(in) :: spacing real(wp), dimension(1:3) :: edge_l real(wp) :: cell_width real(wp), dimension(1:3, 1:3) :: tri_v integer :: i, j !< Loop iterator cell_width = minval(spacing) interpolate = .false. do i = 1, model%ntrs do j = 1, 3 tri_v(1, j) = model%trs(i)%v(1, j) tri_v(2, j) = model%trs(i)%v(2, j) tri_v(3, j) = model%trs(i)%v(3, j) end do edge_l(1) = sqrt((tri_v(1, 2) - tri_v(1, 1))**2 + & (tri_v(2, 2) - tri_v(2, 1))**2 + & (tri_v(3, 2) - tri_v(3, 1))**2) edge_l(2) = sqrt((tri_v(1, 3) - tri_v(1, 2))**2 + & (tri_v(2, 3) - tri_v(2, 2))**2 + & (tri_v(3, 3) - tri_v(3, 2))**2) edge_l(3) = sqrt((tri_v(1, 1) - tri_v(1, 3))**2 + & (tri_v(2, 1) - tri_v(2, 3))**2 + & (tri_v(3, 1) - tri_v(3, 3))**2) if ((edge_l(1) > cell_width) .or. & (edge_l(2) > cell_width) .or. & (edge_l(3) > cell_width)) then interpolate = .true. else interpolate = .false. end if end do end subroutine f_check_interpolation_3D !> This procedure interpolates 2D models. !! @param boundary_v Group of all the boundary vertices of the 2D model without interpolation !! @param boundary_edge_count Output total number of boundary edges !! @param spacing Dimensions of the current levelset cell !! @param interpolated_boundary_v Output all the boundary vertices of the interpolated 2D model !! @param total_vertices Total number of vertices after interpolation pure subroutine f_interpolate_2D(boundary_v, boundary_edge_count, spacing, interpolated_boundary_v, total_vertices) real(wp), intent(in), dimension(:, :, :) :: boundary_v real(wp), dimension(1:3), intent(in) :: spacing real(wp), allocatable, intent(inout), dimension(:, :) :: interpolated_boundary_v integer, intent(inout) :: total_vertices, boundary_edge_count integer :: num_segments integer :: i, j real(wp) :: edge_length, cell_width real(wp), dimension(1:2) :: edge_x, edge_y, edge_del ! Get the number of boundary edges cell_width = minval(spacing(1:2)) num_segments = 0 ! First pass: Calculate the total number of vertices including interpolated ones total_vertices = 1 do i = 1, boundary_edge_count ! Get the coordinates of the two ends of the current edge edge_x(1) = boundary_v(i, 1, 1) edge_y(1) = boundary_v(i, 1, 2) edge_x(2) = boundary_v(i, 2, 1) edge_y(2) = boundary_v(i, 2, 2) ! Compute the length of the edge edge_length = sqrt((edge_x(2) - edge_x(1))**2 + & (edge_y(2) - edge_y(1))**2) ! Determine the number of segments if (edge_length > cell_width) then num_segments = Ifactor_2D*ceiling(edge_length/cell_width) else num_segments = 1 end if ! Each edge contributes num_segments vertices total_vertices = total_vertices + num_segments end do ! Allocate memory for the new boundary vertices array allocate (interpolated_boundary_v(1:total_vertices, 1:3)) ! Fill the new boundary vertices array with original and interpolated vertices total_vertices = 1 do i = 1, boundary_edge_count ! Get the coordinates of the two ends of the current edge edge_x(1) = boundary_v(i, 1, 1) edge_y(1) = boundary_v(i, 1, 2) edge_x(2) = boundary_v(i, 2, 1) edge_y(2) = boundary_v(i, 2, 2) ! Compute the length of the edge edge_length = sqrt((edge_x(2) - edge_x(1))**2 + & (edge_y(2) - edge_y(1))**2) ! Determine the number of segments and interpolation step if (edge_length > cell_width) then num_segments = Ifactor_2D*ceiling(edge_length/cell_width) edge_del(1) = (edge_x(2) - edge_x(1))/num_segments edge_del(2) = (edge_y(2) - edge_y(1))/num_segments else num_segments = 1 edge_del(1) = 0._wp edge_del(2) = 0._wp end if interpolated_boundary_v(1, 1) = edge_x(1) interpolated_boundary_v(1, 2) = edge_y(1) interpolated_boundary_v(1, 3) = 0._wp ! Add original and interpolated vertices to the output array do j = 1, num_segments - 1 total_vertices = total_vertices + 1 interpolated_boundary_v(total_vertices, 1) = edge_x(1) + j*edge_del(1) interpolated_boundary_v(total_vertices, 2) = edge_y(1) + j*edge_del(2) end do ! Add the last vertex of the edge if (num_segments > 0) then total_vertices = total_vertices + 1 interpolated_boundary_v(total_vertices, 1) = edge_x(2) interpolated_boundary_v(total_vertices, 2) = edge_y(2) end if end do end subroutine f_interpolate_2D !> This procedure interpolates 3D models. !! @param model Model to search in. !! @param spacing Dimensions of the current levelset cell !! @param interpolated_boundary_v Output all the boundary vertices of the interpolated 3D model !! @param total_vertices Total number of vertices after interpolation impure subroutine f_interpolate_3D(model, spacing, interpolated_boundary_v, total_vertices) real(wp), dimension(1:3), intent(in) :: spacing type(t_model), intent(in) :: model real(wp), allocatable, intent(inout), dimension(:, :) :: interpolated_boundary_v integer, intent(out) :: total_vertices integer :: i, j, k, num_triangles, num_segments, num_inner_vertices real(wp), dimension(1:3, 1:3) :: tri real(wp), dimension(1:3) :: edge_del, cell_area real(wp), dimension(1:3) :: bary_coord !< Barycentric coordinates real(wp) :: edge_length, cell_width, cell_area_min, tri_area ! Number of triangles in the model num_triangles = model%ntrs cell_width = minval(spacing) ! Find the minimum surface area cell_area(1) = spacing(1)*spacing(2) cell_area(2) = spacing(1)*spacing(3) cell_area(3) = spacing(2)*spacing(3) cell_area_min = minval(cell_area) num_inner_vertices = 0 ! Calculate the total number of vertices including interpolated ones total_vertices = 0 do i = 1, num_triangles do j = 1, 3 ! Get the coordinates of the two vertices of the current edge tri(1, 1) = model%trs(i)%v(j, 1) tri(1, 2) = model%trs(i)%v(j, 2) tri(1, 3) = model%trs(i)%v(j, 3) ! Next vertex in the triangle (cyclic) tri(2, 1) = model%trs(i)%v(mod(j, 3) + 1, 1) tri(2, 2) = model%trs(i)%v(mod(j, 3) + 1, 2) tri(2, 3) = model%trs(i)%v(mod(j, 3) + 1, 3) ! Compute the length of the edge edge_length = sqrt((tri(2, 1) - tri(1, 1))**2 + & (tri(2, 2) - tri(1, 2))**2 + & (tri(2, 3) - tri(1, 3))**2) ! Determine the number of segments if (edge_length > cell_width) then num_segments = Ifactor_3D*ceiling(edge_length/cell_width) else num_segments = 1 end if ! Each edge contributes num_segments vertices total_vertices = total_vertices + num_segments + 1 end do ! Add vertices inside the triangle do k = 1, 3 tri(k, 1) = model%trs(i)%v(k, 1) tri(k, 2) = model%trs(i)%v(k, 2) tri(k, 3) = model%trs(i)%v(k, 3) end do call f_tri_area(tri, tri_area) if (tri_area > threshold_bary*cell_area_min) then num_inner_vertices = Ifactor_bary_3D*ceiling(tri_area/cell_area_min) total_vertices = total_vertices + num_inner_vertices end if end do ! Allocate memory for the new boundary vertices array allocate (interpolated_boundary_v(1:total_vertices, 1:3)) ! Fill the new boundary vertices array with original and interpolated vertices total_vertices = 0 do i = 1, num_triangles ! Loop through the 3 edges of each triangle do j = 1, 3 ! Get the coordinates of the two vertices of the current edge tri(1, 1) = model%trs(i)%v(j, 1) tri(1, 2) = model%trs(i)%v(j, 2) tri(1, 3) = model%trs(i)%v(j, 3) ! Next vertex in the triangle (cyclic) tri(2, 1) = model%trs(i)%v(mod(j, 3) + 1, 1) tri(2, 2) = model%trs(i)%v(mod(j, 3) + 1, 2) tri(2, 3) = model%trs(i)%v(mod(j, 3) + 1, 3) ! Compute the length of the edge edge_length = sqrt((tri(2, 1) - tri(1, 1))**2 + & (tri(2, 2) - tri(1, 2))**2 + & (tri(2, 3) - tri(1, 3))**2) ! Determine the number of segments and interpolation step if (edge_length > cell_width) then num_segments = Ifactor_3D*ceiling(edge_length/cell_width) edge_del(1) = (tri(2, 1) - tri(1, 1))/num_segments edge_del(2) = (tri(2, 2) - tri(1, 2))/num_segments edge_del(3) = (tri(2, 3) - tri(1, 3))/num_segments else num_segments = 1 edge_del = 0._wp end if ! Add original and interpolated vertices to the output array do k = 0, num_segments - 1 total_vertices = total_vertices + 1 interpolated_boundary_v(total_vertices, 1) = tri(1, 1) + k*edge_del(1) interpolated_boundary_v(total_vertices, 2) = tri(1, 2) + k*edge_del(2) interpolated_boundary_v(total_vertices, 3) = tri(1, 3) + k*edge_del(3) end do ! Add the last vertex of the edge total_vertices = total_vertices + 1 interpolated_boundary_v(total_vertices, 1) = tri(2, 1) interpolated_boundary_v(total_vertices, 2) = tri(2, 2) interpolated_boundary_v(total_vertices, 3) = tri(2, 3) end do ! Interpolate verties that are not on edges do k = 1, 3 tri(k, 1) = model%trs(i)%v(k, 1) tri(k, 2) = model%trs(i)%v(k, 2) tri(k, 3) = model%trs(i)%v(k, 3) end do call f_tri_area(tri, tri_area) if (tri_area > threshold_bary*cell_area_min) then num_inner_vertices = Ifactor_bary_3D*ceiling(tri_area/cell_area_min) !Use barycentric coordinates for randomly distributed points do k = 1, num_inner_vertices call random_number(bary_coord(1)) call random_number(bary_coord(2)) if ((bary_coord(1) + bary_coord(2)) >= 1._wp) then bary_coord(1) = 1._wp - bary_coord(1) bary_coord(2) = 1._wp - bary_coord(2) end if bary_coord(3) = 1._wp - bary_coord(1) - bary_coord(2) total_vertices = total_vertices + 1 interpolated_boundary_v(total_vertices, 1) = dot_product(bary_coord, tri(1:3, 1)) interpolated_boundary_v(total_vertices, 2) = dot_product(bary_coord, tri(1:3, 2)) interpolated_boundary_v(total_vertices, 3) = dot_product(bary_coord, tri(1:3, 3)) end do end if end do end subroutine f_interpolate_3D !> This procedure determines the levelset distance and normals of the 3D models without interpolation. !! @param model Model to search in. !! @param point The cell centers of the current level cell !! @param normals The output levelset normals !! @param distance The output levelset distance pure subroutine f_distance_normals_3D(model, point, normals, distance) type(t_model), intent(IN) :: model real(wp), dimension(1:3), intent(in) :: point real(wp), dimension(1:3), intent(out) :: normals real(wp), intent(out) :: distance real(wp), dimension(1:3, 1:3) :: tri real(wp) :: dist_min, dist_t_min real(wp) :: dist_min_normal, dist_buffer_normal real(wp), dimension(1:3) :: midp !< Centers of the triangle facets real(wp), dimension(1:3) :: dist_buffer !< Distance between the cell center and the vertices integer :: i, j, tri_idx !< Iterator dist_min = 1.e12_wp dist_min_normal = 1.e12_wp distance = 0._wp tri_idx = 0 do i = 1, model%ntrs do j = 1, 3 tri(j, 1) = model%trs(i)%v(j, 1) tri(j, 2) = model%trs(i)%v(j, 2) tri(j, 3) = model%trs(i)%v(j, 3) dist_buffer(j) = sqrt((point(1) - tri(j, 1))**2 + & (point(2) - tri(j, 2))**2 + & (point(3) - tri(j, 3))**2) end do ! Get the surface center of each triangle facet do j = 1, 3 midp(j) = (tri(1, j) + tri(2, j) + tri(3, j))/3 end do dist_t_min = minval(dist_buffer(1:3)) dist_buffer_normal = sqrt((point(1) - midp(1))**2 + & (point(2) - midp(2))**2 + & (point(3) - midp(3))**2) if (dist_t_min < dist_min) then dist_min = dist_t_min end if if (dist_buffer_normal < dist_min_normal) then dist_min_normal = dist_buffer_normal tri_idx = i end if end do normals(1) = model%trs(tri_idx)%n(1) normals(2) = model%trs(tri_idx)%n(2) normals(3) = model%trs(tri_idx)%n(3) distance = dist_min end subroutine f_distance_normals_3D !> This procedure determines the levelset distance of 2D models without interpolation. !! @param boundary_v Group of all the boundary vertices of the 2D model without interpolation !! @param boundary_vertex_count Output the total number of boundary vertices !! @param boundary_edge_count Output the total number of boundary edges !! @param point The cell centers of the current levelset cell !! @param spacing Dimensions of the current levelset cell !! @return Distance which the levelset distance without interpolation pure function f_distance(boundary_v, boundary_edge_count, point) result(distance) integer, intent(in) :: boundary_edge_count real(wp), intent(in), dimension(1:boundary_edge_count, 1:3, 1:2) :: boundary_v real(wp), dimension(1:3), intent(in) :: point integer :: i real(wp) :: dist_buffer1, dist_buffer2 real(wp), dimension(1:boundary_edge_count) :: dist_buffer real(wp) :: distance distance = 0._wp do i = 1, boundary_edge_count dist_buffer1 = sqrt((point(1) - boundary_v(i, 1, 1))**2 + & & (point(2) - boundary_v(i, 1, 2))**2) dist_buffer2 = sqrt((point(1) - boundary_v(i, 2, 1))**2 + & & (point(2) - boundary_v(i, 2, 2))**2) dist_buffer(i) = minval((/dist_buffer1, dist_buffer2/)) end do distance = minval(dist_buffer) end function f_distance !> This procedure determines the levelset normals of 2D models without interpolation. !! @param boundary_v Group of all the boundary vertices of the 2D model without interpolation !! @param boundary_edge_count Output the total number of boundary edges !! @param point The cell centers of the current levelset cell !! @param normals Output levelset normals without interpolation pure subroutine f_normals(boundary_v, boundary_edge_count, point, normals) integer, intent(in) :: boundary_edge_count real(wp), intent(in), dimension(1:boundary_edge_count, 1:3, 1:2) :: boundary_v real(wp), dimension(1:3), intent(in) :: point real(wp), dimension(1:3), intent(out) :: normals integer :: i, idx_buffer real(wp) :: dist_min, dist_buffer real(wp) :: midp(1:3) dist_buffer = 0._wp dist_min = initial_distance_buffer idx_buffer = 0 do i = 1, boundary_edge_count midp(1) = (boundary_v(i, 2, 1) + boundary_v(i, 1, 1))/2 midp(2) = (boundary_v(i, 2, 2) + boundary_v(i, 1, 2))/2 midp(3) = 0._wp dist_buffer = sqrt((point(1) - midp(1))**2 + & & (point(2) - midp(2))**2) if (dist_buffer < dist_min) then dist_min = dist_buffer idx_buffer = i end if end do normals(1) = boundary_v(idx_buffer, 3, 1) normals(2) = boundary_v(idx_buffer, 3, 2) normals(3) = 0._wp end subroutine f_normals !> This procedure calculates the barycentric facet area pure subroutine f_tri_area(tri, tri_area) real(wp), dimension(1:3, 1:3), intent(in) :: tri real(wp), intent(out) :: tri_area real(wp), dimension(1:3) :: AB, AC, cross integer :: i !< Loop iterator do i = 1, 3 AB(i) = tri(2, i) - tri(1, i) AC(i) = tri(3, i) - tri(1, i) end do cross(1) = AB(2)*AC(3) - AB(3)*AC(2) cross(2) = AB(3)*AC(1) - AB(1)*AC(3) cross(3) = AB(1)*AC(2) - AB(2)*AC(1) tri_area = 0.5_wp*sqrt(cross(1)**2 + cross(2)**2 + cross(3)**2) end subroutine f_tri_area !> This procedure determines the levelset of interpolated 2D models. !! @param interpolated_boundary_v Group of all the boundary vertices of the interpolated 2D model !! @param total_vertices Total number of vertices after interpolation !! @param point The cell centers of the current levelset cell !! @return Distance which the levelset distance without interpolation pure function f_interpolated_distance(interpolated_boundary_v, total_vertices, point) result(distance) integer, intent(in) :: total_vertices real(wp), intent(in), dimension(1:total_vertices, 1:3) :: interpolated_boundary_v real(wp), dimension(1:3), intent(in) :: point

[qodo-code-review]

PR Code Suggestions ✨

No code suggestions found for the PR.

[codecov]

Codecov Report

All modified and coverable lines are covered by tests ✅

Project coverage is 43.71%. Comparing base (adcc0dd) to head (b4cc8f8).
Report is 5 commits behind head on master.

Additional details and impacted files

@@           Coverage Diff           @@
##           master     #941   +/-   ##
=======================================
  Coverage   43.71%   43.71%           
=======================================
  Files          68       68           
  Lines       18360    18360           
  Branches     2292     2292           
=======================================
  Hits         8026     8026           
  Misses       8945     8945           
  Partials     1389     1389           

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.