diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 1a10554..6031307 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -15,7 +15,7 @@ jobs: fail-fast: false matrix: version: - - '1.9' + - '1.8' os: - ubuntu-20.04 arch: diff --git a/src/FESpaces.jl b/src/FESpaces.jl index dc37615..62d4708 100644 --- a/src/FESpaces.jl +++ b/src/FESpaces.jl @@ -1,28 +1,102 @@ 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 +132,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 +296,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 +576,66 @@ 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 get_face_dofs_permutations( + reffe::Gridap.ReferenceFEs.GenericRefFE{Gridap.ReferenceFEs.RaviartThomas, Dc},Df::Integer) where Dc + first_face = get_offset(get_polytope(reffe),Df) + order = length(get_face_dofs(reffe)[first_face])-1 + nfaces=num_faces(reffe,Df) + if (Df==Dc-1) + facet_polytope = Dc == 2 ? SEGMENT : QUAD + nodes, _ = Gridap.ReferenceFEs.compute_nodes(facet_polytope, [order for i = 1:Dc-1]) + Fill(Gridap.ReferenceFEs._compute_node_permutations(facet_polytope, nodes),nfaces) + elseif (Dc == 3 && Df==1) + nodes, _ = Gridap.ReferenceFEs.compute_nodes(SEGMENT, [order for i = 1:Dc-2]) + Fill(Gridap.ReferenceFEs._compute_node_permutations(SEGMENT, nodes),nfaces) + end +end + function generate_constraints(dmodel::OctreeDistributedDiscreteModel{Dc}, spaces_wo_constraints, reffe, @@ -575,11 +699,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 +733,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) @@ -642,22 +764,19 @@ function generate_constraints(dmodel::OctreeDistributedDiscreteModel{Dc}, gridap_cell_faces[1], gridap_cell_faces[2]) 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) + + face_dofs_permutations = Vector{Vector{Vector{Int}}}(undef, Dc-1) + face_dofs_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 + face_dofs_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, @@ -673,7 +792,7 @@ function generate_constraints(dmodel::OctreeDistributedDiscreteModel{Dc}, face_own_dofs, subface_own_dofs, cell_dof_ids, - node_permutations, + face_dofs_permutations, owner_faces_pindex, owner_faces_lids, ref_constraints, @@ -697,7 +816,7 @@ function generate_constraints(dmodel::OctreeDistributedDiscreteModel{Dc}, face_own_dofs, subface_own_dofs, cell_dof_ids, - node_permutations, + face_dofs_permutations, owner_faces_pindex, owner_faces_lids, ref_constraints, @@ -720,7 +839,7 @@ function generate_constraints(dmodel::OctreeDistributedDiscreteModel{Dc}, face_own_dofs, subface_own_dofs, cell_dof_ids, - node_permutations, + face_dofs_permutations, owner_faces_pindex, owner_faces_lids, ref_constraints, diff --git a/src/OctreeDistributedDiscreteModels.jl b/src/OctreeDistributedDiscreteModels.jl index fab9db6..3a2c9ee 100644 --- a/src/OctreeDistributedDiscreteModels.jl +++ b/src/OctreeDistributedDiscreteModels.jl @@ -63,7 +63,7 @@ mutable struct OctreeDistributedDiscreteModel{Dc,Dp,A,B,C,D,E,F} <: GridapDistri if (isa(dmodel,GridapDistributed.DistributedDiscreteModel)) Gridap.Helpers.@check Dc == Gridap.Geometry.num_cell_dims(dmodel) - Gridap.Helpers.@check Dc == Gridap.Geometry.num_point_dims(dmodel) + Gridap.Helpers.@check Dp == Gridap.Geometry.num_point_dims(dmodel) end A = typeof(parts) @@ -1230,11 +1230,9 @@ function Gridap.Adaptivity.refine(model::OctreeDistributedDiscreteModel{Dc,Dp}; end end -function Gridap.Adaptivity.adapt(model::OctreeDistributedDiscreteModel{Dc,Dp}, - refinement_and_coarsening_flags::MPIArray{<:Vector}; - parts=nothing) where {Dc,Dp} - Gridap.Helpers.@notimplementedif parts!=nothing +function _refine_coarsen_balance!(model::OctreeDistributedDiscreteModel{Dc,Dp}, + refinement_and_coarsening_flags::MPIArray{<:Vector}) where {Dc,Dp} # Variables which are updated accross calls to init_fn_callback_2d current_quadrant_index_within_tree = Cint(0) @@ -1430,6 +1428,16 @@ function Gridap.Adaptivity.adapt(model::OctreeDistributedDiscreteModel{Dc,Dp}, pXest_coarsen!(Val{Dc}, ptr_new_pXest, coarsen_fn_callback_c) pXest_balance!(Val{Dc}, ptr_new_pXest) pXest_update_flags!(Val{Dc},model.ptr_pXest,ptr_new_pXest) + ptr_new_pXest +end + +function Gridap.Adaptivity.adapt(model::OctreeDistributedDiscreteModel{Dc,Dp}, + refinement_and_coarsening_flags::MPIArray{<:Vector}; + parts=nothing) where {Dc,Dp} + + Gridap.Helpers.@notimplementedif parts!=nothing + + ptr_new_pXest = _refine_coarsen_balance!(model, refinement_and_coarsening_flags) # Extract ghost and lnodes ptr_pXest_ghost = setup_pXest_ghost(Val{Dc}, ptr_new_pXest) @@ -2140,7 +2148,6 @@ function setup_non_conforming_distributed_discrete_model(::Type{Val{Dc}}, face_labeling=generate_face_labeling(parts, cell_prange, coarse_discrete_model, - grid, topology, ptr_pXest, ptr_pXest_ghost) diff --git a/src/UniformlyRefinedForestOfOctreesDiscreteModels.jl b/src/UniformlyRefinedForestOfOctreesDiscreteModels.jl index 0753543..dce8aac 100644 --- a/src/UniformlyRefinedForestOfOctreesDiscreteModels.jl +++ b/src/UniformlyRefinedForestOfOctreesDiscreteModels.jl @@ -473,7 +473,6 @@ const ITERATOR_RESTRICT_TO_INTERIOR=Cint(101) function generate_face_labeling(parts, cell_prange, coarse_discrete_model::DiscreteModel{Dc,Dp}, - grid, topology, ptr_pXest, ptr_pXest_ghost) where {Dc,Dp} @@ -501,7 +500,7 @@ function generate_face_labeling(parts, owned_trees_offset[itree+1]=owned_trees_offset[itree]+tree.quadrants.elem_count end - faces_to_entity=map(grid,topology) do grid, topology + faces_to_entity=map(topology) do topology # Iterate over corners num_vertices=Gridap.Geometry.num_faces(topology,0) vertex_to_entity=zeros(Int,num_vertices) @@ -776,8 +775,8 @@ function generate_face_labeling(parts, facet_to_entity = map(x->x[Dc] , faces_to_entity) cell_to_entity = map(x->x[Dc+1], faces_to_entity) - function cell_to_faces(grid,topology,cell_dim,face_dim) - map(grid,topology) do grid,topology + function cell_to_faces(topology,cell_dim,face_dim) + map(topology) do topology Gridap.Geometry.get_faces(topology,cell_dim,face_dim) end end @@ -787,13 +786,13 @@ function generate_face_labeling(parts, update_face_to_entity_with_ghost_data!(vertex_to_entity, cell_prange, num_faces(polytope,0), - cell_to_faces(grid,topology,Dc,0)) + cell_to_faces(topology,Dc,0)) if Dc==3 update_face_to_entity_with_ghost_data!(edget_to_entity, cell_prange, num_faces(polytope,1), - cell_to_faces(grid,topology,Dc,1)) + cell_to_faces(topology,Dc,1)) # map(edget_to_entity) do edget_to_entity # @assert all(edget_to_entity .!= 0) # end @@ -805,12 +804,12 @@ function generate_face_labeling(parts, update_face_to_entity_with_ghost_data!(facet_to_entity, cell_prange, num_faces(polytope,Dc-1), - cell_to_faces(grid,topology,Dc,Dc-1)) + cell_to_faces(topology,Dc,Dc-1)) update_face_to_entity_with_ghost_data!(cell_to_entity, cell_prange, num_faces(polytope,Dc), - cell_to_faces(grid,topology,Dc,Dc)) + cell_to_faces(topology,Dc,Dc)) # map(vertex_to_entity,facet_to_entity,cell_to_entity) do vertex_to_entity,facet_to_entity,cell_to_entity # @assert all(vertex_to_entity .!= 0) @@ -940,7 +939,6 @@ function setup_distributed_discrete_model(::Type{Val{Dc}}, face_labeling=generate_face_labeling(parts, cell_prange, coarse_discrete_model, - grid, topology, ptr_pXest, ptr_pXest_ghost) diff --git a/test/PoissonNonConformingOctreeModelsTests.jl b/test/PoissonNonConformingOctreeModelsTests.jl index deb1548..4b32b05 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,30 +147,33 @@ 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 )) + uhred = solve(op) + e = u - uhred + 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 )) - println("[REDISTRIBUTE SOLUTION] 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 )) + println("[REDISTRIBUTE SOLUTION] el2 < tol: $(el2) < $(tol)") + @assert el2 < tol - fmodel_red + fmodel_red end 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-5 @assert el2 < tol @assert eh1 < tol @@ -193,66 +215,75 @@ 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;num_amr_steps=5) 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) + for i=1:num_amr_steps + rdmodel=test_transfer_ops_and_redistribute(ranks,rdmodel,order,T) end end - function test_3d(ranks,order) + function test_3d(ranks,order,T::Type;num_amr_steps=5) 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) + for i=1:num_amr_steps + 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,) + test(ranks,Val{Dc},perm,order,_field_type(Val{Dc}(),scalar_or_vector)) + end + for Dc=2:3, perm in (1,2), order in (1,4), scalar_or_vector in (: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=2:2, scalar_or_vector in (:scalar,:vector) + test_2d(ranks,order,_field_type(Val{2}(),scalar_or_vector), num_amr_steps=5) + test_3d(ranks,order,_field_type(Val{3}(),scalar_or_vector), num_amr_steps=4) end end - end + diff --git a/test/runtests.jl b/test/runtests.jl index 31acee4..8a3f799 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -37,12 +37,14 @@ function run_tests(testdir) np = [4] extra_args = "-s 2 2 -r 2" elseif f in ["OctreeDistributedDiscreteModelsTests.jl", - "OctreeDistributedDiscreteModelsNoEnvTests.jl"] + "OctreeDistributedDiscreteModelsNoEnvTests.jl", + "AdaptivityFlagsMarkingStrategiesTests.jl"] np = [4] extra_args = "" - elseif f in ["PoissonNonConformingOctreeModelsTests.jl", - "AdaptivityFlagsMarkingStrategiesTests.jl", - "DarcyNonConformingOctreeModelsTests.jl"] + elseif f in ["DarcyNonConformingOctreeModelsTests.jl"] + np = [1,4] + extra_args = "" + elseif f in ["PoissonNonConformingOctreeModelsTests.jl"] np = [1,2,4] extra_args = "" else