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