Skip to content
Snippets Groups Projects
Commit 8f72feb2 authored by Max Kahnt's avatar Max Kahnt
Browse files

Refactor typenames, typedefs. Fix warnings.

parent 47f007f2
No related branches found
No related tags found
No related merge requests found
// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=8 sw=4 sts=4:
template<class MatrixType, class DiscFuncType>
template<class MatrixType, class VectorType>
inline
void ProjectedBlockGSStep<MatrixType, DiscFuncType>::iterate()
void ProjectedBlockGSStep<MatrixType, VectorType>::iterate()
{
if (hasObstacle_->size()!= (unsigned int)this->x_->size())
DUNE_THROW(SolverError, "Size of hasObstacle (" << hasObstacle_->size()
......@@ -18,7 +18,7 @@ void ProjectedBlockGSStep<MatrixType, DiscFuncType>::iterate()
continue;
bool zeroDiagonal = false;
for (int j=0; j<BlockSize; j++) {
for (size_t j=0; j<BlockSize; j++) {
// When using this solver as part of a truncated multigrid solver,
// the diagonal entries of the matrix may get completely truncated
// away. In this case we just do nothing here.
......@@ -45,19 +45,19 @@ void ProjectedBlockGSStep<MatrixType, DiscFuncType>::iterate()
// Solve the local constraint minimization problem
// We use a projected Gauss-Seidel, for lack of anything better
BoxConstraint<field_type,BlockSize> defectObstacle = (*obstacles_)[i];
Obstacle defectObstacle = (*obstacles_)[i];
defectObstacle -= x;
// Initial correction
v = 0;
for (int j=0; j< ((BlockSize==1) ? 1 : 20); j++) {
for (size_t j=0; j< ((BlockSize==1) ? 1 : 20); j++) {
for (int k=0; k<BlockSize; k++) {
for (size_t k=0; k<BlockSize; k++) {
// Compute residual
field_type sr = 0;
for (int l=0; l<BlockSize; l++)
Field sr = 0;
for (size_t l=0; l<BlockSize; l++)
sr += mat[i][i][k][l] * v[l];
sr = r[k] - sr;
......
......@@ -8,32 +8,30 @@
#include "blockgsstep.hh"
#include <dune/solvers/common/boxconstraint.hh>
template<class MatrixType, class DiscFuncType>
class ProjectedBlockGSStep : public BlockGSStep<MatrixType, DiscFuncType>
template<class MatrixType, class VectorType>
class ProjectedBlockGSStep : public BlockGSStep<MatrixType, VectorType>
{
typedef typename DiscFuncType::block_type VectorBlock;
typedef typename DiscFuncType::field_type field_type;
using VectorBlock = typename VectorType::block_type;
using Field = typename VectorType::field_type;
enum {BlockSize = VectorBlock::dimension};
public:
//! Default constructor. Doesn't init anything
ProjectedBlockGSStep() {}
//! Constructor with a linear problem
ProjectedBlockGSStep(const MatrixType& mat, DiscFuncType& x, const DiscFuncType& rhs)
: BlockGSStep<MatrixType,DiscFuncType>(mat, x, rhs)
ProjectedBlockGSStep(const MatrixType& mat, VectorType& x, const VectorType& rhs)
: BlockGSStep<MatrixType,VectorType>(mat, x, rhs)
{}
//! Perform one iteration
virtual void iterate();
const Dune::BitSetVector<BlockSize>* hasObstacle_;
const std::vector<BoxConstraint<field_type,BlockSize> >* obstacles_;
using HasObstacle = Dune::BitSetVector<BlockSize>;
using Obstacle = BoxConstraint<Field, BlockSize>;
const HasObstacle* hasObstacle_;
const std::vector<Obstacle>* obstacles_;
};
#include "projectedblockgsstep.cc"
......
......@@ -14,16 +14,16 @@ void TrustRegionGSStep<MatrixType, VectorType>::iterate()
{
const MatrixType& mat = *this->mat_;
VectorType& x = *this->x_;
const std::vector<BoxConstraint<field_type,blocksize> >& obstacles = *this->obstacles_;
const std::vector<BoxConstraint<Field,BlockSize> >& obstacles = *this->obstacles_;
for (size_t i=0; i<x.size(); i++) {
// Dirichlet nodes in this block
std::bitset<blocksize> ignoreNodes(0);
std::bitset<BlockSize> ignoreNodes(0);
if (this->ignoreNodes_)
ignoreNodes = (*this->ignoreNodes_)[i];
if (ignoreNodes.count() == blocksize)
if (ignoreNodes.count() == BlockSize)
continue;
VectorBlock blockResidual;
......@@ -32,16 +32,16 @@ void TrustRegionGSStep<MatrixType, VectorType>::iterate()
mat[i][i].umv(x[i], blockResidual);
// scalar gauss-seidel
for (int j=0; j<blocksize; j++) {
for (size_t j=0; j<BlockSize; j++) {
if (ignoreNodes[j])
continue;
// Compute residual
field_type r = blockResidual[j] - mat[i][i][j] * x[i];
Field r = blockResidual[j] - mat[i][i][j] * x[i];
r += mat[i][i][j][j] * x[i][j];
const field_type& diag = mat[i][i][j][j];
const Field& diag = mat[i][i][j][j];
// Find local minimum
if (diag > 0) {
......@@ -63,11 +63,11 @@ void TrustRegionGSStep<MatrixType, VectorType>::iterate()
// 1d problem is concave or linear.
// Minimum is attained at one of the boundaries
field_type lBound = obstacles[i].lower(j);
field_type uBound = obstacles[i].upper(j);
Field lBound = obstacles[i].lower(j);
Field uBound = obstacles[i].upper(j);
field_type lValue = 0.5*diag*lBound*lBound - r*lBound;
field_type uValue = 0.5*diag*uBound*uBound - r*uBound;
Field lValue = 0.5*diag*lBound*lBound - r*lBound;
Field uValue = 0.5*diag*uBound*uBound - r*uBound;
x[i][j] = (lValue < uValue) ? lBound : uBound;
......
......@@ -17,15 +17,12 @@
template<class MatrixType, class VectorType>
class TrustRegionGSStep : public ProjectedBlockGSStep<MatrixType, VectorType>
{
using VectorBlock = typename VectorType::block_type;
using Field = VectorType::field_type;
typedef typename VectorType::block_type VectorBlock;
enum {blocksize = VectorBlock::dimension};
typedef typename VectorType::field_type field_type;
enum {BlockSize = VectorBlock::dimension};
public:
//! Default constructor. Doesn't init anything
TrustRegionGSStep() {}
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment