Actual source code: meshpylith.c

  1: #include "src/dm/mesh/meshpylith.h"   /*I      "petscmesh.h"   I*/

  3: #include<list>

  5: namespace ALE {
  6:   namespace PyLith {
  7:     using ALE::Mesh;
  8:     //
  9:     // Builder methods
 10:     //
 11:     inline void Builder::ignoreComments(char *buf, PetscInt bufSize, FILE *f) {
 12:       while((fgets(buf, bufSize, f) != NULL) && ((buf[0] == '#') || (buf[0] == '\0'))) {}
 13:     };
 14:     void Builder::readConnectivity(MPI_Comm comm, const std::string& filename, int& corners, const bool useZeroBase, int& numElements, int *vertices[], int *materials[]) {
 15:       PetscViewer    viewer;
 16:       FILE          *f;
 17:       PetscInt       maxCells = 1024, cellCount = 0;
 18:       PetscInt      *verts;
 19:       PetscInt      *mats;
 20:       char           buf[2048];
 21:       PetscInt       c;
 22:       PetscInt       commRank;

 25:       MPI_Comm_rank(comm, &commRank);
 26:       if (commRank != 0) return;
 27:       PetscViewerCreate(PETSC_COMM_SELF, &viewer);
 28:       PetscViewerSetType(viewer, PETSC_VIEWER_ASCII);
 29:       PetscViewerFileSetMode(viewer, FILE_MODE_READ);
 30:       PetscViewerFileSetName(viewer, filename.c_str());
 31:       PetscViewerASCIIGetPointer(viewer, &f);
 32:       /* Ignore comments */
 33:       ignoreComments(buf, 2048, f);
 34:       do {
 35:         const char *v = strtok(buf, " ");
 36:         int         elementType;

 38:         if (cellCount == maxCells) {
 39:           PetscInt *vtmp, *mtmp;

 41:           vtmp = verts;
 42:           mtmp = mats;
 43:           PetscMalloc2(maxCells*2*corners,PetscInt,&verts,maxCells*2,PetscInt,&mats);
 44:           PetscMemcpy(verts, vtmp, maxCells*corners * sizeof(PetscInt));
 45:           PetscMemcpy(mats,  mtmp, maxCells         * sizeof(PetscInt));
 46:           PetscFree2(vtmp,mtmp);
 47:           maxCells *= 2;
 48:         }
 49:         /* Ignore cell number */
 50:         v = strtok(NULL, " ");
 51:         /* Get element type */
 52:         elementType = atoi(v);
 53:         if (elementType == 1) {
 54:           corners = 8;
 55:         } else if (elementType == 5) {
 56:           corners = 4;
 57:         } else {
 58:           ostringstream msg;

 60:           msg << "We do not accept element type " << elementType << " right now";
 61:           throw ALE::Exception(msg.str().c_str());
 62:         }
 63:         if (cellCount == 0) {
 64:           PetscMalloc2(maxCells*corners,PetscInt,&verts,maxCells,PetscInt,&mats);
 65:         }
 66:         v = strtok(NULL, " ");
 67:         /* Store material type */
 68:         mats[cellCount] = atoi(v);
 69:         v = strtok(NULL, " ");
 70:         /* Ignore infinite domain element code */
 71:         v = strtok(NULL, " ");
 72:         for(c = 0; c < corners; c++) {
 73:           int vertex = atoi(v);
 74: 
 75:           if (!useZeroBase) vertex -= 1;
 76:           verts[cellCount*corners+c] = vertex;
 77:           v = strtok(NULL, " ");
 78:         }
 79:         cellCount++;
 80:       } while(fgets(buf, 2048, f) != NULL);
 81:       PetscViewerDestroy(viewer);
 82:       numElements = cellCount;
 83:       *vertices   = verts;
 84:       *materials  = mats;
 85:     };
 86:     void Builder::readCoordinates(MPI_Comm comm, const std::string& filename, const int dim, int& numVertices, double *coordinates[]) {
 87:       PetscViewer    viewer;
 88:       FILE          *f;
 89:       PetscInt       maxVerts = 1024, vertexCount = 0;
 90:       PetscScalar   *coords;
 91:       double         scaleFactor = 1.0;
 92:       char           buf[2048];
 93:       PetscInt       c;
 94:       PetscInt       commRank;

 97:       MPI_Comm_rank(comm, &commRank);
 98:       if (commRank == 0) {
 99:         PetscViewerCreate(PETSC_COMM_SELF, &viewer);
100:         PetscViewerSetType(viewer, PETSC_VIEWER_ASCII);
101:         PetscViewerFileSetMode(viewer, FILE_MODE_READ);
102:         PetscViewerFileSetName(viewer, filename.c_str());
103:         PetscViewerASCIIGetPointer(viewer, &f);
104:         /* Ignore comments */
105:         ignoreComments(buf, 2048, f);
106:         PetscMalloc(maxVerts*dim * sizeof(PetscScalar), &coords);
107:         /* Read units */
108:         const char *units = strtok(buf, " ");
109:         if (strcmp(units, "coord_units")) {
110:           throw ALE::Exception("Invalid coordinate units line");
111:         }
112:         units = strtok(NULL, " ");
113:         if (strcmp(units, "=")) {
114:           throw ALE::Exception("Invalid coordinate units line");
115:         }
116:         units = strtok(NULL, " ");
117:         if (!strcmp(units, "km")) {
118:           /* Should use Pythia to do units conversion */
119:           scaleFactor = 1000.0;
120:         }
121:         /* Ignore comments */
122:         ignoreComments(buf, 2048, f);
123:         do {
124:           const char *x = strtok(buf, " ");

126:           if (vertexCount == maxVerts) {
127:             PetscScalar *ctmp;

129:             ctmp = coords;
130:             PetscMalloc(maxVerts*2*dim * sizeof(PetscScalar), &coords);
131:             PetscMemcpy(coords, ctmp, maxVerts*dim * sizeof(PetscScalar));
132:             PetscFree(ctmp);
133:             maxVerts *= 2;
134:           }
135:           /* Ignore vertex number */
136:           x = strtok(NULL, " ");
137:           for(c = 0; c < dim; c++) {
138:             coords[vertexCount*dim+c] = atof(x)*scaleFactor;
139:             x = strtok(NULL, " ");
140:           }
141:           vertexCount++;
142:         } while(fgets(buf, 2048, f) != NULL);
143:         PetscViewerDestroy(viewer);
144:         numVertices = vertexCount;
145:         *coordinates = coords;
146:       }
147:     };
148:     // numSplit is the number of split nodes
149:     // splitInd[] is an array of numSplit pairs, <element, vertex>
150:     // splitValues[] is an array of numSplit*dim displacements
151:     void Builder::readSplit(MPI_Comm comm, const std::string& filename, const int dim, const bool useZeroBase, int& numSplit, int *splitInd[], double *splitValues[]) {
152:       PetscViewer    viewer;
153:       FILE          *f;
154:       PetscInt       maxSplit = 1024, splitCount = 0;
155:       PetscInt      *splitId;
156:       PetscScalar   *splitVal;
157:       char           buf[2048];
158:       PetscInt       c;
159:       PetscInt       commRank;

162:       MPI_Comm_rank(comm, &commRank);
163:       if (dim != 3) {
164:         throw ALE::Exception("PyLith only works in 3D");
165:       }
166:       if (commRank != 0) return;
167:       PetscViewerCreate(PETSC_COMM_SELF, &viewer);
168:       PetscViewerSetType(viewer, PETSC_VIEWER_ASCII);
169:       PetscViewerFileSetMode(viewer, FILE_MODE_READ);
170:       PetscExceptionTry1(PetscViewerFileSetName(viewer, filename.c_str()), PETSC_ERR_FILE_OPEN);
171:       if (PetscExceptionValue(ierr)) {
172:         // this means that a caller above me has also tryed this exception so I don't handle it here, pass it up
173:       } else if (PetscExceptionCaught(ierr,PETSC_ERR_FILE_OPEN)) {
174:         // File does not exist
175:         return;
176:       }
177:       PetscViewerASCIIGetPointer(viewer, &f);
178:       /* Ignore comments */
179:       ignoreComments(buf, 2048, f);
180:       PetscMalloc2(maxSplit*2,PetscInt,&splitId,maxSplit*dim,PetscScalar,&splitVal);
181:       do {
182:         const char *s = strtok(buf, " ");

184:         if (splitCount == maxSplit) {
185:           PetscInt    *sitmp;
186:           PetscScalar *svtmp;

188:           sitmp = splitId;
189:           svtmp = splitVal;
190:           PetscMalloc2(maxSplit*2*2,PetscInt,&splitId,maxSplit*dim*2,PetscScalar,&splitVal);
191:           PetscMemcpy(splitId,  sitmp, maxSplit*2   * sizeof(PetscInt));
192:           PetscMemcpy(splitVal, svtmp, maxSplit*dim * sizeof(PetscScalar));
193:           PetscFree2(sitmp,svtmp);
194:           maxSplit *= 2;
195:         }
196:         /* Get element number */
197:         int elem = atoi(s);
198:         if (!useZeroBase) elem -= 1;
199:         splitId[splitCount*2+0] = elem;
200:         s = strtok(NULL, " ");
201:         /* Get node number */
202:         int node = atoi(s);
203:         if (!useZeroBase) node -= 1;
204:         splitId[splitCount*2+1] = node;
205:         s = strtok(NULL, " ");
206:         /* Ignore load history number */
207:         s = strtok(NULL, " ");
208:         /* Get split values */
209:         for(c = 0; c < dim; c++) {
210:           splitVal[splitCount*dim+c] = atof(s);
211:           s = strtok(NULL, " ");
212:         }
213:         splitCount++;
214:       } while(fgets(buf, 2048, f) != NULL);
215:       PetscViewerDestroy(viewer);
216:       numSplit     = splitCount;
217:       *splitInd    = splitId;
218:       *splitValues = splitVal;
219:     };
220:     void Builder::buildSplit(const Obj<pair_section_type>& splitField, int numCells, int numSplit, int splitInd[], double splitVals[]) {
221:       const pair_section_type::patch_type                     patch = 0;
222:       pair_section_type::value_type                          *values;
223:       std::map<pair_section_type::point_type, std::set<int> > elem2index;
224:       int                                                     numValues = 0;

226:       splitField->setName("split");
227:       for(int e = 0; e < numSplit; e++) {
228:         splitField->addFiberDimension(patch, splitInd[e*2+0], 1);
229:         elem2index[splitInd[e*2+0]].insert(e);
230:       }
231:       splitField->allocate();
232:       for(std::map<pair_section_type::point_type, std::set<int> >::const_iterator e_iter = elem2index.begin(); e_iter != elem2index.end(); ++e_iter) {
233:         numValues = std::max(numValues, (int) e_iter->second.size());
234:       }
235:       values = new pair_section_type::value_type[numValues];
236:       for(std::map<pair_section_type::point_type, std::set<int> >::const_iterator e_iter = elem2index.begin(); e_iter != elem2index.end(); ++e_iter) {
237:         const pair_section_type::point_type& e = e_iter->first;
238:         int                                  k = 0;

240:         for(std::set<int>::const_iterator i_iter = e_iter->second.begin(); i_iter != e_iter->second.end(); ++i_iter, ++k) {
241:           const int& i = *i_iter;

243:           if (k >= numValues) {throw ALE::Exception("Invalid split node input");}
244:           values[k].first    = splitInd[i*2+1] + numCells;
245:           values[k].second.x = splitVals[i*3+0];
246:           values[k].second.y = splitVals[i*3+1];
247:           values[k].second.z = splitVals[i*3+2];
248:         }
249:         splitField->update(patch, e, values);
250:       }
251:       delete [] values;
252:     };
253:     void Builder::readTractions(MPI_Comm comm, const std::string& filename, const int dim, const int& corners, const bool useZeroBase, int& numTractions, int& vertsPerFace, int *tractionVertices[], double *tractionValues[]) {
254:       PetscViewer    viewer;
255:       FILE          *f;
256:       PetscInt       maxTractions = 1024, tractionCount = 0;
257:       PetscInt      *tractionVerts;
258:       PetscScalar   *tractionVals;
259:       double         scaleFactor = 1.0;
260:       char           buf[2048];
261:       PetscInt       c;
262:       PetscInt       commRank;

265:       MPI_Comm_rank(comm, &commRank);
266:       if (dim != 3) {
267:         throw ALE::Exception("PyLith only works in 3D");
268:       }
269:       if (commRank != 0) return;
270:       PetscViewerCreate(PETSC_COMM_SELF, &viewer);
271:       PetscViewerSetType(viewer, PETSC_VIEWER_ASCII);
272:       PetscViewerFileSetMode(viewer, FILE_MODE_READ);
273:       PetscExceptionTry1(PetscViewerFileSetName(viewer, filename.c_str()), PETSC_ERR_FILE_OPEN);
274:       if (PetscExceptionValue(ierr)) {
275:         // this means that a caller above me has also tryed this exception so I don't handle it here, pass it up
276:       } else if (PetscExceptionCaught(ierr,PETSC_ERR_FILE_OPEN)) {
277:         // File does not exist
278:         return;
279:       }
280:       /* Logic right now is only good for linear tets and hexes, and should be fixed in the future. */
281:       if (corners == 4) {
282:         vertsPerFace = 3;
283:       } else if (corners == 8) {
284:         vertsPerFace = 4;
285:       } else {
286:         throw ALE::Exception("Unrecognized element type");
287:       }

289:       PetscViewerASCIIGetPointer(viewer, &f);
290:       /* Ignore comments */
291:       ignoreComments(buf, 2048, f);
292:       /* Read units */
293:       const char *units = strtok(buf, " ");
294:       if (strcmp(units, "traction_units")) {
295:         throw ALE::Exception("Invalid traction units line");
296:       }
297:       units = strtok(NULL, " ");
298:       if (strcmp(units, "=")) {
299:         throw ALE::Exception("Invalid traction units line");
300:       }
301:       units = strtok(NULL, " ");
302:       if (!strcmp(units, "MPa")) {
303:         /* Should use Pythia to do units conversion */
304:         scaleFactor = 1.0e6;
305:       }
306:       /* Ignore comments */
307:       ignoreComments(buf, 2048, f);
308:       // Allocate memory.
309:       PetscMalloc2(maxTractions*vertsPerFace,PetscInt,&tractionVerts,maxTractions*dim,PetscScalar,&tractionVals);
310:       do {
311:         const char *s = strtok(buf, " ");

313:         if (tractionCount == maxTractions) {
314:           PetscInt    *titmp;
315:           PetscScalar *tvtmp;

317:           titmp = tractionVerts;
318:           tvtmp = tractionVals;
319:           PetscMalloc2(maxTractions*vertsPerFace*2,PetscInt,&tractionVerts,maxTractions*dim*2,PetscScalar,&tractionVals);
320:           PetscMemcpy(tractionVerts,  titmp, maxTractions*vertsPerFace   * sizeof(PetscInt));
321:           PetscMemcpy(tractionVals, tvtmp, maxTractions*dim * sizeof(PetscScalar));
322:           PetscFree2(titmp,tvtmp);
323:           maxTractions *= 2;
324:         }
325:         /* Get vertices */
326:         int v1 = atoi(s);
327:         if (!useZeroBase) v1 -= 1;
328:         tractionVerts[tractionCount*vertsPerFace+0] = v1;
329:         s = strtok(NULL, " ");
330:         int v2 = atoi(s);
331:         if (!useZeroBase) v2 -= 1;
332:         tractionVerts[tractionCount*vertsPerFace+1] = v2;
333:         s = strtok(NULL, " ");
334:         int v3 = atoi(s);
335:         if (!useZeroBase) v3 -= 1;
336:         tractionVerts[tractionCount*vertsPerFace+2] = v3;
337:         s = strtok(NULL, " ");
338:         /* Get traction values */
339:         for(c = 0; c < dim; c++) {
340:           tractionVals[tractionCount*dim+c] = atof(s);
341:           s = strtok(NULL, " ");
342:         }
343:         tractionCount++;
344:       } while(fgets(buf, 2048, f) != NULL);
345:       PetscViewerDestroy(viewer);
346:       numTractions     = tractionCount;
347:       *tractionVertices    = tractionVerts;
348:       *tractionValues = tractionVals;
349:     };
350:     void Builder::buildTractions(const Obj<real_section_type>& tractionField, const Obj<topology_type>& boundaryTopology, int numCells, int numTractions, int vertsPerFace, int tractionVertices[], double tractionValues[]) {
351:       const real_section_type::patch_type                  patch = 0;
352:       real_section_type::value_type                        values[3];
353:       // Make boundary topology
354:       Obj<sieve_type> boundarySieve = new sieve_type(tractionField->comm(), tractionField->debug());

356:       ALE::New::SieveBuilder<sieve_type>::buildTopology(boundarySieve, 2, numTractions, tractionVertices, 0, false, vertsPerFace, numCells);
357:       boundaryTopology->setPatch(patch, boundarySieve);
358:       boundaryTopology->stratify();
359:       // Make traction field
360:       tractionField->setName("traction");
361:       tractionField->setTopology(boundaryTopology);
362:       tractionField->setFiberDimensionByHeight(patch, 0, 3);
363:       tractionField->allocate();
364:       const Obj<topology_type::label_sequence>& faces = boundaryTopology->heightStratum(patch, 0);
365:       int k = 0;

367:       for(topology_type::label_sequence::iterator f_iter = faces->begin(); f_iter != faces->end(); ++f_iter) {
368:         const topology_type::point_type& face = *f_iter;

370:         values[0] = tractionValues[k*3+0];
371:         values[1] = tractionValues[k*3+1];
372:         values[2] = tractionValues[k*3+2];
373:         k++;
374:         tractionField->update(patch, face, values);
375:       }
376:     };
377:     void Builder::buildMaterials(const Obj<int_section_type>& matField, const int materials[]) {
378:       const int_section_type::patch_type        patch    = 0;
379:       const Obj<topology_type::label_sequence>& elements = matField->getTopology()->heightStratum(patch, 0);

381:       matField->setName("material");
382:       matField->setFiberDimensionByHeight(patch, 0, 1);
383:       matField->allocate();
384:       for(topology_type::label_sequence::iterator e_iter = elements->begin(); e_iter != elements->end(); ++e_iter) {
385:         matField->update(patch, *e_iter, &materials[*e_iter]);
386:       }
387:     };
388:     Obj<Mesh> Builder::readMesh(MPI_Comm comm, const int dim, const std::string& basename, const bool useZeroBase = false, const bool interpolate = false, const int debug = 0) {
389:       Obj<Mesh>             mesh     = new Mesh(comm, dim, debug);
390:       Obj<sieve_type>       sieve    = new sieve_type(comm, debug);
391:       Obj<topology_type>    topology = new topology_type(comm, debug);
392:       int    *cells, *materials;
393:       double *coordinates;
394:       int     numCells = 0, numVertices = 0, numCorners = dim+1;

396:       ALE::PyLith::Builder::readConnectivity(comm, basename+".connect", numCorners, useZeroBase, numCells, &cells, &materials);
397:       ALE::PyLith::Builder::readCoordinates(comm, basename+".coord", dim, numVertices, &coordinates);
398:       ALE::New::SieveBuilder<sieve_type>::buildTopology(sieve, dim, numCells, cells, numVertices, interpolate, numCorners);
399:       sieve->stratify();
400:       topology->setPatch(0, sieve);
401:       topology->stratify();
402:       mesh->setTopology(topology);
403:       ALE::New::SieveBuilder<sieve_type>::buildCoordinates(mesh->getRealSection("coordinates"), dim, coordinates);
404:       Obj<int_section_type> material = mesh->getIntSection("material");
405:       buildMaterials(material, materials);
406:       Obj<ALE::Mesh::pair_section_type> split = createSplit(mesh, basename, useZeroBase);
407:       if (!split.isNull()) {mesh->setPairSection("split", split);}
408:       Obj<ALE::Mesh::real_section_type> traction = createTraction(mesh, basename, useZeroBase);
409:       if (!traction.isNull()) {mesh->setRealSection("traction", traction);}
410:       return mesh;
411:     };
412:     Obj<Builder::pair_section_type> Builder::createSplit(const Obj<Mesh>& mesh, const std::string& basename, const bool useZeroBase = false) {
413:       Obj<pair_section_type> split = NULL;
414:       MPI_Comm comm     = mesh->comm();
415:       int      dim      = mesh->getDimension();
416:       int      numCells = mesh->getTopology()->heightStratum(0, 0)->size();
417:       int     *splitInd;
418:       double  *splitValues;
419:       int      numSplit = 0, hasSplit;

421:       ALE::PyLith::Builder::readSplit(comm, basename+".split", dim, useZeroBase, numSplit, &splitInd, &splitValues);
422:       MPI_Allreduce(&numSplit, &hasSplit, 1, MPI_INT, MPI_MAX, comm);
423:       if (hasSplit) {
424:         split = new pair_section_type(mesh->getTopology());
425:         buildSplit(split, numCells, numSplit, splitInd, splitValues);
426:       }
427:       return split;
428:     };
429:     Obj<Mesh::real_section_type> Builder::createTraction(const Obj<Mesh>& mesh, const std::string& basename, const bool useZeroBase = false) {
430:       Obj<real_section_type> traction = NULL;
431:       MPI_Comm comm       = mesh->comm();
432:       int      debug      = mesh->debug();
433:       int      dim        = mesh->getDimension();
434:       int      numCells   = mesh->getTopology()->heightStratum(0, 0)->size();
435:       int      numCorners = mesh->getTopology()->getPatch(0)->cone(*mesh->getTopology()->heightStratum(0, 0)->begin())->size();
436:       int     *tractionVertices;
437:       double  *tractionValues;
438:       int      numTractions = 0, vertsPerFace = 0, hasTractions;

440:       ALE::PyLith::Builder::readTractions(comm, basename+".traction", dim, numCorners, useZeroBase, numTractions, vertsPerFace, &tractionVertices, &tractionValues);
441:       MPI_Allreduce(&numTractions, &hasTractions, 1, MPI_INT, MPI_MAX, comm);
442:       if (hasTractions) {
443:         Obj<topology_type> boundaryTopology = new topology_type(mesh->comm(), debug);

445:         traction = new real_section_type(boundaryTopology);
446:         buildTractions(traction, boundaryTopology, numCells, numTractions, vertsPerFace, tractionVertices, tractionValues);
447:       }
448:       return traction;
449:     };
450:     //
451:     // Viewer methods
452:     //
455:     PetscErrorCode Viewer::writeVertices(const Obj<Mesh>& mesh, PetscViewer viewer) {
456:       Obj<Builder::real_section_type> coordinates  = mesh->getRealSection("coordinates");
457:       //Mesh::section_type::patch_type patch;
458:       //const double  *array = coordinates->restrict(Mesh::section_type::patch_type());
459:       //int            dim = mesh->getDimension();
460:       //int            numVertices;

464: #if 0
465:       //FIX:
466:       if (vertexBundle->getGlobalOffsets()) {
467:         numVertices = vertexBundle->getGlobalOffsets()[mesh->commSize()];
468:       } else {
469:         numVertices = mesh->getTopology()->depthStratum(0)->size();
470:       }
471:       PetscViewerASCIIPrintf(viewer,"#\n");
472:       PetscViewerASCIIPrintf(viewer,"coord_units = m\n");
473:       PetscViewerASCIIPrintf(viewer,"#\n");
474:       PetscViewerASCIIPrintf(viewer,"#\n");
475:       PetscViewerASCIIPrintf(viewer,"#  Node      X-coord           Y-coord           Z-coord\n");
476:       PetscViewerASCIIPrintf(viewer,"#\n");
477:       if (mesh->commRank() == 0) {
478:         int numLocalVertices = mesh->getTopology()->depthStratum(0)->size();
479:         int vertexCount = 1;

481:         for(int v = 0; v < numLocalVertices; v++) {
482:           PetscViewerASCIIPrintf(viewer,"%7D ", vertexCount++);
483:           for(int d = 0; d < dim; d++) {
484:             if (d > 0) {
485:               PetscViewerASCIIPrintf(viewer," ");
486:             }
487:             PetscViewerASCIIPrintf(viewer,"% 16.8E", array[v*dim+d]);
488:           }
489:           PetscViewerASCIIPrintf(viewer,"\n");
490:         }
491:         for(int p = 1; p < mesh->commSize(); p++) {
492:           double    *remoteCoords;
493:           MPI_Status status;

495:           MPI_Recv(&numLocalVertices, 1, MPI_INT, p, 1, mesh->comm(), &status);
496:           PetscMalloc(numLocalVertices*dim * sizeof(double), &remoteCoords);
497:           MPI_Recv(remoteCoords, numLocalVertices*dim, MPI_DOUBLE, p, 1, mesh->comm(), &status);
498:           for(int v = 0; v < numLocalVertices; v++) {
499:             PetscViewerASCIIPrintf(viewer,"%7D   ", vertexCount++);
500:             for(int d = 0; d < dim; d++) {
501:               if (d > 0) {
502:                 PetscViewerASCIIPrintf(viewer, " ");
503:               }
504:               PetscViewerASCIIPrintf(viewer, "% 16.8E", remoteCoords[v*dim+d]);
505:             }
506:             PetscViewerASCIIPrintf(viewer, "\n");
507:           }
508:         }
509:       } else {
510:         Obj<Mesh::bundle_type> globalOrder = coordinates->getGlobalOrder();
511:         Obj<Mesh::field_type::order_type::coneSequence> cone = globalOrder->getPatch(patch);
512:         const int *offsets = coordinates->getGlobalOffsets();
513:         int        numLocalVertices = (offsets[mesh->commRank()+1] - offsets[mesh->commRank()])/dim;
514:         double    *localCoords;
515:         int        k = 0;

517:         PetscMalloc(numLocalVertices*dim * sizeof(double), &localCoords);
518:         for(Mesh::field_type::order_type::coneSequence::iterator p_iter = cone->begin(); p_iter != cone->end(); ++p_iter) {
519:           int dim = globalOrder->getFiberDimension(patch, *p_iter);

521:           if (dim > 0) {
522:             int offset = coordinates->getFiberOffset(patch, *p_iter);

524:             for(int i = offset; i < offset+dim; ++i) {
525:               localCoords[k++] = array[i];
526:             }
527:           }
528:         }
529:         if (k != numLocalVertices*dim) {
530:           SETERRQ2(PETSC_ERR_PLIB, "Invalid number of coordinates to send %d should be %d", k, numLocalVertices*dim);
531:         }
532:         MPI_Send(&numLocalVertices, 1, MPI_INT, 0, 1, mesh->comm());
533:         MPI_Send(localCoords, numLocalVertices*dim, MPI_DOUBLE, 0, 1, mesh->comm());
534:         PetscFree(localCoords);
535:       }
536: #endif
537:       return(0);
538:     };
541:     PetscErrorCode Viewer::writeElements(const Obj<Mesh>& mesh, const Obj<Builder::int_section_type>& materialField, PetscViewer viewer) {
542:       Obj<Mesh::topology_type> topology = mesh->getTopology();
543: #if 0
544:       Obj<Mesh::sieve_type::traits::heightSequence> elements = topology->heightStratum(0);
545:       Obj<Mesh::bundle_type> elementBundle = mesh->getBundle(topology->depth());
546:       Obj<Mesh::bundle_type> vertexBundle = mesh->getBundle(0);
547:       Obj<Mesh::bundle_type> globalVertex = vertexBundle->getGlobalOrder();
548:       Obj<Mesh::bundle_type> globalElement = elementBundle->getGlobalOrder();
549:       Mesh::bundle_type::patch_type patch;
550:       std::string    orderName("element");
551:       bool           hasMaterial  = !materialField.isNull();
552:       int            dim  = mesh->getDimension();
553:       int            corners = topology->nCone(*elements->begin(), topology->depth())->size();
554:       int            elementType = -1;

558:       if (dim != 3) {
559:         SETERRQ(PETSC_ERR_SUP, "PyLith only supports 3D meshes.");
560:       }
561:       if (corners == 4) {
562:         // Linear tetrahedron
563:         elementType = 5;
564:       } else if (corners == 8) {
565:         // Linear hexahedron
566:         elementType = 1;
567:       } else {
568:         SETERRQ1(PETSC_ERR_SUP, "PyLith Error: Unsupported number of elements vertices: %d", corners);
569:       }
570:       PetscViewerASCIIPrintf(viewer,"#\n");
571:       PetscViewerASCIIPrintf(viewer,"#     N ETP MAT INF     N1     N2     N3     N4     N5     N6     N7     N8\n");
572:       PetscViewerASCIIPrintf(viewer,"#\n");
573:       if (mesh->commRank() == 0) {
574:         int elementCount = 1;

576:         for(Mesh::sieve_type::traits::heightSequence::iterator e_itor = elements->begin(); e_itor != elements->end(); ++e_itor) {
577:           Obj<Mesh::bundle_type::order_type::coneSequence> cone = vertexBundle->getPatch(orderName, *e_itor);

579:           PetscViewerASCIIPrintf(viewer, "%7d %3d", elementCount++, elementType);
580:           if (hasMaterial) {
581:             // No infinite elements
582:             PetscViewerASCIIPrintf(viewer, " %3d %3d", (int) materialField->restrict(patch, *e_itor)[0], 0);
583:           } else {
584:             // No infinite elements
585:             PetscViewerASCIIPrintf(viewer, " %3d %3d", 1, 0);
586:           }
587:           for(Mesh::bundle_type::order_type::coneSequence::iterator c_itor = cone->begin(); c_itor != cone->end(); ++c_itor) {
588:             PetscViewerASCIIPrintf(viewer, " %6d", globalVertex->getIndex(patch, *c_itor).prefix+1);
589:           }
590:           PetscViewerASCIIPrintf(viewer, "\n");
591:         }
592:         for(int p = 1; p < mesh->commSize(); p++) {
593:           int         numLocalElements;
594:           int        *remoteVertices;
595:           MPI_Status  status;

597:           MPI_Recv(&numLocalElements, 1, MPI_INT, p, 1, mesh->comm(), &status);
598:           PetscMalloc(numLocalElements*(corners+1) * sizeof(int), &remoteVertices);
599:           MPI_Recv(remoteVertices, numLocalElements*(corners+1), MPI_INT, p, 1, mesh->comm(), &status);
600:           for(int e = 0; e < numLocalElements; e++) {
601:             // Only linear tetrahedra, material, no infinite elements
602:             int mat = remoteVertices[e*(corners+1)+corners];

604:             PetscViewerASCIIPrintf(viewer, "%7d %3d %3d %3d", elementCount++, elementType, mat, 0);
605:             for(int c = 0; c < corners; c++) {
606:               PetscViewerASCIIPrintf(viewer, " %6d", remoteVertices[e*(corners+1)+c]);
607:             }
608:             PetscViewerASCIIPrintf(viewer, "\n");
609:           }
610:           PetscFree(remoteVertices);
611:         }
612:       } else {
613:         const int *offsets = elementBundle->getGlobalOffsets();
614:         int        numLocalElements = offsets[mesh->commRank()+1] - offsets[mesh->commRank()];
615:         int       *localVertices;
616:         int        k = 0;

618:         PetscMalloc(numLocalElements*(corners+1) * sizeof(int), &localVertices);
619:         for(Mesh::sieve_type::traits::heightSequence::iterator e_itor = elements->begin(); e_itor != elements->end(); ++e_itor) {
620:           Obj<Mesh::bundle_type::order_type::coneSequence> cone = vertexBundle->getPatch(orderName, *e_itor);

622:           if (globalElement->getFiberDimension(patch, *e_itor) > 0) {
623:             for(Mesh::bundle_type::order_type::coneSequence::iterator c_itor = cone->begin(); c_itor != cone->end(); ++c_itor) {
624:               localVertices[k++] = globalVertex->getIndex(patch, *c_itor).prefix;
625:             }
626:             if (hasMaterial) {
627:               localVertices[k++] = (int) materialField->restrict(patch, *e_itor)[0];
628:             } else {
629:               localVertices[k++] = 1;
630:             }
631:           }
632:         }
633:         if (k != numLocalElements*corners) {
634:           SETERRQ2(PETSC_ERR_PLIB, "Invalid number of vertices to send %d should be %d", k, numLocalElements*corners);
635:         }
636:         MPI_Send(&numLocalElements, 1, MPI_INT, 0, 1, mesh->comm());
637:         MPI_Send(localVertices, numLocalElements*(corners+1), MPI_INT, 0, 1, mesh->comm());
638:         PetscFree(localVertices);
639:       }
640: #endif
641:       return(0);
642:     };
645:     PetscErrorCode Viewer::writeVerticesLocal(const Obj<Mesh>& mesh, PetscViewer viewer) {
646:       const Builder::real_section_type::patch_type       patch       = 0;
647:       const Obj<Builder::real_section_type>&             coordinates = mesh->getRealSection("coordinates");
648:       const Obj<Builder::topology_type>&                 topology    = mesh->getTopology();
649:       const Obj<Builder::topology_type::label_sequence>& vertices    = topology->depthStratum(patch, 0);
650:       const Obj<Mesh::numbering_type>&                   vNumbering  = mesh->getFactory()->getLocalNumbering(topology, patch, 0);
651:       int            embedDim = coordinates->getFiberDimension(patch, *vertices->begin());

655:       PetscViewerASCIIPrintf(viewer,"#\n");
656:       PetscViewerASCIIPrintf(viewer,"coord_units = m\n");
657:       PetscViewerASCIIPrintf(viewer,"#\n");
658:       PetscViewerASCIIPrintf(viewer,"#\n");
659:       PetscViewerASCIIPrintf(viewer,"#  Node      X-coord           Y-coord           Z-coord\n");
660:       PetscViewerASCIIPrintf(viewer,"#\n");

662:       for(Mesh::topology_type::label_sequence::iterator v_iter = vertices->begin(); v_iter != vertices->end(); ++v_iter) {
663:         const Builder::real_section_type::value_type *array = coordinates->restrict(patch, *v_iter);

665:         PetscViewerASCIIPrintf(viewer, "%7D ", vNumbering->getIndex(*v_iter)+1);
666:         for(int d = 0; d < embedDim; d++) {
667:           if (d > 0) {
668:             PetscViewerASCIIPrintf(viewer, " ");
669:           }
670:           PetscViewerASCIIPrintf(viewer, "% 16.8E", array[d]);
671:         }
672:         PetscViewerASCIIPrintf(viewer, "\n");
673:       }
674:       return(0);
675:     };
678:     PetscErrorCode Viewer::writeElementsLocal(const Obj<Mesh>& mesh, const Obj<Builder::int_section_type>& materialField, PetscViewer viewer) {
679:       const Mesh::topology_type::patch_type           patch      = 0;
680:       const Obj<Mesh::topology_type>&                 topology   = mesh->getTopology();
681:       const Obj<Mesh::sieve_type>&                    sieve      = topology->getPatch(patch);
682:       const Obj<Mesh::topology_type::label_sequence>& elements   = topology->heightStratum(patch, 0);
683:       const Obj<Mesh::numbering_type>&                eNumbering = mesh->getFactory()->getLocalNumbering(topology, patch, topology->depth());
684:       const Obj<Mesh::numbering_type>&                vNumbering = mesh->getFactory()->getLocalNumbering(topology, patch, 0);
685:       int            dim          = mesh->getDimension();
686:       //int            corners      = sieve->nCone(*elements->begin(), topology->depth())->size();
687:       int            corners      = sieve->cone(*elements->begin())->size();
688:       bool           hasMaterial  = !materialField.isNull();
689:       int            elementType  = -1;

693:       if (dim != 3) {
694:         SETERRQ(PETSC_ERR_SUP, "PyLith only supports 3D meshes.");
695:       }
696:       if (corners == 4) {
697:         // Linear tetrahedron
698:         elementType = 5;
699:       } else if (corners == 8) {
700:         // Linear hexahedron
701:         elementType = 1;
702:       } else {
703:         SETERRQ1(PETSC_ERR_SUP, "PyLith Error: Unsupported number of elements vertices: %d", corners);
704:       }
705:       PetscViewerASCIIPrintf(viewer,"#\n");
706:       PetscViewerASCIIPrintf(viewer,"#     N ETP MAT INF     N1     N2     N3     N4     N5     N6     N7     N8\n");
707:       PetscViewerASCIIPrintf(viewer,"#\n");
708:       for(Mesh::topology_type::label_sequence::iterator e_iter = elements->begin(); e_iter != elements->end(); ++e_iter) {
709:         const Obj<Mesh::sieve_type::traits::coneSequence> cone  = sieve->cone(*e_iter);
710:         Mesh::sieve_type::traits::coneSequence::iterator  begin = cone->begin();
711:         Mesh::sieve_type::traits::coneSequence::iterator  end   = cone->end();

713:         PetscViewerASCIIPrintf(viewer, "%7d %3d", eNumbering->getIndex(*e_iter)+1, elementType);
714:         if (hasMaterial) {
715:           // No infinite elements
716:           PetscViewerASCIIPrintf(viewer, " %3d %3d", (int) materialField->restrict(patch, *e_iter)[0], 0);
717:         } else {
718:           // No infinite elements
719:           PetscViewerASCIIPrintf(viewer, " %3d %3d", 1, 0);
720:         }
721:         for(Mesh::sieve_type::traits::coneSequence::iterator c_iter = begin; c_iter != end; ++c_iter) {
722:           //FIX: Need a global ordering here
723:           PetscViewerASCIIPrintf(viewer, " %6d", vNumbering->getIndex(*c_iter)+1);
724:         }
725:         PetscViewerASCIIPrintf(viewer, "\n");
726:       }
727:       return(0);
728:     };
731:     // The elements seem to be implicitly numbered by appearance, which makes it impossible to
732:     //   number here by bundle, but we can fix it by traversing the elements like the vertices
733:     PetscErrorCode Viewer::writeSplitLocal(const Obj<Mesh>& mesh, const Obj<Builder::pair_section_type>& splitField, PetscViewer viewer) {
734:       const Obj<Mesh::topology_type>&         topology   = mesh->getTopology();
735:       Builder::pair_section_type::patch_type patch      = 0;
736:       const Obj<Mesh::numbering_type>&        eNumbering = mesh->getFactory()->getLocalNumbering(topology, patch, topology->depth());
737:       const Obj<Mesh::numbering_type>&        vNumbering = mesh->getFactory()->getLocalNumbering(topology, patch, 0);

741:       const Builder::pair_section_type::atlas_type::chart_type& chart = splitField->getPatch(patch);

743:       for(Builder::pair_section_type::atlas_type::chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
744:         const Builder::pair_section_type::point_type& e      = *c_iter;
745:         const Builder::pair_section_type::value_type *values = splitField->restrict(patch, e);
746:         const int                                      size   = splitField->getFiberDimension(patch, e);

748:         for(int i = 0; i < size; i++) {
749:           const Builder::pair_section_type::point_type& v      = values[i].first;
750:           const ALE::Mesh::base_type::split_value&      split  = values[i].second;

752:           // No time history
753:           PetscViewerASCIIPrintf(viewer, "%6d %6d 0 %15.9g %15.9g %15.9g\n", eNumbering->getIndex(e)+1, vNumbering->getIndex(v)+1, split.x, split.y, split.z);
754:         }
755:       }
756:       return(0);
757:     };
760:     PetscErrorCode Viewer::writeTractionsLocal(const Obj<Mesh>& mesh, const Obj<Builder::real_section_type>& tractionField, PetscViewer viewer) {
761:       typedef Builder::topology_type     topology_type;
762:       typedef Builder::real_section_type section_type;
763:       const section_type::patch_type            patch      = 0;
764:       const Obj<topology_type>&         boundaryTopology   = tractionField->getTopology();
765:       const Obj<topology_type::sieve_type>&     sieve      = boundaryTopology->getPatch(patch);
766:       const Obj<topology_type::label_sequence>& faces      = boundaryTopology->heightStratum(patch, 0);
767:       const Obj<Mesh::topology_type>&           topology   = mesh->getTopology();
768:       const Obj<Mesh::numbering_type>&          vNumbering = mesh->getFactory()->getLocalNumbering(topology, patch, 0);

772:       PetscViewerASCIIPrintf(viewer,"#\n");
773:       PetscViewerASCIIPrintf(viewer,"traction_units = Pa\n");
774:       PetscViewerASCIIPrintf(viewer,"#\n");
775:       PetscViewerASCIIPrintf(viewer,"#\n");
776:       for(topology_type::label_sequence::iterator f_iter = faces->begin(); f_iter != faces->end(); ++f_iter) {
777:         const topology_type::point_type& face = *f_iter;
778:         const Obj<topology_type::sieve_type::traits::coneSequence>& cone = sieve->cone(face);

780:         for(topology_type::sieve_type::traits::coneSequence::iterator c_iter = cone->begin(); c_iter != cone->end(); ++c_iter) {
781:           const topology_type::point_type& vertex = *c_iter;

783:           PetscViewerASCIIPrintf(viewer, "%6d", vNumbering->getIndex(vertex)+1);
784:           std::cout << vNumbering->getIndex(vertex) << " ("<<vertex<<") ";
785:         }
786:         const section_type::value_type *values = tractionField->restrict(patch, face);

788:         for(int i = 0; i < mesh->getDimension(); ++i) {
789:           if (i > 0) {
790:             PetscViewerASCIIPrintf(viewer, " ");
791:             std::cout << " ";
792:           }
793:           PetscViewerASCIIPrintf(viewer, "%15.9g", values[i]);
794:           std::cout << values[i];
795:         }
796:         PetscViewerASCIIPrintf(viewer,"\n");
797:         std::cout << std::endl;
798:       }
799:       return(0);
800:     };
801:   };
802: };