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.
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.
std::vector< fe::QuadData >	getQuadDatas (const std::vector< util::Point > &nodes)
	Get quadrature data mapped to the physical element (N, dN/dx, x, w, J). Shared isoparametric map; types only fill reference N, dN/dξ, and w.
std::vector< fe::QuadData >	getQuadPoints (const std::vector< util::Point > &nodes)
	Same as getQuadDatas without transforming dN/dξ to dN/dx.

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_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,

   \[ \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 derived as follows:

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 70 of file tetElem.cpp.

    : fe::BaseElem(order, util::vtk_type_tetra) {
 
  if (d_quadOrder > 3) {
    throw std::runtime_error(
        util::io::Msg()
        << "Error: For linear tet element, we only support upto 3 quad "
        "order approximation.\n");
  }
 
  // 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 84 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 125 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 follows 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 98 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 176 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:

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 120 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 93 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 209 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 142 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