ellipticoperator.hh

#ifndef __ELLIPTIC_OPERATOR_HH__
#define __ELLIPTIC_OPERATOR_HH__
 
#include <optional>
#include <dune/common/fmatrix.hh>
 
#include <dune/fem/operator/common/differentiableoperator.hh>
#include <dune/fem/operator/common/temporarylocalmatrix.hh>
#include <dune/fem/operator/common/stencil.hh>
 
#include "constraints/emptyblockconstraints.hh"
#include "../models/modelfacade.hh"
#include "../common/quadrature.hh"
#include "../algorithms/operatorselector.hh"
 
namespace Dune {
 
  namespace ACFem {
 
    template<class Model,
             class DomainFunction,
             class RangeFunction = DomainFunction,
             class Constraints = EmptyBlockConstraints<typename RangeFunction::DiscreteFuncionSpaceType>,
             template<class GridPart> class QuadratureTraits = DefaultQuadratureTraits>
    class EllipticOperator
      : public virtual std::conditional<EffectiveModelTraits<Model>::isLinear,
                                        Fem::LinearOperator<DomainFunction, RangeFunction>,
                                        Fem::Operator<DomainFunction, RangeFunction> >::type
    {
      typedef typename std::conditional<
        EffectiveModelTraits<Model>::isLinear,
        Fem::LinearOperator<DomainFunction, RangeFunction>,
        Fem::Operator<DomainFunction, RangeFunction> >::type
      BaseType;
     public:
      using typename BaseType::DomainFunctionType;
 
      using typename BaseType::RangeFunctionType;
 
      using typename BaseType::DomainFieldType;
 
      using typename BaseType::RangeFieldType;
 
      using DiscreteFunctionType = RangeFunctionType;
 
      static_assert(std::is_base_of<Fem::IsDiscreteFunction, RangeFunctionType>::value,
                    "The RangeFunctionType has to be a discrete function");
 
      static_assert(std::is_base_of<Fem::HasLocalFunction, DomainFunctionType>::value,
                    "The DomainFunctionType has to be provide a local function");
 
      using ModelTraits = EffectiveModelTraits<Model>;
      using ModelType = EffectiveModel<Model>;
      using ModelClosure = ModelFacade<Model>;
      using ConstraintsOperatorType = std::decay_t<Constraints>;
     protected:
      typedef typename DomainFunctionType::FunctionSpaceType DomainFunctionSpaceType;
      typedef typename RangeFunctionType::FunctionSpaceType RangeFunctionSpaceType;
      enum {
        dimRange = RangeFunctionSpaceType::dimRange,
        dimDomain = RangeFunctionSpaceType::dimDomain
      };
      typedef typename RangeFunctionSpaceType::DomainType DomainType;
      typedef typename RangeFunctionSpaceType::RangeType RangeType;
      typedef typename RangeFunctionSpaceType::JacobianRangeType JacobianRangeType;
 
      // how's about
      static_assert(std::is_same<DomainType, typename DomainFunctionSpaceType::DomainType>::value,
                    "Domain and range functions need to live on the same domain");
 
      // how's about divergence constraints?
      static_assert(std::is_same<DomainFunctionSpaceType, RangeFunctionSpaceType>::value,
                    "Domain and range function spaces have to be the same");
 
      typedef typename RangeFunctionType::DiscreteFunctionSpaceType DiscreteFunctionSpaceType;
      typedef typename RangeFunctionType::LocalFunctionType LocalFunctionType;
 
      typedef typename Fem::ConstLocalFunction<DomainFunctionType> LocalDomainFunctionType;
 
      typedef typename DiscreteFunctionSpaceType::GridPartType GridPartType;
      typedef typename GridPartType::IntersectionIteratorType IntersectionIteratorType;
      typedef typename IntersectionIteratorType::Intersection IntersectionType;
 
      typedef QuadratureTraits<GridPartType> QuadratureTraitsType;
      typedef typename QuadratureTraitsType::BulkQuadratureType QuadratureType;
      typedef typename QuadratureTraitsType::BulkMassQuadratureType MassQuadratureType;
      typedef typename QuadratureTraitsType::FaceQuadratureType FaceQuadratureType;
      typedef typename QuadratureTraitsType::FaceMassQuadratureType FaceMassQuadratureType;
 
      template<bool conforming>
      using IntersectionQuadrature = typename QuadratureTraitsType::template IntersectionQuadrature<conforming>;
 
      typedef IntersectionQuadrature<true> BndryQuadratureType;
      typedef typename QuadratureTraitsType::template IntersectionMassQuadrature<true> BndryMassQuadratureType;
      static constexpr bool hasFlux = ModelTraits::template Exists<PDEModel::flux>::value;
      static constexpr bool hasSources = ModelTraits::template Exists<PDEModel::source>::value;
      static constexpr bool hasRobin = ModelTraits::template Exists<PDEModel::robinFlux>::value;
      static constexpr bool hasSingular = ModelTraits:: template Exists<PDEModel::singularFlux>::value;
      static constexpr bool isAffineLinear = ModelTraits::isLinear;
      static constexpr bool hasMassQuadrature = QuadratureTraitsType::hasMassQuadrature;
      static constexpr bool hasFaceMassQuadrature = QuadratureTraitsType::hasMassQuadrature;
 
      using DGSelectorType = OperatorTypeSelector<DiscreteFunctionSpaceType, ModelType>;
      static constexpr bool useDG = DGSelectorType::useDG;
      static constexpr DGFlavour dgFlavour = DGSelectorType::dgFlavour;
      static constexpr int sFormDG = DGSelectorType::sFormDG;
 
     public:
      template<class ConstraintsArg,
               std::enable_if_t<std::is_constructible<Constraints, ConstraintsArg>::value, int> = 0>
      EllipticOperator(const ModelType& model, ConstraintsArg&& constraints)
        : model_(model)
        , constraintsOperator_(constraints)
        , penalty_(DGSelectorType::dgPenalty())
      {}
 
      EllipticOperator(const ModelType &model)
        : model_(model)
        , penalty_(DGSelectorType::dgPenalty())
      {}
 
      EllipticOperator(const EllipticOperator& other)
        : model_(other.model_)
        , constraintsOperator_(std::move(other.constraintsOperator_))
        , penalty_(other.penalty_)
      {}
 
      EllipticOperator(EllipticOperator&& other)
        : model_(std::move(other.model_))
        , constraintsOperator_(std::move(other.constraintsOperator_))
        , penalty_(std::move(other.penalty_))
      {}
 
      virtual void operator() (const DomainFunctionType &u, RangeFunctionType &w) const;
 
      virtual bool symmetric() const
      {
        return ExpressionTraits<ModelType>::isSymmetric;
      }
 
      virtual bool positiveDefinite() const
      {
        return ExpressionTraits<ModelType>::Definiteness::isSemiPositive;
      }
 
     protected:
      const ConstraintsOperatorType& constraints() const {
        return constraintsOperator_;
      }
      const ModelType& model() const {
        return model_;
      }
      ModelType& model() {
        return model_;
      }
      const ModelClosure& closure() const {
        return model_;
      }
      ModelClosure& closure() {
        return model_;
      }
      const double penalty() const {
        return penalty_;
      }
 
     private:
      template<bool conforming>
      void computeInteriorSkeletonContributions(
        const GridPartType& gridPart,
        const IntersectionType& intersection,
        const RangeFieldType& penalty,
        const LocalDomainFunctionType& uLocal,
        const LocalDomainFunctionType& uLocalNb,
        const ModelClosure& model,
        const ModelClosure& modelNb,
        LocalFunctionType& wLocal) const;
 
     protected:
      ModelClosure model_;
      Constraints constraintsOperator_;
      const double penalty_;
    };
 
    template<class JacobianOperator,
             class Model,
             class Constraints = EmptyBlockConstraints<typename JacobianOperator::RangeFunctionType::DiscreteFunctionSpaceType>,
        template<class> class QuadratureTraits = DefaultQuadratureTraits>
    class DifferentiableEllipticOperator
      : public EllipticOperator<Model,
            typename JacobianOperator::DomainFunctionType,
            typename JacobianOperator::RangeFunctionType,
            Constraints, QuadratureTraits>,
        public virtual Fem::DifferentiableOperator<JacobianOperator>,
        public virtual std::conditional<EffectiveModelTraits<Model>::isLinear,
            Fem::AssembledOperator<typename JacobianOperator::DomainFunctionType,
                typename JacobianOperator::RangeFunctionType>,
            Fem::Operator<typename JacobianOperator::DomainFunctionType,
                typename JacobianOperator::RangeFunctionType> >::type
    {
      typedef
      EllipticOperator<Model,
          typename JacobianOperator::DomainFunctionType,
          typename JacobianOperator::RangeFunctionType,
          Constraints, QuadratureTraits>
      BaseType;
     public:
      using JacobianOperatorType = JacobianOperator;
 
      using typename BaseType::ModelType;
      using typename BaseType::ModelClosure;
      using typename BaseType::ModelTraits;
      using typename BaseType::ConstraintsOperatorType;
 
      using typename BaseType::DomainFunctionType;
      using typename BaseType::RangeFunctionType;
 
      using typename BaseType::DomainFieldType;
      using typename BaseType::RangeFieldType;
 
      static_assert(std::is_base_of<Fem::IsDiscreteFunction, DomainFunctionType>::value,
                    "For matrix assembly the DomainFunctionType has to be a discrete function");
 
      // TODO: allow for different discrete function space types.
 
     protected:
      using BaseType::dimDomain; // slighlty misleading name ...
      using BaseType::dimRange; // slighlty misleading name ...
 
      using typename BaseType::DomainFunctionSpaceType;
      using typename BaseType::LocalDomainFunctionType;
 
      using typename BaseType::DiscreteFunctionSpaceType;
      using typename BaseType::LocalFunctionType;
 
      using typename BaseType::GridPartType;
      using typename BaseType::IntersectionType;
 
      using typename BaseType::RangeType;
      using typename BaseType::JacobianRangeType;
 
      using typename BaseType::QuadratureType;
      using typename BaseType::MassQuadratureType;
      using typename BaseType::BndryQuadratureType;
      using typename BaseType::BndryMassQuadratureType;
 
      template<bool conforming>
      using IntersectionQuadrature = typename BaseType::template IntersectionQuadrature<conforming>;
 
      using BaseType::hasFlux;
      using BaseType::hasSources;
      using BaseType::hasRobin;
      using BaseType::hasSingular;
      using BaseType::isAffineLinear;
      using BaseType::hasMassQuadrature;
      using BaseType::hasFaceMassQuadrature;
      using BaseType::useDG;
      using BaseType::dgFlavour;
      using BaseType::sFormDG;
 
      // new typedefs
      using DiscreteDomainFunctionSpaceType = typename DomainFunctionType::DiscreteFunctionSpaceType;
 
      static_assert(std::is_same<DiscreteDomainFunctionSpaceType, DiscreteFunctionSpaceType>::value,
                    "Unfortunately the discrete function space type of domain and range need to coincide. Sorry for that.");
 
      using  ElementMatrixType = Fem::TemporaryLocalMatrix<DiscreteDomainFunctionSpaceType, DiscreteFunctionSpaceType>;
 
     public:
//      using BaseType::EllipticOperator;
      //    using EllipticOperator::EllipticOperator;
      using EllipticOperator<
        Model,
        typename JacobianOperator::DomainFunctionType,
        typename JacobianOperator::RangeFunctionType,
        Constraints, QuadratureTraits>::EllipticOperator;
 
      using BaseType::operator();
      using BaseType::symmetric;
      using BaseType::positiveDefinite;
 
      void jacobian(const DomainFunctionType &u, JacobianOperatorType &jOp) const;
 
     private:
      // storage for basis functions, mutable in order to keep jacobian() const
      mutable std::vector<RangeType> phi_;
      mutable std::vector<JacobianRangeType> dphi_;
 
      mutable std::vector<RangeType> phiNb_;
      mutable std::vector<JacobianRangeType> dphiNb_;
 
      void reserveAssemblyStorage(size_t maxNumBasisFunctions) const
      {
        phi_.resize(maxNumBasisFunctions);
        dphi_.resize(maxNumBasisFunctions);
 
        phiNb_.resize(maxNumBasisFunctions);
        dphiNb_.resize(maxNumBasisFunctions);
      }
 
      void releaseAssemblyStorage() const
      {
        phi_.resize(0);
        dphi_.resize(0);
 
        phiNb_.resize(0);
        dphiNb_.resize(0);
      }
 
      template<bool conforming>
      void computeInteriorSkeletonContributions(
        const GridPartType& gridPart,
        const IntersectionType& intersection,
        const RangeFieldType& penalty,
        const LocalDomainFunctionType& uLocal,
        const LocalDomainFunctionType& uLocalNb,
        const ModelClosure& model,
        const ModelClosure& modelNb,
        ElementMatrixType& elementMatrix,
        ElementMatrixType& elementMatrixNb) const;
 
     protected:
      using BaseType::model_;
      using BaseType::model;
      using BaseType::closure;
      using BaseType::constraints;
      using BaseType::penalty;
    };
 
    // Implementation of application operator.
    template<class Model,
        class DomainFunction, class RangeFunction,
        class Constraints,
        template<class> class QuadratureTraits>
    void EllipticOperator<Model, DomainFunction, RangeFunction, Constraints, QuadratureTraits>
    ::operator()(const DomainFunctionType &u, RangeFunctionType &w) const
    {
      // clear destination
      w.clear();
 
      // get discrete function space
      const auto& dfSpace = w.space();
      const auto& gridPart = dfSpace.gridPart();
 
      // possibly update
      // dirichlet constraints are done weakly in DG
      // as we are not (necessarily) in a lagrange space
      // if that is the case, one could switch to setting them again hard
      constraints().update();
 
      // obtain local functions, note that unfortunately
      // not all HasLocalFunction objects have a localFunction(void)
      // method, hence not "auto"
      LocalDomainFunctionType uLocal(u);
      LocalDomainFunctionType uLocalNb(u);
      LocalFunctionType wLocal(w);
 
      // mutable copy of the model-closure
      auto model = closure();
 
      // For DG we need a copy of the model which lives on the neighbour
      auto modelNb = closure();
 
      // iterate over grid
      const auto end = dfSpace.end();
      for (auto it = dfSpace.begin(); it != end; ++it) {
 
        // get entity (here element)
        const auto& entity = *it;
        // get elements geometry
        const auto& geometry = entity.geometry();
 
        // get local representation of the discrete functions
        uLocal.bind(entity);
 
        wLocal.init(entity);
 
        // call per-element initializer for the integral contributions
        model.bind(entity);
 
        // element integrals, compute all bulk contributions
        if constexpr((hasSources && !hasMassQuadrature) || hasFlux) {
          // obtain quadrature order
          const int quadOrder = uLocal.order() + wLocal.order();
 
          // 2nd and 0 order in one quad-loop
          QuadratureType quadrature(entity, quadOrder);
          const size_t numQuadraturePoints = quadrature.nop();
          for (size_t pt = 0; pt < numQuadraturePoints; ++pt) {
 
            const typename QuadratureType::CoordinateType &x = quadrature.point(pt);
            const double weight = quadrature.weight(pt) * geometry.integrationElement(x);
 
            RangeType phiKern(0);
            JacobianRangeType dPhiKern(0);
 
            // compute the second order contribution
            JacobianRangeType du;
            uLocal.jacobian(quadrature[pt], du);
 
            // compute the source contribution
            RangeType vu;
            uLocal.evaluate(quadrature[pt], vu);
 
            if constexpr (hasSources && !hasMassQuadrature && hasFlux) {
              RangeType phiKern = weight * model.source(quadrature[pt], vu, du);
              JacobianRangeType dPhiKern = weight * model.flux(quadrature[pt], vu, du);
              wLocal.axpy(quadrature[pt], phiKern, dPhiKern);
            } else if constexpr (hasFlux) {
              JacobianRangeType dPhiKern = weight * model.flux(quadrature[pt], vu, du);
              wLocal.axpy(quadrature[pt], dPhiKern);
            } else {
              RangeType phiKern = weight * model.source(quadrature[pt], vu, du);
              wLocal.axpy(quadrature[pt], phiKern);
            }
          }  // end quad loop
        } // end ordinary bulk contributions
 
        // Special code-path to allow e.g. for mass-lumping
        if constexpr (hasSources && hasMassQuadrature) {
          // obtain quadrature order
          const int quadOrder = uLocal.order() + wLocal.order();
 
          MassQuadratureType quadrature(entity, quadOrder);
 
          const size_t numQuadraturePoints = quadrature.nop();
          for (size_t pt = 0; pt < numQuadraturePoints; ++pt) {
 
            const typename MassQuadratureType::CoordinateType &x = quadrature.point(pt);
            const double weight = quadrature.weight(pt) * geometry.integrationElement(x);
 
            // compute the second order contribution
            JacobianRangeType du;
            uLocal.jacobian(quadrature[pt], du);
 
            // compute the source contribution
            RangeType vu;
            uLocal.evaluate(quadrature[pt], vu);
 
            RangeType phiKern = weight * model.source(quadrature[pt], vu, du);
 
            // add to local function
            wLocal.axpy(quadrature[pt], phiKern);
          }
        } // end of mass-lumping
 
        // end bulk block
 
        if constexpr (useDG || hasRobin || hasSingular) { // add skeleton integrals as necessary
 
          const auto area = useDG ? entity.geometry().volume() : 0.0;
 
          const auto end = gridPart.iend(entity);
          for (auto it = gridPart.ibegin(entity); it != end; ++it ) { // looping over intersections
 
            const auto& intersection = *it;
            if (useDG && intersection.neighbor()) {
              const auto neighbor = intersection.outside() ;
              const auto& intersectionGeometry = intersection.geometry();
 
              // compute penalty factor
              const auto intersectionArea = intersectionGeometry.volume();
              //beta is still missing the diffusion coefficient!!
              const auto beta = penalty() * intersectionArea / std::min(area, neighbor.geometry().volume());
 
              modelNb.bind(neighbor);   // model in outside element
              uLocalNb.bind(neighbor); // local u on outisde element
 
              if (intersection.conforming()) {
                computeInteriorSkeletonContributions<true>(
                  gridPart, intersection, beta, uLocal, uLocalNb, model, modelNb, wLocal);
              } else {
                computeInteriorSkeletonContributions<false>(
                  gridPart, intersection, beta, uLocal, uLocalNb, model, modelNb, wLocal);
              }
 
              // unbind neighbour things
              uLocalNb.unbind();
              modelNb.unbind();
 
            } else if (intersection.boundary() && (hasRobin || hasSingular)) {
 
              const auto bdCond = model.classifyBoundary(intersection);
 
              if (bdCond.first && !bdCond.second.all()) {
 
                int quadOrder = 3*dfSpace.order();
                BndryMassQuadratureType intersectionQuad(gridPart, intersection, quadOrder);
                const auto &bndryQuad = intersectionQuad.inside();
                const int numQuadraturePoints = bndryQuad.nop();
 
                for (int pt = 0; pt < numQuadraturePoints; ++pt) {
                  // quadrature weight
                  const double weight = bndryQuad.weight(pt);
 
                  // Get the un-normalized outer normal
                  DomainType integrationNormal
                    = intersection.integrationOuterNormal(bndryQuad.localPoint(pt));
 
                  const double integrationElement = integrationNormal.two_norm();
 
                  integrationNormal /= integrationElement;
 
                  RangeType vu;
                  uLocal.evaluate(bndryQuad[pt], vu);
                  JacobianRangeType du;
                  uLocal.jacobian(bndryQuad[pt], du);
 
                  if constexpr (hasRobin && hasSingular) {
                    RangeType phiKern = weight * integrationElement * model.robinFlux(bndryQuad[pt], integrationNormal, vu, du);
                    JacobianRangeType dphiKern = weight * integrationElement * model.singularFlux(bndryQuad[pt], integrationNormal, vu, du);
                    wLocal.axpy(bndryQuad[pt], phiKern, dphiKern);
                  } else if constexpr (hasRobin) {
                    RangeType phiKern = weight *integrationElement * model.robinFlux(bndryQuad[pt], integrationNormal, vu, du);
                    wLocal.axpy(bndryQuad[pt], phiKern);
                  } else if constexpr (hasSingular) {
                    JacobianRangeType dphiKern = weight * integrationElement * model.singularFlux(bndryQuad[pt], integrationNormal, vu, du);
                    wLocal.axpy(bndryQuad[pt], dphiKern);
                  }
                } // quadrature loop
              } // Robin branch
            } // domain boundary
          } // intersection loop
        } // end of face integrals branch
 
        uLocal.unbind();
        model.unbind();
 
      } // end mesh loop
 
      // communicate data (in parallel runs)
      w.communicate();
 
      // Apply constraints
      constraints()(u, w);
    }
 
    // helper method in order to distinguish between conforming and
    // non-conforming intersections.
    template<class Model,
             class DomainFunction, class RangeFunction,
             class Constraints,
             template<class> class QuadratureTraits>
    template<bool conforming>
    void EllipticOperator<Model, DomainFunction, RangeFunction, Constraints, QuadratureTraits>
    ::computeInteriorSkeletonContributions(
      const GridPartType& gridPart,
      const IntersectionType& intersection,
      const RangeFieldType& penalty,
      const LocalDomainFunctionType& uLocal,
      const LocalDomainFunctionType& uLocalNb,
      const ModelClosure& model,
      const ModelClosure& modelNb,
      LocalFunctionType& wLocal) const
    {
      const int quadOrder = std::max(uLocal.order(),uLocalNb.order()) + wLocal.order();
      IntersectionQuadrature<conforming> intersectionQuad(gridPart, intersection, quadOrder);
      const auto& quadInside = intersectionQuad.inside();
      const auto& quadOutside = intersectionQuad.outside();
 
      const size_t numQuadraturePoints = quadInside.nop();
 
      for (size_t pt = 0; pt < numQuadraturePoints; ++pt) {
        // get coordinate of quadrature point on the reference element of the intersection
        const auto& x = quadInside.localPoint(pt);
        const auto integrationNormal = intersection.integrationOuterNormal(x);
        const auto integrationElement = integrationNormal.two_norm();
        const auto weight = quadInside.weight( pt );
 
        RangeType value;
        JacobianRangeType dvalue;
 
        RangeType vuIn, vuOut, jump;
        JacobianRangeType duIn, duOut;
        uLocal.evaluate(quadInside[pt], vuIn);
        uLocal.jacobian(quadInside[pt], duIn);
        auto aduIn = model.flux(quadInside[pt], vuIn, duIn);
 
        uLocalNb.evaluate(quadOutside[ pt ], vuOut);
        uLocalNb.jacobian(quadOutside[ pt ], duOut);
        auto aduOut = modelNb.flux(quadOutside[pt], vuOut, duOut);
 
 
        jump = vuIn - vuOut;
        // penalty term : beta [u] [phi] = beta (u+ - u-)(phi+ - phi-)=beta (u+ - u-)phi+
        //value = jump;
        value = penalty * integrationElement * jump;
        // {A grad u}.[phi] = {A grad u}.phi+ n_+ = 0.5*(grad u+ + grad u-).n_+ phi+
        //aduIn += aduOut;
        //aduIn *= -0.5;
        value += contractInner(-(aduIn + aduOut)/2_f, integrationNormal);
        //  [ u ] * { A grad phi_en } = -normal(u+ - u-) * 0.5 grad phi_en
        // here we need a dyadic product of u x n
#if 0
        dvalue = -0.5 * sFormDG * outer(jump, integrationNormal);
#else
        for (int r = 0; r < dimRange; ++r) {
          for (int d = 0; d < dimDomain; ++d) {
            dvalue[r][d] = -0.5 * sFormDG * integrationNormal[d] * jump[r];
          }
        }
#endif
 
        auto adValue = model.flux(quadInside[pt], jump, dvalue);
 
        wLocal.axpy(quadInside[pt], (RangeType)(weight * value), (JacobianRangeType)(weight * adValue));
      } // end quadrature loop
    }
 
    // implementation of Jacobian
    template<class JacobianOperator,
             class Model, class Constraints,
             template<class> class QuadratureTraits>
    void DifferentiableEllipticOperator<JacobianOperator, Model, Constraints, QuadratureTraits>
    ::jacobian (const DomainFunctionType &u, JacobianOperator &jOp) const
    {
      // If we use DG we need Diagonal and Neighbor stencil for entity-entity-pairing and entity-neighbor-pairing
      //Otherwise the diagonal Stencil suffices
      typedef typename std::conditional<!(useDG),
                                        Dune::Fem::DiagonalStencil<DiscreteFunctionSpaceType,
                                                                   DiscreteFunctionSpaceType>,
                                        Dune::Fem::DiagonalAndNeighborStencil<DiscreteFunctionSpaceType,
                                                                              DiscreteFunctionSpaceType>
          >::type StencilType;
      typedef Dune::Fem::TemporaryLocalMatrix<DiscreteFunctionSpaceType, DiscreteFunctionSpaceType> ElementMatrixType;
 
      const DiscreteFunctionSpaceType & dfSpace = u.space();
      const GridPartType& gridPart = dfSpace.gridPart();
 
      // Reserve Stencil storage
      jOp.reserve(StencilType(dfSpace, dfSpace));
      jOp.clear();
 
      // BIG FAT NOTE: we need maxNumDofs * blockSize many items
      const size_t maxNumBasisFunctions =
        DiscreteFunctionSpaceType::localBlockSize * dfSpace.blockMapper().maxNumDofs();
      // allocate space for phi_, dphi_, phiNb_, dphiNb_ ...
      reserveAssemblyStorage(maxNumBasisFunctions);
 
      ElementMatrixType elementMatrix(dfSpace, dfSpace);
      ElementMatrixType elementMatrixNb(dfSpace, dfSpace); // for DG
 
      // storage for point of linearization, will never be initialized
      // for affine-linear models
      LocalDomainFunctionType uLocal(u);
      RangeType u0(0);
      JacobianRangeType Du0(0);
 
      // For DG:
      LocalDomainFunctionType uLocalNb(u);
      RangeType u0Nb(0);
      JacobianRangeType Du0Nb(0);
 
      // mutable copy of the model-closure
      auto model = closure();
 
      // For DG we need a copy of the model which lives on the neighbour
      auto modelNb = closure();
 
      // For DG the dirichlet constraints are applied differently
      //via the boundary terms.
      constraints().update();
      //bool totallyConstrained = false;
 
      const auto end = dfSpace.end();
      for (auto it = dfSpace.begin(); it != end; ++it) {
        const auto &entity = *it;
 
        // totallyConstrained = constraints-get-information-on-entity(entity) or so
 
        //if (totallyConstrained) {
          // if all DOFs are fixed then there is nothing to compute ...
        //  continue;
        //}
 
        const auto& geometry = entity.geometry();
 
        model.bind(entity); // call per-element initializer for the operator
 
        if constexpr (!isAffineLinear) { // get values for linearization
          uLocal.bind(entity);
        }
 
        // use a temporary matrix
        elementMatrix.init(entity, entity);
        elementMatrix.clear();
 
        const auto& baseSet = elementMatrix.domainBasisFunctionSet();
        const auto numBaseFunctions = baseSet.size();
 
        if constexpr ((hasSources && !hasMassQuadrature) || hasFlux) {
          // use one qudrature for 0- and 2nd-order term
 
          QuadratureType quadrature(entity, 2*dfSpace.order());
          const size_t numQuadraturePoints = quadrature.nop();
 
          for (size_t pt = 0; pt < numQuadraturePoints; ++pt) {
 
            const auto& x = quadrature.point(pt);
            const auto weight = quadrature.weight(pt) * geometry.integrationElement(x);
 
            // evaluate all basis functions at given quadrature point
            baseSet.evaluateAll(quadrature[pt], phi_);
 
            // evaluate jacobians of all basis functions at given quadrature point
            baseSet.jacobianAll(quadrature[pt], dphi_);
 
            if constexpr (!isAffineLinear) { // get values for linearization
              uLocal.evaluate(quadrature[pt], u0);
              uLocal.jacobian(quadrature[pt], Du0);
            }
 
            for (size_t localCol = 0; localCol < numBaseFunctions; ++localCol) {
              if constexpr (hasSources && !hasMassQuadrature && hasFlux) {
                auto aphi = model.linearizedSource(u0, Du0, quadrature[pt], phi_[localCol], dphi_[localCol]);
                auto adphi = model.linearizedFlux(u0, Du0, quadrature[pt], phi_[localCol], dphi_[localCol]);
                elementMatrix.column(localCol).axpy(phi_, dphi_, aphi, adphi, weight);
              } else if constexpr (hasFlux) {
                auto adphi = model.linearizedFlux(u0, Du0, quadrature[pt], phi_[localCol], dphi_[localCol]);
                elementMatrix.column(localCol).axpy(dphi_, adphi, weight);
              } else {
                auto aphi = model.linearizedSource(u0, Du0, quadrature[pt], phi_[localCol], dphi_[localCol]);
                elementMatrix.column(localCol).axpy(phi_, aphi, weight);
              }
            }
          }
        }
 
        // Special case for e.g. mass-lumping using a "lumping"-quadrature
        if constexpr (hasSources && hasMassQuadrature) {
 
          MassQuadratureType quadrature(entity, 2*dfSpace.order());
          const size_t numQuadraturePoints = quadrature.nop();
 
          for (size_t pt = 0; pt < numQuadraturePoints; ++pt) {
 
            const typename MassQuadratureType::CoordinateType &x = quadrature.point(pt);
            const double weight = quadrature.weight(pt) * geometry.integrationElement(x);
 
            //const typename ElementGeometryType::Jacobian &gjit = geometry.jacobianInverseTransposed(x);
 
            // evaluate all basis functions at given quadrature point
            baseSet.evaluateAll(quadrature[pt], phi_);
 
            // evaluate jacobians of all basis functions at given quadrature point
            baseSet.jacobianAll(quadrature[pt], dphi_);
 
            if constexpr (!isAffineLinear) { // get value for linearization
              uLocal.evaluate(quadrature[pt], u0);
              uLocal.jacobian(quadrature[pt], Du0);
            }
 
            for (size_t localCol = 0; localCol < numBaseFunctions; ++localCol) {
 
              auto aphi = model.linearizedSource(u0, Du0, quadrature[pt], phi_[localCol], dphi_[localCol]);
 
              // get column object and call axpy method
              elementMatrix.column(localCol).axpy(phi_, aphi, weight);
            }
          }
        } // mass lumping part
 
        // end bulk integration
 
        if constexpr (useDG || hasRobin) { // skeleton contributions
 
          const auto area = useDG ? geometry.volume() : 0.0; // not needed for Robin
 
          const auto end = gridPart.iend(entity);
          for (auto it = gridPart.ibegin(entity); it != end; ++it) {
            const auto& intersection = *it;
 
            if (useDG && intersection.neighbor()) {
              const auto neighbor = intersection.outside();
              const auto& intersectionGeometry = intersection.geometry();
 
              // compute penalty factor
              const auto intersectionArea = intersectionGeometry.volume();
              //beta is still missing the diffusion coefficient!!
              const auto beta = penalty() * intersectionArea / std::min(area, neighbor.geometry().volume());
 
              // bind neighbour per-entity objects
              modelNb.bind(neighbor);
              if constexpr (!isAffineLinear) { // get values for linearization
                uLocalNb.bind(neighbor);
              }
 
              // use a temporary matrix for entity-neighbour contributions
              elementMatrixNb.init(neighbor, entity);
              elementMatrixNb.clear();
 
              if (intersection.conforming()) {
                computeInteriorSkeletonContributions<true>(
                  gridPart, intersection, beta, uLocal, uLocalNb, model, modelNb,
                  elementMatrix, elementMatrixNb);
              } else {
                computeInteriorSkeletonContributions<false>(
                  gridPart, intersection, beta, uLocal, uLocalNb, model, modelNb,
                  elementMatrix, elementMatrixNb);
              }
 
              // finished one DG-intersection, add element matrix to global operator.
              jOp.addLocalMatrix(neighbor, entity, elementMatrixNb);
 
              if constexpr (!isAffineLinear) {
                uLocalNb.unbind();
              }
              // unbind neighbour model
              modelNb.unbind();
 
            } else if (intersection.boundary() && (hasRobin || hasSingular)) {
 
              auto bdCond = model.classifyBoundary(intersection);
 
              if (bdCond.first && !bdCond.second.all()) {
 
                // TODO: make more realistic assumptions on the
                // quad-order. HOWTO?
                const auto quadOrder = 3*dfSpace.order();
                BndryMassQuadratureType intersectionQuad(gridPart, intersection, quadOrder);
                const auto& bndryQuad = intersectionQuad.inside();
                const auto numQuadraturePoints = bndryQuad.nop();
 
                for (size_t pt = 0; pt < numQuadraturePoints; ++pt) {
                  // One quad-loop for all contributions
 
                  // quadrature weight
                  const auto weight = bndryQuad.weight(pt);
                  (void)weight;
 
                  // Get the un-normalized outer normal
                  auto normal = intersection.integrationOuterNormal(bndryQuad.localPoint(pt));
                  const auto integrationElement = normal.two_norm();
                  normal /= integrationElement;
 
                  // Get values of all basis functions at the
                  // quad-points. Gradients are not needed.
                  baseSet.evaluateAll(bndryQuad[pt], phi_);
                  baseSet.jacobianAll(bndryQuad[pt], dphi_);
 
                  if constexpr (!isAffineLinear) {
                    // get value for linearization
                    uLocal.evaluate(bndryQuad[pt], u0);
                    uLocal.jacobian(bndryQuad[pt], Du0);
                  }
 
                  for (size_t localCol = 0; localCol < numBaseFunctions; ++localCol) {
                    if  constexpr (hasRobin && hasSingular) {
                      auto alphaphi = integrationElement * model.linearizedRobinFlux(u0, Du0, bndryQuad[pt], normal, phi_[localCol], dphi_[localCol]);
                      auto alphadphi = integrationElement * model.linearizedSingularFlux(u0, Du0, bndryQuad[pt], normal, phi_[localCol], dphi_[localCol]);
                      elementMatrix.column(localCol).axpy(phi_, dphi_, alphaphi, alphadphi, weight);
                    } else if  constexpr (hasRobin) {
                      auto alphaphi = integrationElement * model.linearizedRobinFlux(u0, Du0, bndryQuad[pt], normal, phi_[localCol], dphi_[localCol]);
                      elementMatrix.column(localCol).axpy(phi_, alphaphi, weight);
                    } else if constexpr (hasSingular) {
                      auto alphadphi = integrationElement * model.linearizedSingularFlux(u0, Du0, bndryQuad[pt], normal, phi_[localCol], dphi_[localCol]);
                      elementMatrix.column(localCol).axpy(dphi_, alphadphi, weight);
                    }
                  } // basis function loop
                } // quad loop
              } // Robin branch
            } // boundary contributions
          } // intersection loop
        } // end of skeleton part
 
        jOp.addLocalMatrix(entity, entity, elementMatrix);
 
        if constexpr (!isAffineLinear) {
          uLocal.unbind();
        }
        model.unbind(); // release the model
 
      } // end grid traversal
 
      releaseAssemblyStorage();
 
      constraints().jacobian(u, jOp);
 
      jOp.finalize();
 
      //jOp.print(std::cerr);
    }
 
    // helper method in order to distinguish between conforming and
    // non-conforming intersections.
    template<class JacobianOperator,
             class Model, class Constraints,
             template<class> class QuadratureTraits>
    template<bool conforming>
    void DifferentiableEllipticOperator<JacobianOperator, Model, Constraints, QuadratureTraits>
    ::computeInteriorSkeletonContributions(
      const GridPartType& gridPart,
      const IntersectionType& intersection,
      const RangeFieldType& penalty,
      const LocalDomainFunctionType& uLocal,
      const LocalDomainFunctionType& uLocalNb,
      const ModelClosure& model,
      const ModelClosure& modelNb,
      ElementMatrixType& elementMatrix,
      ElementMatrixType& elementMatrixNb) const
    {
      // FIXME: is this reasonable?
      const auto quadOrder = (!isAffineLinear ? std::max(uLocal.order(), uLocalNb.order()) : 0) + elementMatrix.rangeSpace().order() + 1;
      IntersectionQuadrature<conforming> intersectionQuad(gridPart, intersection, quadOrder);
      const auto& quadInside = intersectionQuad.inside();
      const auto& quadOutside = intersectionQuad.outside();
 
      // get neighbor's base function set
      const auto& baseSetNb = elementMatrixNb.domainBasisFunctionSet();
 
      // refetch to avoid yet another argument to the helper function
      const auto& baseSet = elementMatrix.domainBasisFunctionSet();
      const auto numBaseFunctions = baseSet.size();
 
      // storage for point of linearization
      RangeType u0, u0Nb;
      JacobianRangeType Du0, Du0Nb;
 
      const auto numFaceQuadPoints = quadInside.nop();
      for (size_t pt = 0; pt < numFaceQuadPoints; ++pt) {
        const auto& x = quadInside.localPoint(pt);
        auto normal = intersection.integrationOuterNormal(x);
        const auto faceVol = normal.two_norm();
        normal /= faceVol; // make it into a unit normal
 
        const auto quadWeight = quadInside.weight(pt);
        const auto weight = quadWeight * faceVol;
 
        if constexpr (!isAffineLinear) { // get values for linearization
          uLocal.evaluate(quadInside[pt], u0);
          uLocal.jacobian(quadInside[pt], Du0);
          uLocalNb.evaluate(quadInside[pt], u0Nb);
          uLocalNb.jacobian(quadInside[pt], Du0Nb);
        }
 
        // evaluate basis function of face inside E^- (entity)
 
        baseSet.evaluateAll(quadInside[pt], phi_);
        baseSet.jacobianAll(quadInside[pt], dphi_);
 
        // evaluate basis function of face inside E^+ (neighbor)
 
        baseSetNb.evaluateAll(quadOutside[pt], phiNb_);
        baseSetNb.jacobianAll(quadOutside[pt], dphiNb_);
 
        for (size_t i = 0; i < numBaseFunctions; ++i) {
          auto adphiEn = dphi_[i];
          auto adphiNb = dphiNb_[i];
          dphi_[i] = model.linearizedFlux(u0, Du0, quadInside[ pt ], phi_[i], adphiEn);
          dphiNb_[i] = modelNb.linearizedFlux(u0Nb, Du0Nb, quadOutside[ pt ], phiNb_[i], adphiNb);
        }
 
        // At this point dphi and dphiNb contain the value of
        // the "diffusive" flux applied to the test functions.
 
        // [Assemble skeleton terms: compute factors for axpy method]
        for (size_t localCol = 0; localCol < numBaseFunctions; ++localCol) {
          RangeType valueEn(0), valueNb(0); // NB: usmv() only adds!
          JacobianRangeType dvalueEn(0), dvalueNb(0); // NB: usmv() only adds!
 
          //  -{ A grad u } * [ phi_en ]
          dphi_[localCol].usmv(-0.5, normal, valueEn);
 
          //  -{ A grad u } * [ phi_en ]
          dphiNb_[localCol].usmv(-0.5, normal, valueNb);
 
          //  [ u ] * [ phi_en ] = u^- * phi_en^-
          valueEn.axpy(penalty, phi_[localCol]);
 
          valueNb.axpy(-penalty, phiNb_[localCol] );
          // here we need a dyadic product of u x n
#if 0
          dvalueEn = -0.5 * sFormDG * outer(phi_[localCol], normal);
          dvalueNb = 0.5 * sFormDG * outer(phiNb_[localCol] , normal);
#else
          for (int r = 0; r< dimRange; ++r) {
            for (int d = 0; d< dimDomain; ++d) {
              //  [ u ] * { grad phi_en }
              dvalueEn[r][d] = - 0.5 * sFormDG * phi_[localCol][r] * normal[d];
 
              //  [ u ] * { grad phi_en }
              dvalueNb[r][d] = 0.5 * sFormDG *  phiNb_[localCol][r] * normal[d];
            }
          }
#endif
 
          elementMatrix.column(localCol).axpy(phi_, dphi_, valueEn, dvalueEn, weight);
          elementMatrixNb.column(localCol).axpy(phi_, dphi_, valueNb, dvalueNb, weight);
 
        } // end loop over DoFs
 
      } // end face quadrature loop
 
    }
 
 
 
  } // namespace ACFem
 
} // namespace Dune
 
#endif // #ifndef __ELLIPTIC_OPERATOR_HH__