diff --git a/aggregates_bounding_boxes_tools.jl b/aggregates_bounding_boxes_tools.jl index fff83b5..6f3e2bd 100644 --- a/aggregates_bounding_boxes_tools.jl +++ b/aggregates_bounding_boxes_tools.jl @@ -36,6 +36,48 @@ function setup_aggregate_to_cells(aggregates) end end aggregate_to_cells +end + +""" + Creates an array of arrays with as many entries + as interior cells that are not part of any aggegrate. + For each interior cell, the array + contains the global cell IDs of that cells in the background + model that belong to the same interior cell + + TO-DO: with efficiency in mind we may want to store this + array of arrays as a Gridap.Arrays.Table. +""" +function setup_int_nonagg_cell_to_cells(aggregates) + + size_aggregates=Dict{Int,Int}() + for (i,agg) in enumerate(aggregates) + if agg>0 + if !haskey(size_aggregates,agg) + size_aggregates[agg]=1 + else + size_aggregates[agg]+=1 + end + end + end + + touched=Dict{Int,Int}() + interior_cell_to_cells=Vector{Vector{Int}}() + current_interior_cell=1 + for (i,agg) in enumerate(aggregates) + if agg>0 + if (size_aggregates[agg]==1) + if !haskey(touched,agg) + push!(interior_cell_to_cells,[i]) + touched[agg]=current_interior_cell + current_interior_cell+=1 + else + push!(interior_cell_to_cells[touched[agg]],i) + end + end + end + end + interior_cell_to_cells end function setup_aggregates_bounding_box_model(bgmodel, aggregate_to_cells) @@ -97,6 +139,52 @@ function setup_agg_cells(aggregate_to_cells) end return agg_cells end + +function setup_int_nonagg_cells(int_nonagg_cell_to_cells) + # Generate an array with the global IDs of the cells + # that belong to the interior, yet do not belong to the + # aggregate. We will use "int_nonagg_cells" to refer to those + # cells of the background model that belong to the interior + # but are not part of any of the aggregates. Thus, all interior + # cells but not including the root cells of the aggregates. + int_nonagg_cells=Vector{Int}() + for cell in int_nonagg_cell_to_cells + append!(int_nonagg_cells,cell) + end + return int_nonagg_cells +end + +function setup_root_cells(int_cells, int_nonagg_cells) + # Generate an array with the global IDs of the cells + # that are the root cells of the aggregrates. + root_cells=Vector{Int}() + tester_int_nonagg_cells=Vector{Int}() + for cell in int_cells + if cell ∈ int_nonagg_cells + append!(tester_int_nonagg_cells,cell) + else + append!(root_cells,cell) + end + end + @assert(tester_int_nonagg_cells==int_nonagg_cells) + return root_cells +end + +function setup_cut_cells(agg_cells, root_cells) + # Generate an array with the global IDs of the cells + # that are the cut cells of the aggregrates. + cut_cells=Vector{Int}() + tester_root_cells=Vector{Int}() + for cell in agg_cells + if cell ∈ root_cells + append!(tester_root_cells,cell) + else + append!(cut_cells,cell) + end + end + @assert(tester_root_cells==root_cells) + cut_cells +end function setup_agg_cells_to_aggregate(aggregate_to_cells) # Generate an array that given the local ID of an "agg_cell" @@ -109,8 +197,52 @@ function setup_agg_cells_to_aggregate(aggregate_to_cells) end end agg_cells_to_aggregate +end + +function setup_cut_cells_in_agg_cells_to_aggregate(aggregate_to_cut_cells) + # Generate an array that given the local ID of an "agg_cell" + # returns the ID of the aggregate to which it belongs + # (i.e., flattened version of aggregate_to_cut_cells) + cut_cells_in_agg_cells_to_aggregate=Vector{Int}() + for (i,cells) in enumerate(aggregate_to_cut_cells) + for _ in cells + push!(cut_cells_in_agg_cells_to_aggregate,i) + end + end + cut_cells_in_agg_cells_to_aggregate end - + +function setup_aggregate_to_cut_cells(aggregates, root_cells) + size_aggregates=Dict{Int,Int}() + for (i,agg) in enumerate(aggregates) + if agg>0 + if !haskey(size_aggregates,agg) + size_aggregates[agg]=1 + else + size_aggregates[agg]+=1 + end + end + end + + touched=Dict{Int,Int}() + aggregate_to_cut_cells=Vector{Vector{Int}}() + current_aggregate=1 + for (i,agg) in enumerate(aggregates) + if agg>0 + if (size_aggregates[agg]>1) && i ∉ root_cells + if !haskey(touched,agg) + push!(aggregate_to_cut_cells,[i]) + touched[agg]=current_aggregate + current_aggregate+=1 + else + push!(aggregate_to_cut_cells[touched[agg]],i) + end + end + end + end + aggregate_to_cut_cells +end + function setup_aggregate_to_local_cells(aggregate_to_cells) aggregate_to_local_cells=deepcopy(aggregate_to_cells) current_local_cell=1 diff --git a/bulk_ghost_penalty_canvas_div_on_cut_cells.jl b/bulk_ghost_penalty_canvas_div_on_cut_cells.jl new file mode 100644 index 0000000..66872cd --- /dev/null +++ b/bulk_ghost_penalty_canvas_div_on_cut_cells.jl @@ -0,0 +1,248 @@ +using Gridap +using GridapEmbedded +using FillArrays +using LinearAlgebra + +include("aggregates_bounding_boxes_tools.jl") +include("bulk_ghost_penalty_stab_tools.jl") +include("fields_and_blocks_tools.jl") +include("BulkGhostPenaltyAssembleMaps.jl") + +# Manufactured solution +order = 0 +uex(x) = -VectorValue(2*x[1],2*x[2]) +pex(x) = (x[1]^2 + x[2]^2) +divuex(x) = -4.0 + +# Select geometry +R = 0.2 +geom = disk(R, x0=Point(0.5,0.5)) + +# Setup background model +n=16 +partition = (n,n) +box = get_metadata(geom) +bgmodel = CartesianDiscreteModel((0,1,0,1),partition) +dp = box.pmax - box.pmin +h = dp[1]/n + +# Cut the background model with the mesh +cutgeo = cut(bgmodel,geom) + +# Compute mapping among background model +# cut cells and interior cells +strategy = AggregateAllCutCells() +aggregates = aggregate(strategy,cutgeo) + +aggregate_to_cells=setup_aggregate_to_cells(aggregates) +aggregates_bounding_box_model= + setup_aggregates_bounding_box_model(bgmodel,aggregate_to_cells) + +# Triangulations +Ωbg = Triangulation(bgmodel) +Ωact = Triangulation(cutgeo,ACTIVE) + +# Physical domain +Ω = Triangulation(cutgeo,PHYSICAL) +degree=2*(order+1) +dΩ = Measure(Ω,degree) + +# Set up global spaces +V = FESpace(Ωact, ReferenceFE(raviart_thomas,Float64,order),conformity=:HDiv) +Q = FESpace(Ωact, ReferenceFE(lagrangian,Float64,order), conformity=:L2) +U = TrialFESpace(V) +P = TrialFESpace(Q) +Y = MultiFieldFESpace([V, Q]) +X = MultiFieldFESpace([U, P]) +dx = get_trial_fe_basis(X) +dy = get_fe_basis(Y) +du,dp = dx +dv,dq = dy + +# Aggregate cells and triangulation +agg_cells =setup_agg_cells(aggregate_to_cells) +Ωbg_agg_cells=view(Ωbg,agg_cells) + +# ref_agg_cell_to_agg_cell_map: \hat{K} -> K +ref_agg_cell_to_agg_cell_map=get_cell_map(Ωbg_agg_cells) +agg_cells_to_aggregate =setup_agg_cells_to_aggregate(aggregate_to_cells) +ref_agg_cell_to_ref_bb_map =setup_ref_agg_cell_to_ref_bb_map(aggregates_bounding_box_model, + agg_cells_to_aggregate) + +# Spaces on bounding boxes +reffeₚ_bb =ReferenceFE(lagrangian,Float64,order) # Here we MUST use a Q space (not a P space!) +Qbb=FESpace(aggregates_bounding_box_model,reffeₚ_bb,conformity=:L2) # We need a DG space to represent the L2 projection +Pbb=TrialFESpace(Qbb) +pbb=get_trial_fe_basis(Pbb) +qbb=get_fe_basis(Qbb) +reffeᵤ_bb=ReferenceFE(raviart_thomas,Float64,order) +Vbb=FESpace(aggregates_bounding_box_model,reffeᵤ_bb,conformity=:L2) +Ubb=TrialFESpace(Vbb) +ubb=get_trial_fe_basis(Ubb) +vbb=get_fe_basis(Vbb) + +# Numerical integration (Measures) +degree=2*(order+1) +dΩbg_agg_cells = Measure(Ωbg_agg_cells,degree) + +# LHS of L2 projection on bounding boxes. +aggregate_to_local_cells=setup_aggregate_to_local_cells(aggregate_to_cells) +p_lhs=set_up_bulk_ghost_penalty_lhs(aggregate_to_local_cells, + agg_cells_to_aggregate, + ref_agg_cell_to_ref_bb_map, + dΩbg_agg_cells, + pbb, + qbb) + +u_lhs=set_up_bulk_ghost_penalty_lhs(aggregate_to_local_cells, + agg_cells_to_aggregate, + ref_agg_cell_to_ref_bb_map, + dΩbg_agg_cells, + ubb, + vbb) + +# Selecting relevant global dofs ids of aggregate cells (from background mesh) +Ωbg_agg_cell_dof_ids = get_cell_dof_ids(X,Ωbg_agg_cells) +U_Ωbg_agg_cell_dof_ids = _restrict_to_block(Ωbg_agg_cell_dof_ids, 1) +P_Ωbg_agg_cell_dof_ids = _restrict_to_block(Ωbg_agg_cell_dof_ids, 2) + +# Computing local (per aggregate) dof ids +U_agg_cells_local_dof_ids= + compute_agg_cells_local_dof_ids(U_Ωbg_agg_cell_dof_ids, aggregate_to_local_cells) +P_agg_cells_local_dof_ids= + compute_agg_cells_local_dof_ids(P_Ωbg_agg_cell_dof_ids, aggregate_to_local_cells) + +# Compute global dofs ids per aggregate and reindex these +U_aggregate_dof_ids=compute_aggregate_dof_ids(U_Ωbg_agg_cell_dof_ids,aggregate_to_cells) +U_agg_cells_to_aggregate_dof_ids=lazy_map(Reindex(U_aggregate_dof_ids),agg_cells_to_aggregate) +P_aggregate_dof_ids=compute_aggregate_dof_ids(P_Ωbg_agg_cell_dof_ids,aggregate_to_cells) +P_agg_cells_to_aggregate_dof_ids=lazy_map(Reindex(P_aggregate_dof_ids),agg_cells_to_aggregate) +P_Ωbg_agg_cell_dof_ids + +# parameters +γ = 10.0 # Interior bulk-penalty stabilization parameter +h_U = 1.0 + +########################################### +### STABILIZATION ON Ωagg\Troot ### +########################################### + +## (1) FIND ALL INTERIOR CELLS WHICH ARE NOT PART OF AN AGGREGATE +int_nonagg_cell_to_cells = setup_int_nonagg_cell_to_cells(aggregates) +int_nonagg_cells = setup_int_nonagg_cells(int_nonagg_cell_to_cells) + +## (2) FIND THE INTERIOR CELLS +Ω_agg_cells = view(Ω,agg_cells) # cells in aggregates +int_cells = Ω_agg_cells.parent.b.tface_to_mface + +## (3) FIND THE INTERIOR AGGREGATE (/ROOT) CELLS AND CUT CELLS +root_cells = setup_root_cells(int_cells, int_nonagg_cells) +cut_cells = setup_cut_cells(agg_cells, root_cells) + +## (4) CUT CELLS AND MEASURE +Ωbg_cut_cells = view(Ωbg,cut_cells) # cut cells (part of aggregate) +dΩbg_cut_cells = Measure(Ωbg_cut_cells,degree) + +# (5) DOF IDS related to cut cells +Ωbg_cut_cell_dof_ids = get_cell_dof_ids(X,Ωbg_cut_cells) +U_Ωbg_cut_cell_dof_ids = _restrict_to_block(Ωbg_cut_cell_dof_ids, 1) +P_Ωbg_cut_cell_dof_ids = _restrict_to_block(Ωbg_cut_cell_dof_ids, 2) + +# (6) Defining cut_cells_to_aggregate_dof_ids +aggregate_to_cut_cells = setup_aggregate_to_cut_cells(aggregates, root_cells) +cut_cells_in_agg_cells_to_aggregate = setup_cut_cells_in_agg_cells_to_aggregate(aggregate_to_cut_cells) +U_cut_cells_to_aggregate_dof_ids=lazy_map(Reindex(U_aggregate_dof_ids),cut_cells_in_agg_cells_to_aggregate) +P_cut_cells_to_aggregate_dof_ids=lazy_map(Reindex(P_aggregate_dof_ids),cut_cells_in_agg_cells_to_aggregate) + +# (7) Compute stabilization terms for u and p +wu,ru,cu=bulk_ghost_penalty_stabilization_collect_cell_matrix_on_cut_cells(agg_cells_to_aggregate, + ref_agg_cell_to_ref_bb_map, + dΩbg_agg_cells, + dv, # Test basis + du, # Trial basis (to project) + vbb, # Bounding box space test basis + u_lhs, + U_Ωbg_cut_cell_dof_ids, + U_agg_cells_local_dof_ids, + U_cut_cells_to_aggregate_dof_ids, + γ, + dΩbg_cut_cells) + +wp,rp,cp=bulk_ghost_penalty_stabilization_collect_cell_matrix_on_cut_cells(agg_cells_to_aggregate, + ref_agg_cell_to_ref_bb_map, + dΩbg_agg_cells, + dq, # Test basis + dp, # Trial basis (to project) + qbb, # Bounding box space test basis + p_lhs, + P_Ωbg_cut_cell_dof_ids, + P_agg_cells_local_dof_ids, + P_cut_cells_to_aggregate_dof_ids, + γ, + dΩbg_cut_cells) + +## TEST WITH MANUFACTURED SOLUTIONS +a((u,p),(v,q))=∫(v⋅u+q*p+(∇⋅u)*(∇⋅v))dΩ +wrc=Gridap.FESpaces.collect_cell_matrix(X,Y,a(dx,dy)) +assem=SparseMatrixAssembler(X,Y) +A_nostab = assemble_matrix(assem, wrc) +l((v,q))=∫(v⋅uex+q*pex+(∇⋅v)*divuex)dΩ +b = assemble_vector(l, Y) + +push!(wrc[1], wu...) +push!(wrc[2], ru...) +push!(wrc[3], cu...) + +push!(wrc[1], wp...) +push!(wrc[2], rp...) +push!(wrc[3], cp...) + +A_stabup = assemble_matrix(assem, wrc) +cond_nostab = cond(Array(A_nostab)) +cond_stabup = cond(Array(A_stabup)) + +# CHECK errors when only the p and u stabilization terms are included. +global_l2_proj_dofs_upstab = A_stabup\b +xh_upstab = FEFunction(X, global_l2_proj_dofs_upstab) +uh_upstab,ph_upstab = xh_upstab +euh_upstab = uex-uh_upstab +eph_upstab = pex-ph_upstab +edivuh_upstab = divuex-(∇⋅uh_upstab) +norm_euh_upstab = sum(∫(euh_upstab⋅euh_upstab)*dΩ) +norm_eph_upstab = sum(∫(eph_upstab*eph_upstab)*dΩ) +norm_edivuh_upstab = sum(∫(edivuh_upstab⋅edivuh_upstab)*dΩ) +# No stabilization (n=16) (k=0) κ = 2.3e9 vs (k=2) κ = 2.8e29 +# k=0 (n=16) yields euh = 2.0e-28, eph = 0.00015, edivuh = 1.9e-30, κ = 9.2e5 +# k=2 (n=16) yields euh = 2.0e-18, eph = 4.0e-29, edivuh = 1.7e-19, κ = 4.4e15 + +#DIV stabilization part +wdiv, rdiv, cdiv = div_penalty_stabilization_collect_cell_matrix_on_cut_cells(agg_cells_to_aggregate, + ref_agg_cell_to_ref_bb_map, + dΩbg_agg_cells, + dv, # Test basis + du, # Trial basis (to project) + qbb, # Bounding box space test basis + p_lhs, + U_Ωbg_cut_cell_dof_ids, + U_agg_cells_local_dof_ids, + U_cut_cells_to_aggregate_dof_ids, + γ, + dΩbg_cut_cells) +push!(wrc[1], wdiv...) +push!(wrc[2], rdiv...) +push!(wrc[3], cdiv...) + +A_stab = assemble_matrix(assem, wrc) +cond_stab = cond(Array(A_stab)) + +global_l2_proj_dofs_stab = A_stab\b +xh_stab = FEFunction(X, global_l2_proj_dofs_stab) +uh_stab,ph_stab = xh_stab +euh_stab = uex-uh_stab +eph_stab = pex-ph_stab +edivuh_stab = divuex-(∇⋅uh_stab) +norm_euh_stab = sum(∫(euh_stab⋅euh_stab)*dΩ) +norm_eph_stab = sum(∫(eph_stab*eph_stab)*dΩ) +norm_edivuh_stab = sum(∫(edivuh_stab⋅edivuh_stab)*dΩ) +# k=0 (n=16) yields euh = 1.6e-26, eph = 0.00015, edivuh = 4.2e-29, κ = 7.6e6 +# k=2 (n=16) yields euh = 3.1e-14, eph = 4.1e-29, edivuh = 1.8e-17, κ = 4.7e16 \ No newline at end of file diff --git a/bulk_ghost_penalty_canvas_selecting_cut_cells.jl b/bulk_ghost_penalty_canvas_selecting_cut_cells.jl deleted file mode 100644 index 1d4d723..0000000 --- a/bulk_ghost_penalty_canvas_selecting_cut_cells.jl +++ /dev/null @@ -1,349 +0,0 @@ -using Gridap -using GridapEmbedded -using FillArrays -using LinearAlgebra - -include("BulkGhostPenaltyAssembleMaps.jl") - -# Manufactured solution -order = 0 -uex(x) = -VectorValue(2*x[1],2*x[2]) -divuex(x) = -4.0 - -# Select geometry -R = 0.2 -geom = disk(R, x0=Point(0.5,0.5)) - -# Setup background model -n=16 -partition = (n,n) -box = get_metadata(geom) -bgmodel = CartesianDiscreteModel((0,1,0,1),partition) -dp = box.pmax - box.pmin -h = dp[1]/n - -# Cut the background model with the mesh -cutgeo = cut(bgmodel,geom) - -# Compute mapping among background model -# cut cells and interior cells -strategy = AggregateAllCutCells() -aggregates = aggregate(strategy,cutgeo) - -#INVESTIGATING: Number of aggregates, cells and bg cells. -@assert length(aggregates) == n*n -num_cell = length(aggregates) -using StatsBase -counts = countmap(aggregates,alg=:dict) -agg_counts = filter(((k,v),) -> v > 1 && k != 0, counts) -num_aggregates = length(agg_counts) -num_cells_in_aggregates = 0 -for (k,v) in agg_counts - println(v) - num_cells_in_aggregates +=v -end -num_cells_not_in_aggregates = counts[0] -non_agg_counts = filter(((k,v),) -> v <= 1 || k==0, counts) -num_cells_not_in_aggregates = 0 -for (k,v) in non_agg_counts - println(v) - num_cells_not_in_aggregates +=v -end -num_ext_cells = counts[0] -@assert num_cells_in_aggregates + num_cells_not_in_aggregates == num_cell -println("Out of the $num_cell cells, $num_cells_in_aggregates belong to aggregates and $num_cells_not_in_aggregates are not part of aggregates. Note that we have $num_ext_cells exterior cells") -println("Assuming that there is one root per aggregate, we have $(num_cells_in_aggregates-num_aggregates) cut cells and $num_aggregates roots in the $num_aggregates aggregates cotntaining $num_cells_in_aggregates in total.") - -""" - Creates an array of arrays with as many entries - as aggregates. For each aggregate, the array - contains the global cell IDs of that cells in the background - model that belong to the same aggregate - - TO-DO: with efficiency in mind we may want to store this - array of arrays as a Gridap.Arrays.Table. -""" -function setup_aggregate_to_cells(aggregates) - - size_aggregates=Dict{Int,Int}() - for (i,agg) in enumerate(aggregates) - if agg>0 - if !haskey(size_aggregates,agg) - size_aggregates[agg]=1 - else - size_aggregates[agg]+=1 - end - end - end - - touched=Dict{Int,Int}() - aggregate_to_cells=Vector{Vector{Int}}() - current_aggregate=1 - for (i,agg) in enumerate(aggregates) - if agg>0 - if (size_aggregates[agg]>1) - if !haskey(touched,agg) - push!(aggregate_to_cells,[i]) - touched[agg]=current_aggregate - current_aggregate+=1 - else - push!(aggregate_to_cells[touched[agg]],i) - end - end - end - end - aggregate_to_cells -end - -""" - Creates an array of arrays with as many entries - as interior cells that are not part of any aggegrate. - For each interior cell, the array - contains the global cell IDs of that cells in the background - model that belong to the same interior cell - - TO-DO: with efficiency in mind we may want to store this - array of arrays as a Gridap.Arrays.Table. -""" -function setup_int_nonagg_cell_to_cells(aggregates) - - size_aggregates=Dict{Int,Int}() - for (i,agg) in enumerate(aggregates) - if agg>0 - if !haskey(size_aggregates,agg) - size_aggregates[agg]=1 - else - size_aggregates[agg]+=1 - end - end - end - - touched=Dict{Int,Int}() - interior_cell_to_cells=Vector{Vector{Int}}() - current_interior_cell=1 - for (i,agg) in enumerate(aggregates) - if agg>0 - if (size_aggregates[agg]==1) - if !haskey(touched,agg) - push!(interior_cell_to_cells,[i]) - touched[agg]=current_interior_cell - current_interior_cell+=1 - else - push!(interior_cell_to_cells[touched[agg]],i) - end - end - end - end - interior_cell_to_cells -end - - -function setup_aggregates_bounding_box_model(bgmodel, aggregate_to_cells) - g=get_grid(bgmodel) - cell_coords=get_cell_coordinates(g) - D=num_dims(bgmodel) - xmin=Vector{Float64}(undef,D) - xmax=Vector{Float64}(undef,D) - - # Compute coordinates of the nodes defining the bounding boxes - bounding_box_node_coords= - Vector{Point{D,Float64}}(undef,length(aggregate_to_cells)*2^D) - ptr = [ (((i-1)*2^D)+1) for i in 1:length(aggregate_to_cells)+1 ] - data = collect(1:length(bounding_box_node_coords)) - bounding_box_node_ids = Gridap.Arrays.Table(data,ptr) - for (agg,cells) in enumerate(aggregate_to_cells) - p=first(cell_coords[cells[1]]) - for i in 1:D - xmin[i]=p[i] - xmax[i]=p[i] - end - for cell in cells - for p in cell_coords[cell] - for i in 1:D - xmin[i]=min(xmin[i],p[i]) - xmax[i]=max(xmax[i],p[i]) - end - end - end - bounds = [(xmin[i], xmax[i]) for i in 1:D] - point_iterator = Iterators.product(bounds...) - bounding_box_node_coords[bounding_box_node_ids[agg]] = - reshape([Point(p...) for p in point_iterator],2^D) - end - - # Set up the discrete model of bounding boxes - HEX_AXIS=1 - polytope=Polytope(Fill(HEX_AXIS,D)...) - scalar_reffe=ReferenceFE(polytope,lagrangian,Float64,1) - cell_types=fill(1,length(bounding_box_node_ids)) - cell_reffes=[scalar_reffe] - grid = Gridap.Geometry.UnstructuredGrid(bounding_box_node_coords, - bounding_box_node_ids, - cell_reffes, - cell_types, - Gridap.Geometry.Oriented()) - Gridap.Geometry.UnstructuredDiscreteModel(grid) -end - -int_nonagg_cell_to_cells = setup_int_nonagg_cell_to_cells(aggregates) -aggregate_to_cells=setup_aggregate_to_cells(aggregates) -aggregates_bounding_box_model= - setup_aggregates_bounding_box_model(bgmodel,aggregate_to_cells) - -# colors = color_aggregates(aggregates,bgmodel) -# writevtk(Triangulation(bgmodel),"trian",celldata=["cellin"=>aggregates,"color"=>colors]) -# writevtk(aggregates_bounding_box_model, "bb_model") -# writevtk(bgmodel, "bg_model") - -""" - Changes the domain of a trial/test basis defined on - the reference space of bounding boxes to the reference - space of the agg cells - - TO-DO: in the future, for system of PDEs (MultiField) we should - also take care of blocks (BlockMap) -""" -function change_domain_bb_to_agg_cells(basis_bb, - ref_agg_cell_to_ref_bb_map, - Ω_agg_cells, - agg_cells_to_aggregate) - @assert num_cells(Ω_agg_cells)==length(ref_agg_cell_to_ref_bb_map) - @assert Gridap.CellData.DomainStyle(basis_bb)==ReferenceDomain() - bb_basis_style = Gridap.FESpaces.BasisStyle(basis_bb) - bb_basis_array = Gridap.CellData.get_data(basis_bb) - if (bb_basis_style==Gridap.FESpaces.TrialBasis()) - # Remove transpose map; we will add it later - @assert isa(bb_basis_array,Gridap.Arrays.LazyArray) - @assert isa(bb_basis_array.maps,Fill) - @assert isa(bb_basis_array.maps.value,typeof(transpose)) - bb_basis_array=bb_basis_array.args[1] - end - - bb_basis_array_to_Ω_agg_cells_array = lazy_map(Reindex(bb_basis_array),agg_cells_to_aggregate) - bb_basis_array_to_Ω_agg_cells_array = lazy_map(Broadcasting(∘), - bb_basis_array_to_Ω_agg_cells_array, - ref_agg_cell_to_ref_bb_map) - if (bb_basis_style==Gridap.FESpaces.TrialBasis()) - # Add transpose - bb_basis_array_to_Ω_agg_cells_array=lazy_map(transpose, bb_basis_array_to_Ω_agg_cells_array) - end - - Gridap.CellData.GenericCellField(bb_basis_array_to_Ω_agg_cells_array, - Ω_agg_cells, - ReferenceDomain()) -end - -# Generate an array with the global IDs of the cells -# that belong to an aggregrate. From now on, we will -# use the terminology "agg_cells" to refer to those -# cells of the background model that belong to an aggregate -# (i.e., they can be either cut or interior cells) -agg_cells=Vector{Int}() -for cells in aggregate_to_cells - append!(agg_cells,cells) -end - -# Generate an array with the global IDs of the cells -# that belong to the interior, yet do not belong to the -# aggregate. We will use "int_nonagg_cells" to refer to those -# cells of the background model that belong to the interior -# but are not part of any of the aggregates. -int_nonagg_cells=Vector{Int}() -for cell in int_nonagg_cell_to_cells - append!(int_nonagg_cells,cell) -end - -#INSPECTION -# agg_cells -# int_nonagg_cells -# aggregates - -## Triangulation -Ωbg = Triangulation(bgmodel) -Ω = Triangulation(cutgeo,PHYSICAL) -Ωact = Triangulation(cutgeo,ACTIVE) -Ωcut = Triangulation(cutgeo) - -# Interior cells -Ωbg_agg_cells = view(Ωbg,agg_cells) # cells in aggregates -int_cells = Ω_agg_cells.parent.b.tface_to_mface -# construct the interior aggregate cells, in other words the root cells -root_cells=Vector{Int}() -tester_int_nonagg_cells=Vector{Int}() -for cell in int_cells - if cell ∈ int_nonagg_cells - append!(tester_int_nonagg_cells,cell) - else - append!(root_cells,cell) - end -end -@assert(tester_int_nonagg_cells==int_nonagg_cells) -# Cut cells -cut_cells=Vector{Int}() -tester_root_cells=Vector{Int}() -for cell in agg_cells - if cell ∈ root_cells - append!(tester_root_cells,cell) - else - append!(cut_cells,cell) - end -end -@assert(tester_root_cells==root_cells) - - -# TODO Potentially useful?: -# Ω_agg_cells.parent.a.subgrid.cell_node_ids.data == Ω.a.subgrid.cell_node_ids.data # vector containing cell ids - -# Subdomains in background mesh (aggegrate cells are already defined) -Ωbg_int_cells = view(Ωbg,int_cells) # cells in interior -Ωbg_int_nonagg_cells = view(Ωbg,int_nonagg_cells) # cells in interior but not in aggregate -Ωbg_root_cells = view(Ωbg,root_cells) # root cells: in interior and part of aggregate -Ωbg_cut_cells = view(Ωbg,cut_cells) # cut cells (part of aggregate) - -# Subdomains in physical mesh -Ω_agg_cells = view(Ω,agg_cells) # cells in aggregates -Ω_int_cells = view(Ω,int_cells) # cells in interior -Ω_int_nonagg_cells = view(Ω,int_nonagg_cells) # cells in interior but not in aggregate -Ω_root_cells = view(Ω,root_cells) # root cells: in interior and part of aggregate -Ω_cut_cells = view(Ω,cut_cells) # cut cells (part of aggregate) - -## VISUALISATION: -writevtk(Ωbg,"trian_bg") -writevtk(Ωbg_agg_cells,"trian_bg_agg_cells") -writevtk(Ωbg_int_cells,"trian_bg_int_cells") -writevtk(Ωbg_int_nonagg_cells,"trian_bg_int_nonagg_cells") -writevtk(Ωbg_root_cells,"trian_bg_root_cells") -writevtk(Ωbg_cut_cells,"trian_bg_cut_cells") - -writevtk(Ωact,"trian_act") -writevtk(Ω,"trian_phys") -#TODO: WHY DOESN'T THIS WORK? -# writevtk(Ω_agg_cells,"trian_phys_agg_cells") -# writevtk(Ω_int_cells,"trian_phys_int_cells") -# writevtk(Ω_int_nonagg_cells,"trian_phys_int_nonagg_cells") -# writevtk(Ω_root_cells,"trian_phys_root_cells") -# writevtk(Ω_cut_cells,"trian_phys_cut_coot_cells") -degree = 2*2*order -dΩ = Measure(Ω,degree) -dΩact = Measure(Ωact,degree) -dΩbg_agg_cells = Measure(Ωbg_agg_cells,degree) -dΩbg_cut_cells = Measure(Ωbg_cut_cells,degree) -dΩbg_root_cells = Measure(Ωbg_root_cells,degree) - - -dΩ_agg_cells = Measure(Ω_agg_cells,degree) -dΩcut_cells = Measure(Ω_cut_cells,degree) -dΩroot_cells = Measure(Ω_root_cells,degree) - -# TESTING -area_phys_domain = sum(∫(1.0)dΩ) -area_act_domain = sum(∫(1.0)dΩact) -area_agg_cells = sum(∫(1.0)dΩbg_agg_cells) -area_cut_cells = sum(∫(1.0)dΩbg_cut_cells) -area_root_cells = sum(∫(1.0)dΩbg_root_cells) -@assert(area_agg_cells ≈ area_cut_cells + area_root_cells) - -# TODO: not possible to integrate over the physical part of the domain, e.g. -area_phys_agg_cells = sum(∫(1.0)dΩ_agg_cells) -area_phys_cut_cells = sum(∫(1.0)dΩ_cut_cells) -area_phys_root_cells = sum(∫(1.0)dΩ_root_cells) \ No newline at end of file diff --git a/bulk_ghost_penalty_stab_tools.jl b/bulk_ghost_penalty_stab_tools.jl index c56248e..a4e0e65 100644 --- a/bulk_ghost_penalty_stab_tools.jl +++ b/bulk_ghost_penalty_stab_tools.jl @@ -100,6 +100,108 @@ function bulk_ghost_penalty_stabilization_collect_cell_matrix(agg_cells_to_aggre w, r, c end +## Create collect function on cut cells only +""" + dv, du: Test and trial basis functions. They may be components of a MultiFieldCellField + + # Compute and assemble the bulk penalty stabilization term + # ∫( (dv-dv_l2_proj_agg_cells)*(du-du_l2_proj_agg_cells))*dΩ_cut_cells (long version) + # ∫( (dv)*(du-du_l2_proj_agg_cells))*dΩ_cut_cells (simplified, equivalent version) +""" +function bulk_ghost_penalty_stabilization_collect_cell_matrix_on_cut_cells(agg_cells_to_aggregate, + ref_agg_cell_to_ref_bb_map, + dΩbg_agg_cells, + dv, # Test basis + du, # Trial basis (to project) + dvbb, # Bounding box space test basis + lhs, + Ωbg_cut_cell_dof_ids, + agg_cells_local_dof_ids, + cut_cells_to_aggregate_dof_ids, + γ, + dΩbg_cut_cells) + + + Ωbg_agg_cells=dΩbg_agg_cells.quad.trian + + # Change domain of vbb (test) from Ωbb to Ωagg_cells + dvbb_Ωbg_agg_cells=change_domain_bb_to_agg_cells(dvbb, + ref_agg_cell_to_ref_bb_map, + Ωbg_agg_cells, + agg_cells_to_aggregate) + + du_single_field=_get_single_field_fe_basis(du) + agg_cells_rhs_contribs=get_array(∫(dvbb_Ωbg_agg_cells⋅du_single_field)dΩbg_agg_cells) + ass_rhs_map=BulkGhostPenaltyAssembleRhsMap(agg_cells_local_dof_ids,agg_cells_rhs_contribs) + rhs=lazy_map(ass_rhs_map,aggregate_to_local_cells) + + # TO-DO: optimize using our own optimized version Gridap.Fields.Map + # of backslash that re-uses storage for lu factors among cells, etc. + dv_l2_proj_bb_dofs=lazy_map(\,lhs,rhs) + + # Generate bb-wise array of fields. For each aggregate's bounding box, + # it provides the l2 projection of all basis functions in du + # restricted to the cells included in the bounding box of the aggregate + dv_l2_proj_bb_array=lazy_map(Gridap.Fields.linear_combination, + dv_l2_proj_bb_dofs, + Gridap.CellData.get_data(dvbb)) + + # # Change domain of dv_l2_proj_bb_array from bb to agg_cells + dv_l2_proj_bb_array_agg_cells=lazy_map(Broadcasting(∘), + lazy_map(Reindex(dv_l2_proj_bb_array),agg_cells_to_aggregate), + ref_agg_cell_to_ref_bb_map) + + du_l2_proj_bb_array_agg_cells=lazy_map(transpose, dv_l2_proj_bb_array_agg_cells) + + if (_is_multifield_fe_basis_component(du)) + @assert _is_multifield_fe_basis_component(dv) + @assert _nfields(du)==_nfields(dv) + nfields=_nfields(du) + fieldid=_fieldid(du) + du_l2_proj_bb_array_agg_cells=lazy_map( + Gridap.Fields.BlockMap((1,nfields),fieldid), + du_l2_proj_bb_array_agg_cells) + end + + du_l2_proj_agg_cells = Gridap.CellData.GenericCellField(du_l2_proj_bb_array_agg_cells, Ωbg_agg_cells,ReferenceDomain()) + + # Manually set up the arrays that collect_cell_matrix would return automatically + w = [] + r = [] + c = [] + + dv_du_mat_contribs=get_array(∫(γ*dv⋅du)*dΩbg_cut_cells) + + if (_is_multifield_fe_basis_component(dv)) + nfields=_nfields(dv) + fieldid=_fieldid(dv) + Ωbg_cut_cell_dof_ids=lazy_map(Gridap.Fields.BlockMap(nfields,fieldid),Ωbg_cut_cell_dof_ids) + end + push!(w, dv_du_mat_contribs) + push!(r, Ωbg_cut_cell_dof_ids) + push!(c, Ωbg_cut_cell_dof_ids) + + # In the MultiField case, I have had to add this change domain + # call before setting up the term right below. Otherwise, we get an error + # when trying to multiply the fields. Not sure why this is happening + dvΩagg_cells=Gridap.CellData.change_domain(dv,Ωbg_agg_cells,ReferenceDomain()) + + proj_dv_du_mat_contribs=get_array(∫(γ*(-1.0)*dvΩagg_cells⋅(du_l2_proj_agg_cells))*dΩbg_cut_cells) + + if (_is_multifield_fe_basis_component(du)) + nfields=_nfields(du) + fieldid=_fieldid(du) + cut_cells_to_aggregate_dof_ids= + lazy_map(Gridap.Fields.BlockMap(nfields,fieldid),cut_cells_to_aggregate_dof_ids) + end + + push!(w, proj_dv_du_mat_contribs) + push!(r, Ωbg_cut_cell_dof_ids) + push!(c, cut_cells_to_aggregate_dof_ids) + + w, r, c +end + function div_penalty_stabilization_collect_cell_matrix(agg_cells_to_aggregate, ref_agg_cell_to_ref_bb_map, dΩagg_cells, @@ -214,6 +316,89 @@ function div_penalty_stabilization_collect_cell_matrix(agg_cells_to_aggregate, w,r,c end + +function div_penalty_stabilization_collect_cell_matrix_on_cut_cells(agg_cells_to_aggregate, + ref_agg_cell_to_ref_bb_map, + dΩbg_agg_cells, + dv, # Test basis + du, # Trial basis (to project) + qbb, # Bounding box space test basis + p_lhs, + U_Ωbg_cut_cell_dof_ids, + U_agg_cells_local_dof_ids, + U_cut_cells_to_aggregate_dof_ids, + γ, + dΩbg_cut_cells) + + Ωbg_agg_cells=dΩbg_agg_cells.quad.trian + + ## Compute ∫( (div_dv)*(div_du-div_du_l2_proj_agg_cells))*dΩbg_cut_cells + # Manually set up the arrays that collect_cell_matrix would return automatically + w = [] + r = [] + c = [] + + if (_is_multifield_fe_basis_component(dv)) + nfields=_nfields(dv) + fieldid=_fieldid(dv) + U_Ωbg_cut_cell_dof_ids=lazy_map(Gridap.Fields.BlockMap(nfields,fieldid),U_Ωbg_cut_cell_dof_ids) + end + + ## (I) Let's set up the div_dv_div_du_mat_contribs + # (first part of the integral, that is ∫( (div_dv)*(div_du)*dΩbg_cut_cells + div_dv_div_du_mat_contribs=get_array(∫(γ*(∇⋅dv)⋅(∇⋅du))*dΩbg_cut_cells) + push!(w, div_dv_div_du_mat_contribs) + push!(r, U_Ωbg_cut_cell_dof_ids) + push!(c, U_Ωbg_cut_cell_dof_ids) + + ## (II) Set up ∫( (div_dv)*(-div_du_l2_proj_agg_cells))*dΩbg_cut_cells + div_du=∇⋅(_get_single_field_fe_basis(du)) + + # Compute Π_Q_bb(div_du) + dqbb_Ωagg_cells =change_domain_bb_to_agg_cells(qbb,ref_agg_cell_to_ref_bb_map,Ωbg_agg_cells,agg_cells_to_aggregate) + agg_cells_rhs_contribs=get_array(∫(dqbb_Ωagg_cells⋅div_du)dΩbg_agg_cells) + ass_rhs_map=BulkGhostPenaltyAssembleRhsMap(U_agg_cells_local_dof_ids,agg_cells_rhs_contribs) + rhs=lazy_map(ass_rhs_map,aggregate_to_local_cells) + div_dv_l2_proj_bb_dofs=lazy_map(\,p_lhs,rhs) + + # Generate bb-wise array of fields. For each aggregate's bounding box, + # it provides the l2 projection of all basis functions in dp + # restricted to the cells included in the bounding box of the aggregate + div_dv_l2_proj_bb_array=lazy_map(Gridap.Fields.linear_combination, div_dv_l2_proj_bb_dofs, Gridap.CellData.get_data(qbb)) + + div_dv_l2_proj_bb_array_agg_cells=lazy_map(Broadcasting(∘),lazy_map(Reindex(div_dv_l2_proj_bb_array),agg_cells_to_aggregate),ref_agg_cell_to_ref_bb_map) + div_du_l2_proj_bb_array_agg_cells=lazy_map(transpose, div_dv_l2_proj_bb_array_agg_cells) + + if (_is_multifield_fe_basis_component(du)) + @assert _is_multifield_fe_basis_component(du) + @assert _nfields(du)==_nfields(dv) + nfields=_nfields(du) + fieldid=_fieldid(dv) + div_du_l2_proj_bb_array_agg_cells=lazy_map(Gridap.Fields.BlockMap((fieldid,nfields),fieldid),div_du_l2_proj_bb_array_agg_cells) + end + + div_du_l2_proj_agg_cells = Gridap.CellData.GenericCellField(div_du_l2_proj_bb_array_agg_cells,Ωbg_agg_cells,ReferenceDomain()) + + # In the MultiField case, I have had to add this change domain + # call before setting up the term right below. Otherwise, we get an error + # when trying to multiply the fields. Not sure why this is happening + # dv_Ωagg_cells=Gridap.CellData.change_domain(dv,Ωagg_cells,ReferenceDomain()) + div_dv_Ωagg_cells=Gridap.CellData.change_domain(∇⋅dv,Ωbg_agg_cells,ReferenceDomain()) + + proj_div_dv_div_du_mat_contribs=get_array(∫(γ*(-1.0)*div_dv_Ωagg_cells*(div_du_l2_proj_agg_cells))*dΩbg_cut_cells) + + if (_is_multifield_fe_basis_component(du)) + nfields=_nfields(du) + fieldid=_fieldid(du) + U_cut_cells_to_aggregate_dof_ids=lazy_map(Gridap.Fields.BlockMap(nfields,fieldid),U_cut_cells_to_aggregate_dof_ids) + end + + push!(w, proj_div_dv_div_du_mat_contribs) + push!(r, U_Ωbg_cut_cell_dof_ids) # test + push!(c, U_cut_cells_to_aggregate_dof_ids) # trial + w,r,c +end + function set_up_bulk_ghost_penalty_lhs(aggregate_to_local_cells, agg_cells_to_aggregate, ref_agg_cell_to_ref_bb_map,