burgers_carreau_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_CARREAU_ASSEMBLY_JOB_HPP
#define KERNEL_ASSEMBLY_BURGERS_CARREAU_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 ConvVector_>
    class BurgersCarreauAssemblyJobBase
    {
    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;
 
      bool carreau;
 
      bool frechet_carreau;
 
      DataType_ mu_inf, mu_0, lambda, n, a, reg_eps;
 
    public:
      explicit BurgersCarreauAssemblyJobBase(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),
        carreau(false),
        frechet_carreau(false),
        mu_inf(DataType_(0.)),
        mu_0(DataType_(0.)),
        lambda(DataType_(0.)),
        n(DataType_(0.)),
        a(DataType_(0.)),
        reg_eps(DataType_(1e-100))
      {
      }
 
      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 BurgersCarreauAssemblyJobBase<...>
 
    template<typename Job_, typename DataType_>
    class BurgersCarreauAssemblyTaskBase
    {
    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 bool need_diff, need_conv, need_conv_frechet, need_reac, need_streamdiff, need_carreau, need_carreau_frechet;
 
      DataType_ mu_inf, mu_0, lambda, n, a, reg_eps;
 
      DataType_ gamma_dot, 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;
 
    public:
      explicit BurgersCarreauAssemblyTaskBase(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),
        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_carreau(job_.carreau),
        need_carreau_frechet(job_.frechet_carreau),
        mu_inf(job_.mu_inf),
        mu_0(job_.mu_0),
        lambda(job_.lambda),
        n(job_.n),
        a(job_.a),
        reg_eps(job_.reg_eps),
        gamma_dot(),
        nu_loc(),
        need_loc_v(need_conv),
        need_loc_grad_v(need_conv_frechet || need_carreau),
        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))
      {
        // 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 || need_carreau)
        {
          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...
        gamma_dot = DataType_(1.);
        //default to nu
        nu_loc = nu;
        if(need_carreau)
        {
          strain_rate_tensor_2.set_transpose(loc_grad_v);
          strain_rate_tensor_2.axpy(DataType(1.0), loc_grad_v);
          gamma_dot = Math::sqrt(0.5)*strain_rate_tensor_2.norm_frobenius();
          nu_loc = this->mu_inf + (this->mu_0 - this->mu_inf) * Math::pow((DataType_(1.0)
                                          + Math::pow(this->lambda * this->gamma_dot, this->a)), (this->n-DataType(1.0)) / this->a);
        }
      }
 
      void finish()
      {
        // finish evaluators
        space_eval.finish();
        trafo_eval.finish();
 
        // finish dof mapping
        dof_mapping.finish();
      }
 
      void combine()
      {
        // nothing to do here
      }
    }; // class BurgersCarreauAssemblyTaskBase<...>
 
    template<typename Job_, typename DataType_, int block_size_>
    class BurgersCarreauBlockedAssemblyTaskBase :
      public BurgersCarreauAssemblyTaskBase<Job_, DataType_>
    {
    public:
      typedef BurgersCarreauAssemblyTaskBase<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 BurgersCarreauBlockedAssemblyTaskBase(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_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_carreau_frechet)
            {
              const DataType_ mu_prime = (this->mu_0 - this->mu_inf) * ((this->n-DataType_(1.0)) / this->a) *
                Math::pow( DataType(1.0) + Math::pow(this->lambda * this->gamma_dot, this->a), ((this->n-DataType_(1.0) - this->a) / this->a)) *
                this->a * this->lambda * Math::pow(this->lambda * this->gamma_dot, this->a-DataType(1.0));
              const DataType fac =  mu_prime / (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, this->conv_dim> 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, this->conv_dim> 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 BurgersCarreauBlockedAssemblyTaskBase<...>
 
    template<typename Matrix_, typename Space_, typename ConvVector_ = typename Matrix_::VectorTypeR>
    class BurgersCarreauBlockedMatrixAssemblyJob :
      public BurgersCarreauAssemblyJobBase<typename Matrix_::DataType, Space_, ConvVector_>
    {
    public:
      typedef Matrix_ MatrixType;
      typedef typename MatrixType::DataType DataType;
 
      typedef BurgersCarreauAssemblyJobBase<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 BurgersCarreauBlockedMatrixAssemblyJob(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 BurgersCarreauBlockedAssemblyTaskBase<BurgersCarreauBlockedMatrixAssemblyJob, DataType, block_size>
      {
      public:
        typedef BurgersCarreauBlockedAssemblyTaskBase<BurgersCarreauBlockedMatrixAssemblyJob, DataType, block_size> BaseClass;
 
      protected:
        typename MatrixType::ScatterAxpy scatter_matrix;
 
      public:
        explicit Task(const BurgersCarreauBlockedMatrixAssemblyJob& 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);
        }
      }; // class BurgersCarreauBlockedMatrixAssemblyJob::Task
    }; // class BurgersCarreauBlockedMatrixAssemblyJob
 
    template<typename Vector_, typename Space_, typename ConvVector_ = Vector_>
    class BurgersCarreauBlockedVectorAssemblyJob :
      public BurgersCarreauAssemblyJobBase<typename Vector_::DataType, Space_, ConvVector_>
    {
    public:
      typedef Vector_ VectorType;
      typedef typename VectorType::DataType DataType;
 
      typedef BurgersCarreauAssemblyJobBase<DataType, Space_, ConvVector_> BaseClass;
 
      static constexpr int block_size = Vector_::BlockSize;
 
      VectorType& vector_rhs;
 
      const VectorType& vector_sol;
 
    public:
      explicit BurgersCarreauBlockedVectorAssemblyJob(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 BurgersCarreauBlockedAssemblyTaskBase<BurgersCarreauBlockedVectorAssemblyJob, DataType, block_size>
      {
      public:
        typedef BurgersCarreauBlockedAssemblyTaskBase<BurgersCarreauBlockedVectorAssemblyJob, 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 BurgersCarreauBlockedVectorAssemblyJob& 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);
        }
      }; // class BurgersCarreauBlockedVectorAssemblyJob::Task
    }; // class BurgersCarreauBlockedVectorAssemblyJob<...>
 
    template<typename Job_, typename DataType_>
    class BurgersCarreauScalarAssemblyTaskBase :
      public BurgersCarreauAssemblyTaskBase<Job_, DataType_>
    {
    public:
      typedef BurgersCarreauAssemblyTaskBase<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 BurgersCarreauScalarAssemblyTaskBase(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_loc * 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 BurgersCarreauScalarAssemblyTaskBase<...>
 
    template<typename Matrix_, typename Space_, typename ConvVector_>
    class BurgersCarreauScalarMatrixAssemblyJob :
      public BurgersCarreauAssemblyJobBase<typename Matrix_::DataType, Space_, ConvVector_>
    {
    public:
      typedef Matrix_ MatrixType;
      typedef typename MatrixType::DataType DataType;
 
      typedef BurgersCarreauAssemblyJobBase<DataType, Space_, ConvVector_> BaseClass;
 
      MatrixType& matrix;
 
    public:
      explicit BurgersCarreauScalarMatrixAssemblyJob(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 BurgersCarreauScalarAssemblyTaskBase<BurgersCarreauScalarMatrixAssemblyJob, DataType>
      {
      public:
        typedef BurgersCarreauScalarAssemblyTaskBase<BurgersCarreauScalarMatrixAssemblyJob, DataType> BaseClass;
 
      protected:
        typename MatrixType::ScatterAxpy scatter_matrix;
 
      public:
        explicit Task(const BurgersCarreauScalarMatrixAssemblyJob& 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);
        }
      }; // class BurgersCarreauScalarMatrixAssemblyJob::Task
    }; // class BurgersCarreauScalarMatrixAssemblyJob
 
    template<typename Vector_, typename Space_, typename ConvVector_>
    class BurgersCarreauScalarVectorAssemblyJob :
      public BurgersCarreauAssemblyJobBase<typename Vector_::DataType, Space_, ConvVector_>
    {
    public:
      typedef Vector_ VectorType;
      typedef typename VectorType::DataType DataType;
 
      typedef BurgersCarreauAssemblyJobBase<DataType, Space_, ConvVector_> BaseClass;
 
      VectorType& vector_rhs;
 
      const VectorType& vector_sol;
 
    public:
      explicit BurgersCarreauScalarVectorAssemblyJob(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 BurgersCarreauScalarAssemblyTaskBase<BurgersCarreauScalarVectorAssemblyJob, DataType>
      {
      public:
        typedef BurgersCarreauScalarAssemblyTaskBase<BurgersCarreauScalarVectorAssemblyJob, 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 BurgersCarreauScalarVectorAssemblyJob& 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);
        }
      }; // class BurgersCarreauScalarVectorAssemblyJob::Task
    }; // class BurgersCarreauScalarVectorAssemblyJob<...>
  } // namespace Assembly
} // namespace FEAT
 
#endif // KERNEL_ASSEMBLY_BURGERS_CARREAU_ASSEMBLY_JOB_HPP