Skip to content
Snippets Groups Projects
Commit 411da0d6 authored by Elias Pipping's avatar Elias Pipping Committed by Elias Pipping
Browse files

Inline setupBoundary

parent 06f567ac
No related branches found
No related tags found
No related merge requests found
* Setup
#+name: assembleMassMatrix
#+begin_src c++
{
......@@ -54,6 +55,58 @@
}
#+end_src c++
#+name: setupBoundary
#+begin_src c++
{
typedef typename GridView::template Codim<dim>::Iterator VertexLeafIterator;
Dune::MultipleCodimMultipleGeomTypeMapper
<GridView, Dune::MCMGVertexLayout> const myVertexMapper(leafView);
for (auto it = leafView.template begin<dim>();
it != leafView.template end<dim>();
++it) {
assert(it->geometry().corners() == 1);
Dune::FieldVector<double, dim> const coordinates = it->geometry().corner(0);
size_t const id = myVertexMapper.map(*it);
// Find the center of the lower face
switch (dim) {
case 3:
if (coordinates[2] != lowerLeft[2])
break;
// fallthrough
case 2:
if (coordinates[1] == lowerLeft[1]
&& std::abs(coordinates[0] - lowerLeft[0]) < 1e-8)
specialNode = id;
break;
default:
assert(false);
}
// lower face
if (coordinates[1] == lowerLeft[1]) {
frictionalNodes[id] = true;
ignoreNodes[id][1] = true;
}
// upper face
else if (coordinates[1] == upperRight[1])
ignoreNodes[id] = true;
// right face (except for both corners)
else if (coordinates[0] == upperRight[0])
;
// left face (except for both corners)
else if (coordinates[0] == lowerLeft[0])
;
}
}
#+end_src
* Main
#+begin_src c++ :tangle one-body-sample.cc :noweb yes
#ifdef HAVE_CONFIG_H
#include "config.h"
......@@ -137,63 +190,6 @@
int const dim = DIM;
template <class GridView, class GridCorner>
void
setup_boundary(GridView const &gridView,
Dune::BitSetVector<dim> &ignoreNodes,
Dune::BitSetVector<1> &neumannNodes,
Dune::BitSetVector<1> &frictionalNodes,
GridCorner const &lowerLeft,
GridCorner const &upperRight,
size_t &specialNode)
{
typedef typename GridView::template Codim<dim>::Iterator VertexLeafIterator;
Dune::MultipleCodimMultipleGeomTypeMapper
<GridView, Dune::MCMGVertexLayout> const myVertexMapper(gridView);
for (auto it = gridView.template begin<dim>();
it != gridView.template end<dim>();
++it) {
assert(it->geometry().corners() == 1);
Dune::FieldVector<double, dim> const coordinates = it->geometry().corner(0);
size_t const id = myVertexMapper.map(*it);
// Find the center of the lower face
switch (dim) {
case 3:
if (coordinates[2] != lowerLeft[2])
break;
// fallthrough
case 2:
if (coordinates[1] == lowerLeft[1]
&& std::abs(coordinates[0] - lowerLeft[0]) < 1e-8)
specialNode = id;
break;
default:
assert(false);
}
// lower face
if (coordinates[1] == lowerLeft[1]) {
frictionalNodes[id] = true;
ignoreNodes[id][1] = true;
}
// upper face
else if (coordinates[1] == upperRight[1])
ignoreNodes[id] = true;
// right face (except for both corners)
else if (coordinates[0] == upperRight[0])
;
// left face (except for both corners)
else if (coordinates[0] == lowerLeft[0])
;
}
}
template<class VectorType, class MatrixType, class FunctionType, int dim>
Dune::shared_ptr
<TimeSteppingScheme
......@@ -339,8 +335,7 @@
Dune::BitSetVector<dim> ignoreNodes(finestSize, false);
Dune::BitSetVector<1> neumannNodes(finestSize, false);
Dune::BitSetVector<1> frictionalNodes(finestSize, false);
setup_boundary(leafView, ignoreNodes, neumannNodes, frictionalNodes,
lowerLeft, upperRight, specialNode);
<<setupBoundary>>
assert(specialNode != finestSize);
assert(frictionalNodes[specialNode][0]);
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment