From 013757f05e5f9d4810039254c45b7005e4bbe580 Mon Sep 17 00:00:00 2001
From: Oliver Sander <oliver.sander@tu-dresden.de>
Date: Fri, 22 Mar 2024 06:10:36 +0100
Subject: [PATCH] Allow to use MixedRiemannianPNSolver, where applicable
---
problems/cosserat-continuum-cantilever.parset | 3 +
problems/film-on-substrate.parset | 2 +-
src/cosserat-continuum.cc | 83 ++++++++++++++-----
src/film-on-substrate.cc | 14 +++-
test/filmonsubstratetest.cc | 14 +++-
5 files changed, 90 insertions(+), 26 deletions(-)
diff --git a/problems/cosserat-continuum-cantilever.parset b/problems/cosserat-continuum-cantilever.parset
index bfa21a84..ac1a0765 100644
--- a/problems/cosserat-continuum-cantilever.parset
+++ b/problems/cosserat-continuum-cantilever.parset
@@ -20,6 +20,9 @@ numLevels = 1
# Number of homotopy steps for the Dirichlet boundary conditions
numHomotopySteps = 1
+# Solver type: trustRegion or proximalNewton
+solvertype = trustRegion
+
# Tolerance of the trust region solver
tolerance = 1e-3
diff --git a/problems/film-on-substrate.parset b/problems/film-on-substrate.parset
index 1ce5db4b..8bfd2012 100644
--- a/problems/film-on-substrate.parset
+++ b/problems/film-on-substrate.parset
@@ -61,7 +61,7 @@ initialDeformation = "x"
# Solver parameters
#############################################
-# Inner solver, cholmod or multigrid
+# Solver type: trustRegion or proximalNewton
solvertype = trustRegion
# Number of homotopy steps for the Dirichlet boundary conditions
diff --git a/src/cosserat-continuum.cc b/src/cosserat-continuum.cc
index ee4f3bd2..be9c9288 100644
--- a/src/cosserat-continuum.cc
+++ b/src/cosserat-continuum.cc
@@ -43,6 +43,8 @@
#include <dune/foamgrid/foamgrid.hh>
#endif
+#include <dune/istl/multitypeblockvector.hh>
+
#include <dune/functions/functionspacebases/lagrangebasis.hh>
#include <dune/functions/functionspacebases/compositebasis.hh>
#include <dune/functions/functionspacebases/powerbasis.hh>
@@ -66,6 +68,7 @@
#include <dune/gfe/assemblers/mixedgfeassembler.hh>
#if MIXED_SPACE
+#include <dune/gfe/mixedriemannianpnsolver.hh>
#include <dune/gfe/mixedriemanniantrsolver.hh>
#else
#include <dune/gfe/assemblers/geodesicfeassemblerwrapper.hh>
@@ -504,31 +507,65 @@ int main (int argc, char *argv[]) try
RealTuple<double,3>,
Rotation<double,3> > mixedAssembler(compositeBasis, &localGFEADOLCStiffness);
#if MIXED_SPACE
- MixedRiemannianTrustRegionSolver<GridType,
- CompositeBasis,
- DeformationFEBasis, RealTuple<double,3>,
- OrientationFEBasis, Rotation<double,3> > solver;
- solver.setup(*grid,
- &mixedAssembler,
- deformationFEBasis,
- orientationFEBasis,
- x,
- deformationDirichletDofs,
- orientationDirichletDofs, tolerance,
- maxSolverSteps,
- initialTrustRegionRadius,
- multigridIterations,
- mgTolerance,
- mu, nu1, nu2,
- baseIterations,
- baseTolerance,
- instrumented);
+ if (parameterSet.get<std::string>("solvertype", "trustRegion") == "trustRegion")
+ {
- solver.setScaling(parameterSet.get<FieldVector<double,6> >("solverScaling"));
+ MixedRiemannianTrustRegionSolver<GridType,
+ CompositeBasis,
+ DeformationFEBasis, RealTuple<double,3>,
+ OrientationFEBasis, Rotation<double,3> > solver;
+ solver.setup(*grid,
+ &mixedAssembler,
+ deformationFEBasis,
+ orientationFEBasis,
+ x,
+ deformationDirichletDofs,
+ orientationDirichletDofs, tolerance,
+ maxSolverSteps,
+ initialTrustRegionRadius,
+ multigridIterations,
+ mgTolerance,
+ mu, nu1, nu2,
+ baseIterations,
+ baseTolerance,
+ instrumented);
- solver.setInitialIterate(x);
- solver.solve();
- x = solver.getSol();
+ solver.setScaling(parameterSet.get<FieldVector<double,6> >("solverScaling"));
+
+ solver.setInitialIterate(x);
+ solver.solve();
+ x = solver.getSol();
+ }
+ else
+ {
+ //Create BitVector matching the tangential space
+ const int dimRotationTangent = Rotation<double,3>::TangentVector::dimension;
+ using VectorForBit = MultiTypeBlockVector<std::vector<FieldVector<double,3> >, std::vector<FieldVector<double,dimRotationTangent> > >;
+ using BitVector = Solvers::DefaultBitVector_t<VectorForBit>;
+ BitVector dirichletDofs;
+ dirichletDofs[_0].resize(compositeBasis.size({0}));
+ dirichletDofs[_1].resize(compositeBasis.size({1}));
+ for (size_t i = 0; i < compositeBasis.size({0}); i++) {
+ for (size_t j = 0; j < 3; j++)
+ dirichletDofs[_0][i][j] = deformationDirichletDofs[i][j];
+ }
+ for (size_t i = 0; i < compositeBasis.size({1}); i++) {
+ for (int j = 0; j < dimRotationTangent; j++)
+ dirichletDofs[_1][i][j] = orientationDirichletDofs[i][j];
+ }
+ GFE::MixedRiemannianProximalNewtonSolver<CompositeBasis, DeformationFEBasis, RealTuple<double,3>, OrientationFEBasis, Rotation<double,3>, BitVector> solver;
+ solver.setup(*grid,
+ &mixedAssembler,
+ x,
+ dirichletDofs,
+ tolerance,
+ maxSolverSteps,
+ initialRegularization,
+ instrumented);
+ solver.setInitialIterate(x);
+ solver.solve();
+ x = solver.getSol();
+ }
#else
//The MixedRiemannianTrustRegionSolver can treat the Displacement and Orientation Space as separate ones
//The RiemannianTrustRegionSolver can only treat the Displacement and Rotation together in a ProductManifold.
diff --git a/src/film-on-substrate.cc b/src/film-on-substrate.cc
index 0630bf58..97fc8ef8 100644
--- a/src/film-on-substrate.cc
+++ b/src/film-on-substrate.cc
@@ -64,6 +64,7 @@
#include <dune/gfe/assemblers/sumenergy.hh>
#if MIXED_SPACE
+#include <dune/gfe/mixedriemannianpnsolver.hh>
#include <dune/gfe/mixedriemanniantrsolver.hh>
#else
#include <dune/gfe/assemblers/geodesicfeassemblerwrapper.hh>
@@ -650,7 +651,18 @@ int main (int argc, char *argv[]) try
} else { //parameterSet.get<std::string>("solvertype") == "proximalNewton"
#if MIXED_SPACE
- DUNE_THROW(Exception, "Error: There is no MixedRiemannianProximalNewtonSolver!");
+ GFE::MixedRiemannianProximalNewtonSolver<CompositeBasis, DeformationFEBasis, RealTuple<double,dim>, OrientationFEBasis, Rotation<double,dim>, BitVector> solver;
+ solver.setup(*grid,
+ &mixedAssembler,
+ x,
+ dirichletDofs,
+ tolerance,
+ maxSolverSteps,
+ initialRegularization,
+ instrumented);
+ solver.setInitialIterate(x);
+ solver.solve();
+ x = solver.getSol();
#else
RiemannianProximalNewtonSolver<DeformationFEBasis, RBM, GFEAssemblerWrapper> solver;
solver.setup(*grid,
diff --git a/test/filmonsubstratetest.cc b/test/filmonsubstratetest.cc
index 876871df..c9b40dda 100644
--- a/test/filmonsubstratetest.cc
+++ b/test/filmonsubstratetest.cc
@@ -47,6 +47,7 @@
#include <dune/gfe/assemblers/sumenergy.hh>
#if MIXED_SPACE
+#include <dune/gfe/mixedriemannianpnsolver.hh>
#include <dune/gfe/mixedriemanniantrsolver.hh>
#else
#include <dune/gfe/assemblers/geodesicfeassemblerwrapper.hh>
@@ -458,7 +459,18 @@ int main (int argc, char *argv[])
} else { // solverType == "proximalNewton"
#if MIXED_SPACE
- DUNE_THROW(Exception, "Error: There is no MixedRiemannianProximalNewtonSolver!");
+ GFE::MixedRiemannianProximalNewtonSolver<CompositeBasis, DeformationFEBasis, RealTuple<double,dim>, OrientationFEBasis, Rotation<double,dim>, BitVector> solver;
+ solver.setup(*grid,
+ &mixedAssembler,
+ x,
+ dirichletDofs,
+ tolerance,
+ maxSolverSteps,
+ 1.0 /* initialRegularization */,
+ false /* instrumented */);
+ solver.setInitialIterate(x);
+ solver.solve();
+ x = solver.getSol();
#else
RiemannianProximalNewtonSolver<DeformationFEBasis, RBM, GFEAssemblerWrapper> solver;
solver.setup(*grid,
--
GitLab