From 9f1b91ebf9e20a1a10304dd8fb8e8c728bdb4950 Mon Sep 17 00:00:00 2001 From: Anne Boschman Date: Fri, 8 Nov 2024 17:43:30 +1100 Subject: [PATCH] selecting cut cells --- ...host_penalty_canvas_selecting_cut_cells.jl | 349 ++++++++++++++++++ 1 file changed, 349 insertions(+) create mode 100644 bulk_ghost_penalty_canvas_selecting_cut_cells.jl diff --git a/bulk_ghost_penalty_canvas_selecting_cut_cells.jl b/bulk_ghost_penalty_canvas_selecting_cut_cells.jl new file mode 100644 index 0000000..1d4d723 --- /dev/null +++ b/bulk_ghost_penalty_canvas_selecting_cut_cells.jl @@ -0,0 +1,349 @@ +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