1 #define PETSCDM_DLL 2 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/ 3 4 /*@C 5 DMPlexCreateGmshFromFile - Create a DMPlex mesh from a Gmsh file 6 7 + comm - The MPI communicator 8 . filename - Name of the Gmsh file 9 - interpolate - Create faces and edges in the mesh 10 11 Output Parameter: 12 . dm - The DM object representing the mesh 13 14 Level: beginner 15 16 .seealso: DMPlexCreateFromFile(), DMPlexCreateGmsh(), DMPlexCreate() 17 @*/ 18 PetscErrorCode DMPlexCreateGmshFromFile(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm) 19 { 20 PetscViewer viewer, vheader; 21 PetscMPIInt rank; 22 PetscViewerType vtype; 23 char line[PETSC_MAX_PATH_LEN]; 24 int snum; 25 PetscBool match; 26 int fT; 27 PetscInt fileType; 28 float version; 29 PetscErrorCode ierr; 30 31 PetscFunctionBegin; 32 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 33 /* Determine Gmsh file type (ASCII or binary) from file header */ 34 ierr = PetscViewerCreate(comm, &vheader);CHKERRQ(ierr); 35 ierr = PetscViewerSetType(vheader, PETSCVIEWERASCII);CHKERRQ(ierr); 36 ierr = PetscViewerFileSetMode(vheader, FILE_MODE_READ);CHKERRQ(ierr); 37 ierr = PetscViewerFileSetName(vheader, filename);CHKERRQ(ierr); 38 if (!rank) { 39 /* Read only the first two lines of the Gmsh file */ 40 ierr = PetscViewerRead(vheader, line, 1, NULL, PETSC_STRING);CHKERRQ(ierr); 41 ierr = PetscStrncmp(line, "$MeshFormat", PETSC_MAX_PATH_LEN, &match);CHKERRQ(ierr); 42 if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file"); 43 ierr = PetscViewerRead(vheader, line, 2, NULL, PETSC_STRING);CHKERRQ(ierr); 44 snum = sscanf(line, "%f %d", &version, &fT); 45 fileType = (PetscInt) fT; 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 static PetscErrorCode DMPlexCreateGmsh_ReadElement(PetscViewer viewer, PetscInt numCells, PetscBool binary, PetscBool byteSwap, GmshElement **gmsh_elems) 64 { 65 PetscInt c, p; 66 GmshElement *elements; 67 int i, cellType, dim, numNodes, numNodesIgnore, numElem, numTags; 68 PetscInt pibuf[64]; 69 int ibuf[16]; 70 PetscErrorCode ierr; 71 72 PetscFunctionBegin; 73 ierr = PetscMalloc1(numCells, &elements);CHKERRQ(ierr); 74 for (c = 0; c < numCells;) { 75 ierr = PetscViewerRead(viewer, &ibuf, 3, NULL, PETSC_ENUM);CHKERRQ(ierr); 76 if (byteSwap) ierr = PetscByteSwap(&ibuf, PETSC_ENUM, 3);CHKERRQ(ierr); 77 if (binary) { 78 cellType = ibuf[0]; 79 numElem = ibuf[1]; 80 numTags = ibuf[2]; 81 } else { 82 elements[c].id = ibuf[0]; 83 cellType = ibuf[1]; 84 numTags = ibuf[2]; 85 numElem = 1; 86 } 87 /* http://gmsh.info/doc/texinfo/gmsh.html#MSH-ASCII-file-format */ 88 numNodesIgnore = 0; 89 switch (cellType) { 90 case 1: /* 2-node line */ 91 dim = 1; 92 numNodes = 2; 93 break; 94 case 2: /* 3-node triangle */ 95 dim = 2; 96 numNodes = 3; 97 break; 98 case 3: /* 4-node quadrangle */ 99 dim = 2; 100 numNodes = 4; 101 break; 102 case 4: /* 4-node tetrahedron */ 103 dim = 3; 104 numNodes = 4; 105 break; 106 case 5: /* 8-node hexahedron */ 107 dim = 3; 108 numNodes = 8; 109 break; 110 case 8: /* 3-node 2nd order line */ 111 dim = 1; 112 numNodes = 2; 113 numNodesIgnore = 1; 114 break; 115 case 9: /* 6-node 2nd order triangle */ 116 dim = 2; 117 numNodes = 3; 118 numNodesIgnore = 3; 119 break; 120 case 15: /* 1-node vertex */ 121 dim = 0; 122 numNodes = 1; 123 break; 124 case 6: /* 6-node prism */ 125 case 7: /* 5-node pyramid */ 126 case 10: /* 9-node 2nd order quadrangle */ 127 case 11: /* 10-node 2nd order tetrahedron */ 128 case 12: /* 27-node 2nd order hexhedron */ 129 case 13: /* 19-node 2nd order prism */ 130 case 14: /* 14-node 2nd order pyramid */ 131 default: 132 SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unsupported Gmsh element type %d", cellType); 133 } 134 if (binary) { 135 const PetscInt nint = numNodes + numTags + 1 + numNodesIgnore; 136 for (i = 0; i < numElem; ++i, ++c) { 137 /* Loop over inner binary element block */ 138 elements[c].dim = dim; 139 elements[c].numNodes = numNodes; 140 elements[c].numTags = numTags; 141 142 ierr = PetscViewerRead(viewer, &ibuf, nint, NULL, PETSC_ENUM);CHKERRQ(ierr); 143 if (byteSwap) ierr = PetscByteSwap( &ibuf, PETSC_ENUM, nint);CHKERRQ(ierr); 144 elements[c].id = ibuf[0]; 145 for (p = 0; p < numTags; p++) elements[c].tags[p] = ibuf[1 + p]; 146 for (p = 0; p < numNodes; p++) elements[c].nodes[p] = ibuf[1 + numTags + p]; 147 } 148 } else { 149 elements[c].dim = dim; 150 elements[c].numNodes = numNodes; 151 elements[c].numTags = numTags; 152 ierr = PetscViewerRead(viewer, elements[c].tags, elements[c].numTags, NULL, PETSC_ENUM);CHKERRQ(ierr); 153 ierr = PetscViewerRead(viewer, elements[c].nodes, elements[c].numNodes, NULL, PETSC_ENUM);CHKERRQ(ierr); 154 ierr = PetscViewerRead(viewer, pibuf, numNodesIgnore, NULL, PETSC_ENUM);CHKERRQ(ierr); 155 c++; 156 } 157 } 158 *gmsh_elems = elements; 159 PetscFunctionReturn(0); 160 } 161 162 /*@ 163 DMPlexCreateGmsh - Create a DMPlex mesh from a Gmsh file viewer 164 165 Collective on comm 166 167 Input Parameters: 168 + comm - The MPI communicator 169 . viewer - The Viewer associated with a Gmsh file 170 - interpolate - Create faces and edges in the mesh 171 172 Output Parameter: 173 . dm - The DM object representing the mesh 174 175 Note: http://www.geuz.org/gmsh/doc/texinfo/#MSH-ASCII-file-format 176 and http://www.geuz.org/gmsh/doc/texinfo/#MSH-binary-file-format 177 178 Level: beginner 179 180 .keywords: mesh,Gmsh 181 .seealso: DMPLEX, DMCreate() 182 @*/ 183 PetscErrorCode DMPlexCreateGmsh(MPI_Comm comm, PetscViewer viewer, PetscBool interpolate, DM *dm) 184 { 185 PetscViewerType vtype; 186 GmshElement *gmsh_elem = NULL; 187 PetscSection coordSection; 188 Vec coordinates; 189 PetscScalar *coords, *coordsIn = NULL; 190 PetscInt dim = 0, coordSize, c, v, d, r, cell, *periodicMap = NULL, *periodicMapI = NULL; 191 int i, numVertices = 0, numCells = 0, trueNumCells = 0, numRegions = 0, snum; 192 PetscMPIInt num_proc, rank; 193 char line[PETSC_MAX_PATH_LEN]; 194 PetscBool match, binary, bswap = PETSC_FALSE, periodic = PETSC_FALSE; 195 PetscErrorCode ierr; 196 197 PetscFunctionBegin; 198 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 199 ierr = MPI_Comm_size(comm, &num_proc);CHKERRQ(ierr); 200 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 201 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 202 ierr = PetscLogEventBegin(DMPLEX_CreateGmsh,*dm,0,0,0);CHKERRQ(ierr); 203 ierr = PetscViewerGetType(viewer, &vtype);CHKERRQ(ierr); 204 ierr = PetscStrcmp(vtype, PETSCVIEWERBINARY, &binary);CHKERRQ(ierr); 205 ierr = PetscOptionsGetBool(NULL,NULL,"-dm_plex_periodic_gmsh",&periodic,NULL);CHKERRQ(ierr); 206 if (!rank || binary) { 207 PetscBool match; 208 int fileType, dataSize; 209 float version; 210 211 /* Read header */ 212 ierr = PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);CHKERRQ(ierr); 213 ierr = PetscStrncmp(line, "$MeshFormat", PETSC_MAX_PATH_LEN, &match);CHKERRQ(ierr); 214 if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file"); 215 ierr = PetscViewerRead(viewer, line, 3, NULL, PETSC_STRING);CHKERRQ(ierr); 216 snum = sscanf(line, "%f %d %d", &version, &fileType, &dataSize); 217 if (snum != 3) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse Gmsh file header: %s", line); 218 if (version < 2.0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Gmsh file must be at least version 2.0"); 219 if (dataSize != sizeof(double)) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Data size %d is not valid for a Gmsh file", dataSize); 220 if (binary) { 221 int checkInt; 222 ierr = PetscViewerRead(viewer, &checkInt, 1, NULL, PETSC_ENUM);CHKERRQ(ierr); 223 if (checkInt != 1) { 224 ierr = PetscByteSwap(&checkInt, PETSC_ENUM, 1);CHKERRQ(ierr); 225 if (checkInt == 1) bswap = PETSC_TRUE; 226 else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File type %d is not a valid Gmsh binary file", fileType); 227 } 228 } else if (fileType) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File type %d is not a valid Gmsh ASCII file", fileType); 229 ierr = PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);CHKERRQ(ierr); 230 ierr = PetscStrncmp(line, "$EndMeshFormat", PETSC_MAX_PATH_LEN, &match);CHKERRQ(ierr); 231 if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file"); 232 /* OPTIONAL Read physical names */ 233 ierr = PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);CHKERRQ(ierr); 234 ierr = PetscStrncmp(line, "$PhysicalNames", PETSC_MAX_PATH_LEN, &match);CHKERRQ(ierr); 235 if (match) { 236 ierr = PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);CHKERRQ(ierr); 237 snum = sscanf(line, "%d", &numRegions); 238 if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file"); 239 for (r = 0; r < numRegions; ++r) {ierr = PetscViewerRead(viewer, line, 3, NULL, PETSC_STRING);CHKERRQ(ierr);} 240 ierr = PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);CHKERRQ(ierr); 241 ierr = PetscStrncmp(line, "$EndPhysicalNames", PETSC_MAX_PATH_LEN, &match);CHKERRQ(ierr); 242 if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file"); 243 /* Initial read for vertex section */ 244 ierr = PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);CHKERRQ(ierr); 245 } 246 /* Read vertices */ 247 ierr = PetscStrncmp(line, "$Nodes", PETSC_MAX_PATH_LEN, &match);CHKERRQ(ierr); 248 if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file"); 249 ierr = PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);CHKERRQ(ierr); 250 snum = sscanf(line, "%d", &numVertices); 251 if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file"); 252 ierr = PetscMalloc1(numVertices*3, &coordsIn);CHKERRQ(ierr); 253 if (binary) { 254 size_t doubleSize, intSize; 255 PetscInt elementSize; 256 char *buffer; 257 PetscScalar *baseptr; 258 ierr = PetscDataTypeGetSize(PETSC_ENUM, &intSize);CHKERRQ(ierr); 259 ierr = PetscDataTypeGetSize(PETSC_DOUBLE, &doubleSize);CHKERRQ(ierr); 260 elementSize = (intSize + 3*doubleSize); 261 ierr = PetscMalloc1(elementSize*numVertices, &buffer);CHKERRQ(ierr); 262 ierr = PetscViewerRead(viewer, buffer, elementSize*numVertices, NULL, PETSC_CHAR);CHKERRQ(ierr); 263 if (bswap) ierr = PetscByteSwap(buffer, PETSC_CHAR, elementSize*numVertices);CHKERRQ(ierr); 264 for (v = 0; v < numVertices; ++v) { 265 baseptr = ((PetscScalar*)(buffer+v*elementSize+intSize)); 266 coordsIn[v*3+0] = baseptr[0]; 267 coordsIn[v*3+1] = baseptr[1]; 268 coordsIn[v*3+2] = baseptr[2]; 269 } 270 ierr = PetscFree(buffer);CHKERRQ(ierr); 271 } else { 272 for (v = 0; v < numVertices; ++v) { 273 ierr = PetscViewerRead(viewer, &i, 1, NULL, PETSC_ENUM);CHKERRQ(ierr); 274 ierr = PetscViewerRead(viewer, &(coordsIn[v*3]), 3, NULL, PETSC_DOUBLE);CHKERRQ(ierr); 275 if (i != (int)v+1) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid node number %d should be %d", i, v+1); 276 } 277 } 278 ierr = PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);CHKERRQ(ierr);; 279 ierr = PetscStrncmp(line, "$EndNodes", PETSC_MAX_PATH_LEN, &match);CHKERRQ(ierr); 280 if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file"); 281 /* Read cells */ 282 ierr = PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);CHKERRQ(ierr);; 283 ierr = PetscStrncmp(line, "$Elements", PETSC_MAX_PATH_LEN, &match);CHKERRQ(ierr); 284 if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file"); 285 ierr = PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);CHKERRQ(ierr);; 286 snum = sscanf(line, "%d", &numCells); 287 if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file"); 288 } 289 290 if (!rank || binary) { 291 /* Gmsh elements can be of any dimension/co-dimension, so we need to traverse the 292 file contents multiple times to figure out the true number of cells and facets 293 in the given mesh. To make this more efficient we read the file contents only 294 once and store them in memory, while determining the true number of cells. */ 295 ierr = DMPlexCreateGmsh_ReadElement(viewer, numCells, binary, bswap, &gmsh_elem);CHKERRQ(ierr); 296 for (trueNumCells=0, c = 0; c < numCells; ++c) { 297 if (gmsh_elem[c].dim > dim) {dim = gmsh_elem[c].dim; trueNumCells = 0;} 298 if (gmsh_elem[c].dim == dim) trueNumCells++; 299 } 300 ierr = PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);CHKERRQ(ierr); 301 ierr = PetscStrncmp(line, "$EndElements", 12, &match);CHKERRQ(ierr); 302 if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file"); 303 if (periodic) { 304 PetscInt pVert, *periodicMapT, *aux; 305 int numPeriodic; 306 307 ierr = PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);CHKERRQ(ierr); 308 ierr = PetscStrncmp(line, "$Periodic", 9, &match);CHKERRQ(ierr); 309 if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file"); 310 ierr = PetscMalloc1(numVertices, &periodicMapT);CHKERRQ(ierr); 311 for (i = 0; i < numVertices; i++) periodicMapT[i] = i; 312 ierr = PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);CHKERRQ(ierr); 313 snum = sscanf(line, "%d", &numPeriodic); 314 if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file"); 315 for (i = 0; i < numPeriodic; i++) { 316 int j, edim = -1, slaveTag = -1, masterTag = -1, nNodes; 317 318 ierr = PetscViewerRead(viewer, line, 3, NULL, PETSC_STRING);CHKERRQ(ierr); 319 snum = sscanf(line, "%d %d %d", &edim, &slaveTag, &masterTag); /* slaveTag and masterTag are unused */ 320 if (snum != 3) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file"); 321 ierr = PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);CHKERRQ(ierr); 322 snum = sscanf(line, "%d", &nNodes); 323 if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file"); 324 for (j = 0; j < nNodes; j++) { 325 PetscInt slaveNode, masterNode; 326 327 ierr = PetscViewerRead(viewer, line, 2, NULL, PETSC_STRING);CHKERRQ(ierr); 328 snum = sscanf(line, "%d %d", &slaveNode, &masterNode); 329 if (snum != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file"); 330 periodicMapT[slaveNode - 1] = masterNode - 1; 331 } 332 } 333 ierr = PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);CHKERRQ(ierr); 334 ierr = PetscStrncmp(line, "$EndPeriodic", 12, &match);CHKERRQ(ierr); 335 if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file"); 336 /* we may have slaves of slaves */ 337 for (i = 0; i < numVertices; i++) { 338 while (periodicMapT[periodicMapT[i]] != periodicMapT[i]) { 339 periodicMapT[i] = periodicMapT[periodicMapT[i]]; 340 } 341 } 342 /* periodicMap : from old to new numbering (periodic vertices excluded) 343 periodicMapI: from new to old numbering */ 344 ierr = PetscMalloc1(numVertices, &periodicMap);CHKERRQ(ierr); 345 ierr = PetscMalloc1(numVertices, &periodicMapI);CHKERRQ(ierr); 346 ierr = PetscMalloc1(numVertices, &aux);CHKERRQ(ierr); 347 for (i = 0, pVert = 0; i < numVertices; i++) { 348 if (periodicMapT[i] != i) { 349 pVert++; 350 } else { 351 aux[i] = i - pVert; 352 periodicMapI[i - pVert] = i; 353 } 354 } 355 for (i = 0 ; i < numVertices; i++) { 356 periodicMap[i] = aux[periodicMapT[i]]; 357 } 358 ierr = PetscFree(periodicMapT);CHKERRQ(ierr); 359 ierr = PetscFree(aux);CHKERRQ(ierr); 360 /* remove periodic vertices */ 361 numVertices = numVertices - pVert; 362 } 363 } 364 /* For binary we read on all ranks, but only build the plex on rank 0 */ 365 if (binary && rank) {trueNumCells = 0; numVertices = 0;}; 366 /* Allocate the cell-vertex mesh */ 367 ierr = DMPlexSetChart(*dm, 0, trueNumCells+numVertices);CHKERRQ(ierr); 368 if (!rank) { 369 for (cell = 0, c = 0; c < numCells; ++c) { 370 if (gmsh_elem[c].dim == dim) { 371 ierr = DMPlexSetConeSize(*dm, cell, gmsh_elem[c].numNodes);CHKERRQ(ierr); 372 cell++; 373 } 374 } 375 } 376 ierr = DMSetUp(*dm);CHKERRQ(ierr); 377 /* Add cell-vertex connections */ 378 if (!rank) { 379 PetscInt pcone[8], corner; 380 for (cell = 0, c = 0; c < numCells; ++c) { 381 if (gmsh_elem[c].dim == dim) { 382 for (corner = 0; corner < gmsh_elem[c].numNodes; ++corner) { 383 pcone[corner] = periodicMap ? periodicMap[gmsh_elem[c].nodes[corner] - 1] + trueNumCells 384 : gmsh_elem[c].nodes[corner] -1 + trueNumCells; 385 } 386 if (dim == 3) { 387 /* Tetrahedra are inverted */ 388 if (gmsh_elem[c].numNodes == 4) { 389 PetscInt tmp = pcone[0]; 390 pcone[0] = pcone[1]; 391 pcone[1] = tmp; 392 } 393 /* Hexahedra are inverted */ 394 if (gmsh_elem[c].numNodes == 8) { 395 PetscInt tmp = pcone[1]; 396 pcone[1] = pcone[3]; 397 pcone[3] = tmp; 398 } 399 } 400 ierr = DMPlexSetCone(*dm, cell, pcone);CHKERRQ(ierr); 401 cell++; 402 } 403 } 404 } 405 ierr = MPI_Bcast(&dim, 1, MPIU_INT, 0, comm);CHKERRQ(ierr); 406 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 407 ierr = DMPlexSymmetrize(*dm);CHKERRQ(ierr); 408 ierr = DMPlexStratify(*dm);CHKERRQ(ierr); 409 if (interpolate) { 410 DM idm = NULL; 411 412 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 413 ierr = DMDestroy(dm);CHKERRQ(ierr); 414 *dm = idm; 415 } 416 417 if (!rank) { 418 /* Apply boundary IDs by finding the relevant facets with vertex joins */ 419 PetscInt pcone[8], corner, vStart, vEnd; 420 421 ierr = DMPlexGetDepthStratum(*dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 422 for (c = 0; c < numCells; ++c) { 423 if (gmsh_elem[c].dim == dim-1) { 424 PetscInt joinSize; 425 const PetscInt *join; 426 for (corner = 0; corner < gmsh_elem[c].numNodes; ++corner) { 427 pcone[corner] = periodicMap ? periodicMap[gmsh_elem[c].nodes[corner] - 1] + vStart 428 : gmsh_elem[c].nodes[corner] - 1 + vStart; 429 } 430 ierr = DMPlexGetFullJoin(*dm, gmsh_elem[c].numNodes, (const PetscInt *) pcone, &joinSize, &join);CHKERRQ(ierr); 431 if (joinSize != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not determine Plex facet for element %d", gmsh_elem[c].id); 432 ierr = DMSetLabelValue(*dm, "Face Sets", join[0], gmsh_elem[c].tags[0]);CHKERRQ(ierr); 433 ierr = DMPlexRestoreJoin(*dm, gmsh_elem[c].numNodes, (const PetscInt *) pcone, &joinSize, &join);CHKERRQ(ierr); 434 } 435 } 436 437 /* Create cell sets */ 438 for (cell = 0, c = 0; c < numCells; ++c) { 439 if (gmsh_elem[c].dim == dim) { 440 if (gmsh_elem[c].numTags > 0) { 441 ierr = DMSetLabelValue(*dm, "Cell Sets", cell, gmsh_elem[c].tags[0]);CHKERRQ(ierr); 442 cell++; 443 } 444 } 445 } 446 447 /* Create vertex sets */ 448 for (c = 0; c < numCells; ++c) { 449 if (gmsh_elem[c].dim == 0) { 450 if (gmsh_elem[c].numTags > 0) { 451 PetscInt vid = periodicMap ? periodicMap[gmsh_elem[c].nodes[0] - 1] + vStart 452 : gmsh_elem[c].nodes[0] - 1 + vStart; 453 ierr = DMSetLabelValue(*dm, "Vertex Sets", vid, gmsh_elem[c].tags[0]);CHKERRQ(ierr); 454 } 455 } 456 } 457 } 458 459 /* Read coordinates */ 460 ierr = DMGetCoordinateSection(*dm, &coordSection);CHKERRQ(ierr); 461 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 462 ierr = PetscSectionSetFieldComponents(coordSection, 0, dim);CHKERRQ(ierr); 463 if (periodicMap) { /* we need to localize coordinates on cells */ 464 ierr = PetscSectionSetChart(coordSection, 0, trueNumCells + numVertices);CHKERRQ(ierr); 465 } else { 466 ierr = PetscSectionSetChart(coordSection, trueNumCells, trueNumCells + numVertices);CHKERRQ(ierr); 467 } 468 for (v = trueNumCells; v < trueNumCells+numVertices; ++v) { 469 ierr = PetscSectionSetDof(coordSection, v, dim);CHKERRQ(ierr); 470 ierr = PetscSectionSetFieldDof(coordSection, v, 0, dim);CHKERRQ(ierr); 471 } 472 if (periodicMap) { 473 for (cell = 0, c = 0; c < numCells; ++c) { 474 if (gmsh_elem[c].dim == dim) { 475 ierr = PetscSectionSetDof(coordSection, cell, gmsh_elem[c].numNodes*dim);CHKERRQ(ierr); 476 ierr = PetscSectionSetFieldDof(coordSection, cell, 0, gmsh_elem[c].numNodes*dim);CHKERRQ(ierr); 477 cell++; 478 } 479 } 480 } 481 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 482 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 483 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr); 484 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 485 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 486 ierr = VecSetBlockSize(coordinates, dim);CHKERRQ(ierr); 487 ierr = VecSetType(coordinates, VECSTANDARD);CHKERRQ(ierr); 488 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 489 if (!rank) { 490 491 if (periodicMap) { 492 PetscInt off; 493 494 for (cell = 0, c = 0; c < numCells; ++c) { 495 PetscInt pcone[8], corner; 496 if (gmsh_elem[c].dim == dim) { 497 for (corner = 0; corner < gmsh_elem[c].numNodes; ++corner) { 498 pcone[corner] = gmsh_elem[c].nodes[corner] - 1; 499 } 500 if (dim == 3) { 501 /* Tetrahedra are inverted */ 502 if (gmsh_elem[c].numNodes == 4) { 503 PetscInt tmp = pcone[0]; 504 pcone[0] = pcone[1]; 505 pcone[1] = tmp; 506 } 507 /* Hexahedra are inverted */ 508 if (gmsh_elem[c].numNodes == 8) { 509 PetscInt tmp = pcone[1]; 510 pcone[1] = pcone[3]; 511 pcone[3] = tmp; 512 } 513 } 514 ierr = PetscSectionGetOffset(coordSection, cell, &off);CHKERRQ(ierr); 515 for (corner = 0; corner < gmsh_elem[c].numNodes; ++corner) { 516 PetscInt v = pcone[corner]; 517 for (d = 0; d < dim; ++d) { 518 coords[off++] = coordsIn[v*3+d]; 519 } 520 } 521 cell++; 522 } 523 } 524 for (v = 0; v < numVertices; ++v) { 525 ierr = PetscSectionGetOffset(coordSection, v + trueNumCells, &off);CHKERRQ(ierr); 526 for (d = 0; d < dim; ++d) { 527 coords[off+d] = coordsIn[periodicMapI[v]*3+d]; 528 } 529 } 530 } else { 531 for (v = 0; v < numVertices; ++v) { 532 for (d = 0; d < dim; ++d) { 533 coords[v*dim+d] = coordsIn[v*3+d]; 534 } 535 } 536 } 537 } 538 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 539 ierr = PetscFree(coordsIn);CHKERRQ(ierr); 540 ierr = DMSetCoordinatesLocal(*dm, coordinates);CHKERRQ(ierr); 541 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 542 ierr = PetscFree(periodicMap);CHKERRQ(ierr); 543 /* Clean up intermediate storage */ 544 if (!rank || binary) ierr = PetscFree(gmsh_elem);CHKERRQ(ierr); 545 ierr = PetscLogEventEnd(DMPLEX_CreateGmsh,*dm,0,0,0);CHKERRQ(ierr); 546 PetscFunctionReturn(0); 547 } 548