Skip to content
Snippets Groups Projects
Commit 539eca82 authored by Hörl, Maximilian's avatar Hörl, Maximilian
Browse files

remove debug output

parent 83d2b223
No related branches found
No related tags found
No related merge requests found
......@@ -40,8 +40,6 @@ public:
//NOTE: what would be an appropiate solver here?
A.solve(d, b);
std::cout << b << std::endl;
for (int i = 0; i < dof; i++)
for (int j = 0; j < i; j++)
/*assert*/if(std::abs(A[i][j] - A[j][i]) >
......@@ -130,17 +128,6 @@ private:
//NOTE:
//const int order = 3; //order of the quadrature rule
for (const auto& elem : elements(gridView_))
{
const auto& geo = elem.geometry();
std::cout << mapper_.index(elem) << "\t";
for (int k = 0; k < geo.corners(); k++)
{
std::cout << geo.corner(k) << "\t";
}
std::cout << "\n\n";
}
//we use the basis
// phi_elem,0 (x) = indicator(elem);
// phi_elem,i (x) = x[i]*indicator(elem);
......@@ -152,13 +139,6 @@ private:
const double elemVol = geo.volume();
const auto& center = geo.center();
std::cout << "===============================\nElement " << elemIdx;
for (int k = 0; k < geo.corners(); k++)
{
std::cout << "\t" << geo.corner(k);
}
std::cout << "\n\n";
//in the system of linear equations (SLE) Ad = b,
//the index elemIdxSLE refers to the basis function phi_elem,0
//and the indices elemIdxSLE + i + 1, i = 1,...,dim, refer to the
......@@ -210,13 +190,6 @@ private:
const double intersctVol = intersctGeo.volume();
const auto& intersctCenter = intersctGeo.center();
std::cout << "+++++++++++++++++++++++++++\nIntersection ";
for (int k = 0; k < intersctGeo.corners(); k++)
{
std::cout << "\t" << intersctGeo.corner(k);
}
std::cout << "\tnormal: " << normal << "\n\n";
//TODO: quadrature rule cannot be used for dim = 1!
// const Dune::QuadratureRule<double,dim-1>& secondOrderRule =
// Dune::QuadratureRules<double,dim-1>::rule(
......@@ -227,7 +200,7 @@ private:
//quadraticIntregrals are used
Dune::FieldVector<Scalar, dim> linearIntegrals(0.0);
Dune::FieldMatrix<Scalar, dim, dim> quadraticIntregrals(0.0);
//NOTE: is there a better type for symmetric matrix?
//NOTE: is there a better type for a symmetric matrix?
for (int i = 0; i < dim; i++)
{
......@@ -239,22 +212,11 @@ private:
{
const auto& leftCorner = intersctGeo.corner(0);
const auto& rightCorner = intersctGeo.corner(1);
std::cout << "left corner: " << leftCorner << "\tright corner: " << rightCorner << "\n";
quadraticIntregrals[i][j] = intersctVol / 3 *
( leftCorner[i] * leftCorner[j] +
rightCorner[i] * rightCorner[j]
( leftCorner[i] * leftCorner[j] + rightCorner[i] * rightCorner[j]
+ 0.5 * (leftCorner[i] * rightCorner[j] +
leftCorner[j] * rightCorner[i]) );
std::cout << i << ", " << j << "\t" << leftCorner[i] * leftCorner[j] << "\t" <<
rightCorner[i] * rightCorner[j] << "\t" <<
0.5 * (leftCorner[i] * rightCorner[j] +
leftCorner[j] * rightCorner[i]) << "\t" << intersctVol / 3 *
( leftCorner[i] * leftCorner[j] +
rightCorner[i] * rightCorner[j]
+ 0.5 * (leftCorner[i] * rightCorner[j] +
leftCorner[j] * rightCorner[i]) ) << "\t" << quadraticIntregrals[i][j] << "\n\n";
/*
//use second order quadrature rule for exact evaluation of
// int_intersct x_i*x_j ds
......@@ -270,25 +232,13 @@ private:
quadraticIntregrals[j][i] = quadraticIntregrals[i][j];
}
}
std::cout << "linearIntegrals:\n" << linearIntegrals << "\n\n";
std::cout << "quadraticIntregrals:\n" << quadraticIntregrals << "\n\n";
//exact evaluation of
// int_intersct mu*jump(phi_elem,0)*jump(phi_elem,0) ds
A[elemIdxSLE][elemIdxSLE] += mu * intersctVol;
std::cout << elemIdxSLE << ", " << elemIdxSLE <<"\t" << intersctCenter<< ":\n"
<< mu * intersctVol << "\t" << A[elemIdx][elemIdx] << "\n\n";
if (intersct.neighbor()) //intersct has neighboring element
{
std::cout << "------------------------------\nNeighbor\t";
for (int k = 0; k < intersct.outside().geometry().corners(); k++)
{
std::cout << "\t" << intersct.outside().geometry().corner(k);
}
std::cout << "\n\n";
//index of the neighboring element
const int neighborIdx = mapper_.index(intersct.outside());
const int neighborIdxSLE = (dim + 1)*neighborIdx;
......@@ -302,10 +252,6 @@ private:
// = 0.5 * K * normal[i] * vol(intersct)
A[elemIdxSLE + i + 1][elemIdxSLE] +=
mu * linearIntegrals[i] - 0.5 * K * normal[i] * intersctVol;
std::cout << elemIdxSLE + i + 1 << ", " << elemIdxSLE << ":\n"
<< mu * linearIntegrals[i] << "\t" <<
- 0.5 * K * normal[i] * intersctVol << "\t"
<< A[elemIdxSLE + i + 1][elemIdxSLE] << "\n\n";
for (int j = 0; j <= i; j++)
{
......@@ -319,11 +265,6 @@ private:
+= mu * quadraticIntregrals[i][j]
- 0.5 * K * (normal[i] * linearIntegrals[j]
+ normal[j] * linearIntegrals[i]);
std::cout << elemIdxSLE + i + 1 << ", " << elemIdxSLE + j + 1 << ":\n"
<< mu * quadraticIntregrals[i][j] << "\t" <<
- 0.5 * K * (normal[i] * linearIntegrals[j]
+ normal[j] * linearIntegrals[i]) << "\t" <<
A[elemIdxSLE + i + 1][elemIdxSLE + j + 1] << "\n\n";
}
}
......@@ -351,10 +292,6 @@ private:
// = 0.5 * K * normal[i] * vol(intersct)
A[elemIdxSLE + i + 1][neighborIdxSLE] +=
-mu * linearIntegrals[i] + 0.5 * K * normal[i] * intersctVol;
std::cout << elemIdxSLE + i + 1 << ", " << neighborIdxSLE << ":\n"
<< -mu * linearIntegrals[i] << "\t" <<
0.5 * K * normal[i] * intersctVol << "\t" <<
A[elemIdxSLE + i + 1][neighborIdxSLE] << "\n\n";
//we use the relations
// int_intersct mu*jump(phi_elem,0)
......@@ -366,10 +303,6 @@ private:
// = 0.5 * K * normal[i] * vol(intersct)
A[elemIdxSLE][neighborIdxSLE + i + 1] +=
-mu * linearIntegrals[i] - 0.5 * K * normal[i] * intersctVol;
std::cout << elemIdxSLE << ", " << neighborIdxSLE + i + 1<< ":\n"
<< -mu * linearIntegrals[i] << "\t" <<
-0.5 * K * normal[i] * intersctVol << "\t" <<
A[elemIdxSLE][neighborIdxSLE + i + 1] << "\n\n";
//stiffness matrix A is symmetric
A[neighborIdxSLE][elemIdxSLE + i + 1] +=
......@@ -395,11 +328,6 @@ private:
-mu * quadraticIntregrals[i][j]
- 0.5 * K * (normal[j] * linearIntegrals[i]
- normal[i] * linearIntegrals[j]);
std::cout << elemIdxSLE + i + 1 << ", " << neighborIdxSLE + j + 1<< ":\n"
<< -mu * quadraticIntregrals[i][j] << "\t" <<
- 0.5 * K * (normal[j] * linearIntegrals[i]
- normal[i] * linearIntegrals[j]) << "\t" <<
A[elemIdxSLE + i + 1][neighborIdxSLE + j + 1] << "\n\n";
//stiffness matrix A is symmetric
A[neighborIdxSLE + j + 1][elemIdxSLE + i + 1] +=
......@@ -422,11 +350,6 @@ private:
-mu * quadraticIntregrals[i][j]
- 0.5 * K * (normal[i] * linearIntegrals[j]
- normal[j] * linearIntegrals[i]);
std::cout << elemIdxSLE + j + 1 << ", " << neighborIdxSLE + i + 1<< ":\n"
<< -mu * quadraticIntregrals[i][j] << "\t" <<
- 0.5 * K * (normal[i] * linearIntegrals[j]
- normal[j] * linearIntegrals[i]) << "\t" <<
A[elemIdxSLE + j + 1][neighborIdxSLE + i + 1] << "\n\n";
//stiffness matrix A is symmetric
A[neighborIdxSLE + i + 1][elemIdxSLE + j + 1] +=
......@@ -437,7 +360,6 @@ private:
}
else //boundary facet
{
std::cout << "------------------------------\nBoundary\n\n";
for (int i = 0; i < dim; i++)
{ //we use the relations
// int_intersct mu*jump(phi_elem,0)*jump(phi_elem,i) ds
......@@ -448,11 +370,6 @@ private:
// = K * normal[i] * vol(intersct)
A[elemIdxSLE + i + 1][elemIdxSLE] +=
mu * linearIntegrals[i] - 0.5 * K * normal[i] * intersctVol;
std::cout << elemIdxSLE + i + 1 << ", " << elemIdxSLE << "\t"
<< intersctCenter << ":\n" << mu * linearIntegrals[i] << "\t" <<
- 0.5 * K * normal[i] * intersctVol << "\t"
<< A[elemIdxSLE + i + 1][elemIdxSLE] << "\n\n";
for (int j = 0; j <= i; j++)
{
......@@ -467,11 +384,6 @@ private:
mu * quadraticIntregrals[i][j]
- 0.5 * K * (normal[i] * linearIntegrals[j]
+ normal[j] * linearIntegrals[i]);
std::cout << elemIdxSLE + i + 1 << ", " << elemIdxSLE + j+ 1 << "\t"
<< intersctCenter << ":\n" << mu * quadraticIntregrals[i][j] << "\t" <<
- 0.5 * K * (normal[i] * linearIntegrals[j]
+ normal[j] * linearIntegrals[i]) << "\t"
<< A[elemIdxSLE + i + 1][elemIdxSLE + j + 1] << "\n\n";
}
}
}
......@@ -489,15 +401,12 @@ private:
}
}
}
std::cout << A << std::endl;
}
Matrix A; //stiffness matrix
Vector b; //load vector
Vector d; //solution vector
};
#endif
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment