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