diff --git a/dune/gfe/coupling/Makefile.am b/dune/gfe/coupling/Makefile.am
index 27bfbff2fd4736a24ebc1ba7d3380e39188e7c96..ac4add804aca983a5b2e1acb4a89108cfd0c2aed 100644
--- a/dune/gfe/coupling/Makefile.am
+++ b/dune/gfe/coupling/Makefile.am
@@ -5,6 +5,7 @@ SUBDIRS =
srcincludedir = $(includedir)/dune/common/coupling
srcinclude_HEADERS = rodcontinuumcomplex.hh \
+ rodcontinuumddstep.hh \
rodcontinuumfixedpointstep.hh \
rodcontinuumsteklovpoincarestep.hh
diff --git a/dune/gfe/coupling/rodcontinuumddstep.hh b/dune/gfe/coupling/rodcontinuumddstep.hh
new file mode 100644
index 0000000000000000000000000000000000000000..cc7de07f657f8c853e6c210767ac6cd9061cd339
--- /dev/null
+++ b/dune/gfe/coupling/rodcontinuumddstep.hh
@@ -0,0 +1,288 @@
+#ifndef ROD_CONTINUUM_DD_STEP_HH
+#define ROD_CONTINUUM_DD_STEP_HH
+
+/** \file
+ * \brief Base class for rod-continuum domain decomposition algorithms
+ */
+
+#include <vector>
+
+#include <dune/fufem/functionspacebases/p1nodalbasis.hh>
+#include <dune/fufem/assemblers/boundaryfunctionalassembler.hh>
+#include <dune/fufem/assemblers/localassemblers/neumannboundaryassembler.hh>
+
+#include <dune/gfe/coupling/rodcontinuumcomplex.hh>
+
+
+/** \brief Base class for rod-continua domain decomposition problems
+ * \tparam RodGridType Dune grid type used for rod grids (must be 1d)
+ * \tparam ContinuumGridType Dune grid type used for the continua grids
+ */
+template <class RodGridType, class ContinuumGridType>
+class RodContinuumDDStep
+{
+ static const int dim = ContinuumGridType::dimension;
+
+ // The type used for rod configurations
+ typedef std::vector<RigidBodyMotion<dim> > RodConfigurationType;
+
+ // The type used for continuum configurations
+ typedef Dune::BlockVector<Dune::FieldVector<double,dim> > VectorType;
+ typedef Dune::BlockVector<Dune::FieldVector<double,dim> > ContinuumConfigurationType;
+
+ typedef Dune::BlockVector<Dune::FieldVector<double,6> > RodCorrectionType;
+
+ typedef Dune::BCRSMatrix<Dune::FieldMatrix<double,3,3> > MatrixType;
+
+ typedef typename P1NodalBasis<typename ContinuumGridType::LeafGridView,double>::LocalFiniteElement ContinuumLocalFiniteElement;
+
+public:
+
+ /** \brief Constructor for a complex with one rod and one continuum */
+ RodContinuumDDStep(const RodContinuumComplex<RodGridType,ContinuumGridType>& complex/*,
+ RodAssembler<typename RodGridType::LeafGridView,3>* rodAssembler,
+ RiemannianTrustRegionSolver<RodGridType,RigidBodyMotion<3> >* rodSolver,
+ const MatrixType* stiffnessMatrix3d,
+ const Dune::shared_ptr< ::LoopSolver<VectorType> > solver*/)
+ : complex_(complex)
+ {
+#if 0
+ rods_["rod"].assembler_ = rodAssembler;
+ rods_["rod"].solver_ = rodSolver;
+
+ continua_["continuum"].stiffnessMatrix_ = stiffnessMatrix3d;
+ continua_["continuum"].solver_ = solver;
+
+ mergeRodDirichletAndCouplingBoundaries();
+ mergeContinuumDirichletAndCouplingBoundaries();
+#endif
+ }
+
+ /** \brief Do one domain decomposition step
+ * \param[in,out] lambda The old and new iterate
+ */
+ virtual void iterate(std::map<std::pair<std::string,std::string>, RigidBodyMotion<3> >& lambda) = 0;
+
+protected:
+
+ std::set<std::string> rodsPerContinuum(const std::string& name) const;
+
+ std::set<std::string> continuaPerRod(const std::string& name) const;
+
+ /** \brief Add the content of one map to another, aborting rather than overwriting stuff
+ */
+ template <class X, class Y>
+ static void insert(std::map<X,Y>& map1, const std::map<X,Y>& map2)
+ {
+ int oldSize = map1.size();
+ map1.insert(map2.begin(), map2.end());
+ assert(map1.size() == oldSize + map2.size());
+ }
+
+ //////////////////////////////////////////////////////////////////
+ // Data members related to the coupled problem
+ //////////////////////////////////////////////////////////////////
+ const RodContinuumComplex<RodGridType,ContinuumGridType>& complex_;
+#if 0
+protected:
+
+ //////////////////////////////////////////////////////////////////
+ // Data members related to the rod problems
+ //////////////////////////////////////////////////////////////////
+
+ struct RodData
+ {
+ Dune::BitSetVector<6> dirichletAndCouplingNodes_;
+
+ RodAssembler<typename RodGridType::LeafGridView,3>* assembler_;
+
+ RodLocalStiffness<typename RodGridType::LeafGridView,double>* localStiffness_;
+
+ RiemannianTrustRegionSolver<RodGridType,RigidBodyMotion<3> >* solver_;
+ };
+
+ /** \brief Simple const access to rods */
+ const RodData& rod(const std::string& name) const
+ {
+ assert(rods_.find(name) != rods_.end());
+ return rods_.find(name)->second;
+ }
+
+ std::map<std::string, RodData> rods_;
+
+ typedef typename std::map<std::string, RodData>::iterator RodIterator;
+#endif
+public:
+ /** \todo Should be part of RodData, too */
+ mutable std::map<std::string, RodConfigurationType> rodSubdomainSolutions_;
+#if 0
+protected:
+ //////////////////////////////////////////////////////////////////
+ // Data members related to the continuum problems
+ //////////////////////////////////////////////////////////////////
+
+ struct ContinuumData
+ {
+ const MatrixType* stiffnessMatrix_;
+
+ Dune::shared_ptr< ::LoopSolver<VectorType> > solver_;
+
+ Dune::BitSetVector<dim> dirichletAndCouplingNodes_;
+
+ LinearLocalAssembler<ContinuumGridType,
+ ContinuumLocalFiniteElement,
+ ContinuumLocalFiniteElement,
+ Dune::FieldMatrix<double,dim,dim> >* localAssembler_;
+ };
+
+ /** \brief Simple const access to continua */
+ const ContinuumData& continuum(const std::string& name) const
+ {
+ assert(continua_.find(name) != continua_.end());
+ return continua_.find(name)->second;
+ }
+
+ std::map<std::string, ContinuumData> continua_;
+
+ typedef typename std::map<std::string, ContinuumData>::iterator ContinuumIterator;
+#endif
+public:
+ /** \todo Should be part of ContinuumData, too */
+ mutable std::map<std::string, ContinuumConfigurationType> continuumSubdomainSolutions_;
+
+};
+
+#if 0
+template <class RodGridType, class ContinuumGridType>
+void RodContinuumFixedPointStep<RodGridType,ContinuumGridType>::
+mergeRodDirichletAndCouplingBoundaries()
+{
+ ////////////////////////////////////////////////////////////////////////////////////
+ // For each rod, merge the Dirichlet boundary with all interface boundaries
+ //
+ // Currently, we can really only solve rod problems with complete Dirichlet
+ // boundary. Hence there are more direct ways to construct the
+ // dirichletAndCouplingNodes field. Yet like to keep the analogy to the continuum
+ // problem. And maybe one day we have a more flexible rod solver, too.
+ ////////////////////////////////////////////////////////////////////////////////////
+
+ for (RodIterator rIt = rods_.begin(); rIt != rods_.end(); ++rIt) {
+
+ // name of the current rod
+ const std::string& name = rIt->first;
+
+ // short-cut to avoid frequent map look-up
+ Dune::BitSetVector<6>& dirichletAndCouplingNodes = rods_[name].dirichletAndCouplingNodes_;
+
+ dirichletAndCouplingNodes.resize(complex_.rodGrid(name)->size(1));
+
+ // first copy the true Dirichlet boundary
+ const LeafBoundaryPatch<RodGridType>& dirichletBoundary = complex_.rods_.find(name)->second.dirichletBoundary_;
+
+ for (int i=0; i<dirichletAndCouplingNodes.size(); i++)
+ dirichletAndCouplingNodes[i] = dirichletBoundary.containsVertex(i);
+
+ // get the names of all the continua that we couple with
+ std::set<std::string> continuumNames = continuaPerRod(name);
+
+ for (std::set<std::string>::const_iterator cIt = continuumNames.begin();
+ cIt != continuumNames.end();
+ ++cIt) {
+
+ const LeafBoundaryPatch<RodGridType>& rodInterfaceBoundary
+ = complex_.coupling(std::make_pair(name,*cIt)).rodInterfaceBoundary_;
+
+ /** \todo Use the BoundaryPatch iterator here, for increased efficiency */
+ for (int i=0; i<dirichletAndCouplingNodes.size(); i++) {
+ bool v = rodInterfaceBoundary.containsVertex(i);
+ for (int j=0; j<6; j++)
+ dirichletAndCouplingNodes[i][j] = dirichletAndCouplingNodes[i][j] or v;
+ }
+
+ }
+
+ // We can only handle rod problems with a full Dirichlet boundary
+ assert(dirichletAndCouplingNodes.count()==12);
+
+ }
+
+}
+
+
+template <class RodGridType, class ContinuumGridType>
+void RodContinuumFixedPointStep<RodGridType,ContinuumGridType>::
+mergeContinuumDirichletAndCouplingBoundaries()
+{
+ ////////////////////////////////////////////////////////////////////////////////////
+ // For each continuum, merge the Dirichlet boundary with all interface boundaries
+ ////////////////////////////////////////////////////////////////////////////////////
+
+ for (ContinuumIterator cIt = continua_.begin(); cIt != continua_.end(); ++cIt) {
+
+ // name of the current continuum
+ const std::string& name = cIt->first;
+
+ // short-cut to avoid frequent map look-up
+ Dune::BitSetVector<dim>& dirichletAndCouplingNodes = continua_[name].dirichletAndCouplingNodes_;
+
+ dirichletAndCouplingNodes.resize(complex_.continuumGrid(name)->size(dim));
+
+ // first copy the true Dirichlet boundary
+ const LeafBoundaryPatch<ContinuumGridType>& dirichletBoundary = complex_.continua_.find(name)->second.dirichletBoundary_;
+
+ for (int i=0; i<dirichletAndCouplingNodes.size(); i++)
+ dirichletAndCouplingNodes[i] = dirichletBoundary.containsVertex(i);
+
+ // get the names of all the rods that we couple with
+ std::set<std::string> rodNames = rodsPerContinuum(name);
+
+ for (std::set<std::string>::const_iterator rIt = rodNames.begin();
+ rIt != rodNames.end();
+ ++rIt) {
+
+ const LeafBoundaryPatch<ContinuumGridType>& continuumInterfaceBoundary
+ = complex_.coupling(std::make_pair(*rIt,name)).continuumInterfaceBoundary_;
+
+ /** \todo Use the BoundaryPatch iterator here, for increased efficiency */
+ for (int i=0; i<dirichletAndCouplingNodes.size(); i++) {
+ bool v = continuumInterfaceBoundary.containsVertex(i);
+ for (int j=0; j<dim; j++)
+ dirichletAndCouplingNodes[i][j] = dirichletAndCouplingNodes[i][j] or v;
+ }
+
+ }
+
+ }
+
+}
+#endif
+
+template <class RodGridType, class ContinuumGridType>
+std::set<std::string> RodContinuumDDStep<RodGridType,ContinuumGridType>::
+rodsPerContinuum(const std::string& name) const
+{
+ std::set<std::string> result;
+
+ for (typename RodContinuumComplex<RodGridType,ContinuumGridType>::ConstCouplingIterator it = complex_.couplings_.begin();
+ it!=complex_.couplings_.end(); ++it)
+ if (it->first.second == name)
+ result.insert(it->first.first);
+
+ return result;
+}
+
+template <class RodGridType, class ContinuumGridType>
+std::set<std::string> RodContinuumDDStep<RodGridType,ContinuumGridType>::
+continuaPerRod(const std::string& name) const
+{
+ std::set<std::string> result;
+
+ for (typename RodContinuumComplex<RodGridType,ContinuumGridType>::ConstCouplingIterator it = complex_.couplings_.begin();
+ it!=complex_.couplings_.end(); ++it)
+ if (it->first.first == name)
+ result.insert(it->first.second);
+
+ return result;
+}
+
+#endif
diff --git a/dune/gfe/coupling/rodcontinuumfixedpointstep.hh b/dune/gfe/coupling/rodcontinuumfixedpointstep.hh
index 5ab435afb0f795ba214782b504ad2db36b7b5feb..637619fb328cd31f67ac28a4643c7dc1e72c0782 100644
--- a/dune/gfe/coupling/rodcontinuumfixedpointstep.hh
+++ b/dune/gfe/coupling/rodcontinuumfixedpointstep.hh
@@ -11,14 +11,17 @@
// #include <dune/istl/bcrsmatrix.hh>
// #include <dune/istl/bvector.hh>
//
+#include <dune/fufem/functionspacebases/p1nodalbasis.hh>
#include <dune/fufem/assemblers/boundaryfunctionalassembler.hh>
#include <dune/fufem/assemblers/localassemblers/neumannboundaryassembler.hh>
#include <dune/gfe/coupling/rodcontinuumcomplex.hh>
+#include <dune/gfe/coupling/rodcontinuumddstep.hh>
template <class RodGridType, class ContinuumGridType>
class RodContinuumFixedPointStep
+: public RodContinuumDDStep<RodGridType,ContinuumGridType>
{
static const int dim = ContinuumGridType::dimension;
@@ -44,7 +47,7 @@ public:
RiemannianTrustRegionSolver<RodGridType,RigidBodyMotion<3> >* rodSolver,
const MatrixType* stiffnessMatrix3d,
const Dune::shared_ptr< ::LoopSolver<VectorType> > solver)
- : complex_(complex),
+ : RodContinuumDDStep<RodGridType,ContinuumGridType>(complex),
damping_(damping)
{
rods_["rod"].assembler_ = rodAssembler;
@@ -64,7 +67,7 @@ public:
const std::map<std::string,RiemannianTrustRegionSolver<RodGridType,RigidBodyMotion<3> >*>& rodSolver,
const std::map<std::string,const MatrixType*>& stiffnessMatrix3d,
const std::map<std::string, const Dune::shared_ptr< ::LoopSolver<VectorType> > >& solver)
- : complex_(complex),
+ : RodContinuumDDStep<RodGridType,ContinuumGridType>(complex),
damping_(damping)
{
///////////////////////////////////////////////////////////////////////////////////
@@ -118,24 +121,9 @@ protected:
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;
-
- /** \brief Add the content of one map to another, aborting rather than overwriting stuff
- */
- template <class X, class Y>
- static void insert(std::map<X,Y>& map1, const std::map<X,Y>& map2)
- {
- int oldSize = map1.size();
- map1.insert(map2.begin(), map2.end());
- assert(map1.size() == oldSize + map2.size());
- }
-
//////////////////////////////////////////////////////////////////
// Data members related to the coupled problem
//////////////////////////////////////////////////////////////////
- const RodContinuumComplex<RodGridType,ContinuumGridType>& complex_;
/** \brief Damping factor */
double damping_;
@@ -168,9 +156,6 @@ protected:
typedef typename std::map<std::string, RodData>::iterator RodIterator;
-public:
- /** \todo Should be part of RodData, too */
- mutable std::map<std::string, RodConfigurationType> rodSubdomainSolutions_;
protected:
//////////////////////////////////////////////////////////////////
// Data members related to the continuum problems
@@ -201,10 +186,6 @@ protected:
typedef typename std::map<std::string, ContinuumData>::iterator ContinuumIterator;
-public:
- /** \todo Should be part of ContinuumData, too */
- mutable std::map<std::string, ContinuumConfigurationType> continuumSubdomainSolutions_;
-
};
@@ -229,23 +210,23 @@ mergeRodDirichletAndCouplingBoundaries()
// short-cut to avoid frequent map look-up
Dune::BitSetVector<6>& dirichletAndCouplingNodes = rods_[name].dirichletAndCouplingNodes_;
- dirichletAndCouplingNodes.resize(complex_.rodGrid(name)->size(1));
+ dirichletAndCouplingNodes.resize(this->complex_.rodGrid(name)->size(1));
// first copy the true Dirichlet boundary
- const LeafBoundaryPatch<RodGridType>& dirichletBoundary = complex_.rods_.find(name)->second.dirichletBoundary_;
+ const LeafBoundaryPatch<RodGridType>& dirichletBoundary = this->complex_.rods_.find(name)->second.dirichletBoundary_;
for (int i=0; i<dirichletAndCouplingNodes.size(); i++)
dirichletAndCouplingNodes[i] = dirichletBoundary.containsVertex(i);
// get the names of all the continua that we couple with
- std::set<std::string> continuumNames = continuaPerRod(name);
+ std::set<std::string> continuumNames = this->continuaPerRod(name);
for (std::set<std::string>::const_iterator cIt = continuumNames.begin();
cIt != continuumNames.end();
++cIt) {
const LeafBoundaryPatch<RodGridType>& rodInterfaceBoundary
- = complex_.coupling(std::make_pair(name,*cIt)).rodInterfaceBoundary_;
+ = this->complex_.coupling(std::make_pair(name,*cIt)).rodInterfaceBoundary_;
/** \todo Use the BoundaryPatch iterator here, for increased efficiency */
for (int i=0; i<dirichletAndCouplingNodes.size(); i++) {
@@ -280,23 +261,23 @@ mergeContinuumDirichletAndCouplingBoundaries()
// short-cut to avoid frequent map look-up
Dune::BitSetVector<dim>& dirichletAndCouplingNodes = continua_[name].dirichletAndCouplingNodes_;
- dirichletAndCouplingNodes.resize(complex_.continuumGrid(name)->size(dim));
+ dirichletAndCouplingNodes.resize(this->complex_.continuumGrid(name)->size(dim));
// first copy the true Dirichlet boundary
- const LeafBoundaryPatch<ContinuumGridType>& dirichletBoundary = complex_.continua_.find(name)->second.dirichletBoundary_;
+ const LeafBoundaryPatch<ContinuumGridType>& dirichletBoundary = this->complex_.continua_.find(name)->second.dirichletBoundary_;
for (int i=0; i<dirichletAndCouplingNodes.size(); i++)
dirichletAndCouplingNodes[i] = dirichletBoundary.containsVertex(i);
// get the names of all the rods that we couple with
- std::set<std::string> rodNames = rodsPerContinuum(name);
+ std::set<std::string> rodNames = this->rodsPerContinuum(name);
for (std::set<std::string>::const_iterator rIt = rodNames.begin();
rIt != rodNames.end();
++rIt) {
const LeafBoundaryPatch<ContinuumGridType>& continuumInterfaceBoundary
- = complex_.coupling(std::make_pair(*rIt,name)).continuumInterfaceBoundary_;
+ = this->complex_.coupling(std::make_pair(*rIt,name)).continuumInterfaceBoundary_;
/** \todo Use the BoundaryPatch iterator here, for increased efficiency */
for (int i=0; i<dirichletAndCouplingNodes.size(); i++) {
@@ -312,48 +293,20 @@ mergeContinuumDirichletAndCouplingBoundaries()
}
-template <class RodGridType, class ContinuumGridType>
-std::set<std::string> RodContinuumFixedPointStep<RodGridType,ContinuumGridType>::
-rodsPerContinuum(const std::string& name) const
-{
- std::set<std::string> result;
-
- for (typename RodContinuumComplex<RodGridType,ContinuumGridType>::ConstCouplingIterator it = complex_.couplings_.begin();
- it!=complex_.couplings_.end(); ++it)
- if (it->first.second == name)
- result.insert(it->first.first);
-
- return result;
-}
-
-template <class RodGridType, class ContinuumGridType>
-std::set<std::string> RodContinuumFixedPointStep<RodGridType,ContinuumGridType>::
-continuaPerRod(const std::string& name) const
-{
- std::set<std::string> result;
-
- for (typename RodContinuumComplex<RodGridType,ContinuumGridType>::ConstCouplingIterator it = complex_.couplings_.begin();
- it!=complex_.couplings_.end(); ++it)
- if (it->first.first == name)
- result.insert(it->first.second);
-
- return result;
-}
-
template <class RodGridType, class ContinuumGridType>
std::map<std::pair<std::string,std::string>,RigidBodyMotion<3>::TangentVector> RodContinuumFixedPointStep<RodGridType,ContinuumGridType>::
rodDirichletToNeumannMap(const std::string& rodName,
const std::map<std::pair<std::string,std::string>,RigidBodyMotion<3> >& lambda) const
{
// container for the subdomain solution
- RodConfigurationType& rodX = rodSubdomainSolutions_[rodName];
- rodX.resize(complex_.rodGrid(rodName)->size(1));
+ RodConfigurationType& rodX = this->rodSubdomainSolutions_[rodName];
+ rodX.resize(this->complex_.rodGrid(rodName)->size(1));
///////////////////////////////////////////////////////////
// Set the complete set of Dirichlet values
///////////////////////////////////////////////////////////
- const LeafBoundaryPatch<RodGridType>& dirichletBoundary = complex_.rod(rodName).dirichletBoundary_;
- const RodConfigurationType& dirichletValues = complex_.rod(rodName).dirichletValues_;
+ const LeafBoundaryPatch<RodGridType>& dirichletBoundary = this->complex_.rod(rodName).dirichletBoundary_;
+ const RodConfigurationType& dirichletValues = this->complex_.rod(rodName).dirichletValues_;
for (size_t i=0; i<rodX.size(); i++)
if (dirichletBoundary.containsVertex(i))
@@ -368,7 +321,7 @@ rodDirichletToNeumannMap(const std::string& rodName,
continue;
// Use \lambda as a Dirichlet value for the rod
- const LeafBoundaryPatch<RodGridType>& interfaceBoundary = complex_.coupling(couplingName).rodInterfaceBoundary_;
+ const LeafBoundaryPatch<RodGridType>& interfaceBoundary = this->complex_.coupling(couplingName).rodInterfaceBoundary_;
/** \todo Use the BoundaryPatch iterator, which will be a lot faster
* once we use EntitySeed for its implementation
@@ -386,7 +339,7 @@ rodDirichletToNeumannMap(const std::string& rodName,
////////////////////////////////////////////////////////////////////////////////
// Set initial iterate by interpolating between the Dirichlet values
- RodFactory<typename RodGridType::LeafGridView> rodFactory(complex_.rodGrid(rodName)->leafView());
+ RodFactory<typename RodGridType::LeafGridView> rodFactory(this->complex_.rodGrid(rodName)->leafView());
rodFactory.create(rodX);
rod(rodName).solver_->setInitialSolution(rodX);
@@ -406,7 +359,7 @@ rodDirichletToNeumannMap(const std::string& rodName,
if (couplingName.first != rodName)
continue;
- const LeafBoundaryPatch<RodGridType>& couplingBoundary = complex_.coupling(couplingName).rodInterfaceBoundary_;
+ const LeafBoundaryPatch<RodGridType>& couplingBoundary = this->complex_.coupling(couplingName).rodInterfaceBoundary_;
result[couplingName] = rod(rodName).assembler_->getResultantForce(couplingBoundary, rodX);
}
@@ -426,7 +379,7 @@ continuumNeumannToDirichletMap(const std::string& continuumName,
////////////////////////////////////////////////////
typedef P1NodalBasis<typename ContinuumGridType::LeafGridView,double> P1Basis;
- const LeafBoundaryPatch<ContinuumGridType>& dirichletBoundary = complex_.continuum(continuumName).dirichletBoundary_;
+ const LeafBoundaryPatch<ContinuumGridType>& dirichletBoundary = this->complex_.continuum(continuumName).dirichletBoundary_;
P1Basis basis(dirichletBoundary.gridView());
const MatrixType& stiffnessMatrix = *continuum(continuumName).stiffnessMatrix_;
@@ -437,7 +390,7 @@ continuumNeumannToDirichletMap(const std::string& continuumName,
// The weak form of the volume and Neumann data
/** \brief Not implemented yet! */
- VectorType rhs(complex_.continuumGrid(continuumName)->size(dim));
+ VectorType rhs(this->complex_.continuumGrid(continuumName)->size(dim));
rhs = 0;
typedef typename std::map<std::pair<std::string,std::string>,RigidBodyMotion<3>::TangentVector>::const_iterator ForceIterator;
@@ -450,7 +403,7 @@ continuumNeumannToDirichletMap(const std::string& continuumName,
continue;
// Use 'forceTorque' as a Neumann value for the rod
- const LeafBoundaryPatch<ContinuumGridType>& interfaceBoundary = complex_.coupling(couplingName).continuumInterfaceBoundary_;
+ const LeafBoundaryPatch<ContinuumGridType>& interfaceBoundary = this->complex_.coupling(couplingName).continuumInterfaceBoundary_;
//
VectorType localNeumannValues;
@@ -500,7 +453,7 @@ continuumNeumannToDirichletMap(const std::string& continuumName,
continue;
// Use 'forceTorque' as a Neumann value for the rod
- const LeafBoundaryPatch<ContinuumGridType>& interfaceBoundary = complex_.coupling(couplingName).continuumInterfaceBoundary_;
+ const LeafBoundaryPatch<ContinuumGridType>& interfaceBoundary = this->complex_.coupling(couplingName).continuumInterfaceBoundary_;
computeAverageInterface(interfaceBoundary, x, averageInterface[couplingName]);
@@ -528,7 +481,7 @@ iterate(std::map<std::pair<std::string,std::string>, RigidBodyMotion<3> >& lambd
std::map<std::pair<std::string,std::string>, RigidBodyMotion<3>::TangentVector> forceTorque = rodDirichletToNeumannMap(rodName, lambda);
- insert(rodForceTorque, forceTorque);
+ this->insert(rodForceTorque, forceTorque);
}
@@ -561,7 +514,7 @@ iterate(std::map<std::pair<std::string,std::string>, RigidBodyMotion<3> >& lambd
rodForceTorque,
lambda);
- insert(averageInterface, localAverageInterface);
+ this->insert(averageInterface, localAverageInterface);
}
@@ -583,7 +536,7 @@ iterate(std::map<std::pair<std::string,std::string>, RigidBodyMotion<3> >& lambd
const std::pair<std::string,std::string>& interfaceName = it->first;
- const RigidBodyMotion<dim>& referenceInterface = complex_.coupling(interfaceName).referenceInterface_;
+ const RigidBodyMotion<dim>& referenceInterface = this->complex_.coupling(interfaceName).referenceInterface_;
for (int j=0; j<dim; j++)
it->second.r[j] = (1-damping_) * it->second.r[j]
diff --git a/dune/gfe/coupling/rodcontinuumsteklovpoincarestep.hh b/dune/gfe/coupling/rodcontinuumsteklovpoincarestep.hh
index 0d80155efc847e0d72efb772b592552f1c8220e2..53c506df54889f141dca84a4bc3c9fa82c768b59 100644
--- a/dune/gfe/coupling/rodcontinuumsteklovpoincarestep.hh
+++ b/dune/gfe/coupling/rodcontinuumsteklovpoincarestep.hh
@@ -19,6 +19,7 @@
template <class RodGridType, class ContinuumGridType>
class RodContinuumSteklovPoincareStep
+: public RodContinuumDDStep<RodGridType,ContinuumGridType>
{
static const int dim = ContinuumGridType::dimension;
@@ -51,7 +52,7 @@ public:
ContinuumLocalFiniteElement,
ContinuumLocalFiniteElement,
Dune::FieldMatrix<double,dim,dim> >* localAssembler)
- : complex_(complex),
+ : RodContinuumDDStep<RodGridType,ContinuumGridType>(complex),
preconditioner_(preconditioner),
alpha_(alpha),
richardsonDamping_(richardsonDamping),
@@ -83,7 +84,7 @@ public:
ContinuumLocalFiniteElement,
ContinuumLocalFiniteElement,
Dune::FieldMatrix<double,dim,dim> >*>& localAssembler)
- : complex_(complex),
+ : RodContinuumDDStep<RodGridType,ContinuumGridType>(complex),
preconditioner_(preconditioner),
alpha_(alpha),
richardsonDamping_(richardsonDamping),
@@ -176,24 +177,9 @@ protected:
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;
-
- /** \brief Add the content of one map to another, aborting rather than overwriting stuff
- */
- template <class X, class Y>
- static void insert(std::map<X,Y>& map1, const std::map<X,Y>& map2)
- {
- int oldSize = map1.size();
- map1.insert(map2.begin(), map2.end());
- assert(map1.size() == oldSize + map2.size());
- }
-
//////////////////////////////////////////////////////////////////
// Data members related to the coupled problem
//////////////////////////////////////////////////////////////////
- const RodContinuumComplex<RodGridType,ContinuumGridType>& complex_;
/** \brief Decides which preconditioner is used */
std::string preconditioner_;
@@ -234,9 +220,6 @@ protected:
typedef typename std::map<std::string, RodData>::iterator RodIterator;
-public:
- /** \todo Should be part of RodData, too */
- mutable std::map<std::string, RodConfigurationType> rodSubdomainSolutions_;
protected:
//////////////////////////////////////////////////////////////////
// Data members related to the continuum problems
@@ -267,10 +250,6 @@ protected:
typedef typename std::map<std::string, ContinuumData>::iterator ContinuumIterator;
-public:
- /** \todo Should be part of ContinuumData, too */
- mutable std::map<std::string, ContinuumConfigurationType> continuumSubdomainSolutions_;
-
};
@@ -295,23 +274,23 @@ mergeRodDirichletAndCouplingBoundaries()
// short-cut to avoid frequent map look-up
Dune::BitSetVector<6>& dirichletAndCouplingNodes = rods_[name].dirichletAndCouplingNodes_;
- dirichletAndCouplingNodes.resize(complex_.rodGrid(name)->size(1));
+ dirichletAndCouplingNodes.resize(this->complex_.rodGrid(name)->size(1));
// first copy the true Dirichlet boundary
- const LeafBoundaryPatch<RodGridType>& dirichletBoundary = complex_.rods_.find(name)->second.dirichletBoundary_;
+ const LeafBoundaryPatch<RodGridType>& dirichletBoundary = this->complex_.rods_.find(name)->second.dirichletBoundary_;
for (int i=0; i<dirichletAndCouplingNodes.size(); i++)
dirichletAndCouplingNodes[i] = dirichletBoundary.containsVertex(i);
// get the names of all the continua that we couple with
- std::set<std::string> continuumNames = continuaPerRod(name);
+ std::set<std::string> continuumNames = this->continuaPerRod(name);
for (std::set<std::string>::const_iterator cIt = continuumNames.begin();
cIt != continuumNames.end();
++cIt) {
const LeafBoundaryPatch<RodGridType>& rodInterfaceBoundary
- = complex_.coupling(std::make_pair(name,*cIt)).rodInterfaceBoundary_;
+ = this->complex_.coupling(std::make_pair(name,*cIt)).rodInterfaceBoundary_;
/** \todo Use the BoundaryPatch iterator here, for increased efficiency */
for (int i=0; i<dirichletAndCouplingNodes.size(); i++) {
@@ -346,23 +325,23 @@ mergeContinuumDirichletAndCouplingBoundaries()
// short-cut to avoid frequent map look-up
Dune::BitSetVector<dim>& dirichletAndCouplingNodes = continua_[name].dirichletAndCouplingNodes_;
- dirichletAndCouplingNodes.resize(complex_.continuumGrid(name)->size(dim));
+ dirichletAndCouplingNodes.resize(this->complex_.continuumGrid(name)->size(dim));
// first copy the true Dirichlet boundary
- const LeafBoundaryPatch<ContinuumGridType>& dirichletBoundary = complex_.continua_.find(name)->second.dirichletBoundary_;
+ const LeafBoundaryPatch<ContinuumGridType>& dirichletBoundary = this->complex_.continua_.find(name)->second.dirichletBoundary_;
for (int i=0; i<dirichletAndCouplingNodes.size(); i++)
dirichletAndCouplingNodes[i] = dirichletBoundary.containsVertex(i);
// get the names of all the rods that we couple with
- std::set<std::string> rodNames = rodsPerContinuum(name);
+ std::set<std::string> rodNames = this->rodsPerContinuum(name);
for (std::set<std::string>::const_iterator rIt = rodNames.begin();
rIt != rodNames.end();
++rIt) {
const LeafBoundaryPatch<ContinuumGridType>& continuumInterfaceBoundary
- = complex_.coupling(std::make_pair(*rIt,name)).continuumInterfaceBoundary_;
+ = this->complex_.coupling(std::make_pair(*rIt,name)).continuumInterfaceBoundary_;
/** \todo Use the BoundaryPatch iterator here, for increased efficiency */
for (int i=0; i<dirichletAndCouplingNodes.size(); i++) {
@@ -378,48 +357,20 @@ mergeContinuumDirichletAndCouplingBoundaries()
}
-template <class RodGridType, class ContinuumGridType>
-std::set<std::string> RodContinuumSteklovPoincareStep<RodGridType,ContinuumGridType>::
-rodsPerContinuum(const std::string& name) const
-{
- std::set<std::string> result;
-
- for (typename RodContinuumComplex<RodGridType,ContinuumGridType>::ConstCouplingIterator it = complex_.couplings_.begin();
- it!=complex_.couplings_.end(); ++it)
- if (it->first.second == name)
- result.insert(it->first.first);
-
- return result;
-}
-
-template <class RodGridType, class ContinuumGridType>
-std::set<std::string> RodContinuumSteklovPoincareStep<RodGridType,ContinuumGridType>::
-continuaPerRod(const std::string& name) const
-{
- std::set<std::string> result;
-
- for (typename RodContinuumComplex<RodGridType,ContinuumGridType>::ConstCouplingIterator it = complex_.couplings_.begin();
- it!=complex_.couplings_.end(); ++it)
- if (it->first.first == name)
- result.insert(it->first.second);
-
- return result;
-}
-
template <class RodGridType, class ContinuumGridType>
std::map<std::pair<std::string,std::string>,RigidBodyMotion<3>::TangentVector> RodContinuumSteklovPoincareStep<RodGridType,ContinuumGridType>::
rodDirichletToNeumannMap(const std::string& rodName,
const std::map<std::pair<std::string,std::string>,RigidBodyMotion<3> >& lambda) const
{
// container for the subdomain solution
- RodConfigurationType& rodX = rodSubdomainSolutions_[rodName];
- rodX.resize(complex_.rodGrid(rodName)->size(1));
+ RodConfigurationType& rodX = this->rodSubdomainSolutions_[rodName];
+ rodX.resize(this->complex_.rodGrid(rodName)->size(1));
///////////////////////////////////////////////////////////
// Set the complete set of Dirichlet values
///////////////////////////////////////////////////////////
- const LeafBoundaryPatch<RodGridType>& dirichletBoundary = complex_.rod(rodName).dirichletBoundary_;
- const RodConfigurationType& dirichletValues = complex_.rod(rodName).dirichletValues_;
+ const LeafBoundaryPatch<RodGridType>& dirichletBoundary = this->complex_.rod(rodName).dirichletBoundary_;
+ const RodConfigurationType& dirichletValues = this->complex_.rod(rodName).dirichletValues_;
for (size_t i=0; i<rodX.size(); i++)
if (dirichletBoundary.containsVertex(i))
@@ -434,7 +385,7 @@ rodDirichletToNeumannMap(const std::string& rodName,
continue;
// Use \lambda as a Dirichlet value for the rod
- const LeafBoundaryPatch<RodGridType>& interfaceBoundary = complex_.coupling(couplingName).rodInterfaceBoundary_;
+ const LeafBoundaryPatch<RodGridType>& interfaceBoundary = this->complex_.coupling(couplingName).rodInterfaceBoundary_;
/** \todo Use the BoundaryPatch iterator, which will be a lot faster
* once we use EntitySeed for its implementation
@@ -452,7 +403,7 @@ rodDirichletToNeumannMap(const std::string& rodName,
////////////////////////////////////////////////////////////////////////////////
// Set initial iterate by interpolating between the Dirichlet values
- RodFactory<typename RodGridType::LeafGridView> rodFactory(complex_.rodGrid(rodName)->leafView());
+ RodFactory<typename RodGridType::LeafGridView> rodFactory(this->complex_.rodGrid(rodName)->leafView());
rodFactory.create(rodX);
rod(rodName).solver_->setInitialSolution(rodX);
@@ -472,7 +423,7 @@ rodDirichletToNeumannMap(const std::string& rodName,
if (couplingName.first != rodName)
continue;
- const LeafBoundaryPatch<RodGridType>& couplingBoundary = complex_.coupling(couplingName).rodInterfaceBoundary_;
+ const LeafBoundaryPatch<RodGridType>& couplingBoundary = this->complex_.coupling(couplingName).rodInterfaceBoundary_;
result[couplingName] = rod(rodName).assembler_->getResultantForce(couplingBoundary, rodX);
}
@@ -487,13 +438,13 @@ continuumDirichletToNeumannMap(const std::string& continuumName,
const std::map<std::pair<std::string,std::string>,RigidBodyMotion<3> >& lambda) const
{
// Set initial iterate
- VectorType& x3d = continuumSubdomainSolutions_[continuumName];
- x3d.resize(complex_.continuumGrid(continuumName)->size(dim));
+ VectorType& x3d = this->continuumSubdomainSolutions_[continuumName];
+ x3d.resize(this->complex_.continuumGrid(continuumName)->size(dim));
x3d = 0;
// Copy the true Dirichlet values into it
- const LeafBoundaryPatch<ContinuumGridType>& dirichletBoundary = complex_.continuum(continuumName).dirichletBoundary_;
- const VectorType& dirichletValues = complex_.continuum(continuumName).dirichletValues_;
+ const LeafBoundaryPatch<ContinuumGridType>& dirichletBoundary = this->complex_.continuum(continuumName).dirichletBoundary_;
+ const VectorType& dirichletValues = this->complex_.continuum(continuumName).dirichletValues_;
for (size_t i=0; i<x3d.size(); i++)
if (dirichletBoundary.containsVertex(i))
@@ -508,8 +459,8 @@ continuumDirichletToNeumannMap(const std::string& continuumName,
continue;
// 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_;
+ const LeafBoundaryPatch<ContinuumGridType>& interfaceBoundary = this->complex_.coupling(couplingName).continuumInterfaceBoundary_;
+ const RigidBodyMotion<dim>& referenceInterface = this->complex_.coupling(couplingName).referenceInterface_;
setRotation(interfaceBoundary, x3d, referenceInterface, it->second);
@@ -553,7 +504,7 @@ continuumDirichletToNeumannMap(const std::string& continuumName,
if (couplingName.second != continuumName)
continue;
- const LeafBoundaryPatch<ContinuumGridType>& interfaceBoundary = complex_.coupling(couplingName).continuumInterfaceBoundary_;
+ const LeafBoundaryPatch<ContinuumGridType>& interfaceBoundary = this->complex_.coupling(couplingName).continuumInterfaceBoundary_;
VectorType neumannForces(residual.size());
neumannForces = 0;
@@ -561,7 +512,7 @@ continuumDirichletToNeumannMap(const std::string& continuumName,
weakToStrongBoundaryStress(interfaceBoundary, residual, neumannForces);
/** \todo Is referenceInterface.r the correct center of rotation? */
- const RigidBodyMotion<dim>& referenceInterface = complex_.coupling(couplingName).referenceInterface_;
+ const RigidBodyMotion<dim>& referenceInterface = this->complex_.coupling(couplingName).referenceInterface_;
computeTotalForceAndTorque(interfaceBoundary,
neumannForces,
@@ -597,7 +548,7 @@ linearizedRodNeumannToDirichletMap(const std::string& rodName,
// Assemble the linearized rod problem
////////////////////////////////////////////////////
- const LeafBoundaryPatch<RodGridType>& dirichletBoundary = complex_.rod(rodName).dirichletBoundary_;
+ const LeafBoundaryPatch<RodGridType>& dirichletBoundary = this->complex_.rod(rodName).dirichletBoundary_;
typedef Dune::BCRSMatrix<Dune::FieldMatrix<double,6,6> > MatrixType;
GeodesicFEAssembler<typename RodGridType::LeafGridView, RigidBodyMotion<3> > assembler(dirichletBoundary.gridView(),
@@ -657,7 +608,7 @@ linearizedRodNeumannToDirichletMap(const std::string& rodName,
}
// Use 'forceTorque' as a Neumann value for the rod
- const LeafBoundaryPatch<RodGridType>& interfaceBoundary = complex_.coupling(couplingName).rodInterfaceBoundary_;
+ const LeafBoundaryPatch<RodGridType>& interfaceBoundary = this->complex_.coupling(couplingName).rodInterfaceBoundary_;
const typename RodGridType::LeafGridView::IndexSet& indexSet = interfaceBoundary.gridView().indexSet();
@@ -728,7 +679,7 @@ linearizedRodNeumannToDirichletMap(const std::string& rodName,
if (couplingName.first != rodName)
continue;
- const LeafBoundaryPatch<RodGridType>& interfaceBoundary = complex_.coupling(couplingName).rodInterfaceBoundary_;
+ const LeafBoundaryPatch<RodGridType>& interfaceBoundary = this->complex_.coupling(couplingName).rodInterfaceBoundary_;
const typename RodGridType::LeafGridView::IndexSet& indexSet = interfaceBoundary.gridView().indexSet();
@@ -765,7 +716,7 @@ linearizedContinuumNeumannToDirichletMap(const std::string& continuumName,
* i.e. the linearization at zero
*/
typedef P1NodalBasis<typename ContinuumGridType::LeafGridView,double> P1Basis;
- const LeafBoundaryPatch<ContinuumGridType>& dirichletBoundary = complex_.continuum(continuumName).dirichletBoundary_;
+ const LeafBoundaryPatch<ContinuumGridType>& dirichletBoundary = this->complex_.continuum(continuumName).dirichletBoundary_;
P1Basis basis(dirichletBoundary.gridView());
OperatorAssembler<P1Basis,P1Basis> assembler(basis, basis);
@@ -779,7 +730,7 @@ linearizedContinuumNeumannToDirichletMap(const std::string& continuumName,
// The weak form of the volume and Neumann data
/** \brief Not implemented yet! */
- VectorType rhs(complex_.continuumGrid(continuumName)->size(dim));
+ VectorType rhs(this->complex_.continuumGrid(continuumName)->size(dim));
rhs = 0;
typedef typename std::map<std::pair<std::string,std::string>,RigidBodyMotion<3>::TangentVector>::const_iterator ForceIterator;
@@ -792,7 +743,7 @@ linearizedContinuumNeumannToDirichletMap(const std::string& continuumName,
continue;
// Use 'forceTorque' as a Neumann value for the rod
- const LeafBoundaryPatch<ContinuumGridType>& interfaceBoundary = complex_.coupling(couplingName).continuumInterfaceBoundary_;
+ const LeafBoundaryPatch<ContinuumGridType>& interfaceBoundary = this->complex_.coupling(couplingName).continuumInterfaceBoundary_;
//
VectorType localNeumannValues;
@@ -842,7 +793,7 @@ linearizedContinuumNeumannToDirichletMap(const std::string& continuumName,
continue;
// Use 'forceTorque' as a Neumann value for the rod
- const LeafBoundaryPatch<ContinuumGridType>& interfaceBoundary = complex_.coupling(couplingName).continuumInterfaceBoundary_;
+ const LeafBoundaryPatch<ContinuumGridType>& interfaceBoundary = this->complex_.coupling(couplingName).continuumInterfaceBoundary_;
RigidBodyMotion<3> averageInterface;
computeAverageInterface(interfaceBoundary, x, averageInterface);
@@ -884,7 +835,7 @@ iterate(std::map<std::pair<std::string,std::string>, RigidBodyMotion<3> >& lambd
std::map<std::pair<std::string,std::string>, RigidBodyMotion<3>::TangentVector> forceTorque = rodDirichletToNeumannMap(rodName, lambda);
- insert(rodForceTorque, forceTorque);
+ this->insert(rodForceTorque, forceTorque);
}
@@ -907,7 +858,7 @@ iterate(std::map<std::pair<std::string,std::string>, RigidBodyMotion<3> >& lambd
std::map<std::pair<std::string,std::string>, RigidBodyMotion<3>::TangentVector> forceTorque
= continuumDirichletToNeumannMap(continuumName, lambda);
- insert(continuumForceTorque, forceTorque);
+ this->insert(continuumForceTorque, forceTorque);
}
@@ -953,11 +904,11 @@ iterate(std::map<std::pair<std::string,std::string>, RigidBodyMotion<3> >& lambd
std::map<std::pair<std::string,std::string>, RigidBodyMotion<3>::TangentVector> continuumInterfaceCorrection
= linearizedContinuumNeumannToDirichletMap(continuumName,
- continuumSubdomainSolutions_[continuumName],
+ this->continuumSubdomainSolutions_[continuumName],
residualForceTorque,
lambda);
- insert(continuumCorrection, continuumInterfaceCorrection);
+ this->insert(continuumCorrection, continuumInterfaceCorrection);
}
@@ -982,11 +933,11 @@ iterate(std::map<std::pair<std::string,std::string>, RigidBodyMotion<3> >& lambd
std::map<std::pair<std::string,std::string>, RigidBodyMotion<3>::TangentVector> rodInterfaceCorrection
= linearizedRodNeumannToDirichletMap(rodName,
- rodSubdomainSolutions_[rodName],
+ this->rodSubdomainSolutions_[rodName],
residualForceTorque,
lambda);
- insert(rodCorrection, rodInterfaceCorrection);
+ this->insert(rodCorrection, rodInterfaceCorrection);
}