shapeset.hpp

Go to the documentation of this file.

// This file is a part of NiHu, a C++ BEM template library.
//
// Copyright (C) 2012-2019  Peter Fiala <fiala@hit.bme.hu>
// Copyright (C) 2012-2019  Peter Rucz <rucz@hit.bme.hu>
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program.  If not, see <http://www.gnu.org/licenses/>.
 
#ifndef SHAPESET_HPP_INCLUDED
#define SHAPESET_HPP_INCLUDED
 
#include "../tmp/algorithm.hpp"
#include "../tmp/vector.hpp"
#include "../util/conditional_precompute.hpp"
#include "../util/materialize_sequence.hpp"
 
#include "domain.hpp"
 
#include <iostream>
#include <stdexcept>
#include <type_traits>
 
namespace NiHu
{
 
namespace shape_derivative_index
{
    enum {
        dXI = 0,    
        dETA = 1,   
        dXIXI = 0,  
        dXIETA = 1, 
        dETAXI = 1, 
        dETAETA = 2 
    };
}
 
template <unsigned d>
struct position_dof : std::integral_constant<unsigned, d> {};
typedef position_dof<0> dof0;
typedef position_dof<1> dof1;
typedef position_dof<2> dof2;
typedef position_dof<3> dof3;
 
template <int index>
struct corner_index : std::integral_constant<int, index> {};
 
typedef corner_index<-1> not_corner;
 
// forward declaration
template <class Derived, unsigned Order>
class shape_function;
 
constexpr unsigned num_derivatives(unsigned order, unsigned dim)
{
    return order == 0 ? 1 : (order == 1 ? dim : dim * (dim + 1) / 2);
}
 
namespace shape_set_traits
{
    template <class Derived>
    struct name
    {
        static const std::string value; 
    };
 
    template <class Derived>
    struct domain;
 
    template <class Derived>
    struct num_nodes;
 
    template <class Derived>
    struct id
    {
        enum {
            value = domain_traits::id<typename domain<Derived>::type>::value * 100 +
            num_nodes<Derived>::value
        };
    };
 
    template <class Derived>
    struct polynomial_order;
 
    template <class Derived>
    struct jacobian_order;
 
    template <class Derived, unsigned Order>
    struct shape_complexity;
 
    template <class Derived, unsigned Order>
    struct shape_value_type
    {
        typedef Eigen::Matrix<
            typename domain<Derived>::type::scalar_t,
            num_nodes<Derived>::value,
            num_derivatives(Order, domain<Derived>::type::dimension)
        > type;
    };
 
    template <class Derived, unsigned Order>
    struct factory_functor : conditional_precompute<
        typename shape_complexity<Derived, Order>::type,
        shape_function<Derived, Order>,
        typename domain<Derived>::type::xi_t
    > {};
 
    template <class Derived, unsigned Order>
    struct shape_return_type
    {
        typedef typename factory_functor<Derived, Order>::type::return_type type;
    };
 
    template <class Derived>
    struct position_dof_vector;
 
    template <class Derived>
    struct position_dof_vector_mat
    {
        typedef typename sequence_materializer<
            typename position_dof_vector<Derived>::type
        >::eigen_vector_t type;
    };
 
    template <class Derived>
    struct corner_index_vector;
 
    template <class Derived>
    struct corner_index_vector_mat
    {
        typedef typename sequence_materializer<
            typename corner_index_vector<Derived>::type
        >::eigen_vector_t type;
    };
 
}
 
 
template <class Derived>
class shape_set_base
{
public:
    typedef typename shape_set_traits::domain<Derived>::type domain_t;
 
    enum {
        num_nodes = shape_set_traits::num_nodes<Derived>::value,
        polynomial_order = shape_set_traits::polynomial_order<Derived>::value,
        jacobian_order = shape_set_traits::jacobian_order<Derived>::value,
        id = shape_set_traits::id<Derived>::value,
    };
 
    typedef typename domain_t::scalar_t scalar_t;
    typedef typename domain_t::xi_t xi_t;
 
    typedef std::array<xi_t, num_nodes> corners_t;
    typedef typename corners_t::const_iterator corner_iterator_t;
 
    template <unsigned Order>
    struct shape_value_type : shape_set_traits::shape_value_type<Derived, Order> {};
 
    typedef typename shape_value_type<0>::type shape_t;
    typedef typename shape_value_type<1>::type dshape_t;
    typedef typename shape_value_type<2>::type ddshape_t;
 
    template <unsigned Order>
    static typename shape_set_traits::shape_return_type<Derived, Order>::type
    eval_shape(xi_t const &xi)
    {
        return shape_set_traits::factory_functor<Derived, Order>::type::eval(xi);
    }
 
    typedef typename shape_set_traits::position_dof_vector<Derived>::type position_dof_vector;
    typedef typename shape_set_traits::corner_index_vector<Derived>::type corner_index_vector;
 
    static corner_iterator_t corner_begin()
    {
        return Derived::corner_begin_impl();
    }
 
    static corner_iterator_t corner_end()
    {
        return Derived::corner_begin() + num_nodes;
    }
 
    static xi_t const &corner_at(size_t idx)
    {
        return *(corner_begin() + idx);
    }
 
    static unsigned node_to_domain_corner(unsigned idx)
    {
        typedef typename shape_set_traits::corner_index_vector_mat<Derived>::type runtime_corner_index_vector_t;
        static runtime_corner_index_vector_t runtime_corner_index_vector = materialize_sequence(corner_index_vector());
        int index = runtime_corner_index_vector[idx];
        if (index < 0)
            throw std::runtime_error("Node " + std::to_string(idx) + " is not a corner node in " + shape_set_traits::name<Derived>::value);
        return index;
    }
 
 
    static unsigned position_dof(unsigned idx)
    {
        typedef typename shape_set_traits::position_dof_vector_mat<Derived>::type runtime_position_dof_vector_t;
        static runtime_position_dof_vector_t runtime_position_dof_vector = materialize_sequence(position_dof_vector());
        return runtime_position_dof_vector[idx];
    }
 
};
 
// Forward declaration
template <class Domain>
class constant_shape_set;
 
namespace shape_set_traits
{
    template <class Domain>
    struct domain<constant_shape_set<Domain> > : Domain {};
 
    template <class Domain>
    struct num_nodes<constant_shape_set<Domain> >
    {
        enum { value = 1 };
    };
 
    template <class Domain>
    struct polynomial_order<constant_shape_set<Domain> >
    {
        enum { value = 0 };
    };
 
    template <class Domain>
    struct jacobian_order<constant_shape_set<Domain> >
    {
        enum { value = 0 };
    };
 
    template <class Domain>
    struct shape_complexity<constant_shape_set<Domain>, 0> : matrix_function_complexity::constant {};
 
    template <class Domain>
    struct shape_complexity<constant_shape_set<Domain>, 1> : matrix_function_complexity::zero {};
 
    template <class Domain>
    struct shape_complexity<constant_shape_set<Domain>, 2> : matrix_function_complexity::zero {};
 
    template <class Domain>
    struct position_dof_vector<constant_shape_set<Domain> >
    {
        typedef tmp::vector<position_dof<Domain::dimension> > type;
    };
 
    template <class Domain>
    struct corner_index_vector<constant_shape_set<Domain>>
    {
        typedef tmp::vector<not_corner>  type;
    };
 
}
 
template <class Domain>
class constant_shape_set : public shape_set_base<constant_shape_set<Domain> >
{
public:
    typedef shape_set_base<constant_shape_set<Domain> > base_t;
    typedef typename base_t::domain_t domain_t;
    typedef typename base_t::xi_t xi_t;
    typedef typename base_t::corners_t corners_t;
    typedef typename base_t::corner_iterator_t corner_iterator_t;
 
    static corner_iterator_t corner_begin_impl()
    {
        // Note: function-local static is used here in order to avoid
        //       ambigous order of initialization of Domain's and ShapeSet's
        //       static members.
        static corners_t const m_corners = { Domain::get_center() };
 
        return m_corners.cbegin();
    }
};
 
 
template <class Domain>
class shape_function<constant_shape_set<Domain>, 0>
{
    typedef typename shape_set_traits::shape_value_type<constant_shape_set<Domain>, 0>::type shape_t;
    typedef typename shape_set_traits::domain<constant_shape_set<Domain> >::type::xi_t xi_t;
public:
    static shape_t eval(xi_t const &)
    {
        return shape_t::Ones();
    }
};
 
template <class Domain>
class shape_function<constant_shape_set<Domain>, 1>
{
    typedef typename shape_set_traits::shape_value_type<constant_shape_set<Domain>, 1>::type shape_t;
    typedef typename shape_set_traits::domain<constant_shape_set<Domain>>::type::xi_t xi_t;
public:
    static shape_t eval(xi_t const &xi)
    {
        return shape_t::Zero();
    }
};
 
template <class Domain>
class shape_function<constant_shape_set<Domain>, 2>
{
    typedef typename shape_set_traits::shape_value_type<constant_shape_set<Domain>, 2>::type shape_t;
    typedef typename shape_set_traits::domain<constant_shape_set<Domain>>::type::xi_t xi_t;
public:
    static shape_t eval(xi_t const &xi)
    {
        return shape_t::Zero();
    }
};
 
 
// Forward declaration
template <class Domain>
class isoparam_shape_set;
 
namespace shape_set_traits
{
    template <class Domain>
    struct domain<isoparam_shape_set<Domain> > : Domain {};
 
    template <class Domain>
    struct num_nodes<isoparam_shape_set<Domain> >
    {
        enum { value = Domain::num_corners };
    };
 
    template <class Domain>
    struct polynomial_order<isoparam_shape_set<Domain> >
    {
        enum { value = 1 };
    };
 
    template <class Domain>
    struct position_dof_vector<isoparam_shape_set<Domain> > : tmp::constant_sequence<
        dof0,
        num_nodes<isoparam_shape_set<Domain> >::value,
        tmp::vector<>
    > {};
 
    template <class Domain>
    struct corner_index_vector<isoparam_shape_set<Domain>> : tmp::arithmetic_sequence<
        std::integral_constant<int, 0>,
        num_nodes<isoparam_shape_set<Domain> >::value,
        tmp::vector<>,
        1
    > {};
}
 
template <class Domain>
class isoparam_shape_set
    : public shape_set_base<isoparam_shape_set<Domain> >
{
public:
    typedef shape_set_base<isoparam_shape_set<Domain> > base_t;
    typedef typename base_t::domain_t domain_t;
 
    typedef typename base_t::xi_t xi_t;
 
    typedef typename base_t::corners_t corners_t;
    typedef typename base_t::corner_iterator_t corner_iterator_t;
 
    static corner_iterator_t corner_begin_impl()
    {
        return domain_t::get_corners().cbegin();
    }
};
 
} // end of namespace NiHu
 
#endif // SHAPESET_HPP_INCLUDED