diff --git a/dune/faultnetworks/preconditioners/cellpatchpreconditioner.hh b/dune/faultnetworks/preconditioners/cellpatchpreconditioner.hh
index dae195f1fe60ef0d4ee0b1791edfcce8fed47ad4..a14f005352084a4c1dd86aa6a47c6fa60e08b2c0 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->matrix_.N());
+ rhs.resize(this->mat_->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->matrix_, rhs, localToGlobal[i], boundaryDofs[i]);
+ this->localProblems_[i] = std::make_shared<LocalProblemType>(*this->mat_, rhs, localToGlobal[i], boundaryDofs[i]);
}
}
};
diff --git a/dune/faultnetworks/preconditioners/levelpatchpreconditioner.hh b/dune/faultnetworks/preconditioners/levelpatchpreconditioner.hh
index d651e55ef7fa3f22337b4cdf2e5668e8601d6acb..3d9a9526ad7f09db9d610e7a35cf8ad96625679a 100644
--- a/dune/faultnetworks/preconditioners/levelpatchpreconditioner.hh
+++ b/dune/faultnetworks/preconditioners/levelpatchpreconditioner.hh
@@ -44,7 +44,6 @@ protected:
size_t patchDepth_;
BoundaryMode boundaryMode_;
- MatrixType matrix_;
std::vector<std::shared_ptr<OscLocalProblem<MatrixType, VectorType>> > localProblems_;
public:
@@ -66,7 +65,6 @@ public:
setPatchDepth();
setBoundaryMode();
setDirection();
- setup();
}
// has to setup localProblems_
@@ -123,35 +121,6 @@ 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 459f3dcea43b9518bb3498fbcb43017242166b10..2f12d043cca61bef96032619c7d7be320341b462 100644
--- a/dune/faultnetworks/preconditioners/multilevelpatchpreconditioner.hh
+++ b/dune/faultnetworks/preconditioners/multilevelpatchpreconditioner.hh
@@ -174,7 +174,8 @@ 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();
@@ -236,35 +237,18 @@ private:
mgTransfer_[i]->restrict(rhsHierarchy_[i+1], rhsHierarchy_[i]);
}
- //print(*matrixHierarchy_[0], "mat 0:");
- //print(rhsHierarchy_[0], "rhs0:");
- coarseGlobalPreconditioner_->setProblem(*matrixHierarchy_[0], *xHierarchy_[0], rhsHierarchy_[0]);
+ // coarseGlobalPreconditioner_->setProblem(*matrixHierarchy_[0], *xHierarchy_[0], rhsHierarchy_[0]);
coarseGlobalPreconditioner_->iterate();
- //print(*xHierarchy_[0], "x0:");
for (size_t i=0; i<levelPatchPreconditioners_.size(); i++) {
- //print(*matrixHierarchy_[i+1], "mat:");
- //print(rhsHierarchy_[i+1], "rhs:");
-
- levelPatchPreconditioners_[i]->setProblem(*matrixHierarchy_[i+1], *xHierarchy_[i+1], rhsHierarchy_[i+1]);
+ // 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 bf57590d9e665ad78752b49c635b2f4070dfadf7..dc4064e81981890575975847ab157993829b3c0f 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->matrix_.N());
+ rhs.resize(this->mat_->N());
rhs = 0;
using LocalProblemType = typename std::remove_reference<decltype(*(this->localProblems_[i]))>::type;
- this->localProblems_[i] = std::make_shared<LocalProblemType>(this->matrix_, rhs, localToGlobal, boundaryDofs);
+ this->localProblems_[i] = std::make_shared<LocalProblemType>(*this->mat_, rhs, localToGlobal, boundaryDofs);
}
}
};
diff --git a/src/cantorfaultnetworks/cantorfaultnetwork.cc b/src/cantorfaultnetworks/cantorfaultnetwork.cc
index 4a65d2f2ded47798476c56b5846c5ad5f94103e5..90a9b2b9a97e37b4e5e625b4c1311e15ab5fc85d 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::max(std::round(std::log2(oscGridN)) - coarseResolution, 0.0);
+ const int oscGridLevelIdx = std::round(std::log2(oscGridN)) - coarseResolution;
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,14 +307,6 @@ 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>;
@@ -325,31 +317,37 @@ 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 981054b93d2ad51fbb9ebdd8cf9e5d65e26c3855..0af54c6260ea293ff2e6e68bb3eb56d8c2c31744 100644
--- a/src/cantorfaultnetworks/cantorfaultnetwork.parset
+++ b/src/cantorfaultnetworks/cantorfaultnetwork.parset
@@ -10,12 +10,12 @@ patchDepth = 1
###########################################
[problem0]
-oscDataFile = oscDataLaplace4.mat
+oscDataFile = oscDataLaplace16.mat
# level resolution in 2^(-...)
-coarseResolution = 2
-fineResolution = 3
-exactResolution = 4
+coarseResolution = 3
+fineResolution = 4
+exactResolution = 5
minCantorResolution = 0
penaltyFactor = 1
diff --git a/src/geofaultnetworks/rockfaultnetwork.cc b/src/geofaultnetworks/rockfaultnetwork.cc
index b6983066fc8fb2e47575f8b963e09e92f96ed790..777b2e11cbfbb32089f955dfd0dee024c8fbd301 100644
--- a/src/geofaultnetworks/rockfaultnetwork.cc
+++ b/src/geofaultnetworks/rockfaultnetwork.cc
@@ -314,14 +314,6 @@ 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>;
@@ -332,11 +324,16 @@ 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);