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