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