transportmodel.hh

#ifndef __DUNE_ACFEM_MODELS_MODULES_TRANSPORTMODEL_HH__
#define __DUNE_ACFEM_MODELS_MODULES_TRANSPORTMODEL_HH__
 
#include <dune/fem/function/localfunction/const.hh>
 
#include "../../expressions/terminal.hh"
 
#include "../modelbase.hh"
 
namespace Dune {
 
  namespace ACFem::PDEModel {
 
    template<class FunctionSpace, class GridFunction>
    class TransportModel
      : public ModelBase<FunctionSpace>
      , public Expressions::SelfExpression<TransportModel<FunctionSpace, GridFunction> >
      , public MPL::UniqueTags<ConditionalType<IsConstantExprArg<GridFunction>::value, ConstantExpression, void>,
                               ConditionalType<IsTypedValue<GridFunction>::value, TypedValueExpression, void> >
    {
      static_assert(HasTag<GridFunction, Fem::HasLocalFunction>::value,
                    "GridFunction must provide a local function");
 
      using ThisType = TransportModel;
      using BaseType = ModelBase<FunctionSpace>;
      using LocalFunctionType = Fem::ConstLocalFunction<std::decay_t<GridFunction> >;
     public:
      using GridFunctionType = GridFunction;
      using GridFunctionDecay = std::decay_t<GridFunctionType>;
      using GridPartType = typename GridFunctionDecay::GridPartType;
 
      using typename BaseType::RangeFieldType;
      using typename BaseType::DomainType;
      using typename BaseType::RangeType;
      using typename BaseType::JacobianRangeType;
      using BaseType::dimDomain;
      using BaseType::dimRange;
 
      enum {
        dimWorld = GridPartType::dimensionworld
      };
 
      static_assert((int)GridFunctionDecay::FunctionSpaceType::dimRange == (int)dimWorld
                    &&
                    (int)dimDomain == (int)dimWorld,
                    "This is meant for dimensionworld vector-fields, "
                    "like fluids, deformations etc.");
 
      // Interface methods that need to be reimplemented
 
      template<class FunctionArg, std::enable_if_t<std::is_constructible<LocalFunctionType, FunctionArg>::value, int> = 0>
      TransportModel(FunctionArg&& function, const std::string& name = "")
        : localFunction_(std::forward<FunctionArg>(function))
        , name_(name == "" ? "(-U_iU_jD_iPhi_j+Bndry)" : name)
      {}
 
      std::string name() const
      {
        return name_;
      }
 
      template<class Entity>
      void bind(const Entity& entity)
      {
        localFunction_.bind(entity);
      }
 
      void unbind()
      {
        localFunction_.unbind();
      }
 
      template<class Intersection>
      auto classifyBoundary(const Intersection& intersection)
      {
        // true is correct as higher-level code has to take care of
        // the Dirichlet-(non-Dirichlet) splitting of the boundary.
        return std::make_pair(true, std::bitset<dimRange>());
      }
 
      template<class Quadrature>
      auto linearizedFlux(const QuadraturePoint<Quadrature> &x, const RangeType& value) const
      {
        // We need to provide the tensor-product of value with the
        // value of function_
        const auto velocity = localFunction_.evaluate(x);
 
        JacobianRangeType flux = 0;
        for (int i = 0; i < dimRange; ++i) {
          for (int j = 0; j < dimWorld; ++j) {
            flux[i][j] = - value[i] * velocity[j];
          }
        }
        return flux;
      }
 
      template<class Quadrature>
      auto fluxDivergence(const QuadraturePoint<Quadrature> &x,
                          const RangeType& value,
                          const JacobianRangeType& jacobian) const
      {
        const auto velocity = localFunction_.evaluate(x);
        const auto velocityJacobian = localFunction_.jacobian(x);
 
        RangeType result;
        jacobian.mv(velocity, result);
 
        RangeFieldType velocityDivergence = 0.;
//        for (int i = 0; i < dimDomain; ++i) {
        for (int i = 0; i < dimDomain; ++i) {
          velocityDivergence += velocityJacobian[i][i];
        }
        result.axpy(velocityDivergence, value);
 
        return result;
      }
 
      template<class Quadrature>
      auto linearizedRobinFlux(const QuadraturePoint<Quadrature> &x,
                               const DomainType& unitOuterNormal,
                               const RangeType& value) const
      {
        const auto velocity = localFunction_.evaluate(x);
        return RangeType(value) *= (velocity * unitOuterNormal); // scalar product
      }
 
     protected:
      LocalFunctionType localFunction_;
      std::string name_;
    };
 
 
    template<class Object, class GridFunction>
    constexpr auto
    transportModel(const Object& object,
                   GridFunction&& velocity,
                   const std::string& name = "")
    {
      return expressionClosure(TransportModel<typename Object::FunctionSpaceType, GridFunction>(std::forward<GridFunction>(velocity), name));
    }
 
    template<class Object, class GridFunction, std::enable_if_t<ExpressionTraits<GridFunction>::isZero, int> = 0>
    constexpr auto
    transportModel(const Object& object,
                   GridFunction&& velocity,
                   const std::string& name = "")
    {
      return zeroModel(object);
    }
 
 
 
 
  } // namespace ACFem::PDEModel
 
  namespace ACFem {
 
    using PDEModel::transportModel;
 
  }
 
}  //Namespace Dune
 
 
#endif // __DUNE_ACFEM_MODELS_MODULES_TRANSPORTMODEL_HH__