diff --git a/NEWS.md b/NEWS.md index 524a03951..a2c41f59b 100644 --- a/NEWS.md +++ b/NEWS.md @@ -11,6 +11,11 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - Added AMR-related methods `mark` and `estimate` to `Adaptivity` module. Implemented Dorfler marking strategy. Since PR[#1063](https://github.com/gridap/Gridap.jl/pull/1063). +### 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 ### Added @@ -25,6 +30,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - Fixed issue where barycentric refinement rule in 3D would not produce oriented meshes. Since PR[#1055](https://github.com/gridap/Gridap.jl/pull/1055). ### Changed + - Optimized MonomialBasis low-level functions. Since PR[#1059](https://github.com/gridap/Gridap.jl/pull/1059). ## [0.18.7] - 2024-10-8 diff --git a/docs/make.jl b/docs/make.jl index a2e8cf3ec..b05b679a5 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -23,6 +23,7 @@ pages = [ "Gridap.Adaptivity" => "Adaptivity.md", "Developper notes" => Any[ "dev-notes/block-assemblers.md", + "dev-notes/pullbacks.md", ], ] diff --git a/docs/src/Adaptivity.md b/docs/src/Adaptivity.md index 5a698528a..04e0b3d62 100644 --- a/docs/src/Adaptivity.md +++ b/docs/src/Adaptivity.md @@ -156,6 +156,4 @@ When a cell is refined, we need to be able to evaluate the fields defined on the ```@docs FineToCoarseField -FineToCoarseDofBasis -FineToCoarseRefFE ``` diff --git a/docs/src/ODEs.md b/docs/src/ODEs.md index 3701457b9..abd79dd57 100644 --- a/docs/src/ODEs.md +++ b/docs/src/ODEs.md @@ -44,10 +44,12 @@ More precisely, we would like ``\boldsymbol{u}_{n}`` to be close to ``\boldsymbo ``` This is a condition on the design of the pair (``\mathcal{I}``, ``\mathcal{F}``). -# Classification of ODEs and numerical schemes +## Classification of ODEs and numerical schemes + Essentially, a numerical scheme converts a (continuous) ODE into (discrete) nonlinear systems of equations. These systems of equations can be linear under special conditions on the nature of the ODE and the numerical scheme. Since numerical methods for linear and nonlinear systems of equations can be quite different in terms of cost and implementation, we are interested in solving linear systems whenever possible. This leads us to perform the following classifications. -## Classification of ODEs +### Classification of ODEs + We define a few nonlinearity types based on the expression of the residual. * **Nonlinear**. Nothing special can be said about the residual. * **Quasilinear**. The residual is linear with respect to the highest-order time derivative and the corresponding linear form may depend on time and lower-order time derivatives, i.e. @@ -83,7 +85,8 @@ In particular, for the global residual to be linear, both the implicit and expli > In the special case where the implicit part is linear and the explicit part is quasilinear or semilinear, we could, in theory, identify two linear forms for the global residual. However, introducing this difference would call for an order-dependent classification of ODEs and this would create (infinitely) many new types. Since numerical schemes can rarely take advantage of this extra structure in practice, we still say that the global residual is semilinear in these cases. -## Classification of numerical schemes +### Classification of numerical schemes + We introduce a classification of numerical schemes based on where they evaluate the residual during the state update. * If it is possible (up to a change of variables) to write the system of equations for the state update as evaluations of the residual at known values (that depend on the solution at the current time) for all but the highest-order derivative, we say that the scheme is explicit. @@ -95,7 +98,8 @@ We introduce a classification of numerical schemes based on where they evaluate > ``` > where ``\boldsymbol{x}`` and the unknown of the state update. The scheme is explicit if it is possible to introduce a change of variables such that ``\boldsymbol{u}_{k}`` does not depend on ``\boldsymbol{x}``. Otherwise, it is implicit. -## Classification of systems of equations +### Classification of systems of equations + It is advantageous to introduce this classification of ODE and numerical schemes because the system of equations arising from the discretisation of the ODE by a numerical scheme will be linear or nonlinear depending on whether the scheme is explicit, implicit, or implicit-explicit, and on the type of the ODE. More precisely, we have the following table. | | Nonlinear | Quasilinear | Semilinear | Linear | @@ -107,7 +111,8 @@ When the system is linear, another important practical consideration is whether * If the linear system comes from an explicit scheme, the matrix of the system is constant if the mass matrix is. This means that the ODE has to be quasilinear. * If the linear system comes from an implicit scheme, all the linear forms must be constant for the system to have a constant matrix. -## Reuse across iterations +### Reuse across iterations + For performance reasons, it is thus important that the ODE be described in the most specific way. In particular, we consider that the mass term of a quasilinear ODE is not constant, because if it is, the ODE is semilinear. We enable the user to specify the following constant annotations: * For nonlinear and quasilinear ODE, no quantity can be described as constant. * For a semilinear ODE, whether the mass term is constant. @@ -115,10 +120,12 @@ For performance reasons, it is thus important that the ODE be described in the m If a linear form is constant, regardless of whether the numerical scheme relies on a linear or nonlinear system, it is always possible to compute the jacobian of the residual with respect to the corresponding time derivative only once and retrieve it in subsequent computations of the jacobian. -# High-level API in Gridap +## High-level API in Gridap + The ODE module of `Gridap` relies on the following structure. -## Finite element spaces +### Finite element spaces + The time-dependent counterpart of `TrialFESpace` is `TransientTrialFESpace`. It is built from a standard `TestFESpace` and is equipped with time-dependent Dirichlet boundary conditions. > By definition, test spaces have zero Dirichlet boundary conditions so they need not be seen as time-dependent objects. @@ -140,7 +147,8 @@ U0 = U(t0) ∂ttU0 = ∂ttU(t0) ``` -## Cell fields +### Cell fields + The time-dependent equivalent of `CellField` is `TransientCellField`. It stores the cell field itself together with its derivatives up to the order of the ODE. For example, the following creates a `TransientCellField` with two time derivatives. @@ -151,7 +159,8 @@ u0 = zero(get_free_dof_values(U0)) u = TransientCellField(u0, (∂tu0, ∂ttu0)) ``` -## Finite element operators +### Finite element operators + The time-dependent analog of `FEOperator` is `TransientFEOperator`. It has the following constructors based on the nonlinearity type of the underlying ODE. * `TransientFEOperator(res, jacs, trial, test)` and `TransientFEOperator(res, trial, test; order)` for the version with automatic jacobians. The residual is expected to have the signature `residual(t, u, v)`. @@ -199,7 +208,8 @@ TransientLinearFEOperator((stiffness, mass), res, U, V, constant_forms=(false, t ``` If ``\kappa`` is constant, the keyword `constant_forms` could be replaced by `(true, true)`. -## The `TimeSpaceFunction` constructor +### The `TimeSpaceFunction` constructor + Apply differential operators on a function that depends on time and space is somewhat cumbersome. Let `f` be a function of time and space, and `g(t) = x -> f(t, x)` (as in the prescription of the boundary conditions `g` above). Applying the operator ``\partial_{t} - \Delta`` to `g` and evaluating at ``(t, x)`` is written `∂t(g)(t)(x) - Δ(g(t))(x)`. The constructor `TimeSpaceFunction` allows for simpler notations: let `h = TimeSpaceFunction(g)`. The object `h` is a functor that supports the notations @@ -209,7 +219,8 @@ The constructor `TimeSpaceFunction` allows for simpler notations: let `h = TimeS for all spatial and temporal differential operator, i.e. `op` in `(time_derivative, gradient, symmetric_gradient, divergence, curl, laplacian)` and their symbolic aliases (`∂t`, `∂tt`, `∇`, ...). The operator above applied to `h` and evaluated at `(t, x)` can be conveniently written `∂t(h)(t, x) - Δ(h)(t, x)`. -## Solver and solution +### Solver and solution + The next step is to choose an ODE solver (see below for a full list) and specify the boundary conditions. The solution can then be iterated over until the final time is reached. For example, to use the ``\theta``-method with a nonlinear solver, one could write @@ -234,32 +245,38 @@ for (tn, un) in enumerate(sol) end ``` -# Low-level implementation +## Low-level implementation + We now briefly describe the low-level implementation of the ODE module in `Gridap`. -## ODE operators +### ODE operators + The `ODEOperator` type represents an ODE according to the description above. It implements the `NonlinearOperator` interface, which enables the computation of residuals and jacobians. The algebraic equivalent of `TransientFEOperator` is an `ODEOpFromTFEOp`, which is a subtype of `ODEOperator`. Conceptually, `ODEOpFromTFEOp` can be thought of as an assembled `TransientFEOperator`, i.e. it deals with vectors of degrees of freedom. This operator comes with a cache (`ODEOpFromTFEOpCache`) that stores the transient space, its evaluation at the current time step, a cache for the `TransientFEOperator` itself (if any), and the constant forms (if any). > For now `TransientFEOperator` does not implement the `FEOperator` interface, i.e. it is not possible to evaluate residuals and jacobians directly on it. Rather, they are meant to be evaluated on the `ODEOpFromFEOp`. This is to cut down on the number of conversions between a `TransientCellField` and its vectors of degrees of freedom (one per time derivative). Indeed, when linear forms are constant, no conversion is needed as the jacobian matrix will be stored. -## ODE solvers +### ODE solvers + An ODE solver has to implement the following interface. * `allocate_odecache(odeslvr, odeop, t0, us0)`. This function allocates a cache that can be reused across the three functions `ode_start`, `ode_march!`, and `ode_finish!`. In particular, it is necessary to call `allocate_odeopcache` within this function, so as to instantiate the `ODEOpFromTFEOpCache` and be able to update the Dirichlet boundary conditions in the subsequent functions. * `ode_start(odeslvr, odeop, t0, us0, odecache)`. This function creates the state vectors from the initial conditions. By default, this is the identity. * `ode_march!(stateF, odeslvr, odeop, t0, state0, odecache)`. This is the update map that evolves the state vectors. * `ode_finish!(uF, odeslvr, odeop, t0, tF, stateF, odecache)`. This function converts the state vectors into the evaluation of the solution at the current time step. By default, this copies the first state vector into `uF`. -## Stage operator +### Stage operator + A `StageOperator` represents the linear or nonlinear operator that a numerical scheme relies on to evolve the state vector. It is essentially a special kind of `NonlinearOperator` but it overwrites the behaviour of nonlinear and linear solvers to take advantage of the matrix of a linear system being constant. The following subtypes of `StageOperator` are the building blocks of all numerical schemes. * `LinearStageOperator` represents the system ``\boldsymbol{J} \boldsymbol{x} + \boldsymbol{r} = \boldsymbol{0}``, and can build ``\boldsymbol{J}`` and ``\boldsymbol{r}`` by evaluating the residual at a given point. * `NonlinearStageOperator` represents ``\boldsymbol{r}(\boldsymbol{t}, \boldsymbol{\ell}_{0}(\boldsymbol{x}), \ldots, \boldsymbol{\ell}_{N}(\boldsymbol{x})) = \boldsymbol{0}``, where it is assumed that all the ``\boldsymbol{\ell}_{k}(\boldsymbol{x})`` are linear in ``\boldsymbol{x}``. -## ODE solution +### ODE solution + This type is a simple wrapper around an `ODEOperator`, an `ODESolver`, and initial conditions that can be iterated on to evolve the ODE. -# Numerical schemes formulation and implementation +## Numerical schemes formulation and implementation + We conclude this note by describing some numerical schemes and their implementation in `Gridap`. Suppose that the scheme has been evolved up to time ``t_{n}`` already and that the state vectors ``\{\boldsymbol{s}\}_{n}`` are known. We are willing to evolve the ODE up to time ``t_{n+1} > t_{n}``, i.e. compute the state vectors ``\{\boldsymbol{s}\}_{n+1}``. Generally speaking, a numerical scheme constructs an approximation of the map ``\{\boldsymbol{s}\}_{n} \to \{\boldsymbol{s}\}_{n+1}`` by solving one or more relationships of the type @@ -272,7 +289,8 @@ We now describe the numerical schemes implemented in `Gridap` using this framewo We also briefly characterise these schemes in terms of their order and linear stability. -## ``\theta``-method +### ``\theta``-method + This scheme is used to solve first-order ODEs and relies on the simple state vector ``\{\boldsymbol{s}(t)\} = \{\boldsymbol{u}(t)\}``. This means that the starting and finishing procedures are simply the identity. The ``\theta``-method relies on the following approximation @@ -312,7 +330,8 @@ By looking at the behaviour of the stability function at infinity, we find that * ``\theta = \frac{1}{2}``. The stability region is the whole left complex plane, so the scheme is ``A``-stable. This case is known as the implicit midpoint scheme. * ``\theta > \frac{1}{2}``. The stability region is the whole complex plane except the circle of radius ``\frac{1}{2 \theta - 1}`` centered at ``\left(\frac{1}{2 \theta - 1}, 0\right)``. In particular, the scheme is ``A``-stable. The special case ``\theta = 1`` is known as the Backward Euler scheme. -## Generalised-``\alpha`` scheme for first-order ODEs +### Generalised-``\alpha`` scheme for first-order ODEs + This scheme relies on the state vector ``\{\boldsymbol{s}(t)\} = \{\boldsymbol{u}(t), \partial_{t} \boldsymbol{u}(t)\}``. In particular, it needs a nontrivial starting procedure that evaluates ``\partial_{t} \boldsymbol{u}(t_{0})`` by enforcing a zero residual at ``t_{0}``. The finaliser can still return the first vector of the state vectors. For convenience, let ``\partial_{t} \boldsymbol{u}_{n}`` denote the approximation ``\partial_{t} \boldsymbol{u}(t_{n})``. > Alternatively, the initial velocity can be provided manually: when calling `solve(odeslvr, tfeop, t0, tF, uhs0)`, set `uhs0 = (u0, v0, a0)` instead of `uhs0 = (u0, v0)`. This is useful when enforcing a zero initial residual would lead to a singular system. @@ -344,6 +363,7 @@ t_{n + \alpha_{F}} &= (1 - \alpha_{F}) t_{n} + \alpha_{F} t_{n+1}, \\ The state vector is updated to ``\{\boldsymbol{s}\}_{n+1} = \{\boldsymbol{u}_{n+1}, \partial_{t} \boldsymbol{u}_{n+1}\}``. ##### Analysis + The amplification matrix for the state vector is ```math \boldsymbol{A}(z) = \frac{1}{\alpha_{M} - \alpha_{F} \gamma z} \begin{bmatrix}\alpha_{M} + (1 - \alpha_{F}) \gamma z & \alpha_{M} - \gamma \\ z & \alpha_{M} - 1 + \alpha_{F} (1 - \gamma) z\end{bmatrix}. @@ -384,7 +404,8 @@ This scheme was originally devised to control the damping of high frequencies. O ``` where ``\rho_{\infty}`` is the spectral radius at infinity. Setting ``\rho_{\infty}`` cuts all the highest frequencies in one step, whereas taking ``\rho_{\infty} = 1`` preserves high frequencies. -## Runge-Kutta +### Runge-Kutta + Runge-Kutta methods are multi-stage, i.e. they build estimates of ``\boldsymbol{u}`` at intermediate times between ``t_{n}`` and ``t_{n+1}``. They can be written as follows ```math \begin{align*} @@ -399,6 +420,7 @@ where ``p`` is the number of stages, ``\boldsymbol{A} = (a_{ij})_{1 \leq i, j \l **Implementation details** It is particularly advantageous to save the factorisation of the matrices of the stage operators for Runge-Kutta methods. This is always possible when the method is explicit and the mass matrix is constant, in which case all the stage matrices are the mass matrix. When the method is diagonally-implicit and the stiffness and mass matrices are constant, the matrices of the stage operators are ``\boldsymbol{M} + a_{ii} h_{n} \boldsymbol{K}``. In particular, if two diagonal coefficients coincide, the corresponding operators will have the same matrix. We implement these reuse strategies by storing them in `CompressedArray`s, and introducing a map `i -> NumericalSetup`. ##### Analysis + The stability function of a Runge-Kutta scheme is ```math \rho(z) = 1 + z \boldsymbol{b}^{T} (\boldsymbol{I} - z \boldsymbol{A})^{-1} \boldsymbol{1}. @@ -445,7 +467,8 @@ The analysis of Runge-Kutta methods is well-established but we only derive order \end{array}. ``` -## Implicit-Explicit Runge-Kutta +### Implicit-Explicit Runge-Kutta + When the residual has an implicit-explicit decomposition, usually because we can identify a stiff part that we want to solve implicitly and a nonstiff part that we want to solve explicitly, the Runge-Kutta method reads as follows ```math \begin{align*} @@ -475,7 +498,8 @@ Many methods can be created by padding a DIRK tableau with zeros to give it an a ``` We note that the first column of the matrix and the first weight are all zero, so the first stage for the implicit part does not need to be solved. -## Generalised-``\alpha`` scheme for second-order ODEs +### Generalised-``\alpha`` scheme for second-order ODEs + This scheme relies on the state vector ``\{\boldsymbol{s}(t)\} = \{\boldsymbol{u}(t), \partial_{t} \boldsymbol{u}(t), \partial_{tt} \boldsymbol{u}(t)\}``. It needs a nontrivial starting procedure that evaluates ``\partial_{tt} \boldsymbol{u}(t_{0})`` by enforcing a zero residual at ``t_{0}``. The finaliser can still return the first vector of the state vectors. For convenience, let ``\partial_{tt} \boldsymbol{u}_{n}`` denote the approximation ``\partial_{tt} \boldsymbol{u}(t_{n})``. > The initial acceleration can alternatively be provided manually: when calling `solve(odeslvr, tfeop, t0, tF, uhs0)`, set `uhs0 = (u0, v0, a0)` instead of `uhs0 = (u0, v0)`. This is useful when enforcing a zero initial residual would lead to a singular system. @@ -496,6 +520,7 @@ t_{n + 1 - \alpha_{F}} &= \alpha_{F} t_{n} + (1 - \alpha_{F}) t_{n+1}, \\ The state vector is then updated to ``\{\boldsymbol{s}\}_{n+1} = \{\boldsymbol{u}_{n+1}, \partial_{t} \boldsymbol{u}_{n+1}, \partial_{tt} \boldsymbol{u}_{n+1}\}``. ##### Analysis + The amplification matrix for the state vector is ```math \boldsymbol{A}(z) = \frac{1}{\overline{\alpha_{M}} + \overline{\alpha_{F}} \beta z^{2}} \begin{bmatrix} @@ -526,7 +551,7 @@ This method was also designed to damp high-frequency perturbations so it is comm * The standard generalised-``\alpha`` method is obtained by setting ``\alpha_{M} = \frac{2 \rho_{\infty - 1}}{\rho_{\infty} + 1}``, ``\alpha_{F} = \frac{\rho_{\infty}}{\rho_{\infty} + 1}``. * The Newmark method corresponds to ``\alpha_{F} = \alpha_{M} = 0``. In this case, the values of ``\beta`` and ``\gamma`` are usually chosen as ``\beta = 0``, ``\gamma = \frac{1}{2}`` (explicit central difference scheme), or ``\beta = \frac{1}{4}`` and ``\gamma = \frac{1}{2}`` (midpoint rule). -# Reference +## Reference ```@autodocs Modules = [ODEs,] diff --git a/docs/src/ReferenceFEs.md b/docs/src/ReferenceFEs.md index ff6bee434..8cc845992 100644 --- a/docs/src/ReferenceFEs.md +++ b/docs/src/ReferenceFEs.md @@ -73,8 +73,18 @@ Pages = ["LagrangianRefFEs.jl","LagrangianDofBases.jl","SerendipityRefFEs.jl", ### Moment-Based ReferenceFEs +#### Framework + +```@autodocs +Modules = [ReferenceFEs,] +Order = [:type, :constant, :macro, :function] +Pages = ["MomentBasedReferenceFEs.jl","Pullbacks.jl"] +``` + +#### Available Moment-Based ReferenceFEs + ```@autodocs Modules = [ReferenceFEs,] Order = [:type, :constant, :macro, :function] -Pages = ["MomentBasedReferenceFEs.jl","RaviartThomasRefFEs.jl","NedelecRefFEs.jl","BDMRefFEs.jl","CRRefFEs.jl"] +Pages = ["RaviartThomasRefFEs.jl","NedelecRefFEs.jl","BDMRefFEs.jl","CRRefFEs.jl"] ``` diff --git a/docs/src/dev-notes/pullbacks.md b/docs/src/dev-notes/pullbacks.md new file mode 100644 index 000000000..2a2a98d13 --- /dev/null +++ b/docs/src/dev-notes/pullbacks.md @@ -0,0 +1,82 @@ + +# FE basis transformations + +## Notation + +Consider a reference polytope ``\hat{K}``, mapped to the physical space by a **geometrical map** ``F``, i.e. ``K = F(\hat{K})``. Consider also a linear function space on the reference polytope ``\hat{V}``, and a set of unisolvent degrees of freedom represented by moments in the dual space ``\hat{V}^*``. + +Throughout this document, we will use the following notation: + +- ``\varphi \in V`` is a **physical field** ``\varphi : K \rightarrow \mathbb{R}^k``. A basis of ``V`` is denoted by ``\Phi = \{\varphi\}``. +- ``\hat{\varphi} \in \hat{V}`` is a **reference field** ``\hat{\varphi} : \hat{K} \rightarrow \mathbb{R}^k``. A basis of ``\hat{V}`` is denoted by ``\hat{\Phi} = \{\hat{\varphi}\}``. +- ``\sigma \in V^*`` is a **physical moment** ``\sigma : V \rightarrow \mathbb{R}``. A basis of ``V^*`` is denoted by ``\Sigma = \{\sigma\}``. +- ``\hat{\sigma} \in \hat{V}^*`` is a **reference moment** ``\hat{\sigma} : \hat{V} \rightarrow \mathbb{R}``. A basis of ``\hat{V}^*`` is denoted by ``\hat{\Sigma} = \{\hat{\sigma}\}``. + +## Pullbacks and Pushforwards + +We define a **pushforward** map as ``F^* : \hat{V} \rightarrow V``, mapping reference fields to physical fields. Given a pushforward ``F^*``, we define: + +- The **pullback** ``F_* : V^* \rightarrow \hat{V}^*``, mapping physical moments to reference moments. Its action on physical dofs is defined in terms of the pushforward map ``F^*`` as ``\hat{\sigma} = F_*(\sigma) := \sigma \circ F^*``. +- The **inverse pushforward** ``(F^*)^{-1} : V \rightarrow \hat{V}``, mapping physical fields to reference fields. +- The **inverse pullback** ``(F_*)^{-1} : \hat{V}^* \rightarrow V^*``, mapping reference moments to physical moments. Its action on reference dofs is defined in terms of the inverse pushforward map ``(F^*)^{-1}`` as ``\sigma = (F_*)^{-1}(\hat{\sigma}) := \hat{\sigma} \circ (F^*)^{-1}``. + +## Change of basis + +In many occasions, we will have that (as a basis) + +```math +\hat{\Sigma} \neq F_*(\Sigma), \quad \text{and} \quad \Phi \neq F^*(\hat{\Phi}) +``` + +To maintain conformity and proper scaling in these cases, we define cell-dependent invertible changes of basis ``P`` and ``M``, such that + +```math +\hat{\Sigma} = P F_*(\Sigma), \quad \text{and} \quad \Phi = M F^*(\hat{\Phi}) +``` + +An important result from [1, Theorem 3.1] is that ``P = M^T``. + +!!! details + [1, Lemma 2.6]: A key ingredient is that given ``M`` a matrix we have ``\Sigma (M \Phi) = \Sigma (\Phi) M^T`` since + ```math + [\Sigma (M \Phi)]_{ij} = \sigma_i (M_{jk} \varphi_k) = M_{jk} \sigma_i (\varphi_k) = [\Sigma (\Phi) M^T]_{ij} + ``` + where we have used that moments are linear. + +We then have the following diagram: + +```math +\hat{V}^* \xleftarrow{P} \hat{V}^* \xleftarrow{F_*} V^* \\ +\hat{V} \xrightarrow{F^*} V \xrightarrow{P^T} V +``` + +!!! details + The above diagram is well defined, since we have + ```math + \hat{\Sigma}(\hat{\Phi}) = P F_* (\Sigma)(F^{-*} (P^{-T} \Phi)) = P \Sigma (F^* (F^{-*} P^{-T} \Phi)) = P \Sigma (P^{-T} \Phi) = P \Sigma (\Phi) P^{-1} = Id \\ + \Sigma(\Phi) = F_*^{-1}(P^{-1}\hat{\Sigma})(P^T F^*(\hat{\Phi})) = P^{-1} \hat{\Sigma} (F^{-*}(P^T F^*(\hat{\Phi}))) = P^{-1} \hat{\Sigma} (P^T \hat{\Phi}) = P^{-1} \hat{\Sigma}(\hat{\Phi}) P = Id + ``` + +From an implementation point of view, it is more natural to build ``P^{-1}`` and then retrieve all other matrices by transposition/inversion. + +## Interpolation + +In each cell ``K`` and for ``C_b^k(K)`` the space of functions defined on ``K`` with at least ``k`` bounded derivatives, we define the interpolation operator ``I_K : C_b^k(K) \rightarrow V`` as + +```math +I_K(g) = \Sigma(g) \Phi \quad, \quad \Sigma(g) = P^{-1} \hat{\Sigma}(F^{-*}(g)) +``` + +## Implementation notes + +!!! note + In [2], Covariant and Contravariant Piola maps preserve exactly (without any sign change) the normal and tangential components of a vector field. + I am quite sure that the discrepancy is coming from the fact that the geometrical information in the reference polytope is globally oriented. + For instance, the normals ``n`` and ``\hat{n}`` both have the same orientation, i.e ``n = (||\hat{e}||/||e||) (det J) J^{-T} \hat{n}``. Therefore ``\hat{n}`` is not fully local. See [2, Equation 2.11]. + In our case, we will be including the sign change in the transformation matrices, which will include all cell-and-dof-dependent information. + +## References + +[1] [Kirby 2017, A general approach to transforming finite elements.](https://arxiv.org/abs/1706.09017) + +[2] [Aznaran et al. 2021, Transformations for Piola-mapped elements.](https://arxiv.org/abs/2110.13224) diff --git a/src/Adaptivity/RefinementRules.jl b/src/Adaptivity/RefinementRules.jl index a93d93d20..d617cb091 100644 --- a/src/Adaptivity/RefinementRules.jl +++ b/src/Adaptivity/RefinementRules.jl @@ -25,7 +25,7 @@ end # - The reason why we are saving both the cell maps and the inverse cell maps is to avoid recomputing # them when needed. This is needed for performance when the RefinementRule is used for MacroFEs. # Also, in the case the ref_grid comes from a CartesianGrid, we save the cell maps as -# AffineMaps, which are more efficient than the default linear_combinations. +# AffineFields, which are more efficient than the default linear_combinations. # - We cannot parametrise the RefinementRule by all it's fields, because we will have different types of # RefinementRules in a single mesh. It's the same reason why we don't parametrise the ReferenceFE type. diff --git a/src/Adaptivity/deprecated/FineToCoarseReferenceFEs.jl b/src/Adaptivity/deprecated/FineToCoarseReferenceFEs.jl new file mode 100644 index 000000000..f36ad5793 --- /dev/null +++ b/src/Adaptivity/deprecated/FineToCoarseReferenceFEs.jl @@ -0,0 +1,132 @@ +""" +""" +struct FineToCoarseDofBasis{T,A,B,C} <: AbstractVector{T} + dof_basis :: A + rrule :: B + child_ids :: C + + function FineToCoarseDofBasis(dof_basis::AbstractVector{T},rrule::RefinementRule) where {T<:Dof} + nodes = get_nodes(dof_basis) + child_ids = map(x -> x_to_cell(rrule,x),nodes) + + A = typeof(dof_basis) + B = typeof(rrule) + C = typeof(child_ids) + new{T,A,B,C}(dof_basis,rrule,child_ids) + end +end + +Base.size(a::FineToCoarseDofBasis) = size(a.dof_basis) +Base.axes(a::FineToCoarseDofBasis) = axes(a.dof_basis) +Base.getindex(a::FineToCoarseDofBasis,i::Integer) = getindex(a.dof_basis,i) +Base.IndexStyle(a::FineToCoarseDofBasis) = IndexStyle(a.dof_basis) + +ReferenceFEs.get_nodes(a::FineToCoarseDofBasis) = get_nodes(a.dof_basis) + +# Default behaviour +Arrays.return_cache(b::FineToCoarseDofBasis,field) = return_cache(b.dof_basis,field) +Arrays.evaluate!(cache,b::FineToCoarseDofBasis,field) = evaluate!(cache,b.dof_basis,field) + +# Spetialized behaviour +function Arrays.return_cache(s::FineToCoarseDofBasis{T,<:LagrangianDofBasis},field::FineToCoarseField) where T + b = s.dof_basis + cf = return_cache(field,b.nodes,s.child_ids) + vals = evaluate!(cf,field,b.nodes,s.child_ids) + ndofs = length(b.dof_to_node) + r = ReferenceFEs._lagr_dof_cache(vals,ndofs) + c = CachedArray(r) + return (c, cf) +end + +function Arrays.evaluate!(cache,s::FineToCoarseDofBasis{T,<:LagrangianDofBasis},field::FineToCoarseField) where T + c, cf = cache + b = s.dof_basis + vals = evaluate!(cf,field,b.nodes,s.child_ids) + ndofs = length(b.dof_to_node) + T2 = eltype(vals) + ncomps = num_indep_components(T2) + @check ncomps == num_indep_components(eltype(b.node_and_comp_to_dof)) """\n + Unable to evaluate LagrangianDofBasis. The number of components of the + given Field does not match with the LagrangianDofBasis. + If you are trying to interpolate a function on a FESpace make sure that + both objects have the same value type. + For instance, trying to interpolate a vector-valued function on a scalar-valued FE space + would raise this error. + """ + ReferenceFEs._evaluate_lagr_dof!(c,vals,b.node_and_comp_to_dof,ndofs,ncomps) +end + +function Arrays.return_cache(s::FineToCoarseDofBasis{T,<:MomentBasedDofBasis},field::FineToCoarseField) where T + b = s.dof_basis + cf = return_cache(field,b.nodes,s.child_ids) + vals = evaluate!(cf,field,b.nodes,s.child_ids) + ndofs = num_dofs(b) + r = ReferenceFEs._moment_dof_basis_cache(vals,ndofs) + c = CachedArray(r) + return (c, cf) +end + +function Arrays.evaluate!(cache,s::FineToCoarseDofBasis{T,<:MomentBasedDofBasis},field::FineToCoarseField) where T + c, cf = cache + b = s.dof_basis + vals = evaluate!(cf,field,b.nodes,s.child_ids) + dofs = c.array + ReferenceFEs._eval_moment_dof_basis!(dofs,vals,b) + dofs +end + + +""" + Wrapper for a ReferenceFE which is specialised for + efficiently evaluating FineToCoarseFields. +""" +struct FineToCoarseRefFE{T,D,A} <: ReferenceFE{D} + reffe :: T + dof_basis :: A + + function FineToCoarseRefFE(reffe::ReferenceFE{D},dof_basis::FineToCoarseDofBasis) where D + T = typeof(reffe) + A = typeof(dof_basis) + new{T,D,A}(reffe,dof_basis) + end +end + +ReferenceFEs.num_dofs(reffe::FineToCoarseRefFE) = num_dofs(reffe.reffe) +ReferenceFEs.get_polytope(reffe::FineToCoarseRefFE) = get_polytope(reffe.reffe) +ReferenceFEs.get_prebasis(reffe::FineToCoarseRefFE) = get_prebasis(reffe.reffe) +ReferenceFEs.get_dof_basis(reffe::FineToCoarseRefFE) = reffe.dof_basis +ReferenceFEs.Conformity(reffe::FineToCoarseRefFE) = Conformity(reffe.reffe) +ReferenceFEs.get_face_dofs(reffe::FineToCoarseRefFE) = get_face_dofs(reffe.reffe) +ReferenceFEs.get_shapefuns(reffe::FineToCoarseRefFE) = get_shapefuns(reffe.reffe) +ReferenceFEs.get_metadata(reffe::FineToCoarseRefFE) = get_metadata(reffe.reffe) +ReferenceFEs.get_orders(reffe::FineToCoarseRefFE) = get_orders(reffe.reffe) +ReferenceFEs.get_order(reffe::FineToCoarseRefFE) = get_order(reffe.reffe) + +ReferenceFEs.Conformity(reffe::FineToCoarseRefFE,sym::Symbol) = Conformity(reffe.reffe,sym) +ReferenceFEs.get_face_own_dofs(reffe::FineToCoarseRefFE,conf::Conformity) = get_face_own_dofs(reffe.reffe,conf) + + +function ReferenceFEs.ReferenceFE(p::Polytope,rrule::RefinementRule,name::ReferenceFEName,order) + FineToCoarseRefFE(p,rrule,name,Float64,order) +end + +function ReferenceFEs.ReferenceFE(p::Polytope,rrule::RefinementRule,name::ReferenceFEName,::Type{T},order) where T + FineToCoarseRefFE(p,rrule,name,T,order) +end + +function FineToCoarseRefFE(p::Polytope,rrule::RefinementRule,name::ReferenceFEName,::Type{T},order) where T + @check p == get_polytope(rrule) + reffe = ReferenceFE(p,name,T,order) + dof_basis = FineToCoarseDofBasis(get_dof_basis(reffe),rrule) + return FineToCoarseRefFE(reffe,dof_basis) +end + +# FESpaces constructors + +function FESpaces.TestFESpace(model::DiscreteModel,rrules::AbstractVector{<:RefinementRule},reffe::Tuple{<:ReferenceFEName,Any,Any};kwargs...) + @check num_cells(model) == length(rrules) + @check all(CompressedArray(get_polytopes(model),get_cell_type(model)) .== lazy_map(get_polytope,rrules)) + basis, reffe_args, reffe_kwargs = reffe + reffes = lazy_map(rr -> ReferenceFE(get_polytope(rr),rr,basis,reffe_args...;reffe_kwargs...),rrules) + return TestFESpace(model,reffes;kwargs...) +end diff --git a/src/Arrays/Arrays.jl b/src/Arrays/Arrays.jl index ce7452453..54f39222a 100644 --- a/src/Arrays/Arrays.jl +++ b/src/Arrays/Arrays.jl @@ -47,9 +47,8 @@ export testargs export inverse_map export Broadcasting - export Operation - +export InverseMap # LazyArray diff --git a/src/Arrays/Maps.jl b/src/Arrays/Maps.jl index 7552b80c8..77b5729cf 100644 --- a/src/Arrays/Maps.jl +++ b/src/Arrays/Maps.jl @@ -296,3 +296,16 @@ function inverse_map(f) Function inverse_map is not implemented yet for objects of type $(typeof(f)) """ end + +struct InverseMap{F} <: Map + original::F +end + +function evaluate!(cache,k::InverseMap,args...) + @notimplemented """\n + The inverse evaluation is not implemented yet for maps of type $(typeof(k.original)) + """ +end + +inverse_map(k::Map) = InverseMap(k) +inverse_map(k::InverseMap) = k.original diff --git a/src/FESpaces/ConstantFESpaces.jl b/src/FESpaces/ConstantFESpaces.jl index f7115f46c..8d1f85a8d 100644 --- a/src/FESpaces/ConstantFESpaces.jl +++ b/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() diff --git a/src/FESpaces/FESpaces.jl b/src/FESpaces/FESpaces.jl index b42054859..5537ed784 100644 --- a/src/FESpaces/FESpaces.jl +++ b/src/FESpaces/FESpaces.jl @@ -219,9 +219,7 @@ include("UnconstrainedFESpaces.jl") include("ConformingFESpaces.jl") -include("DivConformingFESpaces.jl") - -include("CurlConformingFESpaces.jl") +include("Pullbacks.jl") include("FESpaceFactories.jl") @@ -253,8 +251,6 @@ include("CLagrangianFESpaces.jl") include("DirichletFESpaces.jl") -#include("ExtendedFESpaces.jl") - include("FESpacesWithLinearConstraints.jl") include("DiscreteModelWithFEMaps.jl") diff --git a/src/FESpaces/Pullbacks.jl b/src/FESpaces/Pullbacks.jl new file mode 100644 index 000000000..26386fc1d --- /dev/null +++ b/src/FESpaces/Pullbacks.jl @@ -0,0 +1,129 @@ + +# TODO: We probably want to export these from Gridap.ReferenceFEs +using Gridap.ReferenceFEs: PushforwardRefFE, Pushforward, Pullback +using Gridap.ReferenceFEs: ContraVariantPiolaMap, CoVariantPiolaMap + +function get_cell_dof_basis( + model::DiscreteModel, cell_reffe::AbstractArray{<:GenericRefFE{T}}, conformity::Conformity +) where T <: PushforwardRefFE + pushforward, cell_change, cell_args = get_cell_pushforward( + Pushforward(T), model, cell_reffe, conformity + ) + cell_ref_dofs = lazy_map(get_dof_basis, cell_reffe) + cell_phy_dofs = lazy_map(inverse_map(Pullback(pushforward)), cell_ref_dofs, cell_args...) + return lazy_map(linear_combination, cell_change, cell_phy_dofs) # TODO: Inverse and transpose +end + +function get_cell_shapefuns( + model::DiscreteModel, cell_reffe::AbstractArray{<:GenericRefFE{T}}, conformity::Conformity +) where T <: PushforwardRefFE + pushforward, cell_change, cell_args = get_cell_pushforward( + Pushforward(T), model, cell_reffe, conformity + ) + cell_ref_fields = lazy_map(get_shapefuns, cell_reffe) + cell_phy_fields = lazy_map(pushforward, cell_ref_fields, cell_args...) + return lazy_map(linear_combination, cell_change, cell_phy_fields) +end + +function get_cell_pushforward( + ::Pushforward, model::DiscreteModel, cell_reffe, conformity +) + @abstractmethod +end + +# ContraVariantPiolaMap + +function get_cell_pushforward( + ::ContraVariantPiolaMap, model::DiscreteModel, cell_reffe, conformity +) + cell_map = get_cell_map(get_grid(model)) + Jt = lazy_map(Broadcasting(∇),cell_map) + change = get_sign_flip(model, cell_reffe) + return ContraVariantPiolaMap(), change, (Jt,) +end + +# CoVariantPiolaMap + +function get_cell_pushforward( + ::CoVariantPiolaMap, model::DiscreteModel, cell_reffe, conformity +) + cell_map = get_cell_map(get_grid(model)) + Jt = lazy_map(Broadcasting(∇),cell_map) + change = lazy_map(r -> Diagonal(ones(num_dofs(r))), cell_reffe) # TODO: Replace by edge-signs + return CoVariantPiolaMap(), change, (Jt,) +end + +# NormalSignMap + +""" + struct NormalSignMap <: Map + ... + end +""" +struct NormalSignMap{T} <: Map + model::T + facet_owners::Vector{Int32} +end + +function NormalSignMap(model) + facet_owners = compute_facet_owners(model) + NormalSignMap(model,facet_owners) +end + +function return_value(k::NormalSignMap,reffe,facet_own_dofs,cell) + Diagonal(fill(one(Float64), num_dofs(reffe))) +end + +function return_cache(k::NormalSignMap,reffe,facet_own_dofs,cell) + model = k.model + Dc = num_cell_dims(model) + topo = get_grid_topology(model) + + cell_facets = get_faces(topo, Dc, Dc-1) + cell_facets_cache = array_cache(cell_facets) + + return cell_facets, cell_facets_cache, CachedVector(Float64) +end + +function evaluate!(cache,k::NormalSignMap,reffe,facet_own_dofs,cell) + cell_facets,cell_facets_cache,dof_sign_cache = cache + facet_owners = k.facet_owners + + setsize!(dof_sign_cache, (num_dofs(reffe),)) + dof_sign = dof_sign_cache.array + fill!(dof_sign, one(eltype(dof_sign))) + + facets = getindex!(cell_facets_cache,cell_facets,cell) + for (lfacet,facet) in enumerate(facets) + owner = facet_owners[facet] + if owner != cell + for dof in facet_own_dofs[lfacet] + dof_sign[dof] = -1.0 + end + end + end + + return Diagonal(dof_sign) +end + +function get_sign_flip(model::DiscreteModel{Dc}, cell_reffe) where Dc + # Comment: lazy_maps on cell_reffes are very optimised, since they are CompressedArray/FillArray + get_facet_own_dofs(reffe) = view(get_face_own_dofs(reffe),get_dimrange(get_polytope(reffe),Dc-1)) + cell_facet_own_dofs = lazy_map(get_facet_own_dofs, cell_reffe) + cell_ids = IdentityVector(Int32(num_cells(model))) + return lazy_map(NormalSignMap(model), cell_reffe, cell_facet_own_dofs, cell_ids) +end + +function compute_facet_owners(model::DiscreteModel{Dc}) where {Dc} + topo = get_grid_topology(model) + facet_to_cell = get_faces(topo, Dc-1, Dc) + + nfacets = num_faces(topo, Dc-1) + owners = Vector{Int32}(undef, nfacets) + for facet in 1:nfacets + facet_cells = view(facet_to_cell, facet) + owners[facet] = first(facet_cells) + end + + return owners +end diff --git a/src/FESpaces/CurlConformingFESpaces.jl b/src/FESpaces/deprecated/CurlConformingFESpaces.jl similarity index 94% rename from src/FESpaces/CurlConformingFESpaces.jl rename to src/FESpaces/deprecated/CurlConformingFESpaces.jl index 0bf9e45cc..55b3d484a 100644 --- a/src/FESpaces/CurlConformingFESpaces.jl +++ b/src/FESpaces/deprecated/CurlConformingFESpaces.jl @@ -45,7 +45,7 @@ function get_cell_shapefuns( cell_reffe_shapefuns = lazy_map(get_shapefuns,cell_reffe) cell_map = get_cell_map(Triangulation(model)) + cell_Jt = lazy_map(Broadcasting(∇),cell_map) k = ReferenceFEs.CoVariantPiolaMap() - lazy_map(k,cell_reffe_shapefuns,cell_map) + lazy_map(k,cell_reffe_shapefuns,cell_Jt) end - diff --git a/src/FESpaces/DivConformingFESpaces.jl b/src/FESpaces/deprecated/DivConformingFESpaces.jl similarity index 97% rename from src/FESpaces/DivConformingFESpaces.jl rename to src/FESpaces/deprecated/DivConformingFESpaces.jl index db08c8ba2..e72724247 100644 --- a/src/FESpaces/DivConformingFESpaces.jl +++ b/src/FESpaces/deprecated/DivConformingFESpaces.jl @@ -47,9 +47,10 @@ function get_cell_shapefuns( sign_flip = get_sign_flip(model, cell_reffe) ) cell_map = get_cell_map(get_grid(model)) + cell_Jt = lazy_map(Broadcasting(∇),cell_map) cell_shapefuns = lazy_map(get_shapefuns,cell_reffe) k = ContraVariantPiolaMap() - lazy_map(k,cell_shapefuns,cell_map,lazy_map(Broadcasting(constant_field),sign_flip)) + lazy_map(k,cell_shapefuns,cell_Jt,lazy_map(Broadcasting(constant_field),sign_flip)) end struct SignFlipMap{T} <: Map diff --git a/src/Fields/AffineMaps.jl b/src/Fields/AffineMaps.jl index 9ce5725fc..33bbe385f 100644 --- a/src/Fields/AffineMaps.jl +++ b/src/Fields/AffineMaps.jl @@ -1,12 +1,25 @@ +struct AffineMap <: Map end + +function return_cache(::AffineMap,G::TensorValue{D1,D2},y0::Point{D2},x::Point{D1}) where {D1,D2} + nothing +end + +function evaluate!( + cache,::AffineMap,G::TensorValue{D1,D2},y0::Point{D2},x::Point{D1} +) where {D1,D2} + x⋅G + y0 +end + + """ A Field with this form y = x⋅G + y0 """ -struct AffineMap{D1,D2,T,L} <:Field +struct AffineField{D1,D2,T,L} <: Field gradient::TensorValue{D1,D2,T,L} origin::Point{D2,T} - function AffineMap( + function AffineField( gradient::TensorValue{D1,D2,T,L}, origin::Point{D2,T}) where {D1,D2,T,L} @@ -14,21 +27,27 @@ struct AffineMap{D1,D2,T,L} <:Field end end -affine_map(gradient,origin) = AffineMap(gradient,origin) +affine_map(gradient,origin) = AffineField(gradient,origin) -function evaluate!(cache,f::AffineMap,x::Point) +function Base.zero(::Type{<:AffineField{D1,D2,T}}) where {D1,D2,T} + gradient = TensorValue{D1,D2}(tfill(zero(T),Val{D1*D2}())) + origin = Point{D2,T}(tfill(zero(T),Val{D2}())) + AffineField(gradient,origin) +end + +function evaluate!(cache,f::AffineField,x::Point) G = f.gradient y0 = f.origin x⋅G + y0 end -function return_cache(f::AffineMap,x::AbstractVector{<:Point}) +function return_cache(f::AffineField,x::AbstractVector{<:Point}) T = return_type(f,testitem(x)) y = similar(x,T,size(x)) CachedArray(y) end -function evaluate!(cache,f::AffineMap,x::AbstractVector{<:Point}) +function evaluate!(cache,f::AffineField,x::AbstractVector{<:Point}) setsize!(cache,size(x)) y = cache.array G = f.gradient @@ -41,11 +60,11 @@ function evaluate!(cache,f::AffineMap,x::AbstractVector{<:Point}) y end -function gradient(h::AffineMap) +function gradient(h::AffineField) ConstantField(h.gradient) end -function push_∇∇(∇∇a::Field,ϕ::AffineMap) +function push_∇∇(∇∇a::Field,ϕ::AffineField) # Assuming ϕ is affine map Jt = ∇(ϕ) Jt_inv = pinvJt(Jt) @@ -80,7 +99,7 @@ end function lazy_map( k::Broadcasting{typeof(push_∇∇)}, cell_∇∇a::AbstractArray, - cell_map::AbstractArray{<:AffineMap}) + cell_map::AbstractArray{<:AffineField}) cell_Jt = lazy_map(∇,cell_map) cell_invJt = lazy_map(Operation(pinvJt),cell_Jt) lazy_map(Broadcasting(Operation(push_∇∇)),cell_∇∇a,cell_invJt) @@ -89,19 +108,19 @@ end function lazy_map( k::Broadcasting{typeof(push_∇∇)}, cell_∇∇at::LazyArray{<:Fill{typeof(transpose)}}, - cell_map::AbstractArray{<:AffineMap}) + cell_map::AbstractArray{<:AffineField}) cell_∇∇a = cell_∇∇at.args[1] cell_∇∇b = lazy_map(k,cell_∇∇a,cell_map) cell_∇∇bt = lazy_map(transpose,cell_∇∇b) cell_∇∇bt end -function inverse_map(f::AffineMap) +function inverse_map(f::AffineField) Jt = f.gradient y0 = f.origin invJt = pinvJt(Jt) x0 = -y0⋅invJt - AffineMap(invJt,x0) + AffineField(invJt,x0) end function lazy_map(::typeof(∇),a::LazyArray{<:Fill{typeof(affine_map)}}) @@ -109,8 +128,10 @@ function lazy_map(::typeof(∇),a::LazyArray{<:Fill{typeof(affine_map)}}) lazy_map(constant_field,gradients) end -function Base.zero(::Type{<:AffineMap{D1,D2,T}}) where {D1,D2,T} - gradient = TensorValue{D1,D2}(tfill(zero(T),Val{D1*D2}())) - origin = Point{D2,T}(tfill(zero(T),Val{D2}())) - AffineMap(gradient,origin) +function lazy_map( + ::typeof(evaluate),a::LazyArray{<:Fill{typeof(affine_map)}},x::AbstractArray +) + gradients = a.args[1] + origins = a.args[2] + lazy_map(Broadcasting(AffineMap()),gradients,origins,x) end diff --git a/src/Fields/ApplyOptimizations.jl b/src/Fields/ApplyOptimizations.jl index 76dd091fd..95df7c730 100644 --- a/src/Fields/ApplyOptimizations.jl +++ b/src/Fields/ApplyOptimizations.jl @@ -135,14 +135,14 @@ function lazy_map( k::Broadcasting{typeof(∇)}, a::LazyArray{<:Fill{typeof(transpose)}}) i_to_basis = lazy_map(k,a.args[1]) - lazy_map( transpose, i_to_basis) + lazy_map(transpose, i_to_basis) end function lazy_map( k::Broadcasting{typeof(∇∇)}, a::LazyArray{<:Fill{typeof(transpose)}}) i_to_basis = lazy_map(k,a.args[1]) - lazy_map( transpose, i_to_basis) + lazy_map(transpose, i_to_basis) end # Gradient rules diff --git a/src/Fields/FieldArrays.jl b/src/Fields/FieldArrays.jl index 4eba60280..485de481b 100644 --- a/src/Fields/FieldArrays.jl +++ b/src/Fields/FieldArrays.jl @@ -368,6 +368,16 @@ function evaluate!(cache,k::LinearCombinationMap{Colon},v::AbstractMatrix,fx::Ab r end +function evaluate!(cache,k::LinearCombinationMap{Colon},v::LinearAlgebra.Diagonal,fx::AbstractVector) + @check length(fx) == size(v,1) + setsize!(cache,(size(v,2),)) + r = cache.array + @inbounds for j in eachindex(fx) + r[j] = outer(fx[j],v.diag[j]) + end + r +end + function return_cache(k::LinearCombinationMap{Colon},v::AbstractMatrix,fx::AbstractMatrix) vf = testitem(fx) vv = testitem(v) @@ -392,6 +402,18 @@ function evaluate!(cache,k::LinearCombinationMap{Colon},v::AbstractMatrix,fx::Ab r end +function evaluate!(cache,k::LinearCombinationMap{Colon},v::LinearAlgebra.Diagonal,fx::AbstractMatrix) + @check size(fx,2) == size(v,1) + setsize!(cache,(size(fx,1),size(v,2))) + r = cache.array + @inbounds for p in 1:size(fx,1) + for j in 1:size(fx,2) + r[p,j] = outer(fx[p,j],v.diag[j]) + end + end + r +end + # Optimizing transpose testitem(a::Transpose{<:Field}) = testitem(a.parent) evaluate!(cache,k::Broadcasting{typeof(∇)},a::Transpose{<:Field}) = transpose(k(a.parent)) diff --git a/src/Fields/Fields.jl b/src/Fields/Fields.jl index 05e5e7a56..0dd060b4c 100644 --- a/src/Fields/Fields.jl +++ b/src/Fields/Fields.jl @@ -56,7 +56,7 @@ export MockFieldArray export Point export inverse_map -export AffineMap +export AffineField export affine_map export gradient diff --git a/src/Fields/FieldsInterfaces.jl b/src/Fields/FieldsInterfaces.jl index bdad6a111..799bc96cc 100644 --- a/src/Fields/FieldsInterfaces.jl +++ b/src/Fields/FieldsInterfaces.jl @@ -288,6 +288,16 @@ function lazy_map(::Operation{typeof(inv)},a::LazyArray{<:Fill{typeof(constant_f lazy_map(constant_field,vinv) end +struct ConstantMap <: Map end + +return_cache(::ConstantMap,v::Number,x::Point) = nothing +evaluate!(cache,::ConstantMap,v::Number,x::Point) = v + +function lazy_map(::typeof(evaluate),a::LazyArray{<:Fill{typeof(constant_field)}},x::AbstractArray) + values = a.args[1] + lazy_map(Broadcasting(ConstantMap()),values,x) +end + ## Make Function behave like Field return_cache(f::FieldGradient{N,<:Function},x::Point) where N = gradient(f.object,Val(N)) diff --git a/src/Fields/InverseFields.jl b/src/Fields/InverseFields.jl index 2248eda8a..f760f3dc6 100644 --- a/src/Fields/InverseFields.jl +++ b/src/Fields/InverseFields.jl @@ -8,8 +8,8 @@ inverse_map(a::Field) = InverseField(a) inverse_map(a::InverseField) = a.original function return_cache(a::InverseField,x::Point) - y₀ = [zero(x)...] # initial guess and solution - F₀ = [zero(x)...] # error + y₀ = [zero(x)...] # initial guess and solution + F₀ = [zero(x)...] # error return return_cache(a.original,x), return_cache(∇(a.original),x), y₀, F₀ end diff --git a/src/Geometry/CartesianGrids.jl b/src/Geometry/CartesianGrids.jl index ee83694c4..5fd239957 100644 --- a/src/Geometry/CartesianGrids.jl +++ b/src/Geometry/CartesianGrids.jl @@ -237,7 +237,7 @@ end # Cell map -struct CartesianMap{D,T,L} <: AbstractArray{AffineMap{D,D,T,L},D} +struct CartesianMap{D,T,L} <: AbstractArray{AffineField{D,D,T,L},D} data::CartesianDescriptor{D,T,typeof(identity)} function CartesianMap(des::CartesianDescriptor{D,T}) where {D,T} L = D*D @@ -256,9 +256,9 @@ function Base.getindex(a::CartesianMap{D,T},I::Vararg{Integer,D}) where {D,T} @inbounds for d in 1:D p[d] = x0[d] + (I[d]-1)*dx[d] end - origin = Point(p) + origin = Point(p) grad = diagonal_tensor(VectorValue(dx)) - AffineMap(grad,origin) + AffineField(grad,origin) end function lazy_map(::typeof(∇),a::CartesianMap) diff --git a/src/ReferenceFEs/AWRefFEs.jl b/src/ReferenceFEs/AWRefFEs.jl index 2f6850afb..72a9a527f 100644 --- a/src/ReferenceFEs/AWRefFEs.jl +++ b/src/ReferenceFEs/AWRefFEs.jl @@ -1,9 +1,11 @@ -struct ArnoldWinther <: ReferenceFEName end + +struct ArnoldWinther <: PushforwardRefFE end const arnoldwinther = ArnoldWinther() """ -ArnoldWintherRefFE(::Type{T},p::Polytope,order::Integer) where T + struct ArnoldWinther <: PushforwardRefFE end + ArnoldWintherRefFE(::Type{T},p::Polytope,order::Integer) where T Arnold-Winther reference finite element. @@ -55,11 +57,11 @@ function ArnoldWintherRefFE(::Type{T},p::Polytope,order::Integer) where T end function ReferenceFE(p::Polytope,::ArnoldWinther, order) - BDMRefFE(Float64,p,order) + ArnoldWintherRefFE(Float64,p,order) end function ReferenceFE(p::Polytope,::ArnoldWinther,::Type{T}, order) where T - BDMRefFE(T,p,order) + ArnoldWintherRefFE(T,p,order) end function Conformity(reffe::GenericRefFE{ArnoldWinther},sym::Symbol) diff --git a/src/ReferenceFEs/BDMRefFEs.jl b/src/ReferenceFEs/BDMRefFEs.jl index 984bb0bba..8fee9e4eb 100644 --- a/src/ReferenceFEs/BDMRefFEs.jl +++ b/src/ReferenceFEs/BDMRefFEs.jl @@ -1,7 +1,9 @@ -struct BDM <: DivConforming end +struct BDM <: PushforwardRefFE end const bdm = BDM() +Pushforward(::Type{<:BDM}) = ContraVariantPiolaMap() + """ BDMRefFE(::Type{et},p::Polytope,order::Integer) where et diff --git a/src/ReferenceFEs/Dofs.jl b/src/ReferenceFEs/Dofs.jl index cee46b14e..2e8d42b07 100644 --- a/src/ReferenceFEs/Dofs.jl +++ b/src/ReferenceFEs/Dofs.jl @@ -1,5 +1,107 @@ abstract type Dof <: Map end +""" + struct LinearCombinationDofVector{T<:Dof,V,F} <: AbstractVector{T} + values :: V + dofs :: F + end + +Type that implements a dof basis (a) as the linear combination of a dof basis +(b). The dofs are first evaluated at dof basis (b) (field `dofs`) and the +dof values are next mapped to dof basis (a) applying a change of basis (field +`values`). + +Fields: + +- `values::AbstractMatrix{<:Number}` the matrix of the change from dof basis (b) to (a) +- `dofs::AbstractVector{T}` A type representing dof basis (b), with `T<:Dof` +""" +struct LinearCombinationDofVector{T,V,F} <: AbstractVector{T} + values::V + dofs::F + function LinearCombinationDofVector( + values::AbstractMatrix{<:Number}, + dofs::AbstractVector{<:Dof} + ) + T = eltype(dofs) + V = typeof(values) + F = typeof(dofs) + new{T,V,F}(values,dofs) + end +end + +Base.size(a::LinearCombinationDofVector) = size(a.values,2) +Base.IndexStyle(::LinearCombinationDofVector) = IndexLinear() +Base.getindex(::LinearCombinationDofVector{T},::Integer) where T = T() + +function linear_combination(a::AbstractMatrix{<:Number}, b::AbstractVector{<:Dof}) + LinearCombinationDofVector(a,b) +end + +function return_cache(b::LinearCombinationDofVector,field) + k = Fields.LinearCombinationMap(:) + cf = return_cache(b.dofs,field) + fx = evaluate!(cf,b.dofs,field) + ck = return_cache(k,b.values,fx) + return cf, ck +end + +function evaluate!(cache,b::LinearCombinationDofVector,field) + cf, ck = cache + k = Fields.LinearCombinationMap(:) + fx = evaluate!(cf,b.dofs,field) + return evaluate!(ck,k,b.values,fx) +end + +""" + struct MappedDofBasis{T<:Dof,MT,BT} <: AbstractVector{T} + F :: MT + σ :: BT + args + end + +Represents η = σ∘F, evaluated as η(φ) = σ(F(φ,args...)) + + - σ : V* -> R is a dof basis + - F : W -> V is a map between function spaces + +Intended combinations would be: + +- σ : V* -> R dof basis in the physical domain and F* : V̂ -> V is a pushforward map. +- ̂σ : V̂* -> R dof basis in the reference domain and (F*)^-1 : V -> V̂ is an inverse pushforward map. + +""" +struct MappedDofBasis{T<:Dof,MT,BT,A} <: AbstractVector{T} + F :: MT + dofs :: BT + args :: A + + function MappedDofBasis(F::Map, dofs::AbstractVector{<:Dof}, args...) + T = eltype(dofs) + MT = typeof(F) + BT = typeof(dofs) + A = typeof(args) + new{T,MT,BT,A}(F,dofs,args) + end +end + +Base.size(b::MappedDofBasis) = size(b.dofs) +Base.IndexStyle(::MappedDofBasis) = IndexLinear() +Base.getindex(b::MappedDofBasis, i) = T() + +function Arrays.return_cache(b::MappedDofBasis, fields) + f_cache = return_cache(b.F,fields,b.args...) + ffields = evaluate!(f_cache,b.F,fields,b.args...) + dofs_cache = return_cache(b.dofs,ffields) + return f_cache, dofs_cache +end + +function Arrays.evaluate!(cache, b::MappedDofBasis, fields) + f_cache, dofs_cache = cache + ffields = evaluate!(f_cache,b.F,fields,b.args...) + evaluate!(dofs_cache,b.dofs,ffields) +end + """ test_dof(dof,field,v;cmp::Function=(==)) diff --git a/src/ReferenceFEs/HHJRefFEs.jl b/src/ReferenceFEs/HHJRefFEs.jl new file mode 100644 index 000000000..9c0c691e5 --- /dev/null +++ b/src/ReferenceFEs/HHJRefFEs.jl @@ -0,0 +1,70 @@ + +struct HellanHerrmannJhonson <: PushforwardRefFE end + +const hhj = HellanHerrmannJhonson() + +Pushforward(::Type{<:HellanHerrmannJhonson}) = DoubleContraVariantPiolaMap() + +""" + struct HellanHerrmannJhonson <: PushforwardRefFE end + HellanHerrmannJhonsonRefFE(::Type{T},p::Polytope,order::Integer) where T + +Hellan-Herrmann-Jhonson reference finite element. + +References: + +- `The Hellan-Herrmann-Johnson method with curved elements`, Arnold and Walker (2020) + +""" +function HellanHerrmannJhonsonRefFE(::Type{T},p::Polytope,order::Integer) where T + @assert p == TRI "HellanHerrmannJhonson Reference FE only defined for TRIangles" + + VT = SymTensorValue{2,T} + prebasis = MonomialBasis{2}(VT,order,Polynomials._p_filter) + fb = MonomialBasis{D-1}(T,order,Polynomials._p_filter) + cb = MonomialBasis{2}(VT,order-1,Polynomials._p_filter) + + function cmom(φ,μ,ds) # Cell and Node moment function: σ_K(φ,μ) = ∫(φ:μ)dK + Broadcasting(Operation(⊙))(φ,μ) + end + function fmom(φ,μ,ds) # Face moment function (normal) : σ_F(φ,μ) = ∫((n·φ·n)*μ)dF + n = get_facet_normal(ds) + φn = Broadcasting(Operation(⋅))(φ,n) + nφn = Broadcasting(Operation(⋅))(n,φn) + Broadcasting(Operation(*))(nφn,μ) + end + + moments = Tuple[ + (get_dimrange(p,1),fmom,fb), # Face moments + (get_dimrange(p,2),cmom,cb) # Cell moments + ] + + return MomentBasedReferenceFE(HellanHerrmannJhonson(),p,prebasis,moments,DivConformity()) +end + +function ReferenceFE(p::Polytope,::HellanHerrmannJhonson, order) + HellanHerrmannJhonsonRefFE(Float64,p,order) +end + +function ReferenceFE(p::Polytope,::HellanHerrmannJhonson,::Type{T}, order) where T + HellanHerrmannJhonsonRefFE(T,p,order) +end + +function Conformity(reffe::GenericRefFE{HellanHerrmannJhonson},sym::Symbol) + hdiv = (:Hdiv,:HDiv) + if sym == :L2 + L2Conformity() + elseif sym in hdiv + DivConformity() + else + @unreachable """\n + It is not possible to use conformity = $sym on a HellanHerrmannJhonson reference FE. + + Possible values of conformity for this reference fe are $((:L2, hdiv...)). + """ + end +end + +function get_face_own_dofs(reffe::GenericRefFE{HellanHerrmannJhonson}, conf::DivConformity) + get_face_dofs(reffe) +end diff --git a/src/ReferenceFEs/LinearCombinationDofVectors.jl b/src/ReferenceFEs/LinearCombinationDofVectors.jl deleted file mode 100644 index 507a34724..000000000 --- a/src/ReferenceFEs/LinearCombinationDofVectors.jl +++ /dev/null @@ -1,46 +0,0 @@ -""" - struct LinearCombinationDofVector{T} <: AbstractVector{Dof} - change_of_basis::Matrix{T} - dof_basis::AbstractVector{<:Dof} - end - -Type that implements a dof basis (a) as the linear combination of a dof basis -(b). The dofs are first evaluated at dof basis (b) (field `dof_basis`) and the -dof values are next mapped to dof basis (a) applying a change of basis (field -`change_of_basis`). - -Fields: - -- `change_of_basis::Matrix{T}` the matrix of the change from dof basis (b) to (a) -- `dof_basis::AbstractVector{<:Dof}` A type representing dof basis (b) -""" -struct LinearCombinationDofVector{T} <: AbstractVector{Dof} - change_of_basis::Matrix{T} - dof_basis::AbstractVector{<:Dof} -end - -@inline Base.size(a::LinearCombinationDofVector) = size(a.dof_basis) -@inline Base.axes(a::LinearCombinationDofVector) = axes(a.dof_basis) -@inline Base.getindex(a::LinearCombinationDofVector,i::Integer) = getindex(a.dof_basis,i) -@inline Base.IndexStyle(::LinearCombinationDofVector) = IndexLinear() - -function linear_combination(a::AbstractMatrix{<:Number}, - b::AbstractVector{<:Dof}) - LinearCombinationDofVector(a,b) -end - -function linear_combination(a::LinearCombinationDofVector{T}, - b::AbstractVector{<:Dof}) where T - linear_combination(a.change_of_basis,b) -end - -function return_cache(b::LinearCombinationDofVector,field) - c, cf = return_cache(b.dof_basis,field) - c, cf, return_cache(*,b.change_of_basis,c) -end - -@inline function evaluate!(cache,b::LinearCombinationDofVector,field) - c, cf, cc = cache - vals = evaluate!(cache,b.dof_basis,field) - evaluate!(cc,*,b.change_of_basis,vals) -end \ No newline at end of file diff --git a/src/ReferenceFEs/MTWRefFEs.jl b/src/ReferenceFEs/MTWRefFEs.jl index c334c10a1..d2f3ede6e 100644 --- a/src/ReferenceFEs/MTWRefFEs.jl +++ b/src/ReferenceFEs/MTWRefFEs.jl @@ -1,9 +1,11 @@ -struct MardalTaiWinther <: ReferenceFEName end + +struct MardalTaiWinther <: PushforwardRefFE end const mtw = MardalTaiWinther() """ -MardalTaiWintherRefFE(::Type{et},p::Polytope,order::Integer) where et + struct MardalTaiWinther <: PushforwardRefFE end + MardalTaiWintherRefFE(::Type{et},p::Polytope,order::Integer) where et Mardal-Tai-Winther reference finite element. diff --git a/src/ReferenceFEs/MomentBasedReferenceFEs.jl b/src/ReferenceFEs/MomentBasedReferenceFEs.jl index 27e551928..72abb3af9 100644 --- a/src/ReferenceFEs/MomentBasedReferenceFEs.jl +++ b/src/ReferenceFEs/MomentBasedReferenceFEs.jl @@ -36,7 +36,7 @@ struct MomentBasedDofBasis{P,V} <: AbstractVector{Moment} end Base.size(a::MomentBasedDofBasis) = (num_dofs(a),) -Base.axes(a::MomentBasedDofBasis) = (Base.OneTo(num_dofs),) +Base.axes(a::MomentBasedDofBasis) = (Base.OneTo(num_dofs(a)),) Base.getindex(a::MomentBasedDofBasis,i::Integer) = Moment() Base.IndexStyle(::MomentBasedDofBasis) = IndexLinear() diff --git a/src/ReferenceFEs/NedelecRefFEs.jl b/src/ReferenceFEs/NedelecRefFEs.jl index 3faa07361..f8bb0aef4 100644 --- a/src/ReferenceFEs/NedelecRefFEs.jl +++ b/src/ReferenceFEs/NedelecRefFEs.jl @@ -1,9 +1,11 @@ struct CurlConformity <: Conformity end -struct Nedelec <: ReferenceFEName end +struct Nedelec <: PushforwardRefFE end const nedelec = Nedelec() +Pushforward(::Type{<:Nedelec}) = CoVariantPiolaMap() + """ NedelecRefFE(::Type{et},p::Polytope,order::Integer) where et @@ -115,43 +117,3 @@ function get_face_dofs(reffe::GenericRefFE{Nedelec,Dc}) where Dc end face_dofs end - -struct CoVariantPiolaMap <: Map end - -function evaluate!( - cache, - ::Broadcasting{typeof(∇)}, - a::Fields.BroadcastOpFieldArray{CoVariantPiolaMap}) - v, Jt = a.args - # Assuming J comes from an affine map - ∇v = Broadcasting(∇)(v) - k = CoVariantPiolaMap() - Broadcasting(Operation(k))(∇v,Jt) -end - -function lazy_map( - ::Broadcasting{typeof(gradient)}, - a::LazyArray{<:Fill{Broadcasting{Operation{CoVariantPiolaMap}}}}) - v, Jt = a.args - ∇v = lazy_map(Broadcasting(∇),v) - k = CoVariantPiolaMap() - lazy_map(Broadcasting(Operation(k)),∇v,Jt) -end - -function evaluate!(cache,::CoVariantPiolaMap,v::Number,Jt::Number) - v⋅transpose(inv(Jt)) # we multiply by the right side to compute the gradient correctly -end - -function evaluate!(cache,k::CoVariantPiolaMap,v::AbstractVector{<:Field},phi::Field) - Jt = ∇(phi) - Broadcasting(Operation(k))(v,Jt) -end - -function lazy_map( - k::CoVariantPiolaMap, - cell_ref_shapefuns::AbstractArray{<:AbstractArray{<:Field}}, - cell_map::AbstractArray{<:Field}) - - cell_Jt = lazy_map(∇,cell_map) - lazy_map(Broadcasting(Operation(k)),cell_ref_shapefuns,cell_Jt) -end diff --git a/src/ReferenceFEs/Pullbacks.jl b/src/ReferenceFEs/Pullbacks.jl index 9c6477238..1a825b460 100644 --- a/src/ReferenceFEs/Pullbacks.jl +++ b/src/ReferenceFEs/Pullbacks.jl @@ -10,8 +10,13 @@ where """ abstract type Pushforward <: Map end +abstract type PushforwardRefFE <: ReferenceFEName end + +Pushforward(::Type{<:PushforwardRefFE}) = @abstractmethod +Pushforward(name::PushforwardRefFE) = Pushforward(typeof(name)) + function Arrays.lazy_map( - k::Pushforward, ref_cell_fields, pf_args... + k::Pushforward, ref_cell_fields::AbstractArray, pf_args::AbstractArray... ) lazy_map(Broadcasting(Operation(k)), ref_cell_fields, pf_args...) end @@ -22,18 +27,34 @@ function Arrays.evaluate!( @abstractmethod end +function evaluate!( + cache, k::Pushforward, f_ref::AbstractVector{<:Field}, args... +) + Broadcasting(Operation(k))(f_ref,args...) +end + function Arrays.lazy_map( ::Broadcasting{typeof(gradient)}, a::LazyArray{<:Fill{Broadcasting{Operation{<:Pushforward}}}} ) - cell_ref_fields, args = a.args + cell_ref_fields, args... = a.args cell_ref_gradient = lazy_map(Broadcasting(∇),cell_ref_fields) return lazy_map(a.maps.value,cell_ref_gradient,args...) end +function Arrays.evaluate!( + cache, + ::Broadcasting{typeof(gradient)}, + a::Fields.BroadcastOpFieldArray{<:Pushforward} +) + v, pf_args... = a.args + grad_v = Broadcasting(∇)(v) + Broadcasting(Operation(a.op))(grad_v,pf_args...) +end + # InversePushforward """ - struct InversePushforward{PF <: Pushforward} <: Map end + const InversePushforward{PF} = InverseMap{PF} where PF <: Pushforward Represents the inverse of a pushforward map F*, defined as (F*)^-1 : V -> V̂ @@ -41,40 +62,27 @@ where - V̂ is a function space on the reference cell K̂ and - V is a function space on the physical cell K. """ -struct InversePushforward{PF} <: Map - pushforward::PF - function InversePushforward(pushforward::Pushforward) - PF = typeof(pushforward) - new{PF}(pushforward) - end -end - -Arrays.inverse_map(pf::Pushforward) = InversePushforward(pf) -Arrays.inverse_map(ipf::InversePushforward) = ipf.pushforward +const InversePushforward{PF} = InverseMap{PF} where PF <: Pushforward function Arrays.lazy_map( - k::InversePushforward, phys_cell_fields, pf_args... + k::InversePushforward, phys_cell_fields::AbstractArray, pf_args::AbstractArray... ) lazy_map(Broadcasting(Operation(k)), phys_cell_fields, pf_args...) end -function Arrays.return_cache( - k::InversePushforward, v_phys::Number, args... +function evaluate!( + cache, k::InversePushforward, f_phys::AbstractVector{<:Field}, args... ) - v_ref_basis = mock_basis(v_phys) - pf_cache = return_cache(k.pushforward,v_ref_basis,args...) - return v_ref_basis, pf_cache + Broadcasting(Operation(k))(f_phys,args...) end -function Arrays.evaluate!( - cache, k::InversePushforward, v_phys::Number, args... +function evaluate!( + cache, k::InversePushforward, f_phys::Field, args... ) - v_ref_basis, pf_cache = cache - change = evaluate!(pf_cache,k.pushforward,v_ref_basis,args...) - return inv(change)⋅v_phys + Operation(k)(f_phys,args...) end -# Pushforward +# Pullback """ struct Pullback{PF <: Pushforward} <: Map end @@ -87,23 +95,25 @@ where Its action on physical dofs σ : V -> R is defined in terms of the pushforward map F* as ̂σ = F**(σ) := σ∘F* : V̂ -> R """ -struct Pullback{PF} <: Map +struct Pullback{PF <: Pushforward} <: Map pushforward::PF - function Pullback(pushforward::Pushforward) - PF = typeof(pushforward) - new{PF}(pushforward) - end end function Arrays.lazy_map( - ::typeof{evaluate},k::LazyArray{<:Fill{<:Pullback}},ref_cell_fields + ::typeof(evaluate),k::LazyArray{<:Fill{<:Pullback}},ref_cell_fields::AbstractArray ) pb = k.maps.value - phys_cell_dofs, pf_args = k.args + phys_cell_dofs, pf_args... = k.args phys_cell_fields = lazy_map(pb.pushforward,ref_cell_fields,pf_args...) return lazy_map(evaluate,phys_cell_dofs,phys_cell_fields) end +function evaluate!( + cache, k::Pullback, σ_phys::AbstractVector{<:Dof}, args... +) + return MappedDofBasis(k.pushforward,σ_phys,args...) +end + # InversePullback """ @@ -115,38 +125,115 @@ where - V̂* is a dof space on the reference cell K̂ and - V* is a dof space on the physical cell K. Its action on reference dofs ̂σ : V̂ -> R is defined in terms of the pushforward map F* as - σ = F**(̂σ) := ̂σ∘(F*)^-1 : V -> R + σ = (F**)^-1(̂σ) := ̂σ∘(F*)^-1 : V -> R """ -struct InversePullback{PF} <: Map - pushforward::PF - function InversePullback(pushforward::Pushforward) - PF = typeof(pushforward) - new{PF}(pushforward) - end -end - -Arrays.inverse_map(pb::Pullback) = InversePullback(pb.pushforward) -Arrays.inverse_map(ipb::InversePullback) = Pullback(ipb.pushforward) +const InversePullback{PB} = InverseMap{PB} where PB <: Pullback function Arrays.lazy_map( - ::typeof{evaluate},k::LazyArray{<:Fill{<:InversePullback}},phys_cell_fields + ::typeof(evaluate), k::LazyArray{<:Fill{<:InversePullback}}, phys_cell_fields::AbstractArray ) pb = inverse_map(k.maps.value) - ref_cell_dofs, pf_args = k.args - ref_cell_fields = lazy_map(inverse_map(pb.pushforward),phys_cell_fields,pf_args...) + ref_cell_dofs, pf_args... = k.args + ref_cell_fields = lazy_map(inverse_map(pb.pushforward), phys_cell_fields, pf_args...) return lazy_map(evaluate,ref_cell_dofs,ref_cell_fields) end +function evaluate!( + cache, k::InversePullback, σ_ref::AbstractVector{<:Dof}, args... +) + pb = inverse_map(k) + return MappedDofBasis(inverse_map(pb.pushforward),σ_ref,args...) +end + # ContraVariantPiolaMap struct ContraVariantPiolaMap <: Pushforward end function evaluate!( - cache,::ContraVariantPiolaMap, - v::Number, - Jt::Number, - sign_flip::Bool + cache, ::ContraVariantPiolaMap, v_ref::Number, Jt::Number +) + idetJ = 1. / meas(Jt) + return v_ref ⋅ (idetJ * Jt) +end + +function evaluate!( + cache, ::InversePushforward{ContraVariantPiolaMap}, v_phys::Number, Jt::Number +) + detJ = meas(Jt) + return v_phys ⋅ (detJ * pinvJt(Jt)) +end + +# TODO: Should this be here? Probably not... + +function Fields.DIV(f::LazyArray{<:Fill}) + df = Fields.DIV(f.args[1]) + k = f.maps.value + lazy_map(k,df) +end + +function Fields.DIV(a::LazyArray{<:Fill{typeof(linear_combination)}}) + i_to_basis = Fields.DIV(a.args[2]) + i_to_values = a.args[1] + lazy_map(linear_combination,i_to_values,i_to_basis) +end + +function Fields.DIV(f::LazyArray{<:Fill{Broadcasting{Operation{ContraVariantPiolaMap}}}}) + ϕrgₖ = f.args[1] + return lazy_map(Broadcasting(divergence),ϕrgₖ) +end + +function Fields.DIV(f::Fill{<:Fields.BroadcastOpFieldArray{ContraVariantPiolaMap}}) + ϕrgₖ = f.value.args[1] + return Fill(Broadcasting(divergence)(ϕrgₖ),length(f)) +end + +# CoVariantPiolaMap + +struct CoVariantPiolaMap <: Pushforward end + +function evaluate!( + cache, ::CoVariantPiolaMap, v_ref::Number, Jt::Number +) + return v_ref ⋅ transpose(pinvJt(Jt)) +end + +function evaluate!( + cache, ::InversePushforward{CoVariantPiolaMap}, v_phys::Number, Jt::Number +) + return v_phys ⋅ transpose(Jt) +end + +# DoubleContraVariantPiolaMap + +struct DoubleContraVariantPiolaMap <: Pushforward end + +function evaluate!( + cache, ::DoubleContraVariantPiolaMap, v_ref::Number, Jt::Number +) + _Jt = meas(Jt) * Jt + return transpose(_Jt) ⋅ v_ref ⋅ _Jt +end + +function evaluate!( + cache, ::InversePushforward{DoubleContraVariantPiolaMap}, v_phys::Number, Jt::Number +) + iJt = meas(Jt) * pinvJt(Jt) + return transpose(iJt) ⋅ v_phys ⋅ iJt +end + +# DoubleCoVariantPiolaMap + +struct DoubleCoVariantPiolaMap <: Pushforward end + +function evaluate!( + cache, ::DoubleCoVariantPiolaMap, v_ref::Number, Jt::Number +) + iJt = pinvJt(Jt) + return iJt ⋅ v_ref ⋅ transpose(iJt) +end + +function evaluate!( + cache, ::InversePushforward{DoubleCoVariantPiolaMap}, v_phys::Number, Jt::Number ) - idetJ = 1/meas(Jt) - ((-1)^sign_flip*v)⋅(idetJ*Jt) + return Jt ⋅ v_phys ⋅ transpose(Jt) end diff --git a/src/ReferenceFEs/RaviartThomasRefFEs.jl b/src/ReferenceFEs/RaviartThomasRefFEs.jl index 29844461b..4920de0c4 100644 --- a/src/ReferenceFEs/RaviartThomasRefFEs.jl +++ b/src/ReferenceFEs/RaviartThomasRefFEs.jl @@ -1,12 +1,13 @@ struct DivConformity <: Conformity end -abstract type DivConforming <: ReferenceFEName end # RaviartThomas -struct RaviartThomas <: DivConforming end +struct RaviartThomas <: PushforwardRefFE end const raviart_thomas = RaviartThomas() +Pushforward(::Type{<:RaviartThomas}) = ContraVariantPiolaMap() + """ RaviartThomasRefFE(::Type{et},p::Polytope,order::Integer) where et @@ -88,59 +89,3 @@ function compute_jacobi_basis(::Type{T},p::ExtrusionPolytope{D},orders) where {D terms = _monomial_terms(extrusion,orders) JacobiPolynomialBasis{D}(T,orders,terms) end - -# ContraVariantPiolaMap - -struct ContraVariantPiolaMap <: Map end - -function evaluate!( - cache, - ::Broadcasting{typeof(∇)}, - a::Fields.BroadcastOpFieldArray{ContraVariantPiolaMap} -) - v, Jt, sign_flip = a.args - ∇v = Broadcasting(∇)(v) - k = ContraVariantPiolaMap() - Broadcasting(Operation(k))(∇v,Jt,sign_flip) -end - -function lazy_map( - ::Broadcasting{typeof(gradient)}, - a::LazyArray{<:Fill{Broadcasting{Operation{ContraVariantPiolaMap}}}} -) - v, Jt, sign_flip = a.args - ∇v = lazy_map(Broadcasting(∇),v) - k = ContraVariantPiolaMap() - lazy_map(Broadcasting(Operation(k)),∇v,Jt,sign_flip) -end - -function lazy_map( - k::ContraVariantPiolaMap, - cell_ref_shapefuns::AbstractArray{<:AbstractArray{<:Field}}, - cell_map::AbstractArray{<:Field}, - sign_flip::AbstractArray{<:AbstractArray{<:Field}} -) - cell_Jt = lazy_map(∇,cell_map) - lazy_map(Broadcasting(Operation(k)),cell_ref_shapefuns,cell_Jt,sign_flip) -end - -function evaluate!( - cache,::ContraVariantPiolaMap, - v::Number, - Jt::Number, - sign_flip::Bool -) - idetJ = 1/meas(Jt) - ((-1)^sign_flip*v)⋅(idetJ*Jt) -end - -function evaluate!( - cache, - k::ContraVariantPiolaMap, - v::AbstractVector{<:Field}, - phi::Field, - sign_flip::AbstractVector{<:Field} -) - Jt = ∇(phi) - Broadcasting(Operation(k))(v,Jt,sign_flip) -end diff --git a/src/ReferenceFEs/ReferenceFEInterfaces.jl b/src/ReferenceFEs/ReferenceFEInterfaces.jl index ead616c1f..a4b1c0b82 100644 --- a/src/ReferenceFEs/ReferenceFEInterfaces.jl +++ b/src/ReferenceFEs/ReferenceFEInterfaces.jl @@ -390,11 +390,11 @@ associated with the dof basis `predofs` and the basis `shapefuns`. It is equivalent to - change = inv(evaluate(predofs,shapefuns)) + change = transpose(inv(evaluate(predofs,shapefuns))) linear_combination(change,predofs) # i.e. transpose(change)*predofs """ function compute_dofs(predofs,shapefuns) - change = inv(evaluate(predofs,shapefuns)) + change = transpose(inv(evaluate(predofs,shapefuns))) linear_combination(change,predofs) end diff --git a/src/ReferenceFEs/ReferenceFEs.jl b/src/ReferenceFEs/ReferenceFEs.jl index abfc4766d..a433e49ba 100644 --- a/src/ReferenceFEs/ReferenceFEs.jl +++ b/src/ReferenceFEs/ReferenceFEs.jl @@ -186,7 +186,6 @@ export ModalC0RefFE export CRRefFE export Lagrangian -export DivConforming export RaviartThomas export BDM export Nedelec @@ -248,6 +247,8 @@ include("StrangQuadratures.jl") include("XiaoGimbutasQuadratures.jl") +include("Pullbacks.jl") + include("MomentBasedReferenceFEs.jl") include("RaviartThomasRefFEs.jl") @@ -264,6 +265,4 @@ include("BezierRefFEs.jl") include("ModalC0RefFEs.jl") -include("LinearCombinationDofVectors.jl") - end # module diff --git a/test/AdaptivityTests/CartesianRefinementTests.jl b/test/AdaptivityTests/CartesianRefinementTests.jl index 1f91dffce..8c9b6675f 100644 --- a/test/AdaptivityTests/CartesianRefinementTests.jl +++ b/test/AdaptivityTests/CartesianRefinementTests.jl @@ -40,8 +40,6 @@ function test_grid_transfers(D,model,parent,order) U_f = TrialFESpace(V_f,sol) V_c = TestFESpace(parent,reffe;dirichlet_tags="boundary") U_c = TrialFESpace(V_c,sol) - V_c_fast = TestFESpace(parent,rrules,reffe;dirichlet_tags="boundary") - U_c_fast = TrialFESpace(V_c_fast,sol) # CellField: Coarse -> Fine cf_c_phy = CellField(sol,ctrian) @@ -80,14 +78,12 @@ function test_grid_transfers(D,model,parent,order) uh_c_inter = interpolate(uh_f,U_c) uh_c_inter2 = interpolate_everywhere(uh_f,U_c) uh_c_inter3 = interpolate_dirichlet(uh_f,U_c) - uh_c_inter4 = interpolate(uh_f,U_c_fast) @test l2_error(uh_f,uh_f_inter,dΩ_f) < 1.e-8 @test l2_error(uh_f,uh_f_inter2,dΩ_f) < 1.e-8 @test l2_error(uh_c,uh_c_inter,dΩ_c) < 1.e-8 @test l2_error(uh_c,uh_c_inter2,dΩ_c) < 1.e-8 - @test l2_error(uh_c,uh_c_inter4,dΩ_c) < 1.e-8 # Coarse FEFunction -> Fine FEFunction, by projection af(u,v) = ∫(v⋅u)*dΩ_f diff --git a/test/AdaptivityTests/EdgeBasedRefinementTests.jl b/test/AdaptivityTests/EdgeBasedRefinementTests.jl index d68ec1cad..929b26556 100644 --- a/test/AdaptivityTests/EdgeBasedRefinementTests.jl +++ b/test/AdaptivityTests/EdgeBasedRefinementTests.jl @@ -38,8 +38,6 @@ function test_grid_transfers(parent,model,order) U_f = TrialFESpace(V_f,sol) V_c = TestFESpace(parent,reffe;dirichlet_tags="boundary") U_c = TrialFESpace(V_c,sol) - V_c_fast = TestFESpace(parent,rrules,reffe;dirichlet_tags="boundary") - U_c_fast = TrialFESpace(V_c_fast,sol) # CellField: Coarse -> Fine cf_c_phy = CellField(sol,ctrian) @@ -78,14 +76,12 @@ function test_grid_transfers(parent,model,order) uh_c_inter = interpolate(uh_f,U_c) uh_c_inter2 = interpolate_everywhere(uh_f,U_c) uh_c_inter3 = interpolate_dirichlet(uh_f,U_c) - uh_c_inter4 = interpolate(uh_f,U_c_fast) @test l2_error(uh_f,uh_f_inter,dΩ_f) < 1.e-8 @test l2_error(uh_f,uh_f_inter2,dΩ_f) < 1.e-8 @test l2_error(uh_c,uh_c_inter,dΩ_c) < 1.e-8 @test l2_error(uh_c,uh_c_inter2,dΩ_c) < 1.e-8 - @test l2_error(uh_c,uh_c_inter4,dΩ_c) < 1.e-8 # Coarse FEFunction -> Fine FEFunction, by projection af(u,v) = ∫(v⋅u)*dΩ_f diff --git a/test/FESpacesTests/ConstantFESpaceTests.jl b/test/FESpacesTests/ConstantFESpaceTests.jl index 980ee6124..cc81846f7 100644 --- a/test/FESpacesTests/ConstantFESpaceTests.jl +++ b/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 diff --git a/test/FESpacesTests/DivConformingFESpacesTests.jl b/test/FESpacesTests/DivConformingFESpacesTests.jl index 39ac66419..46dbd805d 100644 --- a/test/FESpacesTests/DivConformingFESpacesTests.jl +++ b/test/FESpacesTests/DivConformingFESpacesTests.jl @@ -10,19 +10,19 @@ using Gridap.Fields using Gridap.Io function test_div_v_q_equiv(U,V,P,Q,Ω) - v=get_fe_basis(V) - u=get_trial_fe_basis(U) + v = get_fe_basis(V) + u = get_trial_fe_basis(U) - q=get_fe_basis(Q) - p=get_trial_fe_basis(P) + q = get_fe_basis(Q) + p = get_trial_fe_basis(P) - dΩ=Measure(Ω,1) - dΩᵣ=Measure(Ω,1,integration_domain_style=ReferenceDomain()) + dΩ = Measure(Ω,1) + dΩᵣ = Measure(Ω,1,integration_domain_style = ReferenceDomain()) - a1(p,v)=∫(divergence(v)*p)dΩ - a2(p,v)=∫(DIV(v)*p)dΩᵣ + a1(p,v) = ∫(divergence(v)*p)dΩ + a2(p,v) = ∫(DIV(v)*p)dΩᵣ - tol=1.0e-12 + tol = 1.0e-12 assem = SparseMatrixAssembler(P,V) data = collect_cell_matrix(P,V,a1(p,v)) A1 = assemble_matrix(assem,data) @@ -30,8 +30,8 @@ function test_div_v_q_equiv(U,V,P,Q,Ω) A2 = assemble_matrix(assem,data) @test norm(A1-A2) < tol - a3(u,q)=∫(q*divergence(u))dΩ - a4(u,q)=∫(q*DIV(u))dΩᵣ + a3(u,q) = ∫(q*divergence(u))dΩ + a4(u,q) = ∫(q*DIV(u))dΩᵣ assem = SparseMatrixAssembler(U,Q) data = collect_cell_matrix(U,Q,a3(u,q)) A3 = assemble_matrix(assem,data) @@ -39,36 +39,33 @@ function test_div_v_q_equiv(U,V,P,Q,Ω) A4 = assemble_matrix(assem,data) @test norm(A3-A4) < tol - uh=FEFunction(U,rand(num_free_dofs(U))) - l1(q)=∫(q*divergence(uh))dΩ - l2(q)=∫(q*DIV(uh))dΩᵣ - v1=assemble_vector(l1,Q) - v2=assemble_vector(l2,Q) + uh = FEFunction(U,rand(num_free_dofs(U))) + l1(q) = ∫(q*divergence(uh))dΩ + l2(q) = ∫(q*DIV(uh))dΩᵣ + v1 = assemble_vector(l1,Q) + v2 = assemble_vector(l2,Q) @test norm(v1-v2) < tol end -@testset "Test Raviart-Thomas" begin +#@testset "Test Raviart-Thomas" begin - domain =(0,1,0,1) + domain = (0,1,0,1) partition = (3,3) model = CartesianDiscreteModel(domain,partition) order = 1 - u(x) = x reffe = ReferenceFE(raviart_thomas,order) - V = TestFESpace(model,reffe,dirichlet_tags = [1,6]) test_single_field_fe_space(V) U = TrialFESpace(V,u) - reffe = ReferenceFE(lagrangian,Float64,order) - Q = TestFESpace(model,reffe,conformity=:L2) + reffe_p = ReferenceFE(lagrangian,Float64,order) + Q = TestFESpace(model,reffe_p,conformity=:L2) P = TrialFESpace(Q) uh = interpolate(u,U) - e = u - uh Ω = Triangulation(model) @@ -103,7 +100,6 @@ end v(x) = VectorValue(-0.5*x[1]+1.0,-0.5*x[2],-0.5*x[3]) vh = interpolate(v,V) - e = v - vh Ω = Triangulation(model) @@ -150,7 +146,7 @@ end test_div_v_q_equiv(U,V,P,Q,Ω) -end +#end @testset "Test BDM" begin diff --git a/test/FieldsTests/AffineMapsTests.jl b/test/FieldsTests/AffineMapsTests.jl index 284aaec1d..963e3faf6 100644 --- a/test/FieldsTests/AffineMapsTests.jl +++ b/test/FieldsTests/AffineMapsTests.jl @@ -9,7 +9,7 @@ using Test origin = Point(1,1) g = TensorValue(2,0,0,2) -h = AffineMap(g,origin) +h = AffineField(g,origin) @test isa(∇(h),ConstantField) @test isa(Broadcasting(∇)(h),ConstantField) @@ -39,7 +39,7 @@ cell_to_∇hx = lazy_map(evaluate,cell_to_∇h,cell_to_x) test_array(cell_to_hx,fill(hx,ncells)) test_array(cell_to_∇hx,fill(∇hx,ncells)) -T = AffineMap{3,3,Int} +T = AffineField{3,3,Int} @test isa(zero(T),T) #display(cell_to_hx) diff --git a/test/FieldsTests/InverseFieldsTests.jl b/test/FieldsTests/InverseFieldsTests.jl index 8cf7bb5d2..f8db6750a 100644 --- a/test/FieldsTests/InverseFieldsTests.jl +++ b/test/FieldsTests/InverseFieldsTests.jl @@ -8,22 +8,22 @@ using Test b0 = Point(0,0) m0 = TensorValue(1,0,0,1) -id = AffineMap(m0,b0) +id = AffineField(m0,b0) b1 = Point(1,1) m1 = TensorValue(2,0,0,2) -h1 = AffineMap(m1,b1) +h1 = AffineField(m1,b1) b2 = Point(3,3) m2 = TensorValue(4,0,0,4) -h2 = AffineMap(m2,b2) +h2 = AffineField(m2,b2) h = h2 ∘ h1 # (4x+3) ∘ (2x+1) = 8x+7 b3 = Point(7,7) m3 = TensorValue(8,0,0,8) -h3 = AffineMap(m3,b3) +h3 = AffineField(m3,b3) h1inv = inverse_map(h1) h2inv = inverse_map(h2) diff --git a/test/GeometryTests/BoundaryTriangulationsTests.jl b/test/GeometryTests/BoundaryTriangulationsTests.jl index f876d2a86..b775798ef 100644 --- a/test/GeometryTests/BoundaryTriangulationsTests.jl +++ b/test/GeometryTests/BoundaryTriangulationsTests.jl @@ -78,7 +78,7 @@ glue = get_glue(btrian,Val(1)) glue = get_glue(btrian,Val(2)) @test glue.tface_to_mface === btrian.glue.face_to_cell -@test isa(get_cell_map(btrian)[1],AffineMap) +@test isa(get_cell_map(btrian)[1],AffineField) face_s_q = glue.tface_to_mface_map diff --git a/test/GridapTests/issue_879.jl b/test/GridapTests/issue_879.jl index d0eefea69..bd9613aad 100644 --- a/test/GridapTests/issue_879.jl +++ b/test/GridapTests/issue_879.jl @@ -38,6 +38,6 @@ fΓ = interpolate_everywhere(fa, FESpace(Γ,reffe)) # ERROR: LoadError: DimensionMismatch: matrix is not square: dimensions are (1, 2) # Corrections -# Modified src/Fields/AffineMaps.jl +# Modified src/Fields/AffineFields.jl end \ No newline at end of file diff --git a/test/pullbacks.jl b/test/pullbacks.jl new file mode 100644 index 000000000..716c68658 --- /dev/null +++ b/test/pullbacks.jl @@ -0,0 +1,53 @@ + +using FillArrays +using Gridap +using Gridap.ReferenceFEs, Gridap.FESpaces, Gridap.CellData +using Gridap.Fields + +using Gridap.ReferenceFEs: Pullback + +model = CartesianDiscreteModel((0,2,0,4),(2,2)) +Ω = Triangulation(model) +dΩ = Measure(Ω,2) +pts = CellData.get_data(get_cell_points(dΩ)) + +reffe = RaviartThomasRefFE(Float64,QUAD,1) +V = FESpace(model,reffe) + +u(x) = VectorValue(x[1], -x[2]) +uh = interpolate(u,V) + +φ_phys = get_fe_basis(V).cell_basis +φ_ref = φ_phys.args[1] +Jt, sign = φ_phys.args[2:end] + +σ_phys = get_fe_dof_basis(V).cell_dof +σ_ref = Fill(get_dof_basis(reffe),num_cells(model)) + +App = lazy_map(evaluate,σ_phys,φ_phys)[1] +Arr = lazy_map(evaluate,σ_ref,φ_ref)[1] + +pf = ContraVariantPiolaMap() + +# Inverse Pushforward +ipf = inverse_map(pf) +f_ref = lazy_map(ipf,φ_phys,Jt,sign) +Brr = lazy_map(evaluate,σ_ref,f_ref)[1] +Brr == Arr +φ_ref_x = lazy_map(evaluate,φ_ref,pts)[1] +f_ref_x = lazy_map(evaluate,f_ref,pts)[1] +f_ref_x ≈ φ_ref_x + +# Pullback +pb = Pullback(pf) +θ_ref = lazy_map(pb,σ_phys,Jt,sign) +θ_ref_x = lazy_map(evaluate,θ_ref,φ_ref) +σ_ref_x = lazy_map(evaluate,σ_ref,φ_ref) +θ_ref_x ≈ σ_ref_x + +# Inverse Pullback +ipb = inverse_map(pb) +θ_phys = lazy_map(ipb,σ_ref,Jt,sign) +θ_phys_x = lazy_map(evaluate,θ_phys,φ_phys) +σ_phys_x = lazy_map(evaluate,σ_phys,φ_phys) +θ_phys_x ≈ σ_phys_x