From ec2064185952ec11d113692d8990e2cacfbc630a Mon Sep 17 00:00:00 2001
From: Oliver Sander <sander@igpm.rwth-aachen.de>
Date: Mon, 31 Jan 2011 16:57:51 +0000
Subject: [PATCH] moved the part of the Steklov-Poincare step that can handle
continuum--continuum contact into a separate class and file
[[Imported from SVN: r6918]]
---
.../rodcontinuumsteklovpoincarestep.hh | 442 +-----------------
1 file changed, 4 insertions(+), 438 deletions(-)
diff --git a/dune/gfe/coupling/rodcontinuumsteklovpoincarestep.hh b/dune/gfe/coupling/rodcontinuumsteklovpoincarestep.hh
index bde6316a..df60cdc5 100644
--- a/dune/gfe/coupling/rodcontinuumsteklovpoincarestep.hh
+++ b/dune/gfe/coupling/rodcontinuumsteklovpoincarestep.hh
@@ -14,10 +14,6 @@
#include <dune/fufem/assemblers/boundaryfunctionalassembler.hh>
#include <dune/fufem/assemblers/localassemblers/neumannboundaryassembler.hh>
-#if HAVE_DUNE_CONTACT
-#include <dune/contact/nbodyassembler.hh>
-#endif
-
#include <dune/gfe/coupling/rodcontinuumcomplex.hh>
@@ -146,19 +142,7 @@ public:
*/
void iterate(std::map<std::pair<std::string,std::string>, RigidBodyMotion<3> >& lambda);
-#if HAVE_DUNE_CONTACT
- /** \brief Do one Steklov-Poincare step
- * \param[in,out] lambda The old and new iterate
- */
- void iterateWithContact(std::map<std::pair<std::string,std::string>, RigidBodyMotion<3> >& lambda,
- LoopSolver<VectorType>* contactSolver,
- const MatrixType* totalStiffnessMatrix,
- const NBodyAssembler<ContinuumGridType, VectorType>* contactAssembler,
- const std::vector<std::string>& continuumName
- );
-#endif
-
-private:
+protected:
std::map<std::pair<std::string,std::string>,RigidBodyMotion<3>::TangentVector>
rodDirichletToNeumannMap(const std::string& rodName,
@@ -222,7 +206,7 @@ private:
public:
double neumannRegularization_;
-private:
+protected:
//////////////////////////////////////////////////////////////////
// Data members related to the rod problems
@@ -253,7 +237,7 @@ private:
public:
/** \todo Should be part of RodData, too */
mutable std::map<std::string, RodConfigurationType> rodSubdomainSolutions_;
-private:
+protected:
//////////////////////////////////////////////////////////////////
// Data members related to the continuum problems
//////////////////////////////////////////////////////////////////
@@ -286,7 +270,7 @@ private:
public:
/** \todo Should be part of ContinuumData, too */
mutable std::map<std::string, ContinuumConfigurationType> continuumSubdomainSolutions_;
-private:
+
};
@@ -1060,422 +1044,4 @@ iterate(std::map<std::pair<std::string,std::string>, RigidBodyMotion<3> >& lambd
}
}
-
-#if HAVE_DUNE_CONTACT
-/** \brief One preconditioned Richardson step plus a continuum contact problem
-*/
-template <class RodGridType, class ContinuumGridType>
-void RodContinuumSteklovPoincareStep<RodGridType,ContinuumGridType>::
-iterateWithContact(std::map<std::pair<std::string,std::string>, RigidBodyMotion<3> >& lambda,
- LoopSolver<VectorType>* contactSolver,
- const MatrixType* totalStiffnessMatrix,
- const NBodyAssembler<ContinuumGridType, VectorType>* contactAssembler,
- const std::vector<std::string>& continuumName
-)
-{
- ///////////////////////////////////////////////////////////////////
- // Evaluate the Dirichlet-to-Neumann maps for the rods
- ///////////////////////////////////////////////////////////////////
-
- std::map<std::pair<std::string,std::string>, RigidBodyMotion<3>::TangentVector> rodForceTorque;
-
- for (RodIterator it = rods_.begin(); it != rods_.end(); ++it) {
-
- const std::string& rodName = it->first;
-
- std::map<std::pair<std::string,std::string>, RigidBodyMotion<3>::TangentVector> forceTorque = rodDirichletToNeumannMap(rodName, lambda);
-
- insert(rodForceTorque, forceTorque);
-
- }
-
- std::cout << "resultant rod forces and torques: " << std::endl;
- typedef typename std::map<std::pair<std::string,std::string>, RigidBodyMotion<3>::TangentVector>::iterator ForceIterator;
- for (ForceIterator it = rodForceTorque.begin(); it != rodForceTorque.end(); ++it)
- std::cout << " [" << it->first.first << ", " << it->first.second << "] -- "
- << it->second << std::endl;
-
- ///////////////////////////////////////////////////////////////////
- // Evaluate the Dirichlet-to-Neumann map for the continuum
- ///////////////////////////////////////////////////////////////////
-
- std::map<std::pair<std::string,std::string>, RigidBodyMotion<3>::TangentVector> continuumForceTorque;
-
- LinearIterationStep<MatrixType,VectorType>* multigridStep = dynamic_cast<LinearIterationStep<MatrixType,VectorType>*>(contactSolver->iterationStep_);
-
- // Make initial iterate: we start from zero
- for (typename RodContinuumComplex<RodGridType,ContinuumGridType>::ConstContinuumIterator it = complex_.continua_.begin();
- it != complex_.continua_.end(); ++it) {
-
- // Copy the true Dirichlet values into it
- const LeafBoundaryPatch<ContinuumGridType>& dirichletBoundary = complex_.continuum(it->first).dirichletBoundary_;
- const VectorType& dirichletValues = complex_.continuum(it->first).dirichletValues_;
-
- VectorType& thisX = continuumSubdomainSolutions_[it->first];
-
- thisX.resize(dirichletValues.size());
- thisX = 0;
-
- for (size_t i=0; i<thisX.size(); i++)
- if (dirichletBoundary.containsVertex(i))
- thisX[i] = dirichletValues[i];
-
- }
-
- typename std::map<std::pair<std::string,std::string>,RigidBodyMotion<3> >::const_iterator it = lambda.begin();
- for (; it!=lambda.end(); ++it) {
-
- const std::pair<std::string,std::string>& couplingName = it->first;
-
- // Turn \lambda \in TSE(3) into a Dirichlet value for the continuum
- const LeafBoundaryPatch<ContinuumGridType>& interfaceBoundary = complex_.coupling(couplingName).continuumInterfaceBoundary_;
- const RigidBodyMotion<dim>& referenceInterface = complex_.coupling(couplingName).referenceInterface_;
-
- setRotation(interfaceBoundary, continuumSubdomainSolutions_[couplingName.second], referenceInterface, it->second);
-
- }
-
- VectorType totalX3d;
- std::vector<VectorType> xVector(continuumName.size());
- for (int i=0; i<continuumName.size(); i++)
- xVector[i] = continuumSubdomainSolutions_[continuumName[i]];
- NBodyAssembler<ContinuumGridType, VectorType>::concatenateVectors(xVector, totalX3d);
-
- // Make the set off all Dirichlet nodes
- Dune::BitSetVector<dim> totalDirichletNodes(totalStiffnessMatrix->N());
-
- int offset = 0;
- for (int i=0; i<continuumName.size(); i++) {
- const Dune::BitSetVector<dim>& ignoreNodes = continuum(continuumName[i]).dirichletAndCouplingNodes_;
- for (int j=0; j<ignoreNodes.size(); j++)
- totalDirichletNodes[offset + j] = ignoreNodes[j];
- offset += ignoreNodes.size();
- }
-
-
- // separate, untransformed weak volume and Neumann terms
- /** \todo Not implemented yet */
- std::map<std::string,VectorType> rhs3d;
- for (size_t i=0; i<continuumName.size(); i++) {
- rhs3d[continuumName[i]].resize(continuum(continuumName[i]).stiffnessMatrix_->N());
- rhs3d[continuumName[i]] = 0;
- }
- VectorType totalRhs3d(totalStiffnessMatrix->N());
- totalRhs3d = 0;
-
- multigridStep->setProblem(*totalStiffnessMatrix, totalX3d, totalRhs3d);
- multigridStep->ignoreNodes_ = &totalDirichletNodes;
-
- contactSolver->preprocess();
-
- contactSolver->solve();
-
- totalX3d = multigridStep->getSol();
-
- // Separate 3d solution vector
- std::vector<VectorType> x3d(continuumName.size());
- contactAssembler->postprocess(totalX3d, x3d);
-
- // the subdomain solutions in canonical coordinates, stored in a map
- for (size_t i=0; i<x3d.size(); i++)
- continuumSubdomainSolutions_[continuumName[i]] = x3d[i];
-
- //////////////////////////////////////////////////////////////////////////////
- // Extract the residual stresses
- //////////////////////////////////////////////////////////////////////////////
-
- // compute the residual for each continuum
- std::map<std::string,VectorType> residual = rhs3d;
- for (typename std::map<std::string,VectorType>::iterator it = residual.begin(); it != residual.end(); ++it)
- continuum(it->first).stiffnessMatrix_->mmv(continuumSubdomainSolutions_[it->first], it->second);
-
- // Integrate over the residual on the coupling boundary to obtain
- // an element of T^*SE.
- for (typename std::map<std::pair<std::string,std::string>, RigidBodyMotion<3> >::iterator it = lambda.begin(); it!=lambda.end(); ++it) {
-
- const std::pair<std::string,std::string>& couplingName = it->first;
-
- const LeafBoundaryPatch<ContinuumGridType>& interfaceBoundary = complex_.coupling(couplingName).continuumInterfaceBoundary_;
-
- VectorType neumannForces(residual[it->first.second].size());
- neumannForces = 0;
-
- weakToStrongBoundaryStress(interfaceBoundary, residual[it->first.second], neumannForces);
-
- /** \todo Is referenceInterface.r the correct center of rotation? */
- const RigidBodyMotion<dim>& referenceInterface = complex_.coupling(couplingName).referenceInterface_;
-
- Dune::FieldVector<double,3> continuumForce, continuumTorque;
-
- computeTotalForceAndTorque(interfaceBoundary,
- neumannForces,
- complex_.coupling(couplingName).referenceInterface_.r,
- continuumForce, continuumTorque);
-
- continuumForceTorque[couplingName][0] = continuumForce[0];
- continuumForceTorque[couplingName][1] = continuumForce[1];
- continuumForceTorque[couplingName][2] = continuumForce[2];
- continuumForceTorque[couplingName][3] = continuumTorque[0];
- continuumForceTorque[couplingName][4] = continuumTorque[1];
- continuumForceTorque[couplingName][5] = continuumTorque[2];
-
- }
-
-
- std::cout << "resultant continuum force and torque: " << std::endl;
- for (ForceIterator it = continuumForceTorque.begin(); it != continuumForceTorque.end(); ++it)
- std::cout << " [" << it->first.first << ", " << it->first.second << "] -- "
- << it->second << std::endl;
-
- ///////////////////////////////////////////////////////////////
- // Compute the overall Steklov-Poincare residual
- ///////////////////////////////////////////////////////////////
-
- // Flip orientation of all rod forces, to account for opposing normals.
- for (ForceIterator it = rodForceTorque.begin(); it != rodForceTorque.end(); ++it)
- it->second *= -1;
-
- std::map<std::pair<std::string,std::string>, RigidBodyMotion<3>::TangentVector> residualForceTorque = rodForceTorque;
-
- for (ForceIterator it = residualForceTorque.begin(), it2 = continuumForceTorque.begin();
- it != residualForceTorque.end();
- ++it, ++it2) {
- assert(it->first == it2->first);
- it->second += it2->second;
- }
-
- ///////////////////////////////////////////////////////////////
- // Apply the preconditioner
- ///////////////////////////////////////////////////////////////
-
- std::map<std::pair<std::string,std::string>, RigidBodyMotion<3>::TangentVector> continuumCorrection;
- std::map<std::pair<std::string,std::string>, RigidBodyMotion<3>::TangentVector> rodCorrection;
-
- if (preconditioner_=="DirichletNeumann" or preconditioner_=="NeumannNeumann") {
-
- ////////////////////////////////////////////////////////////////////
- // Preconditioner is the linearized Neumann-to-Dirichlet map
- // of the continuum.
- ////////////////////////////////////////////////////////////////////
-
- // Make initial iterate: we start from zero
- std::map<std::string,VectorType> x;
- for (typename RodContinuumComplex<RodGridType,ContinuumGridType>::ConstContinuumIterator it = complex_.continua_.begin();
- it != complex_.continua_.end(); ++it) {
-
- // Copy the true Dirichlet values into it
- const LeafBoundaryPatch<ContinuumGridType>& dirichletBoundary = complex_.continuum(it->first).dirichletBoundary_;
- const VectorType& dirichletValues = complex_.continuum(it->first).dirichletValues_;
-
- VectorType& thisX = x[it->first];
-
- thisX.resize(dirichletValues.size());
- thisX = 0;
-
- for (size_t i=0; i<thisX.size(); i++)
- if (dirichletBoundary.containsVertex(i))
- thisX[i] = dirichletValues[i];
-
- }
-
- VectorType totalX3d;
- std::vector<VectorType> xVector(continuumName.size());
- for (int i=0; i<continuumName.size(); i++)
- xVector[i] = x[continuumName[i]];
- NBodyAssembler<ContinuumGridType, VectorType>::concatenateVectors(xVector, totalX3d);
-
- // Make the set off all Dirichlet nodes
- totalDirichletNodes.unsetAll();
-
- int offset = 0;
- for (int i=0; i<continuumName.size(); i++) {
- const LeafBoundaryPatch<ContinuumGridType>& dirichletBoundary = complex_.continuum(continuumName[i]).dirichletBoundary_;
- int nGridVertices = dirichletBoundary.gridView().indexSet().size(dim);
- for (int j=0; j<nGridVertices; j++)
- totalDirichletNodes[offset + j] = dirichletBoundary.containsVertex(j);
- offset += nGridVertices;
- }
-
-
- // separate, untransformed weak volume and Neumann terms
- /** \todo Not implemented yet */
- std::map<std::string,VectorType> rhs3d;
- for (size_t i=0; i<continuumName.size(); i++) {
- rhs3d[continuumName[i]].resize(continuum(continuumName[i]).stiffnessMatrix_->N());
- rhs3d[continuumName[i]] = 0;
- }
-
- for (ForceIterator it = residualForceTorque.begin(); it!=residualForceTorque.end(); ++it) {
-
- const std::pair<std::string,std::string>& couplingName = it->first;
-
- // Use 'forceTorque' as a Neumann value for the rod
- const LeafBoundaryPatch<ContinuumGridType>& interfaceBoundary = complex_.coupling(couplingName).continuumInterfaceBoundary_;
-
- //
- VectorType localNeumannValues;
- computeAveragePressure<typename ContinuumGridType::LeafGridView>(it->second,
- interfaceBoundary,
- lambda.find(couplingName)->second.r, // center of torque
- localNeumannValues);
-
- typedef P1NodalBasis<typename ContinuumGridType::LeafGridView,double> P1Basis;
- P1Basis basis(interfaceBoundary.gridView());
-
- BoundaryFunctionalAssembler<P1Basis> boundaryFunctionalAssembler(basis, interfaceBoundary);
- BasisGridFunction<P1Basis,VectorType> neumannValuesFunction(basis,localNeumannValues);
- NeumannBoundaryAssembler<ContinuumGridType, Dune::FieldVector<double,3> > localNeumannAssembler(neumannValuesFunction);
- boundaryFunctionalAssembler.assemble(localNeumannAssembler, rhs3d[couplingName.second], false);
-
- }
-
- std::vector<VectorType> rhsVector(continuumName.size());
- for (size_t i=0; i<continuumName.size(); i++)
- rhsVector[i] = rhs3d[continuumName[i]];
-
- NBodyAssembler<ContinuumGridType,VectorType>::concatenateVectors(rhsVector, totalRhs3d);
-
- multigridStep->setProblem(*totalStiffnessMatrix, totalX3d, totalRhs3d);
-
- multigridStep->ignoreNodes_ = &totalDirichletNodes;
-
- contactSolver->preprocess();
-
- contactSolver->solve();
-
- totalX3d = multigridStep->getSol();
-
- // Separate 3d solution vector
- std::vector<VectorType> x3d(continuumName.size());
- contactAssembler->postprocess(totalX3d, x3d);
-
- // the subdomain solutions in canonical coordinates, stored in a map
- for (size_t i=0; i<x3d.size(); i++)
- x[continuumName[i]] = x3d[i];
-
- /////////////////////////////////////////////////////////////////////////////////
- // Average the continuum displacement on the coupling boundary
- /////////////////////////////////////////////////////////////////////////////////
-
- for (typename RodContinuumComplex<RodGridType,ContinuumGridType>::ConstCouplingIterator it = complex_.couplings_.begin();
- it!=complex_.couplings_.end(); ++it) {
-
- const std::pair<std::string,std::string>& couplingName = it->first;
-
- // Use 'forceTorque' as a Neumann value for the rod
- const LeafBoundaryPatch<ContinuumGridType>& interfaceBoundary = complex_.coupling(couplingName).continuumInterfaceBoundary_;
-
- RigidBodyMotion<3> averageInterface;
- computeAverageInterface(interfaceBoundary, x[couplingName.second], averageInterface);
-
- // Compute the linearization
- /** \todo We could loop directly over interfaceCorrection (and save the name lookup)
- * if we could be sure that interfaceCorrection contains all possible entries already
- */
- continuumCorrection[couplingName][0] = averageInterface.r[0];
- continuumCorrection[couplingName][1] = averageInterface.r[1];
- continuumCorrection[couplingName][2] = averageInterface.r[2];
-
- Dune::FieldVector<double,3> infinitesimalRotation = Rotation<3,double>::expInv(averageInterface.q);
- continuumCorrection[couplingName][3] = infinitesimalRotation[0];
- continuumCorrection[couplingName][4] = infinitesimalRotation[1];
- continuumCorrection[couplingName][5] = infinitesimalRotation[2];
-
- }
-
-
- std::cout << "resultant continuum interface corrections: " << std::endl;
- for (ForceIterator it = continuumCorrection.begin(); it != continuumCorrection.end(); ++it)
- std::cout << " [" << it->first.first << ", " << it->first.second << "] -- "
- << it->second << std::endl;
-
-
- }
-
- if (preconditioner_=="NeumannDirichlet" or preconditioner_=="NeumannNeumann") {
-
- ////////////////////////////////////////////////////////////////////
- // Preconditioner is the linearized Neumann-to-Dirichlet map
- // of the rod.
- ////////////////////////////////////////////////////////////////////
-
- for (RodIterator it = rods_.begin(); it != rods_.end(); ++it) {
-
- const std::string& rodName = it->first;
-
- std::map<std::pair<std::string,std::string>, RigidBodyMotion<3>::TangentVector> rodInterfaceCorrection
- = linearizedRodNeumannToDirichletMap(rodName,
- rodSubdomainSolutions_[rodName],
- residualForceTorque,
- lambda);
-
- insert(rodCorrection, rodInterfaceCorrection);
-
- }
-
- std::cout << "resultant rod interface corrections: " << std::endl;
- for (ForceIterator it = rodCorrection.begin(); it != rodCorrection.end(); ++it)
- std::cout << " [" << it->first.first << ", " << it->first.second << "] -- "
- << it->second << std::endl;
-
- }
-
- std::map<std::pair<std::string,std::string>, RigidBodyMotion<3>::TangentVector> interfaceCorrection;
-
- if (preconditioner_=="DirichletNeumann") {
-
- interfaceCorrection = continuumCorrection;
-
- } else if (preconditioner_=="NeumannDirichlet") {
-
- interfaceCorrection = rodCorrection;
-
- } else if (preconditioner_=="NeumannNeumann") {
-
- ////////////////////////////////////////////////////////////////////
- // Preconditioner is a convex combination of the linearized
- // Neumann-to-Dirichlet map of the continuum and the linearized
- // Neumann-to-Dirichlet map of the rod.
- ////////////////////////////////////////////////////////////////////
-
- std::cout << "resultant interface corrections: " << std::endl;
- for (ForceIterator it = rodCorrection.begin(); it != rodCorrection.end(); ++it) {
-
- const std::pair<std::string,std::string> interfaceName = it->first;
-
- std::cout << " [" << interfaceName.first << ", " << interfaceName.second << "]"
- << " -- " << it->second
- << " -- " << continuumCorrection[interfaceName] << std::endl;
-
- // Compute weighted combination of both
- RigidBodyMotion<3>::TangentVector& correction = interfaceCorrection[interfaceName];
- for (int j=0; j<6; j++)
- correction[j] = (alpha_[0] * continuumCorrection[interfaceName][j] + alpha_[1] * rodCorrection[interfaceName][j])
- / (alpha_[0] + alpha_[1]);
-
- }
-
- } else if (preconditioner_=="RobinRobin") {
-
- DUNE_THROW(Dune::NotImplemented, "Robin-Robin preconditioner not implemented yet");
-
- } else
- DUNE_THROW(Dune::NotImplemented, preconditioner_ << " is not a known preconditioner");
-
- ///////////////////////////////////////////////////////////////////////////////
- // Apply the damped correction to the current interface value
- ///////////////////////////////////////////////////////////////////////////////
-
-
- ForceIterator fIt = interfaceCorrection.begin();
- for (typename std::map<std::pair<std::string,std::string>,RigidBodyMotion<3> >::iterator it = lambda.begin();
- it!=lambda.end();
- ++it, ++fIt) {
- assert(it->first == fIt->first);
- fIt->second *= richardsonDamping_;
- it->second = RigidBodyMotion<3>::exp(it->second, fIt->second);
- }
-}
-#endif // HAVE_DUNE_CONTACT
-
#endif
--
GitLab