trace_assembler_basic_jobs.hpp

// FEAT3: Finite Element Analysis Toolbox, Version 3
// Copyright (C) 2010 by Stefan Turek & the FEAT group
// FEAT3 is released under the GNU General Public License version 3,
// see the file 'copyright.txt' in the top level directory for details.
 
#pragma once
 
#include <kernel/assembly/trace_assembler.hpp>
#include <kernel/assembly/function_integral_info.hpp>
#include <kernel/geometry/intern/face_index_mapping.hpp>
#include <kernel/geometry/intern/face_ref_trafo.hpp>
#include <kernel/geometry/intern/congruency_sampler.hpp>
#include <kernel/geometry/intern/congruency_trafo.hpp>
 
namespace FEAT
{
  namespace Assembly
  {
    template<typename DataType_, typename Space_, TrafoTags trafo_config_, TrafoTags facet_trafo_config_, SpaceTags space_config_>
    class TraceAssemblyBasicTaskBase1
    {
    public:
      // we do not support pairwise assembly
      static constexpr bool assemble_pairwise = false;
      typedef Assembly::AsmTraits1<DataType_, Space_, trafo_config_, space_config_> AsmTraits;
      typedef typename AsmTraits::DataType DataType;
      typedef typename AsmTraits::SpaceType SpaceType;
      typedef typename AsmTraits::TrafoType TrafoType;
      typedef typename TrafoType::ShapeType ShapeType;
      typedef typename Shape::FaceTraits<ShapeType, ShapeType::dimension-1>::ShapeType FacetType;
 
      static constexpr int shape_dim = ShapeType::dimension;
      static constexpr int facet_dim = FacetType::dimension;
 
      typedef typename AsmTraits::TrafoEvaluator TrafoEvaluator;
      typedef typename TrafoType::template Evaluator<FacetType, DataType>::Type TrafoFacetEvaluator;
 
      typedef typename AsmTraits::SpaceEvaluator SpaceEvaluator;
 
      static constexpr TrafoTags facet_trafo_config = facet_trafo_config_ | TrafoTags::jac_mat;
 
      typedef typename AsmTraits::TrafoEvalData TrafoEvalData;
      typedef typename TrafoFacetEvaluator::template ConfigTraits<facet_trafo_config>::EvalDataType TrafoFacetEvalData;
      typedef typename AsmTraits::SpaceEvalData SpaceEvalData;
      typedef typename AsmTraits::SpaceBasisData SpaceBasisData;
 
      typedef typename Intern::CubatureTraits<TrafoFacetEvaluator>::RuleType CubatureRuleType;
 
    protected:
      const SpaceType& space;
      const TrafoType& trafo;
      TrafoEvaluator trafo_eval;
      TrafoFacetEvaluator trafo_facet_eval;
      SpaceEvaluator space_eval;
      TrafoEvalData trafo_data;
      TrafoFacetEvalData trafo_facet_data;
      SpaceEvalData space_data;
      typename AsmTraits::DofMapping dof_mapping;
 
      int cell_facet_ori;
 
      Tiny::Matrix<DataType, shape_dim, facet_dim> face_mat;
      Tiny::Matrix<DataType, facet_dim, facet_dim> ori_mat;
      Tiny::Vector<DataType, shape_dim> face_vec;
      Tiny::Vector<DataType, facet_dim> ori_vec;
 
      Tiny::Vector<DataType, shape_dim> cur_point;
      Tiny::Vector<DataType, facet_dim> cur_point_facet;
 
    public:
      explicit TraceAssemblyBasicTaskBase1(const SpaceType& space_) :
        space(space_),
        trafo(space.get_trafo()),
        trafo_eval(trafo),
        trafo_facet_eval(trafo),
        space_eval(space),
        trafo_data(),
        trafo_facet_data(),
        space_data(),
        dof_mapping(space),
        cell_facet_ori(0)
      {
        face_mat.format();
        ori_mat.format();
        face_vec.format();
        ori_vec.format();
      }
 
      virtual ~TraceAssemblyBasicTaskBase1() = default;
 
      void prepare(Index facet, Index cell, int local_facet, int facet_ori)
      {
        // compute facet trafo
        Geometry::Intern::FaceRefTrafo<ShapeType, facet_dim>::compute(face_mat, face_vec, local_facet);
 
        // compute orientation trafo
        Geometry::Intern::CongruencyTrafo<FacetType>::compute(ori_mat, ori_vec, facet_ori);
 
        // compute orientation of actual cell facet
        cell_facet_ori = Geometry::Intern::CongruencySampler<FacetType>::orientation(facet_ori)
          * Shape::ReferenceCell<ShapeType>::facet_orientation(local_facet);
 
        // prepare evaluators
        trafo_facet_eval.prepare(facet);
        trafo_eval.prepare(cell);
        space_eval.prepare(trafo_eval);
 
        // prepare dof mapping
        dof_mapping.prepare(cell);
      }
 
      void prepare_point(Tiny::Vector<DataType, facet_dim>& pt)
      {
        // set cubature point
        cur_point_facet = pt;
 
        // transform point to local facet on reference cell
        cur_point = (face_mat * ((ori_mat * cur_point_facet) + ori_vec)) + face_vec;
 
        // compute trafo data
        trafo_facet_eval(trafo_facet_data, cur_point_facet);
        trafo_eval(trafo_data, cur_point);
 
        // compute normal vector
        trafo_facet_data.normal = Tiny::orthogonal(trafo_facet_data.jac_mat).normalize();
 
        // ensure that the normal vector "points outside"
        if(cell_facet_ori < 0)
          trafo_facet_data.normal.negate();
 
        // copy normal to cell trafo data
        trafo_data.normal = trafo_facet_data.normal;
 
        // compute space data
        space_eval(space_data, trafo_data);
      }
 
      void finish()
      {
        // finish evaluator
        dof_mapping.finish();
 
        // finish evaluators
        space_eval.finish();
        trafo_eval.finish();
        trafo_facet_eval.finish();
      }
    }; // class TraceAssemblyBasicTaskBase1<...>
 
    template<typename DataType_, typename TestSpace_, typename TrialSpace_,
      TrafoTags trafo_config_, TrafoTags facet_trafo_config_, SpaceTags test_config_, SpaceTags trial_config_>
    class TraceAssemblyBasicTaskBase2
    {
    public:
      // we do not support pairwise assembly
      static constexpr bool assemble_pairwise = false;
      typedef Assembly::AsmTraits2<DataType_, TestSpace_, TrialSpace_, trafo_config_, test_config_, trial_config_> AsmTraits;
      typedef typename AsmTraits::DataType DataType;
      typedef typename AsmTraits::TestSpaceType TestSpaceType;
      typedef typename AsmTraits::TrialSpaceType TrialSpaceType;
      typedef typename AsmTraits::TrafoType TrafoType;
      typedef typename TrafoType::ShapeType ShapeType;
      typedef typename Shape::FaceTraits<ShapeType, ShapeType::dimension-1>::ShapeType FacetType;
 
      static constexpr int shape_dim = ShapeType::dimension;
      static constexpr int facet_dim = FacetType::dimension;
 
      typedef typename AsmTraits::TrafoEvaluator TrafoEvaluator;
      typedef typename TrafoType::template Evaluator<FacetType, DataType>::Type TrafoFacetEvaluator;
 
      typedef typename AsmTraits::TestEvaluator TestEvaluator;
      typedef typename AsmTraits::TrialEvaluator TrialEvaluator;
 
      static constexpr TrafoTags facet_trafo_config = facet_trafo_config_ | TrafoTags::jac_mat;
 
      typedef typename AsmTraits::TrafoEvalData TrafoEvalData;
      typedef typename TrafoFacetEvaluator::template ConfigTraits<facet_trafo_config>::EvalDataType TrafoFacetEvalData;
      typedef typename AsmTraits::TestEvalData TestEvalData;
      typedef typename AsmTraits::TestBasisData TestBasisData;
      typedef typename AsmTraits::TrialEvalData TrialEvalData;
      typedef typename AsmTraits::TrialBasisData TrialBasisData;
 
      typedef typename Intern::CubatureTraits<TrafoFacetEvaluator>::RuleType CubatureRuleType;
 
    protected:
      const TestSpaceType& test_space;
      const TrialSpaceType& trial_space;
      const TrafoType& trafo;
      TrafoEvaluator trafo_eval;
      TrafoFacetEvaluator trafo_facet_eval;
      TestEvaluator test_eval;
      TrialEvaluator trial_eval;
      TrafoEvalData trafo_data;
      TrafoFacetEvalData trafo_facet_data;
      TestEvalData test_data;
      TrialEvalData trial_data;
      typename AsmTraits::TestDofMapping test_dof_mapping;
      typename AsmTraits::TrialDofMapping trial_dof_mapping;
 
      int cell_facet_ori;
 
      Tiny::Matrix<DataType, shape_dim, facet_dim> face_mat;
      Tiny::Matrix<DataType, facet_dim, facet_dim> ori_mat;
      Tiny::Vector<DataType, shape_dim> face_vec;
      Tiny::Vector<DataType, facet_dim> ori_vec;
 
      Tiny::Vector<DataType, shape_dim> cur_point;
      Tiny::Vector<DataType, facet_dim> cur_point_facet;
 
    public:
      explicit TraceAssemblyBasicTaskBase2(const TestSpaceType& test_space_, const TrialSpaceType& trial_space_) :
        test_space(test_space_),
        trial_space(trial_space_),
        trafo(test_space.get_trafo()),
        trafo_eval(trafo),
        trafo_facet_eval(trafo),
        test_eval(test_space),
        trial_eval(trial_space),
        trafo_data(),
        trafo_facet_data(),
        test_data(),
        trial_data(),
        test_dof_mapping(test_space),
        trial_dof_mapping(trial_space),
        cell_facet_ori(0)
      {
        face_mat.format();
        ori_mat.format();
        face_vec.format();
        ori_vec.format();
      }
 
      virtual ~TraceAssemblyBasicTaskBase2() = default;
 
      void prepare(Index facet, Index cell, int local_facet, int facet_ori)
      {
        // compute facet trafo
        Geometry::Intern::FaceRefTrafo<ShapeType, facet_dim>::compute(face_mat, face_vec, local_facet);
 
        // compute orientation trafo
        Geometry::Intern::CongruencyTrafo<FacetType>::compute(ori_mat, ori_vec, facet_ori);
 
        // compute orientation of actual cell facet
        cell_facet_ori = Geometry::Intern::CongruencySampler<FacetType>::orientation(facet_ori)
          * Shape::ReferenceCell<ShapeType>::facet_orientation(local_facet);
 
        // prepare evaluators
        trafo_facet_eval.prepare(facet);
        trafo_eval.prepare(cell);
        test_eval.prepare(trafo_eval);
        trial_eval.prepare(trafo_eval);
 
        // prepare dof mapping
        test_dof_mapping.prepare(cell);
        trial_dof_mapping.prepare(cell);
      }
 
      void prepare_point(Tiny::Vector<DataType, facet_dim>& pt)
      {
        // set cubature point
        cur_point_facet = pt;
 
        // transform point to local facet on reference cell
        cur_point = (face_mat * ((ori_mat * cur_point_facet) + ori_vec)) + face_vec;
 
        // compute trafo data
        trafo_facet_eval(trafo_facet_data, cur_point_facet);
        trafo_eval(trafo_data, cur_point);
 
        // compute normal vector
        trafo_facet_data.normal = Tiny::orthogonal(trafo_facet_data.jac_mat).normalize();
 
        // ensure that the normal vector "points outside"
        if(cell_facet_ori < 0)
          trafo_facet_data.normal.negate();
 
        // copy normal to cell trafo data
        trafo_data.normal = trafo_facet_data.normal;
 
        // compute space data
        test_eval(test_data, trafo_data);
        trial_eval(trial_data, trafo_data);
      }
 
      void finish()
      {
        // finish evaluator
        trial_dof_mapping.finish();
        test_dof_mapping.finish();
 
        // finish evaluators
        trial_eval.finish();
        test_eval.finish();
        trafo_eval.finish();
        trafo_facet_eval.finish();
      }
    }; // class TraceAssemblyBasicTaskBase2<...>
 
    template<typename Derived_, typename Vector_, typename Space_,
      TrafoTags trafo_config_, TrafoTags facet_trafo_config_, SpaceTags space_config_>
    class TraceAssemblyBasicVectorTaskCRTP :
      public TraceAssemblyBasicTaskBase1<typename Vector_::DataType, Space_, trafo_config_, facet_trafo_config_ | TrafoTags::jac_det, space_config_>
    {
    public:
      typedef typename Vector_::DataType DataType;
      typedef typename Vector_::ValueType ValueType;
 
      typedef TraceAssemblyBasicTaskBase1<DataType, Space_, trafo_config_, facet_trafo_config_ | TrafoTags::jac_det, space_config_> BaseClass;
      typedef typename BaseClass::AsmTraits AsmTraits;
 
      static constexpr bool need_scatter = true;
      static constexpr bool need_combine = false;
 
      using typename BaseClass::TrafoEvalData;
      using typename BaseClass::TrafoFacetEvalData;
      using typename BaseClass::SpaceBasisData;
 
    protected:
      Vector_& vector;
      typename BaseClass::CubatureRuleType cubature_rule;
      typename AsmTraits::template TLocalVector<ValueType> local_vector;
      typename Vector_::ScatterAxpy scatter_axpy;
      DataType scatter_alpha;
 
    public:
      explicit TraceAssemblyBasicVectorTaskCRTP(Vector_& vector_, const Space_& space_,
        const Cubature::DynamicFactory& cubature_factory_, DataType alpha_) :
        BaseClass(space_),
        vector(vector_),
        cubature_rule(Cubature::ctor_factory, cubature_factory_),
        scatter_axpy(vector),
        scatter_alpha(alpha_)
      {
      }
 
#ifdef DOXYGEN
      void set_point(TrafoFacetEvalData& tau_f, TrafoEvalData& tau);
 
      void eval(ValueType& val, const DataType& weight, const SpaceBasisData& psi);
#endif // DOXYGEN
 
      void assemble()
      {
        // format local vector
        local_vector.format();
 
        // fetch number of local dofs
        const int num_loc_dofs = this->space_eval.get_num_local_dofs();
 
        // loop over all quadrature points and integrate
        for(int k(0); k < cubature_rule.get_num_points(); ++k)
        {
          // prepare point
          static_cast<Derived_&>(*this).prepare_point(cubature_rule.get_point(k));
 
          // evaluate operator
          static_cast<Derived_&>(*this).set_point(this->trafo_facet_data, this->trafo_data);
 
          // test function loop
          for(int i(0); i < num_loc_dofs; ++i)
          {
            // evaluate functional and integrate
            static_cast<Derived_&>(*this).eval(local_vector(i),
              this->trafo_facet_data.jac_det * cubature_rule.get_weight(k), this->space_data.phi[i]);
          } // continue with next test function
        } // continue with next cubature point
      }
 
      void scatter()
      {
        scatter_axpy(local_vector, this->dof_mapping, scatter_alpha);
      }
 
      void combine()
      {
        // nothing to do here
      }
    }; // class TraceAssemblyBasicVectorTaskCRTP<...>
 
    template<typename LinearFunctional_, typename Vector_, typename Space_>
    class TraceAssemblyLinearFunctionalVectorJob
    {
    public:
      typedef typename Vector_::DataType DataType;
      typedef typename Vector_::ValueType ValueType;
 
      static constexpr TrafoTags trafo_config = LinearFunctional_::trafo_config;
      static constexpr SpaceTags space_config = LinearFunctional_::test_config;
 
      class Task :
        public TraceAssemblyBasicVectorTaskCRTP<Task, Vector_, Space_, trafo_config, TrafoTags::none, space_config>
      {
      protected:
        typedef TraceAssemblyBasicVectorTaskCRTP<Task, Vector_, Space_, trafo_config, TrafoTags::none, space_config> BaseClass;
 
        typedef typename BaseClass::AsmTraits AsmTraits;
 
        typename LinearFunctional_::template Evaluator<AsmTraits> func_eval;
 
      public:
        explicit Task(TraceAssemblyLinearFunctionalVectorJob& job) :
          BaseClass(job.vector, job.space, job.cubature_factory, job.alpha),
          func_eval(job.linear_functional)
        {
        }
 
        void prepare(Index facet, Index cell, int local_facet, int facet_ori)
        {
          BaseClass::prepare(facet, cell, local_facet, facet_ori);
          func_eval.prepare(this->trafo_eval);
        }
 
        void set_point(typename BaseClass::TrafoFacetEvalData& DOXY(tau_f), typename BaseClass::TrafoEvalData& tau)
        {
          func_eval.set_point(tau);
        }
 
        void eval(ValueType& val, const DataType& weight, typename BaseClass::SpaceBasisData& psi)
        {
          Tiny::axpy(val, func_eval.eval(psi), weight);
        }
 
        void finish()
        {
          func_eval.finish();
          BaseClass::finish();
        }
      }; // class Task
 
    protected:
      const LinearFunctional_& linear_functional;
      Vector_& vector;
      const Space_& space;
      Cubature::DynamicFactory cubature_factory;
      DataType alpha;
 
    public:
      explicit TraceAssemblyLinearFunctionalVectorJob(const LinearFunctional_& linear_functional_,
        Vector_& vector_, const Space_& space_, String cubature_, DataType alpha_ = DataType(1)) :
        linear_functional(linear_functional_),
        vector(vector_),
        space(space_),
        cubature_factory(cubature_),
        alpha(alpha_)
      {
      }
    }; // class TraceAssemblyLinearFunctionalVectorJob<...>
 
    template<typename Trafo_, typename Vector_, typename LinFunc_, typename Space_>
    void assemble_linear_functional_vector(TraceAssembler<Trafo_>& trace_asm,
      Vector_& vector, const LinFunc_& linear_functional, const Space_& space,
      const String& cubature, const typename Vector_::DataType alpha = typename Vector_::DataType(1))
    {
      XASSERTM(trace_asm.get_trafo() == space.get_trafo(), "trace assembler and space have different trafos");
      XASSERTM(vector.size() == space.get_num_dofs(), "invalid vector length");
 
      TraceAssemblyLinearFunctionalVectorJob<LinFunc_, Vector_, Space_> job(
        linear_functional, vector, space, cubature, alpha);
      trace_asm.assemble(job);
    }
 
    template<typename Derived_, typename Matrix_, typename Space_,
      TrafoTags trafo_config_, TrafoTags facet_trafo_config_, SpaceTags space_config_>
    class TraceAssemblyBasicMatrixTaskCRTP1 :
      public TraceAssemblyBasicTaskBase1<typename Matrix_::DataType, Space_, trafo_config_, facet_trafo_config_ | TrafoTags::jac_det, space_config_>
    {
    public:
      typedef typename Matrix_::DataType DataType;
      typedef typename Matrix_::ValueType ValueType;
 
      typedef TraceAssemblyBasicTaskBase1<DataType, Space_, trafo_config_, facet_trafo_config_ | TrafoTags::jac_det, space_config_> BaseClass;
      typedef typename BaseClass::AsmTraits AsmTraits;
 
      static constexpr bool need_scatter = true;
      static constexpr bool need_combine = false;
 
      using typename BaseClass::TrafoEvalData;
      using typename BaseClass::TrafoFacetEvalData;
      using typename BaseClass::SpaceBasisData;
 
    protected:
      Matrix_& matrix;
      typename BaseClass::CubatureRuleType cubature_rule;
      typename AsmTraits::template TLocalMatrix<ValueType> local_matrix;
      typename Matrix_::ScatterAxpy scatter_axpy;
      DataType scatter_alpha;
 
    public:
      explicit TraceAssemblyBasicMatrixTaskCRTP1(Matrix_& matrix_, const Space_& space_,
        const Cubature::DynamicFactory& cubature_factory_, DataType alpha_) :
        BaseClass(space_),
        matrix(matrix_),
        cubature_rule(Cubature::ctor_factory, cubature_factory_),
        scatter_axpy(matrix),
        scatter_alpha(alpha_)
      {
      }
 
#ifdef DOXYGEN
      void set_point(TrafoFacetEvalData& tau_f, TrafoEvalData& tau);
 
      void eval(ValueType& val, const DataType& weight, const SpaceBasisData& psi, const SpaceBasisData& phi);
#endif // DOXYGEN
 
      void assemble()
      {
        // format local matrix
        local_matrix.format();
 
        // fetch number of local dofs
        const int num_loc_dofs = this->space_eval.get_num_local_dofs();
 
        // loop over all quadrature points and integrate
        for(int k(0); k < cubature_rule.get_num_points(); ++k)
        {
          // prepare point
          static_cast<Derived_&>(*this).prepare_point(cubature_rule.get_point(k));
 
          // evaluate operator
          static_cast<Derived_&>(*this).set_point(this->trafo_facet_data, this->trafo_data);
 
          // test function loop
          for(int i(0); i < num_loc_dofs; ++i)
          {
            // trial function loop
            for(int j(0); j < num_loc_dofs; ++j)
            {
              // evaluate operator and integrate
              static_cast<Derived_&>(*this).eval(local_matrix(i,j),
                this->trafo_facet_data.jac_det * cubature_rule.get_weight(k),
                this->space_data.phi[j], this->space_data.phi[i]);
            } // continue with next trial function
          } // continue with next test function
        } // continue with next cubature point
      }
 
      void scatter()
      {
        scatter_axpy(local_matrix, this->dof_mapping, this->dof_mapping, scatter_alpha);
      }
 
      void combine()
      {
        // nothing to do here
      }
    }; // class TraceAssemblyBasicMatrixTaskCRTP1<...>
 
    template<typename Derived_, typename Matrix_, typename TestSpace_, typename TrialSpace_,
      TrafoTags trafo_config_, TrafoTags facet_trafo_config_, SpaceTags test_config_, SpaceTags trial_config_>
    class TraceAssemblyBasicMatrixTaskCRTP2 :
      public TraceAssemblyBasicTaskBase2<typename Matrix_::DataType, TestSpace_, TrialSpace_, trafo_config_, facet_trafo_config_ | TrafoTags::jac_det, test_config_, trial_config_>
    {
    public:
      typedef typename Matrix_::DataType DataType;
      typedef typename Matrix_::ValueType ValueType;
 
      typedef TraceAssemblyBasicTaskBase2<typename Matrix_::DataType, TestSpace_, TrialSpace_, trafo_config_, facet_trafo_config_ | TrafoTags::jac_det, test_config_, trial_config_> BaseClass;
      typedef typename BaseClass::AsmTraits AsmTraits;
 
      static constexpr bool need_scatter = true;
      static constexpr bool need_combine = false;
 
      using typename BaseClass::TrafoEvalData;
      using typename BaseClass::TrafoFacetEvalData;
      using typename BaseClass::TestBasisData;
      using typename BaseClass::TrialBasisData;
 
    protected:
      Matrix_& matrix;
      typename BaseClass::CubatureRuleType cubature_rule;
      typename AsmTraits::template TLocalMatrix<ValueType> local_matrix;
      typename Matrix_::ScatterAxpy scatter_axpy;
      DataType scatter_alpha;
 
    public:
      explicit TraceAssemblyBasicMatrixTaskCRTP2(Matrix_& matrix_, const TestSpace_& test_space_,
        const TestSpace_& trial_space_, const Cubature::DynamicFactory& cubature_factory_, DataType alpha_) :
        BaseClass(test_space_, trial_space_),
        matrix(matrix_),
        cubature_rule(Cubature::ctor_factory, cubature_factory_),
        scatter_axpy(matrix),
        scatter_alpha(alpha_)
      {
      }
 
#ifdef DOXYGEN
      void set_point(TrafoFacetEvalData& tau_f, TrafoEvalData& tau);
 
      void eval(ValueType& val, const DataType& weight, const TrialBasisData& phi, const TestBasisData& psi),
#endif // DOXYGEN
 
      void assemble()
      {
        // format local matrix
        local_matrix.format();
 
        // fetch number of local dofs
        const int num_loc_test_dofs = this->test_eval.get_num_local_dofs();
        const int num_loc_trial_dofs = this->trial_eval.get_num_local_dofs();
 
        // loop over all quadrature points and integrate
        for(int k(0); k < cubature_rule.get_num_points(); ++k)
        {
          // prepare point
          static_cast<Derived_&>(*this).prepare_point(cubature_rule.get_point(k));
 
          // evaluate operator
          static_cast<Derived_&>(*this).set_point(this->trafo_facet_data, this->trafo_data);
 
          // test function loop
          for(int i(0); i < num_loc_test_dofs; ++i)
          {
            // trial function loop
            for(int j(0); j < num_loc_trial_dofs; ++j)
            {
              // evaluate operator and integrate
              static_cast<Derived_&>(*this).eval(local_matrix(i,j),
                this->trafo_facet_data.jac_det * cubature_rule.get_weight(k),
                this->trial_data.phi[j], this->test_data.phi[i]);
            } // continue with next trial function
          } // continue with next test function
        } // continue with next cubature point
      }
 
      void scatter()
      {
        scatter_axpy(local_matrix, this->test_dof_mapping, this->trial_dof_mapping, scatter_alpha);
      }
 
      void combine()
      {
        // nothing to do here
      }
    }; // class TraceAssemblyBasicMatrixTaskCRTP2<...>
 
    template<typename BilinearOperator_, typename Matrix_, typename Space_>
    class TraceAssemblyBilinearOperatorMatrixJob1
    {
    public:
      typedef typename Matrix_::DataType DataType;
      typedef typename Matrix_::ValueType ValueType;
 
      static constexpr TrafoTags trafo_config = BilinearOperator_::trafo_config;
      static constexpr SpaceTags space_config = BilinearOperator_::test_config | BilinearOperator_::trial_config;
 
      class Task :
        public TraceAssemblyBasicMatrixTaskCRTP1<Task, Matrix_, Space_, trafo_config, TrafoTags::none, space_config>
      {
      protected:
        typedef TraceAssemblyBasicMatrixTaskCRTP1<Task, Matrix_, Space_, trafo_config, TrafoTags::none, space_config> BaseClass;
 
        typedef typename BaseClass::AsmTraits AsmTraits;
 
        typename BilinearOperator_::template Evaluator<AsmTraits> oper_eval;
 
      public:
        explicit Task(TraceAssemblyBilinearOperatorMatrixJob1& job) :
          BaseClass(job.matrix, job.space, job.cubature_factory, job.alpha),
          oper_eval(job.bilinear_operator)
        {
        }
 
        void prepare(Index facet, Index cell, int local_facet, int facet_ori)
        {
          BaseClass::prepare(facet, cell, local_facet, facet_ori);
          oper_eval.prepare(this->trafo_eval);
        }
 
        void set_point(typename BaseClass::TrafoFacetEvalData& DOXY(tau_f), typename BaseClass::TrafoEvalData& tau)
        {
          oper_eval.set_point(tau);
        }
 
        void eval(ValueType& val, const DataType& weight, typename AsmTraits::TrialBasisData& phi, typename AsmTraits::TestBasisData& psi)
        {
          Tiny::axpy(val, oper_eval.eval(phi, psi), weight);
        }
 
        void finish()
        {
          oper_eval.finish();
          BaseClass::finish();
        }
      }; // class Task
 
    protected:
      const BilinearOperator_& bilinear_operator;
      Matrix_& matrix;
      const Space_& space;
      Cubature::DynamicFactory cubature_factory;
      DataType alpha;
 
    public:
      explicit TraceAssemblyBilinearOperatorMatrixJob1(const BilinearOperator_& bilinear_operator_,
        Matrix_& matrix_, const Space_& space_, String cubature_, DataType alpha_ = DataType(1)) :
        bilinear_operator(bilinear_operator_),
        matrix(matrix_),
        space(space_),
        cubature_factory(cubature_),
        alpha(alpha_)
      {
      }
    }; // class TraceAssemblyBilinearOperatorMatrixJob1<...>
 
    template<typename Trafo_, typename Matrix_, typename BilOp_, typename Space_>
    void assemble_bilinear_operator_matrix_1(TraceAssembler<Trafo_>& trace_asm, Matrix_& matrix,
      const BilOp_& bilinear_operator, const Space_& space, const String& cubature,
      const typename Matrix_::DataType alpha = typename Matrix_::DataType(1))
    {
      XASSERTM(trace_asm.get_trafo() == space.get_trafo(), "trace assembler and space have different trafos");
      XASSERTM(matrix.columns() == space.get_num_dofs(), "invalid matrix column count");
      XASSERTM(matrix.rows() == space.get_num_dofs(), "invalid matrix row count");
 
      TraceAssemblyBilinearOperatorMatrixJob1<BilOp_, Matrix_, Space_> job(
        bilinear_operator, matrix, space, cubature, alpha);
      trace_asm.assemble(job);
    }
 
    template<typename BilinearOperator_, typename Matrix_, typename TestSpace_, typename TrialSpace_>
    class TraceAssemblyBilinearOperatorMatrixJob2
    {
    public:
      typedef typename Matrix_::DataType DataType;
      typedef typename Matrix_::ValueType ValueType;
 
      static constexpr TrafoTags trafo_config = BilinearOperator_::trafo_config;
      static constexpr SpaceTags test_config  = BilinearOperator_::test_config;
      static constexpr SpaceTags trial_config = BilinearOperator_::trial_config;
 
      class Task :
        public TraceAssemblyBasicMatrixTaskCRTP2<Task, Matrix_, TestSpace_, TrialSpace_, trafo_config, TrafoTags::none, test_config, trial_config>
      {
      protected:
        typedef TraceAssemblyBasicMatrixTaskCRTP2<Task, Matrix_, TestSpace_, TrialSpace_, trafo_config, TrafoTags::none, test_config, trial_config> BaseClass;
 
        typedef typename BaseClass::AsmTraits AsmTraits;
 
        typename BilinearOperator_::template Evaluator<AsmTraits> oper_eval;
 
      public:
        explicit Task(TraceAssemblyBilinearOperatorMatrixJob2& job) :
          BaseClass(job.matrix, job.test_space, job.trial_space, job.cubature_factory, job.alpha),
          oper_eval(job.bilinear_operator)
        {
        }
 
        void prepare(Index facet, Index cell, int local_facet, int facet_ori)
        {
          BaseClass::prepare(facet, cell, local_facet, facet_ori);
          oper_eval.prepare(this->trafo_eval);
        }
 
        void set_point(typename BaseClass::TrafoFacetEvalData& DOXY(tau_f), typename BaseClass::TrafoEvalData& tau)
        {
          oper_eval.set_point(tau);
        }
 
        void eval(ValueType& val, const DataType& weight, typename AsmTraits::TrialBasisData& phi, typename AsmTraits::TestBasisData& psi)
        {
          Tiny::axpy(val, oper_eval.eval(phi, psi), weight);
        }
 
        void finish()
        {
          oper_eval.finish();
          BaseClass::finish();
        }
      }; // class Task
 
    protected:
      const BilinearOperator_& bilinear_operator;
      Matrix_& matrix;
      const TestSpace_& test_space;
      const TrialSpace_& trial_space;
      Cubature::DynamicFactory cubature_factory;
      DataType alpha;
 
    public:
      explicit TraceAssemblyBilinearOperatorMatrixJob2(const BilinearOperator_& bilinear_operator_,
        Matrix_& matrix_, const TestSpace_& test_space_, const TrialSpace_& trial_space_,
        String cubature_, DataType alpha_ = DataType(1)) :
        bilinear_operator(bilinear_operator_),
        matrix(matrix_),
        test_space(test_space_),
        trial_space(trial_space_),
        cubature_factory(cubature_),
        alpha(alpha_)
      {
      }
    }; // class TraceAssemblyBilinearOperatorMatrixJob2<...>
 
    template<typename Trafo_, typename Matrix_, typename BilOp_, typename TestSpace_, typename TrialSpace_>
    void assemble_bilinear_operator_matrix_2(TraceAssembler<Trafo_>& trace_asm,
      Matrix_& matrix, const BilOp_& bilinear_operator, const TestSpace_& test_space,
      const TrialSpace_& trial_space, const String& cubature,
      const typename Matrix_::DataType alpha = typename Matrix_::DataType(1))
    {
      XASSERTM(trace_asm.get_trafo() == test_space.get_trafo(), "trace assembler and test space have different trafos");
      XASSERTM(trace_asm.get_trafo() == trial_space.get_trafo(), "trace assembler and trial space have different trafos");
      XASSERTM(matrix.columns() == trial_space.get_num_dofs(), "invalid matrix column count");
      XASSERTM(matrix.rows() == test_space.get_num_dofs(), "invalid matrix row count");
 
      TraceAssemblyBilinearOperatorMatrixJob2<BilOp_, Matrix_, TestSpace_, TrialSpace_> job(
        bilinear_operator, matrix, test_space, trial_space, cubature, alpha);
      trace_asm.assemble(job);
    }
 
    template<typename BilinearOperator_, typename Vector_, typename VectorSol_, typename Space_>
    class TraceAssemblyBilinearOperatorApplyVectorJob1
    {
    public:
      typedef typename Vector_::DataType DataType;
      typedef typename Vector_::ValueType ValueType;
 
      static constexpr TrafoTags trafo_config = BilinearOperator_::trafo_config;
      static constexpr SpaceTags space_config = BilinearOperator_::test_config | BilinearOperator_::trial_config;
 
      class Task :
        public TraceAssemblyBasicVectorTaskCRTP<Task, Vector_, Space_, trafo_config, TrafoTags::none, space_config>
      {
      protected:
        typedef TraceAssemblyBasicVectorTaskCRTP<Task, Vector_, Space_, trafo_config, TrafoTags::none, space_config> BaseClass;
 
        typedef typename BaseClass::AsmTraits AsmTraits;
        typedef typename BilinearOperator_::template Evaluator<AsmTraits>::ValueType MatValType;
 
        typename BilinearOperator_::template Evaluator<AsmTraits> oper_eval;
 
        const VectorSol_& vec_sol;
 
        typename VectorSol_::GatherAxpy sol_gather;
 
        typename AsmTraits::template TLocalMatrix<MatValType> local_matrix;
        typename AsmTraits::template TLocalVector<ValueType> local_vec_sol;
 
 
      public:
        explicit Task(TraceAssemblyBilinearOperatorApplyVectorJob1& job) :
          BaseClass(job.vector, job.space, job.cubature_factory, job.alpha),
          oper_eval(job.bilinear_operator),
          vec_sol(job.vec_sol),
          sol_gather(vec_sol)
        {
        }
 
        void prepare(Index facet, Index cell, int local_facet, int facet_ori)
        {
          BaseClass::prepare(facet, cell , local_facet, facet_ori);
          oper_eval.prepare(this->trafo_eval);
          local_vec_sol.format();
          sol_gather(local_vec_sol, this->dof_mapping);
        }
 
        void set_point(typename BaseClass::TrafoFacetEvalData& DOXY(tau_f), typename BaseClass::TrafoEvalData& tau)
        {
          oper_eval.set_point(tau);
        }
 
        void eval(MatValType& val, const DataType& weight,
          typename AsmTraits::TrialBasisData& phi, typename AsmTraits::TestBasisData& psi)
        {
          val += weight * oper_eval.eval(phi, psi);
        }
 
        void assemble()
        {
          // format local matrix
          local_matrix.format();
          this->local_vector.format();
 
          // fetch number of local dofs
          const int num_loc_dofs = this->space_eval.get_num_local_dofs();
 
          // loop over all quadrature points and integrate
          for(int k(0); k < this->cubature_rule.get_num_points(); ++k)
          {
            // prep point
            this->prepare_point(this->cubature_rule.get_point(k));
            // evaluate operator
            this->set_point(this->trafo_facet_data, this->trafo_data);
 
            // test function loop
            for(int i(0); i < num_loc_dofs; ++i)
            {
              // trial function loop
              for(int j(0); j < num_loc_dofs; ++j)
              {
                // evaluate operator and integrate
                this->eval(local_matrix(i,j),
                  this->trafo_facet_data.jac_det * this->cubature_rule.get_weight(k),
                  this->space_data.phi[j], this->space_data.phi[i]);
                // continue with next trial function
              }
              // continue with next test function
            }
            // continue with next cubature point
          }
 
          // apply to local vector
          for(int i(0); i < num_loc_dofs; ++i)
            for(int j(0); j < num_loc_dofs; ++j)
              this->local_vector[i].add_mat_vec_mult(this->local_matrix(i,j), this->local_vec_sol[j]);
        }
 
        void finish()
        {
          oper_eval.finish();
          BaseClass::finish();
        }
      }; // class Task
 
    protected:
      const BilinearOperator_& bilinear_operator;
      Vector_& vector;
      const VectorSol_& vec_sol;
      const Space_& space;
      Cubature::DynamicFactory cubature_factory;
      DataType alpha;
 
    public:
      explicit TraceAssemblyBilinearOperatorApplyVectorJob1(const BilinearOperator_& bilinear_operator_,
        Vector_& vector_, const VectorSol_& vec_sol_, const Space_& space_, const String& cubature_, DataType alpha_ = DataType(1)) :
        bilinear_operator(bilinear_operator_),
        vector(vector_),
        vec_sol(vec_sol_),
        space(space_),
        cubature_factory(cubature_),
        alpha(alpha_)
      {
      }
    }; // class TraceAssemblyBilinearOperatorApplyVectorJob1<...>
 
    template<typename Trafo_, typename Vector_, typename VectorSol_, typename BilOp_, typename Space_>
    void assemble_bilinear_operator_apply_vector_1(TraceAssembler<Trafo_>& trace_asm, Vector_& vector,
      const VectorSol_& vec_sol, const BilOp_& bilinear_operator, const Space_& space, const String& cubature,
      const typename Vector_::DataType alpha = typename Vector_::DataType(1))
    {
      XASSERTM(trace_asm.get_trafo() == space.get_trafo(), "domain assembler and space have different trafos");
      XASSERTM(vector.size() == space.get_num_dofs(), "invalid vector size");
 
      TraceAssemblyBilinearOperatorApplyVectorJob1<BilOp_, Vector_, VectorSol_, Space_> job(
        bilinear_operator, vector, vec_sol, space, cubature, alpha);
      trace_asm.assemble(job);
    }
 
    template<typename DataType_, typename Function_, typename Trafo_, int max_der_>
    class TraceAssemblyAnalyticFunctionIntegralJob
    {
    public:
      typedef DataType_ DataType;
 
      typedef typename AnalyticFunctionIntegral<DataType, Function_>::Type FunctionIntegralType;
 
      class Task
      {
      public:
        static constexpr TrafoTags trafo_config = TrafoTags::img_point | TrafoTags::jac_det;
 
        typedef typename Trafo_::ShapeType ShapeType;
        typedef typename Shape::FaceTraits<ShapeType, ShapeType::dimension-1>::ShapeType FacetType;
 
        typedef typename Trafo_::template Evaluator<FacetType, DataType>::Type TrafoEvaluator;
        typedef typename TrafoEvaluator::template ConfigTraits<trafo_config>::EvalDataType TrafoEvalData;
 
        typedef Analytic::EvalTraits<DataType, Function_> AnalyticEvalTraits;
 
        // we do not support pairwise assembly
        static constexpr bool assemble_pairwise = false;
        static constexpr bool need_scatter = false;
        static constexpr bool need_combine = true;
 
      protected:
        TrafoEvaluator trafo_eval;
        TrafoEvalData trafo_data;
        typename Assembly::Intern::CubatureTraits<TrafoEvaluator>::RuleType cubature_rule;
        typename Function_::template Evaluator<AnalyticEvalTraits> func_eval;
        FunctionIntegralType loc_integral;
        FunctionIntegralType& job_integral;
 
      public:
        explicit Task(TraceAssemblyAnalyticFunctionIntegralJob& job) :
          trafo_eval(job.trafo),
          trafo_data(),
          cubature_rule(Cubature::ctor_factory, job.cubature_factory),
          func_eval(job.function),
          loc_integral(),
          job_integral(job.integral)
        {
        }
 
        void prepare(Index facet, Index DOXY(cell), int DOXY(local_facet), int DOXY(facet_ori))
        {
          trafo_eval.prepare(facet);
        }
 
        void assemble()
        {
          // loop over all quadrature points and integrate
          for(int k(0); k < cubature_rule.get_num_points(); ++k)
          {
            // compute trafo data
            trafo_eval(trafo_data, cubature_rule.get_point(k));
 
            // call helper to integrate
            Intern::AnaFunIntJobHelper<max_der_>::work(loc_integral,
              cubature_rule.get_weight(k) * trafo_data.jac_det, func_eval, trafo_data.img_point);
          }
        }
 
        void scatter()
        {
          // nothing to do here
        }
 
        void finish()
        {
          trafo_eval.finish();
        }
 
        void combine()
        {
          job_integral.push(loc_integral);
        }
      }; // class Task
 
    protected:
      const Trafo_& trafo;
      const Function_& function;
      Cubature::DynamicFactory cubature_factory;
      FunctionIntegralType integral;
 
    public:
      explicit TraceAssemblyAnalyticFunctionIntegralJob(const Function_& function_, const Trafo_& trafo_, const String& cubature_) :
        trafo(trafo_),
        function(function_),
        cubature_factory(cubature_),
        integral()
      {
        integral.max_der = max_der_;
      }
 
      // \returns The integral of the discrete function
      FunctionIntegralType& result()
      {
        return integral;
      }
    }; // class TraceAssemblyAnalyticFunctionIntegralJob<...>
 
    template<int max_der_, typename DataType_, typename Function_, typename Trafo_>
    typename AnalyticFunctionIntegral<DataType_, Function_>::Type integrate_analytic_function(
      TraceAssembler<Trafo_>& trace_asm, const Function_& function, const String& cubature)
    {
      TraceAssemblyAnalyticFunctionIntegralJob<DataType_, Function_, Trafo_, max_der_> job(function, trace_asm.get_trafo(), cubature);
      trace_asm.assemble(job);
      return job.result();
    }
 
    template<typename Vector_, typename Space_, int max_der_>
    class TraceAssemblyDiscreteFunctionIntegralJob
    {
    public:
      typedef typename Vector_::DataType DataType;
      typedef typename Vector_::ValueType ValueType;
 
      static constexpr TrafoTags trafo_config = TrafoTags::none;
      static constexpr SpaceTags space_config = SpaceTags::value |
        (max_der_ >= 1 ? SpaceTags::grad : SpaceTags::none) |
        (max_der_ >= 2 ? SpaceTags::hess : SpaceTags::none);
      static constexpr TrafoTags facet_trafo_config = TrafoTags::jac_det;
 
      typedef typename DiscreteFunctionIntegral<Vector_, Space_>::Type FunctionIntegralType;
 
      class Task :
        public TraceAssemblyBasicTaskBase1<DataType, Space_, trafo_config, facet_trafo_config, space_config>
      {
      public:
        typedef TraceAssemblyBasicTaskBase1<DataType, Space_, trafo_config, facet_trafo_config, space_config> BaseClass;
 
        typedef typename BaseClass::AsmTraits AsmTraits;
 
        // we do not support pairwise assembly
        static constexpr bool assemble_pairwise = false;
        static constexpr bool need_scatter = false;
        static constexpr bool need_combine = true;
 
        using typename BaseClass::TrafoEvalData;
        using typename BaseClass::TrafoFacetEvalData;
 
      protected:
        const Vector_& vector;
        typename BaseClass::CubatureRuleType cubature_rule;
        typename AsmTraits::template TLocalVector<ValueType> local_vector;
        typename Vector_::GatherAxpy gather_axpy;
        FunctionIntegralType loc_integral;
        FunctionIntegralType& job_integral;
 
      public:
        explicit Task(TraceAssemblyDiscreteFunctionIntegralJob& job) :
          BaseClass(job.space),
          vector(job.vector),
          cubature_rule(Cubature::ctor_factory, job.cubature_factory),
          local_vector(),
          gather_axpy(vector),
          loc_integral(),
          job_integral(job.integral)
        {
        }
 
        void assemble()
        {
          // format local vector
          local_vector.format();
 
          // gather local vector data
          gather_axpy(local_vector, this->dof_mapping);
 
          // fetch number of local dofs
          const int num_loc_dofs = this->space_eval.get_num_local_dofs();
 
          // loop over all quadrature points and integrate
          for(int k(0); k < cubature_rule.get_num_points(); ++k)
          {
            // prepare point
            this->prepare_point(cubature_rule.get_point(k));
 
            // do the dirty work
            Intern::DiscFunIntJobHelper<max_der_>::work(loc_integral,
              cubature_rule.get_weight(k) * this->trafo_facet_data.jac_det, this->space_data, local_vector, num_loc_dofs);
          } // continue with next cubature point
        }
 
        void scatter()
        {
          // nothing to do here
        }
 
        void combine()
        {
          job_integral.push(loc_integral);
        }
      }; // class Task
 
    protected:
      const Vector_& vector;
      const Space_& space;
      Cubature::DynamicFactory cubature_factory;
      FunctionIntegralType integral;
 
    public:
      explicit TraceAssemblyDiscreteFunctionIntegralJob(const Vector_& vector_, const Space_& space_, String cubature_):
        vector(vector_),
        space(space_),
        cubature_factory(cubature_),
        integral()
      {
        integral.max_der = max_der_;
      }
 
      // \returns The integral of the discrete function
      FunctionIntegralType& result()
      {
        return integral;
      }
    }; // class TraceAssemblyDiscreteFunctionIntegralJob<...>
 
    template<int max_der_, typename Vector_, typename Trafo_, typename Space_>
    typename DiscreteFunctionIntegral<Vector_, Space_>::Type integrate_discrete_function(
      TraceAssembler<Trafo_>& trace_asm, const Vector_& vector, const Space_& space, const String& cubature)
    {
      XASSERTM(trace_asm.get_trafo() == space.get_trafo(), "domain assembler and space have different trafos");
      XASSERTM(vector.size() == space.get_num_dofs(), "invalid coefficient vector length");
      TraceAssemblyDiscreteFunctionIntegralJob<Vector_, Space_, max_der_> job(vector, space, cubature);
      trace_asm.assemble(job);
      return job.result();
    }
 
    template<typename Function_, typename Vector_, typename Space_, int max_der_>
    class TraceAssemblyErrorFunctionIntegralJob
    {
    public:
      typedef typename Vector_::DataType DataType;
      typedef typename Vector_::ValueType ValueType;
 
      static_assert(Intern::ErrCompatHelper<Function_, Vector_>::valid, "function and vector are incompatible");
 
      typedef typename DiscreteFunctionIntegral<Vector_, Space_>::Type FunctionIntegralType;
 
      static constexpr TrafoTags trafo_config = TrafoTags::none;
      static constexpr SpaceTags space_config = SpaceTags::value |
        (max_der_ >= 1 ? SpaceTags::grad : SpaceTags::none) |
        (max_der_ >= 2 ? SpaceTags::hess : SpaceTags::none);
      static constexpr TrafoTags facet_trafo_config = TrafoTags::img_point | TrafoTags::jac_det;
 
      class Task :
        public TraceAssemblyBasicTaskBase1<DataType, Space_, trafo_config, facet_trafo_config, space_config>
      {
      public:
        typedef TraceAssemblyBasicTaskBase1<DataType, Space_, trafo_config, facet_trafo_config, space_config> BaseClass;
 
        typedef typename BaseClass::AsmTraits AsmTraits;
 
        typedef Analytic::EvalTraits<DataType, Function_> AnalyticEvalTraits;
 
        // we do not support pairwise assembly
        static constexpr bool assemble_pairwise = false;
        static constexpr bool need_scatter = false;
        static constexpr bool need_combine = true;
 
        using typename BaseClass::TrafoEvalData;
        using typename BaseClass::TrafoFacetEvalData;
 
      protected:
        typename Function_::template Evaluator<AnalyticEvalTraits> func_eval;
        const Vector_& vector;
        typename BaseClass::CubatureRuleType cubature_rule;
        typename AsmTraits::template TLocalVector<ValueType> local_vector;
        typename Vector_::GatherAxpy gather_axpy;
        FunctionIntegralType loc_integral;
        FunctionIntegralType& job_integral;
 
      public:
        explicit Task(TraceAssemblyErrorFunctionIntegralJob& job) :
          BaseClass(job.space),
          func_eval(job.function),
          vector(job.vector),
          cubature_rule(Cubature::ctor_factory, job.cubature_factory),
          local_vector(),
          gather_axpy(vector),
          loc_integral(),
          job_integral(job.integral)
        {
        }
 
        void assemble()
        {
          // format local vector
          local_vector.format();
 
          // gather local vector data
          gather_axpy(local_vector, this->dof_mapping);
 
          // fetch number of local dofs
          const int num_loc_dofs = this->space_eval.get_num_local_dofs();
 
          // loop over all quadrature points and integrate
          for(int k(0); k < cubature_rule.get_num_points(); ++k)
          {
            // prepare point
            this->prepare_point(cubature_rule.get_point(k));
 
            // do the dirty work
            Intern::ErrFunIntJobHelper<max_der_>::work(loc_integral,
              cubature_rule.get_weight(k) * this->trafo_facet_data.jac_det, func_eval,
              this->trafo_facet_data.img_point, this->space_data, local_vector, num_loc_dofs);
 
          } // continue with next cubature point
        }
 
        void scatter()
        {
          // nothing to do here
        }
 
        void combine()
        {
          job_integral.push(loc_integral);
        }
      }; // class Task
 
    protected:
      const Function_& function;
      const Vector_& vector;
      const Space_& space;
      Cubature::DynamicFactory cubature_factory;
      FunctionIntegralType integral;
 
    public:
      explicit TraceAssemblyErrorFunctionIntegralJob(const Function_& function_, const Vector_& vector_, const Space_& space_, String cubature_):
        function(function_),
        vector(vector_),
        space(space_),
        cubature_factory(cubature_),
        integral()
      {
        integral.max_der = max_der_;
      }
 
      // \returns The integral of the discrete function
      FunctionIntegralType& result()
      {
        return integral;
      }
    }; // class TraceAssemblyErrorFunctionIntegralJob<...>
 
    template<int max_der_, typename Function_, typename Vector_, typename Trafo_, typename Space_>
    typename DiscreteFunctionIntegral<Vector_, Space_>::Type integrate_error_function(
      TraceAssembler<Trafo_>& trace_asm, const Function_& function,
      const Vector_& vector, const Space_& space, const String& cubature)
    {
      XASSERTM(trace_asm.get_trafo() == space.get_trafo(), "trace assembler and space have different trafos");
      XASSERTM(vector.size() == space.get_num_dofs(), "invalid coefficient vector length");
      TraceAssemblyErrorFunctionIntegralJob<Function_, Vector_, Space_, max_der_> job(function, vector, space, cubature);
      trace_asm.assemble(job);
      return job.result();
    }
 
    template<typename Derived_, typename VectorVelo_, typename VectorPres_, typename SpaceVelo_, typename SpacePres_,
      TrafoTags trafo_config_, TrafoTags facet_trafo_config_, SpaceTags space_velo_config_, SpaceTags space_pres_config_>
    class TraceAssemblyStokesVectorAnalysisTaskCRTP :
      public TraceAssemblyBasicTaskBase2<typename VectorVelo_::DataType, SpaceVelo_, SpacePres_, trafo_config_, facet_trafo_config_ | TrafoTags::jac_det, space_velo_config_, space_pres_config_>
    {
    public:
      typedef typename VectorVelo_::DataType DataType;
 
      typedef TraceAssemblyBasicTaskBase2<DataType, SpaceVelo_, SpacePres_, trafo_config_, facet_trafo_config_ | TrafoTags::jac_det, space_velo_config_, space_pres_config_> BaseClass;
      typedef typename BaseClass::AsmTraits AsmTraits;
 
      typedef typename VectorVelo_::ValueType VeloValueType;
      typedef typename VectorPres_::ValueType PresValueType;
 
      using typename BaseClass::TrafoEvalData;
      using typename BaseClass::TrafoFacetEvalData;
      typedef typename TrafoEvalData::EvalTraits TrafoEvalTraits;
 
      static constexpr int max_local_velo_dofs = AsmTraits::max_local_test_dofs;
      static constexpr int max_local_pres_dofs = AsmTraits::max_local_trial_dofs;
 
      typedef Space::BasisData<Space::StandardVectorEvalTraits<TrafoEvalTraits, max_local_velo_dofs, DataType>, space_velo_config_> VeloData;
      typedef Space::BasisData<Space::StandardScalarEvalTraits<TrafoEvalTraits, max_local_pres_dofs, DataType>, space_pres_config_> PresData;
 
    protected:
      const VectorVelo_& vector_velo;
      const VectorPres_& vector_pres;
      typename BaseClass::CubatureRuleType cubature_rule;
      typename AsmTraits::template TLocalTestVector<VeloValueType> local_vector_velo;
      typename AsmTraits::template TLocalTrialVector<PresValueType> local_vector_pres;
      typename VectorVelo_::GatherAxpy gather_axpy_velo;
      typename VectorPres_::GatherAxpy gather_axpy_pres;
 
    public:
      explicit TraceAssemblyStokesVectorAnalysisTaskCRTP(
        const VectorVelo_& vector_velo_, const VectorPres_& vector_pres_,
        const SpaceVelo_& space_velo_, const SpacePres_& space_pres_,
        const Cubature::DynamicFactory& cubature_factory_) :
        BaseClass(space_velo_, space_pres_),
        vector_velo(vector_velo_),
        vector_pres(vector_pres_),
        cubature_rule(Cubature::ctor_factory, cubature_factory_),
        gather_axpy_velo(vector_velo),
        gather_axpy_pres(vector_pres)
      {
      }
 
#ifdef DOXYGEN
      void eval(TrafoFacetEvalData& tau_f, TrafoEvalData& tau, const DataType& weight,
        const VeloData& velo, const PresData& pres);
#endif // DOXYGEN
 
      void assemble()
      {
        // gather local vectors
        local_vector_velo.format();
        local_vector_pres.format();
        gather_axpy_velo(local_vector_velo, this->test_dof_mapping);
        gather_axpy_pres(local_vector_pres, this->trial_dof_mapping);
 
        // fetch number of local dofs
        const int num_loc_dofs_velo = this->test_eval.get_num_local_dofs();
        const int num_loc_dofs_pres = this->trial_eval.get_num_local_dofs();
 
        VeloData velo_data;
        PresData pres_data;
 
        // loop over all quadrature points and integrate
        for(int k(0); k < cubature_rule.get_num_points(); ++k)
        {
          // prepare point
          static_cast<Derived_&>(*this).prepare_point(cubature_rule.get_point(k));
 
          velo_data.format();
          pres_data.format();
 
          // compute velocity data
          for(int i(0); i < num_loc_dofs_velo; ++i)
          {
            if constexpr (*(space_velo_config_ & SpaceTags::value))
              velo_data.value.axpy(this->test_data.phi[i].value, local_vector_velo(i));
            if constexpr (*(space_velo_config_ & SpaceTags::ref_value))
              velo_data.ref_value.axpy(this->test_data.phi[i].ref_value, local_vector_velo(i));
            if constexpr (*(space_velo_config_ & SpaceTags::grad))
              velo_data.grad.add_outer_product(local_vector_velo(i), this->test_data.phi[i].grad);
            if constexpr (*(space_velo_config_ & SpaceTags::ref_grad))
              velo_data.ref_grad.add_outer_product(local_vector_velo(i), this->test_data.phi[i].ref_grad);
            if constexpr (*(space_velo_config_ & SpaceTags::hess))
              velo_data.hess.add_vec_mat_outer_product(local_vector_velo(i), this->test_data.phi[i].hess);
            if constexpr (*(space_velo_config_ & SpaceTags::ref_hess))
              velo_data.ref_hess.add_vec_mat_outer_product(local_vector_velo(i), this->test_data.phi[i].ref_hess);
          }
 
          // compute velocity data
          for(int i(0); i < num_loc_dofs_pres; ++i)
          {
            if constexpr (*(space_pres_config_ & SpaceTags::value))
              pres_data.value += local_vector_pres(i) * this->trial_data.phi[i].value;
            if constexpr (*(space_pres_config_ & SpaceTags::ref_value))
              pres_data.ref_value += local_vector_pres(i) * this->trial_data.phi[i].ref_value;
            if constexpr (*(space_pres_config_ & SpaceTags::grad))
              pres_data.grad.axpy(local_vector_pres(i), this->trial_data.phi[i].grad);
            if constexpr (*(space_pres_config_ & SpaceTags::ref_grad))
              pres_data.ref_grad.axpy(local_vector_pres(i), this->trial_data.phi[i].ref_grad);
            if constexpr (*(space_pres_config_ & SpaceTags::hess))
              pres_data.hess.axpy(local_vector_pres(i), this->trial_data.phi[i].hess);
            if constexpr (*(space_pres_config_ & SpaceTags::ref_hess))
              pres_data.ref_hess.axpy(local_vector_pres(i), this->trial_data.phi[i].ref_hess);
          }
 
          // evaluate
          static_cast<Derived_&>(*this).eval(this->trafo_facet_data, this->trafo_data,
            cubature_rule.get_weight(k) * this->trafo_facet_data.jac_det, velo_data, pres_data);
        } // continue with next cubature point
      }
    }; // class TraceAssemblyStokesVectorAnalysisTaskCRTP<...>
 
    template<typename VectorVelo_, typename VectorPres_, typename SpaceVelo_, typename SpacePres_>
    class TraceAssemblyStokesBodyForceAssemblyJob
    {
    public:
      typedef typename VectorVelo_::DataType DataType;
      typedef typename VectorVelo_::ValueType VeloValueType;
      typedef typename VectorPres_::ValueType PresValueType;
 
      static constexpr TrafoTags trafo_config = TrafoTags::none;
      static constexpr SpaceTags space_velo_config = SpaceTags::value | SpaceTags::grad;
      static constexpr SpaceTags space_pres_config = SpaceTags::value;
      static constexpr TrafoTags facet_trafo_config = TrafoTags::jac_det | TrafoTags::normal;
 
      static constexpr int dim = SpaceVelo_::shape_dim;
 
      typedef Tiny::Matrix<DataType, 3, 2> RawForces;
 
      class Task :
        public TraceAssemblyStokesVectorAnalysisTaskCRTP<Task, VectorVelo_, VectorPres_, SpaceVelo_, SpacePres_, trafo_config, facet_trafo_config, space_velo_config, space_pres_config>
      {
      public:
        typedef TraceAssemblyStokesVectorAnalysisTaskCRTP<Task, VectorVelo_, VectorPres_, SpaceVelo_, SpacePres_, trafo_config, facet_trafo_config, space_velo_config, space_pres_config> BaseClass;
 
        typedef typename BaseClass::AsmTraits AsmTraits;
 
        // we do not support pairwise assembly
        static constexpr bool assemble_pairwise = false;
        static constexpr bool need_scatter = false;
        static constexpr bool need_combine = true;
 
        using typename BaseClass::TrafoEvalData;
        using typename BaseClass::TrafoFacetEvalData;
        using typename BaseClass::VeloData;
        using typename BaseClass::PresData;
 
      protected:
        RawForces& job_raw_forces;
        RawForces raw_forces;
 
      public:
        explicit Task(TraceAssemblyStokesBodyForceAssemblyJob& job) :
          BaseClass(job.vector_velo, job.vector_pres, job.space_velo, job.space_pres, job.cubature_factory),
          job_raw_forces(job.raw_forces),
          raw_forces()
        {
          raw_forces.format();
        }
 
        void eval(TrafoFacetEvalData& tau_f, TrafoEvalData& DOXY(tau), DataType weight, const VeloData& velo, const PresData& pres)
        {
          this->_eval(weight, tau_f.normal, velo.grad, pres.value);
        }
 
        void scatter()
        {
          // nothing to do here
        }
 
        void combine()
        {
          job_raw_forces += raw_forces;
        }
 
      protected:
        void _eval(const DataType omega, const Tiny::Vector<DataType, 2, 2>& n,
          const Tiny::Matrix<DataType, 2, 2, 2, 2>& grad_v, const DataType val_p)
        {
          raw_forces(0, 0) -= omega * (DataType(2) * grad_v(0,0) * n[0] + (grad_v(0, 1) + grad_v(1, 0)) * n[1]);
          raw_forces(1, 0) -= omega * (DataType(2) * grad_v(1,1) * n[1] + (grad_v(1, 0) + grad_v(0, 1)) * n[0]);
          raw_forces(0, 1) -= omega * val_p * n[0];
          raw_forces(1, 1) -= omega * val_p * n[1];
        }
 
        void _eval(const DataType omega, const Tiny::Vector<DataType, 3, 3>& n,
          const Tiny::Matrix<DataType, 3, 3, 3, 3>& grad_v, const DataType val_p)
        {
          raw_forces(0, 0) -= omega * (DataType(2) * grad_v(0,0) * n[0] + (grad_v(0, 1) + grad_v(1, 0)) * n[1] + (grad_v(0, 2) + grad_v(2, 0)) * n[2]);
          raw_forces(1, 0) -= omega * (DataType(2) * grad_v(1,1) * n[1] + (grad_v(1, 2) + grad_v(2, 1)) * n[2] + (grad_v(1, 0) + grad_v(0, 1)) * n[0]);
          raw_forces(2, 0) -= omega * (DataType(2) * grad_v(2,2) * n[2] + (grad_v(2, 0) + grad_v(0, 2)) * n[0] + (grad_v(2, 1) + grad_v(1, 2)) * n[1]);
          raw_forces(0, 1) -= omega * val_p * n[0];
          raw_forces(1, 1) -= omega * val_p * n[1];
          raw_forces(2, 1) -= omega * val_p * n[2];
        }
      }; // class Task
 
    protected:
      const VectorVelo_& vector_velo;
      const VectorPres_& vector_pres;
      const SpaceVelo_& space_velo;
      const SpacePres_& space_pres;
      Cubature::DynamicFactory cubature_factory;
      RawForces raw_forces;
 
    public:
      explicit TraceAssemblyStokesBodyForceAssemblyJob(const VectorVelo_& vector_velo_, const VectorPres_& vector_pres_,
        const SpaceVelo_& space_velo_, const SpacePres_& space_pres_, String cubature_):
        vector_velo(vector_velo_),
        vector_pres(vector_pres_),
        space_velo(space_velo_),
        space_pres(space_pres_),
        cubature_factory(cubature_),
        raw_forces()
      {
        raw_forces.format();
      }
 
      DataType drag(DataType nu) const
      {
        return nu * raw_forces(0, 0) - raw_forces(0, 1);
      }
 
      DataType lift(DataType nu) const
      {
        return nu * raw_forces(1, 0) - raw_forces(1, 1);
      }
 
      DataType side(DataType nu) const
      {
        return nu * raw_forces(2, 0) - raw_forces(2, 1);
      }
 
      void sync(const Dist::Comm& comm)
      {
        comm.allreduce(&raw_forces.v[0][0], &raw_forces.v[0][0], std::size_t(6), Dist::op_sum);
      }
    }; // class TraceAssemblyStokesBodyForceAssemblyJob<...>
  } // namespace Assembly
} // namespace FEAT