From 2efc08a8f09f3d83bf7eae1f09eb0bd5c8ac302c Mon Sep 17 00:00:00 2001 From: Hridya Dilip Date: Mon, 31 Aug 2020 19:33:00 +1000 Subject: [PATCH 1/8] Adding space_aggregation --- src/AgFEM/AgFEM.jl | 4 + src/AgFEM/CellAggregation_space.jl | 273 +++++++++++++++++++++++++++++ src/Exports.jl | 3 + 3 files changed, 280 insertions(+) create mode 100644 src/AgFEM/CellAggregation_space.jl diff --git a/src/AgFEM/AgFEM.jl b/src/AgFEM/AgFEM.jl index 323049a7..d7be6e2f 100644 --- a/src/AgFEM/AgFEM.jl +++ b/src/AgFEM/AgFEM.jl @@ -15,12 +15,16 @@ using GridapEmbedded.Interfaces using GridapEmbedded.Interfaces: compute_bgcell_to_inoutcut export aggregate +export aggregatespace export color_aggregates +export color_aggregates_space export AggregateAllCutCells +export AggregateSpaceCutCells export AgFEMSpace include("CellAggregation.jl") include("AgFEMSpaces.jl") +include("CellAggregation_space.jl") end # module diff --git a/src/AgFEM/CellAggregation_space.jl b/src/AgFEM/CellAggregation_space.jl new file mode 100644 index 00000000..9677fa4b --- /dev/null +++ b/src/AgFEM/CellAggregation_space.jl @@ -0,0 +1,273 @@ + +function aggregatespace(strategy,cut::EmbeddedDiscretization) + aggregatespace(strategy,cut,cut.geo) +end + +function aggregatespace(strategy,cut::EmbeddedDiscretization,geo) + aggregatespace(strategy,cut,geo,IN) +end + +function aggregatespace(strategy,cut::EmbeddedDiscretization,name::String,in_or_out) + geo = get_geometry(cut.geo,name) + aggregatespace(strategy,cut,geo,in_or_out) +end + +function aggregatespace(strategy,cut::EmbeddedDiscretization,geo::CSG.Geometry,in_or_out) + cell_to_inoutcut = compute_bgcell_to_inoutcut(cut,geo) + facet_to_inoutcut = compute_bgfacet_to_inoutcut(cut.bgmodel,geo) + aggregatespace(strategy,cut.bgmodel,in_or_out,cell_to_inoutcut,facet_to_inoutcut) +end + +function aggregatespace( + strategy, + cut::EmbeddedDiscretization, + cut_facets::EmbeddedFacetDiscretization, + name::String, + in_or_out) + + geo = get_geometry(cut.geo,name) + aggregatespace(strategy,cut,cut_facets,geo,in_or_out) +end + +function aggregatespace( + strategy, + cut::EmbeddedDiscretization, + cut_facets::EmbeddedFacetDiscretization) + + aggregatespace(strategy,cut,cut_facets,cut.geo,IN) +end + +function aggregatespace( + strategy, + cut::EmbeddedDiscretization, + cut_facets::EmbeddedFacetDiscretization, + geo::CSG.Geometry) + + aggregatespace(strategy,cut,cut_facets,geo,IN) +end + +function aggregatespace( + strategy, + cut::EmbeddedDiscretization, + cut_facets::EmbeddedFacetDiscretization, + geo::CSG.Geometry, + in_or_out) + + cell_to_inoutcut = compute_bgcell_to_inoutcut(cut,geo) + facet_to_inoutcut = compute_bgfacet_to_inoutcut(cut_facets,geo) + aggregatespace(strategy,cut.bgmodel,in_or_out,cell_to_inoutcut,facet_to_inoutcut) +end + +function aggregatespace(strategy,model::DiscreteModel,in_or_out,cell_to_inoutcut,facet_to_inoutcut) + @abstractmethod +end + +struct AggregateSpaceCutCells end + +function aggregatespace( + strategy::AggregateSpaceCutCells,model::DiscreteModel,in_or_out,cell_to_inoutcut,facet_to_inoutcut) + _compute_aggregatesspace(cell_to_inoutcut,facet_to_inoutcut,model,in_or_out) +end + +function _compute_aggregatesspace(cell_to_inoutcut,facet_to_inoutcut,model,loc) + @assert loc in (IN,OUT) + trian = get_triangulation(model) + cell_to_coords = get_cell_coordinates(trian) + topo = get_grid_topology(model) + D = num_cell_dims(model) + n_cell = num_cells(topo) + filter = [false,false,true,true] + filters = fill(filter,n_cell) + ptrs = collect(1:n_cell) + + cell_to_faces = get_faces(topo,D,D-1) + face_to_cells = get_faces(topo,D-1,D) + + comp_array = CompressedArray(cell_to_faces,ptrs) + filt_array = FilteredCellArray(comp_array,filters) + cell_to_faces = filt_array + + + + _compute_aggregates_barrier_space( + cell_to_inoutcut,facet_to_inoutcut,loc,cell_to_coords,cell_to_faces,face_to_cells) +end + +function _compute_aggregates_barrier_space( + cell_to_inoutcut,facet_to_inoutcut,loc,cell_to_coords,cell_to_faces,face_to_cells) + + n_cells = length(cell_to_inoutcut) + cell_to_cellin = zeros(Int32,n_cells) + cell_to_touched = fill(false,n_cells) + + for (cell, inoutcut) in enumerate(cell_to_inoutcut) + if inoutcut == loc + cell_to_cellin[cell] = cell + cell_to_touched[cell] = true + end + end + + c1 = array_cache(cell_to_faces) + c2 = array_cache(face_to_cells) + c3 = array_cache(cell_to_coords) + c4 = array_cache(cell_to_coords) + + max_iters = 20 + + all_aggregated = false + for iter in 1:max_iters + all_aggregated = true + for (cell, inoutcut) in enumerate(cell_to_inoutcut) + if inoutcut == CUT && ! cell_to_touched[cell] + neigh_cell = _find_best_neighbor_space( + c1,c2,c3,c4,cell, + cell_to_faces, + face_to_cells, + cell_to_coords, + cell_to_touched, + cell_to_cellin, + facet_to_inoutcut, + loc) + if neigh_cell > 0 + cellin = cell_to_cellin[neigh_cell] + cell_to_cellin[cell] = cellin + else + all_aggregated = false + end + end + end + if all_aggregated + break + end + _touch_aggregated_cells_space!(cell_to_touched,cell_to_cellin) + end + + @assert all_aggregated + + cell_to_cellin +end + +function _find_best_neighbor_space( + c1,c2,c3,c4,cell, + cell_to_faces, + face_to_cells, + cell_to_coords, + cell_to_touched, + cell_to_cellin, + facet_to_inoutcut, + loc) + + + faces = getindex!(c1,cell_to_faces,cell) + dmin = Inf + T = eltype(eltype(face_to_cells)) + best_neigh_cell = zero(T) + i_to_coords = getindex!(c3,cell_to_coords,cell) + for face in faces + inoutcut = facet_to_inoutcut[face] + if inoutcut != CUT && inoutcut != loc + continue + end + neigh_cells = getindex!(c2,face_to_cells,face) + for neigh_cell in neigh_cells + if neigh_cell != cell && cell_to_touched[neigh_cell] + cellin = cell_to_cellin[neigh_cell] + j_to_coords = getindex!(c4,cell_to_coords,cellin) + d = 0.0 + for p in i_to_coords + for q in j_to_coords + d = max(d,Float64(norm(p-q))) + end + end + if (1.0+1.0e-9)*d < dmin + dmin = d + best_neigh_cell = neigh_cell + end + end + end + end + best_neigh_cell +end + +function _touch_aggregated_cells_space!(cell_to_touched,cell_to_cellin) + for (cell,cellin) in enumerate(cell_to_cellin) + if cellin > 0 + cell_to_touched[cell] = true + end + end +end + +function color_aggregates_space(cell_to_cellin,model::DiscreteModel) + topo = get_grid_topology(model) + D = num_cell_dims(model) + cell_to_faces = get_faces(topo,D,D-1) + face_to_cells = get_faces(topo,D-1,D) + _color_aggregates_barrier_space(cell_to_cellin,cell_to_faces,face_to_cells) +end + +function _color_aggregates_barrier_space(cell_to_cellin,cell_to_faces,face_to_cells) + + n_cells = length(cell_to_cellin) + cell_to_isin = fill(false,n_cells) + for cellin in cell_to_cellin + if cellin>0 + cell_to_isin[cellin] = true + end + end + + n_incell = count(cell_to_isin) + g = SimpleGraph(n_incell) + + incell_to_cell = findall(cell_to_isin) + cell_to_incell = zeros(Int,n_cells) + cell_to_incell[cell_to_isin] .= 1:n_incell + + c1 = array_cache(cell_to_faces) + c2 = array_cache(face_to_cells) + + for cell in 1:n_cells + cellin = cell_to_cellin[cell] + if cellin > 0 + incell = cell_to_incell[cellin] + faces = getindex!(c1,cell_to_faces,cell) + for face in faces + neigh_cells = getindex!(c2,face_to_cells,face) + for neigh_cell in neigh_cells + neigh_cellin = cell_to_cellin[neigh_cell] + if neigh_cellin > 0 + neigh_incell = cell_to_incell[neigh_cellin] + if neigh_incell != incell + add_edge!(g,incell,neigh_incell) + end + end + end + end + end + end + + c = greedy_color(g) + + incell_to_nslaves = zeros(Int,n_incell) + for cell in 1:n_cells + cellin = cell_to_cellin[cell] + if cellin > 0 + incell = cell_to_incell[cellin] + incell_to_nslaves[incell] += 1 + end + end + + incell_to_color = c.colors + cell_to_color = zeros(Int,n_cells) + for cell in 1:n_cells + cellin = cell_to_cellin[cell] + if cellin > 0 + incell = cell_to_incell[cellin] + if incell_to_nslaves[incell] > 1 + color = incell_to_color[incell] + cell_to_color[cell] = color + end + end + end + + cell_to_color +end diff --git a/src/Exports.jl b/src/Exports.jl index 3e9c7261..2576f043 100644 --- a/src/Exports.jl +++ b/src/Exports.jl @@ -33,5 +33,8 @@ end @publish AgFEM AgFEMSpace @publish AgFEM aggregate +@publish AgFEM aggregatespace @publish AgFEM color_aggregates +@publish AgFEM color_aggregates_space @publish AgFEM AggregateAllCutCells +@publish AgFEM AggregateSpaceCutCells From a6a55c8a55b6c9c36b6c265bc36a5deb783b13d7 Mon Sep 17 00:00:00 2001 From: Hridya Dilip Date: Tue, 1 Sep 2020 17:56:25 +1000 Subject: [PATCH 2/8] Space Aggregation Filter Test --- .../CellAggregationSpaceFilterTest.jl | 82 +++++++++++++++++++ 1 file changed, 82 insertions(+) create mode 100644 test/AgFEMTests/CellAggregationSpaceFilterTest.jl diff --git a/test/AgFEMTests/CellAggregationSpaceFilterTest.jl b/test/AgFEMTests/CellAggregationSpaceFilterTest.jl new file mode 100644 index 00000000..07c65b70 --- /dev/null +++ b/test/AgFEMTests/CellAggregationSpaceFilterTest.jl @@ -0,0 +1,82 @@ +using Gridap +using GridapEmbedded +using GridapEmbedded.LevelSetCutters +using Gridap.Arrays +using Gridap.Geometry +using Gridap.ReferenceFEs +using Test + +# Background Mesh +domain = (-0.01,1.01,-0.01,1.01) +n_x = 2 +n_t = 2 +partition = (n_x,n_t) +bgmodel = CartesianDiscreteModel(domain,partition) + +# domain +geo = quadrilateral(x0=Point(0,0),d1=VectorValue(1,0),d2=VectorValue(0,1)) +cutgeo = cut(bgmodel,geo) +model = DiscreteModel(cutgeo) + +# AgFEM +#stategy = AggregateSpaceCutCells() +#aggregates = aggregatespace(strategy,cutgeo) + + +trian = get_triangulation(model) +cell_to_coords = get_cell_coordinates(trian) +topo = get_grid_topology(model) +D = num_cell_dims(model) + +n_cell = num_cells(topo) +filter = [false,false,true,true] +filters = fill(filter,n_cell) + +cell_to_faces = get_faces(topo,D,D-1) +face_to_cells = get_faces(topo,D-1,D) + +# Filter Cell_to_Faces +ptrs = collect(1:n_cell) +a = CompressedArray(cell_to_faces,ptrs) +fa = FilteredCellArray(a,filters) +@test fa == [[3,4],[4,7],[9,10],[10,12]] + +# Extracting only the faces to be aggregated +face = get_faces(topo,D-1,D-1) +n_faces = length(face) +filter_f = [false] +filters_f = fill(filter_f,n_faces) + +for i in 1:n_faces + for j in 1:length(fa) + for k in 1:2 + if face[i][1]==fa[j][k] + filters_f[i] = [true] + end + end + end +end + +ptrs_f = collect(1:n_faces) +b = CompressedArray(face,ptrs_f) +fb = FilteredCellArray(b,filters_f) +@test fb == [[],[],[3],[4],[],[],[7],[],[9],[10],[],[12]] + +# Filter Faces_to_cells +ptrs_fc = collect(1:n_faces) +c = CompressedArray(face_to_cells,ptrs_fc) + +filters_fc = deepcopy(filters_f) + +for i in 1:n_faces + l = length(face_to_cells[i]) + println(l) + for j in 2:l + append!(filters_fc[i],filters_f[i][1]) + println(filters_fc[i]) + end +end + +@test filters_fc == [[false],[false,false],[true],[true,true],[false],[false,false],[true],[false],[true],[true,true],[false],[true]] + +fc = FilteredCellArray(c,filters_fc) From 1eeb53bc88057f0c2790ed9fe0d9fec354c9b030 Mon Sep 17 00:00:00 2001 From: Hridya Dilip Date: Wed, 2 Sep 2020 12:42:19 +1000 Subject: [PATCH 3/8] Adding filter in color_aggregates --- src/AgFEM/CellAggregation_space.jl | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/src/AgFEM/CellAggregation_space.jl b/src/AgFEM/CellAggregation_space.jl index 9677fa4b..399cabe0 100644 --- a/src/AgFEM/CellAggregation_space.jl +++ b/src/AgFEM/CellAggregation_space.jl @@ -200,8 +200,16 @@ end function color_aggregates_space(cell_to_cellin,model::DiscreteModel) topo = get_grid_topology(model) D = num_cell_dims(model) + n_cell = num_cells(topo) + filter = [false,false,true,true] + filters = fill(filter,n_cell) + ptrs = collect(1:n_cell) cell_to_faces = get_faces(topo,D,D-1) face_to_cells = get_faces(topo,D-1,D) + comp_array = CompressedArray(cell_to_faces,ptrs) + filt_array = FilteredCellArray(comp_array,filters) + cell_to_faces = filt_array + _color_aggregates_barrier_space(cell_to_cellin,cell_to_faces,face_to_cells) end From bda68266bcecae9d61fbc80356115a48ad94c6f1 Mon Sep 17 00:00:00 2001 From: Hridya Dilip Date: Thu, 3 Sep 2020 14:49:42 +1000 Subject: [PATCH 4/8] Delete --- .../CellAggregationSpaceFilterTest.jl | 82 ------------------- 1 file changed, 82 deletions(-) delete mode 100644 test/AgFEMTests/CellAggregationSpaceFilterTest.jl diff --git a/test/AgFEMTests/CellAggregationSpaceFilterTest.jl b/test/AgFEMTests/CellAggregationSpaceFilterTest.jl deleted file mode 100644 index 07c65b70..00000000 --- a/test/AgFEMTests/CellAggregationSpaceFilterTest.jl +++ /dev/null @@ -1,82 +0,0 @@ -using Gridap -using GridapEmbedded -using GridapEmbedded.LevelSetCutters -using Gridap.Arrays -using Gridap.Geometry -using Gridap.ReferenceFEs -using Test - -# Background Mesh -domain = (-0.01,1.01,-0.01,1.01) -n_x = 2 -n_t = 2 -partition = (n_x,n_t) -bgmodel = CartesianDiscreteModel(domain,partition) - -# domain -geo = quadrilateral(x0=Point(0,0),d1=VectorValue(1,0),d2=VectorValue(0,1)) -cutgeo = cut(bgmodel,geo) -model = DiscreteModel(cutgeo) - -# AgFEM -#stategy = AggregateSpaceCutCells() -#aggregates = aggregatespace(strategy,cutgeo) - - -trian = get_triangulation(model) -cell_to_coords = get_cell_coordinates(trian) -topo = get_grid_topology(model) -D = num_cell_dims(model) - -n_cell = num_cells(topo) -filter = [false,false,true,true] -filters = fill(filter,n_cell) - -cell_to_faces = get_faces(topo,D,D-1) -face_to_cells = get_faces(topo,D-1,D) - -# Filter Cell_to_Faces -ptrs = collect(1:n_cell) -a = CompressedArray(cell_to_faces,ptrs) -fa = FilteredCellArray(a,filters) -@test fa == [[3,4],[4,7],[9,10],[10,12]] - -# Extracting only the faces to be aggregated -face = get_faces(topo,D-1,D-1) -n_faces = length(face) -filter_f = [false] -filters_f = fill(filter_f,n_faces) - -for i in 1:n_faces - for j in 1:length(fa) - for k in 1:2 - if face[i][1]==fa[j][k] - filters_f[i] = [true] - end - end - end -end - -ptrs_f = collect(1:n_faces) -b = CompressedArray(face,ptrs_f) -fb = FilteredCellArray(b,filters_f) -@test fb == [[],[],[3],[4],[],[],[7],[],[9],[10],[],[12]] - -# Filter Faces_to_cells -ptrs_fc = collect(1:n_faces) -c = CompressedArray(face_to_cells,ptrs_fc) - -filters_fc = deepcopy(filters_f) - -for i in 1:n_faces - l = length(face_to_cells[i]) - println(l) - for j in 2:l - append!(filters_fc[i],filters_f[i][1]) - println(filters_fc[i]) - end -end - -@test filters_fc == [[false],[false,false],[true],[true,true],[false],[false,false],[true],[false],[true],[true,true],[false],[true]] - -fc = FilteredCellArray(c,filters_fc) From 8de5774a2521d65c9f458f77adf8f8d3b4218f5e Mon Sep 17 00:00:00 2001 From: Hridya Dilip Date: Thu, 3 Sep 2020 15:02:31 +1000 Subject: [PATCH 5/8] Adding SpaceAggregationTest --- test/AgFEMTests/CellAggregationSpaceTest.jl | 35 +++++++++++++++++++++ 1 file changed, 35 insertions(+) create mode 100644 test/AgFEMTests/CellAggregationSpaceTest.jl diff --git a/test/AgFEMTests/CellAggregationSpaceTest.jl b/test/AgFEMTests/CellAggregationSpaceTest.jl new file mode 100644 index 00000000..0dc0bd0c --- /dev/null +++ b/test/AgFEMTests/CellAggregationSpaceTest.jl @@ -0,0 +1,35 @@ +module CellAggregationSpaceTests + +using Gridap +using GridapEmbedded +using GridapEmbedded.LevelSetCutters +using Gridap.Arrays +using Gridap.Geometry +using Gridap.ReferenceFEs +using Test + +# Background Mesh +domain = (0,2,0,1) +n_x = 10 +n_t = 10 +partition = (n_x,n_t) +bgmodel = CartesianDiscreteModel(domain,partition) + +# domain +geo = quadrilateral(x0=Point(0,0),d1=VectorValue(1,0),d2=VectorValue(1,1)) +cutgeo = cut(bgmodel,geo) +model = DiscreteModel(cutgeo) + +#AgFEM Strategy +strategy = AggregateSpaceCutCells() +aggregates = aggregatespace(strategy,cutgeo) + +colors = color_aggregates_space(aggregates,bgmodel) + +trian = Triangulation(bgmodel) + +#d = "./" +#fi = joinpath(d,"trian") +#writevtk(trian,fi,celldata=["colors"=>colors,"aggregates"=>aggregates]) + +end # module CellAggregationSpaceTests From 3e97e34f7c3de2cbda550ef50c490b094171c65a Mon Sep 17 00:00:00 2001 From: Hridya Dilip <53897026+diliphridya@users.noreply.github.com> Date: Thu, 3 Sep 2020 15:04:20 +1000 Subject: [PATCH 6/8] Rename CellAggregationSpaceTest.jl to CellAggregationSpaceTests.jl --- .../{CellAggregationSpaceTest.jl => CellAggregationSpaceTests.jl} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename test/AgFEMTests/{CellAggregationSpaceTest.jl => CellAggregationSpaceTests.jl} (100%) diff --git a/test/AgFEMTests/CellAggregationSpaceTest.jl b/test/AgFEMTests/CellAggregationSpaceTests.jl similarity index 100% rename from test/AgFEMTests/CellAggregationSpaceTest.jl rename to test/AgFEMTests/CellAggregationSpaceTests.jl From 3292be7332cb130e5b050acfeecd2f9cd853e81a Mon Sep 17 00:00:00 2001 From: Hridya Dilip <53897026+diliphridya@users.noreply.github.com> Date: Thu, 3 Sep 2020 15:05:34 +1000 Subject: [PATCH 7/8] Update AgFEM.jl --- src/AgFEM/AgFEM.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/AgFEM/AgFEM.jl b/src/AgFEM/AgFEM.jl index d7be6e2f..2cb5a348 100644 --- a/src/AgFEM/AgFEM.jl +++ b/src/AgFEM/AgFEM.jl @@ -26,5 +26,5 @@ include("CellAggregation.jl") include("AgFEMSpaces.jl") -include("CellAggregation_space.jl") +include("CellAggregationSpace.jl") end # module From 793c31b9b4367bcb8cc2142f924313658b8230bd Mon Sep 17 00:00:00 2001 From: Hridya Dilip <53897026+diliphridya@users.noreply.github.com> Date: Thu, 3 Sep 2020 15:06:05 +1000 Subject: [PATCH 8/8] Rename CellAggregation_space.jl to CellAggregationSpace.jl --- src/AgFEM/{CellAggregation_space.jl => CellAggregationSpace.jl} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename src/AgFEM/{CellAggregation_space.jl => CellAggregationSpace.jl} (100%) diff --git a/src/AgFEM/CellAggregation_space.jl b/src/AgFEM/CellAggregationSpace.jl similarity index 100% rename from src/AgFEM/CellAggregation_space.jl rename to src/AgFEM/CellAggregationSpace.jl