Working with elements

Table of Contents

• Introduction
• Volume elements
• Surface elements

Introduction

The purpose of this tutorial is to explain how you can work with finite elements in NiHu.

Volume elements

Elements are coordinate transforms \( {\bf x}({\bf \xi})\) that map from an intrinsic domain to the physical coordinate space. Therefore, elements are defined by an interpolation shape function set \( L_i(\xi) \) and nodal locations \( {\bf x}_i \).

A two-dimensional four noded quadrilateral volume element type can be defined as

    typedef NiHu::volume_element<NiHu::quad_1_shape_set, double> Element;

where NiHu::quad_1_shape_set refers to the four noded linear quadrilateral shape functions. The second template argument defines the scalar type of the physical coordinate space.

After the element type has been defined, an element instance can be constructed by providing the nodal locations to the constructor. The nodal locations need to be defined in a matrix of type NiHu::volume_element::coords_t. In our particular case, this is a 2x4 matrix of double elements, where each column stores a nodal coordinate \( {\bf x}_i \). The matrix is filled using Eigen's comma initialiser syntax.

    Element::coords_t coords;
    coords <<
        0.0, 1.0, 1.0, 0.0,
        0.0, 0.0, 1.0, 1.0;
 
    Element elem(coords);

Apparently, our element is a square volume \( 0 \le x,y \le 1 \).

The local intrinsic coordinate type of the element can be accessed as NiHu::volume_element::xi_t.

    Element::xi_t xi;
    xi << 0.0, 0.0;

The element's physical location can be obtained by calling the member NiHu::volume_element::get_x with an argument \( \xi \). In our case, the location is a 2x1 double vector.

    std::cout << "location:\n" << elem.get_x(xi) << std::endl;

The element location gradient can be obtained by calling member NiHu::volume_element::get_dx. In our case, this gradient is a 2x2 matrix, where the first column contains the derivative with respect to the first local coordinate \( {\bf x},_{\xi} \) and the second colum contains \( {\bf x},_{\eta} \). The element Jacobian can be obtained by evaluating the determinant of the gradient matrix.

    std::cout << "gradient:\n" << elem.get_dx(xi) << std::endl;
    std::cout << "jacobian:\n" << elem.get_dx(xi).determinant() << std::endl;

The second derivative of the coordinate mapping is returned by function NiHu::volume_element::get_ddx. In our case, this is a 2x3 matrix, where the first column contains the second derivative with respect to the fisrt intrinsic coordinate \( {\bf x},_{\xi\xi}\), the second column contains \( {\bf x},_{\xi\eta}\), and the third colum contains \( {\bf x},_{\eta\eta} \).

    std::cout << "2nd derivative:\n" << elem.get_ddx(xi) << std::endl;

Surface elements

Surface elements are handled in a similar fashion. The example below defines a quadrilateral surface element in the 3-dimensional physical space. Note that the physical space's dimensionality is automatically deduced from the dimensionality of the intrinsic domain.

    typedef NiHu::surface_element<NiHu::quad_1_shape_set, double> Element;
 
    Element::coords_t coords;
    coords <<
        0.0, 1.0, 1.0, 0.0,
        0.0, 0.0, 1.0, 1.0,
        0.0, 0.0, 0.0, 0.0;
 
    Element elem(coords);
 
    Element::xi_t xi;
    xi << 0.0, 0.0;
 
    std::cout << "location:\n" << elem.get_x(xi) << std::endl;

For the case of the surface element, the location gradient is a 3x2 matrix. The surface elements provides the function NiHu::surface_element::get_normal to return the normal vector. Note that this normal vector is not the unit normal. In three dimensions, its definition is \( {\bf n} = {\bf x},_{\xi} \times {\bf x},_{\eta}\). The Jacobian can be obtained as the norm of the normal vector.

    std::cout << "gradient:\n" << elem.get_dx(xi) << std::endl;
    std::cout << "normal:\n" << elem.get_normal(xi) << std::endl;
    std::cout << "jacobian:\n" << elem.get_normal(xi).norm() << std::endl;

The second derivative is obtained similarly like for volume elements.

    std::cout << "2nd derivative:\n" << elem.get_ddx(xi) << std::endl;