fe::TriElem Class Reference

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

#include <triElem.h>

Inheritance diagram for fe::TriElem:

Collaboration diagram for fe::TriElem:

Public Member Functions
	TriElem (size_t order)
	Constructor.
double	elemSize (const std::vector< util::Point > &nodes) override
	Returns the area 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 triangle element.

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

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

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

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

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

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

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

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

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

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

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

   and determinant of Jacobian is

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

   For linear triangle element, Jacobian (and so \( det(J) \)) is constant. For linear triangle 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, \]

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

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

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

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

Substituting formula for \( N^0_i \) in above to get

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

or

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

Writing above in matrix form, we have

\[ \left[ {\begin{array}{c} x - v^1_x \\ y - v^1_y \end{array}}\right] = \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] \, \left[ {\begin{array}{c} \xi \\ \eta \end{array}}\right]. \]

Denoting the matrix as \( B \)

\[ 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] \]

then

\[ C := B^{-1} = \frac{1}{det(B)} \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]. \]

With \( C \) we have inverse map given by

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

Note that matrix \( B \) is a transpose of the Jacobian of map \( \Phi\) and hence \( det(B) = det(J) \).

Definition at line 91 of file triElem.h.

Constructor & Destructor Documentation

TriElem()

fe::TriElem::TriElem ( size_t  order )

explicit

Constructor.

Parameters

order

Order of quadrature point approximation

Definition at line 16 of file triElem.cpp.

    : fe::BaseElem(order, util::vtk_type_triangle) {
 
  // compute quad data
  this->init();
}

References init().

Here is the call graph for this function:

Member Function Documentation

elemSize()

double fe::TriElem::elemSize ( const std::vector< util::Point > &  nodes )

overridevirtual

Returns the area of element.

If triangle \( T \) is given by points \( v^1, v^2, v^3 \) then the area is

\[ area(T) = \frac{(v^2_x - v^1_x) (v^3_y - v^1_y) - (v^3_x - v^1_x) (v^2_y - v^1_y)}{2}, \]

where \( v^i_x, v^i_y \) are the x and y component of point \( v^i \).

Note that area and Jacobian of map \( \Phi: T^0 \to T \) are related as

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

Here, \( area(T^0) = 0.5 \).

Parameters

nodes

Vertices of element

Returns

vector Vector of shape functions at point p

Implements fe::BaseElem.

Definition at line 23 of file triElem.cpp.

                                                            {
  return 0.5 * ((nodes[1].d_x - nodes[0].d_x) * (nodes[2].d_y - nodes[0].d_y) -
                (nodes[2].d_x - nodes[0].d_x) * (nodes[1].d_y - nodes[0].d_y));
}

getDerShapes() [1/2]

std::vector< std::vector< double > > fe::TriElem::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 126 of file triElem.cpp.

                                          {
 
  // 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
  std::vector<std::vector<double>> r;
  r.push_back(std::vector<double>{-1., -1.});
  r.push_back(std::vector<double>{1., 0.});
  r.push_back(std::vector<double>{0., 1.});
 
  return r;
}

getDerShapes() [2/2]

std::vector< std::vector< double > > fe::TriElem::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)}{\partial x}\) and \( \frac{\partial N_i(x_p, y_p)}{\partial y}\). By using the map \( \Phi : T^0 \to T\) we have

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

and therefore we can write

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

and

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

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} \end{array}}\right] = \left[ {\begin{array}{cc} \frac{\partial x}{\partial \xi} & \frac{\partial y}{\partial \xi} \\ \frac{\partial x}{\partial \eta} & \frac{\partial y}{\partial \eta} \end{array}}\right] \, \left[ {\begin{array}{c} \frac{\partial N_i}{\partial x} \\ \frac{\partial N_i}{\partial y} \end{array}}\right]. \]

The matrix is the Jacobian matrix \( J \) and can be computed easily if vertices of elements are known. Inverse \( J^{-1} \) is given by

\[ J^{-1} = \frac{1}{det(J)} \left[ {\begin{array}{cc} \frac{\partial y}{\partial \eta} & -\frac{\partial y}{\partial \xi} \\ -\frac{\partial x}{\partial \eta} & \frac{\partial x}{\partial \xi} \end{array}}\right]. \]

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} \end{array}}\right] = J^{-1} \left[ {\begin{array}{c} \frac{\partial N^0_i}{\partial \xi} \\ \frac{\partial N^0_i}{\partial \eta} \end{array}}\right]. \]

Here, derivatives \( \frac{\partial N^0_i}{ \partial \xi} \) and \( \frac{\partial N^0_i} {\partial \eta } \) 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 35 of file triElem.cpp.

                                                             {
  // get derivatives of shape function in reference triangle
  auto ders_ref = getDerShapes(mapPointToRefElem(p, nodes));
 
  // get Jacobian and its determinant
  std::vector<std::vector<double>> J;
  auto detJ = getJacobian(p, nodes, &J);
 
  // to hold derivatives
  auto ders = ders_ref;
 
  for (size_t i=0; i<3; i++) {
    // partial N_i/ partial x
    ders[i][0] = (ders_ref[i][0] * J[1][1] - ders_ref[i][1] * J[0][1]) / detJ;
    // partial N_i/ partial y
    ders[i][1] = (-ders_ref[i][0] * J[1][0] + ders_ref[i][1] * J[0][0]) / detJ;
  }
 
  return ders;
}

getJacobian()

double fe::TriElem::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 170 of file triElem.cpp.

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

getQuadDatas()

std::vector< fe::QuadData > fe::TriElem::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 triangle with vertices \( v^1, v^2, v^3 \) and let \( T^0 \) is the reference triangle.

1. To compute quadrature points, we first compute the quadrature points on reference triangle \( T^0 \), and then we use the map \( \Phi: T^0 \to T\) to map the points on reference triangle to the given triangle \( T \).
2. To compute the quadrature weight, we compute the quadrature weight associated to the quadrature point in reference triangle \( 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::TriElem::getDerShapes.

Parameters

nodes

Vector of vertices of an element

Returns

vector Vector of QuadData

Implements fe::BaseElem.

Definition at line 58 of file triElem.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.d_x = qd.d_shapes[0] * nodes[0].d_x + qd.d_shapes[1] * nodes[1].d_x +
        qd.d_shapes[2] * nodes[2].d_x;
    qd.d_p.d_y = qd.d_shapes[0] * nodes[0].d_y + qd.d_shapes[1] * nodes[1].d_y +
        qd.d_shapes[2] * nodes[2].d_y;
 
    // derivatives of shape function
    std::vector<std::vector<double>> ders;
 
    for (size_t i=0; i<3; i++) {
      // partial N_i/ partial x
      auto d1 = (qd.d_derShapes[i][0] * qd.d_J[1][1] - qd.d_derShapes[i][1] *
          qd.d_J[0][1]) / qd.d_detJ;
      // partial N_i/ partial y
      auto d2 = (-qd.d_derShapes[i][0] * qd.d_J[1][0] + qd.d_derShapes[i][1] *
          qd.d_J[0][0]) / qd.d_detJ;
 
      ders.push_back(std::vector<double>{d1, d2});
    }
    qd.d_derShapes = ders;
  }
 
  return qds;
}

getQuadPoints()

std::vector< fe::QuadData > fe::TriElem::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::TriElem::getQuadDatas.

Parameters

nodes

Vector of vertices of an element

Returns

vector Vector of QuadData

Implements fe::BaseElem.

Definition at line 98 of file triElem.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.d_x = qd.d_shapes[0] * nodes[0].d_x + qd.d_shapes[1] * nodes[1].d_x +
        qd.d_shapes[2] * nodes[2].d_x;
    qd.d_p.d_y = qd.d_shapes[0] * nodes[0].d_y + qd.d_shapes[1] * nodes[1].d_y +
        qd.d_shapes[2] * nodes[2].d_y;
  }
 
  return qds;
}

getShapes() [1/2]

std::vector< double > fe::TriElem::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 119 of file triElem.cpp.

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

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

getShapes() [2/2]

std::vector< double > fe::TriElem::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 triangle \( T^0 \) using fe::TriElem::mapPointToRefElem and then compute shape functions at the mapped point using fe::TriElem::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 29 of file triElem.cpp.

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

init()

void fe::TriElem::init ( )

overrideprivatevirtual

Compute the quadrature points for triangle element.

Implements fe::BaseElem.

Definition at line 188 of file triElem.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);
  }
 
  // 2x2 identity matrix
  std::vector<std::vector<double>> ident_mat;
  ident_mat.push_back(std::vector<double>{1., 0.});
  ident_mat.push_back(std::vector<double>{0., 1.});
 
  //
  // first order quad points for triangle
  //
  if (d_quadOrder == 1) {
    d_quads.clear();
    fe::QuadData qd;
    qd.d_w = 0.5;
    qd.d_p = util::Point(1. / 3., 1. / 3., 0.);
    // 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;
    // point 1
    qd.d_w = 1. / 6.;
    qd.d_p = util::Point(1. / 6., 1. / 6., 0.);
    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 = 1. / 6.;
    qd.d_p = util::Point(2. / 3., 1. / 6., 0.);
    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 = 1. / 6.;
    qd.d_p = util::Point(1. / 6., 2. / 3., 0.);
    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;
    // point 1
    qd.d_w = -27. / 96.;
    qd.d_p = util::Point(1. / 3., 1. / 3., 0.);
    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 = 25. / 96.;
    qd.d_p = util::Point(1. / 5., 3. / 5., 0.);
    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 = 25. / 96.;
    qd.d_p = util::Point(1. / 5., 1. / 5., 0.);
    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 = 25. / 96.;
    qd.d_p = util::Point(3. / 5., 1. / 5., 0.);
    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);
  }
 
  //
  // fourth order quad points for triangle
  //
  if (d_quadOrder == 4) {
    d_quads.clear();
    fe::QuadData qd;
    // point 1
    qd.d_w = 0.5 * 0.22338158967801;
    qd.d_p = util::Point(0.44594849091597, 0.44594849091597, 0.);
    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 = 0.5 * 0.22338158967801;
    qd.d_p = util::Point(0.44594849091597, 0.10810301816807, 0.);
    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 = 0.5 * 0.22338158967801;
    qd.d_p = util::Point(0.10810301816807, 0.44594849091597, 0.);
    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 = 0.5 * 0.10995174365532;
    qd.d_p = util::Point(0.09157621350977, 0.09157621350977, 0.);
    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 = 0.5 * 0.10995174365532;
    qd.d_p = util::Point(0.09157621350977, 0.81684757298046, 0.);
    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 6
    qd.d_w = 0.5 * 0.10995174365532;
    qd.d_p = util::Point(0.81684757298046, 0.09157621350977, 0.);
    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);
  }
 
  //
  // fifth order quad points for triangle
  //
  if (d_quadOrder == 5) {
    d_quads.clear();
    fe::QuadData qd;
    // point 1
    qd.d_w = 0.5 * 0.22500000000000;
    qd.d_p = util::Point(0.33333333333333, 0.33333333333333, 0.);
    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 = 0.5 * 0.13239415278851;
    qd.d_p = util::Point(0.47014206410511, 0.47014206410511, 0.);
    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 = 0.5 * 0.13239415278851;
    qd.d_p = util::Point(0.47014206410511, 0.05971587178977, 0.);
    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 = 0.5 * 0.13239415278851;
    qd.d_p = util::Point(0.05971587178977, 0.47014206410511, 0.);
    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 = 0.5 * 0.12593918054483;
    qd.d_p = util::Point(0.10128650732346, 0.10128650732346, 0.);
    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 6
    qd.d_w = 0.5 * 0.12593918054483;
    qd.d_p = util::Point(0.10128650732346, 0.79742698535309, 0.);
    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 7
    qd.d_w = 0.5 * 0.12593918054483;
    qd.d_p = util::Point(0.79742698535309, 0.10128650732346, 0.);
    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 TriElem().

Here is the caller graph for this function:

mapPointToRefElem()

util::Point fe::TriElem::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 139 of file triElem.cpp.

                                                                  {
  auto detB = 2. * elemSize(nodes);
  auto xi = ((nodes[2].d_y - nodes[0].d_y) * (p.d_x - nodes[0].d_x) -
             (nodes[2].d_x - nodes[0].d_x) * (p.d_y - nodes[0].d_y)) /
      detB;
  auto eta = (-(nodes[1].d_y - nodes[0].d_y) * (p.d_x - nodes[0].d_x) +
              (nodes[1].d_x - nodes[0].d_x) * (p.d_y - nodes[0].d_y)) /
      detB;
 
  if (util::isLess(xi, -1.0E-5) || util::isLess(eta, -1.0E-5) ||
      util::isGreater(xi, 1. + 1.0E-5 - eta)) {
    std::cerr << "Error: Trying to map point p = (" << p.d_x << ", " << p.d_y
              << ") in triangle to reference triangle.\n"
              << "But the point p does not belong to triangle = {("
              << nodes[0].d_x << ", " << nodes[0].d_y << "), (" << nodes[1].d_x
              << "," << nodes[1].d_y << "), (" << nodes[2].d_x << ","
              << nodes[2].d_y << ")}.\n"
              << "Coordinates in reference triangle are: xi = " << xi
              << ", eta = " << eta << "\n";
    exit(1);
  }
 
  if (util::isLess(xi, 0.))
    xi = 0.;
  if (util::isLess(eta, 0.))
    eta = 0.;
 
  return {xi, eta, 0.};
}

References util::Point::d_x, util::Point::d_y, util::isGreater(), and util::isLess().

Here is the call graph for this function:

---

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

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