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 PetscInt n = 1, m = 2; 42 /* Read only the first two lines of the Gmsh file */ 43 ierr = PetscViewerRead(vheader, line, &n, PETSC_STRING);CHKERRQ(ierr); 44 ierr = PetscStrncmp(line, "$MeshFormat", PETSC_MAX_PATH_LEN, &match);CHKERRQ(ierr); 45 if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file"); 46 ierr = PetscViewerRead(vheader, line, &m, PETSC_STRING);CHKERRQ(ierr); 47 snum = sscanf(line, "%f %d", &version, &fT); 48 fileType = (PetscInt) fT; 49 if (snum != 2) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse Gmsh file header: %s", line); 50 if (version < 2.0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Gmsh file must be at least version 2.0"); 51 } 52 ierr = MPI_Bcast(&fileType, 1, MPIU_INT, 0, comm);CHKERRQ(ierr); 53 /* Create appropriate viewer and build plex */ 54 if (fileType == 0) vtype = PETSCVIEWERASCII; 55 else vtype = PETSCVIEWERBINARY; 56 ierr = PetscViewerCreate(comm, &viewer);CHKERRQ(ierr); 57 ierr = PetscViewerSetType(viewer, vtype);CHKERRQ(ierr); 58 ierr = PetscViewerFileSetMode(viewer, FILE_MODE_READ);CHKERRQ(ierr); 59 ierr = PetscViewerFileSetName(viewer, filename);CHKERRQ(ierr); 60 ierr = DMPlexCreateGmsh(comm, viewer, interpolate, dm);CHKERRQ(ierr); 61 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 62 ierr = PetscViewerDestroy(&vheader);CHKERRQ(ierr); 63 PetscFunctionReturn(0); 64 } 65 66 67 #undef __FUNCT__ 68 #define __FUNCT__ "DMPlexCreateGmsh" 69 /*@ 70 DMPlexCreateGmsh - Create a DMPlex mesh from a Gmsh file viewer 71 72 Collective on comm 73 74 Input Parameters: 75 + comm - The MPI communicator 76 . viewer - The Viewer associated with a Gmsh file 77 - interpolate - Create faces and edges in the mesh 78 79 Output Parameter: 80 . dm - The DM object representing the mesh 81 82 Note: http://www.geuz.org/gmsh/doc/texinfo/#MSH-ASCII-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 = 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 PetscInt n = 1, m = 3; 115 116 /* Read header */ 117 ierr = PetscViewerRead(viewer, line, &n, 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, &m, 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, &n, 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, &n, 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, &n, 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, &n, 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, &n, PETSC_ENUM);CHKERRQ(ierr); 148 if (bswap) ierr = PetscByteSwap(&i, PETSC_ENUM, 1);CHKERRQ(ierr); 149 ierr = PetscViewerRead(viewer, &(coordsIn[v*3]), &m, 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, &n, 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, &n, 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, &n, 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 PetscInt n = 1; 167 /* Gmsh elements can be of any dimension/co-dimension, so we need to traverse the 168 file contents multiple times to figure out the true number of cells and facets 169 in the given mesh. To make this more efficient we read the file contents only 170 once and store them in memory, while determining the true number of cells. */ 171 ierr = PetscMalloc1(numCells, &gmsh_elem);CHKERRQ(ierr); 172 for (trueNumCells=0, c = 0; c < numCells; ++c) { 173 ierr = DMPlexCreateGmsh_ReadElement(viewer, binary, bswap, &gmsh_elem[c]);CHKERRQ(ierr); 174 if (gmsh_elem[c].dim > dim) {dim = gmsh_elem[c].dim; trueNumCells = 0;} 175 if (gmsh_elem[c].dim == dim) trueNumCells++; 176 } 177 ierr = PetscViewerRead(viewer, line, &n, PETSC_STRING);CHKERRQ(ierr);; 178 ierr = PetscStrncmp(line, "$EndElements", PETSC_MAX_PATH_LEN, &match);CHKERRQ(ierr); 179 if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file"); 180 } 181 /* For binary we read on all ranks, but only build the plex on rank 0 */ 182 if (binary && rank) {trueNumCells = 0; numVertices = 0;}; 183 /* Allocate the cell-vertex mesh */ 184 ierr = DMPlexSetChart(*dm, 0, trueNumCells+numVertices);CHKERRQ(ierr); 185 if (!rank) { 186 for (cell = 0, c = 0; c < numCells; ++c) { 187 if (gmsh_elem[c].dim == dim) { 188 ierr = DMPlexSetConeSize(*dm, cell, gmsh_elem[c].numNodes);CHKERRQ(ierr); 189 cell++; 190 } 191 } 192 } 193 ierr = DMSetUp(*dm);CHKERRQ(ierr); 194 /* Add cell-vertex connections */ 195 if (!rank) { 196 PetscInt pcone[8], corner; 197 for (cell = 0, c = 0; c < numCells; ++c) { 198 if (gmsh_elem[c].dim == dim) { 199 for (corner = 0; corner < gmsh_elem[c].numNodes; ++corner) { 200 pcone[corner] = gmsh_elem[c].nodes[corner] + trueNumCells-1; 201 } 202 ierr = DMPlexSetCone(*dm, cell, pcone);CHKERRQ(ierr); 203 cell++; 204 } 205 } 206 } 207 ierr = MPI_Bcast(&dim, 1, MPIU_INT, 0, comm);CHKERRQ(ierr); 208 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 209 ierr = DMPlexSymmetrize(*dm);CHKERRQ(ierr); 210 ierr = DMPlexStratify(*dm);CHKERRQ(ierr); 211 if (interpolate) { 212 DM idm = NULL; 213 214 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 215 ierr = DMDestroy(dm);CHKERRQ(ierr); 216 *dm = idm; 217 } 218 219 if (!rank) { 220 /* Apply boundary IDs by finding the relevant facets with vertex joins */ 221 PetscInt pcone[8], corner, vStart, vEnd; 222 223 ierr = DMPlexGetDepthStratum(*dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 224 for (c = 0; c < numCells; ++c) { 225 if (gmsh_elem[c].dim == dim-1) { 226 PetscInt joinSize; 227 const PetscInt *join; 228 for (corner = 0; corner < gmsh_elem[c].numNodes; ++corner) { 229 pcone[corner] = gmsh_elem[c].nodes[corner] + vStart - 1; 230 } 231 ierr = DMPlexGetFullJoin(*dm, gmsh_elem[c].numNodes, (const PetscInt *) pcone, &joinSize, &join);CHKERRQ(ierr); 232 if (joinSize != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not determine Plex facet for element %d", gmsh_elem[c].id); 233 ierr = DMPlexSetLabelValue(*dm, "Face Sets", join[0], gmsh_elem[c].tags[0]);CHKERRQ(ierr); 234 ierr = DMPlexRestoreJoin(*dm, gmsh_elem[c].numNodes, (const PetscInt *) pcone, &joinSize, &join);CHKERRQ(ierr); 235 } 236 } 237 } 238 239 /* Read coordinates */ 240 ierr = DMGetCoordinateSection(*dm, &coordSection);CHKERRQ(ierr); 241 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 242 ierr = PetscSectionSetFieldComponents(coordSection, 0, dim);CHKERRQ(ierr); 243 ierr = PetscSectionSetChart(coordSection, trueNumCells, trueNumCells + numVertices);CHKERRQ(ierr); 244 for (v = trueNumCells; v < trueNumCells+numVertices; ++v) { 245 ierr = PetscSectionSetDof(coordSection, v, dim);CHKERRQ(ierr); 246 ierr = PetscSectionSetFieldDof(coordSection, v, 0, dim);CHKERRQ(ierr); 247 } 248 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 249 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 250 ierr = VecCreate(comm, &coordinates);CHKERRQ(ierr); 251 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 252 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 253 ierr = VecSetType(coordinates, VECSTANDARD);CHKERRQ(ierr); 254 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 255 if (!rank) { 256 for (v = 0; v < numVertices; ++v) { 257 for (d = 0; d < dim; ++d) { 258 coords[v*dim+d] = coordsIn[v*3+d]; 259 } 260 } 261 } 262 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 263 ierr = PetscFree(coordsIn);CHKERRQ(ierr); 264 ierr = DMSetCoordinatesLocal(*dm, coordinates);CHKERRQ(ierr); 265 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 266 /* Clean up intermediate storage */ 267 if (!rank || binary) { 268 for (c = 0; c < numCells; ++c) { 269 ierr = PetscFree(gmsh_elem[c].nodes); 270 ierr = PetscFree(gmsh_elem[c].tags); 271 } 272 ierr = PetscFree(gmsh_elem); 273 } 274 PetscFunctionReturn(0); 275 } 276 277 #undef __FUNCT__ 278 #define __FUNCT__ "DMPlexCreateGmsh_ReadElement" 279 PetscErrorCode DMPlexCreateGmsh_ReadElement(PetscViewer viewer, PetscBool binary, PetscBool byteSwap, GmshElement *ele) 280 { 281 int cellType, numElem; 282 PetscInt n = 1; 283 PetscErrorCode ierr; 284 285 PetscFunctionBegin; 286 if (binary) { 287 ierr = PetscViewerRead(viewer, &cellType, &n, PETSC_ENUM);CHKERRQ(ierr); 288 if (byteSwap) ierr = PetscByteSwap(&cellType, PETSC_ENUM, 1);CHKERRQ(ierr); 289 ierr = PetscViewerRead(viewer, &numElem, &n, PETSC_ENUM);CHKERRQ(ierr); 290 if (byteSwap) ierr = PetscByteSwap(&numElem, PETSC_ENUM, 1);CHKERRQ(ierr); 291 ierr = PetscViewerRead(viewer, &(ele->numTags), &n, PETSC_ENUM);CHKERRQ(ierr); 292 if (byteSwap) ierr = PetscByteSwap(&(ele->numTags), PETSC_ENUM, 1);CHKERRQ(ierr); 293 ierr = PetscViewerRead(viewer, &(ele->id), &n, PETSC_ENUM);CHKERRQ(ierr); 294 if (byteSwap) ierr = PetscByteSwap( &(ele->id), PETSC_ENUM, 1);CHKERRQ(ierr); 295 } else { 296 ierr = PetscViewerRead(viewer, &(ele->id), &n, PETSC_ENUM);CHKERRQ(ierr); 297 ierr = PetscViewerRead(viewer, &cellType, &n, PETSC_ENUM);CHKERRQ(ierr); 298 ierr = PetscViewerRead(viewer, &(ele->numTags), &n, PETSC_ENUM);CHKERRQ(ierr); 299 } 300 ierr = PetscMalloc1(ele->numTags, &(ele->tags));CHKERRQ(ierr); 301 ierr = PetscViewerRead(viewer, ele->tags, &ele->numTags, PETSC_ENUM);CHKERRQ(ierr); 302 if (byteSwap) ierr = PetscByteSwap(ele->tags, PETSC_ENUM, ele->numTags);CHKERRQ(ierr); 303 switch (cellType) { 304 case 1: /* 2-node line */ 305 ele->dim = 1; 306 ele->numNodes = 2; 307 break; 308 case 2: /* 3-node triangle */ 309 ele->dim = 2; 310 ele->numNodes = 3; 311 break; 312 case 3: /* 4-node quadrangle */ 313 ele->dim = 2; 314 ele->numNodes = 4; 315 break; 316 case 4: /* 4-node tetrahedron */ 317 ele->dim = 3; 318 ele->numNodes = 4; 319 break; 320 case 5: /* 8-node hexahedron */ 321 ele->dim = 3; 322 ele->numNodes = 8; 323 break; 324 case 15: /* 1-node vertex */ 325 ele->dim = 0; 326 ele->numNodes = 1; 327 break; 328 default: 329 SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unsupported Gmsh element type %d", cellType); 330 } 331 ierr = PetscMalloc1(ele->numNodes, &(ele->nodes));CHKERRQ(ierr); 332 ierr = PetscViewerRead(viewer, ele->nodes, &ele->numNodes, PETSC_ENUM);CHKERRQ(ierr); 333 if (byteSwap) {ierr = PetscByteSwap(ele->nodes, PETSC_ENUM, ele->numNodes);CHKERRQ(ierr);} 334 PetscFunctionReturn(0); 335 } 336