pmassmodel.hh

#ifndef __DUNE_ACFEM_MODELS_MODULES_PMASSMODEL_HH__
#define __DUNE_ACFEM_MODELS_MODULES_PMASSMODEL_HH__
 
#include "../operatorparts/modeladapter.hh"
 
namespace Dune {
 
  namespace ACFem {
 
    template<class FunctionSpace>
    class P_MassOperatorParts
      : public OperatorPartsExpression<P_MassOperatorParts<FunctionSpace> >
    {
      typedef P_MassOperatorParts                    ThisType;
      typedef DefaultOperatorParts<ThisType>    BaseType;
      typedef OperatorPartsInterface<ThisType>  InterfaceType;
      typedef OperatorPartsTraits<ThisType>     TraitsType;
     public:
      typedef FunctionSpace FunctionSpaceType;
 
      typedef typename FunctionSpaceType::DomainType DomainType;
      typedef typename FunctionSpaceType::RangeType RangeType;
      typedef typename FunctionSpaceType::JacobianRangeType JacobianRangeType;
      typedef typename FunctionSpaceType::HessianRangeType HessianRangeType;
      typedef typename FunctionSpaceType::RangeFieldType RangeFieldType;
 
      enum {
        dimDomain = FunctionSpaceType::dimDomain,
        dimRange = FunctionSpaceType::dimRange
      };
 
      typedef typename FunctionSpaceType::DomainFieldType DomainFieldType;
 
      P_MassOperatorParts(const RangeFieldType& exponent, const std::string& name = "")
        : exponent_(exponent - 2.),
          name_(name == "" ? "(|U|^" + std::to_string(exponent_) + " U)" : name)
      {}
 
      std::string name() const
      {
        return name_;
      }
 
      template<class Entity, class Point>
      void source(const Entity& entity,
                  const Point &x,
                  const RangeType& value,
                  const JacobianRangeType& jacobian,
                  RangeType& result) const
      {
        RangeFieldType factor = std::pow(value.two_norm2(), 0.5*exponent_);
 
        result  = value;
        result *= factor;
      }
 
      template<class Entity, class Point>
      void linearizedSource(const RangeType& uBar,
                            const JacobianRangeType& DuBar,
                            const Entity& entity,
                            const Point &x,
                            const RangeType& value,
                            const JacobianRangeType& jacobian,
                            RangeType& result) const
      {
        RangeFieldType norm = uBar.two_norm2();
 
        result  = value;
 
        if (exponent_ == 0.0) {
          return;
        }
 
        result *= std::pow(norm, 0.5 * exponent_);
 
        RangeType tmp(uBar);
 
        tmp *= exponent_*std::pow(norm, 0.5*exponent_ - 1.) * (uBar * value);
 
        result += tmp;
      }
 
     protected:
      const RangeFieldType exponent_;
      const std::string name_;
    };
 
    template<class FunctionSpace>
    struct OperatorPartsTraits<P_MassOperatorParts<FunctionSpace> >
      : public DefaultOperatorPartsTraits<FunctionSpace>
    {
      enum StructureFlags
      {
        isLinear = false, 
        isSymmetric = true, 
        isSemiDefinite = true 
      };
 
      enum ConstituentFlags {
        hasFlux = false,
        hasSources = true,
        hasRobinFlux = false,
      };
    };
 
 
    template<class Object>
    static inline
    P_MassOperatorParts<typename Object::FunctionSpaceType>
    p_MassOperatorParts(const typename Object::FunctionSpaceType::RangeFieldType& p,
                        const Object& object, const std::string& name = "")
    {
      typedef P_MassOperatorParts<typename Object::FunctionSpaceType> OperatorPartsType;
      return OperatorPartsType(p, name);
    }
 
    template<class Object>
    static inline
    OperatorPartsAdapterModel<P_MassOperatorParts<typename Object::FunctionSpaceType>,
                              typename Object::GridPartType>
    p_MassModel(const typename Object::FunctionSpaceType::RangeFieldType& p,
                const Object& object,
                const std::string& name = "")
    {
      typedef P_MassOperatorParts<typename Object::FunctionSpaceType> OperatorPartsType;
      typedef typename Object::GridPartType GridPartType;
      return OperatorPartsAdapterModel<OperatorPartsType, GridPartType>(OperatorPartsType(p, name));
    }
 
 
 
 
  } // namespace ACFem
 
}  //Namespace Dune
 
 
#endif // __DUNE_ACFEM_MODELS_MODULES_PMASSMODEL_HH__