[bugfix] interface basis is now calculated using a Householder matrix

dune/mmdg/mmdg.hh

@@ -41,7 +41,7 @@ public:
     Base::A->compress();
 
     Dune::InverseOperatorResult result;
-    Dune::UMFPack<Matrix> solver(*Base::A);
+    Dune::UMFPack<Matrix> solver(*Base::A); //TODO: ilu0bicgSTAB, cg (ilo0, ssor), gmres
     solver.apply(Base::d, Base::b, result);
 
     //write solution to a vtk file
@@ -832,36 +832,46 @@ private:
     }
   }
 
-  //returns a basis of the interface coordinate system
-  //2d: tau = [ normal[1],
-  //           -normal[0] ]
-  //3d: tau_1 ~ [ normal[1] + normal[2],     tau_2 ~ [ -normal[1],
-  //             -normal[0],                            normal[0] + normal[2]
-  //             -normal[0]             ],             -normal[1]            ]
-  InterfaceBasis interfaceFrame(const Coordinate& normal) const
+  //returns an orthonormal basis of the interface coordinate system
+  InterfaceBasis interfaceFrame (Coordinate normal) const
   {
-    InterfaceBasis iBasis;
+    //find a standard unit vector e_k with k such that |a[k]| != 1
+    int k = 0;
 
-    for (int i = 0; i < dim - 1; i++)
-    {
-      Coordinate tau(0.0);
-      for (int j = 0; j < dim; j++)
+    for (; k < dim; k++)
     {
-        if (j == i)
+      if (std::abs(normal[k]) < 0.9)
       {
-          continue;
+        normal[k] -= 1;
+        break;
       }
-
-        tau[i] += normal[j];
-        tau[j] = -normal[i];
     }
 
-      if constexpr (dim != 2)
+    //regard the unit vector v := (normal - e_k) / norm(normal - e_k)
+    normal /= normal.two_norm();
+
+    //calculate the Householder matrix H(v) = Id - 2*v*v^T
+    CoordinateMatrix householderMatrix(0.0);
+
+    for (int i = 0; i < dim; i++)
     {
-        tau /= tau.two_norm();
+      householderMatrix[i][i] = 1.0;
     }
 
-      iBasis[i] = tau;
+    householderMatrix -= 2.0 * outerProduct(normal);
+
+    //the Householder transformation fulfills the property
+    // H(v) * e_k = +/- normal
+    //since H(v) is orthonormal, an (orthonormal) interface basis is now
+    //given by the columns H(v) * e_j with j != k
+    InterfaceBasis iBasis;
+
+    for (int i = 0; i < dim; i++)
+    {
+      if (i != k)
+      {
+        iBasis[(k > i) ? i : i-1] = householderMatrix[i];
+      }
     }
 
     return iBasis;
@@ -953,6 +963,22 @@ private:
     return mu * mu0 * 2.0 * (1.0 + dim);
   }
 
+  //returns the outer product v * v^T of the vector v
+  const CoordinateMatrix outerProduct (const Coordinate& v) const
+  {
+    CoordinateMatrix vv;
+
+    for (int i = 0; i < dim; i++)
+    {
+      for (int j = 0; j < dim; j++)
+      {
+        vv[i][j] = v[i] * v[j];
+      }
+    }
+
+    return vv;
+  }
+
 };
 
 #endif