From 34f7d168089f9bf1fa4f1add4e58e6a4f3f57dc2 Mon Sep 17 00:00:00 2001 From: Anne Boschman Date: Fri, 30 Aug 2024 11:15:53 +1000 Subject: [PATCH] Generalization to 3D --- bulk_ghost_penalty_canvas.jl | 34 +++++++++++++++------------------- 1 file changed, 15 insertions(+), 19 deletions(-) diff --git a/bulk_ghost_penalty_canvas.jl b/bulk_ghost_penalty_canvas.jl index 2dbcda7..3015e5b 100644 --- a/bulk_ghost_penalty_canvas.jl +++ b/bulk_ghost_penalty_canvas.jl @@ -77,8 +77,10 @@ function setup_aggregates_bounding_box_model(bgmodel, aggregate_to_cells) # Compute coordinates of the nodes defining the bounding boxes bounding_box_node_coords= - Vector{Point{D,Float64}}(undef,length(aggregate_to_cells)*2^D) - current=1 + 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 @@ -93,24 +95,18 @@ function setup_aggregates_bounding_box_model(bgmodel, aggregate_to_cells) end end end - # TO-DO: this nested loop is only valid for D=2. - # Ideally, we should be able to generalize it for - # arbitrary D - for i in (xmin[2],xmax[2]) - for j in (xmin[1],xmax[1]) - bounding_box_node_coords[current]=Point(j,i) - current+=1 - end - end + bounds = [(xmin[i], xmax[i]) for i in 1:D] + point_iterator = Iterators.product([bounds[i] for i in 1:D]...) + bounding_box_node_coords[bounding_box_node_ids[agg]] = + reduce(vcat,[Point(p...) for p in point_iterator]) end - # Set up the discrete model of bounding boxes - 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) - - # TO-DO: this polytope is only valid for D=2 - polytope=QUAD + # Set up the discrete model of bounding boxes + if D==2 + polytope=QUAD + elseif D==3 + polytope=HEX + end scalar_reffe=ReferenceFE(polytope,lagrangian,Float64,1) cell_types=fill(1,length(bounding_box_node_ids)) cell_reffes=[scalar_reffe] @@ -123,7 +119,7 @@ function setup_aggregates_bounding_box_model(bgmodel, aggregate_to_cells) end aggregate_to_cells=setup_aggregate_to_cells(aggregates) -aggregates_bounding_box_model= +aggregates_bounding_box_model = setup_aggregates_bounding_box_model(bgmodel,aggregate_to_cells) colors = color_aggregates(aggregates,bgmodel)