diff --git a/dune/faultnetworks/faultfactories/CMakeLists.txt b/dune/faultnetworks/faultfactories/CMakeLists.txt
index b1baf5e789515893ec4ae134d84a4ca65d83911e..998ed92d5aba3ad6cd6d299bdabc29cc1be55150 100644
--- a/dune/faultnetworks/faultfactories/CMakeLists.txt
+++ b/dune/faultnetworks/faultfactories/CMakeLists.txt
@@ -4,6 +4,7 @@ add_custom_target(faultnetworks_faultfactories_src SOURCES
geofaultfactory.hh
cantorconvergencefaultfactory.hh
cantorfaultfactory.hh
+ rockfaultfactory.hh
sparsecantorfaultfactory.hh
)
@@ -12,5 +13,6 @@ install(FILES
geofaultfactory.hh
cantorconvergencefaultfactory.hh
cantorfaultfactory.hh
+ rockfaultfactory.hh
sparsecantorfaultfactory.hh
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/faultnetworks/faultfactories)
diff --git a/dune/faultnetworks/faultfactories/rockfaultfactory.hh b/dune/faultnetworks/faultfactories/rockfaultfactory.hh
new file mode 100644
index 0000000000000000000000000000000000000000..e8f216e60b7295bc40f132fe125c9f5a7d3efcad
--- /dev/null
+++ b/dune/faultnetworks/faultfactories/rockfaultfactory.hh
@@ -0,0 +1,617 @@
+#ifndef ROCK_FAULT_FACTORY_HH
+#define ROCK_FAULT_FACTORY_HH
+
+#include <math.h>
+#include <queue>
+
+#include <dune/faultnetworks/facehierarchy.hh>
+#include <dune/faultnetworks/utils/debugutils.hh>
+#include <dune/faultnetworks/faultfactories/oscunitcube.hh>
+#include <dune/faultnetworks/levelinterfacenetwork.hh>
+#include <dune/faultnetworks/interfacenetwork.hh>
+#include <dune/faultnetworks/faultinterface.hh>
+#include <dune/faultnetworks/hierarchicleveliterator.hh>
+
+#include <dune/grid/common/mcmgmapper.hh>
+#include <dune/grid/io/file/vtk/vtkwriter.hh>
+#include <dune/grid/uggrid.hh>
+
+#include <dune/fufem/boundarypatch.hh>
+#include <dune/fufem/boundaryiterator.hh>
+
+template<class ctype = double> class Rock {
+protected:
+ int level_;
+
+ // vertex IDs
+ int left_;
+ int right_;
+ int top_;
+ int bottom_;
+
+ Dune::FieldVector<ctype, 2> center_;
+public:
+ void set(int level, int left, int right, int top, int bottom, const Dune::FieldVector<ctype, 2>& center) {
+ level_ = level;
+
+ left_ = left;
+ right_ = right;
+ top_ = top;
+ bottom_ = bottom;
+
+ center_ = center;
+ }
+
+ auto level() const {
+ return level_;
+ }
+
+ auto left() const {
+ return left_;
+ }
+
+ auto right() const {
+ return right_;
+ }
+
+ auto top() const {
+ return top_;
+ }
+
+ auto bottom() const {
+ return bottom_;
+ }
+
+ const auto& center() {
+ return center_;
+ }
+};
+
+
+template <class GridType>
+class LevelRockFaultFactory
+{
+ //! Parameter for mapper class
+ template<int dim>
+ struct FaceMapperLayout
+ {
+ bool contains (Dune::GeometryType gt)
+ {
+ return gt.dim() == dim-1;
+ }
+ };
+
+protected:
+ static const int dimworld = GridType::dimensionworld;
+ static const int dim = GridType::dimension;
+ using ctype = typename GridType::ctype;
+
+ using Coords = typename Dune::FieldVector<ctype, dimworld>;
+ using GV = typename GridType::LevelGridView;
+ using Intersection = typename GridType::LevelIntersection;
+
+ using Element = typename GridType::template Codim<0>::Entity;
+ static const int dimElement = Element::dimension;
+
+ using FaceMapper = typename Dune::MultipleCodimMultipleGeomTypeMapper<GV, FaceMapperLayout >;
+
+
+ using Rock = Rock<ctype>;
+
+ const std::shared_ptr<GridType> grid_;
+ const int level_;
+ const ctype resolution_;
+ const GV gridView_;
+ const double maxAngle_;
+
+ const typename GV::IndexSet& indexSet_;
+
+ FaceMapper faceMapper_;
+
+ std::vector<Intersection> faces_;
+ std::vector<Coords> vertexPositions_;
+ std::vector<std::vector<size_t>> vertexToFaces_;
+ std::vector<int> coarseToLevelVertex_;
+
+ const LevelRockFaultFactory& coarseLevelFactory_;
+ std::vector<Rock> rocks_;
+ std::vector<std::shared_ptr<FaultInterface<GV>>>& faults_;
+
+private:
+ template <typename T> int sgn(T val) {
+ return (T(0) < val) - (val < T(0));
+ }
+
+ bool intersectionAllowed(const size_t faceID, const size_t vertexID,
+ const std::array<size_t, 2>& centerIDs,
+ const std::set<int>& separatingIDs, const std::set<size_t>& faultDofs,
+ const Coords& direction, size_t dim) const {
+
+ const auto& intersection = faces_[faceID];
+
+ //check if "back" edge, cannot deviate from desiredOrientation more than maxAngle degrees (radian)
+ auto orientation = intersection.geometry().center();
+ orientation -= vertexPositions_[vertexID];
+ orientation /= orientation.two_norm();
+
+ if (std::acos(direction*orientation) > maxAngle_)
+ return false;
+
+ // check if intersection has separating dofs or other fault dofs
+ std::set<size_t> intersectionDofs;
+ computeIntersectionDofs(indexSet_, intersection, intersectionDofs);
+ intersectionDofs.erase(vertexID);
+
+ const size_t otherDim = (dim + 1) % 2;
+ for (const auto& isDof : intersectionDofs) {
+ const auto& vertex = vertexPositions_[isDof];
+ std::array<size_t, 2> IDs = {computeID(vertex, 0), computeID(vertex, 1)};
+
+ bool centerPassed = (sgn(direction[dim])>0) ? (centerIDs[dim]<IDs[dim]) : (centerIDs[dim]>IDs[dim]);
+ if (faultDofs.count(isDof) or separatingIDs.count(IDs[otherDim]) or centerPassed) {
+ return false;
+ }
+ }
+
+ return true;
+ }
+
+ /*
+ ctype distance(const Intersection& isec, const Dune::FieldVector<ctype, dimworld> & vertex) const {
+ Dune::FieldVector<ctype, dimworld> vec(vertex);
+ vec += isec.center();
+
+ return isec.unitOuterNormal()*vec;
+ }*/
+
+ //works only in 2D, vertex(1) = vertex(0) - yid, intersection of line (given by vertex and derivative 1) with y axis
+ int computeID(const Coords& vertex, int dim) const {
+ return (int) (vertex[dim]*resolution_);
+ }
+
+ void computeIntersectionDofs(const Intersection& intersection, std::set<size_t>& intersectionDofs) const {
+ intersectionDofs.clear();
+
+ // loop over all vertices of the intersection
+ const auto& insideElement = intersection.inside();
+ const auto& refElement = Dune::ReferenceElements<double,dimElement>::general(insideElement.type());
+
+ for (int i=0; i<refElement.size(intersection.indexInInside(), 1, dimElement); i++) {
+ size_t idxInElement = refElement.subEntity(intersection.indexInInside(), 1, i, dimElement);
+ size_t globalIdx = indexSet_.subIndex(insideElement, idxInElement, dimElement);
+
+ intersectionDofs.insert(globalIdx);
+ }
+ }
+
+ bool isTargetBoundaryVertex(size_t vertexIdx, const std::array<size_t, 2>& targetIDs) const {
+ const auto& vertex = vertexPositions_[vertexIdx];
+ const std::array<size_t, 2> vertexIDs = {computeID(vertex, 0), computeID(vertex, 1)};
+ return (vertexIDs[0] == targetIDs[0]) and (vertexIDs[1] == targetIDs[1]);
+ }
+
+ typedef std::vector<double>::const_iterator ConstVectorIterator;
+
+ struct ordering {
+ bool operator ()(std::pair<size_t, ConstVectorIterator> const& a, std::pair<size_t, ConstVectorIterator> const& b) {
+ return *(a.second) < *(b.second);
+ }
+ };
+
+ void distanceSort(std::vector<size_t>& admissibleFaces, std::vector<double>& distances) const {
+ std::vector<std::pair<size_t, ConstVectorIterator>> order(distances.size());
+
+ size_t j = 0;
+ ConstVectorIterator it = distances.begin();
+ ConstVectorIterator itEnd = distances.end();
+ for (; it != itEnd; ++it, ++j)
+ order[j] = std::make_pair(j, it);
+
+ sort(order.begin(), order.end(), ordering());
+
+ std::vector<size_t> initialAdmissibleFaces(admissibleFaces);
+ for (size_t i=0; i<admissibleFaces.size(); i++)
+ admissibleFaces[i] = initialAdmissibleFaces[order[i].first];
+ }
+
+ /* void generateFaultSeeds(std::vector<FaultCorridor>& faultCorridors, std::vector<std::vector<size_t>>& faultSeeds) {
+
+ //std::cout << "LevelGeoFaultFactory::generateFaultSeeds() " << std::endl;
+ //std::cout << "------------------------------------- " << std::endl << std::endl;
+
+ faultSeeds.resize(faultCorridors.size());
+
+ std::map<int, std::pair<size_t,size_t>> yIDtoFaultCorridor;
+ for (size_t i=0; i<faultCorridors.size(); i++){
+ FaultCorridor& faultCorridor = faultCorridors[i];
+ const std::vector<int> & faultYs = faultCorridor.faultYs();
+
+ faultSeeds[i].resize(faultYs.size());
+ for (size_t j=0; j<faultYs.size(); j++) {
+ yIDtoFaultCorridor[faultYs[j]] = std::make_pair(i, j);
+ }
+ }
+
+ BoundaryIterator<GV> bIt(gridView_, BoundaryIterator<GV>::begin);
+ BoundaryIterator<GV> bEnd(gridView_, BoundaryIterator<GV>::end);
+ for(; bIt!=bEnd; ++bIt) {
+
+ const Element& insideElement = bIt->inside();
+ const auto& geometry = insideElement.geometry();
+
+ // neglect "upper" part of boundary, where (y==1 and x>0) or (y>0 and x==1)
+ const auto& center = bIt->geometry().center();
+ if (center[0] + center[1] >1)
+ continue;
+
+ const auto& refElement = Dune::ReferenceElements<double,dimElement>::general(insideElement.type());
+
+ for (int i=0; i<refElement.size(bIt->indexInInside(), 1, dimElement); i++) {
+ size_t idxInElement = refElement.subEntity(bIt->indexInInside(), 1, i, dimElement);
+ size_t globalIdx = indexSet_.subIndex(insideElement, idxInElement, dimElement);
+
+ const auto& vertex = geometry.corner(idxInElement);
+
+ const int yID = computeYID(vertex);
+
+ if (yIDtoFaultCorridor.count(yID)) {
+ const std::pair<size_t, size_t>& indices = yIDtoFaultCorridor[yID];
+ faultSeeds[indices.first][indices.second] = globalIdx;
+ }
+ }
+ }
+
+ //std::cout << "------------------------------------- " << std::endl << std::endl;
+ }*/
+
+
+ void generateFault(std::shared_ptr<FaultInterface<GV>> fault, const size_t faultSeedIdx, const Coords& center,
+ const size_t corridor,
+ const std::set<size_t>& constSeparatingDofs) {
+ //std::cout << "LevelGeoFaultFactory::generateFault() " << std::endl;
+ //std::cout << "------------------------------------- " << std::endl << std::endl;
+
+ bool success = false;
+
+ std::set<size_t> separatingDofs = constSeparatingDofs;
+
+ auto direction = center;
+ direction -= vertexPositions_[faultSeedIdx];
+ direction /= direction.two_norm();
+
+ size_t dim = 0;
+ if (std::abs(direction[0]) < std::abs(direction[1])) {
+ dim = 1;
+ }
+
+ const std::array<size_t, 2> centerIDs = {computeID(center, 0), computeID(center, 1)};
+
+ std::set<size_t> separatingIDs;
+
+ auto faultSeedID = computeID(vertexPositions_[faultSeedIdx], dim);
+ separatingIDs.insert(faultSeedID + corridor);
+ separatingIDs.insert(faultSeedID - corridor);
+
+ std::vector<size_t> faultDofs(0);
+ faultDofs.push_back(faultSeedIdx);
+
+ std::vector<size_t> faultFaces;
+
+ std::map<size_t, std::vector<size_t>> vertexToAdmissibleFaces;
+
+ std::queue<size_t> vertexQueue;
+ vertexQueue.push(faultSeedIdx);
+
+ while (!vertexQueue.empty()) {
+ const size_t vertexID = vertexQueue.front();
+ vertexQueue.pop();
+
+ if (isTargetBoundaryVertex(vertexID, centerIDs)) {
+ success = true;
+ break;
+ }
+
+ if (vertexToAdmissibleFaces.count(vertexID)==0) {
+ const std::vector<size_t>& faces = vertexToFaces_[vertexID];
+ std::vector<size_t> admissibleFaces;
+
+ for (size_t i=0; i<faces.size(); i++) {
+ if (intersectionAllowed(faces[i], vertexID, centerIDs, separatingDofs, separatingIDs, faultDofs, direction, dim)) {
+ admissibleFaces.push_back(faces[i]);
+ }
+ }
+
+ std::vector<double> distances(admissibleFaces.size(), 0);
+ for (size_t i=0; i<admissibleFaces.size(); i++) {
+ Coords vec(center);
+ vec -= faces_[admissibleFaces[i]].geometry().center();
+ distances[i] = vec.two_norm();
+ }
+ distanceSort(admissibleFaces, distances);
+
+ vertexToAdmissibleFaces[vertexID] = admissibleFaces;
+ }
+
+ std::vector<size_t>& suitableFaces = vertexToAdmissibleFaces[vertexID];
+
+ if (suitableFaces.size()==0) {
+ faultDofs.pop_back();
+ faultFaces.pop_back();
+ vertexToAdmissibleFaces.erase(vertexID);
+
+ separatingDofs.insert(vertexID);
+ vertexQueue.push(faultDofs.back());
+ } else {
+
+ // generate random number from (0,1)
+ double randomNumber = ((double) std::rand() / (RAND_MAX));
+
+ size_t randSelector = (size_t) std::abs(std::log2(randomNumber));
+ if (randSelector >= suitableFaces.size()) {
+ randSelector = suitableFaces.size() -1;
+ }
+
+ size_t nextFaceID = suitableFaces[randSelector];
+
+ suitableFaces.erase(suitableFaces.begin()+randSelector);
+
+ std::set<size_t> nextFaceDofs;
+ computeIntersectionDofs(indexSet_, faces_[nextFaceID], nextFaceDofs);
+ nextFaceDofs.erase(vertexID);
+
+ faultDofs.push_back(*nextFaceDofs.begin());
+ faultFaces.push_back(nextFaceID);
+
+ size_t nextVertexID = *(nextFaceDofs.begin());
+ vertexQueue.push(nextVertexID);
+ }
+ }
+
+ if (!success) {
+ std::cout << "Generating a fault failed: Unable to reach target boundary! This should not happend!" << std::endl;
+ DUNE_THROW(Dune::Exception, "Generating a fault failed: Unable to reach target boundary!");
+ }
+
+ for (size_t i=0; i<faultFaces.size(); i++) {
+ fault->addFace(faces_[faultFaces[i]]);
+ }
+
+ //std::cout << "------------------------------------- " << std::endl << std::endl;
+ }
+
+ void prolong(const Rock& rock, Rock& newRock) {
+ auto newLeft = coarseToLevelVertex_[rock.left()];
+ auto newRight = coarseToLevelVertex_[rock.right()];
+ auto newTop = coarseToLevelVertex_[rock.top()];
+ auto newBottom = coarseToLevelVertex_[rock.bottom()];
+
+ newRock.set(rock.level(), newLeft, newRight, newTop, newBottom, rock.center());
+ }
+
+ bool randomSplit(const Rock& rock) {
+ const auto& center = rock.center();
+
+ double res = 0.0;
+ if (center[0] < center[1]) {
+ res = 100.0*(1.0 - center[0]);
+ } else {
+ res = 100.0*(1.0 - center[1]);
+ }
+
+ double prob = (100.0 * std::rand() / (RAND_MAX + 1.0)) + 1;
+ return res > prob;
+ }
+
+ void split(const Rock& rock, const std::set<size_t>& separatingDofs) {
+ Rock newRock;
+ prolong(rock, newRock);
+
+ bool toBeSplit = (rock.level() == level_-1) and randomSplit(rock);
+ if (!toBeSplit) {
+ rocks_.push_back(newRock);
+ } else {
+ // split rock into 4 subparts by two new intersecting faults
+ std::shared_ptr<FaultInterface<GV>> xFault = std::make_shared<FaultInterface<GV>>(gridView_, level_);
+ generateFault(xFault, newRock, separatingDofs, 0);
+ faults_.push_back(xFault);
+
+ std::shared_ptr<FaultInterface<GV>> yFault = std::make_shared<FaultInterface<GV>>(gridView_, level_);
+ generateFault(yFault, newRock, separatingDofs, 1);
+ faults_.push_back(yFault);
+ }
+ }
+
+ public:
+ LevelRockFaultFactory(const std::shared_ptr<GridType> grid, const int level, const ctype resolution) :
+ grid_(grid),
+ level_(level),
+ resolution_(resolution),
+ gridView_(grid_->levelGridView(level_)),
+ indexSet_(gridView_.indexSet()),
+ faceMapper_(gridView_)
+ {
+ // init faces_, vertexPositions_, vertexToFaces_
+ faces_.resize(faceMapper_.size());
+ vertexPositions_.resize(gridView_.size(dim));
+ vertexToFaces_.resize(gridView_.size(dim));
+ coarseToLevelVertex_.resize(coarseLevelFactory_.gridView().size(dim));
+
+ std::vector<bool> faceHandled(faceMapper_.size(), false);
+
+ for (const auto& elem:elements(gridView_)) {
+ for (const auto& isect:intersections(gridView_, elem)) {
+ //const size_t faceID = faceMapper_.index(isect);
+
+ const int faceID = faceMapper_.subIndex(elem, isect.indexInInside(), 1);
+
+ if (isect.boundary())
+ continue;
+
+ if (faceHandled[faceID])
+ continue;
+
+ faceHandled[faceID] = true;
+ faces_[faceID] = isect;
+
+ const auto& refElement = Dune::ReferenceElements<double,dimElement>::general(elem.type());
+
+ for (int i=0; i<refElement.size(isect.indexInInside(), 1, dimElement); i++) {
+ //size_t vertexID = idxSet.subIndex(elem.subEntity<1>(isect), i, dimElement);
+
+ size_t idxInElement = refElement.subEntity(isect.indexInInside(), 1, i, dimElement);
+ const size_t vertexID = indexSet_.subIndex(elem, idxInElement, dimElement);
+
+ vertexPositions_[vertexID] = elem.geometry().corner(idxInElement);
+ vertexToFaces_[vertexID].push_back(faceID);
+ }
+ }
+ }
+
+ for (const auto& vertex : vertices(gridView_)) {
+ size_t coarseVertexID;
+ bool isCoarseVertex = coarseLevelFactory_.indexSet().index(vertex, coarseVertexID);
+
+ if (isCoarseVertex) {
+ coarseToLevelVertex_[coarseVertexID] = indexSet_.index(vertex);
+ }
+ }
+ }
+
+ void build(const std::set<size_t>& separatingDofs) {
+ faults_.resize(0);
+ const auto& coarseRocks = coarseLevelFactory_.rocks();
+
+ for (size_t i=0; i<coarseRocks.size(); i++) {
+ split(coarseRocks[i], separatingDofs);
+ }
+ }
+
+ const auto& indexSet() const {
+ return indexSet_;
+ }
+
+ const auto& rocks() const {
+ return rocks_;
+ }
+};
+
+template <class GridType>
+class RockFaultFactory {
+
+private:
+ const int coarseResolution_;
+ const size_t maxLevel_;
+ const int coarseGridN_;
+
+ std::vector<double> levelResolutions_;
+
+ std::shared_ptr<GridType> grid_;
+ std::shared_ptr<InterfaceNetwork<GridType>> interfaceNetwork_;
+
+ std::vector<std::shared_ptr<LevelRockFaultFactory<GridType>>> levelRockFaultFactories_;
+
+/*typename std::enable_if<!std::numeric_limits<ctype>::is_integer, bool>::type
+ almost_equal(ctype x, ctype y, int ulp) const
+ {
+ return std::abs(x-y) < std::numeric_limits<ctype>::epsilon() * std::abs(x+y) * ulp
+ || std::abs(x-y) < std::numeric_limits<ctype>::min();
+ }
+
+
+
+bool isSeparatingIntersection(const Intersection& intersection, const std::set<int>& separatingYIDs) const {
+ Dune::FieldVector<ctype, dimworld> faceCenter = intersection.geometry().center();
+
+ ctype yID = faceCenter[1]*coarseGridN_;
+ ctype intpart;
+
+ if (almost_equal(std::modf(yID, &intpart), 0.0, 2)) {
+ std::set<int>::iterator it = separatingYIDs.find((int) intpart);
+ if (it!=separatingYIDs.end())
+ return true;
+ else
+ return false;
+ } else
+ return false;
+}*/
+
+
+public:
+ //setup
+ RockFaultFactory(const int coarseResolution, const size_t maxLevel, const double maxAngle = 2) :
+ coarseResolution_(coarseResolution),
+ maxLevel_(maxLevel),
+ coarseGridN_(std::pow(2, coarseResolution_)),
+ interfaceNetwork_(nullptr)
+ {
+ using GridOb = OscUnitCube<GridType, 2>;
+ using GV = typename GridType::LevelGridView;
+
+ Dune::UGGrid<GridType::dimension>::setDefaultHeapSize(4000);
+ GridOb unitCube(coarseGridN_);
+ grid_ = unitCube.grid();
+ grid_->globalRefine(maxLevel_);
+
+ levelResolutions_.resize(maxLevel_+1);
+ levelRockFaultFactories_.resize(maxLevel_+1);
+
+ // init interface network
+ interfaceNetwork_ = std::make_shared<InterfaceNetwork<GridType>>(*grid_);
+
+ for (size_t i=0; i<=maxLevel_; i++) {
+ std::vector<std::shared_ptr<FaultInterface<GV>>> faults;
+
+ levelResolutions_[i] = std::pow(2, coarseResolution_+i);
+
+ //generate faults on level
+ levelRockFaultFactories_[i] = std::make_shared<LevelRockFaultFactory<GridType>>(grid_, i, levelResolutions_[i]);
+ levelRockFaultFactories_[i]->build(faults, maxAngle);
+
+ for (size_t j=0; j<faults.size(); j++) {
+ interfaceNetwork_->addInterface(faults[j]);
+ }
+
+ if (i==maxLevel_)
+ continue;
+
+ interfaceNetwork_->prolongLevel(i, i+1);
+ }
+ }
+
+ /*
+ void prolongToAll() {
+ // prolong all faults to all subsequent levels
+ for (int i=maxLevel_-1; i>=0; i--) {
+ if (interfaceNetwork_->size(i)>0) {
+ std::set<int> toLevels;
+ for (size_t j=i+1; j<=maxLevel_; j++) {
+ toLevels.insert(j);
+ }
+
+ interfaceNetwork_->prolongLevelInterfaces(i, toLevels);
+ }
+ }
+ interfaceNetwork_->build();
+ }*/
+
+ /*void prolongToAll() {
+ // prolong all faults to all subsequent levels
+ for (int i=interfaceNetwork_->size()-1; i>=0; i--) {
+ interfaceNetwork_->prolongLevelInterfaces(i, maxLevel_);
+ }
+ interfaceNetwork_->build();
+ }*/
+
+ const GridType& grid() const {
+ return *grid_;
+ }
+
+ /*const InterfaceNetwork<GridType>& interfaceNetwork() {
+ return *interfaceNetwork_;
+ }*/
+
+ InterfaceNetwork<GridType>& interfaceNetwork() {
+ return *interfaceNetwork_;
+ }
+};
+#endif
diff --git a/dune/faultnetworks/localproblem.hh b/dune/faultnetworks/localproblem.hh
index 5f115f88ad1ef2f3491def2eaf5d653b24c9c984..b4af08e91fa2108a7d5c62109b816465570a81a1 100644
--- a/dune/faultnetworks/localproblem.hh
+++ b/dune/faultnetworks/localproblem.hh
@@ -147,6 +147,10 @@ public:
void updateLocalRhs(const RangeType& rhs){
localRhs = rhs;
+ /*for (size_t i=0; i<localRhs.size(); i++) {
+ if (!boundaryDofs_[i][0])
+ localRhs[i] = rhs[i];
+ }*/
}
void updateRhsAndBoundary(const RangeType& rhs, const RangeType& boundary) {
@@ -160,6 +164,8 @@ public:
}
void solve(DomainType& x){
+ //print(boundaryDofs_, "boundaryDofs:");
+
#if HAVE_UMFPACK
RangeType localRhsCopy(localRhs);
Dune::InverseOperatorResult res;
diff --git a/dune/faultnetworks/preconditioners/levelpatchpreconditioner.hh b/dune/faultnetworks/preconditioners/levelpatchpreconditioner.hh
index c43f777356a033483741e47b55076aeef54081d8..c9fb33fdb8fce847ece8f595310aa741c60733e6 100644
--- a/dune/faultnetworks/preconditioners/levelpatchpreconditioner.hh
+++ b/dune/faultnetworks/preconditioners/levelpatchpreconditioner.hh
@@ -182,19 +182,46 @@ private:
auto& localProblem = *localProblems_[i];
- VectorType v;
+ VectorType v, v1, v2, v3;
localProblem.restrict(*(this->x_), v);
- VectorType r;
- localProblem.restrict(*(this->rhs_), r);
- Dune::MatrixVector::subtractProduct(r, localProblem.getLocalMat(), v);
- localProblem.updateLocalRhs(r);
+ v1.resize(localProblem.getLocalMat().N());
+ localProblem.getLocalMat().mv(v, v1);
+ v3.resize(matrix_.N());
+ matrix_.mv(*(this->x_), v3);
+ localProblem.restrict(v3, v2);
+
+ //print(v1, "v1");
+ //print(v2, "v2");
+
+ VectorType r1;
+ localProblem.restrict(*(this->rhs_), r1);
+ Dune::MatrixVector::subtractProduct(r1, localProblem.getLocalMat(), v);
+
+ /*localProblem.updateLocalRhs(r);
VectorType x;
localProblem.solve(x);
v += x;
- localProblem.updateVector(v, *(this->x_));
+ localProblem.updateVector(v, *(this->x_));*/
+
+ VectorType localR;
+ VectorType r = *(this->rhs_);
+ Dune::MatrixVector::subtractProduct(r, matrix_, *(this->x_));
+ localProblem.restrict(r, localR);
+
+ // print(r1, "r1");
+ // print(localR, "localR");
+
+ localProblem.updateLocalRhs(localR);
+
+ VectorType x;
+ localProblem.solve(x);
+
+ //VectorType v;
+ localProblem.prolong(x, v);
+ *(this->x_) += v;
}
};
diff --git a/src/cantorfaultnetworks/cantorfaultnetwork.cc b/src/cantorfaultnetworks/cantorfaultnetwork.cc
index 1dd5701b135cfe3df444a169319826a0a22c9ef7..e2ed377d029c7229d6f8ff5125a8c3b088b26ecf 100644
--- a/src/cantorfaultnetworks/cantorfaultnetwork.cc
+++ b/src/cantorfaultnetworks/cantorfaultnetwork.cc
@@ -60,6 +60,8 @@
#include <dune/faultnetworks/faultfactories/cantorfaultfactory.hh>
+#include <dune/solvers/solvers/loopsolver.hh>
+
int main(int argc, char** argv) { try
{
MPIHelper::instance(argc, argv);
@@ -86,7 +88,7 @@ int main(int argc, char** argv) { try
typedef OscLocalAssembler<GridType, DGFE, DGFE> LocalOscAssembler;
typedef IntersectionJumpMassAssembler<GV, DGFE, DGFE> LocalInterfaceAssembler;
-
+
// parse parameter file
ParameterTree parameterSet;
if (argc==2)
@@ -110,7 +112,6 @@ int main(int argc, char** argv) { try
const int minCantorResolution = problemParameters.get<int>("minCantorResolution");
const double c = problemParameters.get<double>("penaltyFactor");
- const int patchDepth = problemParameters.get<int>("patchDepth");
const int maxIterations = problemParameters.get<int>("maxIterations");
const double solverTolerance = problemParameters.get<double>("tol");
@@ -331,13 +332,23 @@ int main(int argc, char** argv) { try
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;
+ preconditioner.setProblem(fineGlobalAssembler.matrix(), x, rhs);
- solver.historyBuffer_ = resultPath + "solutions/level_" + std::to_string(fineResolution) + ".vec";
+ Dune::Solvers::LoopSolver<VectorType>
+ solver (preconditioner, maxIterations, solverTolerance,
+ fineEnergyNorm,
+ Solver::FULL,
+ true, &fineSol);
solver.check();
solver.preprocess();
diff --git a/src/cantorfaultnetworks/cantorfaultnetwork.parset b/src/cantorfaultnetworks/cantorfaultnetwork.parset
index 0d8411bb631bf9e36af423302f000bd1d9312dc8..fe8182cf83ad0868dc2354c936de20b271a04bac 100644
--- a/src/cantorfaultnetworks/cantorfaultnetwork.parset
+++ b/src/cantorfaultnetworks/cantorfaultnetwork.parset
@@ -1,6 +1,11 @@
path = ../data/
resultPath = ../cantorfaultnetworks/results/
+[preconditioner]
+patch = SUPPORT # CELL , SUPPORT
+mode = MULTIPLICATIVE # ADDITIVE , MULTIPLICATIVE
+multDirection = SYMMETRIC # SYMMETRIC , FORWARD , BACKWARD
+patchDepth = 1
###########################################
@@ -10,88 +15,10 @@ oscDataFile = oscDataLaplace32.mat
# level resolution in 2^(-...)
coarseResolution = 0
fineResolution = 5
-exactResolution = 10
-minCantorResolution = 0
-
-penaltyFactor = 1
-patchDepth = 1
-
-# cg solver
-nestedIteration = 0
-tol = 1e-12
-maxIterations = 8
-
-###########################################
-
-[problem1]
-oscDataFile = oscDataLaplace32.mat
-
-# level resolution in 2^(-...)
-coarseResolution = 0
-fineResolution = 6
-exactResolution = 10
+exactResolution = 8
minCantorResolution = 0
penaltyFactor = 1
-patchDepth = 1
-
-# cg solver
-nestedIteration = 0
-tol = 1e-12
-maxIterations = 8
-
-###########################################
-
-[problem2]
-oscDataFile = oscDataLaplace32.mat
-
-# level resolution in 2^(-...)
-coarseResolution = 0
-fineResolution = 7
-exactResolution = 10
-minCantorResolution = 0
-
-penaltyFactor = 1
-patchDepth = 1
-
-# cg solver
-nestedIteration = 0
-tol = 1e-12
-maxIterations = 8
-
-###########################################
-
-[problem3]
-oscDataFile = oscDataLaplace32.mat
-
-# level resolution in 2^(-...)
-coarseResolution = 0
-fineResolution = 8
-exactResolution = 10
-minCantorResolution = 0
-
-penaltyFactor = 1
-patchDepth = 1
-
-# cg solver
-nestedIteration = 0
-tol = 1e-12
-maxIterations = 8
-
-###########################################
-
-
-[problem4]
-oscDataFile = oscDataLaplace32.mat
-
-# level resolution in 2^(-...)
-coarseResolution = 0
-fineResolution = 9
-exactResolution = 10
-minCantorResolution = 0
-
-penaltyFactor = 1
-patchDepth = 1
# cg solver
nestedIteration = 0
diff --git a/src/cantorfaultnetworks/sparsecantorfaultnetwork.cc b/src/cantorfaultnetworks/sparsecantorfaultnetwork.cc
index c0918ab598da0701b34eb5ef4f91a428f0692685..e32efc7b023c37e50fd6ea495fb2872c0a37d500 100644
--- a/src/cantorfaultnetworks/sparsecantorfaultnetwork.cc
+++ b/src/cantorfaultnetworks/sparsecantorfaultnetwork.cc
@@ -59,6 +59,9 @@
#include <dune/faultnetworks/dgmgtransfer.hh>
#include <dune/faultnetworks/faultfactories/sparsecantorfaultfactory.hh>
+#include <dune/faultnetworks/faultfactories/rockfaultfactory.hh>
+
+#include <dune/solvers/solvers/loopsolver.hh>
int main(int argc, char** argv) { try
{
@@ -332,6 +335,15 @@ int main(int argc, char** argv) { try
maxIterations, solverTolerance, &exactEnergyNorm,
Solver::FULL, &exactSol, &fineSol, discMGTransfer, 1.0, true); //((oscGridN+0.0)/fineGridN)
+ /* VectorType x = initialX;
+ preconditioner.setProblem(fineGlobalAssembler.matrix(), x, rhs);
+
+ Dune::Solvers::LoopSolver<VectorType>
+ solver (preconditioner, maxIterations, solverTolerance,
+ fineEnergyNorm,
+ Solver::FULL,
+ true, &fineSol);*/
+
//solver.historyBuffer_ = resultPath + "solutions/level_" + std::to_string(fineCantorLevel) + ".vec";
diff --git a/src/cantorfaultnetworks/sparsecantorfaultnetwork.parset b/src/cantorfaultnetworks/sparsecantorfaultnetwork.parset
index 3b3ea31b5fe0c7cd57ccaad38502e50f9d5090bb..456c5bfc2935601ddc4e07e80fd4070fd2f9aff7 100644
--- a/src/cantorfaultnetworks/sparsecantorfaultnetwork.parset
+++ b/src/cantorfaultnetworks/sparsecantorfaultnetwork.parset
@@ -2,7 +2,7 @@ path = ../data/
resultPath = ../cantorfaultnetworks/results/sparse/
[preconditioner]
-patch = CELL # CELL , SUPPORT
+patch = SUPPORT # CELL , SUPPORT
mode = MULTIPLICATIVE # ADDITIVE , MULTIPLICATIVE
multDirection = SYMMETRIC # SYMMETRIC , FORWARD , BACKWARD
patchDepth = 1
@@ -13,9 +13,9 @@ patchDepth = 1
oscDataFile = oscDataLaplace4.mat
# level resolution in 2^(-...)
-coarseCantorLevel = 2
-fineCantorLevel = 3
-maxCantorLevel = 4
+coarseCantorLevel = 1
+fineCantorLevel = 2
+maxCantorLevel = 3
penaltyFactor = 1