localproblem.hh

#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