plaplacianmodel.hh

#ifndef __DUNE_ACFEM_MODELS_MODULES_PLAPLACIANMODEL_HH__
#define __DUNE_ACFEM_MODELS_MODULES_PLAPLACIANMODEL_HH__
 
#include "../operatorparts/modeladapter.hh"
 
namespace Dune {
 
  namespace ACFem {
 
    template<class FunctionSpace>
    class P_LaplacianOperatorParts
      : public OperatorPartsExpression<P_LaplacianOperatorParts<FunctionSpace> >
    {
      typedef P_LaplacianOperatorParts              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;
 
      enum {
        dimDomain = FunctionSpaceType::dimDomain,
        dimRange = FunctionSpaceType::dimRange
      };
 
      typedef typename FunctionSpaceType::DomainFieldType DomainFieldType;
      typedef typename FunctionSpaceType::RangeFieldType RangeFieldType;
 
      // Interface methods that need to be reimplemented
 
      P_LaplacianOperatorParts(const RangeFieldType& exponent, const std::string& name = "")
        : exponent_(exponent - 2.),
          name_(name == "" ? "(|grad(U)|^" + std::to_string(exponent_) + " grad(U))" : name)
      {}
 
      std::string name() const
      {
        return name_;
      }
 
      template<class Entity, class Point>
      void flux(const Entity& entity,
                const Point &x,
                const RangeType& value,
                const JacobianRangeType& jacobian,
                JacobianRangeType& flux) const
      {
        RangeFieldType factor = std::pow(jacobian.frobenius_norm2(), 0.5*exponent_);
 
        flux = jacobian;
        flux *= factor;
      }
 
      template<class Entity, class Point>
      void linearizedFlux(const RangeType& uBar,
                          const JacobianRangeType& DuBar,
                          const Entity& entity,
                          const Point &x,
                          const RangeType& value,
                          const JacobianRangeType& jacobian,
                          JacobianRangeType& flux) const
      {
        RangeFieldType norm = DuBar.frobenius_norm2();
 
        flux = jacobian;
        flux *= std::pow(norm, 0.5 * exponent_);
 
        JacobianRangeType tmp(DuBar);
 
        RangeFieldType factor = exponent_*std::pow(norm, 0.5*exponent_ - 1.);
        RangeFieldType frobscp = 0;
        for (int i = 0; i< dimRange; ++i) {
          frobscp += (DuBar[i] * jacobian[i]);
        }
 
        tmp *= factor * frobscp;
        flux += tmp;
      }
 
      template<class Entity, class Point>
      void fluxDivergence(const Entity& entity,
                          const Point &x,
                          const RangeType& value,
                          const JacobianRangeType& jacobian,
                          const HessianRangeType& hessian,
                          RangeType& result) const
      {
        RangeFieldType norm    = jacobian.frobenius_norm2();
        RangeFieldType factor1 = std::pow(norm, 0.5*exponent_ - 1.);
        RangeFieldType 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_;
 
        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]);
        }
      }
 
     protected:
      const RangeFieldType exponent_;
      const std::string name_;
    };
 
    template<class FunctionSpace>
    struct OperatorPartsTraits<P_LaplacianOperatorParts<FunctionSpace> >
      : public DefaultOperatorPartsTraits<FunctionSpace>
    {
      enum StructureFlags
      {
        isLinear = false,
        isSymmetric = true,
        isSemiDefinite = true
      };
 
      enum ConstituentFlags {
        hasFlux = true,
        hasSources = false,
        hasRobinFlux = false,
      };
    };
 
 
 
    template<class Object>
    static inline
    P_LaplacianOperatorParts<typename Object::FunctionSpaceType>
    p_LaplacianOperatorParts(const typename Object::FunctionSpaceType::RangeFieldType& p,
                             const Object& object,
                             const std::string& name = "")
    {
      typedef P_LaplacianOperatorParts<typename Object::FunctionSpaceType> OperatorPartsType;
      return OperatorPartsType(p, name);
    }
 
    template<class Object>
    static inline
    OperatorPartsAdapterModel<P_LaplacianOperatorParts<typename Object::FunctionSpaceType>,
                              typename Object::GridPartType>
    p_LaplacianModel(const typename Object::FunctionSpaceType::RangeFieldType& p,
                     const Object& object,
                     const std::string& name = "")
    {
      typedef P_LaplacianOperatorParts<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_PLAPLACIANMODEL_HH__