From e544c081ea41bb3ab833d1626856075fd4c7aa14 Mon Sep 17 00:00:00 2001
From: podlesny <podlesny@mi.fu-berlin.de>
Date: Thu, 27 Jun 2019 17:27:20 +0200
Subject: [PATCH] supportpatchfactory works
---
.../levelpatchpreconditioner.hh | 94 +++++++++----------
.../multilevelpatchpreconditioner.hh | 12 +--
.../preconditioners/supportpatchfactory.hh | 69 +++++---------
3 files changed, 74 insertions(+), 101 deletions(-)
diff --git a/src/spatial-solving/preconditioners/levelpatchpreconditioner.hh b/src/spatial-solving/preconditioners/levelpatchpreconditioner.hh
index d905e4a3..1c8257c1 100644
--- a/src/spatial-solving/preconditioners/levelpatchpreconditioner.hh
+++ b/src/spatial-solving/preconditioners/levelpatchpreconditioner.hh
@@ -19,33 +19,33 @@
#include <dune/fufem/functiontools/boundarydofs.hh>
-template <class PatchSolver, class MatrixType, class VectorType>
+template <class LevelContactNetwork, class PatchSolver, class MatrixType, class VectorType>
class LevelPatchPreconditioner : public LinearIterationStep<MatrixType, VectorType> {
public:
enum Mode {additive, multiplicative};
+ static const int dim = LevelContactNetwork::dim;
private:
using Base = LinearIterationStep<MatrixType, VectorType>;
- using PatchFactory = SupportPatchFactory<GridView, GridView>;
- using Patch = typename SupportPatchFactory::Patch;
+ using ctype = typename LevelContactNetwork::ctype;
- typedef typename BasisType::GridView GridView;
- typedef typename GridView::Grid GridType;
-
- static const int dim = GridType::dimension;
+ using PatchFactory = SupportPatchFactory<LevelContactNetwork>;
+ using Patch = typename PatchFactory::Patch;
const Mode mode_;
- const LevelContactNetwork<GridView, dim>& levelContactNetwork_;
+ const LevelContactNetwork& levelContactNetwork_;
+ const LevelContactNetwork& fineContactNetwork_;
+
+ const int level_;
+ PatchFactory patchFactory_;
std::vector<Patch> patches_;
PatchSolver patchSolver_;
- const GridType& grid_;
- const int level_;
- const BasisType basis_;
+
size_t patchDepth_;
@@ -54,30 +54,28 @@ class LevelPatchPreconditioner : public LinearIterationStep<MatrixType, VectorTy
public:
- // for each coarse patch given by patchLevelBasis: set up local problem, compute local correction
- LevelPatchPreconditioner(const LevelContactNetwork<GridView, dim>& levelContactNetwork,
- const BasisType& patchLevelBasis,
- const LocalAssembler& localAssembler,
- const std::vector<std::shared_ptr<LocalInterfaceAssembler>>& localInterfaceAssemblers,
+ // for each coarse patch given by levelContactNetwork: set up local problem, compute local correction
+ LevelPatchPreconditioner(const LevelContactNetwork& levelContactNetwork,
+ const LevelContactNetwork& fineContactNetwork,
const Mode mode = LevelPatchPreconditioner::Mode::additive) :
- levelInterfaceNetwork_(levelInterfaceNetwork),
- patchLevelBasis_(patchLevelBasis),
- localAssembler_(localAssembler),
- localInterfaceAssemblers_(localInterfaceAssemblers),
- mode_(mode),
- grid_(levelInterfaceNetwork_.grid()),
- level_(levelInterfaceNetwork_.level()),
- basis_(levelInterfaceNetwork_)
- {
+ mode_(mode),
+ levelContactNetwork_(levelContactNetwork),
+ fineContactNetwork_(fineContactNetwork),
+ level_(levelContactNetwork_.level()),
+ patchFactory_(levelContactNetwork_, fineContactNetwork_) {
+
setPatchDepth();
this->verbosity_ = QUIET;
}
// build support patches
void build() {
+ size_t totalNVertices = 0;
+ for (size_t i=0; i<levelContactNetwork_.nBodies(); i++) {
+ totalNVertices += levelContactNetwork_.body(i)->nVertices();
+ }
-
- localProblems_.resize(vertexInElements.size());
+ patches_.resize(totalNVertices);
// init local fine level corrections
if (this->verbosity_ == FULL) {
@@ -90,19 +88,29 @@ class LevelPatchPreconditioner : public LinearIterationStep<MatrixType, VectorTy
timer.start();
}
- const int patchLevel = patchLevelBasis_.faultNetwork().level();
- for (size_t i=0; i<vertexInElements.size(); i++) {
- if (this->verbosity_ != QUIET) {
- std::cout << i << std::endl;
- std::cout << "---------------" << std::endl;
- }
+ Dune::BitSetVector<1> vertexVisited(totalNVertices);
+ vertexVisited.unsetAll();
- Patch patchDofs;
- SupportPatchFactory<decltype(coarseGridView), decltype(fineGridView)> patchFactory(coarseGridView, fineGridView, vertexInElementsBody);
- patchFactory.build(seedObject, patchDofs, patchDepth_);
+ for (size_t bodyIdx=0; bodyIdx<levelContactNetwork_.nBodies(); bodyIdx++) {
+ const auto& gridView = levelContactNetwork_.body(bodyIdx)->gridView();
- if ((i+1) % 10 == 0 && this->verbosity_ == FULL) {
- std::cout << '\r' << std::floor((i+1.0)/patchLevelBasis_.size()*100) << " % done. Elapsed time: " << timer.elapsed() << " seconds. Predicted total setup time: " << timer.elapsed()/(i+1.0)*patchLevelBasis_.size() << " seconds." << std::flush;
+ for (const auto& e : elements(gridView)) {
+ const auto& refElement = Dune::ReferenceElements<double, dim>::general(e.type());
+
+ for (size_t i=0; i<refElement.size(dim); i++) {
+ auto globalIdx = patchFactory_.coarseIndices().vertexIndex(bodyIdx, e, i);
+
+ if (!vertexVisited[globalIdx][0]) {
+ vertexVisited[globalIdx][0] = true;
+
+ /* if (this->verbosity_ != QUIET) {
+ std::cout << i << std::endl;
+ std::cout << "---------------" << std::endl;
+ }*/
+
+ patchFactory_.build(bodyIdx, e, i, patches_[globalIdx], patchDepth_);
+ }
+ }
}
}
@@ -192,16 +200,8 @@ class LevelPatchPreconditioner : public LinearIterationStep<MatrixType, VectorTy
}
}
- const GridView& gridView() const {
- return basis_.getGridView();
- }
-
- const LevelInterfaceNetwork<GridView>& levelInterfaceNetwork() const {
- return levelInterfaceNetwork_;
- }
-
size_t size() const {
- return localProblems_.size();
+ return patches_.size();
}
};
diff --git a/src/spatial-solving/preconditioners/multilevelpatchpreconditioner.hh b/src/spatial-solving/preconditioners/multilevelpatchpreconditioner.hh
index 584cca16..93529434 100644
--- a/src/spatial-solving/preconditioners/multilevelpatchpreconditioner.hh
+++ b/src/spatial-solving/preconditioners/multilevelpatchpreconditioner.hh
@@ -32,10 +32,9 @@ class MultilevelPatchPreconditioner : public LinearIterationStep<MatrixType, Vec
typedef typename Dune::BitSetVector<1> BitVector;
- const InterfaceNetwork<GridType>& interfaceNetwork_;
+ const ContactNetwork<GridType>& interfaceNetwork_;
const BitVector& activeLevels_;
- const std::vector<std::shared_ptr<LocalOperatorAssembler>>& localAssemblers_;
- const std::vector<std::vector<std::shared_ptr<LocalInterfaceAssembler>>>& localInterfaceAssemblers_;
+
const Mode mode_;
int itCounter_;
@@ -45,16 +44,11 @@ class MultilevelPatchPreconditioner : public LinearIterationStep<MatrixType, Vec
std::vector<VectorType> levelX_;
std::vector<VectorType> levelRhs_;
- std::shared_ptr<LevelInterfaceNetwork<GridView>> allFaultLevelInterfaceNetwork_;
- std::shared_ptr<BasisType> allFaultBasis_;
-
std::vector<std::shared_ptr<DGMGTransfer<BasisType>>> mgTransfer_;
public:
- MultilevelPatchPreconditioner(const InterfaceNetwork<GridType>& interfaceNetwork,
+ MultilevelPatchPreconditioner(const ContactNetwork<GridType>& interfaceNetwork,
const BitVector& activeLevels,
- const std::vector<std::shared_ptr<LocalOperatorAssembler>>& localAssemblers,
- const std::vector<std::vector<std::shared_ptr<LocalInterfaceAssembler>>>& localInterfaceAssemblers,
const Mode mode = MultilevelPatchPreconditioner::Mode::additive) :
interfaceNetwork_(interfaceNetwork),
activeLevels_(activeLevels),
diff --git a/src/spatial-solving/preconditioners/supportpatchfactory.hh b/src/spatial-solving/preconditioners/supportpatchfactory.hh
index 8b191cf0..07091fee 100644
--- a/src/spatial-solving/preconditioners/supportpatchfactory.hh
+++ b/src/spatial-solving/preconditioners/supportpatchfactory.hh
@@ -23,7 +23,7 @@ template <class LevelContactNetwork>
class NetworkIndexMapper {
private:
static const int dim = LevelContactNetwork::dim; //TODO
- using ctype = typename LevelContactNetwork::GridView::ctype;
+ using ctype = typename LevelContactNetwork::ctype;
using FiniteElementCache = Dune::PQkLocalFiniteElementCache<ctype, ctype, dim, 1>; // P1 finite element cache
using Element = typename LevelContactNetwork::GridView::template Codim<0>::Entity;
@@ -270,16 +270,6 @@ 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();
-
- BodyElement seedElement;
- seedElement.set(bodyID, coarseElement);
-
auto isPatchIntersection = [this](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);
@@ -302,24 +292,26 @@ class SupportPatchFactory
return false;
};
- // construct coarse patch
- std::vector<BodyElement> patchElements;
- patchElements.emplace_back(seedElement);
+ patchDofs.resize(nVertices_);
+ patchDofs.setAll();
+ std::vector<BodyElement> patchElements;
std::set<size_t> visitedElements;
- //visitedElements.insert(coarseIndices_.elementIndex(seedElement.bodyID, seedElement.element));
std::set<size_t> patchVertices;
patchVertices.insert(coarseIndices_.vertexIndex(bodyID, coarseElement, localVertex));
+ std::set<size_t> visitedBodies;
+
+ BodyElement seedElement;
+ seedElement.set(bodyID, coarseElement);
+
std::queue<BodyElement> patchSeeds;
patchSeeds.push(seedElement);
- std::set<size_t> visitedBodies;
-
for (size_t depth=0; depth<=patchDepth; depth++) {
- std::vector<BodyElement> nextSeeds;
std::set<size_t> newPatchVertices;
+ std::vector<BodyElement> nextSeeds;
while (!patchSeeds.empty()) {
const auto& patchSeed = patchSeeds.front();
@@ -332,8 +324,8 @@ class SupportPatchFactory
}
visitedBodies.insert(patchSeed.bodyID);
-
visitedElements.insert(elemIdx);
+ patchElements.emplace_back(patchSeed);
size_t bodyIdx = patchSeed.bodyID;
const auto& elem = patchSeed.element;
@@ -347,12 +339,9 @@ class SupportPatchFactory
size_t neighborIdx = coarseIndices_.elementIndex(neighbor.bodyID, neighbor.element);
- std::cout << "elementID: " << neighborIdx << std::endl;
-
if (visitedElements.count(neighborIdx))
continue;
- patchElements.emplace_back(neighbor);
patchSeeds.push(neighbor);
} else {
size_t faceIdx = coarseIndices_.faceIndex(bodyIdx, it);
@@ -366,24 +355,22 @@ class SupportPatchFactory
BodyElement neighbor;
neighbor.set(coarseContactNetwork_, intersection);
- size_t neighborIdx = coarseIndices_.elementIndex(neighbor.bodyID, neighbor.element);
+ size_t neighborBody = neighbor.bodyID;
+ size_t neighborIdx = coarseIndices_.elementIndex(neighborBody, neighbor.element);
- if (visitedElements.count(neighborIdx) || visitedBodies.count(neighbor.bodyID))
+ if (visitedElements.count(neighborIdx) || visitedBodies.count(neighborBody))
continue;
std::vector<std::set<size_t>> dofs;
- remoteIntersectionDofs(coarseIndices_, rIs, bodyIdx, neighbor.bodyID, dofs, intersection.flip);
+ remoteIntersectionDofs(coarseIndices_, rIs, bodyIdx, neighborBody, dofs, intersection.flip);
if (isRPatchIntersection(dofs, patchVertices)) {
patchVertices.insert(dofs[1].begin(), dofs[1].end());
- patchElements.emplace_back(neighbor);
patchSeeds.push(neighbor);
} else {
newPatchVertices.insert(dofs[1].begin(), dofs[1].end());
nextSeeds.emplace_back(neighbor);
}
-
- print(patchVertices, "patchVertices: ");
}
}
}
@@ -392,7 +379,7 @@ class SupportPatchFactory
BodyElement neighbor;
neighbor.set(bodyIdx, it.outside());
- size_t neighborIdx = coarseIndices_.elementIndex(neighbor.bodyID, neighbor.element);
+ size_t neighborIdx = coarseIndices_.elementIndex(bodyIdx, neighbor.element);
if (visitedElements.count(neighborIdx))
continue;
@@ -407,9 +394,10 @@ class SupportPatchFactory
BodyElement neighbor;
neighbor.set(coarseContactNetwork_, intersections[i]);
+ size_t neighborBody = neighbor.bodyID;
size_t neighborIdx = coarseIndices_.elementIndex(neighbor.bodyID, neighbor.element);
- if (visitedElements.count(neighborIdx) || visitedBodies.count(neighbor.bodyID))
+ if (visitedElements.count(neighborIdx) || visitedBodies.count(neighborBody))
continue;
newPatchVertices.insert(neighborFaceDofs_[faceIdx].begin(), neighborFaceDofs_[faceIdx].end());
@@ -419,33 +407,22 @@ class SupportPatchFactory
}
}
}
-
patchSeeds.pop();
}
- for (size_t i=0; i<nextSeeds.size(); i++) {
- patchSeeds.push(nextSeeds[i]);
+ for (size_t j=0; j<nextSeeds.size(); j++) {
+ patchSeeds.push(nextSeeds[j]);
}
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;
-
- std::cout << "coarse element: " << elemIdx << std::endl;
-
- visited.insert(elemIdx);
const auto& grid = coarseContactNetwork_.body(coarseElem.bodyID)->gridView().grid();
const auto fineLevel = fineContactNetwork_.body(coarseElem.bodyID)->level();
@@ -458,14 +435,16 @@ class SupportPatchFactory
addLocalDofs(coarseElem, fineElem, visitedElements, patchBoundary, patchDofs);
}
}
-
- std::cout << "Success!" << std::endl;
}
const auto& coarseContactNetwork() {
return coarseContactNetwork_;
}
+ const auto& coarseIndices() const {
+ return coarseIndices_;
+ }
+
private:
void setup(const LevelContactNetwork& levelContactNetwork, const NetworkIndexMapper& indices, std::map<size_t, std::vector<RemoteIntersection>>& faceToRIntersections) {
for (size_t i=0; i<levelContactNetwork.nCouplings(); i++) {
--
GitLab