diff --git a/CHANGELOG.md b/CHANGELOG.md index 13d0f86..75cb4ca 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,6 +2,9 @@ ## [Unreleased] +- Implement polyhedron types (`hex8`, `tet4`, etc) +- Removed `wrapindex` function in favour of `mod1` +- Added handling of boundary edges to `subquad` and `subtri` - Moved excess data from assets folder to dedicated repository diff --git a/Manifest.toml b/Manifest.toml index db412b4..62b174b 100644 --- a/Manifest.toml +++ b/Manifest.toml @@ -2,7 +2,7 @@ julia_version = "1.10.2" manifest_format = "2.0" -project_hash = "8f980f252f20c77185d25ab017fcbeabc72c0bbc" +project_hash = "04794bae62f418ff0e74eb0067582ab95ddb0fe4" [[deps.AbstractFFTs]] deps = ["LinearAlgebra"] diff --git a/Project.toml b/Project.toml index 0e9fd82..56efb65 100644 --- a/Project.toml +++ b/Project.toml @@ -14,6 +14,7 @@ Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc" Rotations = "6038ab10-8711-5258-84ad-4b1120ba62dc" +StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" XML = "72c71f33-b9b6-44de-8c94-c961784809e2" diff --git a/examples/demo_rhombicdodecahedron.jl b/examples/demo_rhombicdodecahedron.jl new file mode 100644 index 0000000..18c9464 --- /dev/null +++ b/examples/demo_rhombicdodecahedron.jl @@ -0,0 +1,21 @@ +using Comodo +using GLMakie +using GeometryBasics +using Rotations + +# Creating faces and vertices for a rhombic dodecahedron +F,V = rhombicdodecahedron(5.0) + +## Visualize mesh +markersize = 25 +strokewidth = 2 +strokecolor = :black +fig = Figure(size = (800,800)) + +ax1 = Axis3(fig[1, 1], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "Rhombic dodecahedron") + +hp1 = poly!(ax1, GeometryBasics.Mesh(V,F), color=:white,transparency=false,strokewidth=strokewidth,strokecolor=strokecolor,shading = FastShading) +hp2 = scatter!(ax1, V,markersize=markersize,color=:red) +hp3 = normalplot(ax1,F,V; color = :green) + +fig \ No newline at end of file diff --git a/examples/demo_subhex.jl b/examples/demo_subhex.jl new file mode 100644 index 0000000..148f3ff --- /dev/null +++ b/examples/demo_subhex.jl @@ -0,0 +1,53 @@ +using Comodo +using GLMakie +using GeometryBasics + +#= +This demo shows the use of `subhex` to refine a a hexhedral mesh through splitting. +=# + +boxDim = [2.0,3.0,1.0] # Dimensions for the box in each direction +boxEl = [2,3,1] # Number of elements to use in each direction +E,V,F,Fb,CFb_type = hexbox(boxDim,boxEl) + +Eh0,Vh0 = subhex(E,V,1;direction=0) +Fh0 = element2faces(Eh0) + +Eh1,Vh1 = subhex(E,V,1;direction=1) +Fh1 = element2faces(Eh1) + +Eh2,Vh2 = subhex(E,V,1;direction=2) +Fh2 = element2faces(Eh2) + +Eh3,Vh3 = subhex(E,V,1;direction=3) +Fh3 = element2faces(Eh3) + + +# Visualisation +strokewidth = 3 +linewidth = 4 + +fig = Figure(size=(1600,800)) + +ax0 = Axis3(fig[1, 1], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "Refined hexahedral mesh, direction=0 (all)") +hp1 = poly!(ax0,GeometryBasics.Mesh(Vh0,Fh0), strokewidth=strokewidth,shading=FastShading,strokecolor=:black, color=:white, transparency=false, overdraw=false) +hp2 = wireframe!(ax0,GeometryBasics.Mesh(V,F), linewidth=linewidth,color=:red, overdraw=false) +# hp3 = normalplot(ax0,Fh0,Vh0;color=:black) + +ax1 = Axis3(fig[1, 2], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "Refined hexahedral mesh, direction=1") +hp1 = poly!(ax1,GeometryBasics.Mesh(Vh1,Fh1), strokewidth=strokewidth,shading=FastShading,strokecolor=:black, color=:white, transparency=false, overdraw=false) +hp2 = wireframe!(ax1,GeometryBasics.Mesh(V,F), linewidth=linewidth,color=:red, overdraw=false) +# hp3 = normalplot(ax1,Fh1,Vh1;color=:black) + +ax2 = Axis3(fig[2, 1], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "Refined hexahedral mesh, direction=2") +hp1 = poly!(ax2,GeometryBasics.Mesh(Vh2,Fh2), strokewidth=strokewidth,shading=FastShading,strokecolor=:black, color=:white, transparency=false, overdraw=false) +hp2 = wireframe!(ax2,GeometryBasics.Mesh(V,F), linewidth=linewidth,color=:red, overdraw=false) +# hp3 = normalplot(ax2,Fh2,Vh2;color=:black) + +ax3 = Axis3(fig[2, 2], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "Refined hexahedral mesh, direction=3") +hp1 = poly!(ax3,GeometryBasics.Mesh(Vh3,Fh3), strokewidth=strokewidth,shading=FastShading,strokecolor=:black, color=:white, transparency=false, overdraw=false) +hp2 = wireframe!(ax3,GeometryBasics.Mesh(V,F), linewidth=linewidth,color=:red, overdraw=false) +# hp3 = normalplot(ax3,Fh3,Vh3;color=:black) + +fig + diff --git a/examples/demo_subtri.jl b/examples/demo_subtri.jl index 393fc28..043d87c 100644 --- a/examples/demo_subtri.jl +++ b/examples/demo_subtri.jl @@ -1,6 +1,7 @@ using Comodo using GLMakie using GeometryBasics +using LinearAlgebra #= This demo shows the use of `subtri` to refine triangulated meshes. Each diff --git a/examples/demo_surface_mesh_smoothing.jl b/examples/demo_surface_mesh_smoothing.jl index 4ea96be..4c93a5e 100644 --- a/examples/demo_surface_mesh_smoothing.jl +++ b/examples/demo_surface_mesh_smoothing.jl @@ -10,6 +10,7 @@ if testCase == 1 # Triangle mesh bunny F = faces(M) V = coordinates(M) V = topoints(V) + F = tofaces(F) F,V,ind1,ind2 = mergevertices(F,V; roundVertices=true) elseif testCase == 2 # Quad mesh sphere F,V = quadsphere(4,100) diff --git a/examples/demo_tet2hex.jl b/examples/demo_tet2hex.jl new file mode 100644 index 0000000..24c541e --- /dev/null +++ b/examples/demo_tet2hex.jl @@ -0,0 +1,53 @@ +using Comodo +using GLMakie +using GeometryBasics + +#= +This demo shows the use of `tet2hex` to convert tetrahedral elements to hexahedral elements +domain. +=# + +testCase = 1 + +if testCase == 1 + E = [tet4{Int64}(1,2,3,4),tet4{Int64}(2,3,4,5),tet4{Int64}(6,7,8,9)] + V = [Point{3,Float64}(-1.0,0.0,0.0), + Point{3,Float64}( 1.0,0.0,0.0), + Point{3,Float64}( 0.0,1.0,0.0), + Point{3,Float64}( 0.0,0.5,1.0), + Point{3,Float64}( 1.0,1.0,1.0), + Point{3,Float64}( 2.0,0.0,0.0), + Point{3,Float64}( 4.0,0.0,0.0), + Point{3,Float64}( 3.0,1.0,0.0), + Point{3,Float64}( 3.0,0.5,1.0), + ] +elseif testCase == 2 + E = [hex8{Int64}(1,2,3,4,5,6,7,8)] + V = [Point{3,Float64}(0.0,0.0,0.0), + Point{3,Float64}(1.0,0.0,0.0), + Point{3,Float64}(1.0,1.0,0.0), + Point{3,Float64}(0.0,1.0,0.0), + Point{3,Float64}(0.0,0.0,1.0), + Point{3,Float64}(1.0,0.0,1.0), + Point{3,Float64}(1.0,1.0,1.0), + Point{3,Float64}(0.0,1.0,1.0), + ] +end +F = element2faces(E) + +Eh,Vh = tet2hex(E,V) +Fh = element2faces(Eh) + +# Visualisation +markersize = 25 + +fig = Figure(size=(1600,800)) + +ax1 = Axis3(fig[1, 1], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "Wireframe of hexahedral mesh") +hp1 = poly!(ax1,GeometryBasics.Mesh(V,F), strokewidth=3,shading=FastShading,strokecolor=:black, color=:white, transparency=true, overdraw=false) +hp2 = normalplot(ax1,F,V) + +hp1 = poly!(ax1,GeometryBasics.Mesh(Vh,Fh), strokewidth=2,shading=FastShading,strokecolor=:red, color=:red, transparency=true, overdraw=false) +hp2 = normalplot(ax1,Fh,Vh;color=:red) + +fig diff --git a/src/Comodo.jl b/src/Comodo.jl index 4d957b2..8e41f11 100644 --- a/src/Comodo.jl +++ b/src/Comodo.jl @@ -11,10 +11,20 @@ import Interpolations # E.g. for resampling curves import BSplineKit import QuadGK import Distances -import DelaunayTriangulation +import DelaunayTriangulation # For regiontrimesh +import StaticArrays # For volumetric mesh definitions include("functions.jl") +# Export imported packages +export GeometryBasics + +# Export types +export Nhedron, tet4, tet10, hex8, hex20, penta6 # Volumetric elements (polyhedra) + +# Export types/structs +export ConnectivitySet + # Export functions export comododir, elements2indices, hexbox, mindist, dist, lerp export gridpoints, gridpoints_equilateral, interp_biharmonic_spline, interp_biharmonic, nbezier, gunique @@ -29,20 +39,15 @@ export con_edge_face, con_edge_edge, con_face_edge, con_face_face, con_face_face export con_vertex_face, con_vertex_vertex, con_vertex_simplex, meshconnectivity export simplexcenter, vertex2simplexdata, simplex2vertexdata export circlepoints, loftlinear, trisurfslice -export wrapindex, edgeangles, count_edge_face, boundaryedges, edges2curve +export edgeangles, count_edge_face, boundaryedges, edges2curve export pointspacingmean, extrudecurve, meshgroup, ray_triangle_intersect export distmarch, mesh_curvature_polynomial, curve_length, evenly_sample, evenly_space export invert_faces, kabsch_rot, sweeploft, loftpoints2surf, revolvecurve export batman, tridisc, triplate, regiontrimesh, scalesimplex, subcurve, dualclad +export tet2hex, element2faces, subhex, rhombicdodecahedron # Export functions: Visualization related export slidercontrol, slider2anim, dirplot, normalplot -# Export types/structs -export ConnectivitySet - -# Export imported packages -export GeometryBasics - end # module diff --git a/src/functions.jl b/src/functions.jl index 849fdf5..e587c30 100644 --- a/src/functions.jl +++ b/src/functions.jl @@ -11,6 +11,16 @@ using BSplineKit # E.g. for resampling curves using QuadGK: quadgk # For numerical integration using Distances using DelaunayTriangulation # For triangular meshing +using StaticArrays + +# Define types +abstract type AbstractPolyhedron{N,T} <: StaticVector{N,T} end +abstract type AbstractNhedron{N,T} <: AbstractPolyhedron{N,T} end +GeometryBasics.@fixed_vector Nhedron = AbstractNhedron +const tet4{T} = Nhedron{4,T} where T<:Integer +const tet10{T} = Nhedron{10,T} where T<:Integer +const hex8{T} = Nhedron{8,T} where T<:Integer +const penta6{T} = Nhedron{6,T} where T<:Integer """ ConnectivitySet(E_uni, con_E2F, con_E2E, F, con_F2E, con_F2F, con_V2E, con_V2F, con_V2V, con_V2V_f, con_F2F_v) @@ -34,6 +44,8 @@ struct ConnectivitySet ConnectivitySet(E_uni, con_E2F, con_E2E, F, con_F2E, con_F2F, con_V2E, con_V2F, con_V2V, con_V2V_f, con_F2F_v) = new(E_uni, con_E2F, con_E2E, F, con_F2E, con_F2F, con_V2E, con_V2F, con_V2V, con_V2V_f, con_F2F_v) end + + """ comododir() @@ -898,9 +910,6 @@ function unique_simplices(F,V=nothing) virtualFaceIndices = sub2ind(n.*ones(Int64,length(F[1])),sort.(F)) _, ind1, ind2 = unique_dict(virtualFaceIndices) - # Fn, ind1, ind2 = gunique(F; return_unique=true, return_index=true, return_inverse=true, return_counts=false, sort_entries=true) - # return Fn, ind1, ind2 - return F[ind1], ind1, ind2 end @@ -1016,16 +1025,10 @@ GeometryBasics.Mesh. The convention is such that for a face referring to the nodes 1-2-3-4, the edges are 1-2, 2-3, 3-4, 4-1. """ function meshedges(F::Array{NgonFace{N,T},1}; unique_only=false) where N where T<:Integer - E = LineFace{Int64}[] - for j1 in 1:N # Loop over each node/point for the current simplex - if j1 cut nBelow = sum(lf) # Number of nodes below if isone(nBelow) # 2-above, 1 below - indP = f[wrapindex(findfirst(lf) .+ (0:2),3)] + indP = f[mod1.(findfirst(lf) .+ (0:2),3)] e1 = sort(indP[[1,2]]) if !haskey(D,e1) @@ -2856,7 +2833,7 @@ function trisurfslice(F::Vector{TriangleFace{TF}},V::Vector{Point{ND,TV}},n = Ve end else # 1-above, 2 below - indP = f[wrapindex(findfirst(.!lf) .+ (0:2),3)] + indP = f[mod1.(findfirst(.!lf) .+ (0:2),3)] e1 = sort(indP[[1,2]]) if !haskey(D,e1) @@ -3915,4 +3892,317 @@ function dualclad(F::Vector{NgonFace{N, TF}},V::Vector{Point{ND,TV}},s; connecti end return Fs,Fq,Vs +end + + +""" + tet2hex(E::Vector{tet4{T}},V::Vector{Point{ND,TV}}) where T<:Integer where ND where TV<:Real + +Converts tetrahedra to hexahedra + +# Description +This function converts the input tetrahedra defined by the element set `E` and the +vertex set `V` to a set of hexahedral elements `Eh` with vertices `Vh`. The conversion +involves a splitting of each tetrahedron into 4 hexahedra. +""" +function tet2hex(E::Vector{tet4{T}},V::Vector{Point{ND,TV}}) where T<:Integer where ND where TV<:Real + # Non-unique tet element faces + Ft = element2faces(E) + + # Unique faces and inverse mapping + Ftu,_,inv_Ft = gunique(Ft; return_unique=true, return_index=true, return_inverse=true, sort_entries=true) + + # Non-unique structured (element by element) tet element edges + Et = Vector{LineFace{T}}(undef,length(E)*6) + for (i,e) in enumerate(E) + ii = 1 + (i-1)*6 + Et[ii ] = LineFace{T}(e[1],e[2]) + Et[ii+1] = LineFace{T}(e[2],e[3]) + Et[ii+2] = LineFace{T}(e[3],e[1]) + Et[ii+3] = LineFace{T}(e[1],e[4]) + Et[ii+4] = LineFace{T}(e[2],e[4]) + Et[ii+5] = LineFace{T}(e[3],e[4]) + end + + #Unique edges and inverse mapping + Etu,_,inv_Et = gunique(Et; return_unique=true, return_index=true, return_inverse=true, sort_entries=true) + + # Create new coordinates + Vc = simplexcenter(E,V) # Element centres + Vf = simplexcenter(Ftu,V) # Face centres + Ve = simplexcenter(Etu,V) # Edge centres + Vh = [V;Vc;Vf;Ve] # Append vertices + + # Create the 4 corner hexahedra for each tetrahedron + offset_Vc = length(V) + offset_F = offset_Vc + length(Vc) + offset_E = offset_F + length(Vf) + inv_Ft = inv_Ft .+ offset_F + inv_Et = inv_Et .+ offset_E + Eh = Vector{hex8{T}}(undef,length(E)*4) # Allocate hexahedral element vector, one hex for each of the 4 nodes per element + for i = eachindex(E) + i_e = 1 + (i-1)*6 # index of first edge + i_f = 1 + (i-1)*4 # index of first face (which happens to also be the first hex element for tetrahedra) + i_vc = i+offset_Vc + Eh[i_f ] = hex8{T}( E[i][1], inv_Et[i_e ], inv_Ft[i_f ], inv_Et[i_e+2], + inv_Et[i_e+3], inv_Ft[i_f+1], i_vc, inv_Ft[i_f+3]) + Eh[i_f+1] = hex8{T}( E[i][2], inv_Et[i_e+1], inv_Ft[i_f ], inv_Et[i_e ], + inv_Et[i_e+4], inv_Ft[i_f+2], i_vc, inv_Ft[i_f+1]) + Eh[i_f+2] = hex8{T}( E[i][3], inv_Et[i_e+2], inv_Ft[i_f ], inv_Et[i_e+1], + inv_Et[i_e+5], inv_Ft[i_f+3], i_vc, inv_Ft[i_f+2]) + Eh[i_f+3] = hex8{T}( E[i][4], inv_Et[i_e+4], inv_Ft[i_f+1], inv_Et[i_e+3], + inv_Et[i_e+5], inv_Ft[i_f+2], i_vc, inv_Ft[i_f+3]) + end + return Eh,Vh +end + +""" + element2faces(E::Vector{Nhedron{N,T}}) where N where T + +Returns element faces + +# Description +This function computes the faces for the input elements defined by `E`. The elements +should be Vectors consisting of `tet4`, `hex8` elements. +""" +function element2faces(E::Vector{Nhedron{N,T}}) where N where T + if eltype(E) <: tet4{T} + nf = 4 + F = Vector{TriangleFace{T}}(undef,length(E)*nf) + for i in eachindex(E) + ii = 1 + (i-1)*nf + F[ii ] = TriangleFace{T}(E[i][3],E[i][2],E[i][1]) + F[ii+1] = TriangleFace{T}(E[i][1],E[i][2],E[i][4]) + F[ii+2] = TriangleFace{T}(E[i][2],E[i][3],E[i][4]) + F[ii+3] = TriangleFace{T}(E[i][3],E[i][1],E[i][4]) + end + elseif eltype(E) <: hex8{T} + nf = 6 + F = Vector{QuadFace{T}}(undef,length(E)*nf) + for i in eachindex(E) + e = E[i] + ii = 1 + (i-1)*nf + F[ii ] = QuadFace{T}(e[4],e[3],e[2],e[1]) # Top + F[ii+1] = QuadFace{T}(e[5],e[6],e[7],e[8]) # Bottom + F[ii+2] = QuadFace{T}(e[1],e[2],e[6],e[5]) # Side 1 + F[ii+3] = QuadFace{T}(e[3],e[4],e[8],e[7]) # Side 2 + F[ii+4] = QuadFace{T}(e[2],e[3],e[7],e[6]) # Front + F[ii+5] = QuadFace{T}(e[5],e[8],e[4],e[1]) # Back + end + end + return F +end + +""" + subhex(E::Vector{hex8{T}},V::Vector{Point{ND,TV}},n::Int64; direction=0) where T<:Integer where ND where TV<:Real + +Split hexahedral elements + +# Description +This function splits the hexahedral elements defined by the elements `E` and +vertices `V`. Splitting is done `n` times as requested. By default the splitting +occurs in all direction (corresponding to the default `direction=0`). If instead +`direction` is set to 1, 2, or 3, then the splitting only occur in the first, second +or third local element direction respectively. Note that this direction depends +on node order used. For a hexahedron where by nodes 1:4 are for the bottom, and +nodes 5:8 are for the top of the element then the directions 1, 2, and 3 correspond +to the x-, y-, and z-direction respectively. +""" +function subhex(E::Vector{hex8{T}},V::Vector{Point{ND,TV}},n::Int64; direction=0) where T<:Integer where ND where TV<:Real + if iszero(n) + return E,V + elseif isone(n) + + # Non-unique hex element faces + Ft = element2faces(E) + + # Non-unique structured (element by element) hexbox element edges + if iszero(direction) + Et = Vector{LineFace{T}}(undef,length(E)*12) + for (i,e) in enumerate(E) + ii = 1 + (i-1)*12 + Et[ii ] = LineFace{T}(e[1],e[2]) + Et[ii+1 ] = LineFace{T}(e[2],e[3]) + Et[ii+2 ] = LineFace{T}(e[3],e[4]) + Et[ii+3 ] = LineFace{T}(e[4],e[1]) + Et[ii+4 ] = LineFace{T}(e[5],e[6]) + Et[ii+5 ] = LineFace{T}(e[6],e[7]) + Et[ii+6 ] = LineFace{T}(e[7],e[8]) + Et[ii+7 ] = LineFace{T}(e[8],e[5]) + Et[ii+8 ] = LineFace{T}(e[1],e[5]) + Et[ii+9 ] = LineFace{T}(e[2],e[6]) + Et[ii+10] = LineFace{T}(e[3],e[7]) + Et[ii+11] = LineFace{T}(e[4],e[8]) + end + elseif isone(direction) # Split in 1st-direction + Et = Vector{LineFace{T}}(undef,length(E)*4) + for (i,e) in enumerate(E) + ii = 1 + (i-1)*4 + Et[ii ] = LineFace{T}(e[1],e[2]) + Et[ii+1 ] = LineFace{T}(e[3],e[4]) + Et[ii+2 ] = LineFace{T}(e[5],e[6]) + Et[ii+3 ] = LineFace{T}(e[7],e[8]) + end + elseif direction==2 # Split in 2nd-direction + Et = Vector{LineFace{T}}(undef,length(E)*4) + for (i,e) in enumerate(E) + ii = 1 + (i-1)*4 + Et[ii ] = LineFace{T}(e[2],e[3]) + Et[ii+1 ] = LineFace{T}(e[4],e[1]) + Et[ii+2 ] = LineFace{T}(e[6],e[7]) + Et[ii+3 ] = LineFace{T}(e[8],e[5]) + end + elseif direction==3 # Split in 3rd-direction + Et = Vector{LineFace{T}}(undef,length(E)*4) + for (i,e) in enumerate(E) + ii = 1 + (i-1)*4 + Et[ii ] = LineFace{T}(e[1],e[5]) + Et[ii+1 ] = LineFace{T}(e[2],e[6]) + Et[ii+2 ] = LineFace{T}(e[3],e[7]) + Et[ii+3 ] = LineFace{T}(e[4],e[8]) + end + end + #Unique edges and inverse mapping + Etu,_,inv_Et = gunique(Et; return_unique=true, return_index=true, return_inverse=true, sort_entries=true) + + # Create mid-edge coordinates + Ve = simplexcenter(Etu,V) # Edge centres + + if iszero(direction) + # Unique faces and inverse mapping + Ftu,_,inv_Ft = gunique(Ft; return_unique=true, return_index=true, return_inverse=true, sort_entries=true) + Vf = simplexcenter(Ftu,V) # Face centres + Vc = simplexcenter(E,V) # Element centres + Vh = [V;Vc;Vf;Ve] # Append vertices + + # Create the 8 corner hexahedra for each hexahedron + offset_Vc = length(V) + offset_F = offset_Vc + length(Vc) + offset_E = offset_F + length(Vf) + inv_Ft = inv_Ft .+ offset_F + inv_Et = inv_Et .+ offset_E + Eh = Vector{hex8{T}}(undef,length(E)*8) # Allocate hexahedral element vector, one hex for each of the 4 nodes per element + for i = eachindex(E) + i_e = 1 + (i-1)*12 # index of first edge + i_f = 1 + (i-1)*6 # index of first face + i_vc = i + offset_Vc + ii = 1 + (i-1)*8 + Eh[ii ] = hex8{T}( E[i][1], inv_Et[i_e ], inv_Ft[i_f ], inv_Et[i_e+3], + inv_Et[i_e+8 ], inv_Ft[i_f+2], i_vc, inv_Ft[i_f+5] ) + + Eh[ii+1] = hex8{T}( E[i][2], inv_Et[i_e+1], inv_Ft[i_f ], inv_Et[i_e], + inv_Et[i_e+9 ], inv_Ft[i_f+4], i_vc, inv_Ft[i_f+2] ) + + Eh[ii+2] = hex8{T}( E[i][3], inv_Et[i_e+2], inv_Ft[i_f ], inv_Et[i_e+1], + inv_Et[i_e+10], inv_Ft[i_f+3], i_vc, inv_Ft[i_f+4] ) + + Eh[ii+3] = hex8{T}( E[i][4], inv_Et[i_e+3], inv_Ft[i_f ], inv_Et[i_e+2], + inv_Et[i_e+11], inv_Ft[i_f+5], i_vc, inv_Ft[i_f+3] ) + + Eh[ii+4] = hex8{T}( E[i][5], inv_Et[i_e+7], inv_Ft[i_f+1], inv_Et[i_e+4], + inv_Et[i_e+8 ], inv_Ft[i_f+5], i_vc, inv_Ft[i_f+2] ) + + Eh[ii+5] = hex8{T}( E[i][6], inv_Et[i_e+4], inv_Ft[i_f+1], inv_Et[i_e+5], + inv_Et[i_e+9 ], inv_Ft[i_f+2], i_vc, inv_Ft[i_f+4] ) + + Eh[ii+6] = hex8{T}( E[i][7], inv_Et[i_e+5], inv_Ft[i_f+1], inv_Et[i_e+6], + inv_Et[i_e+10], inv_Ft[i_f+4], i_vc, inv_Ft[i_f+3] ) + + Eh[ii+7] = hex8{T}( E[i][8], inv_Et[i_e+6], inv_Ft[i_f+1], inv_Et[i_e+7], + inv_Et[i_e+11], inv_Ft[i_f+3], i_vc, inv_Ft[i_f+5] ) + end + elseif isone(direction) # Split in 1st-direction + Vh = [V;Ve] # Append vertices + inv_Et = inv_Et .+ length(V) + Eh = Vector{hex8{T}}(undef,length(E)*2) # Allocate hexahedral element vector, one hex for each of the 4 nodes per element + for i = eachindex(E) + i_e = 1 + (i-1)*4 # index of first edge + i_f = 1 + (i-1)*6 # index of first face + ii = 1 + (i-1)*2 + Eh[ii ] = hex8{T}( Ft[i_f+4][4], Ft[i_f+4][3], Ft[i_f+4][2], Ft[i_f+4][1], + inv_Et[i_e+2], inv_Et[i_e+3], inv_Et[i_e+1], inv_Et[i_e+0] ) + + Eh[ii+1] = hex8{T}( Ft[i_f+5][4], Ft[i_f+5][3], Ft[i_f+5][2], Ft[i_f+5][1], + inv_Et[i_e+0], inv_Et[i_e+1], inv_Et[i_e+3], inv_Et[i_e+2] ) + + end + elseif direction == 2 # Split in 2nd-direction + Vh = [V;Ve] # Append vertices + inv_Et = inv_Et .+ length(V) + Eh = Vector{hex8{T}}(undef,length(E)*2) # Allocate hexahedral element vector, one hex for each of the 4 nodes per element + for i = eachindex(E) + i_e = 1 + (i-1)*4 # index of first edge + i_f = 1 + (i-1)*6 # index of first face + ii = 1 + (i-1)*2 + Eh[ii ] = hex8{T}( Ft[i_f+2][4], Ft[i_f+2][3], Ft[i_f+2][2], Ft[i_f+2][1], + inv_Et[i_e+3], inv_Et[i_e+2], inv_Et[i_e+0], inv_Et[i_e+1] ) + + Eh[ii+1] = hex8{T}( Ft[i_f+3][4], Ft[i_f+3][3], Ft[i_f+3][2], Ft[i_f+3][1], + inv_Et[i_e+2], inv_Et[i_e+3], inv_Et[i_e+1], inv_Et[i_e+0] ) + + end + elseif direction == 3 # Split in 3rd-direction + Vh = [V;Ve] # Append vertices + inv_Et = inv_Et .+ length(V) + Eh = Vector{hex8{T}}(undef,length(E)*2) # Allocate hexahedral element vector, one hex for each of the 4 nodes per element + for i = eachindex(E) + i_e = 1 + (i-1)*4 # index of first edge + i_f = 1 + (i-1)*6 # index of first face + ii = 1 + (i-1)*2 + Eh[ii ] = hex8{T}( Ft[i_f][4], Ft[i_f][3], Ft[i_f][2], Ft[i_f][1], + inv_Et[i_e+0], inv_Et[i_e+1], inv_Et[i_e+2], inv_Et[i_e+3] ) + + Eh[ii+1] = hex8{T}( Ft[i_f+1][4], Ft[i_f+1][3], Ft[i_f+1][2], Ft[i_f+1][1], + inv_Et[i_e+3], inv_Et[i_e+2], inv_Et[i_e+1], inv_Et[i_e+0] ) + + end + end + return Eh,Vh + elseif n>1 + for _ = 1:n + E,V = subhex(E,V,1; direction=direction) + end + return E,V + else + throw(ArgumentError("n should be larger than or equal to 0")) + end +end + +function rhombicdodecahedron(r = 1.0) + F = Vector{QuadFace{Int64}}(undef,12) + F[ 1] = [ 1, 10, 5, 9] + F[ 2] = [ 2, 11, 6, 10] + F[ 3] = [ 3, 12, 7, 11] + F[ 4] = [ 4, 9, 8, 12] + F[ 5] = [ 5, 10, 6, 14] + F[ 6] = [ 6, 11, 7, 14] + F[ 7] = [ 7, 12, 8, 14] + F[ 8] = [ 8, 9, 5, 14] + F[ 9] = [ 1, 13, 2, 10] + F[10] = [ 2, 13, 3, 11] + F[11] = [ 3, 13, 4, 12] + F[12] = [ 4, 13, 1, 9] + + + V = Vector{Point{3,Float64}}(undef,14) + V[ 1] = [-0.5, -0.5, -0.5] + V[ 2] = [ 0.5, -0.5, -0.5] + V[ 3] = [ 0.5, 0.5, -0.5] + V[ 4] = [-0.5, 0.5, -0.5] + V[ 5] = [-0.5, -0.5, 0.5] + V[ 6] = [ 0.5, -0.5, 0.5] + V[ 7] = [ 0.5, 0.5, 0.5] + V[ 8] = [-0.5, 0.5, 0.5] + V[ 9] = [-1.0, -0.0, 0.0] + V[10] = [ 0.0, -1.0, 0.0] + V[11] = [ 1.0, 0.0, 0.0] + V[12] = [ 0.0, 1.0, 0.0] + V[13] = [ 0.0, -0.0, -1.0] + V[14] = [ 0.0, -0.0, 1.0] + + if !isone(r) + V.*=r + end + + return F,V end \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index 0b4c722..c2cc0f6 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1607,8 +1607,8 @@ end E,V,F,Fb,CFb_type = hexbox([1.0,1.0,1.0],[1,1,1]) @test E == [[1, 2, 4, 3, 5, 6, 8, 7]] @test V == Point3{Float64}[[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [1.0, 1.0, 0.0], [0.0, 0.0, 1.0], [1.0, 0.0, 1.0], [0.0, 1.0, 1.0], [1.0, 1.0, 1.0]] - @test F == QuadFace{Int64}[[3, 4, 2, 1], [5, 6, 8, 7], [1, 2, 6, 5], [4, 3, 7, 8], [2, 4, 8, 6], [3, 1, 5, 7]] - @test Fb == QuadFace{Int64}[[3, 4, 2, 1], [5, 6, 8, 7], [1, 2, 6, 5], [4, 3, 7, 8], [2, 4, 8, 6], [3, 1, 5, 7]] + @test F == QuadFace{Int64}[[3, 4, 2, 1], [5, 6, 8, 7], [1, 2, 6, 5], [4, 3, 7, 8], [2, 4, 8, 6], [5, 7, 3, 1]] + @test Fb == QuadFace{Int64}[[3, 4, 2, 1], [5, 6, 8, 7], [1, 2, 6, 5], [4, 3, 7, 8], [2, 4, 8, 6], [5, 7, 3, 1]] @test CFb_type == [1, 2, 3, 4, 5, 6] end @testset "2x2x2 hex box" begin @@ -3043,46 +3043,6 @@ end end -@testset "wrapindex" verbose = true begin - eps_level = 1e-4 - - @testset "single value" begin - n = 5 - m = 2 - @test wrapindex(1,n) == 1 - @test wrapindex(2,n) == 2 - @test wrapindex(n,n) == n - @test wrapindex(n+m,n) == m - end - - @testset "Vector" begin - n = 5 - a = 2 - b = 2*n - m = [1,2,n,n+a,n+b] - @test wrapindex(m,n) == [1,2,n,a,n] - end - - @testset "Unit range" begin - n = 5 - a = 2 - b = 2*n - m = 1:(2*n)+2 - r = wrapindex(m,n) - @test r == [1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2] - end - - @testset "Step range" begin - n = 5 - a = 2 - b = 2*n - m = 1:2:(2*n)+2 - r = wrapindex(m,n) - @test r == [1, 3, 5, 2, 4, 1] - end -end - - @testset "edgeangles" begin eps_level = 1e-4 # Regular cube @@ -4668,4 +4628,123 @@ end [0.8194254478692206, 0.375, 0.36319552381099396]],atol=eps_level) end +end + + + +@testset "tet2hex" verbose=true begin + eps_level = 1e-6 + + E = [tet4{Int64}(1,2,3,4),tet4{Int64}(2,3,4,5),tet4{Int64}(6,7,8,9)] + V = [Point{3,Float64}(-1.0,0.0,0.0), + Point{3,Float64}( 1.0,0.0,0.0), + Point{3,Float64}( 0.0,1.0,0.0), + Point{3,Float64}( 0.0,0.5,1.0), + Point{3,Float64}( 1.0,1.0,1.0), + Point{3,Float64}( 2.0,0.0,0.0), + Point{3,Float64}( 4.0,0.0,0.0), + Point{3,Float64}( 3.0,1.0,0.0), + Point{3,Float64}( 3.0,0.5,1.0), + ] + + Eh,Vh = tet2hex(E,V) + ind = round.(Int64,range(1,length(Vh),10)) + + + @test length(Eh) == length(E)*4 + @test isapprox(Vh[ind],Point{3, Float64}[[-1.0, 0.0, 0.0], [1.0, 1.0, 1.0], [3.0, 0.5, 1.0], + [0.0, 0.3333333333333333, 0.0], [0.6666666666666666, 0.6666666666666666, 0.3333333333333333], + [3.3333333333333335, 0.5, 0.3333333333333333], [-0.5, 0.5, 0.0], [1.0, 0.5, 0.5], + [3.5, 0.5, 0.0], [3.0, 0.75, 0.5]],atol=eps_level) + + @test isa(Eh,Vector{hex8{Int64}}) + @test isa(Vh,Vector{eltype(V)}) + +end + + +@testset "element2faces" verbose=true begin + @testset "hexahedra" begin + boxDim = [1.0,2.0,3.0] # Dimensions for the box in each direction + boxEl = [1,1,1] # Number of elements to use in each direction + E,V,F,Fb,CFb_type = hexbox(boxDim,boxEl) + F = element2faces(E) + @test length(F) == length(E)*6 + @test isa(F,Vector{QuadFace{Int64}}) + @test F[1] == [3,4,2,1] + + boxDim = [1.0,2.0,3.0] # Dimensions for the box in each direction + boxEl = [1,2,3] # Number of elements to use in each direction + E,V,F,Fb,CFb_type = hexbox(boxDim,boxEl) + F = element2faces(E) + @test length(F) == length(E)*6 + @test isa(F,Vector{QuadFace{Int64}}) + @test F[1] == [3,4,2,1] + end + + @testset "tetrahedra" begin + E = [tet4{Int64}(1,2,3,4)] + F = element2faces(E) + @test length(F) == length(E)*4 + @test isa(F,Vector{TriangleFace{Int64}}) + @test F[1] == [3,2,1] + + E = [tet4{Int64}(1,2,3,4),tet4{Int64}(2,3,4,5),tet4{Int64}(6,7,8,9)] + F = element2faces(E) + @test length(F) == length(E)*4 + @test isa(F,Vector{TriangleFace{Int64}}) + @test F[1] == [3,2,1] + end +end + + +@testset "subhex" verbose=true begin + boxDim = [1.0,2.0,3.0] # Dimensions for the box in each direction + boxEl = [1,1,1] # Number of elements to use in each direction + E,V,F,Fb,CFb_type = hexbox(boxDim,boxEl) + + Eh0,Vh0 = subhex(E,V,1;direction=0) + Eh1,Vh1 = subhex(E,V,1;direction=1) + Eh2,Vh2 = subhex(E,V,1;direction=2) + Eh3,Vh3 = subhex(E,V,1;direction=3) + + @test length(Eh0) == length(E)*8 + @test length(Eh1) == length(E)*2 + @test length(Eh2) == length(E)*2 + @test length(Eh3) == length(E)*2 + + @test length(Vh0) == 27 + @test length(Vh1) == 12 + @test length(Vh2) == 12 + @test length(Vh3) == 12 + + @test isa(Vh0,Vector{eltype(V)}) + + Eh0,Vh0 = subhex(E,V,2;direction=0) + @test length(Eh0) == length(E)*8*8 + @test isa(Eh0,Vector{hex8{Int64}}) + @test isa(Vh0,Vector{eltype(V)}) +end + + +@testset "rhombicdodecahedron" verbose=true begin + eps_level = 1e-6 + + F,V = rhombicdodecahedron(1.0) + + Vt = Point{3, Float64}[[-0.5, -0.5, -0.5], [0.5, -0.5, -0.5], [0.5, 0.5, -0.5], + [-0.5, 0.5, -0.5], [-0.5, -0.5, 0.5], [0.5, -0.5, 0.5], [0.5, 0.5, 0.5], [-0.5, 0.5, 0.5], + [-1.0, -0.0, 0.0], [0.0, -1.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, -0.0, -1.0], + [0.0, -0.0, 1.0]] + + @test F == QuadFace{Int64}[[1, 10, 5, 9], [2, 11, 6, 10], [3, 12, 7, 11], [4, 9, 8, 12], + [5, 10, 6, 14], [6, 11, 7, 14], [7, 12, 8, 14], [8, 9, 5, 14], [1, 13, 2, 10], [2, 13, 3, 11], + [3, 13, 4, 12], [4, 13, 1, 9]] + + @test isapprox(V,Vt,atol=eps_level) + + r = π + F,V = rhombicdodecahedron(r) + + @test isapprox(V,Vt.*r,atol=eps_level) end \ No newline at end of file