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 83 Level: beginner 84 85 .keywords: mesh,Gmsh 86 .seealso: DMPLEX, DMCreate() 87 @*/ 88 PetscErrorCode DMPlexCreateGmsh(MPI_Comm comm, PetscViewer viewer, PetscBool interpolate, DM *dm) 89 { 90 PetscViewerType vtype; 91 GmshElement *gmsh_elem; 92 PetscSection coordSection; 93 Vec coordinates; 94 PetscScalar *coords, *coordsIn = NULL; 95 PetscInt dim = 0, coordSize, c, v, d, cell; 96 int numVertices = 0, numCells = 0, trueNumCells = 0, snum; 97 PetscMPIInt num_proc, rank; 98 char line[PETSC_MAX_PATH_LEN]; 99 PetscBool match, binary, bswap = PETSC_FALSE; 100 PetscErrorCode ierr; 101 102 PetscFunctionBegin; 103 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 104 ierr = MPI_Comm_size(comm, &num_proc);CHKERRQ(ierr); 105 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 106 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 107 ierr = PetscViewerGetType(viewer, &vtype);CHKERRQ(ierr); 108 ierr = PetscStrcmp(vtype, PETSCVIEWERBINARY, &binary);CHKERRQ(ierr); 109 if (!rank || binary) { 110 PetscBool match; 111 int fileType, dataSize; 112 float version; 113 114 /* Read header */ 115 ierr = PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);CHKERRQ(ierr); 116 ierr = PetscStrncmp(line, "$MeshFormat", PETSC_MAX_PATH_LEN, &match);CHKERRQ(ierr); 117 if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file"); 118 ierr = PetscViewerRead(viewer, line, 3, NULL, PETSC_STRING);CHKERRQ(ierr); 119 snum = sscanf(line, "%f %d %d", &version, &fileType, &dataSize); 120 if (snum != 3) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse Gmsh file header: %s", line); 121 if (version < 2.0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Gmsh file must be at least version 2.0"); 122 if (dataSize != sizeof(double)) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Data size %d is not valid for a Gmsh file", dataSize); 123 if (binary) { 124 int checkInt; 125 ierr = PetscViewerRead(viewer, &checkInt, 1, NULL, PETSC_ENUM);CHKERRQ(ierr); 126 if (checkInt != 1) { 127 ierr = PetscByteSwap(&checkInt, PETSC_ENUM, 1);CHKERRQ(ierr); 128 if (checkInt == 1) bswap = PETSC_TRUE; 129 else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File type %d is not a valid Gmsh binary file", fileType); 130 } 131 } else if (fileType) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File type %d is not a valid Gmsh ASCII file", fileType); 132 ierr = PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);CHKERRQ(ierr); 133 ierr = PetscStrncmp(line, "$EndMeshFormat", PETSC_MAX_PATH_LEN, &match);CHKERRQ(ierr); 134 if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file"); 135 /* Read vertices */ 136 ierr = PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);CHKERRQ(ierr); 137 ierr = PetscStrncmp(line, "$Nodes", PETSC_MAX_PATH_LEN, &match);CHKERRQ(ierr); 138 if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file"); 139 ierr = PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);CHKERRQ(ierr); 140 snum = sscanf(line, "%d", &numVertices); 141 if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file"); 142 ierr = PetscMalloc1(numVertices*3, &coordsIn);CHKERRQ(ierr); 143 for (v = 0; v < numVertices; ++v) { 144 int i; 145 ierr = PetscViewerRead(viewer, &i, 1, NULL, PETSC_ENUM);CHKERRQ(ierr); 146 if (bswap) ierr = PetscByteSwap(&i, PETSC_ENUM, 1);CHKERRQ(ierr); 147 ierr = PetscViewerRead(viewer, &(coordsIn[v*3]), 3, NULL, PETSC_DOUBLE);CHKERRQ(ierr); 148 if (bswap) ierr = PetscByteSwap(&(coordsIn[v*3]), PETSC_DOUBLE, 3);CHKERRQ(ierr); 149 if (i != (int)v+1) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid node number %d should be %d", i, v+1); 150 } 151 ierr = PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);CHKERRQ(ierr);; 152 ierr = PetscStrncmp(line, "$EndNodes", 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 /* Read cells */ 155 ierr = PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);CHKERRQ(ierr);; 156 ierr = PetscStrncmp(line, "$Elements", PETSC_MAX_PATH_LEN, &match);CHKERRQ(ierr); 157 if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file"); 158 ierr = PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);CHKERRQ(ierr);; 159 snum = sscanf(line, "%d", &numCells); 160 if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file"); 161 } 162 163 if (!rank || binary) { 164 /* Gmsh elements can be of any dimension/co-dimension, so we need to traverse the 165 file contents multiple times to figure out the true number of cells and facets 166 in the given mesh. To make this more efficient we read the file contents only 167 once and store them in memory, while determining the true number of cells. */ 168 ierr = PetscMalloc1(numCells, &gmsh_elem);CHKERRQ(ierr); 169 for (trueNumCells=0, c = 0; c < numCells; ++c) { 170 ierr = DMPlexCreateGmsh_ReadElement(viewer, binary, bswap, &gmsh_elem[c]);CHKERRQ(ierr); 171 if (gmsh_elem[c].dim > dim) {dim = gmsh_elem[c].dim; trueNumCells = 0;} 172 if (gmsh_elem[c].dim == dim) trueNumCells++; 173 } 174 ierr = PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);CHKERRQ(ierr);; 175 ierr = PetscStrncmp(line, "$EndElements", 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 } 178 /* For binary we read on all ranks, but only build the plex on rank 0 */ 179 if (binary && rank) {trueNumCells = 0; numVertices = 0;}; 180 /* Allocate the cell-vertex mesh */ 181 ierr = DMPlexSetChart(*dm, 0, trueNumCells+numVertices);CHKERRQ(ierr); 182 if (!rank) { 183 for (cell = 0, c = 0; c < numCells; ++c) { 184 if (gmsh_elem[c].dim == dim) { 185 ierr = DMPlexSetConeSize(*dm, cell, gmsh_elem[c].numNodes);CHKERRQ(ierr); 186 cell++; 187 } 188 } 189 } 190 ierr = DMSetUp(*dm);CHKERRQ(ierr); 191 /* Add cell-vertex connections */ 192 if (!rank) { 193 PetscInt pcone[8], corner; 194 for (cell = 0, c = 0; c < numCells; ++c) { 195 if (gmsh_elem[c].dim == dim) { 196 for (corner = 0; corner < gmsh_elem[c].numNodes; ++corner) { 197 pcone[corner] = gmsh_elem[c].nodes[corner] + trueNumCells-1; 198 } 199 ierr = DMPlexSetCone(*dm, cell, pcone);CHKERRQ(ierr); 200 cell++; 201 } 202 } 203 } 204 ierr = MPI_Bcast(&dim, 1, MPIU_INT, 0, comm);CHKERRQ(ierr); 205 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 206 ierr = DMPlexSymmetrize(*dm);CHKERRQ(ierr); 207 ierr = DMPlexStratify(*dm);CHKERRQ(ierr); 208 if (interpolate) { 209 DM idm = NULL; 210 211 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 212 ierr = DMDestroy(dm);CHKERRQ(ierr); 213 *dm = idm; 214 } 215 216 if (!rank) { 217 /* Apply boundary IDs by finding the relevant facets with vertex joins */ 218 PetscInt pcone[8], corner, vStart, vEnd; 219 220 ierr = DMPlexGetDepthStratum(*dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 221 for (c = 0; c < numCells; ++c) { 222 if (gmsh_elem[c].dim == dim-1) { 223 PetscInt joinSize; 224 const PetscInt *join; 225 for (corner = 0; corner < gmsh_elem[c].numNodes; ++corner) { 226 pcone[corner] = gmsh_elem[c].nodes[corner] + vStart - 1; 227 } 228 ierr = DMPlexGetFullJoin(*dm, gmsh_elem[c].numNodes, (const PetscInt *) pcone, &joinSize, &join);CHKERRQ(ierr); 229 if (joinSize != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not determine Plex facet for element %d", gmsh_elem[c].id); 230 ierr = DMPlexSetLabelValue(*dm, "Face Sets", join[0], gmsh_elem[c].tags[0]);CHKERRQ(ierr); 231 ierr = DMPlexRestoreJoin(*dm, gmsh_elem[c].numNodes, (const PetscInt *) pcone, &joinSize, &join);CHKERRQ(ierr); 232 } 233 } 234 } 235 236 /* Read coordinates */ 237 ierr = DMGetCoordinateSection(*dm, &coordSection);CHKERRQ(ierr); 238 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 239 ierr = PetscSectionSetFieldComponents(coordSection, 0, dim);CHKERRQ(ierr); 240 ierr = PetscSectionSetChart(coordSection, trueNumCells, trueNumCells + numVertices);CHKERRQ(ierr); 241 for (v = trueNumCells; v < trueNumCells+numVertices; ++v) { 242 ierr = PetscSectionSetDof(coordSection, v, dim);CHKERRQ(ierr); 243 ierr = PetscSectionSetFieldDof(coordSection, v, 0, dim);CHKERRQ(ierr); 244 } 245 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 246 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 247 ierr = VecCreate(comm, &coordinates);CHKERRQ(ierr); 248 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 249 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 250 ierr = VecSetType(coordinates, VECSTANDARD);CHKERRQ(ierr); 251 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 252 if (!rank) { 253 for (v = 0; v < numVertices; ++v) { 254 for (d = 0; d < dim; ++d) { 255 coords[v*dim+d] = coordsIn[v*3+d]; 256 } 257 } 258 } 259 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 260 ierr = PetscFree(coordsIn);CHKERRQ(ierr); 261 ierr = DMSetCoordinatesLocal(*dm, coordinates);CHKERRQ(ierr); 262 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 263 /* Clean up intermediate storage */ 264 if (!rank || binary) { 265 for (c = 0; c < numCells; ++c) { 266 ierr = PetscFree(gmsh_elem[c].nodes); 267 ierr = PetscFree(gmsh_elem[c].tags); 268 } 269 ierr = PetscFree(gmsh_elem); 270 } 271 PetscFunctionReturn(0); 272 } 273 274 #undef __FUNCT__ 275 #define __FUNCT__ "DMPlexCreateGmsh_ReadElement" 276 PetscErrorCode DMPlexCreateGmsh_ReadElement(PetscViewer viewer, PetscBool binary, PetscBool byteSwap, GmshElement *ele) 277 { 278 int cellType, numElem; 279 PetscErrorCode ierr; 280 281 PetscFunctionBegin; 282 if (binary) { 283 ierr = PetscViewerRead(viewer, &cellType, 1, NULL, PETSC_ENUM);CHKERRQ(ierr); 284 if (byteSwap) ierr = PetscByteSwap(&cellType, PETSC_ENUM, 1);CHKERRQ(ierr); 285 ierr = PetscViewerRead(viewer, &numElem, 1, NULL, PETSC_ENUM);CHKERRQ(ierr); 286 if (byteSwap) ierr = PetscByteSwap(&numElem, PETSC_ENUM, 1);CHKERRQ(ierr); 287 ierr = PetscViewerRead(viewer, &(ele->numTags), 1, NULL, PETSC_ENUM);CHKERRQ(ierr); 288 if (byteSwap) ierr = PetscByteSwap(&(ele->numTags), PETSC_ENUM, 1);CHKERRQ(ierr); 289 ierr = PetscViewerRead(viewer, &(ele->id), 1, NULL, PETSC_ENUM);CHKERRQ(ierr); 290 if (byteSwap) ierr = PetscByteSwap( &(ele->id), PETSC_ENUM, 1);CHKERRQ(ierr); 291 } else { 292 ierr = PetscViewerRead(viewer, &(ele->id), 1, NULL, PETSC_ENUM);CHKERRQ(ierr); 293 ierr = PetscViewerRead(viewer, &cellType, 1, NULL, PETSC_ENUM);CHKERRQ(ierr); 294 ierr = PetscViewerRead(viewer, &(ele->numTags), 1, NULL, PETSC_ENUM);CHKERRQ(ierr); 295 } 296 ierr = PetscMalloc1(ele->numTags, &(ele->tags));CHKERRQ(ierr); 297 ierr = PetscViewerRead(viewer, ele->tags, ele->numTags, NULL, PETSC_ENUM);CHKERRQ(ierr); 298 if (byteSwap) ierr = PetscByteSwap(ele->tags, PETSC_ENUM, ele->numTags);CHKERRQ(ierr); 299 switch (cellType) { 300 case 1: /* 2-node line */ 301 ele->dim = 1; 302 ele->numNodes = 2; 303 break; 304 case 2: /* 3-node triangle */ 305 ele->dim = 2; 306 ele->numNodes = 3; 307 break; 308 case 3: /* 4-node quadrangle */ 309 ele->dim = 2; 310 ele->numNodes = 4; 311 break; 312 case 4: /* 4-node tetrahedron */ 313 ele->dim = 3; 314 ele->numNodes = 4; 315 break; 316 case 5: /* 8-node hexahedron */ 317 ele->dim = 3; 318 ele->numNodes = 8; 319 break; 320 case 15: /* 1-node vertex */ 321 ele->dim = 0; 322 ele->numNodes = 1; 323 break; 324 default: 325 SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unsupported Gmsh element type %d", cellType); 326 } 327 ierr = PetscMalloc1(ele->numNodes, &(ele->nodes));CHKERRQ(ierr); 328 ierr = PetscViewerRead(viewer, ele->nodes, ele->numNodes, NULL, PETSC_ENUM);CHKERRQ(ierr); 329 if (byteSwap) {ierr = PetscByteSwap(ele->nodes, PETSC_ENUM, ele->numNodes);CHKERRQ(ierr);} 330 PetscFunctionReturn(0); 331 } 332