4
$\begingroup$

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?

$\endgroup$
  • $\begingroup$ If your interpolation is bilinear you can solve it analytically. $\endgroup$ – nicoguaro Apr 21 '15 at 14:03
  • $\begingroup$ @nicoguaro: No you can't. It's still nonlinear. $\endgroup$ – Wolfgang Bangerth Apr 22 '15 at 12:00
  • $\begingroup$ @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. $\endgroup$ – nicoguaro Apr 22 '15 at 22:01
  • 1
    $\begingroup$ @nicoguaro: Fascinating! I was not aware that there exists a closed-form solution. Can you do the same in 3d? $\endgroup$ – Wolfgang Bangerth Apr 23 '15 at 11:33
  • $\begingroup$ @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). $\endgroup$ – nicoguaro Apr 23 '15 at 12:32
7
$\begingroup$

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

$\endgroup$
  • 1
    $\begingroup$ Thanks, this helped. When I was implementing for triangle elements, I hadn't even really thought about it and just computed the inverse mapping since it was easy. I'm adapting my code to quad elements now and fell into following the same formula that was already implemented, completely overlooking the easiest way to do this facepalm Thanks again $\endgroup$ – Justin Dong Apr 21 '15 at 14:46
  • $\begingroup$ Newton iteration is the best method that I have found. I describe such a technique for general hexahedral elements in section 3.2 of this paper. $\endgroup$ – Jeff Irwin Apr 22 '15 at 2:13
  • $\begingroup$ 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(). $\endgroup$ – Wolfgang Bangerth Apr 22 '15 at 2:22
  • $\begingroup$ @WolfgangBangerth I was unable to find those functions in the link you provided. Is there deal.ii documentation for how implement computation of the inverse transformation with a Newton Iteration? $\endgroup$ – user3482876 Apr 7 '18 at 16:50
  • $\begingroup$ @user3482876: These functions have moved to other files in the same directory since I wrote that comment. But if you search for them, you will still find them. $\endgroup$ – Wolfgang Bangerth Apr 8 '18 at 4:03

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.