Code owners
Assign users and groups as approvers for specific file changes. Learn more.
geofaultfactory.hh 28.88 KiB
#ifndef GEO_FAULT_FACTORY_HH
#define GEO_FAULT_FACTORY_HH
#include <math.h>
#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/fufem/boundarypatch.hh>
#include <dune/fufem/boundaryiterator.hh>
class FaultCorridor {
protected:
int level_;
int lowerY_;
int upperY_;
std::vector<int> faultToLevel_;
int numberOfFaults_;
std::vector<int> separatingYIDs_;
std::vector<int> faultYs_;
private:
void setup() {
numberOfFaults_ = 0;
for (size_t i=0; i<faultToLevel_.size(); i++) {
if (faultToLevel_[i]==level_)
numberOfFaults_++;
}
if (numberOfFaults_>0) {
separatingYIDs_.resize(numberOfFaults_-1);
faultYs_.resize(numberOfFaults_);
int yStep = (upperY_-lowerY_)/numberOfFaults_;
if (yStep<2) {
std::cout << "FaultCorridor::setup() : FaultCorridor does not admit " + std::to_string(numberOfFaults_) + " faults on level " + std::to_string(level_) << std::endl;
DUNE_THROW(Dune::Exception, "FaultCorridor::setup() : FaultCorridor does not admit " + std::to_string(numberOfFaults_) + " faults on level " + std::to_string(level_));
}
int offSet = yStep + (upperY_-lowerY_-yStep*numberOfFaults_)/2;
int prevYID = lowerY_;
for (size_t i=0; i<separatingYIDs_.size(); i++) {
int newYID = lowerY_ + offSet + i*yStep;
separatingYIDs_[i] = newYID;
faultYs_[i] = (prevYID+newYID)/2;
prevYID = newYID;
}
faultYs_[faultYs_.size()-1] = (upperY_+prevYID)/2;
} else {
separatingYIDs_.resize(0);
faultYs_.resize(0);
}
}
public:
FaultCorridor() {}
FaultCorridor(const int level, const int lowerY, const int upperY, std::vector<int> faultToLevel) :
level_(level),
lowerY_(lowerY),
upperY_(upperY),
faultToLevel_(faultToLevel)
{
setup();
}
FaultCorridor(const FaultCorridor& faultCorridor) :
level_(faultCorridor.level_),
lowerY_(faultCorridor.lowerY_),
upperY_(faultCorridor.upperY_),
faultToLevel_(faultCorridor.faultToLevel_),
numberOfFaults_(faultCorridor.numberOfFaults_),
separatingYIDs_(faultCorridor.separatingYIDs_),
faultYs_(faultCorridor.faultYs_)
{}
void nextLevel(std::vector<FaultCorridor>& newFaultCorridors, const int resolutionFactor = 2) {
size_t idx = 0;
int prevYID = resolutionFactor*lowerY_;
for (size_t i=0; i<faultYs_.size(); i++) {
std::vector<int> newFaultToLevel;
for (size_t j=idx; j<faultToLevel_.size(); j++) {
idx++;
if (faultToLevel_[j]>level_)
newFaultToLevel.push_back(faultToLevel_[j]);
else
break;
}
newFaultCorridors.push_back(FaultCorridor(level_+1, prevYID, resolutionFactor*faultYs_[i], newFaultToLevel));
prevYID = resolutionFactor*faultYs_[i];
}
std::vector<int> newFaultToLevel;
for (size_t j=idx; j<faultToLevel_.size(); j++) {
idx++;
if (faultToLevel_[j]>level_)
newFaultToLevel.push_back(faultToLevel_[j]);
else
break;
}
newFaultCorridors.push_back(FaultCorridor(level_+1, prevYID, resolutionFactor*upperY_, newFaultToLevel));
}
bool isfaultY(int candidatY, int& idx) const {
for (size_t i=0; i<faultYs_.size(); i++) {
if (faultYs_[i]<candidatY)
continue;
if (faultYs_[i]==candidatY) {
idx = i;
return true;
} else
return false;
}
return false;
}
int upperY() const {
return upperY_;
}
int lowerY() const {
return lowerY_;
}
const std::vector<int>& faultToLevel() const {
return faultToLevel_;
}
const std::vector<int>& separatingYIDs() const {
return separatingYIDs_;
}
const std::vector<int>& faultYs() const {
return faultYs_;
}
};
template <class GridType>
class LevelGeoFaultFactory
{
//! 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;
typedef typename GridType::ctype ctype;
typedef typename GridType::LevelGridView GV;
typedef typename GridType::LevelIntersection Intersection;
typedef typename GridType::template Codim<0>::Entity Element;
static const int dimElement = Element::dimension;
typedef typename Dune::MultipleCodimMultipleGeomTypeMapper<GV, FaceMapperLayout > FaceMapper;
private:
const std::shared_ptr<GridType> grid_;
const int level_;
const ctype resolution_;
const bool allowLevelIntersections_;
const GV gridView_;
const typename GV::IndexSet& indexSet_;
FaceMapper faceMapper_;
std::vector<Intersection> faces_;
std::vector<Dune::FieldVector<ctype, dimworld>> vertexPositions_;
std::vector<std::vector<size_t>> vertexToFaces_;
private:
bool intersectionAllowed(const size_t faceID, const size_t vertexID, const double maxAngle, const std::set<size_t>& separatingDofs, const std::set<int>& separatingYIDs, const std::set<size_t>& faultDofs) const {
const Intersection& intersection = faces_[faceID];
//check if "back" edge, cannot deviate from desiredOrientation more than maxAngle degrees (radian)
Dune::FieldVector<ctype, dimworld> desiredOrientation(1);
desiredOrientation /= desiredOrientation.two_norm();
Dune::FieldVector<ctype, dimworld> orientation(intersection.geometry().center()); orientation -= vertexPositions_[vertexID];
orientation /= orientation.two_norm();
if (std::acos(desiredOrientation*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);
std::set<size_t>::iterator iDofIt = intersectionDofs.begin();
std::set<size_t>::iterator iDofItEnd = intersectionDofs.end();
for (; iDofIt!=iDofItEnd; ++iDofIt) {
int yID = computeYID(vertexPositions_[*iDofIt]);
if (separatingDofs.count(*iDofIt) || faultDofs.count(*iDofIt) || separatingYIDs.count(yID)) {
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 computeYID(Dune::FieldVector<ctype, dimworld> vertex) const {
return (int) ((vertex[1]-vertex[0])*resolution_);
}
void computeIntersectionDofs(const typename GV::IndexSet& idxSet, const Intersection& intersection, std::set<size_t>& intersectionDofs) const {
intersectionDofs.clear();
// loop over all vertices of the intersection
const Element& 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 = idxSet.subIndex(insideElement, idxInElement, dimElement);
intersectionDofs.insert(globalIdx);
}
}
bool isTargetBoundaryVertex(const Dune::FieldVector<ctype, dimworld>& vertex) const {
bool result = false;
for (size_t i=0; i<vertex.size(); i++) {
if (vertex[i]==1.0) {
result = true;
break;
}
}
return result;
}
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 std::set<size_t>& separatingDofs, const std::set<int>& separatingYIDs, const double maxAngle) {
// compute target at upper/right part of boundary
//std::cout << "LevelGeoFaultFactory::generateFault() " << std::endl;
//std::cout << "------------------------------------- " << std::endl << std::endl;
bool success = false;
const auto& vertexPosition = vertexPositions_[faultSeedIdx];
ctype n = vertexPosition[1]-vertexPosition[0];
Dune::FieldVector<ctype, dimworld> target;
if (n>0) {
target[0] = 1-n;
target[1] = 1;
} else {
target[0] = 1;
target[1] = 1+n;
}
std::set<size_t> faultDofs;
faultDofs.insert(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(vertexPositions_[vertexID])) {
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, maxAngle, separatingDofs, separatingYIDs, faultDofs)) {
admissibleFaces.push_back(faces[i]);
}
}
std::vector<double> distances(admissibleFaces.size(), 0);
for (size_t i=0; i<admissibleFaces.size(); i++) {
Dune::FieldVector<ctype, dimworld> vec(target);
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.erase(vertexID);
faultFaces.pop_back();
vertexToAdmissibleFaces.erase(vertexID);
} 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.insert(nextFaceDofs.begin(), nextFaceDofs.end());
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, FaultCorridors should prevent it!" << 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;
}
public:
LevelGeoFaultFactory(const std::shared_ptr<GridType> grid, const int level, const ctype resolution, const bool allowLevelIntersections = false) :
grid_(grid),
level_(level),
resolution_(resolution),
allowLevelIntersections_(allowLevelIntersections),
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));
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);
}
}
}
}
void build(std::vector<std::shared_ptr<FaultInterface<GV>>>& faults, std::vector<FaultCorridor>& faultCorridors, std::set<size_t> separatingDofs, double maxAngle = 2) {
std::vector<std::vector<size_t>> faultSeeds;
generateFaultSeeds(faultCorridors, faultSeeds);
faults.resize(0);
for (size_t i=0; i<faultCorridors.size(); i++) {
const std::vector<int>& allSeparatingYIDs = faultCorridors[i].separatingYIDs();
if (faultSeeds[i].size()==0)
continue;
if (faultSeeds[i].size()==1) {
std::shared_ptr<FaultInterface<GV>> fault = std::make_shared<FaultInterface<GV>>(gridView_, level_);
std::set<int> separatingYIDs;
generateFault(fault, faultSeeds[i][0], separatingDofs, separatingYIDs, maxAngle);
faults.push_back(fault);
} else {
// generate first fault whose limits are given by coarser fault and one separatingYID
std::shared_ptr<FaultInterface<GV>> fault = std::make_shared<FaultInterface<GV>>(gridView_, level_);
std::set<int> separatingYIDs;
separatingYIDs.insert(allSeparatingYIDs[0]);
generateFault(fault, faultSeeds[i][0], separatingDofs, separatingYIDs, maxAngle);
faults.push_back(fault);
// generate all intermediate faults whose limits are given by two separatingYIDs
for (size_t j=1; j<faultSeeds[i].size()-1; j++) {
std::shared_ptr<FaultInterface<GV>> interface = std::make_shared<FaultInterface<GV>>(gridView_, level_);
std::set<int> separatingYIDs;
separatingYIDs.insert(allSeparatingYIDs[j-1]);
separatingYIDs.insert(allSeparatingYIDs[j]);
generateFault(interface, faultSeeds[i][j], separatingDofs, separatingYIDs, maxAngle);
faults.push_back(interface);
}
// generate last fault whose limits are given by coarser fault and one separatingYID
std::shared_ptr<FaultInterface<GV>> lastFault = std::make_shared<FaultInterface<GV>>(gridView_, level_);
std::set<int> lastSeparatingYIDs;
lastSeparatingYIDs.insert(allSeparatingYIDs[allSeparatingYIDs.size()-1]);
generateFault(lastFault, faultSeeds[i][faultSeeds[i].size()-1], separatingDofs, lastSeparatingYIDs, maxAngle);
faults.push_back(lastFault);
}
}
}
};
template <class GridType>
class GeoFaultFactory
{
//! Parameter for mapper class
template<int dim>
struct FaceMapperLayout
{
bool contains (Dune::GeometryType gt)
{
return gt.dim() == dim-1;
}
};
protected:
typedef OscUnitCube<GridType, 2> GridOb;
static const int dimworld = GridType::dimensionworld ;
static const int dim = GridType::dimension;
typedef typename GridType::ctype ctype;
typedef typename GridType::LevelGridView GV;
typedef typename GridType::LevelIntersection Intersection;
typedef typename GridType::template Codim<0>::Entity Element;
static const int dimElement = Element::dimension;
typedef typename Dune::MultipleCodimMultipleGeomTypeMapper<GV, FaceMapperLayout > FaceMapper;
private:
const std::vector<int>& faultToLevel_;
const int coarseResolution_;
const size_t maxLevel_;
const bool allowIntersections_;
const bool allowLevelIntersections_;
const int coarseGridN_;
std::vector<std::set<size_t>> separatingDofs_;
std::vector<int> faultPerLevel_;
std::shared_ptr<GridType> grid_;
std::shared_ptr<InterfaceNetwork<GridType>> interfaceNetwork_;
std::vector<std::shared_ptr<LevelGeoFaultFactory<GridType>>> levelGeoFaultFactories_;
std::vector<std::vector<FaultCorridor>> faultCorridorHierarchy_;
std::vector<double> levelResolutions_;
private:
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;
}
void computeIntersectionDofs(const typename GV::IndexSet& idxSet, const Intersection& intersection, std::set<size_t>& intersectionDofs) {
intersectionDofs.clear();
// loop over all vertices of the intersection
const Element& 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 = idxSet.subIndex(insideElement, idxInElement, dimElement);
intersectionDofs.insert(globalIdx);
}
}
public:
//setup
GeoFaultFactory(const std::vector<int>& faultToLevel, const int coarseResolution, const size_t maxLevel, const bool allowIntersections = false, const bool allowLevelIntersections = false, const double maxAngle = 2) :
faultToLevel_(faultToLevel),
coarseResolution_(coarseResolution),
maxLevel_(maxLevel),
allowIntersections_(allowIntersections),
allowLevelIntersections_(allowLevelIntersections),
coarseGridN_(std::pow(2, coarseResolution_)),
interfaceNetwork_(nullptr)
{
// set faultPerLevel_:
faultPerLevel_.resize(maxLevel_+1, 0);
for (size_t i=0; i<faultToLevel_.size(); i++) {
faultPerLevel_[faultToLevel_[i]]++;
}
print(faultPerLevel_, "faultPerLevel: ");
Dune::UGGrid<dim>::setDefaultHeapSize(4000);
GridOb unitCube(coarseGridN_);
grid_ = unitCube.grid();
grid_->globalRefine(maxLevel_);
separatingDofs_.resize(maxLevel_+1);
levelResolutions_.resize(maxLevel_+1);
faultCorridorHierarchy_.resize(maxLevel_+1);
levelGeoFaultFactories_.resize(maxLevel_+1);
// init interface network
interfaceNetwork_ = std::make_shared<InterfaceNetwork<GridType>>(*grid_);
std::vector<FaultCorridor>& coarsefaultCorridors = faultCorridorHierarchy_[0];
if (faultPerLevel_[0]==0) {
FaultCorridor newFaultCorridor(0, -coarseGridN_, coarseGridN_, faultToLevel_);
coarsefaultCorridors.push_back(newFaultCorridor);
} else {
int yStep = (2*coarseGridN_-2)/(faultPerLevel_[0]);
int offSet = (2*coarseGridN_-yStep*faultPerLevel_[0])/2 - coarseGridN_;
//std::cout << "yStep: " << yStep << " offSet: " << offSet << std::endl;
FaultCorridor newFaultCorridor(0, offSet, faultPerLevel_[0]*yStep + offSet, faultToLevel_);
coarsefaultCorridors.push_back(newFaultCorridor);
}
for (size_t i=0; i<=maxLevel_; i++) {
std::vector<std::shared_ptr<FaultInterface<GV>>> faults;
std::vector<FaultCorridor>& faultCorridors = faultCorridorHierarchy_[i];
levelResolutions_[i] = std::pow(2, coarseResolution_+i);
//generate faults on level
levelGeoFaultFactories_[i] = std::make_shared<LevelGeoFaultFactory<GridType>>(grid_, i, levelResolutions_[i], allowLevelIntersections_);
levelGeoFaultFactories_[i]->build(faults, faultCorridors, separatingDofs_[i], maxAngle);
for (size_t j=0; j<faults.size(); j++) {
interfaceNetwork_->addInterface(faults[j]);
}
if (i<maxLevel_) {
const int newLevel = i+1;
//compute separating dofs for next level given by current levelinterfacenetwork
interfaceNetwork_->prolongLevel(i, newLevel);
if (!allowIntersections_) {
separatingDofs_[newLevel] = interfaceNetwork_->getInterfaceNetworkDofs(newLevel);
}
//compute faultCorridors for next level
std::vector<FaultCorridor>& newLevelFaultCorridors = faultCorridorHierarchy_[newLevel];
newLevelFaultCorridors.resize(0);
for (size_t j=0; j<faultCorridors.size(); j++) {
faultCorridors[j].nextLevel(newLevelFaultCorridors);
}
}
}
}
/*
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