burgers_assembly_job.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/base.hpp>
#include <kernel/assembly/asm_traits.hpp>
#include <kernel/cubature/dynamic_factory.hpp>
#include <kernel/lafem/dense_vector.hpp>
#include <kernel/lafem/sparse_matrix_bcsr.hpp>
#include <kernel/lafem/dense_vector_blocked.hpp>
#include <kernel/global/vector.hpp>
 
namespace FEAT
{
  namespace Assembly
  {
    template<typename DataType_, typename Space_, typename ConvVector_>
    class BurgersAssemblyJobBase
    {
    public:
      typedef DataType_ DataType;
 
      typedef Space_ SpaceType;
 
      typedef ConvVector_ ConvVectorType;
 
      bool deformation;
 
      DataType_ nu;
 
      DataType_ theta;
 
      DataType_ beta;
 
      DataType_ frechet_beta;
 
      DataType_ sd_delta;
 
      DataType_ sd_nu;
 
      DataType_ sd_v_norm;
 
      const ConvVector_& convection_vector;
 
      const Space_& space;
 
      Cubature::DynamicFactory cubature_factory;
 
      DataType_ alpha;
 
    public:
      explicit BurgersAssemblyJobBase(const ConvVector_& conv_vector, const Space_& space_, const String& cubature) :
        deformation(false),
        nu(DataType_(0)),
        theta(DataType_(0)),
        beta(DataType_(0)),
        frechet_beta(DataType_(0)),
        sd_delta(DataType_(0)),
        sd_nu(DataType_(0)),
        sd_v_norm(DataType_(0)),
        convection_vector(conv_vector),
        space(space_),
        cubature_factory(cubature),
        alpha(DataType_(1))
      {
      }
 
      template<typename IndexType_, int conv_dim_>
      static DataType calc_sd_v_norm(const LAFEM::DenseVectorBlocked<DataType_, IndexType_, conv_dim_>& convect)
      {
        const auto* vals = convect.elements();
        DataType_ r = DataType_(0);
        for(Index i(0); i < convect.size(); ++i)
          r = Math::max(r, vals[i].norm_euclid());
        return r;
      }
 
      template<typename LocalVector_, typename Mirror_>
      static DataType calc_sd_v_norm(const Global::Vector<LocalVector_, Mirror_>& convect)
      {
        const auto* gate = convect.get_gate();
        if(gate != nullptr)
          return gate->max(calc_sd_v_norm(convect.local()));
        else
          return calc_sd_v_norm(convect.local());
      }
 
      template<typename VectorType_>
      void set_sd_v_norm(const VectorType_& convect)
      {
        this->sd_v_norm = calc_sd_v_norm(convect);
      }
    }; // class BurgersAssemblyJobBase<...>
 
    template<typename Job_, typename DataType_>
    class BurgersAssemblyTaskBase
    {
    public:
      static constexpr bool need_scatter = true;
      static constexpr bool need_combine = false;
 
      typedef DataType_ DataType;
 
      typedef typename Job_::SpaceType SpaceType;
 
      typedef typename Job_::ConvVectorType ConvVectorType;
 
      static constexpr int shape_dim = SpaceType::shape_dim;
 
      static constexpr int conv_dim = ConvVectorType::BlockSize;
 
      typedef AsmTraits1<DataType, SpaceType, TrafoTags::jac_det, SpaceTags::value|SpaceTags::grad> AsmTraits;
 
      const Job_& job;
 
      const DataType tol_eps;
 
      const bool deformation;
 
      const DataType nu;
 
      const DataType theta;
 
      const DataType beta;
 
      const DataType frechet_beta;
 
      const DataType sd_delta;
 
      const DataType sd_nu;
 
      const DataType sd_v_norm;
 
      const DataType alpha;
 
      const bool need_diff, need_conv, need_conv_frechet, need_reac, need_streamdiff;
 
      bool need_loc_v, need_loc_grad_v, need_mean_v, need_local_h;
 
      const SpaceType& 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 ConvVectorType::GatherAxpy gather_conv;
 
      static constexpr int max_local_dofs = AsmTraits::max_local_test_dofs;
 
      int num_local_dofs;
 
      Tiny::Vector<Tiny::Vector<DataType, conv_dim>, max_local_dofs> local_conv_dofs;
 
      Tiny::Vector<DataType, conv_dim> loc_v, mean_v;
 
      Tiny::Matrix<DataType, conv_dim, conv_dim> loc_grad_v;
 
      Tiny::Vector<DataType, max_local_dofs> streamdiff_coeffs;
 
      Tiny::Vector<DataType, shape_dim> barycenter;
 
      DataType local_h;
 
      DataType local_delta;
 
    public:
      explicit BurgersAssemblyTaskBase(const Job_& job_) :
        job(job_),
        tol_eps(Math::sqrt(Math::eps<DataType>())),
        deformation(job_.deformation),
        nu(job_.nu),
        theta(job_.theta),
        beta(job_.beta),
        frechet_beta(job_.frechet_beta),
        sd_delta(job_.sd_delta),
        sd_nu(job_.sd_nu),
        sd_v_norm(job_.sd_v_norm),
        alpha(job_.alpha),
        need_diff(Math::abs(nu) > DataType(0)),
        need_conv(Math::abs(beta) > DataType(0)),
        need_conv_frechet(Math::abs(frechet_beta) > DataType(0)),
        need_reac(Math::abs(theta) > DataType(0)),
        need_streamdiff((Math::abs(sd_delta) > DataType(0)) && (sd_v_norm > tol_eps)),
        need_loc_v(need_conv),
        need_loc_grad_v(need_conv_frechet),
        need_mean_v(need_streamdiff),
        need_local_h(need_streamdiff),
        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(),
        gather_conv(job_.convection_vector),
        num_local_dofs(0),
        local_conv_dofs(),
        loc_v(),
        mean_v(),
        loc_grad_v(),
        streamdiff_coeffs(),
        barycenter(),
        local_h(DataType(0)),
        local_delta(DataType(0))
      {
        // compute reference element barycenter
        for(int i(0); i < shape_dim; ++i)
          barycenter[i] = Shape::ReferenceCell<typename SpaceType::ShapeType>::template centre<DataType>(i);
      }
 
      void prepare(const Index cell)
      {
        // prepare dof mapping
        dof_mapping.prepare(cell);
 
        // prepare trafo evaluator
        trafo_eval.prepare(cell);
 
        // prepare space evaluator
        space_eval.prepare(trafo_eval);
 
        // fetch number of local dofs
        num_local_dofs = space_eval.get_num_local_dofs();
 
        // gather our local convection dofs
        local_conv_dofs.format();
        gather_conv(local_conv_dofs, dof_mapping);
 
        // compute mesh width if necessary
        if(need_mean_v || need_local_h || need_streamdiff)
        {
          // reset local h and delta
          local_h = local_delta = DataType(0);
 
          // evaluate trafo and space at barycenter
          trafo_eval(trafo_data, barycenter);
          space_eval(space_data, trafo_data);
 
          // compute velocity at barycenter
          mean_v.format();
          for(int i(0); i < num_local_dofs; ++i)
            mean_v.axpy(space_data.phi[i].value, local_conv_dofs[i]);
 
          // compute norm of mean velocity
          const DataType local_norm_v = mean_v.norm_euclid();
 
          // do we have a non-zero velocity?
          if(local_norm_v > tol_eps)
          {
            // compute local mesh width w.r.t. mean velocity
            local_h = trafo_eval.width_directed(mean_v) * local_norm_v;
 
            // compute local Re_T
            const DataType local_re = (local_norm_v * local_h) / this->sd_nu;
 
            // compute local delta
            local_delta = this->sd_delta * (local_h / this->sd_v_norm) * (DataType(2)*local_re) / (DataType(1) + local_re);
          }
        }
      }
 
      void prepare_point(int cubature_point)
      {
        // compute trafo data
        trafo_eval(trafo_data, cubature_rule.get_point(cubature_point));
 
        // compute basis function data
        space_eval(space_data, trafo_data);
 
        // evaluate convection function and its gradient (if required)
        if(need_loc_v || need_conv || need_streamdiff)
        {
          loc_v.format();
          for(int i(0); i < num_local_dofs; ++i)
          {
            // update velocity value
            loc_v.axpy(space_data.phi[i].value, local_conv_dofs[i]);
          }
        }
        if(need_loc_grad_v || need_conv_frechet)
        {
          loc_grad_v.format();
          for(int i(0); i < num_local_dofs; ++i)
          {
            // update velocity gradient
            loc_grad_v.add_outer_product(local_conv_dofs[i], space_data.phi[i].grad);
          }
        }
 
        // evaluate streamline diffusion coefficients
        if(need_streamdiff)
        {
          for(int i(0); i < num_local_dofs; ++i)
          {
            streamdiff_coeffs[i] = Tiny::dot(loc_v, space_data.phi[i].grad);
          }
        }
      }
 
      void finish()
      {
        // finish evaluators
        space_eval.finish();
        trafo_eval.finish();
 
        // finish dof mapping
        dof_mapping.finish();
      }
 
      void combine()
      {
        // nothing to do here
      }
    }; // class BurgersAssemblyTaskBase<...>
 
    template<typename Job_, typename DataType_, int block_size_>
    class BurgersBlockedAssemblyTaskBase :
      public BurgersAssemblyTaskBase<Job_, DataType_>
    {
    public:
      typedef BurgersAssemblyTaskBase<Job_, DataType_> BaseClass;
 
      typedef DataType_ DataType;
 
      using BaseClass::max_local_dofs;
      using BaseClass::num_local_dofs;
 
    protected:
      typedef Tiny::Matrix<DataType_, block_size_, block_size_> MatrixValue;
      typedef Tiny::Matrix<MatrixValue, max_local_dofs, max_local_dofs> LocalMatrixType;
      LocalMatrixType local_matrix;
 
    public:
      explicit BurgersBlockedAssemblyTaskBase(const Job_& job_) :
        BaseClass(job_),
        local_matrix()
      {
      }
 
      void assemble_burgers_point(const DataType weight)
      {
        // assemble diffusion matrix?
        if(this->need_diff)
        {
          // assemble gradient-tensor diffusion
 
          // test function loop
          for(int i(0); i < this->num_local_dofs; ++i)
          {
            // trial function loop
            for(int j(0); j < this->num_local_dofs; ++j)
            {
              // compute scalar value
              const DataType value = this->nu * weight * Tiny::dot(this->space_data.phi[i].grad, this->space_data.phi[j].grad);
 
              // update local matrix
              local_matrix[i][j].add_scalar_main_diag(value);
            }
          }
 
          // assemble deformation tensor?
          if(this->deformation)
          {
            // test function loop
            for(int i(0); i < this->num_local_dofs; ++i)
            {
              // trial function loop
              for(int j(0); j < this->num_local_dofs; ++j)
              {
                // add outer product of grad(phi) and grad(psi)
                local_matrix[i][j].add_outer_product(this->space_data.phi[j].grad, this->space_data.phi[i].grad, this->nu * weight);
              }
            }
          }
        }
 
        // assemble convection?
        if(this->need_conv)
        {
          // test function loop
          for(int i(0); i < this->num_local_dofs; ++i)
          {
            // trial function loop
            for(int j(0); j < this->num_local_dofs; ++j)
            {
              // compute scalar value
              const DataType value = this->beta * weight * this->space_data.phi[i].value * Tiny::dot(this->loc_v, this->space_data.phi[j].grad);
 
              // update local matrix
              local_matrix[i][j].add_scalar_main_diag(value);
            }
          }
        }
 
        // assemble convection Frechet?
        if(this->need_conv_frechet)
        {
          // test function loop
          for(int i(0); i < this->num_local_dofs; ++i)
          {
            // trial function loop
            for(int j(0); j < this->num_local_dofs; ++j)
            {
              // compute scalar value
              const DataType value = this->frechet_beta * weight * this->space_data.phi[i].value * this->space_data.phi[j].value;
 
              // update local matrix
              local_matrix[i][j].axpy(value, this->loc_grad_v);
            }
          }
        }
 
        // assemble reaction?
        if(this->need_reac)
        {
          // test function loop
          for(int i(0); i < this->num_local_dofs; ++i)
          {
            // trial function loop
            for(int j(0); j < this->num_local_dofs; ++j)
            {
              // compute scalar value
              const DataType value = this->theta * weight *  this->space_data.phi[i].value * this->space_data.phi[j].value;
 
              // update local matrix
              local_matrix[i][j].add_scalar_main_diag(value);
            }
          }
        }
      }
 
      void assemble_streamline_diffusion(const DataType weight)
      {
        // assemble streamline diffusion?
        if(this->need_streamdiff && (this->local_delta > this->tol_eps))
        {
          // test function loop
          for(int i(0); i < this->num_local_dofs; ++i)
          {
            // trial function loop
            for(int j(0); j < this->num_local_dofs; ++j)
            {
              // compute scalar value
              const DataType value = this->local_delta * weight * this->streamdiff_coeffs[i] * this->streamdiff_coeffs[j];
 
              // update local matrix
              local_matrix[i][j].add_scalar_main_diag(value);
            }
          }
        }
      }
 
      void assemble()
      {
        local_matrix.format();
 
        // loop over all quadrature points and integrate
        for(int point(0); point < this->cubature_rule.get_num_points(); ++point)
        {
          // prepare trafo and space for cubature point
          this->prepare_point(point);
 
          // precompute cubature weight
          const DataType weight = this->trafo_data.jac_det * this->cubature_rule.get_weight(point);
 
          // assemble the Burgers operator
          this->assemble_burgers_point(weight);
 
          // assemble the streamline diffusion
          this->assemble_streamline_diffusion(weight);
        } // continue with next cubature point
      }
    }; // class BurgersBlockedAssemblyTaskBase<...>
 
    template<typename Matrix_, typename Space_, typename ConvVector_ = typename Matrix_::VectorTypeR>
    class BurgersBlockedMatrixAssemblyJob :
      public BurgersAssemblyJobBase<typename Matrix_::DataType, Space_, ConvVector_>
    {
    public:
      typedef Matrix_ MatrixType;
      typedef typename MatrixType::DataType DataType;
 
      typedef BurgersAssemblyJobBase<DataType, Space_, ConvVector_> BaseClass;
 
      // no nonsense, please
      static_assert(Matrix_::BlockHeight == Matrix_::BlockWidth, "only square matrix blocks are supported here");
 
      static constexpr int block_size = Matrix_::BlockHeight;
 
      MatrixType& matrix;
 
    public:
      explicit BurgersBlockedMatrixAssemblyJob(Matrix_& matrix_, const ConvVector_& conv_vector,
        const Space_& space_, const String& cubature_name) :
        BaseClass(conv_vector, space_, cubature_name),
        matrix(matrix_)
      {
      }
 
    public:
      class Task :
        public BurgersBlockedAssemblyTaskBase<BurgersBlockedMatrixAssemblyJob, DataType, block_size>
      {
      public:
        typedef BurgersBlockedAssemblyTaskBase<BurgersBlockedMatrixAssemblyJob, DataType, block_size> BaseClass;
 
      protected:
        typename MatrixType::ScatterAxpy scatter_matrix;
 
      public:
        explicit Task(const BurgersBlockedMatrixAssemblyJob& job_) :
          BaseClass(job_),
          scatter_matrix(job_.matrix)
        {
        }
 
        // prepare, assemble and finish are already implemented in the base classes
 
        void scatter()
        {
          this->scatter_matrix(this->local_matrix, this->dof_mapping, this->dof_mapping, this->alpha);
        }
      }; // class BurgersBlockedMatrixAssemblyJob::Task
    }; // class BurgersBlockedMatrixAssemblyJob
 
    template<typename Vector_, typename Space_, typename ConvVector_ = Vector_>
    class BurgersBlockedDiagonalAssemblyJob :
      public BurgersAssemblyJobBase<typename Vector_::DataType, Space_, ConvVector_>
    {
    public:
      typedef Vector_ VectorType;
      typedef typename VectorType::DataType DataType;
 
      typedef BurgersAssemblyJobBase<DataType, Space_, ConvVector_> BaseClass;
 
      static constexpr int block_size = Vector_::BlockSize;
 
      VectorType& vector;
 
    public:
      explicit BurgersBlockedDiagonalAssemblyJob(Vector_& vector_, const ConvVector_& conv_vector,
        const Space_& space_, const String& cubature_name) :
        BaseClass(conv_vector, space_, cubature_name),
        vector(vector_)
      {
      }
 
    public:
      class Task :
        public BurgersBlockedAssemblyTaskBase<BurgersBlockedDiagonalAssemblyJob, DataType, block_size>
      {
      public:
        typedef BurgersBlockedAssemblyTaskBase<BurgersBlockedDiagonalAssemblyJob, DataType, block_size> BaseClass;
 
      protected:
        Vector_& vector;
 
      public:
        explicit Task(const BurgersBlockedDiagonalAssemblyJob& job_) :
          BaseClass(job_),
          vector(job_.vector)
        {
        }
 
        // prepare, assemble and finish are already implemented in the base classes
 
        void scatter()
        {
          auto* val_vec = vector.elements();
          for(int loc_i = 0; loc_i < this->dof_mapping.get_num_local_dofs(); ++loc_i)
          {
            const Index cur_i = this->dof_mapping.get_index(loc_i);
            for(int k = 0; k < block_size; ++k)
              val_vec[cur_i][k] += this->local_matrix[loc_i][loc_i][k][k] * this->alpha;
          }
        }
      }; // class BurgersBlockedMatrixAssemblyJob::Task
    }; // class BurgersBlockedMatrixAssemblyJob
 
    template<typename Vector_, typename Space_, typename ConvVector_ = Vector_>
    class BurgersBlockedVectorAssemblyJob :
      public BurgersAssemblyJobBase<typename Vector_::DataType, Space_, ConvVector_>
    {
    public:
      typedef Vector_ VectorType;
      typedef typename VectorType::DataType DataType;
 
      typedef BurgersAssemblyJobBase<DataType, Space_, ConvVector_> BaseClass;
 
      static constexpr int block_size = Vector_::BlockSize;
 
      VectorType& vector_rhs;
 
      const VectorType& vector_sol;
 
    public:
      explicit BurgersBlockedVectorAssemblyJob(Vector_& rhs_vector_, const Vector_& sol_vector_,
        const ConvVector_& conv_vector, const Space_& space_, const String& cubature_name) :
        BaseClass(conv_vector, space_, cubature_name),
        vector_rhs(rhs_vector_),
        vector_sol(sol_vector_)
      {
        XASSERTM(&rhs_vector_ != &sol_vector_, "rhs and solution vectors must not be the same object");
        // compare void addresses to avoid warnings in case the classes are different
        XASSERTM((void*)&rhs_vector_ != (void*)&conv_vector, "rhs and convection vectors must not be the same object");
      }
 
    public:
      class Task :
        public BurgersBlockedAssemblyTaskBase<BurgersBlockedVectorAssemblyJob, DataType, block_size>
      {
      public:
        typedef BurgersBlockedAssemblyTaskBase<BurgersBlockedVectorAssemblyJob, DataType, block_size> BaseClass;
 
        using BaseClass::max_local_dofs;
 
      protected:
        typename VectorType::GatherAxpy gather_vec_sol;
        typename VectorType::ScatterAxpy scatter_vec_rhs;
        const bool sol_equals_conv;
 
        Tiny::Vector<Tiny::Vector<DataType, block_size>, max_local_dofs> local_sol_dofs;
 
        Tiny::Vector<Tiny::Vector<DataType, block_size>, max_local_dofs> local_vector;
 
      public:
        explicit Task(const BurgersBlockedVectorAssemblyJob& job_) :
          BaseClass(job_),
          gather_vec_sol(job_.vector_sol),
          scatter_vec_rhs(job_.vector_rhs),
          sol_equals_conv((void*)&job_.vector_sol == (void*)&job_.convection_vector),
          local_sol_dofs(),
          local_vector()
        {
        }
 
        void prepare(const Index cell)
        {
          BaseClass::prepare(cell);
 
          // gather local solution dofs if required
          if(!sol_equals_conv)
          {
            this->local_sol_dofs.format();
            this->gather_vec_sol(this->local_sol_dofs, this->dof_mapping);
          }
        }
 
        void compute_local_vector()
        {
          local_vector.format();
 
          // compute local vector from local matrix and sol vector dofs
          if(!sol_equals_conv)
          {
            for(int i(0); i < this->num_local_dofs; ++i)
              for(int j(0); j < this->num_local_dofs; ++j)
                this->local_vector[i].add_mat_vec_mult(this->local_matrix(i,j), this->local_sol_dofs[j]);
          }
          else
          {
            for(int i(0); i < this->num_local_dofs; ++i)
              for(int j(0); j < this->num_local_dofs; ++j)
                this->local_vector[i].add_mat_vec_mult(this->local_matrix(i,j), this->local_conv_dofs[j]);
          }
        }
 
        void assemble()
        {
          // assemble the local matrix
          BaseClass::assemble();
 
          // compute the local vector from the matrix
          compute_local_vector();
        }
 
        void scatter()
        {
          this->scatter_vec_rhs(this->local_vector, this->dof_mapping, this->alpha);
        }
      }; // class BurgersBlockedVectorAssemblyJob::Task
    }; // class BurgersBlockedVectorAssemblyJob<...>
 
    template<typename Job_, typename DataType_>
    class BurgersScalarAssemblyTaskBase :
      public BurgersAssemblyTaskBase<Job_, DataType_>
    {
    public:
      typedef BurgersAssemblyTaskBase<Job_, DataType_> BaseClass;
 
      typedef DataType_ DataType;
 
      using BaseClass::max_local_dofs;
      using BaseClass::num_local_dofs;
 
    protected:
      typedef Tiny::Matrix<DataType, max_local_dofs, max_local_dofs> LocalMatrixType;
      LocalMatrixType local_matrix;
 
    public:
      explicit BurgersScalarAssemblyTaskBase(const Job_& job_) :
        BaseClass(job_),
        local_matrix()
      {
      }
 
      void assemble_burgers_point(const DataType weight)
      {
        // assemble diffusion matrix?
        if(this->need_diff)
        {
          // assemble gradient-tensor diffusion
 
          // test function loop
          for(int i(0); i < this->num_local_dofs; ++i)
          {
            // trial function loop
            for(int j(0); j < this->num_local_dofs; ++j)
            {
              local_matrix[i][j] += this->nu * weight
                * Tiny::dot(this->space_data.phi[i].grad, this->space_data.phi[j].grad);
            }
          }
        }
 
        // assemble convection?
        if(this->need_conv)
        {
          // test function loop
          for(int i(0); i < this->num_local_dofs; ++i)
          {
            // trial function loop
            for(int j(0); j < this->num_local_dofs; ++j)
            {
              local_matrix[i][j] += this->beta * weight
                * this->space_data.phi[i].value * Tiny::dot(this->loc_v, this->space_data.phi[j].grad);
            }
          }
        }
 
        // assemble reaction?
        if(this->need_reac)
        {
          // test function loop
          for(int i(0); i < this->num_local_dofs; ++i)
          {
            // trial function loop
            for(int j(0); j < this->num_local_dofs; ++j)
            {
              local_matrix[i][j] += this->theta * weight
                * this->space_data.phi[i].value * this->space_data.phi[j].value;
            }
          }
        }
      }
 
      void assemble_streamline_diffusion(const DataType weight)
      {
        // assemble streamline diffusion?
        if(this->need_streamdiff && (this->local_delta > this->tol_eps))
        {
          // test function loop
          for(int i(0); i < this->num_local_dofs; ++i)
          {
            // trial function loop
            for(int j(0); j < this->num_local_dofs; ++j)
            {
              local_matrix[i][j] += this->local_delta * weight
                * this->streamdiff_coeffs[i] * this->streamdiff_coeffs[j];
            }
          }
        }
      }
 
      void assemble()
      {
        local_matrix.format();
 
        // loop over all quadrature points and integrate
        for(int point(0); point < this->cubature_rule.get_num_points(); ++point)
        {
          // prepare trafo and space for cubature point
          this->prepare_point(point);
 
          // precompute cubature weight
          const DataType weight = this->trafo_data.jac_det * this->cubature_rule.get_weight(point);
 
          // assemble the Burgers operator
          this->assemble_burgers_point(weight);
 
          // assembles the streamline diffusion stabilization
          this->assemble_streamline_diffusion(weight);
        } // continue with next cubature point
      }
    }; // class BurgersScalarAssemblyTaskBase<...>
 
    template<typename Matrix_, typename Space_, typename ConvVector_>
    class BurgersScalarMatrixAssemblyJob :
      public BurgersAssemblyJobBase<typename Matrix_::DataType, Space_, ConvVector_>
    {
    public:
      typedef Matrix_ MatrixType;
      typedef typename MatrixType::DataType DataType;
 
      typedef BurgersAssemblyJobBase<DataType, Space_, ConvVector_> BaseClass;
 
      MatrixType& matrix;
 
    public:
      explicit BurgersScalarMatrixAssemblyJob(Matrix_& matrix_, const ConvVector_& conv_vector,
        const Space_& space_, const String& cubature_name) :
        BaseClass(conv_vector, space_, cubature_name),
        matrix(matrix_)
      {
      }
 
    public:
      class Task :
        public BurgersScalarAssemblyTaskBase<BurgersScalarMatrixAssemblyJob, DataType>
      {
      public:
        typedef BurgersScalarAssemblyTaskBase<BurgersScalarMatrixAssemblyJob, DataType> BaseClass;
 
      protected:
        typename MatrixType::ScatterAxpy scatter_matrix;
 
      public:
        explicit Task(const BurgersScalarMatrixAssemblyJob& job_) :
          BaseClass(job_),
          scatter_matrix(job_.matrix)
        {
        }
 
        // prepare, assemble and finish are already implemented in the base classes
 
        void scatter()
        {
          this->scatter_matrix(this->local_matrix, this->dof_mapping, this->dof_mapping, this->alpha);
        }
      }; // class BurgersScalarMatrixAssemblyJob::Task
    }; // class BurgersScalarMatrixAssemblyJob
 
    template<typename Vector_, typename Space_, typename ConvVector_>
    class BurgersScalarVectorAssemblyJob :
      public BurgersAssemblyJobBase<typename Vector_::DataType, Space_, ConvVector_>
    {
    public:
      typedef Vector_ VectorType;
      typedef typename VectorType::DataType DataType;
 
      typedef BurgersAssemblyJobBase<DataType, Space_, ConvVector_> BaseClass;
 
      VectorType& vector_rhs;
 
      const VectorType& vector_sol;
 
    public:
      explicit BurgersScalarVectorAssemblyJob(Vector_& rhs_vector_, const Vector_& sol_vector_,
        const ConvVector_& conv_vector, const Space_& space_, const String& cubature_name) :
        BaseClass(conv_vector, space_, cubature_name),
        vector_rhs(rhs_vector_),
        vector_sol(sol_vector_)
      {
        XASSERTM(&rhs_vector_ != &sol_vector_, "rhs and solution vectors must not be the same object");
      }
 
    public:
      class Task :
        public BurgersScalarAssemblyTaskBase<BurgersScalarVectorAssemblyJob, DataType>
      {
      public:
        typedef BurgersScalarAssemblyTaskBase<BurgersScalarVectorAssemblyJob, DataType> BaseClass;
 
        using BaseClass::max_local_dofs;
 
      protected:
        typename VectorType::GatherAxpy gather_vec_sol;
        typename VectorType::ScatterAxpy scatter_vec_rhs;
 
        Tiny::Vector<DataType, max_local_dofs> local_sol_dofs;
 
        Tiny::Vector<DataType, max_local_dofs> local_vector;
 
      public:
        explicit Task(const BurgersScalarVectorAssemblyJob& job_) :
          BaseClass(job_),
          gather_vec_sol(job_.vector_sol),
          scatter_vec_rhs(job_.vector_rhs),
          local_sol_dofs(),
          local_vector()
        {
        }
 
        void prepare(const Index cell)
        {
          BaseClass::prepare(cell);
 
          // gather local solution vector dofs
          this->local_sol_dofs.format();
          this->gather_vec_sol(this->local_sol_dofs, this->dof_mapping);
        }
 
        void compute_local_vector()
        {
          // compute local vector from local matrix and sol vector dofs
          this->local_vector.format();
          for(int i(0); i < this->num_local_dofs; ++i)
            for(int j(0); j < this->num_local_dofs; ++j)
              this->local_vector[i] += this->local_matrix(i,j) * this->local_sol_dofs[j];
          //this->local_vector.set_mat_vec_mult(this->local_matrix, this->local_sol_dofs);
        }
 
        void assemble()
        {
          // assemble the local matrix
          BaseClass::assemble();
 
          // compute local vector from the matrix
          compute_local_vector();
        }
 
        void scatter()
        {
          this->scatter_vec_rhs(this->local_vector, this->dof_mapping, this->alpha);
        }
      }; // class BurgersScalarVectorAssemblyJob::Task
    }; // class BurgersScalarVectorAssemblyJob<...>
  } // namespace Assembly
} // namespace FEAT