From ace344af846479ceac23ecc144f98fd23f48d032 Mon Sep 17 00:00:00 2001 From: Francesc Verdugo Date: Mon, 2 Oct 2023 14:05:33 +0200 Subject: [PATCH] cleanup --- .gitignore | 2 + examples/poisson.jl | 7 +- test/runtests.jl | 551 +------------------------------------------- test/unittests.jl | 320 +++++++++++++++++++++++++ 4 files changed, 328 insertions(+), 552 deletions(-) create mode 100644 test/unittests.jl diff --git a/.gitignore b/.gitignore index 20fe29d1..c8f9e02d 100644 --- a/.gitignore +++ b/.gitignore @@ -1,5 +1,7 @@ *.jl.*.cov *.jl.cov *.jl.mem +*.swp +/output /Manifest.toml /docs/build/ diff --git a/examples/poisson.jl b/examples/poisson.jl index 87880d81..2b52ccb6 100644 --- a/examples/poisson.jl +++ b/examples/poisson.jl @@ -18,6 +18,7 @@ function add_default_params(params_in,params_default) end function default_params_poisson() + outdir = mkpath(joinpath(@__DIR__,"..","output")) params = Dict{Symbol,Any}() params[:mesh] = glk.cartesian_mesh((0,1,0,1),(10,10)) params[:u] = (x) -> sum(x) @@ -27,7 +28,7 @@ function default_params_poisson() params[:neumann_tags] = String[] params[:integration_degree] = 2 params[:solve] = \ - params[:example_name] = "poisson" + params[:example_path] = joinpath(outdir,"poisson") params end @@ -45,7 +46,7 @@ function poisson(params_in) dirichlet_tags = params[:dirichlet_tags] degree = params[:integration_degree] solve = params[:solve] - example_name = params[:example_name] + example_path = params[:example_path] neumann_tags = params[:neumann_tags] # Dirichlet values @@ -340,7 +341,7 @@ function poisson(params_in) results[:nfree] = nfree results[:ncells] = ncells - vtk_grid(example_name,glk.vtk_args(mesh)...) do vtk + vtk_grid(example_path,glk.vtk_args(mesh)...) do vtk glk.vtk_physical_groups!(vtk,mesh) vtk["tag"] = node_to_tag vtk["uh"] = uh diff --git a/test/runtests.jl b/test/runtests.jl index c3aa9a62..9046d74f 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,556 +1,9 @@ module GalerkinToolkitTests -import GalerkinToolkit as glk -using StaticArrays -using WriteVTK -using SparseArrays -using LinearAlgebra -using ForwardDiff using Test -using PartitionedArrays -using Metis +@time @testset "unittests.jl" begin include("unittests.jl") end -# TODO -# before more cleanup: -# generate a high order mesh from a linear non oriented mesh -# -# visualization mesh [done] -# -# remove interpolation [done] -# remove gradient! and value! [done] -# swap order in tabulation matrix. Important!! -# topology reference_faces[end][1].boundary -# cleanups -# parametrize with Ti and Tf -# think about order if it needs to be a tuple - -#lagrangian_reference_face(geo) -#lagrangian_reference_element(geo) -#lagrangian_reference_element(geo,shape=:scaler) -#lagrangian_reference_element(geo,shape=()) -#lagrangian_reference_element(geo,shape=(2,)) - -#segment -#num_dims(segment) -#face_reference_id(segment,d) -#face_reference_id(segment.boundary,d) -#face_reference_id(segment.boundary.topology,d) -#shape_functions(segment,d) -#node_coordinates(segment,d) -# -#mesh -#num_dims(mesh) -#face_reference_id(mesh,d) -#face_reference_id(mesh.topology,d) -#physical_groups(mesh) -# -#mesh2 = set_data(mesh;physical_groups,topology) - -msh = joinpath(@__DIR__,"..","assets","demo.msh") -mesh = glk.mesh_from_gmsh(msh;complexify=false) - -np = 4 -ranks = DebugArray(LinearIndices((np,))) -# TODO inclure global "glue" outside pmesh? -pmesh = glk.partition_mesh(Metis.partition,ranks,mesh,via=:nodes) -pmesh = glk.partition_mesh(Metis.partition,ranks,mesh,via=:cells) - -map(pmesh,ranks) do mesh,rank - pvtk_grid("pdebug",glk.vtk_args(mesh)...;part=rank,nparts=np) do vtk - glk.vtk_physical_groups!(vtk,mesh) - vtk["piece"] = fill(rank,sum(glk.num_faces(mesh))) - vtk["owner"] = local_to_owner(mesh.local_nodes) - vtk["interface"] = map(colors->Int(length(colors)!=1),mesh.local_node_colors) - end -end - -domain = (1,2,1,2) -cells = (2,2) -mesh = glk.cartesian_mesh(domain,cells) - -D = glk.num_dims(mesh) -geo = glk.reference_faces(mesh,D)|>first|>glk.geometry -reffe = glk.lagrangian_reference_element(geo,order=3) - -local_dofs = glk.Object(; - num_dims=Val(2), - face_reference_id = [Int[],Int[],fill(1,glk.num_faces(mesh,2))], - reference_faces = ((),(),[reffe]) -) - -dof_glue = glk.dof_glue_from_mesh_and_local_dofs(mesh,local_dofs) - -geo = glk.unit_n_cube(Val(1)) -reffe = glk.lagrangian_reference_element(geo) -A = reffe.shape_functions.tabulation_matrix(glk.value,reffe.node_coordinates) -B = reffe.shape_functions.tabulation_matrix(ForwardDiff.gradient,reffe.node_coordinates) - -reffe = glk.lagrangian_reference_element(geo,shape=(2,),order=2,major=:node) -A = reffe.shape_functions.tabulation_matrix(glk.value,reffe.node_coordinates) -A = reffe.shape_functions.tabulation_matrix(ForwardDiff.jacobian,reffe.node_coordinates) - -reffe = glk.lagrangian_reference_element(geo,shape=(2,),major=:component) -A = reffe.shape_functions.tabulation_matrix(glk.value,reffe.node_coordinates) -A = reffe.shape_functions.tabulation_matrix(ForwardDiff.jacobian,reffe.node_coordinates) -display(reffe.node_to_dofs) - -reffe = glk.lagrangian_reference_element(geo,shape=(),major=:component) -A = reffe.shape_functions.tabulation_matrix(glk.value,reffe.node_coordinates) -A = reffe.shape_functions.tabulation_matrix(ForwardDiff.jacobian,reffe.node_coordinates) -display(reffe.node_to_dofs) - -geo = glk.unit_n_cube(Val(2)) -reffe = glk.lagrangian_reference_element(geo,order=3,shape=(2,)) -display(reffe.face_own_dofs) -display(reffe.face_own_dof_permutations) - -reffe = glk.lagrangian_reference_element(geo,shape=(2,3)) -A = reffe.shape_functions.tabulation_matrix(glk.value,reffe.node_coordinates) - -geo = glk.unit_simplex(Val(1)) -perms = glk.vertex_permutations_from_geometry(geo) -@test length(perms) == 2 -glk.vertex_permutations(geo) == perms - -geo = glk.unit_simplex(Val(3)) -reffe = glk.lagrangian_reference_face(geo) - -@show reffe -display(reffe) - -perms = glk.vertex_permutations_from_geometry(geo) -glk.vertex_permutations(geo) == perms - -geo = glk.unit_n_cube(Val(2)) -perms = glk.vertex_permutations_from_geometry(geo) -@test length(perms) == 8 -glk.vertex_permutations(geo) == perms - -quad1 = glk.lagrangian_reference_face(geo) -node_perms = glk.interior_node_permutations_from_reference_face(quad1) -@test all(map(i->all(i .!= 0),node_perms)) - -quad2 = glk.lagrangian_reference_face(geo,order=2) -node_perms = glk.interior_node_permutations_from_reference_face(quad2) -@test all(map(i->all(i .!= 0),node_perms)) - -quad3 = glk.lagrangian_reference_face(geo,order=3) -node_perms = glk.interior_node_permutations_from_reference_face(quad3) -@test all(map(i->all(i .!= 0),node_perms)) - -domain = (1,2,1,2,1,2) -cells = (2,2,2) -mesh = glk.structured_simplex_mesh_with_boundary(domain,cells) -vtk_grid("debug",glk.vtk_args(mesh)...) do vtk - glk.vtk_physical_groups!(vtk,mesh) -end - -tet = glk.unit_simplex(Val(3)) -refface = glk.lagrangian_reference_face(tet,order=2) -mesh = glk.mesh_from_reference_face(refface) -vtk_grid("debug",glk.vtk_args(mesh)...) do vtk - vtk["nodeid"] = 1:6 - glk.vtk_physical_groups!(vtk,mesh) -end - -vtk_grid("debug_boundary",glk.vtk_args(glk.boundary(refface))...) do vtk - vtk["nodeid"] = 1:6 -end - -hex = glk.unit_n_cube(Val(3)) -tri_hex = glk.simplexify_reference_geometry(hex) -vtk_grid("debug",glk.vtk_args(tri_hex)...) do vtk - glk.vtk_physical_groups!(vtk,tri_hex) -end - -order = 2 -refface = glk.lagrangian_reference_face(hex;order) -mesh = glk.simplexify_reference_face(refface) - -vtk_grid("debug",glk.vtk_args(mesh)...) |> vtk_save - -domain = (1,2,1,2,1,2) -cells = (2,2,2) -mesh = glk.cartesian_mesh(domain,cells,boundary=false,complexify=false) - -vtk_grid("cartesian",glk.vtk_args(mesh)...) do vtk - glk.vtk_physical_groups!(vtk,mesh) -end - -mesh = glk.cartesian_mesh(domain,cells,complexify=false) - -vtk_grid("cartesian",glk.vtk_args(mesh)...) do vtk - glk.vtk_physical_groups!(vtk,mesh) -end - -domain = (1,2,1,2) -cells = (2,2) -mesh = glk.cartesian_mesh(domain,cells,simplexify=false,boundary=false,complexify=false) -vismesh, visglue = glk.visualization_mesh_from_mesh(mesh) -vismesh, visglue = glk.visualization_mesh_from_mesh(mesh,order=2) -vismesh, visglue = glk.visualization_mesh_from_mesh(mesh,resolution=4) - -pointdata = zeros(glk.num_nodes(vismesh)) -for face in 1:glk.num_faces(mesh,2) - fine_nodes = visglue.face_fine_nodes[face] - for fine_node in fine_nodes - pointdata[fine_node] = face - end -end - -vtk_grid("vismesh",glk.vtk_args(vismesh,2)...) do vtk - vtk["cellid"] = visglue.parent_face - vtk["pointdata"] = pointdata -end - -mesh = glk.cartesian_mesh(domain,cells,simplexify=true,boundary=false,complexify=false) -vismesh, visglue = glk.visualization_mesh_from_mesh(mesh) -vismesh, visglue = glk.visualization_mesh_from_mesh(mesh,order=2) -vismesh, visglue = glk.visualization_mesh_from_mesh(mesh,resolution=4) - -domain = (1,2,1,2,1,2) -cells = (2,2,2) -mesh = glk.cartesian_mesh(domain,cells,boundary=false,complexify=false) - -mesh,_ = glk.complexify_mesh(mesh) - -vtk_grid("cartesian",glk.vtk_args(mesh)...) do vtk - glk.vtk_physical_groups!(vtk,mesh) -end - -mesh = glk.cartesian_mesh(domain,cells) - -vtk_grid("cartesian",glk.vtk_args(mesh)...) do vtk - glk.vtk_physical_groups!(vtk,mesh) -end - -mesh = glk.cartesian_mesh(domain,cells,simplexify=true) - -vtk_grid("cartesian",glk.vtk_args(mesh)...) do vtk - glk.vtk_physical_groups!(vtk,mesh) -end - -mesh = glk.cartesian_mesh(domain,cells,simplexify=true,boundary=false) - -vtk_grid("cartesian",glk.vtk_args(mesh)...) do vtk - glk.vtk_physical_groups!(vtk,mesh) -end - -vertex = glk.unit_n_cube(Val(0)) -vertex1 = glk.lagrangian_reference_face(vertex) - -segment = glk.unit_n_cube(Val(1)) -vtk_grid("segment",glk.vtk_args(segment.boundary)...) |> vtk_save - -segment2 = glk.lagrangian_reference_face(segment) -segment3 = glk.lagrangian_reference_face(segment,order=2) -segment4 = glk.lagrangian_reference_face(segment,order=4) - -quad = glk.unit_n_cube(Val(2)) -vtk_grid("quad",glk.vtk_args(quad.boundary)...) |> vtk_save - -quad4 = glk.lagrangian_reference_face(quad,order=1) -quad9 = glk.lagrangian_reference_face(quad,order=2) -vtk_grid("quad9",glk.vtk_args(quad9.boundary)...) |> vtk_save - -mesh_quad4 = glk.mesh_from_reference_face(quad4) - -vtk_grid("mesh_quad4",glk.vtk_args(mesh_quad4)...) |> vtk_save - -function isoparametric_poisson(mesh) - - # Manufactured solution - u(x) = sum(x) - f(x) = 0.0 - - # Dirichlet values - node_to_tag = zeros(glk.num_nodes(mesh)) - tag_to_name = ["boundary"] - glk.classify_mesh_nodes!(node_to_tag,mesh,tag_to_name) - free_and_dirichlet_nodes = glk.partition_from_mask(i->i==0,node_to_tag) - node_to_x = glk.node_coordinates(mesh) - dirichlet_nodes = last(free_and_dirichlet_nodes) - x_dirichlet = view(node_to_x,dirichlet_nodes) - u_dirichlet = u.(x_dirichlet) - - # Reference cell - d = glk.num_dims(mesh) - ref_cells = glk.reference_faces(mesh,d) - face_to_rid = glk.face_reference_id(mesh,d) - - # Integration rule - degree = 4 - integration_rules = map(ref_cells) do ref_cell - glk.quadrature(glk.geometry(ref_cell),degree) - end - w = map(glk.weights,integration_rules) - q = map(glk.coordinates,integration_rules) - - # Shape functions - ab = map(integration_rules,ref_cells) do integration_rule,ref_cell - q = glk.coordinates(integration_rule) - shape_functions = glk.shape_functions(ref_cell) - m = length(q) - n = glk.num_nodes(ref_cell) - a = glk.tabulation_matrix!(shape_functions)(ForwardDiff.gradient,zeros(eltype(q),m,n),q) - b = glk.tabulation_matrix!(shape_functions)(glk.value,zeros(m,n),q) - (a,b) - end - ∇s = map(first,ab) - s = map(last,ab) - - # Count coo entries - n_coo = 0 - face_to_nodes = glk.face_nodes(mesh,d) - ncells = glk.num_faces(mesh,d) - node_to_free_node = glk.permutation(free_and_dirichlet_nodes) - nfree = length(first(free_and_dirichlet_nodes)) - for cell in 1:ncells - nodes = face_to_nodes[cell] - for nj in nodes - if node_to_free_node[nj] > nfree - continue - end - for ni in nodes - if node_to_free_node[ni] > nfree - continue - end - n_coo += 1 - end - end - end - - # Allocate coo values - I_coo = zeros(Int32,n_coo) - J_coo = zeros(Int32,n_coo) - V_coo = zeros(Float64,n_coo) - b = zeros(nfree) - - # Allocate auxiliary buffers - ∇x = map(i->similar(i,size(i,2)),∇s) - Aes = map(i->zeros(size(i,2),size(i,2)),∇s) - fes = map(i->zeros(size(i,2)),∇s) - ues = map(i->zeros(size(i,2)),∇s) - - Tx = SVector{d,Float64} - TJ = typeof(zero(Tx)*zero(Tx)') - - # Fill coo values - i_coo = 0 - for cell in 1:ncells - rid = face_to_rid[cell] - nodes = face_to_nodes[cell] - Ae = Aes[rid] - ue = ues[rid] - fe = fes[rid] - nl = length(nodes) - ∇se = ∇s[rid] - ∇xe = ∇x[rid] - se = s[rid] - we = w[rid] - nq = length(we) - fill!(Ae,zero(eltype(Ae))) - fill!(fe,zero(eltype(Ae))) - fill!(ue,zero(eltype(Ae))) - for k in 1:nl - nk = nodes[k] - gk = node_to_free_node[nk] - if gk <= nfree - continue - end - ue[k] = u_dirichlet[gk-nfree] - end - for iq in 1:nq - Jt = zero(TJ) - xint = zero(Tx) - for k in 1:nl - x = node_to_x[nodes[k]] - ∇sqx = ∇se[iq,k] - sqx = se[iq,k] - Jt += ∇sqx*x' - xint += sqx*x - end - detJt = det(Jt) - invJt = inv(Jt) - dV = abs(detJt)*we[iq] - for k in 1:nl - ∇xe[k] = invJt*∇se[iq,k] - end - for j in 1:nl - for i in 1:nl - Ae[i,j] += ∇xe[i]⋅∇xe[j]*dV - end - end - fx = f(xint) - for k in 1:nl - fe[k] += fx*se[iq,k]*dV - end - end - for i in 1:nl - for j in 1:nl - fe[i] -= Ae[i,j]*ue[j] - end - end - for i in 1:nl - ni = nodes[i] - gi = node_to_free_node[ni] - if gi > nfree - continue - end - b[gi] += fe[i] - end - for j in 1:nl - nj = nodes[j] - gj = node_to_free_node[nj] - if gj > nfree - continue - end - for i in 1:nl - ni = nodes[i] - gi = node_to_free_node[ni] - if gi > nfree - continue - end - i_coo += 1 - I_coo[i_coo] = gi - J_coo[i_coo] = gj - V_coo[i_coo] = Ae[i,j] - end - end - end - - A = sparse(I_coo,J_coo,V_coo,nfree,nfree) - u_free = A\b - nnodes = glk.num_nodes(mesh) - uh = zeros(nnodes) - uh[first(free_and_dirichlet_nodes)] = u_free - uh[last(free_and_dirichlet_nodes)] = u_dirichlet - - eh1 = 0.0 - for cell in 1:ncells - rid = face_to_rid[cell] - nodes = face_to_nodes[cell] - ue = ues[rid] - nl = length(nodes) - ∇se = ∇s[rid] - ∇xe = ∇x[rid] - se = s[rid] - we = w[rid] - nq = length(we) - for k in 1:nl - nk = nodes[k] - ue[k] = uh[nk] - end - for iq in 1:nq - Jt = zero(TJ) - xint = zero(Tx) - for k in 1:nl - x = node_to_x[nodes[k]] - ∇sqx = ∇se[iq,k] - sqx = se[iq,k] - Jt += ∇sqx*x' - xint += sqx*x - end - detJt = det(Jt) - invJt = inv(Jt) - dV = abs(detJt)*we[iq] - for k in 1:nl - ∇xe[k] = invJt*∇se[iq,k] - end - ux = u(xint) - ∇ux = ForwardDiff.gradient(u,xint) - ∇uhx = zero(∇ux) - uhx = zero(ux) - for k in 1:nl - uek = ue[k] - ∇uhx += uek*∇xe[k] - uhx += uek*se[iq,k] - end - ∇ex = ∇ux - ∇uhx - ex = ux - uhx - eh1 += (∇ex⋅∇ex + ex*ex)*dV - end - end - eh1 = sqrt(eh1) - @test eh1 < 1.0e-10 - - vtk_grid("isoparametric_poisson",glk.vtk_args(mesh)...) do vtk - glk.vtk_physical_groups!(vtk,mesh) - vtk["tag"] = node_to_tag - vtk["uh"] = uh - end - -end - -order = 1 -segment = glk.unit_n_cube(Val(1)) -segment2 = glk.lagrangian_reference_face(segment;order) -quad = glk.unit_n_cube(Val(2)) -quad4 = glk.lagrangian_reference_face(quad;order) - -node_coordinates = SVector{2,Float64}[(0,0),(1,0),(2,0),(0,1),(1,1),(2,1),(0,2),(1,2),(2,2)] -face_nodes = [ - Vector{Int}[], - [[1,2],[2,3],[1,4],[4,7],[7,8],[8,9],[3,6],[6,9]], - [[1,2,4,5],[2,3,5,6],[4,5,7,8],[5,6,8,9]] - ] -face_reference_id = [Int[],Int[1,1,1,1,1,1,1,1],[1,1,1,1]] -reference_faces = ([],[segment2],[quad4]) -physical_groups = [ - Dict([]), - Dict(["face_1"=>[1,2],"face_2"=>[3,4],"boundary"=>[1,2,3,4,5,6,7,8]]), - Dict(["domain"=>[1,2,3,4]])] -mesh = (; - num_dims=Val(2),node_coordinates, - face_nodes,face_reference_id, - reference_faces,physical_groups) - -d = 2 -vtk_grid("mesh2",glk.vtk_args(mesh,d)...) do vtk - glk.vtk_physical_groups!(vtk,mesh,d) -end - -d = 1 -vtk_grid("mesh1",glk.vtk_args(mesh,d)...) do vtk - glk.vtk_physical_groups!(vtk,mesh,d) -end - -vtk_grid("mesh",glk.vtk_args(mesh)...) do vtk - glk.vtk_physical_groups!(vtk,mesh) -end - -isoparametric_poisson(mesh) - -new_mesh, old_to_new = glk.complexify_mesh(mesh) - -isoparametric_poisson(new_mesh) - -vtk_grid("new_mesh",glk.vtk_args(new_mesh)...) do vtk - glk.vtk_physical_groups!(vtk,new_mesh) -end - -msh = joinpath(@__DIR__,"..","assets","demo.msh") -mesh = glk.mesh_from_gmsh(msh;complexify=true) -@test glk.periodic_nodes(mesh) == ([] => []) - - -#topoloy = glk.topology(mesh) -#face_to_cells = glk.face_incidence(topology,d-1,d) -#boundary_faces = findall(cells->length(cells)==1,face_to_cells) -#face_groups = gkl.physical_groups(mesh,d-1) -#face_groups["boundary"] = boundary_faces -#isoparametric_poisson(mesh) - - -vtk_grid("mesh",glk.vtk_args(mesh)...) do vtk - glk.vtk_physical_groups!(vtk,mesh) -end - -for d in 1:3 - vtk_grid("mesh_$d",glk.vtk_args(mesh,d)...) do vtk - glk.vtk_physical_groups!(vtk,mesh,d) - end -end +@time @testset "poisson_test.jl" begin include("poisson_test.jl") end end # module diff --git a/test/unittests.jl b/test/unittests.jl new file mode 100644 index 00000000..2b476377 --- /dev/null +++ b/test/unittests.jl @@ -0,0 +1,320 @@ +module UnitTests + +import GalerkinToolkit as glk +using StaticArrays +using WriteVTK +using SparseArrays +using LinearAlgebra +using ForwardDiff +using Test +using PartitionedArrays +using Metis + + +# TODO +# before more cleanup: +# generate a high order mesh from a linear non oriented mesh +# +# visualization mesh [done] +# +# remove interpolation [done] +# remove gradient! and value! [done] +# swap order in tabulation matrix. Important!! +# topology reference_faces[end][1].boundary +# cleanups +# parametrize with Ti and Tf +# think about order if it needs to be a tuple + +#lagrangian_reference_face(geo) +#lagrangian_reference_element(geo) +#lagrangian_reference_element(geo,shape=:scaler) +#lagrangian_reference_element(geo,shape=()) +#lagrangian_reference_element(geo,shape=(2,)) + +#segment +#num_dims(segment) +#face_reference_id(segment,d) +#face_reference_id(segment.boundary,d) +#face_reference_id(segment.boundary.topology,d) +#shape_functions(segment,d) +#node_coordinates(segment,d) +# +#mesh +#num_dims(mesh) +#face_reference_id(mesh,d) +#face_reference_id(mesh.topology,d) +#physical_groups(mesh) +# +#mesh2 = set_data(mesh;physical_groups,topology) + +outdir = mkpath(joinpath(@__DIR__,"..","output")) + +msh = joinpath(@__DIR__,"..","assets","demo.msh") +mesh = glk.mesh_from_gmsh(msh;complexify=false) + +np = 4 +ranks = DebugArray(LinearIndices((np,))) +# TODO inclure global "glue" outside pmesh? +pmesh = glk.partition_mesh(Metis.partition,ranks,mesh,via=:nodes) +pmesh = glk.partition_mesh(Metis.partition,ranks,mesh,via=:cells) + +map(pmesh,ranks) do mesh,rank + pvtk_grid(joinpath(outdir,"pdebug"),glk.vtk_args(mesh)...;part=rank,nparts=np) do vtk + glk.vtk_physical_groups!(vtk,mesh) + vtk["piece"] = fill(rank,sum(glk.num_faces(mesh))) + vtk["owner"] = local_to_owner(mesh.local_nodes) + vtk["interface"] = map(colors->Int(length(colors)!=1),mesh.local_node_colors) + end +end + +domain = (1,2,1,2) +cells = (2,2) +mesh = glk.cartesian_mesh(domain,cells) + +D = glk.num_dims(mesh) +geo = glk.reference_faces(mesh,D)|>first|>glk.geometry +reffe = glk.lagrangian_reference_element(geo,order=3) + +local_dofs = glk.Object(; + num_dims=Val(2), + face_reference_id = [Int[],Int[],fill(1,glk.num_faces(mesh,2))], + reference_faces = ((),(),[reffe]) +) + +dof_glue = glk.dof_glue_from_mesh_and_local_dofs(mesh,local_dofs) + +geo = glk.unit_n_cube(Val(1)) +reffe = glk.lagrangian_reference_element(geo) +A = reffe.shape_functions.tabulation_matrix(glk.value,reffe.node_coordinates) +B = reffe.shape_functions.tabulation_matrix(ForwardDiff.gradient,reffe.node_coordinates) + +reffe = glk.lagrangian_reference_element(geo,shape=(2,),order=2,major=:node) +A = reffe.shape_functions.tabulation_matrix(glk.value,reffe.node_coordinates) +A = reffe.shape_functions.tabulation_matrix(ForwardDiff.jacobian,reffe.node_coordinates) + +reffe = glk.lagrangian_reference_element(geo,shape=(2,),major=:component) +A = reffe.shape_functions.tabulation_matrix(glk.value,reffe.node_coordinates) +A = reffe.shape_functions.tabulation_matrix(ForwardDiff.jacobian,reffe.node_coordinates) +display(reffe.node_to_dofs) + +reffe = glk.lagrangian_reference_element(geo,shape=(),major=:component) +A = reffe.shape_functions.tabulation_matrix(glk.value,reffe.node_coordinates) +A = reffe.shape_functions.tabulation_matrix(ForwardDiff.jacobian,reffe.node_coordinates) +display(reffe.node_to_dofs) + +geo = glk.unit_n_cube(Val(2)) +reffe = glk.lagrangian_reference_element(geo,order=3,shape=(2,)) +display(reffe.face_own_dofs) +display(reffe.face_own_dof_permutations) + +reffe = glk.lagrangian_reference_element(geo,shape=(2,3)) +A = reffe.shape_functions.tabulation_matrix(glk.value,reffe.node_coordinates) + +geo = glk.unit_simplex(Val(1)) +perms = glk.vertex_permutations_from_geometry(geo) +@test length(perms) == 2 +glk.vertex_permutations(geo) == perms + +geo = glk.unit_simplex(Val(3)) +reffe = glk.lagrangian_reference_face(geo) + +@show reffe +display(reffe) + +perms = glk.vertex_permutations_from_geometry(geo) +glk.vertex_permutations(geo) == perms + +geo = glk.unit_n_cube(Val(2)) +perms = glk.vertex_permutations_from_geometry(geo) +@test length(perms) == 8 +glk.vertex_permutations(geo) == perms + +quad1 = glk.lagrangian_reference_face(geo) +node_perms = glk.interior_node_permutations_from_reference_face(quad1) +@test all(map(i->all(i .!= 0),node_perms)) + +quad2 = glk.lagrangian_reference_face(geo,order=2) +node_perms = glk.interior_node_permutations_from_reference_face(quad2) +@test all(map(i->all(i .!= 0),node_perms)) + +quad3 = glk.lagrangian_reference_face(geo,order=3) +node_perms = glk.interior_node_permutations_from_reference_face(quad3) +@test all(map(i->all(i .!= 0),node_perms)) + +domain = (1,2,1,2,1,2) +cells = (2,2,2) +mesh = glk.structured_simplex_mesh_with_boundary(domain,cells) +vtk_grid(joinpath(outdir,"debug"),glk.vtk_args(mesh)...) do vtk + glk.vtk_physical_groups!(vtk,mesh) +end + +tet = glk.unit_simplex(Val(3)) +refface = glk.lagrangian_reference_face(tet,order=2) +mesh = glk.mesh_from_reference_face(refface) +vtk_grid(joinpath(outdir,"debug"),glk.vtk_args(mesh)...) do vtk + vtk["nodeid"] = 1:6 + glk.vtk_physical_groups!(vtk,mesh) +end + +vtk_grid(joinpath(outdir,"debug_boundary"),glk.vtk_args(glk.boundary(refface))...) do vtk + vtk["nodeid"] = 1:6 +end + +hex = glk.unit_n_cube(Val(3)) +tri_hex = glk.simplexify_reference_geometry(hex) +vtk_grid(joinpath(outdir,"debug"),glk.vtk_args(tri_hex)...) do vtk + glk.vtk_physical_groups!(vtk,tri_hex) +end + +order = 2 +refface = glk.lagrangian_reference_face(hex;order) +mesh = glk.simplexify_reference_face(refface) + +vtk_grid(joinpath(outdir,"debug"),glk.vtk_args(mesh)...) |> vtk_save + +domain = (1,2,1,2,1,2) +cells = (2,2,2) +mesh = glk.cartesian_mesh(domain,cells,boundary=false,complexify=false) + +vtk_grid(joinpath(outdir,"cartesian"),glk.vtk_args(mesh)...) do vtk + glk.vtk_physical_groups!(vtk,mesh) +end + +mesh = glk.cartesian_mesh(domain,cells,complexify=false) + +vtk_grid(joinpath(outdir,"cartesian"),glk.vtk_args(mesh)...) do vtk + glk.vtk_physical_groups!(vtk,mesh) +end + +domain = (1,2,1,2) +cells = (2,2) +mesh = glk.cartesian_mesh(domain,cells,simplexify=false,boundary=false,complexify=false) +vismesh, visglue = glk.visualization_mesh_from_mesh(mesh) +vismesh, visglue = glk.visualization_mesh_from_mesh(mesh,order=2) +vismesh, visglue = glk.visualization_mesh_from_mesh(mesh,resolution=4) + +pointdata = zeros(glk.num_nodes(vismesh)) +for face in 1:glk.num_faces(mesh,2) + fine_nodes = visglue.face_fine_nodes[face] + for fine_node in fine_nodes + pointdata[fine_node] = face + end +end + +vtk_grid(joinpath(outdir,"vismesh"),glk.vtk_args(vismesh,2)...) do vtk + vtk["cellid"] = visglue.parent_face + vtk["pointdata"] = pointdata +end + +mesh = glk.cartesian_mesh(domain,cells,simplexify=true,boundary=false,complexify=false) +vismesh, visglue = glk.visualization_mesh_from_mesh(mesh) +vismesh, visglue = glk.visualization_mesh_from_mesh(mesh,order=2) +vismesh, visglue = glk.visualization_mesh_from_mesh(mesh,resolution=4) + +domain = (1,2,1,2,1,2) +cells = (2,2,2) +mesh = glk.cartesian_mesh(domain,cells,boundary=false,complexify=false) + +mesh,_ = glk.complexify_mesh(mesh) + +vtk_grid(joinpath(outdir,"cartesian"),glk.vtk_args(mesh)...) do vtk + glk.vtk_physical_groups!(vtk,mesh) +end + +mesh = glk.cartesian_mesh(domain,cells) + +vtk_grid(joinpath(outdir,"cartesian"),glk.vtk_args(mesh)...) do vtk + glk.vtk_physical_groups!(vtk,mesh) +end + +mesh = glk.cartesian_mesh(domain,cells,simplexify=true) + +vtk_grid(joinpath(outdir,"cartesian"),glk.vtk_args(mesh)...) do vtk + glk.vtk_physical_groups!(vtk,mesh) +end + +mesh = glk.cartesian_mesh(domain,cells,simplexify=true,boundary=false) + +vtk_grid(joinpath(outdir,"cartesian"),glk.vtk_args(mesh)...) do vtk + glk.vtk_physical_groups!(vtk,mesh) +end + +vertex = glk.unit_n_cube(Val(0)) +vertex1 = glk.lagrangian_reference_face(vertex) + +segment = glk.unit_n_cube(Val(1)) +vtk_grid(joinpath(outdir,"segment"),glk.vtk_args(segment.boundary)...) |> vtk_save + +segment2 = glk.lagrangian_reference_face(segment) +segment3 = glk.lagrangian_reference_face(segment,order=2) +segment4 = glk.lagrangian_reference_face(segment,order=4) + +quad = glk.unit_n_cube(Val(2)) +vtk_grid(joinpath(outdir,"quad"),glk.vtk_args(quad.boundary)...) |> vtk_save + +quad4 = glk.lagrangian_reference_face(quad,order=1) +quad9 = glk.lagrangian_reference_face(quad,order=2) +vtk_grid(joinpath(outdir,"quad9"),glk.vtk_args(quad9.boundary)...) |> vtk_save + +mesh_quad4 = glk.mesh_from_reference_face(quad4) + +vtk_grid(joinpath(outdir,"mesh_quad4"),glk.vtk_args(mesh_quad4)...) |> vtk_save + +order = 1 +segment = glk.unit_n_cube(Val(1)) +segment2 = glk.lagrangian_reference_face(segment;order) +quad = glk.unit_n_cube(Val(2)) +quad4 = glk.lagrangian_reference_face(quad;order) + +node_coordinates = SVector{2,Float64}[(0,0),(1,0),(2,0),(0,1),(1,1),(2,1),(0,2),(1,2),(2,2)] +face_nodes = [ + Vector{Int}[], + [[1,2],[2,3],[1,4],[4,7],[7,8],[8,9],[3,6],[6,9]], + [[1,2,4,5],[2,3,5,6],[4,5,7,8],[5,6,8,9]] + ] +face_reference_id = [Int[],Int[1,1,1,1,1,1,1,1],[1,1,1,1]] +reference_faces = ([],[segment2],[quad4]) +physical_groups = [ + Dict([]), + Dict(["face_1"=>[1,2],"face_2"=>[3,4],"boundary"=>[1,2,3,4,5,6,7,8]]), + Dict(["domain"=>[1,2,3,4]])] +mesh = (; + num_dims=Val(2),node_coordinates, + face_nodes,face_reference_id, + reference_faces,physical_groups) + +d = 2 +vtk_grid(joinpath(outdir,"mesh2"),glk.vtk_args(mesh,d)...) do vtk + glk.vtk_physical_groups!(vtk,mesh,d) +end + +d = 1 +vtk_grid(joinpath(outdir,"mesh1"),glk.vtk_args(mesh,d)...) do vtk + glk.vtk_physical_groups!(vtk,mesh,d) +end + +vtk_grid(joinpath(outdir,"mesh"),glk.vtk_args(mesh)...) do vtk + glk.vtk_physical_groups!(vtk,mesh) +end + +new_mesh, old_to_new = glk.complexify_mesh(mesh) + +vtk_grid(joinpath(outdir,"new_mesh"),glk.vtk_args(new_mesh)...) do vtk + glk.vtk_physical_groups!(vtk,new_mesh) +end + +msh = joinpath(@__DIR__,"..","assets","demo.msh") +mesh = glk.mesh_from_gmsh(msh;complexify=true) +@test glk.periodic_nodes(mesh) == ([] => []) + +vtk_grid(joinpath(outdir,"mesh"),glk.vtk_args(mesh)...) do vtk + glk.vtk_physical_groups!(vtk,mesh) +end + +for d in 1:3 + vtk_grid(joinpath(outdir,"mesh_$d"),glk.vtk_args(mesh,d)...) do vtk + glk.vtk_physical_groups!(vtk,mesh,d) + end +end + +end # module