dune/gfe/mixedriemanniantrsolver.cc

@@ -208,7 +208,7 @@ setup(const GridType& grid,
                                          FufemBasis0>(pkToP1TransferMatrix,p1Basis,*basis0_);
 
         mmgStep0->mgTransfer_.back() = std::make_shared<TruncatedCompressedMGTransfer<CorrectionType0>>();
-        Dune::shared_ptr<TransferOperatorType> topTransferOperator = Dune::make_shared<TransferOperatorType>(pkToP1TransferMatrix);
+        std::shared_ptr<TransferOperatorType> topTransferOperator = std::make_shared<TransferOperatorType>(pkToP1TransferMatrix);
         std::dynamic_pointer_cast<TruncatedCompressedMGTransfer<CorrectionType0>>(mmgStep0->mgTransfer_.back())->setMatrix(topTransferOperator);
 
         for (size_t i=0; i<mmgStep0->mgTransfer_.size()-1; i++){

dune/gfe/riemanniantrsolver.hh

@@ -137,7 +137,7 @@ public:
     void setIgnoreNodes(const Dune::BitSetVector<blocksize>& ignoreNodes)
     {
         ignoreNodes_ = &ignoreNodes;
-        Dune::shared_ptr<LoopSolver<CorrectionType> > loopSolver = std::dynamic_pointer_cast<LoopSolver<CorrectionType> >(innerSolver_);
+        std::shared_ptr<LoopSolver<CorrectionType> > loopSolver = std::dynamic_pointer_cast<LoopSolver<CorrectionType> >(innerSolver_);
         assert(loopSolver);
         loopSolver->iterationStep_->ignoreNodes_ = ignoreNodes_;
     }

dune/gfe/rodassembler.hh

@@ -41,12 +41,9 @@ class RodAssembler<Basis,3> : public GeodesicFEAssembler<Basis, RigidBodyMotion<
 public:
         //! ???
     RodAssembler(const Basis& basis,
-                 RodLocalStiffness<GridView,double>* localStiffness)
-        : GeodesicFEAssembler<Basis, RigidBodyMotion<double,3> >(basis,nullptr)
-        , rodEnergy_(localStiffness)
+                 LocalGeodesicFEStiffness<Basis, RigidBodyMotion<double,3> >* localStiffness)
+    : GeodesicFEAssembler<Basis, RigidBodyMotion<double,3> >(basis,localStiffness)
         {
-            this->localStiffness_ = new LocalGeodesicFEFDStiffness<Basis,RigidBodyMotion<double,3>, double>(localStiffness);
-
             std::vector<RigidBodyMotion<double,3> > referenceConfiguration(basis.size());
 
     for (const auto vertex : Dune::vertices(basis.gridView()))
@@ -59,12 +56,19 @@ public:
                 referenceConfiguration[idx].q = Rotation<double,3>::identity();
             }
 
-    rodEnergy_->setReferenceConfiguration(referenceConfiguration);
+    rodEnergy()->setReferenceConfiguration(referenceConfiguration);
         }
 
+    auto rodEnergy()
+    {
+      // TODO: Does not work for other stiffness implementations
+      auto localFDStiffness = dynamic_cast<LocalGeodesicFEFDStiffness<Basis, RigidBodyMotion<double,3> >*>(this->localStiffness_);
+      return const_cast<RodLocalStiffness<GridView,double>*>(dynamic_cast<const RodLocalStiffness<GridView,double>*>(localFDStiffness->localEnergy_));
+    }
+
         std::vector<RigidBodyMotion<double,3> > getRefConfig()
     {
-      return rodEnergy_->referenceConfiguration_;
+      return rodEnergy()->referenceConfiguration_;
         }
 
   virtual void assembleGradient(const std::vector<RigidBodyMotion<double,3> >& sol,
@@ -82,8 +86,6 @@ public:
         template <class PatchGridView>
         Dune::FieldVector<double,6> getResultantForce(const BoundaryPatch<PatchGridView>& boundary,
                                                       const std::vector<RigidBodyMotion<double,3> >& sol) const;
-
-    RodLocalStiffness<GridView,double>* rodEnergy_;
     }; // end class
 
 

dune/gfe/svd.hh

@@ -41,9 +41,8 @@ void svdcmp(Dune::FieldMatrix<T,m,n>& a_, Dune::FieldVector<T,n>& w, Dune::Field
             a[i+1][j+1] = a_[i][j];
 
     int flag,i,its,j,jj,k,l,nm;
-    T anorm,c,f,g,h,s,scale,x,y,z,*rv1;
-    T rv1_c[n];  // 1 too large to accomodate fortran numbering
-    rv1 = rv1_c-1;
+    T anorm,c,f,g,h,s,scale,x,y,z;
+    T rv1[n+1];  // 1 too large to accomodate fortran numbering
 
     //Householder reduction to bidiagonal form.
     g=scale=anorm=0.0;

src/compute-disc-error.cc

@@ -671,7 +671,7 @@ int main (int argc, char *argv[]) try
 
   const int numLevels = parameterSet.get<int>("numLevels");
 
-  shared_ptr<GridType> grid, referenceGrid;
+  std::shared_ptr<GridType> grid, referenceGrid;
 
   FieldVector<double,dimworld> lower(0), upper(1);
 
@@ -706,8 +706,8 @@ int main (int argc, char *argv[]) try
   {
     std::string path                = parameterSet.get<std::string>("path");
     std::string gridFile            = parameterSet.get<std::string>("gridFile");
-    grid          = shared_ptr<GridType>(GmshReader<GridType>::read(path + "/" + gridFile));
-    referenceGrid = shared_ptr<GridType>(GmshReader<GridType>::read(path + "/" + gridFile));
+    grid          = std::shared_ptr<GridType>(GmshReader<GridType>::read(path + "/" + gridFile));
+    referenceGrid = std::shared_ptr<GridType>(GmshReader<GridType>::read(path + "/" + gridFile));
   }
 
   grid->globalRefine(numLevels-1);

src/cosserat-continuum.cc

@@ -177,7 +177,7 @@ int main (int argc, char *argv[]) try
     typedef UGGrid<dim> GridType;
 #endif
 
-    shared_ptr<GridType> grid;
+    std::shared_ptr<GridType> grid;
 
     FieldVector<double,dimworld> lower(0), upper(1);
 
@@ -201,7 +201,7 @@ int main (int argc, char *argv[]) try
         std::string suffix = gridFile.substr(dotPos, gridFile.length()-dotPos);
 
         if (suffix == ".msh")
-            grid = shared_ptr<GridType>(GmshReader<GridType>::read(path + "/" + gridFile));
+            grid = std::shared_ptr<GridType>(GmshReader<GridType>::read(path + "/" + gridFile));
         else if (suffix == ".vtu" or suffix == ".vtp")
             grid = VTKReader<GridType>::read(path + "/" + gridFile);
     }
@@ -266,12 +266,10 @@ int main (int argc, char *argv[]) try
 
     for (auto&& vertex : vertices(gridView))
     {
-        bool isDirichlet;
-        pythonDirichletVertices.evaluate(vertex.geometry().corner(0), isDirichlet);
+        bool isDirichlet = pythonDirichletVertices(vertex.geometry().corner(0));
         dirichletVertices[indexSet.index(vertex)] = isDirichlet;
 
-        bool isNeumann;
-        pythonNeumannVertices.evaluate(vertex.geometry().corner(0), isNeumann);
+        bool isNeumann = pythonNeumannVertices(vertex.geometry().corner(0));
         neumannVertices[indexSet.index(vertex)] = isNeumann;
     }
 
@@ -429,14 +427,14 @@ int main (int argc, char *argv[]) try
     // ////////////////////////////////////////////////////////////
 
     const ParameterTree& materialParameters = parameterSet.sub("materialParameters");
-    shared_ptr<NeumannFunction> neumannFunction;
+    std::shared_ptr<NeumannFunction> neumannFunction;
     if (parameterSet.hasKey("neumannValues"))
-        neumannFunction = make_shared<NeumannFunction>(parameterSet.get<FieldVector<double,3> >("neumannValues"),
+        neumannFunction = std::make_shared<NeumannFunction>(parameterSet.get<FieldVector<double,3> >("neumannValues"),
                                                        homotopyParameter);
 
-    shared_ptr<VolumeLoad> volumeLoad;
+    std::shared_ptr<VolumeLoad> volumeLoad;
     if (parameterSet.hasKey("volumeLoad"))
-        volumeLoad = make_shared<VolumeLoad>(parameterSet.get<FieldVector<double,3> >("volumeLoad"),
+        volumeLoad = std::make_shared<VolumeLoad>(parameterSet.get<FieldVector<double,3> >("volumeLoad"),
                                                                                           homotopyParameter);
 
     if (mpiHelper.rank() == 0) {

src/film-on-substrate.cc

@@ -348,7 +348,7 @@ int main (int argc, char *argv[]) try
     std::shared_ptr<NeumannFunction> neumannFunction;
     if (parameterSet.hasKey("neumannValues"))
     {
-      neumannFunction = make_shared<NeumannFunction>(parameterSet.get<FieldVector<double,dim> >("neumannValues"),
+      neumannFunction = std::make_shared<NeumannFunction>(parameterSet.get<FieldVector<double,dim> >("neumannValues"),
                                                      homotopyParameter);
 
       std::cout << "Neumann values: " << parameterSet.get<FieldVector<double,dim> >("neumannValues") << std::endl;

src/gradient-flow.cc

@@ -164,8 +164,7 @@ int main (int argc, char *argv[]) try
 
   for (auto&& vertex : vertices(grid->leafGridView()))
   {
-    bool isDirichlet;
-    pythonDirichletVertices.evaluate(vertex.geometry().corner(0), isDirichlet);
+    bool isDirichlet = pythonDirichletVertices(vertex.geometry().corner(0));
     dirichletVertices[indexSet.index(vertex)] = isDirichlet;
   }
 

src/harmonicmaps.cc

@@ -172,7 +172,7 @@ int main (int argc, char *argv[])
     typedef std::conditional<dim==1,OneDGrid,UGGrid<dim> >::type GridType;
 #endif
 
-    shared_ptr<GridType> grid;
+    std::shared_ptr<GridType> grid;
     FieldVector<double,dimworld> lower(0), upper(1);
     std::array<unsigned int,dim> elements;
 
@@ -195,7 +195,7 @@ int main (int argc, char *argv[])
         std::string path                = parameterSet.get<std::string>("path");
         std::string gridFile            = parameterSet.get<std::string>("gridFile");
 
-        grid = shared_ptr<GridType>(GmshReader<GridType>::read(path + "/" + gridFile));
+        grid = std::shared_ptr<GridType>(GmshReader<GridType>::read(path + "/" + gridFile));
     }
 
     grid->globalRefine(numLevels-1);
@@ -224,9 +224,7 @@ int main (int argc, char *argv[])
 
     for (auto&& vertex : vertices(gridView))
     {
-      //bool isDirichlet;
       bool isDirichlet = pythonDirichletVertices(vertex.geometry().corner(0));
-      pythonDirichletVertices.evaluate(vertex.geometry().corner(0), isDirichlet);
       dirichletVertices[indexSet.index(vertex)] = isDirichlet;
     }
 

src/rod3d.cc

@@ -128,7 +128,9 @@ int main (int argc, char *argv[]) try
     RodLocalStiffness<GridView,double> localStiffness(gridView,
                                                       A, J1, J2, E, nu);
 
-    RodAssembler<FEBasis,3> rodAssembler(gridView, &localStiffness);
+    LocalGeodesicFEFDStiffness<FEBasis,RigidBodyMotion<double,3> > localFDStiffness(&localStiffness);
+
+    RodAssembler<FEBasis,3> rodAssembler(gridView, &localFDStiffness);
 
     RiemannianTrustRegionSolver<FEBasis,RigidBodyMotion<double,3> > rodSolver;
 

test/adolctest.cc

@@ -35,7 +35,7 @@ typedef double FDType;
 
 #include <dune/functions/functionspacebases/lagrangebasis.hh>
 #include <dune/functions/functionspacebases/interpolate.hh>
-
+#include <dune/functions/functionspacebases/powerbasis.hh>
 
 #include <dune/gfe/rigidbodymotion.hh>
 #include <dune/gfe/localgeodesicfestiffness.hh>
@@ -464,7 +464,15 @@ int main (int argc, char *argv[]) try
     auto identity = [](const FieldVector<double,2>& x) -> FieldVector<double,3> { return {x[0], x[1], 0};};
 
     std::vector<FieldVector<double,3> > v;
-    Functions::interpolate(feBasis, v, identity);
+    using namespace Functions::BasisFactory;
+
+      auto powerBasis = makeBasis(
+        gridView,
+        power<3>(
+          lagrange<1>(),
+          blockedInterleaved()
+      ));
+    Functions::interpolate(powerBasis, v, identity);
 
     for (size_t i=0; i<x.size(); i++)
       x[i].r = v[i];

test/frameinvariancetest.cc

@@ -74,10 +74,12 @@ int main (int argc, char *argv[]) try
         rotatedX[i].q = rotation.mult(x[i].q);
     }
 
-    RodLocalStiffness<GridView,double> localStiffness(gridView,
-                                                      1,1,1,1e6,0.3);
+    RodLocalStiffness<GridView,double> localRodFirstOrderModel(gridView,
+                                                               1,1,1,1e6,0.3);
 
-    RodAssembler<FEBasis,3> assembler(feBasis, &localStiffness);
+    LocalGeodesicFEFDStiffness<FEBasis,RigidBodyMotion<double,3> > localFDStiffness(&localRodFirstOrderModel);
+
+    RodAssembler<FEBasis,3> assembler(feBasis, &localFDStiffness);
 
     if (std::abs(assembler.computeEnergy(x) - assembler.computeEnergy(rotatedX)) > 1e-6)
         DUNE_THROW(Dune::Exception, "Rod energy not invariant under rigid body motions!");

test/localgeodesicfefunctiontest.cc

@@ -76,8 +76,7 @@ void testPermutationInvariance(const std::vector<TargetSpace>& corners)
     PQkLocalFiniteElementCache<double,double,domainDim,1> feCache;
     typedef typename PQkLocalFiniteElementCache<double,double,domainDim,1>::FiniteElementType LocalFiniteElement;
     
-    GeometryType simplex;
-    simplex.makeSimplex(domainDim);
+    GeometryType simplex = GeometryTypes::simplex(domainDim);
 
     //
     std::vector<TargetSpace> cornersRotated1(domainDim+1);

test/localprojectedfefunctiontest.cc

@@ -118,8 +118,7 @@ void testPermutationInvariance(const std::vector<TargetSpace>& corners)
     PQkLocalFiniteElementCache<double,double,domainDim,1> feCache;
     typedef typename PQkLocalFiniteElementCache<double,double,domainDim,1>::FiniteElementType LocalFiniteElement;
 
-    GeometryType simplex;
-    simplex.makeSimplex(domainDim);
+    GeometryType simplex = GeometryTypes::simplex(domainDim);
 
     //
     std::vector<TargetSpace> cornersRotated1(domainDim+1);

test/rodassemblertest.cc

@@ -551,8 +551,9 @@ int main (int argc, char *argv[]) try
     RodLocalStiffness<GridView,double> localStiffness(gridView,
                                                                     0.01, 0.0001, 0.0001, 2.5e5, 0.3);
 
+    LocalGeodesicFEFDStiffness<Basis,RigidBodyMotion<double,3> > localFDStiffness(&localStiffness);
 
-    RodAssembler<Basis,3> rodAssembler(basis, &localStiffness);
+    RodAssembler<Basis,3> rodAssembler(basis, &localFDStiffness);
 
     std::cout << "Energy: " << rodAssembler.computeEnergy(x) << std::endl;