modelinterface.hh

#ifndef __DUNE_ACFEM_MODELINTERFACE_HH__
#define __DUNE_ACFEM_MODELINTERFACE_HH__
 
#include <dune/fem/misc/bartonnackmaninterface.hh>
#include <dune/fem/gridpart/common/gridpart.hh>
#include <dune/fem/space/common/discretefunctionspace.hh>
 
#include "../expressions/expressionoperations.hh"
#include "../models/boundaryindicator.hh" // for EmptyBoundaryIndicator
#include "../operators/functionals/zerofunctional.hh"
#include "../operators/functionals/nonzerofunctional.hh"
#include "../functions/boundarysupportedfunction.hh"
#include "../functions/constantfunction.hh" // for ZeroGridFunction
#include "../functions/gridfunctionexpression.hh"
#include "../functions/boundaryfunctionexpression.hh"
#include "operatorparts/zerooperatorparts.hh"
 
namespace Dune {
 
  namespace ACFem {
 
    template<class Model>
    struct ModelTraits;
 
    template<class FunctionSpace, class GridPart>
    struct DefaultModelTraits
    {
      typedef FunctionSpace FunctionSpaceType;
 
      typedef typename FunctionSpaceType::ScalarFunctionSpaceType ScalarFunctionSpaceType;
 
      typedef GridPart GridPartType;
 
      typedef ZeroOperatorParts<FunctionSpace> OperatorPartsType;
 
      typedef EmptyBoundaryIndicatorType DirichletIndicatorType;
 
      typedef
      BoundarySupportedFunction<ZeroGridFunction<FunctionSpaceType, GridPartType>,
                                DirichletIndicatorType>
      DirichletBoundaryFunctionType;
 
      typedef
      BoundarySupportedFunction<FractionGridFunction<ScalarFunctionSpaceType, GridPartType, 1L, 1UL>,
                                DirichletIndicatorType>
      DirichletWeightFunctionType;
 
      typedef ZeroGridFunction<FunctionSpaceType, GridPartType> BulkForcesFunctionType;
 
      typedef EmptyBoundaryIndicatorType NeumannIndicatorType;
 
      typedef
      BoundarySupportedFunction<ZeroGridFunction<FunctionSpaceType, GridPartType>,
                                NeumannIndicatorType>
      NeumannBoundaryFunctionType;
 
      template<class DiscreteFunctionSpace>
      struct ForcesFunctionalTraits
      {
        typedef DiscreteFunctionSpace DiscreteFunctionSpaceType;
        typedef ZeroFunctional<DiscreteFunctionSpaceType> FunctionalType;
      };
 
    };
 
    template<class ModelType>
    class ModelInterface
      : public Fem::BartonNackmanInterface<ModelInterface<ModelType>, ModelType>
    {
      typedef ModelInterface ThisType;
      typedef ModelType ImplementationType;
      typedef ModelTraits<ImplementationType> TraitsType;
      typedef Fem::BartonNackmanInterface<ThisType, ImplementationType> BaseType;
      using BaseType::asImp;
     protected:
      ModelInterface() {}
      ModelInterface(const ModelInterface& other) {}
      ModelInterface& operator=(const ModelInterface& other) {}
     public:
      typedef typename TraitsType::FunctionSpaceType FunctionSpaceType;
      typedef typename TraitsType::GridPartType GridPartType;
      typedef typename GridPartType::template Codim<0>::EntityType EntityType;
      typedef typename GridPartType::IntersectionType IntersectionType;
 
      typedef typename FunctionSpaceType::DomainType DomainType;
      typedef typename FunctionSpaceType::RangeType RangeType;
      typedef typename FunctionSpaceType::JacobianRangeType JacobianRangeType;
      typedef typename FunctionSpaceType::HessianRangeType HessianRangeType;
 
      typedef typename FunctionSpaceType::DomainFieldType DomainFieldType;
      typedef typename FunctionSpaceType::RangeFieldType RangeFieldType;
 
      enum { dimDomain = FunctionSpaceType::dimDomain,
             dimRange = FunctionSpaceType::dimRange };
 
      typedef typename TraitsType::OperatorPartsType OperatorPartsType;
 
      typedef typename TraitsType::BulkForcesFunctionType BulkForcesFunctionType;
      typedef typename TraitsType::NeumannBoundaryFunctionType NeumannBoundaryFunctionType;
      typedef typename TraitsType::NeumannIndicatorType   NeumannIndicatorType;
 
      template<class DiscreteFunctionSpace>
      struct ForcesFunctionalTraits
        : public TraitsType::template ForcesFunctionalTraits<DiscreteFunctionSpace>
      {};
 
      typedef typename TraitsType::DirichletBoundaryFunctionType DirichletBoundaryFunctionType;
      typedef typename TraitsType::DirichletWeightFunctionType DirichletWeightFunctionType;
      typedef typename TraitsType::DirichletIndicatorType DirichletIndicatorType;
 
      std::string name() const
      {
        return asImp().name();
      }
 
      OperatorPartsType operatorParts() const
      {
        return asImp().operatorParts();
      }
 
      BulkForcesFunctionType bulkForcesFunction(const GridPartType& gridPart) const {
        return asImp().bulkForcesFunction(gridPart);
      }
 
      template<class DiscreteFunctionSpace>
      typename TraitsType::template ForcesFunctionalTraits<DiscreteFunctionSpace>::FunctionalType
      forcesFunctional(const DiscreteFunctionSpace& space) const
      {
        return asImp().forcesFunctional(space);
      }
 
      static_assert(std::is_same<typename DirichletBoundaryFunctionType::IndicatorType,
                                 typename DirichletWeightFunctionType::IndicatorType>::value,
                    "DirichletBoundaryFunctionType and DirichletWeightFunctionType must be "
                    "BoundarySupportedFunction's sub-ordinate to the same BoundaryIndicator type.");
 
      DirichletBoundaryFunctionType dirichletBoundaryFunction(const GridPartType& gridPart) const {
        CHECK_INTERFACE_IMPLEMENTATION(asImp().dirichletBoundaryFunction(gridPart));
        return asImp().dirichletBoundaryFunction(gridPart);
      }
 
      DirichletWeightFunctionType dirichletWeightFunction(const GridPartType& gridPart) const {
        CHECK_INTERFACE_IMPLEMENTATION(asImp().dirichletWeightFunction(gridPart));
        return asImp().dirichletWeightFunction(gridPart);
      }
 
      NeumannBoundaryFunctionType neumannBoundaryFunction(const GridPartType& gridPart) const {
        return asImp().neumannBoundaryFunction(gridPart);
      }
 
      DirichletIndicatorType dirichletIndicator() const {
        return asImp().dirichletIndicator();
      }
 
      NeumannIndicatorType neumannIndicator() const {
        return asImp().neumannIndicator();
      }
    };
 
    template<class ModelType>
    class DefaultModel
      : public ModelInterface<ModelType>
    {
      typedef ModelType ImplementationType;
      typedef ModelInterface<ImplementationType> BaseType;
      typedef DefaultModel<ImplementationType> ThisType;
      typedef ModelTraits<ImplementationType> TraitsType;
      typedef DefaultModelTraits<typename TraitsType::FunctionSpaceType,
                                 typename TraitsType::GridPartType>
      DefaultTraitsType;
     protected:
      DefaultModel() {}
     public:
      typedef typename DefaultTraitsType::FunctionSpaceType FunctionSpaceType;
      typedef typename DefaultTraitsType::GridPartType GridPartType;
 
      typedef typename FunctionSpaceType::DomainType DomainType;
      typedef typename FunctionSpaceType::RangeType RangeType;
      typedef typename FunctionSpaceType::JacobianRangeType JacobianRangeType;
      typedef typename FunctionSpaceType::HessianRangeType HessianRangeType;
 
      typedef typename FunctionSpaceType::DomainFieldType DomainFieldType;
      typedef typename FunctionSpaceType::RangeFieldType RangeFieldType;
 
      enum {
        dimDomain = FunctionSpaceType::dimDomain,
        dimRange = FunctionSpaceType::dimRange
      };
 
      typedef typename BaseType::EntityType EntityType;
      typedef typename BaseType::IntersectionType IntersectionType;
 
      typedef typename DefaultTraitsType::BulkForcesFunctionType BulkForcesFunctionType;
      typedef typename DefaultTraitsType::DirichletBoundaryFunctionType DirichletBoundaryFunctionType;
      typedef typename DefaultTraitsType::DirichletWeightFunctionType DirichletWeightFunctionType;
      typedef typename DefaultTraitsType::NeumannBoundaryFunctionType NeumannBoundaryFunctionType;
 
      typedef typename DefaultTraitsType::DirichletIndicatorType DirichletIndicatorType;
      typedef typename DefaultTraitsType::NeumannIndicatorType NeumannIndicatorType;
      typedef typename TraitsType::OperatorPartsType OperatorPartsType;
 
      template<class DiscreteFunctionSpace>
      struct ForcesFunctionalTraits
        : public TraitsType::template ForcesFunctionalTraits<DiscreteFunctionSpace>
      {};
 
      std::string name() const
      {
        return "DefaultModel";
      }
 
      OperatorPartsType operatorParts() const
      {
        return OperatorPartsType();
      }
 
      BulkForcesFunctionType bulkForcesFunction(const GridPartType& gridPart) const
      {
        return BulkForcesFunctionType(gridPart);
      }
 
      template<class DiscreteFunctionSpace>
      typename TraitsType::template ForcesFunctionalTraits<DiscreteFunctionSpace>::FunctionalType
      forcesFunctional(const DiscreteFunctionSpace& space) const
      {
        return ZeroFunctional<DiscreteFunctionSpace>(space);
      }
 
      DirichletBoundaryFunctionType dirichletBoundaryFunction(const GridPartType& gridPart) const
      {
        return DirichletBoundaryFunctionType(ZeroGridFunction<FunctionSpaceType, GridPartType>(gridPart),
                                             DirichletIndicatorType());
      }
 
      DirichletWeightFunctionType dirichletWeightFunction(const GridPartType& gridPart) const {
        typedef
          FractionGridFunction<typename FunctionSpaceType::ScalarFunctionSpaceType, GridPartType, 1L, 1UL>
          OneFunctionType;
        return DirichletWeightFunctionType(OneFunctionType(gridPart), DirichletIndicatorType());
      }
 
 
      NeumannBoundaryFunctionType neumannBoundaryFunction(const GridPartType& gridPart) const
      {
        return NeumannBoundaryFunctionType(ZeroGridFunction<FunctionSpaceType, GridPartType>(gridPart),
                                           NeumannIndicatorType());
      }
 
      DirichletIndicatorType dirichletIndicator() const {
        return DirichletIndicatorType();
      }
 
      NeumannIndicatorType neumannIndicator() const {
        return NeumannIndicatorType();
      }
    };
 
    template<class Model>
    struct ModelConstituents
    {
      typedef Model ModelType;
      typedef ModelInterface<ModelType> InterfaceType;
     protected:
      typedef
      Fem::DiscreteFunctionSpaceAdapter<typename InterfaceType::FunctionSpaceType,
                                        typename InterfaceType::GridPartType>
      DiscreteAdapterSpaceType;
      typedef
      typename InterfaceType::template ForcesFunctionalTraits<DiscreteAdapterSpaceType>::FunctionalType
      AdaptedForcesFunctionalType;
      enum {
        allDirichlet =
          std::is_same<typename InterfaceType::DirichletIndicatorType,
                       EntireBoundaryIndicatorType>::value
      };
     public:
      enum ConstituentFlags {
        hasFlux = InterfaceType::OperatorPartsType::hasFlux,
        hasSources = InterfaceType::OperatorPartsType::hasSources,
        hasRobinFlux = InterfaceType::OperatorPartsType::hasRobinFlux,
        hasBulkForces        = !ExpressionTraits<typename InterfaceType::BulkForcesFunctionType>::isZero,
        hasDirichletBoundary = !ExpressionTraits<typename InterfaceType::DirichletIndicatorType>::isZero,
        hasNeumannBoundary   = !allDirichlet && !ExpressionTraits<typename InterfaceType::NeumannIndicatorType>::isZero,
        hasRobinBoundary     = hasRobinFlux,
        hasForcesFunctional  = !ExpressionTraits<AdaptedForcesFunctionalType>::isZero,
      };
    };
 
    template<class Model>
    struct ModelConstituents<ModelInterface<Model> >
      : public ModelConstituents<Model>
    {};
 
    template<class Model>
    void printConstituents(const Model& model, std::ostream& os = std::cout)
    {
      typedef ModelConstituents<Model> HasType;
 
      os << "Constituents of model \"" << model.name() << "\"" << std::endl;
      os << "Flux:              " << std::to_string(HasType::hasFlux) << std::endl;
      os << "Sources:           " << std::to_string(HasType::hasSources) << std::endl;
      os << "RobinFlux:         " << std::to_string(HasType::hasRobinFlux) << std::endl;
      os << "BulkForces:        " << std::to_string(HasType::hasBulkForces) << std::endl;
      os << "ForcesFunctional:  " << std::to_string(HasType::hasForcesFunctional) << std::endl;
      os << "DirichletBoundary: " << std::to_string(HasType::hasDirichletBoundary) << std::endl;
      os << "NeumannBoundary:   " << std::to_string(HasType::hasNeumannBoundary) << std::endl;
      os << "RobinBoundary:     " << std::to_string(HasType::hasRobinBoundary) << std::endl;
      os << std::endl;
 
      // also list the structure flags:
      os << "isLinear:          " << std::to_string(Model::OperatorPartsType::isLinear) << std::endl;
      os << "isSymmetric:       " << std::to_string(Model::OperatorPartsType::isSymmetric) << std::endl;
      os << "isSemiDefinite:    " << std::to_string(Model::OperatorPartsType::isSemiDefinite) << std::endl;
 
      os << std::endl;
    }
 
    template<class DiscreteFunctionSpace, class LinearFunctional>
    struct FixedForcesFunctionalTraits
    {
      typedef DiscreteFunctionSpace DiscreteFunctionSpaceType;
     protected:
      typedef LinearFunctional FunctionalImplType;
      typedef typename DiscreteFunctionSpaceType::FunctionSpaceType FunctionSpaceType;
      typedef typename DiscreteFunctionSpaceType::GridPartType GridPartType;
      typedef
      Fem::DiscreteFunctionSpaceAdapter<FunctionSpaceType, GridPartType>
      AdapterSpaceType;
      enum {
        isSameSpace = std::is_same<typename FunctionalImplType::DiscreteFunctionSpaceType,
                                   DiscreteFunctionSpaceType>::value,
        isAdapterSpace = std::is_same<AdapterSpaceType, DiscreteFunctionSpaceType>::value,
        isZeroFunctional = ExpressionTraits<FunctionalImplType>::isZero
      };
      static_assert(isSameSpace || isAdapterSpace || isZeroFunctional,
                    "Behold! Invalid use of dummy classes!");
     public:
      // The hard requirement is that the type we construct here must
      // define DiscreteFunctionSpace as DiscreteFunctionSpaceType
      typedef
      typename std::conditional<isZeroFunctional, // remap to suitable zero
                                ZeroFunctional<DiscreteFunctionSpaceType>,
                                typename std::conditional<isSameSpace, // identity
                                                          FunctionalImplType,
                                                          NonZeroFunctional<DiscreteFunctionSpaceType> >::type >::type
      FunctionalType;
    };
 
 
 
  } // namespace ACFem
 
}  //Namespace Dune
 
#endif  // __DUNE_ACFEM_MODELINTERFACE_HH__