pmassmodel.hh

#ifndef __DUNE_ACFEM_MODELS_MODULES_PMASSMODEL_HH__
#define __DUNE_ACFEM_MODELS_MODULES_PMASSMODEL_HH__
 
#include "../../common/literals.hh"
#include "../../expressions/terminal.hh"
 
#include "../modelbase.hh"
#include "../expressiontraits.hh"
 
namespace Dune {
 
  namespace ACFem::PDEModel {
 
    using namespace Literals;
 
    template<class FunctionSpace, class PField>
    class P_MassModel
      : public ModelBase<FunctionSpace>
      , public Expressions::SelfExpression<P_MassModel<FunctionSpace, PField> >
      , public MPL::UniqueTags<ConditionalType<ExpressionTraits<PField>::isVolatile, VolatileExpression, void>,
                               ConditionalType<IsConstantExprArg<PField>::value, ConstantExpression, void>,
                               ConditionalType<IsTypedValue<PField>::value, TypedValueExpression, void>,
                               PositiveExpression,
                               SymmetricModel>
    {
      using ThisType = P_MassModel;
      using BaseType = ModelBase<FunctionSpace>;
     public:
      using typename BaseType::RangeFieldType;
      using typename BaseType::RangeType;
      using PFieldType = PField;
      using ExponentType = std::decay_t<decltype(std::declval<PFieldType>() - 2_f)>;
 
      template<class P, std::enable_if_t<std::is_constructible<PFieldType, P>::value, int> = 0>
      P_MassModel(P&& exponent, const std::string& name = "")
        : exponent_(PFieldType(exponent) - 2_f),
          name_(name == "" ? "(|U|^" + toString(exponent_) + " U)" : name)
      {}
 
      std::string name() const
      {
        return name_;
      }
 
      RangeType source(const RangeType& value) const
      {
        auto factor = std::pow(value.two_norm2(), 0.5*exponent_);
 
        auto result = value;
        return result *= factor;
      }
 
      RangeType linearizedSource(const RangeType& uBar,
                                 const RangeType& value) const
      {
        auto norm = uBar.two_norm2();
 
        auto result = value;
 
        if (exponent_ == 0.0) {
          return result;
        }
 
        result *= std::pow(norm, 0.5 * exponent_);
 
        auto tmp = uBar;
 
        tmp *= exponent_*std::pow(norm, 0.5*exponent_ - 1.) * (uBar * value);
 
        return result += tmp;
      }
 
      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_MassModel(PField&& p, const Object& object, const std::string& name = "")
    {
      return expressionClosure(P_MassModel<typename Object::FunctionSpaceType, PField>(std::forward<PField>(p), name));
    }
 
 
 
 
  } // namespace ACFem::PDEModel
 
  namespace ACFem {
 
    using PDEModel::p_MassModel;
 
  }
 
}  //Namespace Dune
 
 
#endif // __DUNE_ACFEM_MODELS_MODULES_PMASSMODEL_HH__