domain_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/domain_assembler.hpp>
#include <kernel/assembly/function_integral_info.hpp>
#include <kernel/analytic/function.hpp>
#include <kernel/lafem/dense_vector.hpp>
#include <kernel/lafem/dense_vector_blocked.hpp>
 
namespace FEAT
{
  namespace Assembly
  {
    template<typename Derived_, typename Vector_, typename Space_,
      TrafoTags trafo_config_, SpaceTags space_config_>
    class DomainAssemblyBasicVectorTaskCRTP
    {
    public:
      typedef typename Vector_::DataType DataType;
      typedef typename Vector_::ValueType ValueType;
 
      static constexpr bool need_scatter = true;
      static constexpr bool need_combine = false;
 
    protected:
      typedef Assembly::AsmTraits1<
        DataType,
        Space_,
        trafo_config_ | TrafoTags::jac_det,
        space_config_
      > AsmTraits;
 
      Vector_& vector;
      const Space_& space;
      const typename AsmTraits::TrafoType& trafo;
      typename AsmTraits::TrafoEvaluator trafo_eval;
      typename AsmTraits::SpaceEvaluator space_eval;
      typename AsmTraits::DofMapping dof_mapping;
      typename AsmTraits::CubatureRuleType cubature_rule;
      typename AsmTraits::TrafoEvalData trafo_data;
      typename AsmTraits::SpaceEvalData space_data;
      typename AsmTraits::template TLocalVector<ValueType> local_vector;
      typename Vector_::ScatterAxpy scatter_axpy;
      DataType scatter_alpha;
 
    public:
      explicit DomainAssemblyBasicVectorTaskCRTP(Vector_& vector_, const Space_& space_,
        const Cubature::DynamicFactory& cubature_factory, DataType alpha_) :
        vector(vector_),
        space(space_),
        trafo(space.get_trafo()),
        trafo_eval(trafo),
        space_eval(space),
        dof_mapping(space),
        cubature_rule(Cubature::ctor_factory, cubature_factory),
        scatter_axpy(vector),
        scatter_alpha(alpha_)
      {
      }
 
#ifdef DOXYGEN
      void set_point(typename AsmTraits::TrafoEvalData& tau);
 
      void eval(ValueType& val, const DataType& weight, const typename AsmTraits::SpaceBasisData& psi);
#endif // DOXYGEN
 
      void prepare(Index cell)
      {
        // prepare dof mapping
        dof_mapping.prepare(cell);
 
        // prepare trafo evaluator
        trafo_eval.prepare(cell);
 
        // prepare space evaluator
        space_eval.prepare(trafo_eval);
      }
 
      void assemble()
      {
        // format local vector
        local_vector.format();
 
        // fetch number of local dofs
        const int num_loc_dofs = space_eval.get_num_local_dofs();
 
        // 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));
 
          // compute basis function data
          space_eval(space_data, trafo_data);
 
          // evaluate operator
          static_cast<Derived_&>(*this).set_point(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),
              trafo_data.jac_det * cubature_rule.get_weight(k),
              space_data.phi[i]);
            // continue with next test function
          }
          // continue with next cubature point
        }
      }
 
      void scatter()
      {
        scatter_axpy(local_vector, dof_mapping, scatter_alpha);
      }
 
      void finish()
      {
        // finish evaluators
        space_eval.finish();
        trafo_eval.finish();
 
        // finish dof mapping
        dof_mapping.finish();
      }
 
      void combine()
      {
        // nothing to do here
      }
    }; // class DomainAssemblyBasicVectorTaskCRTP<...>
 
    template<typename Derived_, typename Matrix_, typename Space_,
      TrafoTags trafo_config_, SpaceTags space_config_>
    class DomainAssemblyBasicMatrixTaskCRTP1
    {
    public:
      typedef typename Matrix_::DataType DataType;
      typedef typename Matrix_::ValueType ValueType;
 
      static constexpr bool need_scatter = true;
      static constexpr bool need_combine = false;
 
    protected:
      typedef Assembly::AsmTraits1<
        DataType,
        Space_,
        trafo_config_ | TrafoTags::jac_det,
        space_config_
      > AsmTraits;
 
      Matrix_& matrix;
      const Space_& space;
      const typename AsmTraits::TrafoType& trafo;
      typename AsmTraits::TrafoEvaluator trafo_eval;
      typename AsmTraits::SpaceEvaluator space_eval;
      typename AsmTraits::DofMapping dof_mapping;
      typename AsmTraits::CubatureRuleType cubature_rule;
      typename AsmTraits::TrafoEvalData trafo_data;
      typename AsmTraits::SpaceEvalData space_data;
      typename AsmTraits::template TLocalMatrix<ValueType> local_matrix;
      typename Matrix_::ScatterAxpy scatter_axpy;
      DataType scatter_alpha;
 
    public:
      explicit DomainAssemblyBasicMatrixTaskCRTP1(Matrix_& matrix_, const Space_& space_,
        const Cubature::DynamicFactory& cubature_factory, DataType alpha_) :
        matrix(matrix_),
        space(space_),
        trafo(space.get_trafo()),
        trafo_eval(trafo),
        space_eval(space),
        dof_mapping(space),
        cubature_rule(Cubature::ctor_factory, cubature_factory),
        scatter_axpy(matrix),
        scatter_alpha(alpha_)
      {
      }
 
#ifdef DOXYGEN
      void set_point(typename AsmTraits::TrafoEvalData& tau);
 
      void eval(ValueType& val, const DataType& weight,
        const typename AsmTraits::TrialBasisData& phi, const typename AsmTraits::TestBasisData& psi);
#endif // DOXYGEN
 
      void prepare(Index cell)
      {
        // prepare dof mapping
        dof_mapping.prepare(cell);
 
        // prepare trafo evaluator
        trafo_eval.prepare(cell);
 
        // prepare space evaluator
        space_eval.prepare(trafo_eval);
      }
 
      void assemble()
      {
        // format local matrix
        local_matrix.format();
 
        // fetch number of local dofs
        const int num_loc_dofs = space_eval.get_num_local_dofs();
 
        // 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));
 
          // compute basis function data
          space_eval(space_data, trafo_data);
 
          // evaluate operator
          static_cast<Derived_&>(*this).set_point(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),
                trafo_data.jac_det * cubature_rule.get_weight(k),
                space_data.phi[j], space_data.phi[i]);
              // continue with next trial function
            }
            // continue with next test function
          }
          // continue with next cubature point
        }
      }
 
      void scatter()
      {
        // incorporate local matrix
        scatter_axpy(local_matrix, dof_mapping, dof_mapping, scatter_alpha);
      }
 
      void finish()
      {
        // finish evaluators
        space_eval.finish();
        trafo_eval.finish();
 
        // finish dof mapping
        dof_mapping.finish();
      }
 
      void combine()
      {
        // nothing to do here
      }
    }; // class DomainAssemblyBasicMatrixTaskCRTP1<...>
 
    template<
      typename Derived_,
      typename Matrix_,
      typename TestSpace_,
      typename TrialSpace_,
      TrafoTags trafo_config_,
      SpaceTags test_config_,
      SpaceTags trial_config_>
    class DomainAssemblyBasicMatrixTaskCRTP2
    {
    public:
      typedef typename Matrix_::DataType DataType;
      typedef typename Matrix_::ValueType ValueType;
 
      static constexpr bool need_scatter = true;
      static constexpr bool need_combine = false;
 
    protected:
      typedef Assembly::AsmTraits2<
        DataType,
        TestSpace_,
        TrialSpace_,
        trafo_config_ | TrafoTags::jac_det,
        test_config_,
        trial_config_
      > AsmTraits;
 
      Matrix_& matrix;
      const TestSpace_& test_space;
      const TrialSpace_& trial_space;
      const typename AsmTraits::TrafoType& trafo;
      typename AsmTraits::TrafoEvaluator trafo_eval;
      typename AsmTraits::TestEvaluator test_eval;
      typename AsmTraits::TrialEvaluator trial_eval;
      typename AsmTraits::TestDofMapping test_dof_mapping;
      typename AsmTraits::TrialDofMapping trial_dof_mapping;
      typename AsmTraits::CubatureRuleType cubature_rule;
      typename AsmTraits::TrafoEvalData trafo_data;
      typename AsmTraits::TestEvalData test_data;
      typename AsmTraits::TrialEvalData trial_data;
      typename AsmTraits::template TLocalMatrix<ValueType> local_matrix;
      typename Matrix_::ScatterAxpy scatter_axpy;
      DataType scatter_alpha;
 
    public:
      explicit DomainAssemblyBasicMatrixTaskCRTP2(Matrix_& matrix_,
        const TestSpace_& test_space_, const TrialSpace_& trial_space_,
        const Cubature::DynamicFactory& cubature_factory, DataType alpha_) :
        matrix(matrix_),
        test_space(test_space_),
        trial_space(trial_space_),
        trafo(test_space.get_trafo()),
        trafo_eval(trafo),
        test_eval(test_space),
        trial_eval(trial_space),
        test_dof_mapping(test_space),
        trial_dof_mapping(trial_space),
        cubature_rule(Cubature::ctor_factory, cubature_factory),
        scatter_axpy(matrix),
        scatter_alpha(alpha_)
      {
      }
 
#ifdef DOXYGEN
      void set_point(typename AsmTraits::TrafoEvalData& tau);
 
      void eval(ValueType& val, const DataType& weight,
        const typename AsmTraits::TrialBasisData& phi, const typename AsmTraits::TestBasisData& psi);
#endif // DOXYGEN
 
      void prepare(Index cell)
      {
        // prepare dof mapping
        test_dof_mapping.prepare(cell);
        trial_dof_mapping.prepare(cell);
 
        // prepare trafo evaluator
        trafo_eval.prepare(cell);
 
        // prepare space evaluators
        test_eval.prepare(trafo_eval);
        trial_eval.prepare(trafo_eval);
      }
 
      void assemble()
      {
        // format local matrix
        local_matrix.format();
 
        // fetch number of local dofs
        const int num_loc_test_dofs = test_eval.get_num_local_dofs();
        const int num_loc_trial_dofs = trial_eval.get_num_local_dofs();
 
        // 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));
 
          // compute basis function data
          test_eval(test_data, trafo_data);
          trial_eval(trial_data, trafo_data);
 
          // evaluate operator
          static_cast<Derived_&>(*this).set_point(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),
                trafo_data.jac_det * cubature_rule.get_weight(k),
                trial_data.phi[j], test_data.phi[i]);
              // continue with next trial function
            }
            // continue with next test function
          }
          // continue with next cubature point
        }
      }
 
      void scatter()
      {
        // incorporate local matrix
        scatter_axpy(local_matrix, test_dof_mapping, trial_dof_mapping, scatter_alpha);
      }
 
      void finish()
      {
        // finish evaluators
        trial_eval.finish();
        test_eval.finish();
        trafo_eval.finish();
 
        // finish dof mapping
        trial_dof_mapping.finish();
        test_dof_mapping.finish();
      }
 
      void combine()
      {
        // nothing to do here
      }
    }; // class DomainAssemblyBasicMatrixTaskCRTP2<...>
 
    template<typename LinearFunctional_, typename Vector_, typename Space_>
    class DomainAssemblyLinearFunctionalVectorJob
    {
    public:
      typedef typename Vector_::DataType DataType;
      typedef typename Vector_::ValueType ValueType;
 
      static constexpr TrafoTags trafo_config = LinearFunctional_::trafo_config;
      static constexpr SpaceTags test_config = LinearFunctional_::test_config;
 
      class Task :
        public DomainAssemblyBasicVectorTaskCRTP<Task, Vector_, Space_, trafo_config, test_config>
      {
      protected:
        typedef DomainAssemblyBasicVectorTaskCRTP<Task, Vector_, Space_, trafo_config, test_config> BaseClass;
 
        typedef typename BaseClass::AsmTraits AsmTraits;
 
        typename LinearFunctional_::template Evaluator<AsmTraits> func_eval;
 
      public:
        explicit Task(DomainAssemblyLinearFunctionalVectorJob& job) :
          BaseClass(job.vector, job.space, job.cubature_factory, job.alpha),
          func_eval(job.linear_functional)
        {
        }
 
        void prepare(Index cell)
        {
          BaseClass::prepare(cell);
          func_eval.prepare(this->trafo_eval);
        }
 
        void set_point(typename AsmTraits::TrafoEvalData& tau)
        {
          func_eval.set_point(tau);
        }
 
        void eval(ValueType& val, const DataType& weight, typename AsmTraits::TestBasisData& 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 DomainAssemblyLinearFunctionalVectorJob(const LinearFunctional_& linear_functional_,
        Vector_& vector_, const Space_& space_, const String& cubature_, DataType alpha_ = DataType(1)) :
        linear_functional(linear_functional_),
        vector(vector_),
        space(space_),
        cubature_factory(cubature_),
        alpha(alpha_)
      {
      }
    }; // class DomainAssemblyLinearFunctionalVectorJob<...>
 
    template<typename Trafo_, typename Vector_, typename LinFunc_, typename Space_>
    void assemble_linear_functional_vector(DomainAssembler<Trafo_>& dom_asm,
      Vector_& vector, const LinFunc_& linear_functional, const Space_& space,
      const String& cubature, const typename Vector_::DataType alpha = typename Vector_::DataType(1))
    {
      XASSERTM(dom_asm.get_trafo() == space.get_trafo(), "domain assembler and space have different trafos");
      XASSERTM(vector.size() == space.get_num_dofs(), "invalid vector length");
 
      DomainAssemblyLinearFunctionalVectorJob<LinFunc_, Vector_, Space_> job(
        linear_functional, vector, space, cubature, alpha);
      dom_asm.assemble(job);
    }
 
    template<typename Function_, typename Vector_, typename Space_>
    class DomainAssemblyForceFunctionalVectorJob
    {
    public:
      typedef typename Vector_::DataType DataType;
      typedef typename Vector_::ValueType ValueType;
 
      static constexpr TrafoTags trafo_config = TrafoTags::img_point | TrafoTags::jac_det;
      static constexpr SpaceTags test_config = SpaceTags::value;
 
      class Task :
        public DomainAssemblyBasicVectorTaskCRTP<Task, Vector_, Space_, trafo_config, test_config>
      {
      protected:
        typedef DomainAssemblyBasicVectorTaskCRTP<Task, Vector_, Space_, trafo_config, test_config> BaseClass;
 
        typedef typename BaseClass::AsmTraits AsmTraits;
 
        typedef Analytic::EvalTraits<DataType, Function_> AnalyticEvalTraits;
        typename Function_::template Evaluator<AnalyticEvalTraits> func_eval;
        typename AnalyticEvalTraits::ValueType func_value;
 
      public:
        explicit Task(DomainAssemblyForceFunctionalVectorJob& job) :
          BaseClass(job.vector, job.space, job.cubature_factory, job.alpha),
          func_eval(job.function)
        {
        }
 
        void set_point(typename AsmTraits::TrafoEvalData& tau)
        {
          func_value = func_eval.value(tau.img_point);
        }
 
        void eval(ValueType& val, const DataType& weight, typename AsmTraits::TestBasisData& psi)
        {
          Tiny::axpy(val, func_value, weight * psi.value);
        }
      }; // class Task
 
    protected:
      const Function_& function;
      Vector_& vector;
      const Space_& space;
      Cubature::DynamicFactory cubature_factory;
      DataType alpha;
 
    public:
      explicit DomainAssemblyForceFunctionalVectorJob(const Function_& function_,
        Vector_& vector_, const Space_& space_, const String& cubature_, DataType alpha_ = DataType(1)) :
        function(function_),
        vector(vector_),
        space(space_),
        cubature_factory(cubature_),
        alpha(alpha_)
      {
      }
    }; // class DomainAssemblyForceFunctionalVectorJob<...>
 
    template<typename Trafo_, typename Vector_, typename Function_, typename Space_>
    void assemble_force_function_vector(DomainAssembler<Trafo_>& dom_asm,
      Vector_& vector, const Function_& function, const Space_& space,
      const String& cubature, const typename Vector_::DataType alpha = typename Vector_::DataType(1))
    {
      XASSERTM(dom_asm.get_trafo() == space.get_trafo(), "domain assembler and space have different trafos");
      XASSERTM(vector.size() == space.get_num_dofs(), "invalid vector length");
 
      DomainAssemblyForceFunctionalVectorJob<Function_, Vector_, Space_> job(
        function, vector, space, cubature, alpha);
      dom_asm.assemble(job);
    }
 
    template<typename BilinearOperator_, typename Matrix_, typename Space_>
    class DomainAssemblyBilinearOperatorMatrixJob1
    {
    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 DomainAssemblyBasicMatrixTaskCRTP1<Task, Matrix_, Space_, trafo_config, space_config>
      {
      protected:
        typedef DomainAssemblyBasicMatrixTaskCRTP1<Task, Matrix_, Space_, trafo_config, space_config> BaseClass;
 
        typedef typename BaseClass::AsmTraits AsmTraits;
 
        typename BilinearOperator_::template Evaluator<AsmTraits> oper_eval;
 
      public:
        explicit Task(DomainAssemblyBilinearOperatorMatrixJob1& job) :
          BaseClass(job.matrix, job.space, job.cubature_factory, job.alpha),
          oper_eval(job.bilinear_operator)
        {
        }
 
        void prepare(Index cell)
        {
          BaseClass::prepare(cell);
          oper_eval.prepare(this->trafo_eval);
        }
 
        void set_point(typename AsmTraits::TrafoEvalData& tau)
        {
          oper_eval.set_point(tau);
        }
 
        void eval(ValueType& val, const DataType& weight,
          typename AsmTraits::TrialBasisData& phi, typename AsmTraits::TestBasisData& psi)
        {
          val += weight * oper_eval.eval(phi, psi);
        }
 
        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 DomainAssemblyBilinearOperatorMatrixJob1(const BilinearOperator_& bilinear_operator_,
        Matrix_& matrix_, const Space_& space_, const String& cubature_, DataType alpha_ = DataType(1)) :
        bilinear_operator(bilinear_operator_),
        matrix(matrix_),
        space(space_),
        cubature_factory(cubature_),
        alpha(alpha_)
      {
      }
    }; // class DomainAssemblyBilinearOperatorMatrixJob1<...>
 
    template<typename Trafo_, typename Matrix_, typename BilOp_, typename Space_>
    void assemble_bilinear_operator_matrix_1(DomainAssembler<Trafo_>& dom_asm, Matrix_& matrix,
      const BilOp_& bilinear_operator, const Space_& space, const String& cubature,
      const typename Matrix_::DataType alpha = typename Matrix_::DataType(1))
    {
      XASSERTM(dom_asm.get_trafo() == space.get_trafo(), "domain 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");
 
      DomainAssemblyBilinearOperatorMatrixJob1<BilOp_, Matrix_, Space_> job(
        bilinear_operator, matrix, space, cubature, alpha);
      dom_asm.assemble(job);
    }
 
    template<typename BilinearOperator_, typename Matrix_, typename TestSpace_, typename TrialSpace_>
    class DomainAssemblyBilinearOperatorMatrixJob2
    {
    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 DomainAssemblyBasicMatrixTaskCRTP2<Task, Matrix_, TestSpace_, TrialSpace_, trafo_config, test_config, trial_config>
      {
      protected:
        typedef DomainAssemblyBasicMatrixTaskCRTP2<Task, Matrix_, TestSpace_, TrialSpace_, trafo_config, test_config, trial_config> BaseClass;
 
        typedef typename BaseClass::AsmTraits AsmTraits;
 
        typename BilinearOperator_::template Evaluator<AsmTraits> oper_eval;
 
      public:
        explicit Task(DomainAssemblyBilinearOperatorMatrixJob2& job) :
          BaseClass(job.matrix, job.test_space, job.trial_space, job.cubature_factory, job.alpha),
          oper_eval(job.bilinear_operator)
        {
        }
 
        void prepare(Index cell)
        {
          BaseClass::prepare(cell);
          oper_eval.prepare(this->trafo_eval);
        }
 
        void set_point(typename AsmTraits::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 DomainAssemblyBilinearOperatorMatrixJob2(const BilinearOperator_& bilinear_operator_,
        Matrix_& matrix_, const TestSpace_& test_space_, const TrialSpace_& trial_space_,
        const 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 DomainAssemblyBilinearOperatorMatrixJob2<...>
 
    template<typename Trafo_, typename Matrix_, typename BilOp_, typename TestSpace_, typename TrialSpace_>
    void assemble_bilinear_operator_matrix_2(DomainAssembler<Trafo_>& dom_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(dom_asm.get_trafo() == test_space.get_trafo(), "domain assembler and test space have different trafos");
      XASSERTM(dom_asm.get_trafo() == trial_space.get_trafo(), "domain 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");
 
      DomainAssemblyBilinearOperatorMatrixJob2<BilOp_, Matrix_, TestSpace_, TrialSpace_> job(
        bilinear_operator, matrix, test_space, trial_space, cubature, alpha);
      dom_asm.assemble(job);
    }
 
    template<typename BilinearOperator_, typename Vector_, typename VectorSol_, typename Space_>
    class DomainAssemblyBilinearOperatorApplyVectorJob1
    {
    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 DomainAssemblyBasicVectorTaskCRTP<Task, Vector_, Space_, trafo_config, space_config>
      {
      protected:
        typedef DomainAssemblyBasicVectorTaskCRTP<Task, Vector_, Space_, trafo_config, 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(DomainAssemblyBilinearOperatorApplyVectorJob1& 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 cell)
        {
          BaseClass::prepare(cell);
          oper_eval.prepare(this->trafo_eval);
          local_vec_sol.format();
          sol_gather(local_vec_sol, this->dof_mapping);
        }
 
        void set_point(typename AsmTraits::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)
          {
            // compute trafo data
            this->trafo_eval(this->trafo_data, this->cubature_rule.get_point(k));
 
            // compute basis function data
            this->space_eval(this->space_data, this->trafo_data);
 
            // evaluate operator
            this->set_point(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_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 DomainAssemblyBilinearOperatorApplyVectorJob1(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 DomainAssemblyBilinearOperatorApplyVectorJob1<...>
 
    template<typename Trafo_, typename Vector_, typename VectorSol_, typename BilOp_, typename Space_>
    void assemble_bilinear_operator_apply_vector_1(DomainAssembler<Trafo_>& dom_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(dom_asm.get_trafo() == space.get_trafo(), "domain assembler and space have different trafos");
      XASSERTM(vector.size() == space.get_num_dofs(), "invalid vector size");
 
      DomainAssemblyBilinearOperatorApplyVectorJob1<BilOp_, Vector_, VectorSol_, Space_> job(
        bilinear_operator, vector, vec_sol, space, cubature, alpha);
      dom_asm.assemble(job);
    }
 
    template<typename DataType_, typename Function_, typename Trafo_, int max_der_>
    class DomainAssemblyAnalyticFunctionIntegralJob
    {
    public:
      static_assert(Function_::can_value, "function cannot compute values");
      static_assert(Function_::can_grad || (max_der_ < 1), "function gradients are required but not available");
      static_assert(Function_::can_hess || (max_der_ < 2), "function hessians are required but not available");
 
      typedef DataType_ DataType;
 
      typedef Analytic::EvalTraits<DataType, Function_> AnalyticEvalTraits;
 
      typedef typename AnalyticFunctionIntegral<DataType, Function_>::Type FunctionIntegralType;
 
    public:
      class Task
      {
      public:
        static constexpr bool need_scatter = false;
        static constexpr bool need_combine = true;
 
      protected:
        static constexpr TrafoTags trafo_config = TrafoTags::img_point | TrafoTags::jac_det;
        typedef typename Trafo_::template Evaluator<typename Trafo_::ShapeType, DataType>::Type TrafoEvaluator;
        TrafoEvaluator trafo_eval;
        typename TrafoEvaluator::template ConfigTraits<trafo_config>::EvalDataType 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(DomainAssemblyAnalyticFunctionIntegralJob& 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 cell)
        {
          trafo_eval.prepare(cell);
        }
 
        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 DomainAssemblyAnalyticFunctionIntegralJob(const Function_& function, const Trafo_& trafo, const String& cubature) :
        _trafo(trafo),
        _function(function),
        _cubature_factory(cubature),
        _integral()
      {
        _integral.max_der = max_der_;
      }
 
      FunctionIntegralType& result()
      {
        return _integral;
      }
    }; // class DomainAssemblyAnalyticFunctionIntegralJob<...>
 
    template<int max_der_, typename DataType_, typename Trafo_, typename Function_>
    typename AnalyticFunctionIntegral<DataType_, Function_>::Type integrate_analytic_function(
      DomainAssembler<Trafo_>& dom_asm, const Function_& function, const String& cubature)
    {
      DomainAssemblyAnalyticFunctionIntegralJob<DataType_, Function_, Trafo_, max_der_> job(function, dom_asm.get_trafo(), cubature);
      dom_asm.assemble(job);
      return job.result();
    }
 
    template<typename Vector_, typename Space_, int max_der_>
    class DomainAssemblyDiscreteFunctionIntegralJob
    {
    public:
      typedef typename Vector_::DataType DataType;
      typedef typename Vector_::ValueType ValueType;
 
      typedef typename DiscreteFunctionIntegral<Vector_, Space_>::Type FunctionIntegralType;
 
    public:
      class Task
      {
      public:
        static constexpr bool need_scatter = false;
        static constexpr bool need_combine = true;
 
      protected:
        static constexpr SpaceTags space_tags = SpaceTags::value |
          (max_der_ >= 1 ? SpaceTags::grad : SpaceTags::none) |
          (max_der_ >= 2 ? SpaceTags::hess : SpaceTags::none);
 
        typedef Assembly::AsmTraits1<DataType, Space_, TrafoTags::jac_det, space_tags> AsmTraits;
        const Vector_& vector;
        const Space_& space;
        const typename AsmTraits::TrafoType& trafo;
        typename AsmTraits::TrafoEvaluator trafo_eval;
        typename AsmTraits::SpaceEvaluator space_eval;
        typename AsmTraits::DofMapping dof_mapping;
        typename AsmTraits::CubatureRuleType cubature_rule;
        typename AsmTraits::TrafoEvalData trafo_data;
        typename AsmTraits::SpaceEvalData space_data;
        typename AsmTraits::template TLocalVector<ValueType> local_vector;
        typename Vector_::GatherAxpy gather_axpy;
        FunctionIntegralType loc_integral;
        FunctionIntegralType& job_integral;
 
      public:
        explicit Task(DomainAssemblyDiscreteFunctionIntegralJob& job) :
          vector(job._vector),
          space(job._space),
          trafo(space.get_trafo()),
          trafo_eval(trafo),
          space_eval(space),
          dof_mapping(space),
          cubature_rule(Cubature::ctor_factory, job._cubature_factory),
          trafo_data(),
          space_data(),
          local_vector(),
          gather_axpy(vector),
          loc_integral(),
          job_integral(job._integral)
        {
        }
 
        void prepare(Index cell)
        {
          // prepare dof mapping
          dof_mapping.prepare(cell);
 
          // prepare trafo evaluator
          trafo_eval.prepare(cell);
 
          // prepare space evaluator
          space_eval.prepare(trafo_eval);
        }
 
        void assemble()
        {
          // format local vector
          local_vector.format();
 
          // gather local vector data
          gather_axpy(local_vector, dof_mapping);
 
          // fetch number of local dofs
          const int num_loc_dofs = space_eval.get_num_local_dofs();
 
          // 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));
 
            // compute basis function data
            space_eval(space_data, trafo_data);
 
            // do the dirty work
            Intern::DiscFunIntJobHelper<max_der_>::work(loc_integral,
              cubature_rule.get_weight(k) * trafo_data.jac_det, space_data, local_vector, num_loc_dofs);
          }
        }
 
        void scatter()
        {
          // nothing to do here
        }
 
        void finish()
        {
          trafo_eval.finish();
        }
 
        void combine()
        {
          job_integral.push(loc_integral);
        }
      }; // class Task
 
    protected:
      const Vector_& _vector;
      const Space_& _space;
      Cubature::DynamicFactory _cubature_factory;
      FunctionIntegralType _integral;
 
    public:
      explicit DomainAssemblyDiscreteFunctionIntegralJob(const Vector_& vector, const Space_& space, const 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 DomainAssemblyDiscreteFunctionIntegralJob<...>
 
    template<int max_der_, typename Vector_, typename Trafo_, typename Space_>
    typename DiscreteFunctionIntegral<Vector_, Space_>::Type integrate_discrete_function(
       DomainAssembler<Trafo_>& dom_asm, const Vector_& vector, const Space_& space, const String& cubature)
    {
      XASSERTM(dom_asm.get_trafo() == space.get_trafo(), "domain assembler and space have different trafos");
      XASSERTM(vector.size() == space.get_num_dofs(), "invalid coefficient vector length");
      DomainAssemblyDiscreteFunctionIntegralJob<Vector_, Space_, max_der_> job(vector, space, cubature);
      dom_asm.assemble(job);
      return job.result();
    }
 
    template<typename Function_, typename Vector_, typename Space_, int max_der_>
    class DomainAssemblyErrorFunctionIntegralJob
    {
    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;
 
      typedef Analytic::EvalTraits<DataType, Function_> AnalyticEvalTraits;
 
    public:
      class Task
      {
      public:
        static constexpr bool need_scatter = false;
        static constexpr bool need_combine = true;
 
      protected:
        static constexpr TrafoTags trafo_config = TrafoTags::img_point | TrafoTags::jac_det;
        static constexpr SpaceTags space_tags = SpaceTags::value |
          (max_der_ >= 1 ? SpaceTags::grad : SpaceTags::none) |
          (max_der_ >= 2 ? SpaceTags::hess : SpaceTags::none);
 
        typedef Assembly::AsmTraits1<DataType, Space_, trafo_config, space_tags> AsmTraits;
        const Vector_& vector;
        const Space_& space;
        const typename AsmTraits::TrafoType& trafo;
        typename AsmTraits::TrafoEvaluator trafo_eval;
        typename AsmTraits::SpaceEvaluator space_eval;
        typename AsmTraits::DofMapping dof_mapping;
        typename AsmTraits::CubatureRuleType cubature_rule;
        typename AsmTraits::TrafoEvalData trafo_data;
        typename AsmTraits::SpaceEvalData space_data;
        typename AsmTraits::template TLocalVector<ValueType> local_vector;
        typename Vector_::GatherAxpy gather_axpy;
        typename Function_::template Evaluator<AnalyticEvalTraits> func_eval;
        FunctionIntegralType loc_integral;
        FunctionIntegralType& job_integral;
 
      public:
        explicit Task(DomainAssemblyErrorFunctionIntegralJob& job) :
          vector(job._vector),
          space(job._space),
          trafo(space.get_trafo()),
          trafo_eval(trafo),
          space_eval(space),
          dof_mapping(space),
          cubature_rule(Cubature::ctor_factory, job._cubature_factory),
          trafo_data(),
          space_data(),
          local_vector(),
          gather_axpy(vector),
          func_eval(job._function),
          loc_integral(),
          job_integral(job._integral)
        {
        }
 
        void prepare(Index cell)
        {
          // prepare dof mapping
          dof_mapping.prepare(cell);
 
          // prepare trafo evaluator
          trafo_eval.prepare(cell);
 
          // prepare space evaluator
          space_eval.prepare(trafo_eval);
        }
 
        void assemble()
        {
          // format local vector
          local_vector.format();
 
          // gather local vector data
          gather_axpy(local_vector, dof_mapping);
 
          // fetch number of local dofs
          const int num_loc_dofs = space_eval.get_num_local_dofs();
 
          // 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));
 
            // compute basis function data
            space_eval(space_data, trafo_data);
 
            // do the dirty work
            Intern::ErrFunIntJobHelper<max_der_>::work(loc_integral,
              cubature_rule.get_weight(k) * trafo_data.jac_det, func_eval, trafo_data.img_point,
              space_data, local_vector, num_loc_dofs);
          }
        }
 
        void scatter()
        {
          // nothing to do here
        }
 
        void finish()
        {
          trafo_eval.finish();
        }
 
        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 DomainAssemblyErrorFunctionIntegralJob(const Function_& function, const Vector_& vector, const Space_& space, const 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 DomainAssemblyErrorFunctionIntegralJob<...>
 
    template<int max_der_, typename Function_, typename Vector_, typename Trafo_, typename Space_>
    typename DiscreteFunctionIntegral<Vector_, Space_>::Type integrate_error_function(
       DomainAssembler<Trafo_>& dom_asm, const Function_& function,
       const Vector_& vector, const Space_& space, const String& cubature)
    {
      XASSERTM(dom_asm.get_trafo() == space.get_trafo(), "domain assembler and space have different trafos");
      XASSERTM(vector.size() == space.get_num_dofs(), "invalid coefficient vector length");
      DomainAssemblyErrorFunctionIntegralJob<Function_, Vector_, Space_, max_der_> job(function, vector, space, cubature);
      dom_asm.assemble(job);
      return job.result();
    }
 
    namespace Intern
    {
      template<typename DT_, typename IT_, int max_der_>
      struct OutVectorHelper
      {
        typedef LAFEM::DenseVectorBlocked<DT_, IT_, max_der_+1> VectorType;
      };
 
      template<typename DT_, typename IT_>
      struct OutVectorHelper<DT_, IT_, 0>
      {
        typedef LAFEM::DenseVector<DT_, IT_> VectorType;
      };
    }
 
    template<typename Function_, typename Vector_, typename Space_, int max_der_>
    class CellErrorFunctionIntegralJob
    {
    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;
 
      typedef Analytic::EvalTraits<DataType, Function_> AnalyticEvalTraits;
 
      typedef typename Intern::OutVectorHelper<DataType, Index, max_der_>::VectorType OutVectorType;
 
      typedef FunctionCellIntegralInfo<FunctionIntegralType, OutVectorType> FunctionCellIntegralType;
 
    public:
      class Task
      {
      public:
        static constexpr bool need_scatter = true;
        static constexpr bool need_combine = true;
 
      protected:
        static constexpr TrafoTags trafo_config = TrafoTags::img_point | TrafoTags::jac_det;
        static constexpr SpaceTags space_tags = SpaceTags::value |
          (max_der_ >= 1 ? SpaceTags::grad : SpaceTags::none) |
          (max_der_ >= 2 ? SpaceTags::hess : SpaceTags::none);
 
        typedef Assembly::AsmTraits1<DataType, Space_, trafo_config, space_tags> AsmTraits;
        const Vector_& vector;
        const Space_& space;
        OutVectorType& out_vector;
        const typename AsmTraits::TrafoType& trafo;
        typename AsmTraits::TrafoEvaluator trafo_eval;
        typename AsmTraits::SpaceEvaluator space_eval;
        typename AsmTraits::DofMapping dof_mapping;
        typename AsmTraits::CubatureRuleType cubature_rule;
        typename AsmTraits::TrafoEvalData trafo_data;
        typename AsmTraits::SpaceEvalData space_data;
        typename AsmTraits::template TLocalVector<ValueType> local_vector;
        typename Vector_::GatherAxpy gather_axpy;
        typename Function_::template Evaluator<AnalyticEvalTraits> func_eval;
        FunctionIntegralType loc_integral;
        FunctionIntegralType& job_integral;
        FunctionIntegralType cell_integral;
 
      public:
        explicit Task(CellErrorFunctionIntegralJob& job) :
          vector(job._vector),
          space(job._space),
          out_vector(job._out_vec),
          trafo(space.get_trafo()),
          trafo_eval(trafo),
          space_eval(space),
          dof_mapping(space),
          cubature_rule(Cubature::ctor_factory, job._cubature_factory),
          trafo_data(),
          space_data(),
          local_vector(),
          gather_axpy(vector),
          func_eval(job._function),
          loc_integral(),
          job_integral(job._integral),
          cell_integral()
        {
        }
 
        void prepare(Index cell)
        {
          // prepare dof mapping
          dof_mapping.prepare(cell);
 
          // prepare trafo evaluator
          trafo_eval.prepare(cell);
 
          // prepare space evaluator
          space_eval.prepare(trafo_eval);
        }
 
        void assemble()
        {
          // format local vector
          local_vector.format();
 
          // gather local vector data
          gather_axpy(local_vector, dof_mapping);
 
          // fetch number of local dofs
          const int num_loc_dofs = space_eval.get_num_local_dofs();
 
          // format cell local integral value
          cell_integral.format();
 
          // 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));
 
            // compute basis function data
            space_eval(space_data, trafo_data);
 
            // do the dirty work
            Intern::ErrFunIntJobHelper<max_der_>::work(cell_integral,
              cubature_rule.get_weight(k) * trafo_data.jac_det, func_eval, trafo_data.img_point,
              space_data, local_vector, num_loc_dofs);
          }
          // add our cell value
          loc_integral.push(cell_integral);
        }
 
        void scatter()
        {
          // get cell index
          Index cell = dof_mapping.get_current_cell_index();
          if constexpr(max_der_ == 0)
          {
            out_vector(cell, cell_integral.norm_h0_sqr);
          }
          else
          {
            // construct tiny vector
            Tiny::Vector<DataType, max_der_+1> loc_vec;
            loc_vec[0] = cell_integral.norm_h0_sqr;
            loc_vec[1] = cell_integral.norm_h1_sqr;
            if constexpr(max_der_ > 1)
              loc_vec[2] = cell_integral.norm_h2_sqr;
            //scatter the tiny vector
            out_vector(cell, loc_vec);
          }
        }
 
        void finish()
        {
          trafo_eval.finish();
        }
 
        void combine()
        {
          job_integral.push(loc_integral);
        }
      }; // class Task
 
    protected:
      const Function_& _function;
      const Vector_& _vector;
      const Space_& _space;
      OutVectorType _out_vec;
      Cubature::DynamicFactory _cubature_factory;
      FunctionIntegralType _integral;
 
    public:
      explicit CellErrorFunctionIntegralJob(const Function_& function, const Vector_& vector, const Space_& space, const String& cubature) :
        _function(function),
        _vector(vector),
        _space(space),
        _cubature_factory(cubature),
        _integral()
      {
        _integral.max_der = max_der_;
        _out_vec = OutVectorType(_space.get_mesh().get_num_elements());
        _out_vec.format();
      }
 
      FunctionCellIntegralType result()
      {
        FunctionCellIntegralType output(_integral, std::move(_out_vec));
        _out_vec = OutVectorType(_space.get_mesh().get_num_elements());
        _out_vec.format();
        return output;
      }
    }; // class CellErrorFunctionIntegralJob<...>
  } // namespace Assembly
} // namespace FEAT