Skip to content
Snippets Groups Projects
Commit 0a83b8e0 authored by podlesny's avatar podlesny
Browse files

.

parent e7a695c2
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
src/spatial-solving/preconditioners/supportpatchfactory.hh
+ 283
− 294
View file @ 0a83b8e0
...@@ -16,128 +16,77 @@ ...@@ -16,128 +16,77 @@
//#include "../../data-structures/levelcontactnetwork.hh" //#include "../../data-structures/levelcontactnetwork.hh"
#include "hierarchicleveliterator.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> template <class LevelContactNetwork>
class BodyHelper { class NetworkIndexMapper {
private: private:
using FiniteElementCache = Dune::PQkLocalFiniteElementCache<ctype, ctype, dim, 1>; // P1 finite element cache
using GridType = typename LevelContactNetwork::DeformedGridType; using GridType = typename LevelContactNetwork::DeformedGridType;
static const int dim = GridType::dimension;
using Element = typename GridType::template Codim<0>::Entity; 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 >; static const int dim = LevelContactNetwork::dimension; //TODO
using LeafFaceMapper = typename Dune::MultipleCodimMultipleGeomTypeMapper<LeafGridView, mcmgLayout(Codim<1>()) >; static const FiniteElementCache cache_;
bool levelMode_;
std::unique_ptr<LevelGridView> levelGridView_;
LevelElementMapper levelElementMapper_;
LevelFaceMapper levelFaceMapper_;
std::unique_ptr<LeafGridView> leafGridView_; const LevelContactNetwork& levelContactNetwork_;
LeafElementMapper leafElementMapper_; const size_t nBodies_;
LeafFaceMapper leafFaceMapper_;
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> size_t vertexOffSet = 0;
void init(const GV& gv) { size_t faceOffSet = 0;
// init vertexInElements size_t elementOffSet = 0;
vertexInElements_.resize(gv.size(dim));
for (size_t i=0; i<vertexInElements_.size(); i++) {
vertexInElements_[i].resize(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)) { vertexOffSet_[i] = vertexOffSet + gridView.size(dim);
const auto& lfe = cache_.get(elem.type()); vertexOffSet = vertexOffSet_[i];
for (size_t j=0; j<lfe.localBasis().size(); ++j) { faceOffSet_[i] = faceOffSet + gridView.size(1);
auto localIdx = lfe.localCoefficients().localKey(j).subEntity(); faceOffSet = faceOffSet_[i];
int dofIndex = indexSet.subIndex(elem, localIdx, dim);
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) : size_t vertexIndex(const size_t bodyID, const Element& elem, const size_t localVertex) const {
levelMode_(false), const auto& gridView = levelContactNetwork_.body(bodyID)->gridView();
leafGridView_(std::make_unique<LeafGridView>(gridView)),
leafElementMapper_(gridView),
leafFaceMapper_(gridView) {
init(gridView);
}
const auto& gridView() const { auto localIdx = cache_.get(elem.type()).localCoefficients().localKey(localVertex).subEntity();
return (levelMode_ ? levelGridView_ : leafGridView_); return vertexOffSet_[bodyID] + gridView.indexSet().subIndex(elem, localIdx, dim);
} }
const auto& elementMapper() const { template <class Intersection>
return (levelMode_ ? levelElementMapper_ : leafElementMapper_); 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 { size_t faceIndex(const size_t bodyID, const Element& elem, const size_t indexInInside) const {
return (levelMode_ ? levelFaceMapper_ : leafFaceMapper_); const auto& gridView = levelContactNetwork_.body(bodyID)->gridView();
return faceOffSet_[bodyID] + gridView.indexSet().subIndex(elem, indexInInside, 1);
} }
const auto& vertexInElements() const { size_t elementIndex(const size_t bodyID, const Element& elem) const {
return vertexInElements_; 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_; return cache_;
} }
}; };
/** /**
* constructs BitSetVector whose entries are * constructs BitSetVector whose entries are
* true, if vertex is outside of patch or part of patch boundary * true, if vertex is outside of patch or part of patch boundary
...@@ -150,130 +99,139 @@ class SupportPatchFactory ...@@ -150,130 +99,139 @@ class SupportPatchFactory
using Patch = Dune::BitSetVector<1>; using Patch = Dune::BitSetVector<1>;
private: 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 static const int dim = LevelContactNetwork::dimension; //TODO
using ctype = typename LevelContactNetwork::ctype; using ctype = typename LevelContactNetwork::ctype;
using BodyHelper = BodyHelper<LevelContactNetwork>;
//using Vertex = typename GridType::template Codim<dim>::Entity; //using Vertex = typename GridType::template Codim<dim>::Entity;
using Element = typename LevelContactNetwork::DeformedGridType::template Codim<0>::Entity; using Element = typename LevelContactNetwork::DeformedGridType::template Codim<0>::Entity;
using BodyElement = BodyElement<Element>; 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 size_t nBodies_;
const LevelContactNetwork& coarseContactNetwork_; const LevelContactNetwork& coarseContactNetwork_;
const LevelContactNetwork& fineContactNetwork_; 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::set<size_t> interiorContactDofs_;
std::map<size_t, std::vector<size_t>> neighborElementDofs_; // map faceID to neighbor element dofs
std::vector<size_t> vertexOffSet_; const NetworkIndexMapper coarseIndices_;
std::vector<size_t> faceOffSet_; const NetworkIndexMapper fineIndices_;
std::vector<size_t> elementOffSet_;
FiniteElementCache cache_;
public: public:
SupportPatchFactory(const LevelContactNetwork& coarseContactNetwork, const LevelContactNetwork& fineContactNetwork) : SupportPatchFactory(const LevelContactNetwork& coarseContactNetwork, const LevelContactNetwork& fineContactNetwork) :
nBodies_(coarseContactNetwork_.nBodies()), nBodies_(coarseContactNetwork_.nBodies()),
coarseContactNetwork_(coarseContactNetwork), coarseContactNetwork_(coarseContactNetwork),
fineContactNetwork_(fineContactNetwork), fineContactNetwork_(fineContactNetwork),
bodyHelper_(nBodies_), coarseIndices_(coarseContactNetwork_),
vertexOffSet_(nBodies_), fineIndices_(fineContactNetwork_) {
faceOffSet_(nBodies_),
elementOffSet_(nBodies_){
assert(nBodies_ == fineContactNetwork_.nBodies()); assert(nBodies_ == fineContactNetwork_.nBodies());
size_t vertexOffSet = 0; setup(coarseContactNetwork_, coarseIndices_, coarseIntersectionMapper_);
size_t faceOffSet = 0; setup(fineContactNetwork_, fineIndices_, fineIntersectionMapper_);
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 neighborFaceDofs_
const auto& contactCouplings = coarseContactNetwork_.nBodyAssembler().getContactCouplings();
for (size_t i=0; i<contactCouplings.size(); i++) { for (size_t i=0; i<contactCouplings.size(); i++) {
const auto& coupling = nBodyAssembler.getCoupling(i); const auto& coupling = nBodyAssembler.getCoupling(i);
const auto& contactCoupling = contactCouplings[i]; const auto& glue = *contactCouplings[i].getGlue();
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); const size_t nmBody = coupling.gridIdx_[0];
neighborElements_[mFaceIdx].emplace_back(nmBodyElement); const size_t mBody = coupling.gridIdx_[1];
std::set<size_t> dofs; for (const auto& rIs : intersections(glue)) {
computeIntersectionDofs(nmIndexSet, rIs, dofs); size_t nmFaceIdx = coarseIndices_.faceIndex(nmBody, rIs.inside(), rIs.indexInInside());
neighborElementDofs_[nmFaceIdx].insert(neighborElementDofs_[nmFaceIdx].end(), dofs.begin(), dofs.end()); size_t mFaceIdx = coarseIndices_.faceIndex(mBody, rIs.outside(), rIs.indexInOutside());
computeIntersectionDofs(mIndexSet, rIs, dofs, false); std::array<std::set<size_t>, 2> dofs;
neighborElementDofs_[mFaceIdx].insert(neighborElementDofs_[mFaceIdx].end(), dofs.begin(), dofs.end()); 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 // neighborElementDofs_ contains dofs of connected intersections that are neighbors to a face given by faceID,
// task: remove boundary dofs and redundant multiples of interior dofs // boundary dofs are contained once, interior dofs multiple times
for (auto& it : neighborElementDofs_) { // task: remove boundary dofs and redundant multiples of interior dofs
auto& dofs = it.second; // TODO: only works in 2D
for (auto& it : neighborElementDofs) {
std::sort(dofs.begin(), dofs.end()); auto& dofs = it.second;
auto& dof = dofs.begin();
while (dof != dofs.end()) { std::sort(dofs.begin(), dofs.end());
auto next = dof; auto& dof = dofs.begin();
if (*dof != *(++next)) while (dof != dofs.end()) {
dof = dofs.erase(dof); 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) { void build(const size_t bodyID, const Element& coarseElement, const size_t localVertex, Patch& patchDofs, const int patchDepth = 0) {
BodyElement seedElement(bodyID, coarseElement); 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 // construct coarse patch
std::vector<BodyElement> patchElements; std::vector<BodyElement> patchElements;
patchElements.emplace_back(seedElement); patchElements.emplace_back(seedElement);
std::set<size_t> visitedElements; std::set<size_t> visitedElements;
visitedElements.insert(elementIndex(seedElement)); visitedElements.insert(coarseIndices_.elementIndex(seedElement.bodyID, seedElement.element));
std::set<size_t> patchVertices; std::set<size_t> patchVertices;
patchVertices.insert(vertexIndex(bodyID, coarseElement, localVertex)); patchVertices.insert(coarseIndices_.vertexIndex(bodyID, coarseElement, localVertex));
std::queue<BodyElement> patchSeeds; std::queue<BodyElement> patchSeeds;
patchSeeds.push(seedElement); patchSeeds.push(seedElement);
...@@ -286,36 +244,37 @@ class SupportPatchFactory ...@@ -286,36 +244,37 @@ class SupportPatchFactory
const auto& patchSeed = patchSeeds.front(); const auto& patchSeed = patchSeeds.front();
patchSeeds.pop(); patchSeeds.pop();
size_t elemIdx = elementIndex(patchSeed); size_t elemIdx = coarseIndices_.elementIndex(patchSeed.bodyID, patchSeed.element);
if (visitedElements.count(elemIdx)) if (visitedElements.count(elemIdx))
continue; continue;
visitedElements.insert(elemIdx); visitedElements.insert(elemIdx);
size_t bodyIdx = patchSeed.bodyID(); size_t bodyIdx = patchSeed.bodyID;
const auto& elem = patchSeed.element(); const auto& elem = patchSeed.element;
const auto& gridView = bodyHelper_[bodyIdx]->gridView(); const auto& gridView = coarseContactNetwork_.body(bodyIdx)->gridView();
for (auto& it : intersections(gridView, elem)) { for (auto& it : intersections(gridView, elem)) {
if (isPatchIntersection(gridView.indexSet(), it, patchVertices, newPatchVertices)) { if (isPatchIntersection(it, bodyIdx, patchVertices, newPatchVertices)) {
if (it.neighbor()) { if (it.neighbor()) {
BodyElement neighbor(bodyIdx, it.outside()); BodyElement neighbor(bodyIdx, it.outside());
if (visitedElements.count(elementIndex(neighbor))) if (visitedElements.count(coarseIndices_.elementIndex(neighbor.bodyID, neighbor.element)))
continue; continue;
patchElements.emplace_back(neighbor); patchElements.emplace_back(neighbor);
patchSeeds.push(neighbor); patchSeeds.push(neighbor);
} else { } else {
size_t faceIdx = faceIndex(bodyIdx, it); size_t faceIdx = coarseIndices_.faceIndex(bodyIdx, it);
if (neighborElements_.count(faceIdx)) { if (coarseIntersectionMapper_.count(faceIdx)) {
const auto& neighbors = neighborElements_[faceIdx]; const auto& intersections = coarseIntersectionMapper_[faceIdx];
for (size_t i=0; i<neighbors.size(); i++) { for (size_t i=0; i<intersections.size(); i++) {
const auto& neighbor = neighbors[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; continue;
patchVertices.insert(neighborElementDofs_[faceIdx].begin(), neighborElementDofs_[faceIdx].end()); patchVertices.insert(neighborElementDofs_[faceIdx].begin(), neighborElementDofs_[faceIdx].end());
...@@ -328,19 +287,20 @@ class SupportPatchFactory ...@@ -328,19 +287,20 @@ class SupportPatchFactory
if (it.neighbor()) { if (it.neighbor()) {
BodyElement neighbor(bodyIdx, it.outside()); BodyElement neighbor(bodyIdx, it.outside());
if (visitedElements.count(elementIndex(neighbor))) if (visitedElements.count(coarseIndices_.elementIndex(neighbor.bodyID, neighbor.element)))
continue; continue;
nextSeeds.emplace_back(neighbor); nextSeeds.emplace_back(neighbor);
} else { } else {
size_t faceIdx = faceIndex(bodyIdx, it); size_t faceIdx = coarseIndices_.faceIndex(bodyIdx, it);
if (neighborElements_.count(faceIdx)) { if (coarseIntersectionMapper_.count(faceIdx)) {
const auto& neighbors = neighborElements_[faceIdx]; const auto& intersections = coarseIntersectionMapper_[faceIdx];
for (size_t i=0; i<neighbors.size(); i++) { for (size_t i=0; i<intersections.size(); i++) {
const auto& neighbor = neighbors[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; continue;
newPatchVertices.insert(neighborElementDofs_[faceIdx].begin(), neighborElementDofs_[faceIdx].end()); newPatchVertices.insert(neighborElementDofs_[faceIdx].begin(), neighborElementDofs_[faceIdx].end());
...@@ -359,42 +319,50 @@ class SupportPatchFactory ...@@ -359,42 +319,50 @@ class SupportPatchFactory
} }
// construct fine patch // construct fine patch
using HierarchicLevelIteratorType = HierarchicLevelIterator<typename LevelContactNetwork::DeformedGridType>;
std::set<size_t> patchBoundary;
for (size_t i=0; i<patchElements.size(); ++i) { for (size_t i=0; i<patchElements.size(); ++i) {
addFinePatchElements(patchElements[i], visitedElements, patchDofs); const auto& coarseElem = patchElements[i];
}
}
private: const auto& grid = coarseContactNetwork_.body(coarseElem.bodyID)->gridView().grid();
template <class IndexSet, class Intersection> const auto fineLevel = fineContactNetwork_.body(coarseElem.bodyID)->gridView().level();
bool isPatchIntersection(const IndexSet& indexSet, const Intersection& is, const std::set<size_t>& patchVertices, std::set<size_t>& newPatchVertices) {
std::set<size_t> intersectionDofs; // add fine patch elements
computeIntersectionDofs(indexSet, is, intersectionDofs); HierarchicLevelIteratorType endHierIt(grid, coarseElem.element, HierarchicLevelIteratorType::PositionFlag::end, fineLevel);
HierarchicLevelIteratorType hierIt(grid, coarseElem.element, HierarchicLevelIteratorType::PositionFlag::begin, fineLevel);
for (auto& dof : intersectionDofs) { for (; hierIt!=endHierIt; ++hierIt) {
if (patchVertices.count(dof)) { const Element& fineElem = *hierIt;
newPatchVertices.insert(intersectionDofs.begin(), intersectionDofs.end()); addLocalDofs(coarseElem, fineElem, visitedElements, patchBoundary, patchDofs);
return true;
} }
} }
return false;
} }
size_t vertexIndex(const size_t bodyID, const Element& elem, const size_t localVertex) const { private:
const auto& gridView = bodyHelper_[bodyID]->gridView(); 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(); for (size_t i=0; i<contactCouplings.size(); i++) {
return vertexOffSet_[bodyID] + gridView.indexSet().subIndex(elem, localIdx, dim); const auto& coupling = nBodyAssembler.getCoupling(i);
} const auto& glue = *contactCouplings[i].getGlue();
template <class Intersection> const size_t nmBody = coupling.gridIdx_[0];
size_t faceIndex(const size_t bodyID, const Intersection& is) const { const size_t mBody = coupling.gridIdx_[1];
const auto& gridView = bodyHelper_[bodyID]->gridView();
return faceOffSet_[bodyID] + gridView.indexSet().subIndex(is.inside(), is.indexInInside(), 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 { RemoteIntersection mIntersection(nmBody, glue, rIsID, true);
const auto& gridView = bodyHelper_[bodyElem.bodyID()]->gridView(); faceToRIntersections[mFaceIdx].emplace_back(mIntersection);
return elementOffSet_[bodyElem.bodyID()] + gridView.indexSet().index(bodyElem.element());
rIsID++;
}
}
} }
template <class Intersection> template <class Intersection>
...@@ -402,112 +370,133 @@ class SupportPatchFactory ...@@ -402,112 +370,133 @@ class SupportPatchFactory
return ReferenceElementHelper<double, dim>::subEntityContainsSubEntity(elem.type(), intersection.indexInInside(), 1, subEntity, 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) { void addLocalDofs(const BodyElement& coarseElem, const Element& fineElem, const std::set<size_t>& coarsePatch, std::set<size_t>& patchBoundary, 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 // insert local dofs
const auto& lfe = cache_.get(fineElem.type()); 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(); 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 // search for parts of the global and inner boundary
for (const auto& is : intersections(gridView, fineElem)) { for (const auto& is : intersections(gridView, fineElem)) {
bool isPatchBoundary = false;
if (is.neighbor()) { if (is.neighbor()) {
auto outsideFather = is.outside(); std::set<size_t> isDofs;
intersectionDofs(fineIndices_, is, bodyID, isDofs);
while (outsideFather.level()>coarseLevel) {
outsideFather = outsideFather.father();
}
BodyElement outside(coarseElem.bodyID(), outsideFather);
if (!coarsePatch.count(elementIndex(outside))) { auto outsideFather = coarseFather(is.outside(), coarseLevel);
isPatchBoundary = true;
if (!coarsePatch.count(coarseIndices_.elementIndex(coarseElem.bodyID, outsideFather))) {
addToPatchBoundary(isDofs, patchBoundary, patchDofs);
} else {
addToPatch(isDofs, patchBoundary, patchDofs);
} }
} else { } 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) { std::set<size_t> rIsDofs;
outsideFather = outsideFather.father(); 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 (!coarsePatch.count(coarseIndices_.elementIndex(outsideID, outsideFather))) {
addToPatchBoundary(rIsDofs, patchBoundary, patchDofs);
if (neighborElements_.count(faceIdx)) { } else {
const auto& outside = neighbors[0]; addToPatch(rIsDofs, patchBoundary, patchDofs);
}
if (coarsePatch.count(elementIndex(outside))) { } else {
isPatchBoundary = false; 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++) { template <class RIntersection>
auto entity = localCoefficients.localKey(i).subEntity(); void remoteIntersectionDofs(const NetworkIndexMapper& indices, const RIntersection& is, const size_t insideID, const size_t outsideID, std::array<std::set<size_t>, 2>& dofs) {
auto codim = localCoefficients.localKey(i).codim(); dofs.clear();
if (containsInsideSubentity(fineElem, is, entity, codim)) { const auto& inside = is.inside();
int dofIndex = indexSet.subIndex(fineElem, entity, dim); const auto& insideGeo = inside.geometry();
patchDofs[dofIndex][0] = true; 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> const auto& outside = is.outside();
void computeIntersectionDofs(const IndexSet& idxSet, const Intersection& intersection, std::set<size_t>& intersectionDofs, bool inside = true) { const auto& outsideGeo = outside.geometry();
intersectionDofs.clear(); const auto& outsideRefElement = Dune::ReferenceElements<ctype, dim>::general(outside.type());
Element* elem; for (int i=0; i<outsideRefElement.size(dim); i++) {
size_t intersectionIndex; const auto& localCorner = isGeo.local(outsideGeo.corner(i));
if (inside) { if (isRefElement.checkInside(localCorner)) {
elem = &intersection.inside(); dofs[1].insert(indices.vertexIndex(outsideID, outside, i));
intersectionIndex = intersection.indexInInside(); }
} else {
elem = &intersection.outside();
intersectionIndex = intersection.indexInOutside();
} }
}
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 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++) { for (int i=0; i<refElement.size(intersectionIndex, 1, dimElement); i++) {
size_t idxInElement = refElement.subEntity(intersectionIndex, 1, i, dimElement); size_t idxInElement = refElement.subEntity(intersectionIndex, 1, i, dimElement);
size_t globalIdx = idxSet.subIndex(*elem, idxInElement, dimElement); dofs.insert(indices.vertexIndex(bodyID, elem, idxInElement));
intersectionDofs.insert(globalIdx);
} }
} }
}; };
......
src/spatial-solving/preconditioners/supportpatchfactory.hh.autosave deleted 100644 → 0
+ 0
− 518
View file @ e7a695c2
#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
src/tests/gridgluefrictiontest.cc
+ 19
− 0
View file @ 0a83b8e0
...@@ -334,6 +334,25 @@ int main(int argc, char *argv[]) { try { ...@@ -334,6 +334,25 @@ int main(int argc, char *argv[]) { try {
print(postprocessed, "postprocessed: "); 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 = ; /* using DualBasis = ;
DualBasis dualBasis; DualBasis dualBasis;
BoundaryFunctionalAssembler<DualBasis> bAssembler(dualBasis, nonmortarBoundary); BoundaryFunctionalAssembler<DualBasis> bAssembler(dualBasis, nonmortarBoundary);
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment