Skip to content

Commit

Permalink
revised and merged functions that allow for selection of (aggregated)…
Browse files Browse the repository at this point in the history 
… cells, e.g. cut,interior, exterior.
  • Loading branch information
@amboschman
amboschman committed Jan 30, 2025
1 parent 6d65126 commit 7cb9098
Show file tree
Hide file tree
Showing 2 changed files with 60 additions and 253 deletions.
255 changes: 38 additions & 217 deletions aggregates_bounding_boxes_tools.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -39,136 +39,47 @@ function setup_aggregate_to_cells(aggregates)
end

"""
Creates an array of arrays with as many entries
as aggregates. For each aggregate, the array
contains the global cell IDs of that cut cells in
the background model that belong to the same aggregate
Creates an array containing the cell IDs of a selection of the background cells. The cell type filter IN (-1) selects the interior cells, OUT (1) the exterior cells, and GridapEmbedded.Interfaces.CUT (0) the cut cells.
TO-DO: publish CUT, so that GridapEmbedded.Interfaces.CUT can be shortened to CUT?
TO-DO: be careful with using restrict_cells(cutgeo,GridapEmbedded.Interfaces.CUT) to replace flatten(restrict_aggregate_to_cells(cutgeo,aggregate_to_cells,GridapEmbedded.Interfaces.CUT))
TO-DO: with efficiency in mind we may want to store this
array of arrays as a Gridap.Arrays.Table.
"""
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
function restrict_cells(cutgeo,cell_type_filter)
restricted_cells=Int64[]
cell_to_inoutcut=compute_bgcell_to_inoutcut(cutgeo,cutgeo.geo)
for cell in 1:length(cell_to_inoutcut)
if cell_to_inoutcut[cell] == cell_type_filter
push!(restricted_cells, cell)
end
aggregate_to_cut_cells
end
end
restricted_cells
end

"""
Creates an array of arrays with as many entries
as aggregates. For each aggregate, the array
contains the global cell IDs of that cut 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 revised_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)
# println("i=$i, agg=$agg")
if agg>0
if (size_aggregates[agg]>1)
if !haskey(touched,agg)
if i ∉ root_cells
push!(aggregate_to_cut_cells,[i])
else
push!(aggregate_to_cut_cells,[])
end
touched[agg]=current_aggregate
current_aggregate+=1
else
if i ∉ root_cells
push!(aggregate_to_cut_cells[touched[agg]],i)
end
end
end
end
# println("touched = $touched")
# println("current agg = $current_aggregate, $aggregate_to_cut_cells")
end
aggregate_to_cut_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.
contains a selection of the global cell IDs of
the cells in the background model that belong to
the same aggregate. The cell type filter IN (-1)
selects the interior cells, OUT (1) the exterior
cells, and GridapEmbedded.Interfaces.CUT (0) the
cut cells.
TO-DO: publish CUT, so that GridapEmbedded.Interfaces.CUT can be shortened to CUT?
"""
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
function restrict_aggregate_to_cells(cutgeo,aggregate_to_cells,cell_type_filter)
aggregate_to_restricted_cells=Vector{Int}[]
cell_to_inoutcut=compute_bgcell_to_inoutcut(cutgeo,cutgeo.geo)
for agg in aggregate_to_cells
restricted_cells = Int64[]
for cell in agg
if cell_to_inoutcut[cell] == cell_type_filter
push!(restricted_cells, cell)
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
push!(aggregate_to_restricted_cells,restricted_cells)
end
interior_cell_to_cells
end
aggregate_to_restricted_cells
end

function setup_aggregates_bounding_box_model(bgmodel, aggregate_to_cells)
g=get_grid(bgmodel)
Expand Down Expand Up @@ -218,83 +129,11 @@ function setup_aggregates_bounding_box_model(bgmodel, aggregate_to_cells)
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)
TO-DO: perhaps merge this function with setup_int_nonagg_cells
"""
function setup_agg_cells(aggregate_to_cells)
agg_cells=Vector{Int}()
for cells in aggregate_to_cells
append!(agg_cells,cells)
end
return agg_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. Thus, all interior
cells but not including the root cells of the aggregates.
TO-DO: perhaps merge this function with setup_agg_cells
"""
function setup_int_nonagg_cells(int_nonagg_cell_to_cells)
int_nonagg_cells=Vector{Int}()
for cell in int_nonagg_cell_to_cells
append!(int_nonagg_cells,cell)
end
return int_nonagg_cells
end

"""
Generate an array with the global IDs of the cells
that are the root cells of the aggregrates.
TO-DO: perhaps merge this function with setup_cut_cells,
e.g. as a function that selects cells from list A which
are not part of a list B.
"""
function setup_root_cells(int_cells, int_nonagg_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)
return root_cells
end

"""
Generate an array with the global IDs of the cells
that are the cut cells of the aggregrates.
TO-DO: perhaps merge this function with setup_cut_cells,
e.g. as a function that selects cells from list A which
are not part of a list B.
Flattens an array of arrays
"""
function setup_cut_cells(agg_cells, root_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)
cut_cells
end
function flatten(array_of_arrays)
vcat(array_of_arrays...)
end

"""
Generate an array that given the local ID of an "agg_cell"
Expand All @@ -306,34 +145,16 @@ end
setup_cut_cells_in_agg_cells_to_aggregate
"""
function setup_agg_cells_to_aggregate(aggregate_to_cells)
agg_cells_to_aggregate=Vector{Int}()
function setup_cells_to_aggregate(aggregate_to_cells)
cells_to_aggregate=Vector{Int}()
for (i,cells) in enumerate(aggregate_to_cells)
for _ in cells
push!(agg_cells_to_aggregate,i)
push!(cells_to_aggregate,i)
end
end
agg_cells_to_aggregate
cells_to_aggregate
end

"""
Generate an array that given the local ID of an cut "agg_cell"
returns the ID of the aggregate to which it belongs
(i.e., flattened version of aggregate_to_cut_cells)
TO-DO: perhaps merge this function with ,
setup_agg_cells_to_aggregate
"""
function setup_cut_cells_in_agg_cells_to_aggregate(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

"""
Renumbers the global cell IDs to local cell IDs in the array of arrays
with as many entries as aggregates. For each aggregate, the array
Expand Down
58 changes: 22 additions & 36 deletions bulk_ghost_penalty_canvas.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -57,13 +57,11 @@ function setup_geometry(nint, ε, pmin, pmax)

bgmodel, cutgeo, h
end
bgmodel, cutgeo, h = setup_geometry(nint, ε, pmin, pmax)

# Compute mapping among background model
# cut cells and interior cells
strategy = AggregateAllCutCells()
aggregates = aggregate(strategy,cutgeo)
bgmodel, cutgeo, h= setup_geometry(nint, ε, pmin, pmax)

# Setup aggregates
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)
Expand All @@ -77,20 +75,19 @@ aggregates_bounding_box_model=
degree=2*2*(order+1)
dΩ = Measure(Ω,degree)

## (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)
@assert length(int_nonagg_cells)>0 # otherwise we can not constrain a dof
# Setup agg cells and triangulation
agg_cells =flatten(aggregate_to_cells) #[CLEAN]: replaced setup_agg_cells
Ωbg_agg_cells=view(Ωbg,agg_cells)

# Pressure space
if problem==0
Q = FESpace(Ωact, ReferenceFE(lagrangian,Float64,order), conformity=:L2)
elseif problem==1
# Set up zero-mean pressure space, with fixed interior dof
Qnzm = FESpace(Ωact, ReferenceFE(lagrangian,Float64,order), conformity=:L2)
@assert nint>2
int_nonagg_dof_to_fix = int_nonagg_cells[1] # pick the first interior cell that is not part of an aggregate to be constrained by the zero mean pressure condition.
spaceWithConstantFixed = Gridap.FESpaces.FESpaceWithConstantFixed(Qnzm,true,int_nonagg_dof_to_fix)
int_cells = restrict_cells(cutgeo,IN) #TODO: this dof may belong to an aggregate (root cell)
nonagg_int_cell = int_cells[findfirst(!in(agg_cells),int_cells)]
spaceWithConstantFixed = Gridap.FESpaces.FESpaceWithConstantFixed(Qnzm,true,nonagg_int_cell)
Qzm_vol_i = assemble_vector(v->∫(v)*dΩ,Qnzm)
Qzm_vol = sum(Qzm_vol_i)
Q = Gridap.FESpaces.ZeroMeanFESpace(spaceWithConstantFixed,Qzm_vol_i,Qzm_vol)
Expand All @@ -107,13 +104,9 @@ 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)
agg_cells_to_aggregate =setup_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)

Expand Down Expand Up @@ -173,31 +166,24 @@ h_U = 1.0
### STABILIZATION ON Ωagg\Troot ###
###########################################

## (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)
# Setup cut cells and triangulation
aggregate_to_cut_cells = restrict_aggregate_to_cells(cutgeo,aggregate_to_cells,GridapEmbedded.Interfaces.CUT)
cut_cells = flatten(aggregate_to_cut_cells)
#TO-DO: look into why cut_cells = restrict_cells(cutgeo,GridapEmbedded.Interfaces.CUT) can not be used instead of the above
Ωbg_cut_cells = view(Ωbg,cut_cells)
dΩbg_cut_cells = Measure(Ωbg_cut_cells,degree)

# (5) DOF IDS related to cut cells
# Selecting relevant global dofs ids of cut cells (from background mesh)
Ω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)
aggregate_to_cut_cells = revised_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)
# Compute global dofs ids per aggregate and reindex these
cut_cells_to_aggregate = setup_cells_to_aggregate(aggregate_to_cut_cells)
U_cut_cells_to_aggregate_dof_ids=lazy_map(Reindex(U_aggregate_dof_ids),cut_cells_to_aggregate)
P_cut_cells_to_aggregate_dof_ids=lazy_map(Reindex(P_aggregate_dof_ids),cut_cells_to_aggregate)

# (7) Compute stabilization terms for u and p
# 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,
Expand Down

0 comments on commit 7cb9098

Please sign in to comment.