burgers_velo_material_assembly_job.hpp

 
// FEAT3: Finite Element Analysis Toolbox, Version 3
// Copyright (C) 2010 - 2023 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
#ifndef KERNEL_ASSEMBLY_BURGERS_VELO_MATERIAL_ASSEMBLY_JOB_HPP
#define KERNEL_ASSEMBLY_BURGERS_VELO_MATERIAL_ASSEMBLY_JOB_HPP 1
 
#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 ViscFunc_, typename ViscDerFunc_, typename ConvVector_>
    class BurgersVeloMaterialAssemblyJobBase
    {
    public:
      typedef DataType_ DataType;
 
      typedef Space_ SpaceType;
 
      typedef ConvVector_ ConvVectorType;
 
      typedef ViscFunc_ ViscosityFunctionType;
 
      typedef ViscDerFunc_ ViscosityDerivFunctionType;
 
      ViscosityFunctionType visco_func;
 
      ViscosityDerivFunctionType visco_der_func;
 
      bool deformation;
 
      DataType_ nu;
 
      DataType_ theta;
 
      DataType_ beta;
 
      DataType_ frechet_beta;
 
      DataType aT;
 
      DataType_ sd_delta;
 
      DataType_ sd_nu;
 
      DataType_ sd_v_norm;
 
      DataType_ reg_eps;
 
      const ConvVector_& convection_vector;
 
      const Space_& space;
 
      Cubature::DynamicFactory cubature_factory;
 
      DataType frechet_material;
 
      DataType alpha;
 
    public:
      explicit BurgersVeloMaterialAssemblyJobBase(const ConvVector_& conv_vector, const Space_& space_, const String& cubature,
                                              ViscosityFunctionType _visc_fun, ViscosityDerivFunctionType _visc_der_fun) :
        visco_func(_visc_fun),
        visco_der_func(_visc_der_fun),
        deformation(false),
        nu(DataType_(0)),
        theta(DataType_(0)),
        beta(DataType_(0)),
        frechet_beta(DataType_(0)),
        aT(DataType(1)),
        sd_delta(DataType_(0)),
        sd_nu(DataType_(0)),
        sd_v_norm(DataType_(0)),
        reg_eps(DataType_(1E-100)),
        convection_vector(conv_vector),
        space(space_),
        cubature_factory(cubature),
        frechet_material(DataType(0)),
        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 BurgersVeloMaterialAssemblyJobBase<...>
 
    template<typename Job_, typename DataType_, typename Function_>
    class BurgersVeloMaterialAssemblyTaskBase
    {
    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;
 
      typedef Function_ FunctionType;
 
      const Job_& job;
 
      const FunctionType visco_func;
 
      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 material_frechet;
 
      const DataType aT;
 
      const bool need_diff, need_conv, need_conv_frechet, need_reac, need_streamdiff, need_material_frechet;
 
      DataType_ gamma_dot, reg_eps, nu_loc;
 
      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::Matrix<DataType, conv_dim, conv_dim> strain_rate_tensor_2;
 
      Tiny::Vector<DataType, max_local_dofs> streamdiff_coeffs;
 
      Tiny::Vector<DataType, shape_dim> barycenter;
 
      DataType local_h;
 
      DataType local_delta;
 
      DataType alpha = DataType(1);
 
    public:
      explicit BurgersVeloMaterialAssemblyTaskBase(const Job_& job_) :
        job(job_),
        visco_func(job_.visco_func),
        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),
        material_frechet(job_.frechet_material),
        aT(job_.aT),
        need_diff(true),
        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_material_frechet(material_frechet > DataType(0)),
        gamma_dot(),
        reg_eps(job_.reg_eps),
        nu_loc(nu),
        need_loc_v(need_conv),
        need_loc_grad_v(true),
        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(),
        strain_rate_tensor_2(),
        streamdiff_coeffs(),
        barycenter(),
        local_h(DataType(0)),
        local_delta(DataType(0)),
        alpha(job.alpha)
      {
        // 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);
          }
        }
 
        //default to 1, only assemble if carreau is needed...
        strain_rate_tensor_2.set_transpose(loc_grad_v);
        strain_rate_tensor_2.axpy(DataType(1.0), loc_grad_v);
        gamma_dot = Math::sqrt(DataType(0.5))*strain_rate_tensor_2.norm_frobenius();
 
        //default to nu
        nu_loc = visco_func(gamma_dot, aT);
      }
 
      void finish()
      {
        // finish evaluators
        space_eval.finish();
        trafo_eval.finish();
 
        // finish dof mapping
        dof_mapping.finish();
      }
 
      void combine()
      {
        // nothing to do here
      }
    }; // class BurgersVeloMaterialAssemblyTaskBase<...>
 
    template<typename Job_, typename DataType_, int block_size_, typename ViscFun_, typename ViscDerFun_>
    class BurgersVeloMaterialBlockedAssemblyTaskBase :
      public BurgersVeloMaterialAssemblyTaskBase<Job_, DataType_, ViscFun_>
    {
    public:
      typedef BurgersVeloMaterialAssemblyTaskBase<Job_, DataType_, ViscFun_> BaseClass;
 
      typedef DataType_ DataType;
      typedef ViscDerFun_ ViscosityDerivFunctionType;
 
      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;
      ViscosityDerivFunctionType visco_der_func;
 
    public:
      explicit BurgersVeloMaterialBlockedAssemblyTaskBase(const Job_& job_) :
        BaseClass(job_),
        local_matrix(),
        visco_der_func(job_.visco_der_func)
      {
      }
 
      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_loc * 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_loc * weight);
              }
            }
            // frechet carreau term "only" makes sense in defo formulation
            if(this->need_material_frechet)
            {
              const DataType fac =  this->material_frechet * visco_der_func(this->gamma_dot, this->aT)/(this->gamma_dot + this->reg_eps);
              //std::cout << fac ;
              // test function loop
              for(int i(0); i < this->num_local_dofs; ++i)
              {
                Tiny::Vector<DataType, BaseClass::ConvVectorType::BlockSize> du_grad_phi;
                du_grad_phi.set_mat_vec_mult(this->strain_rate_tensor_2, this->space_data.phi[i].grad);
 
                // trial function loop
                for(int j(0); j < this->num_local_dofs; ++j)
                {
                  Tiny::Vector<DataType, BaseClass::ConvVectorType::BlockSize> du_grad_psi;
                  du_grad_psi.set_mat_vec_mult(this->strain_rate_tensor_2, this->space_data.phi[j].grad);
                  // add outer product of grad(phi) and grad(psi)
                  local_matrix[i][j].add_outer_product(du_grad_phi, du_grad_psi, fac * 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 BurgersVeloMaterialBlockedAssemblyTaskBase<...>
 
    template<typename Matrix_, typename Space_, typename ViscoFunc_, typename ViscoDerFunc_, typename ConvVector_ = typename Matrix_::VectorTypeR>
    class BurgersVeloMaterialBlockedMatrixAssemblyJob :
      public BurgersVeloMaterialAssemblyJobBase<typename Matrix_::DataType, Space_, ViscoFunc_, ViscoDerFunc_, ConvVector_>
    {
    public:
      typedef Matrix_ MatrixType;
      typedef typename MatrixType::DataType DataType;
 
      typedef BurgersVeloMaterialAssemblyJobBase<DataType, Space_, ViscoFunc_, ViscoDerFunc_, 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 BurgersVeloMaterialBlockedMatrixAssemblyJob(Matrix_& matrix_, const ConvVector_& conv_vector,
        const Space_& space_, const String& cubature_name, ViscoFunc_ visc_fun, ViscoDerFunc_ visc_der_func) :
        BaseClass(conv_vector, space_, cubature_name, visc_fun, visc_der_func),
        matrix(matrix_)
      {
      }
 
    public:
      class Task :
        public BurgersVeloMaterialBlockedAssemblyTaskBase<BurgersVeloMaterialBlockedMatrixAssemblyJob, DataType, block_size, ViscoFunc_, ViscoDerFunc_>
      {
      public:
        typedef BurgersVeloMaterialBlockedAssemblyTaskBase<BurgersVeloMaterialBlockedMatrixAssemblyJob, DataType, block_size, ViscoFunc_, ViscoDerFunc_> BaseClass;
 
      protected:
        typename MatrixType::ScatterAxpy scatter_matrix;
 
      public:
        explicit Task(const BurgersVeloMaterialBlockedMatrixAssemblyJob& 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 BurgersVeloMaterialBlockedMatrixAssemblyJob::Task
    }; // class BurgersVeloMaterialBlockedMatrixAssemblyJob
 
    template<typename DiagonalVector_, typename Space_, typename ViscoFunc_, typename ViscoDerFunc_, typename ConvVector_ = DiagonalVector_>
    class BurgersVeloMaterialBlockedDiagonalAssemblyJob :
      public BurgersVeloMaterialAssemblyJobBase<typename DiagonalVector_::DataType, Space_, ViscoFunc_, ViscoDerFunc_, ConvVector_>
    {
    public:
      typedef DiagonalVector_ VectorType;
      typedef typename VectorType::DataType DataType;
 
      typedef BurgersVeloMaterialAssemblyJobBase<DataType, Space_, ViscoFunc_, ViscoDerFunc_, ConvVector_> BaseClass;
 
      static constexpr int block_size = VectorType::BlockSize;
 
      VectorType& vector;
 
    public:
      explicit BurgersVeloMaterialBlockedDiagonalAssemblyJob(DiagonalVector_& vector_, const ConvVector_& conv_vector,
        const Space_& space_, const String& cubature_name, ViscoFunc_ visc_fun, ViscoDerFunc_ visc_der_func) :
        BaseClass(conv_vector, space_, cubature_name, visc_fun, visc_der_func),
        vector(vector_)
      {
      }
 
    public:
      class Task :
        public BurgersVeloMaterialBlockedAssemblyTaskBase<BurgersVeloMaterialBlockedDiagonalAssemblyJob, DataType, block_size, ViscoFunc_, ViscoDerFunc_>
      {
      public:
        typedef BurgersVeloMaterialBlockedAssemblyTaskBase<BurgersVeloMaterialBlockedDiagonalAssemblyJob, DataType, block_size, ViscoFunc_, ViscoDerFunc_> BaseClass;
 
      protected:
        DiagonalVector_& vector;
 
      public:
        explicit Task(const BurgersVeloMaterialBlockedDiagonalAssemblyJob& 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 BurgersVeloMaterialBlockedMatrixAssemblyJob::Task
    }; // class BurgersVeloMaterialBlockedMatrixAssemblyJob
 
    template<typename Vector_, typename Space_, typename ViscFunc_, typename ViscDerFunc_, typename ConvVector_ = Vector_>
    class BurgersVeloMaterialBlockedVectorAssemblyJob :
      public BurgersVeloMaterialAssemblyJobBase<typename Vector_::DataType, Space_, ViscFunc_, ViscDerFunc_, ConvVector_>
    {
    public:
      typedef Vector_ VectorType;
      typedef typename VectorType::DataType DataType;
 
      typedef BurgersVeloMaterialAssemblyJobBase<DataType, Space_, ViscFunc_, ViscDerFunc_, ConvVector_> BaseClass;
 
      static constexpr int block_size = Vector_::BlockSize;
 
      VectorType& vector_rhs;
 
      const VectorType& vector_sol;
 
    public:
      explicit BurgersVeloMaterialBlockedVectorAssemblyJob(Vector_& rhs_vector_, const Vector_& sol_vector_,
        const ConvVector_& conv_vector, const Space_& space_, const String& cubature_name,
        ViscFunc_ visc_fun, ViscDerFunc_ visc_der_func) :
        BaseClass(conv_vector, space_, cubature_name, visc_fun, visc_der_func),
        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 BurgersVeloMaterialBlockedAssemblyTaskBase<BurgersVeloMaterialBlockedVectorAssemblyJob, DataType, block_size, ViscFunc_, ViscDerFunc_>
      {
      public:
        typedef BurgersVeloMaterialBlockedAssemblyTaskBase<BurgersVeloMaterialBlockedVectorAssemblyJob, DataType, block_size, ViscFunc_, ViscDerFunc_> 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 BurgersVeloMaterialBlockedVectorAssemblyJob& 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 BurgersVeloMaterialBlockedVectorAssemblyJob::Task
    }; // class BurgersVeloMaterialBlockedVectorAssemblyJob<...>
 
    template<typename Job_, typename DataType_, typename ViscFunc_, typename ViscDerFunc_>
    class BurgersVeloMaterialScalarAssemblyTaskBase :
      public BurgersVeloMaterialAssemblyTaskBase<Job_, DataType_, ViscFunc_>
    {
    public:
      typedef BurgersVeloMaterialAssemblyTaskBase<Job_, DataType_, ViscFunc_> 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;
      ViscDerFunc_ visco_der_fun;
 
    public:
      explicit BurgersVeloMaterialScalarAssemblyTaskBase(const Job_& job_) :
        BaseClass(job_),
        local_matrix(),
        visco_der_fun(job_.visc_der_fun)
      {
      }
 
      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_loc * weight
                * Tiny::dot(this->space_data.phi[i].grad, this->space_data.phi[j].grad);
            }
          }
 
          // \TODO: additional frechet term?
        }
 
        // 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 BurgersVeloMaterialScalarAssemblyTaskBase<...>
 
    template<typename Matrix_, typename Space_, typename ViscFunc_, typename ViscDerFunc_, typename ConvVector_>
    class BurgersVeloMaterialScalarMatrixAssemblyJob :
      public BurgersVeloMaterialAssemblyJobBase<typename Matrix_::DataType, Space_, ViscFunc_, ViscDerFunc_, ConvVector_>
    {
    public:
      typedef Matrix_ MatrixType;
      typedef typename MatrixType::DataType DataType;
 
      typedef BurgersVeloMaterialAssemblyJobBase<DataType, Space_, ViscFunc_, ViscDerFunc_, ConvVector_> BaseClass;
 
      MatrixType& matrix;
 
    public:
      explicit BurgersVeloMaterialScalarMatrixAssemblyJob(Matrix_& matrix_, const ConvVector_& conv_vector,
        const Space_& space_, const String& cubature_name,
        ViscFunc_ visc_fun, ViscDerFunc_ visc_der_fun) :
        BaseClass(conv_vector, space_, cubature_name, visc_fun, visc_der_fun),
        matrix(matrix_)
      {
      }
 
    public:
      class Task :
        public BurgersVeloMaterialScalarAssemblyTaskBase<BurgersVeloMaterialScalarMatrixAssemblyJob, DataType, ViscFunc_, ViscDerFunc_>
      {
      public:
        typedef BurgersVeloMaterialScalarAssemblyTaskBase<BurgersVeloMaterialScalarMatrixAssemblyJob, DataType, ViscFunc_, ViscDerFunc_> BaseClass;
 
      protected:
        typename MatrixType::ScatterAxpy scatter_matrix;
 
      public:
        explicit Task(const BurgersVeloMaterialScalarMatrixAssemblyJob& 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 BurgersVeloMaterialScalarMatrixAssemblyJob::Task
    }; // class BurgersVeloMaterialScalarMatrixAssemblyJob
 
    template<typename Vector_, typename Space_, typename ViscoFunc_, typename ViscoDerFunc_, typename ConvVector_>
    class BurgersVeloMaterialScalarVectorAssemblyJob :
      public BurgersVeloMaterialAssemblyJobBase<typename Vector_::DataType, Space_, ViscoFunc_, ViscoDerFunc_, ConvVector_>
    {
    public:
      typedef Vector_ VectorType;
      typedef typename VectorType::DataType DataType;
 
      typedef BurgersVeloMaterialAssemblyJobBase<DataType, Space_, ViscoFunc_, ViscoDerFunc_, ConvVector_> BaseClass;
 
      VectorType& vector_rhs;
 
      const VectorType& vector_sol;
 
    public:
      explicit BurgersVeloMaterialScalarVectorAssemblyJob(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 BurgersVeloMaterialScalarAssemblyTaskBase<BurgersVeloMaterialScalarVectorAssemblyJob, DataType, ViscoFunc_, ViscoDerFunc_>
      {
      public:
        typedef BurgersVeloMaterialScalarAssemblyTaskBase<BurgersVeloMaterialScalarVectorAssemblyJob, DataType, ViscoFunc_, ViscoDerFunc_> 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 BurgersVeloMaterialScalarVectorAssemblyJob& 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 BurgersVeloMaterialScalarVectorAssemblyJob::Task
    }; // class BurgersVeloMaterialScalarVectorAssemblyJob<...>
  } // namespace Assembly
} // namespace FEAT
 
#endif // KERNEL_ASSEMBLY_BURGERS_VELO_MATERIAL_ASSEMBLY_JOB_HPP