# Inverse isoparametric mappings for quadrilateral finite elements

I have an isoparametric mapping $F_{E}: \hat{E} \to E$ where $\hat{E}$ is a reference quadrilateral (square) and $E$ is some quadrilateral in the domain $\Omega$. If $E$ is a parallelogram, then the mapping is affine and computing $F_{E}^{-1}$ is straightforward. But for the general case, we have: $$F_{E}(\hat{x},\hat{y}) = \left(\begin{array}{c}\sum_{i=1}^{4}x_{i}\hat{\psi}_{i}(\hat{x},\hat{y})\\ \sum_{i=1}^{4}y_{i}\hat{\psi}_{i}(\hat{x},\hat{y})\end{array}\right) = \left(\begin{array}{c}x\\ y\end{array}\right)$$

where $\hat{\psi}_{i}$ are the basis functions on $\hat{E}$.

For the discontinuous Galerkin method, I need to compute some integrals over the interfaces (2D) of the elements. A simple example would be the following: $$\int_{e}\langle\langle\nabla\psi_{i}\cdot\mathbf{n}\rangle\rangle[[\psi_{j}]]$$

where $\langle\langle\cdot\rangle\rangle$ and $[[\cdot]]$ denote the average and jump operators, respectively.

We have that $\psi_{j}(x,y) = \hat{\psi}_{j}(F_{E}^{-1}(x,y)) = \hat{\psi}_{j}(\hat{x},\hat{y})$. Furthermore, on a given edge $e$, the point $(x,y) \in e$ maps to two values due to the discontinuities of the basis functions across the interfaces, i.e. we have $\psi_{j}^{-}$ and $\psi_{j}^{+}$ along $e$.

To compute the above integral, I use a quadrature rule on the interval $\hat{e} = [-1,1]$. For a quadrature point $\hat{\zeta} \in \hat{e}$, I map this 1D point to $e$ and call it $\zeta$.

My goal is to compute $\hat{\psi}_{j}^{-}(F_{E-}^{-1}(\zeta))$ and $\hat{\psi}_{j}^{+}(F_{E+}^{-1}(\zeta))$, where $F_{E-}$ and $F_{E+}$ are the mappings associated with the neighboring elements to $e$. For triangular elements and parallelograms, this is really straightforward since the mapping is affine and thus, computing $F_{E}^{-1}$ is easy.

But for general quadrilaterals, is this even possible? I found this paper, but it's quite old and the process described seemed a bit complicated. Is there any other literature anyone can recommend on this? Or some insight into a simpler way to do this?

• If your interpolation is bilinear you can solve it analytically. – nicoguaro Apr 21 '15 at 14:03
• @nicoguaro: No you can't. It's still nonlinear. – Wolfgang Bangerth Apr 22 '15 at 12:00
• @WolfgangBangerth, why nonlinear is synonym of unsolvable? I plug the equations into SymPy and got a solution here. If you do it by hand it is prettier (of course), but I don't have much time now. – nicoguaro Apr 22 '15 at 22:01
• @nicoguaro: Fascinating! I was not aware that there exists a closed-form solution. Can you do the same in 3d? – Wolfgang Bangerth Apr 23 '15 at 11:33
• @Wolfgang, I don't think. I haven't tried though, but it's not as interesting. If available the solution for cubics is messy (more). – nicoguaro Apr 23 '15 at 12:32

You rarely compute the inverse transform (though it's not overly difficult: you just need to do a Newton iteration). The thing is that to compute $\hat\Psi_j^\pm(F_{E\pm}^{-1}(\zeta_k))$ at a quadrature point $\zeta_k$ you already have everything you need: the quadrature points $\zeta_k$ are in fact the mapped locations $\zeta_k=F_{E\pm}(\hat\zeta_k)$ of the quadrature points $\hat\zeta_k$ defined on the reference cell!
So, $\hat\Psi_j^\pm(F_{E\pm}^{-1}(\zeta_k)) = \hat\Psi_j^\pm(\hat\zeta_k)$. In other words, the evaluation of shape functions happens on the reference cell, not the real cell. It's really that easy.
(In the interest of expanding the discussion: The only times when it gets more interesting is if someone wants to know what the value of the finite elements solution is at an arbitrary point $\zeta$ that is not the mapped location of a quadrature point. In that case, you really have to find the cell that $\zeta$ is in, and then do the inverse mapping $F_{E\pm}^{-1}(\zeta)$. This is expensive, and such queries are consequently avoided if possible.)
• Yes, all big finite element packages implement Newton's method in one way or another. Our implementation is here: github.com/dealii/dealii/blob/master/source/fe/mapping_q1.cc Search for functions transform_real_to_unit_cell() and transform_real_to_unit_cell_internal(). – Wolfgang Bangerth Apr 22 '15 at 2:22