Code owners
Assign users and groups as approvers for specific file changes. Learn more.
supportpatchfactory.hh 9.71 KiB
#ifndef SUPPORT_PATCH_FACTORY_HH
#define SUPPORT_PATCH_FACTORY_HH
#include <queue>
#include <set>
#include <dune/common/bitsetvector.hh>
#include <dune/common/fvector.hh>
#include <dune/fufem/referenceelementhelper.hh>
#include <dune/grid/common/mcmgmapper.hh>
#include <dune/faultnetworks/hierarchicleveliterator.hh>
#include <dune/faultnetworks/utils/debugutils.hh>
template <class BasisType>
class SupportPatchFactory
{
protected:
typedef typename BasisType::GridView GV;
typedef typename BasisType::LocalFiniteElement LFE;
typedef typename GV::Grid GridType;
typedef typename GridType::ctype ctype;
static const int dim = GridType::dimension;
typedef typename GridType::template Codim<0>::Entity Element;
typedef typename GridType::template Codim<dim>::Entity Vertex;
typedef HierarchicLevelIterator<GridType> HierarchicLevelIteratorType;
typedef typename Dune::LevelMultipleCodimMultipleGeomTypeMapper<GridType, Dune::MCMGElementLayout > ElementMapper;
private:
const BasisType& coarseBasis_;
const BasisType& fineBasis_;
const int coarseLevel_;
const int fineLevel_;
const std::vector< std::vector <Element>>& vertexInElements_;
const int centerVertexIdx_;
const int patchDepth_;
const GridType& grid;
const GV& coarseGridView;
const GV& fineGridView;
std::vector<int> localToGlobal;
Dune::BitSetVector<1> boundaryDofs;
std::vector<Element> fineRegionElements;
ElementMapper coarseMapper;
public:
//setup
SupportPatchFactory(const BasisType& coarseBasis, const BasisType& fineBasis, const int coarseLevel, const int fineLevel, const std::vector< std::vector <Element>>& vertexInElements, const int centerVertexIdx, const int patchDepth) :
coarseBasis_(coarseBasis),
fineBasis_(fineBasis),
coarseLevel_(coarseLevel),
fineLevel_(fineLevel),
vertexInElements_(vertexInElements),
centerVertexIdx_(centerVertexIdx),
patchDepth_(patchDepth),
grid(coarseBasis_.getGridView().grid()),
coarseGridView(coarseBasis_.getGridView()),
fineGridView(fineBasis_.getGridView()),
coarseMapper(grid, coarseLevel_)
{
fineRegionElements.resize(0);
std::vector<Element> coarseElements;
Dune::BitSetVector<1> visited(coarseMapper.size());
Dune::BitSetVector<1> vertexVisited(vertexInElements_.size());
visited.unsetAll();
vertexVisited.unsetAll();
std::set<int> localDofs;
std::set<int> localBoundaryDofs;
//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;
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);
if (!vertexVisited[dofIndex][0]) {
addVertexSupport(dofIndex, coarseElementNeighbors, visited, vertexVisited);
}
}
}
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, coarseBoundaryElements);
}
localToGlobal.resize(localDofs.size());
boundaryDofs.resize(localDofs.size());
boundaryDofs.unsetAll();
std::set<int>::iterator dofIt = localDofs.begin();
std::set<int>::iterator dofEndIt = localDofs.end();
size_t i=0;
for (; dofIt != dofEndIt; ++dofIt) {
int localDof = *dofIt;
localToGlobal[i] = localDof;
if (localBoundaryDofs.count(localDof)) {
boundaryDofs[i][0] = true;
}
i++;
}
}
std::vector<int>& getLocalToGlobal() {
return localToGlobal;
}
std::vector<Element>& getRegionElements() {
return fineRegionElements;
}
Dune::BitSetVector<1>& getBoundaryDofs() {
return boundaryDofs;
}
private:
/* 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];
for (size_t i=0; i<vertexInElement.size(); i++) {
const Element& elem = vertexInElement[i];
const int elemIdx = coarseMapper.index(elem);
if (!visited[elemIdx][0]) {
coarseElements.push_back(elem);
visited[elemIdx][0] = true;
}
}
vertexVisited[vertexIdx][0] = true;
}
bool containsInsideSubentity(const Element& elem, const typename GridType::LevelIntersection& intersection, int subEntity, int codim) {
return ReferenceElementHelper<double, dim>::subEntityContainsSubEntity(elem.type(), intersection.indexInInside(), 1, subEntity, codim);
}
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, coarseBoundaryElements);
}
}
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);
/*
const auto& fineGeometry = fineElem.geometry();
const auto& coarseGeometry = coarseElem.geometry();
const auto& ref = Dune::ReferenceElements<ctype, dim>::general(fineElem.type()); */
// insert local dofs
for (size_t i=0; i<fineLFE.localBasis().size(); ++i) {
int dofIndex = fineBasis_.index(fineElem, i);
localDofs.insert(dofIndex);
}
// search for parts of the global and inner boundary
typename GridType::LevelIntersectionIterator neighborIt = fineGridView.ibegin(fineElem);
typename GridType::LevelIntersectionIterator neighborItEnd = fineGridView.iend(fineElem);
for (; neighborIt!=neighborItEnd; ++neighborIt) {
const typename GridType::LevelIntersection& intersection = *neighborIt;
bool isInnerBoundary = false;
if (intersection.neighbor()) {
Element outsideFather = intersection.outside();
while (outsideFather.level()>coarseLevel_) {
outsideFather = outsideFather.father();
}
auto outsideFatherIdx = coarseMapper.index(outsideFather);
if (!inCoarsePatch[outsideFatherIdx][0] and !coarseBoundaryElements.count(outsideFatherIdx)) {
isInnerBoundary = true;
}
}
if (intersection.boundary() or isInnerBoundary) {
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();
if (containsInsideSubentity(fineElem, intersection, entity, codim))
localBoundaryDofs.insert(fineBasis_.index(intersection.inside(), i));
}
}
}
}
};
#endif