From 81a90c95746af98d0a77661a1f5f4dce0468157b Mon Sep 17 00:00:00 2001
From: Oliver Sander <sander@igpm.rwth-aachen.de>
Date: Mon, 6 Apr 2009 14:24:14 +0000
Subject: [PATCH] moved to new module dune-solvers
[[Imported from SVN: r2345]]
---
.../iterationsteps/projectedblockgsstep.cc | 87 +++++++++++++++++++
.../iterationsteps/projectedblockgsstep.hh | 41 +++++++++
2 files changed, 128 insertions(+)
create mode 100644 dune-solvers/iterationsteps/projectedblockgsstep.cc
create mode 100644 dune-solvers/iterationsteps/projectedblockgsstep.hh
diff --git a/dune-solvers/iterationsteps/projectedblockgsstep.cc b/dune-solvers/iterationsteps/projectedblockgsstep.cc
new file mode 100644
index 00000000..30902ed0
--- /dev/null
+++ b/dune-solvers/iterationsteps/projectedblockgsstep.cc
@@ -0,0 +1,87 @@
+template<class OperatorType, class DiscFuncType>
+inline
+DiscFuncType ProjectedBlockGSStep<OperatorType, DiscFuncType>::getSol()
+{
+ return *(this->x_);
+}
+
+template<class OperatorType, class DiscFuncType>
+inline
+void ProjectedBlockGSStep<OperatorType, DiscFuncType>::iterate()
+{
+ if (hasObstacle_->size()!= (unsigned int)this->x_->size())
+ DUNE_THROW(SolverError, "Size of hasObstacle (" << hasObstacle_->size()
+ << ") doesn't match solution vector (" << this->x_->size() << ")");
+
+ const OperatorType& mat = *this->mat_;
+
+ for (int i=0; i<this->x_->size(); i++) {
+
+ /** \todo Handle more general boundary conditions */
+ if ((*this->ignoreNodes_)[i][0])
+ continue;
+
+ bool zeroDiagonal = false;
+ for (int 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.
+ zeroDiagonal |= (mat[i][i][j][j] < 1e-10);
+ }
+
+ if (zeroDiagonal)
+ continue;
+
+ VectorBlock r;
+ residual(i, r);
+
+ // Compute x_i += A_{i,i}^{-1} r[i]
+ VectorBlock v;
+ VectorBlock& x = (*this->x_)[i];
+
+
+ if ((*hasObstacle_)[i] == false) {
+
+ // No obstacle --> solve linear problem
+ mat[i][i].solve(v, r);
+
+ } else {
+
+ // Solve the local constraint minimization problem
+ // We use a projected Gauss-Seidel, for lack of anything better
+ BoxConstraint<field_type,BlockSize> defectObstacle = (*obstacles_)[i];
+ defectObstacle -= x;
+
+ // Initial correction
+ v = 0;
+
+ for (int j=0; j<20; j++) {
+
+ for (int k=0; k<BlockSize; k++) {
+
+ // Compute residual
+ field_type sr = 0;
+ for (int l=0; l<BlockSize; l++)
+ sr += mat[i][i][k][l] * v[l];
+
+ sr = r[k] - sr;
+ v[k] += sr / mat[i][i][k][k];
+
+ // Project
+ if (v[k] < defectObstacle.lower(k))
+ v[k] = defectObstacle.lower(k);
+ else if (v[k] > defectObstacle.upper(k))
+ v[k] = defectObstacle.upper(k);
+
+ }
+
+ }
+
+ }
+
+ // Add correction
+ x += v;
+
+ }
+
+}
diff --git a/dune-solvers/iterationsteps/projectedblockgsstep.hh b/dune-solvers/iterationsteps/projectedblockgsstep.hh
new file mode 100644
index 00000000..a81e136e
--- /dev/null
+++ b/dune-solvers/iterationsteps/projectedblockgsstep.hh
@@ -0,0 +1,41 @@
+#ifndef DUNE_PROJECTED_BLOCK_GAUSS_SEIDEL_STEP_HH
+#define DUNE_PROJECTED_BLOCK_GAUSS_SEIDEL_STEP_HH
+
+#include <vector>
+
+#include "blockgsstep.hh"
+#include "boxconstraint.hh"
+
+template<class OperatorType, class DiscFuncType>
+class ProjectedBlockGSStep : public BlockGSStep<OperatorType, DiscFuncType>
+{
+
+ typedef typename DiscFuncType::block_type VectorBlock;
+
+ typedef typename DiscFuncType::field_type field_type;
+
+ enum {BlockSize = VectorBlock::dimension};
+
+public:
+
+ //! Default constructor. Doesn't init anything
+ ProjectedBlockGSStep() {}
+
+ //! Constructor with a linear problem
+ ProjectedBlockGSStep(OperatorType& mat, DiscFuncType& x, DiscFuncType& rhs)
+ : LinearIterationStep<OperatorType,DiscFuncType>(mat, x, rhs)
+ {}
+
+ //! Perform one iteration
+ virtual void iterate();
+
+ virtual DiscFuncType getSol();
+
+ Dune::BitSetVector<1>* hasObstacle_;
+
+ const std::vector<BoxConstraint<field_type,BlockSize> >* obstacles_;
+};
+
+#include "projectedblockgsstep.cc"
+
+#endif
--
GitLab