fe::LineElem Class Reference

A class for mapping and quadrature related operations for linear 2-node line element. More...

#include <lineElem.h>

Inheritance diagram for fe::LineElem:

Collaboration diagram for fe::LineElem:

Public Member Functions
	LineElem (size_t order)
	Constructor for line element.
double	elemSize (const std::vector< util::Point > &nodes) override
	Returns the length 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 Jacobian of the map \( \Phi: T^0 \to T\).
void	init () override
	Compute the quadrature points for line 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 2-node line element.

The reference line element \(T^0 \) is given by vertices at \( v^1 = -1, v^2 = 1 \).

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

   \[N^0_1(\xi) = \frac{1 - \xi}{2}, \quad N^0_2(\xi) = \frac{1 + \xi}{2}. \]

2. Derivative of shape functions are constant and are as follows

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

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

   \[ x(\xi) = \sum_{i=1}^2 N^0_i(\xi) v^i_x, \]

   where \( v^1, v^2\) are vertices of element \( T \). For 1-d points, we simply have \( v^i_x = v^i \).
4. The Jacobian of the map \( \Phi: T^0 \to T \) is given by

   \[ J = \frac{dx}{d\xi}. \]

   Since it is a 1-d case, Jacobian and its determinant are same. For line element the formula for Jacobian is as follows

   \[ J = \frac{dx}{d\xi} = \frac{v^2_x - v^1_x}{2} = \frac{length(T) }{length(T^0)}. \]

5. Inverse map \( \Phi^{-1} \) from \( x \in T \) to \( \xi \in T^0\) is given by

   \[ \xi(x) = \frac{2}{v^2_x - v^1_x} (x - \frac{v^2_x + v^1_x}{2}) = \frac{1}{J}(x - \frac{v^2_x + v^1_x}{2}). \]

Definition at line 49 of file lineElem.h.

Constructor & Destructor Documentation

LineElem()

fe::LineElem::LineElem ( size_t  order )

explicit

Constructor for line element.

Parameters

order

Order of quadrature point approximation

Definition at line 16 of file lineElem.cpp.

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

References init().

Here is the call graph for this function:

Member Function Documentation

elemSize()

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

overridevirtual

Returns the length of element.

If line \( T \) is given by points \( v^1, v^2\) then the length is simply

\[ length(T) = v^2_x - v^1_x. \]

Parameters

nodes

Vertices of element

Returns

vector Vector of shape functions at point p

Implements fe::BaseElem.

Definition at line 23 of file lineElem.cpp.

                                                             {
  return (nodes[1].d_x - nodes[0].d_x);
}

getDerShapes() [1/2]

std::vector< std::vector< double > > fe::LineElem::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 104 of file lineElem.cpp.

                                           {
 
  // N1 = (1 - xi)/2
  // --> d N1/d xi = -1/2
  //
  // N2 = (1 + xi)/2
  // --> d N2/d xi = 1/2
  std::vector<std::vector<double>> r;
  r.push_back(std::vector<double>{-0.5});
  r.push_back(std::vector<double>{0.5});
  return r;
}

getDerShapes() [2/2]

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

overridevirtual

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

Let \( x\) is the point on line \( T \) and let \( \xi \) is the point on reference line \( T^0 \). Let shape functions on \( T\) are \( N_1, N_2 \) and on \( T^0 \) are \( N^0_1, N^0_2 \).

We are interested in \( \frac{\partial N_i(x_p)}{\partial x}\). By using the map \( \xi \to x \) we have

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

and therefore we can write

\[ \frac{\partial N^0_i(\xi)}{\partial \xi} = \frac{\partial N_i}{\partial x} \frac{\partial x}{\partial \xi}. \]

Since \( \frac{\partial x}{\partial \xi} = J \) is a Jacobian of map which can be computed easily if \( v^1, v^2 \) are known, we can invert and obtain the formula

\[ \frac{\partial N_i(\xi)}{\partial x} = \frac{1}{J} \frac{\partial N^0_i(\xi)}{\partial \xi}. \]

Parameters

p	Location of point
nodes	Vertices of element

Returns

vector Vector of derivative of shape functions

Reimplemented from fe::BaseElem.

Definition at line 34 of file lineElem.cpp.

                                                             {
  // get derivatives of shape function in reference triangle
  auto ders_ref = getDerShapes(mapPointToRefElem(p, nodes));
 
  // get Jacobian
  auto detJ = getJacobian(p, nodes, nullptr);
 
  // modify derivatives of shape function
  ders_ref[0][0] = ders_ref[0][0] / detJ;
  ders_ref[1][0] = ders_ref[1][0] / detJ;
 
  return ders_ref;
}

getJacobian()

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

overrideprivatevirtual

Computes Jacobian of the 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 (same as Jacobian in 1-d)

Implements fe::BaseElem.

Definition at line 140 of file lineElem.cpp.

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

getQuadDatas()

std::vector< fe::QuadData > fe::LineElem::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 for each quad point:

• 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 line with vertices \( v^1, v^2 \) and let \( T^0 \) is the reference line.

1. To compute quadrature point, we first compute the quadrature points on reference line \( T^0 \), and then we use the map \( \Phi: T^0 \to T\) to map the points on reference line to the current line \( T \).
2. To compute the quadrature weight, we compute the quadrature weight associated to the quadrature point in reference line \( T^0 \). Suppose \( w^0_q \) is the quadrature weight associated to quadrature point \( \xi_q \in T^0 \), then the quadrature point \( w_q \) associated to the mapped point \( x(\xi_q) \in T \) is given by

   \[ w_q = w^0_q * J \]

   where \( J \) is the Jacobian of map \( \Phi \).
3. We compute shape functions \( N_1, N_2\) associated to \( T \) at quadrature point \( x(\xi_q) \) using formula

   \[ N_i(x(\xi_q)) = N^0_i(\xi_q). \]

4. To compute the derivative of shape functions \( \frac{\partial N_1}{\partial x}, \frac{\partial N_2}{\partial x}\) associated to \( T \), we use the relation between derivatives of shape function in \( T \) and \( T^0 \) described in fe::LineElem::getDerShapes.

Parameters

nodes

Vector of vertices of an element

Returns

vector Vector of QuadData

Implements fe::BaseElem.

Definition at line 50 of file lineElem.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 line
    qd.d_p.d_x = qd.d_shapes[0] * nodes[0].d_x + qd.d_shapes[1] * nodes[1]
        .d_x;
 
    // modify derivative of shape function
    for (size_t i=0; i<2; i++)
      qd.d_derShapes[i][0] = qd.d_derShapes[i][0] / qd.d_detJ;
  }
 
  return qds;
}

getQuadPoints()

std::vector< fe::QuadData > fe::LineElem::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 for each quad point:

• quad point
• quad weight
• shape function evaluated at quad point

This function is a lite version of fe::LineElem::getQuadDatas.

Parameters

nodes

Vector of vertices of an element

Returns

vector Vector of QuadData

Implements fe::BaseElem.

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

getShapes() [1/2]

std::vector< double > fe::LineElem::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 96 of file lineElem.cpp.

                                                          {
 
  // N1 = (1 - xi)/2
  // N2 = (1 + xi)/2
  return std::vector<double>{0.5 * (1. - p.d_x), 0.5 * (1. + p.d_x)};
}

References util::Point::d_x.

getShapes() [2/2]

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

overridevirtual

Returns the values of shape function at point p.

Line \( T \) is given by points \( v^1, v^2\). We first map the point p in \( T \) to reference line \( T^0 \) using fe::LineElem::mapPointToRefElem and then compute shape functions at the mapped point using fe::LineElem::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 28 of file lineElem.cpp.

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

init()

void fe::LineElem::init ( )

overrideprivatevirtual

Compute the quadrature points for line element.

Implements fe::BaseElem.

Definition at line 154 of file lineElem.cpp.

                      {
 
  //
  // compute quad data for reference line element with vertex at p1 = -1 and
  // p2 = 1
  //
  // Shape functions are
  // N1 = (1 - xi)/2
  // N2 = (1 + xi)/2
 
  if (!d_quads.empty())
    return;
 
  // no point in zeroth order
  if (d_quadOrder == 0)
    d_quads.resize(0);
 
  // 1x1 identity matrix
  std::vector<std::vector<double>> ident_mat;
  ident_mat.push_back(std::vector<double>{1.});
 
  //
  // first order quad points
  //
  if (d_quadOrder == 1) {
    d_quads.clear();
    // 1-d points are: {0} and weights are: {2}
    fe::QuadData qd;
    qd.d_w = 2.;
    qd.d_p = util::Point();
    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);
  }
 
  //
  // second order quad points
  //
  if (d_quadOrder == 2) {
    d_quads.clear();
    // 1-d points are: {-1/sqrt{3], 1/sqrt{3}} and weights are: {1,1}
    fe::QuadData qd;
    qd.d_w = 1.;
    qd.d_p = util::Point(-1. / std::sqrt(3.), 0., 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);
 
    qd.d_w = 1.;
    qd.d_p = util::Point(1. / std::sqrt(3.), 0., 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();
    // 1-d points are: {-sqrt{3}/sqrt{5}, 0, sqrt{3}/sqrt{5}}
    // weights are: {5/9, 8/9, 5/9}
    fe::QuadData qd;
    qd.d_w = 5. / 9.;
    qd.d_p = util::Point(-std::sqrt(3.) / std::sqrt(5.), 0., 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);
 
    qd.d_w = 8. / 9.;
    qd.d_p = util::Point();
    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);
 
    qd.d_w = 5. / 9.;
    qd.d_p = util::Point(std::sqrt(3.) / std::sqrt(5.), 0., 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;
    qd.d_w = 0.6521451548625461;
    qd.d_p = util::Point(-0.3399810435848563, 0., 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);
 
    qd.d_w = 0.6521451548625461;
    qd.d_p = util::Point(0.3399810435848563, 0., 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);
 
    qd.d_w = 0.3478548451374538;
    qd.d_p = util::Point(-0.8611363115940526, 0., 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);
 
    qd.d_w = 0.3478548451374538;
    qd.d_p = util::Point(0.8611363115940526, 0., 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;
    qd.d_w = 0.5688888888888889;
    qd.d_p = util::Point();
    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);
 
    qd.d_w = 0.4786286704993665;
    qd.d_p = util::Point(-0.5384693101056831, 0., 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);
 
    qd.d_w = 0.4786286704993665;
    qd.d_p = util::Point(0.5384693101056831, 0., 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);
 
    qd.d_w = 0.2369268850561891;
    qd.d_p = util::Point(-0.9061798459386640, 0., 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);
 
    qd.d_w = 0.2369268850561891;
    qd.d_p = util::Point(0.9061798459386640, 0., 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 LineElem().

Here is the caller graph for this function:

mapPointToRefElem()

util::Point fe::LineElem::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\) are vertices of element \( T\) and let \(T^0 \) is the reference element. Map \(\Phi : T^0 \to T\) is given by

\[ \xi(x) = \frac{2}{v^2_x - v^1_x} (x - \frac{v^2_x + v^1_x}{2}) = \frac{1}{J}(x - \frac{v^2_x + v^1_x}{2}). \]

If mapped point \( \xi\) does not satisfy condition

• \[ -1 \leq \xi \leq 1 \]

  then the point \( \xi \) does not belong to reference line \( T^0\) or equivalently point \(x \) does not belong to line \( T \) and the method issues error. Otherwise the method returns point \( \xi\).

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 118 of file lineElem.cpp.

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

References util::Point::d_x, 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/lineElem.h
• /home/user/project/src/fe/lineElem.cpp