From d12bfe57ac0579017fb843d4c21d567a44ec6574 Mon Sep 17 00:00:00 2001
From: Oliver Sander <oliver.sander@tu-dresden.de>
Date: Thu, 3 Dec 2015 08:03:53 +0100
Subject: [PATCH] Implement trust-region scaling for the mixed-order Cosserat
 implementation

---
 dune/gfe/mixedriemanniantrsolver.cc |  2 ++
 dune/gfe/mixedriemanniantrsolver.hh | 18 +++++++++++++++++-
 src/mixed-cosserat-continuum.cc     |  2 ++
 3 files changed, 21 insertions(+), 1 deletion(-)

diff --git a/dune/gfe/mixedriemanniantrsolver.cc b/dune/gfe/mixedriemanniantrsolver.cc
index 6b26eb0a..345919c4 100644
--- a/dune/gfe/mixedriemanniantrsolver.cc
+++ b/dune/gfe/mixedriemanniantrsolver.cc
@@ -311,6 +311,8 @@ void MixedRiemannianTrustRegionSolver<GridType,Basis0,TargetSpace0,Basis1,Target
     // \todo Use global index set instead of basis for parallel computations
     MaxNormTrustRegion<blocksize0> trustRegion0(assembler_->basis0_.indexSet().size(), initialTrustRegionRadius_);
     MaxNormTrustRegion<blocksize1> trustRegion1(assembler_->basis1_.indexSet().size(), initialTrustRegionRadius_);
+    trustRegion0.set(initialTrustRegionRadius_, std::get<0>(scaling_));
+    trustRegion1.set(initialTrustRegionRadius_, std::get<1>(scaling_));
 
     std::vector<BoxConstraint<field_type,blocksize0> > trustRegionObstacles0;
     std::vector<BoxConstraint<field_type,blocksize1> > trustRegionObstacles1;
diff --git a/dune/gfe/mixedriemanniantrsolver.hh b/dune/gfe/mixedriemanniantrsolver.hh
index 69306d4a..4c311e0c 100644
--- a/dune/gfe/mixedriemanniantrsolver.hh
+++ b/dune/gfe/mixedriemanniantrsolver.hh
@@ -57,7 +57,10 @@ public:
         : NumProc(NumProc::FULL),
 
           h1SemiNorm0_(nullptr), h1SemiNorm1_(nullptr)
-    {}
+    {
+        std::fill(std::get<0>(scaling_).begin(), std::get<0>(scaling_).end(), 1.0);
+        std::fill(std::get<1>(scaling_).begin(), std::get<1>(scaling_).end(), 1.0);
+    }
 
     /** \brief Set up the solver using a monotone multigrid method as the inner solver */
     void setup(const GridType& grid,
@@ -77,6 +80,16 @@ public:
                int baseIterations,
                double baseTolerance,
                bool instrumented);
+
+    void setScaling(const Dune::FieldVector<double,blocksize0+blocksize1>& scaling)
+    {
+      for (int i=0; i<3; i++)
+      {
+        std::get<0>(scaling_)[i] = scaling[i];
+        std::get<1>(scaling_)[i] = scaling[i+3];
+      }
+    }
+
 #if 0
     void setIgnoreNodes(const Dune::BitSetVector<blocksize0>& ignoreNodes)
     {
@@ -114,6 +127,9 @@ protected:
     /** \brief The initial trust-region radius in the maximum-norm */
     double initialTrustRegionRadius_;
 
+    /** \brief Trust-region norm scaling */
+    std::tuple<Dune::FieldVector<double,3>, Dune::FieldVector<double,3> > scaling_;
+
     /** \brief Maximum number of trust-region steps */
     int maxTrustRegionSteps_;
 
diff --git a/src/mixed-cosserat-continuum.cc b/src/mixed-cosserat-continuum.cc
index 9a3512eb..dd7230da 100644
--- a/src/mixed-cosserat-continuum.cc
+++ b/src/mixed-cosserat-continuum.cc
@@ -320,6 +320,8 @@ int main (int argc, char *argv[]) try
                  baseTolerance,
                  instrumented);
 
+        solver.setScaling(parameterSet.get<FieldVector<double,6> >("trustRegionScaling"));
+
         ////////////////////////////////////////////////////////
         //   Set Dirichlet values
         ////////////////////////////////////////////////////////
-- 
GitLab