forked from Ferrite-FEM/Ferrite.jl
-
Notifications
You must be signed in to change notification settings - Fork 0
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
GeometryMapping and FunctionValues as sub-components for FEValues (Fe…
…rrite-FEM#764) This patch should be non-breaking for documented usage of FEValues. The main changes are: * Introduce FunctionValues and GeometryMapping as sub-components for use in FEValues * Use FunctionValues and GeometryMapping in CellValues and FaceValues * FaceValues have a vector of these types, one for each face. * Prepare for generalized mapping by having the cell::AbstractCell as (optional) input to reinit * Allow keyword arguments to CellValues, FaceValues, and PointValues to optionally update the gradient * Replace PointValuesInternal with a PointValues that doesn't update jacobians and gradients * Update both geometric and function shape values and gradients during changes of quadrature points in reinit of PointValues Documentation changes: * Describe the theory behind the mapping of finite elements * Show how a SimpleCellValues object can be created to explain via code how this mapping is done in a simple case * Devdocs for FEValues
- Loading branch information
Showing
30 changed files
with
1,194 additions
and
636 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Large diffs are not rendered by default.
Oops, something went wrong.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,39 @@ | ||
| # [FEValues](@id devdocs-fevalues) | ||
|
|
||
| ## Type definitions | ||
| * `AbstractValues` | ||
| * `AbstractCellValues` | ||
| * [`CellValues`](@ref) | ||
| * `AbstractFaceValues` | ||
| * [`FaceValues`](@ref) | ||
| * [`PointValues`](@ref) | ||
| * `PointValuesInternal` (Optimized version of PointValues) | ||
|
|
||
| ## Internal types | ||
| ```@docs | ||
| Ferrite.GeometryMapping | ||
| Ferrite.MappingValues | ||
| Ferrite.FunctionValues | ||
| ``` | ||
|
|
||
| ## Custom FEValues | ||
| Custom FEValues, `fe_v::AbstractValues`, should normally implement the [`reinit!`](@ref) method. Subtypes of `AbstractValues` have default implementations for some functions, but require some lower-level access functions, specifically | ||
|
|
||
| * [`function_value`](@ref), requires | ||
| * [`shape_value`](@ref) | ||
| * [`getnquadpoints`](@ref) | ||
| * [`getnbasefunctions`](@ref) | ||
| * [`function_gradient`](@ref), [`function_divergence`](@ref), [`function_symmetric_gradient`](@ref), and [`function_curl`](@ref) requires | ||
| * [`shape_gradient`](@ref) | ||
| * [`getnquadpoints`](@ref) | ||
| * [`getnbasefunctions`](@ref) | ||
| * [`spatial_coordinate`](@ref), requires | ||
| * [`geometric_value`](@ref) | ||
| * [`getngeobasefunctions`](@ref) | ||
| * [`getnquadpoints`](@ref) | ||
|
|
||
|
|
||
| ### Array bounds | ||
| * Asking for the `n`th quadrature point must be inside array bounds if `1 <= n <= getnquadpoints(fe_v)`. (`checkquadpoint` can, alternatively, be dispatched to check that `n` is inbounds.) | ||
| * Asking for the `i`th shape value or gradient must be inside array bounds if `1 <= i <= getnbasefunctions(fe_v)` | ||
| * Asking for the `i`th geometric value must be inside array bounds if `1 <= i <= getngeobasefunctions(fe_v)` |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,127 @@ | ||
| # [FEValues](@id fevalues_topicguide) | ||
| A key type of object in `Ferrite.jl` is the so-called `FEValues`, where the most common ones are `CellValues` and `FaceValues`. These objects are used inside the element routines and are used to query the integration weights, shape function values and gradients, and much more; see [`CellValues`](@ref) and [`FaceValues`](@ref). For these values to be correct, it is necessary to reinitialize these for the current cell by using the [`reinit!`](@ref) function. This function maps the values from the reference cell to the actual cell, a process described in detail below, see [Mapping of finite elements](@ref mapping-theory). After that, we show an implementation of a [`SimpleCellValues`](@ref SimpleCellValues) type to illustrate how `CellValues` work for the most standard case, excluding the generalizations and optimization that complicates the actual code. | ||
|
|
||
| ## [Mapping of finite elements](@id mapping_theory) | ||
| The shape functions and gradients stored in an `FEValues` object, are reinitialized for each cell by calling the `reinit!` function. | ||
| The main part of this calculation, considers how to map the values and derivatives of the shape functions, | ||
| defined on the reference cell, to the actual cell. | ||
|
|
||
| The geometric mapping of a finite element from the reference coordinates to the real coordinates is shown in the following illustration. | ||
|
|
||
|  | ||
|
|
||
| This mapping is given by the geometric shape functions, $\hat{N}_i^g(\boldsymbol{\xi})$, such that | ||
| ```math | ||
| \begin{align*} | ||
| \boldsymbol{x}(\boldsymbol{\xi}) =& \sum_{\alpha=1}^N \hat{\boldsymbol{x}}_\alpha \hat{N}_\alpha^g(\boldsymbol{\xi}) \\ | ||
| \boldsymbol{J} :=& \frac{\mathrm{d}\boldsymbol{x}}{\mathrm{d}\boldsymbol{\xi}} = \sum_{\alpha=1}^N \hat{\boldsymbol{x}}_\alpha \otimes \frac{\mathrm{d} \hat{N}_\alpha^g}{\mathrm{d}\boldsymbol{\xi}}\\ | ||
| \boldsymbol{\mathcal{H}} :=& | ||
| \frac{\mathrm{d} \boldsymbol{J}}{\mathrm{d} \boldsymbol{\xi}} = \sum_{\alpha=1}^N \hat{\boldsymbol{x}}_\alpha \otimes \frac{\mathrm{d}^2 \hat{N}^g_\alpha}{\mathrm{d} \boldsymbol{\xi}^2} | ||
| \end{align*} | ||
| ``` | ||
| where the defined $\boldsymbol{J}$ is the jacobian of the mapping, and in some cases we will also need the corresponding hessian, $\boldsymbol{\mathcal{H}}$ (3rd order tensor). | ||
|
|
||
| We require that the mapping from reference coordinates to real coordinates is [diffeomorphic](https://en.wikipedia.org/wiki/Diffeomorphism), meaning that we can express $\boldsymbol{x} = \boldsymbol{x}(\boldsymbol{\xi}(\boldsymbol{x}))$, such that | ||
| ```math | ||
| \begin{align*} | ||
| \frac{\mathrm{d}\boldsymbol{x}}{\mathrm{d}\boldsymbol{x}} = \boldsymbol{I} &= \frac{\mathrm{d}\boldsymbol{x}}{\mathrm{d}\boldsymbol{\xi}} \cdot \frac{\mathrm{d}\boldsymbol{\xi}}{\mathrm{d}\boldsymbol{x}} | ||
| \quad\Rightarrow\quad | ||
| \frac{\mathrm{d}\boldsymbol{\xi}}{\mathrm{d}\boldsymbol{x}} = \left[\frac{\mathrm{d}\boldsymbol{x}}{\mathrm{d}\boldsymbol{\xi}}\right]^{-1} = \boldsymbol{J}^{-1} | ||
| \end{align*} | ||
| ``` | ||
| Depending on the function interpolation, we may want different types of mappings to conserve certain properties of the fields. This results in the different mapping types described below. | ||
|
|
||
| ### Identity mapping | ||
| `Ferrite.IdentityMapping` | ||
|
|
||
| For scalar fields, we always use scalar base functions. For tensorial fields (non-scalar, e.g. vector-fields), the base functions can be constructed from scalar base functions, by using e.g. `VectorizedInterpolation`. From the perspective of the mapping, however, each component is mapped as an individual scalar base function. And for scalar base functions, we only require that the value of the base function is invariant to the element shape (real coordinate), and only depends on the reference coordinate, i.e. | ||
| ```math | ||
| \begin{align*} | ||
| N(\boldsymbol{x}) &= \hat{N}(\boldsymbol{\xi}(\boldsymbol{x}))\nonumber \\ | ||
| \mathrm{grad}(N(\boldsymbol{x})) &= \frac{\mathrm{d}\hat{N}}{\mathrm{d}\boldsymbol{\xi}} \cdot \boldsymbol{J}^{-1} | ||
| \end{align*} | ||
| ``` | ||
|
|
||
| ### Covariant Piola mapping, H(curl) | ||
| `Ferrite.CovariantPiolaMapping` | ||
|
|
||
| The covariant Piola mapping of a vectorial base function preserves the tangential components. For the value, the mapping is defined as | ||
| ```math | ||
| \begin{align*} | ||
| \boldsymbol{N}(\boldsymbol{x}) = \boldsymbol{J}^{-\mathrm{T}} \cdot \hat{\boldsymbol{N}}(\boldsymbol{\xi}(\boldsymbol{x})) | ||
| \end{align*} | ||
| ``` | ||
| which yields the gradient, | ||
| ```math | ||
| \begin{align*} | ||
| \mathrm{grad}(\boldsymbol{N}(\boldsymbol{x})) &= \boldsymbol{J}^{-T} \cdot \frac{\mathrm{d} \hat{\boldsymbol{N}}}{\mathrm{d} \boldsymbol{\xi}} \cdot \boldsymbol{J}^{-1} - \boldsymbol{J}^{-T} \cdot \left[\hat{\boldsymbol{N}}(\boldsymbol{\xi}(\boldsymbol{x}))\cdot \boldsymbol{J}^{-1} \cdot \boldsymbol{\mathcal{H}}\cdot \boldsymbol{J}^{-1}\right] | ||
| \end{align*} | ||
| ``` | ||
|
|
||
| !!! details "Derivation" | ||
| Expressing the gradient, $\mathrm{grad}(\boldsymbol{N})$, in index notation, | ||
| ```math | ||
| \begin{align*} | ||
| \frac{\mathrm{d} N_i}{\mathrm{d} x_j} &= \frac{\mathrm{d}}{\mathrm{d} x_j} \left[J^{-\mathrm{T}}_{ik} \hat{N}_k\right] = \frac{\mathrm{d} J^{-\mathrm{T}}_{ik}}{\mathrm{d} x_j} \hat{N}_k + J^{-\mathrm{T}}_{ik} \frac{\mathrm{d} \hat{N}_k}{\mathrm{d} \xi_l} J_{lj}^{-1} | ||
| \end{align*} | ||
| ``` | ||
|
|
||
| Except for a few elements, $\boldsymbol{J}$ varies as a function of $\boldsymbol{x}$. The derivative can be calculated as | ||
| ```math | ||
| \begin{align*} | ||
| \frac{\mathrm{d} J^{-\mathrm{T}}_{ik}}{\mathrm{d} x_j} &= \frac{\mathrm{d} J^{-\mathrm{T}}_{ik}}{\mathrm{d} J_{mn}} \frac{\mathrm{d} J_{mn}}{\mathrm{d} x_j} = - J_{km}^{-1} J_{in}^{-T} \frac{\mathrm{d} J_{mn}}{\mathrm{d} x_j} \nonumber \\ | ||
| \frac{\mathrm{d} J_{mn}}{\mathrm{d} x_j} &= \mathcal{H}_{mno} J_{oj}^{-1} | ||
| \end{align*} | ||
| ``` | ||
|
|
||
| ### Contravariant Piola mapping, H(div) | ||
| `Ferrite.ContravariantPiolaMapping` | ||
|
|
||
| The covariant Piola mapping of a vectorial base function preserves the normal components. For the value, the mapping is defined as | ||
| ```math | ||
| \begin{align*} | ||
| \boldsymbol{N}(\boldsymbol{x}) = \frac{\boldsymbol{J}}{\det(\boldsymbol{J})} \cdot \hat{\boldsymbol{N}}(\boldsymbol{\xi}(\boldsymbol{x})) | ||
| \end{align*} | ||
| ``` | ||
| This gives the gradient | ||
| ```math | ||
| \begin{align*} | ||
| \mathrm{grad}(\boldsymbol{N}(\boldsymbol{x})) = [\boldsymbol{\mathcal{H}}\cdot\boldsymbol{J}^{-1}] : \frac{[\boldsymbol{I} \underline{\otimes} \boldsymbol{I}] \cdot \hat{\boldsymbol{N}}}{\det(\boldsymbol{J})} | ||
| - \left[\frac{\boldsymbol{J} \cdot \hat{\boldsymbol{N}}}{\det(\boldsymbol{J})}\right] \otimes \left[\boldsymbol{J}^{-T} : \boldsymbol{\mathcal{H}} \cdot \boldsymbol{J}^{-1}\right] | ||
| + \boldsymbol{J} \cdot \frac{\mathrm{d} \hat{\boldsymbol{N}}}{\mathrm{d} \boldsymbol{\xi}} \cdot \frac{\boldsymbol{J}^{-1}}{\det(\boldsymbol{J})} | ||
| \end{align*} | ||
| ``` | ||
|
|
||
| !!! details "Derivation" | ||
| Expressing the gradient, $\mathrm{grad}(\boldsymbol{N})$, in index notation, | ||
| ```math | ||
| \begin{align*} | ||
| \frac{\mathrm{d} N_i}{\mathrm{d} x_j} &= \frac{\mathrm{d}}{\mathrm{d} x_j} \left[\frac{J_{ik}}{\det(\boldsymbol{J})} \hat{N}_k\right] =\nonumber\\ | ||
| &= \frac{\mathrm{d} J_{ik}}{\mathrm{d} x_j} \frac{\hat{N}_k}{\det(\boldsymbol{J})} | ||
| - \frac{\mathrm{d} \det(\boldsymbol{J})}{\mathrm{d} x_j} \frac{J_{ik} \hat{N}_k}{\det(\boldsymbol{J})^2} | ||
| + \frac{J_{ik}}{\det(\boldsymbol{J})} \frac{\mathrm{d} \hat{N}_k}{\mathrm{d} \xi_l} J_{lj}^{-1} \\ | ||
| &= \mathcal{H}_{ikl} J^{-1}_{lj} \frac{\hat{N}_k}{\det(\boldsymbol{J})} | ||
| - J^{-T}_{mn} \mathcal{H}_{mnl} J^{-1}_{lj} \frac{J_{ik} \hat{N}_k}{\det(\boldsymbol{J})} | ||
| + \frac{J_{ik}}{\det(\boldsymbol{J})} \frac{\mathrm{d} \hat{N}_k}{\mathrm{d} \xi_l} J_{lj}^{-1} | ||
| \end{align*} | ||
| ``` | ||
|
|
||
| ## [Walkthrough: Creating `SimpleCellValues`](@id SimpleCellValues) | ||
| In the following, we walk through how to create a `SimpleCellValues` type which | ||
| works similar to `Ferrite.jl`'s `CellValues`, but is not performance optimized and not as general. The main purpose is to explain how the `CellValues` works for the standard case of `IdentityMapping` described above. | ||
| Please note that several internal functions are used, and these may change without a major version increment. Please see the [Developer documentation](@ref) for their documentation. | ||
|
|
||
| ```@eval | ||
| # Include the example here, but modify the Literate output to suit being embedded | ||
| using Literate, Markdown | ||
| base_name = "SimpleCellValues_literate" | ||
| Literate.markdown(string(base_name, ".jl"); name = base_name, execute = true, credit = false, documenter=false) | ||
| content = read(string(base_name, ".md"), String) | ||
| rm(string(base_name, ".md")) | ||
| rm(string(base_name, ".jl")) | ||
| Markdown.parse(content) | ||
| ``` | ||
|
|
||
| ## Further reading | ||
| * [defelement.com](https://defelement.com/ciarlet.html#Mapping+finite+elements) | ||
| * Kirby (2017) [Kirby2017](@cite) |
Oops, something went wrong.