transportmodel.hh

#ifndef __DUNE_ACFEM_MODELS_MODULES_TRANSPORTMODEL_HH__
#define __DUNE_ACFEM_MODELS_MODULES_TRANSPORTMODEL_HH__
 
#include "../operatorparts/modeladapter.hh"
 
namespace Dune {
 
  namespace ACFem {
 
    template<class FunctionSpace, class Velocity>
    class TransportOperatorParts
      : public OperatorPartsExpression<TransportOperatorParts<FunctionSpace, Velocity> >
    {
      typedef TransportOperatorParts           ThisType;
      typedef DefaultOperatorParts<ThisType>   BaseType;
      typedef OperatorPartsInterface<ThisType> InterfaceType;
      typedef OperatorPartsTraits<ThisType>    TraitsType;
      typedef typename Velocity::LocalFunctionType LocalVelocityType;
     public:
      typedef Velocity VelocityType;
      typedef typename TraitsType::FunctionSpaceType FunctionSpaceType;
      typedef typename VelocityType::GridPartType GridPartType;
 
      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 = GridPartType::dimensionworld
      };
 
      static_assert((int)VelocityType::FunctionSpaceType::dimRange == (int)dimWorld && dimDomain == dimWorld,
                    "This is meant for dimensionworld vector-fields, "
                    "like fluids, deformations etc.");
 
      typedef typename FunctionSpaceType::DomainFieldType DomainFieldType;
      typedef typename FunctionSpaceType::RangeFieldType RangeFieldType;
 
      // Interface methods that need to be reimplemented
 
      TransportOperatorParts(const Fem::Function<typename VelocityType::FunctionSpaceType, VelocityType>& velocity,
                     const std::string& name = "")
        : velocity_(velocity),
          localVelocity_(velocity_()),
          name_(name == "" ? "(-U_iU_jD_iPhi_j+Bndry)" : name)
      {}
 
      std::string name() const
      {
        return name_;
      }
 
      template<class Entity>
      void setEntity(const Entity& entity) const
      {
        localVelocity_.init(entity);
      }
 
      template<class Intersection>
      bool setIntersection(const Intersection& intersection) const
      {
        // true is correct as higher-level code has to take care of
        // the Dirichlet-(non-Dirichlet) splitting of the boundary.
        return true;
      }
 
      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
      {
        assert(&entity == &localVelocity_.entity());
 
        // We need to provide the tensor-product of value with the
        // value of velocity_
        typename VelocityType::RangeType velocity;
        localVelocity_.evaluate(x, velocity);
 
        flux = 0;
        for (int i = 0; i < dimRange; ++i) {
          for (int j = 0; j < dimWorld; ++j) {
            flux[i][j] = - value[i] * velocity[j];
          }
        }
      }
 
      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
      {
        assert(&entity == &localVelocity_.entity());
 
        typename VelocityType::RangeType velocity;
        localVelocity_.evaluate(x, velocity);
        typename VelocityType::JacobianRangeType velocityJacobian;
        localVelocity_.jacobian(x, velocityJacobian);
 
        jacobian.mv(velocity, result);
 
        RangeFieldType velocityDivergence = 0.;
        for (int i = 0; i < dimDomain; ++i) {
          velocityDivergence += velocityJacobian[i][i];
        }
        result.axpy(velocityDivergence, value);
      }
 
      template<class Intersection, class Point>
      void linearizedRobinFlux(const RangeType& uBar,
                               const Intersection& intersection,
                               const Point& x,
                               const DomainType& unitOuterNormal,
                               const RangeType& value,
                               RangeType& result) const
      {
        typename VelocityType::RangeType velocity;
        localVelocity_.evaluate(x, velocity);
 
        result  = value;
        result *= (velocity * unitOuterNormal); // scalar product
      }
 
     protected:
      ExpressionStorage<VelocityType> velocity_;
      mutable LocalVelocityType localVelocity_;
      const std::string name_;
    };
 
    template<class FunctionSpace, class Velocity>
    struct OperatorPartsTraits<TransportOperatorParts<FunctionSpace, Velocity> >
      : public DefaultOperatorPartsTraits<FunctionSpace>
    {
      enum StructureFlags
      {
        isLinear = true,
        isSymmetric = false,
        isSemiDefinite = false, // Mmmh.
      };
 
      enum ConstituentFlags {
        hasFlux = true,
        hasSources = false,
        hasRobinFlux= true,
      };
    };
 
 
    template<class Object, class Velocity>
    static inline
    TransportOperatorParts<typename Object::FunctionSpaceType, Velocity>
    transportOperatorParts(const Object& object,
                           const Fem::Function<typename Velocity::FunctionSpaceType, Velocity>& velocity,
                           const std::string& name = "")
    {
      return TransportOperatorParts<typename Object::FunctionSpaceType, Velocity>(velocity, name);
    }
 
    template<class Object, class Velocity>
    static inline
    OperatorPartsAdapterModel<TransportOperatorParts<typename Object::FunctionSpaceType, Velocity>, typename Velocity::GridPartType>
    transportModel(const Object& object,
                   const Fem::Function<typename Velocity::FunctionSpaceType, Velocity>& velocity,
                   const std::string& name = "")
    {
      typedef TransportOperatorParts<typename Object::FunctionSpaceType, Velocity> OperatorPartsType;
      typedef typename Velocity::GridPartType GridPartType;
      return OperatorPartsAdapterModel<OperatorPartsType, GridPartType>(OperatorPartsType(velocity, name));
    }
 
 
 
 
  } // namespace ACFem
 
}  //Namespace Dune
 
 
#endif // __DUNE_ACFEM_MODELS_MODULES_TRANSPORTMODEL_HH__