From 5862d8155b224a1e6b97dbce9de4ef09d6ac946f Mon Sep 17 00:00:00 2001
From: Oliver Sander <sander@igpm.rwth-aachen.de>
Date: Tue, 4 Oct 2011 16:11:16 +0000
Subject: [PATCH] Make rod code compile again now that the general assembler
 framework wants a global function space basis.

So far the P1 basis is hardwired for simplicity.

[[Imported from SVN: r7863]]
---
 .../coupling/rodcontinuumsteklovpoincarestep.hh   |  2 +-
 dune/gfe/rodassembler.cc                          | 15 +++++++++------
 dune/gfe/rodassembler.hh                          |  4 ++--
 dune/gfe/rodlocalstiffness.hh                     | 12 +++++++++++-
 4 files changed, 23 insertions(+), 10 deletions(-)

diff --git a/dune/gfe/coupling/rodcontinuumsteklovpoincarestep.hh b/dune/gfe/coupling/rodcontinuumsteklovpoincarestep.hh
index 66b8f85e..ba0830b5 100644
--- a/dune/gfe/coupling/rodcontinuumsteklovpoincarestep.hh
+++ b/dune/gfe/coupling/rodcontinuumsteklovpoincarestep.hh
@@ -559,7 +559,7 @@ linearizedRodNeumannToDirichletMap(const std::string& rodName,
     const RodLeafBoundaryPatch& dirichletBoundary = this->complex_.rod(rodName).dirichletBoundary_;
 
     typedef Dune::BCRSMatrix<Dune::FieldMatrix<double,6,6> > MatrixType;
-    GeodesicFEAssembler<typename RodGridType::LeafGridView, RigidBodyMotion<3> > assembler(dirichletBoundary.gridView(),
+    GeodesicFEAssembler<P1NodalBasis<typename RodGridType::LeafGridView>, RigidBodyMotion<3> > assembler(dirichletBoundary.gridView(),
                                                                     rod(rodName).localStiffness_);
     MatrixType stiffnessMatrix;
     assembler.assembleMatrix(currentIterate, 
diff --git a/dune/gfe/rodassembler.cc b/dune/gfe/rodassembler.cc
index dc356b2e..6ab87a73 100644
--- a/dune/gfe/rodassembler.cc
+++ b/dune/gfe/rodassembler.cc
@@ -18,7 +18,7 @@ assembleGradient(const std::vector<RigidBodyMotion<3> >& sol,
 {
     using namespace Dune;
 
-    const typename GridView::Traits::IndexSet& indexSet = this->gridView_.indexSet();
+    const typename GridView::Traits::IndexSet& indexSet = this->basis_.getGridView().indexSet();
 
     if (sol.size()!=indexSet.size(gridDim))
         DUNE_THROW(Exception, "Solution vector doesn't match the grid!");
@@ -26,8 +26,8 @@ assembleGradient(const std::vector<RigidBodyMotion<3> >& sol,
     grad.resize(sol.size());
     grad = 0;
 
-    ElementIterator it    = this->gridView_.template begin<0>();
-    ElementIterator endIt = this->gridView_.template end<0>();
+    ElementIterator it    = this->basis_.getGridView().template begin<0>();
+    ElementIterator endIt = this->basis_.getGridView().template end<0>();
 
     // Loop over all elements
     for (; it!=endIt; ++it) {
@@ -36,7 +36,7 @@ assembleGradient(const std::vector<RigidBodyMotion<3> >& sol,
         static const int nDofs = 2;
 
         // Extract local solution
-        Dune::array<RigidBodyMotion<3>, nDofs> localSolution;
+        std::vector<RigidBodyMotion<3> > localSolution(nDofs);
         
         for (int i=0; i<nDofs; i++)
             localSolution[i] = sol[indexSet.subIndex(*it,i,gridDim)];
@@ -44,7 +44,10 @@ assembleGradient(const std::vector<RigidBodyMotion<3> >& sol,
         // Assemble local gradient
         std::vector<FieldVector<double,blocksize> > localGradient(nDofs);
 
-        this->localStiffness_->assembleGradient(*it, localSolution, localGradient);
+        this->localStiffness_->assembleGradient(*it, 
+                                                this->basis_.getLocalFiniteElement(*it),
+                                                localSolution, 
+                                                localGradient);
 
         // Add to global gradient
         for (int i=0; i<nDofs; i++)
@@ -150,7 +153,7 @@ getResultantForce(const BoundaryPatch<PatchGridView>& boundary,
     //    if (gridView_ != &boundary.gridView())
     //        DUNE_THROW(Dune::Exception, "The boundary patch has to match the grid view of the assembler!");
 
-    const typename GridView::Traits::IndexSet& indexSet = this->gridView_.indexSet();
+    const typename GridView::Traits::IndexSet& indexSet = this->basis_.getGridView().indexSet();
 
     if (sol.size()!=indexSet.size(gridDim))
         DUNE_THROW(Exception, "Solution vector doesn't match the grid!");
diff --git a/dune/gfe/rodassembler.hh b/dune/gfe/rodassembler.hh
index cd49bd42..77fac7d8 100644
--- a/dune/gfe/rodassembler.hh
+++ b/dune/gfe/rodassembler.hh
@@ -24,7 +24,7 @@ class RodAssembler
 /** \brief The FEM operator for an extensible, shearable rod in 3d
  */
 template <class GridView>
-class RodAssembler<GridView,3> : public GeodesicFEAssembler<GridView, RigidBodyMotion<3> >
+class RodAssembler<GridView,3> : public GeodesicFEAssembler<P1NodalBasis<GridView>, RigidBodyMotion<3> >
 {
         
     //typedef typename GridType::template Codim<0>::Entity EntityType;
@@ -46,7 +46,7 @@ public:
         //! ???
     RodAssembler(const GridView &gridView,
                  RodLocalStiffness<GridView,double>* localStiffness) 
-        : GeodesicFEAssembler<GridView, RigidBodyMotion<3> >(gridView,localStiffness)
+        : GeodesicFEAssembler<P1NodalBasis<GridView>, RigidBodyMotion<3> >(gridView,localStiffness)
         { 
             std::vector<RigidBodyMotion<3> > referenceConfiguration(gridView.size(gridDim));
 
diff --git a/dune/gfe/rodlocalstiffness.hh b/dune/gfe/rodlocalstiffness.hh
index ac64fff6..89296c02 100644
--- a/dune/gfe/rodlocalstiffness.hh
+++ b/dune/gfe/rodlocalstiffness.hh
@@ -11,7 +11,7 @@
 
 template<class GridView, class RT>
 class RodLocalStiffness 
-    : public LocalGeodesicFEStiffness<GridView,RigidBodyMotion<3> >
+    : public LocalGeodesicFEStiffness<GridView, typename P1NodalBasis<GridView>::LocalFiniteElement, RigidBodyMotion<3> >
 {
     typedef RigidBodyMotion<3> TargetSpace;
 
@@ -95,9 +95,19 @@ public:
         referenceConfiguration_ = referenceConfiguration;
     }
     
+    /** \brief Local element energy for a P1 element */
     virtual RT energy (const Entity& e,
                        const Dune::array<RigidBodyMotion<3>, dim+1>& localSolution) const;
 
+    virtual RT energy (const Entity& e,
+                       const typename P1NodalBasis<GridView>::LocalFiniteElement& localFiniteElement,
+                       const std::vector<RigidBodyMotion<3> >& localSolution) const
+    {
+        assert(localSolution.size()==2);
+        Dune::array<RigidBodyMotion<3>, 2> localSolutionArray = {localSolution[0], localSolution[1]};
+        return energy(e,localSolutionArray);
+    }
+
     /** \brief Assemble the element gradient of the energy functional */
     void assembleGradient(const Entity& element,
                           const std::vector<RigidBodyMotion<3> >& solution,
-- 
GitLab