1 #define PETSCDM_DLL 2 #include <petsc-private/dmpleximpl.h> /*I "petscdmplex.h" I*/ 3 4 #undef __FUNCT__ 5 #define __FUNCT__ "DMPlexCreateGmshFromFile" 6 /*@C 7 DMPlexCreateGmshFromFile - Create a DMPlex mesh from a Gmsh file 8 9 + comm - The MPI communicator 10 . filename - Name of the Gmsh file 11 - interpolate - Create faces and edges in the mesh 12 13 Output Parameter: 14 . dm - The DM object representing the mesh 15 16 Level: beginner 17 18 .seealso: DMPlexCreateFromFile(), DMPlexCreateGmsh(), DMPlexCreate() 19 @*/ 20 PetscErrorCode DMPlexCreateGmshFromFile(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm) 21 { 22 PetscViewer viewer, vheader; 23 PetscMPIInt rank; 24 PetscViewerType vtype; 25 char line[PETSC_MAX_PATH_LEN]; 26 int snum; 27 PetscBool match; 28 int fileType; 29 float version; 30 PetscErrorCode ierr; 31 32 PetscFunctionBegin; 33 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 34 /* Determine Gmsh file type (ASCII or binary) from file header */ 35 ierr = PetscViewerCreate(comm, &vheader);CHKERRQ(ierr); 36 ierr = PetscViewerSetType(vheader, PETSCVIEWERASCII);CHKERRQ(ierr); 37 ierr = PetscViewerFileSetMode(vheader, FILE_MODE_READ);CHKERRQ(ierr); 38 ierr = PetscViewerFileSetName(vheader, filename);CHKERRQ(ierr); 39 if (!rank) { 40 /* Read only the first two lines of the Gmsh file */ 41 ierr = PetscViewerRead(vheader, line, 1, PETSC_STRING);CHKERRQ(ierr); 42 ierr = PetscStrncmp(line, "$MeshFormat", PETSC_MAX_PATH_LEN, &match);CHKERRQ(ierr); 43 if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file"); 44 ierr = PetscViewerRead(vheader, line, 2, PETSC_STRING);CHKERRQ(ierr); 45 snum = sscanf(line, "%f %d", &version, &fileType); 46 if (snum != 2) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse Gmsh file header: %s", line); 47 if (version < 2.0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Gmsh file must be at least version 2.0"); 48 } 49 ierr = MPI_Bcast(&fileType, 1, MPIU_INT, 0, comm);CHKERRQ(ierr); 50 /* Create appropriate viewer and build plex */ 51 if (fileType == 0) vtype = PETSCVIEWERASCII; 52 else vtype = PETSCVIEWERBINARY; 53 ierr = PetscViewerCreate(comm, &viewer);CHKERRQ(ierr); 54 ierr = PetscViewerSetType(viewer, vtype);CHKERRQ(ierr); 55 ierr = PetscViewerFileSetMode(viewer, FILE_MODE_READ);CHKERRQ(ierr); 56 ierr = PetscViewerFileSetName(viewer, filename);CHKERRQ(ierr); 57 ierr = DMPlexCreateGmsh(comm, viewer, interpolate, dm);CHKERRQ(ierr); 58 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 59 ierr = PetscViewerDestroy(&vheader);CHKERRQ(ierr); 60 PetscFunctionReturn(0); 61 } 62 63 64 #undef __FUNCT__ 65 #define __FUNCT__ "DMPlexCreateGmsh" 66 /*@ 67 DMPlexCreateGmsh - Create a DMPlex mesh from a Gmsh file viewer 68 69 Collective on comm 70 71 Input Parameters: 72 + comm - The MPI communicator 73 . viewer - The Viewer associated with a Gmsh file 74 - interpolate - Create faces and edges in the mesh 75 76 Output Parameter: 77 . dm - The DM object representing the mesh 78 79 Note: http://www.geuz.org/gmsh/doc/texinfo/#MSH-ASCII-file-format 80 81 Level: beginner 82 83 .keywords: mesh,Gmsh 84 .seealso: DMPLEX, DMCreate() 85 @*/ 86 PetscErrorCode DMPlexCreateGmsh(MPI_Comm comm, PetscViewer viewer, PetscBool interpolate, DM *dm) 87 { 88 PetscViewerType vtype; 89 GmshElement *gmsh_elem; 90 PetscSection coordSection; 91 Vec coordinates; 92 PetscScalar *coords, *coordsIn = NULL; 93 PetscInt dim = 0, coordSize, c, v, d, cell; 94 int numVertices = 0, numCells = 0, trueNumCells = 0, snum; 95 PetscMPIInt num_proc, rank; 96 char line[PETSC_MAX_PATH_LEN]; 97 PetscBool match, binary, bswap = PETSC_FALSE; 98 PetscErrorCode ierr; 99 100 PetscFunctionBegin; 101 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 102 ierr = MPI_Comm_size(comm, &num_proc);CHKERRQ(ierr); 103 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 104 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 105 ierr = PetscViewerGetType(viewer, &vtype);CHKERRQ(ierr); 106 ierr = PetscStrcmp(vtype, PETSCVIEWERBINARY, &binary);CHKERRQ(ierr); 107 if (!rank || binary) { 108 PetscBool match; 109 int fileType, dataSize; 110 float version; 111 112 /* Read header */ 113 ierr = PetscViewerRead(viewer, line, 1, PETSC_STRING);CHKERRQ(ierr); 114 ierr = PetscStrncmp(line, "$MeshFormat", PETSC_MAX_PATH_LEN, &match);CHKERRQ(ierr); 115 if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file"); 116 ierr = PetscViewerRead(viewer, line, 3, PETSC_STRING);CHKERRQ(ierr); 117 snum = sscanf(line, "%f %d %d", &version, &fileType, &dataSize); 118 if (snum != 3) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse Gmsh file header: %s", line); 119 if (version < 2.0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Gmsh file must be at least version 2.0"); 120 if (dataSize != sizeof(double)) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Data size %d is not valid for a Gmsh file", dataSize); 121 if (binary) { 122 int checkInt; 123 ierr = PetscViewerRead(viewer, &checkInt, 1, PETSC_ENUM);CHKERRQ(ierr); 124 if (checkInt != 1) { 125 ierr = PetscByteSwap(&checkInt, PETSC_ENUM, 1);CHKERRQ(ierr); 126 if (checkInt == 1) bswap = PETSC_TRUE; 127 else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File type %d is not a valid Gmsh binary file", fileType); 128 } 129 } else if (fileType) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File type %d is not a valid Gmsh ASCII file", fileType); 130 ierr = PetscViewerRead(viewer, line, 1, PETSC_STRING);CHKERRQ(ierr); 131 ierr = PetscStrncmp(line, "$EndMeshFormat", PETSC_MAX_PATH_LEN, &match);CHKERRQ(ierr); 132 if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file"); 133 /* Read vertices */ 134 ierr = PetscViewerRead(viewer, line, 1, PETSC_STRING);CHKERRQ(ierr); 135 ierr = PetscStrncmp(line, "$Nodes", PETSC_MAX_PATH_LEN, &match);CHKERRQ(ierr); 136 if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file"); 137 ierr = PetscViewerRead(viewer, line, 1, PETSC_STRING);CHKERRQ(ierr); 138 snum = sscanf(line, "%d", &numVertices); 139 if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file"); 140 ierr = PetscMalloc1(numVertices*3, &coordsIn);CHKERRQ(ierr); 141 for (v = 0; v < numVertices; ++v) { 142 int i; 143 ierr = PetscViewerRead(viewer, &i, 1, PETSC_ENUM);CHKERRQ(ierr); 144 if (bswap) ierr = PetscByteSwap(&i, PETSC_ENUM, 1);CHKERRQ(ierr); 145 ierr = PetscViewerRead(viewer, &(coordsIn[v*3]), 3, PETSC_DOUBLE);CHKERRQ(ierr); 146 if (bswap) ierr = PetscByteSwap(&(coordsIn[v*3]), PETSC_DOUBLE, 3);CHKERRQ(ierr); 147 if (i != (int)v+1) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid node number %d should be %d", i, v+1); 148 } 149 ierr = PetscViewerRead(viewer, line, 1, PETSC_STRING);CHKERRQ(ierr);; 150 ierr = PetscStrncmp(line, "$EndNodes", PETSC_MAX_PATH_LEN, &match);CHKERRQ(ierr); 151 if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file"); 152 /* Read cells */ 153 ierr = PetscViewerRead(viewer, line, 1, PETSC_STRING);CHKERRQ(ierr);; 154 ierr = PetscStrncmp(line, "$Elements", PETSC_MAX_PATH_LEN, &match);CHKERRQ(ierr); 155 if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file"); 156 ierr = PetscViewerRead(viewer, line, 1, PETSC_STRING);CHKERRQ(ierr);; 157 snum = sscanf(line, "%d", &numCells); 158 if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file"); 159 } 160 161 if (!rank || binary) { 162 /* Gmsh elements can be of any dimension/co-dimension, so we need to traverse the 163 file contents multiple times to figure out the true number of cells and facets 164 in the given mesh. To make this more efficient we read the file contents only 165 once and store them in memory, while determining the true number of cells. */ 166 ierr = PetscMalloc1(numCells, &gmsh_elem);CHKERRQ(ierr); 167 for (trueNumCells=0, c = 0; c < numCells; ++c) { 168 ierr = DMPlexCreateGmsh_ReadElement(viewer, binary, bswap, &gmsh_elem[c]);CHKERRQ(ierr); 169 if (gmsh_elem[c].dim > dim) {dim = gmsh_elem[c].dim; trueNumCells = 0;} 170 if (gmsh_elem[c].dim == dim) trueNumCells++; 171 } 172 ierr = PetscViewerRead(viewer, line, 1, PETSC_STRING);CHKERRQ(ierr);; 173 ierr = PetscStrncmp(line, "$EndElements", PETSC_MAX_PATH_LEN, &match);CHKERRQ(ierr); 174 if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file"); 175 } 176 /* For binary we read on all ranks, but only build the plex on rank 0 */ 177 if (binary && rank) {trueNumCells = 0; numVertices = 0;}; 178 /* Allocate the cell-vertex mesh */ 179 ierr = DMPlexSetChart(*dm, 0, trueNumCells+numVertices);CHKERRQ(ierr); 180 if (!rank) { 181 for (cell = 0, c = 0; c < numCells; ++c) { 182 if (gmsh_elem[c].dim == dim) { 183 ierr = DMPlexSetConeSize(*dm, cell, gmsh_elem[c].numNodes);CHKERRQ(ierr); 184 cell++; 185 } 186 } 187 } 188 ierr = DMSetUp(*dm);CHKERRQ(ierr); 189 /* Add cell-vertex connections */ 190 if (!rank) { 191 PetscInt pcone[8], corner; 192 for (cell = 0, c = 0; c < numCells; ++c) { 193 if (gmsh_elem[c].dim == dim) { 194 for (corner = 0; corner < gmsh_elem[c].numNodes; ++corner) { 195 pcone[corner] = gmsh_elem[c].nodes[corner] + trueNumCells-1; 196 } 197 ierr = DMPlexSetCone(*dm, cell, pcone);CHKERRQ(ierr); 198 cell++; 199 } 200 } 201 } 202 ierr = MPI_Bcast(&dim, 1, MPIU_INT, 0, comm);CHKERRQ(ierr); 203 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 204 ierr = DMPlexSymmetrize(*dm);CHKERRQ(ierr); 205 ierr = DMPlexStratify(*dm);CHKERRQ(ierr); 206 if (interpolate) { 207 DM idm = NULL; 208 209 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 210 ierr = DMDestroy(dm);CHKERRQ(ierr); 211 *dm = idm; 212 } 213 214 if (!rank) { 215 /* Apply boundary IDs by finding the relevant facets with vertex joins */ 216 PetscInt pcone[8], corner, vStart, vEnd; 217 218 ierr = DMPlexGetDepthStratum(*dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 219 for (c = 0; c < numCells; ++c) { 220 if (gmsh_elem[c].dim == dim-1) { 221 PetscInt joinSize; 222 const PetscInt *join; 223 for (corner = 0; corner < gmsh_elem[c].numNodes; ++corner) { 224 pcone[corner] = gmsh_elem[c].nodes[corner] + vStart - 1; 225 } 226 ierr = DMPlexGetFullJoin(*dm, gmsh_elem[c].numNodes, (const PetscInt *) pcone, &joinSize, &join);CHKERRQ(ierr); 227 if (joinSize != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not determine Plex facet for element %d", gmsh_elem[c].id); 228 ierr = DMPlexSetLabelValue(*dm, "Face Sets", join[0], gmsh_elem[c].tags[0]);CHKERRQ(ierr); 229 ierr = DMPlexRestoreJoin(*dm, gmsh_elem[c].numNodes, (const PetscInt *) pcone, &joinSize, &join);CHKERRQ(ierr); 230 } 231 } 232 } 233 234 /* Read coordinates */ 235 ierr = DMGetCoordinateSection(*dm, &coordSection);CHKERRQ(ierr); 236 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 237 ierr = PetscSectionSetFieldComponents(coordSection, 0, dim);CHKERRQ(ierr); 238 ierr = PetscSectionSetChart(coordSection, trueNumCells, trueNumCells + numVertices);CHKERRQ(ierr); 239 for (v = trueNumCells; v < trueNumCells+numVertices; ++v) { 240 ierr = PetscSectionSetDof(coordSection, v, dim);CHKERRQ(ierr); 241 ierr = PetscSectionSetFieldDof(coordSection, v, 0, dim);CHKERRQ(ierr); 242 } 243 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 244 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 245 ierr = VecCreate(comm, &coordinates);CHKERRQ(ierr); 246 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 247 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 248 ierr = VecSetType(coordinates, VECSTANDARD);CHKERRQ(ierr); 249 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 250 if (!rank) { 251 for (v = 0; v < numVertices; ++v) { 252 for (d = 0; d < dim; ++d) { 253 coords[v*dim+d] = coordsIn[v*3+d]; 254 } 255 } 256 } 257 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 258 ierr = PetscFree(coordsIn);CHKERRQ(ierr); 259 ierr = DMSetCoordinatesLocal(*dm, coordinates);CHKERRQ(ierr); 260 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 261 /* Clean up intermediate storage */ 262 if (!rank || binary) { 263 for (c = 0; c < numCells; ++c) { 264 ierr = PetscFree(gmsh_elem[c].nodes); 265 ierr = PetscFree(gmsh_elem[c].tags); 266 } 267 ierr = PetscFree(gmsh_elem); 268 } 269 PetscFunctionReturn(0); 270 } 271 272 #undef __FUNCT__ 273 #define __FUNCT__ "DMPlexCreateGmsh_ReadElement" 274 PetscErrorCode DMPlexCreateGmsh_ReadElement(PetscViewer viewer, PetscBool binary, PetscBool byteSwap, GmshElement *ele) 275 { 276 int cellType, numElem; 277 PetscErrorCode ierr; 278 279 PetscFunctionBegin; 280 if (binary) { 281 ierr = PetscViewerRead(viewer, &cellType, 1, PETSC_ENUM);CHKERRQ(ierr); 282 if (byteSwap) ierr = PetscByteSwap(&cellType, PETSC_ENUM, 1);CHKERRQ(ierr); 283 ierr = PetscViewerRead(viewer, &numElem, 1, PETSC_ENUM);CHKERRQ(ierr); 284 if (byteSwap) ierr = PetscByteSwap(&numElem, PETSC_ENUM, 1);CHKERRQ(ierr); 285 ierr = PetscViewerRead(viewer, &(ele->numTags), 1, PETSC_ENUM);CHKERRQ(ierr); 286 if (byteSwap) ierr = PetscByteSwap(&(ele->numTags), PETSC_ENUM, 1);CHKERRQ(ierr); 287 ierr = PetscViewerRead(viewer, &(ele->id), 1, PETSC_ENUM);CHKERRQ(ierr); 288 if (byteSwap) ierr = PetscByteSwap( &(ele->id), PETSC_ENUM, 1);CHKERRQ(ierr); 289 } else { 290 ierr = PetscViewerRead(viewer, &(ele->id), 1, PETSC_ENUM);CHKERRQ(ierr); 291 ierr = PetscViewerRead(viewer, &cellType, 1, PETSC_ENUM);CHKERRQ(ierr); 292 ierr = PetscViewerRead(viewer, &(ele->numTags), 1, PETSC_ENUM);CHKERRQ(ierr); 293 } 294 ierr = PetscMalloc1(ele->numTags, &(ele->tags));CHKERRQ(ierr); 295 ierr = PetscViewerRead(viewer, ele->tags, ele->numTags, PETSC_ENUM);CHKERRQ(ierr); 296 if (byteSwap) ierr = PetscByteSwap(ele->tags, PETSC_ENUM, ele->numTags);CHKERRQ(ierr); 297 switch (cellType) { 298 case 1: /* 2-node line */ 299 ele->dim = 1; 300 ele->numNodes = 2; 301 break; 302 case 2: /* 3-node triangle */ 303 ele->dim = 2; 304 ele->numNodes = 3; 305 break; 306 case 3: /* 4-node quadrangle */ 307 ele->dim = 2; 308 ele->numNodes = 4; 309 break; 310 case 4: /* 4-node tetrahedron */ 311 ele->dim = 3; 312 ele->numNodes = 4; 313 break; 314 case 5: /* 8-node hexahedron */ 315 ele->dim = 3; 316 ele->numNodes = 8; 317 break; 318 case 15: /* 1-node vertex */ 319 ele->dim = 0; 320 ele->numNodes = 1; 321 break; 322 default: 323 SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unsupported Gmsh element type %d", cellType); 324 } 325 ierr = PetscMalloc1(ele->numNodes, &(ele->nodes));CHKERRQ(ierr); 326 ierr = PetscViewerRead(viewer, ele->nodes, ele->numNodes, PETSC_ENUM);CHKERRQ(ierr); 327 if (byteSwap) {ierr = PetscByteSwap(ele->nodes, PETSC_ENUM, ele->numNodes);CHKERRQ(ierr);} 328 PetscFunctionReturn(0); 329 } 330