Skip to content
Snippets Groups Projects
Select Git revision
  • d283d4ee96ef01069022ed6c1d8cccf441a64065
  • master default protected
  • 2021-KornhuberPodlesnyYserentant
  • Dissertation2021
  • transfer-operator-based-matrix-hierarchy
  • KornhuberPodlesny2020
  • 2018-HeidaKornhuberPodlesny
7 results

localproblem.hh

  • user avatar
    .
    podlesny authored
    d283d4ee
    History
    Code owners
    Assign users and groups as approvers for specific file changes. Learn more.
    localproblem.hh 5.49 KiB
    #ifndef OSC_LOCAL_PROBLEM_HH
    #define OSC_LOCAL_PROBLEM_HH
    
    #include <math.h>   
    #include <dune/common/fmatrix.hh>
    #include <dune/common/function.hh>
    #include <dune/common/timer.hh>
    
    #include <dune/istl/matrixindexset.hh>
    #include <dune/matrix-vector/axpy.hh>
    
    #include <dune/istl/umfpack.hh>
    
    #include <dune/fufem/assemblers/localoperatorassembler.hh>
    
    #include <dune/faultnetworks/faultinterface.hh>
    
    template <class MatrixType, class DomainType, class RangeType = DomainType>
    class OscLocalProblem {
      
    private:
        typedef typename MatrixType::block_type block_type;
        typedef typename MatrixType::field_type ctype;
        const MatrixType& mat_;
        const RangeType& rhs_;
    
        const std::vector<int> dofs_; // local to global map
        const Dune::BitSetVector<1> boundaryDofs_; // local boundary dofs
    
        MatrixType localMat;
        RangeType localRhs;
    
    public:
        OscLocalProblem(const MatrixType& mat,
                        const RangeType& rhs,
                        const std::vector<int> dofs,
                        const Dune::BitSetVector<1> boundaryDofs) :
          mat_(mat),
    	  rhs_(rhs),
    	  dofs_(dofs),
          boundaryDofs_(boundaryDofs)
    	{
    
    	if (dofs_.size()!=boundaryDofs_.size())
    	    DUNE_THROW(Dune::Exception, "Sizes of dofs and boundaryDofs do not match!");
    	  
    	// construct globalToLocal map
    	std::map<int, int> globalToLocal;
    	for (size_t i=0; i<dofs_.size(); ++i) {
    	    globalToLocal[dofs_[i]] = i;
    	}
    
    	// build local stiffness matrix
            size_t localDim = dofs_.size();
            Dune::MatrixIndexSet localIdxSet(localDim, localDim);
    	
    	typedef typename MatrixType::row_type RowType;
    	typedef typename RowType::ConstIterator ColumnIterator;
    
        for(size_t rowIdx = 0; rowIdx<localDim; rowIdx++) {
    	    const int globalRowIdx = dofs_[rowIdx]; 
    	    const RowType& row = mat_[globalRowIdx];
    	    
    	    ColumnIterator colIt = row.begin();
    	    const ColumnIterator colEndIt = row.end();
                for(; colIt!=colEndIt; ++colIt) {
    		const int globalColIdx = colIt.index(); 
    		if (*colIt!=0 and globalToLocal.count(globalColIdx)!=0) {
    		    localIdxSet.add(rowIdx, globalToLocal[globalColIdx]);
    		}
                }
            } 
    
            localIdxSet.exportIdx(localMat);
    
            for(size_t rowIdx = 0; rowIdx<localMat.N(); rowIdx++) {
                RowType& row = localMat[rowIdx];
    	    const RowType& globalRow = mat_[dofs_[rowIdx]];
    	    
                ColumnIterator cIt    = row.begin();
                ColumnIterator cEndIt = row.end();
    
                for(; cIt!=cEndIt; ++cIt) {	
    		row[cIt.index()] = globalRow[dofs_[cIt.index()]];
                }
    	} 	
    
    	// init local rhs
        localRhs.resize(localDim);
            localRhs = 0;
    
    	// set local rhs and boundary dofs
    	for (size_t i=0; i<boundaryDofs_.size(); i++) {
    
                if(!boundaryDofs_[i][0]) {
    		localRhs[i] = rhs_[dofs_[i]];
    		continue;
    	    }
    			  
                RowType& row = localMat[i];
    
                ColumnIterator cIt    = row.begin();
                ColumnIterator cEndIt = row.end();
    
                for(; cIt!=cEndIt; ++cIt) {
                    row[cIt.index()] = 0;
                }
                
                row[i] = 1;
    	}
        }
    
        MatrixType& getLocalMat() {
    	return localMat;
        }
        
       /* void setRhs(const RangeType& d){
            localRhs.resize(d.size());
    	
            for (size_t i=0; i<localRhs.size(); ++i) {
                localRhs[i] = d[i];
            }
        }*/
        
        void updateBoundary(const RangeType& boundary){
            for (size_t i=0; i<localRhs.size(); i++) {
                if (boundaryDofs_[i][0])
                    localRhs[i] = boundary[dofs_[i]];
            }
        }
    
        void updateVector(const RangeType& update, RangeType& vec){
            assert(update.size()==dofs_.size());
    
            for (size_t i=0; i<localRhs.size(); i++) {
                if (!boundaryDofs_[i][0])
                    vec[dofs_[i]] = update[i];
            }
        }
    
        void updateRhs(const RangeType& rhs){
            for (size_t i=0; i<localRhs.size(); i++) {
                if (!boundaryDofs_[i][0])
                    localRhs[i] = rhs[dofs_[i]];
            }
        }
    
        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) {
            for (size_t i=0; i<localRhs.size(); i++) {
                if (boundaryDofs_[i][0]) {
                    localRhs[i] = boundary[dofs_[i]];
                } else {
                    localRhs[i] = rhs[dofs_[i]];
                }
            }
        }
    
        void solve(DomainType& x){
            //print(boundaryDofs_, "boundaryDofs:");
    
            #if HAVE_UMFPACK
                RangeType localRhsCopy(localRhs);
                Dune::InverseOperatorResult res;
    
                x.resize(localMat.M());
    
                Dune::UMFPack<MatrixType> directSolver(localMat);
                directSolver.apply(x, localRhsCopy, res);
            #else
            #error No UMFPack!
            #endif
        }
        
        void prolong(DomainType& x, DomainType& res){
            res.resize(mat_.M());
            res = 0;
    
            for (size_t i=0; i<x.size(); i++) {
                    res[dofs_[i]] = x[i];
            }
        }
    
        void restrict(const RangeType& x, RangeType& res){
            res.resize(dofs_.size());
            res = 0;
    
            for (size_t i=0; i<res.size(); i++) {
                if (boundaryDofs_[i][0])
                    res[i] = localRhs[i];
                else
                    res[i] = x[dofs_[i]];
    
            }
        }
    
        void localResidual(const DomainType& x, RangeType& res){
            res = localRhs;
            Dune::MatrixVector::subtractProduct(res, localMat, x);
        }
    
    };
    
    #endif