NLMech  0.1.0
Public Member Functions | Private Member Functions
fe::TetElem Class Reference

A class for mapping and quadrature related operations for linear tetrahedron element. More...

#include <tetElem.h>

Inheritance diagram for fe::TetElem:
Inheritance graph
Collaboration diagram for fe::TetElem:
Collaboration graph

Public Member Functions

double elemSize (const std::vector< util::Point3 > &nodes) override
 Returns the volume of element. More...
 
std::vector< std::vector< double > > getDerShapes (const util::Point3 &p, const std::vector< util::Point3 > &nodes) override
 Returns the values of derivative of shape function at point p. More...
 
std::vector< fe::QuadData > getQuadDatas (const std::vector< util::Point3 > &nodes) override
 Get vector of quadrature data. More...
 
std::vector< fe::QuadData > getQuadPoints (const std::vector< util::Point3 > &nodes) override
 Get vector of quadrature data. More...
 
std::vector< double > getShapes (const util::Point3 &p, const std::vector< util::Point3 > &nodes) override
 Returns the values of shape function at point p. More...
 
void print (int nt=0, int lvl=0) const
 Prints the information about the instance of the object. More...
 
std::string printStr (int nt=0, int lvl=0) const
 Returns the string containing information about the instance of the object. More...
 
 TetElem (size_t order)
 Constructor. More...
 
- Public Member Functions inherited from fe::BaseElem
 BaseElem (size_t order, size_t element_type)
 Constructor. More...
 
size_t getElemType ()
 Get element type. More...
 
size_t getNumQuadPoints ()
 Get number of quadrature points in the data. More...
 
size_t getQuadOrder ()
 Get order of quadrature approximation. More...
 
void print (int nt=0, int lvl=0) const
 Prints the information about the instance of the object. More...
 
std::string printStr (int nt=0, int lvl=0) const
 Returns the string containing information about the instance of the object. More...
 

Private Member Functions

std::vector< std::vector< double > > getDerShapes (const util::Point3 &p) override
 Returns the values of derivative of shape function at point p on reference element. More...
 
double getJacobian (const util::Point3 &p, const std::vector< util::Point3 > &nodes, std::vector< std::vector< double >> *J) override
 Computes the Jacobian of map $ \Phi: T^0 \to T $. More...
 
std::vector< double > getShapes (const util::Point3 &p) override
 Returns the values of shape function at point p on reference element. More...
 
void init () override
 Compute the quadrature points for triangle element.
 
util::Point3 mapPointToRefElem (const util::Point3 &p, const std::vector< util::Point3 > &nodes) override
 Maps point p in a given element to the reference element. More...
 

Detailed Description

A class for mapping and quadrature related operations for linear tetrahedron element.

The reference tetrahedron element $ T^0 $ is given by vertices $ (0,0,0), \, (1,0,0), \, (0,1,0), \, (0,0,1) $.

  1. The shape functions at point $ (\xi, \eta, \zeta) \in T^0 $ are

    \[N^0_1(\xi, \eta, \zeta) = 1- \xi - \eta - \zeta, \quad N^0_2(\xi, \eta, \zeta) = \xi, \quad N^0_3(\xi, \eta, \zeta) = \eta, \quad N^0_4(\xi, \eta, \zeta) = \zeta. \]

  2. For linear tetrahedron element, derivative of shape functions are constant and are as follows

    \[\frac{d N^0_1(\xi, \eta, \zeta)}{d\xi} = -1, \, \frac{d N^0_1(\xi, \eta, \zeta) }{d\eta} = -1, \, \frac{d N^0_1(\xi, \eta, \zeta) }{d\zeta} = -1, \]

    \[\frac{d N^0_2(\xi, \eta, \zeta)}{d\xi} = 1, \, \frac{d N^0_2(\xi, \eta, \zeta)}{d\eta} = 0, \, \frac{d N^0_2(\xi, \eta, \zeta)}{d\zeta} = 0, \]

    \[\frac{d N^0_2(\xi, \eta, \zeta)}{d\xi} = 0, \, \frac{d N^0_2(\xi, \eta, \zeta)}{d\eta} = 1, \, \frac{d N^0_2(\xi, \eta, \zeta)}{d\zeta} = 0, \]

    \[\frac{d N^0_2(\xi, \eta, \zeta)}{d\xi} = 0, \, \frac{d N^0_2(\xi, \eta, \zeta)}{d\eta} = 0, \, \frac{d N^0_2(\xi, \eta, \zeta)}{d\zeta} = 1. \]

  3. Map $ \Phi: T^0 \to T$ is given by

    \[ x(\xi, \eta, \zeta) = \sum_{i=1}^4 N^0_i(\xi, \eta, \zeta) v^i_x, \quad y (\xi, \eta, \zeta) = \sum_{i=1}^4 N^0_i(\xi, \eta, \zeta) v^i_y, \quad z (\xi, \eta, \zeta) = \sum_{i=1}^4 N^0_i(\xi, \eta, \zeta) v^i_z \]

    where $ v^1, v^2, v^3, v^4$ are vertices of element $ T $.
  4. The Jacobian of the map $ \Phi: T^0 \to T$ is given by

    \[ J = \left[ { \begin{array}{ccc} \frac{dx}{d\xi} &\frac{dy}{d\xi} &\frac{dz}{d\xi} \\ \frac{dx}{d\eta} & \frac{dy}{d\eta} &\frac{dz}{d\eta} \\ \frac{dx}{d\zeta} & \frac{dy}{d\zeta} &\frac{dz}{d\zeta} \\ \end{array} } \right] \]

    and determinant of Jacobian is

    \[ det(J) = \frac{dx}{d\xi} (\frac{dy}{d\eta} \times \frac{dz}{d\zeta} - \frac{dy}{d\zeta}\times \frac{dz}{d\eta}) - \frac{dy}{d\xi} (\frac{dx}{d\eta} \times \frac{dz}{d\zeta} - \frac{dx}{d\zeta}\times \frac{dz}{d\eta}) + \frac{dz}{d\xi} (\frac{dx}{d\eta} \times \frac{dy}{d\zeta} - \frac{dx}{d\zeta}\times \frac{dy}{d\eta}). \]

    For linear triangle element, Jacobian (and so $ det(J) $) is constant. For linear tetrahedron elements, we simply have

    \[ \frac{dx}{d\xi} = v^2_x - v^1_x, \quad \frac{dx}{d\eta} = v^3_x - v^1_x, \quad \frac{dx}{d\zeta} = v^4_x - v^1_x, \]

    \[ \frac{dy}{d\xi} = v^2_y - v^1_y, \quad \frac{dy}{d\eta} = v^3_y - v^1_y, \quad \frac{dy}{d\zeta} = v^4_y - v^1_y, \]

    \[ \frac{dz}{d\xi} = v^2_z - v^1_z, \quad \frac{dz}{d\eta} = v^3_z - v^1_z, \quad \frac{dz}{d\zeta} = v^4_z - v^1_z \]

    and $ det(J) = \frac{volume(T)}{volume(T^0)} = 6\times volume(T) $.
  5. Inverse map $ \Phi^{-1} : T \to T^0 $ for linear tetrahedron element can be easily derived. The derivation of the map is provided below:

From map $ (\xi, \eta, \zeta)\in T^0 \to (x,y,z) \in T $ we have

\[ x = \sum_{i=1}^4 N^0_i(\xi, \eta, \zeta) v^i_x, \quad y = \sum_{i=1}^4 N^0_i(\xi, \eta, \zeta) v^i_y, \quad z = \sum_{i=1}^4 N^0_i(\xi, \eta, \zeta) v^i_z. \]

Substituting formula for $ N^0_i $ in above to get

\[ x = (1 - \xi - \eta - \zeta) v^1_x + \xi v^2_x + \eta v^3_x + \zeta v^4_x, \quad y = (1 - \xi - \eta - \zeta) v^1_y + \xi v^2_y + \eta v^3_y + \zeta v^4_y, \quad z = (1 - \xi - \eta - \zeta) v^1_z + \xi v^2_z + \eta v^3_z + \zeta v^4_z \]

or

\[ x - v^1_x = \xi (v^2_x - v^1_x) + \eta (v^3_x - v^1_x) + \zeta (v^4_x - v^1_x), \quad y - v^1_y = \xi (v^2_y - v^1_y) + \eta (v^3_y - v^1_y) + \zeta (v^4_y - v^1_y), \quad z - v^1_z = \xi (v^2_z - v^1_z) + \eta (v^3_z - v^1_z) + \zeta (v^4_z - v^1_z). \]

Writing above in matrix form, we have

\[ \left[ {\begin{array}{c} x - v^1_x \\ y - v^1_y \\ z - v^1_z \end{array}}\right] = \left[ {\begin{array}{ccc} v^2_x - v^1_x & v^3_x - v^1_x & v^4_x - v^1_x\\ v^2_y - v^1_y & v^3_y - v^1_y & v^3_y - v^1_y \\ v^2_z - v^1_z & v^3_z - v^1_z & v^3_z - v^1_z \end{array}}\right] \, \left[ {\begin{array}{c} \xi \\ \eta \\ \zeta \end{array}}\right]. \]

Denoting the matrix as $ B $

\[ B = \left[ {\begin{array}{ccc} v^2_x - v^1_x & v^3_x - v^1_x & v^4_x - v^1_x\\ v^2_y - v^1_y & v^3_y - v^1_y & v^4_y - v^1_y \\ v^2_z - v^1_z & v^3_z - v^1_z & v^4_z - v^1_z \end{array}}\right] \]

. Note that $ B $ is transpose of Jacobin of map $\Phi$ therefore $ det(B) = det(J)$. Inverse of $ B$ is

\[ C := B^{-1} = \frac{1}{det(B)} \left[ {\begin{array}{ccc} B_{22}B_{33} - B_{32}B_{23} & B_{13}B_{32} - B_{33}B_{12} & B_{12}B_{23} - B_{22}B_{13} \\ B_{23}B_{31} - B_{33}B_{21} & B_{11}B_{33} - B_{31}B_{13} & B_{13}B_{21} - B_{23}B_{11} \\ B_{21}B_{32} - B_{31}B_{22} & B_{12}B_{31} - B_{32}B_{11} & B_{11}B_{22} - B_{21}B_{12} \\ \end{array}}\right]. \]

With $ C $ we have inverse map given by

\[ \xi(x,y,z) = C_{11} (x - v^1_x) + C_{12} (y - v^1_y) + C_{13} (z - v^1_z), \quad \eta(x,y,z) = C_{21} (x - v^1_x) + C_{22} (y - v^1_y) + C_{23} (z - v^1_z), \quad \zeta (x,y,z) = C_{31} (x - v^1_x) + C_{32} (y - v^1_y) + C_{33} (z - v^1_z). \]

Constructor & Destructor Documentation

◆ TetElem()

fe::TetElem::TetElem ( size_t  order)
explicit

Constructor.

Parameters
orderOrder of quadrature point approximation

Member Function Documentation

◆ elemSize()

double fe::TetElem::elemSize ( const std::vector< util::Point3 > &  nodes)
overridevirtual

Returns the volume of element.

If tetrahedron $ T $ is given by points $ v^1, v^2, v^3, v^4 $ then the volume is

\[ volume(T) =\frac{1}{3!} \left\vert \left[ {\begin{array}{cccc} v^1_x & v^1_y & v^1_z & 1 \\ v^2_x & v^2_y & v^2_z & 1 \\ v^3_x & v^3_y & v^3_z & 1 \\ v^4_x & v^4_y & v^4_z & 1 \\ \end{array}}\right] \right\vert \]

where $ v^i_x, v^i_y, v^i_z $ are the x, y, z component of point $ v^i $.

Note that volume and Jacobian of map $ \Phi: T^0 \to T $ are related as

\[ volume(T) = volume(T^0) \times det(J). \]

Here, $ volume(T^0) = 1/6 $.

Parameters
nodesVertices of element
Returns
vector Vector of shape functions at point p

Implements fe::BaseElem.

◆ getDerShapes() [1/2]

std::vector< std::vector< double > > fe::TetElem::getDerShapes ( const util::Point3 &  p)
overrideprivatevirtual

Returns the values of derivative of shape function at point p on reference element.

Parameters
pLocation of point
Returns
vector Vector of derivative of shape functions

Implements fe::BaseElem.

◆ getDerShapes() [2/2]

std::vector< std::vector< double > > fe::TetElem::getDerShapes ( const util::Point3 &  p,
const std::vector< util::Point3 > &  nodes 
)
overridevirtual

Returns the values of derivative of shape function at point p.

Below, we present the derivation of the formula:

We are interested in $ \frac{\partial N_i(x_p, y_p, z_p)}{\partial x}$, $ \frac{\partial N_i(x_p, y_p, z_p)}{\partial y}$ and $ \frac{\partial N_i(x_p, y_p, z_p)}{\partial z}$. By using the map $ \Phi : T^0 \to T$ we have

\[ N^0_i(\xi, \eta, \zeta) = N_i(x(\xi,\eta,\zeta), y(\xi, \eta,\eta), z(\xi, \eta,\eta)) \]

and therefore we can write

\[ \frac{\partial N^0_i(\xi, \eta, \zeta)}{\partial \xi} = \frac{\partial N_i}{\partial x} \frac{\partial x}{\partial \xi} + \frac{\partial N_i}{\partial y} \frac{\partial y}{\partial \xi} + \frac{\partial N_i}{\partial z} \frac{\partial z}{\partial \xi}, \]

\[ \frac{\partial N^0_i(\xi, \eta, \zeta)}{\partial \eta} = \frac{\partial N_i}{\partial x} \frac{\partial x}{\partial \eta} + \frac{\partial N_i}{\partial y} \frac{\partial y}{\partial \eta} + \frac{\partial N_i}{\partial z} \frac{\partial z}{\partial \eta} \]

and

\[ \frac{\partial N^0_i(\xi, \eta, \zeta)}{\partial \zeta} = \frac{\partial N_i}{\partial x} \frac{\partial x}{\partial \zeta} + \frac{\partial N_i}{\partial y} \frac{\partial y}{\partial \zeta} + \frac{\partial N_i}{\partial z} \frac{\partial z}{\partial \zeta} \]

which can be written in the matrix form as

\[ \left[ {\begin{array}{c} \frac{\partial N^0_i}{\partial \xi} \\ \frac{\partial N^0_i}{\partial \eta} \\ \frac{\partial N^0_i}{\partial \zeta} \end{array}}\right] = \left[ { \begin{array}{ccc} \frac{dx}{d\xi} &\frac{dy}{d\xi} &\frac{dz}{d\xi} \\ \frac{dx}{d\eta} & \frac{dy}{d\eta} &\frac{dz}{d\eta} \\ \frac{dx}{d\zeta} & \frac{dy}{d\zeta} &\frac{dz}{d\zeta} \\ \end{array} } \right] \, \left[ {\begin{array}{c} \frac{\partial N_i}{\partial x} \\ \frac{\partial N_i}{\partial y} \\ \frac{\partial N_i}{\partial z}\end{array}}\right]. \]

The matrix is the Jacobian matrix $ J $ and can be computed easily if vertices of elements are known. Using $ J^{-1} $ we have following formula for derivatives of the shape function

\[ \left[ {\begin{array}{c} \frac{\partial N_i}{\partial x} \\ \frac{\partial N_i}{\partial y} \\ \frac{\partial N_i}{\partial z}\end{array}}\right] = J^{-1} \left[ {\begin{array}{c} \frac{\partial N^0_i}{\partial \xi} \\ \frac{\partial N^0_i}{\partial \eta} \\ \frac{\partial N^0_i}{\partial \zeta}\end{array}}\right]. \]

Here, derivatives $ \frac{\partial N^0_i}{ \partial \xi} $, $ \frac{\partial N^0_i} {\partial \eta } $ and $ \frac{\partial N^0_i} {\partial \zeta } $ correspond to reference element and are easy to compute.

Parameters
pLocation of point
nodesVertices of element
Returns
vector Vector of derivative of shape functions

Reimplemented from fe::BaseElem.

Here is the call graph for this function:

◆ getJacobian()

double fe::TetElem::getJacobian ( const util::Point3 &  p,
const std::vector< util::Point3 > &  nodes,
std::vector< std::vector< double >> *  J 
)
overrideprivatevirtual

Computes the Jacobian of map $ \Phi: T^0 \to T $.

Parameters
pLocation of point in reference element
nodesVertices of element
JMatrix to store the Jacobian (if not nullptr)
Returns
det(J) Determinant of the Jacobain

Implements fe::BaseElem.

Here is the call graph for this function:

◆ getQuadDatas()

std::vector< fe::QuadData > fe::TetElem::getQuadDatas ( const std::vector< util::Point3 > &  nodes)
overridevirtual

Get vector of quadrature data.

Given element vertices, this method returns the list of quadrature point and essential quantities at quadrature points. Here, order of quadrature approximation is set in the constructor. List of data

  • quad point
  • quad weight
  • shape function evaluated at quad point
  • derivative of shape function evaluated at quad point
  • Jacobian matrix
  • determinant of the Jacobian

Let $ T $ is the given tetrahedron with vertices $ v^1, v^2, v^3, v^4 $ and let $ T^0 $ is the reference tetrahedron.

  1. To compute quadrature points, we first compute the quadrature points on reference tetrahedron $ T^0 $, and then we use the map $ \Phi: T^0 \to T$ to map the points on reference tetrahedron to the given triangle $ T $.
  2. To compute the quadrature weight, we compute the quadrature weight associated to the quadrature point in reference tetrahedron $ T^0 $. Suppose $ w^0_q $ is the quadrature weight associated to quadrature point $ (\xi_q, \eta_q) \in T^0 $, then the quadrature point $ w_q $ associated to the mapped point $ (x(\xi_q, \eta_q), y(\xi_q, \eta_q)) \in T $ is given by

    \[ w_q = w^0_q * det(J) \]

    where $ det(J) $ is the determinant of the Jacobian of map $ \Phi $.
  3. We compute shape functions $ N_1, N_2, N_3$ associated to $ T $ at the quadrature point $ (x(\xi_q, \eta_q), y(\xi_q,\eta_q)) $ using formula

    \[ N_i(x(\xi_q, \eta_q), y(\xi_q, \eta_q)) = N^0_i(\xi_q, \eta_q). \]

  4. To compute derivative of shape functions such as $ \frac{\partial N_i}{\partial x}, \frac{\partial N_i}{\partial y}$ associated to $ T $, we use the relation between derivatives of shape function in $ T $ and $ T^0 $ described in fe::TetElem::getDerShapes.
Parameters
nodesVector of vertices of an element
Returns
vector Vector of QuadData

Implements fe::BaseElem.

Here is the call graph for this function:

◆ getQuadPoints()

std::vector< fe::QuadData > fe::TetElem::getQuadPoints ( const std::vector< util::Point3 > &  nodes)
overridevirtual

Get vector of quadrature data.

Given element vertices, this method returns the list of quadrature point and essential quantities at quadrature points. Here, order of quadrature approximation is set in the constructor. List of data

  • quad point
  • quad weight
  • shape function evaluated at quad point

This function is lite version of fe::TetElem::getQuadDatas.

Parameters
nodesVector of vertices of an element
Returns
vector Vector of QuadData

Implements fe::BaseElem.

◆ getShapes() [1/2]

std::vector< double > fe::TetElem::getShapes ( const util::Point3 &  p)
overrideprivatevirtual

Returns the values of shape function at point p on reference element.

Parameters
pLocation of point
Returns
vector Vector of shape functions at point p

Implements fe::BaseElem.

◆ getShapes() [2/2]

std::vector< double > fe::TetElem::getShapes ( const util::Point3 &  p,
const std::vector< util::Point3 > &  nodes 
)
overridevirtual

Returns the values of shape function at point p.

We first map the point p in $ T $ to reference tetrahedron $ T^0 $ using fe::TetElem::mapPointToRefElem and then compute shape functions at the mapped point using fe::TetElem::getShapes(const util::Point3 &).

Parameters
pLocation of point
nodesVertices of element
Returns
vector Vector of shape functions at point p

Reimplemented from fe::BaseElem.

◆ mapPointToRefElem()

util::Point3 fe::TetElem::mapPointToRefElem ( const util::Point3 &  p,
const std::vector< util::Point3 > &  nodes 
)
overrideprivatevirtual

Maps point p in a given element to the reference element.

Let $ v^1, v^2, v^3$ are three vertices of triangle $ T$ and let $T^0 $ is the reference triangle. Following the introduction to fe::TriElem, the map $ (x,y) \in T $ to $ (\xi, \eta) \in T^0 $ is given by

\[ \xi = C_{11} (x - v^1_x) + C_{12} (y - v^1_y), \quad \eta = C_{21} (x - v^1_x) + C_{22} (y - v^1_y). \]

$ C$ is the inverse of matrix

\[ B = \left[ {\begin{array}{cc} v^2_x - v^1_x & v^3_x - v^1_x \\ v^2_y - v^1_y & v^3_y - v^1_y \end{array}}\right], \]

i.e.

\[ C := B^{-1} = \frac{1}{(v^2_x - v^1_x)(v^3_y - v^1_y) - (v^3_x - v^1_x) (v^2_y - v^1_y)} \left[ {\begin{array}{cc} v^3_y - v^1_y & -(v^3_x - v^1_x) \\ -(v^2_y - v^1_y) & v^2_x - v^1_x \end{array}}\right]. \]

If mapped point $ (\xi, \eta)$ does not satisfy following conditions

  • \[ 0\leq \xi, \eta \]

  • \[ \xi \leq 1 - \eta \quad (or\, equivalently) \quad \eta \leq 1 - \xi \]

    then the point $ (\xi,\eta) $ does not belong to the reference triangle or equivalently point $(x,y) $ does not belong to the triangle $ T $ and the method issues error. Otherwise the method returns point $ (\xi, \eta)$.
Parameters
pLocation of point
nodesVertices of element
Returns
vector Vector of shape functions at point p

Reimplemented from fe::BaseElem.

Here is the call graph for this function:

◆ print()

void fe::TetElem::print ( int  nt = 0,
int  lvl = 0 
) const
inline

Prints the information about the instance of the object.

Parameters
ntNumber of tabs to append before each line of string
lvlLevel of information sought (higher level means more information)
Here is the call graph for this function:

◆ printStr()

std::string fe::TetElem::printStr ( int  nt = 0,
int  lvl = 0 
) const

Returns the string containing information about the instance of the object.

Parameters
ntNumber of tabs to append before each line of string
lvlLevel of information sought (higher level means more information)
Returns
string String containing information about this object

Referenced by print().

Here is the call graph for this function:
Here is the caller graph for this function:
The documentation for this class was generated from the following files: