From 7e10a3a3643ecff971174090e21a2baf02d7cd40 Mon Sep 17 00:00:00 2001 From: Wei Li Date: Fri, 31 May 2024 13:02:52 +1000 Subject: [PATCH] add constructor for Linearised Nedelec space on non-conforming mesh --- src/LinearizedFESpaces.jl | 199 ++++++++++++++------------ test/LinearizedFESpacesTests.jl | 245 ++++++++++++++++++-------------- 2 files changed, 246 insertions(+), 198 deletions(-) diff --git a/src/LinearizedFESpaces.jl b/src/LinearizedFESpaces.jl index a840a60..55c6cbb 100644 --- a/src/LinearizedFESpaces.jl +++ b/src/LinearizedFESpaces.jl @@ -1,100 +1,121 @@ -function _create_ref_rule(Dc,ranks::MPIArray, num_uniform_refs) - if Dc==2 - coarse_model=CartesianDiscreteModel((0,1,0,1),(1,1)) - else - @assert Dc==3 - coarse_model=CartesianDiscreteModel((0,1,0,1,0,1),(1,1,1)) - end - new_comm=MPI.Comm_split(ranks.comm,MPI.Comm_rank(ranks.comm),0) - new_ranks=MPIArray(1,new_comm,(1,)) - model_ref=OctreeDistributedDiscreteModel(new_ranks,coarse_model,num_uniform_refs) - ref_rules=map(local_views(model_ref.dmodel)) do model - Gridap.Adaptivity.RefinementRule(Gridap.Adaptivity.GenericRefinement(), - Dc==2 ? QUAD : HEX, - get_grid(model)) - end - ref_rules.item_ref[] -end - -function _get_order(reffe::Tuple{<:Lagrangian,Any,Any}) - reffe[2][2] +function _create_ref_rule(Dc, ranks::MPIArray, num_uniform_refs) + if Dc == 2 + coarse_model = CartesianDiscreteModel((0, 1, 0, 1), (1, 1)) + else + @assert Dc == 3 + coarse_model = CartesianDiscreteModel((0, 1, 0, 1, 0, 1), (1, 1, 1)) + end + new_comm = MPI.Comm_split(ranks.comm, MPI.Comm_rank(ranks.comm), 0) + new_ranks = MPIArray(1, new_comm, (1,)) + model_ref = OctreeDistributedDiscreteModel(new_ranks, coarse_model, num_uniform_refs) + ref_rules = map(local_views(model_ref.dmodel)) do model + Gridap.Adaptivity.RefinementRule(Gridap.Adaptivity.GenericRefinement(), + Dc == 2 ? QUAD : HEX, + get_grid(model)) + end + ref_rules.item_ref[] end function _create_adaptivity_glue(model::OctreeDistributedDiscreteModel{Dc}, - ref_model::OctreeDistributedDiscreteModel{Dc}, - num_uniform_refinements) where {Dc} - ref_rule=_create_ref_rule(Dc,model.parts,num_uniform_refinements) - num_children=Gridap.Adaptivity.num_subcells(ref_rule) - cell_gids_model = get_cell_gids(model.dmodel) - cell_gids_ref_model = get_cell_gids(ref_model.dmodel) + ref_model::OctreeDistributedDiscreteModel{Dc}, + num_uniform_refinements) where {Dc} + ref_rule = _create_ref_rule(Dc, model.parts, num_uniform_refinements) + num_children = Gridap.Adaptivity.num_subcells(ref_rule) + cell_gids_model = get_cell_gids(model.dmodel) + cell_gids_ref_model = get_cell_gids(ref_model.dmodel) - n2o_cell_map,n2o_cell_to_child_id=map(partition(cell_gids_model),partition(cell_gids_ref_model)) do model_partition, ref_model_partition - num_local_cells_model=local_length(model_partition) - num_owned_cells_model=own_length(model_partition) - num_local_cells_ref_model=local_length(ref_model_partition) - n2o_cell_map=Vector{Int}(undef,num_local_cells_ref_model) - n2o_cell_to_child_id=Vector{Int}(undef,num_local_cells_ref_model) - current=1 - for i=1:num_owned_cells_model - for j=1:num_children - n2o_cell_map[current]=i - n2o_cell_to_child_id[current]=j - current+=1 - end - end - n2o_cell_map,n2o_cell_to_child_id - end |> tuple_of_arrays + n2o_cell_map, n2o_cell_to_child_id = map(partition(cell_gids_model), partition(cell_gids_ref_model)) do model_partition, ref_model_partition + num_local_cells_model = local_length(model_partition) + num_owned_cells_model = own_length(model_partition) + num_local_cells_ref_model = local_length(ref_model_partition) + n2o_cell_map = Vector{Int}(undef, num_local_cells_ref_model) + n2o_cell_to_child_id = Vector{Int}(undef, num_local_cells_ref_model) + current = 1 + for i = 1:num_owned_cells_model + for j = 1:num_children + n2o_cell_map[current] = i + n2o_cell_to_child_id[current] = j + current += 1 + end + end + n2o_cell_map, n2o_cell_to_child_id + end |> tuple_of_arrays - cache = fetch_vector_ghost_values_cache(n2o_cell_map,partition(cell_gids_ref_model)) - fetch_vector_ghost_values!(n2o_cell_map,cache) |> wait - fetch_vector_ghost_values!(n2o_cell_to_child_id,cache) |> wait + cache = fetch_vector_ghost_values_cache(n2o_cell_map, partition(cell_gids_ref_model)) + fetch_vector_ghost_values!(n2o_cell_map, cache) |> wait + fetch_vector_ghost_values!(n2o_cell_to_child_id, cache) |> wait - adaptivity_glue=map(n2o_cell_map,n2o_cell_to_child_id) do n2o_cell_map, n2o_cell_to_child_id - n2o_faces_map = [(d==Dc) ? n2o_cell_map : Int[] for d in 0:Dc] - AdaptivityGlue(n2o_faces_map,n2o_cell_to_child_id,ref_rule) - end -end + adaptivity_glue = map(n2o_cell_map, n2o_cell_to_child_id) do n2o_cell_map, n2o_cell_to_child_id + n2o_faces_map = [(d == Dc) ? n2o_cell_map : Int[] for d in 0:Dc] + AdaptivityGlue(n2o_faces_map, n2o_cell_to_child_id, ref_rule) + end +end function _setup_one_level_refined_octree_model(model::OctreeDistributedDiscreteModel{Dc,Dp}, - fmodel::OctreeDistributedDiscreteModel{Dc,Dp}, - adaptivity_glue) where {Dc,Dp} - @assert model.parts === fmodel.parts - adaptive_models = map(local_views(model), - local_views(fmodel), - adaptivity_glue) do model, fmodel, glue - Gridap.Adaptivity.AdaptedDiscreteModel(fmodel.model,model,glue) + fmodel::OctreeDistributedDiscreteModel{Dc,Dp}, + adaptivity_glue) where {Dc,Dp} + @assert model.parts === fmodel.parts + adaptive_models = map(local_views(model), + local_views(fmodel), + adaptivity_glue) do model, fmodel, glue + Gridap.Adaptivity.AdaptedDiscreteModel(fmodel.model, model, glue) + end + new_fmodel = GridapDistributed.GenericDistributedDiscreteModel(adaptive_models, get_cell_gids(fmodel)) + OctreeDistributedDiscreteModel(Dc, Dp, + model.parts, + new_fmodel, + fmodel.non_conforming_glue, + model.coarse_model, + model.ptr_pXest_connectivity, + pXest_copy(fmodel.pXest_type, fmodel.ptr_pXest), + fmodel.pXest_type, + fmodel.pXest_refinement_rule_type, + false, + fmodel) +end + +function Gridap.LinearizedFESpace(model::OctreeDistributedDiscreteModel{Dc}, + reffe::Tuple{Gridap.ReferenceFEs.Lagrangian,Any,Any}; + kwargs...) where {Dc} + lin_reffe = Gridap._linearize_reffe(reffe) + order = reffe[2][2] + @assert floor(log2(order)) == ceil(log2(order)) "The order of the Lagrangian reference FE must be a power of 2" + num_refinements = Int(log2(order)) + ref_model = model + for i = 1:num_refinements + cell_gids = get_cell_gids(ref_model.dmodel) + ref_coarse_flags = map(partition(cell_gids)) do indices + flags = Vector{Cint}(undef, local_length(indices)) + flags .= refine_flag end - new_fmodel = GridapDistributed.GenericDistributedDiscreteModel(adaptive_models,get_cell_gids(fmodel)) - OctreeDistributedDiscreteModel(Dc,Dp, - model.parts, - new_fmodel, - fmodel.non_conforming_glue, - model.coarse_model, - model.ptr_pXest_connectivity, - pXest_copy(fmodel.pXest_type,fmodel.ptr_pXest), - fmodel.pXest_type, - fmodel.pXest_refinement_rule_type, - false, - fmodel) -end + ref_model, glue = Gridap.Adaptivity.adapt(ref_model, ref_coarse_flags) + end + adaptivity_glue = _create_adaptivity_glue(model, ref_model, num_refinements) + one_level_ref_model = _setup_one_level_refined_octree_model(model, ref_model, adaptivity_glue) + Gridap.FESpace(one_level_ref_model, lin_reffe; kwargs...), one_level_ref_model +end -function Gridap.LinearizedFESpace(model::OctreeDistributedDiscreteModel{Dc}, - reffe::Tuple{Gridap.ReferenceFEs.Lagrangian,Any,Any}; - kwargs...) where {Dc} - lin_reffe=Gridap._linearize_reffe(reffe) - order=_get_order(reffe) - @assert floor(log2(order)) == ceil(log2(order)) "The order of the Lagrangian reference FE must be a power of 2" - num_refinements=Int(log2(order)) - ref_model=model - for i=1:num_refinements - cell_gids=get_cell_gids(ref_model.dmodel) - ref_coarse_flags=map(partition(cell_gids)) do indices - flags=Vector{Cint}(undef,local_length(indices)) - flags.=refine_flag - end - ref_model,glue=Gridap.Adaptivity.adapt(ref_model,ref_coarse_flags) +function Gridap.LinearizedFESpace(model::OctreeDistributedDiscreteModel{Dc}, + reffe::Tuple{Gridap.ReferenceFEs.Nedelec,Any,Any}; + kwargs...) where {Dc} + lin_reffe = (reffe[1], (reffe[2][1], 0), reffe[3]) + order = reffe[2][2] + 1 + @assert floor(log2(order)) == ceil(log2(order)) "The order of the Nedelec reference FE plus one must be a power of 2" + _build_linearized_space(model, lin_reffe, order; kwargs...) +end + +function _build_linearized_space(model, lin_reffe, order; kwargs...) + num_refinements = Int(log2(order)) + ref_model = model + for _ in 1:num_refinements + cell_gids = get_cell_gids(ref_model.dmodel) + ref_coarse_flags = map(partition(cell_gids)) do indices + flags = Vector{Cint}(undef, local_length(indices)) + flags .= refine_flag end - adaptivity_glue=_create_adaptivity_glue(model,ref_model,num_refinements) - one_level_ref_model=_setup_one_level_refined_octree_model(model,ref_model,adaptivity_glue) - Gridap.FESpace(one_level_ref_model,lin_reffe; kwargs...), one_level_ref_model -end \ No newline at end of file + ref_model, _ = Gridap.Adaptivity.adapt(ref_model, ref_coarse_flags) + end + adaptivity_glue = _create_adaptivity_glue(model, ref_model, num_refinements) + one_level_ref_model = _setup_one_level_refined_octree_model(model, ref_model, adaptivity_glue) + Gridap.FESpace(one_level_ref_model, lin_reffe; kwargs...), one_level_ref_model +end diff --git a/test/LinearizedFESpacesTests.jl b/test/LinearizedFESpacesTests.jl index 6a3a838..635ae02 100644 --- a/test/LinearizedFESpacesTests.jl +++ b/test/LinearizedFESpacesTests.jl @@ -1,117 +1,144 @@ module LinearizedFESpacesTests - using Gridap - using PartitionedArrays - using GridapDistributed - using GridapP4est - using FillArrays - using Test - using MPI - - function generate_analytical_functions(order) - u(x) = x[1]^order - f(x) = -Δ(u)(x) - u, f +using Gridap +using PartitionedArrays +using GridapDistributed +using GridapP4est +using FillArrays +using Test +using MPI + +function generate_poisson_analytical_functions(order) + u(x) = x[1]^order + f(x) = -Δ(u)(x) + u, f +end + +function generate_maxwell_analytical_functions(order) + u(x) = VectorValue(x[1]x[2]^order, x[2]x[1]^order) + f(x) = ∇(∇ ⋅ u)(x) - Δ(u)(x) + u(x) + u, f +end + +function generate_weak_residual_form(Vh) + function weak_res(uH, dΩ, f) + vh_basis = Gridap.get_fe_basis(Vh) + ∫(∇(uH) ⋅ ∇(vh_basis)) * dΩ - ∫(f * vh_basis) * dΩ end +end - function generate_weak_residual_form(Vh) - function weak_res(uH, dΩ, f) - vh_basis = Gridap.get_fe_basis(Vh) - ∫(∇(uH)⋅∇(vh_basis))*dΩ - ∫(f*vh_basis)*dΩ +function run(distribute) + # debug_logger = ConsoleLgger(stderr, Logging.Debug) + # global_logger(debug_logger); # Enable the debug logger globally + ranks = distribute(LinearIndices((MPI.Comm_size(MPI.COMM_WORLD),))) + + coarse_model = CartesianDiscreteModel((0, 1, 0, 1), (1, 1)) + model = OctreeDistributedDiscreteModel(ranks, coarse_model, 2) + + cell_gids = get_cell_gids(model.dmodel) + ref_coarse_flags = map(partition(cell_gids)) do indices + flags = Vector{Cint}(undef, local_length(indices)) + flags .= nothing_flag + for i = 1:length(indices) + if i % 2 == 0 + flags[i] = refine_flag end + end + flags end + modelH, _ = Gridap.Adaptivity.adapt(model, ref_coarse_flags) + + # Poisson tests + for order in (2, 4, 8) + degree = 2 * order + 1 + I(x) = x[1]^order + u, f = generate_poisson_analytical_functions(order) + + reffeH = ReferenceFE(lagrangian, Float64, order) + VH = TestFESpace(modelH, reffeH; + conformity=:H1, + dirichlet_tags="boundary") + UH = TrialFESpace(VH, u) + + Vh, modelh = Gridap.LinearizedFESpace(modelH, reffeH; + conformity=:H1, + dirichlet_tags="boundary") + Uh = TrialFESpace(Vh, u) + Ωh = Triangulation(modelh) + dΩh = Measure(Ωh, degree) + + # Test L2-projection + aL2dΩh(uH, vh) = ∫(uH * vh)dΩh + lL2Ωh(vh) = ∫(I * vh)dΩh + op = AffineFEOperator(aL2dΩh, lL2Ωh, UH, Vh) + Ih = solve(op) + eh = sum(∫((I - Ih) * (I - Ih))dΩh) + @test eh < 1.0e-12 + + # Test laplacian (Petrov-Galerkin) + adΩHh(UH, vh) = ∫(∇(UH) ⋅ ∇(vh))dΩh + lΩHh(vh) = ∫(f * vh)dΩh + op = AffineFEOperator(adΩHh, lΩHh, UH, Vh) + uh = solve(op) + eh = sum(∫((u - uh) * (u - uh))dΩh) + @test eh < 1.0e-12 + + weak_res = generate_weak_residual_form(Vh) + rh = weak_res(uh, dΩh, f) + rh_vec = assemble_vector(rh, Vh) + @test norm(rh_vec) < 1.0e-12 + + û(x) = sin(3.2 * x[1]^2) * cos(x[1]) + sin(4.6 * x[1]) * cos(5.2 * x[1]) + ŝ(x) = exp(x[1] / 2) + 2 + k̂(x) = 2 + sin(x[1]) + + q(x) = k̂(x) * ∇(û)(x) + f(x) = -(∇ ⋅ q)(x) + ŝ(x) * û(x) + + VH = TestFESpace(modelH, reffeH; + conformity=:H1, + dirichlet_tags="boundary") + UH = TrialFESpace(VH, û) + + Vh, modelh = Gridap.LinearizedFESpace(modelH, reffeH; + conformity=:H1, + dirichlet_tags="boundary") + + Ωh = Triangulation(modelh) + dΩh = Measure(Ωh, degree) + dv = Gridap.get_fe_basis(Vh) + a(u, v) = ∫(∇(v) ⋅ (k̂ * ∇(u)) + ŝ * u * v)dΩh + l(v) = ∫(f * v)dΩh + affine_op = AffineFEOperator(a, l, UH, Vh) + ũh = solve(affine_op) + + r(u, dv) = a(u, dv) - l(dv) + j(u, du, dv) = a(du, dv) + + tol = 1.0e-12 + op = FEOperator(r, j, UH, Vh) + r, A = Gridap.Algebra.residual_and_jacobian(op, ũh) + @test norm(r) < tol - function run(distribute) - # debug_logger = ConsoleLgger(stderr, Logging.Debug) - # global_logger(debug_logger); # Enable the debug logger globally - ranks = distribute(LinearIndices((MPI.Comm_size(MPI.COMM_WORLD),))) - - coarse_model=CartesianDiscreteModel((0,1,0,1),(1,1)) - model=OctreeDistributedDiscreteModel(ranks,coarse_model,2) - - cell_gids=get_cell_gids(model.dmodel) - ref_coarse_flags=map(partition(cell_gids)) do indices - flags=Vector{Cint}(undef,local_length(indices)) - flags.=nothing_flag - for i=1:length(indices) - if i%2==0 - flags[i]=refine_flag - end - end - flags - end - modelH,_=Gridap.Adaptivity.adapt(model,ref_coarse_flags) - - for order in (2,4,8) - degree = 2*order+1 - I(x) = x[1]^order - u,f=generate_analytical_functions(order) - - reffeH = ReferenceFE(lagrangian,Float64,order) - VH = TestFESpace(modelH,reffeH; - conformity=:H1, - dirichlet_tags="boundary") - UH = TrialFESpace(VH,u) - - Vh,modelh=Gridap.LinearizedFESpace(modelH,reffeH; - conformity=:H1, - dirichlet_tags="boundary") - Uh=TrialFESpace(Vh,u) - Ωh=Triangulation(modelh) - dΩh=Measure(Ωh,degree) - - # Test L2-projection - aL2dΩh(uH,vh)=∫(uH*vh)dΩh - lL2Ωh(vh) =∫(I*vh)dΩh - op=AffineFEOperator(aL2dΩh,lL2Ωh,UH,Vh) - Ih=solve(op) - eh=sum(∫((I-Ih)*(I-Ih))dΩh) - @test eh < 1.0e-12 - - # Test laplacian (Petrov-Galerkin) - adΩHh(UH,vh)=∫(∇(UH)⋅∇(vh))dΩh - lΩHh(vh) =∫(f*vh)dΩh - op=AffineFEOperator(adΩHh,lΩHh,UH,Vh) - uh=solve(op) - eh=sum(∫((u-uh)*(u-uh))dΩh) - @test eh < 1.0e-12 - - weak_res = generate_weak_residual_form(Vh) - rh = weak_res(uh, dΩh, f) - rh_vec = assemble_vector(rh,Vh) - @test norm(rh_vec) < 1.0e-12 - - û(x) = sin(3.2 * x[1]^2) * cos(x[1]) + sin(4.6 * x[1]) * cos(5.2 * x[1]) - ŝ(x) = exp(x[1] / 2) + 2 - k̂(x) = 2 + sin(x[1]) - - q(x) = k̂(x) * ∇(û)(x) - f(x) = -(∇ ⋅ q)(x) + ŝ(x) * û(x) - - VH = TestFESpace(modelH,reffeH; - conformity=:H1, - dirichlet_tags="boundary") - UH = TrialFESpace(VH,û) - - Vh,modelh=Gridap.LinearizedFESpace(modelH,reffeH; - conformity=:H1, - dirichlet_tags="boundary") - - Ωh = Triangulation(modelh) - dΩh = Measure(Ωh, degree) - dv = Gridap.get_fe_basis(Vh) - a(u, v) = ∫(∇(v) ⋅ (k̂ * ∇(u)) + ŝ * u * v)dΩh - l(v) = ∫(f * v)dΩh - affine_op=AffineFEOperator(a, l, UH, Vh) - ũh = solve(affine_op) - - r(u,dv)=a(u,dv)-l(dv) - j(u,du,dv)=a(du,dv) - - tol = 1.0e-12 - op = FEOperator(r,j,UH,Vh) - r,A = Gridap.Algebra.residual_and_jacobian(op,ũh) - @test norm(r) < tol - - end end + + # Maxwell tests + ### TO-DO: on 31 May 2024, order 7 is not working + for order in (1, 3) + u, f = generate_maxwell_analytical_functions(order) + space_kwargs = Dict(:conformity => :Hcurl, :dirichlet_tags => "boundary") + reffe = ReferenceFE(nedelec, Float64, order) + VH = TestFESpace(modelH, reffe; space_kwargs...) + UH = TrialFESpace(VH, u) + Vh, modelh = Gridap.LinearizedFESpace(modelH, reffe; space_kwargs...) + dΩH, dΩh = Measure(Triangulation(modelH), 2order + 1), Measure(Triangulation(modelh), order + 1) + a(u, v, dΩ) = ∫((∇ × u) ⋅ (∇ × v))dΩ + ∫(u ⋅ v)dΩ + l(v, dΩ) = ∫(f ⋅ v)dΩ + uH = solve(AffineFEOperator((u, v) -> a(u, v, dΩH), v -> l(v, dΩH), UH, VH)) + uh = solve(AffineFEOperator((u, v) -> a(u, v, dΩh), v -> l(v, dΩh), UH, Vh)) + eH, eh = uH - u, uh - u + ∇xeH, ∇xeh = ∇ × eH, ∇ × eh + @test sqrt(∑(∫(eH ⋅ eH + ∇xeH ⋅ ∇xeH)dΩH)) < 1e-7 + @test sqrt(∑(∫(eh ⋅ eh + ∇xeh ⋅ ∇xeh)dΩH)) < 1e-7 + end +end end