Merge pull request #234 from gridap/fixes_in_extended_fe_space

Fixes in extended fe space

src/FESpaces/FESpaceFactories.jl

@@ -28,27 +28,32 @@ function FESpace(;kwargs...)
 
   @assert fespace != nothing
 
-  if constraint == nothing
+  restricted_at = _get_restricted_triangulation(kwargs)
+  if restricted_at == nothing
     _fespace = fespace
+  else
+    @assert isa(restricted_at,RestrictedTriangulation)
+    _fespace = ExtendedFESpace(fespace,restricted_at)
+  end
+
+  if constraint == nothing
+    return _fespace
 
   elseif constraint == :zeromean
-    model = _get_kwarg(:model,kwargs)
     order = _get_kwarg(:order,kwargs)
-    trian = get_triangulation(model)
+    model = _get_kwarg(:model,kwargs,nothing)
+    if model === nothing
+      trian = _get_kwarg(:triangulation,kwargs)
+    else
+      trian = get_triangulation(model)
+    end
     quad = CellQuadrature(trian,order)
-    _fespace = ZeroMeanFESpace(fespace,trian,quad)
+    return ZeroMeanFESpace(_fespace,trian,quad)
 
   else
     @unreachable "Unknown constraint value $constraint"
   end
 
-  restricted_at = _get_restricted_triangulation(kwargs)
-  if restricted_at == nothing
-    return _fespace
-  else
-    @assert isa(restricted_at,RestrictedTriangulation)
-    return ExtendedFESpace(_fespace,restricted_at)
-  end
 
 end
 

src/FESpaces/ZeroMeanFESpaces.jl

@@ -29,7 +29,8 @@ function _setup_vols(V,trian,quad)
   U = TrialFESpace(V)
   assem = SparseMatrixAssembler(V,U)
   bh = get_cell_basis(V)
-  cellvec = integrate(bh,trian,quad)
+  bh_trian = restrict(bh,trian)
+  cellvec = integrate(bh_trian,trian,quad)
   cellids = get_cell_id(trian)
   vol_i = assemble_vector(assem,[cellvec],[cellids])
   vol = sum(vol_i)

src/Geometry/BoundaryTriangulations.jl

@@ -65,15 +65,20 @@ end
     BoundaryTriangulation(model::DiscreteModel,face_to_mask::Vector{Bool})
     BoundaryTriangulation(model::DiscreteModel)
 """
+function BoundaryTriangulation(model::DiscreteModel,face_to_mask::Vector{Bool},icell_around::Integer)
+  GenericBoundaryTriangulation(model,face_to_mask,icell_around)
+end
+
 function BoundaryTriangulation(model::DiscreteModel,face_to_mask::Vector{Bool})
-  GenericBoundaryTriangulation(model,face_to_mask)
+  icell_around = 1
+  BoundaryTriangulation(model,face_to_mask,icell_around)
 end
 
 function BoundaryTriangulation(model::DiscreteModel)
   topo = get_grid_topology(model)
   D = num_cell_dims(model)
   face_to_mask = collect(Bool,get_isboundary_face(topo,D-1))
-  GenericBoundaryTriangulation(model,face_to_mask)
+  BoundaryTriangulation(model,face_to_mask)
 end
 
 """

src/Geometry/DiscreteModels.jl

@@ -385,12 +385,12 @@ end
 
 """
 """
-function Triangulation(model,cell_to_oldcell::AbstractVector{<:Integer})
+function Triangulation(model::DiscreteModel,cell_to_oldcell::AbstractVector{<:Integer})
   oldtrian = Triangulation(model)
   RestrictedTriangulation(oldtrian,cell_to_oldcell)
 end
 
-function Triangulation(model,cell_to_mask::AbstractVector{Bool})
+function Triangulation(model::DiscreteModel,cell_to_mask::AbstractVector{Bool})
   oldtrian = Triangulation(model)
   RestrictedTriangulation(oldtrian,cell_to_mask)
 end

src/Geometry/Geometry.jl

@@ -135,6 +135,7 @@ export test_triangulation
 export restrict
 export get_physical_coordinate
 export get_cell_id
+export cell_measure
 
 export Grid
 export get_cell_nodes
@@ -239,8 +240,6 @@ include("DiscreteModels.jl")
 
 include("DiscreteModelPortions.jl")
 
-include("RestrictedDiscreteModels.jl")
-
 include("DiscreteModelMocks.jl")
 
 include("UnstructuredDiscreteModels.jl")
@@ -259,4 +258,6 @@ include("QPointCellFields.jl")
 
 include("AppendedTriangulations.jl")
 
+include("RestrictedDiscreteModels.jl")
+
 end # module

src/Geometry/RestrictedDiscreteModels.jl

@@ -41,3 +41,32 @@ function get_triangulation(model::RestrictedDiscreteModel)
   RestrictedTriangulation(model)
 end
 
+function Triangulation(::Type{ReferenceFE{d}},model::RestrictedDiscreteModel) where d
+  @notimplemented
+end
+
+function Triangulation(model::RestrictedDiscreteModel,cell_to_oldcell::AbstractVector{<:Integer})
+  @notimplemented
+end
+
+function Triangulation(model::RestrictedDiscreteModel,cell_to_mask::AbstractVector{Bool})
+  @notimplemented
+end
+
+function BoundaryTriangulation(model::RestrictedDiscreteModel,face_to_mask::Vector{Bool},icell_around::Integer)
+  d = num_cell_dims(model)-1
+  face_to_oldface = get_face_to_oldface(model,d)
+  oldmodel = get_oldmodel(model)
+  num_oldfaces = num_faces(oldmodel,d)
+  oldface_to_mask = fill(false,num_oldfaces)
+  oldface_to_mask[face_to_oldface] .= face_to_mask
+  BoundaryTriangulation(oldmodel,oldface_to_mask)
+end
+
+function InterfaceTriangulation(model_in::RestrictedDiscreteModel,model_out::RestrictedDiscreteModel)
+  cells_in = get_cell_to_oldcell(model_in)
+  cells_out = get_cell_to_oldcell(model_out)
+  model = get_oldmodel(model_in)
+  InterfaceTriangulation(model,cells_in,cells_out)
+end
+

src/Geometry/SkeletonTriangulations.jl

@@ -23,9 +23,9 @@ end
 """
 function SkeletonTriangulation(model::DiscreteModel,face_to_mask::Vector{Bool})
   left_cell_around = 1
-  left = GenericBoundaryTriangulation(model,face_to_mask,left_cell_around)
+  left = BoundaryTriangulation(model,face_to_mask,left_cell_around)
   right_cell_around = 2
-  right = GenericBoundaryTriangulation(model,face_to_mask,right_cell_around)
+  right = BoundaryTriangulation(model,face_to_mask,right_cell_around)
   SkeletonTriangulation(left,right)
 end
 
@@ -54,13 +54,6 @@ function InterfaceTriangulation(model::DiscreteModel,cells_in,cells_out)
   InterfaceTriangulation(model,cell_to_inout)
 end
 
-function InterfaceTriangulation(model_in::RestrictedDiscreteModel,model_out::RestrictedDiscreteModel)
-  cells_in = get_cell_to_oldcell(model_in)
-  cells_out = get_cell_to_oldcell(model_out)
-  model = get_oldmodel(model_in)
-  InterfaceTriangulation(model,cells_in,cells_out)
-end
-
 function InterfaceTriangulation(model::DiscreteModel,cell_to_inout::AbstractVector{<:Integer})
 
   D = num_cell_dims(model)

src/Geometry/Triangulations.jl

@@ -253,3 +253,24 @@ function _restrict_cell_field(r::AbstractArray,trian)
   GenericCellField(r,cm)
 end
 
+"""
+"""
+function cell_measure(trian_Ω1::Triangulation,n_bgcells::Integer)
+  trian_cut_1 = trian_Ω1
+  quad_cut_1 = CellQuadrature(trian_cut_1,0)
+  subcell1_to_dV = integrate(1,trian_cut_1,quad_cut_1)
+  subcell1_to_bgcell = get_cell_id(trian_cut_1)
+  bgcell_to_dV = zeros(n_bgcells)
+  _meas_K_fill!(bgcell_to_dV,subcell1_to_dV,subcell1_to_bgcell)
+  bgcell_to_dV
+end
+
+function cell_measure(trian_Ω1::Triangulation,trian::Triangulation)
+  cell_measure(trian_Ω1,num_cells(trian))
+end
+
+function _meas_K_fill!(bgcell_to_dV,subcell1_to_dV,subcell1_to_bgcell)
+  for (subcell1, bgcell) in enumerate(subcell1_to_bgcell)
+    bgcell_to_dV[bgcell] += subcell1_to_dV[subcell1]
+  end
+end

test/FESpacesTests/ExtendedFESpacesTests.jl

@@ -109,16 +109,25 @@ collect(evaluate(jump(ε(uh_Γ)),q_Γ))
 V = TestFESpace(model=model_in,valuetype=Float64,reffe=:Lagrangian,order=2,conformity=:H1)
 @test isa(V,ExtendedFESpace)
 
+V = TestFESpace(model=model_in,valuetype=Float64,reffe=:Lagrangian,order=2,conformity=:H1,constraint=:zeromean)
+uh = FEFunction(V,rand(num_free_dofs(V)))
+uh_in = restrict(uh,trian_in)
+@test sum(integrate(uh_in,trian_in,quad_in)) + 1 ≈ 1
+
 V = TestFESpace(model=model,valuetype=Float64,reffe=:Lagrangian,order=2,conformity=:H1)
 @test !isa(V,ExtendedFESpace)
 
 V = TestFESpace(triangulation=trian_in,valuetype=Float64,reffe=:Lagrangian,order=2,conformity=:L2)
 @test isa(V,ExtendedFESpace)
 
+V = TestFESpace(triangulation=trian_in,valuetype=Float64,reffe=:Lagrangian,order=2,conformity=:L2,constraint=:zeromean)
+uh = FEFunction(V,rand(num_free_dofs(V)))
+uh_in = restrict(uh,trian_in)
+@test sum(integrate(uh_in,trian_in,quad_in)) + 1 ≈ 1
+
 V = TestFESpace(triangulation=trian,valuetype=Float64,reffe=:Lagrangian,order=2,conformity=:L2)
 @test !isa(V,ExtendedFESpace)
 
-
 #using Gridap.Visualization
 #writevtk(trian,"trian",cellfields=["uh"=>uh])
 #writevtk(trian_in,"trian_in",cellfields=["uh"=>uh_in])

test/GeometryTests/RestrictedDiscreteModelsTests.jl

@@ -68,4 +68,25 @@ itrian = InterfaceTriangulation(model_fluid,model_solid)
 #writevtk(model,"model")
 #writevtk(trian,"trian")
 
+oldtrian = Triangulation(oldmodel)
+trian_s = SkeletonTriangulation(model)
+trian_b = BoundaryTriangulation(model)
+
+meas_K_b = cell_measure(trian_b,oldtrian)
+meas_K_sl = cell_measure(trian_s.left,oldtrian)
+meas_K_sr = cell_measure(trian_s.right,oldtrian)
+
+oldcell_to_cell = zeros(Int,num_cells(oldmodel))
+oldcell_to_cell[cell_to_oldcell] .= 1:length(cell_to_oldcell)
+
+@test all( oldcell_to_cell[findall(meas_K_b .!= 0)] .!= 0 )
+@test all( oldcell_to_cell[findall(meas_K_sl .!= 0)] .!= 0 )
+@test all( oldcell_to_cell[findall(meas_K_sr .!= 0)] .!= 0 )
+
+#using Gridap.Visualization
+#writevtk(trian_b,"trian_b")
+#writevtk(trian_s,"trian_s")
+#writevtk(trian,"trian")
+#writevtk(oldtrian,"oldtrian",celldata=["meas_K_b"=>meas_K_b,"meas_K_sl"=>meas_K_sl,"meas_K_sr"=>meas_K_sr])
+
 end # module