plaplacianmodel.hh

#ifndef __DUNE_ACFEM_MODELS_MODULES_PLAPLACIANMODEL_HH__
#define __DUNE_ACFEM_MODELS_MODULES_PLAPLACIANMODEL_HH__
 
#include "../../common/literals.hh"
#include "../../expressions/terminal.hh"
#include "../../tensors/tensor.hh"
 
#include "../modelbase.hh"
#include "../expressiontraits.hh"
 
namespace Dune {
 
  namespace ACFem::PDEModel {
 
    using namespace Literals;
 
    template<class FunctionSpace, class PField>
    class P_LaplacianModel
      : public ModelBase<FunctionSpace>
      , public Expressions::SelfExpression<P_LaplacianModel<FunctionSpace, PField> >
      , public MPL::UniqueTags<ConditionalType<ExpressionTraits<PField>::isVolatile, VolatileExpression, void>,
                               ConditionalType<IsConstantExprArg<PField>::value, ConstantExpression, void>,
                               ConditionalType<IsTypedValue<PField>::value, TypedValueExpression, void>,
                               SemiPositiveExpression,
                               SymmetricModel>
    {
      using ThisType = P_LaplacianModel;
      using BaseType = ModelBase<FunctionSpace>;
     public:
      using typename BaseType::DomainType;
      using typename BaseType::RangeFieldType;
      using typename BaseType::RangeType;
      using typename BaseType::JacobianRangeType;
      using typename BaseType::HessianRangeType;
      using PFieldType = PField;
      using ExponentType = std::decay_t<decltype(std::declval<PFieldType>() - 2_f)>;
      using BaseType::dimDomain;
      using BaseType::dimRange;
 
      // Interface methods that need to be reimplemented
 
      P_LaplacianModel(PFieldType&& exponent, const std::string& name = "")
        : exponent_(exponent - 2_f),
          name_(name == "" ? "(|grad(U)|^" + toString(exponent_) + " grad(U))" : name)
      {}
 
      std::string name() const
      {
        return name_;
      }
 
      auto flux(const JacobianRangeType& jacobian) const
      {
        const auto factor = std::pow(jacobian.frobenius_norm2(), 0.5*exponent_);
 
        auto flux = jacobian;
        return flux *= factor;
      }
 
      auto linearizedFlux(const JacobianRangeType& DuBar,
                          const JacobianRangeType& jacobian) const
      {
        const auto norm = DuBar.frobenius_norm2();
 
        auto flux = jacobian;
        flux *= std::pow(norm, 0.5 * exponent_);
 
        auto tmp = DuBar;
 
        auto factor = exponent_*std::pow(norm, 0.5*(exponent_ - 1.));
        auto frobscp = contractInner<2>(DuBar, jacobian);
 
        tmp *= factor * frobscp;
        return flux += tmp;
      }
 
      auto fluxDivergence(const JacobianRangeType& jacobian,
                          const HessianRangeType& hessian) const
      {
        const auto norm    = jacobian.frobenius_norm2();
        const auto factor1 = std::pow(norm, 0.5*exponent_ - 1.);
        const auto factor2 = factor1 * norm;
 
        JacobianRangeType jac_hess = 0;
 
        for (int alpha = 0; alpha < dimRange; ++alpha) {
          DomainType tmp;
          hessian[alpha].mv(jacobian[alpha], tmp);
          jac_hess[alpha] += tmp;
        }
        jac_hess *= factor1 * exponent_;
 
        RangeType result;
        for (int beta = 0; beta < dimRange; ++beta) {
          RangeFieldType laplacian = 0;
          for (int i = 0; i < dimDomain; ++i) {
            laplacian += hessian[beta][i][i];
          }
          result[beta]  = -factor2 * laplacian;  // weighted Laplacian
          result[beta] -= (jac_hess[beta] * jacobian[beta]);
        }
        return result;
      }
 
      const ExponentType& exponent() const { return exponent_; }
 
      constexpr auto p() const { return exponent_ + 2_f; }
 
     protected:
      ExponentType exponent_;
      std::string name_;
    };
 
 
 
    template<class Object, class PField>
    constexpr auto p_LaplacianModel(PField&& p, const Object& object, const std::string& name = "")
    {
      using ModelType = P_LaplacianModel<typename Object::FunctionSpaceType, PField>;
      return expressionClosure(ModelType(std::forward<PField>(p), name));
    }
 
 
 
  } // namespace ACFem::PDEModel
 
  namespace ACFem {
 
    using PDEModel::p_LaplacianModel;
 
  }
 
}  //Namespace Dune
 
 
#endif // __DUNE_ACFEM_MODELS_MODULES_PLAPLACIANMODEL_HH__