1 #define PETSCDM_DLL 2 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/ 3 4 /*@C 5 DMPlexCreatePLYFromFile - Create a DMPlex mesh from a PLY file. 6 7 + comm - The MPI communicator 8 . filename - Name of the .med file 9 - interpolate - Create faces and edges in the mesh 10 11 Output Parameter: 12 . dm - The DM object representing the mesh 13 14 Note: https://en.wikipedia.org/wiki/PLY_(file_format) 15 16 Level: beginner 17 18 .seealso: DMPlexCreateFromFile(), DMPlexCreateMedFromFile(), DMPlexCreateGmsh(), DMPlexCreate() 19 @*/ 20 PetscErrorCode DMPlexCreatePLYFromFile(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm) 21 { 22 PetscViewer viewer; 23 Vec coordinates; 24 PetscSection coordSection; 25 PetscScalar *coords; 26 char line[PETSC_MAX_PATH_LEN], ntype[16], itype[16], name[1024], vtype[16]; 27 PetscBool match, byteSwap = PETSC_FALSE; 28 PetscInt dim = 2, cdim = 3, Nvp = 0, coordSize, xi = -1, yi = -1, zi = -1, v, c, p; 29 PetscMPIInt rank; 30 int snum, Nv, Nc; 31 PetscErrorCode ierr; 32 33 PetscFunctionBegin; 34 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 35 ierr = PetscViewerCreate(comm, &viewer);CHKERRQ(ierr); 36 ierr = PetscViewerSetType(viewer, PETSCVIEWERBINARY);CHKERRQ(ierr); 37 ierr = PetscViewerFileSetMode(viewer, FILE_MODE_READ);CHKERRQ(ierr); 38 ierr = PetscViewerFileSetName(viewer, filename);CHKERRQ(ierr); 39 if (!rank) { 40 PetscBool isAscii, isBinaryBig, isBinaryLittle; 41 42 /* Check for PLY file */ 43 ierr = PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);CHKERRQ(ierr); 44 ierr = PetscStrncmp(line, "ply", PETSC_MAX_PATH_LEN, &match);CHKERRQ(ierr); 45 if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid PLY file"); 46 /* Check PLY format */ 47 ierr = PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);CHKERRQ(ierr); 48 ierr = PetscStrncmp(line, "format", PETSC_MAX_PATH_LEN, &match);CHKERRQ(ierr); 49 if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid PLY file"); 50 ierr = PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);CHKERRQ(ierr); 51 ierr = PetscStrncmp(line, "ascii", PETSC_MAX_PATH_LEN, &isAscii);CHKERRQ(ierr); 52 ierr = PetscStrncmp(line, "binary_big_endian", PETSC_MAX_PATH_LEN, &isBinaryBig);CHKERRQ(ierr); 53 ierr = PetscStrncmp(line, "binary_little_endian", PETSC_MAX_PATH_LEN, &isBinaryLittle);CHKERRQ(ierr); 54 if (isAscii) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "PLY ascii format not yet supported"); 55 else if (isBinaryLittle) byteSwap = PETSC_TRUE; 56 ierr = PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);CHKERRQ(ierr); 57 ierr = PetscStrncmp(line, "1.0", PETSC_MAX_PATH_LEN, &match);CHKERRQ(ierr); 58 if (!match) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid version of PLY file, %s", line); 59 /* Ignore comments */ 60 match = PETSC_TRUE; 61 while (match) { 62 char buf = '\0'; 63 ierr = PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);CHKERRQ(ierr); 64 ierr = PetscStrncmp(line, "comment", PETSC_MAX_PATH_LEN, &match);CHKERRQ(ierr); 65 if (match) while (buf != '\n') {ierr = PetscViewerRead(viewer, &buf, 1, NULL, PETSC_CHAR);CHKERRQ(ierr);} 66 } 67 /* Read vertex information */ 68 ierr = PetscStrncmp(line, "element", PETSC_MAX_PATH_LEN, &match);CHKERRQ(ierr); 69 if (!match) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse PLY file header: %s", line); 70 ierr = PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);CHKERRQ(ierr); 71 ierr = PetscStrncmp(line, "vertex", PETSC_MAX_PATH_LEN, &match);CHKERRQ(ierr); 72 if (!match) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse PLY file header: %s", line); 73 ierr = PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);CHKERRQ(ierr); 74 snum = sscanf(line, "%d", &Nv); 75 if (snum != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse PLY file header: %s", line); 76 match = PETSC_TRUE; 77 while (match) { 78 ierr = PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);CHKERRQ(ierr); 79 ierr = PetscStrncmp(line, "property", PETSC_MAX_PATH_LEN, &match);CHKERRQ(ierr); 80 if (match) { 81 PetscBool matchB; 82 83 ierr = PetscViewerRead(viewer, line, 2, NULL, PETSC_STRING);CHKERRQ(ierr); 84 if (Nvp >= 16) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cannot handle more than 16 property statements in PLY file header: %s", line); 85 snum = sscanf(line, "%s %s", ntype, name); 86 if (snum != 2) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse PLY file header: %s", line); 87 ierr = PetscStrncmp(ntype, "float32", 16, &matchB);CHKERRQ(ierr); 88 if (matchB) { 89 vtype[Nvp] = 'f'; 90 } else { 91 ierr = PetscStrncmp(ntype, "int32", 16, &matchB);CHKERRQ(ierr); 92 if (matchB) { 93 vtype[Nvp] = 'd'; 94 } else { 95 ierr = PetscStrncmp(ntype, "uint8", 16, &matchB);CHKERRQ(ierr); 96 if (matchB) { 97 vtype[Nvp] = 'c'; 98 } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse type in PLY file header: %s", line); 99 } 100 } 101 ierr = PetscStrncmp(name, "x", 16, &matchB);CHKERRQ(ierr); 102 if (matchB) {xi = Nvp;} 103 ierr = PetscStrncmp(name, "y", 16, &matchB);CHKERRQ(ierr); 104 if (matchB) {yi = Nvp;} 105 ierr = PetscStrncmp(name, "z", 16, &matchB);CHKERRQ(ierr); 106 if (matchB) {zi = Nvp;} 107 ++Nvp; 108 } 109 } 110 /* Read cell information */ 111 ierr = PetscStrncmp(line, "element", PETSC_MAX_PATH_LEN, &match);CHKERRQ(ierr); 112 if (!match) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse PLY file header: %s", line); 113 ierr = PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);CHKERRQ(ierr); 114 ierr = PetscStrncmp(line, "face", PETSC_MAX_PATH_LEN, &match);CHKERRQ(ierr); 115 if (!match) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse PLY file header: %s", line); 116 ierr = PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);CHKERRQ(ierr); 117 snum = sscanf(line, "%d", &Nc); 118 if (snum != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse PLY file header: %s", line); 119 ierr = PetscViewerRead(viewer, line, 5, NULL, PETSC_STRING);CHKERRQ(ierr); 120 snum = sscanf(line, "property list %s %s %s", ntype, itype, name); 121 if (snum != 3) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse PLY file header: %s", line); 122 ierr = PetscStrncmp(ntype, "uint8", 1024, &match);CHKERRQ(ierr); 123 if (!match) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid size type in PLY file header: %s", line); 124 ierr = PetscStrncmp(name, "vertex_indices", 1024, &match);CHKERRQ(ierr); 125 if (!match) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid property in PLY file header: %s", line); 126 /* Header should terminate */ 127 ierr = PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);CHKERRQ(ierr); 128 ierr = PetscStrncmp(line, "end_header", PETSC_MAX_PATH_LEN, &match);CHKERRQ(ierr); 129 if (!match) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid PLY file: %s", line); 130 } else { 131 Nc = Nv = 0; 132 } 133 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 134 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 135 ierr = DMPlexSetChart(*dm, 0, Nc+Nv);CHKERRQ(ierr); 136 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 137 ierr = DMSetCoordinateDim(*dm, cdim);CHKERRQ(ierr); 138 /* Read coordinates */ 139 ierr = DMGetCoordinateSection(*dm, &coordSection);CHKERRQ(ierr); 140 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 141 ierr = PetscSectionSetFieldComponents(coordSection, 0, cdim);CHKERRQ(ierr); 142 ierr = PetscSectionSetChart(coordSection, Nc, Nc + Nv);CHKERRQ(ierr); 143 for (v = Nc; v < Nc+Nv; ++v) { 144 ierr = PetscSectionSetDof(coordSection, v, cdim);CHKERRQ(ierr); 145 ierr = PetscSectionSetFieldDof(coordSection, v, 0, cdim);CHKERRQ(ierr); 146 } 147 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 148 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 149 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr); 150 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 151 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 152 ierr = VecSetBlockSize(coordinates, cdim);CHKERRQ(ierr); 153 ierr = VecSetType(coordinates, VECSTANDARD);CHKERRQ(ierr); 154 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 155 if (!rank) { 156 float rbuf[1]; 157 int ibuf[1]; 158 159 for (v = 0; v < Nv; ++v) { 160 for (p = 0; p < Nvp; ++p) { 161 if (vtype[p] == 'f') { 162 ierr = PetscViewerRead(viewer, &rbuf, 1, NULL, PETSC_FLOAT);CHKERRQ(ierr); 163 if (byteSwap) ierr = PetscByteSwap(&rbuf, PETSC_FLOAT, 1);CHKERRQ(ierr); 164 if (p == xi) coords[v*cdim+0] = rbuf[0]; 165 else if (p == yi) coords[v*cdim+1] = rbuf[0]; 166 else if (p == zi) coords[v*cdim+2] = rbuf[0]; 167 } else if (vtype[p] == 'd') { 168 ierr = PetscViewerRead(viewer, &ibuf, 1, NULL, PETSC_INT);CHKERRQ(ierr); 169 if (byteSwap) ierr = PetscByteSwap(&ibuf, PETSC_INT, 1);CHKERRQ(ierr); 170 } else if (vtype[p] == 'c') { 171 ierr = PetscViewerRead(viewer, &ibuf, 1, NULL, PETSC_CHAR);CHKERRQ(ierr); 172 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid vertex property type in PLY file"); 173 } 174 } 175 } 176 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 177 ierr = DMSetCoordinatesLocal(*dm, coordinates);CHKERRQ(ierr); 178 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 179 /* Read topology */ 180 if (!rank) { 181 char ibuf[1]; 182 PetscInt vbuf[16], corners; 183 184 /* Assume same cells */ 185 ierr = PetscViewerRead(viewer, &ibuf, 1, NULL, PETSC_CHAR);CHKERRQ(ierr); 186 corners = ibuf[0]; 187 for (c = 0; c < Nc; ++c) {ierr = DMPlexSetConeSize(*dm, c, corners);CHKERRQ(ierr);} 188 ierr = DMSetUp(*dm);CHKERRQ(ierr); 189 for (c = 0; c < Nc; ++c) { 190 if (c > 0) { 191 ierr = PetscViewerRead(viewer, &ibuf, 1, NULL, PETSC_CHAR);CHKERRQ(ierr); 192 } 193 if (ibuf[0] != corners) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "All cells must have the same number of vertices in PLY file: %D != %D", ibuf[0], corners); 194 ierr = PetscViewerRead(viewer, &vbuf, ibuf[0], NULL, PETSC_INT);CHKERRQ(ierr); 195 if (byteSwap) ierr = PetscByteSwap(&vbuf, PETSC_INT, ibuf[0]);CHKERRQ(ierr); 196 for (v = 0; v < ibuf[0]; ++v) vbuf[v] += Nc; 197 ierr = DMPlexSetCone(*dm, c, vbuf);CHKERRQ(ierr); 198 } 199 } 200 ierr = DMPlexSymmetrize(*dm);CHKERRQ(ierr); 201 ierr = DMPlexStratify(*dm);CHKERRQ(ierr); 202 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 203 if (interpolate) { 204 DM idm; 205 206 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 207 ierr = DMDestroy(dm);CHKERRQ(ierr); 208 *dm = idm; 209 } 210 PetscFunctionReturn(0); 211 } 212