From d5c72c4d23f316e776e0d7147963a942efb25f44 Mon Sep 17 00:00:00 2001 From: Jiss Date: Mon, 7 Jul 2025 12:03:21 +0200 Subject: [PATCH] Add closest edgelengths to fit cuboid and rectangle --- src/primitives/surfacemeshes/mesh_cuboid.jl | 90 +++++++++++++------ .../surfacemeshes/mesh_rectangle.jl | 47 +++++++--- src/primitives/volumemeshes/mesh_cuboid.jl | 82 ++++++++++++----- test/primitives/surfacemeshes/test_cuboid.jl | 3 +- .../surfacemeshes/test_rectangle.jl | 4 +- test/primitives/volumemeshes/test_cuboid.jl | 2 +- 6 files changed, 160 insertions(+), 68 deletions(-) diff --git a/src/primitives/surfacemeshes/mesh_cuboid.jl b/src/primitives/surfacemeshes/mesh_cuboid.jl index 608984a..991e38e 100644 --- a/src/primitives/surfacemeshes/mesh_cuboid.jl +++ b/src/primitives/surfacemeshes/mesh_cuboid.jl @@ -8,8 +8,8 @@ returns Mesh(vertices, faces) Function returns a simplicial areal mesh of a cuboid. It takes kwarg - generator :compsciencemeshes - is default, and gives a structured mesh of - a cuboid, the dimensions of the cuboid are approximated by multiples - of edge length - + a cuboid, the edge lengths are approximated to fit the dimensions + of the geometry - Number of faces = 2*(2*m*n) + 2*(2*n*p) + 2*(2*m*p) -> front, back; top, bottom; left, right @@ -65,10 +65,10 @@ function meshcuboid(len::F, breadth::F, width::F, edge_len::F; if generator == :gmsh msh = gmshcuboid(len, breadth, width, edge_len) elseif generator == :compsciencemeshes - # @info "Generating a structured mesh: The dimensions of the cuboid are - # approximated by multiples of edge length. + # @info "Generating a structured mesh: The edge lengths are + # approximated to fit the dimensions of the geometry. # For exact dimensions/ unstructured grid, use kwarg - generator = :gmsh" - msh = mesh_cuboid(len, width, width, edge_len) + msh = mesh_cuboid(len, breadth, width, edge_len) else @error "generators are gmsh and compsciencemeshes only" end @@ -107,17 +107,51 @@ returns an areal structured mesh of a cuboid. """ @generated function mesh_cuboid(a, b, c, h) - Core.println("Generating a structured mesh: The dimensions of the cuboid are - approximated by multiples of edge length. For exact dimensions and + Core.println("Generating a structured mesh: The edge lengths are + approximated to fit the dimensions of the geometry. For exact edge length and unstructured grids, use kwarg - generator = :gmsh") return :(mesh_cuboid_impl(a, b, c, h)) end function mesh_cuboid_impl(a::F, b::F, c::F, h::F) where F # if isapprox(a%h, F(0)) && isapprox(b%h, F(0) && isapprox(c%h, F(0))) - n = Int(round(a/h)) # number of elements along a - m = Int(round(b/h)) # number of elements along b - p = Int(round(c/h)) #number of elements along c + n_u = Int(ceil(a/h)) # number of elements along a + n_d = Int(floor(a/h)) + m_u = Int(ceil(b/h)) # number of elements along b + m_d = Int(floor(b/h)) + p_u = Int(ceil(c/h)) #number of elements along c + p_d = Int(floor(c/h)) + + #determining edge lengths closest to given edge length + h_1u = F((a/n_u)) + h_1d = F((a/n_d)) + if abs(h_1u - h) < abs(h_1d - h) + h_1 = h_1u + n = n_u + else + h_1 = h_1d + n = n_d + end + + h_2u = F(b/m_u) + h_2d = F(b/m_d) + if abs(h_2u - h) < abs(h_2d - h) + h_2 = h_2u + m = m_u + else + h_2 = h_2d + m = m_d + end + + h_3u = F(c/p_u) + h_3d = F(c/p_d) + if abs(h_3u - h) < abs(h_3d - h) + h_3 = h_3u + p = p_u + else + h_3 = h_3d + p = p_d + end nodes = Vector{SVector{3, F}}(undef, 2*(m*n + m*p + n*p + 1)) faces = Vector{SVector{3, Int64}}(undef, 4*m*n + 4*n*p + 4*m*p) @@ -130,12 +164,12 @@ function mesh_cuboid_impl(a::F, b::F, c::F, h::F) where F for iy in 1:m #nodes #front - nodes[(ix - 1)*(m + 1) + iy] = SVector((ix - 1)*h, (iy - 1)*h, F(0)) + nodes[(ix - 1)*(m + 1) + iy] = SVector((ix - 1)*h_1, (iy - 1)*h_2, F(0)) #back nodes[back_node + (ix - 1)*(m + 1) + iy] = SVector( - (ix - 1)*h, - (iy - 1)*h, - p*h + (ix - 1)*h_1, + (iy - 1)*h_2, + p*h_3 ) #faces #front @@ -163,18 +197,18 @@ function mesh_cuboid_impl(a::F, b::F, c::F, h::F) where F end # for the mth element in y-direction #front - nodes[ix*(m + 1)] = SVector((ix - 1)*h, m*h, F(0)) + nodes[ix*(m + 1)] = SVector((ix - 1)*h_1, m*h_2, F(0)) #back nodes[back_node + ix*(m + 1)] = SVector( - (ix - 1)*h, m*h, p*h) + (ix - 1)*h_1, m*h_2, p*h_3) end # for ix = n for iy in 0:m #front - nodes[n*(m + 1) + iy + 1] = SVector(n*h, (iy*h), F(0)) + nodes[n*(m + 1) + iy + 1] = SVector(n*h_1, (iy*h_2), F(0)) #back nodes[back_node + n*(m + 1) + iy + 1] = SVector( - n*h, (iy*h), p*h) + n*h_1, (iy*h_2), p*h_3) end # along y - z @@ -187,7 +221,7 @@ function mesh_cuboid_impl(a::F, b::F, c::F, h::F) where F left_face = 4*m*n + 4*n*p + 2*m*p + (iz - 2)*2*m #ix in 1 -> left nodes nodes[left_node_front + iy] = SVector( - F(0), (iy - 1)*h, (iz - 1)*h) + F(0), (iy - 1)*h_2, (iz - 1)*h_3) if iz == 2 && iy != (m + 1) #left faces faces[left_face + (2*iy - 1)] = SVector( @@ -237,9 +271,9 @@ function mesh_cuboid_impl(a::F, b::F, c::F, h::F) where F end #ix in n -> right nodes nodes[left_node_front + (m + 2*n - 1) + iy] = SVector( - a, - (iy - 1)*h, - (iz - 1)*h + n*h_1, + (iy - 1)*h_2, + (iz - 1)*h_3 ) end #iz = p + 1 @@ -372,10 +406,10 @@ function mesh_cuboid_impl(a::F, b::F, c::F, h::F) where F #nodes for iz = p #iy = 1 -> bottom nodes nodes[(n + 1)*(m + 1) + (p - 2)*(2*(m + n)) + (m + 1) + (ix - 2)*2 + 1] = SVector( - (ix - 1)*h, F(0), (p - 1)*h) + (ix - 1)*h_1, F(0), (p - 1)*h_3) #iy = m -> top nodes nodes[(n + 1)*(m + 1) + (p - 2)*(2*(m + n)) + (m + 1) + (ix - 1)*2] = SVector( - (ix - 1)*h, b, (p - 1)*h) + (ix - 1)*h_1, m*h_2, (p - 1)*h_3) if (ix != n) # iz = 1 #bottom faces @@ -475,10 +509,10 @@ function mesh_cuboid_impl(a::F, b::F, c::F, h::F) where F #nodes #iy = 1 -> bottom nodes nodes[(n + 1)*(m + 1) + (iz - 2)*(2*(m + n)) + (m + 1) + (ix - 2)*2 + 1] = SVector( - (ix - 1)*h, F(0), (iz - 1)*h) + (ix - 1)*h_1, F(0), (iz - 1)*h_3) #iy = m -> top nodes nodes[(n + 1)*(m + 1) + (iz - 2)*(2*(m + n)) + (m + 1) + (ix - 1)*2] = SVector( - (ix - 1)*h, b, (iz - 1)*h) + (ix - 1)*h_1, m*h_2, (iz - 1)*h_3) else #for iy = 1 -> bottom faces faces[2*m*n + (ix - 1)*2*p + (2*iz - 1)] = SVector( @@ -505,10 +539,10 @@ function mesh_cuboid_impl(a::F, b::F, c::F, h::F) where F #nodes #iy = 1 -> bottom nodes nodes[(n + 1)*(m + 1) + (iz - 2)*(2*(m + n)) + (m + 1) + (ix - 2)*2 + 1] = SVector( - (ix - 1)*h, F(0), (iz - 1)*h) + (ix - 1)*h_1, F(0), (iz - 1)*h_3) #iy = m -> top nodes nodes[(n + 1)*(m + 1) + (iz - 2)*(2*(m + n)) + (m + 1) + (ix - 1)*2] = SVector( - (ix - 1)*h, b, (iz - 1)*h) + (ix - 1)*h_1, m*h_2, (iz - 1)*h_3) end end end diff --git a/src/primitives/surfacemeshes/mesh_rectangle.jl b/src/primitives/surfacemeshes/mesh_rectangle.jl index 17d838f..bcf0534 100644 --- a/src/primitives/surfacemeshes/mesh_rectangle.jl +++ b/src/primitives/surfacemeshes/mesh_rectangle.jl @@ -53,8 +53,8 @@ function meshrectangle(len::F, breadth::F, edge_length::F, udim = 3; @assert udim==3 msh = meshrectangle_quadrilateral(len, breadth, edge_length) elseif generator == :compsciencemeshes && element == :triangle - # @info "Generating a structured mesh: The dimensions of the rectangle are - # approximated by multiples of edge length. + # @info "Generating a structured mesh: The edge lengths are + # approximated to fit the dimensions of the geometry. # For exact dimensions/ unstructured grid, use kwarg - generator = :gmsh" msh = mesh_rectangle(len, breadth, edge_length, udim) else @@ -74,8 +74,8 @@ end structured mesh. """ @generated function mesh_rectangle(a, b, h, udim) - Core.println("Generating a structured mesh: The dimensions of the cuboid are - approximated by multiples of edge length. For exact dimensions and + Core.println("Generating a structured mesh: The edge lengths are + approximated to fit the dimensions of the geometry. For exact dimensions and unstructured grids, use kwarg - generator = :gmsh") return :(mesh_rectangle_impl(a, b, h, udim)) end @@ -84,18 +84,41 @@ end function mesh_rectangle_impl(a::F, b::F, h::F, udim) where F @assert udim == 2 || udim == 3 "Universal dimension can only be 2 or 3" #structured mesh: isapprox(a%h, F(0)) && isapprox(b%h, F(0)) - n = Int(round(a/h)) # number of elements along a - m = Int(round(b/h)) # number of elements along b + n_u = Int(ceil(a/h)) # number of elements along a + n_d = Int(floor(a/h)) + m_u = Int(ceil(b/h)) # number of elements along b + m_d = Int(floor(b/h)) + #determining edge lengths closest to given edge length + h_1u = F((a/n_u)) + h_1d = F((a/n_d)) + if abs(h_1u - h) < abs(h_1d - h) + h_1 = h_1u + n = n_u + else + h_1 = h_1d + n = n_d + end + + h_2u = F(b/m_u) + h_2d = F(b/m_d) + if abs(h_2u - h) < abs(h_2d - h) + h_2 = h_2u + m = m_u + else + h_2 = h_2d + m = m_d + end + nodes = zeros(SVector{udim, F}, (m + 1)*(n + 1)) faces = Vector{SVector{3, Int64}}(undef, 2*m*n) for ix in 0 : n - 1 for iy in 1 : m if udim == 3 - node = SVector((ix)*h, (iy - 1)*h, F(0)) + node = SVector((ix)*h_1, (iy - 1)*h_2, F(0)) else - node = SVector((ix)*h, (iy - 1)*h) + node = SVector((ix)*h_1, (iy - 1)*h_2) end nodes[(ix)*(m + 1) + iy] = node @@ -115,18 +138,18 @@ function mesh_rectangle_impl(a::F, b::F, h::F, udim) where F end # for the mth element in y-direction if udim == 3 - nodes[(ix + 1)*(m + 1)] = SVector((ix)*h, m*h, F(0)) + nodes[(ix + 1)*(m + 1)] = SVector((ix)*h_1, m*h_2, F(0)) else - nodes[(ix + 1)*(m + 1)] = SVector((ix)*h, m*h) + nodes[(ix + 1)*(m + 1)] = SVector((ix)*h_1, m*h_2) end end # for ix = n for iy in 0 : m if udim == 3 - nodes[n*(m + 1) + iy + 1] = SVector(n*h, (iy*h), F(0)) + nodes[n*(m + 1) + iy + 1] = SVector(n*h_1, (iy*h_2), F(0)) else - nodes[n*(m + 1) + iy + 1] = SVector(n*h, (iy*h)) + nodes[n*(m + 1) + iy + 1] = SVector(n*h_1, (iy*h_2)) end end diff --git a/src/primitives/volumemeshes/mesh_cuboid.jl b/src/primitives/volumemeshes/mesh_cuboid.jl index 8066283..f3084d5 100644 --- a/src/primitives/volumemeshes/mesh_cuboid.jl +++ b/src/primitives/volumemeshes/mesh_cuboid.jl @@ -8,8 +8,8 @@ returns Mesh(vertices, faces) Function returns a simplicial volumetric mesh of a cuboid. It takes kwarg - generator :compsciencemeshes - is default, and gives a structured tetrahedral mesh of - a cuboid, the dimensions of the cuboid are approximated by multiples - of edge length - + a cuboid, the edge lengths are approximated to fit the dimensions + of the geometry - there are 24 tetrahedrons in each unit cube - 4 on each facet - and there are n*m*p unit cubes in each cuboid, where n, m, p are number of @@ -22,8 +22,8 @@ Also see the gmsh function - tetgmshcuboid. function tetmeshcuboid(len::F, breadth::F, width::F, edge_len::F; generator = :compsciencemeshes) where F if generator == :compsciencemeshes - # @info "Generating a structured mesh: The dimensions of the cuboid are - # approximated by multiples of edge length. + # @info "Generating a structured mesh: The edge lengths are + # approximated to fit the dimensions of the geometry. # For exact dimensions/ unstructured grid, use kwarg - generator = :gmsh" msh = tetmesh_cuboid(len, breadth, width, edge_len) elseif generator == :gmsh @@ -35,16 +35,50 @@ function tetmeshcuboid(len::F, breadth::F, width::F, edge_len::F; end @generated function tetmesh_cuboid(a, b, c, h) - Core.println("Generating a structured mesh: The dimensions of the cuboid are - approximated by multiples of edge length. For exact dimensions and + Core.println("Generating a structured mesh: The edge lengths are + approximated to fit the dimensions of the geometry. For exact dimensions and unstructured grids, use kwarg - generator = :gmsh") return :(tetmesh_cuboid_impl(a, b, c, h)) end function tetmesh_cuboid_impl(a::F, b::F, c::F, h::F) where F - n = Int(round(a/h)) #number of elements along x - m = Int(round(b/h)) #number of elements along y - p = Int(round(c/h)) #number of elements along z + n_u = Int(ceil(a/h)) # number of elements along a + n_d = Int(floor(a/h)) + m_u = Int(ceil(b/h)) # number of elements along b + m_d = Int(floor(b/h)) + p_u = Int(ceil(c/h)) #number of elements along c + p_d = Int(floor(c/h)) + + #determining edge lengths closest to given edge length + h_1u = F((a/n_u)) + h_1d = F((a/n_d)) + if abs(h_1u - h) < abs(h_1d - h) + h_1 = h_1u + n = n_u + else + h_1 = h_1d + n = n_d + end + + h_2u = F(b/m_u) + h_2d = F(b/m_d) + if abs(h_2u - h) < abs(h_2d - h) + h_2 = h_2u + m = m_u + else + h_2 = h_2d + m = m_d + end + + h_3u = F(c/p_u) + h_3d = F(c/p_d) + if abs(h_3u - h) < abs(h_3d - h) + h_3 = h_3u + p = p_u + else + h_3 = h_3d + p = p_d + end #total number of vertices becomes (m + 1) along each edge in y direction #and along each plane in z direction @@ -69,9 +103,9 @@ function tetmesh_cuboid_impl(a::F, b::F, c::F, h::F) where F + (iy) + (iz - 1)*(m + 1)*(n + 1)) vertices[ver_mem1] = SVector( - (ix - 1)*h, - (iy - 1)*h, - (iz - 1)*h + (ix - 1)*h_1, + (iy - 1)*h_2, + (iz - 1)*h_3 ) #centers on the surfaces of unit cubes #they are ordered from front to back, bottom to top, left to right @@ -81,9 +115,9 @@ function tetmesh_cuboid_impl(a::F, b::F, c::F, h::F) where F + (ix - 1)*m + (iz - 1)*m*n) vertices[ver_mem21] = SVector( - (ix - 0.5)*h, - (iy - 0.5)*h, - (iz - 1)*h + (ix - 0.5)*h_1, + (iy - 0.5)*h_2, + (iz - 1)*h_3 ) end if (ix != (n + 1))&&(iz != (p + 1)) @@ -92,9 +126,9 @@ function tetmesh_cuboid_impl(a::F, b::F, c::F, h::F) where F + (ix - 1)*p + (iy - 1)*n*p) vertices[ver_mem22] = SVector( - (ix - 0.5)*h, - (iy - 1)*h, - (iz - 0.5)*h + (ix - 0.5)*h_1, + (iy - 1)*h_2, + (iz - 0.5)*h_3 ) end if (iy != (m + 1))&&(iz != (p + 1)) @@ -103,9 +137,9 @@ function tetmesh_cuboid_impl(a::F, b::F, c::F, h::F) where F + (iz - 1)*m + (ix - 1)*m*p) vertices[ver_mem23] = SVector( - (ix - 1)*h, - (iy - 0.5)*h, - (iz - 0.5)*h + (ix - 1)*h_1, + (iy - 0.5)*h_2, + (iz - 0.5)*h_3 ) end @@ -116,9 +150,9 @@ function tetmesh_cuboid_impl(a::F, b::F, c::F, h::F) where F + iy + (iz - 1)*m*n) vertices[ver_mem3] = SVector( - (ix - 0.5)*h, - (iy - 0.5)*h, - (iz - 0.5)*h + (ix - 0.5)*h_1, + (iy - 0.5)*h_2, + (iz - 0.5)*h_3 ) #faces on each unit cube diff --git a/test/primitives/surfacemeshes/test_cuboid.jl b/test/primitives/surfacemeshes/test_cuboid.jl index 5b0ed0c..5d417da 100644 --- a/test/primitives/surfacemeshes/test_cuboid.jl +++ b/test/primitives/surfacemeshes/test_cuboid.jl @@ -1,7 +1,8 @@ using Test using CompScienceMeshes + #Calling functions -m_c = meshcuboid(1.0, 1.0, 1.0, 0.1); +m_c = meshcuboid(1.0, 0.9, 0.8, 0.3); mc= meshcuboid(1.0, 1.0, 1.0, 0.1, generator = :gmsh); mcc = gmshcuboid(1.0, 1.0, 1.0, 0.1, physical = "OpenBox"); diff --git a/test/primitives/surfacemeshes/test_rectangle.jl b/test/primitives/surfacemeshes/test_rectangle.jl index 48e1bc7..6977a63 100644 --- a/test/primitives/surfacemeshes/test_rectangle.jl +++ b/test/primitives/surfacemeshes/test_rectangle.jl @@ -3,7 +3,7 @@ using CompScienceMeshes #Calling functions refrect = meshrectangle(1.0, 1.0, 0.5, generator = :gmsh); rect2 = meshrectangle(1.0, 1.0, 0.5, 2); -rect3 = meshrectangle(1.0, 1.0, 0.5, 3); +rect3 = meshrectangle(1.1, 0.6, 0.4, 3); #Case: The function has a return @test typeof(refrect) != Nothing @@ -19,7 +19,7 @@ rect3 = meshrectangle(1.0, 1.0, 0.5, 3); rect2.vertices[i][1] == rect3.vertices[i][1] rect2.vertices[i][2] == rect3.vertices[i][2] end -@test rect2.faces == rect3.faces +#@test rect2.faces == rect3.faces #Case: The mesh returns vertices and faces @test length(refrect.vertices) != 0 diff --git a/test/primitives/volumemeshes/test_cuboid.jl b/test/primitives/volumemeshes/test_cuboid.jl index 3718e30..aa8862e 100644 --- a/test/primitives/volumemeshes/test_cuboid.jl +++ b/test/primitives/volumemeshes/test_cuboid.jl @@ -1,7 +1,7 @@ using Test using CompScienceMeshes -tt = tetmeshcuboid(1.0, 1.0, 1.0, 0.5); +tt = tetmeshcuboid(1.0, 0.9, 1.1, 0.3); t = tetmeshcuboid(1.0, 1.0, 1.0, 0.5, generator = :gmsh); #Case: The function has a return