From edb6204b55506778acf0af62b06b4a1d64930654 Mon Sep 17 00:00:00 2001
From: Oliver Sander <sander@igpm.rwth-aachen.de>
Date: Sun, 23 Jan 2011 11:06:53 +0000
Subject: [PATCH] Change interface of the Steklov-Poincare step class to allow
for more than one coupling
[[Imported from SVN: r6831]]
---
dirneucoupling.cc | 21 +++++----
.../rodcontinuumsteklovpoincarestep.hh | 45 +++++++++++--------
2 files changed, 39 insertions(+), 27 deletions(-)
diff --git a/dirneucoupling.cc b/dirneucoupling.cc
index cbcb5e53..357c2457 100644
--- a/dirneucoupling.cc
+++ b/dirneucoupling.cc
@@ -342,10 +342,13 @@ int main (int argc, char *argv[]) try
// Dirichlet-Neumann Solver
// /////////////////////////////////////////////////////
- // Init interface value
RigidBodyMotion<3> referenceInterface = rodX[0];
- complex.couplings_[std::make_pair("rod","continuum")].referenceInterface_ = referenceInterface;
- RigidBodyMotion<3> lambda = referenceInterface;
+ complex.couplings_[interfaceName].referenceInterface_ = referenceInterface;
+
+ // Init interface value
+ std::map<std::pair<std::string,std::string>, RigidBodyMotion<3> > lambda;
+ lambda[interfaceName] = referenceInterface;
+
FieldVector<double,3> lambdaForce(0);
FieldVector<double,3> lambdaTorque(0);
@@ -359,7 +362,7 @@ int main (int argc, char *argv[]) try
std::cout << "----------------------------------------------------" << std::endl;
// Backup of the current iterate for the error computation later on
- RigidBodyMotion<3> oldLambda = lambda;
+ std::map<std::pair<std::string,std::string>, RigidBodyMotion<3> > oldLambda = lambda;
if (ddType=="FixedPointIteration") {
@@ -367,7 +370,7 @@ int main (int argc, char *argv[]) try
// Dirichlet step for the rod
// //////////////////////////////////////////////////
- rodX[0] = lambda;
+ rodX[0] = lambda[interfaceName];
rodSolver.setInitialSolution(rodX);
rodSolver.solve();
@@ -436,10 +439,10 @@ int main (int argc, char *argv[]) try
// Compute new damped interface value
//////////////////////////////////////////////////////////////
for (int j=0; j<dim; j++)
- lambda.r[j] = (1-damping) * lambda.r[j]
+ lambda[interfaceName].r[j] = (1-damping) * lambda[interfaceName].r[j]
+ damping * (referenceInterface.r[j] + averageInterface.r[j]);
- lambda.q = Rotation<3,double>::interpolate(lambda.q,
+ lambda[interfaceName].q = Rotation<3,double>::interpolate(lambda[interfaceName].q,
referenceInterface.q.mult(averageInterface.q),
damping);
@@ -469,7 +472,7 @@ int main (int argc, char *argv[]) try
} else
DUNE_THROW(NotImplemented, ddType << " is not a known domain decomposition algorithm");
- std::cout << "Lambda: " << lambda << std::endl;
+ std::cout << "Lambda: " << lambda[interfaceName] << std::endl;
// ////////////////////////////////////////////////////////////////////////
// Write the two iterates to disk for later convergence rate measurement
@@ -507,7 +510,7 @@ int main (int argc, char *argv[]) try
// Compute error in the energy norm
// ////////////////////////////////////////////
- double lengthOfCorrection = RigidBodyMotion<3>::distance(oldLambda, lambda);
+ double lengthOfCorrection = RigidBodyMotion<3>::distance(oldLambda[interfaceName], lambda[interfaceName]);
double convRate = lengthOfCorrection / normOfOldCorrection;
diff --git a/dune/gfe/coupling/rodcontinuumsteklovpoincarestep.hh b/dune/gfe/coupling/rodcontinuumsteklovpoincarestep.hh
index 98bf409e..7c5cfb73 100644
--- a/dune/gfe/coupling/rodcontinuumsteklovpoincarestep.hh
+++ b/dune/gfe/coupling/rodcontinuumsteklovpoincarestep.hh
@@ -312,7 +312,7 @@ public:
/** \brief Do one Steklov-Poincare step
* \param[in,out] lambda The old and new iterate
*/
- void iterate(RigidBodyMotion<3>& lambda);
+ void iterate(std::map<std::pair<std::string,std::string>, RigidBodyMotion<3> >& lambda);
private:
@@ -599,13 +599,17 @@ continuumDirichletToNeumannMap(const std::string& continuumName,
/** \brief One preconditioned Richardson step
*/
template <class RodGridType, class ContinuumGridType>
-void RodContinuumSteklovPoincareStep<RodGridType,ContinuumGridType>::iterate(RigidBodyMotion<3>& lambda)
+void RodContinuumSteklovPoincareStep<RodGridType,ContinuumGridType>::
+iterate(std::map<std::pair<std::string,std::string>, RigidBodyMotion<3> >& lambda)
{
+ // temporary
+ std::pair<std::string,std::string> interfaceName = std::make_pair("rod","continuum");
+
///////////////////////////////////////////////////////////////////
// Evaluate the Dirichlet-to-Neumann map for the rod
///////////////////////////////////////////////////////////////////
- RigidBodyMotion<3>::TangentVector rodForceTorque = rodDirichletToNeumannMap(lambda);
+ RigidBodyMotion<3>::TangentVector rodForceTorque = rodDirichletToNeumannMap(lambda[interfaceName]);
std::cout << "resultant rod force and torque: " << rodForceTorque << std::endl;
@@ -613,10 +617,7 @@ void RodContinuumSteklovPoincareStep<RodGridType,ContinuumGridType>::iterate(Rig
// Evaluate the Dirichlet-to-Neumann map for the continuum
///////////////////////////////////////////////////////////////////
- std::map<std::pair<std::string,std::string>, RigidBodyMotion<3> > tmpLambda;
- tmpLambda[std::make_pair("rod","continuum")] = lambda;
-
- std::map<std::pair<std::string,std::string>, RigidBodyMotion<3>::TangentVector > tmpContinuumForceTorque = continuumDirichletToNeumannMap("continuum", tmpLambda);
+ std::map<std::pair<std::string,std::string>, RigidBodyMotion<3>::TangentVector > tmpContinuumForceTorque = continuumDirichletToNeumannMap("continuum", lambda);
RigidBodyMotion<3>::TangentVector continuumForceTorque = tmpContinuumForceTorque[std::make_pair("rod","continuum")];
std::cout << "resultant continuum force and torque: " << continuumForceTorque << std::endl;
@@ -636,7 +637,6 @@ void RodContinuumSteklovPoincareStep<RodGridType,ContinuumGridType>::iterate(Rig
///////////////////////////////////////////////////////////////
Dune::FieldVector<double,6> interfaceCorrection;
- std::pair<std::string,std::string> couplingName = std::make_pair("rod","continuum");
if (preconditioner_=="DirichletNeumann") {
@@ -650,13 +650,15 @@ void RodContinuumSteklovPoincareStep<RodGridType,ContinuumGridType>::iterate(Rig
rhs3d = 0;
LinearizedContinuumNeumannToDirichletMap<typename ContinuumGridType::LeafGridView,MatrixType,VectorType>
- linContNtDMap(complex_.coupling(couplingName).continuumInterfaceBoundary_,
+ linContNtDMap(complex_.coupling(interfaceName).continuumInterfaceBoundary_,
rhs3d,
*continuum("continuum").dirichletValues_,
continuum("continuum").localAssembler_,
continuum("continuum").solver_);
- interfaceCorrection = linContNtDMap.apply(continuumSubdomainSolutions_["continuum"], residualForceTorque, lambda.r);
+ interfaceCorrection = linContNtDMap.apply(continuumSubdomainSolutions_["continuum"],
+ residualForceTorque,
+ lambda[std::make_pair("rod","continuum")].r);
} else if (preconditioner_=="NeumannDirichlet") {
@@ -665,10 +667,12 @@ void RodContinuumSteklovPoincareStep<RodGridType,ContinuumGridType>::iterate(Rig
// of the rod.
////////////////////////////////////////////////////////////////////
- LinearizedRodNeumannToDirichletMap<typename RodGridType::LeafGridView,RodCorrectionType> linRodNtDMap(complex_.coupling(couplingName).rodInterfaceBoundary_,
+ LinearizedRodNeumannToDirichletMap<typename RodGridType::LeafGridView,RodCorrectionType> linRodNtDMap(complex_.coupling(interfaceName).rodInterfaceBoundary_,
rods_["rod"].localStiffness_);
- interfaceCorrection = linRodNtDMap.apply(rodSubdomainSolutions_["rod"], residualForceTorque, lambda.r);
+ interfaceCorrection = linRodNtDMap.apply(rodSubdomainSolutions_["rod"],
+ residualForceTorque,
+ lambda[std::make_pair("rod","continuum")].r);
} else if (preconditioner_=="NeumannNeumann") {
@@ -684,18 +688,21 @@ void RodContinuumSteklovPoincareStep<RodGridType,ContinuumGridType>::iterate(Rig
rhs3d = 0;
LinearizedContinuumNeumannToDirichletMap<typename ContinuumGridType::LeafGridView,MatrixType,VectorType>
- linContNtDMap(complex_.coupling(couplingName).continuumInterfaceBoundary_,
+ linContNtDMap(complex_.coupling(interfaceName).continuumInterfaceBoundary_,
rhs3d,
*continuum("continuum").dirichletValues_,
continuum("continuum").localAssembler_,
continuum("continuum").solver_);
- LinearizedRodNeumannToDirichletMap<typename RodGridType::LeafGridView,RodCorrectionType> linRodNtDMap(complex_.coupling(couplingName).rodInterfaceBoundary_,
+ LinearizedRodNeumannToDirichletMap<typename RodGridType::LeafGridView,RodCorrectionType> linRodNtDMap(complex_.coupling(interfaceName).rodInterfaceBoundary_,
rods_["rod"].localStiffness_);
- Dune::FieldVector<double,6> continuumCorrection = linContNtDMap.apply(continuumSubdomainSolutions_["continuum"], residualForceTorque, lambda.r);
+ Dune::FieldVector<double,6> continuumCorrection = linContNtDMap.apply(continuumSubdomainSolutions_["continuum"],
+ residualForceTorque,
+ lambda[std::make_pair("rod","continuum")].r);
Dune::FieldVector<double,6> rodCorrection = linRodNtDMap.apply(rodSubdomainSolutions_["rod"],
- residualForceTorque, lambda.r);
+ residualForceTorque,
+ lambda[std::make_pair("rod","continuum")].r);
for (int j=0; j<6; j++)
interfaceCorrection[j] = (alpha_[0] * continuumCorrection[j] + alpha_[1] * rodCorrection[j])
@@ -713,8 +720,10 @@ void RodContinuumSteklovPoincareStep<RodGridType,ContinuumGridType>::iterate(Rig
///////////////////////////////////////////////////////////////////////////////
interfaceCorrection *= richardsonDamping_;
- lambda = RigidBodyMotion<3>::exp(lambda, interfaceCorrection);
-
+ typename std::map<std::pair<std::string,std::string>,RigidBodyMotion<3> >::iterator it = lambda.begin();
+ for (; it!=lambda.end(); ++it) {
+ it->second = RigidBodyMotion<3>::exp(it->second, interfaceCorrection);
+ }
}
#endif
\ No newline at end of file
--
GitLab