weakdivergencemodel.hh

#ifndef __DUNE_ACFEM_MODELS_MODULES_WEAKDIVERGENCEMODEL_HH__
#define __DUNE_ACFEM_MODELS_MODULES_WEAKDIVERGENCEMODEL_HH__
 
#include "../operatorparts/modeladapter.hh"
 
namespace Dune {
 
  namespace ACFem {
 
    template<class GridFunction>
    class WeakDivergenceOperatorParts
      : public OperatorPartsExpression<WeakDivergenceOperatorParts<GridFunction> >
    {
      typedef WeakDivergenceOperatorParts       ThisType;
      typedef DefaultOperatorParts<ThisType>    BaseType;
      typedef OperatorPartsInterface<ThisType>  InterfaceType;
      typedef OperatorPartsTraits<ThisType>     TraitsType;
      typedef typename GridFunction::LocalFunctionType LocalFunctionType;
     public:
      typedef GridFunction WeakDivergenceFunctionType;
      typedef typename WeakDivergenceFunctionType::DiscreteFunctionSpaceType::FunctionSpaceType WeakDivergenceFunctionSpaceType;
      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,
        dimWorld = dimDomain,
      };
 
      static_assert(dimRange == 1 && dimDomain == dimWorld,
                    "This is only meant for dimRange == 1 and dimDomain == dimWorld.");
 
      WeakDivergenceOperatorParts(const Fem::Function<WeakDivergenceFunctionSpaceType, WeakDivergenceFunctionType>& function,
                                  const std::string& name = "")
        : function_(function),
          localFunction_(function_()),
          name_(name == "" ? "(div " + function_().name() + ")" : name)
      {}
 
      std::string name() const
      {
        return name_;
      }
 
      // Interface methods that need to be reimplemented
 
      template<class Entity>
      void setEntity(const Entity& entity) const
      {
        localFunction_.init(entity);
      }
 
      template<class Entity, class Point>
      void flux(const Entity& entity,
                const Point &x,
                const RangeType& value,
                const JacobianRangeType& jacobian,
                JacobianRangeType& flux) const
      {
        typename WeakDivergenceFunctionType::RangeType dataValue;
        localFunction_.evaluate(x, dataValue);
 
        flux = 0.;
 
        // value * IdentityMatrix
        for (int i = 0; i < FunctionSpaceType::dimRange; ++i) {
          flux[i] -= dataValue;
        }
      }
 
      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 = 0.;
      }
 
      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
      {
        typename WeakDivergenceFunctionType::JacobianRangeType divergenceData;
        localFunction_.jacobian(x, divergenceData);
 
        result[0] = 0.;
        for (int i = 0; i < dimWorld; ++i) {
          result[0] += divergenceData[i][i];
        }
      }
 
     protected:
      ExpressionStorage<WeakDivergenceFunctionType> function_;
      mutable LocalFunctionType localFunction_;
      const std::string name_;
    };
 
    template<class GridFunction>
    struct OperatorPartsTraits<WeakDivergenceOperatorParts<GridFunction> >
      : public DefaultOperatorPartsTraits<typename GridFunction::FunctionSpaceType::ScalarFunctionSpaceType>
    {
      enum StructureFlags
      {
        isLinear = true, // it is even constant
        isSymmetric = true, // it is even constant
        isSemiDefinite = true // it is even constant
      };
 
      enum ConstituentFlags {
        hasFlux = true, // tested by divergence of test functions
        hasSources = false, // none, after integration by parts
        hasRobinFlux = false,
      };
    };
 
 
 
    template<class GridFunction>
    static inline
    WeakDivergenceOperatorParts<GridFunction>
    weakDivergenceOperatorParts(const Fem::Function<typename GridFunction::FunctionSpaceType, GridFunction>& f,
                                const std::string& name = "")
    {
      typedef WeakDivergenceOperatorParts<GridFunction> OperatorPartsType;
      return OperatorPartsType(f, name);
    }
 
 
    template<class GridFunction>
    static inline
    OperatorPartsAdapterModel<WeakDivergenceOperatorParts<GridFunction>,
                              typename GridFunction::GridPartType>
    weakDivergenceModel(const Fem::Function<typename GridFunction::FunctionSpaceType, GridFunction>& f,
                        const std::string& name = "")
    {
      typedef WeakDivergenceOperatorParts<GridFunction> OperatorPartsType;
      typedef typename GridFunction::GridPartType GridPartType;
      return OperatorPartsAdapterModel<OperatorPartsType, GridPartType>(OperatorPartsType(f, name));
    }
 
 
 
 
  } // namespace ACFem
 
}  //Namespace Dune
 
 
#endif // __DUNE_ACFEM_MODELS_MODULES_WEAKDIVERGENCEMODEL_HH__