From 626e141fa0e67a94978e76331bfc69ecea59971e Mon Sep 17 00:00:00 2001
From: Oliver Sander <sander@igpm.rwth-aachen.de>
Date: Wed, 26 Jan 2011 13:50:58 +0000
Subject: [PATCH] add a scaled identity to the rod stiffness matrix if there is
no Dirichlet boundary: this regularizes the problem
[[Imported from SVN: r6886]]
---
.../gfe/coupling/rodcontinuumsteklovpoincarestep.hh | 13 ++++++++++++-
1 file changed, 12 insertions(+), 1 deletion(-)
diff --git a/dune/gfe/coupling/rodcontinuumsteklovpoincarestep.hh b/dune/gfe/coupling/rodcontinuumsteklovpoincarestep.hh
index b8a70a4d..7522f41a 100644
--- a/dune/gfe/coupling/rodcontinuumsteklovpoincarestep.hh
+++ b/dune/gfe/coupling/rodcontinuumsteklovpoincarestep.hh
@@ -212,6 +212,10 @@ private:
/** \brief Damping factor for the Richardson iteration */
double richardsonDamping_;
+public:
+ double neumannRegularization_;
+private:
+
//////////////////////////////////////////////////////////////////
// Data members related to the rod problems
//////////////////////////////////////////////////////////////////
@@ -609,9 +613,16 @@ linearizedRodNeumannToDirichletMap(const std::string& rodName,
RodCorrectionType rhs(currentIterate.size());
rhs = 0;
assembler.assembleGradient(currentIterate, rhs);
-
+
// The right hand side is the _negative_ gradient
rhs *= -1;
+
+ // If we do not have a Dirichlet boundary at all we add a scaled
+ // identity matrix. That way the matrix gets elliptic again.
+ // Seems to be a common hack in this situation
+ for (size_t i=0; i<stiffnessMatrix.N(); i++)
+ for (int j=0; j<6; j++)
+ stiffnessMatrix[i][i][j][j] += neumannRegularization_;
/////////////////////////////////////////////////////////////////////
// Turn the input force and torque into a Neumann boundary field
--
GitLab