Merge pull request #1069 from gridap/constant-fespaces

ConstantFESpaces on sub-triangulations

NEWS.md

@@ -14,6 +14,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 ### Changed
 
 - Low level optimisations to reduce allocations. `AffineMap` renamed to `AffineField`. New `AffineMap <: Map`, doing the same as `AffineField` without struct allocation. New `ConstantMap <: Map`, doing the same as `ConstantField` without struct allocation. Since PR[#1043](https://github.com/gridap/Gridap.jl/pull/1043).
+- `ConstantFESpaces` can now be built on triangulations. Since PR[#1069](https://github.com/gridap/Gridap.jl/pull/1069)
 
 ## [0.18.8] - 2024-12-2
 

src/FESpaces/ConstantFESpaces.jl

@@ -1,49 +1,63 @@
 
 """
-struct ConstantFESpace <: SingleFieldFESpace
-  # private fields
-end
+    struct ConstantFESpace <: SingleFieldFESpace
+      # private fields
+    end
+
+    ConstantFESpace(model::DiscreteModel; vector_type=Vector{Float64}, field_type=Float64)
+    ConstantFESpace(trian::Triangulation; vector_type=Vector{Float64}, field_type=Float64)
+
+FESpace that is constant over the provided model/triangulation. Typically used as  
+lagrange multipliers. The kwargs `vector_type` and `field_type` are used to specify the
+types of the dof-vector and dof-value respectively.
 """
 struct ConstantFESpace{V,T,A,B,C} <: SingleFieldFESpace
-  model::DiscreteModel
+  trian::Triangulation
   cell_basis::A
   cell_dof_basis::B
   cell_dof_ids::C
 
-  function ConstantFESpace(model;
-                         vector_type::Type{V}=Vector{Float64},
-                         field_type::Type{T}=Float64) where {V,T}
-    function setup_cell_reffe(model::DiscreteModel,
-        reffe::Tuple{<:ReferenceFEName,Any,Any}; kwargs...)
-      basis, reffe_args,reffe_kwargs = reffe
-      cell_reffe = ReferenceFE(model,basis,reffe_args...;reffe_kwargs...)
-    end
+  function ConstantFESpace(
+    model::DiscreteModel{Dc},
+    trian::Triangulation{Dc};
+    vector_type::Type{V}=Vector{Float64},
+    field_type::Type{T}=Float64
+  ) where {Dc,V,T}
+    @assert num_cells(model) == num_cells(trian)
 
-    reffe = ReferenceFE(lagrangian,T,0)
-    cell_reffe = setup_cell_reffe(model,reffe)
+    basis, reffe_args, reffe_kwargs = ReferenceFE(lagrangian,T,0)
+    cell_reffe = ReferenceFE(model,basis,reffe_args...;reffe_kwargs...)
     cell_basis_array = lazy_map(get_shapefuns,cell_reffe)
 
     cell_basis = SingleFieldFEBasis(
-      cell_basis_array,
-      Triangulation(model),
-      TestBasis(),
-      ReferenceDomain())
+      cell_basis_array, trian, TestBasis(), ReferenceDomain()
+    )
+    cell_dof_basis = CellDof(
+      lazy_map(get_dof_basis,cell_reffe),trian,ReferenceDomain()
+    )
+    cell_dof_ids = Fill(Int32(1):Int32(num_indep_components(field_type)),num_cells(trian))
 
-    cell_dof_basis_array = lazy_map(get_dof_basis,cell_reffe)
-    cell_dof_basis = CellDof(cell_dof_basis_array,Triangulation(model),ReferenceDomain())
-
-    cell_dof_ids = Fill(Int32(1):Int32(num_indep_components(field_type)),num_cells(model))
     A = typeof(cell_basis)
     B = typeof(cell_dof_basis)
     C = typeof(cell_dof_ids)
-    new{V,T,A,B,C}(model, cell_basis, cell_dof_basis, cell_dof_ids)
+    new{V,T,A,B,C}(trian, cell_basis, cell_dof_basis, cell_dof_ids)
   end
 end
 
+function ConstantFESpace(model::DiscreteModel; kwargs...)
+  trian = Triangulation(model)
+  ConstantFESpace(model,trian; kwargs...)
+end
+
+function ConstantFESpace(trian::Triangulation; kwargs...)
+  model = get_active_model(trian)
+  ConstantFESpace(model,trian; kwargs...)
+end
+
 TrialFESpace(f::ConstantFESpace) = f
 
 # Delegated functions
-get_triangulation(f::ConstantFESpace) = Triangulation(f.model)
+get_triangulation(f::ConstantFESpace) = f.trian
 
 ConstraintStyle(::Type{<:ConstantFESpace}) = UnConstrained()
 

test/FESpacesTests/ConstantFESpaceTests.jl

@@ -6,9 +6,9 @@ using Test
 domain = (0,1,0,1)
 partition = (4,4)
 model = CartesianDiscreteModel(domain,partition)
-Λ=ConstantFESpace(model)
+Λ = ConstantFESpace(model)
 Gridap.FESpaces.test_fe_space(Λ)
-M=TrialFESpace(Λ)
+M = TrialFESpace(Λ)
 
 order = 2
 u((x,y)) = (x+y)^order
@@ -29,13 +29,17 @@ uh = solve(op)
 @assert sum(∫((uh[1]-u)*(uh[1]-u))dΩ) < 1.0e-14
 abs(sum(∫(uh[2])dΩ)) < 1.0e-12
 
-Λ2=ConstantFESpace(model,field_type=VectorValue{2,Float64})
+Λ2 = ConstantFESpace(model,field_type=VectorValue{2,Float64})
 Gridap.FESpaces.test_fe_space(Λ2)
-M2=TrialFESpace(Λ2)
+M2 = TrialFESpace(Λ2)
 a2(μ,λ) = ∫(λ⋅μ)dΩ
 l2(λ) = ∫(VectorValue(0.0,0.0)⋅λ)dΩ
 op2 = AffineFEOperator(a2,l2,M2,Λ2)
 μ2h = solve(op2)
 @assert sum(∫(μ2h⋅μ2h)dΩ) < 1.0e-12
 
+trian = Triangulation(model,[1,2,3,4])
+Λ3 = ConstantFESpace(trian,field_type=VectorValue{2,Float64})
+Gridap.FESpaces.test_fe_space(Λ3)
+
 end # module