add separate EOCs for bulk and interface

dune/mmdg/mmdg.hh

@@ -48,9 +48,11 @@ public:
     writeVTKoutput();
   }
 
-  //compute L2 error using quadrature rules
+  //compute L2 errors using quadrature rules,
+  //returns {error(bulk), error(interface), error(total)},
+  //the total error is calculated using the norm
   // norm = sqrt( ||.||^2_L2(bulk) + ||.||^2_L2(interface) )
-  Scalar computeL2error(const int quadratureOrder = 5)
+  const std::array<Scalar, 3> computeL2error(const int quadratureOrder = 5)
   {
     const Scalar bulkError = Base::computeL2error(quadratureOrder);
     Scalar interfaceError(0.0);
@@ -87,11 +89,8 @@ public:
       Base::applyQuadratureRule(iGeo, qrule, qlambda);
     }
 
-    //print bulk and interface error
-    std::cout << "ERROR bulk: " << bulkError << "\t ERROR interface: " <<
-      sqrt(interfaceError) << "\n";
-
-    return sqrt(bulkError * bulkError + interfaceError);
+    return {bulkError, sqrt(interfaceError),
+      sqrt(bulkError * bulkError + interfaceError)};
   }
 
   const InterfaceGridView& iGridView_;

src/mmdg.cc

@@ -131,23 +131,33 @@ int main(int argc, char** argv)
     MMDG mmdg(gridView, mapper, iGridView, iMapper, *problem);
     mmdg(mu0, xi);
 
-    //error and total run time
-    const double error = mmdg.computeL2error();
+    //errors and total run time
+    const std::array<double, 3> errors = mmdg.computeL2error();
     const double timeTotal = ClockHelper::elapsedRuntime(timeInitial);
 
-    //determine maximum edge length
-    double largestEdgeLength = 0.;
+    //determine maximum edge lengths
+    double hMax = 0.0; //maximum edge length (bulk)
+    double hMaxInterface = 0.0; //maximum edge length (interface)
+
     for (const auto& edge : edges(gridView))
     {
-      largestEdgeLength = std::max(largestEdgeLength, edge.geometry().volume());
+      if (grid.isInterface(edge))
+      {
+        hMaxInterface = std::max(hMaxInterface, edge.geometry().volume());
+      }
+      else
+      {
+        hMax = std::max(hMax, edge.geometry().volume());
+      }
     }
 
-    //write error, runtime and maximum edge length
+    //write errors, runtime and maximum edge lengths
     std::ofstream errorFile("convergenceData", std::ios::out | std::ios::trunc);
 
     if (errorFile.is_open())
     {
-      errorFile << error << "\n" << timeTotal << "\n" << largestEdgeLength;
+      errorFile << errors[0] << "\n" << errors[1] << "\n" << errors[2] << "\n"
+        << hMax << "\n" << hMaxInterface << "\n" << timeTotal;
       errorFile.close();
     }
     else
@@ -155,9 +165,11 @@ int main(int argc, char** argv)
       std::cout << "Unable to write error data to file.\n";
     }
 
-    //print total run time and error
+    //print total run time and errors
     std::cout << "Finished with a total run time of " << timeTotal <<
-      " Seconds.\nL2-error: " << error << ".\n\n";
+      " Seconds.\nL2-error (bulk): \t" << errors[0] <<
+      "\nL2-error (interface): \t" << errors[1] << "\nL2-error (total): \t"
+      << errors[2] << ".\n\n";
 
     return 0;
   }

src/mmdgAnalysis.py

@@ -41,14 +41,36 @@ def readData():
     file = open("convergenceData", "r")
 
     # read data
-    error = float(file.readline())
+    errorBulk = float(file.readline())
+    errorInterface = float(file.readline())
+    errorTotal = float(file.readline())
+    hMax = float(file.readline())
+    hMaxInterface = float(file.readline())
     timeTotal = float(file.readline())
-    numberOfElements = float(file.readline())
 
     file.close()
 
     # return data
-    return np.array([error, timeTotal, numberOfElements])
+    return np.array([errorBulk, errorInterface, errorTotal, hMax, \
+        hMaxInterface, timeTotal])
+
+
+# determines experimental order of convergence (EOC),
+# data is assumed to be a numpy array with
+# first column: error,
+# second column: maximum edge length
+def EOCloop (data):
+    EOC = np.empty(data.shape[0] - 1)
+
+    for i in range(len(EOC)):
+        if data[i, 0] != 0:
+            EOC[i] = log(data[i, 0] / data[i+1, 0]) / \
+                log(data[i, 1] / data[i+1, 1])
+        else:
+            EOC[i] = -1.0*float('inf')
+    # end for i
+
+    return EOC
 
 
 # === main ===
@@ -72,7 +94,7 @@ dim = 2
 # column 1: total run time
 # column 2: largest edge length
 # different rows correspond to different numbers of grid elements
-data = np.zeros((len(gridfiles), 3))
+data = np.zeros((len(gridfiles), 6))
 
 for i in range(len(gridfiles)):
     # write dgf file
@@ -92,25 +114,25 @@ for i in range(len(gridfiles)):
 data /= numberOfMeans
 
 # determine experimental order of convergence (EOC)
-EOC = np.empty(len(gridfiles) - 1)
-
-for i in range(len(gridfiles) - 1):
-    if data[i, 0] != 0:
-        EOC[i] = log(data[i, 0] / data[i+1, 0]) / \
-            log(data[i, 2] / data[i+1, 2])
-    else:
-        EOC[i] = -1.0*float('inf')
-# end for i
+EOCbulk = EOCloop(data[:, [0,3]])
+EOCinterface = EOCloop(data[:, [1,4]])
+EOCtotal = \
+    EOCloop(np.column_stack([data[:,2], np.maximum(data[:,3], data[:,4])]))
 
 # print EOCs
 print("\n-------------------------------------------------------------\n")
-print("EOCs (" + str(dim) + "d):\t")
-print(EOC)
+print("EOCs (" + str(dim) + "d):")
+print("\nbulk:")
+print(EOCbulk)
+print("\ninterface:")
+print(EOCinterface)
+print("\ntotal:")
+print(EOCtotal)
 print("\n=============================================================\n\n")
 
 # write data to file
 np.savetxt(join("plots", "mmdgAnalysis_" + str(dim) + "d"), data)
-np.savetxt(join("plots", "mmdgEOC_" + str(dim) + "d"), EOC)
+np.savetxt(join("plots", "mmdgEOC_" + str(dim) + "d"), EOCtotal)
 
 # remove output files
 remove("convergenceData")