Commits

Matt Knepley  committed 037a463 Merge

Merge branch 'knepley/feature-plex-gmsh'

* knepley/feature-plex-gmsh:
DMPlex: Fix for orientation in Gmsh files - Fixed test output
DMPlex GMSH: fix MPIU_INT
DMPlex: Slience bogus compiler warning
DMPlex Gmsh: Fixed integer size problem
DMPlex Gmsh: Fixed interface for coordinates
DMPlex ex1: Added test for Gmsh
DMPlex: Added Gmsh input

  • Participants
  • Parent commits b009755, 7083e7f

Comments (5)

  1. Jed Brown

    This still has DMPlexGetCoordinateSection crap in it. You really need to build when you merge. 'master' is currently broken --with-cgns, along with several examples.

    $ git grep DMPlexGetCoordinateSection
    src/dm/impls/plex/examples/tests/ex8.c:  ierr = DMPlexGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr);
    src/dm/impls/plex/plexcgns.c:  ierr = DMPlexGetCoordinateSection(*dm, &coordSection);CHKERRQ(ierr);
    src/dm/impls/plex/plexinterpolate.c:.seealso: DMCopyLabels(), DMGetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM(), DMPlexGetCoordinateSection()
    src/dm/impls/plex/plexinterpolate.c:.seealso: DMCopyCoordinates(), DMGetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM(), DMPlexGetCoordinateSection()
    src/snes/examples/tutorials/ex31.c:  ierr = DMPlexGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr);
    src/snes/examples/tutorials/ex72.c:  ierr = DMPlexGetCoordinateSection(dm, &cSection);CHKERRQ(ierr);
    src/ts/examples/tutorials/ex11.c:  ierr = DMPlexGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr);
    src/ts/examples/tutorials/ex11.c:  ierr = DMPlexGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr);
    src/ts/examples/tutorials/ex11.c:  ierr = DMPlexGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr);
    src/ts/examples/tutorials/ex11.c:  ierr = DMPlexGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr);
    
    1. Jed Brown

      c2166f769c2210509fa086b82dd4549a2bd4a814 in 'knepley/feature-dmda-section', which has not been merged because it does other things. This was broken in the branch and broken in 'master' after the merge. You are supposed to check that the branch-to-merge is functional on its own (otherwise how do you test those changes while developing and how can you bisect later?) and is functional after you merge. Testing an entirely different branch ('next') and inferring that anything contained there is also okay significantly weakens the model.

Files changed (8)

File config/builder.py

                                                                  {'numProcs': 2, 'args': '-dim 2 -cell_simplex 0 -refinement_uniform 1 -interpolate 1 -dm_view ::ascii_latex'},
                                                                  # CGNS tests 10-11 (need to find smaller test meshes)
                                                                  {'numProcs': 1, 'args': '-filename %(meshes)s/tut21.cgns -interpolate 1 -dm_view', 'requires': ['CGNS']},
-                                                                 {'numProcs': 1, 'args': '-filename %(meshes)s/StaticMixer.cgns -interpolate 1 -dm_view', 'requires': ['CGNS']}],
+                                                                 {'numProcs': 1, 'args': '-filename %(meshes)s/StaticMixer.cgns -interpolate 1 -dm_view', 'requires': ['CGNS']},
+                                                                 # Gmsh tests 12-
+                                                                 {'numProcs': 1, 'args': '-filename %(meshes)s/doublet-tet.msh -interpolate 1 -dm_view'}],
                         'src/dm/impls/plex/examples/tests/ex3': [# 0-2 2D P_1 on a triangle
                                                                  {'numProcs': 1, 'args': '-petscspace_order 1 -num_comp 2 -qorder 1'},
                                                                  {'numProcs': 1, 'args': '-petscspace_order 1 -num_comp 2 -qorder 1 -porder 1'},
                                                                      {'numProcs': 1, 'args': '-dim 3'},],
                         'src/dm/impls/plex/examples/tutorials/ex1f90': [{'numProcs': 1, 'args': ''},
                                                                         {'numProcs': 1, 'args': '-dim 3'},],
+                        'src/dm/impls/plex/examples/tutorials/ex2': [# CGNS meshes 0-1
+                                                                     {'numProcs': 1, 'args': '-filename %(meshes)s/tut21.cgns -interpolate 1', 'requires': ['CGNS']},
+                                                                     {'numProcs': 1, 'args': '-filename %(meshes)s/StaticMixer.cgns -interpolate 1', 'requires': ['CGNS']},
+                                                                     # Gmsh meshes 2-2
+                                                                     {'numProcs': 1, 'args': '-filename %(meshes)s/doublet-tet.msh -interpolate 1'},
+                                                                     # Exodus meshes 3-7
+                                                                     {'numProcs': 1, 'args': '-filename %(meshes)s/sevenside-quad.exo -interpolate 1', 'requires': ['exodusii']},
+                                                                     {'numProcs': 1, 'args': '-filename %(meshes)s/sevenside-quad-15.exo -interpolate 1', 'requires': ['exodusii']},
+                                                                     {'numProcs': 1, 'args': '-filename %(meshes)s/squaremotor-30.exo -interpolate 1', 'requires': ['exodusii']},
+                                                                     {'numProcs': 1, 'args': '-filename %(meshes)s/blockcylinder-50.exo -interpolate 1', 'requires': ['exodusii']},
+                                                                     {'numProcs': 1, 'args': '-filename %(meshes)s/simpleblock-100.exo -interpolate 1', 'requires': ['exodusii']},
+                                                                     ],
                         'src/snes/examples/tutorials/ex12':   [# 2D serial P1 test 0-4
                                                                {'numProcs': 1, 'args': '-run_type test -refinement_limit 0.0    -bc_type dirichlet -interpolate 0 -petscspace_order 1 -show_initial -dm_plex_print_fem 1'},
                                                                {'numProcs': 1, 'args': '-run_type test -refinement_limit 0.0    -bc_type dirichlet -interpolate 1 -petscspace_order 1 -show_initial -dm_plex_print_fem 1'},
                                                                {'numProcs': 1, 'args': '-run_type test -dim 3 -refinement_limit 0.0125 -variable_coefficient field    -interpolate 1 -petscspace_order 2 -show_initial -dm_plex_print_fem 1'},
                                                                # Using ExodusII mesh 37-38 BROKEN
                                                                {'numProcs': 1, 'args': '-run_type test -f %(meshes)s/sevenside.exo -refinement_limit 0.0    -bc_type dirichlet -interpolate 1 -show_initial -dm_plex_print_fem 1 -dm_view',
-                                                                'requires': ['Exodus', 'Broken']},
+                                                                'requires': ['exodusii', 'Broken']},
                                                                {'numProcs': 1, 'args': '-run_type test -dim 3 -f /Users/knepley/Downloads/kis_modell_tet2.exo -refinement_limit 0.0    -bc_type dirichlet -interpolate 1 -show_initial -dm_plex_print_fem 1 -dm_view',
-                                                                'requires': ['Exodus', 'Broken']}],
+                                                                'requires': ['exodusii', 'Broken']}],
                         'src/snes/examples/tutorials/ex31':   [# Decoupled field Dirichlet tests 0-6
                                                                {'numProcs': 1, 'args': '-run_type test -refinement_limit 0.0     -forcing_type constant -bc_type dirichlet -interpolate 1 -show_initial -dm_plex_print_fem 1',
                                                                 'setup': './bin/pythonscripts/PetscGenerateFEMQuadrature.py 2 2 2 1 laplacian 2 1 1 1 gradient 2 1 1 1 identity src/snes/examples/tutorials/ex31.h'},

File include/petscdmplex.h

 
 PETSC_EXTERN PetscErrorCode DMPlexCreateExodus(MPI_Comm, PetscInt, PetscBool, DM *);
 PETSC_EXTERN PetscErrorCode DMPlexCreateCGNS(MPI_Comm, PetscInt, PetscBool, DM *);
+PETSC_EXTERN PetscErrorCode DMPlexCreateGmsh(MPI_Comm, PetscViewer, PetscBool, DM *);
 
 PETSC_EXTERN PetscErrorCode DMPlexConstructGhostCells(DM, const char [], PetscInt *, DM *);
 PETSC_EXTERN PetscErrorCode DMPlexConstructCohesiveCells(DM, DMLabel, DM *);

File share/petsc/datafiles/meshes/doublet-tet.msh

+$MeshFormat
+2.2 0 8
+$EndMeshFormat
+$Nodes
+5
+1 -1.0  0.0  0.0
+2  0.0 -1.0  0.0
+3  0.0  0.0 -1.0
+4  0.0  1.0  0.0
+5  1.0  0.0  0.0
+$EndNodes
+$Elements
+2
+1 4 0 3 2 4 1
+2 4 0 4 2 3 5
+$EndElements

File src/dm/impls/plex/examples/tests/ex1.c

   ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr);
   ierr = PetscStrlen(filename, &len);CHKERRQ(ierr);
   if (len) {
+    const char *extGmsh = ".msh";
+    size_t      slen;
+    PetscBool   isGmsh;
+
+    ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extGmsh, 4, &isGmsh);CHKERRQ(ierr);
+    if (isGmsh) {
+      PetscViewer viewer;
+
+      ierr = PetscViewerCreate(comm, &viewer);CHKERRQ(ierr);
+      ierr = PetscViewerSetType(viewer, PETSCVIEWERASCII);CHKERRQ(ierr);
+      ierr = PetscViewerFileSetMode(viewer, FILE_MODE_READ);CHKERRQ(ierr);
+      ierr = PetscViewerFileSetName(viewer, filename);CHKERRQ(ierr);
+      ierr = DMPlexCreateGmsh(comm, viewer, interpolate, dm);CHKERRQ(ierr);
+      ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
+    } else {
 #if defined(PETSC_HAVE_CGNS)
-    int cgid = -1;
+      int cgid = -1;
 
-    if (!rank) {
-      ierr = cg_open(filename, CG_MODE_READ, &cgid);CHKERRQ(ierr);
-      if (cgid <= 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_LIB, "cg_open(\"%s\",...) did not return a valid file ID", filename);
-    }
-    ierr = DMPlexCreateCGNS(comm, cgid, interpolate, dm);CHKERRQ(ierr);
-    if (!rank) {ierr = cg_close(cgid);CHKERRQ(ierr);}
+      if (!rank) {
+        ierr = cg_open(filename, CG_MODE_READ, &cgid);CHKERRQ(ierr);
+        if (cgid <= 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_LIB, "cg_open(\"%s\",...) did not return a valid file ID", filename);
+      }
+      ierr = DMPlexCreateCGNS(comm, cgid, interpolate, dm);CHKERRQ(ierr);
+      if (!rank) {ierr = cg_close(cgid);CHKERRQ(ierr);}
 #else
-    SETERRQ(comm, PETSC_ERR_SUP, "Loading meshes requires CGNS support. Reconfigure using --with-cgns-dir");
+      SETERRQ(comm, PETSC_ERR_SUP, "Loading meshes requires CGNS support. Reconfigure using --with-cgns-dir");
 #endif
+    }
   } else if (cellSimplex) {
     ierr = DMPlexCreateBoxMesh(comm, dim, interpolate, dm);CHKERRQ(ierr);
   } else {

File src/dm/impls/plex/examples/tests/output/ex1_12.out

+Simplical Mesh in 3 dimensions:
+  0-cells: 5
+  1-cells: 9
+  2-cells: 7
+  3-cells: 2
+Labels:
+  depth: 4 strata of sizes (5, 9, 7, 2)

File src/dm/impls/plex/examples/tutorials/output/ex2_7.out

Empty file added.

File src/dm/impls/plex/makefile

 CPPFLAGS =
 CFLAGS   =
 FFLAGS   =
-SOURCEC  = plexcreate.c plex.c plexrefine.c plexinterpolate.c plexpreallocate.c plexreorder.c plexgeometry.c plexlabel.c plexsubmesh.c plexexodusii.c plexcgns.c plexvtk.c plexpoint.c plexvtu.c plexfem.c
+SOURCEC  = plexcreate.c plex.c plexrefine.c plexinterpolate.c plexpreallocate.c plexreorder.c plexgeometry.c plexlabel.c plexsubmesh.c plexexodusii.c plexgmsh.c plexcgns.c plexvtk.c plexpoint.c plexvtu.c plexfem.c
 SOURCEF  =
 SOURCEH  =
 DIRS     = examples

File src/dm/impls/plex/plexgmsh.c

+#define PETSCDM_DLL
+#include <petsc-private/dmpleximpl.h>    /*I   "petscdmplex.h"   I*/
+
+#undef __FUNCT__
+#define __FUNCT__ "DMPlexCreateGmsh"
+/*@
+  DMPlexCreateGmsh - Create a DMPlex mesh from a Gmsh file.
+
+  Collective on comm
+
+  Input Parameters:
++ comm  - The MPI communicator
+. viewer - The Viewer associated with a Gmsh file
+- interpolate - Create faces and edges in the mesh
+
+  Output Parameter:
+. dm  - The DM object representing the mesh
+
+  Note: http://www.geuz.org/gmsh/doc/texinfo/#MSH-ASCII-file-format
+
+  Level: beginner
+
+.keywords: mesh,Gmsh
+.seealso: DMPLEX, DMCreate()
+@*/
+PetscErrorCode DMPlexCreateGmsh(MPI_Comm comm, PetscViewer viewer, PetscBool interpolate, DM *dm)
+{
+  FILE          *fd;
+  PetscSection   coordSection;
+  Vec            coordinates;
+  PetscScalar   *coords, *coordsIn = NULL;
+  PetscInt       dim = 0, coordSize, c, v, d;
+  int            numVertices = 0, numCells = 0, snum;
+  long           fpos = 0;
+  PetscMPIInt    num_proc, rank;
+  char           line[PETSC_MAX_PATH_LEN];
+  PetscBool      match;
+  PetscErrorCode ierr;
+
+  PetscFunctionBegin;
+  ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr);
+  ierr = MPI_Comm_size(comm, &num_proc);CHKERRQ(ierr);
+  ierr = DMCreate(comm, dm);CHKERRQ(ierr);
+  ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr);
+  if (!rank) {
+    PetscBool match;
+    int       fileType, dataSize;
+
+    ierr = PetscViewerASCIIGetPointer(viewer, &fd);CHKERRQ(ierr);
+    /* Read header */
+    fgets(line, PETSC_MAX_PATH_LEN, fd);
+    ierr = PetscStrncmp(line, "$MeshFormat\n", PETSC_MAX_PATH_LEN, &match);CHKERRQ(ierr);
+    if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
+    snum = fscanf(fd, "2.2 %d %d\n", &fileType, &dataSize);CHKERRQ(snum != 2);
+    if (fileType) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File type %d is not a valid Gmsh ASCII file", fileType);
+    if (dataSize != sizeof(double)) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Data size %d is not valid for a Gmsh file", dataSize);
+    fgets(line, PETSC_MAX_PATH_LEN, fd);
+    ierr = PetscStrncmp(line, "$EndMeshFormat\n", PETSC_MAX_PATH_LEN, &match);CHKERRQ(ierr);
+    if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
+    /* Read vertices */
+    fgets(line, PETSC_MAX_PATH_LEN, fd);
+    ierr = PetscStrncmp(line, "$Nodes\n", PETSC_MAX_PATH_LEN, &match);CHKERRQ(ierr);
+    if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
+    snum = fscanf(fd, "%d\n", &numVertices);CHKERRQ(snum != 1);
+    ierr = PetscMalloc(numVertices*3 * sizeof(PetscScalar), &coordsIn);CHKERRQ(ierr);
+    for (v = 0; v < numVertices; ++v) {
+      double x, y, z;
+      int    i;
+
+      snum = fscanf(fd, "%d %lg %lg %lg\n", &i, &x, &y, &z);CHKERRQ(snum != 4);
+      coordsIn[v*3+0] = x; coordsIn[v*3+1] = y; coordsIn[v*3+2] = z;
+      if (i != v+1) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid node number %d should be %d", i, v+1);
+    }
+    fgets(line, PETSC_MAX_PATH_LEN, fd);
+    ierr = PetscStrncmp(line, "$EndNodes\n", PETSC_MAX_PATH_LEN, &match);CHKERRQ(ierr);
+    if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
+    /* Read cells */
+    fgets(line, PETSC_MAX_PATH_LEN, fd);
+    ierr = PetscStrncmp(line, "$Elements\n", PETSC_MAX_PATH_LEN, &match);CHKERRQ(ierr);
+    if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
+    snum = fscanf(fd, "%d\n", &numCells);CHKERRQ(snum != 1);
+  }
+  ierr = DMPlexSetChart(*dm, 0, numCells+numVertices);CHKERRQ(ierr);
+  if (!rank) {
+    fpos = ftell(fd);
+    for (c = 0; c < numCells; ++c) {
+      PetscInt numCorners, t;
+      int      cone[8], i, cellType, numTags, tag;
+
+      snum = fscanf(fd, "%d %d %d", &i, &cellType, &numTags);CHKERRQ(snum != 3);
+      if (numTags) for (t = 0; t < numTags; ++t) {snum = fscanf(fd, "%d", &tag);CHKERRQ(snum != 1);}
+      switch (cellType) {
+      case 1: /* 2-node line */
+        dim = 1;
+        numCorners = 2;
+        snum = fscanf(fd, "%d %d\n", &cone[0], &cone[1]);CHKERRQ(snum != numCorners);
+        break;
+      case 2: /* 3-node triangle */
+        dim = 2;
+        numCorners = 3;
+        snum = fscanf(fd, "%d %d %d\n", &cone[0], &cone[1], &cone[2]);CHKERRQ(snum != numCorners);
+        break;
+      case 3: /* 4-node quadrangle */
+        dim = 2;
+        numCorners = 4;
+        snum = fscanf(fd, "%d %d %d %d\n", &cone[0], &cone[1], &cone[2], &cone[3]);CHKERRQ(snum != numCorners);
+        break;
+      case 4: /* 4-node tetrahedron */
+        dim  = 3;
+        numCorners = 4;
+        snum = fscanf(fd, "%d %d %d %d\n", &cone[0], &cone[1], &cone[2], &cone[3]);CHKERRQ(snum != numCorners);
+        break;
+      case 5: /* 8-node hexahedron */
+        dim = 3;
+        numCorners = 8;
+        snum = fscanf(fd, "%d %d %d %d %d %d %d %d\n", &cone[0], &cone[1], &cone[2], &cone[3], &cone[4], &cone[5], &cone[6], &cone[7]);CHKERRQ(snum != numCorners);
+        break;
+      default:
+        SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unsupported Gmsh element type %d", cellType);
+      }
+      ierr = DMPlexSetConeSize(*dm, c, numCorners);CHKERRQ(ierr);
+      if (i != c+1) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid cell number %d should be %d", i, c+1);
+    }
+  }
+  ierr = DMSetUp(*dm);CHKERRQ(ierr);
+  if (!rank) {
+    ierr = fseek(fd, fpos, SEEK_SET);CHKERRQ(ierr);
+    for (c = 0; c < numCells; ++c) {
+      PetscInt pcone[8], numCorners, corner, t;
+      int      cone[8], i, cellType, numTags, tag;
+
+      snum = fscanf(fd, "%d %d %d", &i, &cellType, &numTags);CHKERRQ(snum != 3);
+      if (numTags) for (t = 0; t < numTags; ++t) {snum = fscanf(fd, "%d", &tag);CHKERRQ(snum != 1);}
+      switch (cellType) {
+      case 1: /* 2-node line */
+        dim = 1;
+        numCorners = 2;
+        snum = fscanf(fd, "%d %d\n", &cone[0], &cone[1]);CHKERRQ(snum != numCorners);
+        break;
+      case 2: /* 3-node triangle */
+        dim = 2;
+        numCorners = 3;
+        snum = fscanf(fd, "%d %d %d\n", &cone[0], &cone[1], &cone[2]);CHKERRQ(snum != numCorners);
+        break;
+      case 3: /* 4-node quadrangle */
+        dim = 2;
+        numCorners = 4;
+        snum = fscanf(fd, "%d %d %d %d\n", &cone[0], &cone[1], &cone[2], &cone[3]);CHKERRQ(snum != numCorners);
+        break;
+      case 4: /* 4-node tetrahedron */
+        dim  = 3;
+        numCorners = 4;
+        snum = fscanf(fd, "%d %d %d %d\n", &cone[0], &cone[1], &cone[2], &cone[3]);CHKERRQ(snum != numCorners);
+        ierr = DMPlexInvertCell(dim, numCorners, cone);CHKERRQ(ierr);
+        break;
+      case 5: /* 8-node hexahedron */
+        dim = 3;
+        numCorners = 8;
+        snum = fscanf(fd, "%d %d %d %d %d %d %d %d\n", &cone[0], &cone[1], &cone[2], &cone[3], &cone[4], &cone[5], &cone[6], &cone[7]);CHKERRQ(snum != numCorners);
+        ierr = DMPlexInvertCell(dim, numCorners, cone);CHKERRQ(ierr);
+        break;
+      default:
+        SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unsupported Gmsh element type %d", cellType);
+      }
+      for (corner = 0; corner < numCorners; ++corner) pcone[corner] = cone[corner] + numCells-1;
+      ierr = DMPlexSetCone(*dm, c, pcone);CHKERRQ(ierr);
+      if (i != c+1) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid cell number %d should be %d", i, c+1);
+    }
+    fgets(line, PETSC_MAX_PATH_LEN, fd);
+    ierr = PetscStrncmp(line, "$EndElements\n", PETSC_MAX_PATH_LEN, &match);CHKERRQ(ierr);
+    if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
+  }
+  ierr = MPI_Bcast(&dim, 1, MPIU_INT, 0, comm);CHKERRQ(ierr);
+  ierr = DMPlexSetDimension(*dm, dim);CHKERRQ(ierr);
+  ierr = DMPlexSymmetrize(*dm);CHKERRQ(ierr);
+  ierr = DMPlexStratify(*dm);CHKERRQ(ierr);
+  if (interpolate) {
+    DM idm;
+
+    ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr);
+    ierr = DMDestroy(dm);CHKERRQ(ierr);
+    *dm  = idm;
+  }
+  /* Read coordinates */
+  ierr = DMGetCoordinateSection(*dm, &coordSection);CHKERRQ(ierr);
+  ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr);
+  ierr = PetscSectionSetFieldComponents(coordSection, 0, dim);CHKERRQ(ierr);
+  ierr = PetscSectionSetChart(coordSection, numCells, numCells + numVertices);CHKERRQ(ierr);
+  for (v = numCells; v < numCells+numVertices; ++v) {
+    ierr = PetscSectionSetDof(coordSection, v, dim);CHKERRQ(ierr);
+    ierr = PetscSectionSetFieldDof(coordSection, v, 0, dim);CHKERRQ(ierr);
+  }
+  ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr);
+  ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr);
+  ierr = VecCreate(comm, &coordinates);CHKERRQ(ierr);
+  ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr);
+  ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr);
+  ierr = VecSetType(coordinates, VECSTANDARD);CHKERRQ(ierr);
+  ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr);
+  if (!rank) {
+    for (v = 0; v < numVertices; ++v) {
+      for (d = 0; d < dim; ++d) {
+        coords[v*dim+d] = coordsIn[v*3+d];
+      }
+    }
+  }
+  ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr);
+  ierr = PetscFree(coordsIn);CHKERRQ(ierr);
+  ierr = DMSetCoordinatesLocal(*dm, coordinates);CHKERRQ(ierr);
+  ierr = VecDestroy(&coordinates);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}