[WIP] continue implementing mmdg

dune/mmdg/mmdg.hh

@@ -28,9 +28,9 @@ public:
       (1 + dim) * gridView.size(0) + dim * iGridView.size(0)),
     iGridView_(iGridView), iMapper_(iMapper) {}
 
-  const void operator() (const Scalar mu)
+  const void operator() (const Scalar mu, const Scalar xi)
   {
-    assembleSLE(mu);
+    assembleSLE(mu, xi);
     Base::A->compress();
   }
 
@@ -39,11 +39,28 @@ public:
 
 private:
   //assemble stiffness matrix A and load vector b
-  void assembleSLE (const Scalar mu)
+  void assembleSLE (const Scalar mu, const Scalar xi)
   {
+    //modell parameters //NOTE: what should be the capture of alpha?
+    const auto alpha = [this](const Coordinate& pos) -> Scalar
+    {
+      return Base::problem_.Kperp(pos) / Base::problem_.aperture(pos);
+    };
+
+    const auto beta = [&alpha, &xi](const Coordinate& pos) -> Scalar
+    {
+      return 4 * alpha(pos) / (2 * xi - 1);
+    };
+
     //quadrature rules
     const QRuleInterface& rule_qInterface = QRulesInterface::rule(iSimplex,
       Base::problem_.quadratureOrder_qInterface());
+    const QRuleInterface& rule_Kpar = QRulesInterface::rule(iSimplex,
+      Base::problem_.quadratureOrder_Kparallel());
+    const QRuleInterface& rule_Kperp = QRulesInterface::rule(iSimplex,
+      Base::problem_.quadratureOrder_Kperp());
+    const QRuleInterface& rule_KperpPlus1 = QRulesInterface::rule(iSimplex,
+      Base::problem_.quadratureOrder_Kperp() + 1);
 
     //we use the bulk basis
     // phi_elem,0 (x) = indicator(elem);
@@ -75,10 +92,11 @@ private:
     {
       const int iElemIdx = iMapper_.index(iElem);
       const auto& iGeo = iElem.geometry();
+      const auto& intersct = Base::grid_.asIntersection(iElem);
 
       //get basis of the interface coordinate system for evaluation of
       //basis function phi_iElem,i with i=1,...,dim-1
-      const std::array<Coordinate, dim-1> iFrame = interfaceFrame(
+      const InterfaceBasis iFrame = interfaceFrame(
         Base::grid_.asIntersection(iElem).centerUnitOuterNormal());
 
       //in the total system of linear equations (SLE) Ad = b,
@@ -111,11 +129,80 @@ private:
 
       //TODO: fill stiffness matrix A with interface entries
 
-      //fill stiffness matrix A with coupling entries
-      for (const auto& elem : elements(Base::gridView_))
+      //evaluate
+      // int_iElem d * (K_par*grad_par(phi_iElem,i)) * grad_par(phi_iElem,j) ds
+      // = int_iElem d * (K_par * tau_i) * tau_j ds
+      //using a quadrature rule for i,j = 1,...,dim-1
+      for (const auto& ip : rule_Kpar)
+      {
+        const auto& qpLocal = ip.position();
+        const auto& qpGlobal = iGeo.global(qpLocal);
+
+        const Scalar weight(ip.weight());
+        const Scalar integrationElem(iGeo.integrationElement(qpLocal));
+
+        for (int i = 0; i < dim-1 ; i++)
+        {
+          Coordinate K_parDotTau_i;
+          Base::problem_.Kparallel(qpGlobal).mv(iFrame[i], K_parDotTau_i);
+
+          for (int j = 0; j <= i; j++) //TODO: symmetrize
+          {
+            Base::A->entry(iElemIdxSLE + i + 1, iElemIdxSLE + j + 1) +=
+              weight * integrationElem * Base::problem_.aperture(qpGlobal)
+              * ( K_parDotTau_i * iFrame[j] );
+          }
+        }
+      }
+
+      //TODO: check interface bilinear form
+
+      //evaluate the coupling term
+      // int_iElem beta * phi_iElem,i * phi_iElem,j ds
+      //using a quadrature rule for i = 0 and j = 0,...,dim-1 //TODO or vice versa (symmetrize)
+      for (const auto& ip : rule_Kperp)
+      {
+        const auto& qpLocal = ip.position();
+        const auto& qpGlobal = iGeo.global(qpLocal);
+
+        const Scalar weight(ip.weight());
+        const Scalar integrationElem(iGeo.integrationElement(qpLocal));
+
+        Base::A->entry(iElemIdxSLE, iElemIdxSLE) +=
+          weight * integrationElem * beta(qpGloabal);
+
+        for (int i = 0; i < dim - 1; i++)
+        {
+          Base::A->entry(iElemIdxSLE + i + 1, iElemIdxSLE + j + 1) +=
+            weight * integrationElem * beta(qpGlobal) * (iFrame[i] * qpGlobal)
+            * (iFrame[j] * qpGlobal);
+        }
+      }
+
+      //evaluate the coupling term
+      // int_iElem beta * phi_iElem,i * phi_iElem,j ds
+      //using a quadrature rule for i,j = 1,...,dim-1 //TODO symmetrize
+      for (const auto& ip : rule_KperpPlus1)
       {
-        //TODO
+        const auto& qpLocal = ip.position();
+        const auto& qpGlobal = iGeo.global(qpLocal);
+
+        const Scalar weight(ip.weight());
+        const Scalar integrationElem(iGeo.integrationElement(qpLocal));
+
+        for (int i = 0; i < dim - 1; i++)
+        {
+          for (int j = 0; j <= i; j++)
+          {
+            Base::A->entry(iElemIdxSLE + i + 1, iElemIdxSLE) +=
+              weight * integrationElem * beta(qpGlobal) * (iFrame[i] * qpGlobal) ;
+          }
+        }
       }
+
+      const auto elem1 = intersct.inside();
+      const auto elem2 = intersct.ouside();
+      //TODO coupling entries
     }
   }
 

dune/mmdg/problems/coupleddgproblem.hh

@@ -20,12 +20,18 @@ class CoupledDGProblem : public DGProblem<Vector, Scalar>
       return Scalar(0.0);
     }
 
-    //aperture of the fracture at position pos
+    //aperture d of the fracture at position pos
     virtual Scalar aperture (const Vector& pos) const
     {
       return Scalar(1.0);
     }
 
+    //(tangential) gradient of the aperture at position pos
+    virtual Vector gradAperture (const Vector& pos) const
+    {
+      return Vector(0.0);
+    }
+
     //tangential permeability tensor of the interface at position pos
     virtual Matrix Kparallel (const Vector& pos) const
     {
@@ -53,6 +59,20 @@ class CoupledDGProblem : public DGProblem<Vector, Scalar>
     {
       return 1;
     }
+
+    //returns the recommended quadrature order to compute an integral
+    //over d(x) * Kparallel(x)
+    virtual int quadratureOrder_Kparallel () const
+    {
+      return 1;
+    }
+
+    //returns the recommended quadrature order to compute an integral
+    //over x * Kperp(x) / d(x)
+    virtual int quadratureOrder_Kperp () const
+    {
+      return 1;
+    }
 };
 
 #endif

src/mmdg.cc

@@ -48,8 +48,9 @@ int main(int argc, char** argv)
     Dune::ParameterTree pt;
     Dune::ParameterTreeParser::readINITree("parameterMMDG.ini", pt);
 
-    //assume constant penalty parameter
+    //assume constant penalty parameter //TODO
     const double mu = std::stod(pt["mu"]);
+    const double xi = std::stod(pt["xi"]); //coupling parameter
 
     //create a grid from .dgf file
     GridFactory gridFactory( "grids/cube" + std::to_string(dim) + "d.dgf" );
@@ -77,7 +78,7 @@ int main(int argc, char** argv)
     }
 
     MMDG mmdg(grid, gridView, mapper, iGridView, iMapper, *problem);
-    mmdg(mu);
+    mmdg(mu, xi);
 
     //print total run time
     const double timeTotal = ClockHelper::elapsedRuntime(timeInitial);

src/parameterMMDG.ini

 mu = 1000
+xi = 0.75
 problem = poisson