ellipticoperator.hh

#ifndef __ELLIPTIC_OPERATOR_HH__
#define __ELLIPTIC_OPERATOR_HH__
 
#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/operatorparts/operatorparts.hh"
#include "../common/quadrature.hh"
#include "../functions/constantfunction.hh"
 
namespace Dune {
 
  namespace ACFem {
 
    template<class OperatorParts,
             class DomainFunction,
             class RangeFunction = DomainFunction,
             class Constraints = EmptyBlockConstraints<RangeFunction>,
             template<class> class QuadratureTraits = DefaultQuadratureTraits>
    class EllipticOperator
      : public virtual std::conditional<OperatorParts::isLinear,
                                        Fem::LinearOperator<DomainFunction, RangeFunction>,
                                        Fem::Operator<DomainFunction, RangeFunction> >::type
    {
     public:
      typedef DomainFunction DomainFunctionType;
      typedef RangeFunction  RangeFunctionType;
      typedef RangeFunctionType DiscreteFunctionType;
      static_assert(std::is_base_of<Fem::IsDiscreteFunction, DiscreteFunctionType>::value,
                    "The RangeFunctionType has to be a discrete function");
      static_assert(std::is_base_of<Fem::HasLocalFunction, DiscreteFunctionType>::value,
                    "The DomainFunctionType has to be provide a local function");
      typedef OperatorParts OperatorPartsType;
      typedef Constraints ConstraintsOperatorType;
     protected:
      typedef typename std::conditional<
        std::is_same<ConstraintsOperatorType, EmptyBlockConstraints<RangeFunction> >::value,
        EmptyBlockConstraints<RangeFunction>,
        const ConstraintsOperatorType&>::type
      ConstraintsOperatorStorageType;
     protected:
      typedef typename ConstraintsOperatorType::LocalOperatorType LocalConstraintsOperatorType;
 
      typedef typename DomainFunctionType::FunctionSpaceType DomainFunctionSpaceType;
      typedef typename RangeFunctionType::FunctionSpaceType RangeFunctionSpaceType;
 
      static_assert(std::is_same<DomainFunctionSpaceType, RangeFunctionSpaceType>::value,
                    "Domain and range function spaces have to be the same");
 
      typedef typename DiscreteFunctionType::DiscreteFunctionSpaceType DiscreteFunctionSpaceType;
      typedef typename DiscreteFunctionType::LocalFunctionType LocalFunctionType;
      typedef typename LocalFunctionType::RangeType RangeType;
      typedef typename LocalFunctionType::JacobianRangeType JacobianRangeType;
 
      typedef typename DomainFunctionType::LocalFunctionType LocalDomainFunctionType;
 
      typedef typename DiscreteFunctionSpaceType::IteratorType IteratorType;
      typedef typename IteratorType::Entity       EntityType;
      typedef typename EntityType::Geometry       ElementGeometryType;
 
      typedef typename DiscreteFunctionSpaceType::DomainType DomainType;
 
      typedef typename DiscreteFunctionSpaceType::GridPartType  GridPartType;
      typedef typename GridPartType::IntersectionIteratorType IntersectionIteratorType;
      typedef typename IntersectionIteratorType::Intersection IntersectionType;
      typedef typename IntersectionType::Geometry  IntersectionGeometryType;
 
      // This need much more work. Autsch.
      typedef QuadratureTraits<GridPartType> QuadratureTraitsType;
      typedef typename QuadratureTraitsType::BulkQuadratureType QuadratureType;
      typedef typename QuadratureTraitsType::BulkMassQuadratureType MassQuadratureType;
      typedef typename QuadratureTraitsType::FaceQuadratureType FaceQuadratureType;
      typedef typename QuadratureTraitsType::FaceMassQuadratureType FaceMassQuadratureType;
 
      // use IntersectionQuadrature to create appropriate face quadratures
      typedef Dune::Fem::IntersectionQuadrature<FaceQuadratureType, true> IntersectionQuadratureType;
      // In principle we know that BndryQuadratureType equals FaceQuadratureType, but stick to the redirection.
      typedef typename IntersectionQuadratureType::FaceQuadratureType     BndryQuadratureType;
 
      typedef Dune::Fem::IntersectionQuadrature<FaceMassQuadratureType, true> IntersectionMassQuadratureType;
      // In principle we know that BndryQuadratureType equals FaceQuadratureType, but stick to the redirection.
      typedef typename IntersectionMassQuadratureType::FaceQuadratureType     BndryMassQuadratureType;
 
     public:
      EllipticOperator(const OperatorPartsType &parts,
                       const ConstraintsOperatorType& constraints = ConstraintsOperatorType())
        : operatorParts_(parts), constraintsOperator_(constraints)
      {}
 
      virtual void operator()(const DomainFunctionType &u, RangeFunctionType &w) const;
 
      virtual bool symmetric() const
      {
        return OperatorPartsType::isSymmetric;
      }
 
      virtual bool positiveDefinite() const
      {
        return OperatorPartsType::isSemiDefinite;
      }
 
     protected:
      const ConstraintsOperatorType& constraints() const { return constraintsOperator_; }
      const OperatorPartsType& operatorParts() const { return operatorParts_; }
      OperatorPartsType& operatorParts() { return operatorParts_; }
     private:
      OperatorPartsType operatorParts_;
      ConstraintsOperatorStorageType constraintsOperator_;
    };
 
    template<class JacobianOperator,
             class OperatorParts,
             class Constraints = DifferentiableEmptyBlockConstraints<JacobianOperator>,
             template<class> class QuadratureTraits = DefaultQuadratureTraits>
    class DifferentiableEllipticOperator
      : public EllipticOperator<OperatorParts,
                                typename JacobianOperator::DomainFunctionType,
                                typename JacobianOperator::RangeFunctionType,
                                Constraints, QuadratureTraits>,
        public Fem::DifferentiableOperator<JacobianOperator>,
        public virtual std::conditional<OperatorParts::isLinear,
                                        Fem::AssembledOperator<typename JacobianOperator::DomainFunctionType,
                                                               typename JacobianOperator::RangeFunctionType>,
                                        Fem::Operator<typename JacobianOperator::DomainFunctionType,
                                                      typename JacobianOperator::RangeFunctionType> >::type
    {
      typedef
      EllipticOperator<OperatorParts,
                       typename JacobianOperator::DomainFunctionType,
                       typename JacobianOperator::RangeFunctionType,
                       Constraints, QuadratureTraits>
      BaseType;
     public:
      typedef JacobianOperator JacobianOperatorType;
      typedef typename JacobianOperatorType::DomainFunctionType DomainFunctionType;
      typedef typename JacobianOperatorType::RangeFunctionType RangeFunctionType;
      typedef DomainFunctionType DiscreteFunctionType;
      typedef typename DomainFunctionType::LocalFunctionType LocalDomainFunctionType;
      typedef OperatorParts OperatorPartsType;
      typedef Constraints ConstraintsOperatorType;
     protected:
      typedef typename ConstraintsOperatorType::LocalOperatorType LocalConstraintsOperatorType;
     protected:
      typedef typename DiscreteFunctionType::DiscreteFunctionSpaceType DiscreteFunctionSpaceType;
      typedef typename DiscreteFunctionType::LocalFunctionType LocalFunctionType;
      typedef typename LocalFunctionType::RangeType RangeType;
      typedef typename LocalFunctionType::JacobianRangeType JacobianRangeType;
 
      typedef typename DiscreteFunctionSpaceType::IteratorType IteratorType;
      typedef typename IteratorType::Entity       EntityType;
      typedef typename EntityType::Geometry       ElementGeometryType;
 
      typedef typename DiscreteFunctionSpaceType::DomainType DomainType;
 
      typedef typename DiscreteFunctionSpaceType::GridPartType  GridPartType;
      typedef typename GridPartType::IntersectionIteratorType IntersectionIteratorType;
      typedef typename IntersectionIteratorType::Intersection IntersectionType;
      typedef typename IntersectionType::Geometry  IntersectionGeometryType;
 
      // This need much more work. Autsch.
      typedef QuadratureTraits<GridPartType> QuadratureTraitsType;
      typedef typename QuadratureTraitsType::BulkQuadratureType QuadratureType;
      typedef typename QuadratureTraitsType::BulkMassQuadratureType MassQuadratureType;
      typedef typename QuadratureTraitsType::FaceQuadratureType FaceQuadratureType;
      typedef typename QuadratureTraitsType::FaceMassQuadratureType FaceMassQuadratureType;
 
      // use IntersectionQuadrature to create appropriate face quadratures
      typedef Dune::Fem::IntersectionQuadrature<FaceQuadratureType, true> IntersectionQuadratureType;
      // In principle we know that BndryQuadratureType equals FaceQuadratureType, but stick to the redirection.
      typedef typename IntersectionQuadratureType::FaceQuadratureType     BndryQuadratureType;
 
      typedef Dune::Fem::IntersectionQuadrature<FaceMassQuadratureType, true> IntersectionMassQuadratureType;
      // In principle we know that BndryQuadratureType equals FaceQuadratureType, but stick to the redirection.
      typedef typename IntersectionMassQuadratureType::FaceQuadratureType     BndryMassQuadratureType;
 
     public:
      DifferentiableEllipticOperator(const OperatorPartsType &parts,
                                     const ConstraintsOperatorType& constraints = ConstraintsOperatorType())
        : BaseType(parts, constraints)
      {}
 
      using BaseType::operator();
      using BaseType::symmetric;
      using BaseType::positiveDefinite;
 
      void jacobian(const DiscreteFunctionType &u, JacobianOperatorType &jOp) const;
 
     protected:
      using BaseType::operatorParts;
      using BaseType::constraints;
    };
 
    // Implementation of Application operator
    template<class OperatorParts,
             class DomainFunction, class RangeFunction,
             class Constraints,
             template<class> class QuadratureTraits>
    void EllipticOperator<OperatorParts, DomainFunction, RangeFunction, Constraints, QuadratureTraits>
    ::operator()(const DomainFunction &u, RangeFunction &w) const
    {
      // clear destination
      w.clear();
 
      // get discrete function space
      const DiscreteFunctionSpaceType &dfSpace = w.space();
      const GridPartType& gridPart = dfSpace.gridPart();
 
      const bool hasFlux    = OperatorPartsType::hasFlux;
      const bool hasSources = OperatorPartsType::hasSources;
      const bool hasRobin   = OperatorPartsType::hasRobinFlux;
      const bool hasMassQuadrature = !std::is_same<QuadratureType, MassQuadratureType>::value;
 
      // possibly update
      constraints().rebuild();
      const bool doConstraints = constraints().size() > 0;
 
      // iterate over grid
      const IteratorType end = dfSpace.end();
      for (IteratorType it = dfSpace.begin(); it != end; ++it) {
 
        // get entity (here element)
        const EntityType &entity = *it;
        // get elements geometry
        const ElementGeometryType &geometry = entity.geometry();
 
        // get local representation of the discrete functions
        const LocalDomainFunctionType uLocal = u.localFunction(entity);
 
        LocalFunctionType wLocal = w.localFunction(entity);
 
        // call per-element initializer for the integral contributions
        operatorParts().setEntity(entity);
 
        // obtain quadrature order
        const int quadOrder = uLocal.order() + wLocal.order();
 
        if (hasSources || hasFlux) { // element integral, this computes all bulk contributions
          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 (hasSources && !hasMassQuadrature) {
              RangeType avu(0);
              operatorParts().source(entity, quadrature[pt], vu, du, avu);
              phiKern += avu;
            }
 
            if (hasFlux) {
              JacobianRangeType adu(0);
              // apply flux (i.e. 2nd and first order terms)
              operatorParts().flux(entity, quadrature[pt], vu, du, adu);
              dPhiKern += adu;
            }
 
            // add to local function
            phiKern *= weight;
            wLocal.axpy(quadrature[pt], phiKern);
 
            // add to local function
            dPhiKern *= weight;
            wLocal.axpy(quadrature[pt], dPhiKern);
 
          }  // end quad loop
 
          // Special code-path to allow e.g. for mass-lumping
          if (hasSources && hasMassQuadrature) {
            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;
              operatorParts().source(entity, quadrature[pt], vu, du, phiKern);
 
              // add to local function
              phiKern *= weight;
              wLocal.axpy(quadrature[pt], phiKern);
            }
          }
 
        } // end bulk block
 
        if (hasRobin) { // boundary integral
 
          IntersectionIteratorType iend = gridPart.iend(entity);
          for (IntersectionIteratorType it = gridPart.ibegin(entity); it != iend; ++it) {
 
            const IntersectionType &intersection = *it;
 
            assert(intersection.conforming() == true);
 
            if (intersection.neighbor() || !intersection.boundary()) {
              // Although one might want to augment the operator with
              // stabilization terms later ("edge oriented FEM"). Maybe.
              continue;
            }
 
            if (!operatorParts().setIntersection(intersection)) {
              continue;
            }
 
            if (std::is_same<FaceQuadratureType, FaceMassQuadratureType>::value) {
              int quadOrder = 3*dfSpace.order();
 
              IntersectionQuadratureType intersectionQuad(gridPart, intersection, quadOrder);
              const BndryQuadratureType &bndryQuad = intersectionQuad.inside();
              const int numQuadraturePoints = bndryQuad.nop();
 
              for (int pt = 0; pt < numQuadraturePoints; ++pt) {
 
                // One quad-loop for all contributions
                RangeType phiKern;
 
                // 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);
 
                RangeType alphau;
                operatorParts().robinFlux(intersection, bndryQuad[pt], integrationNormal, vu, phiKern);
                phiKern *= weight * integrationElement;
                wLocal.axpy(bndryQuad[pt], phiKern);
              }
 
            } else { // potential mass-lumping
 
              IntersectionMassQuadratureType intersectionQuad(gridPart, intersection, quadOrder);
              const BndryMassQuadratureType &bndryQuad = intersectionQuad.inside();
              const int numQuadraturePoints = bndryQuad.nop();
 
              for (int pt = 0; pt < numQuadraturePoints; ++pt) {
                // One quad-loop for all contributions
 
                // 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);
 
                RangeType phiKern;
                operatorParts().robinFlux(intersection, bndryQuad[pt], integrationNormal, vu, phiKern);
                phiKern *= weight *integrationElement;
                wLocal.axpy(bndryQuad[pt], phiKern);
 
              } // end boundary mass-quadrature loop
 
            }
 
          } // end Intersection loop
 
        } // end boundary integral block
 
        if (false && doConstraints) {
          // It is probably cheaper to apply constraints globally
          LocalConstraintsOperatorType localConstraints(constraints().localOperator(entity));
 
          localConstraints(uLocal, wLocal);
        }
 
      } // end mesh loop
 
      // communicate data (in parallel runs)
      w.communicate();
 
      // Apply constraints
      if (doConstraints) {
        constraints()(u, w);
      }
    }
 
    // implementation of Jacobian
    template<class JacobianOperator, class OperatorParts, class Constraints, template<class> class QuadratureTraits>
    void DifferentiableEllipticOperator<JacobianOperator, OperatorParts, Constraints, QuadratureTraits>
    ::jacobian (const DiscreteFunctionType &u, JacobianOperator &jOp) const
    {
      typedef typename JacobianOperator::LocalMatrixType LocalMatrixType;
      typedef typename DiscreteFunctionSpaceType::BasisFunctionSetType BasisFunctionSetType;
      // Diagonal stencil for entity-entity-pairing
      typedef Dune::Fem::DiagonalStencil<DiscreteFunctionSpaceType, DiscreteFunctionSpaceType> StencilType;
      typedef Dune::Fem::TemporaryLocalMatrix<DiscreteFunctionSpaceType, DiscreteFunctionSpaceType> ElementMatrixType;
 
      const DiscreteFunctionSpaceType &dfSpace = u.space();
      const GridPartType& gridPart = dfSpace.gridPart();
 
      // DiagonalStencil
 
      jOp.reserve(StencilType(dfSpace, dfSpace));
      jOp.clear();
 
      // BIG FAT NOTE: we need maxNumDofs * blockSize many items
      const unsigned maxNumBasisFunctions =
             DiscreteFunctionSpaceType::localBlockSize * dfSpace.blockMapper().maxNumDofs();
      std::vector<typename LocalFunctionType::RangeType> phi(maxNumBasisFunctions);
      std::vector<typename LocalFunctionType::JacobianRangeType> dphi(maxNumBasisFunctions);
 
      const bool hasFlux        = OperatorPartsType::hasFlux;
      const bool hasSources     = OperatorPartsType::hasSources;
      const bool hasRobin       = OperatorPartsType::hasRobinFlux;
      const bool isAffineLinear = OperatorPartsType::isLinear;
      const bool hasMassQuadrature = !std::is_same<QuadratureType, MassQuadratureType>::value;
 
      // possibly update
      constraints().rebuild();
      const bool doConstraints  = constraints().size() > 0;
 
      RangeType u0(0);
      JacobianRangeType Du0(0);
 
      LocalConstraintsOperatorType localConstraints(constraints().localOperator());
      bool totallyConstrained = false;
 
      ElementMatrixType elementMatrix(dfSpace, dfSpace);
 
      const IteratorType end = dfSpace.end();
      for (IteratorType it = dfSpace.begin(); it != end; ++it) {
        const EntityType &entity = *it;
        const ElementGeometryType &geometry = entity.geometry();
 
        // call per-element initializer for the operator
        operatorParts().setEntity(entity);
 
        const LocalDomainFunctionType uLocal = u.localFunction(entity);
        //auto elementMatrix = jOp.localMatrix(entity, entity);
        elementMatrix.init(entity, entity);
        elementMatrix.clear();
 
        const BasisFunctionSetType &baseSet = elementMatrix.domainBasisFunctionSet();
        const unsigned int numBaseFunctions = baseSet.size();
 
        if (doConstraints) {
          localConstraints.init(entity);
          totallyConstrained = localConstraints.totallyConstrained();
        }
 
        if (!totallyConstrained && (hasFlux || hasSources)) { // bulk integration
 
          QuadratureType quadrature(entity, 2*dfSpace.order());
          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);
 
            //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], /*gjit,*/ dphi);
 
            if (!isAffineLinear) { // get values for linearization
              uLocal.evaluate(quadrature[pt], u0);
              uLocal.jacobian(quadrature[pt], Du0);
            }
 
            RangeType aphi(0);          // Zero order
            JacobianRangeType adphi(0); // Second order
 
            for (unsigned int localCol = 0; localCol < numBaseFunctions; ++localCol) {
 
              if (hasSources && !hasMassQuadrature) {
                operatorParts().linearizedSource(u0, Du0, entity, quadrature[pt], phi[localCol], dphi[localCol], aphi);
              } else {
                aphi = 0;
              }
 
              if (hasFlux) {
                operatorParts().linearizedFlux(u0, Du0, entity, quadrature[pt], phi[localCol], dphi[localCol], adphi);
              } else {
                adphi = 0;
              }
 
              // get column object and call axpy method
              elementMatrix.column(localCol).axpy(phi, dphi, aphi, adphi, weight);
            }
 
          }
 
          // Special case for e.g. mass-lumping using a "lumping"-quadrature
          if (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);
 
              if (!isAffineLinear) { // get value for linearization
                uLocal.evaluate(quadrature[pt], u0);
                uLocal.jacobian(quadrature[pt], Du0);
              }
 
              RangeType aphi(0);          // Zero order
              JacobianRangeType adphi(0); // Second order
 
              for (unsigned int localCol = 0; localCol < numBaseFunctions; ++localCol) {
 
                operatorParts().linearizedSource(u0, Du0, entity, quadrature[pt], phi[localCol], dphi[localCol], aphi);
 
                adphi = 0;
 
                // get column object and call axpy method
                elementMatrix.column(localCol).axpy(phi, dphi, aphi, adphi, weight);
              }
            }
          } // mass lumping part
 
        } // end bulk integration
 
        if (!totallyConstrained && hasRobin) { // boundary contribution
 
          IntersectionIteratorType iend = dfSpace.gridPart().iend(entity);
          for (IntersectionIteratorType it = dfSpace.gridPart().ibegin(entity); it != iend; ++it) {
 
            const IntersectionType &intersection = *it;
 
            assert(intersection.conforming() == true);
 
            if (intersection.neighbor() || !intersection.boundary()) {
              // Although one might want to augment the operator with
              // stabilization terms later ("edge oriented FEM"). Maybe.
              continue;
            }
 
            if (!operatorParts().setIntersection(intersection)) {
              continue;
            }
 
 
            // TODO: make more realistic assumptions on the
            // quad-order. HOWTO?
            int quadOrder = 3*dfSpace.order();
 
            IntersectionMassQuadratureType intersectionQuad(gridPart, intersection, quadOrder);
            const BndryMassQuadratureType &bndryQuad = intersectionQuad.inside();
            const int numQuadraturePoints = bndryQuad.nop();
 
            for (int pt = 0; pt < numQuadraturePoints; ++pt) {
              // One quad-loop for all contributions
 
              // quadrature weight
              const double weight = bndryQuad.weight(pt);
 
              // Get the un-normalized outer normal
              DomainType integrationNormal
                = intersection.integrationOuterNormal(bndryQuad.localPoint(pt));
 
              // Only needed for Robin
              double integrationElement = integrationNormal.two_norm();
 
              integrationNormal /= integrationElement;
 
              // Get values of all basis functions at the
              // quad-points. Gradients are not needed.
              baseSet.evaluateAll(bndryQuad[pt], phi);
 
              if (!isAffineLinear) {
                // get value for linearization
                uLocal.evaluate(bndryQuad[pt], u0);
              }
 
              for (unsigned int localCol = 0; localCol < numBaseFunctions; ++localCol) {
 
                RangeType alphaphi;
                operatorParts().linearizedRobinFlux(u0, intersection, bndryQuad[pt], integrationNormal, phi[localCol], alphaphi);
                alphaphi *= integrationElement;
 
                // get column object and call axpy method
                elementMatrix.column(localCol).axpy(phi, alphaphi, weight);
              } // basis function loop
            } // quad loop
 
          } // intersection loop
 
        } // boundary contribution
 
        if (doConstraints) {
          // As we can access matrix entries only via the LocalMatrix
          // objects applying constraints globally would require yet
          // another mesh-loop. So we do it locally.
          localConstraints.jacobian(uLocal /* unused */, elementMatrix);
        }
 
        jOp.addLocalMatrix(entity, entity, elementMatrix);
 
      } // end grid traversal
 
      jOp.communicate();
    }
 
 
 
  } // namespace ACFem
 
} // namespace Dune
 
#endif // #ifndef ELLIPTIC_HH