laplacianmodel.hh

#ifndef __DUNE_ACFEM_MODELS_MODULES_LAPLACIANMODEL_HH__
#define __DUNE_ACFEM_MODELS_MODULES_LAPLACIANMODEL_HH__
 
#include "../operatorparts/modeladapter.hh"
 
namespace Dune {
 
  namespace ACFem {
 
    template<class FunctionSpace>
    class LaplacianOperatorParts
      : public OperatorPartsExpression<LaplacianOperatorParts<FunctionSpace> >
    {
      typedef 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
      LaplacianOperatorParts(const std::string& name = "")
        : name_(name == "" ? "Laplace" : name)
      {}
 
      std::string name() const
      {
        return name_;
      }
 
      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
      {
        flux = jacobian;
      }
 
      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
      {
        // This is then the trace of the Hessian
        for (int i = 0; i < dimRange; ++i) {
          result[i] = 0;
          for (int j = 0; j < dimDomain; ++j) {
            result[i] -= hessian[i][j][j];
          }
        }
      }
 
     protected:
      const std::string name_;
    };
 
    template<class FunctionSpace>
    struct OperatorPartsTraits<LaplacianOperatorParts<FunctionSpace> >
      : public DefaultOperatorPartsTraits<FunctionSpace>
    {
      enum StructureFlags
      {
        isLinear = true,
        isSymmetric = true,
        isSemiDefinite = true
      };
 
      enum ConstituentFlags {
        hasFlux = true,
        hasSources = false,
        hasRobinFlux = false,
      };
    };
 
 
    template<class Object>
    static inline
    LaplacianOperatorParts<typename Object::FunctionSpaceType>
    laplacianOperatorParts(const Object& object, const std::string& name = "")
    {
      typedef LaplacianOperatorParts<typename Object::FunctionSpaceType> OperatorPartsType;
      return OperatorPartsType(name);
    }
 
    template<class Object>
    static inline
    OperatorPartsAdapterModel<LaplacianOperatorParts<typename Object::FunctionSpaceType>,
                              typename Object::GridPartType>
    laplacianModel(const Object& object, const std::string& name = "")
    {
      typedef LaplacianOperatorParts<typename Object::FunctionSpaceType> OperatorPartsType;
      typedef typename Object::GridPartType GridPartType;
      return OperatorPartsAdapterModel<OperatorPartsType, GridPartType>(OperatorPartsType(name));
    }
 
 
 
 
  } // namespace ACFem
 
}  //Namespace Dune
 
 
#endif // __DUNE_ACFEM_MODELS_MODULES_LAPLACIANMODEL_HH__