Actual source code: meshlagrit.c
petsc-3.4.5 2014-06-29
1: #include <petscdmmesh_formats.hh> /*I "petscmesh.h" I*/
3: #if 0
5: void FlipCellOrientation(pylith::int_array *const cells, const int numCells, const int numCorners, const int meshDim)
6: {
7: if (3 == meshDim && 4 == numCorners) {
8: for (int iCell = 0; iCell < numCells; ++iCell) {
9: const int i1 = iCell*numCorners+1;
10: const int i2 = iCell*numCorners+2;
11: const int tmp = (*cells)[i1];
12: (*cells)[i1] = (*cells)[i2];
13: (*cells)[i2] = tmp;
14: }
15: }
16: }
18: //Obj<PETSC_MESH_TYPE> m = ALE::LaGriT::Builder::readMesh(comm, 3, options->baseFilename, options->interpolate, options->debug);'
19: Obj<PETSC_MESH_TYPE> m = new PETSC_MESH_TYPE(comm, options->dim, options->debug);
20: Obj<PETSC_MESH_TYPE::sieve_type> sieve = new PETSC_MESH_TYPE::sieve_type(comm, options->debug);
21: bool flipEndian = false;
22: int dim;
23: pylith::int_array cells;
24: pylith::double_array coordinates;
25: pylith::int_array materialIds;
26: int numCells = 0, numVertices = 0, numCorners = 0;
28: if (!m->commRank()) {
29: if (pylith::meshio::GMVFile::isAscii(options->baseFilename)) {
30: pylith::meshio::GMVFileAscii filein(options->baseFilename);
31: filein.read(&coordinates, &cells, &materialIds, &dim, &dim, &numVertices, &numCells, &numCorners);
32: if (options->interpolate) {
33: FlipCellOrientation(&cells, numCells, numCorners, dim);
34: }
35: } else {
36: pylith::meshio::GMVFileBinary filein(options->baseFilename, flipEndian);
37: filein.read(&coordinates, &cells, &materialIds, &dim, &dim, &numVertices, &numCells, &numCorners);
38: if (!options->interpolate) {
39: FlipCellOrientation(&cells, numCells, numCorners, dim);
40: }
41: } // if/else
42: }
43: ALE::SieveBuilder<PETSC_MESH_TYPE>::buildTopology(sieve, dim, numCells, const_cast<int*>(&cells[0]), numVertices, options->interpolate, numCorners, -1, m->getArrowSection("orientation"));
44: m->setSieve(sieve);
45: m->stratify();
46: ALE::SieveBuilder<PETSC_MESH_TYPE>::buildCoordinates(m, dim, const_cast<double*>(&coordinates[0]));
48: DMMeshCreate(comm, &mesh);
49: DMMeshSetMesh(mesh, m);
50: DMMeshIDBoundary(mesh);
51: #endif
53: namespace ALE {
54: namespace LaGriT {
55: void Builder::readInpFile(MPI_Comm comm, const std::string& filename, const int dim, const int numCorners, int& numElements, int *vertices[], int& numVertices, PetscReal *coordinates[])
56: {
57: PetscViewer viewer;
58: FILE *f;
59: PetscReal *coords;
60: PetscInt *verts;
61: PetscInt commRank;
62: char buf[2048];
65: MPI_Comm_rank(comm, &commRank);CHKERRXX(ierr);
66: if (commRank != 0) return;
67: PetscViewerCreate(PETSC_COMM_SELF, &viewer);CHKERRXX(ierr);
68: PetscViewerSetType(viewer, PETSCVIEWERASCII);CHKERRXX(ierr);
69: PetscViewerFileSetMode(viewer, FILE_MODE_READ);CHKERRXX(ierr);
70: PetscViewerFileSetName(viewer, filename.c_str());CHKERRXX(ierr);
71: PetscViewerASCIIGetPointer(viewer, &f);CHKERRXX(ierr);
72: // Read header
73: if (!fgets(buf, 2048, f)) throw ALE::Exception("File ended prematurely");
74: // Number of vertices
75: const char *x = strtok(buf, " ");
76: numVertices = atoi(x);
77: // Number of elements
78: x = strtok(NULL, " "); numElements = atoi(x);
79: // Element type
80: x = strtok(NULL, " ");
81: // ???
82: x = strtok(NULL, " ");
83: // ???
84: x = strtok(NULL, " ");
85: PetscMalloc(numVertices*dim * sizeof(PetscReal), &coords);CHKERRXX(ierr);
86: for (PetscInt i = 0; i < numVertices; ++i) {
87: if (!fgets(buf, 2048, f)) throw ALE::Exception("File ended prematurely");
88: x = strtok(buf, " ");
89: // Ignore vertex number
90: x = strtok(NULL, " ");
91: for (int c = 0; c < dim; c++) {
92: coords[i*dim+c] = atof(x);
94: x = strtok(NULL, " ");
95: }
96: }
97: *coordinates = coords;
98: PetscMalloc(numElements*numCorners * sizeof(PetscInt), &verts);CHKERRXX(ierr);
99: for (PetscInt i = 0; i < numElements; ++i) {
100: if (!fgets(buf, 2048, f)) throw ALE::Exception("File ended prematurely");
101: x = strtok(buf, " ");
102: // Ignore element number
103: x = strtok(NULL, " ");
104: // Ignore ??? (material)
105: x = strtok(NULL, " ");
106: // Ignore element type
107: x = strtok(NULL, " ");
108: for (int c = 0; c < numCorners; c++) {
109: verts[i*numCorners+c] = atoi(x) - 1;
111: x = strtok(NULL, " ");
112: }
113: }
114: *vertices = verts;
115: PetscViewerDestroy(&viewer);CHKERRXX(ierr);
116: };
117: Obj<Builder::Mesh> Builder::readMesh(MPI_Comm comm, const int dim, const std::string& filename, const bool interpolate = false, const int debug = 0)
118: {
119: typedef ALE::Mesh<PetscInt,PetscScalar> FlexMesh;
120: Obj<Mesh> mesh = new Mesh(comm, dim, debug);
121: Obj<sieve_type> sieve = new sieve_type(comm, debug);
122: Obj<FlexMesh> m = new FlexMesh(comm, dim, debug);
123: Obj<FlexMesh::sieve_type> s = new FlexMesh::sieve_type(comm, debug);
124: std::map<Mesh::point_type,Mesh::point_type> renumbering;
125: int *cells;
126: PetscReal *coordinates;
127: int numCells = 0, numVertices = 0, numCorners = dim+1;
128: PetscErrorCode ierr;
130: mesh->setSieve(sieve);
131: Builder::readInpFile(comm, filename, dim, numCorners, numCells, &cells, numVertices, &coordinates);
132: ALE::SieveBuilder<FlexMesh>::buildTopology(s, dim, numCells, cells, numVertices, interpolate, numCorners, -1, m->getArrowSection("orientation"));
133: m->setSieve(s);
134: m->stratify();
135: ALE::SieveBuilder<FlexMesh>::buildCoordinates(m, dim+1, coordinates);
136: PetscFree(cells);CHKERRXX(ierr);
137: PetscFree(coordinates);CHKERRXX(ierr);
138: ALE::ISieveConverter::convertMesh(*m, *mesh, renumbering, false);
139: return mesh;
140: };
141: void Builder::readFault(Obj<Builder::Mesh> mesh, const std::string& filename)
142: {
143: throw ALE::Exception("Not implemented for optimized sieves");
144: };
145: #if 0
146: void Builder::readFault(Obj<Builder::Mesh> mesh, const std::string& filename)
147: {
148: const int numCells = mesh->heightStratum(0)->size();
149: PetscViewer viewer;
150: FILE *f;
151: char buf[2048];
152: PetscInt numPsets;
155: if (mesh->commRank() != 0) return;
156: PetscViewerCreate(PETSC_COMM_SELF, &viewer);CHKERRXX(ierr);
157: PetscViewerSetType(viewer, PETSCVIEWERASCII);CHKERRXX(ierr);
158: PetscViewerFileSetMode(viewer, FILE_MODE_READ);CHKERRXX(ierr);
159: PetscViewerFileSetName(viewer, filename.c_str());CHKERRXX(ierr);
160: PetscViewerASCIIGetPointer(viewer, &f);CHKERRXX(ierr);
161: // Read header
162: if (!fgets(buf, 2048, f)) throw ALE::Exception("File ended prematurely");
163: // Check file type
164: const char *x = strtok(buf, " ");
165: std::string fileType("pset");
166: if (fileType != x) throw ALE::Exception("Invalid file type");
167: // Ignore set type
168: x = strtok(NULL, " ");
169: // Number of psets
170: x = strtok(NULL, " ");
171: numPsets = atoi(x);
172: for (PetscInt p = 0; p < numPsets; ++p) {
173: if (!fgets(buf, 2048, f)) throw ALE::Exception("File ended prematurely");
174: // Read name
175: x = strtok(buf, " ");
176: const Obj<Mesh::int_section_type>& fault = mesh->getIntSection(x);
177: // Vertices per line
178: x = strtok(NULL, " ");
179: const PetscInt vertsPerLine = atoi(x);
180: // Total vertices
181: x = strtok(NULL, " ");
182: const PetscInt totalVerts = atoi(x);
184: for (PetscInt v = 0; v < totalVerts; ++v) {
185: if (v%vertsPerLine == 0) {
186: if (!fgets(buf, 2048, f)) throw ALE::Exception("File ended prematurely");
187: x = strtok(buf, " ");
188: } else {
189: x = strtok(NULL, " ");
190: }
191: const int vv = atoi(x) + numCells - 1;
193: fault->setFiberDimension(vv, 1);
194: }
195: fault->allocatePoint();
196: }
197: PetscViewerDestroy(&viewer);CHKERRXX(ierr);
198: };
199: #endif
200: }
201: }