diff --git a/cmake/FindParmetis.cmake b/cmake/FindParmetis.cmake
index b00a11d41..e9a544418 100644
--- a/cmake/FindParmetis.cmake
+++ b/cmake/FindParmetis.cmake
@@ -28,7 +28,13 @@ if(NOT EXISTS "${METIS_LIBRARY}")
   message(FATAL_ERROR "metis library not found")
 endif()
 
-set(PARMETIS_LIBRARIES ${PARMETIS_LIBRARY} ${METIS_LIBRARY})
+find_library(GK_LIBRARY GKlib PATHS "${PARMETIS_PREFIX}/lib")
+if(EXISTS "${GK_LIBRARY}")
+	set(PARMETIS_LIBRARIES ${PARMETIS_LIBRARY} ${METIS_LIBRARY} ${GK_LIBRARY})
+else()
+	set(PARMETIS_LIBRARIES ${PARMETIS_LIBRARY} ${METIS_LIBRARY})
+endif()
+
 set(PARMETIS_INCLUDE_DIRS ${PARMETIS_INCLUDE_DIR} ${METIS_INCLUDE_DIR})
 
 include(FindPackageHandleStandardArgs)
diff --git a/mds/apfMDS.h b/mds/apfMDS.h
index ef745e576..fbc8749d8 100644
--- a/mds/apfMDS.h
+++ b/mds/apfMDS.h
@@ -192,8 +192,12 @@ int getMdsIndex(Mesh2* in, MeshEntity* e);
   so call apf::reorderMdsMesh after any mesh modification. */
 MeshEntity* getMdsEntity(Mesh2* in, int dimension, int index);
 
+int gmshMajorVersion(const char* filename);
+
 Mesh2* loadMdsFromGmsh(gmi_model* g, const char* filename);
 
+Mesh2* loadMdsDmgFromGmsh(const char* fnameDmg, const char* filename);
+
 Mesh2* loadMdsFromUgrid(gmi_model* g, const char* filename);
 
 void printUgridPtnStats(gmi_model* g, const char* ugridfile, const char* ptnfile,
diff --git a/mds/mdsGmsh.cc b/mds/mdsGmsh.cc
index 63e5ce463..7d4a26567 100644
--- a/mds/mdsGmsh.cc
+++ b/mds/mdsGmsh.cc
@@ -52,6 +52,8 @@ struct Reader {
   char* line;
   char* word;
   size_t linecap;
+  int major_version;
+  int minor_version;
   bool isQuadratic;
   std::map<long, Node> nodeMap;
   std::map<long, apf::MeshEntity*> entMap[4];
@@ -60,19 +62,6 @@ struct Reader {
   std::vector<int> physicalType[4];
 };
 
-void initReader(Reader* r, apf::Mesh2* m, const char* filename)
-{
-  r->mesh = m;
-  r->file = fopen(filename, "r");
-  if (!r->file) {
-    lion_eprint(1,"couldn't open Gmsh file \"%s\"\n",filename);
-    abort();
-  }
-  r->line = static_cast<char*>(malloc(1));
-  r->line[0] = '\0';
-  r->linecap = 1;
-  r->isQuadratic = false;
-}
 
 void freeReader(Reader* r)
 {
@@ -97,6 +86,16 @@ long getLong(Reader* r)
   return x;
 }
 
+double getDouble(Reader* r)
+{
+  double x;
+  int nchars;
+  int ret = sscanf(r->word, "%lf%n", &x, &nchars);
+  PCU_ALWAYS_ASSERT(ret == 1);
+  r->word += nchars;
+  return x;
+}
+
 bool startsWith(char const* prefix, char const* s)
 {
   int ls = strlen(s);
@@ -118,18 +117,114 @@ void checkMarker(Reader* r, char const* marker)
   PCU_ALWAYS_ASSERT(startsWith(marker, r->line));
 }
 
-void readNode(Reader* r)
+void initReader(Reader* r, apf::Mesh2* m, const char* filename)
+{
+  r->mesh = m;
+  r->file = fopen(filename, "r");
+  if (!r->file) {
+    lion_eprint(1,"couldn't open Gmsh file \"%s\"\n",filename);
+    abort();
+  }
+  r->line = static_cast<char*>(malloc(1));
+  r->line[0] = '\0';
+  r->linecap = 1;
+  r->isQuadratic = false;
+  seekMarker(r, "$MeshFormat");
+  int fileType, dataSize;
+  int ret = sscanf(r->line, "%d.%d %d %d\n",
+      &r->major_version, &r->minor_version, &fileType, &dataSize);
+  PCU_ALWAYS_ASSERT(ret==4);
+}
+
+void readNode(Reader* r, int bm=-1)
 {
-  long id;
   Node n;
   apf::Vector3& p = n.point;
-  sscanf(r->line, "%ld %lf %lf %lf", &id, &p[0], &p[1], &p[2]);
-  r->nodeMap[id] = n;
+  if(r->major_version == 4) {
+    sscanf(r->line, "%lf %lf %lf", &p[0], &p[1], &p[2]);
+    r->nodeMap[bm] = n;
+  } else if(r->major_version == 2) {
+    long id;
+    sscanf(r->line, "%ld %lf %lf %lf", &id, &p[0], &p[1], &p[2]);
+    r->nodeMap[id] = n;
+  }
   getLine(r);
 }
 
-void readNodes(Reader* r)
+void readEntities(Reader* r,const char* fnameDmg) 
 {
+  seekMarker(r, "$Entities");
+  long nlde,ilde,iud,tag,isign,nMV,nME,nMF,nMR;
+  double x,y,z;
+  FILE* f = fopen(fnameDmg, "w");
+  sscanf(r->line, "%ld %ld %ld %ld", &nMV, &nME, &nMF, &nMR);
+  fprintf(f, "%ld %ld %ld %ld \n", nMR, nMF, nME, nMV); // just reverse order 
+  fprintf(f, "%f %f %f \n ", 0.0, 0.0, 0.0); // Probaby model bounding box?
+  fprintf(f, "%f %f %f \n", 0.0, 0.0, 0.0); // 
+   
+  getLine(r); // because readNode gets the next line we need this outside  for Nodes_Block
+  for (long i = 0; i < nMV; ++i){  
+    sscanf(r->line, "%ld %lf %lf %lf %ld ", &tag, &x, &y, &z, &iud);
+    fprintf(f, "%ld %lf %lf %lf \n",tag,x,y,z);
+    getLine(r); 
+  }
+  for (long i = 0; i < nME; ++i){
+    tag = getLong(r); 
+    fprintf(f, "%ld", tag);
+    for (int i=0; i< 6; ++i) x=getDouble(r);  // read past min maxes
+    iud = getLong(r); 
+    for(long j =0; j < iud; ++j) isign=getLong(r); // read past iud user tags
+    nlde=getLong(r);  // 2 in straight edged models but...
+    for(long j =0; j < nlde; ++j) {
+      ilde=getLong(r);
+      fprintf(f, " %ld", std::abs(ilde)); // modVerts started from 1
+    }
+    fprintf(f, "\n");
+    getLine(r); 
+  }
+  for (long i = 0; i < nMF; ++i){
+    tag = getLong(r); 
+    fprintf(f, "%ld %d\n", tag, 1);
+    for (int i=0; i< 6; ++i) x=getDouble(r);  // read past min maxes
+    iud = getLong(r); 
+    for(long j =0; j < iud; ++j) isign=getLong(r); // read past iud user tags
+    nlde=getLong(r); 
+    fprintf(f, "  %ld \n", nlde);
+    for(long j =0; j < nlde; ++j) {
+      ilde=getLong(r); 
+      if(ilde > 0 ) 
+        isign=1;
+      else
+        isign=0;
+      fprintf(f, "    %ld %ld \n", std::abs(ilde),isign); 
+    }
+    getLine(r); 
+  }   
+  for (long i = 0; i < nMR; ++i){ 
+    tag = getLong(r); 
+    fprintf(f, "%ld %d \n", tag, 1);
+    for (int i=0; i< 6; ++i) x=getDouble(r);  // read past min maxes
+    iud = getLong(r); 
+    for(long j =0; j < iud; ++j) getLong(r); // read past iud user tags
+    nlde=getLong(r); 
+    fprintf(f, "%ld \n", nlde);
+    for(long j =0; j < nlde; ++j) {
+      ilde=getLong(r); 
+      if(ilde > 0 ) 
+        isign=1;
+      else
+        isign=0;
+      fprintf(f, "%ld %ld \n", std::abs(ilde),isign); 
+    }
+    getLine(r); 
+  }   
+  checkMarker(r, "$EndEntities");
+  fclose(f);
+}
+
+void readNodesV2(Reader* r)
+{
+  PCU_ALWAYS_ASSERT(r->major_version == 2);
   seekMarker(r, "$Nodes");
   long n = getLong(r);
   getLine(r);
@@ -138,6 +233,27 @@ void readNodes(Reader* r)
   checkMarker(r, "$EndNodes");
 }
 
+void readNodesV4(Reader* r)
+{
+  PCU_ALWAYS_ASSERT(r->major_version == 4);
+  seekMarker(r, "$Nodes");
+  long Num_EntityBlocks,Num_Nodes,Nodes_Block,edim,etag,junk1,junk2,junk3;
+  sscanf(r->line, "%ld %ld %ld %ld", &Num_EntityBlocks, &Num_Nodes, &junk1, &junk2);
+  getLine(r); // because readNode gets the next line we need this outside  for Nodes_Block
+  for (long i = 0; i < Num_EntityBlocks; ++i){
+    sscanf(r->line, "%ld %ld %ld %ld", &edim, &etag, &junk3, &Nodes_Block);
+    long blockMap[Nodes_Block];
+    for (long j = 0; j < Nodes_Block; ++j){
+      getLine(r);
+      sscanf(r->line, "%ld", &blockMap[j]);
+    }
+    getLine(r);
+    for (long j = 0; j < Nodes_Block; ++j)
+      readNode(r,blockMap[j]);  // has a genLine at end
+  }    
+  checkMarker(r, "$EndNodes");
+}
+
 apf::MeshEntity* lookupVert(Reader* r, long nodeId, apf::ModelEntity* g)
 {
   PCU_ALWAYS_ASSERT(r->nodeMap.count(nodeId));
@@ -149,34 +265,38 @@ apf::MeshEntity* lookupVert(Reader* r, long nodeId, apf::ModelEntity* g)
   return n.entity;
 }
 
-void readElement(Reader* r)
+void readElement(Reader* r, long gmshType=-1, long gtag=-1)
 {
   long id = getLong(r);
-  long gmshType = getLong(r);
+  if(r->major_version == 2) {
+    gmshType = getLong(r);
+  }
   if (isQuadratic(gmshType))
     r->isQuadratic = true;
   int apfType = apfFromGmsh(gmshType);
   PCU_ALWAYS_ASSERT(0 <= apfType);
   int nverts = apf::Mesh::adjacentCount[apfType][0];
   int dim = apf::Mesh::typeDimension[apfType];
-  long ntags = getLong(r);
-  /* The Gmsh 4.9 documentation on the legacy 2.* format states:
-   * "By default, the first tag is the tag of the physical entity to which the
-   * element belongs; the second is the tag of the elementary model entity to
-   * which the element belongs; the third is the number of mesh partitions to
-   * which the element belongs, followed by the partition ids (negative
-   * partition ids indicate ghost cells). A zero tag is equivalent to no tag.
-   * Gmsh and most codes using the MSH 2 format require at least the first two
-   * tags (physical and elementary tags)."
-   * A physical entity is a user defined grouping of elementary model entities.
-   * An elementary model entity is a geometric model entity. */
-  PCU_ALWAYS_ASSERT(ntags >= 2);
-  const int physType = static_cast<int>(getLong(r));
-  PCU_ALWAYS_ASSERT(dim>=0 && dim<4);
-  r->physicalType[dim].push_back(physType);
-  long gtag = getLong(r);
-  for (long i = 2; i < ntags; ++i)
-    getLong(r); /* discard all other element tags */
+  if(r->major_version == 2) {
+    long ntags = getLong(r);
+    /* The Gmsh 4.9 documentation on the legacy 2.* format states:
+     * "By default, the first tag is the tag of the physical entity to which the
+     * element belongs; the second is the tag of the elementary model entity to
+     * which the element belongs; the third is the number of mesh partitions to
+     * which the element belongs, followed by the partition ids (negative
+     * partition ids indicate ghost cells). A zero tag is equivalent to no tag.
+     * Gmsh and most codes using the MSH 2 format require at least the first two
+     * tags (physical and elementary tags)."
+     * A physical entity is a user defined grouping of elementary model entities.
+     * An elementary model entity is a geometric model entity. */
+    PCU_ALWAYS_ASSERT(ntags >= 2);
+    const int physType = static_cast<int>(getLong(r));
+    PCU_ALWAYS_ASSERT(dim>=0 && dim<4);
+    r->physicalType[dim].push_back(physType);
+    gtag = getLong(r);
+    for (long i = 2; i < ntags; ++i)
+      getLong(r); /* discard all other element tags */
+  }
   apf::ModelEntity* g = r->mesh->findModelEntity(dim, gtag);
   apf::Downward verts;
   for (int i = 0; i < nverts; ++i) {
@@ -193,8 +313,9 @@ void readElement(Reader* r)
   getLine(r);
 }
 
-void readElements(Reader* r)
+void readElementsV2(Reader* r)
 {
+  PCU_ALWAYS_ASSERT(r->major_version == 2);
   seekMarker(r, "$Elements");
   long n = getLong(r);
   getLine(r);
@@ -203,6 +324,23 @@ void readElements(Reader* r)
   checkMarker(r, "$EndElements");
 }
 
+void readElementsV4(Reader* r)
+{
+  PCU_ALWAYS_ASSERT(r->major_version == 4);
+  seekMarker(r, "$Elements");
+  long Num_EntityBlocks,Num_Elements,Elements_Block,Edim,gtag,gmshType,junk1,junk2;
+  sscanf(r->line, "%ld %ld %ld %ld", &Num_EntityBlocks, &Num_Elements, &junk1, &junk2);
+  getLine(r);
+  for (long i = 0; i < Num_EntityBlocks; ++i){
+    sscanf(r->line, "%ld %ld %ld %ld", &Edim, &gtag, &gmshType, &Elements_Block);
+    getLine(r);
+    for (long j = 0; j < Elements_Block; ++j) {
+       readElement(r,gmshType,gtag);
+    }
+  }
+  checkMarker(r, "$EndElements");
+}
+
 void setElmPhysicalType(Reader* r, apf::Mesh2* m) {
   apf::MeshEntity* e;
   apf::MeshTag* tag = m->createIntTag("gmsh_physical_entity", 1);
@@ -281,18 +419,44 @@ void readGmsh(apf::Mesh2* m, const char* filename)
 {
   Reader r;
   initReader(&r, m, filename);
-  readNodes(&r);
-  readElements(&r);
-  m->acceptChanges();
-  setElmPhysicalType(&r,m);
-  freeReader(&r);
+  if(r.major_version == 4) {
+    readNodesV4(&r);
+    readElementsV4(&r);
+    m->acceptChanges();
+  } else if(r.major_version == 2) {
+    readNodesV2(&r);
+    readElementsV2(&r);
+    m->acceptChanges();
+    setElmPhysicalType(&r,m);
+  }
   if (r.isQuadratic)
     readQuadratic(&r, m, filename);
+  freeReader(&r);
 }
+}  // closes original namespace 
+
+namespace apf {
 
+int gmshMajorVersion(const char* filename) {
+  Reader r;
+  Mesh2* m=NULL;
+  initReader(&r, m,  filename);
+  int version = r.major_version;
+  freeReader(&r);
+  return version;
+}
+
+void gmshFindDmg(const char* fnameDmg, const char* filename)
+{
+  Reader r;
+  
+  Mesh2* m=NULL;
+  initReader(&r, m,  filename);
+  PCU_ALWAYS_ASSERT(r.major_version==4);
+  readEntities(&r, fnameDmg);
+  freeReader(&r);
 }
 
-namespace apf {
 
 Mesh2* loadMdsFromGmsh(gmi_model* g, const char* filename)
 {
@@ -301,4 +465,12 @@ Mesh2* loadMdsFromGmsh(gmi_model* g, const char* filename)
   return m;
 }
 
+Mesh2* loadMdsDmgFromGmsh(const char*fnameDmg, const char* filename)
+{
+  gmshFindDmg(fnameDmg, filename);  // new function that scans $Entities and writes a dmg 
+  Mesh2* m = makeEmptyMdsMesh(gmi_load(fnameDmg), 0, false);
+  readGmsh(m, filename);
+  return m;
+}
+
 }
diff --git a/pumi-meshes b/pumi-meshes
index 98f48ca2c..d03989485 160000
--- a/pumi-meshes
+++ b/pumi-meshes
@@ -1 +1 @@
-Subproject commit 98f48ca2cb8f1f3ad3db0c5fe4ad1c8df1f20811
+Subproject commit d0398948585a3a0487f0e6cc61fe591573531fc4
diff --git a/test/gmsh.cc b/test/gmsh.cc
index 307fab07e..b2f3911ec 100644
--- a/test/gmsh.cc
+++ b/test/gmsh.cc
@@ -13,20 +13,45 @@ int main(int argc, char** argv)
   MPI_Init(&argc,&argv);
   PCU_Comm_Init();
   lion_set_verbosity(1);
-  if ( argc != 4 ) {
+  if ( argc != 5 ) {
     if ( !PCU_Comm_Self() )
-      printf("Usage: %s <model> <in .msh> <out .smb>\n", argv[0]);
+      printf("Usage: %s <in .dmg> <in .msh> <out .smb> <out .dmg>\n"
+             "The input .msh and output .smb file names are required. \n"
+             "If 'none' is specified as the input model file name then \n"
+             "a output model (.dmg) will be written to the specified filename. \n"
+             "When a **gmsh v2** .msh is passed in, a minimal model will be created from "
+             "the mesh.\n"
+             "When a **gmsh v4** .msh is passed in, a topological model will be created "
+             "from the geometric model entities defined in the gmsh input file.\n", argv[0]);
     MPI_Finalize();
     exit(EXIT_FAILURE);
   }
   gmi_register_null();
   gmi_register_mesh();
-  apf::Mesh2* m = apf::loadMdsFromGmsh(gmi_load(argv[1]), argv[2]);
-  // if input model is null derive a basic model for verify to pass.
-  if (std::string(argv[1]).compare(".null") == 0)
-    apf::deriveMdsModel(m);
+
+  std::string model(argv[1]);
+  std::string gmsh(argv[2]);
+  std::string outMesh(argv[3]);
+  std::string outModel(argv[4]);
+
+  const int gmshVersion = apf::gmshMajorVersion(gmsh.c_str());
+  fprintf(stderr, "version %d\n", gmshVersion);
+  apf::Mesh2* m = NULL;
+  if (gmshVersion == 2) {
+    if (model.compare("none") == 0) {
+      m = apf::loadMdsFromGmsh(gmi_load(".null"), gmsh.c_str());
+      apf::deriveMdsModel(m);
+      gmi_write_dmg(m->getModel(),outModel.c_str());
+    } else {
+      m = apf::loadMdsFromGmsh(gmi_load(model.c_str()), gmsh.c_str());
+    }
+  } else if (gmshVersion == 4) {
+    if (model.compare("none") == 0) {
+      m = apf::loadMdsDmgFromGmsh(outModel.c_str(), gmsh.c_str());
+    }
+  }
   m->verify();
-  m->writeNative(argv[3]);
+  m->writeNative(outMesh.c_str());
   m->destroyNative();
   apf::destroyMesh(m);
   PCU_Comm_Free();
diff --git a/test/testing.cmake b/test/testing.cmake
index de1120f82..b9d7011e7 100644
--- a/test/testing.cmake
+++ b/test/testing.cmake
@@ -191,11 +191,23 @@ mpi_test(create_misSquare 1
   mis_test)
 
 set(MDIR ${MESHES}/gmsh)
-mpi_test(twoQuads 1
+mpi_test(gmshv2TwoQuads 1
   ./from_gmsh
-  ".null"
+  "none"
   "${MDIR}/twoQuads.msh"
-  "${MDIR}/twoQuads.smb")
+  "${MDIR}/twoQuads.smb"
+  "${MDIR}/twoQuads.dmg")
+
+set(MDIR ${MESHES}/gmsh/v4)
+mpi_test(gmshV4AirFoil 1
+  ./from_gmsh
+  "none"
+  "${MDIR}/AirfoilDemo.msh"
+  "${MDIR}/AirfoilDemo.smb"
+  "${MDIR}/AirfoilDemo.dmg")
+add_test(NAME gmshV4AirFoil_dmgDiff
+  COMMAND diff -r ${MDIR}/AirfoilDemo.dmg AirfoilDemo_gold.dmg
+  WORKING_DIRECTORY ${MDIR})
 
 set(MDIR ${MESHES}/ugrid)
 mpi_test(naca_ugrid 2