Transfinite mapping

By performing a so-called transfinite mapping, all operations on basis functions (p.a. multiplication, differentiation and integration) within a single mesh element (defined by physical coordinates) can be performed on a standard element (defined by natural coordinates, which are within the range -1 and +1). As a result, all these operations become independent of the specific mesh element. Therefore, this approach is more efficient and easier to implement correctly.

Interpolation within a standard element

A standard element is defined as

(1)\[\Omega_{st} = \left\{ \left( \xi, \eta, \zeta \right) | -1 \leq \xi,\eta,\zeta \leq 1 \right\},\]

where the number of natural coordinates (\(\xi\), \(\eta\), \(\zeta\)) matches the dimension of the mesh. In the figure below, the 2D transfinite mapping from the physical to the natural domain is visualized.

Figure made with TikZ

Fig. 1 2D transfinite mapping

Each natural location within the standard element corresponds to a (unique) physical location in the mesh element. The transfinite mapping is defined as the mapping from the physical coordinate to the natural coordinate.

Consider a linear edge \(\mathbf{e}_1\) to be defined by two vertices \(\mathbf{v}_1\) and \(\mathbf{v}_2\). This edge can be parametrized by a natural coordinate \(\xi\) (which has a range of -1 to +1) as

(2)\[\begin{split}\begin{split} \mathbf{e}_1\left( \xi \right) &= \frac{\mathbf{v}_1+\mathbf{v}_2}{2} + \xi \frac{\mathbf{v}_2-\mathbf{v}_1}{2},\\ &= \mathbf{v}_1 \frac{1-\xi}{2} + \mathbf{v}_2 \frac{1+\xi}{2}. \end{split}\end{split}\]

Similar relations for the other edges can be determined. Next, consider a flat face \(\mathbf{f}\) to be defined by four edges \(\mathbf{e}_1\), \(\mathbf{e}_2\), \(\mathbf{e}_3\), and \(\mathbf{e}_4\). This face can be parametrized by the natural coordinate \(\xi\) and \(\eta\) (which both have ranges of -1 to +1) as

(3)\[\mathbf{f}_5\left( \xi, \eta \right) = \mathbf{e}_1 \frac{1-\eta}{2} + \mathbf{e}_2 \frac{1+\xi}{2} + \mathbf{e}_2 \frac{1+\eta}{2} + \mathbf{e}_4 \frac{1-\xi}{2} - \left\{ \mathbf{v}_1 \frac{1-\xi}{2} \frac{1-\eta}{2} + \mathbf{v}_2 \frac{1+\xi}{2} \frac{1-\eta}{2} + \mathbf{v}_3 \frac{1+\xi}{2} \frac{1+\eta}{2} + \mathbf{v}_4 \frac{1-\xi}{2} \frac{1+\eta}{2} \right\}.\]

In this last relation, the need for the subtraction of the vertices can be understood after realizing that these are added twice by the edge contributions. Again, similar relations for the other faces can be determined.

Finally, a hexahedron (volume) is defined by 6 faces, and its parametrization can be written as

(4)\[\begin{split}\begin{split} \mathbf{h}\left( \xi, \eta, \zeta \right) = &\mathbf{f}_1 \frac{1-\xi}{2} + \mathbf{f}_2 \frac{1+\xi}{2} + \mathbf{f}_3 \frac{1-\eta}{2} + \mathbf{f}_4 \frac{1+\eta}{2} + \mathbf{f}_5 \frac{1-\zeta}{2} + \mathbf{f}_6 \frac{1+\zeta}{2} \\ &- \mathbf{e}_1 \frac{1-\eta}{2} \frac{1-\zeta}{2} - \mathbf{e}_3 \frac{1+\eta}{2} \frac{1-\zeta}{2} - \mathbf{e}_5 \frac{1-\eta}{2} \frac{1+\zeta}{2} - \mathbf{e}_7 \frac{1+\eta}{2} \frac{1+\zeta}{2} \\ &- \mathbf{e}_2 \frac{1+\xi}{2} \frac{1-\zeta}{2} - \mathbf{e}_4 \frac{1-\xi}{2} \frac{1-\zeta}{2} - \mathbf{e}_6 \frac{1+\xi}{2} \frac{1+\zeta}{2} - \mathbf{e}_8 \frac{1-\xi}{2} \frac{1+\zeta}{2} \\ &- \mathbf{e}_9 \frac{1-\xi}{2} \frac{1-\eta}{2} - \mathbf{e}_{10} \frac{1+\xi}{2} \frac{1-\eta}{2} - \mathbf{e}_{11} \frac{1+\xi}{2} \frac{1+\eta}{2} - \mathbf{e}_{12} \frac{1-\xi}{2} \frac{1+\eta}{2}\\ &+ \left\{ \left( \mathbf{v}_1 \frac{1-\xi}{2} \frac{1-\eta}{2} + \mathbf{v}_2 \frac{1+\xi}{2} \frac{1-\eta}{2} + \mathbf{v}_3 \frac{1+\xi}{2} \frac{1+\eta}{2} + \mathbf{v}_4 \frac{1-\xi}{2} \frac{1+\eta}{2} \right) \frac{1-\zeta}{2} + \left( \mathbf{v}_5 \frac{1-\xi}{2} \frac{1-\eta}{2} + \mathbf{v}_6 \frac{1+\xi}{2} \frac{1-\eta}{2} + \mathbf{v}_7 \frac{1+\xi}{2} \frac{1+\eta}{2} + \mathbf{v}_8 \frac{1-\xi}{2} \frac{1+\eta}{2} \right) \frac{1+\zeta}{2} \right\}, \end{split}\end{split}\]

where the first part is the contribution of the faces, the second part a correction to prevent duplication of edges, and the third part to ensure a correct contribution of vertices (the first and second part cancel out the vertices exactly).

Note that the entire parametrization depends on the definition of position/direction of vertices, edges and faces. The following definitions are assumed here:

Vertex id

\(\xi\) coordinate

\(\eta\) coordinate

\(\zeta\) coordinate

1

-1

-1

-1

2

+1

-1

-1

3

+1

+1

-1

4

-1

+1

-1

5

-1

-1

+1

6

+1

-1

+1

7

+1

+1

+1

8

-1

+1

+1

Edge id

vertex1 id

vertex2 id

direction/variable

constant values

1

1

2

\(+\xi\)

\(\eta\equiv-1\), \(\zeta\equiv-1\)

2

2

3

\(+\eta\)

\(\xi \equiv+1\), \(\zeta\equiv-1\)

3

4

3

\(+\xi\)

\(\eta\equiv+1\), \(\zeta\equiv-1\)

4

1

4

\(+\eta\)

\(\xi \equiv-1\), \(\zeta\equiv-1\)

5

5

6

\(+\xi\)

\(\eta\equiv-1\), \(\zeta\equiv+1\)

6

6

7

\(+\eta\)

\(\xi \equiv+1\), \(\zeta\equiv+1\)

7

8

7

\(+\xi\)

\(\eta\equiv+1\), \(\zeta\equiv+1\)

8

5

8

\(+\eta\)

\(\xi \equiv-1\), \(\zeta\equiv+1\)

9

1

5

\(+\zeta\)

\(\xi \equiv-1\), \(\eta\equiv-1\)

10

2

6

\(+\zeta\)

\(\xi \equiv+1\), \(\eta\equiv-1\)

11

3

7

\(+\zeta\)

\(\xi \equiv+1\), \(\eta\equiv+1\)

12

4

8

\(+\zeta\)

\(\xi \equiv-1\), \(\eta\equiv+1\)

Face id

edge1 id

edge2 id

edge3 id

edge4 id

normal direction

constant value

1

4

12

-8

-9

\(+\xi\)

\(\xi\equiv-1\)

2

2

11

-6

-10

\(+\xi\)

\(\xi\equiv+1\)

3

1

10

-5

-9

\(+\eta\)

\(\eta\equiv-1\)

4

3

11

-7

-12

\(+\eta\)

\(\eta\equiv+1\)

5

1

2

-3

-4

\(+\zeta\)

\(\zeta\equiv-1\)

6

5

6

-7

-8

\(+\zeta\)

\(\zeta\equiv+1\)

Note here that a minus sign indicates an opposite edge direction.

Differentiation within a standard element

In partial differential equations (PDEs) the differentials are defined in the physical domain. In order to use the differentiation within a standard element, a (back) transformation from the natural to the physical domain needs to be made.

1D element

For a 1D element, the first derivative is given by the following chain rule

(5)\[\frac{\partial}{\partial x} = \frac{\partial \xi}{\partial x} \frac{\partial}{\partial \xi},\]

where the partial derivative \(\frac{\partial \xi}{\partial x}\) is just the inverse of \(\frac{\partial x}{\partial \xi}\). The second derivative is given by

(6)\[\begin{split}\begin{split} \frac{\partial^2}{\partial x^2} &= \frac{\partial \xi}{\partial x} \frac{\partial}{\partial \xi} \left( \frac{\partial \xi}{\partial x} \frac{\partial}{\partial \xi} \right),\\ &= \left( \frac{\partial \xi}{\partial x} \right)^2 \frac{\partial^2}{\partial \xi^2} + \frac{\partial^2 \xi}{\partial x^2} \frac{\partial}{\partial \xi} \end{split}\end{split}\]

where the partial derivatives \(\frac{\partial \xi}{\partial x}\) and \(\frac{\partial^2 \xi}{\partial x^2}\) are just the inverses of \(\frac{\partial x}{\partial \xi}\) and \(\frac{\partial^2 x}{\partial \xi^2}\), respectively.

2D element

For a 2D element the chain rule gives

(7)\[\begin{split}\nabla = \begin{bmatrix} \frac{\partial}{\partial x} \\ \frac{\partial}{\partial y} \end{bmatrix} = \begin{bmatrix} \frac{\partial \xi}{\partial x} \frac{\partial}{\partial \xi} + \frac{\partial \eta}{\partial x} \frac{\partial}{\partial \eta} \\ \frac{\partial \xi}{\partial y} \frac{\partial}{\partial \xi} + \frac{\partial \eta}{\partial y} \frac{\partial}{\partial \eta} \end{bmatrix},\end{split}\]

where \(\nabla\) is the nable operator. The partial derivative with respect to \(x\) and \(y\) can be obtained analytically for a linear mapping, but a more general technique is required for curvilinear mappings. The partial derivatives with respect to \(\xi\) and \(\eta\) can be obtained directly from the transfinite mapping. The following chain rule can be used for the total change of a general function \(u(\xi,\eta)\)

(8)\[\textrm{d}u(\xi,\eta) = \frac{\partial u}{\partial \xi} \textrm{d}\xi + \frac{\partial u}{\partial \eta} \textrm{d}\eta.\]

If the function \(u\) is replaced by the mapping \(\mathcal{X}(\xi,\eta)\), the following system is obtained

(9)\[\begin{split} \begin{bmatrix} \textrm{d}x\\ \textrm{d}y \end{bmatrix} = \begin{bmatrix} \frac{\partial x}{\partial \xi} & \frac{\partial x}{\partial \eta} \\ \frac{\partial y}{\partial \xi} & \frac{\partial y}{\partial \eta} \end{bmatrix} \begin{bmatrix} \textrm{d}\xi\\ \textrm{d}\eta \end{bmatrix}.\end{split}\]

This can be inverted to obtain

(10)\[\begin{split} \begin{bmatrix} \textrm{d}\xi\\ \textrm{d}\eta \end{bmatrix} = \frac{1}{\left| J_{\textrm{2D}} \right|} \begin{bmatrix} \frac{\partial y}{\partial \eta} & -\frac{\partial x}{\partial \eta} \\ -\frac{\partial y}{\partial \xi} & \frac{\partial x}{\partial \xi} \end{bmatrix} \begin{bmatrix} \textrm{d}x\\ \textrm{d}y \end{bmatrix},\end{split}\]

where \(\left| J_{\textrm{2D}} \right|\) is the 2D Jacobian defined by

(11)\[\begin{split}\left| J_{\textrm{2D}} \right| = \begin{vmatrix} \frac{\partial x}{\partial \xi} & \frac{\partial x}{\partial \eta} \\ \frac{\partial y}{\partial \xi} & \frac{\partial y}{\partial \eta} \end{vmatrix} = \frac{\partial x}{\partial \xi}\frac{\partial y}{\partial \eta} - \frac{\partial x}{\partial \eta}\frac{\partial y}{\partial \xi}.\end{split}\]

As the transfinite mapping is assumed to be one-to-one and have an inverse, the chain rule can also be applied directly to \(\xi\) and \(\eta\) to obtain

(12)\[\begin{split} \begin{bmatrix} \textrm{d}\xi\\ \textrm{d}\eta \end{bmatrix} = \begin{bmatrix} \frac{\partial \xi}{\partial x} & \frac{\partial \xi}{\partial y} \\ \frac{\partial \eta}{\partial x} & \frac{\partial \eta}{\partial y} \end{bmatrix} \begin{bmatrix} \textrm{d}x\\ \textrm{d}y \end{bmatrix}.\end{split}\]

Combining these results leads to the following partial derivatives

(13)\[ \frac{\partial \xi}{\partial x} = \frac{1}{\left| J_{\textrm{2D}} \right|} \frac{\partial y}{\partial \eta}, \qquad \frac{\partial \xi}{\partial y} = -\frac{1}{\left| J_{\textrm{2D}} \right|} \frac{\partial x}{\partial \eta}, \qquad \frac{\partial \eta}{\partial x} = -\frac{1}{\left| J_{\textrm{2D}} \right|} \frac{\partial y}{\partial \xi}, \qquad \frac{\partial \eta}{\partial y} = \frac{1}{\left| J_{\textrm{2D}} \right|} \frac{\partial x}{\partial \xi}.\]
(14)\[\begin{split}\nabla = \begin{bmatrix} \frac{\partial}{\partial x} \\ \frac{\partial}{\partial y} \end{bmatrix} = \frac{1}{\left| J_{\textrm{2D}} \right|} \underbrace{ \begin{bmatrix} \frac{\partial y}{\partial \eta} & -\frac{\partial x}{\partial \eta}\\ -\frac{\partial y}{\partial \xi} & \frac{\partial x}{\partial \xi} \end{bmatrix} }_{M} \begin{bmatrix} D_{\xi}\\ D_{\eta} \end{bmatrix},\end{split}\]

where \(D_{\xi}\) and \(D_{\eta}\) are differential operators of the basis functions. Note that the array \(M\) becomes the identity matrix if \(x\) and \(y\) are aligned with \(\xi\) and \(\eta\), respectively.

A similar derivation for the second derivatives can be made, which leads to

(15)\[\begin{split}\begin{split} \frac{\partial^2}{\partial \xi^2} &= \left( \frac{\partial x}{\partial \xi} \right)^2 \frac{\partial^2}{\partial x^2} + \left( \frac{\partial y}{\partial \xi} \right)^2 \frac{\partial^2}{\partial y^2} + 2 \frac{\partial x}{\partial \xi} \frac{\partial y}{\partial \xi} \frac{\partial^2}{\partial x \partial y} + \frac{\partial^2 x}{\partial \xi^2} \frac{\partial}{\partial x} + \frac{\partial^2 y}{\partial \xi^2} \frac{\partial}{\partial y},\\ \frac{\partial^2}{\partial \eta^2} &= \left( \frac{\partial x}{\partial \eta} \right)^2 \frac{\partial^2}{\partial x^2} + \left( \frac{\partial y}{\partial \eta} \right)^2 \frac{\partial^2}{\partial y^2} + 2 \frac{\partial x}{\partial \eta} \frac{\partial y}{\partial \eta} \frac{\partial^2}{\partial x \partial y} + \frac{\partial^2 x}{\partial \eta^2} \frac{\partial}{\partial x} + \frac{\partial^2 y}{\partial \eta^2} \frac{\partial}{\partial y},\\ \frac{\partial^2}{\partial \xi \partial \eta} &= \frac{\partial x}{\partial \xi} \frac{\partial x}{\partial \eta} \frac{\partial^2}{\partial x^2} + \frac{\partial y}{\partial \xi} \frac{\partial y}{\partial \eta} \frac{\partial^2}{\partial y^2} + \left( \frac{\partial x}{\partial \eta} \frac{\partial y}{\partial \xi} + \frac{\partial y}{\partial \eta} \frac{\partial x}{\partial \xi} \right) \frac{\partial^2}{\partial x \partial y} + \frac{\partial^2 x}{\partial \xi \partial \eta} \frac{\partial}{\partial x} + \frac{\partial^2 y}{\partial \xi \partial \eta} \frac{\partial}{\partial y}, \end{split}\end{split}\]

All derivatives of the physical with respect to the natural space can be derived directly from the transfinite mapping. However, as the PDEs require derivative operators with respect to physical coordinates, for practical use this needs to be rewritten as

(16)\[\begin{split}\begin{Bmatrix} \frac{\partial^2}{\partial^2 x} \\ \frac{\partial^2}{\partial^2 y} \\ \frac{\partial^2}{\partial x \partial y} \end{Bmatrix} = \underbrace{ \begin{bmatrix} \left( \frac{\partial x}{\partial \xi} \right)^2 & \left( \frac{\partial y}{\partial \xi} \right)^2 & 2 \frac{\partial x}{\partial \xi} \frac{\partial y}{\partial \xi}\\ \left( \frac{\partial x}{\partial \eta} \right)^2 & \left( \frac{\partial y}{\partial \eta} \right)^2 & 2 \frac{\partial x}{\partial \eta} \frac{\partial y}{\partial \eta}\\ \frac{\partial x}{\partial \xi} \frac{\partial x}{\partial \eta} & \frac{\partial y}{\partial \xi} \frac{\partial y}{\partial \eta} & \frac{\partial x}{\partial \eta} \frac{\partial y}{\partial \xi} + \frac{\partial y}{\partial \eta} \frac{\partial x}{\partial \xi} \end{bmatrix}^{-1}}_{\hat{J_2}} \left( \begin{Bmatrix} \frac{\partial^2}{\partial^2 \xi} \\ \frac{\partial^2}{\partial^2 \eta} \\ \frac{\partial^2}{\partial \xi \partial \eta} \end{Bmatrix} - \underbrace{ \begin{bmatrix} \frac{\partial^2 x}{\partial \xi^2} & \frac{\partial^2 y}{\partial \xi^2}\\ \frac{\partial^2 x}{\partial \eta^2} & \frac{\partial^2 y}{\partial \eta^2}\\ \frac{\partial^2 x}{\partial \xi \partial \eta} & \frac{\partial^2 y}{\partial \xi \partial \eta} \end{bmatrix}}_{\bar{J_2}} \begin{Bmatrix} \frac{\partial}{\partial x}\\ \frac{\partial}{\partial y} \end{Bmatrix} \right).\end{split}\]

This implementation is used in the FEHandler class in Nemesis, see the listing below:

Listing 1 Simplified implementation of the Dxx operator (see the getDxx method in the FEHandler.m class)
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
 function out = getDxx( self, i_element, use_scaling )
   assert( nargin>1, 'Please provide an element number, or use getDxx_array' )
   if ( nargin == 2 )
     use_scaling = true;
   end

   finite_element = self.getFiniteElement( i_element );                                                    % Get the finite element of element i_element
   J_inverse_10 = self.quadrature.Jacobian{ i_element }{ 2 }.inverse_xx_10;                                % \hat{J_2}(1,1)
   J_inverse_01 = self.quadrature.Jacobian{ i_element }{ 2 }.inverse_xx_01;                                % \hat{J_2}(1,2)
   J_inverse_11 = self.quadrature.Jacobian{ i_element }{ 2 }.inverse_xx_11;                                % \hat{J_2}(1,3)

   J_2nd_11 = self.quadrature.Jacobian{ i_element }{ 2 }.derivative_11;                                    % \bar{J_2}(1,1)
   J_2nd_12 = self.quadrature.Jacobian{ i_element }{ 2 }.derivative_12;                                    % \bar{J_2}(1,2)
   J_2nd_21 = self.quadrature.Jacobian{ i_element }{ 2 }.derivative_21;                                    % \bar{J_2}(2,1)
   J_2nd_22 = self.quadrature.Jacobian{ i_element }{ 2 }.derivative_22;                                    % \bar{J_2}(1,2)
   J_2nd_31 = self.quadrature.Jacobian{ i_element }{ 2 }.derivative_31;                                    % \bar{J_2}(3,1)
   J_2nd_32 = self.quadrature.Jacobian{ i_element }{ 2 }.derivative_32;                                    % \bar{J_2}(3,2)

   % NB: Do not use dofs scaling, as this is done in the end of this function
   Dx = self.getDx( i_element, false );                                                                    % Dx operator (unscaled!)
   Dy = self.getDy( i_element, false );                                                                    % Dy operator (unscaled!)

   out = finite_element.zero_matrix;                                                                       % Initialize as a zero matrix
   out = out + ( finite_element.DD{1,1} - ( Dx .* J_2nd_11(:) + Dy .* J_2nd_12(:) ) ) .* J_inverse_10(:);  % DD{1,1} = \frac{\partial^2}{\partial^2 \xi}
   out = out + ( finite_element.DD{2,2} - ( Dx .* J_2nd_21(:) + Dy .* J_2nd_22(:) ) ) .* J_inverse_01(:);  % DD{2,2} = \frac{\partial^2}{\partial^2 \eta}
   out = out + ( finite_element.DD{1,2} - ( Dx .* J_2nd_31(:) + Dy .* J_2nd_32(:) ) ) .* J_inverse_11(:);  % DD{1,2} = \frac{\partial^2}{\partial \xi \partial \eta}

   if ( use_scaling )
     out = out * self.dofs.scaling.elemental{ i_element };                                                 % Scaling of dofs, which necessary for AMR
   end
 end