diff --git a/examples/demo_geosphere.jl b/examples/demo_geosphere.jl index d094af5..276f83a 100644 --- a/examples/demo_geosphere.jl +++ b/examples/demo_geosphere.jl @@ -16,18 +16,29 @@ n2 = 1 # Number of refinement iterations F2,V2 = geosphere(n2,r) n3 = 3 # Number of refinement iterations -F3,V3 = geosphere(n3,r) +method3 = :linear +F3,V3 = geosphere(n3,r; method=method3) + +n4 = 3 # Number of refinement iterations +method4 = :Loop +F4,V4 = geosphere(n4,r; method=method4) + #Visualize mesh +lineWidth = 1 + fig = Figure(size=(1600,800)) -ax1=Axis3(fig[1, 1], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "A geodesic sphere n=0") -hp1=poly!(ax1,V1,F1, strokewidth=2,color=:white, shading = FastShading) +ax1=Axis3(fig[1, 1], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "A geodesic sphere n=$n1") +hp1=poly!(ax1,V1,F1, strokewidth=lineWidth,color=:white, shading = FastShading) + +ax2=Axis3(fig[1, 2], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "A geodesic sphere n=$n2") +hp2=poly!(ax2,V2,F2, strokewidth=lineWidth,color=:white, shading = FastShading) -ax2=Axis3(fig[1, 2], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "A geodesic sphere n=1") -hp2=poly!(ax2,V2,F2, strokewidth=2,color=:white, shading = FastShading) +ax3=Axis3(fig[2, 1], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "A geodesic sphere n=$n3, method=$method3") +hp3=poly!(ax3,V3,F3, strokewidth=lineWidth,color=:white, shading = FastShading) -ax3=Axis3(fig[1, 3], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "A geodesic sphere n=3") -hp3=poly!(ax3,V3,F3, strokewidth=2,color=:white, shading = FastShading) +ax4=Axis3(fig[2, 2], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "A geodesic sphere n=$n4, method=$method4") +hp4=poly!(ax4,V4,F4, strokewidth=lineWidth,color=:white, shading = FastShading) fig \ No newline at end of file diff --git a/examples/demo_hemisphere.jl b/examples/demo_hemisphere.jl index a37098c..cb48ba9 100644 --- a/examples/demo_hemisphere.jl +++ b/examples/demo_hemisphere.jl @@ -6,36 +6,44 @@ using Statistics using LinearAlgebra -# Demo for building a hemissphere from a isosaheddron . +# Demo for building a hemissphere from an isosaheddron. # Example geometry for a sphere that is cut so some edges are boundary edges n = 1# Number of refinement steps of the geodesic sphere -r = 1.0 # Sphere radius +r = 1.0 # hemisphere radius -r = 1.0 -F1,V1 = hemisphere(1,r; face_type=:tri) -F2,V2 = hemisphere(2,r; face_type=:tri) -F3,V3 = hemisphere(3,r; face_type=:tri) +F1,V1 = hemisphere(0,r; face_type=:tri,closed=true) +F2,V2 = hemisphere(1,r; face_type=:tri) +F3,V3 = hemisphere(2,r; face_type=:tri,closed=true) +F4,V4 = hemisphere(3,r; face_type=:tri) -F4,V4 = hemisphere(1,r; face_type=:quad) -F5,V5 = hemisphere(2,r; face_type=:quad) -F6,V6 = hemisphere(3,r; face_type=:quad) +F5,V5 = hemisphere(0,r; face_type=:quad,closed=true) +F6,V6 = hemisphere(1,r; face_type=:quad) +F7,V7 = hemisphere(2,r; face_type=:quad,closed=true) +F8,V8 = hemisphere(2,r; face_type=:quad) # Visualization fig = Figure(size=(1200,800)) ax1 = Axis3(fig[1, 1], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "face_type=:tri, n=0") hp1 = poly!(ax1,GeometryBasics.Mesh(V1,F1), strokewidth=1,color=:white, strokecolor=:black, shading = FastShading, transparency=false) -ax2 = Axis3(fig[1, 2], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "face_type=:tri, n=2") -hp2 = poly!(ax2,GeometryBasics.Mesh(V2,F2), strokewidth=1,color=:white, strokecolor=:black, shading = FastShading, transparency=false) -ax3 = Axis3(fig[1, 3], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "face_type=:tri, n=3") -hp4 = poly!(ax3,GeometryBasics.Mesh(V3,F3), strokewidth=1,color=:white, strokecolor=:black, shading = FastShading, transparency=false) +# normalplot(ax1,F1,V1) -ax4 = Axis3(fig[2, 1], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "face_type=:quad, n=0") +ax2 = Axis3(fig[1, 2], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "face_type=:tri, n=1") +hp2 = poly!(ax2,GeometryBasics.Mesh(V2,F2), strokewidth=1,color=:white, strokecolor=:black, shading = FastShading, transparency=false) +ax3 = Axis3(fig[1, 3], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "face_type=:tri, n=2") +hp3 = poly!(ax3,GeometryBasics.Mesh(V3,F3), strokewidth=1,color=:white, strokecolor=:black, shading = FastShading, transparency=false) +ax4 = Axis3(fig[1, 4], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "face_type=:tri, n=3") hp4 = poly!(ax4,GeometryBasics.Mesh(V4,F4), strokewidth=1,color=:white, strokecolor=:black, shading = FastShading, transparency=false) -ax5 = Axis3(fig[2, 2], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "face_type=:quad, n=2") + +ax5 = Axis3(fig[2, 1], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "face_type=:quad, n=0") hp5 = poly!(ax5,GeometryBasics.Mesh(V5,F5), strokewidth=1,color=:white, strokecolor=:black, shading = FastShading, transparency=false) -ax6 = Axis3(fig[2, 3], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "face_type=:quad, n=3") +# normalplot(ax5,F5,V5) +ax6 = Axis3(fig[2, 2], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "face_type=:quad, n=1") hp6 = poly!(ax6,GeometryBasics.Mesh(V6,F6), strokewidth=1,color=:white, strokecolor=:black, shading = FastShading, transparency=false) +ax7 = Axis3(fig[2, 3], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "face_type=:quad, n=2") +hp7 = poly!(ax7,GeometryBasics.Mesh(V7,F7), strokewidth=1,color=:white, strokecolor=:black, shading = FastShading, transparency=false) +ax8 = Axis3(fig[2, 4], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "face_type=:quad, n=3") +hp8 = poly!(ax8,GeometryBasics.Mesh(V8,F8), strokewidth=1,color=:white, strokecolor=:black, shading = FastShading, transparency=false) fig \ No newline at end of file diff --git a/examples/demo_perlin.jl b/examples/demo_perlin.jl new file mode 100644 index 0000000..2d1061f --- /dev/null +++ b/examples/demo_perlin.jl @@ -0,0 +1,130 @@ +using GLMakie +using Random + +function randangle(siz) + A = Matrix{Float64}(undef,siz) + for i in eachindex(A) + A[i] = rand()*pi*rand((-1,1)) + end + return A +end + +Random.seed!(1) + +# Define grid +size_grid = (25,25) # grid size +sampleFactor = 35 # Pixel sample factor wrt grid +pixelSize = 1/sampleFactor # Pixel size assuming grid has unit steps + +# Create grid vectors +A = randangle(size_grid) # Random angles + +# for w = 0:1:19 +# A .+= (2*pi/20) + +# A = pi/2 .* ones(size_grid) +# for i in eachindex(A) +# A[i] *= rand((-1,1)) +# # println(i) +# end +# A[7]*=-1 + +Ux = cos.(A) # Unit vector x-component +Uy = sin.(A) # Unit vector y-component + +# Initialise image +size_image = (size_grid .- 1) .* sampleFactor # image size +M = Matrix{Float64}(undef,size_image) # Start as undef Float64 + +# Pre-compute grid cell quantities +xy = range(0+pixelSize/2,1-pixelSize/2,sampleFactor) # x or y coordinates within a grid cell + +X = [x for x in 0:sampleFactor:(size_grid[2]-1)*sampleFactor, j in 0:sampleFactor:(size_grid[1]-1)*sampleFactor]' +Y = [y for i in 0:sampleFactor:(size_grid[2]-1)*sampleFactor, y in 0:sampleFactor:(size_grid[1]-1)*sampleFactor]' + +# ig = ceil(Int,i./sampleFactor) # Grid row index + # ip = mod1(i,sampleFactor) # Cell row index + + # jg = ceil(Int,j./sampleFactor) # Grid column index + # jp = mod1(j,sampleFactor) # Cell column index + +# ceil.(collect(1:1:25)./sampleFactor) + +function fade_perlin(a,b,t) + # q = 3.0*t^3 - 2.0*t^3 # Smoothstep + q = 0.5 .- 0.5.*cos.(t*pi) # Cosine + # q = t # Linear + + # Fade function q = 6t⁵-15t⁴+10t³ + # q = t^3 * (t * (6.0 * t - 15.0) + 10.0) # Perlin fade function + + return (1.0-q)*a + q*b +end + +xc = [0,1,1,0] +yc = [0,0,1,1] +for ip in 1:sampleFactor # For each cell row + for jp in 1:sampleFactor # For each cell column + for ig in 1:size_grid[1]-1 # For each grid row + for jg in 1:size_grid[2]-1 # For each grid column + i = (ig-1)*sampleFactor + ip # Pixel row index + j = (jg-1)*sampleFactor + jp # Pixel column index + + # Current pixel cell coordinates + px = xy[jp] + py = xy[ip] + + # Offset vector components + xc1 = px # -xc[1] Offset vector 1 x + xc2 = px-xc[2] # Offset vector 2 x + xc3 = px-xc[3] # Offset vector 3 x + xc4 = px-xc[4] # Offset vector 4 x + + yc1 = py # -yc[2] Offset vector 1 y + yc2 = py-yc[2] # Offset vector 2 y + yc3 = py-yc[3] # Offset vector 3 y + yc4 = py-yc[4] # Offset vector 4 y + + u1x = Ux[ig ,jg] + u2x = Ux[ig ,jg+1] + u3x = Ux[ig+1,jg+1] + u4x = Ux[ig+1,jg] + + u1y = Uy[ig ,jg] + u2y = Uy[ig ,jg+1] + u3y = Uy[ig+1,jg+1] + u4y = Uy[ig+1,jg] + + d1 = xc1.*u1x + yc1.*u1y + d2 = xc2.*u2x + yc2.*u2y + d3 = xc3.*u3x + yc3.*u3y + d4 = xc4.*u4x + yc4.*u4y + + # Interpolation + # Fade function 6t⁵-15t⁴+10t³ + d12 = fade_perlin(d1,d2,px) + d34 = fade_perlin(d4,d3,px) + d = fade_perlin(d12,d34,py) + + M[i,j] = d + + end + end + end +end + +fig = Figure(size=(1500,1500)) +# ax1 = Axis3(fig[1, 1], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "Edge angles") +ax1 = Axis(fig[1, 1], aspect = DataAspect(), title = "Perlin noise",limits=(-sampleFactor,size_image[2]+sampleFactor,-sampleFactor,size_image[1]+sampleFactor) ) +hm = image!(ax1, M',interpolate=false,colormap = Makie.Reverse(:Spectral),colorrange=(-0.5,0.5)) # + +arrows!(ax1,reduce(vcat,X), reduce(vcat,Y), reduce(vcat,Ux), reduce(vcat,Uy), arrowsize = 10, lengthscale = sampleFactor/3, color = :black, linewidth=1) + +Colorbar(fig[1, 2], hm) + +fig + +# save("/home/kevin/DATA/Julia/Comodo.jl/assets/img/perlin_noise_c"*string(w)*".jpg",fig,px_per_unit = 1) +# save("/home/kevin/DATA/Julia/Comodo.jl/assets/img/perlin_noise_large.jpg",fig,px_per_unit = 2) + +# end \ No newline at end of file diff --git a/examples/demo_quaddisc.jl b/examples/demo_quaddisc.jl new file mode 100644 index 0000000..529786c --- /dev/null +++ b/examples/demo_quaddisc.jl @@ -0,0 +1,48 @@ +using Comodo +using GLMakie +using GeometryBasics +using LinearAlgebra + +#= +This is a demonstration of the capabilities of the `quaddisc` function which +generates the faces `F` and vertices `V` for a quadrangulated disc (circle). +=# + +# Define input parameters + + +r = 1.0 # Radius + +n1 = 0 +method1 = :linear +F1,V1 = quaddisc(r,n1; method=method1) + +n2 = 1 +method2 = :Catmull_Clark +F2,V2 = quaddisc(r,n2; method=method2) + +n3 = 2 +method3 = :Catmull_Clark +F3,V3 = quaddisc(r,n3; method=method3) + +n4 = 3 +method4 = :Catmull_Clark +F4,V4 = quaddisc(r,n4; method=method4, orientation=:down) + +# Visualization +fig = Figure(size=(800,800)) + +ax1 = Axis3(fig[1, 1], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "Quandrangulated disc, n=$n1, method=$method1") +hp1 = poly!(ax1,GeometryBasics.Mesh(V1,F1), strokewidth=1,color=:white,shading=FastShading,transparency=false) + +ax2 = Axis3(fig[1, 2], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "Quandrangulated disc, n=$n2, method=$method2") +hp2 = poly!(ax2,GeometryBasics.Mesh(V2,F2), strokewidth=1,color=:white,shading=FastShading,transparency=false) + +ax3 = Axis3(fig[2, 1], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "Quandrangulated disc, n=$n3, method=$method3") +hp3 = poly!(ax3,GeometryBasics.Mesh(V3,F3), strokewidth=1,color=:white,shading=FastShading,transparency=false) + +ax4 = Axis3(fig[2, 2], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "Quandrangulated disc, n=$n4, method=$method4") +hp4 = poly!(ax4,GeometryBasics.Mesh(V4,F4), strokewidth=1,color=:white,shading=FastShading,transparency=false) + + +fig diff --git a/examples/demo_quadplate.jl b/examples/demo_quadplate.jl index fa59d1b..9df4a4a 100644 --- a/examples/demo_quadplate.jl +++ b/examples/demo_quadplate.jl @@ -2,18 +2,15 @@ using Comodo using GLMakie using GeometryBasics -plateDim = [20.0,24.0] -plateElem = [11,16] - -F,V = quadplate(plateDim,plateElem) -N = facenormal(F,V) - +plateDim1 = [20.0,24.0] +plateElem1 = [11,16] +orientation1 = :up +F1,V1 = quadplate(plateDim1,plateElem1; orientation=orientation1) ## Visualize mesh -fig = Figure(size=(1000,1000)) +fig = Figure(size=(1200,800)) + +ax1 = Axis3(fig[1, 1], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "Quadrilateral mesh plate") +hp2 = poly!(ax1,GeometryBasics.Mesh(V1,F1), strokewidth=3,color=:white, shading = FastShading) -ax1=Axis3(fig[1, 1], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "A quadrilateral mesh of a plate",azimuth=-pi/2,elevation=pi/2) -hp2=poly!(ax1,GeometryBasics.Mesh(V,F), strokewidth=3,color=:white, shading = FastShading) -# hp2 = scatter!(ax1, V,markersize=15,color=:orange) -# normalplot(ax1,F,V) fig diff --git a/examples/demo_subquad.jl b/examples/demo_subquad.jl index ac36881..84840f4 100644 --- a/examples/demo_subquad.jl +++ b/examples/demo_subquad.jl @@ -51,11 +51,11 @@ ax1 = Axis3(fig[1, 1], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", wireframe!(ax1,M, linewidth=lineWidth,color=:red, overdraw=false) poly!(ax1,GeometryBasics.Mesh(Vn1,Fn1), strokewidth=strokewidth1,color=:white,shading=FastShading,transparency=false) -ax2 = Axis3(fig[1, 2], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "Linearly, n=2") +ax2 = Axis3(fig[1, 2], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "Linear, n=2") wireframe!(ax2,M, linewidth=lineWidth,color=:red, overdraw=false) poly!(ax2,GeometryBasics.Mesh(Vn2,Fn2), strokewidth=strokewidth1,color=:white,shading=FastShading,transparency=false) -ax3 = Axis3(fig[1, 3], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "Linearly, n=3") +ax3 = Axis3(fig[1, 3], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "Linear, n=3") hp1 = wireframe!(ax3,M, linewidth=lineWidth,color=:red, overdraw=false) hp2 = poly!(ax3,GeometryBasics.Mesh(Vn3,Fn3), strokewidth=strokewidth1,color=:white,shading=FastShading,transparency=false) Legend(fig[1, 4],[hp1,hp2],["Initial","Refined"]) diff --git a/examples/demo_tridisc.jl b/examples/demo_tridisc.jl index 7017618..6480cd6 100644 --- a/examples/demo_tridisc.jl +++ b/examples/demo_tridisc.jl @@ -9,17 +9,40 @@ generates the faces `F` and vertices `V` for a triangulated disc (circle). =# # Define input parameters -r = 2.0 # Radius -n = 3 # Number of refinement iterations +r = 1.0 # Radius # Create disc mesh -F,V = tridisc(r,n) +n1 = 0 +F1,V1 = tridisc(r,n1) -## Visualization +n2 = 2 +F2,V2 = tridisc(r,n2) + +n3 = 3 +ngon3 = 5 +method3 = :linear +F3,V3 = tridisc(r,n3; ngon=ngon3, method=method3) + +n4 = 3 +ngon4 = 6 +method4 = :Loop +F4,V4 = tridisc(r,n4; ngon=ngon4, method=method4) + +# Visualization fig = Figure(size=(800,800)) -ax1 = Axis3(fig[1, 1], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "A triangulated mesh of a circular domain (disc)") -# scatter!(ax1,V,markersize=10,color = :black) -hp2 = poly!(ax1,GeometryBasics.Mesh(V,F), strokewidth=1,color=:white,shading=FastShading,transparency=false) -# normalplot(ax1,F,V) +ax1 = Axis3(fig[1, 1], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "Triangulated disc, n=$n1") +hp1 = poly!(ax1,GeometryBasics.Mesh(V1,F1), strokewidth=1,color=:white,shading=FastShading,transparency=false) +# normalplot(ax1,F1,V1) + +ax2 = Axis3(fig[1, 2], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "Triangulated disc, n=$n2") +hp2 = poly!(ax2,GeometryBasics.Mesh(V2,F2), strokewidth=1,color=:white,shading=FastShading,transparency=false) + +ax3 = Axis3(fig[2, 1], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "Triangulated disc, n=$n3, ngon=$ngon3, method=$method3") +hp3 = poly!(ax3,GeometryBasics.Mesh(V3,F3), strokewidth=1,color=:white,shading=FastShading,transparency=false) + +ax4 = Axis3(fig[2, 2], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "Triangulated disc, n=$n4, ngon=$ngon4, method=$method4") +hp4 = poly!(ax4,GeometryBasics.Mesh(V4,F4), strokewidth=1,color=:white,shading=FastShading,transparency=false) + + fig diff --git a/examples/demo_triplate.jl b/examples/demo_triplate.jl new file mode 100644 index 0000000..4f5d467 --- /dev/null +++ b/examples/demo_triplate.jl @@ -0,0 +1,18 @@ +using Comodo +using GLMakie +using GeometryBasics + +plateDim1 = [20.0,24.0] +pointSpacing1 = 2.0 + +orientation1 = :up +F1,V1 = triplate(plateDim1,pointSpacing1; orientation=orientation1) + +## Visualization +linewidth = 1 +fig = Figure(size=(1200,800)) + +ax1 = Axis3(fig[1, 1], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "Triangulated mesh plate") +hp2 = poly!(ax1,GeometryBasics.Mesh(V1,F1), strokewidth=linewidth,color=:white, shading = FastShading) +# normalplot(ax1,F1,V1) +fig diff --git a/src/Comodo.jl b/src/Comodo.jl index 2c02024..f1c1d72 100644 --- a/src/Comodo.jl +++ b/src/Comodo.jl @@ -34,7 +34,7 @@ export gunique, unique_simplices, unique_dict, occursonce export subtri, subquad, geosphere,hemisphere,quad2tri export icosahedron, tetrahedron, cube, dodecahedron, octahedron, platonicsolid export tofaces, topoints, togeometrybasics_mesh -export quadplate, quadsphere, smoothmesh_laplacian, smoothmesh_hc +export triplate, quadplate, quadsphere, smoothmesh_laplacian, smoothmesh_hc export mergevertices, separate_vertices, remove_unused_vertices export edgecrossproduct, facearea, facenormal, vertexnormal, normalizevector export con_edge_face, con_edge_edge, con_face_edge, con_face_face, con_face_face_v, con_vertex_edge, con_vertex_vertex_f @@ -45,7 +45,7 @@ export edgeangles, count_edge_face, boundaryedges, boundaryfaces, boundaryfacein 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 batman, tridisc, quaddisc, regiontrimesh, scalesimplex, subcurve, dualclad export tet2hex, element2faces, subhex, rhombicdodecahedron, tri2quad export tetgenmesh, surfacevolume, tetvolume, extrudefaces, filletcurve export squircle, circlerange, edgefaceangles,faceanglesegment diff --git a/src/functions.jl b/src/functions.jl index bc67adf..9e0e256 100644 --- a/src/functions.jl +++ b/src/functions.jl @@ -1829,12 +1829,12 @@ function subquad(F::Vector{NgonFace{4,TF}},V::Vector{Point{ND,TV}},n::Int; metho Fn = Vector{QuadFace{TF}}(undef,length(F)*4) nv = length(V) ne = length(Eu) - for q in eachindex(F) - i = 1 + (q-1)*4 + for q in eachindex(F) for ii = 0:3 - Fn[i+ii] = QuadFace{TF}([F[q][ii+1], con_F2E[q][ii+1]+nv, q+nv+ne, con_F2E[q][1+mod(3+ii,4)]+nv]) + Fn[q+ii*length(F)] = QuadFace{TF}([F[q][ii+1], con_F2E[q][ii+1]+nv, q+nv+ne, con_F2E[q][1+mod(3+ii,4)]+nv]) end end + return Fn,Vn elseif n>1 for _ =1:n @@ -1853,7 +1853,7 @@ end """ - geosphere(n::Int,r::T) where T <: Real + geosphere(n::Int,r::T; method=:linear) where T <: Real Returns a geodesic sphere triangulation @@ -1864,17 +1864,28 @@ refinement iterations `n` and the radius `r`. Geodesic spheres (aka Buckminster- spheres) are triangulations of a sphere that have near uniform edge lenghts. The algorithm starts with a regular icosahedron. Next this icosahedron is refined `n` times, while nodes are pushed to a sphere surface with radius `r` at each -iteration. +iteration. Two methods are available, i.e. `:linear` (default) and `:Loop` +(see also `subtri`). The former features simply linear splitting while the latter +features the Loop method which may produce a smoother result. """ -function geosphere(n::Int,r::T) where T <: Real +function geosphere(n::Int,r::T; method=:linear) where T <: Real M = platonicsolid(4,r) V = coordinates(M) - F = faces(M) - for _ = 1:n - numPointsNow = length(V) - F,V = subtri(F,V,1) - @inbounds for i in numPointsNow+1:length(V) - V[i] = pushtoradius_(V[i],r) # Overwrite points + F = faces(M) + for _ = 1:n + # Set push start + if method == :linear + s = length(V) # push from here as linear leaves original in place + elseif method == :Loop + s = 1 # push all as Loop alters coordinates + end + + # Sub-divide sphere + F,V = subtri(F,V,1; method=method) + + # Push altered/new points to sphere + @inbounds for i in s:length(V) + V[i] = pushtoradius_(V[i],r) end end return F,V @@ -1894,7 +1905,7 @@ The algorithm starts with a regular icosahedron that is adjusted to generate a h Next this icosahedron is refined `n` times, while nodes are pushed to a sphere surface with radius `r` at each iteration. """ -function hemisphere(n::Int,r::T; face_type=:tri) where T <: Real +function hemisphere(n::Int,r::T; face_type=:tri, closed=false) where T <: Real if n<0 throw(ArgumentError("n should be >= 0")) end @@ -1904,7 +1915,7 @@ function hemisphere(n::Int,r::T; face_type=:tri) where T <: Real F,V = geosphere(1,r) # A once refined icosahedron has an "equator" line # Rotating sphere so a vertex "points north" and equator is in xy-plane - Q = RotXYZ(0.0,atan(Base.MathConstants.golden),0.0) # Rotation matrix + Q = RotXYZ(0.0,atan(Base.MathConstants.golden),0.0) # Rotation matrix V = [Point{3, Float64}(Q*v) for v in V] # Rotate vertices elseif face_type == :quad F,V = quadsphere(1,r) # A once refined cube has an "equator" line @@ -1917,18 +1928,51 @@ function hemisphere(n::Int,r::T; face_type=:tri) where T <: Real F = [F[i] for i in findall(map(f-> mean([V[j][3] for j in f])>=-searchTol,F))] # Remove faces below equator F,V = remove_unused_vertices(F,V) # Cleanup/remove unused vertices after faces were removed + if closed==true + if face_type == :tri + Fb,Vb = tridisc(r,1; ngon=5, method=:linear, orientation=:down) + Q = RotXYZ(0.0,0.0,pi/2.0) # Rotation matrix + Vb = [Point{3, Float64}(Q*v) for v in Vb] # Rotate vertices + elseif face_type == :quad + Fb,Vb = quaddisc(r,0; method=:linear, orientation=:down) + end + + # Add base and merge nodes + indTopFaces = 1:length(F) # Indices of top + append!(F, [f .+ length(V) for f in Fb]) + append!(V,Vb) + F,V,_,_ = mergevertices(F,V) + end + # Refine sphere if needed if n>0 # If larger then refining is needed - for _ = 1:n # Refine once n times - numPointsNow = length(V) # Number of points prior to refining - # Change refine behaviour based on face_type - if face_type == :tri - F,V = subtri(F,V,1) + @inbounds for _ = 1:n # Refine once n times + nv = length(V) # Number of points prior to refining + nf = length(F) + + # Change refine behaviour based on face_type + if face_type == :tri + F,V = subtri(F,V,1) elseif face_type == :quad - F,V = subquad(F,V,1) + F,V = subquad(F,V,1) end + + if closed == true + indTopFaces = [j .+ i*nf for i = 0:3 for j in indTopFaces] + indPush = Vector{Int64}() + for f in F[indTopFaces] + for i in f + if i>nv + push!(indPush,i) + end + end + end + else + indPush = nv+1:length(V) + end + # Now push newly introduced points to the sphere - @inbounds for i in numPointsNow+1:length(V) + @inbounds for i in indPush V[i] = pushtoradius_(V[i],r) # Overwrite points end end @@ -2372,16 +2416,17 @@ point set. function mergevertices(F::Vector{NgonFace{N,TF}},V::Vector{Point{ND,TV}}; roundVertices = true, numDigitsMerge=nothing) where N where TF<:Integer where ND where TV<:Real m = length(V) - if roundVertices + if roundVertices == true if isnothing(numDigitsMerge) E = meshedges(F) d = [sqrt( sum((V[e[1]] .- V[e[2]]).^2) ) for e in E] pointSpacing = mean(d) numDigitsMerge = 6-round(Int,log10(pointSpacing)) end - + # Create rounded coordinates to help obtain unique set - VR = [round.(v,digits = numDigitsMerge) for v in V] + # Note -0.0+0.0 = 0.0 so addition of zero points helps avoid 0.0 and -0.0 being seen as unique + VR = [round.(v,digits = numDigitsMerge)+Point{ND,TV}(0.0,0.0,0.0) for v in V] # Get unique indices and reverse for rounded vertices _,ind1,ind2 = gunique(VR; return_index=true, return_inverse=true,sort_entries=false) @@ -2392,8 +2437,8 @@ function mergevertices(F::Vector{NgonFace{N,TF}},V::Vector{Point{ND,TV}}; roundV if length(V) != m # If the length has changed # Correct indices for faces - for q in eachindex(F) - F[q] = ind2[F[q]] + for i in eachindex(F) + F[i] = ind2[F[i]] end end @@ -2416,7 +2461,7 @@ in the range (0,1). If `λ=0` then no smoothing occurs. If `λ=1` then pure Laplacian mean based smoothing occurs. For intermediate values a linear blending between the two occurs. """ -function smoothmesh_laplacian(F::Vector{NgonFace{N,TF}},V::Vector{Point{ND,TV}}, n=1, λ=0.5; con_V2V=nothing, constrained_points=nothing) where N where TF<:Integer where ND where TV<:Real +function smoothmesh_laplacian(F::Vector{NgonFace{N,TF}},V::Vector{Point{ND,TV}}, n=1, λ=0.5; con_V2V=nothing, tolDist=nothing, constrained_points=nothing) where N where TF<:Integer where ND where TV<:Real if λ>1.0 || λ<0.0 throw(ArgumentError("λ should be in the range 0-1")) @@ -2431,7 +2476,9 @@ function smoothmesh_laplacian(F::Vector{NgonFace{N,TF}},V::Vector{Point{ND,TV}}, E_uni = meshedges(F;unique_only=true) con_V2V = con_vertex_vertex(E_uni) end - for _ = 1:n + c = 0 + while c1 for _ in 1:(growsteps-1) ind = unique(reduce(vcat,con_V2V[ind])) @@ -3487,14 +3549,7 @@ function mesh_curvature_polynomial(F::Vector{NgonFace{N,TF}},V::Vector{Point{ND, end vr = [Q*(v-V[q]) for v in V[ind]] # Rotate point set to a 2D problem - - # if !isnothing(indBoundary) - # if in(q,indBoundary) - # # println("wow") - # push!(vr,normalizevector(mean(-vr))) - # end - # end - + # Set up polynomial fit T = Matrix{Float64}(undef,(length(ind),5)) w = Vector{Float64}(undef,length(ind)) @@ -3529,9 +3584,6 @@ function mesh_curvature_polynomial(F::Vector{NgonFace{N,TF}},V::Vector{Point{ND, U1[q] = Vec3{Float64}(NaN,NaN,NaN) U2[q] = Vec3{Float64}(NaN,NaN,NaN) end - - - end H = 0.5 * (K1.+K2) # Mean curvature @@ -3979,42 +4031,123 @@ function batman(n::Int; symmetric = false, dir=:acw) end """ - tridisc(r=1.0,n=0) + tridisc(r=1.0,n=0; ngon=6, method = :linear, orientation=:up) # Description Generates the faces `F` and vertices `V` for a triangulated disc (circle). The algorithm starts with a triangulated hexagon (obtained if `n=0`) and uses -iterative subtriangulation circular coordinate correction to obtain the final -mesh. -""" -function tridisc(r=1.0,n=0) +iterative subtriangulation, and uses iterative subdivision (and pushing of +boundary points to circular boundary) to obtain the final mesh. The subdivision +`method` is an optional input, and is either `:Loop` (default) or `:linear`. +Lastly the optional input `orientation`, which can be `:up` or `:down` sets the +face normal direction. +""" +function tridisc(r=1.0,n=0; ngon=6, method=:Loop, orientation=:up) + if !in(orientation,(:up,:down)) + throw(ArgumentError("Orientation not supported. Use :up or :down")) + end + # Create a triangulated hexagon - V = circlepoints(r,6) + V = circlepoints(r,ngon) push!(V,Point{3,Float64}(0.0,0.0,0.0)) - F = [TriangleFace{Int}(1,2,7), TriangleFace{Int}(2,3,7), TriangleFace{Int}(3,4,7), TriangleFace{Int}(4,5,7), TriangleFace{Int}(5,6,7), TriangleFace{Int}(6,1,7)] + if orientation==:up + F = [TriangleFace{Int64}(i,mod1(i+1,ngon),ngon+1) for i in 1:ngon] + else#if orientation==:down + F = [TriangleFace{Int64}(i,ngon+1,mod1(i+1,ngon)) for i in 1:ngon] + end # Refine mesh n times if n>0 @inbounds for _ = 1:n - F,V = subtri(F,V,1; method = :linear) - indB = elements2indices(boundaryedges(F)) - V[indB] = r.*V[indB]./norm.(V[indB]) # Push to conform to circle - end - indB = elements2indices(boundaryedges(F)) + nv = length(V) + F,V = subtri(F,V,1; method = method) + indBoundary = unique(reduce(vcat,boundaryedges(F))) + @inbounds for i in indBoundary + if i>nv || method == :Loop + V[i] *= r/norm(V[i]) + end + end + end end return F,V end -# """ -# triplate(xSpan,ySpan,pointSpacing::T) where T <: Real +""" + triplate(xSpan,ySpan,pointSpacing::T) where T <: Real -# # Description -# Generates a triangulated mesh for a plate. -# """ -# function triplate(xSpan,ySpan,pointSpacing::T) where T <: Real -# return gridpoints_equilateral(xSpan,ySpan,pointSpacing; return_faces = true, rectangular=true) -# end +# Description +Generates a triangulated mesh for a plate. +""" +function triplate(plateDim,pointSpacing::T; orientation=:up) where T <: Real + if !in(orientation,(:up,:down)) + throw(ArgumentError("Orientation not supported. Use :up or :down")) + end + + F, V = gridpoints_equilateral((-plateDim[1]/2,plateDim[1]/2),(-plateDim[2]/2,plateDim[2]/2),pointSpacing; return_faces = true, rectangular=true) + if orientation == :down + return invert_faces(F), V + else + return F,V + end +end + +""" + quaddisc(r,n; method = :Catmull_Clark, orientation=:up) + +# Description + +Generates the faces `F` and vertices `V` for a quadrangulated (circle). The +algorithm starts with 12 quadrilateral faces and an octagon boundary +(obtained if `n=0`), and uses iterative subdivision (and pushing of boundary +points to circular boundary) to obtain the final mesh. The subdivision +`method` is an optional input, and is either `:Catmull_Clark` (default) or +`:linear`. Lastly the optional input `orientation`, which can be `:up` or +`:down` sets the face normal direction. +""" +function quaddisc(r,n; method = :Catmull_Clark, orientation=:up) + + if !in(orientation,(:up,:down)) + throw(ArgumentError("Orientation not supported. Use :up or :down")) + end + + # Create disc mesh + V = circlepoints(r,8) + V = append!(V, circlepoints(r/2,8)) + V = push!(V, Point{3,Float64}(0.0,0.0,0.0)) + F = Vector{QuadFace{Int64}}(undef,12) + + # Add outer ring + for i = 1:8 + F[i] = QuadFace{Int64}(i, mod1(i+1,8), mod1(i+1+8,8)+8, mod1(i+8,16)) + end + + # Add core + for i = 1:4 + j = 9 + (i-1)*2 + F[i+8] = QuadFace{Int64}(j, mod1(j+1+8,8)+8, mod1(j+2+8,8)+8, 17) + end + + # Invert if needed + if orientation==:down + F = invert_faces(F) + end + + # Refine + for _ = 1:1:n + nv = length(V) + F,V = subquad(F,V,1; method=method) + indBoundary = unique(reduce(vcat,boundaryedges(F))) + @inbounds for i in indBoundary + if i>nv || method == :Catmull_Clark + V[i] *= r/norm(V[i]) + end + end + # e = (2.0 .* pi .* r) / length(indBoundary)#(8.0 .* 2.0 .^(n.-1)) + # V = smoothmesh_laplacian(F,V, 1e3; constrained_points=indBoundary, tolDist=e/100) + end + return F,V +end """ @@ -5171,8 +5304,7 @@ function faceanglesegment(F::Vector{NgonFace{NF,TF}},V::Vector{Point{ND,TV}}; de if !isnothing(indEdgeFound) # If we found an edge if abs(A[indEdgeFound])<=angleThreshold # Is angle at this edge low enough for indFace in con_E2F[indEdgeFound] # then add the other face - if indFace!=indFaceNow # If other face - # println(indFace) + if indFace!=indFaceNow # If other face push!(indToCheckNew,indFace) # Add to check for next round G[indFace] = g # Assign group label end diff --git a/test/runtests.jl b/test/runtests.jl index bd6dbc4..ed2ca86 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1599,10 +1599,7 @@ end @test eltype(F) == eltype(Fn) @test length(Fn) == 4^n*length(F) @test eltype(V) == eltype(Vn) - @test isapprox(Vn[ind], Point{3, Float64}[[-0.5773502691896258, -0.5773502691896258, -0.5773502691896258], - [0.2886751345948129, -0.2886751345948129, -0.5773502691896258], [0.4330127018922194, 0.5773502691896258, -0.5773502691896258], - [-0.5773502691896258, 0.2886751345948129, 0.14433756729740646], [0.14433756729740646, -0.14433756729740646, 0.5773502691896258], - [0.14433756729740646, -0.5773502691896258, -0.43301270189221935]], atol=eps_level) + @test isapprox(Vn[ind], Point{3, Float64}[[-0.5773502691896258, -0.5773502691896258, -0.5773502691896258], [-0.2886751345948129, 0.5773502691896258, 0.2886751345948129], [-0.5773502691896258, 0.0, 0.14433756729740646], [-0.2886751345948129, 0.14433756729740646, 0.5773502691896258], [0.4330127018922194, -0.43301270189221935, -0.5773502691896258], [0.14433756729740646, -0.5773502691896258, -0.43301270189221935]], atol=eps_level) end @testset "Catmull_Clark" begin @@ -1615,10 +1612,7 @@ end @test eltype(F) == eltype(Fn) @test length(Fn) == 4^n*length(F) @test eltype(V) == eltype(Vn) - @test isapprox(Vn[ind], Point{3, Float64}[[-0.2895661072324513, -0.2895661072324513, -0.2895661072324513], - [0.18414681640860342, -0.18414681640860342, -0.4257509268592806], [0.24244200756986245, 0.31179169810728824, -0.31179169810728824], - [-0.44889359308574917, 0.19141642225574024, 0.09422035643025145], [0.0977285611909523, -0.0977285611909523, 0.47360764269461497], - [0.0907121516695506, -0.407828803431474, -0.2770228830681994]], atol=eps_level) + @test isapprox(Vn[ind], Point{3, Float64}[[-0.2895661072324513, -0.2895661072324513, -0.2895661072324513], [-0.18414681640860342, 0.4257509268592806, 0.18414681640860342], [-0.48187698248769556, 0.0, 0.09948266357130273], [-0.19141642225574024, 0.09422035643025145, 0.44889359308574917], [0.25066958302055375, -0.25066958302055375, -0.3566674840045866], [0.0907121516695506, -0.407828803431474, -0.2770228830681994]], atol=eps_level) end @testset "Catmull_Clark, unconstrained boundary" begin @@ -1638,9 +1632,7 @@ end @test eltype(F) == eltype(Fn) @test length(Fn) == 4^n*length(F) @test eltype(V) == eltype(Vn) - @test isapprox(Vn[ind], Point{3, Float64}[[0.5078125, 3.2959746043559335e-17, 0.0], [0.4228515625000001, 0.21143198334581026, 1.7320508075688776], - [-0.4296875, -0.13531646934131847, 0.8660254037844388], [-0.39453125, -0.26048420348203816, 0.6495190528383292], - [-0.14843749999999978, 0.47360764269461497, 0.2165063509461097], [0.09765624999999989, -0.4397785253592852, 1.5155444566227678]], atol=eps_level) + @test isapprox(Vn[ind], Point{3, Float64}[[0.5078125, 3.2959746043559335e-17, 0.0], [-0.39453125, -0.2604842034820381, 0.0], [0.09765624999999989, -0.4397785253592852, 0.8660254037844388], [0.4228515625000001, 0.21143198334581026, 1.0825317547305486], [-0.33593750000000006, -0.36535446722155995, 0.2165063509461097], [0.09765624999999989, -0.4397785253592852, 1.5155444566227678]], atol=eps_level) end @testset "Catmull_Clark, constrained boundary" begin @@ -1657,12 +1649,11 @@ end Fn, Vn = subquad(F, V, n; method=:Catmull_Clark, constrain_boundary=true) ind = round.(Int,range(1,length(Vn),6)) # Sample indices + V_true = Point{3, Float64}[[1.0, 0.0, 0.0], [-0.5, -0.43301270189221924, 0.0], [0.08666992187499989, -0.47149332286115675, 0.8660254037844388], [0.4898681640625001, 0.21333487119592254, 1.0825317547305486], [-0.4472656250000001, -0.5581804360329389, 0.2165063509461097], [0.07128906249999988, -0.5158940393637769, 1.5155444566227678]] @test eltype(F) == eltype(Fn) @test length(Fn) == 4^n*length(F) @test eltype(V) == eltype(Vn) - @test isapprox(Vn[ind], Point{3, Float64}[[1.0, 0.0, 0.0], [0.6250000000000001, 0.2165063509461097, 1.7320508075688776], - [-0.45166015625, -0.1606883073428157, 0.8660254037844388], [-0.4296875, -0.3175708389854069, 0.6495190528383292], - [-0.2597656249999997, 0.6664336115059939, 0.2165063509461097], [0.07128906249999988, -0.5158940393637769, 1.5155444566227678]], atol=eps_level) + @test isapprox(Vn[ind], V_true, atol=eps_level) end @testset "errors" begin @@ -1691,6 +1682,16 @@ end [0.21302286564912973, -0.5712516591357086, -0.7926492292592814], [0.840177885327139, -0.5192584897281833, 0.15643446504023087], [-0.7838430424199713, 0.08108629344330351, -0.6156420208736807]], atol=eps_level) @test isapprox(norm.(V),fill(r,length(V))) # Test radii + + n = 4 + F, V = geosphere(n, r; method=:Loop) + ind = round.(Int,range(1,length(V),6)) # Sample indices + @test F isa Vector{TriangleFace{Int}} + @test length(F) == 4^n*20 # 20 since the icosahedron is the start geometry + @test V isa Vector{Point{3,Float64}} + V_true = Point{3, Float64}[[2.545820645877024e-18, -0.5257311121191336, -0.8506508083520399], [0.9876761238145257, 0.13653930829468125, -0.07650419424530518], [-0.2650128038205533, -0.9439540325795487, -0.196771436412852], [0.19531467334213626, -0.8698675244963651, 0.45297093527490284], [-0.22030556728004222, 0.665396606253752, 0.7132410626228752], [-0.8270273152493076, 0.0310097456068568, -0.561305812823028]] + @test isapprox(V[ind],V_true, atol=eps_level) + @test isapprox(norm.(V),fill(r,length(V))) # Test radii end @@ -2508,16 +2509,16 @@ end @testset "quadsphere" begin eps_level = 1e-4 - r = 1.0 F, V = quadsphere(3, r) - + ind = round.(Int,range(1,length(V),6)) + V_true = Point{3, Float64}[[-0.5773502691896258, -0.5773502691896258, -0.5773502691896258], [-0.36700100107065714, 0.8547634353587377, 0.36700100107065714], [-0.9807852804032304, 0.0, 0.1950903220161283], [-0.3815651304866988, 0.18679164007781993, 0.9052717461589679], [0.4942387316114103, -0.4942387316114103, -0.7151616267322294], [0.1731315037005652, -0.8165807660291813, -0.550655368608694]] @test V isa Vector{Point3{Float64}} @test length(V) == 386 - @test isapprox(V[1], [-0.5773502691896258, -0.5773502691896258, -0.5773502691896258], atol=eps_level) + @test isapprox(V[ind], V_true, atol=eps_level) @test F isa Vector{QuadFace{Int}} @test length(F) == 384 - @test F[1] == [1, 99, 291, 187] + @test F[1] == [1, 99, 291, 116] end @testset "simplex2vertexdata" verbose = true begin @@ -2556,34 +2557,41 @@ end # Quads DF = collect(Float64,1:length(Fq)) # Face data (here face numbers) DV = simplex2vertexdata(Fq,DF,Vq; weighting=:none) - @test isapprox(DV,[12.0, 9.666666666666666, 12.666666666666666, 15.666666666666666, 13.0, 10.0, 12.333333333333334, 14.666666666666666, 6.5, 11.5, 8.5, 10.0, 14.0, 9.0, 12.5, 16.5, 20.5, 16.5, 11.5, 13.0, 2.5, 6.5, 10.5, 14.5, 18.5, 22.5], atol=eps_level) + D_true = [13.333333333333334, 14.666666666666666, 17.333333333333332, 20.0, 11.666666666666666, 9.0, 7.666666666666667, 6.333333333333333, 11.0, 6.5, 11.5, 9.0, 10.0, 14.5, 12.5, 13.5, 14.5, 13.5, 18.0, 15.5, 10.0, 11.0, 12.0, 13.0, 14.0, 15.0] + @test isapprox(DV,D_true, atol=eps_level) # Triangles DF = collect(Float64,1:length(Ft)) # Face data (here face numbers) DV = simplex2vertexdata(Ft,DF,Vt; weighting=:none) - @test isapprox(DV,[10.333333333333334, 11.333333333333334, 13.333333333333334, 7.0, 6.833333333333333, 7.166666666666667, 8.166666666666666, 7.666666666666667, 8.666666666666666, 8.5], atol=eps_level) + D_true = [10.333333333333334, 11.333333333333334, 13.333333333333334, 7.0, 6.833333333333333, 7.166666666666667, 8.166666666666666, 7.666666666666667, 8.666666666666666, 8.5] + @test isapprox(DV,D_true, atol=eps_level) # Hexahedral elements DF = collect(Float64,1:length(Eh)) # Face data (here face numbers) DV = simplex2vertexdata(Eh,DF,Vh; weighting=:none) - @test isapprox(DV,[1.0, 1.5, 2.0, 2.0, 2.5, 3.0, 3.0, 3.5, 4.0, 3.0, 3.5, 4.0, 4.0, 4.5, 5.0, 5.0, 5.5, 6.0, 5.0, 5.5, 6.0, 6.0, 6.5, 7.0, 7.0, 7.5, 8.0], atol=eps_level) + D_true = [1.0, 1.5, 2.0, 2.0, 2.5, 3.0, 3.0, 3.5, 4.0, 3.0, 3.5, 4.0, 4.0, 4.5, 5.0, 5.0, 5.5, 6.0, 5.0, 5.5, 6.0, 6.0, 6.5, 7.0, 7.0, 7.5, 8.0] + @test isapprox(DV,D_true, atol=eps_level) end @testset "Vector Float64 data, weighting=:area" begin # Single element DF = [1.0] # Face data DV = simplex2vertexdata(F1,DF,V1; weighting=:area) - @test isapprox(DV,ones(Float64,length(F1[1])), atol=eps_level) + D_true = [1.0, 1.0, 1.0, 1.0] + @test isapprox(DV,D_true, atol=eps_level) + # Quads DF = collect(Float64,1:length(Fq)) # Face data (here face numbers) DV = simplex2vertexdata(Fq,DF,Vq; weighting=:area) - @test isapprox(DV,[12.0, 9.666666666666668, 12.666666666666666, 15.666666666666664, 13.0, 10.0, 12.333333333333332, 14.666666666666666, 6.500000000000001, 11.5, 8.5, 10.0, 14.0, 9.0, 12.5, 16.5, 20.5, 16.5, 11.5, 13.000000000000002, 2.5, 6.500000000000001, 10.5, 14.5, 18.5, 22.5], atol=eps_level) + D_true = [13.333333333333332, 14.666666666666668, 17.333333333333336, 20.0, 11.666666666666666, 9.0, 7.666666666666667, 6.333333333333333, 11.0, 6.500000000000001, 11.5, 9.0, 10.0, 14.5, 12.500000000000002, 13.5, 14.5, 13.5, 18.0, 15.5, 10.0, 10.999999999999998, 12.000000000000002, 13.000000000000002, 14.0, 15.0] + @test isapprox(DV,D_true, atol=eps_level) # Triangles DF = collect(Float64,1:length(Ft)) # Face data (here face numbers) DV = simplex2vertexdata(Ft,DF,Vt; weighting=:area) - @test isapprox(DV,[10.333333333333334, 11.333333333333334, 13.333333333333336, 6.999999999999999, 5.095917942265425, 5.646428199482246, 6.646428199482247, 6.146428199482247, 6.49489742783178, 6.545407685048602], atol=eps_level) + D_true = [10.333333333333334, 11.333333333333334, 13.333333333333336, 6.999999999999999, 5.095917942265425, 5.646428199482246, 6.646428199482247, 6.146428199482247, 6.49489742783178, 6.545407685048602] + @test isapprox(DV,D_true, atol=eps_level) end @testset "Vector of Vector Float64 data, weighting=:none" begin @@ -2595,17 +2603,20 @@ end # Quads DF = [i.*[1.0,2.0,3.0] for i in eachindex(Fq)] # Vector data for each face DV = simplex2vertexdata(Fq,DF,Vq; weighting=:none) - @test isapprox(DV,[[12.0, 24.0, 36.0], [9.666666666666666, 19.333333333333332, 29.0], [12.666666666666666, 25.333333333333332, 38.0], [15.666666666666666, 31.333333333333332, 47.0], [13.0, 26.0, 39.0], [10.0, 20.0, 30.0], [12.333333333333334, 24.666666666666668, 37.0], [14.666666666666666, 29.333333333333332, 44.0], [6.5, 13.0, 19.5], [11.5, 23.0, 34.5], [8.5, 17.0, 25.5], [10.0, 20.0, 30.0], [14.0, 28.0, 42.0], [9.0, 18.0, 27.0], [12.5, 25.0, 37.5], [16.5, 33.0, 49.5], [20.5, 41.0, 61.5], [16.5, 33.0, 49.5], [11.5, 23.0, 34.5], [13.0, 26.0, 39.0], [2.5, 5.0, 7.5], [6.5, 13.0, 19.5], [10.5, 21.0, 31.5], [14.5, 29.0, 43.5], [18.5, 37.0, 55.5], [22.5, 45.0, 67.5]], atol=eps_level) + D_true = [[13.333333333333334, 26.666666666666668, 40.0], [14.666666666666666, 29.333333333333332, 44.0], [17.333333333333332, 34.666666666666664, 52.0], [20.0, 40.0, 60.0], [11.666666666666666, 23.333333333333332, 35.0], [9.0, 18.0, 27.0], [7.666666666666667, 15.333333333333334, 23.0], [6.333333333333333, 12.666666666666666, 19.0], [11.0, 22.0, 33.0], [6.5, 13.0, 19.5], [11.5, 23.0, 34.5], [9.0, 18.0, 27.0], [10.0, 20.0, 30.0], [14.5, 29.0, 43.5], [12.5, 25.0, 37.5], [13.5, 27.0, 40.5], [14.5, 29.0, 43.5], [13.5, 27.0, 40.5], [18.0, 36.0, 54.0], [15.5, 31.0, 46.5], [10.0, 20.0, 30.0], [11.0, 22.0, 33.0], [12.0, 24.0, 36.0], [13.0, 26.0, 39.0], [14.0, 28.0, 42.0], [15.0, 30.0, 45.0]] + @test isapprox(DV,D_true, atol=eps_level) # Triangles DF = [i.*[1.0,2.0,3.0] for i in eachindex(Ft)] # Vector data for each face DV = simplex2vertexdata(Ft,DF,Vt; weighting=:none) - @test isapprox(DV,[[10.333333333333334, 20.666666666666668, 31.0], [11.333333333333334, 22.666666666666668, 34.0], [13.333333333333334, 26.666666666666668, 40.0], [7.0, 14.0, 21.0], [6.833333333333333, 13.666666666666666, 20.5], [7.166666666666667, 14.333333333333334, 21.5], [8.166666666666666, 16.333333333333332, 24.5], [7.666666666666667, 15.333333333333334, 23.0], [8.666666666666666, 17.333333333333332, 26.0], [8.5, 17.0, 25.5]], atol=eps_level) - + D_true = [[10.333333333333334, 20.666666666666668, 31.0], [11.333333333333334, 22.666666666666668, 34.0], [13.333333333333334, 26.666666666666668, 40.0], [7.0, 14.0, 21.0], [6.833333333333333, 13.666666666666666, 20.5], [7.166666666666667, 14.333333333333334, 21.5], [8.166666666666666, 16.333333333333332, 24.5], [7.666666666666667, 15.333333333333334, 23.0], [8.666666666666666, 17.333333333333332, 26.0], [8.5, 17.0, 25.5]] + @test isapprox(DV,D_true, atol=eps_level) + # Hexahedral elements DF = [i.*[1.0,2.0,3.0] for i in eachindex(Eh)] # Vector data for each face DV = simplex2vertexdata(Eh,DF,Vh; weighting=:none) - @test isapprox(DV,[[1.0, 2.0, 3.0], [1.5, 3.0, 4.5], [2.0, 4.0, 6.0], [2.0, 4.0, 6.0], [2.5, 5.0, 7.5], [3.0, 6.0, 9.0], [3.0, 6.0, 9.0], [3.5, 7.0, 10.5], [4.0, 8.0, 12.0], [3.0, 6.0, 9.0], [3.5, 7.0, 10.5], [4.0, 8.0, 12.0], [4.0, 8.0, 12.0], [4.5, 9.0, 13.5], [5.0, 10.0, 15.0], [5.0, 10.0, 15.0], [5.5, 11.0, 16.5], [6.0, 12.0, 18.0], [5.0, 10.0, 15.0], [5.5, 11.0, 16.5], [6.0, 12.0, 18.0], [6.0, 12.0, 18.0], [6.5, 13.0, 19.5], [7.0, 14.0, 21.0], [7.0, 14.0, 21.0], [7.5, 15.0, 22.5], [8.0, 16.0, 24.0]], atol=eps_level) + D_true = [[1.0, 2.0, 3.0], [1.5, 3.0, 4.5], [2.0, 4.0, 6.0], [2.0, 4.0, 6.0], [2.5, 5.0, 7.5], [3.0, 6.0, 9.0], [3.0, 6.0, 9.0], [3.5, 7.0, 10.5], [4.0, 8.0, 12.0], [3.0, 6.0, 9.0], [3.5, 7.0, 10.5], [4.0, 8.0, 12.0], [4.0, 8.0, 12.0], [4.5, 9.0, 13.5], [5.0, 10.0, 15.0], [5.0, 10.0, 15.0], [5.5, 11.0, 16.5], [6.0, 12.0, 18.0], [5.0, 10.0, 15.0], [5.5, 11.0, 16.5], [6.0, 12.0, 18.0], [6.0, 12.0, 18.0], [6.5, 13.0, 19.5], [7.0, 14.0, 21.0], [7.0, 14.0, 21.0], [7.5, 15.0, 22.5], [8.0, 16.0, 24.0]] + @test isapprox(DV,D_true, atol=eps_level) end @testset "Vector of Matrix Float64 data, weighting=:none" begin @@ -2617,17 +2628,23 @@ end # Quads DF = [i.*[1.0 2.0 3.0; 4.0 5.0 6.0] for i in eachindex(Fq)] # Matrix data for each face DV = simplex2vertexdata(Fq,DF,Vq; weighting=:area) - @test isapprox(DV,[[12.0 24.0 36.0; 48.0 60.0 72.0], [9.666666666666668 19.333333333333336 28.999999999999996; 38.66666666666667 48.33333333333333 57.99999999999999], [12.666666666666666 25.333333333333332 38.0; 50.666666666666664 63.33333333333333 76.0], [15.666666666666664 31.33333333333333 47.0; 62.66666666666666 78.33333333333333 94.0], [13.0 26.0 39.0; 52.0 65.0 78.0], [10.0 20.0 30.0; 40.0 50.0 60.0], [12.333333333333332 24.666666666666664 37.0; 49.33333333333333 61.666666666666664 74.0], [14.666666666666666 29.333333333333332 44.0; 58.666666666666664 73.33333333333334 88.0], [6.500000000000001 13.000000000000002 19.5; 26.000000000000004 32.5 39.0], [11.5 23.0 34.5; 46.0 57.5 69.0], [8.5 17.0 25.500000000000004; 34.0 42.5 51.00000000000001], [10.0 20.0 30.0; 40.0 50.0 60.0], [14.0 28.0 42.0; 56.0 70.0 84.0], [9.0 18.0 27.0; 36.0 44.99999999999999 54.0], [12.5 25.0 37.50000000000001; 50.0 62.5 75.00000000000001], [16.5 33.0 49.5; 66.0 82.5 99.0], [20.5 41.0 61.5; 82.0 102.50000000000001 123.0], [16.5 33.0 49.5; 66.0 82.5 99.0], [11.5 23.0 34.5; 46.0 57.5 69.0], [13.000000000000002 26.000000000000004 39.0; 52.00000000000001 65.0 78.0], [2.5 5.0 7.5; 10.0 12.5 15.0], [6.500000000000001 13.000000000000002 19.5; 26.000000000000004 32.5 39.0], [10.5 21.0 31.5; 42.0 52.49999999999999 63.0], [14.5 29.0 43.5; 58.0 72.49999999999999 87.0], [18.5 37.0 55.5; 74.0 92.5 111.0], [22.5 45.0 67.5; 90.0 112.50000000000001 135.0]], atol=eps_level) - + ind = round.(Int,range(1,length(DV),6)) + D_true = [[13.333333333333332 26.666666666666664 40.0; 53.33333333333333 66.66666666666667 80.0], [9.0 18.0 27.0; 36.0 45.0 54.0], [11.5 23.0 34.5; 46.0 57.5 69.0], [13.5 27.0 40.5; 54.0 67.5 81.0], [10.0 20.0 30.0; 40.0 49.99999999999999 60.0], [15.0 30.0 45.0; 60.0 75.0 90.0]] + @test isapprox(DV[ind],D_true, atol=eps_level) + # Triangles DF = [i.*[1.0 2.0 3.0; 4.0 5.0 6.0] for i in eachindex(Ft)] # Matrix data for each face DV = simplex2vertexdata(Ft,DF,Vt; weighting=:none) - @test isapprox(DV,[[10.333333333333334 20.666666666666668 31.0; 41.333333333333336 51.666666666666664 62.0], [11.333333333333334 22.666666666666668 34.0; 45.333333333333336 56.666666666666664 68.0], [13.333333333333334 26.666666666666668 40.0; 53.333333333333336 66.66666666666667 80.0], [7.0 14.0 21.0; 28.0 35.0 42.0], [6.833333333333333 13.666666666666666 20.5; 27.333333333333332 34.166666666666664 41.0], [7.166666666666667 14.333333333333334 21.5; 28.666666666666668 35.833333333333336 43.0], [8.166666666666666 16.333333333333332 24.5; 32.666666666666664 40.833333333333336 49.0], [7.666666666666667 15.333333333333334 23.0; 30.666666666666668 38.333333333333336 46.0], [8.666666666666666 17.333333333333332 26.0; 34.666666666666664 43.333333333333336 52.0], [8.5 17.0 25.5; 34.0 42.5 51.0]], atol=eps_level) + ind = round.(Int,range(1,length(DV),6)) + D_true = [[10.333333333333334 20.666666666666668 31.0; 41.333333333333336 51.666666666666664 62.0], [13.333333333333334 26.666666666666668 40.0; 53.333333333333336 66.66666666666667 80.0], [6.833333333333333 13.666666666666666 20.5; 27.333333333333332 34.166666666666664 41.0], [7.166666666666667 14.333333333333334 21.5; 28.666666666666668 35.833333333333336 43.0], [7.666666666666667 15.333333333333334 23.0; 30.666666666666668 38.333333333333336 46.0], [8.5 17.0 25.5; 34.0 42.5 51.0]] + @test isapprox(DV[ind],D_true, atol=eps_level) # Hexahedral elements DF = [i.*[1.0 2.0 3.0; 4.0 5.0 6.0] for i in eachindex(Eh)] # Matrix data for each face DV = simplex2vertexdata(Eh,DF,Vh; weighting=:none) - @test isapprox(DV,[[1.0 2.0 3.0; 4.0 5.0 6.0], [1.5 3.0 4.5; 6.0 7.5 9.0], [2.0 4.0 6.0; 8.0 10.0 12.0], [2.0 4.0 6.0; 8.0 10.0 12.0], [2.5 5.0 7.5; 10.0 12.5 15.0], [3.0 6.0 9.0; 12.0 15.0 18.0], [3.0 6.0 9.0; 12.0 15.0 18.0], [3.5 7.0 10.5; 14.0 17.5 21.0], [4.0 8.0 12.0; 16.0 20.0 24.0], [3.0 6.0 9.0; 12.0 15.0 18.0], [3.5 7.0 10.5; 14.0 17.5 21.0], [4.0 8.0 12.0; 16.0 20.0 24.0], [4.0 8.0 12.0; 16.0 20.0 24.0], [4.5 9.0 13.5; 18.0 22.5 27.0], [5.0 10.0 15.0; 20.0 25.0 30.0], [5.0 10.0 15.0; 20.0 25.0 30.0], [5.5 11.0 16.5; 22.0 27.5 33.0], [6.0 12.0 18.0; 24.0 30.0 36.0], [5.0 10.0 15.0; 20.0 25.0 30.0], [5.5 11.0 16.5; 22.0 27.5 33.0], [6.0 12.0 18.0; 24.0 30.0 36.0], [6.0 12.0 18.0; 24.0 30.0 36.0], [6.5 13.0 19.5; 26.0 32.5 39.0], [7.0 14.0 21.0; 28.0 35.0 42.0], [7.0 14.0 21.0; 28.0 35.0 42.0], [7.5 15.0 22.5; 30.0 37.5 45.0], [8.0 16.0 24.0; 32.0 40.0 48.0]], atol=eps_level) + ind = round.(Int,range(1,length(DV),6)) + D_true = [[1.0 2.0 3.0; 4.0 5.0 6.0], [3.0 6.0 9.0; 12.0 15.0 18.0], [3.5 7.0 10.5; 14.0 17.5 21.0], [5.5 11.0 16.5; 22.0 27.5 33.0], [6.0 12.0 18.0; 24.0 30.0 36.0], [8.0 16.0 24.0; 32.0 40.0 48.0]] + @test isapprox(DV[ind],D_true, atol=eps_level) end end @@ -2663,17 +2680,23 @@ end # Quads DV = collect(Float64,1:length(Vq)) # Vertex data (here node number) DF = vertex2simplexdata(Fq,DV) - @test isapprox(DF,[12.75, 11.5, 14.25, 16.0, 13.25, 12.75, 12.75, 12.75, 14.25, 13.75, 12.25, 12.75, 14.25, 14.75, 14.25, 13.75, 14.5, 15.0, 16.25, 15.75, 16.0, 15.5, 16.25, 16.75], atol=eps_level) + ind = round.(Int,range(1,length(DF),6)) + D_true = [12.75, 16.0, 14.75, 12.25, 16.0, 16.75] + @test isapprox(DF[ind],D_true, atol=eps_level) # Triangles DV = collect(Float64,1:length(Vt)) # Vertex data (here node number) DF = vertex2simplexdata(Ft,DV) - @test isapprox(DF,[8.0, 6.333333333333333, 7.333333333333333, 8.333333333333334, 5.333333333333333, 6.0, 5.666666666666667, 6.333333333333333, 5.333333333333333, 4.333333333333333, 6.333333333333333, 6.333333333333333, 7.333333333333333, 4.666666666666667, 5.666666666666667, 6.666666666666667], atol=eps_level) + ind = round.(Int,range(1,length(DF),6)) + D_true = [8.0, 8.333333333333334, 5.666666666666667, 4.333333333333333, 7.333333333333333, 6.666666666666667] + @test isapprox(DF[ind],D_true, atol=eps_level) # Hexahedral elements DV = collect(Float64,1:length(Vh)) # Vertex data (here node number) DF = vertex2simplexdata(Eh,DV) - @test isapprox(DF,[7.5, 8.5, 10.5, 11.5, 16.5, 17.5, 19.5, 20.5], atol=eps_level) + ind = round.(Int,range(1,length(DF),6)) + D_true = [7.5, 8.5, 11.5, 16.5, 19.5, 20.5] + @test isapprox(DF[ind],D_true, atol=eps_level) end @testset "Vector of Vector Float64 data" begin @@ -2686,6 +2709,7 @@ end # Quads DV = [i.*[1.0,2.0,π] for i in eachindex(Vq)] # Vector data for each vertex DF = vertex2simplexdata(Fq,DV) + @testset "simplexcenter" begin eps_level = 1e-4 F, V = geosphere(2, 1.0) @@ -2695,17 +2719,23 @@ end @test length(VC) == length(F) @test isapprox(VC[1:30:end], Point3{Float64}[[-0.3504874080794224, 0.0, -0.9175879469807824], [0.9252211181650858, 0.17425248968910703, -0.2898716471939399], [-0.17425248968910703, -0.2898716471939399, -0.9252211181650858], [0.870241439674047, 0.4295322227262335, -0.15715894749713352], [-0.0876218520198556, 0.14524212567637496, -0.9709982913596904], [0.9709982913596904, 0.0876218520198556, 0.14524212567637496], [-0.5759258522984322, -0.1565463433383235, -0.7844827600958122], [0.7844827600958122, 0.5759258522984322, -0.1565463433383235], [-0.4295322227262335, -0.15715894749713352, -0.870241439674047], [0.8917525488507145, 0.08663063766925146, -0.41096206852816675], [-0.08663063766925146, -0.41096206852816675, -0.8917525488507145]], atol=eps_level) end - @test isapprox(DF,[[12.75, 25.5, 40.05530633326987], [11.5, 23.0, 36.12831551628262], [14.25, 28.5, 44.76769531365455], [16.0, 32.0, 50.26548245743669], [13.25, 26.5, 41.62610266006476], [12.75, 25.5, 40.05530633326986], [12.75, 25.5, 40.055306333269854], [12.75, 25.5, 40.05530633326986], [14.25, 28.5, 44.76769531365455], [13.75, 27.5, 43.19689898685965], [12.25, 24.5, 38.48451000647496], [12.75, 25.5, 40.05530633326986], [14.25, 28.5, 44.767695313654556], [14.75, 29.5, 46.33849164044945], [14.25, 28.5, 44.767695313654556], [13.75, 27.5, 43.19689898685965], [14.5, 29.0, 45.553093477052], [15.0, 30.0, 47.1238898038469], [16.25, 32.5, 51.050880620834135], [15.75, 31.5, 49.48008429403924], [16.0, 32.0, 50.26548245743669], [15.5, 31.0, 48.69468613064179], [16.25, 32.5, 51.05088062083414], [16.75, 33.5, 52.62167694762904]], atol=eps_level) + ind = round.(Int,range(1,length(DF),6)) + D_true = [[12.75, 25.5, 40.05530633326987], [16.0, 32.0, 50.26548245743669], [14.75, 29.5, 46.33849164044945], [12.25, 24.5, 38.48451000647496], [16.0, 32.0, 50.26548245743669], [16.75, 33.5, 52.62167694762904]] + @test isapprox(DF[ind],D_true, atol=eps_level) # Triangles DV = [i.*[1.0,2.0,π] for i in eachindex(Vt)] # Vector data for each vertex DF = vertex2simplexdata(Ft,DV) - @test isapprox(DF,[[8.0, 16.0, 25.132741228718345], [6.333333333333333, 12.666666666666666, 19.896753472735355], [7.333333333333333, 14.666666666666666, 23.03834612632515], [8.333333333333334, 16.666666666666668, 26.179938779914945], [5.333333333333333, 10.666666666666666, 16.755160819145562], [6.0, 12.0, 18.84955592153876], [5.666666666666667, 11.333333333333334, 17.802358370342162], [6.333333333333333, 12.666666666666666, 19.896753472735355], [5.333333333333333, 10.666666666666666, 16.755160819145562], [4.333333333333333, 8.666666666666666, 13.613568165555769], [6.333333333333333, 12.666666666666666, 19.896753472735355], [6.333333333333333, 12.666666666666666, 19.896753472735355], [7.333333333333333, 14.666666666666666, 23.03834612632515], [4.666666666666667, 9.333333333333334, 14.66076571675237], [5.666666666666667, 11.333333333333334, 17.80235837034216], [6.666666666666667, 13.333333333333334, 20.943951023931955]], atol=eps_level) + ind = round.(Int,range(1,length(DF),6)) + D_true = [[8.0, 16.0, 25.132741228718345], [8.333333333333334, 16.666666666666668, 26.179938779914945], [5.666666666666667, 11.333333333333334, 17.802358370342162], [4.333333333333333, 8.666666666666666, 13.613568165555769], [7.333333333333333, 14.666666666666666, 23.03834612632515], [6.666666666666667, 13.333333333333334, 20.943951023931955]] + @test isapprox(DF[ind],D_true, atol=eps_level) # Hexahedral elements DV = [i.*[1.0,2.0,π] for i in eachindex(Vh)] # Vector data for each face DF = vertex2simplexdata(Ft,DV) - @test isapprox(DF,[[8.0, 16.0, 25.132741228718345], [6.333333333333333, 12.666666666666666, 19.896753472735355], [7.333333333333333, 14.666666666666666, 23.03834612632515], [8.333333333333334, 16.666666666666668, 26.179938779914945], [5.333333333333333, 10.666666666666666, 16.755160819145562], [6.0, 12.0, 18.84955592153876], [5.666666666666667, 11.333333333333334, 17.802358370342162], [6.333333333333333, 12.666666666666666, 19.896753472735355], [5.333333333333333, 10.666666666666666, 16.755160819145562], [4.333333333333333, 8.666666666666666, 13.613568165555769], [6.333333333333333, 12.666666666666666, 19.896753472735355], [6.333333333333333, 12.666666666666666, 19.896753472735355], [7.333333333333333, 14.666666666666666, 23.03834612632515], [4.666666666666667, 9.333333333333334, 14.66076571675237], [5.666666666666667, 11.333333333333334, 17.80235837034216], [6.666666666666667, 13.333333333333334, 20.943951023931955]], atol=eps_level) + ind = round.(Int,range(1,length(DF),6)) + D_true = [[8.0, 16.0, 25.132741228718345], [8.333333333333334, 16.666666666666668, 26.179938779914945], [5.666666666666667, 11.333333333333334, 17.802358370342162], [4.333333333333333, 8.666666666666666, 13.613568165555769], [7.333333333333333, 14.666666666666666, 23.03834612632515], [6.666666666666667, 13.333333333333334, 20.943951023931955]] + @test isapprox(DF[ind],D_true, atol=eps_level) end @testset "Vector of Matrix Float64 data" begin @@ -2718,17 +2748,25 @@ end # Quads DV = [i.*[1.0 2.0 3.0; 4.0 5.0 6.0] for i in eachindex(Vq)] # Matrix data for each vertex DF = vertex2simplexdata(Fq,DV) - @test isapprox(DF,[[12.75 25.5 38.25; 51.0 63.75 76.5], [11.5 23.0 34.5; 46.0 57.5 69.0], [14.25 28.5 42.75; 57.0 71.25 85.5], [16.0 32.0 48.0; 64.0 80.0 96.0], [13.25 26.5 39.75; 53.0 66.25 79.5], [12.75 25.5 38.25; 51.0 63.75 76.5], [12.75 25.5 38.25; 51.0 63.75 76.5], [12.75 25.5 38.25; 51.0 63.75 76.5], [14.25 28.5 42.75; 57.0 71.25 85.5], [13.75 27.5 41.25; 55.0 68.75 82.5], [12.25 24.5 36.75; 49.0 61.25 73.5], [12.75 25.5 38.25; 51.0 63.75 76.5], [14.25 28.5 42.75; 57.0 71.25 85.5], [14.75 29.5 44.25; 59.0 73.75 88.5], [14.25 28.5 42.75; 57.0 71.25 85.5], [13.75 27.5 41.25; 55.0 68.75 82.5], [14.5 29.0 43.5; 58.0 72.5 87.0], [15.0 30.0 45.0; 60.0 75.0 90.0], [16.25 32.5 48.75; 65.0 81.25 97.5], [15.75 31.5 47.25; 63.0 78.75 94.5], [16.0 32.0 48.0; 64.0 80.0 96.0], [15.5 31.0 46.5; 62.0 77.5 93.0], [16.25 32.5 48.75; 65.0 81.25 97.5], [16.75 33.5 50.25; 67.0 83.75 100.5]], atol=eps_level) + ind = round.(Int,range(1,length(DF),6)) + D_true = [[12.75 25.5 38.25; 51.0 63.75 76.5], [16.0 32.0 48.0; 64.0 80.0 96.0], [14.75 29.5 44.25; 59.0 73.75 88.5], [12.25 24.5 36.75; 49.0 61.25 73.5], [16.0 32.0 48.0; 64.0 80.0 96.0], [16.75 33.5 50.25; 67.0 83.75 100.5]] + @test isapprox(DF[ind],D_true, atol=eps_level) # Triangles DV = [i.*[1.0 2.0 3.0; 4.0 5.0 6.0] for i in eachindex(Vt)] # Matrix data for each vertex DF = vertex2simplexdata(Ft,DV) - @test isapprox(DF,[[8.0 16.0 24.0; 32.0 40.0 48.0], [6.333333333333333 12.666666666666666 19.0; 25.333333333333332 31.666666666666668 38.0], [7.333333333333333 14.666666666666666 22.0; 29.333333333333332 36.666666666666664 44.0], [8.333333333333334 16.666666666666668 25.0; 33.333333333333336 41.666666666666664 50.0], [5.333333333333333 10.666666666666666 16.0; 21.333333333333332 26.666666666666668 32.0], [6.0 12.0 18.0; 24.0 30.0 36.0], [5.666666666666667 11.333333333333334 17.0; 22.666666666666668 28.333333333333332 34.0], [6.333333333333333 12.666666666666666 19.0; 25.333333333333332 31.666666666666668 38.0], [5.333333333333333 10.666666666666666 16.0; 21.333333333333332 26.666666666666668 32.0], [4.333333333333333 8.666666666666666 13.0; 17.333333333333332 21.666666666666668 26.0], [6.333333333333333 12.666666666666666 19.0; 25.333333333333332 31.666666666666668 38.0], [6.333333333333333 12.666666666666666 19.0; 25.333333333333332 31.666666666666668 38.0], [7.333333333333333 14.666666666666666 22.0; 29.333333333333332 36.666666666666664 44.0], [4.666666666666667 9.333333333333334 14.0; 18.666666666666668 23.333333333333332 28.0], [5.666666666666667 11.333333333333334 17.0; 22.666666666666668 28.333333333333332 34.0], [6.666666666666667 13.333333333333334 20.0; 26.666666666666668 33.333333333333336 40.0]], atol=eps_level) + ind = round.(Int,range(1,length(DF),6)) + D_true = [[8.0 16.0 24.0; 32.0 40.0 48.0], [8.333333333333334 16.666666666666668 25.0; 33.333333333333336 41.666666666666664 50.0], [5.666666666666667 11.333333333333334 17.0; 22.666666666666668 28.333333333333332 34.0], [4.333333333333333 8.666666666666666 13.0; 17.333333333333332 21.666666666666668 26.0], [7.333333333333333 14.666666666666666 22.0; 29.333333333333332 36.666666666666664 44.0], [6.666666666666667 13.333333333333334 20.0; 26.666666666666668 33.333333333333336 40.0]] + @test isapprox(DF[ind],D_true, atol=eps_level) + # Hexahedral elements DV = [i.*[1.0 2.0 3.0; 4.0 5.0 6.0] for i in eachindex(Vh)] # Matrix data for each face DF = vertex2simplexdata(Eh,DV) - @test isapprox(DF,[[7.5 15.0 22.5; 30.0 37.5 45.0], [8.5 17.0 25.5; 34.0 42.5 51.0], [10.5 21.0 31.5; 42.0 52.5 63.0], [11.5 23.0 34.5; 46.0 57.5 69.0], [16.5 33.0 49.5; 66.0 82.5 99.0], [17.5 35.0 52.5; 70.0 87.5 105.0], [19.5 39.0 58.5; 78.0 97.5 117.0], [20.5 41.0 61.5; 82.0 102.5 123.0]], atol=eps_level) + ind = round.(Int,range(1,length(DF),6)) + D_true = [[7.5 15.0 22.5; 30.0 37.5 45.0], [8.5 17.0 25.5; 34.0 42.5 51.0], [11.5 23.0 34.5; 46.0 57.5 69.0], [16.5 33.0 49.5; 66.0 82.5 99.0], [19.5 39.0 58.5; 78.0 97.5 117.0], [20.5 41.0 61.5; 82.0 102.5 123.0]] + @test isapprox(DF[ind],D_true, atol=eps_level) + end end @@ -2743,7 +2781,8 @@ end @test VC isa typeof(V) @test length(VC) == length(F) - @test isapprox(VC[ind], Point3{Float64}[[-0.3504874080794224, 0.0, -0.9175879469807824], [-0.4295322227262335, 0.15715894749713352, -0.870241439674047], [-0.7866254783422401, 0.3950724390612581, -0.4435636037400384], [-0.6300791350039167, 0.29832147806366355, -0.6968609080759566], [-0.805121911181463, 0.14017131621592582, -0.5511333847440926]], atol=eps_level) + VC_true = Point{3, Float64}[[-0.3504874080794224, 0.0, -0.9175879469807824], [-0.4295322227262335, 0.15715894749713352, -0.870241439674047], [-0.7866254783422401, 0.3950724390612581, -0.4435636037400384], [-0.6300791350039167, 0.29832147806366355, -0.6968609080759566], [-0.805121911181463, 0.14017131621592582, -0.5511333847440926]] + @test isapprox(VC[ind], VC_true, atol=eps_level) end @testset "Quadrilaterals" begin @@ -2752,7 +2791,8 @@ end ind = round.(Int,range(1,length(VC),5)) @test VC isa typeof(V) @test length(VC) == length(F) - @test isapprox(VC[ind], Point3{Float64}[[-0.4802860138667546, -0.4802860138667546, -0.6949720954766154], [-0.4802860138667546, 0.4802860138667546, 0.6949720954766154], [-0.7899092719339054, -0.16747661189958585, -0.5326696383253359], [0.7899092719339054, -0.16747661189958585, 0.5326696383253359], [0.16747661189958585, -0.7899092719339054, -0.5326696383253359]], atol=eps_level) + VC_true = Point{3, Float64}[[-0.4802860138667546, -0.4802860138667546, -0.6949720954766154], [-0.532669638325336, -0.16747661189958585, -0.7899092719339054], [0.532669638325336, -0.7899092719339054, -0.16747661189958585], [0.18742110835893674, -0.9256306250953279, -0.18742110835893674], [0.16747661189958585, -0.7899092719339054, -0.532669638325336]] + @test isapprox(VC[ind], VC_true, atol=eps_level) end end