make style

examples/MeshPart/makefile

@@ -0,0 +1,61 @@
+# Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
+# at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+# LICENSE and NOTICE for details. LLNL-CODE-806117.
+#
+# This file is part of the MFEM library. For more information and source code
+# availability visit https://mfem.org.
+#
+# MFEM is free software; you can redistribute it and/or modify it under the
+# terms of the BSD-3 license. We welcome feedback and contributions, see file
+# CONTRIBUTING.md for details.
+
+# Use the MFEM build directory
+MFEM_DIR ?= ../..
+MFEM_BUILD_DIR ?= ../..
+SRC = $(if $(MFEM_DIR:../..=),$(MFEM_DIR)/examples/MeshPart/,)
+CONFIG_MK = $(MFEM_BUILD_DIR)/config/config.mk
+
+MFEM_LIB_FILE = mfem_is_not_built
+-include $(CONFIG_MK)
+
+SEQ_EXAMPLES = test_mesh_partition
+PAR_EXAMPLES =
+ifeq ($(MFEM_USE_MPI),NO)
+   EXAMPLES = $(SEQ_EXAMPLES)
+else
+   EXAMPLES = $(PAR_EXAMPLES) $(SEQ_EXAMPLES)
+endif
+
+.SUFFIXES:
+.SUFFIXES: .o .cpp .mk
+.PHONY: all clean
+.PRECIOUS: %.o
+
+COMMON_O= mesh_partition.o
+
+# Remove built-in rules
+%: %.cpp
+%.o: %.cpp
+
+all: $(EXAMPLES)
+
+# Rules for building the EXAMPLES
+
+%: $(SRC)%.cpp $(COMMON_O) $(MFEM_LIB_FILE) $(CONFIG_MK)
+	$(MFEM_CXX) $(MFEM_FLAGS) $< -o $@ $(COMMON_O) $(MFEM_LIBS)
+
+# Rules for compiling miniapp dependencies
+$(COMMON_O) $($(EXAMPLES)): \
+%.o: $(SRC)%.cpp $(SRC)%.hpp $(CONFIG_MK)
+	$(MFEM_CXX) $(MFEM_FLAGS) -c $(<) -o $(@)
+
+# Generate an error message if the MFEM library is not built and exit
+$(MFEM_LIB_FILE):
+	$(error The MFEM library is not built)
+
+clean:
+	rm -f *.o *~ $(SEQ_EXAMPLES) $(PAR_EXAMPLES)
+	rm -rf *.dSYM *.TVD.*breakpoints
+	rm output/*
+
+

examples/MeshPart/mesh_partition.cpp

@@ -0,0 +1,301 @@
+
+#include "mesh_partition.hpp"
+
+
+Subdomain::Subdomain(const Mesh & mesh0_) 
+: mesh0(&mesh0_), dim(mesh0->Dimension()), sdim(mesh0->SpaceDimension()) 
+{ 
+   // MFEM_VERIFY(dim == 3, "Only 3D domains for now are supported");
+   if(mesh0->NURBSext)
+   {
+      MFEM_ABORT("Nurbs meshes are not supported yet");
+   }
+}
+
+void Subdomain::BuildSubMesh(const Array<int> & elems, const entity_type & etype)
+{
+   // if mesh nodes are defined we use them for the vertices
+   // otherwise we use the vertices them selfs
+   cout << "entity type " << etype << endl;
+
+   int nv = mesh0->GetNV();
+   int subelems = elems.Size();
+   Array<int> vmarker(nv); vmarker = 0;
+   int numvertices = 0;
+   for (int ie = 0; ie<subelems; ie++)
+   {
+      int el = elems[ie];
+      Array<int> vertices;
+      switch (etype)
+      {
+      case 0: mesh0->GetElementVertices(el,vertices); break;
+      case 1: mesh0->GetFaceVertices(el,vertices); break;
+      default:
+         MFEM_ABORT("Wrong entity type choice");
+         break;
+      }
+      for (int iv=0; iv<vertices.Size(); iv++)
+      {
+         int v = vertices[iv];
+         if (vmarker[v]) continue;
+         vmarker[v] = 1;
+         numvertices++;
+      }
+   }
+   cout << "Num of new vertices: " << numvertices << endl;
+
+   // Construct new mesh
+   Mesh * meshptr = nullptr;
+   switch (etype)
+   {
+   case 0: 
+      mesh = new Mesh(dim,numvertices, subelems); 
+      element_map = elems;
+      meshptr = mesh;
+      break;
+   case 1: 
+      surface_mesh = new Mesh(dim-1,numvertices, subelems,0,sdim);
+      surface_element_map = elems;
+      meshptr = surface_mesh;
+      break;
+      default:
+         MFEM_ABORT("Wrong entity type choice");
+         break;
+   }
+
+   Vector values;
+   const GridFunction * nodes0 = mesh0->GetNodes();
+   int vk = 0;
+   // if (nodes0) // this is NOT NEEDED here
+   // {
+   //    vcoords.SetSize(sdim, mesh0->GetNV());
+   //    for (int i = 0; i< sdim; i++)
+   //    {
+   //       nodes0->GetNodalValues(values,i+1);
+   //       vcoords.SetRow(i,values);
+   //       cout << "values size = " << values.Size() << endl;
+   //    }
+   //    for (int iv = 0; iv<mesh0->GetNV(); ++iv)
+   //    {
+   //       if (!vmarker[iv]) continue;
+   //       meshptr->AddVertex(vcoords.GetColumn(iv));
+   //       vmarker[iv] = ++vk;
+   //    }
+   // }
+   // else
+   {
+      for (int iv = 0; iv<mesh0->GetNV(); ++iv)
+      {
+         if (!vmarker[iv]) continue;
+         meshptr->AddVertex(mesh0->GetVertex(iv));
+         vmarker[iv] = ++vk;
+      }
+   }
+   // Add elements
+   for (int ie = 0; ie<subelems; ie++)
+   {
+      const Element * el = nullptr;
+      switch (etype)
+      {
+         case 0: el = mesh0->GetElement(elems[ie]); break;
+         case 1: el = mesh0->GetFace(elems[ie]); break;
+         default: MFEM_ABORT("Wrong entity type choice"); break;
+      }
+      Element * nel = meshptr->NewElement(el->GetGeometryType());
+      int nv0 = el->GetNVertices();
+      const int * v0 = el->GetVertices();
+      Array<int> v1(nv0);
+      for (int i=0; i<nv0; i++)
+      {
+         v1[i] = vmarker[v0[i]]-1;
+      }   
+      nel->SetVertices(v1.GetData());
+      meshptr->AddElement(nel);
+   }
+
+   meshptr->FinalizeTopology();
+   
+   if (nodes0)
+   {
+      cout << "nodes not null" << endl;
+      // Extract Nodes GridFunction and determine its type
+      const FiniteElementSpace * fes0 = nodes0->FESpace();
+
+      Ordering::Type ordering = fes0->GetOrdering();
+      int order = fes0->FEColl()->GetOrder();
+      bool discont = fes0->IsDGSpace();   
+
+      cout << "discont = " << discont << endl;
+      // Set curvature of the same type as original mesh
+      meshptr->SetCurvature(order, discont, sdim, ordering);
+
+      const FiniteElementSpace * fes1 = meshptr->GetNodalFESpace();
+      GridFunction * nodes = meshptr->GetNodes();
+
+      Array<int> vdofs0;
+      Array<int> vdofs;
+      Vector loc_vec;
+      
+      // Copy nodes to submesh
+      for (int e = 0; e < elems.Size(); e++)
+      {
+         fes1->GetElementVDofs(e, vdofs);
+         switch (etype)
+         {
+         case 0: 
+            fes0->GetElementVDofs(elems[e], vdofs0); 
+            nodes0->GetSubVector(vdofs0, loc_vec);
+            break;
+         case 1: 
+            if (!discont)
+            {
+               fes0->GetFaceVDofs(elems[e], vdofs0); 
+               nodes0->GetSubVector(vdofs0, loc_vec);
+            }
+            else
+            {
+               const FiniteElement * el = fes1->GetFE(e);
+               const IntegrationRule & ir = el->GetNodes(); 
+               int np = ir.GetNPoints();
+               FaceElementTransformations * Tr = 
+                  const_cast<Mesh *>(mesh0)->GetFaceElementTransformations(elems[e]);
+               int el1 = Tr->Elem1No;
+               loc_vec.SetSize(vdofs.Size());
+               for (int i = 0; i<np; i++)
+               {
+                  Tr->SetAllIntPoints(&ir[i]);
+                  const IntegrationPoint & ip = Tr->GetElement1IntPoint();
+                  Vector val;
+                  nodes0->GetVectorValue(el1,ip,val);
+                  for (int j = 0; j<val.Size(); j++)
+                  {
+                     loc_vec[i+j*np] = val[j];
+                  }
+               }
+            }
+            break;
+         default:
+            MFEM_ABORT("Wrong entity type choice");
+            break;
+         }
+         nodes->SetSubVector(vdofs, loc_vec);
+      }
+   }
+   
+   meshptr->Finalize();
+}
+
+
+void Subdomain::BuildDofMap(const entity_type & etype)
+{
+   Array<int> elems;
+   FiniteElementSpace * fesptr = nullptr;
+   const FiniteElementCollection *fec = fes0->FEColl();
+   switch(etype)
+   {
+   case 0: 
+      fesptr = new FiniteElementSpace(mesh,fec);
+      elems = element_map;
+      break;
+   case 1: 
+      fesptr = new FiniteElementSpace(surface_mesh,fec);
+      elems = surface_element_map;
+      break;
+   default:
+      MFEM_ABORT("Wrong entity type choice");
+      break;
+   }
+   
+   Array<int> dofs(fesptr->GetVSize());
+   for (int iel = 0; iel<elems.Size(); ++iel)
+   {
+      // index in the global mesh
+      int iel_idx = elems[iel];
+      // get the dofs of this element
+      Array<int> ldofs;
+      Array<int> gdofs;
+      switch(etype)
+      {
+         case 0: fes0->GetElementVDofs(iel_idx,gdofs); break;
+         case 1: fes0->GetFaceVDofs(iel_idx,gdofs); break;
+         default:  MFEM_ABORT("Wrong entity type");  break;
+      }
+      fesptr->GetElementDofs(iel,ldofs);
+         // the sizes have to match
+      MFEM_VERIFY(gdofs.Size() == ldofs.Size(),
+                  "Size inconsistency");
+     // loop through the dofs and take into account the signs;
+      int ndof = ldofs.Size();
+      for (int i = 0; i<ndof; ++i)
+      {
+         int ldof_ = ldofs[i];
+         int gdof_ = gdofs[i];
+         int ldof = (ldof_ >= 0) ? ldof_ : abs(ldof_) - 1;
+         int gdof = (gdof_ >= 0) ? gdof_ : abs(gdof_) - 1;
+         dofs[ldof] = gdof;
+      }
+   }
+   switch(etype)
+   {
+   case 0: 
+      dof_map = dofs; 
+      fes = fesptr;
+      break;
+   case 1: 
+      surface_dof_map = dofs; 
+      surface_fes = fesptr;
+      break;
+   default:  
+      MFEM_ABORT("Wrong entity type");  break;
+   }
+}
+
+void Subdomain::BuildProlongationMatrix(const entity_type & etype)
+{
+   Array<int> dofs;
+   SparseMatrix * Ptr = nullptr;
+   switch (etype)
+   {
+   case 0: 
+      if (!dof_map.Size()) BuildDofMap(etype);
+      dofs = dof_map; 
+      Ptr = P;
+      break;
+   case 1: 
+      if (!surface_dof_map.Size()) BuildDofMap(etype);
+      dofs = surface_dof_map; 
+      Ptr = Pf;
+      break;
+   default: 
+      MFEM_ABORT("Wrong entity type"); 
+      break;
+   }
+
+   int height = fes0->GetVSize();
+   int width = dofs.Size();
+   Ptr = new SparseMatrix(height,width);
+   for (int i = 0; i< dofs.Size(); i++)
+   {
+      int j = dofs[i];
+      Ptr->Set(j,i,1.);
+   }
+   Ptr->Finalize();
+
+   switch (etype)
+   {
+      case 0: P = Ptr; break;
+      case 1: Pf = Ptr; break;
+      default: MFEM_ABORT("Wrong entity type"); break;
+   }
+
+}
+
+Subdomain::~Subdomain()
+{
+   delete mesh;
+   delete surface_mesh;
+   delete fes; 
+   delete surface_fes;
+   delete P; 
+   delete Pf; 
+}

examples/MeshPart/mesh_partition.hpp

@@ -0,0 +1,106 @@
+#pragma once
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+using namespace std;
+using namespace mfem;
+
+
+class Subdomain
+{
+public:
+   enum entity_type
+   {
+      volume, 
+      surface
+   };
+private:
+   const Mesh *mesh0=nullptr;
+   int dim, sdim;
+   const FiniteElementSpace *fes0=nullptr;
+   DenseMatrix vcoords;
+   Mesh *mesh=nullptr; // Submesh 
+   Mesh *surface_mesh=nullptr; //Surface mesh 
+   FiniteElementSpace *fes=nullptr; // FE Space on the submesh
+   FiniteElementSpace *surface_fes=nullptr; // FE Space on the submesh
+
+   Array<int> element_map, surface_element_map;
+   Array<int> dof_map, surface_dof_map;
+   SparseMatrix * P=nullptr; 
+   SparseMatrix * Pf=nullptr; 
+   void BuildDofMap(const entity_type & etype);
+   void BuildProlongationMatrix(const entity_type & etype);
+   void BuildSubMesh(const Array<int> & elems, const entity_type & etype);
+
+public:
+   Subdomain(const Mesh & mesh_);
+   Mesh * GetSubMesh(const Array<int> & elems)
+   {
+      if(!mesh) BuildSubMesh(elems, entity_type::volume);
+      return mesh; 
+   }
+   Mesh * GetSurfaceMesh(const Array<int> & surface_elems)
+   {
+      if (!surface_mesh) BuildSubMesh(surface_elems, entity_type::surface);
+      return surface_mesh;
+   }
+   void SetFESpace(const FiniteElementSpace & fes0_)
+   {
+      fes0 = &fes0_;
+   }
+   void GetElementMap(Array<int> & element_map_)
+   {
+      element_map_ = element_map;
+   }
+   void GetFaceElementMap(Array<int> & surface_element_map_)
+   {
+      surface_element_map_ = surface_element_map;
+   }
+   void GetDofMap(Array<int> & dof_map_)
+   {
+      if (!dof_map.Size()) BuildDofMap(entity_type::volume);
+      dof_map_ = dof_map;
+   }
+   void GetSurfaceDofMap(Array<int> & surface_dof_map_)
+   {
+      if (!surface_dof_map.Size()) BuildDofMap(entity_type::surface);
+      surface_dof_map_ = surface_dof_map;
+   }
+   SparseMatrix * GetProlonationMatrix()
+   {
+      if (!P) BuildProlongationMatrix(entity_type::volume);
+      return P;
+   }
+   SparseMatrix * GetSurfaceProlonationMatrix()
+   {
+      if (!Pf) BuildProlongationMatrix(entity_type::surface);
+      return Pf;
+   }
+   FiniteElementSpace * GetSubFESpace(const entity_type & etype)
+   {
+      switch (etype)
+      {
+      case 0: 
+         if (!fes)
+         {
+            MFEM_VERIFY(mesh, "Volume mesh not built");
+            BuildDofMap(etype);
+         }
+         return fes; 
+         break;
+      case 1: 
+         if (!surface_fes)
+         {
+            MFEM_VERIFY(surface_mesh, "Surface mesh not built");
+            BuildDofMap(etype);
+         }
+         return surface_fes; 
+         break;
+      default: 
+         MFEM_ABORT("Wrong entity type"); 
+         return nullptr; 
+         break;
+      }
+   }
+   ~Subdomain();
+};

examples/MeshPart/test1.cpp

@@ -0,0 +1,96 @@
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+int main(int argc, char *argv[])
+{
+   // 1. Parse command line options
+   const char *mesh_file = "../../data/periodic-annulus-sector.msh";
+   int order = 1;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh", "Mesh file to use.");
+   args.AddOption(&order, "-o", "--order", "Finite element polynomial degree");
+   args.ParseCheck();
+
+   // 2. Read the mesh from the given mesh file, and refine once uniformly.
+   Mesh orig_mesh(mesh_file);
+   orig_mesh.CheckElementOrientation(true);
+   orig_mesh.CheckBdrElementOrientation(true);
+
+   // mesh.EnsureNodes();
+   // mesh.UniformRefinement();
+   // Array<int> elems({1,3,21,10,20,2,0});
+   Array<int> elems({0,1,2,3});
+   // int nel = orig_mesh.GetNE();
+   int nel = elems.Size();
+   // Array<int> elems(nel);
+   // for (int i = 0; i<nel; i++)
+   // {
+      // elems[i] = i;
+   // }
+   // elems.Print();
+
+   Mesh new_mesh = Mesh::ExtractMesh(orig_mesh,elems);
+   new_mesh.CheckElementOrientation(true);
+   new_mesh.CheckBdrElementOrientation(true);
+   {
+      char vishost[] = "localhost";
+      int  visport   = 19916;
+      socketstream mesh0_sock(vishost, visport);
+      mesh0_sock.precision(8);
+      mesh0_sock << "mesh\n" << orig_mesh << "keys n \n" << flush;
+
+      socketstream mesh1_sock(vishost, visport);
+      mesh1_sock.precision(8);
+      mesh1_sock << "mesh\n" << new_mesh << "keys n \n" << flush;
+   }
+
+   Array<int> faces;
+   for (int i = 0; i<orig_mesh.GetNBE(); i++)
+   {
+      if (orig_mesh.GetBdrAttribute(i) >= 1)
+         faces.Append(orig_mesh.GetBdrFace(i));
+   }
+      // faces.Append(orig_mesh.GetBdrFace(1));
+      // faces.Append(orig_mesh.GetBdrFace(2));
+
+   Mesh surface_mesh = Mesh::ExtractSurfaceMesh(orig_mesh,faces);
+   surface_mesh.CheckElementOrientation(true);
+   surface_mesh.CheckBdrElementOrientation(true);
+
+   surface_mesh.Print();
+
+   if (surface_mesh.Dimension() > 1)
+   {
+      char vishost[] = "localhost";
+      int  visport   = 19916;
+
+      socketstream mesh2_sock(vishost, visport);
+      mesh2_sock.precision(8);
+      mesh2_sock << "mesh\n" << surface_mesh << "keys n \n" << flush;
+   }
+   else
+   {
+      ParaViewDataCollection paraview_dc("surf_mesh", &surface_mesh);
+      paraview_dc.SetPrefixPath("ParaView");
+      paraview_dc.SetLevelsOfDetail(3);
+      paraview_dc.SetCycle(0);
+      paraview_dc.SetDataFormat(VTKFormat::BINARY);
+      paraview_dc.SetHighOrderOutput(true);
+      paraview_dc.SetTime(0.0); // set the time
+      H1_FECollection fec(order,surface_mesh.Dimension());
+      FiniteElementSpace fespace(&surface_mesh,&fec);
+      GridFunction gf(&fespace);
+      gf.Randomize();
+      paraview_dc.RegisterField("solution",&gf);
+      paraview_dc.Save();
+   }
+
+   return 0;
+}
+

examples/MeshPart/test2.cpp

@@ -0,0 +1,180 @@
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+#include "mesh_partition.hpp"
+
+using namespace std;
+using namespace mfem;
+
+double sin_func(const Vector & x);
+
+int main(int argc, char *argv[])
+{
+   // 1. Parse command line options
+   const char *mesh_file = "../../data/periodic-annulus-sector.msh";
+   int order = 1;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh", "Mesh file to use.");
+   args.AddOption(&order, "-o", "--order", "Finite element polynomial degree");
+   args.ParseCheck();
+
+   // 2. Read the mesh from the given mesh file, and refine once uniformly.
+   Mesh mesh(mesh_file);
+   // mesh.EnsureNodes();
+   // mesh.UniformRefinement();
+   // Array<int> elems0({0,1,2,3});
+   int nel = mesh.GetNE();
+   // int nel = 5;
+   Array<int> elems0(nel/2);
+   for (int i = 0; i<nel/2; i++)
+   {
+      elems0[i] = i;
+   }
+   elems0.Print();
+   // elems0.Append(24);
+   // elems0.Append(23);
+   // elems0.Append(26);
+   Subdomain subdomain0(mesh);
+
+   Mesh * submesh = subdomain0.GetSubMesh(elems0);
+
+
+   // cout << "number of boundary elements = " << mesh.GetNBE() << endl;
+   Array<int> faces(mesh.GetNBE()/2);
+   for (int i = 0; i<mesh.GetNBE()/2; i++)
+   {
+      faces[i] = mesh.GetBdrFace(i);
+   }
+   
+   Mesh * surfmesh = subdomain0.GetSurfaceMesh(faces);
+
+   H1_FECollection fec(order, mesh.Dimension());
+   FiniteElementSpace fespace(&mesh, &fec);
+   FunctionCoefficient coeff(sin_func);
+   GridFunction gf(&fespace);
+   gf.ProjectCoefficient(coeff);
+   {
+      char vishost[] = "localhost";
+      int  visport   = 19916;
+      socketstream mesh_sock(vishost, visport);
+      mesh_sock.precision(8);
+      // mesh_sock << "mesh\n" << mesh << "keys n \n" << flush;
+      mesh_sock << "solution\n" << mesh << gf <<  "keys jnmR \n" 
+      << "valuerange 0 1.0 \n" << flush;
+         // << flush;
+   }
+
+   subdomain0.SetFESpace(fespace);
+
+   SparseMatrix * P = subdomain0.GetProlonationMatrix();
+   FiniteElementSpace * elem_fes = 
+               subdomain0.GetSubFESpace(Subdomain::entity_type::volume);
+   GridFunction gf_e(elem_fes);
+   cout << "Size P = " << P->Height() << " x " << P->Width() << endl;
+   cout << "gf_e.Size = " << gf_e.Size() << endl;
+   cout << "gf.Size = " << gf.Size() << endl;
+   P->MultTranspose(gf,gf_e);
+
+   SparseMatrix * Pb = subdomain0.GetSurfaceProlonationMatrix();
+   FiniteElementSpace * bdr_elem_fes = 
+               subdomain0.GetSubFESpace(Subdomain::entity_type::surface);
+   GridFunction gf_b(bdr_elem_fes);
+   Pb->MultTranspose(gf,gf_b);
+
+   {
+      char vishost[] = "localhost";
+      int  visport   = 19916;
+      if (submesh)
+      {
+         socketstream mesh0_sock(vishost, visport);
+         mesh0_sock.precision(8);
+         // mesh0_sock << "mesh\n" << *submesh << "keys n \n" << flush;
+         mesh0_sock << "solution\n" << *submesh << gf_e << "keys nmR \n"
+         << "valuerange 0 1.0 \n" << flush;
+         // << flush;
+      }
+      if (surfmesh && mesh.Dimension() == 3)
+      {
+         socketstream mesh1_sock(vishost, visport);
+         mesh1_sock.precision(8);
+         // mesh1_sock << "mesh\n" << *bdrmesh0 << "keys n \n" << flush;
+         mesh1_sock << "solution\n" << *surfmesh << gf_b
+         << "valuerange 0 1.0 \n" << flush;
+         // << flush;
+      }
+   }
+
+   // ParaViewDataCollection paraview_dc("mesh_partition", surfmesh);
+   // paraview_dc.SetPrefixPath("ParaView");
+   // const FiniteElementSpace * fes_ = surfmesh->GetNodalFESpace();
+   // int ord = (fes_) ? fes_->GetOrder(0) : order;
+   // paraview_dc.SetLevelsOfDetail(ord);
+   // paraview_dc.SetCycle(0);
+   // paraview_dc.SetDataFormat(VTKFormat::BINARY);
+   // paraview_dc.SetHighOrderOutput(true);
+   // paraview_dc.SetTime(0.0); // set the time
+   // paraview_dc.RegisterField("solution",&gf_b);
+   // paraview_dc.Save();
+
+
+
+   // // ---------------------------------------------------------
+   // FiniteElementCollection *fec1 = new H1_FECollection(order, submesh->Dimension());
+   // FiniteElementSpace fespace1(submesh, fec1);
+   // Array<int> ess_tdof_list;
+   // if (submesh->bdr_attributes.Size())
+   // {
+   //    Array<int> ess_bdr(mesh.bdr_attributes.Max());
+   //    ess_bdr = 1;
+   //    fespace1.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+   // }
+
+   // LinearForm b(&fespace1);
+   // ConstantCoefficient one(1.0);
+   // b.AddDomainIntegrator(new DomainLFIntegrator(one));
+   // b.Assemble();
+
+   // GridFunction x(&fespace1);
+   // x = 0.0;
+
+   // // 9. Set up the bilinear form a(.,.) on the finite element space
+   // //    corresponding to the Laplacian operator -Delta, by adding the Diffusion
+   // //    domain integrator.
+   // BilinearForm a(&fespace1);
+   // a.AddDomainIntegrator(new DiffusionIntegrator(one));
+   // a.Assemble();
+   // OperatorPtr A;
+   // Vector B, X;
+   // a.FormLinearSystem(ess_tdof_list, x, b, A, X, B);
+   // cout << "Size of linear system: " << A->Height() << endl;
+   // // Use a simple symmetric Gauss-Seidel preconditioner with PCG.
+   // GSSmoother M((SparseMatrix&)(*A));
+   // PCG(*A, M, B, X, 1, 200, 1e-12, 0.0);
+   // // 12. Recover the solution as a finite element grid function.
+   // a.RecoverFEMSolution(X, b, x);
+   // {
+   //    char vishost[] = "localhost";
+   //    int  visport   = 19916;
+   //    socketstream sol_sock2(vishost, visport);
+   //    sol_sock2.precision(8);
+   //    sol_sock2 << "solution\n" << *submesh << x << flush;
+   // }
+
+   return 0;
+}
+
+double sin_func(const Vector & x)
+{
+   Vector c(x.Size());
+   c.Randomize();
+   // double dotp = c*x;
+   // return (sin(10.0*M_PI*dotp));
+   // return sin(2.*M_PI*x[0]);
+   // return 1.-x[1]*x[1]/4.0;
+   // return (0.5-x[1])*(0.5-x[1]);
+   return x[1];
+   // double r = sqrt(x[0]*x[0] + x[1]*x[1] + x[2]*x[2]);
+   // return r;
+}

examples/MeshPart/test_mesh_partition.cpp

@@ -0,0 +1,178 @@
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+#include "mesh_partition.hpp"
+
+using namespace std;
+using namespace mfem;
+
+double sin_func(const Vector & x);
+
+int main(int argc, char *argv[])
+{
+   // 1. Parse command line options
+   const char *mesh_file = "../data/star.mesh";
+   int order = 1;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh", "Mesh file to use.");
+   args.AddOption(&order, "-o", "--order", "Finite element polynomial degree");
+   args.ParseCheck();
+
+   // 2. Read the mesh from the given mesh file, and refine once uniformly.
+   Mesh mesh(mesh_file);
+   // mesh.Print(cout);
+   mesh.UniformRefinement();
+
+   // mesh.EnsureNodes();
+   // Array<int> elems0({0,4,8,12,16});
+   // Array<int> elems0({0,4,8,12,16});
+   // Array<int> elems0({6,7,8});
+   Array<int> elems0({0,1,2,3,4});
+   // Array<int> elems0({7,6,17,20,21,22});
+   // Array<int> elems0({0,1,2,3,4});
+   // Array<int> elems0({104,103,86,109});
+   // elems0.Print(cout, elems0.Size());
+   Subdomain subdomain0(mesh);
+   Mesh * submesh0 = subdomain0.GetSubMesh(elems0);
+
+   // Array<int> bdrelems0({8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23});
+   cout << "number of boundary elements = " << mesh.GetNBE() << endl;
+   Array<int> bdrelems0({0,1});
+   Array<int> faces(bdrelems0.Size());
+   for (int i = 0; i<bdrelems0.Size(); i++)
+   {
+      faces[i] = mesh.GetBdrFace(bdrelems0[i]);
+   }
+   Mesh * surfmesh0 = subdomain0.GetSurfaceMesh(faces);
+
+   H1_FECollection fec(order, mesh.Dimension());
+   FiniteElementSpace fespace(&mesh, &fec);
+   FunctionCoefficient coeff(sin_func);
+   GridFunction gf(&fespace);
+   gf.ProjectCoefficient(coeff);
+   {
+      char vishost[] = "localhost";
+      int  visport   = 19916;
+      socketstream mesh_sock(vishost, visport);
+      mesh_sock.precision(8);
+      mesh_sock << "solution\n" << mesh << gf 
+      << "valuerange -1.0 1.0 \n" << flush;
+   }
+
+   subdomain0.SetFESpace(fespace);
+
+   SparseMatrix * P = subdomain0.GetProlonationMatrix();
+   FiniteElementSpace * elem_fes = 
+               subdomain0.GetSubFESpace(Subdomain::entity_type::volume);
+   GridFunction gf_e(elem_fes);
+   P->MultTranspose(gf,gf_e);
+
+   SparseMatrix * Pf = subdomain0.GetSurfaceProlonationMatrix();
+   FiniteElementSpace * face_elem_fes = 
+               subdomain0.GetSubFESpace(Subdomain::entity_type::surface);
+   GridFunction gf_f(face_elem_fes);
+   Pf->MultTranspose(gf,gf_f);
+   {
+      char vishost[] = "localhost";
+      int  visport   = 19916;
+      if (submesh0)
+      {
+         socketstream mesh0_sock(vishost, visport);
+         mesh0_sock.precision(8);
+         // mesh0_sock << "mesh\n" << *submesh0 << "keys n \n" << flush;
+         mesh0_sock << "solution\n" << *submesh0 << gf_e 
+         << "valuerange -1.0 1.0 \n" << flush;
+      }
+      if (surfmesh0 && mesh.Dimension()==3)
+      {
+         socketstream mesh1_sock(vishost, visport);
+         mesh1_sock.precision(8);
+         // mesh1_sock << "mesh\n" << *bdrmesh0 << "keys n \n" << flush;
+         mesh1_sock << "solution\n" << *surfmesh0 << gf_f 
+         << "valuerange -1.0 1.0 \n" << flush;
+      }
+   }
+
+   // Array<int> bdr_faces;
+   // for (int i =0; i<mesh.GetNBE(); i++)
+   // {
+   //    int attr = mesh.GetBdrAttribute(i);
+   //    if (attr == 4)
+   //    {
+   //       bdr_faces.Append(mesh.GetBdrFace(i));
+   //    }
+   // }
+
+   // H1_FECollection fec(order, mesh.Dimension());
+   // FiniteElementSpace fespace(&mesh, &fec);
+   // FunctionCoefficient coeff(sin_func);
+   // GridFunction gf(&fespace);
+   // gf.ProjectCoefficient(coeff);
+   // {
+   //    char vishost[] = "localhost";
+   //    int  visport   = 19916;
+   //    socketstream mesh_sock(vishost, visport);
+   //    mesh_sock.precision(8);
+   //    // mesh_sock << "mesh\n" << mesh << "keys n \n" << flush;
+   //    mesh_sock << "solution\n" << mesh << gf << flush; 
+   //    // << "valuerange -5000.0 5000.0 \n" << flush;
+   //    // << "valuerange -1.0 1.0 \n" << flush;
+   // }
+
+   // Subdomain subdomain1(mesh);
+   // subdomain1.SetFESpace(fespace);
+   // Mesh * bdrmesh0 = subdomain1.GetBdrSurfaceMesh(bdr_faces);
+   // SparseMatrix * Pb = subdomain1.GetBdrProlonationMatrix();
+   // FiniteElementSpace * bdr_elem_fes = 
+   //             subdomain1.GetSubFESpace(Subdomain::entity_type::bdr);
+   // GridFunction gf_f(bdr_elem_fes);
+   // Pb->MultTranspose(gf,gf_f);
+
+
+
+   // // gf.Print();
+   // // gf_f.Print();
+   // // bdrmesh0->Print(cout);
+
+   // if (bdrmesh0)
+   // {
+   //    char vishost[] = "localhost";
+   //    int  visport   = 19916;
+   //    socketstream mesh1_sock(vishost, visport);
+   //    mesh1_sock.precision(8);
+   //    // mesh1_sock << "mesh\n" << *bdrmesh0 << "keys n \n" << flush;
+   //    mesh1_sock << "solution\n" << *bdrmesh0 << gf_f << flush;
+   //    // << "valuerange -5000.0 5000.0 \n" << flush;
+   // }
+
+   ParaViewDataCollection paraview_dc("mesh_partition", surfmesh0);
+   paraview_dc.SetPrefixPath("ParaView");
+   const FiniteElementSpace * fes_ = surfmesh0->GetNodalFESpace();
+   int ord = (fes_) ? fes_->GetOrder(0) : order;
+   paraview_dc.SetLevelsOfDetail(5);
+   paraview_dc.SetCycle(0);
+   paraview_dc.SetDataFormat(VTKFormat::BINARY);
+   paraview_dc.SetHighOrderOutput(true);
+   paraview_dc.SetTime(0.0); // set the time
+   paraview_dc.RegisterField("solution",&gf_f);
+   paraview_dc.Save();
+
+
+   return 0;
+}
+
+double sin_func(const Vector & x)
+{
+   Vector c(x.Size());
+   c.Randomize();
+   // double dotp = c*x;
+   double dotp = x.Sum();
+   // return (sin(10.0*M_PI*dotp));
+   // return sin(2.*M_PI*x[0]);
+   // return 1.-x[1]*x[1]/4.0;
+   return (0.5-x[1])*(0.5-x[1]);
+   // double r = sqrt(x[0]*x[0] + x[1]*x[1] + x[2]*x[2]);
+   // return r;
+}

mesh/mesh.cpp

@@ -4440,6 +4440,311 @@ void Mesh::MakeSimplicial_(const Mesh &orig_mesh, int *vglobal)
    MFEM_ASSERT(CheckBdrElementOrientation(false) == 0, "");
 }
 
+Mesh Mesh::ExtractMesh(const Mesh &orig_mesh, const Array<int> & elems)
+{
+   Mesh mesh;
+   mesh.ExtractMesh_(orig_mesh, elems);
+   return mesh;
+}
+
+void Mesh::ExtractMesh_(const Mesh &orig_mesh, const Array<int> & elems)
+{
+   int dim = orig_mesh.Dimension();
+   int sdim = orig_mesh.SpaceDimension();
+
+   int nv = orig_mesh.GetNV();
+   int nf = orig_mesh.GetNumFaces();
+   // vertex marker
+   Array<int> vmarker(nv); vmarker = 0;
+   int new_nv = 0;
+   int new_ne = elems.Size();
+
+   // Count and mark the vertices to be added to the new mesh
+   Array<int> vertices;
+   for (int iel=0; iel<new_ne; ++iel)
+   {
+      int el = elems[iel];
+      orig_mesh.GetElementVertices(el,vertices);
+      for (int iv=0; iv<vertices.Size(); ++iv)
+      {
+         int v = vertices[iv];
+         if (vmarker[v]) { continue; }
+         vmarker[v] = 1;
+         new_nv++;
+      }
+   }
+
+   // Count the bdry elements to be added to the new mesh
+   Array<int> BoundaryMarker(nf); BoundaryMarker = 0;
+   int new_nbe = 0;
+   for (int i = 0; i<new_ne; i++)
+   {
+      int el = elems[i];
+      Array<int> faces, ori;
+      switch (dim)
+      {
+         case 1: orig_mesh.GetElementVertices(el,faces); break;
+         case 2: orig_mesh.GetElementEdges(el,faces,ori); break;
+         default: orig_mesh.GetElementFaces(el,faces,ori); break;
+      }
+      for (int f=0; f<faces.Size(); ++f) { BoundaryMarker[faces[f]]++; }
+   }
+
+   for (int f=0; f<nf; ++f) { if (BoundaryMarker[f] == 1) { new_nbe++; } }
+
+   InitMesh(dim,sdim,new_nv,new_ne,new_nbe);
+
+   // Add vertices
+   int vk = 0;
+   for (int iv = 0; iv<nv; ++iv)
+   {
+      if (!vmarker[iv]) { continue; }
+      AddVertex(orig_mesh.GetVertex(iv));
+      // save the vertex unique id to the marker
+      vmarker[iv] = ++vk;
+   }
+
+   // add BdrElements
+   int attr = 0;
+   for (int f=0; f<nf; ++f)
+   {
+      if (BoundaryMarker[f] != 1) { continue; }
+      const Element * el = orig_mesh.GetFace(f);
+      Element * new_bel = nullptr;
+      if (dim == 1)
+      {
+         int pt = vmarker[f]-1;
+         new_bel = (Element*)new Point(&pt);
+         new_bel->SetAttribute(++attr);
+      }
+      else
+      {
+         new_bel = NewElement(el->GetGeometryType());
+         int nv0 = el->GetNVertices();
+         const int * v0 = el->GetVertices();
+         Array<int> v1(nv0);
+         for (int i=0; i<nv0; i++)
+         {
+            v1[i] = vmarker[v0[i]]-1;
+         }
+         new_bel->SetVertices(v1.GetData());
+         new_bel->SetAttribute(++attr);
+      }
+      AddBdrElement(new_bel);
+   }
+
+   // add elements
+   for (int e=0; e<new_ne; ++e)
+   {
+      const Element * el = orig_mesh.GetElement(elems[e]);
+      Element * new_el = NewElement(el->GetGeometryType());
+      int nv0 = el->GetNVertices();
+      const int * v0 = el->GetVertices();
+      Array<int> v1(nv0);
+      for (int i=0; i<nv0; i++) { v1[i] = vmarker[v0[i]]-1; }
+      new_el->SetVertices(v1.GetData());
+      new_el->SetAttribute(el->GetAttribute());
+      AddElement(new_el);
+   }
+
+   FinalizeTopology();
+
+   const GridFunction * orig_nodes = orig_mesh.GetNodes();
+   if (orig_nodes)
+   {
+      const FiniteElementSpace * orig_fes = orig_nodes->FESpace();
+
+      Ordering::Type ordering = orig_fes->GetOrdering();
+      int order = orig_fes->FEColl()->GetOrder();
+      bool discont = orig_fes->IsDGSpace();
+
+      SetCurvature(order, discont, sdim, ordering);
+
+      const FiniteElementSpace * new_fes = GetNodalFESpace();
+      GridFunction * new_nodes = GetNodes();
+
+      Array<int> orig_vdofs;
+      Array<int> new_vdofs;
+      Vector vec;
+
+      // Copy nodes to submesh
+      for (int e = 0; e < new_ne; e++)
+      {
+         new_fes->GetElementVDofs(e, new_vdofs);
+         orig_fes->GetElementVDofs(elems[e], orig_vdofs);
+         orig_nodes->GetSubVector(orig_vdofs, vec);
+         new_nodes->SetSubVector(new_vdofs, vec);
+      }
+   }
+   Finalize();
+}
+
+Mesh Mesh::ExtractSurfaceMesh(const Mesh &orig_mesh, const Array<int> & faces)
+{
+   Mesh mesh;
+   mesh.ExtractSurfaceMesh_(orig_mesh, faces);
+   return mesh;
+}
+
+void Mesh::ExtractSurfaceMesh_(const Mesh &orig_mesh, const Array<int> & faces)
+{
+   int dim = orig_mesh.Dimension()-1;
+
+   MFEM_VERIFY(dim > 0, "Only dim > 1 is supported");
+
+   int sdim = orig_mesh.SpaceDimension();
+
+   int nv = orig_mesh.GetNV();
+   int nf = (dim == 2) ? orig_mesh.GetNEdges() : nv;
+   // vertex marker
+   Array<int> vmarker(nv); vmarker = 0;
+   int new_nv = 0;
+   int new_ne = faces.Size();
+
+   // Count and mark the vertices to be added to the new mesh
+   Array<int> vertices;
+   for (int f=0; f<new_ne; ++f)
+   {
+      int el = faces[f];
+      orig_mesh.GetFaceVertices(el,vertices);
+      for (int iv=0; iv<vertices.Size(); ++iv)
+      {
+         int v = vertices[iv];
+         if (vmarker[v]) { continue; }
+         vmarker[v] = 1;
+         new_nv++;
+      }
+   }
+
+   // Count the bdry elements to be added to the new mesh
+   Array<int> BoundaryMarker(nf); BoundaryMarker = 0;
+   int new_nbe = 0;
+   for (int i = 0; i<new_ne; i++)
+   {
+      int el = faces[i];
+      Array<int> edges, ori;
+      switch (dim)
+      {
+         case 1: orig_mesh.GetEdgeVertices(el,edges); break; // these are vertices
+         case 2: orig_mesh.GetFaceEdges(el,edges,ori); break;
+         default: MFEM_ABORT("Unreachable"); break;
+      }
+      for (int f=0; f<edges.Size(); ++f) { BoundaryMarker[edges[f]]++; }
+   }
+
+   for (int f=0; f<nf; ++f) { if (BoundaryMarker[f] == 1) { new_nbe++; } }
+
+   InitMesh(dim,sdim,new_nv,new_ne,new_nbe);
+
+   // Add vertices
+   int vk = 0;
+   for (int iv = 0; iv<nv; ++iv)
+   {
+      if (!vmarker[iv]) { continue; }
+      AddVertex(orig_mesh.GetVertex(iv));
+      // save the vertex unique id to the marker
+      vmarker[iv] = ++vk;
+   }
+
+   // add BdrElements
+   int attr = 0;
+   for (int f=0; f<nf; ++f)
+   {
+      if (BoundaryMarker[f] != 1) { continue; }
+      Element * new_bel = nullptr;
+      if (dim == 1)
+      {
+         int pt = vmarker[f]-1;
+         new_bel = (Element*)new Point(&pt);
+         new_bel->SetAttribute(++attr);
+      }
+      else
+      {
+         new_bel = NewElement(mfem::Geometry::SEGMENT);
+         int nv0 = 2;
+         Array<int> vert;
+         orig_mesh.GetEdgeVertices(f,vert);
+         Array<int> v1(nv0);
+         for (int i=0; i<nv0; i++)
+         {
+            v1[i] = vmarker[vert[i]]-1;
+         }
+         new_bel->SetVertices(v1.GetData());
+         new_bel->SetAttribute(++attr);
+      }
+      AddBdrElement(new_bel);
+   }
+
+   // add elements
+   for (int e=0; e<new_ne; ++e)
+   {
+      const Element * el = orig_mesh.GetFace(faces[e]);
+      Element * new_el = NewElement(el->GetGeometryType());
+      int nv0 = el->GetNVertices();
+      const int * v0 = el->GetVertices();
+      Array<int> v1(nv0);
+      for (int i=0; i<nv0; i++) { v1[i] = vmarker[v0[i]]-1; }
+      new_el->SetVertices(v1.GetData());
+      new_el->SetAttribute(el->GetAttribute());
+      AddElement(new_el);
+   }
+
+   FinalizeTopology();
+
+   const GridFunction * orig_nodes = orig_mesh.GetNodes();
+   if (orig_nodes)
+   {
+      const FiniteElementSpace * orig_fes = orig_nodes->FESpace();
+
+      Ordering::Type ordering = orig_fes->GetOrdering();
+      int order = orig_fes->FEColl()->GetOrder();
+      bool discont = orig_fes->IsDGSpace();
+
+      SetCurvature(order, discont, sdim, ordering);
+
+      const FiniteElementSpace * new_fes = GetNodalFESpace();
+      GridFunction * new_nodes = GetNodes();
+
+      Array<int> orig_vdofs;
+      Array<int> new_vdofs;
+      Vector vec;
+
+      // Copy nodes to submesh
+      for (int e = 0; e < new_ne; e++)
+      {
+         new_fes->GetElementVDofs(e, new_vdofs);
+         if (!discont)
+         {
+            orig_fes->GetFaceVDofs(faces[e], orig_vdofs);
+            orig_nodes->GetSubVector(orig_vdofs, vec);
+         }
+         else
+         {
+            const FiniteElement * el = new_fes->GetFE(e);
+            const IntegrationRule & ir = el->GetNodes();
+            int np = ir.GetNPoints();
+            FaceElementTransformations * Tr =
+               const_cast<Mesh *>(&orig_mesh)->GetFaceElementTransformations(faces[e]);
+            int el1 = Tr->Elem1No;
+            vec.SetSize(new_vdofs.Size());
+            for (int i = 0; i<np; i++)
+            {
+               Tr->SetAllIntPoints(&ir[i]);
+               const IntegrationPoint & ip = Tr->GetElement1IntPoint();
+               Vector val;
+               orig_nodes->GetVectorValue(el1,ip,val);
+               for (int j = 0; j<val.Size(); j++)
+               {
+                  vec[i+j*np] = val[j];
+               }
+            }
+         }
+         new_nodes->SetSubVector(new_vdofs, vec);
+      }
+   }
+   Finalize();
+}
+
 Mesh Mesh::MakePeriodic(const Mesh &orig_mesh, const std::vector<int> &v2v)
 {
    Mesh periodic_mesh(orig_mesh, true); // Make a copy of the original mesh

mesh/mesh.hpp

@@ -486,6 +486,12 @@ protected:
    // Internal helper used in MakeSimplicial (and ParMesh::MakeSimplicial).
    void MakeSimplicial_(const Mesh &orig_mesh, int *vglobal);
 
+   /// Internal helper used in ExtractMesh
+   void ExtractMesh_(const Mesh &orig_mesh, const Array<int> & elems);
+
+   /// Internal helper used in ExtractMesh
+   void ExtractSurfaceMesh_(const Mesh &orig_mesh, const Array<int> & faces);
+
 public:
 
    Mesh() { SetEmpty(); }
@@ -579,6 +585,16 @@ public:
        new mesh. Periodic meshes are not supported. */
    static Mesh MakeSimplicial(const Mesh &orig_mesh);
 
+   /** Create a mesh by extracting from @a orig_mesh the given list of
+       @a elems. Periodic and high-order meshes are also supported.
+       @warning NonConforming meshes are not supported. */
+   static Mesh ExtractMesh(const Mesh &orig_mesh, const Array<int> & elems);
+
+   /** Create a surface mesh by extracting from @a orig_mesh the given list of
+       @a faces. Periodic and high-order meshes are also supported.
+       @warning NonConforming meshes are not supported. */
+   static Mesh ExtractSurfaceMesh(const Mesh &orig_mesh, const Array<int> & faces);
+
    /// Create a periodic mesh by identifying vertices of @a orig_mesh.
    /** Each vertex @a i will be mapped to vertex @a v2v[i], such that all
        vertices that are coincident under the periodic mapping get mapped to