Forked from
agnumpde / dune-fufem
362 commits behind the upstream repository.
-
Jonathan Youett authoredJonathan Youett authored
Code owners
Assign users and groups as approvers for specific file changes. Learn more.
boundarypatch.hh 31.25 KiB
// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set ts=8 sw=4 et sts=4:
#ifndef BOUNDARY_PATCH_HH
#define BOUNDARY_PATCH_HH
#include <vector>
#include <set>
#include <memory>
#include <dune/common/bitsetvector.hh>
#include <dune/common/version.hh>
#include <dune/geometry/referenceelements.hh>
#include <dune/grid/common/mcmgmapper.hh>
#include <dune/grid/common/gridview.hh>
#include <dune/grid/common/gridfactory.hh>
#include <dune/fufem/globalintersectioniterator.hh>
#include <dune/fufem/boundaryiterator.hh>
// Forward declaration
template <class GridView> class BoundaryPatchIterator;
/** \brief Insertion property for segment insertion indices
*
* You can use a BoundaryPatchInsertionIndexProperty
* with BoundaryPatch::insertFacesByProperty()
* to insert all faces contained a set boundary
* segments described by insertion indices.
*/
template <class GridView>
class BoundaryPatchInsertionIndexProperty
{
typedef typename Dune::GridFactory<typename GridView::Grid> GridFactory;
public:
typedef typename GridView::Intersection Intersection;
/** \brief Create property from GridFactory and container of insertion indices
*/
template<class T>
BoundaryPatchInsertionIndexProperty(const GridFactory& factory, const T& t) : factory_(factory)
{
indices_.insert(t.begin(), t.end());
}
/** \brief Create property from GridFactory and a single insertion index
*/
BoundaryPatchInsertionIndexProperty(const GridFactory& factory, int i) : factory_(factory)
{
indices_.insert(i);
}
/** \brief Check if intersection is contained in one of the boundary segments
*/
bool operator() (const Intersection& i) const
{
return (factory_.wasInserted(i) and (indices_.find(factory_.insertionIndex(i))!=indices_.end()));
}
private:
const GridFactory& factory_;
std::set<int> indices_;
};
/** \brief Insertion property for segment indices
*
* You can use a BoundaryPatchSegmentIndexProperty
* with BoundaryPatch::insertFacesByProperty()
* to insert all faces contained a set boundary
* segments described by boundary segment indices.
*/
template <class GridView>
class BoundaryPatchSegmentIndexProperty
{
public:
typedef typename GridView::Intersection Intersection;
/** \brief Create property from a container of boundary segment indices
*/
template<class T>
BoundaryPatchSegmentIndexProperty(const T& t)
{
indices_.insert(t.begin(), t.end());
}
/** \brief Create property from a single boundary segment index
*/
BoundaryPatchSegmentIndexProperty(int i)
{
indices_.insert(i);
}
/** \brief Check if intersection is contained in one of the boundary segments
*/
bool operator() (const Intersection& i) const
{
return (indices_.find(i.boundarySegmentIndex())!=indices_.end());
}
private:
std::set<int> indices_;
};
/** \brief Insertion property for vertex vector
*
* You can use a BoundaryPatchEnclosingVerticesProperty
* with BoundaryPatch::insertFacesByProperty()
* to insert all faces that are enclosed by vertices
* marked in a vector.
*/
template <class GridView, int ncomp = 1>
class BoundaryPatchEnclosingVerticesProperty
{
typedef typename GridView::IndexSet IndexSet;
static const int dim = GridView::dimension;
public:
typedef typename GridView::Intersection Intersection;
/** \brief Create property from a marker vector
*/
BoundaryPatchEnclosingVerticesProperty(const GridView& gridView, const Dune::BitSetVector<ncomp>& vertices) :
indexSet_(gridView.indexSet()),
vertices_(vertices)
{}
/** \brief Check if intersection is enclosed by vertices in the vector
*/
bool operator() (const Intersection& i) const
{
const auto inside = i.inside();
int localFaceIndex = i.indexInInside();
const auto& refElement = Dune::ReferenceElements<double, dim>::general(inside.type());
// Get global node ids
int n = refElement.size(localFaceIndex, 1, dim);
// Using ReferenceElement::subEntity is OK here, because we loop
// over all subEntities (i.e. sub-vertices) and just return false
// if _any_ of them is not marked.
for (int i=0; i<n; i++) {
int localVertexIndex = refElement.subEntity(localFaceIndex, 1, i, dim);
if (not(vertices_[indexSet_.subIndex(inside, localVertexIndex, dim)].any()))
return false;
}
return true;
}
private:
const IndexSet& indexSet_;
const Dune::BitSetVector<ncomp>& vertices_;
};
/** \brief Insertion property for vertex vector
*
* You can use a BoundaryPatchTouchingVerticesProperty
* with BoundaryPatch::insertFacesByProperty()
* to insert all faces that are touched by any vertices
* marked in a vector.
*/
template <class GridView>
class BoundaryPatchTouchingVerticesProperty
{
typedef typename GridView::IndexSet IndexSet;
static const int dim = GridView::dimension;
public:
typedef typename GridView::Intersection Intersection;
/** \brief Create property from a marker vector
*/
BoundaryPatchTouchingVerticesProperty(const GridView& gridView, const Dune::BitSetVector<1>& vertices) :
indexSet_(gridView.indexSet()),
vertices_(vertices)
{}
/** \brief Check if intersection is enclosed by vertices in the vector
*/
bool operator() (const Intersection& i) const
{
const auto inside = i.inside();
int localFaceIndex = i.indexInInside();
const auto& refElement = Dune::ReferenceElements<double, dim>::general(inside.type());
// Get global node ids
int n = refElement.size(localFaceIndex, 1, dim);
// Using ReferenceElement::subEntity is OK here, because we loop
// over all subEntities (i.e. sub-vertices) and return true as soon
// as _any_ of them is marked.
for (int i=0; i<n; i++) {
int localVertexIndex = refElement.subEntity(localFaceIndex, 1, i, dim);
if (vertices_[indexSet_.subIndex(inside, localVertexIndex, dim)][0])
return true;
}
return false;
}
private:
const IndexSet& indexSet_;
const Dune::BitSetVector<1>& vertices_;
};
/** \brief Encapsulate a part of a grid boundary
\tparam GridView The grid view on which this boundary patch lives
*/
template <class GV>
class BoundaryPatch
{
public:
typedef GV GridView;
protected:
static const int dim = GridView::dimension;
typedef typename GridView::IndexSet IndexSetType;
typedef typename GridView::IntersectionIterator IntersectionIterator;
typedef typename GridView::Intersection Intersection;
typedef typename GridView::template Codim<0>::Entity EntityType;
typedef typename GridView::template Codim<0>::Iterator ElementIterator;
typedef typename Dune::MultipleCodimMultipleGeomTypeMapper<GridView> MapperType;
public:
/** \brief Export the appropriate iterator type */
typedef BoundaryPatchIterator<GridView> iterator;
BoundaryPatch() : verticesAreUpToDate_(false)
{}
/** \brief Constructor for a given grid view
\param initialValue If true: BoundaryPatch contains all boundary faces
*/
BoundaryPatch(const GridView& gridView, bool initialValue = false)
: verticesAreUpToDate_(false)
{
setup(gridView, initialValue);
}
BoundaryPatch(const BoundaryPatch& other) {
vertices_ = other.vertices_;
verticesAreUpToDate_ = other.verticesAreUpToDate_;
faces_ = other.faces_;
if (other.gridView_.get() != NULL)
gridView_ = std::unique_ptr<GridView>(new GridView(*other.gridView_));
if (other.mapper_.get() != NULL)
mapper_ = std::unique_ptr<MapperType>(new MapperType(*gridView_,Dune::mcmgLayout(Dune::Codim<1>())));
}
/** \brief Construct boundary patch from given vertices
\param gridView The grid view
\param vertices A vector of which any bit component is set if the corresponding vertex belongs to the patch.
*/
template <int blocksize>
BoundaryPatch(const GridView& gridView, const Dune::BitSetVector<blocksize>& vertices)
{
setup(gridView, vertices);
}
/** \brief Setup for a given grid view
*
* \param initialValue If true: BoundaryPatch contains all boundary faces
*/ void setup(const GridView& gridView, bool initialValue = false)
{
gridView_ = std::unique_ptr<GridView>(new GridView(gridView));
mapper_ = std::unique_ptr<MapperType>(new MapperType(*gridView_,Dune::mcmgLayout(Dune::Codim<1>())));
faces_.resize(gridView_->size(1));
faces_.unsetAll();
verticesAreUpToDate_ = false;
// include all boundary faces, if requested
if (initialValue)
{
BoundaryIterator<GridView> bIt(*gridView_, BoundaryIterator<GridView>::begin);
BoundaryIterator<GridView> bEnd(*gridView_, BoundaryIterator<GridView>::end);
for(; bIt!=bEnd; ++bIt)
addFace(bIt->inside(), bIt->indexInInside());
}
}
/** \brief Setup for a given grid view and a bitfield
*
* \param gridView The grid view
* \param vertices A vector of which any bit component is set if the corresponding vertex belongs to the patch.
*/
template <int blocksize>
void setup(const GridView& gridView, const Dune::BitSetVector<blocksize>& vertices)
{
// setup empty BoundaryPatch
setup(gridView);
// create property for insertion by vertex vector
BoundaryPatchEnclosingVerticesProperty<GridView, blocksize> prop(*gridView_, vertices);
// insert by property
insertFacesByProperty(prop);
}
/** \brief Return iterator on first boundary face */
iterator begin() const
{
return iterator(gridView_->template begin<0>(),
gridView_->template end<0>(), *this);
}
/** \brief Return iterator after last boundary patch face */
iterator end() const
{
return iterator(gridView_->template end<0>(),
gridView_->template end<0>(), *this);
}
/** \brief Clear the boundary patch, but don't disconnect it from the grid and level */
void clear() {
vertices_.unsetAll();
verticesAreUpToDate_ = false;
faces_.unsetAll();
}
void getFaces(Dune::BitSetVector<1>& faces) const {
faces = faces_;
}
/** \brief Returns the number of faces in the surface */
int numFaces() const {
return faces_.count();
}
/** \brief Add a boundary face to the patch */
void addFace(const IntersectionIterator& nIt) {
this->addFace(*nIt);
}
/** \brief Add a boundary face to the patch */
void addFace(const Intersection& nIt) {
this->faces_[mapper_->subIndex(nIt.inside(), nIt.indexInInside(),1)] = true;
}
/** \brief Add a boundary face to the patch */
void addFace(const EntityType& en, int fIdx) {
this->faces_[mapper_->subIndex(en, fIdx,1)] = true;
}
/** \brief Insert boundary faces filtered by an insertion property
*
* \tparam InsertionProperty a class providing bool operator()(Intersection)
* \param property All boundary intersections with property(intersection)==true are inserted
*/
template<class InsertionProperty>
void insertFacesByProperty(const InsertionProperty& property)
{
BoundaryIterator<GridView> bIt(*gridView_, BoundaryIterator<GridView>::begin);
BoundaryIterator<GridView> bEnd(*gridView_, BoundaryIterator<GridView>::end);
for(; bIt!=bEnd; ++bIt)
{
if (property(*bIt))
addFace(bIt->inside(), bIt->indexInInside());
}
this->verticesAreUpToDate_ = false;
}
/** \brief Return true if the BoundaryPatch contains a certain boundary face */
template<class IntersectionType>
bool contains(const IntersectionType& intersection) const {
return this->faces_[mapper_->subIndex(intersection.inside(), intersection.indexInInside(),1)][0];
}
/** \brief Return true if the BoundaryPatch contains a certain boundary face */
bool contains(const IntersectionIterator& nIt) const {
return this->faces_[mapper_->subIndex(nIt->inside(), nIt->indexInInside(),1)][0];
}
/** \brief Return true if the BoundaryPatch contains a certain boundary face */
bool contains(const EntityType& en, int faceIdx) const {
return this->faces_[mapper_->subIndex(en, faceIdx,1)][0];
}
/** \brief Return true if the BoundaryPatch contains a certain subentity in some face */
bool contains(const EntityType& en, int subEntity, int codim) const {
if (codim==0) // disregard element dofs
return false;
if (codim==1)
return this->faces_[mapper_->subIndex(en, subEntity, 1)][0];
const auto& refElement = Dune::ReferenceElements<double, GridView::dimension>::general(en.type());
for (int faceIdx = 0; faceIdx<refElement.size(1); ++faceIdx)
{
if (this->faces_[mapper_->subIndex(en, faceIdx, 1)][0])
if (isSubEntity(refElement, faceIdx, 1, subEntity, codim))
return true;
}
return false;
}
/** \brief Return true if the BoundaryPatch contains a face of given element
*
*/
bool containsFaceOf(const EntityType& e) const
{
int numFaces = e.subEntities(1);
for (int faceIdx = 0; faceIdx<numFaces; ++faceIdx)
{
if (contains(e, faceIdx))
return true;
}
return false;
}
/** \brief Return reference to the carrier grid view */
const GridView& gridView() const {return *gridView_;}
/** \brief Check if BoundaryPatch is initialized */
bool isInitialized() const
{
return gridView_.get();
}
/** \brief Partial deep copy */
BoundaryPatch& operator=(const BoundaryPatch& other) {
vertices_ = other.vertices_;
verticesAreUpToDate_ = other.verticesAreUpToDate_;
faces_ = other.faces_;
gridView_ = std::unique_ptr<GridView>(new GridView(*other.gridView_));
mapper_ = std::unique_ptr<MapperType>(new MapperType(*gridView_,Dune::mcmgLayout(Dune::Codim<1>())));
return *this;
}
bool operator==(const BoundaryPatch& other) {
size_t s = faces_.size();
if (s != other.faces_.size())
return false;
// There is no BitSetVector::operator==
for (size_t i = 0; i < s; ++i)
if (faces_[i] != other.faces_[i])
return false;
return true;
}
bool containsVertex(size_t v) const {
if (!verticesAreUpToDate_)
makeVertices();
#ifdef DUNE_RANGE_CHECKING
if (v<0 || v>=vertices_.size())
DUNE_THROW(Dune::RangeError, "Trying to check for vertex " << v
<< ", but the vertices are only in the range [0..." << int(vertices_.size())-1
<< "]!");
#endif
return vertices_[v][0];
}
/** \brief Returns the area of the boundary patch */
typename GridView::Grid::ctype area() const {
typename GridView::Grid::ctype patchArea = 0;
// loop over all elements
iterator it = begin();
iterator endIt = end();
for (; it!=endIt; ++it)
patchArea += it->geometry().volume();
return patchArea;
}
void getVertices(Dune::BitSetVector<1>& vertices) const {
if (!this->verticesAreUpToDate_) this->makeVertices();
vertices = this->vertices_;
}
const Dune::BitSetVector<1>* getVertices() const {
if (!this->verticesAreUpToDate_)
this->makeVertices();
return &this->vertices_;
}
/**
* \brief Get the vertices which are on the boundary of the boundary patch
*
* ATTENTION: Don't use this method, until you're absolutely sure that the internals
* work correctly for your grid implementation! It does not work for general
* grids. See the internal documentation for details.
*/
void getPatchBoundaryVertices(Dune::BitSetVector<1>& patchBoundary) const {
const IndexSetType& indexSet = gridView_->indexSet();
patchBoundary.resize(indexSet.size(dim));
patchBoundary.unsetAll();
// loop over all elements
for (const auto& e : elements(*gridView_)) {
const auto& refElement = Dune::ReferenceElements<double, dim>::general(e.type());
// Loop over all neighbors
for (const auto& is : intersections(*gridView_,e)) {
// if the element is a boundary element, but _not_ a patch element
if (is.boundary() && !contains(is)) {
// Get global node ids
int n = refElement.size(is.indexInInside(),1,dim);
// ATTENTION: This implementation is broken!
//
// It will not work in general, because there's no guarantee that:
//
// subIndex(e, refElement.subEntity(face,1,i,dim)) == subIndex(e.subEntity<1>(face), i, dim)
//
// The latter is what you really want because it uses the correct embedding
// of the face into the element (which is not the same as the one in the
// reference element)
for (int i=0; i<n; i++) {
int faceIdxi = refElement.subEntity(is.indexInInside(), 1, i, dim);
int globalIdx = indexSet.subIndex(e, faceIdxi,dim);
if (this->containsVertex(globalIdx))
patchBoundary[globalIdx] = true;
}
}
}
}
}
/**
* \brief Get the edges which are on the boundary of the boundary patch
*
* \todo This only works when dim==3
*
* ATTENTION: Don't use this method, until you're absolutely sure that the internals * work correctly for your grid implementation! It does not work for general
* grids. See the internal documentation for details.
*/
void getPatchBoundaryEdges(Dune::BitSetVector<1>& patchBoundaryEdges) const {
assert(dim==3);
const IndexSetType& indexSet = gridView_->indexSet();
patchBoundaryEdges.resize(indexSet.size(dim-1));
patchBoundaryEdges.unsetAll();
// loop over all elements
for (const auto& e : elements(*gridView_)) {
const auto& refElement = Dune::ReferenceElements<double, dim>::general(e.type());
// Loop over all neighbors
for (const auto& is : intersections(*gridView_,e)) {
// if the element is a boundary element, but _not_ a patch element
if (is.boundary() && !contains(is)) {
// Loop over the face edges
int n = refElement.size(is.indexInInside(),1,dim-1);
for (int i=0; i<n; i++) {
// ATTENTION: This implementation is broken!
//
// It will not work in general, because there's no guarantee that:
//
// subIndex(e, refElement.subEntity(face,1,i,dim-1)) == subIndex(e.subEntity<1>(face), i, dim-1)
//
// The latter is what you really want because it uses the correct embedding
// of the face into the element (which is not the same as the one in the
// reference element).
int edge = refElement.subEntity(is.indexInInside(), 1, i, dim-1);
// If 'edge' would be correct, the following use of ReferenceElement::subEntity
// would be OK, because we just do a check for all vertices we pass, i.e. 0 and 1.
int v0 = refElement.subEntity(edge, dim-1, 0, dim);
int v1 = refElement.subEntity(edge, dim-1, 1, dim);
if (containsVertex(indexSet.subIndex(e, v0,dim))
&& containsVertex(indexSet.subIndex(e, v1,dim)))
patchBoundaryEdges[indexSet.subIndex(e, edge,dim-1)] = true;
}
}
}
}
}
/** \brief Returns the number of vertices in the surface */
int numVertices() const {
if (!this->verticesAreUpToDate_)
this->makeVertices();
return this->vertices_.count();
}
/**
* \brief Return vertex unit normals
*
* ATTENTION: This will only work correctly if unitOuterNormal
* is constant along intersections. See the internal documentation for details.
* */ void getNormals(std::vector<Dune::FieldVector<double, dim> >& normals) const {
normals = getNormals();
}
/** \brief Return vertex unit normals */
std::vector<Dune::FieldVector<double,dim> > getNormals() const {
if (!this->verticesAreUpToDate_)
this->makeVertices();
typedef Dune::FieldVector<double, dim> Vector;
std::vector<Vector> normals(this->vertices_.size(), Vector(0.0));
// Loop over all boundary segments
iterator it = begin();
iterator endIt = end();
for (; it!=endIt; ++it) {
const auto inside = it->inside();
const auto& refElement = Dune::ReferenceElements<double, dim>::general(inside.type());
// Compute actual surface normals
int n = refElement.size(it->indexInInside(),1,dim);
for (int i=0; i<n; i++) {
// ATTENTION: This implementation does only work if unitOuterNormal is constant on the intersection.
// It will not work in general, because there's no guarantee that:
//
// subIndex(e, refElement.subEntity(inside,1,i,dim)) == subIndex(e.subEntity<1>(inside), i, dim)
//
// The latter is what you really want because it uses the correct embedding
// of the face into the element (which is not the same as the one in the
// reference element).
//
// If unitOuterNormal is constant we add the same faceNormal to all vertices
// of this face, hence the proper order does not matter. Otherwise this
// implementation may compute the faceNormal at the wrong posInBoundary
// because posInElement and vertexIdx don't correspond to each other.
int faceIdxi = refElement.subEntity(it->indexInInside(), 1, i, dim);
// local position of vertex in element
const Dune::FieldVector<double, dim>& posInElement = refElement.position(faceIdxi, dim);
// local position in boundary segment
Dune::FieldVector<double, dim-1> posInBoundary = it->geometryInInside().local(posInElement);
// Get normal there
const Dune::FieldVector<double, dim>& faceNormal = it->unitOuterNormal(posInBoundary);
int vertexIdx = gridView_->indexSet().subIndex(inside, faceIdxi,dim);
normals[vertexIdx] += faceNormal;
}
}
// Normalize vertex normals
for (size_t i=0; i<normals.size(); i++)
if (this->vertices_[i][0])
normals[i] /= normals[i].two_norm();
return normals;
}
/** \brief Return vertex unit normals corr. to local indices */
std::vector<Dune::FieldVector<double,dim> > getLocalNormals() const {
if (!this->verticesAreUpToDate_)
this->makeVertices();
auto globalToLocal = makeGlobalToLocal();
using Vector = Dune::FieldVector<double, dim>;
std::vector<Vector> normals(this->vertices_.count(), Vector(0.0));
// Loop over all boundary segments
for (const auto& it : *this) {
const auto inside = it.inside();
const auto& refElement = Dune::ReferenceElements<double, dim>::general(inside.type());
// Compute actual surface normals
int n = refElement.size(it.indexInInside(),1,dim);
const auto& geomInInside = it.geometryInInside();
for (int i=0; i<n; i++) {
// See the Warning in the method above!
int faceIdxi = refElement.subEntity(it.indexInInside(), 1, i, dim);
// local position of vertex in element
auto& posInElement = refElement.position(faceIdxi, dim);
// local position in boundary segment
auto posInBoundary = geomInInside.local(posInElement);
// Get normal there
const auto& faceNormal = it.unitOuterNormal(posInBoundary);
int vertexIdx = gridView_->indexSet().subIndex(inside, faceIdxi,dim);
normals[globalToLocal[vertexIdx]] += faceNormal;
}
}
// Normalize vertex normals
for (size_t i=0; i<normals.size(); i++)
normals[i] /= normals[i].two_norm();
return normals;
}
void makeGlobalToLocal(std::vector<int>& globalToLocal) const {
globalToLocal = makeGlobalToLocal();
}
std::vector<int> makeGlobalToLocal() const {
if (!this->verticesAreUpToDate_)
this->makeVertices();
std::vector<int> globalToLocal(this->vertices_.size());
int idx = 0;
for (size_t i=0; i< globalToLocal.size(); i++)
globalToLocal[i] = (this->vertices_[i][0]) ? idx++ : -1;
return globalToLocal;
}
/** \brief Adds another patch in the sense that afterwards this patch
* contains all the vertices included in any one of the two patches
*
* \param patch Patch to be added*/
BoundaryPatch<GridView> addPatch(const BoundaryPatch<GridView>& patch) {
/* test whether both patches live on the same grid
* for LevelBoundaryPatches it is not tested whether both patches refer to the same level*/
if (&this->gridView().grid() != &patch.gridView().grid())
DUNE_THROW(Dune::NotImplemented, "addPatch only works for patches defined on the same grid");
Dune::BitSetVector<1> vertices, additionalVertices;
this->getVertices(vertices);
patch.getVertices(additionalVertices);
// create combined BitSetVector
for (size_t i=0; i<vertices.size(); i++)
if (additionalVertices[i][0])
vertices[i] = true;
GridView gridView = this->gridView();
this->setup(gridView, vertices);
return *this;
}
/** \brief Determines the overlap of this boundary patch and another
* by reducing this patch to the common vertices.
*
* \param patch Patch to be intersected with*/
BoundaryPatch<GridView> intersect(const BoundaryPatch<GridView>& patch) {
/* test whether both patches live on the same grid
* for LevelBoundaryPatches it is not tested whether both patches refer to the same level*/
if (&this->gridView().grid() != &patch.gridView().grid())
DUNE_THROW(Dune::NotImplemented, "addPatch only works for patches defined on the same grid");
Dune::BitSetVector<1> vertices, intersectionVertices;
this->getVertices(vertices);
patch.getVertices(intersectionVertices);
// create intersection BitSetVector
for (size_t i=0; i<vertices.size(); i++)
if (vertices[i][0] && intersectionVertices[i][0])
vertices[i][0] = true;
else
vertices[i][0] = false;
GridView gridView = this->gridView();
this->setup(gridView, vertices);
return *this;
}
protected:
// Given a reference element for an entity E this
// checks if subentity (j, jCodim) of E is a subentity of subentity (i,iCodim) of E
// e.g. if the j-th vertex of E is also a vertex of the i-th edge of E
template<class RefElement>
static bool isSubEntity(const RefElement& referenceElement, int i, int iCodim, int j, int jCodim)
{
for (int jj = 0; jj<referenceElement.size(i, iCodim, jCodim); jj++)
{
// This use of ReferenceElement::subEntity is OK, because
// we just check if any of the returned indices matches j.
if (referenceElement.subEntity(i, iCodim, jj, jCodim) == j)
return true;
}
return false;
}
void makeVertices() const {
const typename GridView::IndexSet& indexSet = gridView_->indexSet();
this->vertices_.resize(gridView_->size(dim)); this->vertices_.unsetAll();
// loop over all elements
iterator it = begin();
iterator endIt = end();
for (; it!=endIt; ++it) {
const auto inside = it->inside();
const auto& refElement = Dune::ReferenceElements<double, dim>::general(inside.type());
// Get global node ids
int n = refElement.size(it->indexInInside(),1,dim);
// This use of ReferenceElement::subEntity is OK, because
// we mark the vertices corresponding to all returned indices.
for (int i=0; i<n; i++) {
int faceIdxi = refElement.subEntity(it->indexInInside(), 1, i, dim);
vertices_[indexSet.template subIndex(inside,faceIdxi,dim)] = true;
}
}
verticesAreUpToDate_ = true;
}
mutable Dune::BitSetVector<1> vertices_;
mutable bool verticesAreUpToDate_;
Dune::BitSetVector<1> faces_;
std::unique_ptr<GridView> gridView_;
std::unique_ptr<MapperType> mapper_;
};
/** \brief Iterator of boundary intersections of grid view
*
* \tparam GridView The grid view on which this boundary patch lives
*/
template <class GridView>
class BoundaryPatchIterator
: public GlobalIntersectionIterator<GridView, BoundaryPatchIterator<GridView> >
{
// Only BoundaryPatch can call this class' constructors
friend class BoundaryPatch<GridView>;
typedef GlobalIntersectionIterator<GridView, BoundaryPatchIterator<GridView> > Base;
enum {dim=GridView::dimension};
typedef typename GridView::template Codim<0>::Iterator ElementIterator;
typedef typename GridView::template Codim<0>::Entity Element;
public:
/** \brief The type of objects we are iterating over */
typedef typename Base::Intersection Intersection;
/** \brief Return true if the given element contains no boundary patch faces (and hence can be skipped) */
bool skipElement(const Element& e)
{
return not(boundaryPatchPtr_->containsFaceOf(e));
}
bool skipIntersection(const Intersection& i)
{
return not(boundaryPatchPtr_->contains(i.inside(), i.indexInInside()));
}
protected:
/** \brief Create iterator for given grid view
*
* \param eIt Element iterator to the first element to consider
* \param endEIt Element iterator after the last element to consider
*/
BoundaryPatchIterator(const ElementIterator& eIt, const ElementIterator& endEIt,
const BoundaryPatch<GridView>& boundaryPatch) :
Base(boundaryPatch.gridView(), eIt, endEIt),
boundaryPatchPtr_(&boundaryPatch)
{
Base::initialIncrement();
}
const BoundaryPatch<GridView>* boundaryPatchPtr_;
};
#endif