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 PetscSection coordSection; 30 Vec coordinates; 31 PetscScalar *coords, *coordsIn = NULL; 32 PetscInt dim = 0, tdim = 0, coordSize, c, v, d; 33 int numVertices = 0, numCells = 0, trueNumCells = 0, numFacets = 0, snum; 34 long fpos = 0; 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 fpos = ftell(fd); 86 /* The Gmsh format disguises facets as elements, so we have to run through all "element" entries 87 to get the correct numCells and decide the topological dimension of the mesh */ 88 trueNumCells = 0; 89 for (c = 0; c < numCells; ++c) { 90 PetscInt numCorners, t; 91 int cone[8], i, cellType, numTags, tag; 92 snum = fscanf(fd, "%d %d %d", &i, &cellType, &numTags);CHKERRQ(snum != 3); 93 if (numTags) for (t = 0; t < numTags; ++t) {snum = fscanf(fd, "%d", &tag);CHKERRQ(snum != 1);} 94 switch (cellType) { 95 case 1: /* 2-node line */ 96 dim = 1; 97 numCorners = 2; 98 snum = fscanf(fd, "%d %d\n", &cone[0], &cone[1]);CHKERRQ(snum != numCorners); 99 break; 100 case 2: /* 3-node triangle */ 101 dim = 2; 102 numCorners = 3; 103 snum = fscanf(fd, "%d %d %d\n", &cone[0], &cone[1], &cone[2]);CHKERRQ(snum != numCorners); 104 break; 105 case 3: /* 4-node quadrangle */ 106 dim = 2; 107 numCorners = 4; 108 snum = fscanf(fd, "%d %d %d %d\n", &cone[0], &cone[1], &cone[2], &cone[3]);CHKERRQ(snum != numCorners); 109 break; 110 case 4: /* 4-node tetrahedron */ 111 dim = 3; 112 numCorners = 4; 113 snum = fscanf(fd, "%d %d %d %d\n", &cone[0], &cone[1], &cone[2], &cone[3]);CHKERRQ(snum != numCorners); 114 break; 115 case 5: /* 8-node hexahedron */ 116 dim = 3; 117 numCorners = 8; 118 snum = fscanf(fd, "%d %d %d %d %d %d %d %d\n", &cone[0], &cone[1], &cone[2], &cone[3], &cone[4], &cone[5], &cone[6], &cone[7]);CHKERRQ(snum != numCorners); 119 break; 120 default: 121 SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unsupported Gmsh element type %d", cellType); 122 } 123 if (dim > tdim) { 124 tdim = dim; 125 trueNumCells = 0; 126 } 127 trueNumCells++; 128 } 129 } 130 ierr = DMPlexSetChart(*dm, 0, trueNumCells+numVertices);CHKERRQ(ierr); 131 numFacets = numCells - trueNumCells; 132 if (!rank) { 133 ierr = fseek(fd, fpos, SEEK_SET);CHKERRQ(ierr); 134 for (c = 0; c < numCells; ++c) { 135 PetscInt numCorners, t; 136 int cone[8], i, cellType, numTags, tag; 137 138 snum = fscanf(fd, "%d %d %d", &i, &cellType, &numTags);CHKERRQ(snum != 3); 139 if (numTags) for (t = 0; t < numTags; ++t) {snum = fscanf(fd, "%d", &tag);CHKERRQ(snum != 1);} 140 switch (cellType) { 141 case 1: /* 2-node line */ 142 dim = 1; 143 numCorners = 2; 144 snum = fscanf(fd, "%d %d\n", &cone[0], &cone[1]);CHKERRQ(snum != numCorners); 145 break; 146 case 2: /* 3-node triangle */ 147 dim = 2; 148 numCorners = 3; 149 snum = fscanf(fd, "%d %d %d\n", &cone[0], &cone[1], &cone[2]);CHKERRQ(snum != numCorners); 150 break; 151 case 3: /* 4-node quadrangle */ 152 dim = 2; 153 numCorners = 4; 154 snum = fscanf(fd, "%d %d %d %d\n", &cone[0], &cone[1], &cone[2], &cone[3]);CHKERRQ(snum != numCorners); 155 break; 156 case 4: /* 4-node tetrahedron */ 157 dim = 3; 158 numCorners = 4; 159 snum = fscanf(fd, "%d %d %d %d\n", &cone[0], &cone[1], &cone[2], &cone[3]);CHKERRQ(snum != numCorners); 160 break; 161 case 5: /* 8-node hexahedron */ 162 dim = 3; 163 numCorners = 8; 164 snum = fscanf(fd, "%d %d %d %d %d %d %d %d\n", &cone[0], &cone[1], &cone[2], &cone[3], &cone[4], &cone[5], &cone[6], &cone[7]);CHKERRQ(snum != numCorners); 165 break; 166 default: 167 SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unsupported Gmsh element type %d", cellType); 168 } 169 if (dim == tdim) { 170 ierr = DMPlexSetConeSize(*dm, c-numFacets, numCorners);CHKERRQ(ierr); 171 } 172 if (i != c+1) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid cell number %d should be %d", i, c+1); 173 } 174 } 175 ierr = DMSetUp(*dm);CHKERRQ(ierr); 176 if (!rank) { 177 ierr = fseek(fd, fpos, SEEK_SET);CHKERRQ(ierr); 178 for (c = 0; c < numCells; ++c) { 179 PetscInt pcone[8], numCorners, corner, t; 180 int cone[8], i, cellType, numTags, tag; 181 182 snum = fscanf(fd, "%d %d %d", &i, &cellType, &numTags);CHKERRQ(snum != 3); 183 if (numTags) for (t = 0; t < numTags; ++t) {snum = fscanf(fd, "%d", &tag);CHKERRQ(snum != 1);} 184 switch (cellType) { 185 case 1: /* 2-node line */ 186 dim = 1; 187 numCorners = 2; 188 snum = fscanf(fd, "%d %d\n", &cone[0], &cone[1]);CHKERRQ(snum != numCorners); 189 break; 190 case 2: /* 3-node triangle */ 191 dim = 2; 192 numCorners = 3; 193 snum = fscanf(fd, "%d %d %d\n", &cone[0], &cone[1], &cone[2]);CHKERRQ(snum != numCorners); 194 break; 195 case 3: /* 4-node quadrangle */ 196 dim = 2; 197 numCorners = 4; 198 snum = fscanf(fd, "%d %d %d %d\n", &cone[0], &cone[1], &cone[2], &cone[3]);CHKERRQ(snum != numCorners); 199 break; 200 case 4: /* 4-node tetrahedron */ 201 dim = 3; 202 numCorners = 4; 203 snum = fscanf(fd, "%d %d %d %d\n", &cone[0], &cone[1], &cone[2], &cone[3]);CHKERRQ(snum != numCorners); 204 ierr = DMPlexInvertCell(dim, numCorners, cone);CHKERRQ(ierr); 205 break; 206 case 5: /* 8-node hexahedron */ 207 dim = 3; 208 numCorners = 8; 209 snum = fscanf(fd, "%d %d %d %d %d %d %d %d\n", &cone[0], &cone[1], &cone[2], &cone[3], &cone[4], &cone[5], &cone[6], &cone[7]);CHKERRQ(snum != numCorners); 210 ierr = DMPlexInvertCell(dim, numCorners, cone);CHKERRQ(ierr); 211 break; 212 default: 213 SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unsupported Gmsh element type %d", cellType); 214 } 215 if (dim == tdim) { 216 for (corner = 0; corner < numCorners; ++corner) pcone[corner] = cone[corner] + trueNumCells-1; 217 ierr = DMPlexSetCone(*dm, c-numFacets, pcone);CHKERRQ(ierr); 218 } 219 if (i != c+1) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid cell number %d should be %d", i, c+1); 220 } 221 fgets(line, PETSC_MAX_PATH_LEN, fd); 222 ierr = PetscStrncmp(line, "$EndElements\n", PETSC_MAX_PATH_LEN, &match);CHKERRQ(ierr); 223 if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file"); 224 } 225 ierr = MPI_Bcast(&dim, 1, MPIU_INT, 0, comm);CHKERRQ(ierr); 226 ierr = DMPlexSetDimension(*dm, dim);CHKERRQ(ierr); 227 ierr = DMPlexSymmetrize(*dm);CHKERRQ(ierr); 228 ierr = DMPlexStratify(*dm);CHKERRQ(ierr); 229 if (interpolate) { 230 DM idm; 231 232 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 233 ierr = DMDestroy(dm);CHKERRQ(ierr); 234 *dm = idm; 235 } 236 /* Read coordinates */ 237 ierr = DMGetCoordinateSection(*dm, &coordSection);CHKERRQ(ierr); 238 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 239 ierr = PetscSectionSetFieldComponents(coordSection, 0, dim);CHKERRQ(ierr); 240 ierr = PetscSectionSetChart(coordSection, numCells, numCells + numVertices);CHKERRQ(ierr); 241 for (v = numCells; v < numCells+numVertices; ++v) { 242 ierr = PetscSectionSetDof(coordSection, v, dim);CHKERRQ(ierr); 243 ierr = PetscSectionSetFieldDof(coordSection, v, 0, dim);CHKERRQ(ierr); 244 } 245 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 246 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 247 ierr = VecCreate(comm, &coordinates);CHKERRQ(ierr); 248 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 249 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 250 ierr = VecSetType(coordinates, VECSTANDARD);CHKERRQ(ierr); 251 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 252 if (!rank) { 253 for (v = 0; v < numVertices; ++v) { 254 for (d = 0; d < dim; ++d) { 255 coords[v*dim+d] = coordsIn[v*3+d]; 256 } 257 } 258 } 259 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 260 ierr = PetscFree(coordsIn);CHKERRQ(ierr); 261 ierr = DMSetCoordinatesLocal(*dm, coordinates);CHKERRQ(ierr); 262 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 263 PetscFunctionReturn(0); 264 } 265 266 #undef __FUNCT__ 267 #define __FUNCT__ "DMPlexView_Ascii" 268 PetscErrorCode DMPlexView_Ascii(DM dm, PetscViewer viewer) 269 { 270 271 } 272