WIP: New vector interpolations and mappings

benchmark/helper.jl

@@ -66,17 +66,17 @@ function _generalized_ritz_galerkin_assemble_local_matrix(grid::Ferrite.Abstract
 end
 
 # Minimal Petrov-Galerkin type local assembly loop. We assume that both function spaces share the same integration rule. Test is applied from the left.
-function _generalized_petrov_galerkin_assemble_local_matrix(grid::Ferrite.AbstractGrid, cellvalues_shape::CellValues{<: Ferrite.InterpolationByDim{dim}}, f_shape, cellvalues_test::CellValues{<: Ferrite.InterpolationByDim{dim}}, f_test, op) where {dim}
+function _generalized_petrov_galerkin_assemble_local_matrix(grid::Ferrite.AbstractGrid, cellvalues_shape::CellValues, f_shape, cellvalues_test::CellValues, f_test, op)
     n_basefuncs_shape = getnbasefunctions(cellvalues_shape)
     n_basefuncs_test = getnbasefunctions(cellvalues_test)
     Ke = zeros(n_basefuncs_test, n_basefuncs_shape)
 
     #implicit assumption: Same geometry!
-    X_shape = zeros(Vec{dim,Float64}, Ferrite.getngeobasefunctions(cellvalues_shape))
+    X_shape = zeros(get_coordinate_type(grid), Ferrite.getngeobasefunctions(cellvalues_shape))
     getcoordinates!(X_shape, grid, 1)
     reinit!(cellvalues_shape, X_shape)
 
-    X_test = zeros(Vec{dim,Float64}, Ferrite.getngeobasefunctions(cellvalues_test))
+    X_test = zeros(get_coordinate_type(grid), Ferrite.getngeobasefunctions(cellvalues_test))
     getcoordinates!(X_test, grid, 1)
     reinit!(cellvalues_test, X_test)
 

docs/Project.toml

@@ -16,6 +16,7 @@ Optim = "429524aa-4258-5aef-a3af-852621145aeb"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Tensors = "48a634ad-e948-5137-8d70-aa71f2a747f4"
 TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
 UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"

docs/make.jl

@@ -47,6 +47,8 @@ bibtex_plugin = CitationBibliography(
         "Tutorials" => [
             "Tutorials overview" => "tutorials/index.md",
             "tutorials/heat_equation.md",
+            "tutorials/heat_equation_rt.md",
+            "tutorials/heat_equation_triangle.md",
             "tutorials/linear_elasticity.md",
             "tutorials/incompressible_elasticity.md",
             "tutorials/hyperelasticity.md",
@@ -61,6 +63,7 @@ bibtex_plugin = CitationBibliography(
         "Topic guides" => [
             "Topic guide overview" => "topics/index.md",
             "topics/fe_intro.md",
+            "topics/FEValues.md",
             "topics/degrees_of_freedom.md",
             "topics/assembly.md",
             "topics/boundary_conditions.md",

docs/src/assets/references.bib

@@ -60,6 +60,15 @@ @article{Scroggs2022
 @phdthesis{Cenanovic2017,
   author={Mirza Cenanovic},
   title={Finite element methods for surface problems},
-  school={J\o{}pk\o{}ping University, School of Engineering},
+  school={Jönköping University, School of Engineering},
   year=2017
 }
+@misc{Kirby2017,
+      title={A general approach to transforming finite elements}, 
+      author={Robert C. Kirby},
+      year={2017},
+      eprint={1706.09017},
+      doi={10.48550/arXiv.1706.09017},
+      archivePrefix={arXiv},
+      primaryClass={math.NA}
+}

docs/src/devdocs/FEValues.md

@@ -0,0 +1,42 @@
+# [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
+```
+
+## How CellValues and FaceValues works
+* The function interpolation, `ip_fun`, determines how it should be mapped, by defining `get_mapping_type(ip_fun)` for its type.
+* The mapping type, e.g. `IdentityMapping`, decides the requirements for `GeometryValues`, specifically if the `hessian` $\partial^2M/\partial\xi^2$,
+  of the geometric shape functions, $M(\xi)$, on the reference cell should be precalculated or not. 
+  ***Note:*** *This should also in the future be determined by the required order of derivatives to be mapped in `FunctionValues`*
+* As the first part of `reinit!`, the `MappingValues` are calculated based on the cell's coordinates. If the `GeometricMapping` contains the hessian 
+  on the reference cell, the `hessian` on the actual cell, $\partial^2M/\partial x^2$, is calculated and returned in `MappingValues`. Otherwise, only
+  the jacobian, $\partial M/\partial x$, is calculated. 
+* In the second part of `reinit!`, The `MappingValues`, containing sufficient information for the current quadrature point, is given to `apply_mapping!`. 
+  This allows the shape values and gradients stored in `FunctionValues`, to be mapped to the current cell by calling `apply_mapping!`.
+
+## Custom FEValues
+Custom FEValues, `fe_v`, should normally implement the `reinit!` method.
+Additionally, for normal functionality the `getnquadpoints` should be implemented.
+Note that 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.)
+
+Supporting `function_value`, `function_gradient`, `function_symmetric_gradient`, `function_divergence`, and `function_curl`,
+requires implementing `getnbasefunctions`, `shape_value`, and `shape_gradient`. 
+Note that asking for the `i`th shape value or gradient must be inside array bounds if `1<=i<:getnbasefunctions(fe_v)`
+
+Supporting `spatial_coordinate` requires implementing `getngeobasefunctions` and `geometric_value`.
+Note that asking for the `i`th geometric value must be inside array bounds if `1<=i<:getngeobasefunctions(fe_v)`

docs/src/devdocs/index.md

@@ -5,5 +5,5 @@ developing the library.
 
 ```@contents
 Depth = 1
-Pages = ["reference_cells.md", "interpolations.md", "elements.md", "dofhandler.md", "performance.md"]
+Pages = ["reference_cells.md", "interpolations.md", "elements.md", "FEValues.md", "dofhandler.md", "performance.md"]
 ```

docs/src/devdocs/interpolations.md

@@ -40,6 +40,7 @@ Ferrite.reference_coordinates(::Interpolation)
 Ferrite.face_to_element_transformation
 Ferrite.is_discontinuous(::Interpolation)
 Ferrite.adjust_dofs_during_distribution(::Interpolation)
+Ferrite.get_mapping_type
 ```
 
 for all entities which exist on that reference element. The dof functions default to having no

docs/src/literate-howto/postprocessing.jl

@@ -46,7 +46,7 @@ include("../tutorials/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::CellValues{<:ScalarInterpolation}, dh::DofHandler, a::AbstractVector{T}) where T
+function compute_heat_fluxes(cellvalues::CellValues, dh::DofHandler, a::AbstractVector{T}) where T
 
     n = getnbasefunctions(cellvalues)
     cell_dofs = zeros(Int, n)