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