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.

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 triangle element.

The reference triangle element \( T^0 \) is given by vertices \( (0,0), \, (1,0), \, (0,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. \]
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. \]
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 \).
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,
\[ \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) \).
Inverse map \( \Phi^{-1} : T \to T^9 \) for linear triangle element can be derived as follows:

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 18 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 25 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 66 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 follows 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 37 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 111 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);
}

◆ 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 59 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 31 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 129 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 79 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)) {
    throw std::runtime_error(
        util::io::Msg()
        << "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");
  }
 
  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

Public Member Functions

Private Member Functions

Additional Inherited Members

Detailed Description

Constructor & Destructor Documentation

◆ TriElem()

Member Function Documentation

◆ elemSize()

◆ getDerShapes() [1/2]

◆ getDerShapes() [2/2]

◆ getJacobian()

◆ getShapes() [1/2]

◆ getShapes() [2/2]

◆ init()

◆ mapPointToRefElem()