292 lines
20 KiB
TeX
292 lines
20 KiB
TeX
\documentclass{article}[twocolumn]
|
|
\usepackage{graphicx} % Required for inserting images
|
|
\usepackage{booktabs}
|
|
\usepackage{amsmath}
|
|
\usepackage{amssymb}
|
|
\usepackage{geometry}
|
|
\geometry{
|
|
letterpaper,
|
|
left=20mm,
|
|
top=20mm,
|
|
right=20mm,
|
|
bottom=30mm
|
|
}
|
|
|
|
\title{Deriving The Full Lane Emden Weak Form}
|
|
\author{Emily M. Boudreaux}
|
|
\date{April 2025}
|
|
|
|
\begin{document}
|
|
|
|
\maketitle
|
|
|
|
\section{Continuous Variational Form}
|
|
We start with the strong form of the Lane-Emden equation in three dimensions.
|
|
\begin{align}
|
|
\Delta \theta + \theta^{n} = 0
|
|
\end{align}
|
|
We put this into weak form by multiplying by some scalar test function $v^{\theta}$ which lives in the Sobolev space $H^{1}(\Omega)$
|
|
\begin{align}
|
|
\int_{\Omega}v^{\theta}\Delta\theta dV + \int_{\Omega}v^{\theta}\theta^{n}dV = 0
|
|
\end{align}
|
|
Applying Green's first identity
|
|
\begin{align}
|
|
\oint_{\partial \Omega}v^{\theta}\left(\nabla\theta\cdot\hat{n}\right)dA - \int_{\Omega}\nabla v^{\theta}\cdot\nabla\theta dV + \int_{\Omega}v^{\theta}\theta^{n} dV = 0
|
|
\end{align}
|
|
We let the surface integral go to zero as the value of $\theta$, and therefore $v_{\theta}$, at the surface of the domain is physically constrained to equal 0.
|
|
\begin{align}
|
|
- \int_{\Omega}\nabla v^{\theta}\cdot\nabla\theta dV + \int_{\Omega}v^{\theta}\theta^{n} dV = 0
|
|
\end{align}
|
|
Now we define a new variable $\phi\equiv\nabla\theta$ so that we can eventually apply essential boundary conditions to both $\theta$ and $\nabla\theta$. We must also then find the variational form of this expression. For that we multiply by some vector test function, $\vec{v}^{\phi}$ which will live in some vector space (In MFEM we will eventually use a Raviart-Thomas space, denoted $RT^{0}(\Omega)$)
|
|
\begin{align}
|
|
\int_{\Omega}\vec{v}^{\phi}\cdot\vec{\phi}dV - \int_{\Omega}\vec{v}^{\phi}\cdot\nabla\theta = 0
|
|
\end{align}
|
|
So then the final, continuous variational system of equations which we have is
|
|
\begin{align}
|
|
-\int_{\Omega}\nabla v^{\theta}\cdot\phi dV + \int_{\Omega}v^{\theta}\theta^{n}dV &= 0 \\
|
|
\int_{\Omega}\vec{v}^{\phi}\cdot\vec{\phi}dV - \int_{\Omega}\vec{v}^{\phi}\cdot\nabla\theta &= 0
|
|
\end{align}
|
|
\section{Discritized Variational Form}
|
|
In order to work with this in FEM we need to discritize this. First, Let $\theta_{h}$ be some discrete approximation of $\theta$ which lives on $v^{\theta}_{h}\subset v^{\theta}$ such that
|
|
\begin{align}
|
|
\theta_{h} = \sum_{i=1}^{N^{\theta}_{dof}}\theta_{i}N_{i}^{\theta}
|
|
\end{align}
|
|
Where $\{N_{i}^{\theta}\}_{i=1}^{N^{\theta}_{dof}}$ is a set of basis functions which span $v_{h}^{\theta}$ and $\theta_{i}$ are scalar degrees of freedom associated to each basis function. Similarly we can discritize $\phi$ by first letting $\phi_{h}$ be a discrete approximation of $\phi$ which live on $v_{h}^{\phi}\subset v^{\phi}$
|
|
\begin{align}
|
|
\phi_{h} = \sum_{j=1}^{N^{\phi}_{dof}}\phi_{j}N_{j}^{\phi}
|
|
\end{align}
|
|
where $\{N_{j}^{\phi}\}_{j=1}^{N_{dof}^{\phi}}$ is a set of vector basis functions which span $v_{h}^{\phi}$. Let us further let the column vectors $\bar{\theta}$ and $\bar{\phi}$ be
|
|
\begin{align}
|
|
\bar\theta &\equiv \left[\theta_{1},\ ...\ ,\theta_{N^{\theta}_{dof}}\right]^{T}\\
|
|
\bar\phi &\equiv \left[\phi_{1},\ ...\ ,\phi_{N^{\phi}_{dof}}\right]^{T}
|
|
\end{align}
|
|
|
|
In order to discritize the weak form we need to adopt a method of selecting test functions for $\theta$ and $\phi$. In the Galerkin method the test functions are chosen from the same finite dimensional subspaces which the approximate solutions are defined on. This is typically done by selecting each basis function to be a test function. This means then that we approximate
|
|
\begin{align}
|
|
v^{\theta}_{h} &= N_{k}^{\theta} \ \ \ \forall \ \ \ k=1,\ ...\ ,N^{\theta}_{dof} \\
|
|
v^{\phi}_{h} &= N_{\ell}^{\phi} \ \ \ \forall \ \ \ \ell=1,\ ...\ ,N^{\phi}_{dof}
|
|
\end{align}
|
|
We can now substitute these discritzed expressions for $\theta$, $\phi$, $v^{\theta}$, and $v^{\phi}$ back into the weak form...
|
|
\begin{align}
|
|
-\int_{\Omega}\nabla N_{k}^{\theta}\left(\sum_{j=1}^{N^{\phi}_{dof}}\phi_{j}\vec{N_{j}^{\phi}}\right)dV + \int_{\Omega}N_{k}^{\theta}\left(\sum_{i=1}^{N_{dof}^{\theta}}\theta_{i}N_{i}^{\theta}\right)^{n}dV &= 0 \\
|
|
\int_{\Omega}\vec{N}_{\ell}^{\phi}\cdot\left(\sum_{j=1}^{N^{\phi}_{dof}}\phi_{j}\vec{N_{j}^{\phi}}\right)dV - \int_{\Omega}\vec{N}_{\ell}^{\phi}\cdot\nabla\left(\sum_{i=1}^{N^{\theta}_{dof}}\theta_{i}N_{i}^{\theta}\right)dV &= 0
|
|
\end{align}
|
|
|
|
I want to pause here and make a note of a point of symbolics which might become confusing latter. We are going to be substituting the basis function, $N^{a}_{b}$, into various places in these equations. However, depending on if we substitute them in for the test functions, $v^{a}$, or the trial functions, $\theta$ and $\phi$, the semantic meaning of those basis functions changes. Any basis function set, $N^{a}_{b}$, used to represented a test function will eventually represent the \textbf{range} of the operator; whereas, any basis function set used to represent a trial function will eventually represent the \textbf{domain} of the operator. This becomes confusing since we use the same symbolics for them. Therefore, for the rest of this derivation I will use $N^{a}_{b}$ to represent the trial function basis set and $\psi^{a}_{b}$ to represent the test function basis set. Using this new symbology we can rewrite the previous two equations as the equivalent forms
|
|
\begin{align}
|
|
-\int_{\Omega}\nabla \psi_{k}^{\theta}\left(\sum_{j=1}^{N^{\phi}_{dof}}\phi_{j}\vec{N_{j}^{\phi}}\right)dV + \int_{\Omega}\psi_{k}^{\theta}\left(\sum_{i=1}^{N_{dof}^{\theta}}\theta_{i}N_{i}^{\theta}\right)^{n}dV &= 0 \\
|
|
\int_{\Omega}\vec{\psi}_{\ell}^{\phi}\cdot\left(\sum_{j=1}^{N^{\phi}_{dof}}\phi_{j}\vec{N_{j}^{\phi}}\right)dV - \int_{\Omega}\vec{\psi}_{\ell}^{\phi}\cdot\nabla\left(\sum_{i=1}^{N^{\theta}_{dof}}\theta_{i}N_{i}^{\theta}\right)dV &= 0
|
|
\end{align}
|
|
Now we exploit the linearity of summation and integration to move the sums out of the integrals
|
|
\begin{align}
|
|
-\sum_{j=1}^{N_{dof}^{\phi}}\phi_{j}\int_{\Omega}\nabla \psi_{k}^{\theta}\cdot\vec{N}_{j}^{\phi}dV + \int_{\Omega}\psi_{k}^{\theta}\left(\theta_{h}\right)^{n} &= 0 \\
|
|
\sum_{j=1}^{N_{dof}^{\phi}}\phi_{j}\int_{\Omega}\vec{\psi}_{\ell}^{\phi}\cdot \vec{N}_{j}^{\phi}dV - \sum_{i=1}^{N_{dof}^{\theta}}\theta_{i}\int_{\Omega}\vec{\psi}_{\ell}^{\phi}\cdot\nabla N_{i}^{\theta} dV &=0
|
|
\end{align}
|
|
We will now define $M_{kj}$, $D_{\ell j}$, and $Q_{\ell i}$ such that
|
|
\begin{align}
|
|
M_{kj} &\equiv \int_{\Omega}\nabla \psi_{k}^{\theta}\cdot \vec{N}_{j}^{\phi}dV \\
|
|
D_{\ell j} &\equiv \int_{\Omega}\vec{\psi}_{\ell}^{\phi}\cdot\vec{N}_{j}^{\phi}dV \\
|
|
Q_{\ell i} &\equiv \int_{\Omega}\vec{\psi}_{\ell}^{\phi}\cdot\nabla N_{i}^{\theta} dV
|
|
\end{align}
|
|
Further we will define $\mathbf{M},\ \mathbf{D},\ \text{and} \ \mathbf{Q}$ such that they are the matrices associated to $M_{kj}$, $D_{\ell j}$, and $Q_{\ell j}$.
|
|
|
|
Note that we do not define a matrix for the non-linear term. This is because we need to treat it as a separate term in computational FEM software, so it is useful for us to split it out now. Instead, let us define $f(\theta)$ to handle the non linear term such that
|
|
\begin{align}
|
|
f(\bar{\theta}) \equiv \int_{\Omega}\psi_{k}^{\theta}\left(\theta_{h}\right)^{n}dV
|
|
\end{align}
|
|
We can write the variational form of our system of equations as
|
|
\begin{align}
|
|
-\sum_{j=1}^{N_{dof}^{\phi}}\phi_{j}M_{kj} + f(\bar{\theta)} &= 0 \\
|
|
\sum_{j=1}^{N^{\phi}_{dof}}\phi_{j}D_{\ell j} - \sum_{i=1}^{N_{dof}^{\theta}}\theta_{i}Q_{\ell i} &= 0
|
|
\end{align}
|
|
Or using the notation we defined
|
|
\begin{align}
|
|
-\mathbf{M}\bar{\phi} + f(\bar{\theta}) &= 0 \\
|
|
\mathbf{D}\bar{\phi} - \mathbf{Q}\bar{\theta} &= 0
|
|
\end{align}
|
|
We can then set this up as a matrix operation
|
|
\begin{align}
|
|
\begin{bmatrix}
|
|
0 & -\mathbf{M} \\
|
|
-\mathbf{Q} & \mathbf{D}
|
|
\end{bmatrix}
|
|
\begin{bmatrix}
|
|
\bar{\theta} \\
|
|
\bar{\phi}
|
|
\end{bmatrix} + \begin{bmatrix}
|
|
f(\bar{\theta}) \\
|
|
0
|
|
\end{bmatrix} = \begin{bmatrix}
|
|
0 \\
|
|
0
|
|
\end{bmatrix}
|
|
\end{align}
|
|
|
|
From this form we can easily see that the residual matrix is
|
|
|
|
\begin{align}
|
|
R &= \begin{bmatrix}
|
|
f(\bar{\theta}) - M\bar{\phi} \\
|
|
D\bar{\phi} - Q\bar{\theta}
|
|
\end{bmatrix}
|
|
\end{align}
|
|
|
|
\subsection{The Jacobian}
|
|
We need to define the Jacobian of this system of equations so that we can use it
|
|
in our Newton-Raphson method. Generally the Jacobian is the matrix of partial derivitives wrt. the state vector. We will let our state vector, $X$, be
|
|
\begin{align}
|
|
X = \begin{bmatrix}
|
|
\bar{\theta} \\
|
|
\bar{\phi}
|
|
\end{bmatrix}
|
|
\end{align}
|
|
|
|
So then the Jacobian is
|
|
\begin{align}
|
|
J &= \begin{bmatrix}
|
|
\frac{\partial}{\partial \theta}\left(f(\theta) - M\phi\right) & \frac{\partial}{\partial \phi}\left(f(\theta) - M\phi\right) \\
|
|
\frac{\partial}{\partial \theta}\left(D\phi - Q\theta\right) & \frac{\partial}{\partial \phi}\left(D\phi - Q\theta\right)
|
|
\end{bmatrix} \\
|
|
J &= \begin{bmatrix}
|
|
\frac{df}{d\theta} - \phi\frac{\partial M}{\partial \theta} & -M-\phi\frac{\partial M}{\partial \phi} \\
|
|
-Q - \theta\frac{\partial Q}{\partial \theta} & D + \phi\frac{\partial D}{\partial \phi} - \theta\frac{\partial Q}{\partial \phi}
|
|
\end{bmatrix}
|
|
\end{align}
|
|
|
|
Finally, we know that the matrices $M$, $D$, and $Q$ are constant with respect to $\theta$ and $\phi$. Therefore, we can drop the partial derivatives with respect to $\theta$ and $\phi$ from the Jacobian. This gives us
|
|
\begin{align}
|
|
\mathbf{J} &= \begin{bmatrix}
|
|
\frac{df}{d\theta} & -M \\
|
|
-Q & D
|
|
\end{bmatrix}
|
|
\end{align}
|
|
|
|
\noindent In a fully assembled, distritized form this will look like
|
|
|
|
\begin{align}
|
|
\mathbf{J} = \begin{bmatrix} \frac{df}{d\theta}_{00} & \dots & \frac{df}{d\theta}_{0n_{\theta}} & -M_{00} & \dots & -M_{0n_{\phi}} \\
|
|
\vdots & \ddots & & \vdots & \ddots & \\
|
|
\frac{df}{d\theta}_{n_{\theta}0} & & \frac{df}{d\theta}_{n_{\theta}n_{\theta}} & -M_{n_{\theta}0} & & -M_{n_{\theta}n_{\phi}} \\
|
|
-Q_{00} & \dots & -Q_{0n_{\theta}} & D_{00} & \dots & D_{0n_{\phi}} \\
|
|
\vdots & \ddots & & \vdots & \ddots & \\
|
|
-Q_{n_{\phi}0} & & -Q_{n_{\phi}n_{\theta}} & D_{n_{\phi}0} & & D_{n_{\phi}n_{\phi}}
|
|
\end{bmatrix}
|
|
\end{align}
|
|
|
|
\noindent Where $N_{dof}^{\theta} = n_{\theta}$ is the number of degrees of freedom on $\theta$, and $N_{dof}^{\phi} = n_{\phi}$ is the number of degrees of freedom on $\phi$. Note how the Jacobian is a matrix of size $\left(N_{dof}^{\theta} + N_{dof}^{\phi} \times N_{dof}^{\theta} + N_{dof}^{\phi}\right)$
|
|
|
|
\subsection{Preconditioner}
|
|
Due to the eventual size of these matrices we would like to be able to solve
|
|
each step in this using a memory efficiet approach. Krylov solvers, such as
|
|
GMRES, allow for matrix free iterative solutions (as long the concept of
|
|
multiplication is defined). However, for these systems to be well formed for
|
|
such solvers it is useful for us to use a preconditioner. However, this is a
|
|
somewhat strongly coupled system where we cannot simply use the inverse
|
|
diagonals of the matrix as a preconditioner. Instead, to encode the coupling we
|
|
will use Schur's Compliment. Each Newton iteration we solve the equation
|
|
\begin{align}
|
|
\mathbf{J}\Delta \vec{x} = \vec{b}
|
|
\end{align}
|
|
If we expand this out
|
|
\begin{align}
|
|
\begin{bmatrix} \mathbf{\dot{f}} & -\mathbf{M} \\
|
|
-\mathbf{Q} & \mathbf{D}
|
|
\end{bmatrix}\begin{bmatrix} \theta \\
|
|
\phi
|
|
\end{bmatrix} = \begin{bmatrix} b_{0} \\
|
|
b_{1}\end{bmatrix}
|
|
\end{align}
|
|
We can pull out the first equation from this system
|
|
\begin{align}
|
|
\mathbf{\dot{f}}\theta - \mathbf{M}\phi &= b_{0} \\
|
|
\theta &= \mathbf{\dot{f}}^{-1}b_{0} + \mathbf{\dot{f}}^{-1}\mathbf{M}\phi
|
|
\end{align}
|
|
Then if we pull out the second equation from the system
|
|
\begin{align}
|
|
-\mathbf{Q}\theta + \mathbf{D}\phi &= b_{1} \\
|
|
-\mathbf{Q}\left(\mathbf{\dot{f}}^{-1}b_{0} + \mathbf{\dot{f}}^{-1}\mathbf{M}\phi\right) + \mathbf{D}\phi &= b_{1}
|
|
\end{align}
|
|
rearanging terms a bit
|
|
\begin{align}
|
|
-\mathbf{Q}\mathbf{\dot{f}}^{-1}b_{0}-\mathbf{Q}\mathbf{\dot{f}}^{-1}\mathbf{M}\phi + \mathbf{D}\phi &= b_{1} \\
|
|
\left(\mathbf{D} - \mathbf{Q}\mathbf{\dot{f}}^{-1}\mathbf{M}\right)\phi &= b_{1} + \mathbf{Q}\mathbf{\dot{f}}^{-1}b_{0}
|
|
\end{align}
|
|
The term $\mathbf{D}-\mathbf{Q}\mathbf{\dot{f}}^{-1}\mathbf{M}$ is Schur's Compliment for this system, and we will represent this by the symbol $\mathbf{\tilde{S}}$. We can use Schur's Compilment to precondition our equation if we let the preconditioner be of the form
|
|
\begin{align}
|
|
\mathbf{P} = \begin{bmatrix} \mathbf{\dot{f}}^{-1} & 0 \\
|
|
0 & \mathbf{\tilde{S}}^{-1}
|
|
\end{bmatrix}
|
|
\end{align}
|
|
So then the preconditioned equation which can be more easily solved by some Krylov solver (such as GMRES) is
|
|
\begin{align}
|
|
\mathbf{P}\mathbf{J}\Delta \vec{x} = \mathbf{P}\vec{b}
|
|
\end{align}
|
|
|
|
It is easy to see here that for this system to be solvable / well defined both
|
|
$\mathbf{\tilde{S}}$ and $\mathbf{\dot{f}}$ need to be invertable
|
|
matrices. They are both easily shown to be square (with $\mathbf{\tilde{S}}$
|
|
having a size $\left(N_{dof}^{\phi}\times N_{dof}^{\phi}\right)$ and
|
|
$\mathbf{\dot{f}}$ having a size $\left(N_{dof}^{\theta}\times
|
|
N_{dof}^{\theta}\right)$).
|
|
|
|
|
|
\subsection{A Few Quick Notes}
|
|
A few notes on the dimensions of $\mathbf{M}$, $\mathbf{Q}$, $\mathbf{D}$, and $f(\bar{\theta})$.
|
|
\begin{itemize}
|
|
\item $\mathbf{M}$ is a matrix of size $\left(N_{dof}^{\theta}\ \times\ N_{dof}^{\phi}\right)$.
|
|
\item $\mathbf{Q}$ is a matrix of size $\left(N_{dof}^{\phi}\ \times\ N_{dof}^{\theta}\right)$.
|
|
\item $\mathbf{D}$ is a matrix of size $\left(N_{dof}^{\phi}\ \times\ N_{dof}^{\phi}\right)$.
|
|
\item $f(\bar{\theta})$ is a vector of size $N_{dof}^{\theta}$.
|
|
\end{itemize}
|
|
|
|
\section{Representation in FEM}
|
|
We will make use of the MFEM library\footnote{https://mfem.org/} to encode this system of equations. This document is not intended to be a comprehensive guide to using MFEM; rather, here we will provide an explanation for how to translate $\mathbf{M}$, $\mathbf{D}$, and $\mathbf{Q}$ into pre-existing MFEM integrators. The non linear term must be written as a custom integrator and an explanation of this process is outside the scope of this document.
|
|
|
|
\subsection{MFEM Integrators}
|
|
MFEM provides an extensive set of integrators. Of interest here are the BilinearFormIntegrators and MixedBilinearFormIntegrators. We will explain how to use these by following the process of deciding how $\mathbf{M}$ should be represented. Recall that
|
|
\begin{align}
|
|
\mathbf{M} &= \left[M_{kj}\right] \\
|
|
M_{kj} &= \int_{\Omega}\nabla\psi_{k}^{\theta}\cdot\vec{N}_{j}^{\phi}dV
|
|
\end{align}
|
|
Also recall that $\psi$ denotes the test space while $N$ denotes the trial space. MFEM provides a robust set of integrators. Because $\mathbf{M}$ is composed of terms from the $\theta$ and $\psi$ spaces it is what is called a Mixed form. Therefore, we will look at the mixed form integrators provided by MFEM.
|
|
\begin{table}[htbp]
|
|
\scriptsize
|
|
\centering
|
|
\caption{Selection of MFEM Mixed Bilinear Form Integrators}
|
|
\label{tab:mfem_mixed_integrators}
|
|
% Adjust column specifications (l, c, r, p{}) as needed for alignment and width
|
|
\begin{tabular}{@{}lllllll@{}}
|
|
\toprule
|
|
Class Name & Domain & Range & Coef. & Operator & Continuous Op. & Dimension \\
|
|
\midrule
|
|
MixedDotProductIntegrator & ND, RT & H1, L2 & V & $(\vec{\lambda} \cdot \vec{u}, v)$ & $\vec{\lambda} \cdot \vec{u}$ & 2D, 3D \\
|
|
MixedScalarCrossProductIntegrator & ND, RT & H1, L2 & V & $(\vec{\lambda} \times \vec{u}, v)$ & $\vec{\lambda} \times \vec{u}$ & 2D \\
|
|
MixedVectorWeakDivergenceIntegrator & ND, RT & H1 & S, D, M & $(-\lambda \vec{u}, \nabla v)$ & $\nabla \cdot (\lambda \vec{u})$ & 2D, 3D \\
|
|
MixedWeakDivCrossIntegrator & ND, RT & H1 & V & $(-\vec{\lambda} \times \vec{u}, \nabla v)$ & $\nabla \cdot (\vec{\lambda} \times \vec{u})$ & 3D \\
|
|
MixedVectorMassIntegrator & ND, RT & ND, RT & S, D, M & $(\lambda \vec{u}, \vec{v})$ & $\lambda \vec{u}$ & 2D, 3D \\
|
|
MixedCrossProductIntegrator & ND, RT & ND, RT & V & $(\vec{\lambda} \times \vec{u}, \vec{v})$ & $\vec{\lambda} \times \vec{u}$ & 3D \\
|
|
MixedVectorWeakCurlIntegrator & ND, RT & ND & S, D, M & $(\lambda \vec{u}, \nabla \times \vec{v})$ & $\nabla \times (\lambda \vec{u})$ & 3D \\
|
|
MixedWeakCurlCrossIntegrator & ND, RT & ND & V & $(\vec{\lambda} \times \vec{u}, \nabla \times \vec{v})$ & $\nabla \times (\vec{\lambda} \times \vec{u})$ & 3D \\
|
|
MixedScalarWeakCurlCrossIntegrator & ND, RT & ND & V & $(\vec{\lambda} \times \vec{u}, \nabla \times \vec{v})$ & $\nabla \times (\vec{\lambda} \times \vec{u})$ & 2D \\ % Note: Same Op as previous? Verify source.
|
|
MixedWeakGradDotIntegrator & ND, RT & RT & V & $(-\vec{\lambda} \cdot \vec{u}, \nabla \cdot \vec{v})$ & $\nabla(\vec{\lambda} \cdot \vec{u})$ & 2D, 3D \\
|
|
MixedScalarCurlIntegrator & ND & H1, L2 & S & $(\lambda \nabla \times \vec{u}, v)$ & $\lambda \nabla \times \vec{u}$ & 2D \\
|
|
MixedCrossCurlGradIntegrator & ND & H1 & V & $(\vec{\lambda} \times (\nabla \times \vec{u}), \nabla v)$ & $-\nabla \cdot (\vec{\lambda} \times (\nabla \times \vec{u}))$ & 3D \\
|
|
MixedVectorCurlIntegrator & ND & ND, RT & S, D, M & $(\lambda \nabla \times \vec{u}, \vec{v})$ & $\lambda \nabla \times \vec{u}$ & 3D \\
|
|
MixedCrossCurlIntegrator & ND & ND, RT & V & $(\vec{\lambda} \times (\nabla \times \vec{u}), \vec{v})$ & $\vec{\lambda} \times (\nabla \times \vec{u})$ & 3D \\
|
|
MixedScalarCrossCurlIntegrator & ND & ND, RT & V & $(\vec{\lambda} \times \hat{z} (\nabla \times \vec{u}), \vec{v})$ & $\vec{\lambda} \times \hat{z} (\nabla \times \vec{u})$ & 2D \\
|
|
MixedCurlCurlIntegrator & ND & ND & S, D, M & $(\lambda \nabla \times \vec{u}, \nabla \times \vec{v})$ & $\nabla \times (\lambda \nabla \times \vec{u})$ & 3D \\
|
|
MixedCrossCurlCurlIntegrator & ND & ND & V & $(\vec{\lambda} \times (\nabla \times \vec{u}), \nabla \times \vec{v})$ & $\nabla \times (\vec{\lambda} \times (\nabla \times \vec{u}))$ & 3D \\
|
|
MixedScalarDivergenceIntegrator & RT & H1, L2 & S & $(\lambda \nabla \cdot \vec{u}, v)$ & $\lambda \nabla \cdot \vec{u}$ & 2D, 3D \\
|
|
MixedDivGradIntegrator & RT & H1 & V & $(\vec{\lambda} (\nabla \cdot \vec{u}), \nabla v)$ & $-\nabla \cdot (\vec{\lambda} (\nabla \cdot \vec{u}))$ & 2D, 3D \\
|
|
MixedVectorDivergenceIntegrator & RT & ND, RT & V & $(\vec{\lambda} (\nabla \cdot \vec{u}), \vec{v})$ & $\vec{\lambda} (\nabla \cdot \vec{u})$ & 2D, 3D \\
|
|
\bottomrule
|
|
\end{tabular}
|
|
\end{table}
|
|
|
|
There is a lot of information in this table so we will break it down. First we need to identify what spaces the domain and range of our operator exist in. The range of the operator is that which contains the test function while the domain is the space containing the trial function. For $\mathbf{M}$ the test function is in the $\theta$ space, or H1, while the trial function is in the $\phi$ space, or RT.
|
|
|
|
Next we look at the Operator column. These define the operation within the integral where $(a,b)$ is the inner product of $a$ and $b$. More specifically, the MFEM documentation provides that $u$ is the trial function and $v$ is the test function. So then we are looking for an integrator which has the operator of the inner product of the trial function and the gradient of the test function while also satisfying the domain and range constraints. Upon investigation of this table we can see that the \texttt{MixedVectorWeakDivegenceIntegrator} has a range defined on H1 and domain defined on RT, just like we need. Further, its operator is given as $(-\lambda \vec{u},\nabla v)$. This is the same form as we have. $\left(-\lambda\vec{u} \rightarrow 1\times\vec{N}_{j}^{\phi}, \nabla v\rightarrow \nabla \psi_{k}^{\theta}\right)$. Therefore, \texttt{MixedVectorWeakDivergenceIntegrator} will compute the matrix $\textbf{M}$ over our domain without any modifications (as long as we are careful about the sign on the coefficient $\lambda$).
|
|
|
|
The other integrators map to \texttt{VectorFEMassIntegrator} and \texttt{MixedVectorGradientIntegrator} for $\mathbf{D}$ and $\mathbf{Q}$ respectively. From here one would assemble these, along with the non-linear term into a block form; however, that is beyond the scope of this document.
|
|
\end{document}
|