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:

Collaboration diagram for fe::TetElem:

Public Member Functions
	TetElem (size_t order)
	Constructor.
double	elemSize (const std::vector< util::Point > &nodes) override
	Returns the volume of element.
std::vector< double >	getShapes (const util::Point &p, const std::vector< util::Point > &nodes) override
	Returns the values of shape function at point p.
std::vector< std::vector< double > >	getDerShapes (const util::Point &p, const std::vector< util::Point > &nodes) override
	Returns the values of derivative of shape function at point p.
std::vector< fe::QuadData >	getQuadDatas (const std::vector< util::Point > &nodes) override
	Get vector of quadrature data.
std::vector< fe::QuadData >	getQuadPoints (const std::vector< util::Point > &nodes) override
	Get vector of quadrature data.
Public Member Functions inherited from fe::BaseElem
	BaseElem (size_t order, size_t element_type)
	Constructor.
size_t	getElemType ()
	Get element type.
size_t	getQuadOrder ()
	Get order of quadrature approximation.
size_t	getNumQuadPoints ()
	Get number of quadrature points in the data.

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

Additional Inherited Members
Protected Attributes inherited from fe::BaseElem
size_t	d_quadOrder
	Order of quadrature point integration approximation.
size_t	d_numQuadPts
	Number of quadrature points for order d_quadOrder.
size_t	d_elemType
	Element type.
std::vector< fe::QuadData >	d_quads
	Quadrature data collection.

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 jacobian 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). \]

Definition at line 141 of file tetElem.h.

Constructor & Destructor Documentation

TetElem()

fe::TetElem::TetElem ( size_t  order )

explicit

Constructor.

Parameters

order

Order of quadrature point approximation

Definition at line 68 of file tetElem.cpp.

    : fe::BaseElem(order, util::vtk_type_triangle) {
 
  if (d_quadOrder > 3) {
    std::cout << "Error: For linear tet element, we only support upto 3 quad "
                 "order approximation.\n";
    exit(1);
  }
 
  // compute quad data
  this->init();
}

References fe::BaseElem::d_quadOrder, and init().

Here is the call graph for this function:

Member Function Documentation

elemSize()

double fe::TetElem::elemSize ( const std::vector< util::Point > &  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

nodes

Vertices of element

Returns

vector Vector of shape functions at point p

Implements fe::BaseElem.

Definition at line 81 of file tetElem.cpp.

                                                            {
  // volume of tet element is (1/6) a * (b x c),
  // where a = v2 - v1, b = v3 - v1, c = v4 - v1
  auto a = nodes[1] - nodes[0];
  auto b = nodes[2] - nodes[0];
  auto c = nodes[3] - nodes[0];
  return (1. / 6.) * a * (b.cross(c));
}

getDerShapes() [1/2]

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

overrideprivatevirtual

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

Parameters

p	Location of point

Returns

vector Vector of derivative of shape functions

Implements fe::BaseElem.

Definition at line 186 of file tetElem.cpp.

                        {
 
  // d N1/d xi = -1, d N1/d eta = -1, d N1/d zeta = -1,
  // d N2/ d xi = 1, d N2/d eta = 0, d N2/d zeta = 0,
  // d N3/ d xi = 0, d N3/d eta = 1, d N3/d zeta = 0,
  // d N4/ d xi = 0, d N4/d eta = 0, d N4/d zeta = 1,
  std::vector<std::vector<double>> r;
  r.push_back(std::vector<double>{-1., -1., -1.});
  r.push_back(std::vector<double>{1., 0., 0.});
  r.push_back(std::vector<double>{0., 1., 0.});
  r.push_back(std::vector<double>{0., 0., 1.});
 
  return r;
}

getDerShapes() [2/2]

std::vector< std::vector< double > > fe::TetElem::getDerShapes ( const util::Point &  p, const std::vector< util::Point > &  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

p	Location of point
nodes	Vertices of element

Returns

vector Vector of derivative of shape functions

Reimplemented from fe::BaseElem.

Definition at line 95 of file tetElem.cpp.

                                                           {
 
  // get derivatives of shape function in reference tet element
  auto ders_ref = getDerShapes(mapPointToRefElem(p, nodes));
 
  // get Jacobian and its determinant
  std::vector<std::vector<double>> J;
  auto detJ = getJacobian(p, nodes, &J);
 
  auto J_inv = util::inv(J);
 
  // to hold derivatives
  std::vector<std::vector<double>> ders(ders_ref.size(),
                                        std::vector<double>(3, 0.));
 
  // grad N_i = J_inv * grad N_i^ref
  for (size_t i = 0; i < 4; i++)
    ders[i] = util::dot(J_inv, ders_ref[i]);
 
  return ders;
}

References util::dot(), and util::inv().

Here is the call graph for this function:

getJacobian()

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

overrideprivatevirtual

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

Parameters

p	Location of point in reference element
nodes	Vertices of element
J	Matrix to store the Jacobian (if not nullptr)

Returns

det(J) Determinant of the Jacobain

Implements fe::BaseElem.

Definition at line 237 of file tetElem.cpp.

                                                                 {
  if (J != nullptr) {
    J->resize(3);
    (*J)[0] = std::vector<double>{nodes[1].d_x - nodes[0].d_x,
                                  nodes[1].d_y - nodes[0].d_y,
                                  nodes[1].d_z - nodes[0].d_z};
    (*J)[1] = std::vector<double>{nodes[2].d_x - nodes[0].d_x,
                                  nodes[2].d_y - nodes[0].d_y,
                                  nodes[2].d_z - nodes[0].d_z};
    (*J)[2] = std::vector<double>{nodes[3].d_x - nodes[0].d_x,
                                  nodes[3].d_y - nodes[0].d_y,
                                  nodes[3].d_z - nodes[0].d_z};
 
    return util::det(*J);
  } else {
    std::vector<std::vector<double>> J_local;
    J_local.resize(3);
    J_local[0] = std::vector<double>{nodes[1].d_x - nodes[0].d_x,
                                  nodes[1].d_y - nodes[0].d_y,
                                  nodes[1].d_z - nodes[0].d_z};
    J_local[1] = std::vector<double>{nodes[2].d_x - nodes[0].d_x,
                                  nodes[2].d_y - nodes[0].d_y,
                                  nodes[2].d_z - nodes[0].d_z};
    J_local[2] = std::vector<double>{nodes[3].d_x - nodes[0].d_x,
                                  nodes[3].d_y - nodes[0].d_y,
                                  nodes[3].d_z - nodes[0].d_z};
 
    return util::det(J_local);
  }
}

References util::det().

Here is the call graph for this function:

getQuadDatas()

std::vector< fe::QuadData > fe::TetElem::getQuadDatas ( const std::vector< util::Point > &  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

nodes

Vector of vertices of an element

Returns

vector Vector of QuadData

Implements fe::BaseElem.

Definition at line 118 of file tetElem.cpp.

                                       {
 
  // copy quad data associated to reference element
  auto qds = d_quads;
 
  // modify data
  for (auto &qd : qds) {
 
    // get Jacobian and determinant
    qd.d_detJ = getJacobian(qd.d_p, nodes, &(qd.d_J));
 
    // transform quad weight
    qd.d_w *= qd.d_detJ;
 
    // map point to triangle
    qd.d_p = util::Point();
    for (size_t i=0; i<4; i++) {
      qd.d_p.d_x += qd.d_shapes[i] * nodes[i].d_x;
      qd.d_p.d_y += qd.d_shapes[i] * nodes[i].d_y;
      qd.d_p.d_z += qd.d_shapes[i] * nodes[i].d_z;
    }
 
    // get inverse of Jacobian
    auto J_inv = util::inv(qd.d_J);
 
    // to hold derivatives
    std::vector<std::vector<double>> ders(qd.d_derShapes.size(),
                                          std::vector<double>(3, 0.));
 
    // grad N_i = J_inv * grad N_i^ref
    for (size_t i = 0; i < 4; i++)
      ders[i] = util::dot(J_inv, qd.d_derShapes[i]);
 
    qd.d_derShapes = ders;
  }
 
  return qds;
}

References util::dot(), and util::inv().

Here is the call graph for this function:

getQuadPoints()

std::vector< fe::QuadData > fe::TetElem::getQuadPoints ( const std::vector< util::Point > &  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

nodes

Vector of vertices of an element

Returns

vector Vector of QuadData

Implements fe::BaseElem.

Definition at line 158 of file tetElem.cpp.

                                       {
  // copy quad data associated to reference element
  auto qds = d_quads;
 
  // modify data
  for (auto &qd : qds) {
 
    // transform quad weight
    qd.d_w *= getJacobian(qd.d_p, nodes, nullptr);
 
    // map point to triangle
    qd.d_p = util::Point();
    for (size_t i=0; i<4; i++) {
      qd.d_p.d_x += qd.d_shapes[i] * nodes[i].d_x;
      qd.d_p.d_y += qd.d_shapes[i] * nodes[i].d_y;
      qd.d_p.d_z += qd.d_shapes[i] * nodes[i].d_z;
    }
  }
 
  return qds;
}

getShapes() [1/2]

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

overrideprivatevirtual

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

Parameters

p	Location of point

Returns

vector Vector of shape functions at point p

Implements fe::BaseElem.

Definition at line 181 of file tetElem.cpp.

                                                         {
  // N1 = 1 - xi - eta - zeta, N2 = xi, N3 = eta, N4 = zeta
  return std::vector<double>{1. - p.d_x - p.d_y - p.d_z, p.d_x, p.d_y, p.d_z};
}

References util::Point::d_x, util::Point::d_y, and util::Point::d_z.

getShapes() [2/2]

std::vector< double > fe::TetElem::getShapes ( const util::Point &  p, const std::vector< util::Point > &  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::Point &).

Parameters

p	Location of point
nodes	Vertices of element

Returns

vector Vector of shape functions at point p

Reimplemented from fe::BaseElem.

Definition at line 90 of file tetElem.cpp.

                                                           {
  return getShapes(mapPointToRefElem(p, nodes));
}

init()

void fe::TetElem::init ( )

overrideprivatevirtual

Compute the quadrature points for triangle element.

Implements fe::BaseElem.

Definition at line 270 of file tetElem.cpp.

                     {
  //
  // compute quad data for reference triangle with vertex at
  // (0,0), (1,0), (0,1)
  //
 
  if (!d_quads.empty()) return;
 
  // no point in zeroth order
  if (d_quadOrder == 0) {
    d_quads.resize(0);
  }
 
  // 3x3 identity matrix
  std::vector<std::vector<double>> ident_mat;
  ident_mat.push_back(std::vector<double>{1., 0., 0.});
  ident_mat.push_back(std::vector<double>{0., 1., 0.});
  ident_mat.push_back(std::vector<double>{0., 0., 1.});
 
  //
  // These datas are from LibMesh code
  // See: https://libmesh.github.io/doxygen/quadrature__gauss__3D_8C_source.html
 
  //
  // first order quad points for triangle
  //
  if (d_quadOrder == 1) {
    d_quads.clear();
    fe::QuadData qd;
    qd.d_w = 1. / 6.;
    qd.d_p = util::Point(1. / 4., 1. / 4., 1. / 4.);
    // N1 = 1 - xi - eta, N2 = xi, N3 = eta
    qd.d_shapes = getShapes(qd.d_p);
    // d N1/d xi = -1, d N1/d eta = -1, d N2/ d xi = 1, d N2/d eta = 0,
    // d N3/ d xi = 0, d N3/d eta = 1
    qd.d_derShapes = getDerShapes(qd.d_p);
    qd.d_J = ident_mat;
    qd.d_detJ = 1.;
    d_quads.push_back(qd);
  }
 
  //
  // second order quad points for triangle
  //
  if (d_quadOrder == 2) {
    d_quads.clear();
    fe::QuadData qd;
 
    double w = 1. / 24.;
    double a = 0.585410196624969;
    double b = 0.138196601125011;
    // point 1
    qd.d_w = w;
    qd.d_p = util::Point(a, b, b);
    qd.d_shapes = getShapes(qd.d_p);
    qd.d_derShapes = getDerShapes(qd.d_p);
    qd.d_J = ident_mat;
    qd.d_detJ = 1.;
    d_quads.push_back(qd);
    // point 2
    qd.d_w = w;
    qd.d_p = util::Point(b, a, b);
    qd.d_shapes = getShapes(qd.d_p);
    qd.d_derShapes = getDerShapes(qd.d_p);
    qd.d_J = ident_mat;
    qd.d_detJ = 1.;
    d_quads.push_back(qd);
    // point 3
    qd.d_w = w;
    qd.d_p = util::Point(b, b, a);
    qd.d_shapes = getShapes(qd.d_p);
    qd.d_derShapes = getDerShapes(qd.d_p);
    qd.d_J = ident_mat;
    qd.d_detJ = 1.;
    d_quads.push_back(qd);
    // point 4
    qd.d_w = w;
    qd.d_p = util::Point(b, b, b);
    qd.d_shapes = getShapes(qd.d_p);
    qd.d_derShapes = getDerShapes(qd.d_p);
    qd.d_J = ident_mat;
    qd.d_detJ = 1.;
    d_quads.push_back(qd);
  }
 
  //
  // third order quad points for triangle
  //
  if (d_quadOrder == 3) {
    d_quads.clear();
    fe::QuadData qd;
 
    double w1 = -2. / 15.;
    double w2 = 0.075;
 
    double a = 0.25;
    double b = 0.5;
    double c = 1. / 6.;
 
    // point 1
    qd.d_w = w1;
    qd.d_p = util::Point(a, a, a);
    qd.d_shapes = getShapes(qd.d_p);
    qd.d_derShapes = getDerShapes(qd.d_p);
    qd.d_J = ident_mat;
    qd.d_detJ = 1.;
    d_quads.push_back(qd);
    // point 2
    qd.d_w = w2;
    qd.d_p = util::Point(b, c, c);
    qd.d_shapes = getShapes(qd.d_p);
    qd.d_derShapes = getDerShapes(qd.d_p);
    qd.d_J = ident_mat;
    qd.d_detJ = 1.;
    d_quads.push_back(qd);
    // point 3
    qd.d_w = w2;
    qd.d_p = util::Point(c, b, c);
    qd.d_shapes = getShapes(qd.d_p);
    qd.d_derShapes = getDerShapes(qd.d_p);
    qd.d_J = ident_mat;
    qd.d_detJ = 1.;
    d_quads.push_back(qd);
    // point 4
    qd.d_w = w2;
    qd.d_p = util::Point(c, c, b);
    qd.d_shapes = getShapes(qd.d_p);
    qd.d_derShapes = getDerShapes(qd.d_p);
    qd.d_J = ident_mat;
    qd.d_detJ = 1.;
    d_quads.push_back(qd);
    // point 5
    qd.d_w = w2;
    qd.d_p = util::Point(c, c, c);
    qd.d_shapes = getShapes(qd.d_p);
    qd.d_derShapes = getDerShapes(qd.d_p);
    qd.d_J = ident_mat;
    qd.d_detJ = 1.;
    d_quads.push_back(qd);
  }
}

References fe::QuadData::d_derShapes, fe::QuadData::d_detJ, fe::QuadData::d_J, fe::QuadData::d_p, fe::QuadData::d_shapes, and fe::QuadData::d_w.

Referenced by TetElem().

Here is the caller graph for this function:

mapPointToRefElem()

util::Point fe::TetElem::mapPointToRefElem ( const util::Point &  p, const std::vector< util::Point > &  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

p	Location of point
nodes	Vertices of element

Returns

vector Vector of shape functions at point p

Reimplemented from fe::BaseElem.

Definition at line 203 of file tetElem.cpp.

                                                           {
 
  // get Jacobian matrix and compute its transpose
  std::vector<std::vector<double>> J(3, std::vector<double>(3, 0.));
  auto detJ = getJacobian(p, nodes, &J);
 
  // get transpose of Jacobian
  auto B = util::transpose(J);
  auto detB = detJ;
 
  // get inverse of B
  auto B_inv = util::inv(B);
 
  // get vector from first vertex to point p
  std::vector<double> vec_p = {p.d_x - nodes[0].d_x,
                               p.d_y - nodes[0].d_y,
                               p.d_z - nodes[0].d_z};
  // multiply B_inv to vector to transform point
  auto p_ref = util::dot(B_inv, vec_p);
 
  // check point
  checkPoint(p_ref, nodes);
 
  if (util::isLess(p_ref[0], 0.)) p_ref[0] = 0.;
  if (util::isLess(p_ref[1], 0.)) p_ref[1] = 0.;
  if (util::isLess(p_ref[2], 0.)) p_ref[2] = 0.;
  if (util::isGreater(p_ref[0], 1.)) p_ref[0] = 1.;
  if (util::isGreater(p_ref[1], 1.)) p_ref[1] = 1.;
  if (util::isGreater(p_ref[2], 1.)) p_ref[2] = 1.;
 
  return util::Point(p_ref);
}

References util::Point::d_x, util::Point::d_y, util::Point::d_z, util::dot(), util::inv(), util::isGreater(), util::isLess(), and util::transpose().

Here is the call graph for this function:

---

The documentation for this class was generated from the following files:

• /home/user/project/src/fe/tetElem.h
• /home/user/project/src/fe/tetElem.cpp