diff --git a/src/one-body-sample.cc b/src/one-body-sample.cc index 7c08e4aa1b060670ceadf01d0fac6646647ab10b..f00312f20b6e75dd1a9dde155c7cb6e4ecd9579f 100644 --- a/src/one-body-sample.cc +++ b/src/one-body-sample.cc @@ -25,6 +25,8 @@ #include <dune/solvers/norms/energynorm.hh> #include <dune/solvers/common/numproc.hh> // Solver::FULL +#include <dune/grid/common/mcmgmapper.hh> + #include <exception> int const dim = 2; @@ -83,7 +85,41 @@ int main() { BlockGSStep<OperatorType, VectorType> blockGSStep(stiffnessMatrix, u, f); EnergyNorm<OperatorType, VectorType> energyNorm(blockGSStep); - Dune::BitSetVector<VectorType::block_type::dimension> ignoreNodes(false); + Dune::BitSetVector<VectorType::block_type::dimension> ignoreNodes( + grid.size(grid.maxLevel(), dim), false); + + { // Play around with the boundary + typedef GridType::LeafGridView GridView; + GridView leafView = grid.leafView(); + typedef GridView::Codim<dim>::Iterator VertexLeafIterator; + + typedef Dune::MultipleCodimMultipleGeomTypeMapper< + GridView, Dune::MCMGVertexLayout> VertexMapper; + VertexMapper myVertexMapper(leafView); + + size_t bounding_nodes = 0; + for (VertexLeafIterator it = leafView.begin<dim>(); + it != leafView.end<dim>(); ++it) { + Dune::GeometryType gt = it->type(); + assert(it->geometry().corners() == 1); + SmallVector coordinates = it->geometry().corner(0); + bool bounding(false); + for (int i = 0; i < dim; ++i) { + if (coordinates[i] == end_points[i] || coordinates[i] == 0) { + bounding = true; + break; + } + } + if (bounding) { + ++bounding_nodes; + size_t id = myVertexMapper.map(*it); + std::cout << "Ignoring id #" << id << std::endl; + ignoreNodes[id] = true; + } + } + std::cout << "Number of bounding nodes: " << bounding_nodes << std::endl; + } + blockGSStep.ignoreNodes_ = &ignoreNodes; LoopSolver<VectorType> solver(&blockGSStep, solver_maxIterations,