meancurvaturemodel.hh

#ifndef __DUNE_ACFEM_MODELS_MODULES_MEANCURVATUREMODEL_HH__
#define __DUNE_ACFEM_MODELS_MODULES_MEANCURVATUREMODEL_HH__
 
#include "../operatorparts/modeladapter.hh"
#include "../../expressions/parameterinterface.hh"
 
namespace Dune {
 
  namespace ACFem {
 
    template<class FunctionSpace, class Parameter = TrivialParameter<typename FunctionSpace::RangeType> >
    class MeanCurvatureOperatorParts
      : public OperatorPartsExpression<MeanCurvatureOperatorParts<FunctionSpace, Parameter> >
    {
      typedef MeanCurvatureOperatorParts       ThisType;
      typedef DefaultOperatorParts<ThisType>   BaseType;
      typedef OperatorPartsInterface<ThisType> InterfaceType;
      typedef OperatorPartsTraits<ThisType>    TraitsType;
     public:
      typedef Parameter ParameterType;
 
      typedef typename TraitsType::FunctionSpaceType 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
      MeanCurvatureOperatorParts(const ParameterInterface<ParameterType>& eps = 1.0,
                                 const std::string& name = "")
        : regularization_(eps),
          name_(name == "" ? "MC" : 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 eps = regularization_().value();
        RangeFieldType factor = std::sqrt(eps*eps + jacobian.frobenius_norm2());
 
        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 eps = regularization_().value();
        RangeFieldType factor = std::sqrt(eps*eps + DuBar.frobenius_norm2());
 
        flux  = jacobian;
        flux /= factor;
 
        JacobianRangeType correction;
 
        correction  = DuBar;
        correction *= (DuBar[0] * jacobian[0]) / factor / factor / factor;
 
        flux -= correction;
      }
 
      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 eps = regularization_().value();
        RangeFieldType factor = std::sqrt(eps*eps + jacobian.frobenius_norm2());
 
        // This is then the trace of the Hessian
        result = 0.;
        for (int j = 0; j < dimDomain; ++j) {
          result -= hessian[0][j][j];
        }
 
        // + the correction from the non-linearity
        RangeFieldType correction = 0;
        for (int i = 0; i < dimDomain; ++i) {
          correction += jacobian[0][i] * hessian[0][i][i] * jacobian[0][i];
          for (int j = i+1; j < dimDomain; ++j) {
            correction += 2.0*jacobian[0][i] * hessian[0][i][j] * jacobian[0][j];
          }
        }
        result += correction/factor/factor;
        result /= factor;
      }
 
     protected:
      ExpressionStorage<ParameterType> regularization_;
      const std::string name_;
    };
 
    template<class FunctionSpace, class Parameter>
    struct OperatorPartsTraits<MeanCurvatureOperatorParts<FunctionSpace, Parameter> >
      : public DefaultOperatorPartsTraits<FunctionSpace>
    {
      enum StructureFlags
      {
        isLinear = false,
        isSymmetric = true,
        isSemiDefinite = true
      };
 
      enum ConstituentFlags {
        hasFlux = true,
        hasSources = false,
        hasRobinFlux = false
      };
    };
 
 
    template<class Object, class Parameter = TrivialParameter<typename Object::FunctionSpaceType::RangeType> >
    static inline
    MeanCurvatureOperatorParts<typename Object::FunctionSpaceType, Parameter>
    meanCurvatureOperatorParts(const Parameter& regularization,
                               const Object& object,
                               const std::string& name = "")
    {
      return MeanCurvatureOperatorParts<typename Object::FunctionSpaceType, Parameter>(regularization, name);
    }
 
    template<class Object, class Parameter = TrivialParameter<typename Object::FunctionSpaceType::RangeType> >
    static inline
    OperatorPartsAdapterModel<MeanCurvatureOperatorParts<typename Object::FunctionSpaceType, Parameter>,
                              typename Object::GridPartType>
    meanCurvatureModel(const Parameter& regularization,
                       const Object& object,
                       const std::string& name = "")
    {
      typedef MeanCurvatureOperatorParts<typename Object::FunctionSpaceType, Parameter> OperatorPartsType;
      typedef typename Object::GridPartType GridPartType;
      return OperatorPartsAdapterModel<OperatorPartsType, GridPartType>(OperatorPartsType(regularization, name));
    }
 
 
 
 
 
  } // namespace ACFem
 
}  //Namespace Dune
 
 
#endif // __DUNE_ACFEM_MODELS_MODULES_MEANCURVATUREMODEL_HH__