add dg.hh and dgproblem.hh

dune/mmdg/dg.hh

+#ifndef DUNE_MMDG_DG_HH
+#define DUNE_MMDG_DG_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;
+
+  //constructor
+  DG (const GridView& gridView, const Mapper& mapper,
+    const Problem& problem) :
+    gridView_(gridView), mapper_(mapper), problem_(problem) {}
+
+  const void operator() (const Scalar K, const Scalar mu) const
+  {
+    //NOTE: what is an appropriate sparse matrix type?
+    const int dof = (1 + dim)*gridView_.size(0); //degrees of freedom
+    using Matrix = Dune::DynamicMatrix<Scalar>;
+    using Vector = Dune::DynamicVector<Scalar>;
+    Matrix A(dof, dof, 0.0); //stiffness matrix
+    Vector b(dof, 0.0); //load vector
+    Vector d;
+
+    //NOTE:
+//    const int order = 3; //order of the quadrature rule
+
+    //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();
+
+  /*    //TODO: can be done outside of the loop?
+      const Dune::QuadratureRule<double,dim>& rule =
+        Dune::QuadratureRules<double,dim>::rule(geo.type(),order);
+
+      //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 auto weight = ip.weight();
+
+        //quadrature for int_elem q*phi_elem,0 dV
+        b[(1 + dim)*elemIdx] += weight * q(qp);
+
+        //quadrature for int_elem q*phi_elem,i dV
+        for (int i = 0; i < dim; i++)
+        {
+          b[(1 + dim)*elemIdx + i + 1] += weight * qp[i] * q(qp);
+        }
+      }
+  */
+
+      //NOTE: makeshift solution for source term q = -1
+      b[(1 + dim)*elemIdx] = -elemVol;
+      for (int i = 0; i < dim; i++)
+      {
+        b[(1 + dim)*elemIdx + i + 1] = -elemVol * geo.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[(1 + dim)*elemIdx + i + 1][(1 + dim)*elemIdx + 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 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;
+            }
+*/
+          }
+        }
+
+        //exact evaluation of
+        // int_intersct mu*jump(phi_elem,0)*jump(phi_elem,0) ds
+        A[(1 + dim)*elemIdx][(1 + dim)*elemIdx] +=
+          mu * intersctGeo.volume();
+
+        if (intersct.neighbor()) //intersct has neighboring element
+        {
+          //index of the neighboring element
+          const int neighborIdx = mapper_.index(intersct.outside());
+
+          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[(1 + dim)*elemIdx + i + 1][(1 + dim)*elemIdx] +=
+              mu * linearIntegrals[i]
+              - 0.5 * K * normal[i] * intersctVol;
+
+            //stiffness matrix A is symmetric
+            A[(1 + dim)*elemIdx][(1 + dim)*elemIdx + i + 1] +=
+              A[(1 + dim)*elemIdx + i + 1][(1 + dim)*elemIdx];
+
+            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[(1 + dim)*elemIdx + i + 1][(1 + dim)*elemIdx + j + 1]
+                += mu * quadraticIntregrals[i][j]
+                  - 0.5 * K * (normal[i] * linearIntegrals[j]
+                    + normal[j] * linearIntegrals[i]);
+
+              //stiffness matrix A is symmetric
+              A[(1 + dim)*elemIdx + j + 1][(1 + dim)*elemIdx + i + 1]
+                += A[(1 + dim)*elemIdx + i + 1]
+                  [(1 + dim)*elemIdx + j + 1];
+            }
+          }
+
+          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[(1 + dim)*elemIdx][(1 + dim)*neighborIdx] +=
+            -mu * intersctVol;
+
+          //stiffness matrix A is symmetric
+          A[(1 + dim)*neighborIdx][(1 + dim)*elemIdx] +=
+            A[(1 + dim)*elemIdx][(1 + dim)*neighborIdx];
+
+          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[(1 + dim)*elemIdx + i + 1][(1 + dim)*neighborIdx] +=
+              -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[(1 + dim)*elemIdx][(1 + dim)*neighborIdx + i + 1] +=
+              -mu * linearIntegrals[i]
+              - 0.5 * K * normal[i] * intersctVol;
+
+            //stiffness matrix A is symmetric
+            A[(1 + dim)*neighborIdx][(1 + dim)*elemIdx + i + 1] +=
+              A[(1 + dim)*elemIdx + i + 1][(1 + dim)*neighborIdx];
+            A[(1 + dim)*neighborIdx + i + 1][(1 + dim)*elemIdx] +=
+              A[(1 + dim)*elemIdx][(1 + dim)*neighborIdx + 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[(1 + dim)*elemIdx + i + 1]
+                [(1 + dim)*neighborIdx + j + 1] +=
+                -mu * quadraticIntregrals[i][j]
+                - 0.5 * K * (normal[j] * linearIntegrals[i]
+                  - normal[i] * linearIntegrals[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[(1 + dim)*elemIdx + j + 1]
+                [(1 + dim)*neighborIdx + i + 1] +=
+                -mu * quadraticIntregrals[i][j]
+                - 0.5 * K * (normal[i] * linearIntegrals[j]
+                  - normal[j] * linearIntegrals[i]);
+
+              //stiffness matrix A is symmetric
+              A[(1 + dim)*neighborIdx + j + 1]
+                [(1 + dim)*elemIdx + i + 1] +=
+                A[(1 + dim)*elemIdx + i + 1]
+                  [(1 + dim)*neighborIdx + j + 1];
+              A[(1 + dim)*neighborIdx + i + 1]
+                [(1 + dim)*elemIdx + j + 1] +=
+                A[(1 + dim)*elemIdx + j + 1]
+                  [(1 + dim)*neighborIdx + 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[(1 + dim)*elemIdx + i + 1][(1 + dim)*elemIdx] +=
+              mu * linearIntegrals[i]
+              - 0.5 * K * normal[i] * intersctVol;
+
+            //stiffness matrix A is symmetric
+            A[(1 + dim)*elemIdx][(1 + dim)*elemIdx + i + 1] +=
+              A[(1 + dim)*elemIdx + i + 1][(1 + dim)*elemIdx];
+
+            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[(1 + dim)*elemIdx + i + 1]
+                [(1 + dim)*elemIdx + j + 1] +=
+                mu * quadraticIntregrals[i][j]
+                - 0.5 * K * (normal[i] * linearIntegrals[j]
+                  + normal[j] * linearIntegrals[i]);
+
+              //stiffness matrix A is symmetric
+              A[(1 + dim)*elemIdx + j + 1][(1 + dim)*elemIdx + i + 1]
+                += A[(1 + dim)*elemIdx + i + 1]
+                  [(1 + dim)*elemIdx + j + 1];
+            }
+          }
+        }
+      }
+    }
+
+    std::cout << A << std::endl;
+
+    //NOTE: what would be an appropiate solver here?
+    A.solve(d, b);
+
+    //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((dim + 1)*gridView_.size(0), 0.0);
+
+    for (const auto& elem : elements(gridView_))
+    {
+      const int elemIdx = mapper_.index(elem);
+      const auto& geo = elem.geometry();
+
+      for (int k = 0; k < geo.corners(); k++)
+      {
+        //contribution of the basis function
+        // phi_elem,0 (x) = indicator(elem);
+        //at the kth corner of elem
+        pressure[(dim + 1)*elemIdx + k] = d[(1 + dim)*elemIdx];
+
+        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[(dim + 1)*elemIdx + k] +=
+            d[(1 + dim)*elemIdx + i + 1] * geo.corner(k)[i];
+        }
+      }
+    }
+
+    //create output directory if necessary
+    mkdir("data", S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
+
+    auto vtkWriter =
+      std::make_shared<Dune::VTKWriter<GridView>>(gridView_);
+    Dune::VTKSequenceWriter<GridView>
+      vtkSqWriter(vtkWriter, "pressure", "data", "data");
+
+    using VTKFunction = Dune::VTKFunction<GridView>;
+    using P1Function = Dune::NonconformingP1VTKFunction<GridView,
+      Dune::DynamicVector<double>>;
+    VTKFunction* vtkOperator =
+      new P1Function(gridView_, pressure, "pressure");
+
+    //NOTE: why does this only work with a sequence writer?
+    vtkSqWriter.addVertexData(
+      std::shared_ptr<const VTKFunction>(vtkOperator) );
+    vtkSqWriter.write(0.0);
+  }
+
+  const GridView& gridView_;
+
+  const Mapper& mapper_; //element mapper for gridView_
+
+  const Problem& problem_; //the DG problem to be solved
+};
+
+#endif

dune/mmdg/problems/dgproblem.hh

+#ifndef __DUNE_MMDG_DGPROBLEM_HH__
+#define __DUNE_MMDG_DGPROBLEM_HH__
+
+//abstract class providing an interface for a problem that is to be solved with
+//a discontinuous Galerkin scheme
+template<class Coordinate, class ctype>
+class DGProblem
+{
+  public:
+    using Scalar = ctype;
+    using Vector = Coordinate;
+
+    //the exact solution at position pos
+    virtual Scalar exactSolution (const Vector& pos) const
+    {
+      return Scalar(0.0);
+    }
+
+    //indicates whether an exact solution is implemented for the problem
+    virtual bool hasExactSolution() const
+    {
+      return false;
+    };
+
+    //source term
+    virtual Scalar q (const Vector& pos) const //= 0;
+    {
+      return Scalar(0.0);
+    }
+
+    //boundary condition at position pos
+    virtual Scalar boundary (const Vector& pos) const //= 0;
+    {
+      return Scalar(0.0);
+    }
+};
+
+#endif

src/dg.cc

@@ -4,6 +4,7 @@
 
 #include <sys/stat.h>
 #include <sys/types.h>
+#include <type_traits>
 
 #include <dune/common/parallel/mpihelper.hh>
 #include <dune/common/exceptions.hh>
@@ -16,13 +17,15 @@
 #include <dune/grid/io/file/dgfparser/dgfyasp.hh>
 #include <dune/grid/io/file/dgfparser/dgfug.hh>
 #include <dune/grid/io/file/vtk.hh>
+#include <dune/grid/yaspgrid.hh>
 
 #include <dune/geometry/quadraturerules.hh>
 
 #include <dune/mmesh/mmesh.hh>
 
+#include <dune/mmdg/dg.hh>
 #include <dune/mmdg/clockhelper.hh>
-#include <dune/mmdg/nonconformingp1vtkfunction.hh>
+#include <dune/mmdg/problems/dgproblem.hh>
 
 int main(int argc, char** argv)
 {
@@ -35,7 +38,8 @@ int main(int argc, char** argv)
 
     static constexpr int dim = GRIDDIM;
 
-    using Grid = Dune::MovingMesh<dim>;
+    using Grid = std::conditional<dim == 1, Dune::YaspGrid<dim>,
+      Dune::MovingMesh<dim>>::type;
     using GridFactory = Dune::DGFGridFactory<Grid>;
     using GridView = typename Grid::LeafGridView;
     using Mapper = typename Dune::MultipleCodimMultipleGeomTypeMapper<GridView>;
@@ -63,303 +67,16 @@ int main(int argc, char** argv)
     //element mapper for the leaf grid view
     const Mapper mapper(gridView, Dune::mcmgElementLayout());
 
-    const int dof = (1 + dim)*gridView.size(0); //degrees of freedom
-    using Matrix = Dune::DynamicMatrix<double>; //Dune::FieldMatrix<double, dof, dof>;
-    using Vector = Dune::DynamicVector<double>; //Dune::FieldVector<double, dof>;
-    Matrix A(dof, dof, 0.0); //stiffness matrix
-    Vector b(dof, 0.0);
-    Vector d;
+    const DGProblem<Dune::FieldVector<double, dim>, double> dummyProblem;
 
-    const auto q = [](const Coordinate& x) {return x.two_norm();}; //source term
+    const DG dg(gridView, mapper, dummyProblem);
 
-    const int order = 3; //order of the quadrature rule
-
-    //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();
-
-      //TODO: can be done outside of the loop?
-      const Dune::QuadratureRule<double,dim>& rule =
-        Dune::QuadratureRules<double,dim>::rule(geo.type(),order);
-
-      for (const auto& ip : rule)
-      {
-        const auto qp = ip.position();
-        //NOTE: is the volume of the element taken into account automatically?
-        const auto weight = ip.weight();
-
-        //quadrature for int_elem q*phi_elem,0 dV
-        b[(1 + dim)*elemIdx] += weight * q(qp);
-
-        //quadrature for int_elem q*phi_elem,i dV
-        for (int i = 0; i < dim; i++)
-        {
-          b[(1 + dim)*elemIdx + i + 1] += weight * qp[i] * q(qp);
-        }
-      }
-
-      for (int i = 0; i < dim; i++)
-      {
-        //exact evaluation of
-        // int_elem K*grad(phi_elem,i)*grad(phi_elem,i) dV
-        A[(1 + dim)*elemIdx + i + 1][(1 + dim)*elemIdx + 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<double, dim> linearIntegrals(0.0);
-        Dune::FieldMatrix<double, dim, dim> quadraticIntregrals(0.0);
-
-        for (int i = 0; i < dim; i++)
-        {
-          //use midpoint rule for exact evaluation of
-          // int_intersct x_i ds
-          linearIntegrals[i] = mu * intersctVol * 0.5 * intersctCenter[i];
-
-          for (int j = 0; j <= i; j++)
-          {
-            //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;
-            }
-          }
-        }
-
-        //exact evaluation of
-        // int_intersct mu*jump(phi_elem,0)*jump(phi_elem,0) ds
-        A[(1 + dim)*elemIdx][(1 + dim)*elemIdx] += mu * intersctGeo.volume();
-
-        if (intersct.neighbor()) //intersct has neighboring element
-        {
-          //index of the neighboring element
-          const int neighborIdx = mapper.index(intersct.outside());
-
-          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[(1 + dim)*elemIdx + i + 1][(1 + dim)*elemIdx] =
-              mu * linearIntegrals[i] - 0.5 * K * normal[i] * intersctVol;
-
-            //stiffness matrix A is symmetric
-            A[(1 + dim)*elemIdx][(1 + dim)*elemIdx + i + 1] =
-              A[(1 + dim)*elemIdx + i + 1][(1 + dim)*elemIdx];
-
-            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[(1 + dim)*elemIdx + i + 1][(1 + dim)*elemIdx + j + 1] =
-                mu * quadraticIntregrals[i][j] - 0.5 * K * (normal[i] *
-                  linearIntegrals[j] + normal[j] * linearIntegrals[i]);
-
-              //stiffness matrix A is symmetric
-              A[(1 + dim)*elemIdx + j + 1][(1 + dim)*elemIdx + i + 1] =
-                A[(1 + dim)*elemIdx + i + 1][(1 + dim)*elemIdx + j + 1];
-            }
-          }
-
-          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[(1 + dim)*elemIdx][(1 + dim)*neighborIdx] = -mu * intersctVol;
-
-          //stiffness matrix A is symmetric
-          A[(1 + dim)*neighborIdx][(1 + dim)*elemIdx] =
-            A[(1 + dim)*elemIdx][(1 + dim)*neighborIdx];
-
-          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[(1 + dim)*elemIdx + i + 1][(1 + dim)*neighborIdx] =
-              -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[(1 + dim)*elemIdx][(1 + dim)*neighborIdx + i + 1] =
-              -mu * linearIntegrals[i] - 0.5 * K * normal[i] * intersctVol;
-
-            //stiffness matrix A is symmetric
-            A[(1 + dim)*neighborIdx][(1 + dim)*elemIdx + i + 1] =
-              A[(1 + dim)*elemIdx + i + 1][(1 + dim)*neighborIdx];
-            A[(1 + dim)*neighborIdx + i + 1][(1 + dim)*elemIdx] =
-              A[(1 + dim)*elemIdx][(1 + dim)*neighborIdx + 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[(1 + dim)*elemIdx + i + 1][(1 + dim)*neighborIdx + j + 1] =
-                -mu * quadraticIntregrals[i][j]
-                - 0.5 * K * (normal[j] * linearIntegrals[i]
-                  - normal[i] * linearIntegrals[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[(1 + dim)*elemIdx + j + 1][(1 + dim)*neighborIdx + i + 1] =
-                -mu * quadraticIntregrals[i][j]
-                - 0.5 * K * (normal[i] * linearIntegrals[j]
-                  - normal[j] * linearIntegrals[i]);
-
-              //stiffness matrix A is symmetric
-              A[(1 + dim)*neighborIdx + j + 1][(1 + dim)*elemIdx + i + 1] =
-                A[(1 + dim)*elemIdx + i + 1][(1 + dim)*neighborIdx + j + 1];
-              A[(1 + dim)*neighborIdx + i + 1][(1 + dim)*elemIdx + j + 1] =
-                A[(1 + dim)*elemIdx + j + 1][(1 + dim)*neighborIdx + 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[(1 + dim)*elemIdx + i + 1][(1 + dim)*elemIdx] =
-              mu * linearIntegrals[i] - K * normal[i] * intersctVol;
-
-            //stiffness matrix A is symmetric
-            A[(1 + dim)*elemIdx][(1 + dim)*elemIdx + i + 1] =
-              A[(1 + dim)*elemIdx + i + 1][(1 + dim)*elemIdx];
-
-            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[(1 + dim)*elemIdx + i + 1][(1 + dim)*elemIdx + j + 1] =
-                mu * quadraticIntregrals[i][j] - K * (normal[i] *
-                  linearIntegrals[j] + normal[j] * linearIntegrals[i]);
-
-              //stiffness matrix A is symmetric
-              A[(1 + dim)*elemIdx + j + 1][(1 + dim)*elemIdx + i + 1] =
-                A[(1 + dim)*elemIdx + i + 1][(1 + dim)*elemIdx + j + 1];
-            }
-          }
-        }
-      }
-    }
-
-    A.solve(d, b);
-
-    //storage for pressure data, for each element we store the pressure at the
-    //corners of the element, the VTKOperatorP1 will create the output using a
-    //linear interpolation for each element
-    Dune::DynamicVector<double> pressure(3*gridView.size(0), 0.0);
-
-    for (const auto& elem : elements(gridView))
-    {
-      const int elemIdx = mapper.index(elem);
-      const auto& geo = elem.geometry();
-
-      for (int k = 0; k < geo.corners(); k++)
-      {
-        //contribution of the basis function
-        // phi_elem,0 (x) = indicator(elem);
-        //at the kth corner of elem
-        pressure[3*elemIdx + k] = d[(1 + dim)*elemIdx];
-
-        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[3*elemIdx + k] +=
-            d[(1 + dim)*elemIdx + i + 1] * geo.corner(k)[i];
-        }
-      }
-    }
-
-/*
-    typedef Dune::VTKFunction< GridView > VTKFunction;
-    typedef Dune::NonconformingP1VTKFunction<GridView,
-    Dune::DynamicVector<double>> P1Function;
-    VTKFunction* p1 = new P1Function(gridView, lp, "linearpressure");
-    vtkwriter.addVertexData( std::shared_ptr< const VTKFunction >(p1) );
-  */
-
-    //create output directory if necessary
-    mkdir("data", S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
-
-    auto vtkWriter = std::make_shared<Dune::VTKWriter<GridView>>(gridView);
-    Dune::VTKSequenceWriter<GridView>
-      vtkSqWriter(vtkWriter, "pressure", "data", "data");
-
-    using VTKOperatorP1 = Dune::VTKFunction<GridView>;
-    using P1Function = Dune::NonconformingP1VTKFunction<GridView,
-      Dune::DynamicVector<double>>;
-    VTKOperatorP1* vtkOperator = new P1Function(gridView, pressure, "pressure");
-
-
-    vtkSqWriter.addVertexData( std::shared_ptr<const VTKOperatorP1>(vtkOperator) );
-    vtkSqWriter.write(0.0);
+    dg(K, mu);
 
     //print total run time
     const double timeTotal = ClockHelper::elapsedRuntime(timeInitial);
-    std::cout << "Finished with a total run time of " << timeTotal << ".\n";
+    std::cout << "Finished with a total run time of " << timeTotal <<
+      " Seconds.\n";
 
     return 0;
   }

src/grids/cube1d.dgf

@@ -3,7 +3,7 @@ DGF
 Interval
 0	
 1	
-5	
+20	
 #
 
 Simplex

src/grids/cube2d.dgf

@@ -3,7 +3,7 @@ DGF
 Interval
 0	0	
 1	1	
-2	2	
+20	20	
 #
 
 Simplex