[BREAKING] Embedded elements

CHANGELOG.md

@@ -8,10 +8,16 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 ## [Unreleased]
 ### Added
 - `DofHandler` now supports fields on subdomains and mixed grids. ([#667][github-667])
+- Support for embedded elements in `CellValues`. ([#651][github-651])
+- Support for vectorized interpolations with dimension different from the spatial 
+  dimension in `CellValues`. ([#651][github-651])
+
 ### Removed
 - **BREAKING**: `MixedDofHandler` has been renamed to `DofHandler`. Update by replacing 
 `MixedDofHandler` to `DofHandler`. ([#667][github-667])
-
+- **BREAKING**: All `CellValues` have a finer granularity on the different types of
+  dimensions used. Custom `CellValues` must be updated to respect the new type 
+  parameters. ([#651][github-651])
 
 ## [0.3.14] - 2023-04-03
 ### Added
@@ -414,6 +420,7 @@ poking into Ferrite internals:
 [github-645]: https://github.com/Ferrite-FEM/Ferrite.jl/pull/645
 [github-647]: https://github.com/Ferrite-FEM/Ferrite.jl/pull/647
 [github-650]: https://github.com/Ferrite-FEM/Ferrite.jl/pull/650
+[github-651]: https://github.com/Ferrite-FEM/Ferrite.jl/pull/651
 [github-653]: https://github.com/Ferrite-FEM/Ferrite.jl/pull/653
 [github-656]: https://github.com/Ferrite-FEM/Ferrite.jl/pull/656
 [github-658]: https://github.com/Ferrite-FEM/Ferrite.jl/pull/658

Project.toml

@@ -9,6 +9,7 @@ NearestNeighbors = "b8a86587-4115-5ab1-83bc-aa920d37bbce"
 Preferences = "21216c6a-2e73-6563-6e65-726566657250"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Tensors = "48a634ad-e948-5137-8d70-aa71f2a747f4"
 WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192"
 

docs/Manifest.toml

@@ -329,7 +329,7 @@ uuid = "29a986be-02c6-4525-aec4-84b980013641"
 version = "1.2.9"
 
 [[deps.Ferrite]]
-deps = ["EnumX", "LinearAlgebra", "NearestNeighbors", "Preferences", "Reexport", "SparseArrays", "Tensors", "WriteVTK"]
+deps = ["EnumX", "LinearAlgebra", "NearestNeighbors", "Preferences", "Reexport", "SparseArrays", "StaticArrays", "Tensors", "WriteVTK"]
 path = ".."
 uuid = "c061ca5d-56c9-439f-9c0e-210fe06d3992"
 version = "0.3.14"

docs/src/literate/landau.jl

@@ -60,7 +60,7 @@ struct ThreadCache{CV, T, DIM, F <: Function, GC <: GradientConfig, HC <: Hessia
     gradconf         ::GC
     hessconf         ::HC
 end
-function ThreadCache(dpc::Int, nodespercell, cvP::CellValues{DIM, T}, modelparams, elpotential) where {DIM, T}
+function ThreadCache(dpc::Int, nodespercell, cvP::CellValues{DIM, DIM, T}, modelparams, elpotential) where {DIM, T}
     element_indices  = zeros(Int, dpc)
     element_dofs     = zeros(dpc)
     element_gradient = zeros(dpc)

docs/src/literate/postprocessing.jl

@@ -46,7 +46,7 @@ include("heat_equation.jl");
 # Next we define a function that computes the heat flux for each integration point in the domain.
 # Fourier's law is adopted, where the conductivity tensor is assumed to be isotropic with unit
 # conductivity ``\lambda = 1 ⇒ q = - \nabla u``, where ``u`` is the temperature.
-function compute_heat_fluxes(cellvalues::CellScalarValues{dim,T}, dh::DofHandler, a) where {dim,T}
+function compute_heat_fluxes(cellvalues::CellScalarValues{dim,dim,T}, dh::DofHandler, a) where {dim,T}
 
     n = getnbasefunctions(cellvalues)
     cell_dofs = zeros(Int, n)

src/FEValues/cell_values.jl

@@ -1,7 +1,7 @@
 # Defines CellScalarValues and CellVectorValues and common methods
 """
-    CellScalarValues([::Type{T}], quad_rule::QuadratureRule, func_interpol::Interpolation, [geom_interpol::Interpolation])
-    CellVectorValues([::Type{T}], quad_rule::QuadratureRule, func_interpol::Interpolation, [geom_interpol::Interpolation])
+    CellScalarValues([::Type{T}], quad_rule::QuadratureRule, func_interpol::Interpolation, [geo_interpol::Interpolation])
+    CellVectorValues([::Type{T}], quad_rule::QuadratureRule, func_interpol::Interpolation, [geo_interpol::Interpolation])
 
 A `CellValues` object facilitates the process of evaluating values of shape functions, gradients of shape functions,
 values of nodal functions, gradients and divergences of nodal functions etc. in the finite element cell. There are
@@ -13,7 +13,7 @@ utilizes scalar shape functions and `CellVectorValues` utilizes vectorial shape
 * `T`: an optional argument (default to `Float64`) to determine the type the internal data is stored as.
 * `quad_rule`: an instance of a [`QuadratureRule`](@ref)
 * `func_interpol`: an instance of an [`Interpolation`](@ref) used to interpolate the approximated function
-* `geom_interpol`: an optional instance of a [`Interpolation`](@ref) which is used to interpolate the geometry
+* `geo_interpol`: an optional instance of a [`Interpolation`](@ref) which is used to interpolate the geometry
 
 **Common methods:**
 * [`reinit!`](@ref)
@@ -34,137 +34,222 @@ utilizes scalar shape functions and `CellVectorValues` utilizes vectorial shape
 CellValues, CellScalarValues, CellVectorValues
 
 # CellScalarValues
-struct CellScalarValues{dim,T<:Real,refshape<:AbstractRefShape} <: CellValues{dim,T,refshape}
+#   sdim = spatial dimension
+#   rdim = reference element dimension
+struct CellScalarValues{sdim,rdim,T<:Real,refshape<:AbstractRefShape} <: CellValues{sdim,rdim,T,refshape}
     N::Matrix{T}
-    dNdx::Matrix{Vec{dim,T}}
-    dNdξ::Matrix{Vec{dim,T}}
+    dNdx::Matrix{SVector{sdim,T}}
+    dNdξ::Matrix{SVector{rdim,T}}
     detJdV::Vector{T}
     M::Matrix{T}
-    dMdξ::Matrix{Vec{dim,T}}
-    qr::QuadratureRule{dim,refshape,T}
+    dMdξ::Matrix{SVector{rdim,T}}
+    qr::QuadratureRule{rdim,refshape,T}
     # The following fields are deliberately abstract -- they are never used in
     # performance critical code, just stored here for convenience.
-    func_interp::Interpolation{dim,refshape}
-    geo_interp::Interpolation{dim,refshape}
+    func_interp::Interpolation{rdim,refshape}
+    geo_interp::Interpolation{rdim,refshape}
 end
 
+# FIXME sdim should be something like `getdim(value(geo_interpol))``
 function CellScalarValues(quad_rule::QuadratureRule, func_interpol::Interpolation,
-        geom_interpol::Interpolation=func_interpol)
-    CellScalarValues(Float64, quad_rule, func_interpol, geom_interpol)
+        geo_interpol::Interpolation=func_interpol; sdim::Int=getdim(func_interpol))
+    CellScalarValues(Float64, quad_rule, func_interpol, geo_interpol, Val(sdim))
 end
 
-function CellScalarValues(::Type{T}, quad_rule::QuadratureRule{dim,shape}, func_interpol::Interpolation,
-        geom_interpol::Interpolation=func_interpol) where {dim,T,shape<:AbstractRefShape}
+# FIXME sdim should be something like `length(value(geo_interpol))`
+function CellScalarValues(valtype::Type{T}, quad_rule::QuadratureRule, func_interpol::Interpolation,
+        geo_interpol::Interpolation=func_interpol; sdim::Int=getdim(func_interpol))  where {T}
+    CellScalarValues(valtype, quad_rule, func_interpol, geo_interpol, Val(sdim))
+end
+
+function CellScalarValues(::Type{T}, quad_rule::QuadratureRule{rdim,shape}, func_interpol::Interpolation{rdim,shape},
+        geo_interpol::Interpolation{rdim,shape}, ::Val{sdim}) where {rdim,T,shape<:AbstractRefShape,sdim}
 
-    @assert getdim(func_interpol) == getdim(geom_interpol)
-    @assert getrefshape(func_interpol) == getrefshape(geom_interpol) == shape
     n_qpoints = length(getweights(quad_rule))
 
     # Function interpolation
     n_func_basefuncs = getnbasefunctions(func_interpol)
     N    = fill(zero(T)          * T(NaN), n_func_basefuncs, n_qpoints)
-    dNdx = fill(zero(Vec{dim,T}) * T(NaN), n_func_basefuncs, n_qpoints)
-    dNdξ = fill(zero(Vec{dim,T}) * T(NaN), n_func_basefuncs, n_qpoints)
+    dNdx = fill(zero(SVector{sdim,T}) * T(NaN), n_func_basefuncs, n_qpoints)
+    dNdξ = fill(zero(SVector{rdim,T}) * T(NaN), n_func_basefuncs, n_qpoints)
 
     # Geometry interpolation
-    n_geom_basefuncs = getnbasefunctions(geom_interpol)
+    n_geom_basefuncs = getnbasefunctions(geo_interpol)
     M    = fill(zero(T)          * T(NaN), n_geom_basefuncs, n_qpoints)
-    dMdξ = fill(zero(Vec{dim,T}) * T(NaN), n_geom_basefuncs, n_qpoints)
+    dMdξ = fill(zero(SVector{rdim,T}) * T(NaN), n_geom_basefuncs, n_qpoints)
 
     for (qp, ξ) in enumerate(quad_rule.points)
         for i in 1:n_func_basefuncs
             dNdξ[i, qp], N[i, qp] = gradient(ξ -> value(func_interpol, i, ξ), ξ, :all)
         end
         for i in 1:n_geom_basefuncs
-            dMdξ[i, qp], M[i, qp] = gradient(ξ -> value(geom_interpol, i, ξ), ξ, :all)
+            dMdξ[i, qp], M[i, qp] = gradient(ξ -> value(geo_interpol, i, ξ), ξ, :all)
         end
     end
 
     detJdV = fill(T(NaN), n_qpoints)
 
-    CellScalarValues{dim,T,shape}(N, dNdx, dNdξ, detJdV, M, dMdξ, quad_rule, func_interpol, geom_interpol)
+    CellScalarValues{sdim,rdim,T,shape}(N, dNdx, dNdξ, detJdV, M, dMdξ, quad_rule, func_interpol, geo_interpol)
 end
 
 # CellVectorValues
-struct CellVectorValues{dim,T<:Real,refshape<:AbstractRefShape,M} <: CellValues{dim,T,refshape}
-    N::Matrix{Vec{dim,T}}
-    dNdx::Matrix{Tensor{2,dim,T,M}}
-    dNdξ::Matrix{Tensor{2,dim,T,M}}
+#   vdim = vector dimension (i.e. dimension of evaluation of what `value` should return)
+#   sdim = spatial dimension
+#   rdim = reference element dimension
+#   M1   = number of elements in the matrix dNdx (should be vdim × sdim)
+#   M2   = number of elements in the matrix dNdξ (should be vdim × rdim)
+struct CellVectorValues{vdim,sdim,rdim,T<:Real,refshape<:AbstractRefShape,M1,M2} <: CellValues{sdim,rdim,T,refshape}
+    N::Matrix{SVector{vdim,T}} # vdim
+    dNdx::Matrix{SMatrix{vdim,sdim,T,M1}} # vdim × sdim
+    dNdξ::Matrix{SMatrix{vdim,rdim,T,M2}} # vdim × rdim
     detJdV::Vector{T}
     M::Matrix{T}
-    dMdξ::Matrix{Vec{dim,T}}
-    qr::QuadratureRule{dim,refshape,T}
+    dMdξ::Matrix{SVector{rdim,T}} # rdim
+    qr::QuadratureRule{rdim,refshape,T} #rdim
     # The following fields are deliberately abstract -- they are never used in
     # performance critical code, just stored here for convenience.
-    func_interp::Interpolation{dim,refshape}
-    geo_interp::Interpolation{dim,refshape}
+    func_interp::Interpolation{rdim,refshape} # rdim
+    geo_interp::Interpolation{rdim,refshape} # rdim
 end
 
-function CellVectorValues(quad_rule::QuadratureRule, func_interpol::Interpolation, geom_interpol::Interpolation=func_interpol)
-    CellVectorValues(Float64, quad_rule, func_interpol, geom_interpol)
+# FIXME sdim should be something like `length(value(geo_interpol))`
+# FIXME vdim should be something like `length(value(func_interpol))`
+function CellVectorValues(quad_rule::QuadratureRule, func_interpol::Interpolation, 
+        geo_interpol::Interpolation=func_interpol; sdim::Int=getdim(geo_interpol), vdim::Int=getdim(func_interpol))
+    CellVectorValues(Float64, quad_rule, func_interpol, geo_interpol, Val(sdim), Val(vdim))
 end
 
-function CellVectorValues(::Type{T}, quad_rule::QuadratureRule{dim,shape}, func_interpol::Interpolation,
-        geom_interpol::Interpolation=func_interpol) where {dim,T,shape<:AbstractRefShape}
+# FIXME sdim should be something like `length(value(geo_interpol))`
+# FIXME vdim should be something like `length(value(func_interpol))`
+function CellVectorValues(valuetype::Type{T}, quad_rule::QuadratureRule, func_interpol::Interpolation, 
+        geo_interpol::Interpolation=func_interpol; sdim::Int=getdim(geo_interpol), vdim::Int=getdim(func_interpol)) where {T}
+    CellVectorValues(valuetype, quad_rule, func_interpol, geo_interpol, Val(sdim), Val(vdim))
+end
 
-    @assert getdim(func_interpol) == getdim(geom_interpol)
-    @assert getrefshape(func_interpol) == getrefshape(geom_interpol) == shape
+function CellVectorValues(::Type{T}, quad_rule::QuadratureRule{rdim,shape}, func_interpol::Interpolation,
+        geo_interpol::Interpolation, ::Val{sdim}, ::Val{vdim}) where {rdim,T,shape<:AbstractRefShape,sdim,vdim}
+    @assert getrefshape(func_interpol) == getrefshape(geo_interpol) == shape
     n_qpoints = length(getweights(quad_rule))
 
+    M1 = vdim*sdim
+    M2 = vdim*rdim
+
     # Function interpolation
-    n_func_basefuncs = getnbasefunctions(func_interpol) * dim
-    N    = fill(zero(Vec{dim,T})      * T(NaN), n_func_basefuncs, n_qpoints)
-    dNdx = fill(zero(Tensor{2,dim,T}) * T(NaN), n_func_basefuncs, n_qpoints)
-    dNdξ = fill(zero(Tensor{2,dim,T}) * T(NaN), n_func_basefuncs, n_qpoints)
+    n_func_basefuncs = getnbasefunctions(func_interpol) * vdim
+    N    = fill(zero(SVector{vdim,T})      * T(NaN), n_func_basefuncs, n_qpoints)
+    dNdx = fill(zero(SMatrix{vdim,sdim,T,M1}) * T(NaN), n_func_basefuncs, n_qpoints)
+    dNdξ = fill(zero(SMatrix{vdim,rdim,T,M2}) * T(NaN), n_func_basefuncs, n_qpoints)
 
     # Geometry interpolation
-    n_geom_basefuncs = getnbasefunctions(geom_interpol)
+    n_geom_basefuncs = getnbasefunctions(geo_interpol)
     M    = fill(zero(T)          * T(NaN), n_geom_basefuncs, n_qpoints)
-    dMdξ = fill(zero(Vec{dim,T}) * T(NaN), n_geom_basefuncs, n_qpoints)
+    dMdξ = fill(zero(SVector{rdim,T}) * T(NaN), n_geom_basefuncs, n_qpoints)
 
     for (qp, ξ) in enumerate(quad_rule.points)
         basefunc_count = 1
         for basefunc in 1:getnbasefunctions(func_interpol)
             dNdξ_temp, N_temp = gradient(ξ -> value(func_interpol, basefunc, ξ), ξ, :all)
-            for comp in 1:dim
-                N_comp = zeros(T, dim)
+            for comp in 1:vdim
+                N_comp = zeros(T, vdim)
                 N_comp[comp] = N_temp
-                N[basefunc_count, qp] = Vec{dim,T}((N_comp...,))
+                N[basefunc_count, qp] = SVector{vdim,T}((N_comp...,))
 
-                dN_comp = zeros(T, dim, dim)
+                dN_comp = zeros(T, vdim, rdim)
                 dN_comp[comp, :] = dNdξ_temp
-                dNdξ[basefunc_count, qp] = Tensor{2,dim,T}((dN_comp...,))
+                dNdξ[basefunc_count, qp] = SMatrix{vdim,rdim,T}((dN_comp...,))
                 basefunc_count += 1
             end
         end
         for basefunc in 1:n_geom_basefuncs
-            dMdξ[basefunc, qp], M[basefunc, qp] = gradient(ξ -> value(geom_interpol, basefunc, ξ), ξ, :all)
+            dMdξ[basefunc, qp], M[basefunc, qp] = gradient(ξ -> value(geo_interpol, basefunc, ξ), ξ, :all)
         end
     end
 
     detJdV = fill(T(NaN), n_qpoints)
-    MM = Tensors.n_components(Tensors.get_base(eltype(dNdx)))
 
-    CellVectorValues{dim,T,shape,MM}(N, dNdx, dNdξ, detJdV, M, dMdξ, quad_rule, func_interpol, geom_interpol)
+    CellVectorValues{vdim,sdim,rdim,T,shape,M1,M2}(N, dNdx, dNdξ, detJdV, M, dMdξ, quad_rule, func_interpol, geo_interpol)
 end
 
-function reinit!(cv::CellValues{dim}, x::AbstractVector{Vec{dim,T}}) where {dim,T}
+"""
+Reinit for volumetric elements, i.e. elements whose reference dimension match the spatial dimension.
+"""
+function reinit!(cv::CellValues{sdim,sdim}, x::AbstractVector{Vec{sdim,T}}) where {sdim,T}
     n_geom_basefuncs = getngeobasefunctions(cv)
     n_func_basefuncs = getnbasefunctions(cv)
     length(x) == n_geom_basefuncs || throw_incompatible_coord_length(length(x), n_geom_basefuncs)
+    TdNdx = typeof(cv.dNdx[1, 1])
 
     @inbounds for i in 1:length(cv.qr.weights)
         w = cv.qr.weights[i]
-        fecv_J = zero(Tensor{2,dim})
+        fecv_J = zero(Tensor{2,sdim,T})
         for j in 1:n_geom_basefuncs
-            fecv_J += x[j] ⊗ cv.dMdξ[j, i]
+            fecv_J += x[j] ⊗ tensor_cast(cv.dMdξ[j, i]) # TODO cleanup?
         end
         detJ = det(fecv_J)
         detJ > 0.0 || throw_detJ_not_pos(detJ)
         cv.detJdV[i] = detJ * w
         Jinv = inv(fecv_J)
         for j in 1:n_func_basefuncs
-            cv.dNdx[j, i] = cv.dNdξ[j, i] ⋅ Jinv
+            cv.dNdx[j, i] = TdNdx((tensor_cast(cv.dNdξ[j, i]) ⋅ Jinv).data) # TODO cleanup?
+        end
+    end
+end
+
+# Hotfix to get the dots right for embedded elements.
+@inline dothelper(x::V,A::M) where {V<:SVector,M<:Union{SMatrix,MMatrix}} = A'*x
+@inline dothelper(B::M1,A::M2) where {M1<:SMatrix,M2<:Union{SMatrix,MMatrix}} = B*A
+
+"""
+Embedding determinant for surfaces in 3D.
+
+TLDR: "det(J) =" ||∂x/∂ξ₁ × ∂x/∂ξ₂||₂
+
+The transformation theorem for some function f on a 2D surface in 3D space leads to
+  ∫ f ⋅ dS = ∫ f ⋅ (∂x/∂ξ₁ × ∂x/∂ξ₂) dξ₁dξ₂ = ∫ f ⋅ n ||∂x/∂ξ₁ × ∂x/∂ξ₂||₂ dξ₁dξ₂
+where ||∂x/∂ξ₁ × ∂x/∂ξ₂||₂ is "detJ" and n is the unit normal.
+See e.g. https://scicomp.stackexchange.com/questions/41741/integration-of-d-1-dimensional-functions-on-finite-element-surfaces for simple explanation.
+For more details see e.g. the doctoral thesis by Mirza Cenanovic **Finite element methods for surface problems* (2017), Ch. 2 **Trangential Calculus**.
+"""
+edet(J::MMatrix{3,2,T}) where {T} = norm(J[:,1] × J[:,2])
+
+"""
+Embedding determinant for curves in 2D and 3D.
+
+TLDR: "det(J) =" ||∂x/∂ξ||₂
+
+The transformation theorem for some function f on a 1D curve in 2D and 3D space leads to
+  ∫ f ⋅ dE = ∫ f ⋅ ∂x/∂ξ dξ = ∫ f ⋅ t ||∂x/∂ξ||₂ dξ
+where ||∂x/∂ξ||₂ is "detJ" and t is "the unit tangent".
+See e.g. https://scicomp.stackexchange.com/questions/41741/integration-of-d-1-dimensional-functions-on-finite-element-surfaces for simple explanation.
+"""
+edet(J::Union{MMatrix{2,1,T},MMatrix{3,1,T}}) where {T} = norm(J)
+
+"""
+Reinit for embedded elements, i.e. elements whose reference dimension is smaller than the spatial dimension.
+"""
+function reinit!(cv::CellValues{sdim,rdim}, x::AbstractVector{Vec{sdim,T}}) where {sdim,rdim,T}
+    @assert sdim > rdim "This reinit only works for embedded elements. Maybe you swapped the reference and spatial dimensions?"
+    n_geom_basefuncs = getngeobasefunctions(cv)
+    n_func_basefuncs = getnbasefunctions(cv)
+    length(x) == n_geom_basefuncs || throw_incompatible_coord_length(length(x), n_geom_basefuncs)
+
+    @inbounds for i in 1:length(cv.qr.weights)
+        w = cv.qr.weights[i]
+        fecv_J = zero(MMatrix{sdim,rdim,T}) # TODO replace with MixedTensor (see https://github.com/Ferrite-FEM/Tensors.jl/pull/188)
+        for j in 1:n_geom_basefuncs
+            #fecv_J += x[j] ⊗ cv.dMdξ[j, i] # TODO via Tensors.jl
+            for k in 1:sdim, l in 1:rdim
+                fecv_J[k,l] += x[j][k] * cv.dMdξ[j, i][l]
+            end
+        end
+        detJ = edet(fecv_J)
+        detJ > 0.0 || throw_detJ_not_pos(detJ)
+        cv.detJdV[i] = detJ * w
+        # Compute "left inverse" of J
+        Jinv = pinv(fecv_J)
+        for j in 1:n_func_basefuncs
+            #cv.dNdx[j, i] = cv.dNdξ[j, i] ⋅ Jinv # TODO via Tensors.jl
+            cv.dNdx[j, i] = dothelper(cv.dNdξ[j, i], Jinv)
         end
     end
 end