assembly.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_ASSEMBLY_HPP_INCLUDED
#define NIHU_ASSEMBLY_HPP_INCLUDED
 
#include "function_space.hpp"
#include "integral_operator.hpp"
#include "../util/matrix_block.hpp"
 
namespace NiHu
{
 
template <class TestSpace, class TrialSpace, class OnSameMesh>
class assembly
{
private:
    template <class Operator, class TestField, class TrialField, bool isTrivial =
        integral_operator_traits<Operator>::is_local &&
        !std::is_same<typename TestField::elem_t, typename TrialField::elem_t>::value
        >
    struct eval_on_impl
    { struct type {
        template <class result_t>
        void operator() (result_t &result,
            integral_operator_base<Operator> const &op,
            function_space_base<TestSpace> const &test_space,
            function_space_base<TrialSpace> const &trial_space)
        {
            if (!integral_operator_traits<Operator>::is_local)
            {
                for (auto test_it = test_space.template field_begin<TestField>();
                    test_it != test_space.template field_end<TestField>(); ++test_it)
                    for (auto trial_it = trial_space.template field_begin<TrialField>();
                        trial_it != trial_space.template field_end<TrialField>(); ++trial_it)
                    {
                        block(result, test_it->get_dofs(), trial_it->get_dofs())
                            += op.eval_on_fields(*test_it, *trial_it, OnSameMesh());
                    }
            }
            else    // local operator
            {
                // we can be sure that the element types are the same and element id's are monotonic
                auto test_it = test_space.template field_begin<TestField>();
                auto test_end = test_space.template field_end<TestField>();
                auto trial_it = trial_space.template field_begin<TrialField>();
                auto trial_end = trial_space.template field_end<TrialField>();
 
                for (; test_it != test_end; ++test_it)
                {
                    auto test_id = test_it->get_elem().get_id();
                    for (; trial_it != trial_end && trial_it->get_elem().get_id() < test_id; ++trial_it);
                    if (trial_it->get_elem().get_id() == test_id)
                    {
                        block(result, test_it->get_dofs(), trial_it->get_dofs())
                            += op.eval_on_fields(*test_it, *trial_it, OnSameMesh());
                    }
                }
            }
        }
    };};
 
 
    template <class Operator, class TestField, class TrialField>
    struct eval_on_impl<Operator, TestField, TrialField, true>
    { struct type {
        template <class result_t>
        void operator() (
            result_t &,
            integral_operator_base<Operator> const &,
            function_space_base<TestSpace> const &,
            function_space_base<TrialSpace> const &)
        {
        }
    };};
 
    template <class Operator, class TestField, class TrialField>
    struct eval_on : eval_on_impl<Operator, TestField, TrialField>
    {
    };
 
 
public:
    template <class result_t, class Operator>
    static result_t &eval_into(
        result_t &result,
        integral_operator_base<Operator> const &op,
        function_space_base<TestSpace> const &test_space,
        function_space_base<TrialSpace> const &trial_space)
    {
        if(!OnSameMesh::value && integral_operator_traits<Operator>::is_local)
            return result;
        
        tmp::d_call_each<
            typename function_space_traits<TestSpace>::field_type_vector_t,
            typename function_space_traits<TrialSpace>::field_type_vector_t,
            eval_on<Operator, tmp::_1, tmp::_2>,
            result_t &,
            integral_operator_base<Operator> const &,
            function_space_base<TestSpace> const &,
            function_space_base<TrialSpace> const &
        >(result, op, test_space, trial_space);
 
        return result;
    }
};
 
}
 
#endif /* NIHU_ASSEMBLY_HPP_INCLUDED */