added dynamic determination of orientation in 2D

fem/fespace.cpp

@@ -1298,17 +1298,89 @@ void FiniteElementSpace::MakeVDimMatrix(SparseMatrix &mat) const
    delete vmat;
 }
 
+void FiniteElementSpace::MakePeriodic(const std::vector<int> &v2v)
+{
+   Array<int> source_vertices, target_vertices, dofs, tDofs, lDofs;
+   for (int k=0; k < mesh->GetNBE(); k++)
+   {
+      mesh->GetBdrElementVertices(k, source_vertices);
+      if (std::all_of(source_vertices.begin(), source_vertices.end(), [&](int v){ return v2v[v] != v; }))
+      {
+         GetBdrElementDofs(k, dofs);
+         lDofs.Append(dofs);
+         int index = -1;
+         for (int l=0; l < mesh->GetNBE(); l++)
+         {
+            mesh->GetBdrElementVertices(l, target_vertices);
+            if (std::all_of(source_vertices.begin(), source_vertices.end(), [&](int v){ return target_vertices.Find(v2v[v]) != -1; }))
+            {
+               index = l;
+               break;
+            }
+         }
+         MFEM_ASSERT(index != -1, "could not find boundary element containing mapped vertex.")
+         GetBdrElementDofs(index, dofs);
+         // Append in the order consistent with the orientation of the target element vertices / dofs
+         if (target_vertices[0] == v2v[source_vertices[0]])
+            tDofs.Append(dofs);
+         else
+            for (int m=0; m < dofs.Size(); m++)
+               tDofs.Append(dofs[dofs.Size()-1-m]);
+      }
+   }
+//   std::cout << "tdofs" << std::endl;
+//   tDofs.Print();
+//   std::cout << "ldofs" << std::endl;
+//   lDofs.Print();
+
+   Array<int> local2true(ndofs);
+   int count = 0;
+   for (int k=0; k < ndofs; k++)
+   {
+      local2true[k] = count;
+      if (lDofs.Find(k) == -1)
+         count++;
+   }
+   for (int k=0; k < lDofs.Size(); k++)
+      local2true[lDofs[k]] = local2true[tDofs[k]];
+
+//   std::cout << "map" << std::endl;
+//   local2true.Print();
+
+
+   periodicProlongationMatrix =
+      std::make_unique<SparseMatrix>(ndofs, count);
+   for (int i=0; i < ndofs; i++)
+      periodicProlongationMatrix->Add(i,local2true[i],1.0);
+   periodicProlongationMatrix->Finalize();
+//   std::cout << "Prolongation" << std::endl;
+//   periodicProlongationMatrix->ToDenseMatrix()->Print(std::cout);
+
+   periodicRestrictionMatrix =
+      std::make_unique<SparseMatrix>(count, ndofs);
+   int i = 0;
+   for (int j=0; j < ndofs; j++)
+   {
+      if (lDofs.Find(j) != -1)
+         continue;
+      periodicRestrictionMatrix->Add(i,j,1.0);
+      i++;
+   }
+   periodicRestrictionMatrix->Finalize();
+ //  std::cout << "Restriction:" << std::endl;
+ //  periodicRestrictionMatrix->ToDenseMatrix()->Print(std::cout);
+}
 
 const SparseMatrix* FiniteElementSpace::GetConformingProlongation() const
 {
-   if (Conforming()) { return NULL; }
+   if (Conforming()) { return periodicProlongationMatrix.get(); }
    if (!cP_is_set) { BuildConformingInterpolation(); }
    return cP.get();
 }
 
 const SparseMatrix* FiniteElementSpace::GetConformingRestriction() const
 {
-   if (Conforming()) { return NULL; }
+   if (Conforming()) { return periodicRestrictionMatrix.get(); }
    if (!cP_is_set) { BuildConformingInterpolation(); }
    if (cR && !R_transpose) { R_transpose.reset(new TransposeOperator(*cR)); }
    return cR.get();

fem/fespace.hpp

@@ -223,6 +223,9 @@ class FiniteElementSpace
    friend void Mesh::Swap(Mesh &, bool);
    friend class LORBase;
 
+   std::unique_ptr<SparseMatrix> periodicProlongationMatrix;
+   std::unique_ptr<SparseMatrix> periodicRestrictionMatrix;
+
 protected:
    /// The mesh that FE space lives on (not owned).
    Mesh *mesh;
@@ -593,6 +596,9 @@ public:
    /// The returned SparseMatrix is owned by the FiniteElementSpace.
    const SparseMatrix *GetHpConformingRestriction() const;
 
+void MakePeriodic(const std::vector<int> &v2v);
+
+
    /// The returned Operator is owned by the FiniteElementSpace.
    virtual const Operator *GetProlongationMatrix() const
    { return GetConformingProlongation(); }