diff --git a/dune/mmdg/dg.hh b/dune/mmdg/dg.hh index 73835628eec65a9616fc96423d9c4b57f6b56b9e..7a0043a3a18c44cf13d124a5bf977e241e2ae12b 100644 --- a/dune/mmdg/dg.hh +++ b/dune/mmdg/dg.hh @@ -145,6 +145,8 @@ public: //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 { @@ -161,6 +163,10 @@ public: // = 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++) { @@ -174,6 +180,11 @@ public: += 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"; } } @@ -201,6 +212,10 @@ public: // = 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) @@ -212,6 +227,10 @@ public: // = 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] += @@ -237,6 +256,11 @@ public: -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] += @@ -259,6 +283,11 @@ public: -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] += @@ -279,6 +308,11 @@ public: // = 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++) { @@ -293,6 +327,11 @@ public: 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"; } } } @@ -316,17 +355,42 @@ 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]) > std::numeric_limits<Scalar>::epsilon()) std::cout << i << ", " << j << std::endl; + //NOTE: dim = 1 + auto exactSolution = [](const auto& pos) {return 0.5*pos[0]*(pos[0] - 1);}; + + //NOTE: dim = 2 +/* auto exactSolution = [](const auto& pos) + { + double solution = 0.0; + for (int m = 1; m <= 21; m = m + 2) + { + for (int n = 1; n <= 21; n = n + 2) + { + solution += sin(M_PI * m * pos[0]) * sin(M_PI * n * pos[1]) + / (m * n * (m*m + n*n)); + } + } + solution *= -16 / pow(M_PI, 4); + return solution; + }; +*/ + + + //storage for pressure data, for each element we store the //pressure at the corners of the element, the VTKFunction will //create the output using a linear interpolation for each element - Dune::DynamicVector<Scalar> - pressure((dim + 1)*gridView_.size(0), 0.0); + Dune::DynamicVector<Scalar> pressure((dim + 1)*gridView_.size(0), 0.0); + Dune::DynamicVector<double> + exactPressure((dim + 1)*gridView_.size(0), 0.0); for (const auto& elem : elements(gridView_)) { @@ -335,6 +399,8 @@ public: for (int k = 0; k < geo.corners(); k++) { + exactPressure[elemIdxSLE + k] = exactSolution(geo.corner(k)); + //contribution of the basis function // phi_elem,0 (x) = indicator(elem); //at the kth corner of elem @@ -350,71 +416,20 @@ public: } } - Dune::DynamicVector<Scalar> - pressureTEST(gridView_.size(0), 0.0); - - for (const auto& elem : elements(gridView_)) - { - const int elemIdx = mapper_.index(elem); - const int elemIdxSLE = (dim + 1)*elemIdx; - const auto& center = elem.geometry().center(); - pressureTEST[elemIdx] = d[elemIdxSLE]; - - for (int i = 0; i < dim; i++) - { - pressureTEST[elemIdx] += d[elemIdxSLE + i + 1] * center[i]; - } - } - - - //NOTE: dim = 1 -// auto exactSolution = [](const auto& pos) {return 0.5*pos[0]*(pos[0] - 1);}; - - //NOTE: dim = 2 - auto exactSolution = [](const auto& pos) - { - double solution = 0.0; - for (int m = 1; m <= 21; m = m + 2) - { - for (int n = 1; n <= 21; n = n + 2) - { - solution += sin(M_PI * m * pos[0]) * sin(M_PI * n * pos[1]) - / (m * n * (m*m + n*n)); - } - } - solution *= -16 / pow(M_PI, 4); - return solution; - }; - - - Dune::DynamicVector<double> exactPressure(gridView_.size(0)); - for (const auto& elem : elements(gridView_)) - { - exactPressure[mapper_.index(elem)] = - exactSolution(elem.geometry().center()); - } - - //create output directory if necessary mkdir("data", S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH); Dune::VTKWriter<GridView> vtkWriter(gridView_, Dune::VTK::nonconforming); - // std::make_shared<Dune::VTKWriter<GridView>>(gridView_); -// Dune::VTKSequenceWriter<GridView> -// vtkSqWriter(vtkWriter, "pressure", "data", "data"); using VTKFunction = Dune::VTKFunction<GridView>; using P1Function = Dune::NonconformingP1VTKFunction<GridView, Dune::DynamicVector<double>>; - VTKFunction* vtkOperator = - new P1Function(gridView_, pressure, "pressure"); - - //NOTE: why does this only work with a sequence writer? - vtkWriter.addVertexData( - std::shared_ptr<const VTKFunction>(vtkOperator) ); - vtkWriter.addCellData(exactPressure, "exactPressure"); - vtkWriter.addCellData(pressureTEST, "pressureTEST"); - vtkWriter.write("linearpressure"); + + vtkWriter.addVertexData( std::shared_ptr<const VTKFunction>( + new P1Function(gridView_, pressure, "pressure"))); + vtkWriter.addVertexData( std::shared_ptr<const VTKFunction>( + new P1Function(gridView_, exactPressure, "exactPressure"))); + vtkWriter.pwrite("pressureData", "", "data"); } const GridView& gridView_;