gradientmodel.hh

#ifndef __DUNE_ACFEM_MODELS_MODULES_GRADIENTMODEL_HH__
#define __DUNE_ACFEM_MODELS_MODULES_GRADIENTMODEL_HH__
 
#include <dune/fem/gridpart/common/gridpart.hh>
 
#include "../../common/gridfunctionspace.hh"
#include "../../expressions/terminal.hh"
 
#include "../modelbase.hh"
 
namespace Dune {
 
  namespace ACFem::PDEModel {
 
    template<class FunctionSpace>
    class GradientModel
      : public ModelBase<FunctionSpace,
                         typename Fem::ToNewDimRangeFunctionSpace<FunctionSpace, FunctionSpace::dimDomain>::Type>
      , public Expressions::SelfExpression<GradientModel<FunctionSpace> >
      , public TypedValueExpression
    {
      static_assert(FunctionSpace::dimRange == 1, "Domain function space must be scalar");
      using ThisType = GradientModel;
      using BaseType = ModelBase<FunctionSpace,
                                 typename Fem::ToNewDimRangeFunctionSpace<FunctionSpace, FunctionSpace::dimDomain>::Type>;
     public:
      using typename BaseType::DomainType;
      using typename BaseType::DomainRangeType;
      using typename BaseType::DomainJacobianRangeType;
      using typename BaseType::RangeRangeType;
      using typename BaseType::RangeJacobianRangeType;
      using BaseType::rangeDimRange;
 
      GradientModel(const std::string& name)
        : name_(name == "" ? "(grad p)" : name)
      {}
 
      std::string name() const
      {
        return name_;
      }
 
      auto linearizedFlux(const DomainRangeType& value) const
      {
        RangeJacobianRangeType result(0);
        for (int i = 0; i < rangeDimRange; ++i) {
          result[i][i] = -value;
        }
        return result;
      }
 
      auto fluxDivergence(const DomainJacobianRangeType& jacobian) const
      {
        return jacobian[0];
      }
 
      auto linearizedRobinFlux(const DomainType& unitOuterNormal,
                               const DomainRangeType& value) const
      {
        return RangeRangeType(unitOuterNormal) *= value;
      }
 
      template<class Intersection>
      auto classifyBoundary(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 std::make_pair(true, std::bitset<rangeDimRange>());
      }
 
     protected:
      const std::string name_;
    };
 
 
 
    template<class Object,
             std::enable_if_t<Object::FunctionSpaceType::ScalarFunctionSpaceType::dimRange == 1,
                              int> = 0>
    static inline auto
    gradientModel(const Object& object, const std::string& name = "")
    {
      return expressionClosure(GradientModel<typename Object::FunctionSpaceType::ScalarFunctionSpaceType>(name));
    }
 
    template<class GridPartTraits>
    static inline auto
    gradientModel(const Fem::GridPartInterface<GridPartTraits>& gridPart,
                  const std::string& name = "")
    {
      typedef typename Fem::GridPartInterface<GridPartTraits>::GridPartType GridPartType;
      return gradientModel(ScalarGridFunctionSpace<GridPartType>{}, name);
    }
 
 
 
 
  } // namespace ACFem::PDEModel
 
  namespace ACFem {
 
    using PDEModel::gradientModel;
 
  }
 
}  //Namespace Dune
 
 
#endif // __DUNE_ACFEM_MODELS_MODULES_GRADIENTMODEL_HH__