Skip to content
Snippets Groups Projects
Commit edb6204b authored by Oliver Sander's avatar Oliver Sander Committed by sander@FU-BERLIN.DE
Browse files

Change interface of the Steklov-Poincare step class to allow for more than one coupling 

[[Imported from SVN: r6831]]
parent 0f8a4d32
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
dirneucoupling.cc
+ 12
− 9
View file @ edb6204b
......@@ -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;
......
dune/gfe/coupling/rodcontinuumsteklovpoincarestep.hh
+ 27
− 18
View file @ edb6204b
......@@ -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
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment