From d6067ef5ebf4181ac77238ac6db3da5dde360a1a Mon Sep 17 00:00:00 2001 From: "Alberto F. Martin" Date: Sat, 21 Oct 2023 23:35:59 +1100 Subject: [PATCH] H1-conforming vector-values FEs work with 1 MPI task Did not test yet with more than one --- src/FESpaces.jl | 227 ++++++++++++++---- test/PoissonNonConformingOctreeModelsTests.jl | 124 ++++++---- 2 files changed, 257 insertions(+), 94 deletions(-) diff --git a/src/FESpaces.jl b/src/FESpaces.jl index dc37615..bd3e0fe 100644 --- a/src/FESpaces.jl +++ b/src/FESpaces.jl @@ -1,28 +1,124 @@ +# function _build_constraint_coefficients_matrix_in_ref_space(Dc, reffe::Tuple{<:Lagrangian,Any,Any}) +# function _h_refined_reffe(reffe::Tuple{<:Lagrangian,Any,Any}) +# (reffe[1], (reffe[2][1], 2 * reffe[2][2]), reffe[3]) +# end +# cell_polytope = Dc == 2 ? QUAD : HEX +# basis, reffe_args, reffe_kwargs = reffe +# cell_reffe = ReferenceFE(cell_polytope, basis, reffe_args...; reffe_kwargs...) +# h_refined_reffe = _h_refined_reffe(reffe) +# basis, reffe_args, reffe_kwargs = h_refined_reffe +# cell_reffe_h_refined = ReferenceFE(cell_polytope, basis, reffe_args...; reffe_kwargs...) +# dof_basis_h_refined = Gridap.CellData.get_dof_basis(cell_reffe_h_refined) +# coarse_shape_funs = Gridap.ReferenceFEs.get_shapefuns(cell_reffe) +# ref_constraints = evaluate(dof_basis_h_refined, coarse_shape_funs) +# end + +# function _generate_face_subface_ldof_to_cell_ldof(Df,Dc,reffe::Tuple{<:Lagrangian,Any,Any}) +# cell_polytope = (Dc == 2) ? QUAD : HEX +# rr = Gridap.Adaptivity.RedRefinementRule(cell_polytope) +# order=reffe[2][2] +# Gridap.Adaptivity.get_face_subface_ldof_to_cell_ldof(rr, Tuple(fill(order,Dc)), Df) +# end + function _build_constraint_coefficients_matrix_in_ref_space(Dc, reffe::Tuple{<:Lagrangian,Any,Any}) - function _h_refined_reffe(reffe::Tuple{<:Lagrangian,Any,Any}) - (reffe[1], (reffe[2][1], 2 * reffe[2][2]), reffe[3]) - end cell_polytope = Dc == 2 ? QUAD : HEX basis, reffe_args, reffe_kwargs = reffe cell_reffe = ReferenceFE(cell_polytope, basis, reffe_args...; reffe_kwargs...) - h_refined_reffe = _h_refined_reffe(reffe) - basis, reffe_args, reffe_kwargs = h_refined_reffe - cell_reffe_h_refined = ReferenceFE(cell_polytope, basis, reffe_args...; reffe_kwargs...) - dof_basis_h_refined = Gridap.CellData.get_dof_basis(cell_reffe_h_refined) - coarse_shape_funs = Gridap.ReferenceFEs.get_shapefuns(cell_reffe) - ref_constraints = evaluate(dof_basis_h_refined, coarse_shape_funs) + + # TO-DO: How can modelH be created such that it is tailored to cell_polytope? + modelH= _generate_unit_hypercube_model(Dc) + modelh=refine(modelH,2) + + VH=TestFESpace(modelH,cell_reffe) + Vh=TestFESpace(modelh,cell_reffe) + + uH=get_fe_basis(VH) + uHh=change_domain(uH,get_triangulation(modelh),ReferenceDomain()) + σRTh=Gridap.FESpaces.get_fe_dof_basis(Vh) + ref_constraints_contribs=σRTh(uHh) + + ref_constraints = Matrix{Float64}(undef,num_free_dofs(Vh),num_free_dofs(VH)) + cell_dof_ids = get_cell_dof_ids(Vh) + cache_cell_dof_ids = array_cache(cell_dof_ids) + cache_ref_constraints_contribs = array_cache(ref_constraints_contribs) + for cell=1:length(cell_dof_ids) + current_cell_dof_ids=getindex!(cache_cell_dof_ids,cell_dof_ids,cell) + current_ref_constraints_contribs=getindex!(cache_ref_constraints_contribs,ref_constraints_contribs,cell) + ref_constraints[current_cell_dof_ids,:]=current_ref_constraints_contribs + end + ref_constraints end +function _fill_face_subface_ldof_to_cell_ldof!(face_subface_ldof_to_cell_ldof, + num_faces, + coarse_faces_to_child_ids, + face_dofs, + cells_dof_ids, + first_face) + for coarse_face_id=1:num_faces + for (subface,child_id) in enumerate(coarse_faces_to_child_ids[coarse_face_id,:]) + @debug "coarse_face_id: $(coarse_face_id), subface: $(subface), child_id: $(child_id)" + cell_dof_ids=cells_dof_ids[child_id] + for (i,dof) in enumerate(cell_dof_ids[face_dofs[first_face+coarse_face_id]]) + @debug "i: $(i), dof: $(dof)" + face_subface_ldof_to_cell_ldof[coarse_face_id][subface][i]=dof + end + end + end +end + + function _generate_face_subface_ldof_to_cell_ldof(Df,Dc,reffe::Tuple{<:Lagrangian,Any,Any}) cell_polytope = (Dc == 2) ? QUAD : HEX - rr = Gridap.Adaptivity.RedRefinementRule(cell_polytope) - order=reffe[2][2] - Gridap.Adaptivity.get_face_subface_ldof_to_cell_ldof(rr, Tuple(fill(order,Dc)), Df) + coarse_faces_to_child_ids = (Dc == 2) ? _coarse_faces_to_child_ids_2D : + _coarse_faces_to_child_ids_3D + + basis, reffe_args, reffe_kwargs = reffe + cell_reffe = ReferenceFE(cell_polytope, basis, reffe_args...; reffe_kwargs...) + + # TO-DO: How can modelH be created such that it is tailored to cell_polytope? + modelH= _generate_unit_hypercube_model(Dc) + modelh=refine(modelH,2) + Vh=TestFESpace(modelh,reffe) + cells_dof_ids=get_cell_dof_ids(Vh) + + first_face = get_offset(get_polytope(cell_reffe),Df) + face_dofs=get_face_dofs(cell_reffe) + num_dofs_x_face = length(face_dofs[first_face+1]) + + if (Df==Dc-1) # Facets + num_faces = 2*Dc + num_subfaces = 2^(Dc-1) + face_subface_ldof_to_cell_ldof=_allocate_face_subface_ldof_to_cell_ldof(num_faces,num_subfaces,num_dofs_x_face) + _fill_face_subface_ldof_to_cell_ldof!(face_subface_ldof_to_cell_ldof, + num_faces, + coarse_faces_to_child_ids, + face_dofs, + cells_dof_ids, + first_face) + face_subface_ldof_to_cell_ldof + else + @assert Df==1 && Dc==3 # Edges in 3D + num_edges=12 + num_subedges=2 + edge_subedge_ldof_to_cell_ldof=_allocate_face_subface_ldof_to_cell_ldof(num_edges,num_subedges,num_dofs_x_face) + _fill_face_subface_ldof_to_cell_ldof!(edge_subedge_ldof_to_cell_ldof, + num_edges, + _coarse_edges_to_child_ids_3D, + face_dofs, + cells_dof_ids, + first_face) + edge_subedge_ldof_to_cell_ldof + end end + + const _coarse_faces_to_child_ids_2D=[1 2; 3 4; 1 3; 2 4] const _coarse_faces_to_child_ids_3D=[1 2 3 4; 5 6 7 8; 1 2 5 6; 3 4 7 8; 1 3 5 7; 2 4 6 8; ] +const _coarse_edges_to_child_ids_3D=[1 2; 3 4; 5 6; 7 8; 1 3; 2 4; 5 7; 6 8; 1 5; 2 6; 3 7; 4 8; ] + function _generate_unit_hypercube_model(Dc) @assert Dc==2 || Dc==3 @@ -58,27 +154,19 @@ function _generate_face_subface_ldof_to_cell_ldof(Df,Dc,reffe::Tuple{<:RaviartTh modelH= _generate_unit_hypercube_model(Dc) modelh=refine(modelH,2) RTh=TestFESpace(modelh,reffe) - num_faces = 2*Dc num_subfaces = 2^(Dc-1) first_face = get_offset(get_polytope(cell_reffe),Df) - face_own_dofs=get_face_own_dofs(cell_reffe) num_dofs_x_face = length(face_own_dofs[first_face+1]) - face_subface_ldof_to_cell_ldof=_allocate_face_subface_ldof_to_cell_ldof(num_faces,num_subfaces,num_dofs_x_face) - RTh_cell_dof_ids=get_cell_dof_ids(RTh) - for coarse_face_id=1:num_faces - for (subface,child_id) in enumerate(coarse_faces_to_child_ids[coarse_face_id,:]) - @debug "coarse_face_id: $(coarse_face_id), subface: $(subface), child_id: $(child_id)" - cell_dof_ids=RTh_cell_dof_ids[child_id] - for (i,dof) in enumerate(cell_dof_ids[face_own_dofs[first_face+coarse_face_id]]) - @debug "i: $(i), dof: $(dof)" - face_subface_ldof_to_cell_ldof[coarse_face_id][subface][i]=dof - end - end - end + _fill_face_subface_ldof_to_cell_ldof!(face_subface_ldof_to_cell_ldof, + num_faces, + coarse_faces_to_child_ids, + face_own_dofs, + RTh_cell_dof_ids, + first_face) else @assert Df==1 && Dc==3 # Edges in 3D num_edges=12 @@ -230,16 +318,18 @@ function _restrict_face_dofs_to_face_dim(cell_reffe,Df) first_face_faces = Gridap.ReferenceFEs.get_faces(polytope)[first_face] face_dofs_to_face_dim = face_own_dofs[first_face_faces] touched=Dict{Int,Int}() + dofs=Int64[] current=1 for (i,face_dofs) in enumerate(face_dofs_to_face_dim) for (j,dof) in enumerate(face_dofs) - if haskey(touched,dof) - face_dofs[j]=touched[dof] - else - touched[dof]=current - face_dofs[j]=touched[dof] - current=current+1 - end + push!(dofs,dof) + end + end + sort!(dofs) + touched = Dict( dofs[i]=>i for i=1:length(dofs)) + for (i,face_dofs) in enumerate(face_dofs_to_face_dim) + for (j,dof) in enumerate(face_dofs) + face_dofs[j]=touched[dof] end end face_dofs_to_face_dim @@ -508,10 +598,51 @@ function _compute_owner_edges_pindex_and_lids( owner_edge_lid, _ = owner_edges_lids[owner_edge] owner_edges_pindex[owner_edge_lid]=pindex end - end + end owner_edges_pindex, owner_edges_lids end +using Gridap.ReferenceFEs + +function get_nodes_permutations(reffe::GenericLagrangianRefFE{GradConformity}) + p = get_polytope(reffe) + face_nodes = get_face_nodes(reffe) + dofs = get_dof_basis(reffe) + interior_nodes = dofs.nodes[face_nodes[end]] + Gridap.ReferenceFEs._compute_node_permutations(p,interior_nodes) +end + +function get_nodes_permutations(reffe::ReferenceFE,d::Integer) + p = get_polytope(reffe) + range = get_dimrange(p,d) + get_face_nodes_permutations(reffe)[range] +end + +function get_face_nodes_permutations(reffe::GenericLagrangianRefFE{GradConformity}) + nodes_permutations = get_nodes_permutations(reffe) + reffaces = reffe.reffe.metadata + _reffaces = vcat(reffaces...) + face_nodes_permutations = map(get_nodes_permutations,_reffaces) + push!(face_nodes_permutations,nodes_permutations) + face_nodes_permutations +end + +function get_face_dofs_permutations(reffe::LagrangianRefFE) + dofs = get_dof_basis(reffe) + face_nodes_permutations = get_face_nodes_permutations(reffe) + face_nodes = get_face_nodes(reffe) + face_dofs = get_face_dofs(reffe) + face_dofs_permutations = Gridap.ReferenceFEs._generate_face_own_dofs_permutations( + face_nodes_permutations, dofs.node_and_comp_to_dof, face_nodes, face_dofs) + face_dofs_permutations +end + +function get_face_dofs_permutations(reffe::ReferenceFE,d::Integer) + p = get_polytope(reffe) + range = get_dimrange(p,d) + get_face_dofs_permutations(reffe)[range] +end + function generate_constraints(dmodel::OctreeDistributedDiscreteModel{Dc}, spaces_wo_constraints, reffe, @@ -575,11 +706,10 @@ function generate_constraints(dmodel::OctreeDistributedDiscreteModel{Dc}, basis, reffe_args, reffe_kwargs = reffe cell_reffe = ReferenceFE(Dc == 2 ? QUAD : HEX, basis, reffe_args...; reffe_kwargs...) - reffe_cell = cell_reffe cell_dof_ids = get_cell_dof_ids(V) - face_own_dofs = Gridap.ReferenceFEs.get_face_own_dofs(reffe_cell) - face_dofs = Gridap.ReferenceFEs.get_face_dofs(reffe_cell) + face_own_dofs = Gridap.ReferenceFEs.get_face_own_dofs(cell_reffe) + face_dofs = Gridap.ReferenceFEs.get_face_dofs(cell_reffe) hanging_faces_owner_face_dofs = Vector{Vector{Vector{Int}}}(undef, Dc) @@ -610,7 +740,6 @@ function generate_constraints(dmodel::OctreeDistributedDiscreteModel{Dc}, end basis, reffe_args, reffe_kwargs = reffe - face_reffe = ReferenceFE(facet_polytope, basis, reffe_args...; reffe_kwargs...) pindex_to_cfvertex_to_fvertex = Gridap.ReferenceFEs.get_vertex_permutations(facet_polytope) if (Dc == 3) @@ -644,20 +773,26 @@ function generate_constraints(dmodel::OctreeDistributedDiscreteModel{Dc}, end - node_permutations = Vector{Vector{Vector{Int}}}(undef, Dc - 1) - nodes, _ = Gridap.ReferenceFEs.compute_nodes(facet_polytope, [reffe_args[2] for i = 1:Dc-1]) - node_permutations[Dc-1] = Gridap.ReferenceFEs._compute_node_permutations(facet_polytope, nodes) + # node_permutations = Vector{Vector{Vector{Int}}}(undef, Dc - 1) + # nodes, _ = Gridap.ReferenceFEs.compute_nodes(facet_polytope, [reffe_args[2] for i = 1:Dc-1]) + # node_permutations[Dc-1] = Gridap.ReferenceFEs._compute_node_permutations(facet_polytope, nodes) + # if (Dc == 3) + # nodes, _ = Gridap.ReferenceFEs.compute_nodes(edget_polytope, [reffe_args[2] for i = 1:Dc-2]) + # node_permutations[1] = Gridap.ReferenceFEs._compute_node_permutations(edget_polytope, nodes) + # end + + node_permutations = Vector{Vector{Vector{Int}}}(undef, Dc-1) + node_permutations[Dc-1] = + first(get_face_dofs_permutations(cell_reffe, Dc-1)) if (Dc == 3) - nodes, _ = Gridap.ReferenceFEs.compute_nodes(edget_polytope, [reffe_args[2] for i = 1:Dc-2]) - node_permutations[1] = Gridap.ReferenceFEs._compute_node_permutations(edget_polytope, nodes) - end + node_permutations[1] = + first(get_face_dofs_permutations(cell_reffe, 1)) + end subface_own_dofs = Vector{Vector{Vector{Int}}}(undef, Dc - 1) subface_own_dofs[Dc-1] = _restrict_face_dofs_to_face_dim(cell_reffe,Dc-1) - #Gridap.ReferenceFEs.get_face_own_dofs(face_reffe) if (Dc == 3) subface_own_dofs[1] = _restrict_face_dofs_to_face_dim(cell_reffe,1) - #Gridap.ReferenceFEs.get_face_own_dofs(edge_reffe) end if (_num_face_own_dofs(cell_reffe,0)>0) _generate_constraints!(0, diff --git a/test/PoissonNonConformingOctreeModelsTests.jl b/test/PoissonNonConformingOctreeModelsTests.jl index deb1548..9680ddb 100644 --- a/test/PoissonNonConformingOctreeModelsTests.jl +++ b/test/PoissonNonConformingOctreeModelsTests.jl @@ -11,15 +11,35 @@ module PoissonNonConformingOctreeModelsTests include("CoarseDiscreteModelsTools.jl") - function test_transfer_ops_and_redistribute(ranks,dmodel,order) + function generate_analytical_problem_functions(T::Type{Float64},order) # Define manufactured functions u(x) = x[1]+x[2]^order f(x) = -Δ(u)(x) + u,f + end + + function generate_analytical_problem_functions(T::Type{VectorValue{2,Float64}},order) + # Define manufactured functions + u(x) = VectorValue(x[1]+x[2]^order,x[1]^order+x[2]) + f(x) = -Δ(u)(x) + u,f + end + + function generate_analytical_problem_functions(T::Type{VectorValue{3,Float64}},order) + # Define manufactured functions + u(x) = VectorValue(x[1]+x[2]^order+x[3],x[1]^order+x[2]+x[3],x[1]+x[2]+x[3]^order) + f(x) = -Δ(u)(x) + u,f + end + + function test_transfer_ops_and_redistribute(ranks,dmodel,order,T::Type) + # Define manufactured functions + u,f = generate_analytical_problem_functions(T,order) degree = 2*order+1 - reffe=ReferenceFE(lagrangian,Float64,order) + reffe=ReferenceFE(lagrangian,T,order) VH=FESpace(dmodel,reffe;dirichlet_tags="boundary") UH=TrialFESpace(VH,u) - ref_coarse_flags=map(ranks,partition(get_cell_gids(dmodel.dmodel))) do rank,indices + ref_coarse_flags=map(ranks,partition(get_cell_gids(dmodel.dmodel))) do rank,indices flags=zeros(Cint,length(indices)) flags.=nothing_flag @@ -45,17 +65,17 @@ module PoissonNonConformingOctreeModelsTests dΩH = Measure(ΩH,degree) aH(u,v) = ∫( ∇(v)⊙∇(u) )*dΩH - bH(v) = ∫(v*f)*dΩH + bH(v) = ∫(v⋅f)*dΩH op = AffineFEOperator(aH,bH,UH,VH) uH = solve(op) e = u - uH # # Compute errors - el2 = sqrt(sum( ∫( e*e )*dΩH )) - eh1 = sqrt(sum( ∫( e*e + ∇(e)⋅∇(e) )*dΩH )) + el2 = sqrt(sum( ∫( e⋅e )*dΩH )) + eh1 = sqrt(sum( ∫( e⋅e + ∇(e)⊙∇(e) )*dΩH )) - tol=1e-6 + tol=1e-5 println("[SOLVE COARSE] el2 < tol: $(el2) < $(tol)") println("[SOLVE COARSE] eh1 < tol: $(eh1) < $(tol)") @assert el2 < tol @@ -66,7 +86,7 @@ module PoissonNonConformingOctreeModelsTests dΩh = Measure(Ωh,degree) ah(u,v) = ∫( ∇(v)⊙∇(u) )*dΩh - bh(v) = ∫(v*f)*dΩh + bh(v) = ∫(v⋅f)*dΩh op = AffineFEOperator(ah,bh,Uh,Vh) uh = solve(op) @@ -77,8 +97,8 @@ module PoissonNonConformingOctreeModelsTests # # Compute errors - el2 = sqrt(sum( ∫( e*e )*dΩh )) - eh1 = sqrt(sum( ∫( e*e + ∇(e)⋅∇(e) )*dΩh )) + el2 = sqrt(sum( ∫( e⋅e )*dΩh )) + eh1 = sqrt(sum( ∫( e⋅e + ∇(e)⊙∇(e) )*dΩh )) println("[SOLVE FINE] el2 < tol: $(el2) < $(tol)") println("[SOLVE FINE] eh1 < tol: $(eh1) < $(tol)") @@ -88,8 +108,7 @@ module PoissonNonConformingOctreeModelsTests # prolongation via interpolation uHh=interpolate(uH,Uh) e = uh - uHh - el2 = sqrt(sum( ∫( e*e )*dΩh )) - tol=1e-6 + el2 = sqrt(sum( ∫( e⋅e )*dΩh )) println("[INTERPOLATION] el2 < tol: $(el2) < $(tol)") @assert el2 < tol @@ -100,14 +119,14 @@ module PoissonNonConformingOctreeModelsTests oph = AffineFEOperator(ahp,lhp,Uh,Vh) uHh = solve(oph) e = uh - uHh - el2 = sqrt(sum( ∫( e*e )*dΩh )) + el2 = sqrt(sum( ∫( e⋅e )*dΩh )) println("[L2 PROJECTION] el2 < tol: $(el2) < $(tol)") @assert el2 < tol # restriction via interpolation uhH=interpolate(uh,UH) e = uH - uhH - el2 = sqrt(sum( ∫( e*e )*dΩh )) + el2 = sqrt(sum( ∫( e⋅e )*dΩh )) println("[INTERPOLATION] el2 < tol: $(el2) < $(tol)") @assert el2 < tol @@ -118,7 +137,7 @@ module PoissonNonConformingOctreeModelsTests oph = AffineFEOperator(aHp,lHp,UH,VH) uhH = solve(oph) e = uH - uhH - el2 = sqrt(sum( ∫( e*e )*dΩH )) + el2 = sqrt(sum( ∫( e⋅e )*dΩH )) fmodel_red, red_glue=GridapDistributed.redistribute(fmodel); Vhred=FESpace(fmodel_red,reffe,conformity=:H1;dirichlet_tags="boundary") @@ -128,19 +147,19 @@ module PoissonNonConformingOctreeModelsTests dΩhred = Measure(Ωhred,degree) ahred(u,v) = ∫( ∇(v)⊙∇(u) )*dΩhred - bhred(v) = ∫(v*f)*dΩhred + bhred(v) = ∫(v⋅f)*dΩhred op = AffineFEOperator(ahred,bhred,Uhred,Vhred) uhred = solve(op) e = u - uhred - el2 = sqrt(sum( ∫( e*e )*dΩhred )) + el2 = sqrt(sum( ∫( e⋅e )*dΩhred )) println("[SOLVE FINE REDISTRIBUTED] el2 < tol: $(el2) < $(tol)") @assert el2 < tol uhred2 = GridapDistributed.redistribute_fe_function(uh,Vhred,fmodel_red,red_glue) e = u - uhred2 - el2 = sqrt(sum( ∫( e*e )*dΩhred )) + el2 = sqrt(sum( ∫( e⋅e )*dΩhred )) println("[REDISTRIBUTE SOLUTION] el2 < tol: $(el2) < $(tol)") @assert el2 < tol @@ -149,9 +168,12 @@ module PoissonNonConformingOctreeModelsTests function test_refine_and_coarsen_at_once(ranks, dmodel::OctreeDistributedDiscreteModel{Dc}, - order) where Dc - u(x) = x[1]+x[2]^order - f(x) = -Δ(u)(x) + order, + T::Type) where Dc + + # Define manufactured functions + u,f = generate_analytical_problem_functions(T,order) + degree = 2*order+1 ref_coarse_flags=map(ranks,partition(get_cell_gids(dmodel.dmodel))) do rank,indices flags=zeros(Cint,length(indices)) @@ -165,7 +187,7 @@ module PoissonNonConformingOctreeModelsTests end fmodel,glue=adapt(dmodel,ref_coarse_flags); - reffe=ReferenceFE(lagrangian,Float64,order) + reffe=ReferenceFE(lagrangian,T,order) VH=FESpace(dmodel,reffe,conformity=:H1;dirichlet_tags="boundary") UH=TrialFESpace(VH,u) @@ -175,17 +197,17 @@ module PoissonNonConformingOctreeModelsTests dΩH = Measure(ΩH,degree) aH(u,v) = ∫( ∇(v)⊙∇(u) )*dΩH - bH(v) = ∫(v*f)*dΩH + bH(v) = ∫(v⋅f)*dΩH op = AffineFEOperator(aH,bH,UH,VH) uH = solve(op) e = u - uH # # Compute errors - el2 = sqrt(sum( ∫( e*e )*dΩH )) - eh1 = sqrt(sum( ∫( e*e + ∇(e)⋅∇(e) )*dΩH )) + el2 = sqrt(sum( ∫( e⋅e )*dΩH )) + eh1 = sqrt(sum( ∫( e⋅e + ∇(e)⊙∇(e) )*dΩH )) - tol=1e-7 + tol=1e-6 @assert el2 < tol @assert eh1 < tol @@ -193,66 +215,72 @@ module PoissonNonConformingOctreeModelsTests dΩh = Measure(Ωh,degree) ah(u,v) = ∫( ∇(v)⊙∇(u) )*dΩh - bh(v) = ∫(v*f)*dΩh + bh(v) = ∫(v⋅f)*dΩh op = AffineFEOperator(ah,bh,Uh,Vh) uh = solve(op) e = u - uh # # Compute errors - el2 = sqrt(sum( ∫( e*e )*dΩh )) - eh1 = sqrt(sum( ∫( e*e + ∇(e)⋅∇(e) )*dΩh )) + el2 = sqrt(sum( ∫( e⋅e )*dΩh )) + eh1 = sqrt(sum( ∫( e⋅e + ∇(e)⊙∇(e) )*dΩh )) - tol=1e-7 @assert el2 < tol @assert eh1 < tol # prolongation via interpolation uHh=interpolate(uH,Uh) e = uh - uHh - el2 = sqrt(sum( ∫( e*e )*dΩh )) - tol=1e-7 + el2 = sqrt(sum( ∫( e⋅e )*dΩh )) @assert el2 < tol end - function test_2d(ranks,order) + function test_2d(ranks,order,T::Type) coarse_model=CartesianDiscreteModel((0,1,0,1),(1,1)) dmodel=OctreeDistributedDiscreteModel(ranks,coarse_model,2) - test_refine_and_coarsen_at_once(ranks,dmodel,order) + test_refine_and_coarsen_at_once(ranks,dmodel,order,T) rdmodel=dmodel for i=1:5 - rdmodel=test_transfer_ops_and_redistribute(ranks,rdmodel,order) + rdmodel=test_transfer_ops_and_redistribute(ranks,rdmodel,order,T) end end - function test_3d(ranks,order) + function test_3d(ranks,order,T::Type) coarse_model=CartesianDiscreteModel((0,1,0,1,0,1),(1,1,1)) dmodel=OctreeDistributedDiscreteModel(ranks,coarse_model,2) - test_refine_and_coarsen_at_once(ranks,dmodel,order) + test_refine_and_coarsen_at_once(ranks,dmodel,order,T) rdmodel=dmodel for i=1:5 - rdmodel=test_transfer_ops_and_redistribute(ranks,rdmodel,order) + rdmodel=test_transfer_ops_and_redistribute(ranks,rdmodel,order,T) end end - function test(ranks,TVDc::Type{Val{Dc}}, perm, order) where Dc + function test(ranks,TVDc::Type{Val{Dc}}, perm, order,T::Type) where Dc coarse_model = setup_model(TVDc,perm) model = OctreeDistributedDiscreteModel(ranks, coarse_model, 1) - test_transfer_ops_and_redistribute(ranks,model,order) + test_transfer_ops_and_redistribute(ranks,model,order,T) end + + function _field_type(::Val{Dc}, scalar_or_vector::Symbol) where Dc + if scalar_or_vector==:scalar + Float64 + else + @assert scalar_or_vector==:vector + VectorValue{Dc,Float64} + end + end function run(distribute) - # debug_logger = ConsoleLogger(stderr, Logging.Debug) + # debug_logger = ConsoleLgger(stderr, Logging.Debug) # global_logger(debug_logger); # Enable the debug logger globally - ranks = distribute(LinearIndices((MPI.Comm_size(MPI.COMM_WORLD),))) - for Dc=2:3, perm=1:4, order=1:4 - test(ranks,Val{Dc},perm,order) + for Dc=2:3, perm=1:4, order=1:4, scalar_or_vector in (:scalar,:vector) + test(ranks,Val{Dc},perm,order,_field_type(Val{Dc}(),scalar_or_vector)) end - for order=1:2 - test_2d(ranks,order) - test_3d(ranks,order) + for order=4:4,scalar_or_vector in (:scalar,:vector) + test_2d(ranks,order,_field_type(Val{2}(),scalar_or_vector)) + test_3d(ranks,order,_field_type(Val{3}(),scalar_or_vector)) end end - end +