dg.hh

#ifndef DUNE_MMDG_DG_HH
#define DUNE_MMDG_DG_HH

#include <cmath>

#include <dune/common/dynmatrix.hh>
#include <dune/common/dynvector.hh>

#include <dune/grid/io/file/vtk.hh>

#include <dune/geometry/quadraturerules.hh>

#include <dune/mmdg/nonconformingp1vtkfunction.hh>

template<class GridView, class Mapper, class Problem>
class DG
{
public:
  using Scalar = typename GridView::ctype;
  static constexpr int dim = GridView::dimension;
  using Matrix = Dune::DynamicMatrix<Scalar>; //NOTE: what is an appropriate sparse matrix type? -> BCRS
  using Vector = Dune::DynamicVector<Scalar>;
  using VTKFunction = Dune::VTKFunction<GridView>;
  using P1Function = Dune::NonconformingP1VTKFunction<GridView,
    Dune::DynamicVector<double>>;

  //constructor
  DG (const GridView& gridView, const Mapper& mapper,
    const Problem& problem) :
    gridView_(gridView), mapper_(mapper), problem_(problem),
    dof((1 + dim) * gridView.size(0))
    {
      A = Matrix(dof, dof, 0.0); //initialize stiffness matrix
      b = Vector(dof, 0.0); //initialize load vector
      d = Vector(dof, 0.0); //initialize solution vector
    }

  const void operator() (const Scalar K, const Scalar mu)
  {
    //assemble stiffness matrix A and load vector b
    assembleSLE(K, mu);

    //NOTE: what would be an appropiate solver here?
    A.solve(d, b);

    //write solution to a vtk file
    writeVTKoutput();
  }

  const GridView& gridView_;
  const Mapper& mapper_; //element mapper for gridView_
  const Problem& problem_; //the DG problem to be solved

  const int dof; //degrees of freedom

private:
  //assemble stiffness matrix A and load vector b
  void assembleSLE (const Scalar K, const Scalar mu)
  {
    //we use the basis
    // phi_elem,0 (x) = indicator(elem);
    // phi_elem,i (x) = x[i]*indicator(elem);
    //for elem in elements(gridView) and i = 1,...,dim
    for (const auto& elem : elements(gridView_))
    {
      const int elemIdx = mapper_.index(elem);
      const auto& geo = elem.geometry();
      const double elemVol = geo.volume();
      const auto& center = geo.center();

      //in the system of linear equations (SLE) Ad = b,
      //the index elemIdxSLE refers to the basis function phi_elem,0
      //and the indices elemIdxSLE + i + 1, i = 1,...,dim, refer to the
      //basis function phi_elem,i
      const int elemIdxSLE = (dim + 1)*elemIdx;

/*      //TODO: can be done outside of the loop?
      const Dune::QuadratureRule<double,dim>& rule =
        Dune::QuadratureRules<double,dim>::rule(geo.type(), problem_.quadratureOrder());
      Dune::FieldVector<double, dim+1> update(0.0);

      //NOTE: how are quadrature rules in Dune applied correctly?
      for (const auto& ip : rule)
      {
        const auto& qp = ip.position();
        //NOTE: is the volume of the element taken into account
        //automatically?
        const double weight = ip.weight();

        //quadrature for int_elem q*phi_elem,0 dV
        b[elemIdxSLE] += weight * problem_.q(qp);

        //quadrature for int_elem q*phi_elem,i dV
        for (int i = 0; i < dim; i++)
        {
          b[elemIdxSLE + i + 1] += weight * qp[i] * problem_.q(qp);
        }
      }
*/
      //NOTE: makeshift solution for source term q = -1
      b[elemIdxSLE] += -elemVol;
      for (int i = 0; i < dim; i++)
      {
        b[elemIdxSLE + i + 1] += -elemVol * center[i];
      }

      for (int i = 0; i < dim; i++)
      {
        //exact evaluation of
        // int_elem K*grad(phi_elem,i)*grad(phi_elem,i) dV
        A[elemIdxSLE + i + 1][elemIdxSLE + i + 1] += K * elemVol;
      }

      //iterate over all intersection with the boundary of elem
      for (const auto& intersct : intersections(gridView_, elem))
      {
        const auto& normal = intersct.centerUnitOuterNormal();
        const auto& intersctGeo = intersct.geometry();
        const double intersctVol = intersctGeo.volume();
        const auto& intersctCenter = intersctGeo.center();

        //TODO: quadrature rule cannot be used for dim = 1!
  //        const Dune::QuadratureRule<double,dim-1>& secondOrderRule =
  //          Dune::QuadratureRules<double,dim-1>::rule(
  //          intersct.type(), 2);

        //storage for multiply used integral values,
        //note that only the lower diagonal and diagonal entries of
        //quadraticIntregrals are used
        Dune::FieldVector<Scalar, dim> linearIntegrals(0.0);
        Dune::FieldMatrix<Scalar, dim, dim> quadraticIntregrals(0.0);
        //NOTE: is there a better type for a symmetric matrix?

        for (int i = 0; i < dim; i++)
        {
          //use midpoint rule for exact evaluation of
          // int_intersct x_i ds
          linearIntegrals[i] = intersctVol * intersctCenter[i];

          for (int j = 0; j <= i; j++)
          {
            const auto& leftCorner = intersctGeo.corner(0);
            const auto& rightCorner = intersctGeo.corner(1);

            quadraticIntregrals[i][j] = intersctVol / 3 *
              ( leftCorner[i] * leftCorner[j] + rightCorner[i] * rightCorner[j]
              + 0.5 * (leftCorner[i] * rightCorner[j] +
                leftCorner[j] * rightCorner[i]) );
  /*
            //use second order quadrature rule for exact evaluation of
            // int_intersct x_i*x_j ds
            //NOTE: use Simpson's rule instead manually?
            for (const auto& ip : secondOrderRule)
            {
              const auto& qp = ip.position();
              quadraticIntregrals[i][j] += //NOTE: volume necessary?
                ip.weight() * qp[i] * qp[j] * intersctVol;
            }
  */
            //NOTE: probably unnecessary
            quadraticIntregrals[j][i] = quadraticIntregrals[i][j];
          }
        }

        //exact evaluation of
        // int_intersct mu*jump(phi_elem,0)*jump(phi_elem,0) ds
        A[elemIdxSLE][elemIdxSLE] += mu * intersctVol;

        if (intersct.neighbor()) //intersct has neighboring element
        {
          //index of the neighboring element
          const int neighborIdx = mapper_.index(intersct.outside());
          const int neighborIdxSLE = (dim + 1)*neighborIdx;

          for (int i = 0; i < dim; i++)
          { //we use the relations
            // int_intersct mu*jump(phi_elem,0)*jump(phi_elem,i) ds
            // = mu * int_intersct x_i ds
            //and
            // int_intersct avg(K*grad(phi_elem,i))*jump(phi_elem,0) ds
            // = 0.5 * K * normal[i] * vol(intersct)
            A[elemIdxSLE + i + 1][elemIdxSLE] +=
              mu * linearIntegrals[i] - 0.5 * K * normal[i] * intersctVol;

            for (int j = 0; j <= i; j++)
            {
              //we use the relations
              // int_intersct mu*jump(phi_elem,i)*jump(phi_elem,j) ds
              // = mu * int_intersct x_i*x_j ds
              //and
              // int_intersct avg(K*grad(phi_elem,i))*jump(phi_elem,j) ds
              // = 0.5 * K * normal[i] * int_intersct x_j ds
              A[elemIdxSLE + i + 1][elemIdxSLE + j + 1]
                += mu * quadraticIntregrals[i][j]
                  - 0.5 * K * (normal[i] * linearIntegrals[j]
                    + normal[j] * linearIntegrals[i]);
            }
          }

          if (neighborIdx > elemIdx)
          { //make sure that each facet is considered only once
            continue;
          }

          //exact evaluation of
          // int_intersct mu*jump(phi_elem,0)*jump(phi_neighbor,0) ds
          A[elemIdxSLE][neighborIdxSLE] += -mu * intersctVol;

          //stiffness matrix A is symmetric
          A[neighborIdxSLE][elemIdxSLE] += A[elemIdxSLE][neighborIdxSLE];

          for (int i = 0; i < dim; i++)
          {
            //we use the relations
            // int_intersct mu*jump(phi_elem,i)
            //  *jump(phi_neighbor,0) ds
            // = -mu*int_intersct x_i ds
            //and
            // int_intersct avg(K*grad(phi_neighbor,i))
            //  *jump(phi_elem,0) ds
            // = 0.5 * K * normal[i] * vol(intersct)
            A[elemIdxSLE + i + 1][neighborIdxSLE] +=
              -mu * linearIntegrals[i] + 0.5 * K * normal[i] * intersctVol;

            //we use the relations
            // int_intersct mu*jump(phi_elem,0)
            //  *jump(phi_neighbor,i) ds
            // = -mu*int_intersct x_i ds
            //and
            // int_intersct avg(K*grad(phi_neighbor,i))
            //  *jump(phi_elem,0) ds
            // = 0.5 * K * normal[i] * vol(intersct)
            A[elemIdxSLE][neighborIdxSLE + i + 1] +=
              -mu * linearIntegrals[i] - 0.5 * K * normal[i] * intersctVol;

            //stiffness matrix A is symmetric
            A[neighborIdxSLE][elemIdxSLE + i + 1] +=
              A[elemIdxSLE + i + 1][neighborIdxSLE];
            A[neighborIdxSLE + i + 1][elemIdxSLE] +=
              A[elemIdxSLE][neighborIdxSLE + i + 1];

            for (int j = 0; j <= i; j++)
            {
              //we use the relations
              // int_intersct mu*jump(phi_elem,i)
              //  *jump(phi_neighbor,j) ds
              // = -mu*int_intersct x_i*x_j ds
              //and
              // int_intersct avg(K*grad(phi_neighbor,j))
              //  *jump(phi_elem,i) ds
              // = 0.5 * K * normal[j] * int_intersct x_i ds
              //as well as
              // int_intersct avg(K*grad(phi_elem,i))
              //  *jump(phi_neighbor,j) ds
              // = -0.5 * K * normal[i] * int_intersct x_j ds
              A[elemIdxSLE + i + 1][neighborIdxSLE + j + 1] +=
                -mu * quadraticIntregrals[i][j]
                - 0.5 * K * (normal[j] * linearIntegrals[i]
                  - normal[i] * linearIntegrals[j]);

              //stiffness matrix A is symmetric
              A[neighborIdxSLE + j + 1][elemIdxSLE + i + 1] +=
                A[elemIdxSLE + i + 1][neighborIdxSLE + j + 1];

              if (i != j)
              {
                // int_intersct mu*jump(phi_elem,j)
                //  *jump(phi_neighbor,i) ds
                // = -mu*int_intersct x_i*x_j ds
                //and
                // int_intersct avg(K*grad(phi_neighbor,i))
                //  *jump(phi_elem,j) ds
                // = 0.5 * K * normal[i] * int_intersct x_j ds
                //as well as
                // int_intersct avg(K*grad(phi_elem,j))
                //  *jump(phi_neighbor,i) ds
                // = -0.5 * K * normal[j] * int_intersct x_i ds
                A[elemIdxSLE + j + 1][neighborIdxSLE + i + 1] +=
                  -mu * quadraticIntregrals[i][j]
                  - 0.5 * K * (normal[i] * linearIntegrals[j]
                    - normal[j] * linearIntegrals[i]);

                //stiffness matrix A is symmetric
                A[neighborIdxSLE + i + 1][elemIdxSLE + j + 1] +=
                  A[elemIdxSLE + j + 1][neighborIdxSLE + i + 1];
              }
            }
          }
        }
        else //boundary facet
        {
          for (int i = 0; i < dim; i++)
          { //we use the relations
            // int_intersct mu*jump(phi_elem,0)*jump(phi_elem,i) ds
            // = mu * int_intersct x_i ds
            //and for boundary facets
            // int_intersct avg(K*grad(phi_elem,i))
            //  *jump(phi_elem,0) ds
            // = K * normal[i] * vol(intersct)
            A[elemIdxSLE + i + 1][elemIdxSLE] +=
              mu * linearIntegrals[i] - 0.5 * K * normal[i] * intersctVol;

            for (int j = 0; j <= i; j++)
            {
              //we use the relations
              // int_intersct mu*jump(phi_elem,i)*jump(phi_elem,j) ds
              // = mu * int_intersct x_i*x_j ds
              //and for boundary facets
              // int_intersct avg(K*grad(phi_elem,i))
              //  *jump(phi_elem,j) ds
              // = 0.5 * K * normal[i] * int_intersct x_j ds
              A[elemIdxSLE + i + 1][elemIdxSLE + j + 1] +=
                mu * quadraticIntregrals[i][j]
                - 0.5 * K * (normal[i] * linearIntegrals[j]
                  + normal[j] * linearIntegrals[i]);
            }
          }
        }
      }

      //stiffness matrix A is symmetric
      for (int i = 0; i < dim; i++)
      {
        A[elemIdxSLE][elemIdxSLE + i + 1] = A[elemIdxSLE + i + 1][elemIdxSLE];

        for(int j = 0; j < i; j++)
        {
          A[elemIdxSLE + j + 1][elemIdxSLE + i + 1] =
            A[elemIdxSLE + i + 1][elemIdxSLE + j + 1];
        }
      }
    }

    //NOTE: check if A is symmetric
    for (int i = 0; i < dof; i++)
      for (int j = 0; j < i; j++)
        /*assert*/if(std::abs(A[i][j] - A[j][i]) >
          std::numeric_limits<Scalar>::epsilon())
            std::cout << i << ", " << j << std::endl;
  }

  //writes the solution to a vtk file
  void writeVTKoutput () const
  {
    //storage for pressure data, for each element we store the
    //pressure at the corners of the element, the VTKFunction will
    //create the output using a linear interpolation for each element
    Dune::DynamicVector<Scalar> pressure(dof, 0.0);
    Dune::DynamicVector<Scalar> exactPressure(dof, 0.0);

    for (const auto& elem : elements(gridView_))
    {
      const int elemIdxSLE = (dim + 1)*mapper_.index(elem);
      const auto& geo = elem.geometry();

      for (int k = 0; k < geo.corners(); k++)
      {
        if (problem_.hasExactSolution())
        {
          exactPressure[elemIdxSLE + k] = problem_.exactSolution(geo.corner(k));
        }

        //contribution of the basis function
        // phi_elem,0 (x) = indicator(elem);
        //at the kth corner of elem
        pressure[elemIdxSLE + k] = d[elemIdxSLE];

        for (int i = 0; i < dim; i++)
        {
          //contribution of the basis function
          // phi_elem,i (x) = x[i]*indicator(elem);
          //at the kth corner of elem
          pressure[elemIdxSLE + k] += d[elemIdxSLE + i + 1] * geo.corner(k)[i];
        }
      }
    }

    //create output directory if necessary
    mkdir("data", S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);

    Dune::VTKWriter<GridView> vtkWriter(gridView_, Dune::VTK::nonconforming);
    vtkWriter.addVertexData( std::shared_ptr<const VTKFunction>(
      new P1Function(gridView_, pressure, "pressure")));

    if (problem_.hasExactSolution())
    {
      vtkWriter.addVertexData( std::shared_ptr<const VTKFunction>(
        new P1Function(gridView_, exactPressure, "exactPressure")));
    }

    vtkWriter.pwrite("pressureData", "", "data");
  }


  Matrix A; //stiffness matrix
  Vector b; //load vector
  Vector d; //solution vector
};

#endif