element.hpp

Go to the documentation of this file.

// This file is a part of NiHu, a C++ BEM template library.
//
// Copyright (C) 2012-2014  Peter Fiala <fiala@hit.bme.hu>
// Copyright (C) 2012-2014  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 ELEMENT_HPP_INCLUDED
#define ELEMENT_HPP_INCLUDED
 
#include <boost/math/constants/constants.hpp>
 
#include "../util/crtp_base.hpp"
#include "../tmp/bool.hpp"
#include "shapeset.hpp"
 
#include <iostream>
 
namespace NiHu
{
 
class element_overlapping
{
public:
    unsigned get_num(void) const
    {
        return m_num;
    }
 
    unsigned get_ind1(void) const
    {
        return m_ind1;
    }
 
    unsigned get_ind2(void) const
    {
        return m_ind2;
    }
 
    element_overlapping(
        unsigned num = 0, unsigned ind1 = 0, unsigned ind2 = 0) :
        m_num(num), m_ind1(ind1), m_ind2(ind2)
    {
    }
    
    std::ostream &print_debug(std::ostream &os) const
    {
        return os << "Element overlapping info: " << std::endl
            << "Number of coincident nodes: " << m_num << std::endl
            << "Start index on element 1 : " << m_ind1 << " and on element 2 " << m_ind2;
    }
 
private:
    unsigned m_num;
    unsigned m_ind1, m_ind2;
};
 
 
template <class Derived, unsigned Order>
class location_impl;
 
namespace element_traits
{
    template <class Derived>
    struct name
    {
        static const std::string value;
    };
 
    template <class Derived>
    struct space_type;
 
    template <class Derived>
    struct lset;
 
    template <class Derived>
    struct is_surface_element;
 
    template <class Derived>
    struct id
    {
        enum {
            value = element_traits::space_type<Derived>::type::dimension * 10000 +
            shape_set_traits::id<typename element_traits::lset<Derived>::type>::value
        };
    };
 
    template <class Derived, unsigned Order>
    struct location_complexity : shape_set_traits::shape_complexity<
        typename lset<Derived>::type,
        Order
    > {};
 
    template <class Derived>
    struct coords_type
    {
        typedef Eigen::Matrix<
            typename space_type<Derived>::type::scalar_t,
            space_type<Derived>::type::dimension,
            lset<Derived>::type::num_nodes
        > type;
    };
 
    template <class Derived, unsigned Order>
    struct location_value_type
    {
        typedef Eigen::Matrix<
            typename space_type<Derived>::type::scalar_t,
            space_type<Derived>::type::dimension,
            num_derivatives(Order, shape_set_traits::domain<typename lset<Derived>::type>::type::dimension)
        > type;
    };
 
    template <class Derived, unsigned Order>
    struct location_factory_functor : conditional_precompute_instance<
        typename location_complexity<Derived, Order>::type,
        location_impl<Derived, Order>,
        typename coords_type<Derived>::type,
        typename shape_set_traits::domain<
            typename lset<Derived>::type
        >::type::xi_t
    > {};
 
    template <class Derived, unsigned Order>
    struct location_return_type
    {
        typedef typename location_factory_functor<
            Derived,
            Order
        >::type::return_type type;
    };
}
 
typedef Eigen::Matrix<unsigned, 1, 1> elem_id_t;
 
 
template <class Derived, unsigned Order>
class location_impl
{
public:
    typedef typename element_traits::lset<Derived>::type lset_t;
    typedef typename element_traits::coords_type<Derived>::type coords_t;
    typedef typename element_traits::location_value_type<Derived, Order>::type ret_t;
    typedef typename shape_set_traits::domain<lset_t>::type::xi_t xi_t;
 
    static ret_t eval(coords_t const &coords, xi_t const &xi)
    {
        return coords * lset_t::template eval_shape<Order>(xi);
    }
};
 
 
template <class Derived>
class element_base
{
public:
    typedef Derived type;
 
    typedef typename element_traits::space_type<Derived>::type space_t;
    typedef typename element_traits::lset<Derived>::type lset_t;
    typedef typename space_t::scalar_t scalar_t;
 
    typedef typename lset_t::domain_t domain_t;
    typedef typename domain_t::xi_t xi_t;
 
    enum {
        id = element_traits::id<Derived>::value,
        num_nodes = lset_t::num_nodes,
    };
 
    typedef typename element_traits::location_value_type<Derived, 0>::type x_t;
    typedef typename element_traits::location_return_type<Derived, 0>::type x_return_type;
    typedef typename element_traits::location_value_type<Derived, 1>::type dx_t;
    typedef typename element_traits::location_return_type<Derived, 1>::type dx_return_type;
    typedef typename element_traits::location_value_type<Derived, 2>::type ddx_t;
    typedef typename element_traits::location_return_type<Derived, 2>::type ddx_return_type;
 
    typedef Eigen::Matrix<unsigned, num_nodes, 1> nodes_t;
    typedef typename element_traits::coords_type<Derived>::type coords_t;
    typedef elem_id_t id_t;
 
protected:
    id_t m_id;
    nodes_t m_nodes;
    coords_t m_coords;
    x_t m_center;
    scalar_t m_linear_size_estimate;
 
    typename element_traits::location_factory_functor<Derived, 0>::type x_computer;
    typename element_traits::location_factory_functor<Derived, 1>::type dx_computer;
    typename element_traits::location_factory_functor<Derived, 2>::type ddx_computer;
 
public:
    NIHU_CRTP_HELPERS
 
    element_base() {}
 
    element_base(coords_t const &coords, unsigned id = 0, nodes_t const &nodes = nodes_t())
        :m_id((id_t()<<id).finished())
        , m_nodes(nodes)
        , m_coords(coords)
        , m_center(get_x(domain_t::get_center()))
        , x_computer(m_coords, xi_t())
        , dx_computer(m_coords, xi_t())
        , ddx_computer(m_coords, xi_t())
    {
        // compute linear size estimate (average edge length)
        scalar_t d = 0.0;
        for (unsigned e = 0; e < domain_t::num_edges; ++e)
        {
            auto const &edge = domain_t::get_edge(e);
            x_t x0 = get_x(domain_t::get_corner(edge[0]));
            x_t x1 = get_x(domain_t::get_corner(edge[1]));
            d += (x1-x0).norm();
        }
        d /= scalar_t(domain_t::num_edges);
        m_linear_size_estimate = d;
    }
 
    id_t const &get_id(void) const
    {
        return m_id;
    }
 
    nodes_t const &get_nodes(void) const
    {
        return m_nodes;
    }
 
    coords_t const &get_coords(void) const
    {
        return m_coords;
    }
 
    x_return_type get_x(xi_t const &xi) const
    {
        return x_computer(m_coords, xi);
    }
 
    dx_return_type get_dx(xi_t const &xi) const
    {
        return dx_computer(m_coords, xi);
    }
 
    ddx_return_type get_ddx(xi_t const &xi) const
    {
        return ddx_computer(m_coords, xi);
    }
 
    x_t const &get_center(void) const
    {
        return m_center;
    }
 
    scalar_t const &get_linear_size_estimate(void) const
    {
        return m_linear_size_estimate;
    }
    
    template <class OtherElem>
    element_overlapping get_overlapping(OtherElem const &other) const
    {
        unsigned num_coinc = 0; // number of coincident nodes
        unsigned start_ind1 = 0, start_ind2 = 0;
 
        auto const &otherNodes = other.get_nodes();
        for (unsigned i = 0; i < num_nodes; ++i)
        {
            for (unsigned j = 0; j < OtherElem::num_nodes; ++j)
            {
                if (m_nodes[i] == otherNodes[j])
                {
                    if (num_coinc == 0)
                    {
                        start_ind1 = i;
                        start_ind2 = j;
                    }
                    else
                    {
                        if (i < start_ind1
                            || (i==num_nodes-1
                            && start_ind1 == 0
                            && num_coinc == 1))
                            start_ind1 = i;
                        if (j < start_ind2 ||
                            (j==OtherElem::num_nodes-1 &&
                            start_ind2 == 0 &&
                            num_coinc == 1))
                            start_ind2 = j;
                    }
                    num_coinc++;
                }
            }
        }
        if (num_coinc == 0)
            return element_overlapping();
        return element_overlapping(num_coinc, start_ind1, start_ind2);
    }
    
    Derived get_linearized_elem(xi_t const &xi0) const
    {
        // coordinates of the linearized elem
        coords_t coords;
 
        // phyisical location and its derivative at the reference point
        x_t x0 = get_x(xi0);
        auto dx0 = get_dx(xi0);
 
        // compute corners of linearized elem
        for (size_t i = 0; i < num_nodes; ++i)
        {
            xi_t const &xi = lset_t::corner_at(i);
            xi_t dxi = xi - xi0;
            coords.col(i) = x0 + dx0 * dxi;
        }
 
        // construct linearized elem
        return Derived(coords);
    }
};
 
inline bool operator<(elem_id_t const &lhs, elem_id_t const &rhs)
{
    return lhs(0,0) < rhs(0,0);
}
 
 
// forward declaration
template <class LSet, class scalar_t>
class surface_element;
 
// forward declaration
template <class LSet, class scalar_t>
class volume_element;
 
 
namespace element_traits
{
    /* the space dimension is computed from the Lset domain's dimension */
    template <class LSet, class Scalar>
    struct space_type<volume_element<LSet, Scalar> >
        : space<Scalar, LSet::domain_t::dimension> {};
 
    template <class LSet, class Scalar>
    struct lset<volume_element<LSet, Scalar> >
    {
        typedef LSet type;
    };
 
    template <class LSet, class Scalar>
    struct is_surface_element<volume_element<LSet, Scalar> > : std::false_type {};
}
 
 
 
template<class Derived, class enable = void>
class normal_impl;
 
namespace element_traits
{
    /* the space dimension is computed from the Lset domain's dimension (dL+1) */
    template <class LSet, class Scalar>
    struct space_type<surface_element<LSet, Scalar> >
        : space<Scalar, LSet::domain_t::dimension + 1> {};
 
    template <class LSet, class Scalar>
    struct lset<surface_element<LSet, Scalar> >
    {
        typedef LSet type;
    };
 
    template <class LSet, class Scalar>
    struct is_surface_element<surface_element<LSet, Scalar> > : std::true_type {};
}
 
template <class Derived>
class normal_impl<Derived, typename std::enable_if<
    element_traits::space_type<Derived>::type::dimension == 3
>::type>
{
    typedef typename element_traits::location_value_type<Derived, 0>::type x_t;
    typedef typename element_traits::location_value_type<Derived, 1>::type dx_t;
public:
    static x_t eval(dx_t const &m)
    {
        return m.col(shape_derivative_index::dXI).cross(m.col(shape_derivative_index::dETA));
    }
};
 
template <class Derived>
class normal_impl<Derived, typename std::enable_if<
    element_traits::space_type<Derived>::type::dimension == 2
>::type>
{
    typedef typename element_traits::location_value_type<Derived, 0>::type x_t;
    typedef typename element_traits::location_value_type<Derived, 1>::type dx_t;
public:
    static x_t eval(dx_t const &m)
    {
        using namespace boost::math::double_constants;
        return Eigen::Rotation2D<typename Derived::scalar_t>(-pi/2.0) * m;
    }
};
 
namespace element_traits
{
    template <class Derived>
    struct normal_factory_functor : conditional_precompute_instance<
        typename location_complexity<Derived, 1>::type,
        normal_impl<Derived>,
        typename location_value_type<Derived, 1>::type
    > {};
 
    template <class Derived>
    struct normal_return_type
    {
        typedef typename normal_factory_functor<
            Derived
        >::type::return_type type;
    };
}
 
template <class LSet, class scalar_t>
class surface_element :
    public element_base<surface_element<LSet, scalar_t> >
{
public:
    typedef element_base<surface_element<LSet, scalar_t> > base_t;
 
    using typename base_t::coords_t;
    using typename base_t::nodes_t;
    using typename base_t::xi_t;
    using typename base_t::x_t;
    using typename base_t::dx_t;
 
    using domain_t = typename base_t::domain_t;
 
    using base_t::get_dx;
 
    typedef typename element_traits::normal_return_type<surface_element>::type normal_return_type;
 
protected:
    typename element_traits::normal_factory_functor<surface_element>::type normal_computer;
 
public:
    surface_element(
        coords_t const &coords,
        unsigned id = 0,
        nodes_t const &nodes = nodes_t())
        : base_t(coords, id, nodes),
        normal_computer(get_dx(xi_t()))
    {
    }
    
    void flip_inplace(void)
    {
        *this = flip();
    }
    
    surface_element flip(void) const
    {
        return surface_element(this->get_coords().colwise().reverse(), this->get_id()(0), this->get_nodes().reverse());
    }
 
    normal_return_type get_normal(xi_t const &xi = domain_t::get_center()) const
    {
        return normal_computer(get_dx(xi));
    }
};
 
 
 
template <class LSet, class scalar_t>
class volume_element :
    public element_base<volume_element<LSet, scalar_t> >
{
public:
    typedef element_base<volume_element<LSet, scalar_t> > base_t;
 
    using typename base_t::coords_t;
    using typename base_t::nodes_t;
    using typename base_t::xi_t;
    using typename base_t::x_t;
    using typename base_t::dx_t;
 
    using domain_t = typename base_t::domain_t;
 
    using base_t::get_dx;
 
    volume_element(
        coords_t const &coords,
        unsigned id = 0,
        nodes_t const &nodes = nodes_t())
        : base_t(coords, id, nodes)
    {
    }
};
 
} // namespace NiHu
 
 
#endif