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
double	elemSize (const std::vector< util::Point3 > &nodes) override
	Returns the area of element. More...

std::vector< std::vector< double > >	getDerShapes (const util::Point3 &p, const std::vector< util::Point3 > &nodes) override
	Returns the values of derivative of shape function at point p. More...

std::vector< fe::QuadData >	getQuadDatas (const std::vector< util::Point3 > &nodes) override
	Get vector of quadrature data. More...

std::vector< fe::QuadData >	getQuadPoints (const std::vector< util::Point3 > &nodes) override
	Get vector of quadrature data. More...

std::vector< double >	getShapes (const util::Point3 &p, const std::vector< util::Point3 > &nodes) override
	Returns the values of shape function at point p. More...

void	print (int nt=0, int lvl=0) const
	Prints the information about the instance of the object. More...

std::string	printStr (int nt=0, int lvl=0) const
	Returns the string containing information about the instance of the object. More...

	TriElem (size_t order)
	Constructor. More...

Public Member Functions inherited from fe::BaseElem
	BaseElem (size_t order, size_t element_type)
	Constructor. More...

size_t	getElemType ()
	Get element type. More...

size_t	getNumQuadPoints ()
	Get number of quadrature points in the data. More...

size_t	getQuadOrder ()
	Get order of quadrature approximation. More...

void	print (int nt=0, int lvl=0) const
	Prints the information about the instance of the object. More...

std::string	printStr (int nt=0, int lvl=0) const
	Returns the string containing information about the instance of the object. More...

Private Member Functions
std::vector< std::vector< double > >	getDerShapes (const util::Point3 &p) override
	Returns the values of derivative of shape function at point p on reference element. More...

double	getJacobian (const util::Point3 &p, const std::vector< util::Point3 > &nodes, std::vector< std::vector< double >> *J) override
	Computes the Jacobian of map $\Phi: T^0 \to T$ . More...

std::vector< double >	getShapes (const util::Point3 &p) override
	Returns the values of shape function at point p on reference element. More...

void	init () override
	Compute the quadrature points for triangle element.

util::Point3	mapPointToRefElem (const util::Point3 &p, const std::vector< util::Point3 > &nodes) override
	Maps point p in a given element to the reference element. More...

Detailed Description

A class for mapping and quadrature related operations for linear 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 are vertices of element .
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 ) 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)$ .
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) $ .

Constructor & Destructor Documentation

◆ TriElem()

fe::TriElem::TriElem ( size_t order )

explicit

Constructor.

Parameters

order Order of quadrature point approximation

Member Function Documentation

◆ elemSize()

double fe::TriElem::elemSize ( const std::vector< util::Point3 > & 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.

◆ getDerShapes() [1/2]

std::vector< std::vector< double > > fe::TriElem::getDerShapes ( const util::Point3 & 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.

◆ getDerShapes() [2/2]

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

overridevirtual

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

Below, we present the derivation of the formula:

We are interested in $\frac{\partial N_i(x_p, y_p)}{\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.

◆ getJacobian()

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

overrideprivatevirtual

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

Parameters

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.

◆ getQuadDatas()

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

overridevirtual

Get vector of quadrature data.

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

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

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

To compute quadrature points, we first compute the quadrature points on reference triangle , and then we use the map $\Phi: T^0 \to T$ to map the points on reference triangle to the given triangle .
To compute the quadrature weight, we compute the quadrature weight associated to the quadrature point in reference triangle . Suppose is the quadrature weight associated to quadrature point $(\xi_q, \eta_q) \in T^0$ , then the quadrature point 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 is the determinant of the Jacobian of map $\Phi$ .
We compute shape functions associated to 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).$
To compute derivative of shape functions such as $\frac{\partial N_i}{\partial x}, \frac{\partial N_i}{\partial y}$ associated to , we use the relation between derivatives of shape function in and described in fe::TriElem::getDerShapes.

Parameters

nodes Vector of vertices of an element

Returns: vector Vector of QuadData

Implements fe::BaseElem.

◆ getQuadPoints()

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

overridevirtual

Get vector of quadrature data.

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

quad point
quad weight
shape function evaluated at quad point

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

Parameters

nodes Vector of vertices of an element

Returns: vector Vector of QuadData

Implements fe::BaseElem.

◆ getShapes() [1/2]

std::vector< double > fe::TriElem::getShapes ( const util::Point3 & 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.

◆ getShapes() [2/2]

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

overridevirtual

Returns the values of shape function at point p.

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

Parameters

p	Location of point
nodes	Vertices of element

Returns: vector Vector of shape functions at point p

Reimplemented from fe::BaseElem.

◆ mapPointToRefElem()

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

overrideprivatevirtual

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

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

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

$ C$ is the inverse of matrix

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

i.e.

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

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

$0\leq \xi, \eta$
$\xi \leq 1 - \eta \quad (or\, equivalently) \quad \eta \leq 1 - \xi$
then the point $(\xi,\eta)$ does not belong to the reference triangle or equivalently point does not belong to the triangle and the method issues error. Otherwise the method returns point $(\xi, \eta)$ .