From 0ec44c7d75b880a902d4af9092328da3f79fa6b0 Mon Sep 17 00:00:00 2001
From: Oliver Sander <sander@igpm.rwth-aachen.de>
Date: Sat, 2 Apr 2011 17:01:22 +0000
Subject: [PATCH] merge the continuum Neumann-to-Dirichlet map from the
Steklov-Poincare step class
[[Imported from SVN: r7063]]
---
.../coupling/rodcontinuumfixedpointstep.hh | 155 ++++++++++--------
1 file changed, 85 insertions(+), 70 deletions(-)
diff --git a/dune/gfe/coupling/rodcontinuumfixedpointstep.hh b/dune/gfe/coupling/rodcontinuumfixedpointstep.hh
index f850ef3f..99ab514b 100644
--- a/dune/gfe/coupling/rodcontinuumfixedpointstep.hh
+++ b/dune/gfe/coupling/rodcontinuumfixedpointstep.hh
@@ -113,10 +113,11 @@ protected:
rodDirichletToNeumannMap(const std::string& rodName,
const std::map<std::pair<std::string,std::string>,RigidBodyMotion<3> >& lambda) const;
- std::map<std::pair<std::string,std::string>,RigidBodyMotion<3>::TangentVector>
- continuumDirichletToNeumannMap(const std::string& continuumName,
- const std::map<std::pair<std::string,std::string>,RigidBodyMotion<3> >& lambda) const;
-
+ std::map<std::pair<std::string,std::string>,RigidBodyMotion<3>::TangentVector>
+ continuumNeumannToDirichletMap(const std::string& continuumName,
+ 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::set<std::string> rodsPerContinuum(const std::string& name) const;
std::set<std::string> continuaPerRod(const std::string& name) const;
@@ -413,106 +414,120 @@ rodDirichletToNeumannMap(const std::string& rodName,
return result;
}
-#if 0
template <class RodGridType, class ContinuumGridType>
-std::map<std::pair<std::string,std::string>,RigidBodyMotion<3>::TangentVector> RodContinuumFixedPointStep<RodGridType,ContinuumGridType>::
-continuumDirichletToNeumannMap(const std::string& continuumName,
- const std::map<std::pair<std::string,std::string>,RigidBodyMotion<3> >& lambda) const
+std::map<std::pair<std::string,std::string>,RigidBodyMotion<3>::TangentVector>
+RodContinuumFixedPointStep<RodGridType,ContinuumGridType>::
+continuumNeumannToDirichletMap(const std::string& continuumName,
+ 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
{
- // Set initial iterate
- VectorType& x3d = continuumSubdomainSolutions_[continuumName];
- x3d.resize(complex_.continuumGrid(continuumName)->size(dim));
- x3d = 0;
+ ////////////////////////////////////////////////////
+ // Assemble the linearized problem
+ ////////////////////////////////////////////////////
- // Copy the true Dirichlet values into it
+ /** \todo We are actually assembling standard linear elasticity,
+ * i.e. the linearization at zero
+ */
+ typedef P1NodalBasis<typename ContinuumGridType::LeafGridView,double> P1Basis;
const LeafBoundaryPatch<ContinuumGridType>& dirichletBoundary = complex_.continuum(continuumName).dirichletBoundary_;
- const VectorType& dirichletValues = complex_.continuum(continuumName).dirichletValues_;
+ P1Basis basis(dirichletBoundary.gridView());
+ OperatorAssembler<P1Basis,P1Basis> assembler(basis, basis);
+
+ MatrixType stiffnessMatrix;
+ assembler.assemble(*continuum(continuumName).localAssembler_, stiffnessMatrix);
- for (size_t i=0; i<x3d.size(); i++)
- if (dirichletBoundary.containsVertex(i))
- x3d[i] = dirichletValues[i];
- typename std::map<std::pair<std::string,std::string>,RigidBodyMotion<3> >::const_iterator it = lambda.begin();
- for (; it!=lambda.end(); ++it) {
+ /////////////////////////////////////////////////////////////////////
+ // Turn the input force and torque into a Neumann boundary field
+ /////////////////////////////////////////////////////////////////////
+
+ // The weak form of the volume and Neumann data
+ /** \brief Not implemented yet! */
+ VectorType rhs(complex_.continuumGrid(continuumName)->size(dim));
+ rhs = 0;
+
+ 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.second != continuumName)
continue;
-
- // Turn \lambda \in TSE(3) into a Dirichlet value for the continuum
+
+ // Use 'forceTorque' as a Neumann value for the rod
const LeafBoundaryPatch<ContinuumGridType>& interfaceBoundary = complex_.coupling(couplingName).continuumInterfaceBoundary_;
- const RigidBodyMotion<dim>& referenceInterface = complex_.coupling(couplingName).referenceInterface_;
- setRotation(interfaceBoundary, x3d, referenceInterface, it->second);
-
+ //
+ VectorType localNeumannValues;
+ computeAveragePressure<typename ContinuumGridType::LeafGridView>(forceTorque.find(couplingName)->second,
+ interfaceBoundary,
+ centerOfTorque.find(couplingName)->second.r,
+ localNeumannValues);
+
+ BoundaryFunctionalAssembler<P1Basis> boundaryFunctionalAssembler(basis, interfaceBoundary);
+ BasisGridFunction<P1Basis,VectorType> neumannValuesFunction(basis,localNeumannValues);
+ NeumannBoundaryAssembler<ContinuumGridType, Dune::FieldVector<double,3> > localNeumannAssembler(neumannValuesFunction);
+ boundaryFunctionalAssembler.assemble(localNeumannAssembler, rhs, false);
+
}
// Set the correct Dirichlet nodes
+ /** \brief Don't write this BitSetVector at each iteration */
+ Dune::BitSetVector<dim> dirichletNodes(rhs.size(),false);
+ for (size_t i=0; i<dirichletNodes.size(); i++)
+ dirichletNodes[i] = dirichletBoundary.containsVertex(i);
dynamic_cast<IterationStep<VectorType>* >(continuum(continuumName).solver_->iterationStep_)->ignoreNodes_
- = &continuum(continuumName).dirichletAndCouplingNodes_;
-
- // Right hand side vector: currently without Neumann and volume terms
- VectorType rhs3d(x3d.size());
- rhs3d = 0;
+ = &dirichletNodes;
+
+ // Initial iterate is 0, all Dirichlet values are 0 (because we solve a correction problem
+ VectorType x(dirichletNodes.size());
+ x = 0;
- // Solve the Dirichlet problem
assert( (dynamic_cast<LinearIterationStep<MatrixType,VectorType>* >(continuum(continuumName).solver_->iterationStep_)) );
- dynamic_cast<LinearIterationStep<MatrixType,VectorType>* >(continuum(continuumName).solver_->iterationStep_)->setProblem(*continuum(continuumName).stiffnessMatrix_, x3d, rhs3d);
+ dynamic_cast<LinearIterationStep<MatrixType,VectorType>* >(continuum(continuumName).solver_->iterationStep_)->setProblem(stiffnessMatrix, x, rhs);
- continuum(continuumName).solver_->preprocess();
+ //solver.preprocess();
continuum(continuumName).solver_->solve();
- x3d = dynamic_cast<IterationStep<VectorType>* >(continuum(continuumName).solver_->iterationStep_)->getSol();
-
- // Integrate over the residual on the coupling boundary to obtain
- // an element of T^*SE.
- Dune::FieldVector<double,3> continuumForce, continuumTorque;
-
- VectorType residual = rhs3d;
- continuum(continuumName).stiffnessMatrix_->mmv(x3d,residual);
-
- //////////////////////////////////////////////////////////////////////////////
- // Extract the residual stresses
- //////////////////////////////////////////////////////////////////////////////
-
- std::map<std::pair<std::string,std::string>, RigidBodyMotion<3>::TangentVector > result;
+ x = continuum(continuumName).solver_->iterationStep_->getSol();
- for (it = lambda.begin(); it!=lambda.end(); ++it) {
+ /////////////////////////////////////////////////////////////////////////////////
+ // Average the continuum displacement on the coupling boundary
+ /////////////////////////////////////////////////////////////////////////////////
+ std::map<std::pair<std::string,std::string>,RigidBodyMotion<3>::TangentVector> interfaceCorrection;
+
+ for (ForceIterator it = forceTorque.begin(); it!=forceTorque.end(); ++it) {
const std::pair<std::string,std::string>& couplingName = it->first;
-
+
if (couplingName.second != continuumName)
continue;
+ // Use 'forceTorque' as a Neumann value for the rod
const LeafBoundaryPatch<ContinuumGridType>& interfaceBoundary = complex_.coupling(couplingName).continuumInterfaceBoundary_;
-
- VectorType neumannForces(residual.size());
- neumannForces = 0;
- weakToStrongBoundaryStress(interfaceBoundary, residual, neumannForces);
+ RigidBodyMotion<3> averageInterface;
+ computeAverageInterface(interfaceBoundary, x, 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
+ */
+ interfaceCorrection[couplingName][0] = averageInterface.r[0];
+ interfaceCorrection[couplingName][1] = averageInterface.r[1];
+ interfaceCorrection[couplingName][2] = averageInterface.r[2];
+
+ Dune::FieldVector<double,3> infinitesimalRotation = Rotation<3,double>::expInv(averageInterface.q);
+ interfaceCorrection[couplingName][3] = infinitesimalRotation[0];
+ interfaceCorrection[couplingName][4] = infinitesimalRotation[1];
+ interfaceCorrection[couplingName][5] = infinitesimalRotation[2];
- /** \todo Is referenceInterface.r the correct center of rotation? */
- const RigidBodyMotion<dim>& referenceInterface = complex_.coupling(couplingName).referenceInterface_;
-
- computeTotalForceAndTorque(interfaceBoundary,
- neumannForces,
- referenceInterface.r,
- continuumForce, continuumTorque);
-
- result[couplingName][0] = continuumForce[0];
- result[couplingName][1] = continuumForce[1];
- result[couplingName][2] = continuumForce[2];
- result[couplingName][3] = continuumTorque[0];
- result[couplingName][4] = continuumTorque[1];
- result[couplingName][5] = continuumTorque[2];
-
}
- return result;
+ return interfaceCorrection;
}
-#endif
/** \brief One preconditioned Richardson step
--
GitLab