diff --git a/dune/faultnetworks/preconditioners/cellpatchpreconditioner.hh b/dune/faultnetworks/preconditioners/cellpatchpreconditioner.hh
index a14f005352084a4c1dd86aa6a47c6fa60e08b2c0..dae195f1fe60ef0d4ee0b1791edfcce8fed47ad4 100644
--- a/dune/faultnetworks/preconditioners/cellpatchpreconditioner.hh
+++ b/dune/faultnetworks/preconditioners/cellpatchpreconditioner.hh
@@ -28,7 +28,7 @@ public:
const auto& boundaryDofs = patchFactory.getBoundaryDofs();
VectorType rhs;
- rhs.resize(this->mat_->N());
+ rhs.resize(this->matrix_.N());
rhs = 0;
std::cout << "CellPatchPreconditioner::build() level: " << this->level_ << std::endl;
@@ -38,7 +38,7 @@ public:
for (size_t i=0; i<cellCount; i++) {
using LocalProblemType = typename std::remove_reference<decltype(*(this->localProblems_[i]))>::type;
- this->localProblems_[i] = std::make_shared<LocalProblemType>(*this->mat_, rhs, localToGlobal[i], boundaryDofs[i]);
+ this->localProblems_[i] = std::make_shared<LocalProblemType>(this->matrix_, rhs, localToGlobal[i], boundaryDofs[i]);
}
}
};
diff --git a/dune/faultnetworks/preconditioners/levelpatchpreconditioner.hh b/dune/faultnetworks/preconditioners/levelpatchpreconditioner.hh
index 3d9a9526ad7f09db9d610e7a35cf8ad96625679a..d651e55ef7fa3f22337b4cdf2e5668e8601d6acb 100644
--- a/dune/faultnetworks/preconditioners/levelpatchpreconditioner.hh
+++ b/dune/faultnetworks/preconditioners/levelpatchpreconditioner.hh
@@ -44,6 +44,7 @@ protected:
size_t patchDepth_;
BoundaryMode boundaryMode_;
+ MatrixType matrix_;
std::vector<std::shared_ptr<OscLocalProblem<MatrixType, VectorType>> > localProblems_;
public:
@@ -65,6 +66,7 @@ public:
setPatchDepth();
setBoundaryMode();
setDirection();
+ setup();
}
// has to setup localProblems_
@@ -121,6 +123,35 @@ public:
}
private:
+ void setup() {
+ // assemble stiffness matrix for entire level including all faults
+ GlobalFaultAssembler<BasisType, BasisType> globalFaultAssembler(basis_, basis_, levelInterfaceNetwork_);
+ globalFaultAssembler.assembleOperator(localAssembler_, localInterfaceAssemblers_, matrix_);
+
+ // set boundary conditions
+ Dune::BitSetVector<1> globalBoundaryDofs;
+ BoundaryPatch<GridView> boundaryPatch(levelInterfaceNetwork_.levelGridView(), true);
+ constructBoundaryDofs(boundaryPatch, basis_, globalBoundaryDofs);
+
+ typedef typename MatrixType::row_type RowType;
+ typedef typename RowType::ConstIterator ColumnIterator;
+
+ for(size_t i=0; i<globalBoundaryDofs.size(); i++) {
+ if(!globalBoundaryDofs[i][0])
+ continue;
+
+ RowType& row = matrix_[i];
+
+ ColumnIterator cIt = row.begin();
+ ColumnIterator cEndIt = row.end();
+
+ for(; cIt!=cEndIt; ++cIt) {
+ row[cIt.index()] = 0;
+ }
+ row[i] = 1;
+ }
+ }
+
void iterateAdd() {
VectorType it, x;
for (size_t i=0; i<localProblems_.size(); i++) {
diff --git a/dune/faultnetworks/preconditioners/multilevelpatchpreconditioner.hh b/dune/faultnetworks/preconditioners/multilevelpatchpreconditioner.hh
index 2f12d043cca61bef96032619c7d7be320341b462..459f3dcea43b9518bb3498fbcb43017242166b10 100644
--- a/dune/faultnetworks/preconditioners/multilevelpatchpreconditioner.hh
+++ b/dune/faultnetworks/preconditioners/multilevelpatchpreconditioner.hh
@@ -174,8 +174,7 @@ public:
ignoreNodesHierarchy_[i] = new BitVector();
mgTransfer_[i]->restrictToFathers(*ignoreNodesHierarchy_[i+1], *ignoreNodesHierarchy_[i]);
}
- } else
- DUNE_THROW(SolverError, "We need a set of nodes to ignore");
+ }
coarseGlobalPreconditioner_->setProblem(*matrixHierarchy_[0], *xHierarchy_[0], rhsHierarchy_[0]);
coarseGlobalPreconditioner_->preprocess();
@@ -237,18 +236,35 @@ private:
mgTransfer_[i]->restrict(rhsHierarchy_[i+1], rhsHierarchy_[i]);
}
- // coarseGlobalPreconditioner_->setProblem(*matrixHierarchy_[0], *xHierarchy_[0], rhsHierarchy_[0]);
+ //print(*matrixHierarchy_[0], "mat 0:");
+ //print(rhsHierarchy_[0], "rhs0:");
+ coarseGlobalPreconditioner_->setProblem(*matrixHierarchy_[0], *xHierarchy_[0], rhsHierarchy_[0]);
coarseGlobalPreconditioner_->iterate();
+ //print(*xHierarchy_[0], "x0:");
for (size_t i=0; i<levelPatchPreconditioners_.size(); i++) {
- // levelPatchPreconditioners_[i]->setProblem(*matrixHierarchy_[i+1], *xHierarchy_[i+1], rhsHierarchy_[i+1]);
+ //print(*matrixHierarchy_[i+1], "mat:");
+ //print(rhsHierarchy_[i+1], "rhs:");
+
+ levelPatchPreconditioners_[i]->setProblem(*matrixHierarchy_[i+1], *xHierarchy_[i+1], rhsHierarchy_[i+1]);
levelPatchPreconditioners_[i]->iterate();
+ //print(*xHierarchy_[i+1], "xi+1:");
+
VectorType x(xHierarchy_[i+1]->size());
mgTransfer_[i]->prolong(*xHierarchy_[i], x);
+ //print(x, "xi:");
+
*xHierarchy_[i+1] += x;
+ //print(*xHierarchy_[i+1], "xi+1:");
}
+
+
+ //print(*xHierarchy_[levelPatchPreconditioners_.size()] , "xback:");
+ //print(*xHierarchy_.back() , "xback:");
+ *(this->x_) = *xHierarchy_.back();
+ //print(*(this->x_) , "x:");
}
diff --git a/dune/faultnetworks/preconditioners/supportpatchpreconditioner.hh b/dune/faultnetworks/preconditioners/supportpatchpreconditioner.hh
index dc4064e81981890575975847ab157993829b3c0f..bf57590d9e665ad78752b49c635b2f4070dfadf7 100644
--- a/dune/faultnetworks/preconditioners/supportpatchpreconditioner.hh
+++ b/dune/faultnetworks/preconditioners/supportpatchpreconditioner.hh
@@ -60,11 +60,11 @@ public:
const auto& boundaryDofs = patchFactory.getBoundaryDofs();
VectorType rhs;
- rhs.resize(this->mat_->N());
+ rhs.resize(this->matrix_.N());
rhs = 0;
using LocalProblemType = typename std::remove_reference<decltype(*(this->localProblems_[i]))>::type;
- this->localProblems_[i] = std::make_shared<LocalProblemType>(*this->mat_, rhs, localToGlobal, boundaryDofs);
+ this->localProblems_[i] = std::make_shared<LocalProblemType>(this->matrix_, rhs, localToGlobal, boundaryDofs);
}
}
};
diff --git a/src/cantorfaultnetworks/cantorfaultnetwork.cc b/src/cantorfaultnetworks/cantorfaultnetwork.cc
index 90a9b2b9a97e37b4e5e625b4c1311e15ab5fc85d..4a65d2f2ded47798476c56b5846c5ad5f94103e5 100644
--- a/src/cantorfaultnetworks/cantorfaultnetwork.cc
+++ b/src/cantorfaultnetworks/cantorfaultnetwork.cc
@@ -141,7 +141,7 @@ int main(int argc, char** argv) { try
matrixReader.read(matrix);
const int oscGridN = matrix.N();
- const int oscGridLevelIdx = std::round(std::log2(oscGridN)) - coarseResolution;
+ const int oscGridLevelIdx = std::max(std::round(std::log2(oscGridN)) - coarseResolution, 0.0);
if (oscGridN>fineGridN)
DUNE_THROW(Dune::Exception, "Provided oscData too large!");
@@ -161,7 +161,7 @@ int main(int argc, char** argv) { try
CantorFaultFactory<GridType> faultFactory(levelToCantorLevel, coarseResolution, exactLevelIdx, maxCantorLevel);
const InterfaceNetwork<GridType>& interfaceNetwork = faultFactory.interfaceNetwork();
- if (problemCount==0) {
+ /*if (problemCount==0) {
std::vector<int> writeLevels {1, 2, 8};
for (size_t i=0; i<writeLevels.size(); i++) {
int writeLevel = writeLevels[i];
@@ -170,7 +170,7 @@ int main(int argc, char** argv) { try
LevelInterfaceNetworkWriter networkWriter(resultPath + "levelinterfacenetwork_" + std::to_string(writeLevel) + ".tikz");
networkWriter.write(interfaceNetwork.levelInterfaceNetwork(writeLevel), writeGrid);
}
- }
+ }*/
const GridType& grid = faultFactory.grid();
@@ -307,6 +307,14 @@ int main(int argc, char** argv) { try
}
}
+ // set initial iterate
+ VectorType initialX;
+ DGMGTransfer<DGBasis> coarseFineTransfer(*coarseBasis, *fineBasis);
+ coarseFineTransfer.prolong(coarseSol, initialX);
+
+ VectorType rhs(fineRhs);
+
+ // setup preconditioner
const auto& preconditionerParset = parameterSet.sub("preconditioner");
using MultilevelPatchPreconditioner = MultilevelPatchPreconditioner<DGBasis, LocalOscAssembler, LocalInterfaceAssembler, MatrixType, VectorType>;
@@ -317,37 +325,31 @@ int main(int argc, char** argv) { try
levelFaultPreconditioner.setPatchDepth(patchDepth);
levelFaultPreconditioner.setBoundaryMode(LevelPatchPreconditioner::BoundaryMode::homogeneous);
}*/
+ preconditioner.setProblem(fineGlobalAssembler.matrix(), initialX, rhs);
std::cout << "Setup complete, starting preconditioned cg iteration!" << std::endl;
std::cout << std::endl << std::endl;
- // set initial iterate
- VectorType initialX;
- DGMGTransfer<DGBasis> coarseFineTransfer(*coarseBasis, *fineBasis);
- coarseFineTransfer.prolong(coarseSol, initialX);
-
- VectorType rhs(fineRhs);
-
DGMGTransfer<DGBasis> discMGTransfer(*fineBasis, *exactBasis);
// solve
- /*OscCGSolver<MatrixType, VectorType, DGBasis>
+ OscCGSolver<MatrixType, VectorType, DGBasis>
solver(*fineBasis, &fineGlobalAssembler.matrix(), &initialX, &rhs, &preconditioner,
maxIterations, solverTolerance, &exactEnergyNorm,
- Solver::FULL, &exactSol, &fineSol, discMGTransfer, 1.0, true); //((oscGridN+0.0)/fineGridN) */
+ Solver::FULL, &exactSol, &fineSol, discMGTransfer, 1.0, true); //((oscGridN+0.0)/fineGridN)
//solver.historyBuffer_ = resultPath + "solutions/level_" + std::to_string(fineResolution) + ".vec";
- VectorType x = initialX;
+ /*VectorType x = initialX;
preconditioner.setProblem(fineGlobalAssembler.matrix(), x, rhs);
Dune::Solvers::LoopSolver<VectorType>
solver (preconditioner, maxIterations, solverTolerance,
fineEnergyNorm,
Solver::FULL,
- true, &fineSol);
+ true, &fineSol); */
solver.check();
solver.preprocess();
diff --git a/src/cantorfaultnetworks/cantorfaultnetwork.parset b/src/cantorfaultnetworks/cantorfaultnetwork.parset
index 0af54c6260ea293ff2e6e68bb3eb56d8c2c31744..981054b93d2ad51fbb9ebdd8cf9e5d65e26c3855 100644
--- a/src/cantorfaultnetworks/cantorfaultnetwork.parset
+++ b/src/cantorfaultnetworks/cantorfaultnetwork.parset
@@ -10,12 +10,12 @@ patchDepth = 1
###########################################
[problem0]
-oscDataFile = oscDataLaplace16.mat
+oscDataFile = oscDataLaplace4.mat
# level resolution in 2^(-...)
-coarseResolution = 3
-fineResolution = 4
-exactResolution = 5
+coarseResolution = 2
+fineResolution = 3
+exactResolution = 4
minCantorResolution = 0
penaltyFactor = 1
diff --git a/src/geofaultnetworks/rockfaultnetwork.cc b/src/geofaultnetworks/rockfaultnetwork.cc
index 777b2e11cbfbb32089f955dfd0dee024c8fbd301..b6983066fc8fb2e47575f8b963e09e92f96ed790 100644
--- a/src/geofaultnetworks/rockfaultnetwork.cc
+++ b/src/geofaultnetworks/rockfaultnetwork.cc
@@ -314,6 +314,14 @@ int main(int argc, char** argv) { try
}
}
+ // set initial iterate
+ VectorType initialX;
+ DGMGTransfer<DGBasis> coarseFineTransfer(*coarseBasis, *fineBasis);
+ coarseFineTransfer.prolong(coarseSol, initialX);
+
+ VectorType rhs(fineRhs);
+
+ // setup preconditioner
const auto& preconditionerParset = parameterSet.sub("preconditioner");
using MultilevelPatchPreconditioner = MultilevelPatchPreconditioner<DGBasis, LocalOscAssembler, LocalInterfaceAssembler, MatrixType, VectorType>;
@@ -324,16 +332,11 @@ int main(int argc, char** argv) { try
levelFaultPreconditioner.setPatchDepth(patchDepth);
levelFaultPreconditioner.setBoundaryMode(LevelPatchPreconditioner::BoundaryMode::homogeneous);
}*/
+ preconditioner.setProblem(fineGlobalAssembler.matrix(), initialX, rhs);
std::cout << "Setup complete, starting preconditioned cg iteration!" << std::endl;
std::cout << std::endl << std::endl;
- // set initial iterate
- VectorType initialX;
- DGMGTransfer<DGBasis> coarseFineTransfer(*coarseBasis, *fineBasis);
- coarseFineTransfer.prolong(coarseSol, initialX);
-
- VectorType rhs(fineRhs);
DGMGTransfer<DGBasis> discMGTransfer(*fineBasis, *exactBasis);