bugfix at ref space, MomentFitting tests completes gridap#22

momentfiting does not work in simplex bgmeshes

gridap/GridapEmbedded.jl#74

src/Embedded.jl

@@ -134,6 +134,7 @@ end
 function _cut_stl(model::DiscreteModel,geom::STLGeometry;kwargs...)
   subcell_grid, subface_grid, bsubface_grid, labels = subtriangulate(model,geom;kwargs...)
 
+  D = num_dims(model)
   inout_dict = Dict{Int8,Int8}(
     FACE_IN => IN, FACE_OUT => OUT, FACE_CUT => CUT, UNSET => OUT )
 
@@ -154,18 +155,15 @@ function _cut_stl(model::DiscreteModel,geom::STLGeometry;kwargs...)
   subfacets = SubFacetData(data...)
 
   bface_to_bgcell = labels.bface_to_bgcell
-  bface_to_points = get_cell_node_ids( bsubface_grid )
-  point_to_coords = get_node_coordinates( bsubface_grid )
-
-
-  bface_to_bgface = boundary_facet_to_background_facet(
-    model,bface_to_bgcell,labels.bface_to_lbgface)
-  ftrian = Grid(ReferenceFE{num_dims(model)-1},model)
-  point_to_rcoords = send_to_ref_space(ftrian,bface_to_bgface,bsubface_grid)
+  bface_to_lbgface = labels.bface_to_lbgface
+  bface_to_bgface = boundary_facet_to_background_facet(model,bface_to_bgcell,bface_to_lbgface)
   newbfaces = unique_boundary_facets(bface_to_bgcell,bface_to_bgface)
-  bface_to_points = bface_to_points[newbfaces]
+  bsubface_grid = isolate_cell_coordinates(bsubface_grid,newbfaces)
   bface_to_bgface = bface_to_bgface[newbfaces]
   bface_to_io = labels.bface_to_io[newbfaces]
+  bface_to_points = get_cell_node_ids( bsubface_grid )
+  point_to_coords = get_node_coordinates( bsubface_grid )
+  point_to_rcoords = send_to_ref_space(Val{D-1}(),model,bface_to_bgface,bsubface_grid)
   bface_to_io = [ replace( bface_to_io, inout_dict... ) ]
   data = bface_to_points,bface_to_bgface,point_to_coords,point_to_rcoords
   bsubfacets = SubCellData(data...)
@@ -180,6 +178,18 @@ function _cut_stl(model::DiscreteModel,geom::STLGeometry;kwargs...)
   (bgface_to_ioc,bsubfacets,bface_to_io,oid_to_ls)
 end
 
+
+function isolate_cell_coordinates(grid::Grid,cells=1:num_cells(grid))
+  cell_coords = Table(get_cell_coordinates(grid)[cells])
+  coords = cell_coords.data
+  data = 1:length(cell_coords.data)
+  ptrs = cell_coords.ptrs
+  cell_nodes = Table(data,ptrs)
+  reffes = get_reffes(grid)
+  cell_types = get_cell_type(grid)[cells]
+  UnstructuredGrid(coords,cell_nodes,reffes,cell_types)
+end
+
 function boundary_facet_to_background_facet(model::DiscreteModel,f_to_bgc,f_to_lbgf)
   topo = get_grid_topology(model)
   D = num_dims(model)
@@ -213,6 +223,16 @@ function send_to_ref_space(
   send_to_ref_space(get_grid(model),cell_to_bgcell,subgrid)
 end
 
+function send_to_ref_space(
+  ::Val{D},
+  model::DiscreteModel,
+  cell_to_bgcell::Vector,
+  subgrid::Grid) where D
+
+  grid = Grid(ReferenceFE{D},model)
+  send_to_ref_space(grid,cell_to_bgcell,subgrid)
+end
+
 function send_to_ref_space(grid::Grid,cell_to_bgcell::Vector,subgrid::Grid)
   bgcell_map = get_cell_map(grid)
   bgcell_invmap = lazy_map(inverse_map,bgcell_map)

test/CutFacetsTests.jl

@@ -4,6 +4,8 @@ using Gridap
 using GridapEmbedded
 using STLCutters
 using Test
+using Gridap.ReferenceFEs
+
 
 geo = STLGeometry(joinpath(@__DIR__,"data/cube.stl"))
 geoₐ = cube(x0=Point(0.0,0.0,0.5))
@@ -103,7 +105,6 @@ dΛoutₐ⁻ = Measure(Λoutₐ.⁻,degree)
 @test sum(∫(1)dΛout⁻) ≈ sum(∫(1)dΛoutₐ⁻)
 
 
-using Gridap.ReferenceFEs
 cutgeo = cut(model, geo)
 Ω = Triangulation(cutgeo,PHYSICAL)
 Ω_act_in = Triangulation(cutgeo,ACTIVE_IN,geo)
@@ -113,8 +114,6 @@ dΩᵐ_out = Measure(Ω_act_out,Quadrature(momentfitted,cutgeo,degree,in_or_out=
 
 dΩ = Measure(Ω,degree)
 f = x -> x[1] + 1
-f = 1
-# @test ∑(∫(f)dΩᵐ_in) ≈ ∑(∫(f)dΩ)
-
+@test ∑(∫(f)dΩᵐ_in) ≈ ∑(∫(f)dΩ)
 
 end # module

test/Integration.jl

@@ -1,4 +1,4 @@
-#module IntegrationTests
+module IntegrationTests
 
 using STLCutters
 using Gridap
@@ -77,8 +77,7 @@ dΩᵐ_in = Measure(Ω_act_in,Quadrature(momentfitted,cutgeo,degree,in_or_out=IN
 dΩᵐ_out = Measure(Ω_act_out,Quadrature(momentfitted,cutgeo,degree,in_or_out=OUT))
 
 f = x -> x[1] + 1
-f = 1
-# @test ∑(∫(f)dΩᵐ_in) ≈ ∑(∫(f)dΩ)
+@test ∑(∫(f)dΩᵐ_in) ≈ ∑(∫(f)dΩ)
 
 # Simplex background
 #
@@ -130,4 +129,7 @@ b = sum( ∫( v⋅n_Γd⋅∇(u) ) * dΓd )
 
 # Moment fitted
 
-#end
+# Not working until GridapEmbedded #74 is fixed
+# https://github.com/gridap/GridapEmbedded.jl/issues/74
+
+end