Added comments

src/dune-simon.cc

@@ -25,6 +25,7 @@
 using namespace Dune;
 
 
+// Generate the index set Z_{eta,>}
 void indexSet(const int& eta, int& c, std::vector<int>& k0, std::vector<int>& k1)
 {
   c = 0;
@@ -51,6 +52,7 @@ void indexSet(const int& eta, int& c, std::vector<int>& k0, std::vector<int>& k1
 }
 
 
+// Generate a random graph function
 template <class VectorType>
 void generateSurfaceFunction(int& c, std::vector<int>& k0, std::vector<int>& k1, VectorType& xk,
   VectorType& yk, const double& kappa, const double& sigma)
@@ -73,6 +75,7 @@ void generateSurfaceFunction(int& c, std::vector<int>& k0, std::vector<int>& k1,
 }
 
 
+// Evaluate the first fundamental form of the surface at a given point x
 template <class CoordinateType, class VectorType, class MatrixType>
 void getFundamentalMatrix(CoordinateType& x, MatrixType& fundamentalmatrix, int& c,
   std::vector<int>& k0, std::vector<int>& k1, VectorType& xk, VectorType& yk, double& pi)
@@ -99,6 +102,7 @@ void getFundamentalMatrix(CoordinateType& x, MatrixType& fundamentalmatrix, int&
 }
 
 
+// Assemble the local stiffness and mass matrices
 template <class LocalView, class VectorType, class MatrixType>
 void getLocalMatrices(const LocalView& localview, MatrixType& localstiffnessmatrix,
   MatrixType& localmassmatrix, int& c, std::vector<int>& k0, std::vector<int>& k1,
@@ -177,6 +181,7 @@ void getLocalMatrices(const LocalView& localview, MatrixType& localstiffnessmatr
 }
 
 
+// Return true/false depending on if the index is on the right or top boundary of the domain
 template <class IndexType>
 bool isOnBoundary(IndexType& index, const int& elementnumber, const int& feorder)
 {
@@ -205,6 +210,7 @@ bool isOnBoundary(IndexType& index, const int& elementnumber, const int& feorder
 }
 
 
+// Shift index, if it is on the right/top boundary, to the identified one on the left/bottom boundary
 template <class IndexType>
 void indexCorrection(IndexType& index, const int& elementnumber, const int& feorder)
 {
@@ -276,6 +282,7 @@ void getOccupationPattern(const FEBasis& febasis, MatrixIndexSet& matrixindexset
 
   }
 
+  // Add nonzero entries on the diagonal of the boundary nodes
   for(int i = 0; i < n; i++){
     if (isOnBoundary(i, elementnumber, feorder)){
       matrixindexset.add(i,i);
@@ -352,6 +359,7 @@ void assembleMatrices(const FEBasis& febasis, MatrixType& stiffnessmatrix, Matri
 
   }
 
+  // Make the stiffness matrix to the identity for all boundary nodes to ensure that the final matrix is invertible
   for(int i = 0; i < n; i++){
     if (isOnBoundary(i, elementnumber, feorder)){
       stiffnessmatrix[i][i] = 1;
@@ -361,6 +369,7 @@ void assembleMatrices(const FEBasis& febasis, MatrixType& stiffnessmatrix, Matri
 }
 
 
+// Copy all values of a vector to the identified boundary points
 template <class VectorType>
 void copyToBoundary(VectorType& vector, const int& n, const int& elementnumber, const int& feorder)
 {
@@ -374,6 +383,7 @@ void copyToBoundary(VectorType& vector, const int& n, const int& elementnumber,
 }
 
 
+// Interpolate the surface function in the grid view
 template <class FEBasis, class GridView, class VectorType>
 void interpolateSurfaceFunction(const FEBasis& febasis, GridView& gridview, VectorType& surface,
   std::vector<char>& nonboundarynodes, const int& elementnumber, const int& feorder, int& c, std::vector<int>& k0,
@@ -400,6 +410,7 @@ void interpolateSurfaceFunction(const FEBasis& febasis, GridView& gridview, Vect
 }
 
 
+// Compute one implicit Euler step
 template <class FEBasis, class Solver, class VectorType, class MatrixType>
 void implicitEuler(const FEBasis& febasis, Solver& solver, VectorType& iter,
   MatrixType& massmatrix)