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, cell; 97 int i, numVertices = 0, numCells = 0, trueNumCells = 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 /* Read vertices */ 138 ierr = PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);CHKERRQ(ierr); 139 ierr = PetscStrncmp(line, "$Nodes", PETSC_MAX_PATH_LEN, &match);CHKERRQ(ierr); 140 if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file"); 141 ierr = PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);CHKERRQ(ierr); 142 snum = sscanf(line, "%d", &numVertices); 143 if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file"); 144 ierr = PetscMalloc1(numVertices*3, &coordsIn);CHKERRQ(ierr); 145 if (binary) { 146 size_t doubleSize, intSize; 147 PetscInt elementSize; 148 char *buffer; 149 PetscScalar *baseptr; 150 ierr = PetscDataTypeGetSize(PETSC_ENUM, &intSize);CHKERRQ(ierr); 151 ierr = PetscDataTypeGetSize(PETSC_DOUBLE, &doubleSize);CHKERRQ(ierr); 152 elementSize = (intSize + 3*doubleSize); 153 ierr = PetscMalloc1(elementSize*numVertices, &buffer);CHKERRQ(ierr); 154 ierr = PetscViewerRead(viewer, buffer, elementSize*numVertices, NULL, PETSC_CHAR);CHKERRQ(ierr); 155 if (bswap) ierr = PetscByteSwap(buffer, PETSC_CHAR, elementSize*numVertices);CHKERRQ(ierr); 156 for (v = 0; v < numVertices; ++v) { 157 baseptr = ((PetscScalar*)(buffer+v*elementSize+intSize)); 158 coordsIn[v*3+0] = baseptr[0]; 159 coordsIn[v*3+1] = baseptr[1]; 160 coordsIn[v*3+2] = baseptr[2]; 161 } 162 ierr = PetscFree(buffer);CHKERRQ(ierr); 163 } else { 164 for (v = 0; v < numVertices; ++v) { 165 ierr = PetscViewerRead(viewer, &i, 1, NULL, PETSC_ENUM);CHKERRQ(ierr); 166 ierr = PetscViewerRead(viewer, &(coordsIn[v*3]), 3, NULL, PETSC_DOUBLE);CHKERRQ(ierr); 167 if (i != (int)v+1) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid node number %d should be %d", i, v+1); 168 } 169 } 170 ierr = PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);CHKERRQ(ierr);; 171 ierr = PetscStrncmp(line, "$EndNodes", PETSC_MAX_PATH_LEN, &match);CHKERRQ(ierr); 172 if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file"); 173 /* Read cells */ 174 ierr = PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);CHKERRQ(ierr);; 175 ierr = PetscStrncmp(line, "$Elements", PETSC_MAX_PATH_LEN, &match);CHKERRQ(ierr); 176 if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file"); 177 ierr = PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);CHKERRQ(ierr);; 178 snum = sscanf(line, "%d", &numCells); 179 if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file"); 180 } 181 182 if (!rank || binary) { 183 /* Gmsh elements can be of any dimension/co-dimension, so we need to traverse the 184 file contents multiple times to figure out the true number of cells and facets 185 in the given mesh. To make this more efficient we read the file contents only 186 once and store them in memory, while determining the true number of cells. */ 187 ierr = DMPlexCreateGmsh_ReadElement(viewer, numCells, binary, bswap, &gmsh_elem);CHKERRQ(ierr); 188 for (trueNumCells=0, c = 0; c < numCells; ++c) { 189 if (gmsh_elem[c].dim > dim) {dim = gmsh_elem[c].dim; trueNumCells = 0;} 190 if (gmsh_elem[c].dim == dim) trueNumCells++; 191 } 192 ierr = PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);CHKERRQ(ierr);; 193 ierr = PetscStrncmp(line, "$EndElements", PETSC_MAX_PATH_LEN, &match);CHKERRQ(ierr); 194 if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file"); 195 } 196 /* For binary we read on all ranks, but only build the plex on rank 0 */ 197 if (binary && rank) {trueNumCells = 0; numVertices = 0;}; 198 /* Allocate the cell-vertex mesh */ 199 ierr = DMPlexSetChart(*dm, 0, trueNumCells+numVertices);CHKERRQ(ierr); 200 if (!rank) { 201 for (cell = 0, c = 0; c < numCells; ++c) { 202 if (gmsh_elem[c].dim == dim) { 203 ierr = DMPlexSetConeSize(*dm, cell, gmsh_elem[c].numNodes);CHKERRQ(ierr); 204 cell++; 205 } 206 } 207 } 208 ierr = DMSetUp(*dm);CHKERRQ(ierr); 209 /* Add cell-vertex connections */ 210 if (!rank) { 211 PetscInt pcone[8], corner; 212 for (cell = 0, c = 0; c < numCells; ++c) { 213 if (gmsh_elem[c].dim == dim) { 214 for (corner = 0; corner < gmsh_elem[c].numNodes; ++corner) { 215 pcone[corner] = gmsh_elem[c].nodes[corner] + trueNumCells-1; 216 } 217 ierr = DMPlexSetCone(*dm, cell, pcone);CHKERRQ(ierr); 218 cell++; 219 } 220 } 221 } 222 ierr = MPI_Bcast(&dim, 1, MPIU_INT, 0, comm);CHKERRQ(ierr); 223 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 224 ierr = DMPlexSymmetrize(*dm);CHKERRQ(ierr); 225 ierr = DMPlexStratify(*dm);CHKERRQ(ierr); 226 if (interpolate) { 227 DM idm = NULL; 228 229 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 230 ierr = DMDestroy(dm);CHKERRQ(ierr); 231 *dm = idm; 232 } 233 234 if (!rank) { 235 /* Apply boundary IDs by finding the relevant facets with vertex joins */ 236 PetscInt pcone[8], corner, vStart, vEnd; 237 238 ierr = DMPlexGetDepthStratum(*dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 239 for (c = 0; c < numCells; ++c) { 240 if (gmsh_elem[c].dim == dim-1) { 241 PetscInt joinSize; 242 const PetscInt *join; 243 for (corner = 0; corner < gmsh_elem[c].numNodes; ++corner) { 244 pcone[corner] = gmsh_elem[c].nodes[corner] + vStart - 1; 245 } 246 ierr = DMPlexGetFullJoin(*dm, gmsh_elem[c].numNodes, (const PetscInt *) pcone, &joinSize, &join);CHKERRQ(ierr); 247 if (joinSize != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not determine Plex facet for element %d", gmsh_elem[c].id); 248 ierr = DMSetLabelValue(*dm, "Face Sets", join[0], gmsh_elem[c].tags[0]);CHKERRQ(ierr); 249 ierr = DMPlexRestoreJoin(*dm, gmsh_elem[c].numNodes, (const PetscInt *) pcone, &joinSize, &join);CHKERRQ(ierr); 250 } 251 } 252 253 /* Create cell sets */ 254 for (cell = 0, c = 0; c < numCells; ++c) { 255 if (gmsh_elem[c].dim == dim) { 256 if (gmsh_elem[c].numTags > 0) { 257 ierr = DMSetLabelValue(*dm, "Cell Sets", cell, gmsh_elem[c].tags[0]);CHKERRQ(ierr); 258 cell++; 259 } 260 } 261 } 262 } 263 264 /* Read coordinates */ 265 ierr = DMGetCoordinateSection(*dm, &coordSection);CHKERRQ(ierr); 266 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 267 ierr = PetscSectionSetFieldComponents(coordSection, 0, dim);CHKERRQ(ierr); 268 ierr = PetscSectionSetChart(coordSection, trueNumCells, trueNumCells + numVertices);CHKERRQ(ierr); 269 for (v = trueNumCells; v < trueNumCells+numVertices; ++v) { 270 ierr = PetscSectionSetDof(coordSection, v, dim);CHKERRQ(ierr); 271 ierr = PetscSectionSetFieldDof(coordSection, v, 0, dim);CHKERRQ(ierr); 272 } 273 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 274 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 275 ierr = VecCreate(comm, &coordinates);CHKERRQ(ierr); 276 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 277 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 278 ierr = VecSetType(coordinates, VECSTANDARD);CHKERRQ(ierr); 279 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 280 if (!rank) { 281 for (v = 0; v < numVertices; ++v) { 282 for (d = 0; d < dim; ++d) { 283 coords[v*dim+d] = coordsIn[v*3+d]; 284 } 285 } 286 } 287 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 288 ierr = PetscFree(coordsIn);CHKERRQ(ierr); 289 ierr = DMSetCoordinatesLocal(*dm, coordinates);CHKERRQ(ierr); 290 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 291 /* Clean up intermediate storage */ 292 if (!rank || binary) ierr = PetscFree(gmsh_elem);CHKERRQ(ierr); 293 ierr = PetscLogEventEnd(DMPLEX_CreateGmsh,*dm,0,0,0);CHKERRQ(ierr); 294 PetscFunctionReturn(0); 295 } 296 297 #undef __FUNCT__ 298 #define __FUNCT__ "DMPlexCreateGmsh_ReadElement" 299 PetscErrorCode DMPlexCreateGmsh_ReadElement(PetscViewer viewer, PetscInt numCells, PetscBool binary, PetscBool byteSwap, GmshElement **gmsh_elems) 300 { 301 PetscInt c, p; 302 GmshElement *elements; 303 int i, cellType, dim, numNodes, numElem, numTags; 304 int ibuf[16]; 305 PetscErrorCode ierr; 306 307 PetscFunctionBegin; 308 ierr = PetscMalloc1(numCells, &elements);CHKERRQ(ierr); 309 for (c = 0; c < numCells;) { 310 ierr = PetscViewerRead(viewer, &ibuf, 3, NULL, PETSC_ENUM);CHKERRQ(ierr); 311 if (byteSwap) ierr = PetscByteSwap(&ibuf, PETSC_ENUM, 3);CHKERRQ(ierr); 312 if (binary) { 313 cellType = ibuf[0]; 314 numElem = ibuf[1]; 315 numTags = ibuf[2]; 316 } else { 317 elements[c].id = ibuf[0]; 318 cellType = ibuf[1]; 319 numTags = ibuf[2]; 320 numElem = 1; 321 } 322 switch (cellType) { 323 case 1: /* 2-node line */ 324 dim = 1; 325 numNodes = 2; 326 break; 327 case 2: /* 3-node triangle */ 328 dim = 2; 329 numNodes = 3; 330 break; 331 case 3: /* 4-node quadrangle */ 332 dim = 2; 333 numNodes = 4; 334 break; 335 case 4: /* 4-node tetrahedron */ 336 dim = 3; 337 numNodes = 4; 338 break; 339 case 5: /* 8-node hexahedron */ 340 dim = 3; 341 numNodes = 8; 342 break; 343 case 15: /* 1-node vertex */ 344 dim = 0; 345 numNodes = 1; 346 break; 347 default: 348 SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unsupported Gmsh element type %d", cellType); 349 } 350 if (binary) { 351 const PetscInt nint = numNodes + numTags + 1; 352 for (i = 0; i < numElem; ++i, ++c) { 353 /* Loop over inner binary element block */ 354 elements[c].dim = dim; 355 elements[c].numNodes = numNodes; 356 elements[c].numTags = numTags; 357 358 ierr = PetscViewerRead(viewer, &ibuf, nint, NULL, PETSC_ENUM);CHKERRQ(ierr); 359 if (byteSwap) ierr = PetscByteSwap( &ibuf, PETSC_ENUM, nint);CHKERRQ(ierr); 360 elements[c].id = ibuf[0]; 361 for (p = 0; p < numTags; p++) elements[c].tags[p] = ibuf[1 + p]; 362 for (p = 0; p < numNodes; p++) elements[c].nodes[p] = ibuf[1 + numTags + p]; 363 } 364 } else { 365 elements[c].dim = dim; 366 elements[c].numNodes = numNodes; 367 elements[c].numTags = numTags; 368 ierr = PetscViewerRead(viewer, elements[c].tags, elements[c].numTags, NULL, PETSC_ENUM);CHKERRQ(ierr); 369 ierr = PetscViewerRead(viewer, elements[c].nodes, elements[c].numNodes, NULL, PETSC_ENUM);CHKERRQ(ierr); 370 c++; 371 } 372 } 373 *gmsh_elems = elements; 374 PetscFunctionReturn(0); 375 } 376