modeldirichletvalues.hh

#ifndef __DUNE_ACFEM_FUNCTIONS_MODULES_MODELDIRICHLETVALUES_HH__
#define __DUNE_ACFEM_FUNCTIONS_MODULES_MODELDIRICHLETVALUES_HH__
 
#include <dune/fem/function/localfunction/bindable.hh>
#include "../../models/modelfacade.hh"
#include "../../models/modelboundaryindicators.hh"
 
#include "zero.hh"
 
namespace Dune {
 
  namespace ACFem {
 
    namespace GridFunction {
 
      template<class Model, class GridPart>
      class ModelDirichletValues
        : public Fem::BindableGridFunction<GridPart, typename ModelTraits<Model>::RangeType>
      {
        static_assert(IsPDEModel<Model>::value,
                      "Model better would be a proper PDE-model ...");
        using BaseType = Fem::BindableGridFunction<GridPart, typename ModelTraits<Model>::RangeType>;
       public:
        using GridPartType = GridPart;
        using ModelType = std::decay_t<Model>;
        using ModelTraits = ModelIntrospection::Traits<Model>;
        using typename BaseType::EntityType;
        using typename BaseType::RangeFieldType;
        using RealType = typename Dune::FieldTraits<RangeFieldType>::real_type;
        using typename BaseType::DomainType;
        using typename BaseType::RangeType;
        using typename BaseType::JacobianRangeType;
        using HessianRangeType = typename ModelType::HessianRangeType; // NOT from BaseType
        static constexpr std::size_t dimRange = ModelType::dimRange;
        static constexpr std::size_t dimDomain = ModelType::dimDomain;
        static constexpr std::size_t dirichlet = ModelIntrospection::dirichlet;
        static constexpr std::size_t linearizedDirichlet = ModelIntrospection::linearizedDirichlet;
        static constexpr bool precomputeJacobian =
          ModelTraits::template IsPiecewiseConstant<linearizedDirichlet>::value
          &&
          ModelTraits::template IsAffineLinear<dirichlet>::value;
        using DirichletJacobianRangeType = FieldMatrix<RangeFieldType, dimRange, dimRange>;
 
        ModelDirichletValues(Model&& model,
                             const GridPartType& gridPart,
                             const std::string& name = "")
          : BaseType(gridPart)
          , model_(model)
          , dirichlet_(model_)
          , linearizedDirichlet_(model_)
          , name_(name == "" ? "Dirichlet("+model_.name()+")" : name)
          , tolerance_(Fem::Parameter::getValue<RealType>("acfem.dirichlet.newton.tolerance", 2.0*std::numeric_limits<RangeFieldType>::epsilon()))
          , maxIter_(Fem::Parameter::getValue<RealType>("acfem.dirichlet.newton.iterationLimit", 100))
        {}
 
        ModelDirichletValues(const ModelDirichletValues& other)
          : BaseType(other)
          , model_(other.model_)
          , dirichlet_(model_)
          , linearizedDirichlet_(model_)
          , name_(other.name_)
          , tolerance_(other.tolerance_)
          , maxIter_(other.maxIter_)
        {}
 
        ModelDirichletValues(ModelDirichletValues&& other)
          : BaseType(other)
          , model_(std::move(other.model_))
          , dirichlet_(model_)
          , linearizedDirichlet_(model_)
          , name_(std::move(other.name_))
          , tolerance_(other.tolerance_)
          , maxIter_(other.maxIter_)
        {}
 
        using BaseType::entity;
        using BaseType::gridPart;
 
        const std::string& name() const { return name_; }
 
        // New  Bindable interface
        void bind(const EntityType& entity)
        {
          BaseType::bind(entity);
          model_.bind(BaseType::entity());
        }
 
        // New  Bindable interface
        void unbind()
        {
          model_.unbind();
          BaseType::unbind();
        }
 
        template<class Intersection>
        auto classifyIntersection(const Intersection& intersection)
        {
          auto supported = model_.classifyBoundary(intersection);
          if constexpr (precomputeJacobian) {
            // FIXME: will be wrong when working on a manifold
            jacobian_ = computeJacobian(RangeType(), quadraturePoint(DomainType()));
          }
          return supported;
        }
 
        int order() const
        {
          return std::numeric_limits<int>::max();
        }
 
        template<class Point,
                 std::enable_if_t<(IsQuadraturePoint<Point>::value || IsFieldVector<Point>::value), int> = 0>
        auto evaluate(Point& x) const
        {
          /*
            G(u) = u*|u|^2 - g
            G'(u)v = 2 (uv) u + |u|^2  v,
 
            G(u) = alpha u - g
            G'(u)v = alpha v
          */
 
          // Some start value. In case of non-convergence one could try
          // to modify it ...
          auto values = RangeType(0.5);
 
          for (int k = 0; k <= maxIter_; ++k) {
 
            assert(k < maxIter_); // rather throw or something
 
            auto G = dirichlet_(x, values);
 
            RangeType d;
            if constexpr (precomputeJacobian) {
                jacobian_.solve(d, G);
              } else {
              computeJacobian(values, x).solve(d, G);
            }
 
#if 0
            static_assert(ModelTraits::template IsAffineLinear<dirichlet>::value,
                          "non-linear Dirichlet values");
#endif
 
            if constexpr (ModelTraits::template IsAffineLinear<dirichlet>::value) {
              values -= d;
              break;
            }
 
            if (d.two_norm2() < tolerance_ * tolerance_) {
              break;
            }
 
            values -= d;
 
          }
          return values;
        }
 
        template<class Point,
                 std::enable_if_t<(IsQuadraturePoint<Point>::value || IsFieldVector<Point>::value), int> = 0>
        void evaluate(const Point& point, RangeType& result) const
        {
          result = evaluate(point);
        }
 
        template<class Point,
                 std::enable_if_t<(IsQuadraturePoint<Point>::value || IsFieldVector<Point>::value), int> = 0>
        auto jacobian(const Point& x) const
        {
#if !__DUNE_ACFEM_MAKE_CHECK__
          static_assert(sizeof(Point) != 0, "It does not make sense to differentiate Dirichlet values.");
#endif
          return JacobianRangeType(0);
        }
 
        template<class Point,
                 std::enable_if_t<(IsQuadraturePoint<Point>::value || IsFieldVector<Point>::value), int> = 0>
        void jacobian(const Point& x, JacobianRangeType& result) const
        {
#if !__DUNE_ACFEM_MAKE_CHECK__
          static_assert(sizeof(Point) != 0, "It does not make sense to differentiate Dirichlet values.");
#endif
          result = jacobian(x);
        }
 
        // hessian of local function
        template<class Point,
                 std::enable_if_t<(IsQuadraturePoint<Point>::value || IsFieldVector<Point>::value), int> = 0>
        auto hessian(const Point& x) const
        {
#if !__DUNE_ACFEM_MAKE_CHECK__
          static_assert(sizeof(Point) != 0, "It does not make sense to differentiate Dirichlet values.");
#endif
          return HessianRangeType(0);
        }
 
        // hessian of local function
        template<class Point,
                 std::enable_if_t<(IsQuadraturePoint<Point>::value || IsFieldVector<Point>::value), int> = 0>
        void hessian(const Point& x, HessianRangeType& ret) const
        {
#if !__DUNE_ACFEM_MAKE_CHECK__
          static_assert(sizeof(Point) != 0, "It does not make sense to differentiate Dirichlet values.");
#endif
          ret = hessian(x);
        }
 
        const ModelType& model() const { return model_; }
        ModelType& model() { return model_; }
 
       protected:
        template<class Quadrature>
        auto computeJacobian(const RangeType& values, const QuadraturePoint<Quadrature>& x) const
        {
          DirichletJacobianRangeType jacobian;
 
          forLoop<dimRange>([&](auto j) {
              using J = decltype(j);
              auto ej = kroneckerDelta<J::value>(IndexSequence<dimRange>{});
              auto tmp = linearizedDirichlet_(values, x, ej);
              for (unsigned i = 0; i < dimRange; ++i) {
                jacobian[i][J::value] = tmp[i];
              }
            });
          return jacobian;
        }
 
       protected:
        Model model_;
        ModelClosureMethod<ModelType, dirichlet> dirichlet_;
        ModelClosureMethod<ModelType, linearizedDirichlet> linearizedDirichlet_;
        std::string name_;
        RealType tolerance_;
        int maxIter_;
        DirichletJacobianRangeType jacobian_;
      };
 
      template<class Model, class GridPart,
               std::enable_if_t<EffectiveModelTraits<Model>::template Exists<PDEModel::dirichlet>::value, int> = 0>
      auto modelDirichletValues(Model&& model,
                                const GridPart& gridPart,
                                const std::string& name = "")
      {
        static_assert(IsPDEModel<Model>::value, "Model better would be a proper PDE-model ...");
        return ModelDirichletValues<Model, GridPart>(std::forward<Model>(model), gridPart, name);
      }
 
      template<class Model, class GridPart,
               std::enable_if_t<!EffectiveModelTraits<Model>::template Exists<PDEModel::dirichlet>::value, int> = 0>
      auto modelDirichletValues(Model&& model,
                                const GridPart& gridPart,
                                const std::string& name = "")
      {
        static_assert(IsPDEModel<Model>::value, "Model better would be a proper PDE-model ...");
        return GridFunction::zeros<typename std::decay_t<Model>::FunctionSpaceType>(gridPart);
      }
 
 
 
    } // NS GridFunction::
 
    using GridFunction::modelDirichletValues;
 
  } // ACFem::
 
} // Dune::
 
#endif // __DUNE_ACFEM_FUNCTIONS_MODULES_MODELDIRICHLETVALUES_HH__