diff --git a/src/dune-simon.cc b/src/dune-simon.cc
index 8394fd256738781b541fcdd7549d1ab96e739311..c55627efddef344e57d0712533be409597b8d7f9 100644
--- a/src/dune-simon.cc
+++ b/src/dune-simon.cc
@@ -26,9 +26,10 @@ using namespace Dune;
// Generate the index set Z_{eta,>}
-void indexSet(const int& eta, int& c, std::vector<int>& k0, std::vector<int>& k1)
+void indexSet(const int& eta, int& modenumber, std::vector<int>& k0,
+ std::vector<int>& k1)
{
- c = 0;
+ modenumber = 0;
for(int i = 0; i <= eta; i++){
int s;
@@ -43,7 +44,7 @@ void indexSet(const int& eta, int& c, std::vector<int>& k0, std::vector<int>& k1
auto norm = std::sqrt(i*i + j*j);
if (norm <= eta){
- c += 1;
+ modenumber += 1;
k0.push_back(i);
k1.push_back(j);
}
@@ -54,19 +55,19 @@ 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)
+void generateSurfaceFunction(int& modenumber, std::vector<int>& k0, std::vector<int>& k1,
+ VectorType& xk, VectorType& yk, const double& kappa, const double& sigma)
{
double pi = std::acos(-1.0);
- xk.resize(c);
- yk.resize(c);
+ xk.resize(modenumber);
+ yk.resize(modenumber);
std::random_device rd{};
std::mt19937 gen{rd()};
std::normal_distribution<> d;
- for(int i = 0; i < c; i++){
+ for(int i = 0; i < modenumber; i++){
auto normsq = k0[i]*k0[i] + k1[i]*k1[i];
double factor = std::sqrt(2*normsq*(4*pi*pi*kappa*normsq + sigma));
xk[i] = d(gen)/factor;
@@ -77,17 +78,19 @@ 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)
+void getFundamentalMatrix(CoordinateType& x, MatrixType& fundamentalmatrix,
+ int& modenumber, std::vector<int>& k0, std::vector<int>& k1, VectorType& xk,
+ VectorType& yk, double& pi)
{
double sp;
double factor;
std::vector<double> gradient;
+
gradient.resize(2);
gradient[0] = 0;
gradient[1] = 0;
- for(int i=0; i < c; i++){
+ for(int i=0; i < modenumber; i++){
sp = 2*pi*(x[0]*k0[i] + x[1]*k1[i]);
factor = std::cos(sp)*xk[i] + std::sin(sp)*yk[i];
@@ -105,8 +108,8 @@ 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,
- VectorType& xk, VectorType& yk, double& pi)
+ MatrixType& localmassmatrix, int& modenumber, std::vector<int>& k0,
+ std::vector<int>& k1, VectorType& xk, VectorType& yk, double& pi)
{
// Get the grid element from the local FE basis view
typedef typename LocalView::Element Element;
@@ -159,7 +162,7 @@ void getLocalMatrices(const LocalView& localview, MatrixType& localstiffnessmatr
localFiniteElement.localBasis().evaluateFunction(quadPos, values);
FieldMatrix<double,dim,dim> fundamentalMatrix;
- getFundamentalMatrix(globalQuadPos, fundamentalMatrix, c, k0, k1, xk, yk, pi);
+ getFundamentalMatrix(globalQuadPos, fundamentalMatrix, modenumber, k0, k1, xk, yk, pi);
auto transformationFactor = std::sqrt(fundamentalMatrix.determinant());
auto transformationMatrix = fundamentalMatrix;
transformationMatrix.invert();
@@ -245,8 +248,8 @@ void indexCorrection(IndexType& index, const int& elementnumber, const int& feor
// Get the occupation pattern of the stiffness/mass matrix
template <class FEBasis>
-void getOccupationPattern(const FEBasis& febasis, MatrixIndexSet& matrixindexset, const int& elementnumber,
- const int& feorder)
+void getOccupationPattern(const FEBasis& febasis, MatrixIndexSet& matrixindexset,
+ const int& elementnumber, const int& feorder)
{
// Total number of grid vertices
const int n = febasis.size();
@@ -294,9 +297,10 @@ void getOccupationPattern(const FEBasis& febasis, MatrixIndexSet& matrixindexset
// Assemble the Laplace stiffness and mass matrix on the given grid view
template <class FEBasis, class VectorType, class MatrixType>
-void assembleMatrices(const FEBasis& febasis, MatrixType& stiffnessmatrix, MatrixType& massmatrix,
- const int& elementnumber, const int& feorder, int& c, std::vector<int>& k0, std::vector<int>& k1, VectorType& xk,
- VectorType& yk, double& pi)
+void assembleMatrices(const FEBasis& febasis, MatrixType& stiffnessmatrix,
+ MatrixType& massmatrix, const int& elementnumber, const int& feorder, int& modenumber,
+ std::vector<int>& k0, std::vector<int>& k1, VectorType& xk, VectorType& yk,
+ double& pi)
{
const int n = febasis.size();
stiffnessmatrix.setSize(n,n);
@@ -332,10 +336,10 @@ void assembleMatrices(const FEBasis& febasis, MatrixType& stiffnessmatrix, Matri
// Now let's get the element stiffness matrix
// A dense matrix is used for the element stiffness matrix
Matrix<FieldMatrix<double,1,1>> localStiffnessMatrix;
- // getLocalStiffnessMatrix(localView, localstiffnessmatrix, c, k0, k1, xk, yk);
+ // getLocalStiffnessMatrix(localView, localstiffnessmatrix, modenumber, k0, k1, xk, yk);
Matrix<FieldMatrix<double,1,1>> localMassMatrix;
- // getLocalMassMatrix(localView, localmassmatrix, c, k0, k1, xk, yk);
- getLocalMatrices(localView, localStiffnessMatrix, localMassMatrix, c, k0, k1, xk, yk, pi);
+ // getLocalMassMatrix(localView, localmassmatrix, modenumber, k0, k1, xk, yk);
+ getLocalMatrices(localView, localStiffnessMatrix, localMassMatrix, modenumber, k0, k1, xk, yk, pi);
// Add element stiffness matrix onto the global stiffness matrix
for (size_t i=0; i<localStiffnessMatrix.N(); i++) {
@@ -369,9 +373,10 @@ void assembleMatrices(const FEBasis& febasis, MatrixType& stiffnessmatrix, Matri
}
-// Copy all values of a vector to the identified boundary points
+// Copy all values on the bottom/left to the identified boundary points on the top/right
template <class VectorType>
-void copyToBoundary(VectorType& vector, const int& n, const int& elementnumber, const int& feorder)
+void copyToBoundary(VectorType& vector, const int& n, const int& elementnumber,
+ const int& feorder)
{
for(int i = 0; i < n; i++){
if (isOnBoundary(i, elementnumber, feorder)){
@@ -383,19 +388,20 @@ void copyToBoundary(VectorType& vector, const int& n, const int& elementnumber,
}
-// Interpolate the surface function in the grid view
+// Interpolate the surface function on a given 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,
- std::vector<int>& k1, VectorType& xk, VectorType& yk, double& pi)
+void interpolateSurfaceFunction(const FEBasis& febasis, GridView& gridview,
+ VectorType& surface, std::vector<char>& nonboundarynodes, const int& elementnumber,
+ const int& feorder, int& modenumber, std::vector<int>& k0, std::vector<int>& k1,
+ VectorType& xk, VectorType& yk, double& pi)
{
const int n = febasis.size();
- auto graphFunction = [c, k0, k1, xk, yk, pi](const auto x) {
+ auto graphFunction = [modenumber, k0, k1, xk, yk, pi](const auto x) {
double value = 0;
double sp;
- for(int i=0; i < c; i++){
+ for(int i=0; i < modenumber; i++){
sp = 2*pi*(x[0]*k0[i] + x[1]*k1[i]);
value += std::sin(sp)*xk[i] - std::cos(sp)*yk[i];
}
@@ -452,19 +458,19 @@ int main (int argc, char *argv[]) try
// Setting up the random surface
- int c;
+ int modeNumber;
std::vector<int> k0;
std::vector<int> k1;
- indexSet(eta, c, k0, k1);
+ indexSet(eta, modeNumber, k0, k1);
- std::cout << "Number of Fourier modes is " << 2*c << std::endl;
+ std::cout << "Number of Fourier modes is " << 2*modeNumber << std::endl;
VectorType Xk;
VectorType Yk;
if (answer1 == "y"){
- generateSurfaceFunction(c, k0, k1, Xk, Yk, kappa, sigma);
+ generateSurfaceFunction(modeNumber, k0, k1, Xk, Yk, kappa, sigma);
if (answer2 == "y"){
storeMatrixMarket(Xk, "X" + name1);
@@ -499,7 +505,7 @@ int main (int argc, char *argv[]) try
if (answer3 == "y"){
timer.reset();
- assembleMatrices(feBasis, stiffnessMatrix, massMatrix, elementNumber, feOrder, c, k0, k1, Xk, Yk, pi);
+ assembleMatrices(feBasis, stiffnessMatrix, massMatrix, elementNumber, feOrder, modeNumber, k0, k1, Xk, Yk, pi);
std::cout << "Assembling the matrices took " << timer.elapsed() << "s" << std::endl;
if (answer4 == "y"){
@@ -531,7 +537,7 @@ int main (int argc, char *argv[]) try
VectorType surface;
if (answer5 == "y"){
timer.reset();
- interpolateSurfaceFunction(feBasis, gridView, surface, nonBoundaryNodes, elementNumber, feOrder, c, k0, k1, Xk, Yk, pi);
+ interpolateSurfaceFunction(feBasis, gridView, surface, nonBoundaryNodes, elementNumber, feOrder, modeNumber, k0, k1, Xk, Yk, pi);
std::cout << "Interpolating the surface function took " << timer.elapsed() << "s" << std::endl;
}