Rebase to master

In particular this means
* No more dune-fufem bases
* Functioning periodicity support
* Get rid of deprecated PythonFunction

src/quasiconvexity-test.cc

@@ -24,12 +24,11 @@
 #include <dune/functions/functionspacebases/interpolate.hh>
 #include <dune/functions/functionspacebases/lagrangebasis.hh>
 #include <dune/functions/gridfunctions/discreteglobalbasisfunction.hh>
+#include <dune/functions/functionspacebases/periodicbasis.hh>
 #include <dune/functions/functionspacebases/powerbasis.hh>
 
 #include <dune/fufem/boundarypatch.hh>
 #include <dune/fufem/functiontools/boundarydofs.hh>
-#include <dune/fufem/functionspacebases/dunefunctionsbasis.hh>
-#include <dune/fufem/functionspacebases/periodicbasis.hh>
 #include <dune/fufem/dunepython.hh>
 
 #include <dune/solvers/solvers/iterativesolver.hh>
@@ -38,16 +37,16 @@
 #include <dune/elasticity/common/trustregionsolver.hh>
 #include <dune/elasticity/assemblers/localadolcstiffness.hh>
 #include <dune/elasticity/assemblers/feassembler.hh>
-#include <dune/elasticity/materials/neffenergy.hh>
-#include <dune/elasticity/materials/coshenergy.hh>
-#include <dune/elasticity/materials/dacorognaenergy.hh>
-#include <dune/elasticity/materials/expacoshenergy.hh>
-#include <dune/elasticity/materials/fungenergy.hh>
+//#include <dune/elasticity/materials/neffenergy.hh>
+//#include <dune/elasticity/materials/coshenergy.hh>
+//#include <dune/elasticity/materials/dacorognaenergy.hh>
+//#include <dune/elasticity/materials/expacoshenergy.hh>
+//#include <dune/elasticity/materials/fungenergy.hh>
 #include <dune/elasticity/materials/localintegralenergy.hh>
-#include <dune/elasticity/materials/nefflogenergy.hh>
-#include <dune/elasticity/materials/nefflogenergy2.hh>
-#include <dune/elasticity/materials/neffreducedenergy.hh>
-#include <dune/elasticity/materials/neffw2energy.hh>
+//#include <dune/elasticity/materials/nefflogenergy.hh>
+//#include <dune/elasticity/materials/nefflogenergy2.hh>
+//#include <dune/elasticity/materials/neffreducedenergy.hh>
+//#include <dune/elasticity/materials/neffw2energy.hh>
 
 // grid dimension
 const int dim = 2;
@@ -233,7 +232,7 @@ void writeDisplacement(const Basis& basis, const Vector& periodicX, const FieldM
   // Add the underlying affine displacement.  For this we need
   // a representation in a non-periodic basis.
   BCRSMatrix<double> matrix;
-  assembleGlobalBasisTransferMatrix(matrix, basis, nonperiodicBasis);
+  assembleGlobalBasisTransferMatrix(matrix, basis, nonperiodicPowerBasis);
 
   Vector nonperiodicX(nonperiodicBasis.size());
   matrix.mv(periodicX, nonperiodicX);
@@ -397,25 +396,6 @@ int main (int argc, char *argv[]) try
   GridView gridView = grid->leafGridView();
 #endif
 
-  // FE basis spanning the FE space that we are working in
-  //using FEBasis = Functions::LagrangeBasis<GridView, order>;
-  static_assert(order==1, "Periodic grids are only implemented for order==1");
-  using FEBasis = Fufem::PeriodicBasis<GridView>;
-  FEBasis feBasis(gridView);
-
-#if 1
-  // This is needed for interpolation, but PeriodicBasis
-  // doesn't quite support being put into power bases yet.
-  using namespace Dune::Fufem::BasisFactory;
-  using namespace Dune::Functions::BasisFactory;
-  auto periodicPowerBasis = makeBasis(
-    gridView,
-    power<dim>(
-      periodic(),
-      blockedInterleaved()
-  ));
-#endif
-
   // Check whether two points are equal on R/Z x R/Z
   auto equivalent = [](const FieldVector<double,2>& x, const FieldVector<double,2>& y)
   {
@@ -440,51 +420,37 @@ int main (int argc, char *argv[]) try
   std::cout << "Grid has " << boundaryVertices.size() << " boundary vertices" << std::endl;
 
   bool periodic = parameterSet.get<bool>("periodic", false);
+
+  //using namespace Dune::Fufem::BasisFactory;
+  using namespace Dune::Functions::BasisFactory;
+  Functions::BasisFactory::PeriodicIndexSet periodicIndices;
   if (periodic)
   {
-#if 0
-    for (const auto& v1 : vertices(gridView))
-      for (const auto& v2 : vertices(gridView))
-        if (equivalent(v1.geometry().corner(0), v2.geometry().corner(0)))
-        {
-          //std::cout << "unifying " << v1.geometry().corner(0) << " with " << v2.geometry().corner(0) << std::endl;
-          feBasis.preBasis().unifyIndexPair({gridView.indexSet().index(v1)}, {gridView.indexSet().index(v2)});
-        }
-#else
+    if (order!=1)
+      DUNE_THROW(NotImplemented, "Periodicity detection only implemented for order==1");
+
     for (const auto& v1 : boundaryVertices)
       for (const auto& v2 : boundaryVertices)
         if (equivalent(v1.second, v2.second))
-        {
-          //std::cout << "unifying " << v1.geometry().corner(0) << " with " << v2.geometry().corner(0) << std::endl;
-          feBasis.preBasis().unifyIndexPair({v1.first}, {v2.first});
-        }
-#endif
+          periodicIndices.unifyIndexPair(v1.first, v2.first);
   }
 
-  feBasis.preBasis().initializeIndices();
+  // FE basis spanning the FE space that we are working in
+  auto feBasis = makeBasis(
+      gridView,
+      power<dim>(
+        Functions::BasisFactory::periodic(lagrange<order>(), periodicIndices),
+        blockedInterleaved()
+    ));
+  using FEBasis = decltype(feBasis);
 
   std::cout << "Basis has " << feBasis.size() << " degrees of freedom" << std::endl;
 
   // Get positions of the Lagrange nodes
   std::vector<FieldVector<double,dim> > lagrangeNodes(feBasis.size());
 
-#if 0  // Use this once periodicPowerBasis works
   const auto identityFct = [](const FieldVector<double,dim>& p) {return p;};
-  Functions::interpolate(periodicPowerBasis, lagrangeNodes, identityFct);
-#else
-  const auto identityX = [](const FieldVector<double,dim>& p) {return p[0];};
-  const auto identityY = [](const FieldVector<double,dim>& p) {return p[1];};
-
-  std::vector<double> lagrangeNodesX(feBasis.size());
-  std::vector<double> lagrangeNodesY(feBasis.size());
-  Functions::interpolate(feBasis, lagrangeNodesX, identityX);
-  Functions::interpolate(feBasis, lagrangeNodesY, identityY);
-  for (std::size_t i=0; i<lagrangeNodes.size(); i++)
-  {
-    lagrangeNodes[i][0] = lagrangeNodesX[i];
-    lagrangeNodes[i][1] = lagrangeNodesY[i];
-  }
-#endif
+  Functions::interpolate(feBasis, lagrangeNodes, identityFct);
 
   ///////////////////////////////////////////
   //   Set Dirichlet values
@@ -510,7 +476,8 @@ int main (int argc, char *argv[]) try
     dirichletVertices.setAll();
     BoundaryPatch<GridView> dirichletBoundary(gridView, dirichletVertices);
 
-    BitSetVector<1> dirichletNodes(feBasis.size(), false);
+    // TODO: Do we still need to distinguish between dirichletNodes and dirichletDofs?
+    BitSetVector<dim> dirichletNodes(feBasis.size(), false);
     constructBoundaryDofs(dirichletBoundary,feBasis,dirichletNodes);
 
     for (size_t i=0; i<feBasis.size(); i++)
@@ -561,10 +528,10 @@ int main (int argc, char *argv[]) try
     // The python class that implements the Dirichlet values
     Python::Module initialIterateModule = Python::import(parameterSet.get<std::string>("initialIteratePython"));
 
-    auto f = Python::Callable(initialIterateModule.get("initialIterate"));
-    PythonFunction<FieldVector<double,dim>, FieldVector<double,dim> > initialIteratePython(f);
+    // Method of that python class
+    auto pythonInitialIterate = Python::makeFunction<FieldVector<double,dim>(const FieldVector<double,dim>&)>(initialIterateModule.get("f"));
 
-    Dune::Functions::interpolate(powerBasis, x, initialIteratePython);
+    Functions::interpolate(powerBasis, x, pythonInitialIterate);
   }
 
   if (parameterSet.hasKey("initialIterateFile"))
@@ -658,15 +625,12 @@ int main (int argc, char *argv[]) try
 
     // Assembler using ADOL-C
     std::cout << "Selected energy is: " << parameterSet.get<std::string>("energy") << std::endl;
-    std::shared_ptr<Elasticity::LocalEnergy<GridView,
-                                     FEBasis::LocalView::Tree::FiniteElement,
-                                     std::vector<Dune::FieldVector<adouble, dim> > > > elasticEnergy;
+    std::shared_ptr<Elasticity::LocalEnergy<FEBasis::LocalView,adouble> > elasticEnergy;
 
     if (parameterSet.get<std::string>("energy") == "kpower")
     {
       auto density = std::make_shared<KPowerDensity<dim,adouble> >(F0, materialParameters);
-      elasticEnergy = std::make_shared<Elasticity::LocalIntegralEnergy<GridView,
-                                                                       FEBasis::LocalView::Tree::FiniteElement,
+      elasticEnergy = std::make_shared<Elasticity::LocalIntegralEnergy<FEBasis::LocalView,
                                                                        adouble> >(density);
     }
 
@@ -704,8 +668,7 @@ int main (int argc, char *argv[]) try
     if (parameterSet.get<std::string>("energy") == "expacosh")
     {
       auto density = std::make_shared<ExpACosHDensity<dim,adouble> >(F0,materialParameters);
-      elasticEnergy = std::make_shared<Elasticity::LocalIntegralEnergy<GridView,
-                                                                       FEBasis::LocalView::Tree::FiniteElement,
+      elasticEnergy = std::make_shared<Elasticity::LocalIntegralEnergy<FEBasis::LocalView,
                                                                        adouble> >(density);
     }
 #endif
@@ -717,16 +680,14 @@ int main (int argc, char *argv[]) try
     if (parameterSet.get<std::string>("energy") == "magic")
     {
       auto density = std::make_shared<MagicDensity<dim,adouble> >(F0,materialParameters);
-      elasticEnergy = std::make_shared<Elasticity::LocalIntegralEnergy<GridView,
-                                                                       FEBasis::LocalView::Tree::FiniteElement,
+      elasticEnergy = std::make_shared<Elasticity::LocalIntegralEnergy<FEBasis::LocalView,
                                                                        adouble> >(density);
     }
 
     if (parameterSet.get<std::string>("energy") == "burkholder")
     {
       auto density = std::make_shared<BurkholderDensity<dim,adouble> >(F0,materialParameters);
-      elasticEnergy = std::make_shared<Elasticity::LocalIntegralEnergy<GridView,
-                                                                       FEBasis::LocalView::Tree::FiniteElement,
+      elasticEnergy = std::make_shared<Elasticity::LocalIntegralEnergy<FEBasis::LocalView,
                                                                        adouble> >(density);
     }
 
@@ -743,23 +704,16 @@ int main (int argc, char *argv[]) try
     if (parameterSet.get<std::string>("energy") == "voss")
     {
       auto density = std::make_shared<VossDensity<dim,adouble> >(F0,materialParameters);
-      elasticEnergy = std::make_shared<Elasticity::LocalIntegralEnergy<GridView,
-                                                                       FEBasis::LocalView::Tree::FiniteElement,
+      elasticEnergy = std::make_shared<Elasticity::LocalIntegralEnergy<FEBasis::LocalView,
                                                                        adouble> >(density);
     }
 
     if(!elasticEnergy)
       DUNE_THROW(Exception, "Error: Selected energy not available!");
 
-    LocalADOLCStiffness<GridView,
-                        FEBasis::LocalView::Tree::FiniteElement,
-                        SolutionType> localADOLCStiffness(elasticEnergy.get());
+    auto localADOLCStiffness = std::make_shared<Elasticity::LocalADOLCStiffness<FEBasis::LocalView> >(elasticEnergy);
 
-    // dune-fufem-style FE basis for the transition from dune-fufem to dune-functions
-    typedef DuneFunctionsBasis<FEBasis> FufemFEBasis;
-    FufemFEBasis fufemFEBasis(feBasis);
-
-    FEAssembler<FufemFEBasis,SolutionType> assembler(fufemFEBasis, &localADOLCStiffness);
+    Elasticity::FEAssembler<FEBasis,SolutionType> assembler(feBasis, localADOLCStiffness);
 
     // /////////////////////////////////////////////////
     //   Create a Riemannian trust-region solver