common_operators.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
 
// includes, FEAT
#include <kernel/assembly/bilinear_operator.hpp>
#include <iostream>
 
namespace FEAT
{
  namespace Assembly
  {
    namespace Common
    {
      class LaplaceOperator :
        public BilinearOperator
      {
      public:
        static constexpr TrafoTags trafo_config = TrafoTags::none;
        static constexpr SpaceTags test_config = SpaceTags::grad;
        static constexpr SpaceTags trial_config = SpaceTags::grad;
 
        template<typename AsmTraits_>
        class Evaluator :
          public BilinearOperator::Evaluator<AsmTraits_>
        {
        public:
          typedef typename AsmTraits_::DataType DataType;
          typedef typename AsmTraits_::TrafoData TrafoData;
          typedef typename AsmTraits_::TestBasisData TestBasisData;
          typedef typename AsmTraits_::TrialBasisData TrialBasisData;
          typedef DataType ValueType;
 
        public:
          explicit Evaluator(const LaplaceOperator& DOXY(operat))
          {
          }
 
          ValueType eval(const TrialBasisData& phi, const TestBasisData& psi)
          {
            return dot(phi.grad, psi.grad);
          }
        }; // class LaplaceOperator::Evaluator<...>
      }; // class LaplaceOperator
 
      template<int dimension_>
      class LaplaceOperatorBlocked :
        public BilinearOperator
      {
      public:
        static constexpr int BlockHeight = dimension_;
        static constexpr int BlockWidth = dimension_;
 
        static constexpr TrafoTags trafo_config = TrafoTags::none;
        static constexpr SpaceTags test_config = SpaceTags::grad;
        static constexpr SpaceTags trial_config = SpaceTags::grad;
 
        template<typename AsmTraits_>
        class Evaluator :
          public BilinearOperator::Evaluator<AsmTraits_>
        {
        public:
          typedef typename AsmTraits_::DataType DataType;
          typedef typename AsmTraits_::TrafoData TrafoData;
          typedef typename AsmTraits_::TestBasisData TestBasisData;
          typedef typename AsmTraits_::TrialBasisData TrialBasisData;
          typedef Tiny::Matrix<DataType, dimension_, dimension_> ValueType;
 
        public:
          explicit Evaluator(const LaplaceOperatorBlocked& DOXY(operat))
          {
          }
 
          ValueType eval(const TrialBasisData& phi, const TestBasisData& psi)
          {
            ValueType r(DataType(0));
            r.add_scalar_main_diag(Tiny::dot(phi.grad, psi.grad));
            return r;
          }
        }; // class LaplaceOperatorBlocked::Evaluator<...>
      }; // class LaplaceOperatorBlocked
 
      class LaplaceBeltramiOperator :
        public Assembly::BilinearOperator
      {
      public:
        static constexpr TrafoTags trafo_config = TrafoTags::jac_det;
        static constexpr SpaceTags test_config = SpaceTags::ref_grad;
        static constexpr SpaceTags trial_config = SpaceTags::ref_grad;
 
        template<typename AsmTraits_>
        class Evaluator :
          public Assembly::BilinearOperator::Evaluator<AsmTraits_>
        {
        public:
          typedef typename AsmTraits_::DataType DataType;
          typedef typename AsmTraits_::TrafoEvaluator TrafoEvaluator;
          typedef typename AsmTraits_::TrafoData TrafoData;
          typedef typename AsmTraits_::TestBasisData TestBasisData;
          typedef typename AsmTraits_::TrialBasisData TrialBasisData;
          typedef DataType ValueType;
 
        private:
          Tiny::Matrix<DataType, TrafoEvaluator::domain_dim, TrafoEvaluator::domain_dim> gram_mat, gram_inv;
 
        public:
          explicit Evaluator(const LaplaceBeltramiOperator& DOXY(operat))
          {
          }
 
          void set_point(const TrafoData& tau)
          {
            // compute the Jacobi-Gram matrix
            gram_mat.set_gram(tau.jac_mat);
            // and invert it
            gram_inv.set_inverse(gram_mat);
          }
 
          ValueType eval(const TrialBasisData& phi, const TestBasisData& psi)
          {
            return gram_inv.scalar_product(phi.ref_grad, psi.ref_grad);
          }
        }; // class LaplaceBeltramiOperator::Evaluator<...>
      }; // class LaplaceBeltramiOperator
 
      class BiharmonicOperator :
        public BilinearOperator
      {
      public:
        static constexpr TrafoTags trafo_config = TrafoTags::none;
        static constexpr SpaceTags test_config = SpaceTags::hess;
        static constexpr SpaceTags trial_config = SpaceTags::hess;
 
        template<typename AsmTraits_>
        class Evaluator :
          public BilinearOperator::Evaluator<AsmTraits_>
        {
        public:
          typedef typename AsmTraits_::DataType DataType;
          typedef typename AsmTraits_::TrafoData TrafoData;
          typedef typename AsmTraits_::TestBasisData TestBasisData;
          typedef typename AsmTraits_::TrialBasisData TrialBasisData;
          typedef DataType ValueType;
 
        public:
          explicit Evaluator(const BiharmonicOperator& DOXY(operat))
          {
          }
 
          ValueType eval(const TrialBasisData& phi, const TestBasisData& psi)
          {
            return phi.hess.trace() * psi.hess.trace();
          }
        }; // class BiharmonicOperator::Evaluator<...>
      }; // class BiharmonicOperator
 
      class IdentityOperator :
        public BilinearOperator
      {
      public:
        static constexpr TrafoTags trafo_config = TrafoTags::none;
        static constexpr SpaceTags test_config = SpaceTags::value;
        static constexpr SpaceTags trial_config = SpaceTags::value;
 
        template<typename AsmTraits_>
        class Evaluator :
          public BilinearOperator::Evaluator<AsmTraits_>
        {
        public:
          typedef typename AsmTraits_::DataType DataType;
          typedef typename AsmTraits_::TrafoData TrafoData;
          typedef typename AsmTraits_::TestBasisData TestBasisData;
          typedef typename AsmTraits_::TrialBasisData TrialBasisData;
          typedef DataType ValueType;
 
        public:
          explicit Evaluator(const IdentityOperator& DOXY(operat))
          {
          }
 
          ValueType eval(const TrialBasisData& phi, const TestBasisData& psi)
          {
            return phi.value * psi.value;
          }
        }; // class IdentityOperator::Evaluator<...>
      }; // class IdentityOperator
 
      template<int dimension_>
      class IdentityOperatorBlocked :
        public BilinearOperator
      {
      public:
        static constexpr int BlockHeight = dimension_;
        static constexpr int BlockWidth = dimension_;
 
        static constexpr TrafoTags trafo_config = TrafoTags::none;
        static constexpr SpaceTags test_config = SpaceTags::value;
        static constexpr SpaceTags trial_config = SpaceTags::value;
 
        template<typename AsmTraits_>
        class Evaluator :
          public BilinearOperator::Evaluator<AsmTraits_>
        {
        public:
          typedef typename AsmTraits_::DataType DataType;
          typedef Tiny::Matrix<DataType, dimension_, dimension_> ValueType;
          typedef typename AsmTraits_::TrafoData TrafoData;
          typedef typename AsmTraits_::TestBasisData TestBasisData;
          typedef typename AsmTraits_::TrialBasisData TrialBasisData;
 
        public:
          explicit Evaluator(const IdentityOperatorBlocked& DOXY(operat))
          {
          }
 
          ValueType eval(const TrialBasisData& phi, const TestBasisData& psi)
          {
            ValueType r(DataType(0));
            r.add_scalar_main_diag(phi.value * psi.value);
            return r;
          }
        }; // class IdentityOperatorBlocked::Evaluator<...>
      }; // class IdentityOperatorBlocked
 
      class TrialDerivativeOperator :
        public BilinearOperator
      {
      public:
        int deriv;
 
        static constexpr TrafoTags trafo_config = TrafoTags::none;
        static constexpr SpaceTags test_config = SpaceTags::grad;
        static constexpr SpaceTags trial_config = SpaceTags::value;
 
        template<typename AsmTraits_>
        class Evaluator :
          public BilinearOperator::Evaluator<AsmTraits_>
        {
        public:
          typedef typename AsmTraits_::DataType DataType;
          typedef typename AsmTraits_::TrafoData TrafoData;
          typedef typename AsmTraits_::TestBasisData TestBasisData;
          typedef typename AsmTraits_::TrialBasisData TrialBasisData;
          typedef DataType ValueType;
 
        protected:
          const int deriv;
 
        public:
          explicit Evaluator(const TrialDerivativeOperator& operat) :
            deriv(operat.deriv)
          {
          }
 
          // copy pasted since Doxygen does not like the operator part in
          // \copydoc BilinearOperator::Evaluator::operator()
          ValueType eval(const TrialBasisData& phi, const TestBasisData& psi)
          {
            return phi.value * psi.grad[deriv];
          }
        }; // class TrialDerivativeOperator::Evaluator<...>
 
        /*
         *
         * \param[in] derivative
         * The index of the derivative for this operator:
         *  - 0: X-derivative
         *  - 1: Y-derivative
         *  - 2: Z-derivative
         *  - ...
        */
        explicit TrialDerivativeOperator(int derivative) :
          deriv(derivative)
        {
        }
      }; // class TrialDerivativeOperator
 
      class TestDerivativeOperator :
        public BilinearOperator
      {
      public:
        int deriv;
 
        static constexpr TrafoTags trafo_config = TrafoTags::none;
        static constexpr SpaceTags test_config = SpaceTags::grad;
        static constexpr SpaceTags trial_config = SpaceTags::value;
 
        template<typename AsmTraits_>
        class Evaluator :
          public BilinearOperator::Evaluator<AsmTraits_>
        {
        public:
          typedef typename AsmTraits_::DataType DataType;
          typedef typename AsmTraits_::TrafoData TrafoData;
          typedef typename AsmTraits_::TestBasisData TestBasisData;
          typedef typename AsmTraits_::TrialBasisData TrialBasisData;
          typedef DataType ValueType;
 
        protected:
          const int deriv;
 
        public:
          explicit Evaluator(const TestDerivativeOperator& operat) :
            deriv(operat.deriv)
          {
          }
 
          // copy pasted since Doxygen does not like the operator part in
          // \copydoc BilinearOperator::Evaluator::operator()
          ValueType eval(const TrialBasisData& phi, const TestBasisData& psi)
          {
            return phi.value * psi.grad[deriv];
          }
        }; // class TestDerivativeOperator::Evaluator<...>
 
        /*
         *
         * \param[in] derivative
         * The index of the derivative for this operator:
         *  - 0: X-derivative
         *  - 1: Y-derivative
         *  - 2: Z-derivative
         *  - ...
        */
        explicit TestDerivativeOperator(int derivative) :
          deriv(derivative)
        {
        }
      }; // class TestDerivativeOperator
 
      class DivDivOperator :
        public Assembly::BilinearOperator
      {
      public:
        int ir;
        int ic;
 
        static constexpr TrafoTags trafo_config = TrafoTags::none;
        static constexpr SpaceTags test_config = SpaceTags::grad;
        static constexpr SpaceTags trial_config = SpaceTags::grad;
 
        template<typename AsmTraits_>
        class Evaluator :
          public Assembly::BilinearOperator::Evaluator<AsmTraits_>
        {
        public:
          typedef typename AsmTraits_::DataType DataType;
          typedef typename AsmTraits_::TrafoData TrafoData;
          typedef typename AsmTraits_::TestBasisData TestBasisData;
          typedef typename AsmTraits_::TrialBasisData TrialBasisData;
          typedef DataType ValueType;
 
        protected:
          const int ir;
          const int ic;
 
        public:
          explicit Evaluator(const DivDivOperator& operat) :
            ir(operat.ir), ic(operat.ic)
          {
          }
 
          ValueType eval(const TrialBasisData& phi, const TestBasisData& psi)
          {
            return phi.grad[ic] * psi.grad[ir];
          }
        }; // class DivDivOperator::Evaluator<...>
 
        explicit DivDivOperator(int ir_, int ic_) :
          ir(ir_), ic(ic_)
        {
        }
      }; // class DivDivOperator
 
      class DuDvOperator :
        public BilinearOperator
      {
      public:
        int ir;
        int ic;
 
        static constexpr TrafoTags trafo_config = TrafoTags::none;
        static constexpr SpaceTags test_config = SpaceTags::grad;
        static constexpr SpaceTags trial_config = SpaceTags::grad;
 
        template<typename AsmTraits_>
        class Evaluator :
          public BilinearOperator::Evaluator<AsmTraits_>
        {
        public:
          typedef typename AsmTraits_::DataType DataType;
          typedef typename AsmTraits_::TrafoData TrafoData;
          typedef typename AsmTraits_::TestBasisData TestBasisData;
          typedef typename AsmTraits_::TrialBasisData TrialBasisData;
          typedef DataType ValueType;
 
        protected:
          const int ir;
          const int ic;
 
        public:
          explicit Evaluator(const DuDvOperator& operat) :
            ir(operat.ir), ic(operat.ic)
          {
          }
 
          ValueType eval(const TrialBasisData& phi, const TestBasisData& psi)
          {
            return ((ir==ic) ? dot(phi.grad,psi.grad) : DataType(0)) + phi.grad[ir] * psi.grad[ic];
          }
        }; // class DuDvOperator::Evaluator<...>
 
        explicit DuDvOperator(int ir_, int ic_) :
          ir(ir_), ic(ic_)
        {
        }
      }; // class DuDvOperator
 
      template<int dimension_>
      class DuDvOperatorBlocked :
        public BilinearOperator
      {
      public:
        static constexpr int BlockHeight = dimension_;
        static constexpr int BlockWidth = dimension_;
 
        static constexpr TrafoTags trafo_config = TrafoTags::none;
        static constexpr SpaceTags test_config = SpaceTags::grad;
        static constexpr SpaceTags trial_config = SpaceTags::grad;
 
        template<typename AsmTraits_>
        class Evaluator :
          public BilinearOperator::Evaluator<AsmTraits_>
        {
        public:
          typedef typename AsmTraits_::DataType DataType;
          typedef Tiny::Matrix<DataType, dimension_, dimension_> ValueType;
          typedef typename AsmTraits_::TrafoData TrafoData;
          typedef typename AsmTraits_::TestBasisData TestBasisData;
          typedef typename AsmTraits_::TrialBasisData TrialBasisData;
 
        public:
          explicit Evaluator(const DuDvOperatorBlocked& DOXY(operat))
          {
          }
 
          ValueType eval(const TrialBasisData& phi, const TestBasisData& psi)
          {
            ValueType r(DataType(0));
            r.set_outer_product(phi.grad, psi.grad);
            r.add_scalar_main_diag(Tiny::dot(phi.grad, psi.grad));
            return r;
          }
 
          /*OperatorValueType operator()(const TrialBasisData& phi, const TestBasisData& psi)
          {
            OperatorValueType r(DataType(0));
            for(int i(0); i < OperatorValueType::m; ++i)
            {
              r(i,i) = dot(phi.grad, psi.grad) + phi.grad[i]*psi.grad[i];
              for(int j(0); j < i; ++j)
              {
                r(i,j) = phi.grad[i] * psi.grad[j];
                r(j,i) = phi.grad[j] * psi.grad[i];
              }
            }
            return r;
          }*/
        }; // class DuDvOperatorBlocked::Evaluator<...>
      }; // class DuDVOperatorBlocked
 
      template<int dimension_>
      class NormalGradientTrialOperatorBlocked :
        public BilinearOperator
      {
      public:
        static constexpr int BlockHeight = dimension_;
        static constexpr int BlockWidth = dimension_;
 
        static constexpr TrafoTags trafo_config = TrafoTags::normal;
        static constexpr SpaceTags test_config = SpaceTags::value;
        static constexpr SpaceTags trial_config = SpaceTags::grad;
 
        template<typename AsmTraits_>
        class Evaluator :
          public BilinearOperator::Evaluator<AsmTraits_>
        {
        public:
          typedef typename AsmTraits_::DataType DataType;
          typedef Tiny::Matrix<DataType, dimension_, dimension_> ValueType;
          typedef typename AsmTraits_::TrafoData TrafoData;
          typedef typename AsmTraits_::TestBasisData TestBasisData;
          typedef typename AsmTraits_::TrialBasisData TrialBasisData;
 
        protected:
          Tiny::Vector<DataType, dimension_> _normal;
 
        public:
          explicit Evaluator(const NormalGradientTrialOperatorBlocked&) : _normal(DataType(0))
          {
          }
 
          void set_point(const TrafoData& tau)
          {
            _normal = tau.normal;
          }
 
          ValueType eval(const TrialBasisData& phi, const TestBasisData& psi)
          {
            ValueType r(DataType(0));
            r.add_scalar_main_diag(Tiny::dot(_normal, phi.grad) * psi.value);
 
            return r;
          }
        }; // class NormalGradientTrialOperatorBlocked::Evaluator<...>
      }; // class NormalGradientTrialOperatorBlocked
 
      template<int dimension_>
      class NormalTransposedGradientTrialOperatorBlocked :
        public BilinearOperator
      {
      public:
        static constexpr int BlockHeight = dimension_;
        static constexpr int BlockWidth = dimension_;
 
        static constexpr TrafoTags trafo_config = TrafoTags::normal;
        static constexpr SpaceTags test_config = SpaceTags::value;
        static constexpr SpaceTags trial_config = SpaceTags::grad;
 
        template<typename AsmTraits_>
        class Evaluator :
          public BilinearOperator::Evaluator<AsmTraits_>
        {
        public:
          typedef typename AsmTraits_::DataType DataType;
          typedef Tiny::Matrix<DataType, dimension_, dimension_> ValueType;
          typedef typename AsmTraits_::TrafoData TrafoData;
          typedef typename AsmTraits_::TestBasisData TestBasisData;
          typedef typename AsmTraits_::TrialBasisData TrialBasisData;
 
        protected:
          Tiny::Vector<DataType, dimension_> _normal;
 
        public:
          explicit Evaluator(const NormalTransposedGradientTrialOperatorBlocked&) : _normal(DataType(0))
          {
          }
 
          void set_point(const TrafoData& tau)
          {
            _normal = tau.normal;
          }
 
          ValueType eval(const TrialBasisData& phi, const TestBasisData& psi)
          {
            ValueType r(DataType(0));
            r.add_outer_product(phi.grad, _normal, psi.value);
            return r;
          }
        }; // class NormalTransposedGradientTrialOperatorBlocked::Evaluator<...>
 
      }; // class NormalTransposedGradientTrialOperatorBlocked
 
      template<int dimension_>
      class GradientTrialOperatorBlocked :
        public BilinearOperator
      {
      public:
        static constexpr int BlockHeight = dimension_;
        static constexpr int BlockWidth = 1;
 
        static constexpr TrafoTags trafo_config = TrafoTags::none;
        static constexpr SpaceTags test_config = SpaceTags::value;
        static constexpr SpaceTags trial_config = SpaceTags::grad;
 
        template<typename AsmTraits_>
        class Evaluator :
          public BilinearOperator::Evaluator<AsmTraits_>
        {
        public:
          typedef typename AsmTraits_::DataType DataType;
          typedef Tiny::Matrix<DataType, dimension_, 1> ValueType;
          typedef typename AsmTraits_::TrafoData TrafoData;
          typedef typename AsmTraits_::TestBasisData TestBasisData;
          typedef typename AsmTraits_::TrialBasisData TrialBasisData;
 
        public:
          explicit Evaluator(const GradientTrialOperatorBlocked&)
          {
          }
 
          ValueType eval(const TrialBasisData& phi, const TestBasisData& psi)
          {
            ValueType r(DataType(0));
            for(int i(0); i < dimension_; ++i)
              r(i,0) = phi.grad[i] * psi.value;
            return r;
          }
        }; // class GradientTrialOperatorBlocked::Evaluator<...>
      }; // class GradientTrialOperatorBlocked
 
      template<int dimension_>
      class GradientTestOperatorBlocked :
        public BilinearOperator
      {
      public:
        static constexpr int BlockHeight = dimension_;
        static constexpr int BlockWidth = 1;
 
        static constexpr TrafoTags trafo_config = TrafoTags::none;
        static constexpr SpaceTags test_config = SpaceTags::grad;
        static constexpr SpaceTags trial_config = SpaceTags::value;
 
        template<typename AsmTraits_>
        class Evaluator :
          public BilinearOperator::Evaluator<AsmTraits_>
        {
        public:
          typedef typename AsmTraits_::DataType DataType;
          typedef Tiny::Matrix<DataType, dimension_, 1> ValueType;
          typedef typename AsmTraits_::TrafoData TrafoData;
          typedef typename AsmTraits_::TestBasisData TestBasisData;
          typedef typename AsmTraits_::TrialBasisData TrialBasisData;
 
        public:
          explicit Evaluator(const GradientTestOperatorBlocked&)
          {
          }
 
          ValueType eval(const TrialBasisData& phi, const TestBasisData& psi)
          {
            ValueType r(DataType(0));
            for(int i(0); i < dimension_; ++i)
              r(i,0) = phi.value * psi.grad[i];
            return r;
          }
        }; // class GradientTestOperatorBlocked::Evaluator<...>
      }; // class GradientTestOperatorBlocked
 
      template<int dim_, int nsc_>
      class StressDivergenceOperator :
        public Assembly::BilinearOperator
      {
      public:
        static constexpr int BlockHeight = dim_;
        static constexpr int BlockWidth = nsc_;
 
        static constexpr TrafoTags trafo_config = TrafoTags::none;
        static constexpr SpaceTags test_config = SpaceTags::value;
        static constexpr SpaceTags trial_config = SpaceTags::grad;
 
        template<typename AsmTraits_>
        class Evaluator :
          public BilinearOperator::Evaluator<AsmTraits_>
        {
        public:
          typedef Tiny::Matrix<typename AsmTraits_::DataType, dim_, nsc_> ValueType;
          typedef typename AsmTraits_::TestBasisData TestBasisData;
          typedef typename AsmTraits_::TrialBasisData TrialBasisData;
 
        public:
          explicit Evaluator(const StressDivergenceOperator& DOXY(operat)) {}
 
          template<typename T_>
          static void eval(Tiny::Matrix<T_, 2, 4, 2, 4>& R, const Tiny::Vector<T_, 2, 2>& der, const T_ u)
          {
            // u_1 = dx sigma_11 + dy sigma_12 = dx sigma_1 + dy sigma_2
            //     = (dx, dy, 0, 0) : (sigma_11, sigma_12, sigma_21, sigma_22)
            R(0,0) = u * der(0);
            R(0,1) = u * der(1);
            R(0,2) = T_(0);
            R(0,3) = T_(0);
 
            // u_2 = dx sigma_21 + dy sigma_22 = dx sigma_3 + dy sigma_2
            //     = (0, 0, dx, dy) : (sigma_11, sigma_12, sigma_21, sigma_22)
            R(1,0) = T_(0);
            R(1,1) = T_(0);
            R(1,2) = u * der(0);
            R(1,3) = u * der(1);
          }
 
          template<typename T_>
          static void eval(Tiny::Matrix<T_, 2, 3, 2, 3>& R, const Tiny::Vector<T_, 2, 2>& der, const T_ u)
          {
            // we have:
            //  [ sigma_11 sigma_12 ]   [ sigma_1 sigma_3 ]
            //  [ sigma_21 sigma_22 ] = [ sigma_3 sigma_2 ]
 
            // u_1 = dx sigma_11 + dy sigma_12 = dx sigma_1 + dy sigma_3
            //     = (dx, 0, dy) : (sigma_11, sigma_22, sigma_12)
            R(0,0) = u * der(0);
            R(0,1) = T_(0);
            R(0,2) = u * der(1);
 
            // u_2 = dx sigma_21 + dy sigma_22 = dx sigma_3 + dy sigma_2
            //     = (0, dy, dx) : (sigma_11, sigma_22, sigma_12)
            R(1,0) = T_(0);
            R(1,1) = u * der(1);
            R(1,2) = u * der(0);
          }
 
          template<typename T_>
          static void eval(Tiny::Matrix<T_, 3, 9, 3, 9>& R, const Tiny::Vector<T_, 3, 3>& der, const T_ u)
          {
            // u_1 = dx sigma_11 + dy sigma_12 + dz sigma_13 = dx sigma_1 + dy sigma_2 + dz sigma_3
            R(0,0) = u * der(0);
            R(0,1) = u * der(1);
            R(0,2) = u * der(2);
            R(0,3) = T_(0);
            R(0,4) = T_(0);
            R(0,5) = T_(0);
            R(0,6) = T_(0);
            R(0,7) = T_(0);
            R(0,8) = T_(0);
 
            // u_2 = dx sigma_21 + dy sigma_22 + dz sigma_23 = dx sigma_4 + dy sigma_5 + dz sigma_6
            R(1,0) = T_(0);
            R(1,1) = T_(0);
            R(1,2) = T_(0);
            R(1,3) = u * der(0);
            R(1,4) = u * der(1);
            R(1,5) = u * der(2);
            R(1,6) = T_(0);
            R(1,7) = T_(0);
            R(1,8) = T_(0);
 
            // u_3 = dx sigma_31 + dy sigma_32 + dz sigma_33 = dx sigma_7 + dy sigma_8 + dz sigma_9
            R(2,0) = T_(0);
            R(2,1) = T_(0);
            R(2,2) = T_(0);
            R(2,3) = T_(0);
            R(2,4) = T_(0);
            R(2,5) = T_(0);
            R(2,6) = u * der(0);
            R(2,7) = u * der(1);
            R(2,8) = u * der(2);
          }
 
          template<typename T_>
          static void eval(Tiny::Matrix<T_, 3, 6, 3, 6>& R, const Tiny::Vector<T_, 3, 3>& der, const T_ u)
          {
            // we have:
            //  [ sigma_11 sigma_12 sigma_13 ]   [ sigma_1 sigma_4 sigma_6 ]
            //  [ sigma_21 sigma_22 sigma_23 ] = [ sigma_4 sigma_2 sigma_5 ]
            //  [ sigma_22 sigma_23 sigma_33 ]   [ sigma_6 sigma_5 sigma_3 ]
 
            // u_1 = dx sigma_11 + dy sigma_12 + dz sigma_13 = dx sigma_1 + dy sigma_4 + dz sigma_6
            R(0,0) = u * der(0);
            R(0,1) = T_(0);
            R(0,2) = T_(0);
            R(0,3) = u * der(1);
            R(0,4) = T_(0);
            R(0,5) = u * der(2);
 
            // u_2 = dx sigma_21 + dy sigma_22 + dz sigma_23 = dx sigma_4 + dy sigma_2 + dz sigma_5
            R(1,0) = T_(0);
            R(1,1) = u * der(1);
            R(1,2) = T_(0);
            R(1,3) = u * der(0);
            R(1,4) = u * der(2);
            R(1,5) = T_(0);
 
            // u_3 = dx sigma_31 + dy sigma_32 + dz sigma_33 = dx sigma_6 + dy sigma_5 + dz sigma_3
            R(2,0) = T_(0);
            R(2,1) = T_(0);
            R(2,2) = u * der(2);
            R(2,3) = T_(0);
            R(2,4) = u * der(1);
            R(2,5) = u * der(0);
          }
 
          ValueType eval(const TrialBasisData& phi, const TestBasisData& psi)
          {
            ValueType R;
            eval(R, phi.grad, psi.value);
            return R;
          }
        }; // class StressDivergenceOperator::Evaluator
      }; // class StressDivergenceOperator
 
      template<int dim_, int nsc_>
      class StrainRateTensorOperator :
        public Assembly::BilinearOperator
      {
      public:
        static constexpr int BlockHeight = nsc_;
        static constexpr int BlockWidth = dim_;
 
        static constexpr TrafoTags trafo_config = TrafoTags::none;
        static constexpr SpaceTags test_config = SpaceTags::value;
        static constexpr SpaceTags trial_config = SpaceTags::grad;
 
        template<typename AsmTraits_>
        class Evaluator :
          public BilinearOperator::Evaluator<AsmTraits_>
        {
        public:
          typedef Tiny::Matrix<typename AsmTraits_::DataType, nsc_, dim_> ValueType;
          typedef typename AsmTraits_::TestBasisData TestBasisData;
          typedef typename AsmTraits_::TrialBasisData TrialBasisData;
 
        public:
          explicit Evaluator(const StrainRateTensorOperator& DOXY(operat)) {}
 
          template<typename T_>
          static void eval(Tiny::Matrix<T_, 4, 2, 4, 2>& K, const Tiny::Vector<T_, 2, 2>& der, const T_ s)
          {
            // sigma_1 [11] = dx u_1 = (dx, 0) : (u_1, u_2)
            K(0,0) = s * der(0);
            K(0,1) = T_(0);
 
            // sigma_2 [12] = 1/2 * (dy u_1 + dx u_2) = 1/2 * (dy, dx) : (u_1, u_2)
            K(1,0) = s * der(1) / T_(2);
            K(1,1) = s * der(0) / T_(2);
 
            // sigma_3 [21] = 1/2 * (dy u_1 + dx u_2) = 1/2 * (dy, dx) : (u_1, u_2)
            K(2,0) = s * der(1) / T_(2);
            K(2,1) = s * der(0) / T_(2);
 
            // sigma_4 [22] = dy u_2 = (0, dy) : (u_1, u_2)
            K(3,0) = T_(0);
            K(3,1) = s * der(1);
          }
 
          template<typename T_>
          static void eval(Tiny::Matrix<T_, 3, 2, 3, 2>& K, const Tiny::Vector<T_, 2, 2>& der, const T_ s)
          {
            // we have:
            //  [ sigma_11 sigma_12 ]   [ sigma_1 sigma_3 ]
            //  [ sigma_21 sigma_22 ] = [ sigma_3 sigma_2 ]
 
            // sigma_1 [11] = dx u_1 = (dx, 0) : (u_1, u_2)
            K(0,0) = s * der(0);
            K(0,1) = T_(0);
 
            // sigma_2 [22] = dy u_2 = (0, dy) : (u_1, u_2)
            K(1,0) = T_(0);
            K(1,1) = s * der(1);
 
            // sigma_3 [12] = 1/2 * (dy u_1 + dx u_2) = 1/2 * (dy, dx) : (u_1, u_2)
            K(2,0) = s * der(1) / T_(2);
            K(2,1) = s * der(0) / T_(2);
          }
 
          template<typename T_>
          static void eval(Tiny::Matrix<T_, 9, 3, 9, 3>& K, const Tiny::Vector<T_, 3, 3>& der, const T_ s)
          {
            // sigma_1 [11] = dx u_1 = (dx, 0, 0) : (u_1, u_2, u_3)
            K(0,0) = s * der(0);
            K(0,1) = T_(0);
            K(0,2) = T_(0);
 
            // sigma_2 [12] = 1/2 * (dy u_1 + dx u_2) = 1/2 * (dy, dx, 0) : (u_1, u_2, u_3)
            K(1,0) = s * der(1) / T_(2);
            K(1,1) = s * der(0) / T_(2);
            K(1,2) = T_(0);
 
            // sigma_3 [13] = 1/2 * (dz u_1 + dx u_3) = 1/2 * (dz, 0, dx) : (u_1, u_2, u_3)
            K(2,0) = s * der(2) / T_(2);
            K(2,1) = T_(0);
            K(2,2) = s * der(0) / T_(2);
 
            // sigma_4 [21] = 1/2 * (dx u_2 + dy u_1) = 1/2 * (dy, dx, 0) : (u_1, u_2, u_3)
            K(3,0) = s * der(1) / T_(2);
            K(3,1) = s * der(0) / T_(2);
            K(3,2) = T_(0);
 
            // sigma_5 [22] = dy u_2 = (0, dy, 0) : (u_1, u_2, u_3)
            K(4,0) = T_(0);
            K(4,1) = s * der(1);
            K(4,2) = T_(0);
 
            // sigma_6 [23] = 1/2 * (dz u_2 + dy u_3) = 1/2 * (0, dz, dy) : (u_1, u_2, u_3)
            K(5,0) = T_(0);
            K(5,1) = s * der(2) / T_(2);
            K(5,2) = s * der(1) / T_(2);
 
            // sigma_7 [31] = 1/2 * (dx u_3 + dz u_1) = 1/2 * (dz, 0, dx) : (u_1, u_2, u_3)
            K(6,0) = s * der(2) / T_(2);
            K(6,1) = T_(0);
            K(6,1) = s * der(0) / T_(2);
 
            // sigma_8 [32] = 1/2 * (dy u_3 + dz u_2) = 1/2 * (0, dz, dy) : (u_1, u_2, u_3)
            K(7,0) = T_(0);
            K(7,1) = s * der(2) / T_(2);
            K(7,2) = s * der(1) / T_(2);
 
            // sigma_9 [33] = dz u3 = (0, 0, dz) : (u_1, u_2, u_3)
            K(8,0) = T_(0);
            K(8,1) = T_(0);
            K(8,2) = s * der(2);
          }
 
          template<typename T_>
          static void eval(Tiny::Matrix<T_, 6, 3, 6, 3>& K, const Tiny::Vector<T_, 3, 3>& der, const T_ s)
          {
            // we have:
            //  [ sigma_11 sigma_12 sigma_13 ]   [ sigma_1 sigma_4 sigma_6 ]
            //  [ sigma_21 sigma_22 sigma_23 ] = [ sigma_4 sigma_2 sigma_5 ]
            //  [ sigma_22 sigma_23 sigma_33 ]   [ sigma_6 sigma_5 sigma_3 ]
 
            // sigma_1 [11] = dx u_1 = (dx, 0, 0) : (u_1, u_2, u_3)
            K(0,0) = s * der(0);
            K(0,1) = T_(0);
            K(0,2) = T_(0);
 
            // sigma_2 [22] = dy u_2 = (0, dy, 0) : (u_1, u_2, u_3)
            K(1,0) = T_(0);
            K(1,1) = s * der(1);
            K(1,2) = T_(0);
 
            // sigma_3 [33] = dz u3 = (0, 0, dz) : (u_1, u_2, u_3)
            K(2,0) = T_(0);
            K(2,1) = T_(0);
            K(2,2) = s * der(2);
 
            // sigma_4 [12] = 1/2 * (dy u_1 + dx u_2) = 1/2 * (dy, dx, 0) : (u_1, u_2, u_3)
            K(3,0) = s * der(1) / T_(2);
            K(3,1) = s * der(0) / T_(2);
            K(3,2) = T_(0);
 
            // sigma_5 [23] = 1/2 * (dz u_2 + dy u_3) = 1/2 * (0, dz, dy) : (u_1, u_2, u_3)
            K(4,0) = T_(0);
            K(4,1) = s * der(2) / T_(2);
            K(4,2) = s * der(1) / T_(2);
 
            // sigma_6 [13] = 1/2 * (dz u_1 + dx u_3) = 1/2 * (dz, 0, dx) : (u_1, u_2, u_3)
            K(5,0) = s * der(2) / T_(2);
            K(5,1) = T_(0);
            K(5,2) = s * der(0) / T_(2);
          }
 
          ValueType eval(const TrialBasisData& phi, const TestBasisData& psi)
          {
            ValueType K;
            eval(K, phi.grad, psi.value);
            return K;
          }
        }; // class StrainRateTensorOperator::Evaluator
      }; // class StrainRateTensorOperator
    } // namespace Common
  } // namespace Assembly
} // namespace FEAT