diff --git a/dune/faultnetworks/patchfactories/cellpatchfactory.hh b/dune/faultnetworks/patchfactories/cellpatchfactory.hh
index 7b2379e46ce7eb606f25a34c582dd54f9114ad60..c2572ac0bb05cdc7918f701ff285cbbd8d9ed1fa 100644
--- a/dune/faultnetworks/patchfactories/cellpatchfactory.hh
+++ b/dune/faultnetworks/patchfactories/cellpatchfactory.hh
@@ -15,7 +15,7 @@
template <class BasisType>
class CellPatchFactory
{
- protected:
+protected:
typedef typename BasisType::GridView GV;
typedef typename BasisType::LocalFiniteElement LFE;
@@ -26,42 +26,35 @@ class CellPatchFactory
typedef typename Dune::LevelMultipleCodimMultipleGeomTypeMapper<GridType, Dune::MCMGElementLayout > ElementMapper;
- private:
- const BasisType& basis_;
- const int level_;
+ const BasisType& coarseBasis_;
+ const BasisType& fineBasis_;
+
+ const int coarseLevel_;
+ const int fineLevel_;
+
+ const GridType& grid;
+ const GV& coarseGridView;
+ const GV& fineGridView;
- const GridType& grid;
- const GV& gridView;
-
std::vector<std::vector<int>> cellLocalToGlobal;
std::vector<Dune::BitSetVector<1>> cellBoundaryDofs;
std::vector<std::vector<Element>> cellElements;
ElementMapper mapper;
-private:
- /*void print(const Dune::BitSetVector<1>& x, std::string message){
- std::cout << message << std::endl;
- for (size_t i=0; i<x.size(); i++) {
- std::cout << x[i][0] << " ";
- }
- std::cout << std::endl << std::endl;
- }*/
-
- bool containsInsideSubentity(const Element& elem, const typename GridType::LevelIntersection& intersection, int subEntity, int codim) {
- return ReferenceElementHelper<double, dim>::subEntityContainsSubEntity(elem.type(), intersection.indexInInside(), 1, subEntity, codim);
- }
-
public:
//setup
- CellPatchFactory(const BasisType& basis, const int level) :
- basis_(basis),
- level_(level),
- grid(basis_.getGridView().grid()),
- gridView(basis_.getGridView()),
- mapper(grid, level_) {
-
- const auto& interfaces = basis_.faultNetwork();
+ CellPatchFactory(const BasisType& coarseBasis, const BasisType& fineBasis, const int coarseLevel, const int fineLevel) :
+ coarseBasis_(coarseBasis),
+ fineBasis_(fineBasis),
+ coarseLevel_(coarseLevel),
+ fineLevel_(fineLevel),
+ grid(coarseBasis_.getGridView().grid()),
+ coarseGridView(coarseBasis_.getGridView()),
+ fineGridView(fineBasis_.getGridView()),
+ mapper(grid, coarseLevel_) {
+
+ const auto& interfaces = coarseBasis_.faultNetwork();
cellElements.resize(0);
@@ -69,7 +62,7 @@ public:
visited.unsetAll();
std::queue<Element> cellSeeds;
- cellSeeds.push(*gridView. template begin <0>());
+ cellSeeds.push(*coarseGridView. template begin <0>());
size_t cellCount = 0;
while (!cellSeeds.empty()) {
@@ -81,13 +74,11 @@ public:
cellCount++;
- cellElements.resize(cellCount);
- auto& elements = cellElements.back();
+ std::vector<Element> coarseElements;
std::queue<Element> elementQueue;
elementQueue.push(cellSeed);
-
cellLocalToGlobal.resize(cellCount);
cellBoundaryDofs.resize(cellCount);
std::set<int> dofs;
@@ -101,17 +92,10 @@ public:
if (visited[elemID][0] == true)
continue;
- elements.push_back(elem);
+ coarseElements.push_back(elem);
visited[elemID][0] = true;
- const LFE& lfe = basis_.getLocalFiniteElement(elem);
- // insert dofs
- for (size_t i=0; i<lfe.localBasis().size(); ++i) {
- int dofIndex = basis_.index(elem, i);
- dofs.insert(dofIndex);
- }
-
- for (const auto& is : intersections(gridView, elem)) {
+ for (const auto& is : intersections(coarseGridView, elem)) {
if (is.neighbor()) {
const auto neighbor = is.outside();
const auto neighborId = mapper.index(neighbor);
@@ -124,20 +108,16 @@ public:
}
}
}
+ }
+ }
- // set boundaryDofs
- if (is.boundary()) {
- const auto& localCoefficients = lfe.localCoefficients();
-
- for (size_t i=0; i<localCoefficients.size(); i++) {
- auto entity = localCoefficients.localKey(i).subEntity();
- auto codim = localCoefficients.localKey(i).codim();
+ cellElements.resize(cellCount);
+ auto& elements = cellElements.back();
- if (containsInsideSubentity(elem, is, entity, codim))
- boundaryDofs.insert(basis_.index(elem, i));
- }
- }
- }
+ // construct fine patch
+ for (size_t i=0; i<coarseElements.size(); i++) {
+ const Element& elem = coarseElements[i];
+ addFinePatchElements(elem, dofs, boundaryDofs, elements);
}
auto& localToGlobal = cellLocalToGlobal.back();
@@ -170,5 +150,54 @@ public:
auto& getBoundaryDofs() {
return cellBoundaryDofs;
}
+
+private:
+ bool containsInsideSubentity(const Element& elem, const typename GridType::LevelIntersection& intersection, int subEntity, int codim) {
+ return ReferenceElementHelper<double, dim>::subEntityContainsSubEntity(elem.type(), intersection.indexInInside(), 1, subEntity, codim);
+ }
+
+ void addFinePatchElements(const Element& coarseElem, std::set<int>& localDofs, std::set<int>& localBoundaryDofs, std::vector<Element>& fineElements) {
+ using HierarchicLevelIteratorType = HierarchicLevelIterator<GridType>;
+ HierarchicLevelIteratorType endHierIt(grid, coarseElem, HierarchicLevelIteratorType::PositionFlag::end, fineLevel_);
+ for (HierarchicLevelIteratorType hierIt(grid, coarseElem, HierarchicLevelIteratorType::PositionFlag::begin, fineLevel_); hierIt!=endHierIt; ++hierIt) {
+
+ const Element& fineElem = *hierIt;
+ fineElements.push_back(fineElem);
+
+ addLocalDofs(coarseElem, fineElem, localDofs, localBoundaryDofs);
+ }
+ }
+
+ void addLocalDofs(const Element& coarseElem, const Element& fineElem, std::set<int>& localDofs, std::set<int>& localBoundaryDofs) {
+ const LFE& fineLFE = fineBasis_.getLocalFiniteElement(fineElem);
+
+ /*
+ const auto& fineGeometry = fineElem.geometry();
+ const auto& coarseGeometry = coarseElem.geometry();
+
+ const auto& ref = Dune::ReferenceElements<ctype, dim>::general(fineElem.type()); */
+
+ // insert dofs
+ for (size_t i=0; i<fineLFE.localBasis().size(); ++i) {
+ int dofIndex = fineBasis_.index(fineElem, i);
+ localDofs.insert(dofIndex);
+ }
+
+ for (const auto& is : intersections(fineGridView, fineElem)) {
+
+ // set boundaryDofs
+ if (is.boundary()) {
+ const auto& localCoefficients = fineLFE.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))
+ localBoundaryDofs.insert(fineBasis_.index(fineElem, i));
+ }
+ }
+ }
+ }
};
#endif
diff --git a/dune/faultnetworks/patchfactories/supportpatchfactory.hh b/dune/faultnetworks/patchfactories/supportpatchfactory.hh
index 07a297d13cd0151395e9a15ba0e726aeeafc47fd..6630e305eff35ae51fb7f2b4f001c51c6749e0e8 100644
--- a/dune/faultnetworks/patchfactories/supportpatchfactory.hh
+++ b/dune/faultnetworks/patchfactories/supportpatchfactory.hh
@@ -203,7 +203,7 @@ private:
const typename GridType::LevelIntersection& intersection = *neighborIt;
bool isInnerBoundary = false;
- if (intersection.neighbor()) {
+ if (intersection.neighbor()) {
Element outsideFather = intersection.outside();
while (outsideFather.level()>coarseLevel_) {
@@ -213,7 +213,7 @@ private:
if (!inCoarsePatch[coarseMapper.index(outsideFather)][0]) {
isInnerBoundary = true;
}
- }
+ }
if (intersection.boundary() or isInnerBoundary) {
typedef typename LFE::Traits::LocalCoefficientsType LocalCoefficients;
diff --git a/dune/faultnetworks/preconditioners/cellpatchpreconditioner.hh b/dune/faultnetworks/preconditioners/cellpatchpreconditioner.hh
index 8fc86fc57494db238a27461a1e3ca282ecc07833..c4c2b04fa72c2ee04a0c0e48140516f913575d7c 100644
--- a/dune/faultnetworks/preconditioners/cellpatchpreconditioner.hh
+++ b/dune/faultnetworks/preconditioners/cellpatchpreconditioner.hh
@@ -15,13 +15,15 @@ private:
public:
CellPatchPreconditioner(const LevelInterfaceNetwork<GridView>& levelInterfaceNetwork,
+ const BasisType& patchLevelBasis,
const LocalAssembler& localAssembler,
const std::vector<std::shared_ptr<LocalInterfaceAssembler>>& localInterfaceAssemblers,
const Mode mode = Mode::additive) :
- Base(levelInterfaceNetwork, localAssembler, localInterfaceAssemblers, mode) {}
+ Base(levelInterfaceNetwork, patchLevelBasis, localAssembler, localInterfaceAssemblers, mode) {}
void build() {
- CellPatchFactory<BasisType> patchFactory(this->basis_, this->level_);
+ const int patchLevel = this->patchLevelBasis_.faultNetwork().level();
+ CellPatchFactory<BasisType> patchFactory(this->patchLevelBasis_, this->basis_, patchLevel, this->level_);
const auto& localToGlobal = patchFactory.getLocalToGlobal();
const auto& boundaryDofs = patchFactory.getBoundaryDofs();
diff --git a/dune/faultnetworks/preconditioners/levelpatchpreconditioner.hh b/dune/faultnetworks/preconditioners/levelpatchpreconditioner.hh
index d8c3090d5c77782b7ef9c4838b43b1222bc35f96..0cfd9163d4c5c3c0756c479caa2d0fdfcca73091 100644
--- a/dune/faultnetworks/preconditioners/levelpatchpreconditioner.hh
+++ b/dune/faultnetworks/preconditioners/levelpatchpreconditioner.hh
@@ -30,6 +30,7 @@ protected:
typedef typename GridView::Grid GridType;
const LevelInterfaceNetwork<GridView>& levelInterfaceNetwork_;
+ const BasisType& patchLevelBasis_;
const LocalAssembler& localAssembler_;
const std::vector<std::shared_ptr<LocalInterfaceAssembler>>& localInterfaceAssemblers_;
@@ -48,10 +49,12 @@ protected:
public:
LevelPatchPreconditioner(const LevelInterfaceNetwork<GridView>& levelInterfaceNetwork,
+ const BasisType& patchLevelBasis,
const LocalAssembler& localAssembler,
const std::vector<std::shared_ptr<LocalInterfaceAssembler>>& localInterfaceAssemblers,
const Mode mode = LevelPatchPreconditioner::Mode::ADDITIVE) :
levelInterfaceNetwork_(levelInterfaceNetwork),
+ patchLevelBasis_(patchLevelBasis),
localAssembler_(localAssembler),
localInterfaceAssemblers_(localInterfaceAssemblers),
mode_(mode),
diff --git a/dune/faultnetworks/preconditioners/multilevelpatchpreconditioner.hh b/dune/faultnetworks/preconditioners/multilevelpatchpreconditioner.hh
index d238009eed605b788e6fc358365dad5dbef727e7..e4bf038e631ba11dc8f103e985b782ead0f786b6 100644
--- a/dune/faultnetworks/preconditioners/multilevelpatchpreconditioner.hh
+++ b/dune/faultnetworks/preconditioners/multilevelpatchpreconditioner.hh
@@ -283,8 +283,14 @@ private:
}
if (activeLevels_[i][0]) {
- levelPatchPreconditioners_.push_back(std::make_shared<CellPreconditioner>(interfaceNetwork_.levelInterfaceNetwork(i), *localAssemblers_[i], localInterfaceAssemblers_[i], mode_));
- levelPatchPreconditioners_.back()->setBoundaryMode(LevelPatchPreconditionerType::BoundaryMode::homogeneous);
+ const LevelInterfaceNetwork<GridView>& levelNetwork = interfaceNetwork_.levelInterfaceNetwork(i);
+
+ if (levelPatchPreconditioners_.size() == 0) {
+ levelPatchPreconditioners_.push_back(std::make_shared<CellPreconditioner>(levelNetwork, coarseGlobalPreconditioner_->basis(), *localAssemblers_[i], localInterfaceAssemblers_[i], mode_));
+ } else {
+ levelPatchPreconditioners_.push_back(std::make_shared<CellPreconditioner>(levelNetwork, levelPatchPreconditioners_.back()->basis(), *localAssemblers_[i], localInterfaceAssemblers_[i], mode_));
+ levelPatchPreconditioners_.back()->setBoundaryMode(LevelPatchPreconditionerType::BoundaryMode::homogeneous);
+ }
maxLevel = i;
}
diff --git a/dune/faultnetworks/preconditioners/supportpatchpreconditioner.hh b/dune/faultnetworks/preconditioners/supportpatchpreconditioner.hh
index 7acb2c1323aa3f370667ec046ef1fa336947f265..0711cd198a9a7c6a805f465546b41cedd4b72462 100644
--- a/dune/faultnetworks/preconditioners/supportpatchpreconditioner.hh
+++ b/dune/faultnetworks/preconditioners/supportpatchpreconditioner.hh
@@ -14,8 +14,6 @@ private:
using GridType = typename Base::GridType;
using Mode = typename Base::Mode;
- const BasisType& patchLevelBasis_;
-
public:
// for each coarse patch given by patchLevelBasis: set up local problem, compute local correction
SupportPatchPreconditioner(const LevelInterfaceNetwork<GridView>& levelInterfaceNetwork,
@@ -23,8 +21,8 @@ public:
const LocalAssembler& localAssembler,
const std::vector<std::shared_ptr<LocalInterfaceAssembler>>& localInterfaceAssemblers,
const Mode mode = Mode::ADDITIVE) :
- Base(levelInterfaceNetwork, localAssembler, localInterfaceAssemblers, mode),
- patchLevelBasis_(patchLevelBasis) {}
+ Base(levelInterfaceNetwork, patchLevelBasis, localAssembler, localInterfaceAssemblers, mode)
+ {}
void build() {
// init vertexInElements
@@ -32,7 +30,7 @@ public:
using Element = typename GridType::template Codim<0>::Entity;
std::vector<std::vector<Element>> vertexInElements;
- const GridView& patchLevelGridView = patchLevelBasis_.getGridView();
+ const GridView& patchLevelGridView = this->patchLevelBasis_.getGridView();
vertexInElements.resize(patchLevelGridView.size(dim));
for (size_t i=0; i<vertexInElements.size(); i++) {
vertexInElements[i].resize(0);
@@ -40,10 +38,10 @@ public:
for (const auto& elem : elements(patchLevelGridView)) {
// compute coarseGrid vertexInElements
- const auto& coarseFE = patchLevelBasis_.getLocalFiniteElement(elem);
+ const auto& coarseFE = this->patchLevelBasis_.getLocalFiniteElement(elem);
for (size_t i=0; i<coarseFE.localBasis().size(); ++i) {
- int dofIndex = patchLevelBasis_.indexInGridView(elem, i);
+ int dofIndex = this->patchLevelBasis_.indexInGridView(elem, i);
vertexInElements[dofIndex].push_back(elem);
}
}
@@ -51,9 +49,9 @@ public:
std::cout << "SupportPatchPreconditioner::build() level: " << this->level_ << std::endl;
this->localProblems_.resize(vertexInElements.size());
- const int patchLevel = patchLevelBasis_.faultNetwork().level();
+ const int patchLevel = this->patchLevelBasis_.faultNetwork().level();
for (size_t i=0; i<vertexInElements.size(); i++) {
- SupportPatchFactory<BasisType> patchFactory(patchLevelBasis_, this->basis_, patchLevel, this->level_, vertexInElements, i, this->patchDepth_);
+ SupportPatchFactory<BasisType> patchFactory(this->patchLevelBasis_, this->basis_, patchLevel, this->level_, vertexInElements, i, this->patchDepth_);
const auto& localToGlobal = patchFactory.getLocalToGlobal();
const auto& boundaryDofs = patchFactory.getBoundaryDofs();
diff --git a/src/cantorfaultnetworks/sparsecantorfaultnetwork.parset b/src/cantorfaultnetworks/sparsecantorfaultnetwork.parset
index f5370c0508d456d281c9930a30b35e02bdc47ce1..aad8acdf9af41f44e0ce44d1a5f012114f7bde43 100644
--- a/src/cantorfaultnetworks/sparsecantorfaultnetwork.parset
+++ b/src/cantorfaultnetworks/sparsecantorfaultnetwork.parset
@@ -2,8 +2,8 @@ path = ../data/
resultPath = ../cantorfaultnetworks/results/sparse/
[preconditioner]
-patch = CELL # CELL , SUPPORT
-mode = ADDITIVE # ADDITIVE , MULTIPLICATIVE
+patch = SUPPORT # CELL , SUPPORT
+mode = MULTIPLICATIVE # ADDITIVE , MULTIPLICATIVE
multDirection = SYMMETRIC # SYMMETRIC , FORWARD , BACKWARD
patchDepth = 1