From e5201c44e5914bef924ca4a0d75b9f8b954e321b Mon Sep 17 00:00:00 2001
From: Oliver Sander <sander@igpm.rwth-aachen.de>
Date: Fri, 17 Apr 2009 20:00:18 +0000
Subject: [PATCH] RodSolver is now called RiemannianTrustRegionSolver, and can
be parametrized with the TargetSpace. Some remains of RigidBodyMotion are
still there, though.
[[Imported from SVN: r4041]]
---
rod3d.cc | 4 +-
src/riemanniantrsolver.cc | 84 ++++++++++++++++++---------------------
src/riemanniantrsolver.hh | 23 ++++++-----
3 files changed, 54 insertions(+), 57 deletions(-)
diff --git a/rod3d.cc b/rod3d.cc
index 6251722c..f7dfc756 100644
--- a/rod3d.cc
+++ b/rod3d.cc
@@ -16,7 +16,7 @@
#include "src/rodwriter.hh"
#include "src/rotation.hh"
#include "src/rodassembler.hh"
-#include "src/rodsolver.hh"
+#include "src/riemanniantrsolver.hh"
typedef RigidBodyMotion<3> TargetSpace;
@@ -120,7 +120,7 @@ int main (int argc, char *argv[]) try
RodAssembler<GridType> rodAssembler(grid);
rodAssembler.setShapeAndMaterial(A, J1, J2, E, nu);
- RodSolver<GridType> rodSolver;
+ RiemannianTrustRegionSolver<GridType,RigidBodyMotion<3> > rodSolver;
rodSolver.setup(grid,
&rodAssembler,
x,
diff --git a/src/riemanniantrsolver.cc b/src/riemanniantrsolver.cc
index 33705f26..31b3fefd 100644
--- a/src/riemanniantrsolver.cc
+++ b/src/riemanniantrsolver.cc
@@ -15,39 +15,33 @@
#include <dune-solvers/norms/energynorm.hh>
#include <dune-solvers/norms/h1seminorm.hh>
-#include "rigidbodymotion.hh"
-#include "quaternion.hh"
#include "maxnormtrustregion.hh"
// for debugging
#include "../test/fdcheck.hh"
-// Number of degrees of freedom:
-// 7 (x, y, z, q_1, q_2, q_3, q_4) for a spatial rod
-
-template <class GridType>
-void RodSolver<GridType>::setup(const GridType& grid,
- const GeodesicFEAssembler<typename GridType::LeafGridView, RigidBodyMotion<3> >* assembler,
- const SolutionType& x,
- const Dune::BitSetVector<blocksize>& dirichletNodes,
- double tolerance,
- int maxTrustRegionSteps,
- double initialTrustRegionRadius,
- int multigridIterations,
- double mgTolerance,
- int mu,
- int nu1,
- int nu2,
- int baseIterations,
- double baseTolerance,
- bool instrumented)
+template <class GridType, class TargetSpace>
+void RiemannianTrustRegionSolver<GridType,TargetSpace>::
+setup(const GridType& grid,
+ const GeodesicFEAssembler<typename GridType::LeafGridView,TargetSpace>* assembler,
+ const SolutionType& x,
+ const Dune::BitSetVector<blocksize>& dirichletNodes,
+ double tolerance,
+ int maxTrustRegionSteps,
+ double initialTrustRegionRadius,
+ int multigridIterations,
+ double mgTolerance,
+ int mu,
+ int nu1,
+ int nu2,
+ int baseIterations,
+ double baseTolerance,
+ bool instrumented)
{
- using namespace Dune;
-
- grid_ = &grid;
- assembler_ = assembler;
+ grid_ = &grid;
+ assembler_ = assembler;
x_ = x;
- tolerance_ = tolerance;
+ this->tolerance_ = tolerance;
maxTrustRegionSteps_ = maxTrustRegionSteps;
initialTrustRegionRadius_ = initialTrustRegionRadius;
multigridIterations_ = multigridIterations;
@@ -94,12 +88,12 @@ void RodSolver<GridType>::setup(const GridType& grid,
// //////////////////////////////////////////////////////////////////////////////////////
// Assemble a Laplace matrix to create a norm that's equivalent to the H1-norm
// //////////////////////////////////////////////////////////////////////////////////////
- LeafP1Function<GridType,double> u(grid),f(grid);
- LaplaceLocalStiffness<typename GridType::LeafGridView,double> laplaceStiffness;
- LeafP1OperatorAssembler<GridType,double,1>* A = new LeafP1OperatorAssembler<GridType,double,1>(grid);
+ Dune::LeafP1Function<GridType,double> u(grid),f(grid);
+ Dune::LaplaceLocalStiffness<typename GridType::LeafGridView,double> laplaceStiffness;
+ Dune::LeafP1OperatorAssembler<GridType,double,1>* A = new Dune::LeafP1OperatorAssembler<GridType,double,1>(grid);
A->assemble(laplaceStiffness,u,f);
- typedef typename LeafP1OperatorAssembler<GridType,double,1>::RepresentationType LaplaceMatrixType;
+ typedef typename Dune::LeafP1OperatorAssembler<GridType,double,1>::RepresentationType LaplaceMatrixType;
if (h1SemiNorm_)
delete h1SemiNorm_;
@@ -124,7 +118,7 @@ void RodSolver<GridType>::setup(const GridType& grid,
delete hessianMatrix_;
hessianMatrix_ = new MatrixType;
- MatrixIndexSet indices(grid_->size(1), grid_->size(1));
+ Dune::MatrixIndexSet indices(grid_->size(1), grid_->size(1));
assembler_->getNeighborsPerVertex(indices);
indices.exportIdx(*hessianMatrix_);
@@ -134,7 +128,7 @@ void RodSolver<GridType>::setup(const GridType& grid,
hasObstacle_.resize(numLevels);
for (int i=0; i<hasObstacle_.size(); i++)
- hasObstacle_[i].resize(grid_->size(i, 1),true);
+ hasObstacle_[i].resize(grid_->size(i, gridDim),true);
// ////////////////////////////////////
// Create the transfer operators
@@ -152,8 +146,8 @@ void RodSolver<GridType>::setup(const GridType& grid,
}
-template <class GridType>
-void RodSolver<GridType>::solve()
+template <class GridType, class TargetSpace>
+void RiemannianTrustRegionSolver<GridType,TargetSpace>::solve()
{
MaxNormTrustRegion<blocksize> trustRegion(x_.size(), initialTrustRegionRadius_);
@@ -176,7 +170,7 @@ void RodSolver<GridType>::solve()
// /////////////////////////////////////////////////////
for (int i=0; i<maxTrustRegionSteps_; i++) {
- if (this->verbosity_ == FULL) {
+ if (this->verbosity_ == Solver::FULL) {
std::cout << "----------------------------------------------------" << std::endl;
std::cout << " Trust-Region Step Number: " << i
<< ", radius: " << trustRegion.radius()
@@ -284,7 +278,7 @@ void RodSolver<GridType>::solve()
}
- if (this->verbosity_ == FULL) {
+ if (this->verbosity_ == NumProc::FULL) {
double translationMax = 0;
double rotationMax = 0;
for (size_t j=0; j<corr.size(); j++) {
@@ -296,11 +290,11 @@ void RodSolver<GridType>::solve()
printf("infinity norm of the correction: %g %g\n", translationMax, rotationMax);
}
- if (corr.infinity_norm() < tolerance_) {
- if (this->verbosity_ == FULL)
+ if (corr.infinity_norm() < this->tolerance_) {
+ if (this->verbosity_ == NumProc::FULL)
std::cout << "CORRECTION IS SMALL ENOUGH" << std::endl;
- if (this->verbosity_ != QUIET)
+ if (this->verbosity_ != NumProc::QUIET)
std::cout << i+1 << " trust-region steps were taken." << std::endl;
break;
}
@@ -334,7 +328,7 @@ void RodSolver<GridType>::solve()
hessianMatrix_->umv(corr, tmp);
double modelDecrease = (rhs*corr) - 0.5 * (corr*tmp);
- if (this->verbosity_ == FULL) {
+ if (this->verbosity_ == NumProc::FULL) {
std::cout << "Absolute model decrease: " << modelDecrease
<< ", functional decrease: " << oldEnergy - energy << std::endl;
std::cout << "Relative model decrease: " << modelDecrease / energy
@@ -344,16 +338,16 @@ void RodSolver<GridType>::solve()
assert(modelDecrease >= 0);
if (energy >= oldEnergy) {
- if (this->verbosity_ == FULL)
+ if (this->verbosity_ == NumProc::FULL)
printf("Richtung ist keine Abstiegsrichtung!\n");
}
if (energy >= oldEnergy &&
(std::abs(oldEnergy-energy)/energy < 1e-9 || modelDecrease/energy < 1e-9)) {
- if (this->verbosity_ == FULL)
+ if (this->verbosity_ == NumProc::FULL)
std::cout << "Suspecting rounding problems" << std::endl;
- if (this->verbosity_ != QUIET)
+ if (this->verbosity_ != NumProc::QUIET)
std::cout << i+1 << " trust-region steps were taken." << std::endl;
x_ = newIterate;
@@ -377,12 +371,12 @@ void RodSolver<GridType>::solve()
} else {
// unsuccessful iteration
trustRegion.scale(0.5);
- if (this->verbosity_ == FULL)
+ if (this->verbosity_ == NumProc::FULL)
std::cout << "Unsuccessful iteration!" << std::endl;
}
// Write current energy
- if (this->verbosity_ == FULL)
+ if (this->verbosity_ == NumProc::FULL)
std::cout << "--- Current energy: " << energy << " ---" << std::endl;
// /////////////////////////////////////////////////////////////////////
diff --git a/src/riemanniantrsolver.hh b/src/riemanniantrsolver.hh
index 6c00fa49..2d27153c 100644
--- a/src/riemanniantrsolver.hh
+++ b/src/riemanniantrsolver.hh
@@ -16,11 +16,14 @@
#include "rigidbodymotion.hh"
-/** \brief Riemannian trust-region solver for 3d Cosserat rod problems */
-template <class GridType>
-class RodSolver : public IterativeSolver<std::vector<RigidBodyMotion<3> >, Dune::BitSetVector<6> >
+/** \brief Riemannian trust-region solver for geodesic finite-element problems */
+template <class GridType, class TargetSpace>
+class RiemannianTrustRegionSolver
+ : public IterativeSolver<std::vector<TargetSpace>, Dune::BitSetVector<6> >
{
- const static int blocksize = 6;
+ const static int blocksize = TargetSpace::TangentVector::size;
+
+ const static int gridDim = GridType::dimension;
// Centralize the field type here
typedef double field_type;
@@ -28,17 +31,17 @@ class RodSolver : public IterativeSolver<std::vector<RigidBodyMotion<3> >, Dune:
// Some types that I need
typedef Dune::BCRSMatrix<Dune::FieldMatrix<field_type, blocksize, blocksize> > MatrixType;
typedef Dune::BlockVector<Dune::FieldVector<field_type, blocksize> > CorrectionType;
- typedef std::vector<RigidBodyMotion<3> > SolutionType;
+ typedef std::vector<TargetSpace> SolutionType;
public:
- RodSolver()
- : IterativeSolver<std::vector<RigidBodyMotion<3> >, Dune::BitSetVector<6> >(0,100,NumProc::FULL),
+ RiemannianTrustRegionSolver()
+ : IterativeSolver<std::vector<TargetSpace>, Dune::BitSetVector<blocksize> >(0,100,NumProc::FULL),
hessianMatrix_(NULL), h1SemiNorm_(NULL)
{}
void setup(const GridType& grid,
- const GeodesicFEAssembler<typename GridType::LeafGridView, RigidBodyMotion<3> >* rodAssembler,
+ const GeodesicFEAssembler<typename GridType::LeafGridView, TargetSpace>* rodAssembler,
const SolutionType& x,
const Dune::BitSetVector<blocksize>& dirichletNodes,
double tolerance,
@@ -99,7 +102,7 @@ protected:
MatrixType* hessianMatrix_;
/** \brief The assembler for the material law */
- const GeodesicFEAssembler<typename GridType::LeafGridView, RigidBodyMotion<3> >* assembler_;
+ const GeodesicFEAssembler<typename GridType::LeafGridView, TargetSpace>* assembler_;
/** \brief The multigrid solver */
LoopSolver<CorrectionType>* mmgSolver_;
@@ -116,6 +119,6 @@ protected:
};
-#include "rodsolver.cc"
+#include "riemanniantrsolver.cc"
#endif
--
GitLab