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