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.
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 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_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, \( 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 18 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

\[ 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 25 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 59 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 follows from \( v^1, v^2 \), 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 36 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 96 of file lineElem.cpp.

                                                                 {
  // Reference line is [-1, 1], dN/dξ = ±1/2, so J = dx/dξ = L/2.
  const double Jxi = 0.5 * (nodes[1].d_x - nodes[0].d_x);
 
  if (J != nullptr) {
    J->resize(1);
    (*J)[0] = std::vector<double>{Jxi};
  }
 
  return Jxi;
}

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 51 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 30 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 110 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 73 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) ) {
    throw std::runtime_error(
        util::io::Msg()
        << "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");
  }
 
  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