Inverse isoparametric mappings for quadrilateral finite elements

Question

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?

Answer 1

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.)