field.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 NIHU_FIELD_HPP_INCLUDED
#define NIHU_FIELD_HPP_INCLUDED
 
#include "../util/crtp_base.hpp"
#include "corner_angle.hpp"
#include "element.hpp"
#include "field_dimension.hpp"
 
namespace NiHu
{
 
namespace field_option
{
struct isoparametric {};
struct constant {};
}
 
namespace field_traits
{
template <class Derived>
struct elem_type;
 
template <class Derived>
struct nset_type;
 
template <class Derived>
struct quantity_dimension;
 
template <class Derived>
struct is_dof_vector_stored : std::false_type {};
 
template <class Derived>
struct is_dirac : std::false_type {};
 
template <class Derived>
struct id
{
    enum {
        value =
        elem_type<Derived>::type::id * 100 + nset_type<Derived>::type::num_nodes
    };
};
 
template <class Derived>
struct dof_vector_type
{
    typedef Eigen::Matrix<
        unsigned,
        quantity_dimension<Derived>::value * nset_type<Derived>::type::num_nodes,
        1
    > type;
};
 
template <class Derived>
struct dof_vector_return_type : std::conditional<
    is_dof_vector_stored<Derived>::value,
    typename dof_vector_type<Derived>::type const &,
    typename dof_vector_type<Derived>::type
> {};
 
template <class Derived>
struct eval_return_type : shape_set_traits::shape_return_type<typename nset_type<Derived>::type, 0> {};
} // namespace NiHu::field_traits
 
 
template <class Derived>
class field_base
{
public:
    NIHU_CRTP_HELPERS
 
    typedef Derived type;
 
    typedef typename field_traits::elem_type<Derived>::type elem_t;
    typedef typename elem_t::x_t x_t;
    typedef typename elem_t::xi_t xi_t;
    typedef typename field_traits::nset_type<Derived>::type nset_t;
    typedef typename field_traits::dof_vector_type<Derived>::type dofs_t;
    typedef typename field_traits::dof_vector_return_type<Derived>::type dofs_return_t;
 
    typedef typename field_traits::eval_return_type<Derived>::type eval_return_t;
 
    enum {
        id = field_traits::id<Derived>::value,
        quantity_dimension = field_traits::quantity_dimension<Derived>::value,
        num_dofs = quantity_dimension * nset_t::num_nodes
    };
 
    eval_return_t eval(xi_t const &xi) const
    {
        return derived().eval(xi);
    }
 
    elem_t const &get_elem() const
    {
        return derived().get_elem();
    }
 
    dofs_return_t get_dofs() const
    {
        return derived().get_dofs();
    }
 
    double get_corner_angle(unsigned nset_corner_idx) const
    {
        return corner_angle_computer<elem_t, nset_t>::eval(get_elem(), nset_corner_idx);
    }
};
 
 
 
template <class Derived>
class field_impl;
 
 
// forward declaration
template <class Field>
class dirac_field;
 
namespace field_traits
{
template <class Derived>
struct elem_type<dirac_field<Derived> > : elem_type<Derived> {};
 
template <class Derived>
struct nset_type<dirac_field<Derived> > : nset_type<Derived> {};
 
template <class Derived>
struct quantity_dimension<dirac_field<Derived> > : quantity_dimension<Derived> {};
 
template <class Derived>
struct is_dirac<dirac_field<Derived> > : std::true_type {};
}
 
 
template <class Field>
class dirac_field :
    public field_base<dirac_field<Field> >,
    public field_impl<Field>
{
public:
    typedef dirac_field type;
 
    typedef field_impl<Field> impl_t;
 
    typedef typename field_traits::elem_type<dirac_field<Field> >::type elem_t;
    typedef Field field_t;
 
    typedef typename field_base<dirac_field<Field>>::nset_t nset_t;
 
    using impl_t::get_elem;
    using impl_t::get_dofs;
};
 
 
// forward declaration
template <class ElemType, class FieldOption, class Dimension = field_dimension::_1d>
class field_view;
 
namespace field_traits
{
template <class ElemType, class FieldOption, class Dimension>
struct elem_type<field_view<ElemType, FieldOption, Dimension> > : ElemType {};
 
template <class ElemType, class Dimension>
struct nset_type<field_view<ElemType, field_option::isoparametric, Dimension> >
{
    typedef  typename ElemType::lset_t type;
};
 
template <class ElemType, class Dimension>
struct nset_type<field_view<ElemType, field_option::constant, Dimension> >
{
    typedef constant_shape_set<typename ElemType::domain_t> type;
};
 
template <class ElemType, class FieldOption, class Dimension>
struct quantity_dimension<field_view<ElemType, FieldOption, Dimension> > : Dimension {};
}
 
 
template <class ElemType, class Dimension>
class field_impl<field_view<ElemType, field_option::isoparametric, Dimension> > :
    public ElemType
{
    typedef field_view<ElemType, field_option::isoparametric, Dimension> Derived;
    typedef typename field_traits::dof_vector_type<Derived>::type dofs_t;
    typedef typename field_traits::dof_vector_return_type<Derived>::type dofs_return_t;
    typedef typename field_traits::eval_return_type<Derived>::type eval_return_t;
    typedef typename field_traits::nset_type<Derived>::type nset_t;
    typedef typename ElemType::xi_t xi_t;
public:
 
    eval_return_t eval(xi_t const &xi) const
    {
        return nset_t::template eval_shape<0>(xi);
    }
 
 
    ElemType const &get_elem() const
    {
        return *this;   // static cast
    }
 
    dofs_return_t get_dofs() const
    {
        enum { D = Dimension::value };
        dofs_t dofs;
        for (unsigned n = 0; n < ElemType::num_nodes; ++n)
            for (unsigned d = 0; d < D; ++d)
                dofs(n*D + d) = ElemType::get_nodes()(n)*D + d;
        return dofs;
    }
};
 
 
 
 
template <class ElemType, class Dimension>
class field_impl<field_view< ElemType, field_option::constant, Dimension> > :
    public ElemType
{
private:
    typedef field_view<ElemType, field_option::constant, Dimension> Derived;
    typedef typename field_traits::dof_vector_type<Derived>::type dofs_t;
    typedef typename field_traits::dof_vector_return_type<Derived>::type dofs_return_t;
    typedef typename field_traits::eval_return_type<Derived>::type eval_return_t;
    typedef typename field_traits::nset_type<Derived>::type nset_t;
    typedef typename ElemType::xi_t xi_t;
public:
 
    eval_return_t eval(xi_t const &xi) const
    {
        return nset_t::template eval_shape<0>(xi);
    }
 
    ElemType const &get_elem() const
    {
        return *this; // static cast
    }
 
    dofs_return_t get_dofs() const
    {
        enum { D = Dimension::value };
        return dofs_t::LinSpaced(D, ElemType::get_id()(0)*D, ElemType::get_id()(0)*D + D - 1);
    }
};
 
 
template <class ElemType, class Option, class Dimension>
class field_view :
    public field_base<field_view<ElemType, Option, Dimension> >,
    public field_impl<field_view<ElemType, Option, Dimension> >
{
    typedef field_base<field_view<ElemType, Option, Dimension> > crtp_base_t;
    typedef field_impl<field_view<ElemType, Option, Dimension> > impl_t;
 
public:
    typedef field_view type;
 
    typedef typename crtp_base_t::elem_t elem_t;
 
    enum {
        id = crtp_base_t::id
    };
 
    typedef typename crtp_base_t::nset_t nset_t;
 
    using impl_t::eval;
    using impl_t::get_elem;
    using impl_t::get_dofs;
};
 
 
template <class ElemType, class NSet, class Dimension>
class field;
 
 
namespace field_traits
{
template <class ElemType, class NSet, class Dimension>
struct elem_type<field<ElemType, NSet, Dimension> > : ElemType {};
 
template <class ElemType, class NSet, class Dimension>
struct nset_type<field<ElemType, NSet, Dimension> >
{
    typedef NSet type;
};
 
template <class ElemType, class NSet, class Dimension>
struct quantity_dimension<field<ElemType, NSet, Dimension> > : Dimension {};
} // end of namespace field_traits
 
 
template <class ElemType, class NSet, class Dimension>
class field_impl<field<ElemType, NSet, Dimension> >
{
public:
    typedef field<ElemType, NSet, Dimension> Derived;
    typedef ElemType elem_t;
    typedef typename ElemType::xi_t xi_t;
    typedef typename field_traits::dof_vector_type<Derived>::type dofs_t;
    typedef typename field_traits::dof_vector_return_type<Derived>::type dofs_return_t;
 
    typedef typename field_traits::eval_return_type<Derived>::type eval_return_t;
 
    field_impl(element_base<elem_t> const &elem, dofs_t const &dofs) :
        m_elem(elem.derived()), m_dofs(dofs)
    {
    }
 
    eval_return_t eval(xi_t const &xi) const
    {
        return NSet::template eval_shape<0>(xi);
    }
 
    elem_t const &get_elem() const
    {
        return m_elem;
    }
 
    dofs_return_t get_dofs() const
    {
        return m_dofs;
    }
 
protected:
    elem_t m_elem;
    dofs_t m_dofs;
};
 
template <class ElemType, class NSet, class Dimension = field_dimension::_1d>
class field :
    public field_base<field<ElemType, NSet, Dimension> >,
    public field_impl<field<ElemType, NSet, Dimension> >
{
public:
    typedef field type;
    typedef field_base<field<ElemType, NSet, Dimension> > base_t;
    typedef field_impl<field<ElemType, NSet, Dimension> > impl_t;
 
    typedef typename base_t::elem_t elem_t;
    typedef typename base_t::dofs_t dofs_t;
 
    using impl_t::eval;
    using impl_t::get_elem;
    using impl_t::get_dofs;
 
    field(element_base<elem_t> const &elem, dofs_t const &dofs) :
        impl_t(elem, dofs)
    {
    }
};
 
template <class Elem, class Option, class Dimension = field_dimension::_1d>
field_view<Elem, Option, Dimension> const &
create_field_view(element_base<Elem> const & e, Option, Dimension dim = Dimension())
{
    return static_cast<field_view<Elem, Option, Dimension> const &>(e.derived());
}
 
 
 
 
template <class Field>
class directional_derivative_field;
 
 
namespace field_traits
{
template <class Field>
struct elem_type<directional_derivative_field<Field>> : elem_type<Field> {};
 
template <class Field>
struct nset_type<directional_derivative_field<Field>> : nset_type<Field> {};
 
template <class Field>
struct quantity_dimension<directional_derivative_field<Field>> : quantity_dimension<Field> {};
 
template <class Field>
struct eval_return_type<directional_derivative_field<Field>>
    : shape_set_traits::shape_value_type<typename nset_type<directional_derivative_field<Field>>::type, 0> {};
 
} // end of namespace field_traits
 
 
template <class Field>
class directional_derivative_field :
    public field_base<directional_derivative_field<Field>>
{
public:
    typedef directional_derivative_field type;
 
    typedef field_base<directional_derivative_field<Field>> crtp_base_t;
 
    typedef typename crtp_base_t::nset_t nset_t;
    typedef typename crtp_base_t::elem_t elem_t;
    typedef typename crtp_base_t::xi_t xi_t;
    typedef typename crtp_base_t::x_t x_t;
    typedef typename crtp_base_t::dofs_return_t dofs_return_t;
    typedef typename crtp_base_t::eval_return_t eval_return_t;
 
    template <class Direction>
    directional_derivative_field(Field const &field, Eigen::MatrixBase<Direction> const &direction)
        : m_field(field)
        , m_direction(direction.derived())
    {
    }
 
    eval_return_t eval(xi_t const &xi) const
    {
        auto dx = get_elem().get_dx(xi);
        double dx_xi = dx.col(0).norm();        // dx / dxi
        double dx_eta = dx.col(1).norm();       // dx / deta
        x_t s = dx.col(0).normalized();         // tangential unit vector
        x_t t_eta = dx.col(1).normalized();
        x_t n = s.cross(t_eta).normalized();    // unit normal
        x_t t = n.cross(s);                     // second tangential unit vector
 
        x_t const &m = m_direction;
        double a = (s.dot(m) - s.dot(t_eta) / t.dot(t_eta) * t.dot(m)) / dx_xi;
        double b = t.dot(m) / t.dot(t_eta) / dx_eta;
 
        auto dN = nset_t::template eval_shape<1>(xi);
        return a * dN.col(0) + b * dN.col(1);
    }
 
    elem_t const &get_elem() const
    {
        return m_field.get_elem();
    }
 
    dofs_return_t get_dofs() const
    {
        return m_field.get_dofs();
    }
 
private:
    Field const &m_field;
    x_t const &m_direction;
};
 
template <class Field, class Direction>
directional_derivative_field<Field>
directional_derivative(field_base<Field> const &field, Eigen::MatrixBase<Direction> const &direction)
{
    return directional_derivative_field<Field>(field.derived(), direction);
}
 
 
template <class Elem, class Dimension = field_dimension::_1d>
field_view<Elem, field_option::constant, Dimension> const &
constant_view(element_base<Elem> const & e, Dimension dim = Dimension())
{
    return create_field_view(e, field_option::constant(), dim);
}
 
template <class Elem, class Dimension = field_dimension::_1d>
field_view<Elem, field_option::isoparametric, Dimension> const &
isoparametric_view(element_base<Elem> const & e, Dimension dim = Dimension())
{
    return create_field_view(e, field_option::isoparametric(), dim);
}
 
template <class Field>
dirac_field<Field> const &
dirac(field_base<Field> const & f)
{
    return static_cast<dirac_field<Field> const &>(
        static_cast<field_impl<Field> const &>(f.derived()));
}
 
} // end of namespace NiHu
 
#endif // NIHU_FIELD_HPP_INCLUDED