function_space.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 FUNCTION_SPACE_HPP_INCLUDED
#define FUNCTION_SPACE_HPP_INCLUDED
 
#include "../util/casted_iterator.hpp"
#include "../util/crtp_base.hpp"
#include "../util/type2tag.hpp"
#include "field.hpp"
#include "field_dimension.hpp"
#include "mesh.hpp"
 
namespace NiHu
{
 
template <class Derived>
struct function_space_traits;
 
 
template <class Derived>
class function_space_base
{
public:
    NIHU_CRTP_HELPERS
 
    typedef function_space_traits<Derived> traits_t;
 
    typedef typename traits_t::mesh_t mesh_t;
 
    typedef typename traits_t::field_type_vector_t field_type_vector_t;
 
    typedef typename tmp::deref<
        typename tmp::begin<
            typename mesh_t::elem_type_vector_t
        >::type
    >::type::x_t x_t;
 
    enum {
        quantity_dimension = traits_t::quantity_dimension
    };
 
    mesh_t const &get_mesh() const
    {
        return derived().get_mesh();
    }
 
    template <class FieldType>
    typename traits_t::template iterator<FieldType>::type
        field_begin() const
    {
        return derived().template field_begin<FieldType>();
    }
 
    template <class FieldType>
    typename traits_t::template iterator<FieldType>::type
        field_end() const
    {
        return derived().template field_end<FieldType>();
    }
 
    size_t get_num_dofs() const
    {
        return derived().get_num_dofs();
    }
};
 
 
template <class Derived>
class function_space_impl;
 
 
// forward declaration
template<class Mesh, class FieldOption, class Dimension>
class function_space_view;
 
 
template <class Mesh, class FieldOption, class Dimension>
struct function_space_traits<function_space_view<Mesh, FieldOption, Dimension> >
{
    typedef Mesh mesh_t;
 
    enum {
        quantity_dimension = Dimension::value
    };
 
    typedef typename mesh_t::elem_type_vector_t etv_t;
 
    typedef typename tmp::transform<
        etv_t,
        tmp::inserter <
            typename tmp::empty<etv_t>::type,
            tmp::push_back<tmp::_1, tmp::_2>
        >,
        field_view<tmp::_1, FieldOption, Dimension>
    >::type field_type_vector_t;
 
    template <class field_t>
    struct iterator : casted_iterator<
        typename mesh_elem_iterator_t<typename field_t::elem_t>::type,
        field_t
    > {};
};
 
 
namespace function_space_impl_internal
{
    template <class Mesh, class FieldOption, class Dimension>
    struct get_num_dofs;
    
    template <class Mesh, class Dimension>
    struct get_num_dofs<Mesh, field_option::constant, Dimension>
    {
        typedef function_space_traits< function_space_view<Mesh, field_option::constant, Dimension> > traits_t;
        static size_t eval(Mesh const& mesh)
        {
            return mesh.get_num_elements() * traits_t::quantity_dimension;
        }
    };
 
    template <class Mesh, class Dimension>
    struct get_num_dofs<Mesh, field_option::isoparametric, Dimension>
    {
        typedef function_space_traits< function_space_view<Mesh, field_option::isoparametric, Dimension> > traits_t;
        static size_t eval(Mesh const& mesh)
        {
            return mesh.get_num_points() * traits_t::quantity_dimension;;
        }
    };
 
} // namespace NiHu::function_space_impl_internal
 
template<class Mesh, class FieldOption, class Dimension>
class function_space_impl<function_space_view<Mesh, FieldOption, Dimension> > :
    public Mesh
{
public:
    typedef function_space_traits< function_space_view<Mesh, FieldOption, Dimension> > traits_t;
 
#if 0
private:
    size_t get_num_dofs_impl(field_option::constant) const
    {
        return Mesh::get_num_elements() * traits_t::quantity_dimension;
    }
 
    size_t get_num_dofs_impl(field_option::isoparametric) const
    {
        return Mesh::get_num_points() * traits_t::quantity_dimension;
    }
#endif
 
public:
    Mesh const &get_mesh() const
    {
        return static_cast<Mesh const &> (*this);
    }
 
    template <class FieldType>
    typename traits_t::template iterator<FieldType>::type
        field_begin() const
    {
        return Mesh::template begin<typename FieldType::elem_t>();
    }
 
    template <class FieldType>
    typename traits_t::template iterator<FieldType>::type
        field_end() const
    {
        return Mesh::template end<typename FieldType::elem_t>();
    }
 
    size_t get_num_dofs() const
    {
        return function_space_impl_internal::get_num_dofs<Mesh, FieldOption, Dimension>::eval(*this);
        //return get_num_dofs_impl(FieldOption());
    }
};
 
 
template<class Mesh, class FieldOption, class Dimension = field_dimension::_1d>
class function_space_view :
    public function_space_base<function_space_view<Mesh, FieldOption, Dimension> >,
    public function_space_impl<function_space_view<Mesh, FieldOption, Dimension> >
{
public:
    typedef function_space_impl<function_space_view<Mesh, FieldOption, Dimension> > impl_t;
 
    using impl_t::get_mesh;
    using impl_t::field_begin;
    using impl_t::field_end;
    using impl_t::get_num_dofs;
};
 
 
template <class Mesh, class Option, class Dimension = field_dimension::_1d>
function_space_view<Mesh, Option, Dimension> const &
    create_function_space_view(mesh_base<Mesh> const &msh, Option, Dimension dim = Dimension())
{
    return static_cast<function_space_view<Mesh, Option, Dimension> const &>(msh.derived());
}
 
template <class Mesh, class Dimension = field_dimension::_1d>
function_space_view<Mesh, field_option::isoparametric, Dimension> const &
    isoparametric_view(mesh_base<Mesh> const &msh, Dimension dim = Dimension())
{
    return create_function_space_view(msh.derived(), field_option::isoparametric(), dim);
}
 
template <class Mesh, class Dimension = field_dimension::_1d>
function_space_view<Mesh, field_option::constant, Dimension> const &
    constant_view(mesh_base<Mesh> const &msh, Dimension dim = Dimension())
{
    return create_function_space_view(msh.derived(), field_option::constant(), dim);
}
 
 
// forward declaration
template <class FuncSpace>
class dirac_space;
 
template <class FuncSpace>
struct function_space_traits<dirac_space<FuncSpace> >
{
    typedef function_space_traits<FuncSpace> base_traits;
 
    typedef typename base_traits::mesh_t mesh_t;
 
    enum {
        quantity_dimension = base_traits::quantity_dimension
    };
 
    typedef typename tmp::transform<
        typename base_traits::field_type_vector_t,
        tmp::inserter<
            typename tmp::empty<
                typename base_traits::field_type_vector_t
            >::type,
            tmp::push_back<tmp::_1, tmp::_2>
        >,
        dirac_field<tmp::_1>
    >::type field_type_vector_t;
 
    template <class dirac_field_t>
    struct iterator : casted_iterator<
        typename base_traits::template iterator<
            typename dirac_field_t::field_t
        >::type,
        dirac_field_t,
        field_impl<typename dirac_field_t::field_t>
    > {};
};
 
 
template <class FuncSpace>
class dirac_space :
    public function_space_base<dirac_space<FuncSpace> >,
    public function_space_impl<FuncSpace>
{
public:
    typedef function_space_traits<dirac_space<FuncSpace> > traits_t;
    typedef function_space_impl<FuncSpace> impl_t;
    
    // explicitly define mesh_t because it would be ambiguous otherwise
    typedef typename traits_t::mesh_t mesh_t;
 
    // explicitly define field_type_vector_t because it would be ambiguous otherwise
    typedef typename traits_t::field_type_vector_t field_type_vector_t;
 
    using impl_t::get_mesh;
    using impl_t::get_num_dofs;
 
    template <class dirac_field_t>
    typename traits_t::template iterator<dirac_field_t>::type
        field_begin() const
    {
        return impl_t::template field_begin<typename dirac_field_t::field_t>();
    }
 
    template <class dirac_field_t>
    typename traits_t::template iterator<dirac_field_t>::type
        field_end() const
    {
        return impl_t::template field_end<typename dirac_field_t::field_t>();
    }
};
 
 
template <class FuncSpace>
dirac_space<FuncSpace> const &
    dirac(function_space_base<FuncSpace> const &space)
{
    return static_cast<dirac_space<FuncSpace> const &>(
        static_cast<function_space_impl<FuncSpace> const &>(space.derived()));
}
 
 
 
// forward declaration
template <class FieldTypeVector>
class function_space;
 
template <class FieldTypeVector>
struct field_2_elem_type_vector
{
    template <class field_t>
    struct elemize
    {
        typedef typename field_t::elem_t type;
    };
 
    typedef typename tmp::unique<
        typename tmp::transform<
            FieldTypeVector,
            tmp::inserter<
                typename tmp::empty<FieldTypeVector>::type,
                tmp::push_back<tmp::_1, tmp::_2>
            >,
            elemize<tmp::_1>
        >::type
    >::type type;
};
 
 
template <class FieldTypeVector>
struct function_space_traits<function_space<FieldTypeVector> >
{
    typedef mesh<typename field_2_elem_type_vector<FieldTypeVector>::type> mesh_t;
 
    typedef FieldTypeVector field_type_vector_t;
 
    enum {
        quantity_dimension = tmp::deref<typename tmp::begin<field_type_vector_t>::type>::type::quantity_dimension
    };
 
    template <class field_t>
    struct iterator
    {
        typedef typename eigen_std_vector<field_t>::type::const_iterator type;
    };
};
 
 
template <class FieldTypeVector>
class function_space_impl<function_space<FieldTypeVector> > :
    public function_space_traits<function_space<FieldTypeVector> >::mesh_t
{
public:
    typedef function_space_traits<function_space<FieldTypeVector> > traits_t;
 
    typedef typename traits_t::field_type_vector_t field_type_vector_t;
 
    typedef typename traits_t::mesh_t mesh_t;
 
    typedef typename tmp::inherit<
        typename tmp::transform<
            field_type_vector_t,
            tmp::inserter<
                typename tmp::empty<field_type_vector_t>::type,
                tmp::push_back<tmp::_1,tmp::_2>
            >,
            eigen_std_vector<tmp::_1>
        >::type
    >::type field_container_t;
 
protected:
    template <class field_t>
    struct field_adder { struct type {
        bool operator() (unsigned const input[], function_space_impl &fsp)
        {
            typedef typename field_t::elem_t elem_t;
            if (input[0] == field_t::id)
            {
                // construct element and push
                typename elem_t::nodes_t nodes;
                typename elem_t::coords_t coords;
                for (unsigned i = 0; i < elem_t::num_nodes; ++i)
                {
                    nodes[i] = input[i+1];
                    coords.col(i) = fsp.points[nodes[i]];
                }
                elem_t const &elemref = fsp.push_element(elem_t(coords, fsp.m_num_elements++, nodes));
 
                // construct field and push
                typename field_t::dofs_t dofs;
                for (unsigned i = 0; i < field_t::num_dofs; ++i)
                    dofs[i] = input[i+elem_t::num_nodes+1];
                fsp.push_field(field_t(elemref, dofs));
 
                fsp.m_num_dofs = std::max(fsp.m_num_dofs, dofs.maxCoeff()+1);
                return true;
            }
            return false;
        }
    };};
 
public:
    function_space_impl() : m_num_dofs(0)
    {
    }
 
    template <class node_matrix_t, class field_matrix_t>
    function_space_impl(node_matrix_t const &nodes, field_matrix_t const &fields)
        : m_num_dofs(0)
    {
        unsigned const N_MAX_FIELD = 1000;
        double c[mesh_t::nDim];
 
        for (int i = 0; i < nodes.rows(); ++i)
        {
            for (unsigned j = 0; j < mesh_t::nDim; ++j)
                c[j] = nodes(i,j);
            this->add_node(c);
        }
        unsigned e[N_MAX_FIELD];
        for (int i = 0; i < fields.rows(); ++i)
        {
            for (int j = 0; j < fields.cols(); ++j)
                e[j] = (unsigned)fields(i,j);
            add_field(e);
        }
    }
 
    mesh_t const &get_mesh() const
    {
        return static_cast<mesh_t const &>(*this);
    }
 
    template <class FieldType>
    typename traits_t::template iterator<FieldType>::type
        field_begin() const
    {
        return m_fields.eigen_std_vector<FieldType>::type::begin();
    }
 
    template <class FieldType>
    typename traits_t::template iterator<FieldType>::type
        field_end() const
    {
        return m_fields.eigen_std_vector<FieldType>::type::end();
    }
 
    bool add_field(unsigned const input[])
    {
        return tmp::call_until<
            field_type_vector_t,
            field_adder<tmp::_1>,
            unsigned const*,
            function_space_impl &
        >(input, *this);
    }
 
    template <class field_t>
    field_t const &push_field(field_base<field_t> const &f)
    {
        m_fields.eigen_std_vector<field_t>::type::push_back(f.derived());
        return *(m_fields.eigen_std_vector<field_t>::type::rbegin());
    }
 
    unsigned get_num_fields() const
    {
        return this->get_num_elements();
    }
 
    unsigned get_num_dofs() const
    {
        return m_num_dofs;
    }
 
protected:
    field_container_t m_fields;
    unsigned m_num_dofs;
};
 
 
template <class FieldTypeVector>
class function_space :
    public function_space_base<function_space<FieldTypeVector> >,
    public function_space_impl<function_space<FieldTypeVector> >
{
public:
    typedef function_space_impl<function_space<FieldTypeVector> > impl_t;
 
    typedef typename impl_t::mesh_t mesh_t;
 
    typedef typename impl_t::field_type_vector_t field_type_vector_t;
 
    using impl_t::field_begin;
    using impl_t::field_end;
    using impl_t::get_num_dofs;
    using impl_t::get_mesh;
 
    function_space() :
        impl_t()
    {
    }
 
    template <class node_matrix_t, class field_matrix_t>
    function_space(node_matrix_t const &nodes, field_matrix_t const &fields)
        : impl_t(nodes, fields)
    {
    }
};
 
template <class nodes_t, class elements_t, class...fields_t>
function_space<tmp::vector<typename tag2type<fields_t>::type...> >
    create_function_space(nodes_t const &nodes, elements_t const &elements, fields_t const &...fields)
{
    return function_space<tmp::vector<typename tag2type<fields_t>::type...> >(nodes, elements);
}
 
 
template <class DerivedField, class OriginalIterator>
class directional_derivative_field_iterator;
 
template <class Field, class OriginalIterator>
class directional_derivative_field_iterator<directional_derivative_field<Field>, OriginalIterator>
    : public OriginalIterator
{
public:
    typedef Field original_field_t;
    typedef directional_derivative_field<original_field_t> derivative_field_t;
    
    typedef typename derivative_field_t::x_t x_t;
 
    typedef OriginalIterator original_iterator_t;
    
    typedef directional_derivative_field<original_field_t> value_t;
 
    template <class Direction>
    directional_derivative_field_iterator(original_iterator_t const& iter, Eigen::MatrixBase<Direction> const& direction)
        : original_iterator_t(iter)
        , m_direction(direction.derived())
    {
 
    }
 
    value_t operator*() const
    {
        // TODO: check what the original iterator returns upon dereferencing
        //       if it is an object created on-the-fly, error is expected
        //       as the directional_derivative_field stores the original field by reference
 
        return value_t(this->original_iterator_t::operator*(), m_direction);
    }
 
private:
    //original_iterator_t m_iterator;
    x_t const &m_direction;
};
 
 
template <class FunctionSpace>
class directional_derivative_function_space;
 
template <class FunctionSpace>
struct function_space_traits<directional_derivative_function_space<FunctionSpace>>
{
    typedef FunctionSpace original_function_space_t;
 
    typedef function_space_traits<original_function_space_t> base_traits;
 
    typedef typename original_function_space_t::mesh_t mesh_t;
 
    typedef typename tmp::transform<
        typename base_traits::field_type_vector_t,
        tmp::inserter<
            typename tmp::empty<
                typename base_traits::field_type_vector_t
            >::type,
            tmp::push_back<tmp::_1, tmp::_2>
        >,
        directional_derivative_field<tmp::_1>
    >::type field_type_vector_t;
 
    enum {
        quantity_dimension = original_function_space_t::quantity_dimension
    };
 
    template <class field_t>
    struct iterator;
 
    template <class original_field_t>
    struct iterator<directional_derivative_field<original_field_t> >
    {
    private:
        typedef typename function_space_traits<original_function_space_t>::template iterator<original_field_t>::type original_iterator_t;
        typedef directional_derivative_field<original_field_t> directional_field_t;
    public:
 
        //typedef directional_derivative_field_iterator<field_t, FunctionSpace> type;
        typedef directional_derivative_field_iterator<directional_field_t, original_iterator_t> type;
    };
};
 
template <class FunctionSpace>
class directional_derivative_function_space
    : public function_space_base<directional_derivative_function_space<FunctionSpace>>
{
public:
    typedef function_space_base<directional_derivative_function_space<FunctionSpace>> base_t;
    typedef typename base_t::x_t x_t;
    typedef typename base_t::mesh_t mesh_t;
    typedef typename base_t::traits_t traits_t;
 
    template <class Direction>
    directional_derivative_function_space(FunctionSpace const &fs, Eigen::MatrixBase<Direction> const &direction)
        : m_fs(fs)
        , m_derivative_direction(direction)
    {
    }
 
    mesh_t const &get_mesh() const
    {
        return m_fs.get_mesh();
    }
 
    template <class FieldType>
    typename traits_t::template iterator<FieldType>::type
        field_begin() const
    {
        typedef typename traits_t::template iterator<FieldType>::type return_t;
        typedef typename return_t::original_field_t original_field_t;
 
        return return_t(m_fs.template field_begin<original_field_t>(), m_derivative_direction);
    }
 
    template <class FieldType>
    typename traits_t::template iterator<FieldType>::type
        field_end() const
    {
        typedef typename traits_t::template iterator<FieldType>::type return_t;
        typedef typename return_t::original_field_t original_field_t;
 
        return return_t(m_fs.template field_end<original_field_t>(), m_derivative_direction);
    }
 
private:
    FunctionSpace const &m_fs;
    x_t const m_derivative_direction;
};
 
template <class FunctionSpace, class Direction>
directional_derivative_function_space<FunctionSpace>
directional_derivative(function_space_base<FunctionSpace> const &fs, Eigen::MatrixBase<Direction> const &direction)
{
    return directional_derivative_function_space<FunctionSpace>(fs.derived(), direction);
}
 
 
template <class FuncSpace>
struct is_function_space : std::is_base_of<
    function_space_base<typename std::decay<FuncSpace>::type>,
    typename std::decay<FuncSpace>::type
> {};
 
} // namespace NiHu
 
#endif // FUNCTION_SPACE_HPP_INCLUDED