From 726b6952c863d92169e68a225632c48357921eb8 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Maximilian=20H=C3=B6rl?=
<maximilian.hoerl@mathematik.uni-stuttgart.de>
Date: Fri, 17 Jan 2020 17:59:04 +0100
Subject: [PATCH] [WIP] add debug output
---
dune/mmdg/dg.hh | 53 ++++++++++++++++++++++++++++++++++++++++++++++---
1 file changed, 50 insertions(+), 3 deletions(-)
diff --git a/dune/mmdg/dg.hh b/dune/mmdg/dg.hh
index 7a0043a..80d2b02 100644
--- a/dune/mmdg/dg.hh
+++ b/dune/mmdg/dg.hh
@@ -36,6 +36,17 @@ public:
//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);
@@ -48,6 +59,13 @@ public:
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
@@ -99,6 +117,13 @@ public:
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 << "\n\n";
+
//TODO: quadrature rule cannot be used for dim = 1!
// const Dune::QuadratureRule<double,dim-1>& secondOrderRule =
// Dune::QuadratureRules<double,dim-1>::rule(
@@ -121,11 +146,22 @@ public:
{
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]
+ 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
@@ -141,6 +177,8 @@ public:
quadraticIntregrals[i][j] = quadraticIntregrals[j][i];
}
}
+ 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
@@ -150,6 +188,14 @@ public:
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;
@@ -298,6 +344,7 @@ public:
}
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
@@ -364,10 +411,10 @@ public:
std::cout << i << ", " << j << std::endl;
//NOTE: dim = 1
- auto exactSolution = [](const auto& pos) {return 0.5*pos[0]*(pos[0] - 1);};
+// auto exactSolution = [](const auto& pos) {return 0.5*pos[0]*(pos[0] - 1);};
//NOTE: dim = 2
-/* auto exactSolution = [](const auto& pos)
+ auto exactSolution = [](const auto& pos)
{
double solution = 0.0;
for (int m = 1; m <= 21; m = m + 2)
@@ -381,7 +428,7 @@ public:
solution *= -16 / pow(M_PI, 4);
return solution;
};
-*/
+
--
GitLab