diff --git a/src/multi-body-problem.cc b/src/multi-body-problem.cc
index 6ab6c3dd4864db9acc42226d07046974396ab0e0..aa662b5e41d01164eb6b4ed050e9395376b63c2b 100644
--- a/src/multi-body-problem.cc
+++ b/src/multi-body-problem.cc
@@ -185,11 +185,36 @@ int main(int argc, char *argv[]) {
const auto & coarseContactNetwork = *contactNetwork.level(0);
const auto & fineContactNetwork = *contactNetwork.level(1);
SupportPatchFactory<typename ContactNetwork::LevelContactNetwork> supportPatchFactory(coarseContactNetwork, fineContactNetwork);
- size_t bodyID = 0;
+ size_t bodyID = 1;
size_t patchDepth = 0;
- writeToVTK(fineContactNetwork.body(bodyID)->gridView(), "", "body_" + std::to_string(bodyID) + "_fine");
- writeToVTK(coarseContactNetwork.body(bodyID)->gridView(), "", "body_" + std::to_string(bodyID) + "_coarse");
+ std::cout << std::endl;
+
+ // print coarse dofs
+ for (size_t i=0; i<bodyCount; i++) {
+ std::cout << "Coarse dofs body " << i << std::endl;
+ const auto& gv = coarseContactNetwork.body(i)->gridView();
+ printDofLocation(gv);
+
+ ScalarVector dofs(gv.size(dims));
+ for (size_t j=0; j<dofs.size(); j++) {
+ dofs[j] = j;
+ }
+ writeToVTK(gv, dofs, "", "body_" + std::to_string(i) + "_coarse");
+ }
+
+ // print fine dofs
+ for (size_t i=0; i<bodyCount; i++) {
+ std::cout << "Fine dofs body " << i << std::endl;
+ const auto& gv = fineContactNetwork.body(i)->gridView();
+ printDofLocation(gv);
+
+ ScalarVector dofs(gv.size(dims));
+ for (size_t j=0; j<dofs.size(); j++) {
+ dofs[j] = j;
+ }
+ writeToVTK(gv, dofs, "", "body_" + std::to_string(i) + "_fine");
+ }
using Patch = typename SupportPatchFactory<typename ContactNetwork::LevelContactNetwork>::Patch;
Patch patch;
@@ -230,10 +255,10 @@ int main(int argc, char *argv[]) {
writeToVTK(gv, patchVec, "", "patch_" + std::to_string(globalIdx) + "_body_" + std::to_string(j));
}
}
-
-
}
}
+
+ std::cout << "Done! :)" << std::endl;
return 1;
//printMortarBasis<Vector>(contactNetwork.nBodyAssembler());
diff --git a/src/spatial-solving/preconditioners/supportpatchfactory.hh b/src/spatial-solving/preconditioners/supportpatchfactory.hh
index 75db9ad29689aaac6250bfa567b13f3da6febb9c..a7274816d1b90f26900b9ca9e90fafa6ef1a6760 100644
--- a/src/spatial-solving/preconditioners/supportpatchfactory.hh
+++ b/src/spatial-solving/preconditioners/supportpatchfactory.hh
@@ -17,6 +17,8 @@
//#include "../../data-structures/levelcontactnetwork.hh"
#include "hierarchicleveliterator.hh"
+#include "../../utils/debugutils.hh"
+
template <class LevelContactNetwork>
class NetworkIndexMapper {
private:
@@ -50,14 +52,14 @@ class NetworkIndexMapper {
for (size_t i=0; i<nBodies_; i++) {
const auto& gridView = levelContactNetwork_.body(i)->gridView();
- vertexOffSet_[i] = vertexOffSet + gridView.size(dim);
- vertexOffSet = vertexOffSet_[i];
+ vertexOffSet_[i] = vertexOffSet;
+ vertexOffSet = vertexOffSet_[i] + gridView.size(dim);
- faceOffSet_[i] = faceOffSet + gridView.size(1);
- faceOffSet = faceOffSet_[i];
+ faceOffSet_[i] = faceOffSet;
+ faceOffSet = faceOffSet_[i] + gridView.size(1);
- elementOffSet_[i] = elementOffSet + gridView.size(0);
- elementOffSet = elementOffSet_[i];
+ elementOffSet_[i] = elementOffSet;
+ elementOffSet = elementOffSet_[i] + gridView.size(0);
}
}
@@ -108,16 +110,6 @@ class SupportPatchFactory
using Element = typename LevelContactNetwork::GridView::template Codim<0>::Entity;
- 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
size_t contactCouplingID;
@@ -132,9 +124,28 @@ class SupportPatchFactory
flip = _flip;
}
- template <class GridGlue>
- const auto& get(const GridGlue& glue) const {
- return glue.getIntersection(intersectionID);
+ auto get(const LevelContactNetwork& levelContactNetwork) const {
+ return std::move(levelContactNetwork.glue(contactCouplingID)->getIntersection(intersectionID));
+ }
+ };
+
+ struct BodyElement {
+ size_t bodyID;
+ Element element;
+
+ void set(const size_t _bodyID, const Element& _elem) {
+ bodyID = _bodyID;
+ element = _elem;
+ }
+
+ void set(const LevelContactNetwork& levelContactNetwork, const RemoteIntersection& rIs) {
+ const auto& is = rIs.get(levelContactNetwork);
+
+ if (rIs.flip) {
+ set(rIs.bodyID, is.inside());
+ } else {
+ set(rIs.bodyID, is.outside());
+ }
}
};
@@ -187,11 +198,57 @@ class SupportPatchFactory
std::vector<std::set<size_t>> 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());
+
+ if (dofs[1].size()>0) {
+ neighborFaceDofs_[nmFaceIdx].insert(neighborFaceDofs_[nmFaceIdx].end(), dofs[1].begin(), dofs[1].end());
+ }
+
+ if (dofs[0].size()>0) {
+ neighborFaceDofs_[mFaceIdx].insert(neighborFaceDofs_[mFaceIdx].end(), dofs[0].begin(), dofs[0].end());
+ }
}
+ }
+
+ for (size_t i=0; i<nBodies_; i++) {
+ std::cout << "Coarse body" << i << ":" << std::endl;
+
+ for (const auto& e : elements(coarseContactNetwork_.body(i)->gridView())) {
+ std::cout << "elemID: " << coarseIndices_.elementIndex(i, e) << std::endl;
+ std::cout << "vertexIDs: ";
+ const int dimElement = Element::dimension;
+ const auto& refElement = Dune::ReferenceElements<double, dimElement>::general(e.type());
+
+ for (int j=0; j<refElement.size(dim); j++) {
+ std::cout << coarseIndices_.vertexIndex(i, e, j) << " ";
+ }
+ std::cout << std::endl;
+
+ std::cout << "intersectionIDs: ";
+ for (const auto& is : intersections(coarseContactNetwork_.body(i)->gridView(), e)) {
+ std::cout << coarseIndices_.faceIndex(i, e, is.indexInInside()) << " (";
+
+ std::set<size_t> dofs;
+ intersectionDofs(coarseIndices_, is, i, dofs);
+ for (auto& d : dofs) {
+ std::cout << d << " ";
+ }
+ std::cout << "); ";
+
+ }
+ std::cout << std::endl << "--------------" << std::endl;
}
+ }
+
+ std::cout << std::endl;
+ for (auto& kv : neighborFaceDofs_) {
+ std::cout << "faceID: " << kv.first << " dofs: ";
+ const auto& dofs = kv.second;
+ for (auto& d : dofs) {
+ std::cout << d << " ";
+ }
+ std::cout << std::endl;
+ }
/*
// neighborElementDofs_ contains dofs of connected intersections that are neighbors to a face given by faceID,
// boundary dofs are contained once, interior dofs multiple times
@@ -215,6 +272,8 @@ class SupportPatchFactory
void build(const size_t bodyID, const Element& coarseElement, const size_t localVertex, Patch& patchDofs, const size_t patchDepth = 0) {
std::cout << "Building SupportPatch... ";
+ std::cout << "elemID: " << coarseIndices_.elementIndex(bodyID, coarseElement) << " vertexID: " << coarseIndices_.vertexIndex(bodyID, coarseElement, localVertex) << std::endl;
+
patchDofs.resize(nVertices_);
patchDofs.setAll();
@@ -239,7 +298,7 @@ class SupportPatchFactory
patchElements.emplace_back(seedElement);
std::set<size_t> visitedElements;
- visitedElements.insert(coarseIndices_.elementIndex(seedElement.bodyID, seedElement.element));
+ //visitedElements.insert(coarseIndices_.elementIndex(seedElement.bodyID, seedElement.element));
std::set<size_t> patchVertices;
patchVertices.insert(coarseIndices_.vertexIndex(bodyID, coarseElement, localVertex));
@@ -272,7 +331,11 @@ class SupportPatchFactory
BodyElement neighbor;
neighbor.set(bodyIdx, it.outside());
- if (visitedElements.count(coarseIndices_.elementIndex(neighbor.bodyID, neighbor.element)))
+ size_t neighborIdx = coarseIndices_.elementIndex(neighbor.bodyID, neighbor.element);
+
+ std::cout << "elementID: " << neighborIdx << std::endl;
+
+ if (visitedElements.count(neighborIdx))
continue;
patchElements.emplace_back(neighbor);
@@ -284,15 +347,20 @@ class SupportPatchFactory
const auto& intersections = coarseIntersectionMapper_[faceIdx];
for (size_t i=0; i<intersections.size(); i++) {
BodyElement neighbor;
- const auto& glue = *coarseContactNetwork_.glue(intersections[i].contactCouplingID);
- neighbor.set(intersections[i].bodyID, intersections[i].get(glue).outside());
+ neighbor.set(coarseContactNetwork_, intersections[i]);
+
+ size_t neighborIdx = coarseIndices_.elementIndex(neighbor.bodyID, neighbor.element);
- if (visitedElements.count(coarseIndices_.elementIndex(neighbor.bodyID, neighbor.element)))
+ if (visitedElements.count(neighborIdx))
continue;
+ std::cout << "elementID: " << coarseIndices_.elementIndex(neighbor.bodyID, neighbor.element) << std::endl;
+
patchVertices.insert(neighborFaceDofs_[faceIdx].begin(), neighborFaceDofs_[faceIdx].end());
patchElements.emplace_back(neighbor);
patchSeeds.push(neighbor);
+
+ print(patchVertices, "patchVertices: ");
}
}
}
@@ -301,7 +369,9 @@ class SupportPatchFactory
BodyElement neighbor;
neighbor.set(bodyIdx, it.outside());
- if (visitedElements.count(coarseIndices_.elementIndex(neighbor.bodyID, neighbor.element)))
+ size_t neighborIdx = coarseIndices_.elementIndex(neighbor.bodyID, neighbor.element);
+
+ if (visitedElements.count(neighborIdx))
continue;
nextSeeds.emplace_back(neighbor);
@@ -312,10 +382,11 @@ class SupportPatchFactory
const auto& intersections = coarseIntersectionMapper_[faceIdx];
for (size_t i=0; i<intersections.size(); i++) {
BodyElement neighbor;
- const auto& glue = *coarseContactNetwork_.glue(intersections[i].contactCouplingID);
- neighbor.set(intersections[i].bodyID, intersections[i].get(glue).outside());
+ neighbor.set(coarseContactNetwork_, intersections[i]);
+
+ size_t neighborIdx = coarseIndices_.elementIndex(neighbor.bodyID, neighbor.element);
- if (visitedElements.count(coarseIndices_.elementIndex(neighbor.bodyID, neighbor.element)))
+ if (visitedElements.count(neighborIdx))
continue;
newPatchVertices.insert(neighborFaceDofs_[faceIdx].begin(), neighborFaceDofs_[faceIdx].end());
@@ -331,14 +402,24 @@ class SupportPatchFactory
patchSeeds.push(nextSeeds[i]);
}
patchVertices.insert(newPatchVertices.begin(), newPatchVertices.end());
+ print(patchVertices, "patchVertices: ");
}
+ std::cout << "constructing fine patch... " << std::endl;
+
// construct fine patch
using HierarchicLevelIteratorType = HierarchicLevelIterator<typename LevelContactNetwork::GridType>;
+ std::set<size_t> visited;
std::set<size_t> patchBoundary;
for (size_t i=0; i<patchElements.size(); ++i) {
const auto& coarseElem = patchElements[i];
+ size_t elemIdx = coarseIndices_.elementIndex(coarseElem.bodyID, coarseElem.element);
+ if (visited.count(elemIdx))
+ continue;
+
+ visited.insert(elemIdx);
+
const auto& grid = coarseContactNetwork_.body(coarseElem.bodyID)->gridView().grid();
const auto fineLevel = fineContactNetwork_.body(coarseElem.bodyID)->level();
@@ -370,15 +451,18 @@ class SupportPatchFactory
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;
nmIntersection.set(mBody, i, rIsID, false);
faceToRIntersections[nmFaceIdx].emplace_back(nmIntersection);
- RemoteIntersection mIntersection;
- mIntersection.set(nmBody, i, rIsID, true);
- faceToRIntersections[mFaceIdx].emplace_back(mIntersection);
+ if (rIs.neighbor()) {
+ size_t mFaceIdx = indices.faceIndex(mBody, rIs.outside(), rIs.indexInOutside());
+
+ RemoteIntersection mIntersection;
+ mIntersection.set(nmBody, i, rIsID, true);
+ faceToRIntersections[mFaceIdx].emplace_back(mIntersection);
+ }
rIsID++;
}
@@ -393,8 +477,6 @@ class SupportPatchFactory
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& gridView = fineContactNetwork_.body(coarseElem.bodyID)->gridView();
- const auto& indexSet = gridView.indexSet();
-
const size_t bodyID = coarseElem.bodyID;
const auto coarseLevel = coarseElem.element.level();
@@ -413,28 +495,34 @@ class SupportPatchFactory
addToPatch(isDofs, patchBoundary, patchDofs);
}
} else {
- /* size_t faceIdx = fineIndices_.faceIndex(bodyID, is);
+ size_t faceIdx = fineIndices_.faceIndex(bodyID, is);
if (fineIntersectionMapper_.count(faceIdx)) {
const auto& intersections = fineIntersectionMapper_[faceIdx];
for (size_t i=0; i<intersections.size(); i++) {
const auto& intersection = intersections[i];
- const auto& glue = *fineContactNetwork_.glue(intersection.contactCouplingID);
- const auto& rIs = intersection.get(glue);
+ const auto& rIs = intersection.get(fineContactNetwork_);
const size_t outsideID = intersection.bodyID;
std::vector<std::set<size_t>> rIsDofs;
- remoteIntersectionDofs(fineIndices_, rIs, bodyID, outsideID, rIsDofs);
+ if (intersection.flip) {
+ remoteIntersectionDofs(fineIndices_, rIs, outsideID, bodyID, rIsDofs);
+ } else {
+ remoteIntersectionDofs(fineIndices_, rIs, bodyID, outsideID, rIsDofs);
+ }
if (rIs.neighbor()) {
Element outsideFather;
+ const size_t outsideLevel = coarseContactNetwork_.body(outsideID)->level();
+
if (intersection.flip)
- outsideFather = coarseFather(rIs.inside(), fineContactNetwork_.body(outsideID)->level());
+ outsideFather = coarseFather(rIs.inside(), outsideLevel);
else
- outsideFather = coarseFather(rIs.outside(), fineContactNetwork_.body(outsideID)->level());
+ outsideFather = coarseFather(rIs.outside(), outsideLevel);
std::set<size_t> totalRIsDofs(rIsDofs[0]);
totalRIsDofs.insert(rIsDofs[1].begin(), rIsDofs[1].end());
+
if (!coarsePatch.count(coarseIndices_.elementIndex(outsideID, outsideFather))) {
addToPatchBoundary(totalRIsDofs, patchBoundary, patchDofs);
} else {
@@ -449,7 +537,7 @@ class SupportPatchFactory
intersectionDofs(fineIndices_, is, bodyID, isDofs);
addToPatchBoundary(isDofs, patchBoundary, patchDofs);
- } */
+ }
}
}
}
@@ -489,12 +577,15 @@ class SupportPatchFactory
const auto& inside = is.inside();
const auto& insideGeo = inside.geometry();
const auto& insideRefElement = Dune::ReferenceElements<ctype, dim>::general(inside.type());
+ const int dimElement = Dune::ReferenceElements<ctype, dim>::dimension;
- for (int i=0; i<insideRefElement.size(dim); i++) {
- const auto& localCorner = isGeo.local(insideGeo.corner(i));
+ for (int i=0; i<insideRefElement.size(is.indexInInside(), 1, dimElement); i++) {
+ size_t idxInElement = insideRefElement.subEntity(is.indexInInside(), 1, i, dimElement);
+
+ const auto& localCorner = isGeo.local(insideGeo.corner(idxInElement));
if (isRefElement.checkInside(localCorner)) {
- dofs[0].insert(indices.vertexIndex(insideID, inside, i));
+ dofs[0].insert(indices.vertexIndex(insideID, inside, idxInElement));
}
}
@@ -503,11 +594,13 @@ class SupportPatchFactory
const auto& outsideGeo = outside.geometry();
const auto& outsideRefElement = Dune::ReferenceElements<ctype, dim>::general(outside.type());
- for (int i=0; i<outsideRefElement.size(dim); i++) {
- const auto& localCorner = isGeo.local(outsideGeo.corner(i));
+ for (int i=0; i<outsideRefElement.size(is.indexInOutside(), 1, dimElement); i++) {
+ size_t idxInElement = outsideRefElement.subEntity(is.indexInOutside(), 1, i, dimElement);
+
+ const auto& localCorner = isGeo.local(outsideGeo.corner(idxInElement));
if (isRefElement.checkInside(localCorner)) {
- dofs[1].insert(indices.vertexIndex(outsideID, outside, i));
+ dofs[1].insert(indices.vertexIndex(outsideID, outside, idxInElement));
}
}
}
diff --git a/src/utils/debugutils.hh b/src/utils/debugutils.hh
index 32c923f46fc17722f1b77df367b56a59eb0cc201..4cd6024a74fd2904dab51b878371e24b1a01e3bb 100644
--- a/src/utils/debugutils.hh
+++ b/src/utils/debugutils.hh
@@ -129,18 +129,18 @@
std::cout << x[i] << " ";
}
std::cout << std::endl << std::endl;
- }
+ }*/
- void print(const std::set<size_t>& x, std::string message){
+ template <class T>
+ void print(const std::set<T>& x, std::string message){
std::cout << message << std::endl;
- std::set<size_t>::iterator dofIt = x.begin();
- std::set<size_t>::iterator dofEndIt = x.end();
- for (; dofIt != dofEndIt; ++dofIt) {
- std::cout << *dofIt << " ";
+
+ for (const auto& c : x) {
+ std::cout << c << " ";
}
std::cout << std::endl << std::endl;
}
-
+/*
void print(const std::set<int>& x, std::string message){
std::cout << message << std::endl;
std::set<int>::iterator dofIt = x.begin();