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 273 /* Read coordinates */ 274 ierr = DMGetCoordinateSection(*dm, &coordSection);CHKERRQ(ierr); 275 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 276 ierr = PetscSectionSetFieldComponents(coordSection, 0, dim);CHKERRQ(ierr); 277 ierr = PetscSectionSetChart(coordSection, trueNumCells, trueNumCells + numVertices);CHKERRQ(ierr); 278 for (v = trueNumCells; v < trueNumCells+numVertices; ++v) { 279 ierr = PetscSectionSetDof(coordSection, v, dim);CHKERRQ(ierr); 280 ierr = PetscSectionSetFieldDof(coordSection, v, 0, dim);CHKERRQ(ierr); 281 } 282 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 283 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 284 ierr = VecCreate(comm, &coordinates);CHKERRQ(ierr); 285 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 286 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 287 ierr = VecSetType(coordinates, VECSTANDARD);CHKERRQ(ierr); 288 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 289 if (!rank) { 290 for (v = 0; v < numVertices; ++v) { 291 for (d = 0; d < dim; ++d) { 292 coords[v*dim+d] = coordsIn[v*3+d]; 293 } 294 } 295 } 296 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 297 ierr = PetscFree(coordsIn);CHKERRQ(ierr); 298 ierr = DMSetCoordinatesLocal(*dm, coordinates);CHKERRQ(ierr); 299 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 300 /* Clean up intermediate storage */ 301 if (!rank || binary) ierr = PetscFree(gmsh_elem);CHKERRQ(ierr); 302 ierr = PetscLogEventEnd(DMPLEX_CreateGmsh,*dm,0,0,0);CHKERRQ(ierr); 303 PetscFunctionReturn(0); 304 } 305 306 #undef __FUNCT__ 307 #define __FUNCT__ "DMPlexCreateGmsh_ReadElement" 308 PetscErrorCode DMPlexCreateGmsh_ReadElement(PetscViewer viewer, PetscInt numCells, PetscBool binary, PetscBool byteSwap, GmshElement **gmsh_elems) 309 { 310 PetscInt c, p; 311 GmshElement *elements; 312 int i, cellType, dim, numNodes, numElem, numTags; 313 int ibuf[16]; 314 PetscErrorCode ierr; 315 316 PetscFunctionBegin; 317 ierr = PetscMalloc1(numCells, &elements);CHKERRQ(ierr); 318 for (c = 0; c < numCells;) { 319 ierr = PetscViewerRead(viewer, &ibuf, 3, NULL, PETSC_ENUM);CHKERRQ(ierr); 320 if (byteSwap) ierr = PetscByteSwap(&ibuf, PETSC_ENUM, 3);CHKERRQ(ierr); 321 if (binary) { 322 cellType = ibuf[0]; 323 numElem = ibuf[1]; 324 numTags = ibuf[2]; 325 } else { 326 elements[c].id = ibuf[0]; 327 cellType = ibuf[1]; 328 numTags = ibuf[2]; 329 numElem = 1; 330 } 331 switch (cellType) { 332 case 1: /* 2-node line */ 333 dim = 1; 334 numNodes = 2; 335 break; 336 case 2: /* 3-node triangle */ 337 dim = 2; 338 numNodes = 3; 339 break; 340 case 3: /* 4-node quadrangle */ 341 dim = 2; 342 numNodes = 4; 343 break; 344 case 4: /* 4-node tetrahedron */ 345 dim = 3; 346 numNodes = 4; 347 break; 348 case 5: /* 8-node hexahedron */ 349 dim = 3; 350 numNodes = 8; 351 break; 352 case 15: /* 1-node vertex */ 353 dim = 0; 354 numNodes = 1; 355 break; 356 default: 357 SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unsupported Gmsh element type %d", cellType); 358 } 359 if (binary) { 360 const PetscInt nint = numNodes + numTags + 1; 361 for (i = 0; i < numElem; ++i, ++c) { 362 /* Loop over inner binary element block */ 363 elements[c].dim = dim; 364 elements[c].numNodes = numNodes; 365 elements[c].numTags = numTags; 366 367 ierr = PetscViewerRead(viewer, &ibuf, nint, NULL, PETSC_ENUM);CHKERRQ(ierr); 368 if (byteSwap) ierr = PetscByteSwap( &ibuf, PETSC_ENUM, nint);CHKERRQ(ierr); 369 elements[c].id = ibuf[0]; 370 for (p = 0; p < numTags; p++) elements[c].tags[p] = ibuf[1 + p]; 371 for (p = 0; p < numNodes; p++) elements[c].nodes[p] = ibuf[1 + numTags + p]; 372 } 373 } else { 374 elements[c].dim = dim; 375 elements[c].numNodes = numNodes; 376 elements[c].numTags = numTags; 377 ierr = PetscViewerRead(viewer, elements[c].tags, elements[c].numTags, NULL, PETSC_ENUM);CHKERRQ(ierr); 378 ierr = PetscViewerRead(viewer, elements[c].nodes, elements[c].numNodes, NULL, PETSC_ENUM);CHKERRQ(ierr); 379 c++; 380 } 381 } 382 *gmsh_elems = elements; 383 PetscFunctionReturn(0); 384 } 385