From 92b8b9c36da0459700d669f8f8afb0ac6c126aef Mon Sep 17 00:00:00 2001
From: Oliver Sander <sander@igpm.rwth-aachen.de>
Date: Mon, 24 Jan 2011 08:03:56 +0000
Subject: [PATCH] make the linearized NtD map work for rods with zero, one, or
two couplings
[[Imported from SVN: r6855]]
---
.../rodcontinuumsteklovpoincarestep.hh | 62 +++++++++++++++----
1 file changed, 49 insertions(+), 13 deletions(-)
diff --git a/dune/gfe/coupling/rodcontinuumsteklovpoincarestep.hh b/dune/gfe/coupling/rodcontinuumsteklovpoincarestep.hh
index a8dd0923..56b27a43 100644
--- a/dune/gfe/coupling/rodcontinuumsteklovpoincarestep.hh
+++ b/dune/gfe/coupling/rodcontinuumsteklovpoincarestep.hh
@@ -497,16 +497,14 @@ linearizedRodNeumannToDirichletMap(const std::string& rodName,
const std::map<std::pair<std::string,std::string>,RigidBodyMotion<3>::TangentVector>& forceTorque,
const std::map<std::pair<std::string,std::string>,RigidBodyMotion<3> >& centerOfTorque) const
{
- std::pair<std::string,std::string> interfaceName = std::make_pair(rodName,"continuum");
-
////////////////////////////////////////////////////
// Assemble the linearized rod problem
////////////////////////////////////////////////////
- const LeafBoundaryPatch<RodGridType>& interface = complex_.coupling(interfaceName).rodInterfaceBoundary_;
+ const LeafBoundaryPatch<RodGridType>& dirichletBoundary = complex_.rod(rodName).dirichletBoundary_;
typedef Dune::BCRSMatrix<Dune::FieldMatrix<double,6,6> > MatrixType;
- GeodesicFEAssembler<typename RodGridType::LeafGridView, RigidBodyMotion<3> > assembler(interface.gridView(),
+ GeodesicFEAssembler<typename RodGridType::LeafGridView, RigidBodyMotion<3> > assembler(dirichletBoundary.gridView(),
rod(rodName).localStiffness_);
MatrixType stiffnessMatrix;
assembler.assembleMatrix(currentIterate,
@@ -525,17 +523,36 @@ linearizedRodNeumannToDirichletMap(const std::string& rodName,
// Turn the input force and torque into a Neumann boundary field
/////////////////////////////////////////////////////////////////////
- // The weak form of the Neumann data
- rhs[0] += forceTorque.find(interfaceName)->second;
+ typedef typename std::map<std::pair<std::string,std::string>,RigidBodyMotion<3>::TangentVector>::const_iterator ForceIterator;
+
+ for (ForceIterator it = forceTorque.begin(); it!=forceTorque.end(); ++it) {
+
+ const std::pair<std::string,std::string>& couplingName = it->first;
+
+ if (couplingName.first != rodName)
+ continue;
+
+ // Use 'forceTorque' as a Neumann value for the rod
+ const LeafBoundaryPatch<RodGridType>& interfaceBoundary = complex_.coupling(couplingName).rodInterfaceBoundary_;
+
+ /** \todo Use the BoundaryPatch iterator, which will be a lot faster
+ * once we use EntitySeed for its implementation
+ */
+ for (size_t i=0; i<rhs.size(); i++)
+ if (interfaceBoundary.containsVertex(i))
+ rhs[i] += it->second;
+ }
///////////////////////////////////////////////////////////
// Modify matrix and rhs to contain one Dirichlet node
///////////////////////////////////////////////////////////
- int dIdx = rhs.size()-1; // hardwired index of the single Dirichlet node
- rhs[dIdx] = 0;
- stiffnessMatrix[dIdx] = 0;
- stiffnessMatrix[dIdx][dIdx] = Dune::ScaledIdentityMatrix<double,6>(1);
+ for (size_t i=0; i<rhs.size(); i++)
+ if (dirichletBoundary.containsVertex(i)) {
+ rhs[i] = 0;
+ stiffnessMatrix[i] = 0;
+ stiffnessMatrix[i][i] = Dune::ScaledIdentityMatrix<double,6>(1);
+ }
//////////////////////////////////////////////////
// Solve the Neumann problem
@@ -559,10 +576,29 @@ linearizedRodNeumannToDirichletMap(const std::string& rodName,
// Solve!
cg.apply(x, rhs, statistics);
- std::cout << "Linear rod interface correction: " << x[0] << std::endl;
-
+ ///////////////////////////////////////////////////////////////////////////////////////
+ // Extract the solution at the interface boundaries
+ ///////////////////////////////////////////////////////////////////////////////////////
std::map<std::pair<std::string,std::string>,RigidBodyMotion<3>::TangentVector> interfaceCorrection;
- interfaceCorrection[interfaceName] = x[0];
+
+ for (ForceIterator it = forceTorque.begin(); it!=forceTorque.end(); ++it) {
+
+ const std::pair<std::string,std::string>& couplingName = it->first;
+
+ if (couplingName.first != rodName)
+ continue;
+
+ // Use 'forceTorque' as a Neumann value for the rod
+ const LeafBoundaryPatch<RodGridType>& interfaceBoundary = complex_.coupling(couplingName).rodInterfaceBoundary_;
+
+ /** \todo Use the BoundaryPatch iterator, which will be a lot faster
+ * once we use EntitySeed for its implementation
+ * \todo We assume here that a coupling boundary consists of a single point only (not two)
+ */
+ for (size_t i=0; i<rhs.size(); i++)
+ if (interfaceBoundary.containsVertex(i))
+ interfaceCorrection[couplingName] = rhs[i];
+ }
return interfaceCorrection;
}
--
GitLab