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 = PetscLogEventBegin(DMPLEX_CreateGmsh,*dm,0,0,0);CHKERRQ(ierr); 108 ierr = PetscViewerGetType(viewer, &vtype);CHKERRQ(ierr); 109 ierr = PetscStrcmp(vtype, PETSCVIEWERBINARY, &binary);CHKERRQ(ierr); 110 if (!rank || binary) { 111 PetscBool match; 112 int fileType, dataSize; 113 float version; 114 115 /* Read header */ 116 ierr = PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);CHKERRQ(ierr); 117 ierr = PetscStrncmp(line, "$MeshFormat", PETSC_MAX_PATH_LEN, &match);CHKERRQ(ierr); 118 if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file"); 119 ierr = PetscViewerRead(viewer, line, 3, NULL, PETSC_STRING);CHKERRQ(ierr); 120 snum = sscanf(line, "%f %d %d", &version, &fileType, &dataSize); 121 if (snum != 3) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse Gmsh file header: %s", line); 122 if (version < 2.0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Gmsh file must be at least version 2.0"); 123 if (dataSize != sizeof(double)) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Data size %d is not valid for a Gmsh file", dataSize); 124 if (binary) { 125 int checkInt; 126 ierr = PetscViewerRead(viewer, &checkInt, 1, NULL, PETSC_ENUM);CHKERRQ(ierr); 127 if (checkInt != 1) { 128 ierr = PetscByteSwap(&checkInt, PETSC_ENUM, 1);CHKERRQ(ierr); 129 if (checkInt == 1) bswap = PETSC_TRUE; 130 else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File type %d is not a valid Gmsh binary file", fileType); 131 } 132 } else if (fileType) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File type %d is not a valid Gmsh ASCII file", fileType); 133 ierr = PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);CHKERRQ(ierr); 134 ierr = PetscStrncmp(line, "$EndMeshFormat", PETSC_MAX_PATH_LEN, &match);CHKERRQ(ierr); 135 if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file"); 136 /* Read vertices */ 137 ierr = PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);CHKERRQ(ierr); 138 ierr = PetscStrncmp(line, "$Nodes", PETSC_MAX_PATH_LEN, &match);CHKERRQ(ierr); 139 if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file"); 140 ierr = PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);CHKERRQ(ierr); 141 snum = sscanf(line, "%d", &numVertices); 142 if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file"); 143 ierr = PetscMalloc1(numVertices*3, &coordsIn);CHKERRQ(ierr); 144 for (v = 0; v < numVertices; ++v) { 145 int i; 146 ierr = PetscViewerRead(viewer, &i, 1, NULL, PETSC_ENUM);CHKERRQ(ierr); 147 if (bswap) ierr = PetscByteSwap(&i, PETSC_ENUM, 1);CHKERRQ(ierr); 148 ierr = PetscViewerRead(viewer, &(coordsIn[v*3]), 3, NULL, PETSC_DOUBLE);CHKERRQ(ierr); 149 if (bswap) ierr = PetscByteSwap(&(coordsIn[v*3]), PETSC_DOUBLE, 3);CHKERRQ(ierr); 150 if (i != (int)v+1) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid node number %d should be %d", i, v+1); 151 } 152 ierr = PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);CHKERRQ(ierr);; 153 ierr = PetscStrncmp(line, "$EndNodes", PETSC_MAX_PATH_LEN, &match);CHKERRQ(ierr); 154 if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file"); 155 /* Read cells */ 156 ierr = PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);CHKERRQ(ierr);; 157 ierr = PetscStrncmp(line, "$Elements", PETSC_MAX_PATH_LEN, &match);CHKERRQ(ierr); 158 if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file"); 159 ierr = PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);CHKERRQ(ierr);; 160 snum = sscanf(line, "%d", &numCells); 161 if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file"); 162 } 163 164 if (!rank || binary) { 165 /* Gmsh elements can be of any dimension/co-dimension, so we need to traverse the 166 file contents multiple times to figure out the true number of cells and facets 167 in the given mesh. To make this more efficient we read the file contents only 168 once and store them in memory, while determining the true number of cells. */ 169 ierr = PetscMalloc1(numCells, &gmsh_elem);CHKERRQ(ierr); 170 for (trueNumCells=0, c = 0; c < numCells; ++c) { 171 ierr = DMPlexCreateGmsh_ReadElement(viewer, binary, bswap, &gmsh_elem[c]);CHKERRQ(ierr); 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) { 266 for (c = 0; c < numCells; ++c) { 267 ierr = PetscFree(gmsh_elem[c].nodes);CHKERRQ(ierr); 268 ierr = PetscFree(gmsh_elem[c].tags);CHKERRQ(ierr); 269 } 270 ierr = PetscFree(gmsh_elem);CHKERRQ(ierr); 271 } 272 ierr = PetscLogEventEnd(DMPLEX_CreateGmsh,*dm,0,0,0);CHKERRQ(ierr); 273 PetscFunctionReturn(0); 274 } 275 276 #undef __FUNCT__ 277 #define __FUNCT__ "DMPlexCreateGmsh_ReadElement" 278 PetscErrorCode DMPlexCreateGmsh_ReadElement(PetscViewer viewer, PetscBool binary, PetscBool byteSwap, GmshElement *ele) 279 { 280 int cellType, numElem; 281 PetscErrorCode ierr; 282 283 PetscFunctionBegin; 284 if (binary) { 285 ierr = PetscViewerRead(viewer, &cellType, 1, NULL, PETSC_ENUM);CHKERRQ(ierr); 286 if (byteSwap) ierr = PetscByteSwap(&cellType, PETSC_ENUM, 1);CHKERRQ(ierr); 287 ierr = PetscViewerRead(viewer, &numElem, 1, NULL, PETSC_ENUM);CHKERRQ(ierr); 288 if (byteSwap) ierr = PetscByteSwap(&numElem, PETSC_ENUM, 1);CHKERRQ(ierr); 289 ierr = PetscViewerRead(viewer, &(ele->numTags), 1, NULL, PETSC_ENUM);CHKERRQ(ierr); 290 if (byteSwap) ierr = PetscByteSwap(&(ele->numTags), PETSC_ENUM, 1);CHKERRQ(ierr); 291 ierr = PetscViewerRead(viewer, &(ele->id), 1, NULL, PETSC_ENUM);CHKERRQ(ierr); 292 if (byteSwap) ierr = PetscByteSwap( &(ele->id), PETSC_ENUM, 1);CHKERRQ(ierr); 293 } else { 294 ierr = PetscViewerRead(viewer, &(ele->id), 1, NULL, PETSC_ENUM);CHKERRQ(ierr); 295 ierr = PetscViewerRead(viewer, &cellType, 1, NULL, PETSC_ENUM);CHKERRQ(ierr); 296 ierr = PetscViewerRead(viewer, &(ele->numTags), 1, NULL, PETSC_ENUM);CHKERRQ(ierr); 297 } 298 ierr = PetscMalloc1(ele->numTags, &(ele->tags));CHKERRQ(ierr); 299 ierr = PetscViewerRead(viewer, ele->tags, ele->numTags, NULL, PETSC_ENUM);CHKERRQ(ierr); 300 if (byteSwap) ierr = PetscByteSwap(ele->tags, PETSC_ENUM, ele->numTags);CHKERRQ(ierr); 301 switch (cellType) { 302 case 1: /* 2-node line */ 303 ele->dim = 1; 304 ele->numNodes = 2; 305 break; 306 case 2: /* 3-node triangle */ 307 ele->dim = 2; 308 ele->numNodes = 3; 309 break; 310 case 3: /* 4-node quadrangle */ 311 ele->dim = 2; 312 ele->numNodes = 4; 313 break; 314 case 4: /* 4-node tetrahedron */ 315 ele->dim = 3; 316 ele->numNodes = 4; 317 break; 318 case 5: /* 8-node hexahedron */ 319 ele->dim = 3; 320 ele->numNodes = 8; 321 break; 322 case 15: /* 1-node vertex */ 323 ele->dim = 0; 324 ele->numNodes = 1; 325 break; 326 default: 327 SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unsupported Gmsh element type %d", cellType); 328 } 329 ierr = PetscMalloc1(ele->numNodes, &(ele->nodes));CHKERRQ(ierr); 330 ierr = PetscViewerRead(viewer, ele->nodes, ele->numNodes, NULL, PETSC_ENUM);CHKERRQ(ierr); 331 if (byteSwap) {ierr = PetscByteSwap(ele->nodes, PETSC_ENUM, ele->numNodes);CHKERRQ(ierr);} 332 PetscFunctionReturn(0); 333 } 334