deformationtensormodel.hh

#ifndef __DUNE_ACFEM_MODELS_MODULES_DEFORMATIONTENSORMODEL_HH__
#define __DUNE_ACFEM_MODELS_MODULES_DEFORMATIONTENSORMODEL_HH__
 
#include "../../expressions/expressionoperations.hh"
#include "../../expressions/terminal.hh"
 
#include "../modelbase.hh"
#include "../expressiontraits.hh"
 
namespace Dune {
 
  namespace ACFem::PDEModel {
 
    template<class FunctionSpace>
    class DeformationTensorModel
      : public ModelBase<FunctionSpace>
      , public Expressions::SelfExpression<DeformationTensorModel<FunctionSpace> >
      , public MPL::UniqueTags<TypedValueExpression, SemiPositiveExpression, SymmetricModel>
    {
      using ThisType = DeformationTensorModel;
      using BaseType = ModelBase<FunctionSpace>;
     public:
      using typename BaseType::RangeType;
      using typename BaseType::JacobianRangeType;
      using typename BaseType::HessianRangeType;
 
      using BaseType::dimDomain;
      using BaseType::dimRange;
      enum {
        dimWorld = dimDomain
      };
 
      static_assert(dimDomain == dimRange && dimDomain == dimWorld,
                    "This is meant for dimensionworld vector-fields, "
                    "like fluids, deformations etc.");
 
      // Interface methods that need to be reimplemented
 
      DeformationTensorModel(const std::string& name = "")
        : name_(name == "" ? "(D_iU_j+ D_jU_i)" : name)
      {}
 
      std::string name() const
      {
        return name_;
      }
 
      JacobianRangeType linearizedFlux(const JacobianRangeType& jacobian) const
      {
        JacobianRangeType flux;
        // Why is there no transposeMatrix stuff, or an adaper??? anyhow ...
        for (int i = 0; i < dimWorld; ++i) {
          flux[i][i] = 2.*jacobian[i][i];
          for (int j = i+1; j < dimWorld; ++j) {
            flux[i][j] = flux[j][i] = jacobian[i][j] + jacobian[j][i];
          }
        }
        return flux;
      }
 
      template<class Point>
      RangeType fluxDivergence(const HessianRangeType& hessian) const
      {
        RangeType result;
        // just do it in a straight-forward fashion ...
        for (int alpha = 0; alpha < dimWorld; ++alpha) {
          result[alpha] = 0;
          for (int i = 0; i < dimWorld; ++i) {
            // For divergence free vector fields the second term is zero.
            result[alpha] -= hessian[alpha][i][i] + hessian[i][alpha][i];
          }
        }
        return result;
      }
 
     protected:
      std::string name_;
    };
 
 
    template<class Object>
    inline auto deformationTensorModel(const Object& object, const std::string& name = "")
    {
      typedef DeformationTensorModel<typename Object::FunctionSpaceType> ModelType;
      return expressionClosure(ModelType(name));
    }
 
 
 
 
  } // namespace ACFem
 
  namespace ACFem {
 
    using PDEModel::deformationTensorModel;
 
  }
 
}  //Namespace Dune
 
 
#endif // __DUNE_ACFEM_MODELS_MODULES_DEFORMATIONTENSORMODEL_HH__