Introduction

The purpose of this tutorial is to demonstrate how a custom element/field type can be introduced into the NiHu toolbox.

We define a custom 3-dimensional four noded quadrangle surface field. The geometrical representation will be a standard 4-noded linear quadrilateral element (NiHu::quad_1_elem) having its nodes at the four element corners. The field is going to be defined so that its shape function nodes are located at the Gaussian nodes of the element. The new field is shown in the figure below.

Geometrical (red) and field (blue) nodal locations of the four noded quadrilateral element in the reference domain

In the following sections we define the new shape function set, introduce the new field type, and present an example how the new field type is used in a collocational or Galerkin type BEM.

The Gaussian shape set

The Gaussian shape set is defined over the quadrilateral domain NiHu::quad_domain located between coordinates \((-1,-1)\) and \((+1,+1)\).

The four nodal locations are

\( \displaystyle \left(\xi_i,\eta_i\right) = \left(\pm\sqrt{3}/3, \pm\sqrt{3}/3\right)\)

The shape functions are defined as

\( \displaystyle L_1(\xi,\eta) = (1-\sqrt{3}\xi)(1-\sqrt{3}\eta)/4 \\ \displaystyle L_2(\xi,\eta) = (1+\sqrt{3}\xi)(1-\sqrt{3}\eta)/4 \\ \displaystyle L_3(\xi,\eta) = (1+\sqrt{3}\xi)(1+\sqrt{3}\eta)/4 \\ \displaystyle L_4(\xi,\eta) = (1-\sqrt{3}\xi)(1+\sqrt{3}\eta)/4\)

and their derivatives with respect to both variables are

\( \displaystyle L'_{1,\xi}(\xi,\eta) = (-\sqrt{3})(1-\sqrt{3}\eta)/4, \quad L'_{1,\eta}(\xi,\eta) = (1-\sqrt{3}\xi)(-\sqrt{3})/4 \\ \displaystyle L'_{2,\xi}(\xi,\eta) = (+\sqrt{3})(1-\sqrt{3}\eta)/4, \quad L'_{1,\eta}(\xi,\eta) = (1-\sqrt{3}\xi)(-\sqrt{3})/4 \\ \displaystyle L'_{3,\xi}(\xi,\eta) = (+\sqrt{3})(1+\sqrt{3}\eta)/4, \quad L'_{1,\eta}(\xi,\eta) = (1-\sqrt{3}\xi)(-\sqrt{3})/4 \\ \displaystyle L'_{4,\xi}(\xi,\eta) = (-\sqrt{3})(1+\sqrt{3}\eta)/4, \quad L'_{1,\eta}(\xi,\eta) = (1-\sqrt{3}\xi)(-\sqrt{3})/4\)

The new shape function set will be termed quad_1_gauss_shape_set, and is introduced with a forward declaration:

// predefine the new shape set

class quad_1_gauss_shape_set;

Before defining the shape functions, we first define some basic properties of the shape function set by specialising the traits class template NiHu::shape_set_traits to the new shape set as follows:

namespace shape_set_traits
{
    // the shape set is defined over the quad domain
    template <>
    struct domain<quad_1_gauss_shape_set> : quad_domain {};
 
    // it has 4 shapeset nodes
    template <>
    struct num_nodes<quad_1_gauss_shape_set> { enum { value = 4 }; };
 
    // if used as geometrical shape set, Jacobian is 1-st order in xi/eta
    template <>
    struct jacobian_order<quad_1_gauss_shape_set> { enum { value = 1 }; };
 
    // N(xi,eta) is first order
    template <>
    struct polynomial_order<quad_1_gauss_shape_set> { enum { value = 1 }; };
 
    // N(xi, eta) and its derivatives are computed on the fly
    template <unsigned Order>
    struct shape_complexity<quad_1_gauss_shape_set, Order>
    {
        typedef matrix_function_complexity::general type;
    };
 
    // DOF locations are inside the element (2DOF)
    template <>
    struct position_dof_vector<quad_1_gauss_shape_set>
    {
        typedef tmp::vector<dof2, dof2, dof2, dof2> type;
    };
} // end of namespace shape_set_traits

We have defined the shape set's domain as the NiHu::quad_domain.
The number of shape set nodes is 4.
The polynomial order of the shape functions (the highest power of \(\xi\) or \(\eta\) in the definition of \(L_i(\xi,\eta)\)) is 1.
The polynomial order of the Jacobian (the highest power of \(\xi\) or \(\eta\) in the product of the derivatives \(L'_{\xi}\cdot L'_{\eta}\)) is 1 too.

We mention here that the order of the Jacobian is only needed if the shape function set is used as a geometrical interpolation function set. This will not be the case now, but we keep our code consistent with other shape function definitions.

After having the shape set traits defined, we can define the shape function class itself. The new shape function class must be derived from the CRTP base NiHu::shape_set_base. This base class defines a general interface for all shape sets, including convenient type definitions of the variable vector \(\xi\), shape function vector \(L_i(\xi)\) and its gradient matrix \(\nabla_{\bf \xi} L_i(\xi)\).

The new derived class must define three static member functions.

eval_shape evaluates the shape functions.
eval_dshape evaluates the gradient of the shape functions.
corner_begin returns a pointer to the first corner of the shape set.

// define the new shape set class
class quad_1_gauss_shape_set
    : public shape_set_base<quad_1_gauss_shape_set>
{
    typedef shape_set_base<quad_1_gauss_shape_set> base_t;
public:
 
    using base_t::corner_at;
    
    // return iterator to corners
    static xi_t const *corner_begin_impl(void)
    {
        return m_corners;
    }
 
protected:
    static const xi_t m_corners[num_nodes];
};

The functions returning the shape function set and its derivatives are defined as

// define expression for computing the zeroth order derivative N(xi, eta)
template <>
class shape_function<quad_1_gauss_shape_set, 0>
{
    typedef shape_set_traits::shape_value_type<quad_1_gauss_shape_set, 0>::type shape_t;
    typedef shape_set_traits::domain<quad_1_gauss_shape_set>::type::xi_t xi_t;
public:
    static shape_t eval(xi_t const &_xi)
    {
        auto xi = _xi[0], eta = _xi[1];
        return ( quad_1_gauss_shape_set::shape_t() <<
            (1.0 - std::sqrt(3.0)*xi) * (1.0 - std::sqrt(3.0)*eta),
            (1.0 + std::sqrt(3.0)*xi) * (1.0 - std::sqrt(3.0)*eta),
            (1.0 - std::sqrt(3.0)*xi) * (1.0 + std::sqrt(3.0)*eta),
            (1.0 + std::sqrt(3.0)*xi) * (1.0 + std::sqrt(3.0)*eta)
        ).finished() / 4.0;
    }
};
 
 
// define expression for computing the first order derivative N(xi, eta)
template <>
class shape_function<quad_1_gauss_shape_set, 1>
{
    typedef shape_set_traits::shape_value_type<quad_1_gauss_shape_set, 1>::type shape_t;
    typedef shape_set_traits::domain<quad_1_gauss_shape_set>::type::xi_t xi_t;
public:
    static shape_t eval(xi_t const &_xi)
    {
        auto xi = _xi[0], eta = _xi[1];
        return ( quad_1_gauss_shape_set::dshape_t() <<
            -std::sqrt(3.0) * (1.0 - std::sqrt(3.0)*eta) , (1.0 - std::sqrt(3.0)*xi) * -std::sqrt(3.0),
            +std::sqrt(3.0) * (1.0 - std::sqrt(3.0)*eta) , (1.0 + std::sqrt(3.0)*xi) * -std::sqrt(3.0),
            -std::sqrt(3.0) * (1.0 + std::sqrt(3.0)*eta) , (1.0 - std::sqrt(3.0)*xi) * +std::sqrt(3.0),
            +std::sqrt(3.0) * (1.0 + std::sqrt(3.0)*eta) , (1.0 + std::sqrt(3.0)*xi) * +std::sqrt(3.0)
        ).finished() / 4.0;
    }
};

The shape function's nodal locations are stored in the static array m_corners. The corner_begin function returns the address of the array.

// define the corners array
quad_1_gauss_shape_set::xi_t
    const quad_1_gauss_shape_set::m_corners[quad_1_gauss_shape_set::num_nodes] = {
        quad_1_gauss_shape_set::xi_t(-std::sqrt(3.0)/3.0, -std::sqrt(3.0)/3.0),
        quad_1_gauss_shape_set::xi_t(+std::sqrt(3.0)/3.0, -std::sqrt(3.0)/3.0),
        quad_1_gauss_shape_set::xi_t(-std::sqrt(3.0)/3.0, +std::sqrt(3.0)/3.0),
        quad_1_gauss_shape_set::xi_t(+std::sqrt(3.0)/3.0, +std::sqrt(3.0)/3.0)
};

That's all, we have defined the shape function set. From now on, it can be used for geometrical interpolation or field interpolation purposes. Furthermore, when using this shape function set in the collocational BEM context, the shape function nodes defined above are automatically used to generate weakly singular quadratures around the collocation points.

The field

The new field is going to be based on the standard NiHu::quad_1_elem element, extended with our new shape function set. The field will be termed quad_1_gauss_field, and is defined using a simple type definition:

// define the new field type
typedef field<
    quad_1_elem,            // over a linear quad elem
    quad_1_gauss_shape_set, // using the new shape set
    field_dimension::_1d    // scalar field
> quad_1_gauss_field;

Each field type is automatically assigned an integer identifier. However, we can override the field id definition by specialising the template structure (metafunction) NiHu::field_id to our new field type

namespace field_traits
{
    template <>
    struct id<quad_1_gauss_field> { enum {value = 666}; };
} // end of namespace field_traits

The defined field id will be used in the function space definition matrix to distinguish between different kind of fields.

Finally, a field tag is assigned to the field type using the metafunction NiHu::type2tag

// define a tag to the new type

typedef type2tag<quad_1_gauss_field> quad_1_gauss_field_tag;

Our new field type is ready to use in collocational, Galerkin or general BEM methods.

Table of Contents

Introduction

The Gaussian shape set

The field