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