Basis separation for edge/corner dofs.

docker/Dockerfile

@@ -29,8 +29,11 @@ WORKDIR /usr/lib
 RUN arch=$(uname -m) && sudo ln -s /usr/lib/${arch}-linux-gnu/libscalapack-openmpi.so.2.1.0 ./
 
 # re-install mfem with mumps, with cmake
+WORKDIR $LIB_DIR
+RUN sudo rm -rf mfem_parallel
+RUN sudo git clone https://github.com/mfem/mfem.git mfem_parallel
 WORKDIR $LIB_DIR/mfem
-RUN sudo git checkout v4.6
+RUN sudo git checkout v4.8
 WORKDIR $LIB_DIR/mfem/build
 RUN sudo -E cmake .. -DBUILD_SHARED_LIBS=YES -DMFEM_USE_MPI=YES -DMFEM_USE_GSLIB=${MG} -DMFEM_USE_METIS=YES -DMFEM_USE_METIS_5=YES -DMFEM_USE_MUMPS=YES -DGSLIB_DIR="$LIB_DIR/gslib/build" -DMUMPS_DIR="$LIB_DIR/mumps/build/local"
 RUN sudo -E make -j 16

include/component_topology_handler.hpp

@@ -156,6 +156,8 @@ class ComponentTopologyHandler : public TopologyHandler
    virtual void TransferToGlobal(Array<GridFunction*> &us, Array<GridFunction*> &global_u, const int &num_var)
    { mfem_error("ComponentTopologyHandler does not yet support global grid function/mesh!\n"); }
 
+   std::map<DofTag, Array<int>> SplitDofsByBoundary(FiniteElementSpace &fes);
+
 protected:
    // Get vertex orientation of face2 (from mesh2) with respect to face1 (mesh1).
    int GetOrientation(BlockMesh *comp1, const Element::Type &be_type, const Array<int> &vtx1, const Array<int> &vtx2);

include/topology_handler.hpp

@@ -51,6 +51,56 @@ struct TopologyData {
    Array<int> *global_bdr_attributes = NULL;
 };
 
+struct DofTag
+{
+   const Array<int> battr;
+
+   DofTag(const Array<int> &battrs)
+      : battr(battrs)
+   {
+      if (!battr.IsSorted())
+         mfem_error("DofTag: boundary attributes must be sorted beforehand!\n");
+   }
+
+   bool operator==(const DofTag &dtag) const
+   {
+      if (battr.Size() != dtag.battr.Size())
+         return false;
+
+      for (int k = 0; k < battr.Size(); k++)
+         if (battr[k] != dtag.battr[k])
+            return false;
+
+      return true;
+   }
+
+   bool operator<(const DofTag &dtag) const
+   {
+      if (battr.Size() == dtag.battr.Size())
+      {
+         for (int k = 0; k < battr.Size(); k++)
+            if (battr[k] != dtag.battr[k])
+               return (battr[k] < dtag.battr[k]);
+
+         return false;
+      }
+      else
+         return (battr.Size() < dtag.battr.Size());
+   }
+
+   std::string print() const
+   {
+      std::string tag = "(";
+      for (int k = 0 ; k < battr.Size(); k++)
+      {
+         tag += std::to_string(battr[k]) + ",";
+      }
+      tag += ")";
+
+      return tag;
+   }
+};
+
 const TopologyHandlerMode SetTopologyHandlerMode();
 
 class TopologyHandler
@@ -98,6 +148,8 @@ class TopologyHandler
    virtual Mesh* GetMesh(const int k) = 0;
    virtual Mesh* GetGlobalMesh() = 0;
 
+   std::map<DofTag, Array<int>> SplitBoundaryDofs(Mesh* const mesh, const FiniteElementSpace *fes) const;
+
    /*
       Methods only for ComponentTopologyHandler 
    */

sketches/CMakeLists.txt

@@ -28,6 +28,7 @@ add_executable(ns_mms ns_mms.cpp $<TARGET_OBJECTS:scaleupROMObj>)
 add_executable(ns_dg_mms ns_dg_mms.cpp $<TARGET_OBJECTS:scaleupROMObj>)
 add_executable(ns_rom ns_rom.cpp $<TARGET_OBJECTS:scaleupROMObj>)
 add_executable(usns usns.cpp $<TARGET_OBJECTS:scaleupROMObj>)
+add_executable(dof_separation dof_separation.cpp $<TARGET_OBJECTS:scaleupROMObj>)
 
 file(COPY inputs/gen_interface.yml DESTINATION ${CMAKE_BINARY_DIR}/sketches/inputs)
 file(COPY meshes/2x2.mesh DESTINATION ${CMAKE_BINARY_DIR}/sketches/meshes)

sketches/dof_separation.cpp

@@ -0,0 +1,216 @@
+// Copyright 2023 Lawrence Livermore National Security, LLC. See the top-level LICENSE file for details.
+//
+// SPDX-License-Identifier: MIT
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+#include <algorithm>
+#include "etc.hpp"
+
+using namespace std;
+using namespace mfem;
+
+struct DofTag
+{
+   Array<int> battr;
+
+   DofTag(const Array<int> &battrs)
+      : battr(battrs)
+   { battr.Sort(); }
+
+   bool operator==(const DofTag &dtag) const
+   {
+      if (battr.Size() != dtag.battr.Size())
+         return false;
+
+      for (int k = 0; k < battr.Size(); k++)
+         if (battr[k] != dtag.battr[k])
+            return false;
+
+      return true;
+   }
+
+   bool operator<(const DofTag &dtag) const
+   {
+      if (battr.Size() == dtag.battr.Size())
+      {
+         for (int k = 0; k < battr.Size(); k++)
+            if (battr[k] != dtag.battr[k])
+               return (battr[k] < dtag.battr[k]);
+
+         return false;
+      }
+      else
+         return (battr.Size() < dtag.battr.Size());
+   }
+
+   std::string print() const
+   {
+      std::string tag = "(";
+      for (int k = 0 ; k < battr.Size(); k++)
+      {
+         tag += std::to_string(battr[k]) + ",";
+      }
+      tag += ")";
+
+      return tag;
+   }
+};
+
+int main(int argc, char *argv[])
+{
+   StopWatch chrono;
+
+   // 1. Parse command-line options.
+   const char *mesh_file = "../data/star.mesh";
+   int order = 1;
+   int refine = 0;
+   bool pa = false;
+   const char *device_config = "cpu";
+   bool visualization = 1;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree).");
+   args.AddOption(&refine, "-r", "--refine",
+                  "Number of refinements.");
+   args.Parse();
+   if (!args.Good())
+   {
+      args.PrintUsage(cout);
+      return 1;
+   }
+   args.PrintOptions(cout);
+
+   // 3. Read the mesh from the given mesh file. We can handle triangular,
+   //    quadrilateral, tetrahedral, hexahedral, surface and volume meshes with
+   //    the same code.
+   Mesh mesh(mesh_file, 1, 1);
+   int dim = mesh.Dimension();
+
+   // 4. Refine the mesh to increase the resolution. In this example we do
+   //    'ref_levels' of uniform refinement. We choose 'ref_levels' to be the
+   //    largest number that gives a final mesh with no more than 10,000
+   //    elements.
+   for (int l = 0; l < refine; l++)
+   {
+      mesh.UniformRefinement();
+   }
+
+   // 5. Define a finite element space on the mesh. Here we use the
+   //    Raviart-Thomas finite elements of the specified order.
+   // FiniteElementCollection *hdiv_coll(new RT_FECollection(order, dim));
+   L2_FECollection l2_coll(order, dim);
+   H1_FECollection h1_coll(order+1, dim);
+   L2_FECollection pl2_coll(order, dim);
+   H1_FECollection ph1_coll(order, dim);
+
+   FiniteElementSpace fes(&mesh, &h1_coll);
+   FiniteElementSpace ufes(&mesh, &l2_coll, dim);
+   FiniteElementSpace pfes(&mesh, &ph1_coll);
+
+   // 6. Define the BlockStructure of the problem, i.e. define the array of
+   //    offsets for each variable. The last component of the Array is the sum
+   //    of the dimensions of each block.
+   Array<int> block_offsets(dim+2); // number of variables + 1
+   block_offsets[0] = 0;
+   for (int d = 1; d <= dim; d++)
+      block_offsets[d] = ufes.GetNDofs();
+   block_offsets[dim+1] = pfes.GetVSize();
+   block_offsets.PartialSum();
+
+   std::cout << "***********************************************************\n";
+   for (int d = 1; d < block_offsets.Size(); d++)
+      printf("dim(%d) = %d\n", d, block_offsets[d] - block_offsets[d-1]);
+   printf("dim(q) = %d\n", block_offsets.Last());
+   std::cout << "***********************************************************\n";
+
+   printf("ufes vsize: %d\n", ufes.GetVSize());
+   printf("pfes vsize: %d\n", pfes.GetVSize());
+
+   Array<int> udofs(ufes.GetVSize());
+   Array<int> pdofs(pfes.GetVSize());
+   udofs = 0;
+   pdofs = 0;
+
+   int max_battr = -1;
+   for (int i = 0; i < pfes.GetNBE(); i++)
+      max_battr = max(max_battr, mesh.GetBdrAttribute(i));
+
+   Array2D<int> pbattr(pfes.GetVSize(), max_battr);
+   pbattr = 0;
+
+   Array<int> vdofs;
+   FaceElementTransformations *tr;
+   printf("pfes nbe: %d\n", pfes.GetNBE());
+   for (int i = 0; i < pfes.GetNBE(); i++)
+   {
+      const int bdr_attr = mesh.GetBdrAttribute(i);
+
+      pfes.GetBdrElementVDofs(i, vdofs);
+      for (int k = 0 ; k < vdofs.Size(); k++)
+         printf("%d\t", vdofs[k]);
+
+      tr = mesh.GetBdrFaceTransformations(i);
+      if (tr != NULL)
+      {
+         pfes.GetElementVDofs(tr->Elem1No, vdofs);
+         printf(" / ");
+         for (int k = 0 ; k < vdofs.Size(); k++)
+         {
+            printf("%d\t", vdofs[k]);
+            pdofs[vdofs[k]] += 1;
+            pbattr(vdofs[k], bdr_attr-1) = 1;
+         }
+      }
+      printf("\n");
+   }
+
+   for (int k = 0; k < pdofs.Size(); k++)
+   {
+      printf("pdofs[%d] = %d / ", k, pdofs[k]);
+      int nbattr = 0;
+      for (int d = 0; d < max_battr; d++)
+      {
+         printf("%d ", pbattr(k, d));
+         nbattr += pbattr(k, d);
+      }
+      printf(": %d\n", nbattr);
+   }
+
+   std::map<DofTag, Array<int>> pdofbattr;
+   for (int k = 0; k < pfes.GetVSize(); k++)
+   {
+      Array<int> battrs(0);
+      for (int d = 0; d < max_battr; d++)
+      {
+         if (pbattr(k, d) > 0)
+            battrs.Append(d+1);
+      }
+      battrs.Sort();
+
+      DofTag pdtag(battrs);
+      if (!pdofbattr.count(pdtag))
+         pdofbattr[pdtag] = Array<int>(0);
+
+      pdofbattr[pdtag].Append(k);
+   }
+
+   for(std::map<DofTag, Array<int>>::iterator iter = pdofbattr.begin(); iter != pdofbattr.end(); ++iter)
+   {
+      const DofTag *key = &(iter->first);
+      const Array<int> *val = &(iter->second);
+
+      printf("%s: ", key->print().c_str());
+      for (int k = 0 ; k < val->Size(); k++)
+      {
+         printf("%d ", (*val)[k]);
+      }
+      printf("\n");
+   }
+
+   return 0;
+}

sketches/domain_decomposition.cpp

@@ -132,7 +132,7 @@ int main(int argc, char* argv[])
     for (int m = 0; m < 1; m++){
       printf("Mesh %d\n", m);
       for (int i = 0; i < meshes[m]->GetNBE(); i++) {
-        int fn = meshes[m]->GetBdrFace(i);
+        int fn = meshes[m]->GetBdrElementFaceIndex(i);
         printf("Boundary Element %d : Face %d\n", i, fn);
       }
       // for (int i = 0; i < meshes[m]->GetNBE(); i++) {
@@ -146,7 +146,7 @@ int main(int argc, char* argv[])
       //   printf("Boundary element %d - Adjacent element %d\n", i, elem_id);
       //   printf("Face index: %d, face orientation : %d\n", face_info / 64, face_info % 64);
       //
-      //   int fn = meshes[m]->GetBdrFace(i);
+      //   int fn = meshes[m]->GetBdrElementFaceIndex(i);
       //   int face_inf[2];
       //   meshes[m]->GetFaceInfos(fn, &face_inf[0], &face_inf[1]);
       //   printf("From faces_info - Face index: %d, face orientation : %d\n", face_inf[0] / 64, face_inf[0] % 64);

sketches/generate_interface.cpp

@@ -46,7 +46,7 @@ int main(int argc, char *argv[])
          for (int be = 0; be < meshes[m]->GetNBE(); be++)
          {
             int battr = meshes[m]->GetBdrAttribute(be);
-            int face = meshes[m]->GetBdrFace(be);
+            int face = meshes[m]->GetBdrElementFaceIndex(be);
             Array<int> vtx;
             // comp.GetFaceVertices(face, vtx);
             meshes[m]->GetBdrElementVertices(be, vtx);
@@ -115,7 +115,7 @@ int main(int argc, char *argv[])
       for (int be = 0; be < comp.GetNBE(); be++)
       {
          int battr = comp.GetBdrAttribute(be);
-         int face = comp.GetBdrFace(be);
+         int face = comp.GetBdrElementFaceIndex(be);
          Array<int> vtx;
          // comp.GetFaceVertices(face, vtx);
          comp.GetBdrElementVertices(be, vtx);
@@ -168,7 +168,7 @@ int main(int argc, char *argv[])
       // These will be stored in a 'interface pair' file, with the comp/battr pair above.
       Array<std::map<int,int>*> vtx_2to1(n_interface_pairs);
       // Array<std::map<int,int>*> be_2to1(n_interface_pairs);
-      Array<Array<std::pair<int,int>>*> be_pairs(n_interface_pairs);
+      // Array<Array<std::pair<int,int>>*> be_pairs(n_interface_pairs);
       for (int k = 0; k < n_interface_pairs; k++)
       {
          int bidx1 = comp.bdr_attributes.Find(if_battr1[k]);
@@ -177,18 +177,18 @@ int main(int argc, char *argv[])
          assert(bdr_vtx[bidx1]->Size() == bdr_vtx[bidx2]->Size());
 
          vtx_2to1[k] = new std::map<int,int>;
-         be_pairs[k] = new Array<std::pair<int,int>>(0);
+         // be_pairs[k] = new Array<std::pair<int,int>>(0);
       }
 
       (*vtx_2to1[0])[1] = 3;
       (*vtx_2to1[0])[0] = 2;
       // (*be_pairs[0])[0] = 1;
-      be_pairs[0]->Append({1, 0});
+      // be_pairs[0]->Append({1, 0});
 
       (*vtx_2to1[1])[0] = 1;
       (*vtx_2to1[1])[2] = 3;
       // (*be_pairs[0])[2] = 3;
-      be_pairs[1]->Append({3, 2});
+      // be_pairs[1]->Append({3, 2});
 
       // These are the informations that will be stored in the 'global topology' file.
       Array<int> global_comps(topol_data.numSub);
@@ -252,8 +252,8 @@ int main(int argc, char *argv[])
       //          tmp.BE2 = (*be_pair)[be].second;
 
       //          // use the face index from each component mesh.
-      //          int f1 = meshes[tmp.Mesh1]->GetBdrFace(tmp.BE1);
-      //          int f2 = meshes[tmp.Mesh2]->GetBdrFace(tmp.BE2);
+      //          int f1 = meshes[tmp.Mesh1]->GetBdrElementFaceIndex(tmp.BE1);
+      //          int f2 = meshes[tmp.Mesh2]->GetBdrElementFaceIndex(tmp.BE2);
       //          int Inf1, Inf2, dump;
       //          meshes[tmp.Mesh1]->GetFaceInfos(f1, &Inf1, &dump);
       //          meshes[tmp.Mesh2]->GetFaceInfos(f2, &Inf2, &dump);