Merge pull request #68 from pmartorell/moment_fitted

Moment fitted machinery

.github/workflows/ci.yml

@@ -8,7 +8,7 @@ jobs:
       fail-fast: false
       matrix:
         version:
-          - '1.6'
+          - '1.8'
         os:
           - ubuntu-latest
         arch:
@@ -42,7 +42,7 @@ jobs:
       fail-fast: false
       matrix:
         version:
-          - '1.6'
+          - '1.8'
         os:
           - ubuntu-latest
         arch:
@@ -77,7 +77,7 @@ jobs:
       - uses: actions/checkout@v2
       - uses: julia-actions/setup-julia@v1
         with:
-          version: '1.6'
+          version: '1.8'
       - run: |
           julia --project=docs -e '
             using Pkg

.github/workflows/ci_x86.yml

@@ -14,7 +14,7 @@ jobs:
       fail-fast: false
       matrix:
         version:
-          - '1.6'
+          - '1.8'
         os:
           - ubuntu-latest
         arch:

.gitignore

@@ -6,4 +6,5 @@
 /docs/build/
 /docs/site/
 Manifest.toml
-.vscode/
+.vscode/
+*.swp

NEWS.md

@@ -4,6 +4,11 @@ All notable changes to this project will be documented in this file.
 The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
+## [0.8.2] - 2023-08-22
+
+### Added
+- Moment fitted machinery. Since PR [#68](https://github.com/gridap/GridapEmbedded.jl/pull/68).
+
 ## [0.8.0] - 2021-11-03
 
 ### Changed

src/AgFEM/AggregateBoundingBoxes.jl

@@ -1,25 +1,30 @@
+function init_bboxes(cell_to_coords,args...;kwargs...)
+  @abstractmethod
+end
+
 function init_bboxes(cell_to_coords)
   # RMK: Assuming first node is min and last node is max of BBox
-  [ [cell_to_coords[c][1],cell_to_coords[c][end]] for c in 1:length(cell_to_coords) ]
+  map(x-> [ x[1], x[end] ], vec(cell_to_coords) )
 end
 
-function init_bboxes(cell_to_coords,cut::EmbeddedDiscretization)
+function init_bboxes(cell_to_coords,cut::EmbeddedDiscretization;in_or_out=IN)
   bgcell_to_cbboxes = init_bboxes(cell_to_coords)
   cut_bgtrian = Triangulation(cut,CUT,cut.geo)
   cut_bgmodel = get_active_model(cut_bgtrian)
   ccell_to_bgcell = get_cell_to_parent_cell(cut_bgmodel)
-  ccell_to_cbboxes = init_cut_bboxes(cut,ccell_to_bgcell)
+  ccell_to_cbboxes = init_cut_bboxes(cut,ccell_to_bgcell,in_or_out)
   for (cc,cbb) in enumerate(ccell_to_cbboxes)
     bgcell_to_cbboxes[ccell_to_bgcell[cc]] = cbb
   end
   bgcell_to_cbboxes
 end
 
-function init_cut_bboxes(cut,ccell_to_bgcell)
-  subcell_to_inout   = compute_subcell_to_inout(cut,cut.geo) # WARNING! This function has a bug
-  subcell_is_in      = lazy_map(i->i==IN,subcell_to_inout)
+function init_cut_bboxes(cut,ccell_to_bgcell,in_or_out)
+  subcell_to_inout   = compute_subcell_to_inout(cut,cut.geo)
+  subcell_is_in      = lazy_map(i->i==in_or_out,subcell_to_inout)
+  subcell_to_bgcell  = cut.subcells.cell_to_bgcell
   inscell_to_subcell = findall(subcell_is_in)
-  inscell_to_bgcell  = lazy_map(Reindex(cut.subcells.cell_to_bgcell),inscell_to_subcell)
+  inscell_to_bgcell  = lazy_map(Reindex(subcell_to_bgcell),inscell_to_subcell)
   inscell_to_bboxes  = init_subcell_bboxes(cut,inscell_to_subcell)
   lazy_map(ccell_to_bgcell) do bg
     bg_to_scs = findall(map(x->x==bg,inscell_to_bgcell))
@@ -29,12 +34,12 @@ function init_cut_bboxes(cut,ccell_to_bgcell)
 end
 
 function init_subcell_bboxes(cut,inscell_to_subcell)
-  subcell_to_points = lazy_map(Reindex(cut.subcells.cell_to_points),inscell_to_subcell)
+  subcell_to_points = cut.subcells.cell_to_points
   point_to_coords   = cut.subcells.point_to_coords
-  subcell_to_coords = lazy_map(subcell_to_points) do points
-    point_to_coords[points]
-  end
-  lazy_map(compute_subcell_bbox,subcell_to_coords)
+  inscell_to_points = lazy_map(Reindex(subcell_to_points),inscell_to_subcell)
+  inscell_to_coords = lazy_map( #
+    Broadcasting(Reindex(point_to_coords)),inscell_to_points)
+  lazy_map(compute_subcell_bbox,inscell_to_coords)
 end
 
 function compute_subcell_bbox(subcell_to_coords)
@@ -71,35 +76,25 @@ function reset_bboxes_at_cut_cells!(root_to_agg_bbox,
   end
 end
 
-function compute_cell_bboxes(model::DiscreteModelPortion,cell_to_root)
-  compute_cell_bboxes(get_parent_model(model),cell_to_root)
-end
-
-function compute_cell_bboxes(model::DiscreteModel,cell_to_root)
-  trian = Triangulation(model)
-  compute_cell_bboxes(trian,cell_to_root)
+function compute_cell_bboxes(model,cell_to_root,args...;kwargs...)
+  @abstractmethod
 end
 
-function compute_cell_bboxes(trian::Triangulation,cell_to_root)
-  cell_to_coords = get_cell_coordinates(trian)
-  root_to_agg_bbox = init_bboxes(cell_to_coords)
-  compute_bboxes!(root_to_agg_bbox,cell_to_root)
-  reset_bboxes_at_cut_cells!(root_to_agg_bbox,cell_to_root,cell_to_coords)
-  root_to_agg_bbox
+function compute_cell_bboxes(model::DiscreteModelPortion,
+                             cell_to_root,args...;kwargs...)
+  compute_cell_bboxes(get_parent_model(model),cell_to_root,args...;kwargs...)
 end
 
-function compute_cell_bboxes(model::DiscreteModelPortion,cut::EmbeddedDiscretization,cell_to_root)
-  compute_cell_bboxes(get_parent_model(model),cut,cell_to_root)
-end
-
-function compute_cell_bboxes(model::DiscreteModel,cut::EmbeddedDiscretization,cell_to_root)
+function compute_cell_bboxes(model::DiscreteModel,
+                             cell_to_root,args...;kwargs...)
   trian = Triangulation(model)
-  compute_cell_bboxes(trian,cut,cell_to_root)
+  compute_cell_bboxes(trian,cell_to_root,args...;kwargs...)
 end
 
-function compute_cell_bboxes(trian::Triangulation,cut::EmbeddedDiscretization,cell_to_root)
+function compute_cell_bboxes(trian::Triangulation,
+                             cell_to_root,args...;kwargs...)
   cell_to_coords = get_cell_coordinates(trian)
-  root_to_agg_bbox = init_bboxes(cell_to_coords,cut)
+  root_to_agg_bbox = init_bboxes(cell_to_coords,args...;kwargs...)
   compute_bboxes!(root_to_agg_bbox,cell_to_root)
   reset_bboxes_at_cut_cells!(root_to_agg_bbox,cell_to_root,cell_to_coords)
   root_to_agg_bbox
@@ -111,8 +106,9 @@ function compute_bbox_dfaces(model::DiscreteModel,cell_to_agg_bbox)
   D = num_dims(gt)
   for d = 1:D-1
     dface_to_Dfaces = get_faces(gt,d,D)
-    d_bboxes =
-      [ _compute_bbox_dface(dface_to_Dfaces,cell_to_agg_bbox,face) for face in 1:num_faces(gt,d) ]
+    d_bboxes = map(1:num_faces(gt,d)) do face
+      _compute_bbox_dface(dface_to_Dfaces,cell_to_agg_bbox,face)
+    end
     bboxes = push!(bboxes,d_bboxes)
   end
   bboxes = push!(bboxes,cell_to_agg_bbox)
@@ -121,16 +117,18 @@ end
 function _compute_bbox_dface(dface_to_Dfaces,cell_to_agg_bbox,i)
   cells_around_dface_i = getindex(dface_to_Dfaces,i)
   bboxes_around_dface_i = cell_to_agg_bbox[cells_around_dface_i]
-  bbmins = [ bboxes_around_dface_i[i][1].data for i in 1:length(bboxes_around_dface_i) ]
-  bbmaxs = [ bboxes_around_dface_i[i][2].data for i in 1:length(bboxes_around_dface_i) ]
+  bbmins = map( i->i[1].data, bboxes_around_dface_i )
+  bbmaxs = map( i->i[2].data, bboxes_around_dface_i )
   [min.(bbmins...),max.(bbmaxs...)]
 end
 
 function _compute_cell_to_dface_bboxes(model::DiscreteModel,dbboxes)
   gt = get_grid_topology(model)
   trian = Triangulation(model)
   ctc = get_cell_coordinates(trian)
-  bboxes = [ __compute_cell_to_dface_bboxes(gt,ctc,dbboxes,cell) for cell in 1:num_cells(model) ]
+  bboxes = map(1:num_cells(model)) do cell
+    __compute_cell_to_dface_bboxes(gt,ctc,dbboxes,cell)
+  end
   CellPoint(bboxes,trian,PhysicalDomain()).cell_ref_point
 end
 
@@ -150,24 +148,20 @@ function __compute_cell_to_dface_bboxes(gt::GridTopology,ctc,dbboxes,cell::Int)
   cdbboxes = vcat(cdbboxes,[dbboxes[D][cell][1],dbboxes[D][cell][end]])
 end
 
-function compute_cell_to_dface_bboxes(model::DiscreteModel,cell_to_root)
-  cbboxes = compute_cell_bboxes(model,cell_to_root)
-  dbboxes = compute_bbox_dfaces(model,cbboxes)
-  _compute_cell_to_dface_bboxes(model,dbboxes)
-end
-
-function compute_cell_to_dface_bboxes(model::DiscreteModelPortion,cell_to_root)
-  bboxes = compute_cell_to_dface_bboxes(get_parent_model(model),cell_to_root)
-  bboxes[get_cell_to_parent_cell(model)]
+function compute_cell_to_dface_bboxes(model,cell_to_root,args...;kwargs...)
+  @abstractmethod
 end
 
-function compute_cell_to_dface_bboxes(model::DiscreteModel,cut::EmbeddedDiscretization,cell_to_root)
-  cbboxes = compute_cell_bboxes(model,cut,cell_to_root)
+function compute_cell_to_dface_bboxes(model::DiscreteModel,
+                                      cell_to_root,args...;kwargs...)
+  cbboxes = compute_cell_bboxes(model,cell_to_root,args...;kwargs...)
   dbboxes = compute_bbox_dfaces(model,cbboxes)
   _compute_cell_to_dface_bboxes(model,dbboxes)
 end
 
-function compute_cell_to_dface_bboxes(model::DiscreteModelPortion,cut::EmbeddedDiscretization,cell_to_root)
-  bboxes = compute_cell_to_dface_bboxes(get_parent_model(model),cut,cell_to_root)
+function compute_cell_to_dface_bboxes(model::DiscreteModelPortion,
+                                      cell_to_root,args...;kwargs...)
+  pmodel = get_parent_model(model)
+  bboxes = compute_cell_to_dface_bboxes(pmodel,cell_to_root,args...;kwargs...)
   bboxes[get_cell_to_parent_cell(model)]
 end

src/Exports.jl

@@ -47,4 +47,6 @@ end
 @publish AgFEM AggregateCutCellsByThreshold
 @publish AgFEM AggregateAllCutCells
 @publish AgFEM compute_cell_bboxes
-@publish AgFEM compute_cell_to_dface_bboxes
+@publish AgFEM compute_cell_to_dface_bboxes
+
+@publish MomentFittedQuadratures momentfitted

src/GridapEmbedded.jl

@@ -8,6 +8,8 @@ include("LevelSetCutters/LevelSetCutters.jl")
 
 include("AgFEM/AgFEM.jl")
 
+include("MomentFittedQuadratures/MomentFittedQuadratures.jl")
+
 include("Exports.jl")
 
 end # module

src/Interfaces/EmbeddedDiscretizations.jl

@@ -1,5 +1,7 @@
 
-struct EmbeddedDiscretization{Dp,T} <: GridapType
+abstract type AbstractEmbeddedDiscretization <: GridapType end
+
+struct EmbeddedDiscretization{Dp,T} <: AbstractEmbeddedDiscretization
   bgmodel::DiscreteModel
   ls_to_bgcell_to_inoutcut::Vector{Vector{Int8}}
   subcells::SubCellData{Dp,Dp,T}
@@ -46,6 +48,14 @@ function DiscreteModel(cut::EmbeddedDiscretization,geo::CSG.Geometry,in_or_out)
   #DiscreteModel(cut.bgmodel,cell_list)
 end
 
+function get_background_model(cut::EmbeddedDiscretization)
+  cut.bgmodel
+end
+
+function get_geometry(cut::EmbeddedDiscretization)
+  cut.geo
+end
+
 function compute_bgcell_to_inoutcut(cut::EmbeddedDiscretization,geo::CSG.Geometry)
 
   tree = get_tree(geo)

src/Interfaces/EmbeddedFacetDiscretizations.jl

@@ -1,5 +1,5 @@
 
-struct EmbeddedFacetDiscretization{Dc,Dp,T} <: GridapType
+struct EmbeddedFacetDiscretization{Dc,Dp,T} <: AbstractEmbeddedDiscretization
   bgmodel::DiscreteModel{Dp,Dp}
   ls_to_facet_to_inoutcut::Vector{Vector{Int8}}
   subfacets::SubCellData{Dc,Dp,T}
@@ -8,6 +8,14 @@ struct EmbeddedFacetDiscretization{Dc,Dp,T} <: GridapType
   geo::CSG.Geometry
 end
 
+function get_background_model(cut::EmbeddedFacetDiscretization)
+  cut.bgmodel
+end
+
+function get_geometry(cut::EmbeddedFacetDiscretization)
+  cut.geo
+end
+
 function SkeletonTriangulation(cut::EmbeddedFacetDiscretization)
   SkeletonTriangulation(cut,PHYSICAL_IN)
 end
@@ -243,7 +251,7 @@ struct SubFacetBoundaryTriangulation{Dc,Dp,T} <: Triangulation{Dc,Dp}
     subfacet_to_facet_map = _setup_cell_ref_map(subfacets,subgrid)
     face_ref_map = lazy_map(Reindex(glue.tface_to_mface_map),subfacet_to_facet)
     cell_ref_map = lazy_map(∘,face_ref_map,subfacet_to_facet_map)
-    
+
     new{Dc,Dp,T}(
       facets,
       subfacets,

src/Interfaces/Interfaces.jl

@@ -8,6 +8,7 @@ using Gridap.Geometry
 using Gridap.CellData
 using Gridap.Visualization
 
+import GridapEmbedded.CSG: get_geometry
 import Gridap.Geometry: UnstructuredGrid
 import Gridap.Geometry: BoundaryTriangulation
 import Gridap.Geometry: SkeletonTriangulation

src/LevelSetCutters/LevelSetCutters.jl

@@ -98,6 +98,14 @@ function cut_facets(background::DiscreteModel,geom::DiscreteGeometry)
   cut_facets(cutter,background,geom)
 end
 
+function cut_facets(cut::EmbeddedDiscretization,geo::CSG.Geometry)
+  cut_facets(cut.bgmodel,geo)
+end
+
+function cut_facets(cut::EmbeddedFacetDiscretization,args...)
+  cut
+end
+
 function _cut_ls_facets(model::DiscreteModel,geom)
   D = num_cell_dims(model)
   grid = Grid(ReferenceFE{D-1},model)