Extending FEniCS to Work in Higher Dimensions Using Tensor Product Finite Elements
We present a method to extend the finite element library FEniCS to solve problems with domains in dimensions above three by constructing tensor product finite elements. This methodology only requires that the high dimensional domain is structured as a Cartesian product of two lower dimensional subdomains. In this study we consider Dirichlet problems for scalar linear partial differential equations, though the methodology can be extended to non-linear problems. The utilization of tensor product finite elements allows us to construct a global system of linear algebraic equations that only relies on the finite element infrastructure of the lower dimensional subdomains contained in FEniCS. We demonstrate the effectiveness of our methodology in four distinctive test cases. The first test case is a Poisson equation posed in a four dimensional domain which is a Cartesian product of two unit squares solved using the classical Galerkin finite element method. The second test case is the wave equation in space-time, where the computational domain is a Cartesian product of a two dimensional space grid and a one dimensional time interval. In this second case we also employ the Galerkin method. The third test case is an advection dominated advection-diffusion equation where the global domain is a Cartesian product of two one dimensional intervals in which the streamline upwind Petrov-Galerkin method is applied to ensure discrete stability. The final test case uses the Galerkin approach to solve a Poisson problem on a Cartesian product of two intervals with a spatially varying, non-separable diffusivity term. In all cases, a p=1 basis is used and optimal L^2 convergence rates of order h^{p+1} of the errors are achieved with respect to h refinement
\section{Introduction}

In the last decade, open source libraries with high level APIs that automate the process of solving partial differential equations (PDEs) using the finite element method (FEM) have become valuable tools in computational research. Some prominent libraries of this kind include Firedrake~\cite{rathgeber2016firedrake}, deal.II~\cite{bangerth2007deal}, MFEM~\cite{anderson2021mfem}, and  FEniCS~\cite{alnaes2015fenics}. FEniCS is one of the most widely used of these libraries. FEniCS, along with most other FEM libraries, only include capability of handling up to three dimensional problems (GetFEM~\cite{renard2020getfem} is a notable exception). FEniCS specifically supports unstructured meshes with corresponding basis functions in one, two, and three dimensions. Extension of FEniCS to solve problems in higher dimensions is of interest for problems such as spectral wind wave models, spatial population genetics, quantum mechanics, and even 3d space-time problems.

FEniCS does not contain support for fully unstructured meshes in dimensions higher than three. However, at the expense of losing a fully unstructured mesh, a higher dimensional space can be discretized with the available tools in FEniCS if it is the Cartesian product of two lower dimensional (3 or lower) spaces that can themselves be unstructured. It is well known that in a Cartesian product space, a finite element basis can be constructed called the product basis which spans the function space of the full domain~\cite{brenner2008mathematical,quarteroni2010numerical,ern2013theory}. This research seeks to use the infrastructure of an FEM library, such as FEniCS, in the lower dimensional spaces to construct a basis for this high dimensional space. 

 Using FEM to discretize a structured domain via a product basis has been done before without the use of FEM libraries such as FEniCS. In 1978, Banks~\cite{Bank1978} used tensor product finite elements to solve a 2D Poisson problem on a structured grid in order to find a faster solver. In the study, the domain was a Cartesian product between two unit intervals resulting in a uniform square mesh. Banks used the tensor product of two one dimensional quadratic and cubic basis functions to construct the polynomial basis of the full 2D domain. In 1980, Baker performed numerical tests as well as a stability analysis on a tensor product finite method with application to convection-dominated fluid flow problems~\cite{BAKER1981215}. The stability analysis showed the basic algorithm is spatially fourth- order accurate in its most elementary embodiment and the numerical experiments on the convection-dominated model test problems confirmed the basic viability of the developed algorithm, and its tensor product formulation. Recently, Du \eirik{\emph{et al.}~\cite{DU2013181}} used tensor product finite elements to construct a fast solver for an electromagnetics scattering problem of a cavity. In  Firedrake~\cite{rathgeber2016firedrake}, there exists a capability to create tensor product finite elements only up to three dimensions~\cite{mcrae2016automated}. A package built on deal.II called Hyperdeal~\cite{munch2020hyperdeal} has the capability of creating tensor product finite elements in up to six dimensions. deal.II is different to a library such as FEniCS since deal.II uses quadrilateral elements in two dimensions whereas FEniCS uses triangles. \eirik{Here, we present a general software framework built on the components of the FEniCS library as well as other open source Python libraries to construct high dimensional meshes and corresponding FE discretizations.  }


Following this introduction, the general set of problems for which this paper focuses on will be defined and the notation for the product basis will be introduced \eirik{in Section~\ref{sec:model_prob}}. Then, a set of three different model problems will be described in detail and the derivation of the system of equations using tensor product finite elements will be shown for each case. The first model problem will be discussed in Section~\ref{sec:model_prob} which is an N dimensional Poisson problem where the domain is a Cartesian product between two subdomains of dimension three or less. In Section~\ref{sec:prob2} the second model problem is discussed which is a wave equation where the domain is decomposed as a Cartesian product between space and time. The third model problem is discussed in Section~\ref{sec:prob3} which is an advection dominated advection-diffusion equation where the advection aligns with one of the subdomains, a streamlined upwind Petrov-Galerkin method (SUPG) is formulated for this case. After each system of equations is derived, numerical tests were run using FEniCS for each model problem where specific boundary conditions were given. For each case, error and convergence rates are tabulated in Section~\ref{sec:numerical_results}. Lastly, conclusions and recommendations for future work are given in Section~\ref{sec:Conclusions}.


\section{Methods} \label{sec:methods}

To present the proposed methodology and algorithms, we 
consider the following class of problems, i.e., PDEs:
\begin{equation*}
\begin{split}
    \mathscr{L} u = f \quad \textrm{in} \quad \Omega = \Omega_1 \times \Omega_2.
\end{split}
\end{equation*}
Where $\mathscr{L}$ is a linear differential operator, $f$ is a forcing function, and the domain $\Omega$ is defined as a Cartesian product between two lower dimensional Lipschitz domains.  For example, if the global domain $\Omega  \subset \mathbb{R}^2$, then $\Omega$ can be defined by the Cartesian product of 2 intervals $\Omega_1\subset \mathbb{R}$ and  $\Omega_2 \subset \mathbb{R}$. Hence, we  consider general domains of the form $\Omega = \Omega_1 \times \Omega_2 $, see Figure~\ref{fig:cartesian_illustration} for an illustration.
\begin{figure}[h!]
    \centering
    \includegraphics[width=\textwidth]{CartesianProductIllust.png}
    \caption{Example of a global domain $\Omega$ that is a Cartesian product of two lower dimensional subdomains $\Omega_1, \Omega_2$.}
    \label{fig:cartesian_illustration}
\end{figure}
FEniCS can discretize up to three dimensional objects, hence,  in practice this framework can be used to define domains that are Cartesian products of up to six dimensions. Furthermore, this process can be done iteratively to yield even higher dimensional domains as $\Omega_1$ itself could be a product of two other spaces and so on. We note that in this presentation, we consider only symmetric functional settings which admit well posed  Galerkin FE discrtetization. 

It can be shown that if a function $u$ is defined on a Cartesian product domain $\Omega = \Omega_1 \times \Omega_2$ then we can construct a basis of a polynomial space to be used for the FEM \eirik{ by exploiting this underlying geometric structure, see, e.g.,}~\cite{quarteroni2010numerical,brenner2008mathematical,ern2013theory}. 
 \eirik{With} a basis for a polynomial space for the first lower dimensional subdomain $\Omega_1$:
\begin{equation}
    \{\phi_i\}_{i=1}^N,
\end{equation}
and a basis for the second subdomain $\Omega_2$:
\begin{equation}
    \{\psi_j\}_{j=1}^M.
\end{equation}
The basis for the entire domain can subsequently be constructed and any arbitrary function $u$ whose domain is in $\Omega = \Omega_1 \times \Omega_2$ can be approximated via the product basis:
\begin{equation}
    u \approx \sum_{i=1}^N\sum_{j=1}^M u_{i,j} \phi_i\, \psi_j. 
\end{equation}
\eirik{Note that in the following,} we use boldface letters and symbols to denote vector quantities, e.g., $\dX = \dx \dy$. \eirik{To avoid raising complexity in implementation, we also assume in the numerical test cases (Section~\ref{sec:numerical_results}) that all physical parameters in the following model problems are constant and that the Dirichlet boundary conditions are homogeneous.}
\matt{However, it is possible to remove these assumptions by projecting the parameters and boundary conditions to the appropriate function space.}

\subsection{Model Problem 1: N dimensional Poisson Equation} \label{sec:model_prob}
As a first model problem to illustrate the methodology, we consider the Poisson equation:
\begin{equation}\label{eqn:model_prob1}
\begin{split}
        - \Delta \, u &= f \quad \textrm{in} \quad \Omega, \\
        u &= u_D \quad \textrm{on} \quad \partial \Omega, 
\end{split}
\end{equation}
where the source $f$ is in $L^2(\Omega)$ and the source data $u_D$ is assumed to be sufficiently regular.
To define the weak formulation for~\eqref{eqn:model_prob1}, multiply both sides with a test function $v$ in $L^2(\Omega)$ and integrate over $\Omega$:
\begin{equation}\label{eqn:weak_form1}
    \int_{\Omega}  - \Delta \, u\, v \dX = \int_{\Omega} f \, v \dX \quad \forall v \in L^2(\Omega), 
\end{equation}
Integrating by parts on the left side (assuming Dirichlet boundary conditions on entire boundary) gives the following weak formulation\matt{: find $u \in \mathcal{U}(\Omega)$ such that}
\begin{equation} \label{eq:weak_form}
    \int_{\Omega}  \Nabla u \cdot \Nabla v \dX = \int_{\Omega} f\,v \dX \quad \forall \quad v \in \mathcal{V}(\Omega), 
\end{equation}
where the function space $\mathcal{V}(\Omega)$ is the Hilbert space with zero trace on the boundary  $H^1_0(\Omega)$ and $\mathcal{U}(\Omega)$ is $\mathcal{V}(\Omega)$ with a finite energy lift on the boundary so that the Dirichlet condition $u=u_D$ is satisfied.
The weak formulation in~\eqref{eq:weak_form}, and its corresponding discretization is known to be well posed, see, e.g.,~\cite{becker1981finite}. 

With a well posed weak formulation at hand, we can discretize this weak form using a finite element basis. In this case, the functional setting dictates the use of a $C^0$ continuous \eirik{polynomial} basis for the classical FEM.
Hence, we approximate the trial functions with the product basis:
\begin{equation}
    u \approx \sum_{i=1}^N\sum_{j=1}^M u_{i,j} \, \phi_i \, \psi_j,
\end{equation}
and the test functions by its product basis:
\begin{equation}
    v \approx \sum_{k=1}^N\sum_{l=1}^M v_{k,l} \, \gamma_k \, \beta_l.
\end{equation}
The weak form from~\eqref{eqn:weak_form1} can then be discretized by substitution of the product bases:
\begin{equation}
\displaystyle    \int_{\Omega}  (\Nabla \sum_{i=1}^N\sum_{j=1}^M u_{i,j} \, \phi_i \, \psi_j) \cdot ( \Nabla \sum_{k=1}^N\sum_{l=1}^M  v_{k,l}\,\gamma_k\, \beta_l) \dX = \int_{\Omega} f\, (\sum_{k=1}^N\sum_{l=1}^M  v_{k,l}\,\gamma_k\, \beta_l ) \dX.
\end{equation}
Due to arbitrariness of the test function and the bilinearity of the weak form, this implies the following form:
\begin{equation}
    \sum_{i=1}^N\sum_{j=1}^M \int_{\Omega}  (\Nabla   \, \phi_i \, \psi_j) \cdot ( \Nabla \,\gamma_k\, \beta_l) \dX \, u_{i,j} = \int_{\Omega} f\, (\,\gamma_k\, \beta_l ) \dX.
\end{equation}
This left hand side \eirik{results} in a \eirik{product between a} 4 dimensional \eirik{and} a 2 dimensional tensor $u_{i,j}$. \eirik{Consequently,} the 4 dimensional tensor $A_{ijkl}$ will be of the form:
\begin{equation}
        \int_{\Omega}   \Nabla  \phi_i \,\psi_j \cdot \Nabla  \gamma_k \,\beta_l \, \dX.
\end{equation}
By the product rule we have:
\begin{equation} \label{eq:chain_ruled}
        \int_{\Omega}   \Nabla  \phi_i \, \psi_j \cdot \Nabla  \gamma_k \,\beta_l \dX  = \int_{\Omega}\, (\phi_i \, \Nabla \psi_j + \psi_j\, \Nabla \phi_i) \cdot (\gamma_k\, \Nabla \beta_l + \beta_l\, \Nabla \gamma_k) \dX.
\end{equation}
By construction, \eirik{the} $\phi$\eirik{'s}  and$\gamma$\eirik{'s} only vary in the first subdomain $\Omega_1$ whereas \eirik{the} $\psi$\eirik{'s} and $\beta$\eirik{'s} only vary in the second subdomain $\Omega_2$. \eirik{Thus, the integral form~\eqref{eq:chain_ruled} can be simplified. To this end, we use the following notation convention:} if the gradient operator on the entire domain \eirik{$\Omega$} is $\Nabla$, define the gradient on the subdomains \eirik{$\Omega_1$ and $\Omega_2$} \eirik{as} $\Nabla_1$, $\Nabla_2$\eirik{, respectively}. \eirik{Hence, by} construction $\Nabla = (\Nabla_1, \Nabla_2)$.
Rewriting \eirik{ the gradients in~\eqref{eq:chain_ruled} then} gives:
\begin{equation} \label{eq:expansion}
\begin{split}
   & \int_{\Omega} (\phi_i \, (\Nabla_1, \Nabla_2)\,\psi_j + \psi_j \,(\Nabla_1, \Nabla_2)\, \phi_i) \cdot (\gamma_k\, (\Nabla_1, \Nabla_2)\, \beta_l + \beta_l (\Nabla_1, \Nabla_2)\gamma_k)\dX,  \\ \qquad &=
    \int_{\Omega} (  (\mathbf{0}, \phi_i\,\Nabla_2 \psi_j) +  (\psi_j\,\Nabla_1\phi_i, \mathbf{0}) ) \cdot ( (\mathbf{0}, \gamma_k\, \Nabla_2\beta_l)  +  (\beta_l\,\Nabla_1\gamma_k, \mathbf{0}))\dX,  \\ &=
    \int_{\Omega}  (\psi_j\,\Nabla_1\phi_i, \phi_i\,\Nabla_2 \psi_j)  \cdot  (\beta_l\,\Nabla_1\gamma_k, \gamma_k\, \Nabla_2\beta_l)  \dX, \\ \qquad &= 
    \int_{\Omega}  \Nabla_1\phi_i \cdot \Nabla_1\gamma_k\,\psi_j\,\beta_l +  \phi_i \,\gamma_k \,\Nabla_2 \psi_j \cdot \Nabla_2\,\beta_l  \dX. 
\end{split}
\end{equation}
Since \eirik{the} domain $\Omega$ is a Cartesian product of \eirik{the} subdomains $\Omega_1, \Omega_2$ we can rewrite the \eirik{last} integral \eirik{in~\eqref{eq:expansion},
as}:
\begin{equation}
           \int_{\Omega_1} \int_{\Omega_2}  \Nabla_1\phi_i \cdot \Nabla_1\gamma_k\,\psi_j\,\beta_l +  \phi_i\, \gamma_k\, \Nabla_2 \psi_j \cdot \Nabla_2\beta_l  \dy  \dx.
\end{equation}
\eirik{An application of} Fubini's theorem \eirik{gives}:
\begin{equation}
               \int_{\Omega_1} \Nabla_1\phi_i \cdot \Nabla_1\gamma_k \dx \int_{\Omega_2}  \psi_j\,\beta_l \dy + \int_{\Omega_1} \phi_i\, \gamma_k \dx \int_{\Omega_2}  \Nabla_2 \psi_j \cdot \Nabla_2\beta_l   \dy,
\end{equation}
These can be written as Kronecker products of  matrices,  \eirik{e.g., stiffness matrices computed with} FEniCS. \eirik{To continue the discussion, we introduce}  the following notation:
\begin{equation}
    \begin{split}
        K_{11} = \int_{\Omega_1} \Nabla_1\phi_i \cdot \Nabla_1\gamma_k \dx, \quad K_{22} = \int_{\Omega_2}  \psi_j\;\beta_l \dy,\\
        K_{12} = \int_{\Omega_1} \phi_i\; \gamma_k \dx, \quad
        K_{21} = \int_{\Omega_2}  \Nabla_2 \psi_j \cdot \Nabla_2\beta_l   \dy,
    \end{split}
\end{equation}
\eirik{where} the first index indicates subdomain (\eirik{i.e.,} 1 or 2 \eirik{in this case}) and the second indicates the term from the weak form (1 is the first order operator and 2  the \eirik{$0^{th}$} order). \eirik{Hence, each local matrix $K_{ij}$ can be constructed independently, and the global stiffness matrix $A_{ijkl}$ can be defined by the sum of 2 Kronecker products.}
The Kronecker product allows us to represent \eirik{the matrix of a} 4D tensor  \eirik{by smaller 2D matrices}:
\begin{equation}
    A = K_{11} \otimes K_{22} + K_{12} \otimes K_{21}.
\end{equation}
\eirik{A similar reasoning for the right hand side leads to}:
\begin{equation}
    \int_{\Omega}  f_{dof}\;\phi_k\; \beta_l \dX  = f_{dof} F_1 \otimes F_2
\end{equation}
Where $f_{dof}$ is the pointwise value of $f$ at each $k,l$ coordinate in the global space. The forcing vectors are defined as:
\begin{equation}
    F_{1} = \int_{\Omega_1}   \phi_k \dx \quad F_{2} = \int_{\Omega_2}  \beta_l \dy.
\end{equation}
\eirik{Finally,} the entire system \eirik{of equations} becomes:
\begin{equation} \label{eq:final_poisson_form}
    (K_{11} \otimes K_{22} + K_{12} \otimes K_{21})(U) = f_{dof} F_1 \otimes F_2,
\end{equation}
where $U$ are the values at all d.o.f in the Cartesian product space for $u$. 
\eirik{A subsequent application of boundary conditions to~\eqref{eq:final_poisson_form}  leads to the final system of linear algebraic equations.} 


\subsection{Model Problem 2: Arbitrary dimensional Space-Time Wave Equation}\label{sec:prob2}
\eirik{The second model problem we consider is the linear wave equation. We consider this transient problem to highlight the application of this methodology to transient problems where space-time finite elements are employed. Thus, we have the following model problem:}
\begin{equation}\label{eqn:model_prob2}
\begin{split}
        \frac{\partial^2 u}{\partial t}-c^2\Delta u & = 0 \quad \textrm{in} \quad \Omega_T, \\
        \eirik{u} &\eirik{ =u_D \quad \textrm{on} \quad \partial \Omega,} \\
        \eirik{u} &\eirik{ =u_{initial} \quad \textrm{on} \quad \Omega,}
\end{split}
\end{equation}
\eirik{where $c$ denotes the wave speed}.
In this problem the space-time domain $\Omega_T$ is a Cartesian product domain of one, two, or three dimensional \eirik{spatial domain and a one dimensional temporal domain} ($\Omega_T = \Omega \times (0,T)$). The weak form of~\eqref{eqn:model_prob2} is obtained by multiplying by a test function and integrating over the entire space-time domain:
\begin{equation}
    \int_{\Omega_T} \frac{\partial^2 u}{\partial t} v  - c^2\Delta u v \dX  \, {\rm d t}     = 0 \quad \forall \quad v \in L^2(\Omega_T),
\end{equation}
\eirik{and subsequent integration } by parts in \eirik{space and time}  gives the following weak form. . Find $u \in \mathcal{U}(\Omega_T)$:
\begin{equation}\label{eqn:weak_form2}
   \int_{\Omega_T} -\frac{\partial u}{\partial t} \frac{\partial v}{\partial t} + c^2\Nabla u \cdot \Nabla v \dX \, {\rm d t}     + \int_{\partial \Omega_T|_{t=T}}\frac{\partial u}{\partial t}v  \dX = 0 \quad \forall \quad v \in \mathcal{V}(\Omega_T),
\end{equation}
where \eirik{we have applied Dirichlet conditions to the space-time boundary $\partial \Omega_T$, except at the final time boundary.} $\mathcal{V}(\Omega_T)$ is the space of all $H^1_0(\Omega_T)$ functions except on $\partial \Omega_T|_{t=T}$, where the trace is an unkown, and $\mathcal{U}(\Omega_T)$ is the space of all $H^1_0(\Omega_T)$ plus the trace $u=u_D$ on $\partial \Omega_T$ except the aforementioned part of the boundary $\Omega_T|_{t=T}$. 

\eirik{The corresponding discretization of}~\eqref{eqn:weak_form2} using a product basis gives a very similar system of equations to the Poisson problem \eirik{considered in Section~\ref{sec:model_prob}}, with a slight variation. The term $K_{21}$ must include \eirik{additional integrals} and becomes:
\begin{equation}
    K_{21} = \int_0^T  \frac{\partial \psi_j}{\partial t} \frac{\partial \beta_l}{\partial t}   dt - \int_{\partial \Omega_2} \Nabla \psi_j \cdot \mathbf{n} \beta_l  \, {\rm d s}    .
\end{equation}
\eirik{Consequently}, the global system of equations becomes:
\begin{equation} \label{eq:global_wave_sys}
    (c^2 K_{11} \otimes K_{22} - K_{12} \otimes K_{21})U = \mathbf{0}.
\end{equation}
\eirik{Finally, application of boundary and initial conditions to~\eqref{eq:global_wave_sys} results in the final system of linear algebraic equations.} 

\subsection{Model Problem 3: SUPG Stabilized Advection Dominated Advection Diffusion Equation}\label{sec:prob3}

\eirik{To highlight the versatility of our approach to consider non-standard FE techniques, we consider a PDE which is known to lead to stability issues in the Galerkin FE setting. Hence, we consider an advection-diffusion  PDE  in which advection is the dominant:}  
\begin{equation}\label{eqn:model_problem3}
\begin{split}
        -\kappa \, \Delta u + \mathbf{b} \cdot \Nabla u = f, \quad \textrm{on} \quad \Omega \\
         u = u_D \quad \textrm{on} \quad \partial \Omega,,
\end{split}
\end{equation}
where $\Omega$ is defined as a Cartesian Product of two lower dimensional Lipschitz domains: $\Omega = \Omega_1 \times \Omega_2$. \eirik{By following the standard procedure of deriving integral formulations,} we get the corresponding weak form:
\begin{equation}\label{eqn:weak_form3}
        \int_{\Omega} \kappa \Nabla u \cdot \Nabla v + \mathbf{b} \cdot \Nabla u v dx = \int_{\Omega}fv dx \quad \forall v \in \mathcal{V}(\Omega)
\end{equation}
In problems where the advection term dominates the diffusion, that is when Peclet number $Pe = \frac{L \| \mathbf{b} \| }{\kappa} >> 1$, where $L$ is the characteristic length, the standard Galerkin \eirik{method applied to~\eqref{eqn:weak_form3} may result in a}   discretization that is \eirik{unstable. } \eirik{This issue of stability can be overcome by careful design of the FE mesh or through stabilization techniques that ensure satisfaction of the discrete \emph{inf-sup} condition. 
Here, we consider the SUPG method introduced by Brooks and Hughes~\cite{brooks1982streamline} since it is widely used and has well developed criteria for discrete stability.}
\eirik{The SUPG method leads to stable FE discretizations by adjusting the discretized weak form~\eqref{eqn:weak_form3} with a penalized residual, i.e., find $u_h \in \mathcal{U}_h(\Omega)$: }
\begin{equation}\label{eqn:weak_form3b}
   \kappa(\Nabla u_h , \Nabla v_h)_{\Omega} + (\mathbf{b} \cdot \Nabla u_h, v_h)_{\Omega} + (\underbrace{-\kappa\,\Delta u_h + \mathbf{b} \cdot \Nabla u_h -f}_{\eirik{Residual}}, \tau(\mathbf{b} \cdot \Nabla v_h)  )_{\Omega} = (f,v_h)_{\Omega} \quad \forall v_h \in \mathcal{V}_h(\Omega),
\end{equation}
\eirik{where the trial and test spaces consist of standard piecewise polynomials. Note that when the residual is zero, the stabilization term vanishes, i.e., it is consistent with the weak form~\eqref{eqn:weak_form3}.}

\eirik{As for the preceding model problems,} we wish to construct the 4 dimensional tensor $A_{ijkl}$ using the product bases of the test and trial spaces. \eirik{Substitution of} the product basis into (\ref{eqn:weak_form3b}) gives:
\begin{equation}
        \kappa(\Nabla (\phi_i\,\psi_j) , \Nabla (\gamma_k\,\beta_l))_{\Omega} + (\mathbf{b} \cdot \Nabla (\phi_i\, \psi_j), \gamma_k \beta_l)_{\Omega}
        + (-\kappa\Delta (\phi_i\,\psi_j) + \mathbf{b} \cdot \Nabla (\phi_i\,\psi_j) -f, \tau[\mathbf{b} \cdot \Nabla(\gamma_k\,\beta_l)]  )_{\Omega} = (f,\,\gamma_k\, \beta_l)_{\Omega}.
\end{equation}
\eirik{Application of} the product rule gives:
\begin{equation}
 \begin{split}
         \kappa(\psi_j \Nabla \phi_i + \phi_i \Nabla \psi_j , \beta_l \Nabla \gamma_k + \gamma_k \Nabla \beta_l)_{\Omega} + (\mathbf{b} \cdot (\psi_j \Nabla \phi_i +\phi_i \Nabla \psi_j), \gamma_k \beta_l)_{\Omega} +\\
         (-\kappa(\psi_j \Delta (\phi_i) + \phi_i \Delta (\psi_j)) + \mathbf{b} \cdot (\psi_j \Nabla (\phi_i) + \phi_i \Nabla(\psi_j)) -f,
         \tau[\mathbf{b} \cdot (\beta_l \Nabla (\gamma_k) + \gamma_k \Nabla (\beta_l))]  )_{\Omega} = (f \gamma_k  \beta_l)_{\Omega}, 
 \end{split}
 \end{equation}
\eirik{which we expand} using Fubini's theorem:
\begin{equation} \label{eq:fubinisplit}
\begin{split}
    \kappa \left[ (\Nabla \phi_i, \Nabla \gamma_k)_{\Omega_1}(\psi_j,\beta_l)_{\Omega_2} +  (\phi_i,\gamma_k)_{\Omega_1}( \Nabla \psi_j ,  \Nabla \beta_l)_{\Omega_2}\right] + (\mathbf{b} \cdot \Nabla \phi_i, \gamma_k)_{\Omega_1} (\psi_j,\beta_l)_{\Omega_2}
    + (\phi_i,\gamma_k)_{\Omega_1}(\mathbf{b} \cdot \Nabla \psi_j,\beta_l)_{\Omega_2} - \\
    \kappa \tau [(\Delta \phi_i, \mathbf{b} \cdot \Nabla \gamma_k)_{\Omega_1} (\psi_j \beta_l)_{\Omega_2}
    + (\Delta \phi_i , \gamma_k)_{\Omega_1}(\psi_j, \mathbf{b} \cdot  \Nabla \beta_l)_{\Omega_2} +    (\phi_i, \mathbf{b} \cdot \Nabla \gamma_k)_{\Omega_1}(\Delta \psi_j, \beta_l)_{\Omega_2} + (\phi_i,\gamma_k)_{\Omega_1}(\Delta \psi_j, \mathbf{b} \cdot \Nabla \beta_l)_{\Omega_2}] + \\
    \tau [ (\mathbf{b} \cdot \Nabla \phi_i, \mathbf{b} \cdot \Nabla \gamma_k)_{\Omega_1}(\psi_j,\beta_l)_{\Omega_2} + (\mathbf{b} \cdot \Nabla \phi_i, \gamma_k)_{\Omega_1}(\psi_j , \mathbf{b} \cdot \Nabla \beta_l)_{\Omega_2} + \\
    (\phi_i, \mathbf{b} \cdot \Nabla \gamma_k)_{\Omega_1}(\mathbf{b} \cdot \Nabla \psi_j, \beta_l )_{\Omega_2} + (\phi_i,\gamma_k)_{\Omega_1}(\mathbf{b} \cdot \Nabla \psi_j, \mathbf{b} \cdot \Nabla \beta_l)_{\Omega_2}] \\
    = f_{dof}(\gamma_k)_{\Omega_1}(\beta_l)_{\Omega_2} - \tau f_{dof}(\mathbf{b} \cdot \gamma_k)_{\Omega_1}(\beta_l)_{\Omega_2} - \tau f_{dof} (\gamma_k)_{\Omega_1}(\mathbf{b} \cdot \Nabla \beta_l)_{\Omega_2},
\end{split}
\end{equation}
\eirik{where $f_{dof}$ represents the function value $f$ at each corresponding degree of freedom .} The discrete weak form in~\eqref{eq:fubinisplit} can be represented as a global stiffness matrix compromised of the following smaller submatrices:
\begin{equation}
    \begin{split}
        K_{11} = (\Nabla \phi_i, \Nabla \gamma_k)_{\Omega_1} \quad  &
        K_{21} = ( \Nabla \psi_j ,  \Nabla \beta_l)_{\Omega_2}\\
        K_{12} = (\phi_i,\gamma_k)_{\Omega_1} \quad  &
        K_{22} = (\psi_j, \beta_l)_{\Omega_2}\\
        K_{13} = (\mathbf{b} \cdot \Nabla \phi_i,\gamma_k)_{\Omega_1}\quad  &
        K_{23} = (\mathbf{b} \cdot \Nabla \psi_j,\beta_l)_{\Omega_2}\\
        K_{14} = (\Delta \phi_i, \mathbf{b} \cdot \Nabla \gamma_k)_{\Omega_1}\quad  &
        K_{24} = (\Delta \psi_j, \mathbf{b} \cdot \Nabla \beta_l)_{\Omega_2}\\
        K_{15} = (\Delta \phi_i , \gamma_k)_{\Omega_1} \quad  &
        K_{25} = (\Delta \psi_j , \beta_l)_{\Omega_2}  \\
        K_{16} = (\phi_i, \mathbf{b} \cdot  \Nabla \gamma_k)_{\Omega_1}\quad  &
        K_{26} = (\psi_j, \mathbf{b} \cdot  \Nabla \beta_l)_{\Omega_2}\\
        K_{17} = (\mathbf{b} \cdot \Nabla \phi_i, \mathbf{b} \cdot \Nabla \gamma_k)_{\Omega_1}\quad  &
        K_{27} = (\mathbf{b} \cdot \Nabla \psi_j, \mathbf{b} \cdot \Nabla \beta_l)_{\Omega_2}\\
        F_{1} = (\gamma_k)_{\Omega_1} \quad  &
        F_{2} = (\beta_l)_{\Omega_2} \\
        F_{3} = (\mathbf{b} \cdot \Nabla \beta_l)_{\Omega_2} \quad  &
        F_{4} = (\mathbf{b} \cdot \Nabla \gamma_k)_{\Omega_1}.
    \end{split}
\end{equation}
Then the global system of equations to solve would be:
\begin{equation}
\begin{split}
    [ \kappa (K_{11} \otimes K_{22} + K_{12}\otimes K_{21}) + K_{13} \otimes K_{22} + K_{12} \otimes K_{23} - \kappa \tau(K_{14} \otimes K_{22} + K_{15} \otimes K_{26} + K_{16} \otimes K_{25} + K_{12} \otimes K_{24} ) + \\
    \tau ( K_{17} \otimes K_{22} + K_{13} \otimes K_{26} + K_{16} \otimes K_{23} + K_{12} \otimes K_{27} )]\mathbf{U} = f_{dof} F_1 \otimes F_2 - \tau f_{dof}(F_4 \otimes F_2 + F_1 \otimes F_3).
\end{split}
\end{equation}
One advantage of this method is that it is possible to align one subdomain with the velocity vector. \eirik{In this special case} when $\mathbf{b}$ is only non-zero along one subdomain (e.g., $\Omega_2$) the above weak form~\eqref{eq:fubinisplit} reduces to:
\begin{equation}\label{eqn:weak_form3c}
\begin{split}
    \kappa \left[ (\Nabla \phi_i, \Nabla \gamma_k)_{\Omega_1}(\psi_j,\beta_l)_{\Omega_2} +  (\phi_i,\gamma_k)_{\Omega_1}( \Nabla \psi_j ,  \Nabla \beta_l)_{\Omega_2}\right]+ (\phi_i,\gamma_k)_{\Omega_1}(\mathbf{b} \cdot \Nabla \psi_j,\beta_l)_{\Omega_2} - \\
    \kappa \tau\left[ (\Delta \phi_i , \gamma_k)_{\Omega_1}(\psi_j, \mathbf{b} \cdot  \Nabla \beta_l)_{\Omega_2} + (\phi_i,\gamma_k)_{\Omega_1}(\Delta \psi_j, \mathbf{b} \cdot \Nabla \beta_l)_{\Omega_2}\right] \\
    + \tau \left[   (\phi_i,\gamma_k)_{\Omega_1}(\mathbf{b} \cdot \Nabla \psi_j, \mathbf{b} \cdot \Nabla \beta_l)_{\Omega_2}\right] = f_{dof}(\gamma_k)_{\Omega_1}(\beta_l)_{\Omega_2}  - \tau f_{dof} (\gamma_k)_{\Omega_1}(\mathbf{b} \cdot \Nabla \beta_l)_{\Omega_2}.
\end{split}
\end{equation}
\eirik{Hence, the local matrices are defined:}
\begin{equation}
    \begin{split}
        K_{11} = (\Nabla \phi_i, \Nabla \gamma_k)_{\Omega_1} \quad  &
        K_{21} = ( \Nabla \psi_j ,  \Nabla w)_{\Omega_2}\\
        K_{12} = (\phi_i,\gamma_k)_{\Omega_1}\quad  &
        K_{22} = (\psi_j, \beta_l)_{\Omega_2}\\
         K_{13} = (\Delta \phi_i , \gamma_k)_{\Omega_1} \quad  &
        K_{23} = (\mathbf{b} \cdot \Nabla \psi_j,\beta_l)_{\Omega_2}\\
        K_{24} = (\psi_j, \mathbf{b} \cdot  \Nabla \beta_l)_{\Omega_2}\quad  &
        K_{25} = (\Delta \psi_j, \mathbf{b} \cdot \Nabla \beta_l)_{\Omega_2}\\
        K_{26} = (\mathbf{b} \cdot \Nabla \psi_j, \mathbf{b} \cdot \Nabla \beta_l)_{\Omega_2}\quad  &
        F_{1} = (\gamma_k)_{\Omega_1} \\
        F_{2} = (\beta_l)_{\Omega_2} \quad  &
        F_{3} = (\mathbf{b} \cdot \Nabla \beta_l)_{\Omega_2},
    \end{split}
\end{equation}
\eirik{substitution}  into~\eqref{eqn:weak_form3c} yields the global system:
\begin{equation} \label{eqn:discrete_AD_1direction}
    (\kappa K_{11} \otimes K_{22} + \kappa K_{12} \otimes K_{21} + K_{12} \otimes K_{23} -\kappa \tau K_{13} \otimes K_{24} - \kappa \tau K_{12} \otimes K_{25} + \tau K_{12} \otimes K_{26})\mathbf{U} = f_{dof} F_1 \otimes F_2 - \tau f_{dof} F_1 \otimes F_3.
\end{equation}

\section{Numerical Verifications}\label{sec:numerical_results}
For each of the three problem \eirik{introduced in Section~\ref{sec:methods}, we consider and implement a specific test case} in FEniCS. \eirik{To verify the developed framework, we investigate the $h-$convergence properties of the implemented methods by consideration of the  rate of convergence of the FE solutions. The test and trial spaces in all cases  consist of continuous Lagrange polynomials of degree 1.}

\subsection{Case 1: 4-D Poisson Equation }

\eirik{ We first consider the Poisson PDE in the high dimensional space of order four. We define the computational domain as } 
a tensor product between two unit squares, \eirik{i.e.,}
$\Omega = ((0,1)\times(0,1)) \times ((0,1) \times (0,1))$.
\eirik{We select} the forcing function defined as $f = \text{4}\pi^2 u_{exact}$  and the exact solution as $u_{exact} = \Pi^4_{i=1} sin(\pi x_i) $,
 where $x_i$ is a coordinate in the domain $\Omega$.  

 \eirik{In Table~\ref{tab:results_poisson}, the convergence data for the four  dimensional Poisson problem is listed. Note that the convergence of the FE solution is optimal, as  the rate of convergence of the root mean square error (RMSE) and $l_{\infty}$ norm approaches $\mathcal{O}(h^{p+1})$}.
\begin{table}[h!]
\centering
\caption{\label{tab:results_poisson}  Error estimation results for the 4D Poisson problem.}
\begin{tabular}{@{}llllll@{}}
\toprule
{dofs \hspace{6mm}} & { h \hspace{6mm}} & {$l_{\infty}$ \hspace{6mm}} & {  RMSE  \hspace{6mm}} & {RMSE rate} & {$l_{\infty}$ rate} \\
\midrule \midrule

256 &  0.333    & 1.40E-01  & 3.16E-02  & - & - \\
1296 & 0.200    & 7.88E-02  & 1.61E-02  & 1.44 & 1.13 \\
2401 & 0.167    & 6.55E-02  & 1.21E-02  & 1.57 & 1.01  \\
4096 & 0.143    & 4.50E-02  & 9.39E-03  & 1.65 & 2.44 \\
6561 & 0.125    & 3.76E-02  & 7.48E-03  & 1.70 & 1.34\\
1461 & 0.100   & 2.42E-02  & 5.06E-03  & 1.75 & 1.96 \\
20736 & 0.091   & 1.95E-02  & 4.26E-03  & 1.79 & 2.31 \\
28561 & 0.083   & 1.69E-02  & 3.64E-03  & 1.81 & 1.61 \\
38416 & 0.077   & 1.41E-02  & 3.15E-03  & 1.83  & 2.26 \\
50625 & 0.071   & 1.25E-02  & 2.74E-03. & 1.84  & 1.67  \\



\bottomrule
\end{tabular}
\end{table}



\subsection{Case 2: 2D Space-Time Wave Equation}

\eirik{As a verification of the space-time wave model problem}~\eqref{eqn:model_prob2}, \eirik{we select} the space-time domain as the Cartesian product of a unit square \eirik{spatial} domain $\Omega$ and an interval time domain, i.e., $\Omega_T = ((0,1) \times (0,1)) \times (0,T)$.  The wave propagation speed is $c=1$ and \eirik{we consider a}  manufactured solution $u_{exact} = sin(x-ct) + sin(y-ct)$. \eirik{This solution is used to ascertain boundary and initial conditions needed to solve the resulting system of equations.}
\eirik{In Table~\ref{tab:wave_spacetime}, the convergence results are presented along with the time interval element size, denoted by dt, and the space-time CFL number. The RMSE is observed to converge at the expected optimal rate, whereas the  }  $l_{\infty}$ error \eirik{exhibits a reduced rate for the finer meshes}.
\begin{table}[h!]
\centering
\caption{\label{tab:wave_spacetime}  Results for the Wave Space-Time problem.}
\begin{tabular}{@{}llllllll@{}}
\toprule
{dofs \hspace{6mm}} & { h \hspace{6mm}} & { dt \hspace{6mm}} & { CFL \hspace{6mm}} & {$l_{\infty}$ \hspace{6mm}} & {  RMSE  \hspace{15mm}} & {RMSE rate} & {$l_{\infty}$ rate} \\
\midrule \midrule

200 &  0.250 & 0.14 & 0.57  & 2.25E-04  & 5.84E-05 & - & -   \\
1215 & 0.125 & 0.07 & 0.57  & 5.06E-05  & 1.51E-05 & 1.95 & 2.22  \\
3718 & 0.083 & 0.05 & 0.57  & 2.53E-05  & 7.01E-06 & 1.89 & 1.71 \\
8381 & 0.063 & 0.04 & 0.57  & 1.31E-05  & 3.95E-06 & 1.99 & 1.31 \\


\bottomrule
\end{tabular}
\end{table}


\subsection{Case 3: SUPG Stabilized Advection Dominated Advection Diffusion Equation}
\eirik{As a final numerical verification, we consider a special case of advection dominated advection diffusion equation. In particular, we consider the case in which the advection 
acts in a single direction aligned with a coordinate axis, see~\eqref{eqn:discrete_AD_1direction}.    }
 \eirik{We consider a case where the} domain $\Omega$ is a tensor product of 2 unit intervals: $\Omega = (0,1) \times (0,1)$, diffusivity constant $\kappa =\frac{1}{100}$, and \eirik{advection  vector is} $\mathbf{b} = (0,1)$. The analytic solution in this case \eirik{is inspired by the work of Egger and Sch{\"o}berl~\cite{egger2010hybrid}}:
\begin{equation} \label{eq:ADR_exact}
u_{exact} =  (-4(x-0.5)^2+1)\left[ y + \frac{e^{\frac{1}{\kappa} \cdot b_y \cdot y} - 1 }{1 - e^{\frac{1}{\kappa} \cdot b_y}}\right],
\end{equation}
which implies that the forcing term must be $f = -\kappa\,\Delta u_{exact} + \mathbf{b} \cdot \Nabla u_{exact}$ \eirik{and we enforce the corresponding} homogeneous Dirichlet on the boundary $\partial \Omega$.
\eirik{In Figures~\ref{fig:cartesian_ADR}, and~\ref{fig:analytic_ADR}, we show the approximate FE solution and the analytic solution, respectively. As expected, the SUPG stabilization results in a stable solution at this relatively coarse FE mesh.}
\eirik{In Table~\ref{tab:results_ADR}, the convergence results for this final case are presented. Both the RMSE and $l_{\infty}$ error converge at the expected optimal rates.} 
\begin{table}[h]
\centering
\caption{\label{tab:results_ADR}  Error estimation results for advection dominated advection diffusion problem.}
\begin{tabular}{@{}llllll@{}}
\toprule
{dofs \hspace{6mm}} & { h \hspace{6mm}} & {$l_{\infty}$ \hspace{6mm}} & {  RMSE  \hspace{15mm}} & {RMSE rate} & {$l_{\infty}$ rate} \\
\midrule \midrule

9    & 0.50 & 2.50E-01   & 8.33E-02 & -  & - \\
16   & 0.33 & 1.19E-01   & 4.65E-02 & 1.44 & 1.84 \\
25   & 0.25 & 7.38E-02   & 3.01E-02 & 1.51 & 1.65 \\
36   & 0.20 & 4.79E-02   & 2.08E-02 & 1.65 & 1.93  \\
49   & 0.17 & 3.33E-02   & 1.52E-02 & 1.73 & 1.99 \\
256  & 0.07 & 6.41E-03   & 2.85E-03 & 1.83 & 1.80 \\
676  & 0.04 & 2.52E-03   & 1.07E-03 & 1.92 & 1.83 \\
2601 & 0.02 & 6.50E-04   & 2.74E-04 & 1.96 & 1.96 \\





\bottomrule
\end{tabular}
\end{table}
\begin{figure}[h!]
    \centering
    \includegraphics[width=0.5\textwidth]{Analytic_256_dof.png}
    \caption{Analytic solution projected onto FE mesh with 256 dofs.}
    \label{fig:analytic_ADR}
\end{figure}
\begin{figure}[h!]
    \centering
    \includegraphics[width=0.5\textwidth]{Cartesian_Product_256_dof.png}
    \caption{Cartesian product FE solution at 256 dofs.}
    \label{fig:cartesian_ADR}
\end{figure}
\section{Conclusions}\label{sec:Conclusions}

\eirik{In this paper, we have introduced and implemented tensor product FE routines  for high dimensional problems in FEniCS. This methodology allows us to extend the}
FEniCS library to domains with more than three dimensions so long as they are a Cartesian product of subdomains three or lower.
\eirik{To verify the developed methodology, we consider three test cases utilizing classical and stabilized FE methods. } For each test case, \eirik{we observe the expected } convergence to the analytic solutions with respect to grid refinement was demonstrated in both the RMSE and $l_{\infty}$ norms.

\eirik{We consider only linear PDEs here }  since this allowed for the explicit construction of a single linear system of algebraic equations. \eirik{Hence, future studies should investigate potential extensions to nonlinear PDEs.} Additionally, the global system of equations was solved naively by explicitly constructing the global stiffness matrix and inverting. \eirik{However,} the global matrix is sparse with highly structured blocks which should allow for faster solvers that would greatly reduce run time. \eirik{We refer to existing works~\cite{BIALECKI1993369,Gao_thesis,Bank1978}, where  related}  problems  were considered \eirik{and leave the consideration of such solvers for future studies.} 
\eirik{Further extensions to mixed FE methods, } such as those discussed in the book by Brezzi \emph{et al.}~\cite{brezzi2012mixed}. \eirik{Finally, } FE methods utilizing with discontinuous test/trial spaces \eirik{ could be considered due to their extensive use in engineering applications}. Furthermore, full integration of this method into the FEniCS API would be valuable for both simplicity of future implementations as well as for performance.

 

The code used in the numerical studies is publicly available and can be found on GitHub at \\ \url{https://github.com/Markloveland/Cartesian_Product_FEM} 


\section{Acknowledgements}

Author Loveland has been supported by the CSEM Fellowship from the Oden Institute at the University of Texas at Austin. Authors Loveland, Valseth, and Dawson have been supported by the United States National Science Foundation - NSF PREEVENTS Track
2 Program, under NSF Grant Number 1855047 and the Department of Homeland Security Coastal Resilience Center research project "Accurate and Fast Wave Modeling and Coupling with ADCIRC".
Author Lukac has been supported by the University of Oregon.

 
\bibliographystyle{elsarticle-num}
 
\section{Introduction}

In the last decade, open source libraries with high level APIs that automate the process of solving partial differential equations (PDEs) using the finite element method (FEM) have become valuable tools in computational research. Some prominent libraries of this kind include Firedrake~\cite{rathgeber2016firedrake}, deal.II~\cite{bangerth2007deal}, MFEM~\cite{anderson2021mfem}, and  FEniCS~\cite{alnaes2015fenics}. FEniCS is one of the most widely used of these libraries. FEniCS, along with most other FEM libraries, only include capability of handling up to three dimensional problems (GetFEM~\cite{renard2020getfem} is a notable exception). FEniCS specifically supports unstructured meshes with corresponding basis functions in one, two, and three dimensions. Extension of FEniCS to solve problems in higher dimensions is of interest for problems such as spectral wind wave models, spatial population genetics, quantum mechanics, and even 3d space-time problems.

FEniCS does not contain support for fully unstructured meshes in dimensions higher than three. However, at the expense of losing a fully unstructured mesh, a higher dimensional space can be discretized with the available tools in FEniCS if it is the Cartesian product of two lower dimensional (3 or lower) spaces that can themselves be unstructured. It is well known that in a Cartesian product space, a finite element basis can be constructed called the product basis which spans the function space of the full domain~\cite{brenner2008mathematical,quarteroni2010numerical,ern2013theory}. This research seeks to use the infrastructure of an FEM library, such as FEniCS, in the lower dimensional spaces to construct a basis for this high dimensional space. 

 Using FEM to discretize a structured domain via a product basis has been done before without the use of FEM libraries such as FEniCS. In 1978, Banks~\cite{Bank1978} used tensor product finite elements to solve a 2D Poisson problem on a structured grid in order to find a faster solver. In the study, the domain was a Cartesian product between two unit intervals resulting in a uniform square mesh. Banks used the tensor product of two one dimensional quadratic and cubic basis functions to construct the polynomial basis of the full 2D domain. In 1980, Baker performed numerical tests as well as a stability analysis on a tensor product finite method with application to convection-dominated fluid flow problems~\cite{BAKER1981215}. The stability analysis showed the basic algorithm is spatially fourth- order accurate in its most elementary embodiment and the numerical experiments on the convection-dominated model test problems confirmed the basic viability of the developed algorithm, and its tensor product formulation. Recently, Du \eirik{\emph{et al.}~\cite{DU2013181}} used tensor product finite elements to construct a fast solver for an electromagnetics scattering problem of a cavity. In  Firedrake~\cite{rathgeber2016firedrake}, there exists a capability to create tensor product finite elements only up to three dimensions~\cite{mcrae2016automated}. A package built on deal.II called Hyperdeal~\cite{munch2020hyperdeal} has the capability of creating tensor product finite elements in up to six dimensions. deal.II is different to a library such as FEniCS since deal.II uses quadrilateral elements in two dimensions whereas FEniCS uses triangles. \eirik{Here, we present a general software framework built on the components of the FEniCS library as well as other open source Python libraries to construct high dimensional meshes and corresponding FE discretizations.  }


Following this introduction, the general set of problems for which this paper focuses on will be defined and the notation for the product basis will be introduced \eirik{in Section~\ref{sec:model_prob}}. Then, a set of three different model problems will be described in detail and the derivation of the system of equations using tensor product finite elements will be shown for each case. The first model problem will be discussed in Section~\ref{sec:model_prob} which is an N dimensional Poisson problem where the domain is a Cartesian product between two subdomains of dimension three or less. In Section~\ref{sec:prob2} the second model problem is discussed which is a wave equation where the domain is decomposed as a Cartesian product between space and time. The third model problem is discussed in Section~\ref{sec:prob3} which is an advection dominated advection-diffusion equation where the advection aligns with one of the subdomains, a streamlined upwind Petrov-Galerkin method (SUPG) is formulated for this case. After each system of equations is derived, numerical tests were run using FEniCS for each model problem where specific boundary conditions were given. For each case, error and convergence rates are tabulated in Section~\ref{sec:numerical_results}. Lastly, conclusions and recommendations for future work are given in Section~\ref{sec:Conclusions}.


\section{Methods} \label{sec:methods}

To present the proposed methodology and algorithms, we 
consider the following class of problems, i.e., PDEs:
\begin{equation*}
\begin{split}
    \mathscr{L} u = f \quad \textrm{in} \quad \Omega = \Omega_1 \times \Omega_2.
\end{split}
\end{equation*}
Where $\mathscr{L}$ is a linear differential operator, $f$ is a forcing function, and the domain $\Omega$ is defined as a Cartesian product between two lower dimensional Lipschitz domains.  For example, if the global domain $\Omega  \subset \mathbb{R}^2$, then $\Omega$ can be defined by the Cartesian product of 2 intervals $\Omega_1\subset \mathbb{R}$ and  $\Omega_2 \subset \mathbb{R}$. Hence, we  consider general domains of the form $\Omega = \Omega_1 \times \Omega_2 $, see Figure~\ref{fig:cartesian_illustration} for an illustration.
\begin{figure}[h!]
    \centering
    \includegraphics[width=\textwidth]{CartesianProductIllust.png}
    \caption{Example of a global domain $\Omega$ that is a Cartesian product of two lower dimensional subdomains $\Omega_1, \Omega_2$.}
    \label{fig:cartesian_illustration}
\end{figure}
FEniCS can discretize up to three dimensional objects, hence,  in practice this framework can be used to define domains that are Cartesian products of up to six dimensions. Furthermore, this process can be done iteratively to yield even higher dimensional domains as $\Omega_1$ itself could be a product of two other spaces and so on. We note that in this presentation, we consider only symmetric functional settings which admit well posed  Galerkin FE discrtetization. 

It can be shown that if a function $u$ is defined on a Cartesian product domain $\Omega = \Omega_1 \times \Omega_2$ then we can construct a basis of a polynomial space to be used for the FEM \eirik{ by exploiting this underlying geometric structure, see, e.g.,}~\cite{quarteroni2010numerical,brenner2008mathematical,ern2013theory}. 
 \eirik{With} a basis for a polynomial space for the first lower dimensional subdomain $\Omega_1$:
\begin{equation}
    \{\phi_i\}_{i=1}^N,
\end{equation}
and a basis for the second subdomain $\Omega_2$:
\begin{equation}
    \{\psi_j\}_{j=1}^M.
\end{equation}
The basis for the entire domain can subsequently be constructed and any arbitrary function $u$ whose domain is in $\Omega = \Omega_1 \times \Omega_2$ can be approximated via the product basis:
\begin{equation}
    u \approx \sum_{i=1}^N\sum_{j=1}^M u_{i,j} \phi_i\, \psi_j. 
\end{equation}
\eirik{Note that in the following,} we use boldface letters and symbols to denote vector quantities, e.g., $\dX = \dx \dy$. \eirik{To avoid raising complexity in implementation, we also assume in the numerical test cases (Section~\ref{sec:numerical_results}) that all physical parameters in the following model problems are constant and that the Dirichlet boundary conditions are homogeneous.}
\matt{However, it is possible to remove these assumptions by projecting the parameters and boundary conditions to the appropriate function space.}

\subsection{Model Problem 1: N dimensional Poisson Equation} \label{sec:model_prob}
As a first model problem to illustrate the methodology, we consider the Poisson equation:
\begin{equation}\label{eqn:model_prob1}
\begin{split}
        - \Delta \, u &= f \quad \textrm{in} \quad \Omega, \\
        u &= u_D \quad \textrm{on} \quad \partial \Omega, 
\end{split}
\end{equation}
where the source $f$ is in $L^2(\Omega)$ and the source data $u_D$ is assumed to be sufficiently regular.
To define the weak formulation for~\eqref{eqn:model_prob1}, multiply both sides with a test function $v$ in $L^2(\Omega)$ and integrate over $\Omega$:
\begin{equation}\label{eqn:weak_form1}
    \int_{\Omega}  - \Delta \, u\, v \dX = \int_{\Omega} f \, v \dX \quad \forall v \in L^2(\Omega), 
\end{equation}
Integrating by parts on the left side (assuming Dirichlet boundary conditions on entire boundary) gives the following weak formulation\matt{: find $u \in \mathcal{U}(\Omega)$ such that}
\begin{equation} \label{eq:weak_form}
    \int_{\Omega}  \Nabla u \cdot \Nabla v \dX = \int_{\Omega} f\,v \dX \quad \forall \quad v \in \mathcal{V}(\Omega), 
\end{equation}
where the function space $\mathcal{V}(\Omega)$ is the Hilbert space with zero trace on the boundary  $H^1_0(\Omega)$ and $\mathcal{U}(\Omega)$ is $\mathcal{V}(\Omega)$ with a finite energy lift on the boundary so that the Dirichlet condition $u=u_D$ is satisfied.
The weak formulation in~\eqref{eq:weak_form}, and its corresponding discretization is known to be well posed, see, e.g.,~\cite{becker1981finite}. 

With a well posed weak formulation at hand, we can discretize this weak form using a finite element basis. In this case, the functional setting dictates the use of a $C^0$ continuous \eirik{polynomial} basis for the classical FEM.
Hence, we approximate the trial functions with the product basis:
\begin{equation}
    u \approx \sum_{i=1}^N\sum_{j=1}^M u_{i,j} \, \phi_i \, \psi_j,
\end{equation}
and the test functions by its product basis:
\begin{equation}
    v \approx \sum_{k=1}^N\sum_{l=1}^M v_{k,l} \, \gamma_k \, \beta_l.
\end{equation}
The weak form from~\eqref{eqn:weak_form1} can then be discretized by substitution of the product bases:
\begin{equation}
\displaystyle    \int_{\Omega}  (\Nabla \sum_{i=1}^N\sum_{j=1}^M u_{i,j} \, \phi_i \, \psi_j) \cdot ( \Nabla \sum_{k=1}^N\sum_{l=1}^M  v_{k,l}\,\gamma_k\, \beta_l) \dX = \int_{\Omega} f\, (\sum_{k=1}^N\sum_{l=1}^M  v_{k,l}\,\gamma_k\, \beta_l ) \dX.
\end{equation}
Due to arbitrariness of the test function and the bilinearity of the weak form, this implies the following form:
\begin{equation}
    \sum_{i=1}^N\sum_{j=1}^M \int_{\Omega}  (\Nabla   \, \phi_i \, \psi_j) \cdot ( \Nabla \,\gamma_k\, \beta_l) \dX \, u_{i,j} = \int_{\Omega} f\, (\,\gamma_k\, \beta_l ) \dX.
\end{equation}
This left hand side \eirik{results} in a \eirik{product between a} 4 dimensional \eirik{and} a 2 dimensional tensor $u_{i,j}$. \eirik{Consequently,} the 4 dimensional tensor $A_{ijkl}$ will be of the form:
\begin{equation}
        \int_{\Omega}   \Nabla  \phi_i \,\psi_j \cdot \Nabla  \gamma_k \,\beta_l \, \dX.
\end{equation}
By the product rule we have:
\begin{equation} \label{eq:chain_ruled}
        \int_{\Omega}   \Nabla  \phi_i \, \psi_j \cdot \Nabla  \gamma_k \,\beta_l \dX  = \int_{\Omega}\, (\phi_i \, \Nabla \psi_j + \psi_j\, \Nabla \phi_i) \cdot (\gamma_k\, \Nabla \beta_l + \beta_l\, \Nabla \gamma_k) \dX.
\end{equation}
By construction, \eirik{the} $\phi$\eirik{'s}  and$\gamma$\eirik{'s} only vary in the first subdomain $\Omega_1$ whereas \eirik{the} $\psi$\eirik{'s} and $\beta$\eirik{'s} only vary in the second subdomain $\Omega_2$. \eirik{Thus, the integral form~\eqref{eq:chain_ruled} can be simplified. To this end, we use the following notation convention:} if the gradient operator on the entire domain \eirik{$\Omega$} is $\Nabla$, define the gradient on the subdomains \eirik{$\Omega_1$ and $\Omega_2$} \eirik{as} $\Nabla_1$, $\Nabla_2$\eirik{, respectively}. \eirik{Hence, by} construction $\Nabla = (\Nabla_1, \Nabla_2)$.
Rewriting \eirik{ the gradients in~\eqref{eq:chain_ruled} then} gives:
\begin{equation} \label{eq:expansion}
\begin{split}
   & \int_{\Omega} (\phi_i \, (\Nabla_1, \Nabla_2)\,\psi_j + \psi_j \,(\Nabla_1, \Nabla_2)\, \phi_i) \cdot (\gamma_k\, (\Nabla_1, \Nabla_2)\, \beta_l + \beta_l (\Nabla_1, \Nabla_2)\gamma_k)\dX,  \\ \qquad &=
    \int_{\Omega} (  (\mathbf{0}, \phi_i\,\Nabla_2 \psi_j) +  (\psi_j\,\Nabla_1\phi_i, \mathbf{0}) ) \cdot ( (\mathbf{0}, \gamma_k\, \Nabla_2\beta_l)  +  (\beta_l\,\Nabla_1\gamma_k, \mathbf{0}))\dX,  \\ &=
    \int_{\Omega}  (\psi_j\,\Nabla_1\phi_i, \phi_i\,\Nabla_2 \psi_j)  \cdot  (\beta_l\,\Nabla_1\gamma_k, \gamma_k\, \Nabla_2\beta_l)  \dX, \\ \qquad &= 
    \int_{\Omega}  \Nabla_1\phi_i \cdot \Nabla_1\gamma_k\,\psi_j\,\beta_l +  \phi_i \,\gamma_k \,\Nabla_2 \psi_j \cdot \Nabla_2\,\beta_l  \dX. 
\end{split}
\end{equation}
Since \eirik{the} domain $\Omega$ is a Cartesian product of \eirik{the} subdomains $\Omega_1, \Omega_2$ we can rewrite the \eirik{last} integral \eirik{in~\eqref{eq:expansion},
as}:
\begin{equation}
           \int_{\Omega_1} \int_{\Omega_2}  \Nabla_1\phi_i \cdot \Nabla_1\gamma_k\,\psi_j\,\beta_l +  \phi_i\, \gamma_k\, \Nabla_2 \psi_j \cdot \Nabla_2\beta_l  \dy  \dx.
\end{equation}
\eirik{An application of} Fubini's theorem \eirik{gives}:
\begin{equation}
               \int_{\Omega_1} \Nabla_1\phi_i \cdot \Nabla_1\gamma_k \dx \int_{\Omega_2}  \psi_j\,\beta_l \dy + \int_{\Omega_1} \phi_i\, \gamma_k \dx \int_{\Omega_2}  \Nabla_2 \psi_j \cdot \Nabla_2\beta_l   \dy,
\end{equation}
These can be written as Kronecker products of  matrices,  \eirik{e.g., stiffness matrices computed with} FEniCS. \eirik{To continue the discussion, we introduce}  the following notation:
\begin{equation}
    \begin{split}
        K_{11} = \int_{\Omega_1} \Nabla_1\phi_i \cdot \Nabla_1\gamma_k \dx, \quad K_{22} = \int_{\Omega_2}  \psi_j\;\beta_l \dy,\\
        K_{12} = \int_{\Omega_1} \phi_i\; \gamma_k \dx, \quad
        K_{21} = \int_{\Omega_2}  \Nabla_2 \psi_j \cdot \Nabla_2\beta_l   \dy,
    \end{split}
\end{equation}
\eirik{where} the first index indicates subdomain (\eirik{i.e.,} 1 or 2 \eirik{in this case}) and the second indicates the term from the weak form (1 is the first order operator and 2  the \eirik{$0^{th}$} order). \eirik{Hence, each local matrix $K_{ij}$ can be constructed independently, and the global stiffness matrix $A_{ijkl}$ can be defined by the sum of 2 Kronecker products.}
The Kronecker product allows us to represent \eirik{the matrix of a} 4D tensor  \eirik{by smaller 2D matrices}:
\begin{equation}
    A = K_{11} \otimes K_{22} + K_{12} \otimes K_{21}.
\end{equation}
\eirik{A similar reasoning for the right hand side leads to}:
\begin{equation}
    \int_{\Omega}  f_{dof}\;\phi_k\; \beta_l \dX  = f_{dof} F_1 \otimes F_2
\end{equation}
Where $f_{dof}$ is the pointwise value of $f$ at each $k,l$ coordinate in the global space. The forcing vectors are defined as:
\begin{equation}
    F_{1} = \int_{\Omega_1}   \phi_k \dx \quad F_{2} = \int_{\Omega_2}  \beta_l \dy.
\end{equation}
\eirik{Finally,} the entire system \eirik{of equations} becomes:
\begin{equation} \label{eq:final_poisson_form}
    (K_{11} \otimes K_{22} + K_{12} \otimes K_{21})(U) = f_{dof} F_1 \otimes F_2,
\end{equation}
where $U$ are the values at all d.o.f in the Cartesian product space for $u$. 
\eirik{A subsequent application of boundary conditions to~\eqref{eq:final_poisson_form}  leads to the final system of linear algebraic equations.} 


\subsection{Model Problem 2: Arbitrary dimensional Space-Time Wave Equation}\label{sec:prob2}
\eirik{The second model problem we consider is the linear wave equation. We consider this transient problem to highlight the application of this methodology to transient problems where space-time finite elements are employed. Thus, we have the following model problem:}
\begin{equation}\label{eqn:model_prob2}
\begin{split}
        \frac{\partial^2 u}{\partial t}-c^2\Delta u & = 0 \quad \textrm{in} \quad \Omega_T, \\
        \eirik{u} &\eirik{ =u_D \quad \textrm{on} \quad \partial \Omega,} \\
        \eirik{u} &\eirik{ =u_{initial} \quad \textrm{on} \quad \Omega,}
\end{split}
\end{equation}
\eirik{where $c$ denotes the wave speed}.
In this problem the space-time domain $\Omega_T$ is a Cartesian product domain of one, two, or three dimensional \eirik{spatial domain and a one dimensional temporal domain} ($\Omega_T = \Omega \times (0,T)$). The weak form of~\eqref{eqn:model_prob2} is obtained by multiplying by a test function and integrating over the entire space-time domain:
\begin{equation}
    \int_{\Omega_T} \frac{\partial^2 u}{\partial t} v  - c^2\Delta u v \dX  \, {\rm d t}     = 0 \quad \forall \quad v \in L^2(\Omega_T),
\end{equation}
\eirik{and subsequent integration } by parts in \eirik{space and time}  gives the following weak form. . Find $u \in \mathcal{U}(\Omega_T)$:
\begin{equation}\label{eqn:weak_form2}
   \int_{\Omega_T} -\frac{\partial u}{\partial t} \frac{\partial v}{\partial t} + c^2\Nabla u \cdot \Nabla v \dX \, {\rm d t}     + \int_{\partial \Omega_T|_{t=T}}\frac{\partial u}{\partial t}v  \dX = 0 \quad \forall \quad v \in \mathcal{V}(\Omega_T),
\end{equation}
where \eirik{we have applied Dirichlet conditions to the space-time boundary $\partial \Omega_T$, except at the final time boundary.} $\mathcal{V}(\Omega_T)$ is the space of all $H^1_0(\Omega_T)$ functions except on $\partial \Omega_T|_{t=T}$, where the trace is an unkown, and $\mathcal{U}(\Omega_T)$ is the space of all $H^1_0(\Omega_T)$ plus the trace $u=u_D$ on $\partial \Omega_T$ except the aforementioned part of the boundary $\Omega_T|_{t=T}$. 

\eirik{The corresponding discretization of}~\eqref{eqn:weak_form2} using a product basis gives a very similar system of equations to the Poisson problem \eirik{considered in Section~\ref{sec:model_prob}}, with a slight variation. The term $K_{21}$ must include \eirik{additional integrals} and becomes:
\begin{equation}
    K_{21} = \int_0^T  \frac{\partial \psi_j}{\partial t} \frac{\partial \beta_l}{\partial t}   dt - \int_{\partial \Omega_2} \Nabla \psi_j \cdot \mathbf{n} \beta_l  \, {\rm d s}    .
\end{equation}
\eirik{Consequently}, the global system of equations becomes:
\begin{equation} \label{eq:global_wave_sys}
    (c^2 K_{11} \otimes K_{22} - K_{12} \otimes K_{21})U = \mathbf{0}.
\end{equation}
\eirik{Finally, application of boundary and initial conditions to~\eqref{eq:global_wave_sys} results in the final system of linear algebraic equations.} 

\subsection{Model Problem 3: SUPG Stabilized Advection Dominated Advection Diffusion Equation}\label{sec:prob3}

\eirik{To highlight the versatility of our approach to consider non-standard FE techniques, we consider a PDE which is known to lead to stability issues in the Galerkin FE setting. Hence, we consider an advection-diffusion  PDE  in which advection is the dominant:}  
\begin{equation}\label{eqn:model_problem3}
\begin{split}
        -\kappa \, \Delta u + \mathbf{b} \cdot \Nabla u = f, \quad \textrm{on} \quad \Omega \\
         u = u_D \quad \textrm{on} \quad \partial \Omega,,
\end{split}
\end{equation}
where $\Omega$ is defined as a Cartesian Product of two lower dimensional Lipschitz domains: $\Omega = \Omega_1 \times \Omega_2$. \eirik{By following the standard procedure of deriving integral formulations,} we get the corresponding weak form:
\begin{equation}\label{eqn:weak_form3}
        \int_{\Omega} \kappa \Nabla u \cdot \Nabla v + \mathbf{b} \cdot \Nabla u v dx = \int_{\Omega}fv dx \quad \forall v \in \mathcal{V}(\Omega)
\end{equation}
In problems where the advection term dominates the diffusion, that is when Peclet number $Pe = \frac{L \| \mathbf{b} \| }{\kappa} >> 1$, where $L$ is the characteristic length, the standard Galerkin \eirik{method applied to~\eqref{eqn:weak_form3} may result in a}   discretization that is \eirik{unstable. } \eirik{This issue of stability can be overcome by careful design of the FE mesh or through stabilization techniques that ensure satisfaction of the discrete \emph{inf-sup} condition. 
Here, we consider the SUPG method introduced by Brooks and Hughes~\cite{brooks1982streamline} since it is widely used and has well developed criteria for discrete stability.}
\eirik{The SUPG method leads to stable FE discretizations by adjusting the discretized weak form~\eqref{eqn:weak_form3} with a penalized residual, i.e., find $u_h \in \mathcal{U}_h(\Omega)$: }
\begin{equation}\label{eqn:weak_form3b}
   \kappa(\Nabla u_h , \Nabla v_h)_{\Omega} + (\mathbf{b} \cdot \Nabla u_h, v_h)_{\Omega} + (\underbrace{-\kappa\,\Delta u_h + \mathbf{b} \cdot \Nabla u_h -f}_{\eirik{Residual}}, \tau(\mathbf{b} \cdot \Nabla v_h)  )_{\Omega} = (f,v_h)_{\Omega} \quad \forall v_h \in \mathcal{V}_h(\Omega),
\end{equation}
\eirik{where the trial and test spaces consist of standard piecewise polynomials. Note that when the residual is zero, the stabilization term vanishes, i.e., it is consistent with the weak form~\eqref{eqn:weak_form3}.}

\eirik{As for the preceding model problems,} we wish to construct the 4 dimensional tensor $A_{ijkl}$ using the product bases of the test and trial spaces. \eirik{Substitution of} the product basis into (\ref{eqn:weak_form3b}) gives:
\begin{equation}
        \kappa(\Nabla (\phi_i\,\psi_j) , \Nabla (\gamma_k\,\beta_l))_{\Omega} + (\mathbf{b} \cdot \Nabla (\phi_i\, \psi_j), \gamma_k \beta_l)_{\Omega}
        + (-\kappa\Delta (\phi_i\,\psi_j) + \mathbf{b} \cdot \Nabla (\phi_i\,\psi_j) -f, \tau[\mathbf{b} \cdot \Nabla(\gamma_k\,\beta_l)]  )_{\Omega} = (f,\,\gamma_k\, \beta_l)_{\Omega}.
\end{equation}
\eirik{Application of} the product rule gives:
\begin{equation}
 \begin{split}
         \kappa(\psi_j \Nabla \phi_i + \phi_i \Nabla \psi_j , \beta_l \Nabla \gamma_k + \gamma_k \Nabla \beta_l)_{\Omega} + (\mathbf{b} \cdot (\psi_j \Nabla \phi_i +\phi_i \Nabla \psi_j), \gamma_k \beta_l)_{\Omega} +\\
         (-\kappa(\psi_j \Delta (\phi_i) + \phi_i \Delta (\psi_j)) + \mathbf{b} \cdot (\psi_j \Nabla (\phi_i) + \phi_i \Nabla(\psi_j)) -f,
         \tau[\mathbf{b} \cdot (\beta_l \Nabla (\gamma_k) + \gamma_k \Nabla (\beta_l))]  )_{\Omega} = (f \gamma_k  \beta_l)_{\Omega}, 
 \end{split}
 \end{equation}
\eirik{which we expand} using Fubini's theorem:
\begin{equation} \label{eq:fubinisplit}
\begin{split}
    \kappa \left[ (\Nabla \phi_i, \Nabla \gamma_k)_{\Omega_1}(\psi_j,\beta_l)_{\Omega_2} +  (\phi_i,\gamma_k)_{\Omega_1}( \Nabla \psi_j ,  \Nabla \beta_l)_{\Omega_2}\right] + (\mathbf{b} \cdot \Nabla \phi_i, \gamma_k)_{\Omega_1} (\psi_j,\beta_l)_{\Omega_2}
    + (\phi_i,\gamma_k)_{\Omega_1}(\mathbf{b} \cdot \Nabla \psi_j,\beta_l)_{\Omega_2} - \\
    \kappa \tau [(\Delta \phi_i, \mathbf{b} \cdot \Nabla \gamma_k)_{\Omega_1} (\psi_j \beta_l)_{\Omega_2}
    + (\Delta \phi_i , \gamma_k)_{\Omega_1}(\psi_j, \mathbf{b} \cdot  \Nabla \beta_l)_{\Omega_2} +    (\phi_i, \mathbf{b} \cdot \Nabla \gamma_k)_{\Omega_1}(\Delta \psi_j, \beta_l)_{\Omega_2} + (\phi_i,\gamma_k)_{\Omega_1}(\Delta \psi_j, \mathbf{b} \cdot \Nabla \beta_l)_{\Omega_2}] + \\
    \tau [ (\mathbf{b} \cdot \Nabla \phi_i, \mathbf{b} \cdot \Nabla \gamma_k)_{\Omega_1}(\psi_j,\beta_l)_{\Omega_2} + (\mathbf{b} \cdot \Nabla \phi_i, \gamma_k)_{\Omega_1}(\psi_j , \mathbf{b} \cdot \Nabla \beta_l)_{\Omega_2} + \\
    (\phi_i, \mathbf{b} \cdot \Nabla \gamma_k)_{\Omega_1}(\mathbf{b} \cdot \Nabla \psi_j, \beta_l )_{\Omega_2} + (\phi_i,\gamma_k)_{\Omega_1}(\mathbf{b} \cdot \Nabla \psi_j, \mathbf{b} \cdot \Nabla \beta_l)_{\Omega_2}] \\
    = f_{dof}(\gamma_k)_{\Omega_1}(\beta_l)_{\Omega_2} - \tau f_{dof}(\mathbf{b} \cdot \gamma_k)_{\Omega_1}(\beta_l)_{\Omega_2} - \tau f_{dof} (\gamma_k)_{\Omega_1}(\mathbf{b} \cdot \Nabla \beta_l)_{\Omega_2},
\end{split}
\end{equation}
\eirik{where $f_{dof}$ represents the function value $f$ at each corresponding degree of freedom .} The discrete weak form in~\eqref{eq:fubinisplit} can be represented as a global stiffness matrix compromised of the following smaller submatrices:
\begin{equation}
    \begin{split}
        K_{11} = (\Nabla \phi_i, \Nabla \gamma_k)_{\Omega_1} \quad  &
        K_{21} = ( \Nabla \psi_j ,  \Nabla \beta_l)_{\Omega_2}\\
        K_{12} = (\phi_i,\gamma_k)_{\Omega_1} \quad  &
        K_{22} = (\psi_j, \beta_l)_{\Omega_2}\\
        K_{13} = (\mathbf{b} \cdot \Nabla \phi_i,\gamma_k)_{\Omega_1}\quad  &
        K_{23} = (\mathbf{b} \cdot \Nabla \psi_j,\beta_l)_{\Omega_2}\\
        K_{14} = (\Delta \phi_i, \mathbf{b} \cdot \Nabla \gamma_k)_{\Omega_1}\quad  &
        K_{24} = (\Delta \psi_j, \mathbf{b} \cdot \Nabla \beta_l)_{\Omega_2}\\
        K_{15} = (\Delta \phi_i , \gamma_k)_{\Omega_1} \quad  &
        K_{25} = (\Delta \psi_j , \beta_l)_{\Omega_2}  \\
        K_{16} = (\phi_i, \mathbf{b} \cdot  \Nabla \gamma_k)_{\Omega_1}\quad  &
        K_{26} = (\psi_j, \mathbf{b} \cdot  \Nabla \beta_l)_{\Omega_2}\\
        K_{17} = (\mathbf{b} \cdot \Nabla \phi_i, \mathbf{b} \cdot \Nabla \gamma_k)_{\Omega_1}\quad  &
        K_{27} = (\mathbf{b} \cdot \Nabla \psi_j, \mathbf{b} \cdot \Nabla \beta_l)_{\Omega_2}\\
        F_{1} = (\gamma_k)_{\Omega_1} \quad  &
        F_{2} = (\beta_l)_{\Omega_2} \\
        F_{3} = (\mathbf{b} \cdot \Nabla \beta_l)_{\Omega_2} \quad  &
        F_{4} = (\mathbf{b} \cdot \Nabla \gamma_k)_{\Omega_1}.
    \end{split}
\end{equation}
Then the global system of equations to solve would be:
\begin{equation}
\begin{split}
    [ \kappa (K_{11} \otimes K_{22} + K_{12}\otimes K_{21}) + K_{13} \otimes K_{22} + K_{12} \otimes K_{23} - \kappa \tau(K_{14} \otimes K_{22} + K_{15} \otimes K_{26} + K_{16} \otimes K_{25} + K_{12} \otimes K_{24} ) + \\
    \tau ( K_{17} \otimes K_{22} + K_{13} \otimes K_{26} + K_{16} \otimes K_{23} + K_{12} \otimes K_{27} )]\mathbf{U} = f_{dof} F_1 \otimes F_2 - \tau f_{dof}(F_4 \otimes F_2 + F_1 \otimes F_3).
\end{split}
\end{equation}
One advantage of this method is that it is possible to align one subdomain with the velocity vector. \eirik{In this special case} when $\mathbf{b}$ is only non-zero along one subdomain (e.g., $\Omega_2$) the above weak form~\eqref{eq:fubinisplit} reduces to:
\begin{equation}\label{eqn:weak_form3c}
\begin{split}
    \kappa \left[ (\Nabla \phi_i, \Nabla \gamma_k)_{\Omega_1}(\psi_j,\beta_l)_{\Omega_2} +  (\phi_i,\gamma_k)_{\Omega_1}( \Nabla \psi_j ,  \Nabla \beta_l)_{\Omega_2}\right]+ (\phi_i,\gamma_k)_{\Omega_1}(\mathbf{b} \cdot \Nabla \psi_j,\beta_l)_{\Omega_2} - \\
    \kappa \tau\left[ (\Delta \phi_i , \gamma_k)_{\Omega_1}(\psi_j, \mathbf{b} \cdot  \Nabla \beta_l)_{\Omega_2} + (\phi_i,\gamma_k)_{\Omega_1}(\Delta \psi_j, \mathbf{b} \cdot \Nabla \beta_l)_{\Omega_2}\right] \\
    + \tau \left[   (\phi_i,\gamma_k)_{\Omega_1}(\mathbf{b} \cdot \Nabla \psi_j, \mathbf{b} \cdot \Nabla \beta_l)_{\Omega_2}\right] = f_{dof}(\gamma_k)_{\Omega_1}(\beta_l)_{\Omega_2}  - \tau f_{dof} (\gamma_k)_{\Omega_1}(\mathbf{b} \cdot \Nabla \beta_l)_{\Omega_2}.
\end{split}
\end{equation}
\eirik{Hence, the local matrices are defined:}
\begin{equation}
    \begin{split}
        K_{11} = (\Nabla \phi_i, \Nabla \gamma_k)_{\Omega_1} \quad  &
        K_{21} = ( \Nabla \psi_j ,  \Nabla w)_{\Omega_2}\\
        K_{12} = (\phi_i,\gamma_k)_{\Omega_1}\quad  &
        K_{22} = (\psi_j, \beta_l)_{\Omega_2}\\
         K_{13} = (\Delta \phi_i , \gamma_k)_{\Omega_1} \quad  &
        K_{23} = (\mathbf{b} \cdot \Nabla \psi_j,\beta_l)_{\Omega_2}\\
        K_{24} = (\psi_j, \mathbf{b} \cdot  \Nabla \beta_l)_{\Omega_2}\quad  &
        K_{25} = (\Delta \psi_j, \mathbf{b} \cdot \Nabla \beta_l)_{\Omega_2}\\
        K_{26} = (\mathbf{b} \cdot \Nabla \psi_j, \mathbf{b} \cdot \Nabla \beta_l)_{\Omega_2}\quad  &
        F_{1} = (\gamma_k)_{\Omega_1} \\
        F_{2} = (\beta_l)_{\Omega_2} \quad  &
        F_{3} = (\mathbf{b} \cdot \Nabla \beta_l)_{\Omega_2},
    \end{split}
\end{equation}
\eirik{substitution}  into~\eqref{eqn:weak_form3c} yields the global system:
\begin{equation} \label{eqn:discrete_AD_1direction}
    (\kappa K_{11} \otimes K_{22} + \kappa K_{12} \otimes K_{21} + K_{12} \otimes K_{23} -\kappa \tau K_{13} \otimes K_{24} - \kappa \tau K_{12} \otimes K_{25} + \tau K_{12} \otimes K_{26})\mathbf{U} = f_{dof} F_1 \otimes F_2 - \tau f_{dof} F_1 \otimes F_3.
\end{equation}

\section{Numerical Verifications}\label{sec:numerical_results}
For each of the three problem \eirik{introduced in Section~\ref{sec:methods}, we consider and implement a specific test case} in FEniCS. \eirik{To verify the developed framework, we investigate the $h-$convergence properties of the implemented methods by consideration of the  rate of convergence of the FE solutions. The test and trial spaces in all cases  consist of continuous Lagrange polynomials of degree 1.}

\subsection{Case 1: 4-D Poisson Equation }

\eirik{ We first consider the Poisson PDE in the high dimensional space of order four. We define the computational domain as } 
a tensor product between two unit squares, \eirik{i.e.,}
$\Omega = ((0,1)\times(0,1)) \times ((0,1) \times (0,1))$.
\eirik{We select} the forcing function defined as $f = \text{4}\pi^2 u_{exact}$  and the exact solution as $u_{exact} = \Pi^4_{i=1} sin(\pi x_i) $,
 where $x_i$ is a coordinate in the domain $\Omega$.  

 \eirik{In Table~\ref{tab:results_poisson}, the convergence data for the four  dimensional Poisson problem is listed. Note that the convergence of the FE solution is optimal, as  the rate of convergence of the root mean square error (RMSE) and $l_{\infty}$ norm approaches $\mathcal{O}(h^{p+1})$}.
\begin{table}[h!]
\centering
\caption{\label{tab:results_poisson}  Error estimation results for the 4D Poisson problem.}
\begin{tabular}{@{}llllll@{}}
\toprule
{dofs \hspace{6mm}} & { h \hspace{6mm}} & {$l_{\infty}$ \hspace{6mm}} & {  RMSE  \hspace{6mm}} & {RMSE rate} & {$l_{\infty}$ rate} \\
\midrule \midrule

256 &  0.333    & 1.40E-01  & 3.16E-02  & - & - \\
1296 & 0.200    & 7.88E-02  & 1.61E-02  & 1.44 & 1.13 \\
2401 & 0.167    & 6.55E-02  & 1.21E-02  & 1.57 & 1.01  \\
4096 & 0.143    & 4.50E-02  & 9.39E-03  & 1.65 & 2.44 \\
6561 & 0.125    & 3.76E-02  & 7.48E-03  & 1.70 & 1.34\\
1461 & 0.100   & 2.42E-02  & 5.06E-03  & 1.75 & 1.96 \\
20736 & 0.091   & 1.95E-02  & 4.26E-03  & 1.79 & 2.31 \\
28561 & 0.083   & 1.69E-02  & 3.64E-03  & 1.81 & 1.61 \\
38416 & 0.077   & 1.41E-02  & 3.15E-03  & 1.83  & 2.26 \\
50625 & 0.071   & 1.25E-02  & 2.74E-03. & 1.84  & 1.67  \\



\bottomrule
\end{tabular}
\end{table}



\subsection{Case 2: 2D Space-Time Wave Equation}

\eirik{As a verification of the space-time wave model problem}~\eqref{eqn:model_prob2}, \eirik{we select} the space-time domain as the Cartesian product of a unit square \eirik{spatial} domain $\Omega$ and an interval time domain, i.e., $\Omega_T = ((0,1) \times (0,1)) \times (0,T)$.  The wave propagation speed is $c=1$ and \eirik{we consider a}  manufactured solution $u_{exact} = sin(x-ct) + sin(y-ct)$. \eirik{This solution is used to ascertain boundary and initial conditions needed to solve the resulting system of equations.}
\eirik{In Table~\ref{tab:wave_spacetime}, the convergence results are presented along with the time interval element size, denoted by dt, and the space-time CFL number. The RMSE is observed to converge at the expected optimal rate, whereas the  }  $l_{\infty}$ error \eirik{exhibits a reduced rate for the finer meshes}.
\begin{table}[h!]
\centering
\caption{\label{tab:wave_spacetime}  Results for the Wave Space-Time problem.}
\begin{tabular}{@{}llllllll@{}}
\toprule
{dofs \hspace{6mm}} & { h \hspace{6mm}} & { dt \hspace{6mm}} & { CFL \hspace{6mm}} & {$l_{\infty}$ \hspace{6mm}} & {  RMSE  \hspace{15mm}} & {RMSE rate} & {$l_{\infty}$ rate} \\
\midrule \midrule

200 &  0.250 & 0.14 & 0.57  & 2.25E-04  & 5.84E-05 & - & -   \\
1215 & 0.125 & 0.07 & 0.57  & 5.06E-05  & 1.51E-05 & 1.95 & 2.22  \\
3718 & 0.083 & 0.05 & 0.57  & 2.53E-05  & 7.01E-06 & 1.89 & 1.71 \\
8381 & 0.063 & 0.04 & 0.57  & 1.31E-05  & 3.95E-06 & 1.99 & 1.31 \\


\bottomrule
\end{tabular}
\end{table}


\subsection{Case 3: SUPG Stabilized Advection Dominated Advection Diffusion Equation}
\eirik{As a final numerical verification, we consider a special case of advection dominated advection diffusion equation. In particular, we consider the case in which the advection 
acts in a single direction aligned with a coordinate axis, see~\eqref{eqn:discrete_AD_1direction}.    }
 \eirik{We consider a case where the} domain $\Omega$ is a tensor product of 2 unit intervals: $\Omega = (0,1) \times (0,1)$, diffusivity constant $\kappa =\frac{1}{100}$, and \eirik{advection  vector is} $\mathbf{b} = (0,1)$. The analytic solution in this case \eirik{is inspired by the work of Egger and Sch{\"o}berl~\cite{egger2010hybrid}}:
\begin{equation} \label{eq:ADR_exact}
u_{exact} =  (-4(x-0.5)^2+1)\left[ y + \frac{e^{\frac{1}{\kappa} \cdot b_y \cdot y} - 1 }{1 - e^{\frac{1}{\kappa} \cdot b_y}}\right],
\end{equation}
which implies that the forcing term must be $f = -\kappa\,\Delta u_{exact} + \mathbf{b} \cdot \Nabla u_{exact}$ \eirik{and we enforce the corresponding} homogeneous Dirichlet on the boundary $\partial \Omega$.
\eirik{In Figures~\ref{fig:cartesian_ADR}, and~\ref{fig:analytic_ADR}, we show the approximate FE solution and the analytic solution, respectively. As expected, the SUPG stabilization results in a stable solution at this relatively coarse FE mesh.}
\eirik{In Table~\ref{tab:results_ADR}, the convergence results for this final case are presented. Both the RMSE and $l_{\infty}$ error converge at the expected optimal rates.} 
\begin{table}[h]
\centering
\caption{\label{tab:results_ADR}  Error estimation results for advection dominated advection diffusion problem.}
\begin{tabular}{@{}llllll@{}}
\toprule
{dofs \hspace{6mm}} & { h \hspace{6mm}} & {$l_{\infty}$ \hspace{6mm}} & {  RMSE  \hspace{15mm}} & {RMSE rate} & {$l_{\infty}$ rate} \\
\midrule \midrule

9    & 0.50 & 2.50E-01   & 8.33E-02 & -  & - \\
16   & 0.33 & 1.19E-01   & 4.65E-02 & 1.44 & 1.84 \\
25   & 0.25 & 7.38E-02   & 3.01E-02 & 1.51 & 1.65 \\
36   & 0.20 & 4.79E-02   & 2.08E-02 & 1.65 & 1.93  \\
49   & 0.17 & 3.33E-02   & 1.52E-02 & 1.73 & 1.99 \\
256  & 0.07 & 6.41E-03   & 2.85E-03 & 1.83 & 1.80 \\
676  & 0.04 & 2.52E-03   & 1.07E-03 & 1.92 & 1.83 \\
2601 & 0.02 & 6.50E-04   & 2.74E-04 & 1.96 & 1.96 \\





\bottomrule
\end{tabular}
\end{table}
\begin{figure}[h!]
    \centering
    \includegraphics[width=0.5\textwidth]{Analytic_256_dof.png}
    \caption{Analytic solution projected onto FE mesh with 256 dofs.}
    \label{fig:analytic_ADR}
\end{figure}
\begin{figure}[h!]
    \centering
    \includegraphics[width=0.5\textwidth]{Cartesian_Product_256_dof.png}
    \caption{Cartesian product FE solution at 256 dofs.}
    \label{fig:cartesian_ADR}
\end{figure}
\section{Conclusions}\label{sec:Conclusions}

\eirik{In this paper, we have introduced and implemented tensor product FE routines  for high dimensional problems in FEniCS. This methodology allows us to extend the}
FEniCS library to domains with more than three dimensions so long as they are a Cartesian product of subdomains three or lower.
\eirik{To verify the developed methodology, we consider three test cases utilizing classical and stabilized FE methods. } For each test case, \eirik{we observe the expected } convergence to the analytic solutions with respect to grid refinement was demonstrated in both the RMSE and $l_{\infty}$ norms.

\eirik{We consider only linear PDEs here }  since this allowed for the explicit construction of a single linear system of algebraic equations. \eirik{Hence, future studies should investigate potential extensions to nonlinear PDEs.} Additionally, the global system of equations was solved naively by explicitly constructing the global stiffness matrix and inverting. \eirik{However,} the global matrix is sparse with highly structured blocks which should allow for faster solvers that would greatly reduce run time. \eirik{We refer to existing works~\cite{BIALECKI1993369,Gao_thesis,Bank1978}, where  related}  problems  were considered \eirik{and leave the consideration of such solvers for future studies.} 
\eirik{Further extensions to mixed FE methods, } such as those discussed in the book by Brezzi \emph{et al.}~\cite{brezzi2012mixed}. \eirik{Finally, } FE methods utilizing with discontinuous test/trial spaces \eirik{ could be considered due to their extensive use in engineering applications}. Furthermore, full integration of this method into the FEniCS API would be valuable for both simplicity of future implementations as well as for performance.

 

The code used in the numerical studies is publicly available and can be found on GitHub at \\ \url{https://github.com/Markloveland/Cartesian_Product_FEM} 


\section{Acknowledgements}

Author Loveland has been supported by the CSEM Fellowship from the Oden Institute at the University of Texas at Austin. Authors Loveland, Valseth, and Dawson have been supported by the United States National Science Foundation - NSF PREEVENTS Track
2 Program, under NSF Grant Number 1855047 and the Department of Homeland Security Coastal Resilience Center research project "Accurate and Fast Wave Modeling and Coupling with ADCIRC".
Author Lukac has been supported by the University of Oregon.

 
\bibliographystyle{elsarticle-num}
 
2022-02-03T02:04:18.000
2202
2202.00762
en
https://arxiv.org/abs/2202.00762
We present a method to extend the finite element library FEniCS to solve problems with domains in dimensions above three by constructing tensor product finite elements. This methodology only requires that the high dimensional domain is structured as a Cartesian product of two lower dimensional subdomains. In this study we consider Dirichlet problems for scalar linear partial differential equations, though the methodology can be extended to non-linear problems. The utilization of tensor product finite elements allows us to construct a global system of linear algebraic equations that only relies on the finite element infrastructure of the lower dimensional subdomains contained in FEniCS. We demonstrate the effectiveness of our methodology in four distinctive test cases. The first test case is a Poisson equation posed in a four dimensional domain which is a Cartesian product of two unit squares solved using the classical Galerkin finite element method. The second test case is the wave equation in space-time, where the computational domain is a Cartesian product of a two dimensional space grid and a one dimensional time interval. In this second case we also employ the Galerkin method. The third test case is an advection dominated advection-diffusion equation where the global domain is a Cartesian product of two one dimensional intervals in which the streamline upwind Petrov-Galerkin method is applied to ensure discrete stability. The final test case uses the Galerkin approach to solve a Poisson problem on a Cartesian product of two intervals with a spatially varying, non-separable diffusivity term. In all cases, a p=1 basis is used and optimal L^2 convergence rates of order h^{p+1} of the errors are achieved with respect to h refinement
Numerical Analysis (math.NA); Mathematical Software (cs.MS)
Extending FEniCS to Work in Higher Dimensions Using Tensor Product Finite Elements
Qwen/Qwen-72B
1. YES
2. YES