From acfdc49ca65d2d192b7a35f39e9d9f1d49ca0b3f Mon Sep 17 00:00:00 2001
From: podlesny <podlesny@zedat.fu-berlin.de>
Date: Thu, 4 Jun 2020 14:48:11 +0200
Subject: [PATCH] .
---
dune/faultnetworks/localproblem.hh | 4 +
.../patchfactories/supportpatchfactory.hh | 87 ++++++++++++++-----
.../levelpatchpreconditioner.hh | 4 +-
.../supportpatchpreconditioner.hh | 8 +-
dune/faultnetworks/solvers/osccgsolver.cc | 6 +-
.../sparsecantorfaultnetwork.cc | 2 +-
.../sparsecantorfaultnetwork.parset | 4 +-
7 files changed, 80 insertions(+), 35 deletions(-)
diff --git a/dune/faultnetworks/localproblem.hh b/dune/faultnetworks/localproblem.hh
index 7ae7a56..5f115f8 100644
--- a/dune/faultnetworks/localproblem.hh
+++ b/dune/faultnetworks/localproblem.hh
@@ -145,6 +145,10 @@ public:
}
}
+ void updateLocalRhs(const RangeType& rhs){
+ localRhs = rhs;
+ }
+
void updateRhsAndBoundary(const RangeType& rhs, const RangeType& boundary) {
for (size_t i=0; i<localRhs.size(); i++) {
if (boundaryDofs_[i][0]) {
diff --git a/dune/faultnetworks/patchfactories/supportpatchfactory.hh b/dune/faultnetworks/patchfactories/supportpatchfactory.hh
index 6630e30..97618bb 100644
--- a/dune/faultnetworks/patchfactories/supportpatchfactory.hh
+++ b/dune/faultnetworks/patchfactories/supportpatchfactory.hh
@@ -12,6 +12,8 @@
#include <dune/faultnetworks/hierarchicleveliterator.hh>
+#include <dune/faultnetworks/utils/debugutils.hh>
+
template <class BasisType>
class SupportPatchFactory
{
@@ -76,33 +78,68 @@ class SupportPatchFactory
std::set<int> localDofs;
std::set<int> localBoundaryDofs;
- addVertexSupport(centerVertexIdx_, coarseElements, visited, vertexVisited);
+ //std::set<int> centerVertices;
+ //centerVertices.insert(centerVertexIdx_);
+
+ std::set<int> coarseBoundaryElements;
+
+ addVertexSupport(centerVertexIdx_, coarseElements, visited, vertexVisited);
+ /*for (size_t i=0; i<coarseElements.size(); ++i) {
+ const auto elem = coarseElements[i];
+
+ for (const auto& is : intersections(coarseGridView, elem)) {
+ if (!is.neighbor() or is.boundary())
+ continue;
+
+ if (visited[coarseMapper.index(is.outside())][0])
+ continue;
+
+ const auto& localCoefficients = coarseBasis_.getLocalFiniteElement(is.inside()).localCoefficients();
+
+ for (size_t i=0; i<localCoefficients.size(); i++) {
+ unsigned int entity = localCoefficients.localKey(i).subEntity();
+ unsigned int codim = localCoefficients.localKey(i).codim();
+
+ auto idx = coarseBasis_.index(is.inside(), i);
+ auto as = containsInsideSubentity(elem, is, entity, codim);
+ auto asd = centerVertices.count(idx);
+
+ if (containsInsideSubentity(elem, is, entity, codim) and centerVertices.count(coarseBasis_.index(is.inside(), i))) {
+ coarseBoundaryElements.insert(coarseMapper.index(is.outside()));
+ break;
+ }
+ }
+ }
+ }*/
+
// construct coarse patch
for (int depth=1; depth<patchDepth_; depth++) {
- std::vector<Element> coarseElementNeighbors;
+ std::vector<Element> coarseElementNeighbors;
- for (size_t i=0; i<coarseElements.size(); ++i) {
- const Element& elem = coarseElements[i];
- const LFE& coarseLFE = coarseBasis_.getLocalFiniteElement(elem);
+ for (size_t i=0; i<coarseElements.size(); ++i) {
+ const Element& elem = coarseElements[i];
+ const LFE& coarseLFE = coarseBasis_.getLocalFiniteElement(elem);
- for (size_t j=0; j<coarseLFE.localBasis().size(); ++j) {
- int dofIndex = coarseBasis_.indexInGridView(elem, j);
+ for (size_t j=0; j<coarseLFE.localBasis().size(); ++j) {
+ int dofIndex = coarseBasis_.indexInGridView(elem, j);
- if (!vertexVisited[dofIndex][0]) {
- addVertexSupport(dofIndex, coarseElementNeighbors, visited, vertexVisited);
- }
- }
- }
+ if (!vertexVisited[dofIndex][0]) {
+ addVertexSupport(dofIndex, coarseElementNeighbors, visited, vertexVisited);
+ }
+ }
+ }
- coarseElements.insert(coarseElements.begin(), coarseElementNeighbors.begin(), coarseElementNeighbors.end());
+ coarseElements.insert(coarseElements.begin(), coarseElementNeighbors.begin(), coarseElementNeighbors.end());
}
+ //print(coarseBoundaryElements, "coarseBoundaryElements");
+
// construct fine patch
for (size_t i=0; i<coarseElements.size(); ++i) {
const Element& elem = coarseElements[i];
- addFinePatchElements(elem, visited, localDofs, localBoundaryDofs);
+ addFinePatchElements(elem, visited, localDofs, localBoundaryDofs, coarseBoundaryElements);
}
@@ -141,13 +178,13 @@ class SupportPatchFactory
private:
- void print(const Dune::BitSetVector<1>& x, std::string message){
+ /* 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;
- }
+ }*/
void addVertexSupport(const int vertexIdx, std::vector<Element>& coarseElements, Dune::BitSetVector<1>& visited, Dune::BitSetVector<1>& vertexVisited) {
const std::vector<Element>& vertexInElement = vertexInElements_[vertexIdx];
@@ -169,18 +206,19 @@ private:
}
- void addFinePatchElements(const Element& coarseElem, const Dune::BitSetVector<1>& inCoarsePatch, std::set<int>& localDofs, std::set<int>& localBoundaryDofs) {
+ void addFinePatchElements(const Element& coarseElem, const Dune::BitSetVector<1>& inCoarsePatch, std::set<int>& localDofs, std::set<int>& localBoundaryDofs, const std::set<int>& coarseBoundaryElements) {
HierarchicLevelIteratorType endHierIt(grid, coarseElem, HierarchicLevelIteratorType::PositionFlag::end, fineLevel_);
for (HierarchicLevelIteratorType hierIt(grid, coarseElem, HierarchicLevelIteratorType::PositionFlag::begin, fineLevel_); hierIt!=endHierIt; ++hierIt) {
const Element& fineElem = *hierIt;
fineRegionElements.push_back(fineElem);
- addLocalDofs(coarseElem, fineElem, inCoarsePatch, localDofs, localBoundaryDofs);
+ addLocalDofs(coarseElem, fineElem, inCoarsePatch, localDofs, localBoundaryDofs, coarseBoundaryElements);
}
}
- void addLocalDofs(const Element& coarseElem, const Element& fineElem, const Dune::BitSetVector<1>& inCoarsePatch, std::set<int>& localDofs, std::set<int>& localBoundaryDofs) {
+ void addLocalDofs(const Element& coarseElem, const Element& fineElem, const Dune::BitSetVector<1>& inCoarsePatch,
+ std::set<int>& localDofs, std::set<int>& localBoundaryDofs, const std::set<int>& coarseBoundaryElements) {
const LFE& fineLFE = fineBasis_.getLocalFiniteElement(fineElem);
/*
@@ -210,8 +248,9 @@ private:
outsideFather = outsideFather.father();
}
- if (!inCoarsePatch[coarseMapper.index(outsideFather)][0]) {
- isInnerBoundary = true;
+ auto outsideFatherIdx = coarseMapper.index(outsideFather);
+ if (!inCoarsePatch[outsideFatherIdx][0] and !coarseBoundaryElements.count(outsideFatherIdx)) {
+ isInnerBoundary = true;
}
}
@@ -219,9 +258,9 @@ private:
typedef typename LFE::Traits::LocalCoefficientsType LocalCoefficients;
const LocalCoefficients* localCoefficients = &fineBasis_.getLocalFiniteElement(intersection.inside()).localCoefficients();
- for (size_t i=0; i<localCoefficients->size(); i++) {
- unsigned int entity = localCoefficients->localKey(i).subEntity();
- unsigned int codim = localCoefficients->localKey(i).codim();
+ for (size_t i=0; i<localCoefficients->size(); i++) {
+ unsigned int entity = localCoefficients->localKey(i).subEntity();
+ unsigned int codim = localCoefficients->localKey(i).codim();
if (containsInsideSubentity(fineElem, intersection, entity, codim))
localBoundaryDofs.insert(fineBasis_.index(intersection.inside(), i));
diff --git a/dune/faultnetworks/preconditioners/levelpatchpreconditioner.hh b/dune/faultnetworks/preconditioners/levelpatchpreconditioner.hh
index 0cfd916..c43f777 100644
--- a/dune/faultnetworks/preconditioners/levelpatchpreconditioner.hh
+++ b/dune/faultnetworks/preconditioners/levelpatchpreconditioner.hh
@@ -186,7 +186,9 @@ private:
localProblem.restrict(*(this->x_), v);
VectorType r;
- localProblem.localResidual(v, r);
+ localProblem.restrict(*(this->rhs_), r);
+ Dune::MatrixVector::subtractProduct(r, localProblem.getLocalMat(), v);
+ localProblem.updateLocalRhs(r);
VectorType x;
localProblem.solve(x);
diff --git a/dune/faultnetworks/preconditioners/supportpatchpreconditioner.hh b/dune/faultnetworks/preconditioners/supportpatchpreconditioner.hh
index 0711cd1..68f544c 100644
--- a/dune/faultnetworks/preconditioners/supportpatchpreconditioner.hh
+++ b/dune/faultnetworks/preconditioners/supportpatchpreconditioner.hh
@@ -6,6 +6,8 @@
#include <dune/faultnetworks/preconditioners/levelpatchpreconditioner.hh>
#include <dune/faultnetworks/patchfactories/supportpatchfactory.hh>
+#include <dune/faultnetworks/utils/debugutils.hh>
+
template <class BasisType, class LocalAssembler, class LocalInterfaceAssembler, class MatrixType, class VectorType>
class SupportPatchPreconditioner : public LevelPatchPreconditioner<BasisType, LocalAssembler, LocalInterfaceAssembler, MatrixType, VectorType> {
private:
@@ -31,7 +33,7 @@ public:
std::vector<std::vector<Element>> vertexInElements;
const GridView& patchLevelGridView = this->patchLevelBasis_.getGridView();
- vertexInElements.resize(patchLevelGridView.size(dim));
+ vertexInElements.resize(patchLevelGridView.size(dim)); //(this->patchLevelBasis_.size());
for (size_t i=0; i<vertexInElements.size(); i++) {
vertexInElements[i].resize(0);
}
@@ -41,13 +43,15 @@ public:
const auto& coarseFE = this->patchLevelBasis_.getLocalFiniteElement(elem);
for (size_t i=0; i<coarseFE.localBasis().size(); ++i) {
- int dofIndex = this->patchLevelBasis_.indexInGridView(elem, i);
+ int dofIndex = this->patchLevelBasis_.indexInGridView(elem, i); //index
vertexInElements[dofIndex].push_back(elem);
}
}
std::cout << "SupportPatchPreconditioner::build() level: " << this->level_ << std::endl;
+ //printBasisDofLocation(this->patchLevelBasis_);
+
this->localProblems_.resize(vertexInElements.size());
const int patchLevel = this->patchLevelBasis_.faultNetwork().level();
for (size_t i=0; i<vertexInElements.size(); i++) {
diff --git a/dune/faultnetworks/solvers/osccgsolver.cc b/dune/faultnetworks/solvers/osccgsolver.cc
index 4a63d47..772f36e 100644
--- a/dune/faultnetworks/solvers/osccgsolver.cc
+++ b/dune/faultnetworks/solvers/osccgsolver.cc
@@ -232,10 +232,6 @@ void OscCGSolver<MatrixType, VectorType, DGBasisType>::solve()
std::cout << "--------------------\n";
-
- if (this->useRelativeError_)
- std::cout << "Norm of final discretization error: " << (discretizationError*normOfInitialDiscretizationError) << std::endl;
- else
- std::cout << "Norm of final discretization error: " << discretizationError << std::endl;
+ std::cout << "Norm of final discretization error: " << discretizationError << std::endl;
}
}
diff --git a/src/cantorfaultnetworks/sparsecantorfaultnetwork.cc b/src/cantorfaultnetworks/sparsecantorfaultnetwork.cc
index b359981..c0918ab 100644
--- a/src/cantorfaultnetworks/sparsecantorfaultnetwork.cc
+++ b/src/cantorfaultnetworks/sparsecantorfaultnetwork.cc
@@ -333,7 +333,7 @@ int main(int argc, char** argv) { try
Solver::FULL, &exactSol, &fineSol, discMGTransfer, 1.0, true); //((oscGridN+0.0)/fineGridN)
- solver.historyBuffer_ = resultPath + "solutions/level_" + std::to_string(fineCantorLevel) + ".vec";
+ //solver.historyBuffer_ = resultPath + "solutions/level_" + std::to_string(fineCantorLevel) + ".vec";
solver.check();
solver.preprocess();
diff --git a/src/cantorfaultnetworks/sparsecantorfaultnetwork.parset b/src/cantorfaultnetworks/sparsecantorfaultnetwork.parset
index aad8acd..3b3ea31 100644
--- a/src/cantorfaultnetworks/sparsecantorfaultnetwork.parset
+++ b/src/cantorfaultnetworks/sparsecantorfaultnetwork.parset
@@ -2,7 +2,7 @@ path = ../data/
resultPath = ../cantorfaultnetworks/results/sparse/
[preconditioner]
-patch = SUPPORT # CELL , SUPPORT
+patch = CELL # CELL , SUPPORT
mode = MULTIPLICATIVE # ADDITIVE , MULTIPLICATIVE
multDirection = SYMMETRIC # SYMMETRIC , FORWARD , BACKWARD
patchDepth = 1
@@ -13,7 +13,7 @@ patchDepth = 1
oscDataFile = oscDataLaplace4.mat
# level resolution in 2^(-...)
-coarseCantorLevel = 1
+coarseCantorLevel = 2
fineCantorLevel = 3
maxCantorLevel = 4
--
GitLab