ellipticmodel.hh

#ifndef ELLIPTC_MODEL_HH
#define ELLIPTC_MODEL_HH
 
#include <cassert>
#include <cmath>
 
#include <dune/fem/solver/timeprovider.hh>
#include <dune/fem/io/parameter.hh>
 
#include "modelinterface.hh"
#include "probleminterface.hh"
#include "../functions/gridfunctionexpression.hh"
#include "../functions/gridfunctionwrapper.hh"
 
namespace Dune {
 
  namespace ACFem {
 
    // DiffusionModel
    // --------------
 
    template<class FunctionSpace, class GridPart, class Problem = ProblemInterface<FunctionSpace> >
    struct EllipticModel
      : public DefaultModel<EllipticModel<FunctionSpace, GridPart, Problem> >
    {
     private:
      typedef EllipticModel            ThisType;
      typedef DefaultModel<ThisType>   BaseType;
      typedef ModelInterface<ThisType> InterfaceType;
      typedef ModelTraits<ThisType>    TraitsType;
     public:
      typedef typename InterfaceType::FunctionSpaceType FunctionSpaceType;
      typedef typename InterfaceType::GridPartType GridPartType;
      typedef typename InterfaceType::EntityType EntityType;
      typedef typename InterfaceType::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;
 
      typedef Problem   ProblemType ;
 
      typedef typename ProblemType::DirichletBoundaryType DirichletBoundaryType;
      typedef typename ProblemType::NeumannBoundaryType   NeumannBoundaryType;
      typedef typename ProblemType::RightHandSideType     BulkForcesType;
      typedef typename ProblemType::ExactSolutionType     ExactSolutionType; // not part of the ModelInterface
 
      typedef typename InterfaceType::DirichletBoundaryFunctionType DirichletBoundaryFunctionType;
      typedef typename InterfaceType::DirichletWeightFunctionType DirichletWeightFunctionType;
      typedef typename InterfaceType::NeumannBoundaryFunctionType NeumannBoundaryFunctionType;
      typedef typename InterfaceType::BulkForcesFunctionType BulkForcesFunctionType;
      typedef typename TraitsType::ExactSolutionFunctionType ExactSolutionFunctionType; // not part of the ModelInterface
 
      typedef typename TraitsType::DirichletIndicatorType DirichletIndicatorType;
      typedef typename TraitsType::NeumannIndicatorType NeumannIndicatorType;
      typedef typename TraitsType::RobinIndicatorType RobinIndicatorType;
 
      EllipticModel(const ProblemType& problem)
        : problem_(problem),
          dirichletFunction_(problem_),
          neumannFunction_(problem_),
          bulkForcesFunction_(problem_),
          solutionFunction_(problem_),
          dirichletIndicator_(*this),
          neumannIndicator_(*this),
          robinIndicator_(*this),
          lastRobinIndication_(false)
      {}
 
      void setEntity(const EntityType& entity) const
      {
        // just here to remind us ...
      }
 
      bool setIntersection(const IntersectionType& intersection) const
      {
        lastRobinIndication_ = robinIndicator_.applies(intersection);
        return lastRobinIndication_;
      }
 
      template<class Entity, class Point>
      void flux(const Entity& entity,
                const Point &x,
                const RangeType& value,
                const JacobianRangeType& jacobian,
                JacobianRangeType& flux) const
      {
        // Default is a linear model
        linearizedFlux(value, jacobian, entity, x, value, jacobian, flux);
      }
 
      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
      {
        // The default model does not have a non-linearity.
        const DomainType xGlobal = entity.geometry().global(coordinate(x));
        problem_.secondOrderCoefficient(xGlobal, jacobian, flux);
      }
 
      template<class Entity, class Point>
      void source(const Entity& entity,
                  const Point &x,
                  const RangeType& value,
                  const JacobianRangeType& jacobian,
                  RangeType& result) const
      {
        linearizedSource(value, jacobian, entity, x, value, jacobian, result);
      }
 
      template<class Entity, class Point>
      void linearizedSource(const RangeType& uBar,
                            const JacobianRangeType& DuBar,
                            const Entity& entity,
                            const Point &x,
                            const RangeType& value,
                            const JacobianRangeType& jacobian,
                            RangeType& result) const
      {
        const DomainType xGlobal = entity.geometry().global(coordinate(x));
        result = 0;
        if (problem_.has(ProblemType::ZeroOrderCoefficient)) {
          RangeType tmp;
          problem_.zeroOrderCoefficient(xGlobal, value, tmp);
          result += tmp;
        }
        if (problem_.has(ProblemType::FirstOrderCoefficient)) {
          RangeType tmp;
          problem_.firstOrderCoefficient(xGlobal, value, jacobian, tmp);
          result += tmp;
        }
      }
 
      template<class Intersection, class Point>
      void robinFlux(const Intersection& intersection,
                     const Point& x,
                     const DomainType& unitOuterNormal,
                     const RangeType& value,
                     RangeType& result) const
      {
        linearizedRobinFlux(value, intersection, x, unitOuterNormal, value, result);
      }
 
      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
      {
        if (!lastRobinIndication_) {
          result = 0;
          return;
        }
        const DomainType xGlobal = intersection.inside()->geometry().global(coordinate(x));
 
        problem_.robinData(xGlobal, value, result);
      }
 
      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
      {
        const DomainType xGlobal = entity.geometry().global(coordinate(x));
 
        result = 0;
        if (problem_.has(ProblemType::SecondOrderCoefficient)) {
          problem_.secondOrderCoefficient(xGlobal, jacobian, hessian, result);
        }
      }
 
     public:
      BulkForcesFunctionType bulkForcesFunction(const GridPartType& gridPart) const
      {
        typedef BulkForcesFunctionType FunctionType;
 
        return FunctionType("Bulk Forces", bulkForcesFunction_, gridPart);
      }
 
      DirichletBoundaryFunctionType dirichletBoundaryFunction(const GridPartType& gridPart) const
      {
        typedef typename TraitsType::DirichletBoundaryGridFunctionType GridFunctionType;
        typedef DirichletBoundaryFunctionType FunctionType;
 
        auto gridFunction = GridFunctionType("Dirichlet Boundary Values", dirichletFunction_, gridPart);
 
        return FunctionType(gridFunction, dirichletIndicator_);
      }
 
      DirichletWeightFunctionType dirichletWeightFunction(const GridPartType& gridPart) const
      {
        typedef DirichletWeightFunctionType FunctionType;
 
        auto gridFunction = oneFunction(dirichletBoundaryFunction(gridPart));
 
        return FunctionType(gridFunction, dirichletIndicator_);
      }
 
      NeumannBoundaryFunctionType
      neumannBoundaryFunction(const GridPartType& gridPart) const
      {
        typedef typename TraitsType::NeumannBoundaryGridFunctionType GridFunctionType;
        typedef NeumannBoundaryFunctionType FunctionType;
 
        auto gridFunction = GridFunctionType("Neumann Boundary Values", neumannFunction_, gridPart);
 
        return FunctionType(gridFunction, neumannIndicator_);
      }
 
      ExactSolutionFunctionType
      exactSolutionFunction(const GridPartType& gridPart) const
      {
        typedef ExactSolutionFunctionType FunctionType;
 
        return FunctionType("Exact Solution", solutionFunction_, gridPart);
      }
 
      DirichletIndicatorType dirichletIndicator() const {
        return dirichletIndicator_;
      }
 
      NeumannIndicatorType neumannIndicator() const {
        return neumannIndicator_;
      }
 
      RobinIndicatorType robinIndicator() const {
        return robinIndicator_;
      }
 
      const ProblemType& problem() const { return problem_; }
 
     protected:
      const ProblemType& problem_; 
      DirichletBoundaryType dirichletFunction_;
      NeumannBoundaryType   neumannFunction_;
      BulkForcesType        bulkForcesFunction_;
      ExactSolutionType     solutionFunction_;
 
      DirichletIndicatorType dirichletIndicator_;
      NeumannIndicatorType   neumannIndicator_;
      RobinIndicatorType     robinIndicator_;
      mutable bool lastRobinIndication_;
    };
 
    template<class FunctionSpace, class GridPart, class Problem>
    struct ModelTraits<EllipticModel<FunctionSpace, GridPart, Problem> >
      : public DefaultModelTraits<FunctionSpace, GridPart>
    {
     private:
      typedef EllipticModel<FunctionSpace, GridPart, Problem> ModelType;
      typedef Problem   ProblemType;
 
      typedef typename ProblemType::DirichletBoundaryType DirichletBoundaryType;
      typedef typename ProblemType::NeumannBoundaryType   NeumannBoundaryType;
      typedef typename ProblemType::RightHandSideType     BulkForcesType;
      typedef typename ProblemType::ExactSolutionType     ExactSolutionType;
 
      typedef typename FunctionSpace::DomainType DomainType;
     public:
      enum {
        isLinear = true,      // The ProblemInterface only models linear operators
        isSymmetric = true,   // really?
        isSemiDefinite = true // really?
      };
 
      enum ConstituentFlags {
        hasFlux = true,    
        hasSources = true, 
      };
 
      typedef FunctionSpace FunctionSpaceType;
      typedef typename FunctionSpace::ScalarFunctionSpaceType ScalarFunctionSpaceType;
      typedef GridPart GridPartType;
 
      struct DirichletIndicator
        : DefaultBoundaryIndicator<DirichletIndicator>
      {
        DirichletIndicator(const ModelType& model) : model_(model)
        {}
 
        template<class Intersection>
        bool applies(const Intersection& intersection) const
        {
          const int bndId = intersection.boundaryId();
          const DomainType center = intersection.geometry().center();
 
          return model_.problem().isDirichletSegment(bndId, center);
        }
       private:
        const ModelType& model_;
      };
 
      typedef DirichletIndicator DirichletIndicatorType;
 
      struct NeumannIndicator
        : DefaultBoundaryIndicator<NeumannIndicator>
      {
        NeumannIndicator(const ModelType& model) : model_(model)
        {}
 
        template<class Intersection>
        bool applies(const Intersection& intersection) const
        {
          const int bndId = intersection.boundaryId();
          const DomainType center = intersection.geometry().center();
 
          return model_.problem().isNeumannSegment(bndId, center);
        }
       private:
        const ModelType& model_;
      };
      typedef NeumannIndicator NeumannIndicatorType;
 
      struct RobinIndicator
        : DefaultBoundaryIndicator<RobinIndicator>
      {
        RobinIndicator(const ModelType& model) : model_(model)
        {}
 
        template<class Intersection>
        bool applies(const Intersection& intersection) const
        {
          const int bndId = intersection.boundaryId();
          const DomainType center = intersection.geometry().center();
 
          return model_.problem().isRobinSegment(bndId, center);
        }
       private:
        const ModelType& model_;
      };
      typedef RobinIndicator RobinIndicatorType;
 
      typedef GridFunctionAdapter<BulkForcesType, GridPartType> BulkForcesFunctionType;
      typedef GridFunctionAdapter<ExactSolutionType, GridPartType> ExactSolutionFunctionType;
 
      typedef GridFunctionAdapter<DirichletBoundaryType, GridPartType> DirichletBoundaryGridFunctionType;
      typedef
      BoundarySupportedFunction<DirichletBoundaryGridFunctionType, DirichletIndicatorType>
      DirichletBoundaryFunctionType;
      typedef
      BoundarySupportedFunction<FractionGridFunction<ScalarFunctionSpaceType, GridPartType, 1L, 1UL>,
                                DirichletIndicatorType>
      DirichletWeightFunctionType;
 
      typedef GridFunctionAdapter<NeumannBoundaryType, GridPartType> NeumannBoundaryGridFunctionType;
      typedef
      BoundarySupportedFunction<NeumannBoundaryGridFunctionType, NeumannIndicatorType>
      NeumannBoundaryFunctionType;
    };
 
 
 
  } // ACFem
 
} // Dune
 
#endif // #ifndef ELLIPTC_MODEL_HH