1 #define PETSCDM_DLL 2 #include <petsc-private/dmpleximpl.h> /*I "petscdmplex.h" I*/ 3 4 #undef __FUNCT__ 5 #define __FUNCT__ "DMPlexCreateGmsh" 6 /*@ 7 DMPlexCreateGmsh - Create a DMPlex mesh from a Gmsh file. 8 9 Collective on comm 10 11 Input Parameters: 12 + comm - The MPI communicator 13 . viewer - The Viewer associated with a Gmsh file 14 - interpolate - Create faces and edges in the mesh 15 16 Output Parameter: 17 . dm - The DM object representing the mesh 18 19 Note: http://www.geuz.org/gmsh/doc/texinfo/#MSH-ASCII-file-format 20 21 Level: beginner 22 23 .keywords: mesh,Gmsh 24 .seealso: DMPLEX, DMCreate() 25 @*/ 26 PetscErrorCode DMPlexCreateGmsh(MPI_Comm comm, PetscViewer viewer, PetscBool interpolate, DM *dm) 27 { 28 FILE *fd; 29 GmshElement *gmsh_elem; 30 PetscSection coordSection; 31 Vec coordinates; 32 PetscScalar *coords, *coordsIn = NULL; 33 PetscInt dim = 0, coordSize, c, v, d, cell; 34 int numVertices = 0, numCells = 0, trueNumCells = 0, snum; 35 PetscMPIInt num_proc, rank; 36 char line[PETSC_MAX_PATH_LEN]; 37 PetscBool match; 38 PetscErrorCode ierr; 39 40 PetscFunctionBegin; 41 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 42 ierr = MPI_Comm_size(comm, &num_proc);CHKERRQ(ierr); 43 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 44 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 45 if (!rank) { 46 PetscBool match; 47 int fileType, dataSize; 48 49 ierr = PetscViewerASCIIGetPointer(viewer, &fd);CHKERRQ(ierr); 50 /* Read header */ 51 fgets(line, PETSC_MAX_PATH_LEN, fd); 52 ierr = PetscStrncmp(line, "$MeshFormat\n", PETSC_MAX_PATH_LEN, &match);CHKERRQ(ierr); 53 if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file"); 54 snum = fscanf(fd, "2.2 %d %d\n", &fileType, &dataSize);CHKERRQ(snum != 2); 55 if (fileType) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File type %d is not a valid Gmsh ASCII file", fileType); 56 if (dataSize != sizeof(double)) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Data size %d is not valid for a Gmsh file", dataSize); 57 fgets(line, PETSC_MAX_PATH_LEN, fd); 58 ierr = PetscStrncmp(line, "$EndMeshFormat\n", PETSC_MAX_PATH_LEN, &match);CHKERRQ(ierr); 59 if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file"); 60 /* Read vertices */ 61 fgets(line, PETSC_MAX_PATH_LEN, fd); 62 ierr = PetscStrncmp(line, "$Nodes\n", PETSC_MAX_PATH_LEN, &match);CHKERRQ(ierr); 63 if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file"); 64 snum = fscanf(fd, "%d\n", &numVertices);CHKERRQ(snum != 1); 65 ierr = PetscMalloc(numVertices*3 * sizeof(PetscScalar), &coordsIn);CHKERRQ(ierr); 66 for (v = 0; v < numVertices; ++v) { 67 double x, y, z; 68 int i; 69 70 snum = fscanf(fd, "%d %lg %lg %lg\n", &i, &x, &y, &z);CHKERRQ(snum != 4); 71 coordsIn[v*3+0] = x; coordsIn[v*3+1] = y; coordsIn[v*3+2] = z; 72 if (i != v+1) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid node number %d should be %d", i, v+1); 73 } 74 fgets(line, PETSC_MAX_PATH_LEN, fd); 75 ierr = PetscStrncmp(line, "$EndNodes\n", PETSC_MAX_PATH_LEN, &match);CHKERRQ(ierr); 76 if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file"); 77 /* Read cells */ 78 fgets(line, PETSC_MAX_PATH_LEN, fd); 79 ierr = PetscStrncmp(line, "$Elements\n", PETSC_MAX_PATH_LEN, &match);CHKERRQ(ierr); 80 if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file"); 81 snum = fscanf(fd, "%d\n", &numCells);CHKERRQ(snum != 1); 82 } 83 84 if (!rank) { 85 /* Gmsh elements can be of any dimension/co-dimension, so we need to traverse the 86 file contents multiple times to figure out the true number of cells and facets 87 in the given mesh. To make this more efficient we read the file contents only 88 once and store them in memory, while determining the true number of cells. */ 89 ierr = PetscMalloc1(numCells, &gmsh_elem);CHKERRQ(ierr); 90 for (trueNumCells=0, c = 0; c < numCells; ++c) { 91 ierr = DMPlexCreateGmsh_ReadElement(fd, &gmsh_elem[c]);CHKERRQ(ierr); 92 if (gmsh_elem[c].dim > dim) {dim = gmsh_elem[c].dim; trueNumCells = 0;} 93 if (gmsh_elem[c].dim == dim) trueNumCells++; 94 } 95 fgets(line, PETSC_MAX_PATH_LEN, fd); 96 ierr = PetscStrncmp(line, "$EndElements\n", PETSC_MAX_PATH_LEN, &match);CHKERRQ(ierr); 97 if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file"); 98 } 99 /* Allocate the cell-vertex mesh */ 100 ierr = DMPlexSetChart(*dm, 0, trueNumCells+numVertices);CHKERRQ(ierr); 101 if (!rank) { 102 for (cell = 0, c = 0; c < numCells; ++c) { 103 if (gmsh_elem[c].dim == dim) { 104 ierr = DMPlexSetConeSize(*dm, cell, gmsh_elem[c].numNodes);CHKERRQ(ierr); 105 cell++; 106 } 107 } 108 } 109 ierr = DMSetUp(*dm);CHKERRQ(ierr); 110 /* Add cell-vertex connections */ 111 if (!rank) { 112 PetscInt pcone[8], corner; 113 for (cell = 0, c = 0; c < numCells; ++c) { 114 if (gmsh_elem[c].dim == dim) { 115 for (corner = 0; corner < gmsh_elem[c].numNodes; ++corner) { 116 pcone[corner] = gmsh_elem[c].nodes[corner] + trueNumCells-1; 117 } 118 ierr = DMPlexSetCone(*dm, cell, pcone);CHKERRQ(ierr); 119 cell++; 120 } 121 } 122 } 123 ierr = MPI_Bcast(&dim, 1, MPIU_INT, 0, comm);CHKERRQ(ierr); 124 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 125 ierr = DMPlexSymmetrize(*dm);CHKERRQ(ierr); 126 ierr = DMPlexStratify(*dm);CHKERRQ(ierr); 127 if (interpolate) { 128 DM idm = NULL; 129 130 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 131 ierr = DMDestroy(dm);CHKERRQ(ierr); 132 *dm = idm; 133 } 134 135 if (!rank) { 136 /* Apply boundary IDs by finding the relevant facets with vertex joins */ 137 PetscInt pcone[8], corner, vStart, vEnd; 138 139 ierr = DMPlexGetDepthStratum(*dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 140 for (c = 0; c < numCells; ++c) { 141 if (gmsh_elem[c].dim == dim-1) { 142 PetscInt joinSize; 143 const PetscInt *join; 144 for (corner = 0; corner < gmsh_elem[c].numNodes; ++corner) { 145 pcone[corner] = gmsh_elem[c].nodes[corner] + vStart - 1; 146 } 147 ierr = DMPlexGetFullJoin(*dm, gmsh_elem[c].numNodes, (const PetscInt *) pcone, &joinSize, &join);CHKERRQ(ierr); 148 if (joinSize != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not determine Plex facet for element %d", gmsh_elem[c].id); 149 ierr = DMPlexSetLabelValue(*dm, "Face Sets", join[0], gmsh_elem[c].tags[0]);CHKERRQ(ierr); 150 ierr = DMPlexRestoreJoin(*dm, gmsh_elem[c].numNodes, (const PetscInt *) pcone, &joinSize, &join);CHKERRQ(ierr); 151 } 152 } 153 } 154 155 /* Read coordinates */ 156 ierr = DMGetCoordinateSection(*dm, &coordSection);CHKERRQ(ierr); 157 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 158 ierr = PetscSectionSetFieldComponents(coordSection, 0, dim);CHKERRQ(ierr); 159 ierr = PetscSectionSetChart(coordSection, trueNumCells, trueNumCells + numVertices);CHKERRQ(ierr); 160 for (v = trueNumCells; v < trueNumCells+numVertices; ++v) { 161 ierr = PetscSectionSetDof(coordSection, v, dim);CHKERRQ(ierr); 162 ierr = PetscSectionSetFieldDof(coordSection, v, 0, dim);CHKERRQ(ierr); 163 } 164 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 165 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 166 ierr = VecCreate(comm, &coordinates);CHKERRQ(ierr); 167 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 168 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 169 ierr = VecSetType(coordinates, VECSTANDARD);CHKERRQ(ierr); 170 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 171 if (!rank) { 172 for (v = 0; v < numVertices; ++v) { 173 for (d = 0; d < dim; ++d) { 174 coords[v*dim+d] = coordsIn[v*3+d]; 175 } 176 } 177 } 178 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 179 ierr = PetscFree(coordsIn);CHKERRQ(ierr); 180 ierr = DMSetCoordinatesLocal(*dm, coordinates);CHKERRQ(ierr); 181 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 182 /* Clean up intermediate storage */ 183 if (!rank) { 184 for (c = 0; c < numCells; ++c) { 185 ierr = PetscFree(gmsh_elem[c].nodes); 186 ierr = PetscFree(gmsh_elem[c].tags); 187 } 188 ierr = PetscFree(gmsh_elem); 189 } 190 PetscFunctionReturn(0); 191 } 192 193 #undef __FUNCT__ 194 #define __FUNCT__ "DMPlexCreateGmsh_ReadElement" 195 PetscErrorCode DMPlexCreateGmsh_ReadElement(FILE *fd, GmshElement *ele) 196 { 197 int snum, cellType; 198 PetscInt t; 199 PetscErrorCode ierr; 200 201 PetscFunctionBegin; 202 snum = fscanf(fd, "%d %d %d", &(ele->id), &cellType, &(ele->numTags));CHKERRQ(snum != 3); 203 ierr = PetscMalloc1(ele->numTags, &(ele->tags));CHKERRQ(ierr); 204 for (t=0; t<ele->numTags; t++) {snum = fscanf(fd, "%d", &(ele->tags[t]));CHKERRQ(snum != 1);} 205 switch (cellType) { 206 case 1: /* 2-node line */ 207 ele->dim = 1; 208 ele->numNodes = 2; 209 ierr = PetscMalloc1(ele->numNodes, &(ele->nodes));CHKERRQ(ierr); 210 snum = fscanf(fd, "%d %d\n", &(ele->nodes[0]), &(ele->nodes[1]));CHKERRQ(snum != ele->numNodes); 211 break; 212 case 2: /* 3-node triangle */ 213 ele->dim = 2; 214 ele->numNodes = 3; 215 ierr = PetscMalloc1(ele->numNodes, &(ele->nodes));CHKERRQ(ierr); 216 snum = fscanf(fd, "%d %d %d\n", &(ele->nodes[0]), &(ele->nodes[1]), &(ele->nodes[2]));CHKERRQ(snum != ele->numNodes); 217 break; 218 case 3: /* 4-node quadrangle */ 219 ele->dim = 2; 220 ele->numNodes = 4; 221 ierr = PetscMalloc1(ele->numNodes, &(ele->nodes));CHKERRQ(ierr); 222 snum = fscanf(fd, "%d %d %d %d\n", &(ele->nodes[0]), &(ele->nodes[1]), &(ele->nodes[2]), &(ele->nodes[3]));CHKERRQ(snum != ele->numNodes); 223 break; 224 case 4: /* 4-node tetrahedron */ 225 ele->dim = 3; 226 ele->numNodes = 4; 227 ierr = PetscMalloc1(ele->numNodes, &(ele->nodes));CHKERRQ(ierr); 228 snum = fscanf(fd, "%d %d %d %d\n", &(ele->nodes[0]), &(ele->nodes[1]), &(ele->nodes[2]), &(ele->nodes[3]));CHKERRQ(snum != ele->numNodes); 229 break; 230 case 5: /* 8-node hexahedron */ 231 ele->dim = 3; 232 ele->numNodes = 8; 233 ierr = PetscMalloc1(ele->numNodes, &(ele->nodes));CHKERRQ(ierr); 234 snum = fscanf(fd, "%d %d %d %d %d %d %d %d\n", &(ele->nodes[0]), &(ele->nodes[1]), &(ele->nodes[2]), &(ele->nodes[3]), &(ele->nodes[4]), &(ele->nodes[5]), &(ele->nodes[6]), &(ele->nodes[7]));CHKERRQ(snum != ele->numNodes); 235 break; 236 case 15: /* 1-node vertex */ 237 ele->dim = 0; 238 ele->numNodes = 1; 239 ierr = PetscMalloc1(ele->numNodes, &(ele->nodes));CHKERRQ(ierr); 240 snum = fscanf(fd, "%d\n", &(ele->nodes[0]));CHKERRQ(snum != ele->numNodes); 241 break; 242 default: 243 SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unsupported Gmsh element type %d", cellType); 244 } 245 PetscFunctionReturn(0); 246 } 247