interiorpenaltydgassembler.hh

Go to the documentation of this file.

// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
 
// SPDX-FileCopyrightText: Copyright © DUNE-FUFEM Project contributors, see file AUTHORS.md
// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception OR LGPL-3.0-or-later
 
#ifndef IPDG_ASSEMBLER_HH
#define IPDG_ASSEMBLER_HH
 
#ifndef DUNE_FUFEM_DISABLE_HEADER_DEPRECATION
#warning This header is deprecated and will be removed after 2.11. Use Dune::Fufem::Forms instead.
#endif
 
// dune-common includes
#include <dune/common/exceptions.hh>
#include <dune/common/fvector.hh>
#include <dune/common/fmatrix.hh>
#include <dune/common/typetraits.hh>
 
// dune-istl includes
#include <dune/istl/matrix.hh>
#include <dune/istl/scaledidmatrix.hh>
 
// dune-fufem includes
#include "dune/fufem/quadraturerules/quadraturerulecache.hh"
 
namespace Dune::Fufem::Impl {
 
  template <class Matrix>
  requires (not Dune::IsNumber<Matrix>::value)
  static void addToDiagonal(Matrix& x, const typename Matrix::field_type a) {
    for (int i = 0; i < Matrix::rows; ++i)
      x[i][i] += a;
  }
 
  // Specialization for scalars (to be interpreted as 1x1 matrices)
  template <class Matrix>
  requires Dune::IsNumber<Matrix>::value
  static void addToDiagonal(Matrix& x, const Matrix& a)
  {
    x += a;
  }
 
  template <typename FieldType, int n>
  static void addToDiagonal(
      Dune::ScaledIdentityMatrix<FieldType, n>& x,
      const typename Dune::ScaledIdentityMatrix<FieldType, n>::field_type a) {
    x.scalar() += a;
  }
 
} // end namespace Dune::Fufem::Impl
 
 
template <class GridType>
class
[[deprecated("This class is deprecated and will be removed after 2.11. Use Dune::Fufem::Forms instead.")]]
InteriorPenaltyDGAssembler
{
    enum class DGType {SIPG = -1, IIPG = 0, NIPG = 1};
 
    static const int dim      = GridType::dimension -1; // We're on the edges
    static const int dimworld = GridType::dimensionworld;
 
  public:
 
    [[deprecated("This class is deprecated and will be removed after 2.11. Use Dune::Fufem::Forms instead.")]]
    InteriorPenaltyDGAssembler() :
      sigma0_(0),
      dirichlet_(true),
      dgType_((double) DGType::SIPG),
      sigma1_(0) {}
 
    [[deprecated("This class is deprecated and will be removed after 2.11. Use Dune::Fufem::Forms instead.")]]
    InteriorPenaltyDGAssembler(double penalty, bool dirichlet=true, DGType dgType=DGType::SIPG, double gradientPenalty=0) :
      sigma0_(penalty),
      dirichlet_(dirichlet),
      dgType_((double) dgType),
      sigma1_(gradientPenalty) {}
 
    void setPenalty(double penaltyParameter) {
      if(penaltyParameter<0.0)
        DUNE_THROW(Dune::MathError, "Penalty parameter must not be negative");
      sigma0_ = penaltyParameter;
    }
 
    void setGradientPenalty(double gradientPenalty) {
      if(gradientPenalty<0.0)
        DUNE_THROW(Dune::MathError, "Penalty parameter must not be negative");
      sigma1_ = gradientPenalty;
    }
 
    template <class Edge, class BoolMatrix, class TestLocalFE, class AnsatzLocalFE>
    void indices(const Edge& it, BoolMatrix& isNonZero, const TestLocalFE& tFE, const AnsatzLocalFE& aFE) const
    {
      isNonZero = dirichlet_; // For boundary edges, only in case of Dirichlet values we need to assemble a local matrix.
    }
 
    template <class Edge, class BoolMatrix, class TestLocalFE, class AnsatzLocalFE>
    void indices(const Edge& it, BoolMatrix& isNonZero, const TestLocalFE& tFEi, const AnsatzLocalFE& aFEi, const TestLocalFE& tFEo, const AnsatzLocalFE& aFEo) const
    {
      isNonZero = true;
    }
 
 
    template <class Edge, class LocalMatrix, class TestLocalFE, class AnsatzLocalFE>
    void assemble(const Edge& edge, LocalMatrix& localMatrix, const TestLocalFE& tFE, const AnsatzLocalFE& aFE) const
    {
      localMatrix = 0.0;
 
      // For Neumann data, we don't assemble on boundary edges
      // TODO: There is no need that this has to be decided globally. It would be very easy to
      // just decide this edge by edge for mixed b.c..
      if (!dirichlet_)
        return;
 
      using FVdimworld = typename Dune::template FieldVector<double,dimworld>;
 
      using JacobianType = typename TestLocalFE::Traits::LocalBasisType::Traits::JacobianType;
      using RangeType = typename TestLocalFE::Traits::LocalBasisType::Traits::RangeType;
 
      auto rows = localMatrix.N();
      auto cols = localMatrix.M();
 
      // get geometries of the edge and the inside element
      const auto edgeGeometry = edge.geometry();
      const auto insideGeometry = edge.inside().geometry();
 
      // Get the length of the edge
      const auto edgeLength= edgeGeometry.volume();
 
      // get quadrature rule
      QuadratureRuleKey tFEquad(edge.type(), tFE.localBasis().order());
      QuadratureRuleKey quadKey = tFEquad.square();
 
      const auto& quad = QuadratureRuleCache<double, dim>::rule(quadKey);
 
      // store gradients of shape functions and base functions
      std::vector<JacobianType> referenceGradients(tFE.localBasis().size());
      std::vector<FVdimworld> gradients(tFE.localBasis().size());
 
      // store values of shape functions
      std::vector<RangeType> tFEvalues(tFE.localBasis().size());
 
      // We also need a unit normal on the edge pointing away from the inner element
      const auto outerNormal = edge.centerUnitOuterNormal();
 
      // loop over quadrature points
      for (size_t pt=0; pt < quad.size(); ++pt)
      {
        // get quadrature point
        const auto& quadPos = quad[pt].position();
 
        // get transposed inverse of Jacobian of transformation
        const auto& invJacobian = insideGeometry.jacobianInverseTransposed(edge.geometryInInside().global(quadPos));
 
        // get integration factor
        const auto integrationElement = edgeGeometry.integrationElement(quadPos);
 
        // get gradients of shape functions
        tFE.localBasis().evaluateJacobian(edge.geometryInInside().global(quadPos), referenceGradients);
 
        // transform gradients
        for (size_t i=0; i<gradients.size(); ++i)
          invJacobian.mv(referenceGradients[i][0], gradients[i]);
 
        // evaluate basis functions
        tFE.localBasis().evaluateFunction(edge.geometryInInside().global(quadPos), tFEvalues);
 
        // compute matrix entries
        auto z = quad[pt].weight() * integrationElement;
        for (size_t i=0; i<rows; ++i)
        {
          for (size_t j=0; j<cols; ++j)
          {
            auto zij = -z*tFEvalues[i][0]*(gradients[j]*outerNormal);
            zij += dgType_*z*tFEvalues[j]*(gradients[i]*outerNormal);
            zij += sigma0_*z/edgeLength*tFEvalues[i]*tFEvalues[j];
            // The sigma1_ stabilization term does not effect boundary edges
            Dune::Fufem::Impl::addToDiagonal(localMatrix[i][j],zij);
          }
        }
      }
    }
 
    template <class Edge, class LocalMatrix, class TestLocalFE, class AnsatzLocalFE>
    void assemble(const Edge& edge, LocalMatrix& localMatrix, const TestLocalFE& tFEinside, const AnsatzLocalFE& aFEinside, const TestLocalFE& tFEoutside, const AnsatzLocalFE& aFEoutside) const
    {
 
      localMatrix = 0.0;
 
      // to avoid too much duplicate code, we write this method, which was designed for the vintage intersection operator assembler, by calling
      // the newer assembleBlockwise method.
      using MatrixContainer = Dune::Matrix<LocalMatrix>;
      auto mc = MatrixContainer(2,2);
      mc[0][0].setSize(tFEinside.size(), aFEinside.size());
      mc[0][1].setSize(tFEinside.size(), aFEoutside.size());
      mc[1][0].setSize(tFEoutside.size(), aFEinside.size());
      mc[1][1].setSize(tFEoutside.size(), aFEoutside.size());
 
      // compute entries
      assembleBlockwise(edge, mc, tFEinside, aFEinside, tFEoutside, aFEoutside);
 
      // redistribute them into the localMatrix
      auto rows = localMatrix.N();
      auto cols = localMatrix.M();
 
      const auto insideSize = tFEinside.localBasis().size();
 
      for (size_t i=0; i<rows; ++i)
      {
        for (size_t j=0; j<cols; ++j)
        {
          // Depending on which elements i and j belong to, perform different operations.
          // In particular, i (resp. j) < insizeSize corresponds to functions from the inner element
          // while the others correspond to the outer element
          if (i < insideSize)
            localMatrix[i][j] += (j< insideSize) ?  mc[0][0][i][j] : mc[0][1][i][j-insideSize];
          else
            localMatrix[i][j] += (j< insideSize) ? mc[1][0][i-insideSize][j] : mc[1][1][i-insideSize][j-insideSize];
        }
      }
    }
 
    template <class Edge, class MC, class TestLocalFE, class AnsatzLocalFE>
    void assembleBlockwise(const Edge& edge, MC& matrixContainer, const TestLocalFE& tFEinside, const AnsatzLocalFE& aFEinside, const TestLocalFE& tFEoutside, const AnsatzLocalFE& aFEoutside) const
    {
 
      using FVdimworld = typename Dune::template FieldVector<double,dimworld>;
 
      using JacobianType = typename TestLocalFE::Traits::LocalBasisType::Traits::JacobianType;
      using RangeType = typename TestLocalFE::Traits::LocalBasisType::Traits::RangeType;
 
      // get geometry and store it
      const auto edgeGeometry = edge.geometry();
      const auto insideGeometry = edge.inside().geometry();
      const auto outsideGeometry = edge.outside().geometry();
 
      const auto edgeLength = edgeGeometry.volume();
      matrixContainer = 0.0;
 
      // get quadrature rule
      QuadratureRuleKey tFEquad(edge.type(), std::max(tFEinside.localBasis().order(),tFEoutside.localBasis().order()));
      QuadratureRuleKey quadKey = tFEquad.square();
 
      const auto& quad = QuadratureRuleCache<double, dim>::rule(quadKey);
 
      // store gradients for both the inner and outer elements
      std::vector<JacobianType> insideReferenceGradients(tFEinside.localBasis().size());
      std::vector<JacobianType> outsideReferenceGradients(tFEoutside.localBasis().size());
      std::vector<FVdimworld> insideGradients(tFEinside.localBasis().size());
      std::vector<FVdimworld> outsideGradients(tFEoutside.localBasis().size());
 
      // store values of shape functions
      std::vector<RangeType> tFEinsideValues(tFEinside.localBasis().size());
      std::vector<RangeType> tFEoutsideValues(tFEoutside.localBasis().size());
 
      const auto outerNormal = edge.centerUnitOuterNormal();
 
      // loop over quadrature points
      for (size_t pt=0; pt < quad.size(); ++pt)
      {
        // get quadrature point
        const auto& quadPos = quad[pt].position();
 
        // get transposed inverse of Jacobian of transformation
        const auto& invJacobian = insideGeometry.jacobianInverseTransposed(edge.geometryInInside().global(quadPos));
        const auto& outsideInvJacobian = outsideGeometry.jacobianInverseTransposed(edge.geometryInOutside().global(quadPos));
 
        // get integration factor
        const auto integrationElement = edgeGeometry.integrationElement(quadPos);
 
        // get gradients of shape functions on both the inside and outside element
        tFEinside.localBasis().evaluateJacobian(edge.geometryInInside().global(quadPos), insideReferenceGradients);
        tFEoutside.localBasis().evaluateJacobian(edge.geometryInOutside().global(quadPos), outsideReferenceGradients);
 
        // transform gradients
        for (size_t i=0; i<insideGradients.size(); ++i)
          invJacobian.mv(insideReferenceGradients[i][0], insideGradients[i]);
        for (size_t i=0; i<outsideGradients.size(); ++i)
          outsideInvJacobian.mv(outsideReferenceGradients[i][0], outsideGradients[i]);
 
        // evaluate basis functions
        tFEinside.localBasis().evaluateFunction(edge.geometryInInside().global(quadPos), tFEinsideValues);
        tFEoutside.localBasis().evaluateFunction(edge.geometryInOutside().global(quadPos), tFEoutsideValues);
 
        // compute matrix entries
        auto z = quad[pt].weight() * integrationElement;
 
        // Basis functions from inside as test functions
        for (size_t i=0; i<tFEinside.localBasis().size(); ++i)
        {
          // Basis functions from inside as ansatz functions
          for (size_t j=0; j<aFEinside.localBasis().size(); ++j)
          {
            // M11, see Riviere, p. 54f
            auto zij = -0.5*z*tFEinsideValues[i][0]*(insideGradients[j]*outerNormal);
            zij += 0.5*dgType_*z*tFEinsideValues[j]*(insideGradients[i]*outerNormal);
            zij += sigma0_*z/edgeLength*tFEinsideValues[i]*tFEinsideValues[j];
            zij += sigma1_*z/edgeLength*(insideGradients[i]*outerNormal)*(insideGradients[j]*outerNormal);
            Dune::Fufem::Impl::addToDiagonal(matrixContainer[0][0][i][j],zij);
          }
          // Basis functions from outside as ansatz functions
          for (size_t j=0; j<aFEoutside.localBasis().size(); ++j) {
            // M12
            auto zij = -0.5*z*tFEinsideValues[i][0]*(outsideGradients[j]*outerNormal);
            zij -= 0.5*dgType_*z*tFEoutsideValues[j]*(insideGradients[i]*outerNormal);
            zij -= sigma0_*z/edgeLength*tFEinsideValues[i]*tFEoutsideValues[j];
            zij -= sigma1_*z/edgeLength*(insideGradients[i]*outerNormal)*(outsideGradients[j]*outerNormal);
            Dune::Fufem::Impl::addToDiagonal(matrixContainer[0][1][i][j],zij);
          }
        }
        // Basis functions from outside as test functions
        for (size_t i=0; i<tFEoutside.localBasis().size(); ++i) {
          // Basis functions from inside as ansatz functions
          for (size_t j=0; j<aFEinside.localBasis().size(); ++j) {
            // M21, see Riviere, p. 54f
            auto zij = 0.5*z*tFEoutsideValues[i][0]*(insideGradients[j]*outerNormal);
            zij += 0.5*dgType_*z*tFEinsideValues[j]*(outsideGradients[i]*outerNormal);
            zij -= sigma0_*z/edgeLength*tFEoutsideValues[i]*tFEinsideValues[j];
            zij -= sigma1_*z/edgeLength*(outsideGradients[i]*outerNormal)*(insideGradients[j]*outerNormal);
            Dune::Fufem::Impl::addToDiagonal(matrixContainer[1][0][i][j],zij);
          }
          // Basis functions from outside as ansatz functions
          for (size_t j=0; j<aFEoutside.localBasis().size(); ++j) {
            // M22
            auto zij = 0.5*z*tFEoutsideValues[i][0]*(outsideGradients[j]*outerNormal);
            zij -= 0.5*dgType_*z*tFEoutsideValues[j]*(outsideGradients[i]*outerNormal);
            zij += sigma0_*z/edgeLength*tFEoutsideValues[i]*tFEoutsideValues[j];
            zij += sigma1_*z/edgeLength*(outsideGradients[i]*outerNormal)*(outsideGradients[j]*outerNormal);
            Dune::Fufem::Impl::addToDiagonal(matrixContainer[1][1][i][j],zij);
          }
        }
      }
    }
 
  private:
 
 
    double sigma0_; // Penalizes discontinuities
 
 
    bool dirichlet_;
 
    const double dgType_; // -1 leads to the symmetric interior penalty DG variant. {-1,0,1} are possible values.
 
    double sigma1_; // Penalizes jumps of the normal derivative along the non-boundary edges
};
#endif