Skip to content
Snippets Groups Projects
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