diff --git a/src/spatial-solving/preconditioners/supportpatchfactory.hh b/src/spatial-solving/preconditioners/supportpatchfactory.hh
index f518ecc0de92b9fcdbaab92bdfca58b36f26072a..76d64e3761a95ae0196c94e0a11e46bedb856b0f 100644
--- a/src/spatial-solving/preconditioners/supportpatchfactory.hh
+++ b/src/spatial-solving/preconditioners/supportpatchfactory.hh
@@ -16,128 +16,77 @@
//#include "../../data-structures/levelcontactnetwork.hh"
#include "hierarchicleveliterator.hh"
-template <class Element>
-class BodyElement {
-private:
- size_t bodyID_;
- Element element_;
-
-public:
- BodyElement() {}
- BodyElement(size_t bodyID, const Element& element) {
- set(bodyID, element);
- }
-
- void set(size_t bodyID, const Element& element) {
- bodyID_ = bodyID;
- element_ = element;
- }
-
- const auto& bodyID() const {
- return bodyID_;
- }
- const auto& element() const {
- return element_;
- }
-};
-
-
template <class LevelContactNetwork>
-class BodyHelper {
+class NetworkIndexMapper {
private:
+ using FiniteElementCache = Dune::PQkLocalFiniteElementCache<ctype, ctype, dim, 1>; // P1 finite element cache
using GridType = typename LevelContactNetwork::DeformedGridType;
- static const int dim = GridType::dimension;
-
using Element = typename GridType::template Codim<0>::Entity;
- //using Vertex = typename GridType::template Codim<dim>::Entity;
-
- using LevelGridView = typename LevelContactNetwork::DeformedLevelGridView;
- using LeafGridView = typename LevelContactNetwork::DeformedLeafGridView;
-
- using HierarchicLevelIteratorType = HierarchicLevelIterator<GridType>;
- using LevelElementMapper = typename Dune::MultipleCodimMultipleGeomTypeMapper<LevelGridView, Dune::MCMGElementLayout >;
- using LevelFaceMapper = typename Dune::MultipleCodimMultipleGeomTypeMapper<LevelGridView, mcmgLayout(Codim<1>()) >;
- using LeafElementMapper = typename Dune::MultipleCodimMultipleGeomTypeMapper<LeafGridView, Dune::MCMGElementLayout >;
- using LeafFaceMapper = typename Dune::MultipleCodimMultipleGeomTypeMapper<LeafGridView, mcmgLayout(Codim<1>()) >;
-
- bool levelMode_;
-
- std::unique_ptr<LevelGridView> levelGridView_;
- LevelElementMapper levelElementMapper_;
- LevelFaceMapper levelFaceMapper_;
+ static const int dim = LevelContactNetwork::dimension; //TODO
+ static const FiniteElementCache cache_;
- std::unique_ptr<LeafGridView> leafGridView_;
- LeafElementMapper leafElementMapper_;
- LeafFaceMapper leafFaceMapper_;
+ const LevelContactNetwork& levelContactNetwork_;
+ const size_t nBodies_;
- std::vector<std::vector<Element>>& vertexInElements_;
+ std::vector<size_t> vertexOffSet_;
+ std::vector<size_t> faceOffSet_;
+ std::vector<size_t> elementOffSet_;
- //std::vector<size_t> globalToLocal_;
+public:
+ NetworkIndexMapper(const LevelContactNetwork& levelContactNetwork) :
+ levelContactNetwork_(levelContactNetwork),
+ nBodies_(levelContactNetwork_.nBodies()),
+ vertexOffSet_(nBodies_),
+ faceOffSet_(nBodies_),
+ elementOffSet_(nBodies_) {
- template <class GV>
- void init(const GV& gv) {
- // init vertexInElements
- vertexInElements_.resize(gv.size(dim));
- for (size_t i=0; i<vertexInElements_.size(); i++) {
- vertexInElements_[i].resize(0);
- }
+ size_t vertexOffSet = 0;
+ size_t faceOffSet = 0;
+ size_t elementOffSet = 0;
- const auto& indexSet = gv.indexSet();
+ for (size_t i=0; i<nBodies_; i++) {
+ const auto& gridView = levelContactNetwork_.body(i)->gridView();
- for (const auto& elem : elements(gv)) {
- const auto& lfe = cache_.get(elem.type());
+ vertexOffSet_[i] = vertexOffSet + gridView.size(dim);
+ vertexOffSet = vertexOffSet_[i];
- for (size_t j=0; j<lfe.localBasis().size(); ++j) {
- auto localIdx = lfe.localCoefficients().localKey(j).subEntity();
- int dofIndex = indexSet.subIndex(elem, localIdx, dim);
+ faceOffSet_[i] = faceOffSet + gridView.size(1);
+ faceOffSet = faceOffSet_[i];
- vertexInElements_[dofIndex].push_back(elem);
- }
+ elementOffSet_[i] = elementOffSet + gridView.size(0);
+ elementOffSet = elementOffSet_[i];
}
}
-public:
-
- BodyHelper(const LevelGridView& gridView) :
- levelMode_(true),
- levelGridView_(std::make_unique<LevelGridView>(gridView)),
- levelElementMapper_(gridView),
- leafFaceMapper_(gridView) {
-
- init(gridView);
- }
- BodyHelper(const LeafGridView& gridView) :
- levelMode_(false),
- leafGridView_(std::make_unique<LeafGridView>(gridView)),
- leafElementMapper_(gridView),
- leafFaceMapper_(gridView) {
-
- init(gridView);
- }
+ size_t vertexIndex(const size_t bodyID, const Element& elem, const size_t localVertex) const {
+ const auto& gridView = levelContactNetwork_.body(bodyID)->gridView();
- const auto& gridView() const {
- return (levelMode_ ? levelGridView_ : leafGridView_);
+ auto localIdx = cache_.get(elem.type()).localCoefficients().localKey(localVertex).subEntity();
+ return vertexOffSet_[bodyID] + gridView.indexSet().subIndex(elem, localIdx, dim);
}
- const auto& elementMapper() const {
- return (levelMode_ ? levelElementMapper_ : leafElementMapper_);
+ template <class Intersection>
+ size_t faceIndex(const size_t bodyID, const Intersection& is) const {
+ const auto& gridView = levelContactNetwork_.body(bodyID)->gridView();
+ return faceOffSet_[bodyID] + gridView.indexSet().subIndex(is.inside(), is.indexInInside(), 1);
}
- const auto& faceMapper() const {
- return (levelMode_ ? levelFaceMapper_ : leafFaceMapper_);
+ size_t faceIndex(const size_t bodyID, const Element& elem, const size_t indexInInside) const {
+ const auto& gridView = levelContactNetwork_.body(bodyID)->gridView();
+ return faceOffSet_[bodyID] + gridView.indexSet().subIndex(elem, indexInInside, 1);
}
- const auto& vertexInElements() const {
- return vertexInElements_;
+ size_t elementIndex(const size_t bodyID, const Element& elem) const {
+ const auto& gridView = levelContactNetwork_.body(bodyID)->gridView();
+ return elementOffSet_[bodyID] + gridView.indexSet().index(elem);
}
- const auto& cache() const {
+ static const auto& cache() const {
return cache_;
}
};
-
/**
* constructs BitSetVector whose entries are
* true, if vertex is outside of patch or part of patch boundary
@@ -150,130 +99,139 @@ class SupportPatchFactory
using Patch = Dune::BitSetVector<1>;
private:
+ struct BodyElement {
+ size_t bodyID;
+ Element element;
+
+ void set(const size_t _bodyID, const Element& _elem) {
+ bodyID = _bodyID;
+ element = _elem;
+ }
+ };
+
+ struct RemoteIntersection {
+ size_t bodyID; // used to store bodyID of outside elem
+ GridGlue& glue;
+ size_t intersectionID;
+
+ bool flip = false;
+
+ void set(const size_t _bodyID, const GridGlue& _glue, const size_t _intersectionID, bool _flip = false) {
+ bodyID = _bodyID;
+ glue = _glue;
+ intersectionID = _intersectionID;
+ flip = _flip;
+ }
+
+ auto& get() {
+ return (flip) ? glue.getIntersection(intersectionID).flip() : glue.getIntersection(intersectionID);
+ }
+ };
+
static const int dim = LevelContactNetwork::dimension; //TODO
using ctype = typename LevelContactNetwork::ctype;
- using BodyHelper = BodyHelper<LevelContactNetwork>;
//using Vertex = typename GridType::template Codim<dim>::Entity;
using Element = typename LevelContactNetwork::DeformedGridType::template Codim<0>::Entity;
using BodyElement = BodyElement<Element>;
- using FiniteElementCache = Dune::PQkLocalFiniteElementCache<ctype, ctype, dim, 1>; // P1 finite element cache
+
+ using NetworkIndexMapper = NetworkIndexMapper<LevelContactNetwork>;
+
+ using GridGlue = ;
const size_t nBodies_;
const LevelContactNetwork& coarseContactNetwork_;
const LevelContactNetwork& fineContactNetwork_;
- std::vector<std::unique_ptr<BodyHelper>> bodyHelper_;
+ std::map<size_t, std::vector<RemoteIntersection>> coarseIntersectionMapper_; // map faceID to neighbor elements
+ std::map<size_t, std::vector<RemoteIntersection>> fineIntersectionMapper_; // map faceID to neighbor elements
+ std::map<size_t, std::vector<size_t>> neighborFaceDofs_; // map faceID to neighbor element dofs
- std::map<size_t, std::vector<BodyElement>> neighborElements_; // map faceID to neighbor elements
- std::map<size_t, std::vector<size_t>> neighborElementDofs_; // map faceID to neighbor element dofs
+ std::set<size_t> interiorContactDofs_;
- std::vector<size_t> vertexOffSet_;
- std::vector<size_t> faceOffSet_;
- std::vector<size_t> elementOffSet_;
-
- FiniteElementCache cache_;
+ const NetworkIndexMapper coarseIndices_;
+ const NetworkIndexMapper fineIndices_;
public:
SupportPatchFactory(const LevelContactNetwork& coarseContactNetwork, const LevelContactNetwork& fineContactNetwork) :
nBodies_(coarseContactNetwork_.nBodies()),
coarseContactNetwork_(coarseContactNetwork),
fineContactNetwork_(fineContactNetwork),
- bodyHelper_(nBodies_),
- vertexOffSet_(nBodies_),
- faceOffSet_(nBodies_),
- elementOffSet_(nBodies_){
+ coarseIndices_(coarseContactNetwork_),
+ fineIndices_(fineContactNetwork_) {
assert(nBodies_ == fineContactNetwork_.nBodies());
- size_t vertexOffSet = 0;
- size_t faceOffSet = 0;
- size_t elementOffSet = 0;
-
- for (size_t i=0; i<nBodies_; i++) {
- const auto& gridView = coarseContactNetwork_.body(i)->gridView();
- bodyHelper_[i] = std::make_unique<BodyHelper>(gridView);
-
- vertexOffSet_[i] = vertexOffSet + gridView.size(dim);
- vertexOffSet = vertexOffSet_[i];
-
- faceOffSet_[i] = faceOffSet + gridView.size(1);
- faceOffSet = faceOffSet_[i];
-
- elementOffSet_[i] = elementOffSet + gridView.size(0);
- elementOffSet = elementOffSet_[i];
- }
-
- // init boundary mapping/transfer
- const auto& nBodyAssembler = coarseContactNetwork_.nBodyAssembler();
- const auto& contactCouplings = nBodyAssembler.getContactCouplings();
+ setup(coarseContactNetwork_, coarseIndices_, coarseIntersectionMapper_);
+ setup(fineContactNetwork_, fineIndices_, fineIntersectionMapper_);
+ // setup neighborFaceDofs_
+ const auto& contactCouplings = coarseContactNetwork_.nBodyAssembler().getContactCouplings();
for (size_t i=0; i<contactCouplings.size(); i++) {
- const auto& coupling = nBodyAssembler.getCoupling(i);
- const auto& contactCoupling = contactCouplings[i];
-
- const size_t nmBody = coupling.gridIdx_[0];
- const size_t mBody = coupling.gridIdx_[1];
-
- const auto& nmBoundary = contactCoupling.nonmortarBoundary();
- const auto& mBoundary = contactCoupling.mortarBoundary();
-
- const auto& nmIndexSet = nonmortarBoundary.gridView().indexSet();
- const auto& mIndexSet = mortarBoundary.gridView().indexSet();
-
- for (const auto& rIs : intersections(*contactCoupling.getGlue())) {
- const auto& inside = rIs.inside();
- const auto& outside = rIs.outside();
-
- int nmFaceIdx = nmIndexSet.subIndex(inside, rIs.indexInInside(), 1);
- int mFaceIdx = mIndexSet.subIndex(outside, rIs.indexInOutside(), 1);
-
- BodyElement nmBodyElement(nmBody, inside);
- BodyElement mBodyElement(mBody, outside);
+ const auto& coupling = nBodyAssembler.getCoupling(i);
+ const auto& glue = *contactCouplings[i].getGlue();
- neighborElements_[nmFaceIdx].emplace_back(mBodyElement);
- neighborElements_[mFaceIdx].emplace_back(nmBodyElement);
+ const size_t nmBody = coupling.gridIdx_[0];
+ const size_t mBody = coupling.gridIdx_[1];
- std::set<size_t> dofs;
- computeIntersectionDofs(nmIndexSet, rIs, dofs);
- neighborElementDofs_[nmFaceIdx].insert(neighborElementDofs_[nmFaceIdx].end(), dofs.begin(), dofs.end());
+ for (const auto& rIs : intersections(glue)) {
+ size_t nmFaceIdx = coarseIndices_.faceIndex(nmBody, rIs.inside(), rIs.indexInInside());
+ size_t mFaceIdx = coarseIndices_.faceIndex(mBody, rIs.outside(), rIs.indexInOutside());
- computeIntersectionDofs(mIndexSet, rIs, dofs, false);
- neighborElementDofs_[mFaceIdx].insert(neighborElementDofs_[mFaceIdx].end(), dofs.begin(), dofs.end());
+ std::array<std::set<size_t>, 2> dofs;
+ remoteIntersectionDofs(coarseIndices_, rIs, nmBody, mBody, dofs);
+ neighborFaceDofs_[nmFaceIdx].insert(neighborFaceDofs_[nmFaceIdx].end(), dofs[1].begin(), dofs[1].end());
+ neighborFaceDofs_[mFaceIdx].insert(neighborFaceDofs_[mFaceIdx].end(), dofs[0].begin(), dofs[0].end());
+ }
}
- }
- // neighborElementDofs_ contains dofs of connected intersections that are neighbors to a face given by faceID,
- // boundary dofs are contained once, interior dofs multiple times
- // task: remove boundary dofs and redundant multiples of interior dofs
- for (auto& it : neighborElementDofs_) {
- auto& dofs = it.second;
-
- std::sort(dofs.begin(), dofs.end());
- auto& dof = dofs.begin();
- while (dof != dofs.end()) {
- auto next = dof;
- if (*dof != *(++next))
- dof = dofs.erase(dof);
+ /*
+ // neighborElementDofs_ contains dofs of connected intersections that are neighbors to a face given by faceID,
+ // boundary dofs are contained once, interior dofs multiple times
+ // task: remove boundary dofs and redundant multiples of interior dofs
+ // TODO: only works in 2D
+ for (auto& it : neighborElementDofs) {
+ auto& dofs = it.second;
+
+ std::sort(dofs.begin(), dofs.end());
+ auto& dof = dofs.begin();
+ while (dof != dofs.end()) {
+ auto next = dof;
+ if (*dof != *(++next))
+ dof = dofs.erase(dof);
+ }
+ dofs.erase(std::unique(dofs.begin(), dofs.end()), dofs.end());
}
- dofs.erase(std::unique(dofs.begin(), dofs.end()), dofs.end());
- }
+ */
}
void build(const size_t bodyID, const Element& coarseElement, const size_t localVertex, Patch& patchDofs, const int patchDepth = 0) {
BodyElement seedElement(bodyID, coarseElement);
+ auto isPatchIntersection = [](auto& is, const size_t bodyID, const std::set<size_t>& patchVertices, std::set<size_t>& newPatchVertices) {
+ std::set<size_t> dofs;
+ intersectionDofs(coarseIndices_, is, bodyID, dofs);
+
+ for (auto& dof : dofs) {
+ if (patchVertices.count(dof)) {
+ newPatchVertices.insert(dofs.begin(), dofs.end());
+ return true;
+ }
+ }
+ return false;
+ };
+
// construct coarse patch
std::vector<BodyElement> patchElements;
patchElements.emplace_back(seedElement);
std::set<size_t> visitedElements;
- visitedElements.insert(elementIndex(seedElement));
+ visitedElements.insert(coarseIndices_.elementIndex(seedElement.bodyID, seedElement.element));
std::set<size_t> patchVertices;
- patchVertices.insert(vertexIndex(bodyID, coarseElement, localVertex));
+ patchVertices.insert(coarseIndices_.vertexIndex(bodyID, coarseElement, localVertex));
std::queue<BodyElement> patchSeeds;
patchSeeds.push(seedElement);
@@ -286,36 +244,37 @@ class SupportPatchFactory
const auto& patchSeed = patchSeeds.front();
patchSeeds.pop();
- size_t elemIdx = elementIndex(patchSeed);
+ size_t elemIdx = coarseIndices_.elementIndex(patchSeed.bodyID, patchSeed.element);
if (visitedElements.count(elemIdx))
continue;
visitedElements.insert(elemIdx);
- size_t bodyIdx = patchSeed.bodyID();
- const auto& elem = patchSeed.element();
- const auto& gridView = bodyHelper_[bodyIdx]->gridView();
+ size_t bodyIdx = patchSeed.bodyID;
+ const auto& elem = patchSeed.element;
+ const auto& gridView = coarseContactNetwork_.body(bodyIdx)->gridView();
for (auto& it : intersections(gridView, elem)) {
- if (isPatchIntersection(gridView.indexSet(), it, patchVertices, newPatchVertices)) {
+ if (isPatchIntersection(it, bodyIdx, patchVertices, newPatchVertices)) {
if (it.neighbor()) {
BodyElement neighbor(bodyIdx, it.outside());
- if (visitedElements.count(elementIndex(neighbor)))
+ if (visitedElements.count(coarseIndices_.elementIndex(neighbor.bodyID, neighbor.element)))
continue;
patchElements.emplace_back(neighbor);
patchSeeds.push(neighbor);
} else {
- size_t faceIdx = faceIndex(bodyIdx, it);
+ size_t faceIdx = coarseIndices_.faceIndex(bodyIdx, it);
- if (neighborElements_.count(faceIdx)) {
- const auto& neighbors = neighborElements_[faceIdx];
- for (size_t i=0; i<neighbors.size(); i++) {
- const auto& neighbor = neighbors[i];
+ if (coarseIntersectionMapper_.count(faceIdx)) {
+ const auto& intersections = coarseIntersectionMapper_[faceIdx];
+ for (size_t i=0; i<intersections.size(); i++) {
+ BodyElement neighbor;
+ neighbor.set(intersections[i].bodyID, intersections[i].get().outside());
- if (visitedElements.count(elementIndex(neighbor)))
+ if (visitedElements.count(coarseIndices_.elementIndex(neighbor.bodyID, neighbor.element)))
continue;
patchVertices.insert(neighborElementDofs_[faceIdx].begin(), neighborElementDofs_[faceIdx].end());
@@ -328,19 +287,20 @@ class SupportPatchFactory
if (it.neighbor()) {
BodyElement neighbor(bodyIdx, it.outside());
- if (visitedElements.count(elementIndex(neighbor)))
+ if (visitedElements.count(coarseIndices_.elementIndex(neighbor.bodyID, neighbor.element)))
continue;
nextSeeds.emplace_back(neighbor);
} else {
- size_t faceIdx = faceIndex(bodyIdx, it);
+ size_t faceIdx = coarseIndices_.faceIndex(bodyIdx, it);
- if (neighborElements_.count(faceIdx)) {
- const auto& neighbors = neighborElements_[faceIdx];
- for (size_t i=0; i<neighbors.size(); i++) {
- const auto& neighbor = neighbors[i];
+ if (coarseIntersectionMapper_.count(faceIdx)) {
+ const auto& intersections = coarseIntersectionMapper_[faceIdx];
+ for (size_t i=0; i<intersections.size(); i++) {
+ BodyElement neighbor;
+ neighbor.set(intersections[i].bodyID, intersections[i].get().outside());
- if (visitedElements.count(elementIndex(neighbor)))
+ if (visitedElements.count(coarseIndices_.elementIndex(neighbor.bodyID, neighbor.element)))
continue;
newPatchVertices.insert(neighborElementDofs_[faceIdx].begin(), neighborElementDofs_[faceIdx].end());
@@ -359,42 +319,50 @@ class SupportPatchFactory
}
// construct fine patch
+ using HierarchicLevelIteratorType = HierarchicLevelIterator<typename LevelContactNetwork::DeformedGridType>;
+ std::set<size_t> patchBoundary;
for (size_t i=0; i<patchElements.size(); ++i) {
- addFinePatchElements(patchElements[i], visitedElements, patchDofs);
- }
- }
+ const auto& coarseElem = patchElements[i];
-private:
- template <class IndexSet, class Intersection>
- bool isPatchIntersection(const IndexSet& indexSet, const Intersection& is, const std::set<size_t>& patchVertices, std::set<size_t>& newPatchVertices) {
- std::set<size_t> intersectionDofs;
- computeIntersectionDofs(indexSet, is, intersectionDofs);
-
- for (auto& dof : intersectionDofs) {
- if (patchVertices.count(dof)) {
- newPatchVertices.insert(intersectionDofs.begin(), intersectionDofs.end());
- return true;
+ const auto& grid = coarseContactNetwork_.body(coarseElem.bodyID)->gridView().grid();
+ const auto fineLevel = fineContactNetwork_.body(coarseElem.bodyID)->gridView().level();
+
+ // add fine patch elements
+ HierarchicLevelIteratorType endHierIt(grid, coarseElem.element, HierarchicLevelIteratorType::PositionFlag::end, fineLevel);
+ HierarchicLevelIteratorType hierIt(grid, coarseElem.element, HierarchicLevelIteratorType::PositionFlag::begin, fineLevel);
+ for (; hierIt!=endHierIt; ++hierIt) {
+ const Element& fineElem = *hierIt;
+ addLocalDofs(coarseElem, fineElem, visitedElements, patchBoundary, patchDofs);
}
}
- return false;
}
- size_t vertexIndex(const size_t bodyID, const Element& elem, const size_t localVertex) const {
- const auto& gridView = bodyHelper_[bodyID]->gridView();
+private:
+ void setup(const LevelContactNetwork& levelContactNetwork, const NetworkIndexMapper& indices, std::map<size_t, std::vector<RemoteIntersection>>& faceToRIntersections) {
+ const auto& nBodyAssembler = levelContactNetwork.nBodyAssembler();
+ const auto& contactCouplings = nBodyAssembler.getContactCouplings();
- auto localIdx = cache_.get(elem.type()).localCoefficients().localKey(localVertex).subEntity();
- return vertexOffSet_[bodyID] + gridView.indexSet().subIndex(elem, localIdx, dim);
- }
+ for (size_t i=0; i<contactCouplings.size(); i++) {
+ const auto& coupling = nBodyAssembler.getCoupling(i);
+ const auto& glue = *contactCouplings[i].getGlue();
- template <class Intersection>
- size_t faceIndex(const size_t bodyID, const Intersection& is) const {
- const auto& gridView = bodyHelper_[bodyID]->gridView();
- return faceOffSet_[bodyID] + gridView.indexSet().subIndex(is.inside(), is.indexInInside(), 1);
- }
+ const size_t nmBody = coupling.gridIdx_[0];
+ const size_t mBody = coupling.gridIdx_[1];
+
+ size_t rIsID = 0;
+ for (const auto& rIs : intersections(glue)) {
+ size_t nmFaceIdx = indices.faceIndex(nmBody, rIs.inside(), rIs.indexInInside());
+ size_t mFaceIdx = indices.faceIndex(mBody, rIs.outside(), rIs.indexInOutside());
+
+ RemoteIntersection nmIntersection(mBody, glue, rIsID, false);
+ faceToRIntersections[nmFaceIdx].emplace_back(nmIntersection);
- size_t elementIndex(const BodyElement& bodyElem) const {
- const auto& gridView = bodyHelper_[bodyElem.bodyID()]->gridView();
- return elementOffSet_[bodyElem.bodyID()] + gridView.indexSet().index(bodyElem.element());
+ RemoteIntersection mIntersection(nmBody, glue, rIsID, true);
+ faceToRIntersections[mFaceIdx].emplace_back(mIntersection);
+
+ rIsID++;
+ }
+ }
}
template <class Intersection>
@@ -402,112 +370,133 @@ class SupportPatchFactory
return ReferenceElementHelper<double, dim>::subEntityContainsSubEntity(elem.type(), intersection.indexInInside(), 1, subEntity, codim);
}
- void addFinePatchElements(const BodyElement& coarseElem, const std::set<size_t>& coarsePatch, Patch& patchDofs) {
- using HierarchicLevelIteratorType = HierarchicLevelIterator<typename LevelContactNetwork::DeformedGridType>;
-
- const auto& grid = coarseContactNetwork_.body(coarseElem.bodyID())->gridView().grid();
- const auto fineLevel = fineContactNetwork_.body(coarseElem.bodyID())->gridView().level();
-
- HierarchicLevelIteratorType endHierIt(grid, coarseElem.element(), HierarchicLevelIteratorType::PositionFlag::end, fineLevel);
- HierarchicLevelIteratorType hierIt(grid, coarseElem.element(), HierarchicLevelIteratorType::PositionFlag::begin, fineLevel);
- for (; hierIt!=endHierIt; ++hierIt) {
- const Element& fineElem = *hierIt;
- addLocalDofs(coarseElem, fineElem, coarsePatch, patchDofs);
- }
- }
-
- void addLocalDofs(const BodyElement& coarseElem, const Element& fineElem, const std::set<size_t>& coarsePatch, Patch& patchDofs) {
+ void addLocalDofs(const BodyElement& coarseElem, const Element& fineElem, const std::set<size_t>& coarsePatch, std::set<size_t>& patchBoundary, Patch& patchDofs) {
// insert local dofs
const auto& lfe = cache_.get(fineElem.type());
- const auto& gridView = fineContactNetwork_.body(coarseElem.bodyID())->gridView();
+ const auto& gridView = fineContactNetwork_.body(coarseElem.bodyID)->gridView();
const auto& indexSet = gridView.indexSet();
-
- for (size_t i=0; i<lfe.localBasis().size(); ++i) {
- auto localIdx = lfe.localCoefficients().localKey(i).subEntity();
- int dofIndex = indexSet.subIndex(fineElem, localIdx, dim);
- patchDofs[dofIndex][0] = false;
- }
- const auto coarseLevel = coarseElem.element().level();
+ const size_t bodyID = coarseElem.bodyID;
+
+ const auto coarseLevel = coarseElem.element.level();
// search for parts of the global and inner boundary
for (const auto& is : intersections(gridView, fineElem)) {
- bool isPatchBoundary = false;
-
if (is.neighbor()) {
- auto outsideFather = is.outside();
-
- while (outsideFather.level()>coarseLevel) {
- outsideFather = outsideFather.father();
- }
-
- BodyElement outside(coarseElem.bodyID(), outsideFather);
+ std::set<size_t> isDofs;
+ intersectionDofs(fineIndices_, is, bodyID, isDofs);
- if (!coarsePatch.count(elementIndex(outside))) {
- isPatchBoundary = true;
+ auto outsideFather = coarseFather(is.outside(), coarseLevel);
+
+ if (!coarsePatch.count(coarseIndices_.elementIndex(coarseElem.bodyID, outsideFather))) {
+ addToPatchBoundary(isDofs, patchBoundary, patchDofs);
+ } else {
+ addToPatch(isDofs, patchBoundary, patchDofs);
}
} else {
- isPatchBoundary = true;
+ size_t faceIdx = fineIndices_.faceIndex(bodyID, is);
+ if (fineIntersectionMapper_.count(faceIdx)) {
+ const auto& intersections = coarseIntersectionMapper_[faceIdx];
- auto outsideFather = is.outside();
+ for (size_t i=0; i<intersections.size(); i++) {
+ const auto& intersection = intersections[i];
+ const auto& rIs = intersection.get();
+ const size_t outsideID = intersection.bodyID;
- while (outsideFather.level()>coarseLevel) {
- outsideFather = outsideFather.father();
- }
+ std::set<size_t> rIsDofs;
+ remoteIntersectionDofs(fineIndices_, rIs, bodyID, outsideID, rIsDofs);
- BodyElement outside(coarseElem.bodyID(), outsideFather);
+ if (rIs.neighbor()) {
+ auto outsideFather = coarseFather(rIs.outside(), fineContactNetwork_.body(outsideID)->gridView().level());
- size_t faceIdx = faceIndex(coarseElem.bodyID(), is);
-
- if (neighborElements_.count(faceIdx)) {
- const auto& outside = neighbors[0];
-
- if (coarsePatch.count(elementIndex(outside))) {
- isPatchBoundary = false;
+ if (!coarsePatch.count(coarseIndices_.elementIndex(outsideID, outsideFather))) {
+ addToPatchBoundary(rIsDofs, patchBoundary, patchDofs);
+ } else {
+ addToPatch(rIsDofs, patchBoundary, patchDofs);
+ }
+ } else {
+ addToPatchBoundary(rIsDofs, patchBoundary, patchDofs);
+ }
}
+ } else {
+ std::set<size_t> isDofs;
+ intersectionDofs(fineIndices_, is, bodyID, isDofs);
+
+ addToPatchBoundary(isDofs, patchBoundary, patchDofs);
}
}
-
- if (isPatchBoundary) {
- const auto& localCoefficients = cache_.get(fineElem.type()).localCoefficients();
+ }
+ }
+
+ auto&& coarseFather(const Element& fineElem, const size_t coarseLevel) const {
+ Element coarseElem = fineElem;
+ while (coarseElem.level() > coarseLevel) {
+ coarseElem = coarseElem.father();
+ }
+ return coarseElem;
+ }
+
+ void addToPatchBoundary(const std::set<size_t>& dofs, std::set<size_t>& patchBoundary, Patch& patchDofs) {
+ for (auto& dof : dofs) {
+ patchBoundary.insert(dof);
+ patchDofs[dof][0] = true;
+ }
+ }
+
+ template <class Intersection>
+ void addToPatch(const std::set<size_t>& dofs, const std::set<size_t>& patchBoundary, Patch& patchDofs) {
+ for (auto& dof : dofs) {
+ if (!patchBoundary.count(dof)) {
+ patchDofs[dof][0] = false;
+ }
+ }
+ }
- for (size_t i=0; i<localCoefficients.size(); i++) {
- auto entity = localCoefficients.localKey(i).subEntity();
- auto codim = localCoefficients.localKey(i).codim();
+ template <class RIntersection>
+ void remoteIntersectionDofs(const NetworkIndexMapper& indices, const RIntersection& is, const size_t insideID, const size_t outsideID, std::array<std::set<size_t>, 2>& dofs) {
+ dofs.clear();
- if (containsInsideSubentity(fineElem, is, entity, codim)) {
- int dofIndex = indexSet.subIndex(fineElem, entity, dim);
- patchDofs[dofIndex][0] = true;
- }
- }
+ const auto& inside = is.inside();
+ const auto& insideGeo = inside.geometry();
+ const auto& insideRefElement = Dune::ReferenceElements<ctype, dim>::general(inside.type());
+
+ for (int i=0; i<insideRefElement.size(dim); i++) {
+ const auto& localCorner = isGeo.local(insideGeo.corner(i));
+
+ if (isRefElement.checkInside(localCorner)) {
+ dofs[0].insert(indices.vertexIndex(insideID, inside, i));
}
- }
- }
+ }
+ const auto& isGeo = is.geometry();
+ const auto& isRefElement = Dune::ReferenceElements<ctype, dim-1>::general(is.type());
- template <class IndexSet, class Intersection>
- void computeIntersectionDofs(const IndexSet& idxSet, const Intersection& intersection, std::set<size_t>& intersectionDofs, bool inside = true) {
- intersectionDofs.clear();
+ const auto& outside = is.outside();
+ const auto& outsideGeo = outside.geometry();
+ const auto& outsideRefElement = Dune::ReferenceElements<ctype, dim>::general(outside.type());
- Element* elem;
- size_t intersectionIndex;
+ for (int i=0; i<outsideRefElement.size(dim); i++) {
+ const auto& localCorner = isGeo.local(outsideGeo.corner(i));
- if (inside) {
- elem = &intersection.inside();
- intersectionIndex = intersection.indexInInside();
- } else {
- elem = &intersection.outside();
- intersectionIndex = intersection.indexInOutside();
+ if (isRefElement.checkInside(localCorner)) {
+ dofs[1].insert(indices.vertexIndex(outsideID, outside, i));
+ }
}
+ }
+
+ template <class Intersection>
+ void intersectionDofs(const NetworkIndexMapper& indices, const Intersection& is, const size_t bodyID, std::set<size_t>& dofs) {
+ dofs.clear();
+
+ const Element& elem = is.inside();
+ size_t intersectionIndex = is.indexInInside();
const int dimElement = Element::dimension;
- const auto& refElement = Dune::ReferenceElements<double, dimElement>::general(elem->type());
+ const auto& refElement = Dune::ReferenceElements<double, dimElement>::general(elem.type());
for (int i=0; i<refElement.size(intersectionIndex, 1, dimElement); i++) {
size_t idxInElement = refElement.subEntity(intersectionIndex, 1, i, dimElement);
- size_t globalIdx = idxSet.subIndex(*elem, idxInElement, dimElement);
-
- intersectionDofs.insert(globalIdx);
+ dofs.insert(indices.vertexIndex(bodyID, elem, idxInElement));
}
}
};
diff --git a/src/spatial-solving/preconditioners/supportpatchfactory.hh.autosave b/src/spatial-solving/preconditioners/supportpatchfactory.hh.autosave
deleted file mode 100644
index 869cbd83d4f0cf01737f5d7d41f7cd53c392da96..0000000000000000000000000000000000000000
--- a/src/spatial-solving/preconditioners/supportpatchfactory.hh.autosave
+++ /dev/null
@@ -1,518 +0,0 @@
-#ifndef SUPPORT_PATCH_FACTORY_HH
-#define SUPPORT_PATCH_FACTORY_HH
-
-#include <set>
-#include <queue>
-#include <memory>
-
-#include <dune/common/bitsetvector.hh>
-#include <dune/common/fvector.hh>
-
-#include <dune/fufem/referenceelementhelper.hh>
-#include <dune/grid/common/mcmgmapper.hh>
-
-#include <dune/localfunctions/lagrange/pqkfactory.hh>
-
-//#include "../../data-structures/levelcontactnetwork.hh"
-#include "hierarchicleveliterator.hh"
-
-template <class Element>
-class BodyElement {
-private:
- size_t bodyID_;
- Element element_;
-
-public:
- BodyElement() {}
- BodyElement(size_t bodyID, const Element& element) {
- set(bodyID, element);
- }
-
- void set(size_t bodyID, const Element& element) {
- bodyID_ = bodyID;
- element_ = element;
- }
-
- const auto& bodyID() const {
- return bodyID_;
- }
- const auto& element() const {
- return element_;
- }
-};
-
-
-template <class LevelContactNetwork>
-class BodyHelper {
-private:
- using GridType = typename LevelContactNetwork::DeformedGridType;
- static const int dim = GridType::dimension;
-
- using Element = typename GridType::template Codim<0>::Entity;
- //using Vertex = typename GridType::template Codim<dim>::Entity;
-
- using LevelGridView = typename LevelContactNetwork::DeformedLevelGridView;
- using LeafGridView = typename LevelContactNetwork::DeformedLeafGridView;
-
- using HierarchicLevelIteratorType = HierarchicLevelIterator<GridType>;
- using LevelElementMapper = typename Dune::MultipleCodimMultipleGeomTypeMapper<LevelGridView, Dune::MCMGElementLayout >;
- using LevelFaceMapper = typename Dune::MultipleCodimMultipleGeomTypeMapper<LevelGridView, mcmgLayout(Codim<1>()) >;
-
- using LeafElementMapper = typename Dune::MultipleCodimMultipleGeomTypeMapper<LeafGridView, Dune::MCMGElementLayout >;
- using LeafFaceMapper = typename Dune::MultipleCodimMultipleGeomTypeMapper<LeafGridView, mcmgLayout(Codim<1>()) >;
-
- bool levelMode_;
-
- std::unique_ptr<LevelGridView> levelGridView_;
- LevelElementMapper levelElementMapper_;
- LevelFaceMapper levelFaceMapper_;
-
- std::unique_ptr<LeafGridView> leafGridView_;
- LeafElementMapper leafElementMapper_;
- LeafFaceMapper leafFaceMapper_;
-
- std::vector<std::vector<Element>>& vertexInElements_;
-
- //std::vector<size_t> globalToLocal_;
-
- template <class GV>
- void init(const GV& gv) {
- // init vertexInElements
- vertexInElements_.resize(gv.size(dim));
- for (size_t i=0; i<vertexInElements_.size(); i++) {
- vertexInElements_[i].resize(0);
- }
-
- const auto& indexSet = gv.indexSet();
-
- for (const auto& elem : elements(gv)) {
- const auto& lfe = cache_.get(elem.type());
-
- for (size_t j=0; j<lfe.localBasis().size(); ++j) {
- auto localIdx = lfe.localCoefficients().localKey(j).subEntity();
- int dofIndex = indexSet.subIndex(elem, localIdx, dim);
-
- vertexInElements_[dofIndex].push_back(elem);
- }
- }
- }
-public:
-
- BodyHelper(const LevelGridView& gridView) :
- levelMode_(true),
- levelGridView_(std::make_unique<LevelGridView>(gridView)),
- levelElementMapper_(gridView),
- leafFaceMapper_(gridView) {
-
- init(gridView);
- }
-
- BodyHelper(const LeafGridView& gridView) :
- levelMode_(false),
- leafGridView_(std::make_unique<LeafGridView>(gridView)),
- leafElementMapper_(gridView),
- leafFaceMapper_(gridView) {
-
- init(gridView);
- }
-
- const auto& gridView() const {
- return (levelMode_ ? levelGridView_ : leafGridView_);
- }
-
- const auto& elementMapper() const {
- return (levelMode_ ? levelElementMapper_ : leafElementMapper_);
- }
-
- const auto& faceMapper() const {
- return (levelMode_ ? levelFaceMapper_ : leafFaceMapper_);
- }
-
- const auto& vertexInElements() const {
- return vertexInElements_;
- }
-
- const auto& cache() const {
- return cache_;
- }
-};
-
-
-/**
- * constructs BitSetVector whose entries are
- * true, if vertex is outside of patch or part of patch boundary
- * false, if vertex is in interior of patch
- */
-template <class LevelContactNetwork>
-class SupportPatchFactory
-{
-public:
- using Patch = Dune::BitSetVector<1>;
-
-private:
- static const int dim = LevelContactNetwork::dimension; //TODO
- using ctype = typename LevelContactNetwork::ctype;
-
- using BodyHelper = BodyHelper<LevelContactNetwork>;
- //using Vertex = typename GridType::template Codim<dim>::Entity;
-
- using Element = typename LevelContactNetwork::DeformedGridType::template Codim<0>::Entity;
- using BodyElement = BodyElement<Element>;
- using FiniteElementCache = Dune::PQkLocalFiniteElementCache<ctype, ctype, dim, 1>; // P1 finite element cache
-
- const size_t nBodies_;
-
- const LevelContactNetwork& coarseContactNetwork_;
- const LevelContactNetwork& fineContactNetwork_;
-
- std::vector<std::unique_ptr<BodyHelper>> bodyHelper_;
-
- std::map<size_t, std::vector<BodyElement>> neighborElements_; // map faceID to neighbor elements
- std::map<size_t, std::vector<size_t>> neighborElementDofs_; // map faceID to neighbor element dofs
-
- std::vector<size_t> vertexOffSet_;
- std::vector<size_t> faceOffSet_;
- std::vector<size_t> elementOffSet_;
-
- FiniteElementCache cache_;
-
-public:
- SupportPatchFactory(const LevelContactNetwork& coarseContactNetwork, const LevelContactNetwork& fineContactNetwork) :
- nBodies_(coarseContactNetwork_.nBodies()),
- coarseContactNetwork_(coarseContactNetwork),
- fineContactNetwork_(fineContactNetwork),
- bodyHelper_(nBodies_),
- vertexOffSet_(nBodies_),
- faceOffSet_(nBodies_),
- elementOffSet_(nBodies_){
-
- assert(nBodies_ == fineContactNetwork_.nBodies());
-
- size_t vertexOffSet = 0;
- size_t faceOffSet = 0;
- size_t elementOffSet = 0;
-
- for (size_t i=0; i<nBodies_; i++) {
- const auto& gridView = coarseContactNetwork_.body(i)->gridView();
- bodyHelper_[i] = std::make_unique<BodyHelper>(gridView);
-
- vertexOffSet_[i] = vertexOffSet + gridView.size(dim);
- vertexOffSet = vertexOffSet_[i];
-
- faceOffSet_[i] = faceOffSet + gridView.size(1);
- faceOffSet = faceOffSet_[i];
-
- elementOffSet_[i] = elementOffSet + gridView.size(0);
- elementOffSet = elementOffSet_[i];
- }
-
- // init boundary mapping/transfer
- const auto& nBodyAssembler = coarseContactNetwork_.nBodyAssembler();
- const auto& contactCouplings = nBodyAssembler.getContactCouplings();
-
- for (size_t i=0; i<contactCouplings.size(); i++) {
- const auto& coupling = nBodyAssembler.getCoupling(i);
- const auto& contactCoupling = contactCouplings[i];
-
- const size_t nmBody = coupling.gridIdx_[0];
- const size_t mBody = coupling.gridIdx_[1];
-
- const auto& nmBoundary = contactCoupling.nonmortarBoundary();
- const auto& mBoundary = contactCoupling.mortarBoundary();
-
- const auto& nmIndexSet = nonmortarBoundary.gridView().indexSet();
- const auto& mIndexSet = mortarBoundary.gridView().indexSet();
-
- for (const auto& rIs : intersections(*contactCoupling.getGlue())) {
- const auto& inside = rIs.inside();
- const auto& outside = rIs.outside();
-
- int nmFaceIdx = nmIndexSet.subIndex(inside, rIs.indexInInside(), 1);
- int mFaceIdx = mIndexSet.subIndex(outside, rIs.indexInOutside(), 1);
-
- BodyElement nmBodyElement(nmBody, inside);
- BodyElement mBodyElement(mBody, outside);
-
- neighborElements_[nmFaceIdx].emplace_back(mBodyElement);
- neighborElements_[mFaceIdx].emplace_back(nmBodyElement);
-
- std::set<size_t> dofs;
- computeIntersectionDofs(nmIndexSet, rIs, dofs);
- neighborElementDofs_[nmFaceIdx].insert(neighborElementDofs_[nmFaceIdx].end(), dofs.begin(), dofs.end());
-
- computeIntersectionDofs(mIndexSet, rIs, dofs, false);
- neighborElementDofs_[mFaceIdx].insert(neighborElementDofs_[mFaceIdx].end(), dofs.begin(), dofs.end());
- }
- }
-
- // neighborElementDofs_ contains dofs of connected intersections that are neighbors to a face given by faceID,
- // boundary dofs are contained once, interior dofs multiple times
- // task: remove boundary dofs and redundant multiples of interior dofs
- for (auto& it : neighborElementDofs_) {
- auto& dofs = it.second;
-
- std::sort(dofs.begin(), dofs.end());
- auto& dof = dofs.begin();
- while (dof != dofs.end()) {
- auto next = dof;
- if (*dof != *(++next))
- dof = dofs.erase(dof);
- }
- dofs.erase(std::unique(dofs.begin(), dofs.end()), dofs.end());
- }
- }
-
- void build(const size_t bodyID, const Element& coarseElement, const size_t localVertex, Patch& patchDofs, const int patchDepth = 0) {
- BodyElement seedElement(bodyID, coarseElement);
-
- // construct coarse patch
- std::vector<BodyElement> patchElements;
- patchElements.emplace_back(seedElement);
-
- std::set<size_t> visitedElements;
- visitedElements.insert(elementIndex(seedElement));
-
- std::set<size_t> patchVertices;
- patchVertices.insert(vertexIndex(bodyID, coarseElement, localVertex));
-
- std::queue<BodyElement> patchSeeds;
- patchSeeds.push(seedElement);
-
- for (size_t depth=0; depth<=patchDepth; depth++) {
- std::vector<BodyElement> nextSeeds;
- std::set<size_t> newPatchVertices;
-
- while (!patchSeeds.empty()) {
- const auto& patchSeed = patchSeeds.front();
- patchSeeds.pop();
-
- size_t elemIdx = elementIndex(patchSeed);
-
- if (visitedElements.count(elemIdx))
- continue;
-
- visitedElements.insert(elemIdx);
-
- size_t bodyIdx = patchSeed.bodyID();
- const auto& elem = patchSeed.element();
- const auto& gridView = bodyHelper_[bodyIdx]->gridView();
-
- for (auto& it : intersections(gridView, elem)) {
- if (isPatchIntersection(gridView.indexSet(), it, patchVertices, newPatchVertices)) {
- if (it.neighbor()) {
- BodyElement neighbor(bodyIdx, it.outside());
-
- if (visitedElements.count(elementIndex(neighbor)))
- continue;
-
- patchElements.emplace_back(neighbor);
- patchSeeds.push(neighbor);
- } else {
- size_t faceIdx = faceIndex(bodyIdx, it);
-
- if (neighborElements_.count(faceIdx)) {
- const auto& neighbors = neighborElements_[faceIdx];
- for (size_t i=0; i<neighbors.size(); i++) {
- const auto& neighbor = neighbors[i];
-
- if (visitedElements.count(elementIndex(neighbor)))
- continue;
-
- patchVertices.insert(neighborElementDofs_[faceIdx].begin(), neighborElementDofs_[faceIdx].end());
- patchElements.emplace_back(neighbor);
- patchSeeds.push(neighbor);
- }
- }
- }
- } else {
- if (it.neighbor()) {
- BodyElement neighbor(bodyIdx, it.outside());
-
- if (visitedElements.count(elementIndex(neighbor)))
- continue;
-
- nextSeeds.emplace_back(neighbor);
- } else {
- size_t faceIdx = faceIndex(bodyIdx, it);
-
- if (neighborElements_.count(faceIdx)) {
- const auto& neighbors = neighborElements_[faceIdx];
- for (size_t i=0; i<neighbors.size(); i++) {
- const auto& neighbor = neighbors[i];
-
- if (visitedElements.count(elementIndex(neighbor)))
- continue;
-
- newPatchVertices.insert(neighborElementDofs_[faceIdx].begin(), neighborElementDofs_[faceIdx].end());
- nextSeeds.emplace_back(neighbor);
- }
- }
- }
- }
- }
- }
-
- for (size_t i=0; i<nextSeeds.size(); i++) {
- patchSeeds.push(nextSeeds[i]);
- }
- patchVertices.insert(newPatchVertices.begin(), newPatchVertices.end());
- }
-
- // construct fine patch
- for (size_t i=0; i<patchElements.size(); ++i) {
- addFinePatchElements(patchElements[i], visitedElements, patchDofs);
- }
- }
-
-private:
- template <class IndexSet, class Intersection>
- bool isPatchIntersection(const IndexSet& indexSet, const Intersection& is, const std::set<size_t>& patchVertices, std::set<size_t>& newPatchVertices) {
- std::set<size_t> intersectionDofs;
- computeIntersectionDofs(indexSet, is, intersectionDofs);
-
- for (auto& dof : intersectionDofs) {
- if (patchVertices.count(dof)) {
- newPatchVertices.insert(intersectionDofs.begin(), intersectionDofs.end());
- return true;
- }
- }
- return false;
- }
-
- size_t vertexIndex(const size_t bodyID, const Element& elem, const size_t localVertex) const {
- const auto& gridView = bodyHelper_[bodyID]->gridView();
-
- auto localIdx = cache_.get(elem.type()).localCoefficients().localKey(localVertex).subEntity();
- return vertexOffSet_[bodyID] + gridView.indexSet().subIndex(elem, localIdx, dim);
- }
-
- template <class Intersection>
- size_t faceIndex(const size_t bodyID, const Intersection& is) const {
- const auto& gridView = bodyHelper_[bodyID]->gridView();
- return faceOffSet_[bodyID] + gridView.indexSet().subIndex(is.inside(), is.indexInInside(), 1);
- }
-
- size_t elementIndex(const BodyElement& bodyElem) const {
- const auto& gridView = bodyHelper_[bodyElem.bodyID()]->gridView();
- return elementOffSet_[bodyElem.bodyID()] + gridView.indexSet().index(bodyElem.element());
- }
-
- template <class Intersection>
- bool containsInsideSubentity(const Element& elem, const Intersection& intersection, int subEntity, int codim) {
- return ReferenceElementHelper<double, dim>::subEntityContainsSubEntity(elem.type(), intersection.indexInInside(), 1, subEntity, codim);
- }
-
- void addFinePatchElements(const BodyElement& coarseElem, const std::set<size_t>& coarsePatch, Patch& patchDofs) {
- using HierarchicLevelIteratorType = HierarchicLevelIterator<typename LevelContactNetwork::DeformedGridType>;
-
- const auto& grid = coarseContactNetwork_.body(coarseElem.bodyID())->gridView().grid();
- const auto fineLevel = fineContactNetwork_.body(coarseElem.bodyID())->gridView().level();
-
- HierarchicLevelIteratorType endHierIt(grid, coarseElem.element(), HierarchicLevelIteratorType::PositionFlag::end, fineLevel);
- HierarchicLevelIteratorType hierIt(grid, coarseElem.element(), HierarchicLevelIteratorType::PositionFlag::begin, fineLevel);
- for (; hierIt!=endHierIt; ++hierIt) {
- const Element& fineElem = *hierIt;
- addLocalDofs(coarseElem, fineElem, coarsePatch, patchDofs);
- }
- }
-
- void addLocalDofs(const BodyElement& coarseElem, const Element& fineElem, const std::set<size_t>& coarsePatch, Patch& patchDofs) {
- // insert local dofs
- const auto& lfe = cache_.get(fineElem.type());
- const auto& gridView = fineContactNetwork_.body(coarseElem.bodyID())->gridView();
- const auto& indexSet = gridView.indexSet();
-
- for (size_t i=0; i<lfe.localBasis().size(); ++i) {
- auto localIdx = lfe.localCoefficients().localKey(i).subEntity();
- int dofIndex = indexSet.subIndex(fineElem, localIdx, dim);
- patchDofs[dofIndex][0] = false;
- }
-
- const auto coarseLevel = coarseElem.element().level();
-
- // search for parts of the global and inner boundary
- for (const auto& is : intersections(gridView, fineElem)) {
- bool isPatchBoundary = false;
-
- if (is.neighbor()) {
- auto outsideFather = is.outside();
-
- while (outsideFather.level()>coarseLevel) {
- outsideFather = outsideFather.father();
- }
-
- BodyElement outside(coarseElem.bodyID(), outsideFather);
-
- if (!coarsePatch.count(elementIndex(outside))) {
- isPatchBoundary = true;
- }
- } else {
- isPatchBoundary = true;
-
- auto outsideFather = is.outside();
-
- while (outsideFather.level()>coarseLevel) {
- outsideFather = outsideFather.father();
- }
-
- BodyElement outside(coarseElem.bodyID(), outsideFather);
-
- size_t faceIdx = faceIndex(coarseElem.bodyID(), is);
-
- if (neighborElements_.count(faceIdx)) {
- const auto& neighbors = neighborElements_[faceIdx];
- for (size_t i=0; i<neighbors.size(); i++) {
- const auto& neighbor = neighbors[i];
-
- if (coarsePatch.count(elementIndex(outside))) {
- isPatchBoundary = false;
- }
- }
-
- }
- }
-
- if (isPatchBoundary) {
- const auto& localCoefficients = cache_.get(fineElem.type()).localCoefficients();
-
- for (size_t i=0; i<localCoefficients.size(); i++) {
- auto entity = localCoefficients.localKey(i).subEntity();
- auto codim = localCoefficients.localKey(i).codim();
-
- if (containsInsideSubentity(fineElem, is, entity, codim)) {
- int dofIndex = indexSet.subIndex(fineElem, entity, dim);
- patchDofs[dofIndex][0] = true;
- }
- }
- }
- }
- }
-
-
- template <class IndexSet, class Intersection>
- void computeIntersectionDofs(const IndexSet& idxSet, const Intersection& intersection, std::set<size_t>& intersectionDofs, bool inside = true) {
- intersectionDofs.clear();
-
- Element* elem;
- size_t intersectionIndex;
-
- if (inside) {
- elem = &intersection.inside();
- intersectionIndex = intersection.indexInInside();
- } else {
- elem = &intersection.outside();
- intersectionIndex = intersection.indexInOutside();
- }
-
- const int dimElement = Element::dimension;
- const auto& refElement = Dune::ReferenceElements<double, dimElement>::general(elem->type());
-
- for (int i=0; i<refElement.size(intersectionIndex, 1, dimElement); i++) {
- size_t idxInElement = refElement.subEntity(intersectionIndex, 1, i, dimElement);
- size_t globalIdx = idxSet.subIndex(*elem, idxInElement, dimElement);
-
- intersectionDofs.insert(globalIdx);
- }
- }
-};
-#endif
diff --git a/src/tests/gridgluefrictiontest.cc b/src/tests/gridgluefrictiontest.cc
index d31f58d3490ba44c1012a520af030400f1787ef3..6ac022ea9ed5a7867e7f219ac9e422c42facf5f5 100644
--- a/src/tests/gridgluefrictiontest.cc
+++ b/src/tests/gridgluefrictiontest.cc
@@ -334,6 +334,25 @@ int main(int argc, char *argv[]) { try {
print(postprocessed, "postprocessed: ");
+ const auto& contactCouplings = nBodyAssembler.getContactCouplings();
+ const auto& contactCoupling = contactCouplings[0];
+ const auto& glue = *contactCoupling->getGlue();
+
+ size_t isCount = 0;
+ for (const auto& rIs : intersections(glue)) {
+ std::cout << "intersection id: " << isCount
+ << " insideElement: " << gridView0.indexSet().index(rIs.inside())
+ << " outsideElement: " << gridView1.indexSet().index(rIs.outside()) << std::endl;
+ isCount++;
+ }
+
+ for (size_t i=0; i<isCount; i++) {
+ const auto& is = glue.getIntersection(i);
+ std::cout << "intersection id: " << i
+ << " insideElement: " << gridView0.indexSet().index(is.inside())
+ << " outsideElement: " << gridView1.indexSet().index(is.outside()) << std::endl;
+ }
+
/* using DualBasis = ;
DualBasis dualBasis;
BoundaryFunctionalAssembler<DualBasis> bAssembler(dualBasis, nonmortarBoundary);