Implement periodic boundary conditions

src/quasiconvexity-test.cc

@@ -11,6 +11,7 @@
 #include <dune/fufem/utilities/adolcnamespaceinjections.hh>
 
 #include <dune/common/bitsetvector.hh>
+#include <dune/common/float_cmp.hh>
 #include <dune/common/parametertree.hh>
 #include <dune/common/parametertreeparser.hh>
 
@@ -28,6 +29,7 @@
 #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>
@@ -267,41 +269,106 @@ int main (int argc, char *argv[]) try
 #endif
 
   // FE basis spanning the FE space that we are working in
-  typedef Dune::Functions::LagrangeBasis<GridView, order> FEBasis;
+  //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);
 
-  // /////////////////////////////////////////
-  //   Read Dirichlet values
-  // /////////////////////////////////////////
+#if 0
+  // 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
 
-  BitSetVector<1> dirichletVertices(gridView.size(dim), false);
+  // Check whether two points are equal on R/Z x R/Z
+  auto equivalent = [](const FieldVector<double,2>& x, const FieldVector<double,2>& y)
+  {
+    return ( (FloatCmp::eq(x[0],y[0]) or FloatCmp::eq(x[0]+1,y[0]) or FloatCmp::eq(x[0]-1,y[0]))
+           and (FloatCmp::eq(x[1],y[1]) or FloatCmp::eq(x[1]+1,y[1]) or FloatCmp::eq(x[1]-1,y[1])) );
+  };
 
-  const GridView::IndexSet& indexSet = gridView.indexSet();
+  // Check whether two points are equal
+  auto equal = [](const FieldVector<double,2>& x, const FieldVector<double,2>& y)
+  {
+    return FloatCmp::eq(x[0],y[0]) && FloatCmp::eq(x[1],y[1]);
+  };
 
-  // Make Python function that computes which vertices are on the Dirichlet boundary,
-  // based on the vertex positions.
+  bool periodic = parameterSet.get<bool>("periodic");
+  if (periodic)
+  {
+    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)});
+        }
+  }
 
-  std::string lambda = std::string("lambda x: (") + parameterSet.get<std::string>("dirichletVerticesPredicate") + std::string(")");
-  PythonFunction<FieldVector<double,dim>, bool> pythonDirichletVertices(Python::evaluate(lambda));
+  feBasis.preBasis().initializeIndices();
 
-  for (auto&& v : vertices(gridView))
+  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++)
   {
-    bool isDirichlet;
-    pythonDirichletVertices.evaluate(v.geometry().corner(0), isDirichlet);
-    dirichletVertices[indexSet.index(v)] = isDirichlet;
+    lagrangeNodes[i][0] = lagrangeNodesX[i];
+    lagrangeNodes[i][1] = lagrangeNodesY[i];
   }
+#endif
+
+  ///////////////////////////////////////////
+  //   Set Dirichlet values
+  ///////////////////////////////////////////
 
+  BitSetVector<1> dirichletVertices(gridView.size(dim), false);
+
+  BitSetVector<dim> dirichletDofs(feBasis.size(), false);
+
+  if (periodic)
+  {
+    // Only the lower-left corner is Dirichlet boundary
+    for (size_t i=0; i<lagrangeNodes.size(); i++)
+      if (equivalent(lagrangeNodes[i], {0,0}))
+      {
+        dirichletDofs[i] = true;
+        break;
+      }
+  }
+  else
+  {
+    // The entire boundary is Dirichlet boundary
+    dirichletVertices.setAll();
     BoundaryPatch<GridView> dirichletBoundary(gridView, dirichletVertices);
 
     BitSetVector<1> dirichletNodes(feBasis.size(), false);
     constructBoundaryDofs(dirichletBoundary,feBasis,dirichletNodes);
 
-  BitSetVector<dim> dirichletDofs(feBasis.size(), false);
     for (size_t i=0; i<feBasis.size(); i++)
       if (dirichletNodes[i][0])
         for (int j=0; j<dim; j++)
           dirichletDofs[i][j] = true;
-
+  }
 #if 0
   // HACK!
 #warning Teichmueller hack!
@@ -380,7 +447,6 @@ int main (int argc, char *argv[]) try
       if (not inFile)
         DUNE_THROW(IOError, "File " << parameterSet.get<std::string>("initialIterateFile") << " could not be opened.");
 
-      int vertex = 0;
       for (std::string line; std::getline(inFile, line);)
       {
         // The first line starts with 'x' -- skip it
@@ -398,21 +464,17 @@ int main (int argc, char *argv[]) try
         FieldVector<double,2> pos{atof(result[0].c_str()), atof(result[1].c_str())};
         FieldVector<double,2> u{atof(result[2].c_str()), atof(result[3].c_str())};
 
-        //x[vertex] = u - (F0-Id)*pos;
+        // phi = u - (F0-Id)*pos;
         FieldVector<double,2> phi = { u[0] - (F0[0][0]-1)*pos[0] -  F0[0][1]   *pos[1],
                                       u[1] -  F0[1][0]   *pos[0] - (F0[1][1]-1)*pos[1] };
 
-        for (const auto& vert : vertices(gridView))
-          if ( (vert.geometry().corner(0)-pos).two_norm() < 1e-8)
+        for (std::size_t i=0; i<lagrangeNodes.size(); i++)
+          if ( (periodic && equivalent(lagrangeNodes[i], pos)) || (!periodic && equal(lagrangeNodes[i], pos)) )
           {
-            x[gridView.indexSet().index(vert)] = phi;
-            std::cout << "pos: " << pos << ",  vertex: " << gridView.indexSet().index(vert) << std::endl;
+            x[i] = phi;
+            //std::cout << "pos: " << pos << ",  lagrange node: " << i << std::endl;
             break;
           }
-        //x[vertex] = phi;
-
-        //std::cout << "pos: " << pos << ",  phi: " << x[vertex] << std::endl;
-        vertex++;
       }
 
       inFile.close();
@@ -425,13 +487,13 @@ int main (int argc, char *argv[]) try
 
   // Output initial iterate (of homotopy loop)
   SolutionType identity;
-  Dune::Functions::interpolate(powerBasis, identity, [&F0](FieldVector<double,dim> x)
+  Dune::Functions::interpolate(powerBasis, identity, [&F0](FieldVector<double,dim> pos)
   {
-    return FieldVector<double,2>{ (F0[0][0]-1)*x[0] + F0[0][1]*x[1], F0[1][0]*x[0] + (F0[1][1]-1)*x[1] };
+    return FieldVector<double,2>{ (F0[0][0]-1)*pos[0] + F0[0][1]*pos[1], F0[1][0]*pos[0] + (F0[1][1]-1)*pos[1] };
   });
 
   SolutionType displacement = x;
-  displacement -= identity;
+  displacement += identity;
 
   auto displacementFunction = Dune::Functions::makeDiscreteGlobalBasisFunction<FieldVector<double,dim>>(feBasis, displacement);
   auto localDisplacementFunction = localFunction(displacementFunction);
@@ -635,7 +697,7 @@ int main (int argc, char *argv[]) try
 #endif
     // Compute the displacement
     auto displacement = x;
-    displacement -= identity;
+    displacement += identity;
 
     auto displacementFunction = Dune::Functions::makeDiscreteGlobalBasisFunction<FieldVector<double,dim>>(feBasis, displacement);