Actual source code: meshpcice.c

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

  3: namespace ALE {
  4:   namespace PCICE {
  5:     //
  6:     // Builder methods
  7:     //
  8:     void Builder::readConnectivity(MPI_Comm comm, const std::string& filename, int& corners, const bool useZeroBase, int& numElements, int *vertices[]) {
  9:       PetscViewer    viewer;
 10:       FILE          *f;
 11:       PetscInt       numCells, cellCount = 0;
 12:       PetscInt      *verts;
 13:       char           buf[2048];
 14:       PetscInt       c;
 15:       PetscInt       commRank;

 18:       MPI_Comm_rank(comm, &commRank);

 20:       if (commRank != 0) return;
 21:       PetscViewerCreate(PETSC_COMM_SELF, &viewer);
 22:       PetscViewerSetType(viewer, PETSC_VIEWER_ASCII);
 23:       PetscViewerFileSetMode(viewer, FILE_MODE_READ);
 24:       PetscViewerFileSetName(viewer, filename.c_str());
 25:       if (ierr) {
 26:         ostringstream txt;
 27:         txt << "Could not open PCICE connectivity file: " << filename;
 28:         throw ALE::Exception(txt.str().c_str());
 29:       }
 30:       PetscViewerASCIIGetPointer(viewer, &f);
 31:       if (fgets(buf, 2048, f) == NULL) {
 32:         throw ALE::Exception("Invalid connectivity file: Missing number of elements");
 33:       }
 34:       const char *sizes = strtok(buf, " ");
 35:       numCells = atoi(sizes);
 36:       sizes = strtok(NULL, " ");
 37:       if (sizes != NULL) {
 38:         corners = atoi(sizes);
 39:         std::cout << "Reset corners to " << corners << std::endl;
 40:       }
 41:       PetscMalloc(numCells*corners * sizeof(PetscInt), &verts);
 42:       while(fgets(buf, 2048, f) != NULL) {
 43:         const char *v = strtok(buf, " ");
 44: 
 45:         /* Ignore cell number */
 46:         v = strtok(NULL, " ");
 47:         for(c = 0; c < corners; c++) {
 48:           int vertex = atoi(v);
 49: 
 50:           if (!useZeroBase) vertex -= 1;
 51:           verts[cellCount*corners+c] = vertex;
 52:           v = strtok(NULL, " ");
 53:         }
 54:         cellCount++;
 55:       }
 56:       PetscViewerDestroy(viewer);
 57:       numElements = numCells;
 58:       *vertices = verts;
 59:     };
 60:     void Builder::readCoordinates(MPI_Comm comm, const std::string& filename, const int dim, int& numVertices, double *coordinates[]) {
 61:       PetscViewer    viewer;
 62:       FILE          *f;
 63:       PetscInt       numVerts, vertexCount = 0;
 64:       PetscScalar   *coords;
 65:       char           buf[2048];
 66:       PetscInt       c;
 67:       PetscInt       commRank;

 70:       MPI_Comm_rank(comm, &commRank);

 72:       if (commRank != 0) return;
 73:       PetscViewerCreate(PETSC_COMM_SELF, &viewer);
 74:       PetscViewerSetType(viewer, PETSC_VIEWER_ASCII);
 75:       PetscViewerFileSetMode(viewer, FILE_MODE_READ);
 76:       PetscViewerFileSetName(viewer, filename.c_str());
 77:       if (ierr) {
 78:         ostringstream txt;
 79:         txt << "Could not open PCICE coordinate file: " << filename;
 80:         throw ALE::Exception(txt.str().c_str());
 81:       }
 82:       PetscViewerASCIIGetPointer(viewer, &f);
 83:       numVerts = atoi(fgets(buf, 2048, f));
 84:       PetscMalloc(numVerts*dim * sizeof(PetscScalar), &coords);
 85:       while(fgets(buf, 2048, f) != NULL) {
 86:         const char *x = strtok(buf, " ");
 87: 
 88:         /* Ignore vertex number */
 89:         x = strtok(NULL, " ");
 90:         for(c = 0; c < dim; c++) {
 91:           coords[vertexCount*dim+c] = atof(x);
 92:           x = strtok(NULL, " ");
 93:         }
 94:         vertexCount++;
 95:       }
 96:       PetscViewerDestroy(viewer);
 97:       numVertices = numVerts;
 98:       *coordinates = coords;
 99:     };
100:     Obj<Mesh> Builder::readMesh(MPI_Comm comm, const int dim, const std::string& basename, const bool useZeroBase = true, const bool interpolate = true, const int debug = 0) {
101:       return readMesh(comm, dim, basename+".nodes", basename+".lcon", useZeroBase, interpolate, debug);
102:     };
103:     Obj<Mesh> Builder::readMesh(MPI_Comm comm, const int dim, const std::string& coordFilename, const std::string& adjFilename, const bool useZeroBase = true, const bool interpolate = true, const int debug = 0) {
104:       Obj<Mesh>          mesh     = new Mesh(comm, dim, debug);
105:       Obj<sieve_type>    sieve    = new sieve_type(comm, debug);
106:       Obj<topology_type> topology = new topology_type(comm, debug);
107:       int    *cells = NULL;
108:       double *coordinates = NULL;
109:       int     numCells = 0, numVertices = 0, numCorners = dim+1;

112:       ALE::PCICE::Builder::readConnectivity(comm, adjFilename, numCorners, useZeroBase, numCells, &cells);
113:       ALE::PCICE::Builder::readCoordinates(comm, coordFilename, dim, numVertices, &coordinates);
114:       ALE::New::SieveBuilder<sieve_type>::buildTopology(sieve, dim, numCells, cells, numVertices, interpolate, numCorners);
115:       sieve->stratify();
116:       topology->setPatch(0, sieve);
117:       topology->stratify();
118:       mesh->setTopology(topology);
119:       ALE::New::SieveBuilder<sieve_type>::buildCoordinates(mesh->getRealSection("coordinates"), dim, coordinates);
120:       if (cells) {PetscFree(cells);}
121:       if (coordinates) {PetscFree(coordinates);}
122:       return mesh;
123:     };
124:     // Creates boundary sections:
125:     //   IBC[NBFS,2]:     ALL
126:     //     BL[NBFS,1]:
127:     //     BNVEC[NBFS,2]:
128:     //   BCFUNC[NBCF,NV]: ALL
129:     //   IBNDFS[NBN,2]:   STILL NEED 4-5
130:     //     BNNV[NBN,2]
131:     void Builder::readBoundary(const Obj<Mesh>& mesh, const std::string& bcFilename, const int numBdFaces, const int numBdVertices) {
132:       const Mesh::topology_type::patch_type patch = 0;
133:       PetscViewer    viewer;
134:       FILE          *f;
135:       char           buf[2048];

138:       const Obj<Mesh::int_section_type>& ibc    = mesh->getIntSection("IBC");
139:       const Obj<Mesh::int_section_type>& ibndfs = mesh->getIntSection("IBNDFS");
140:       const Obj<Mesh::int_section_type>& ibcnum = mesh->getIntSection("IBCNUM");
141:       const Obj<Mesh::real_section_type>&    bl     = mesh->getRealSection("BL");
142:       const Obj<Mesh::real_section_type>&    bnvec  = mesh->getRealSection("BNVEC");
143:       const Obj<Mesh::real_section_type>&    bnnv   = mesh->getRealSection("BNNV");
144:       const Obj<Mesh::real_section_type>&    bcvec  = mesh->getRealSection("BCVEC");
145:       if (mesh->commRank() != 0) {
146:         mesh->distributeBCValues();
147:         return;
148:       }
149:       PetscViewerCreate(PETSC_COMM_SELF, &viewer);
150:       PetscViewerSetType(viewer, PETSC_VIEWER_ASCII);
151:       PetscViewerFileSetMode(viewer, FILE_MODE_READ);
152:       PetscViewerFileSetName(viewer, bcFilename.c_str());
153:       PetscViewerASCIIGetPointer(viewer, &f);
154:       // Create IBC section
155:       int *tmpIBC = new int[numBdFaces*4];
156:       std::map<int,std::set<int> > elem2Idx;
157:       std::map<int,int> bfReorder;
158:       for(int bf = 0; bf < numBdFaces; bf++) {
159:         const char *x = strtok(fgets(buf, 2048, f), " ");

161:         // Ignore boundary face number
162:         x = strtok(NULL, " ");
163:         tmpIBC[bf*4+0] = atoi(x);
164:         x = strtok(NULL, " ");
165:         tmpIBC[bf*4+1] = atoi(x);
166:         x = strtok(NULL, " ");
167:         tmpIBC[bf*4+2] = atoi(x);
168:         x = strtok(NULL, " ");
169:         tmpIBC[bf*4+3] = atoi(x);
170:         const int elem = tmpIBC[bf*4+0]-1;

172:         ibc->addFiberDimension(patch, elem, 4);
173:         ibcnum->addFiberDimension(patch, elem, 1);
174:         bl->addFiberDimension(patch, elem, 1);
175:         bnvec->addFiberDimension(patch, elem, 2);
176:         bcvec->addFiberDimension(patch, elem, 4);
177:         elem2Idx[elem].insert(bf);
178:       }
179:       ibc->allocate();
180:       ibcnum->allocate();
181:       bl->allocate();
182:       bnvec->allocate();
183:       bcvec->allocate();
184:       const Mesh::int_section_type::chart_type& chart = ibc->getPatch(patch);
185:       int num = 1;

187:       for(Mesh::int_section_type::chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
188:         const int elem = *p_iter;
189:         int bfNum[2];
190:         int k = 0;

192:         for(std::set<int>::const_iterator i_iter = elem2Idx[elem].begin(); i_iter != elem2Idx[elem].end(); ++i_iter) {
193:           bfReorder[(*i_iter)+1] = num;
194:           bfNum[k++] = num;
195:           num++;
196:         }
197:         ibcnum->update(patch, elem, bfNum);
198:       }
199:       for(int bf = 0; bf < numBdFaces; bf++) {
200:         const int elem = tmpIBC[bf*4]-1;

202:         if (elem2Idx[elem].size() > 1) {
203:           if (*elem2Idx[elem].begin() == bf) {
204:             int values[8];
205:             int k = 0;

207:             for(std::set<int>::const_iterator i_iter = elem2Idx[elem].begin(); i_iter != elem2Idx[elem].end(); ++i_iter) {
208:               for(int v = 0; v < 4; ++v) {
209:                 values[k*4+v] = tmpIBC[*i_iter*4+v];
210:               }
211:               k++;
212:             }
213:             ibc->update(patch, elem, values);
214:           }
215:         } else {
216:           ibc->update(patch, elem, &tmpIBC[bf*4]);
217:         }
218:       }
219:       delete [] tmpIBC;
220:       // Create BCFUNC section
221:       int numBcFunc = atoi(strtok(fgets(buf, 2048, f), " "));
222:       for(int bc = 0; bc < numBcFunc; bc++) {
223:         const char *x = strtok(fgets(buf, 2048, f), " ");
224:         Mesh::bc_value_type value;

226:         // Ignore function number
227:         x = strtok(NULL, " ");
228:         value.rho = atof(x);
229:         x = strtok(NULL, " ");
230:         value.u   = atof(x);
231:         x = strtok(NULL, " ");
232:         value.v   = atof(x);
233:         x = strtok(NULL, " ");
234:         value.p   = atof(x);
235:         mesh->setBCValue(bc+1, value);
236:       }
237:       mesh->distributeBCValues();
238:       // Create IBNDFS section
239:       const int numElements = mesh->getTopology()->heightStratum(patch, 0)->size();
240:       int      *tmpIBNDFS   = new int[numBdVertices*3];

242:       for(int bv = 0; bv < numBdVertices; bv++) {
243:         const char *x = strtok(fgets(buf, 2048, f), " ");

245:         // Ignore boundary node number
246:         x = strtok(NULL, " ");
247:         tmpIBNDFS[bv*3+0] = atoi(x);
248:         x = strtok(NULL, " ");
249:         tmpIBNDFS[bv*3+1] = atoi(x);
250:         x = strtok(NULL, " ");
251:         tmpIBNDFS[bv*3+2] = atoi(x);
252:         ibndfs->setFiberDimension(patch, tmpIBNDFS[bv*3+0]-1+numElements, 5);
253:       }
254:       ibndfs->allocate();
255:       for(int bv = 0; bv < numBdVertices; bv++) {
256:         int values[5];

258:         values[0] = tmpIBNDFS[bv*3+0];
259:         // Covert to new boundary face numbers
260:         values[1] = bfReorder[tmpIBNDFS[bv*3+1]];
261:         values[2] = bfReorder[tmpIBNDFS[bv*3+2]];
262:         values[3] = 0;
263:         values[4] = 0;
264:         ibndfs->update(patch, values[0]-1+numElements, values);
265:       }
266:       PetscViewerDestroy(viewer);
267:       // Create BNNV[NBN,2]
268:       const int dim = mesh->getDimension();

270:       for(int bv = 0; bv < numBdVertices; bv++) {
271:         bnnv->setFiberDimension(patch, tmpIBNDFS[bv*3+0]-1+numElements, dim);
272:       }
273:       bnnv->allocate();
274:       delete [] tmpIBNDFS;
275:     };
276:     void Builder::outputVerticesLocal(const Obj<Mesh>& mesh, int *numVertices, int *dim, double *coordinates[], const bool columnMajor) {
277:       const Mesh::real_section_type::patch_type            patch      = 0;
278:       const Obj<Mesh::real_section_type>&                  coordSec   = mesh->getRealSection("coordinates");
279:       if (!coordSec->hasPatch(patch)) {
280:         *numVertices = 0;
281:         *dim         = 0;
282:         *coordinates = NULL;
283:         return;
284:       }
285:       const Obj<Mesh::topology_type>&                 topology   = mesh->getTopology();
286:       const Obj<Mesh::topology_type::label_sequence>& vertices   = topology->depthStratum(patch, 0);
287:       const Obj<Mesh::numbering_type>&                vNumbering = mesh->getFactory()->getLocalNumbering(topology, patch, 0);
288:       int            size     = vertices->size();
289:       int            embedDim = coordSec->getFiberDimension(patch, *vertices->begin());
290:       double        *coords;

293:       PetscMalloc(vertices->size()*embedDim * sizeof(double), &coords);
294:       for(Mesh::topology_type::label_sequence::iterator v_iter = vertices->begin(); v_iter != vertices->end(); ++v_iter) {
295:         const Mesh::real_section_type::value_type *array = coordSec->restrict(patch, *v_iter);
296:         const int                             row   = vNumbering->getIndex(*v_iter);

298:         if (columnMajor) {
299:           for(int d = 0; d < embedDim; d++) {
300:             coords[d*size + row] = array[d];
301:           }
302:         } else {
303:           for(int d = 0; d < embedDim; d++) {
304:             coords[row*embedDim + d] = array[d];
305:           }
306:         }
307:       }
308:       *numVertices = size;
309:       *dim         = embedDim;
310:       *coordinates = coords;
311:     };
312:     void Builder::outputElementsLocal(const Obj<Mesh>& mesh, int *numElements, int *numCorners, int *vertices[], const bool columnMajor) {
313:       const Mesh::topology_type::patch_type           patch      = 0;
314:       const Obj<Mesh::topology_type>&                 topology   = mesh->getTopology();
315:       if (!topology->hasPatch(patch)) {
316:         *numElements = 0;
317:         *numCorners  = 0;
318:         *vertices    = NULL;
319:         return;
320:       }
321:       const Obj<Mesh::sieve_type>&                    sieve      = topology->getPatch(patch);
322:       const Obj<Mesh::topology_type::label_sequence>& elements   = topology->heightStratum(patch, 0);
323:       const Obj<Mesh::numbering_type>&                eNumbering = mesh->getFactory()->getLocalNumbering(topology, patch, topology->depth());
324:       const Obj<Mesh::numbering_type>&                vNumbering = mesh->getFactory()->getLocalNumbering(topology, patch, 0);
325:       int            size         = elements->size();
326:       //int            corners      = sieve->nCone(*elements->begin(), topology->depth())->size();
327:       int            corners      = sieve->cone(*elements->begin())->size();
328:       int           *v;

331:       PetscMalloc(elements->size()*corners * sizeof(int), &v);
332:       for(Mesh::topology_type::label_sequence::iterator e_iter = elements->begin(); e_iter != elements->end(); ++e_iter) {
333:         const Obj<Mesh::sieve_type::traits::coneSequence> cone  = sieve->cone(*e_iter);
334:         Mesh::sieve_type::traits::coneSequence::iterator  begin = cone->begin();
335:         Mesh::sieve_type::traits::coneSequence::iterator  end   = cone->end();

337:         const int row = eNumbering->getIndex(*e_iter);
338:         int       c   = -1;
339:         if (columnMajor) {
340:           for(Mesh::sieve_type::traits::coneSequence::iterator c_iter = begin; c_iter != end; ++c_iter) {
341:             v[(++c)*size + row] = vNumbering->getIndex(*c_iter)+1;
342:           }
343:         } else {
344:           for(Mesh::sieve_type::traits::coneSequence::iterator c_iter = begin; c_iter != end; ++c_iter) {
345:             v[row*corners + ++c] = vNumbering->getIndex(*c_iter)+1;
346:           }
347:         }
348:       }
349:       *numElements = size;
350:       *numCorners  = corners;
351:       *vertices    = v;
352:     };
355:     PetscErrorCode Viewer::writeVertices(const ALE::Obj<ALE::Mesh>& mesh, PetscViewer viewer) {
356:       ALE::Obj<ALE::Mesh::real_section_type> coordinates = mesh->getRealSection("coordinates");
357: #if 0
358:       ALE::Mesh::field_type::patch_type patch;
359:       const double  *array = coordinates->restrict(patch);
360:       int            numVertices;

364:       //FIX:
365:       if (vertexBundle->getGlobalOffsets()) {
366:         numVertices = vertexBundle->getGlobalOffsets()[mesh->commSize()];
367:       } else {
368:         numVertices = mesh->getTopology()->depthStratum(0)->size();
369:       }
370:       PetscViewerASCIIPrintf(viewer, "%D\n", numVertices);
371:       if (mesh->commRank() == 0) {
372:         int numLocalVertices = mesh->getTopology()->depthStratum(0)->size();
373:         int embedDim = coordinates->getFiberDimension(patch, *mesh->getTopology()->depthStratum(0)->begin());
374:         int vertexCount = 1;

376:         for(int v = 0; v < numLocalVertices; v++) {
377:           PetscViewerASCIIPrintf(viewer, "%7D   ", vertexCount++);
378:           for(int d = 0; d < embedDim; d++) {
379:             if (d > 0) {
380:               PetscViewerASCIIPrintf(viewer, " ");
381:             }
382:             PetscViewerASCIIPrintf(viewer, "% 12.5E", array[v*embedDim+d]);
383:           }
384:           PetscViewerASCIIPrintf(viewer, "\n");
385:         }
386:         for(int p = 1; p < mesh->commSize(); p++) {
387:           double    *remoteCoords;
388:           MPI_Status status;

390:           MPI_Recv(&numLocalVertices, 1, MPI_INT, p, 1, mesh->comm(), &status);
391:           PetscMalloc(numLocalVertices*embedDim * sizeof(double), &remoteCoords);
392:           MPI_Recv(remoteCoords, numLocalVertices*embedDim, MPI_DOUBLE, p, 1, mesh->comm(), &status);
393:           for(int v = 0; v < numLocalVertices; v++) {
394:             PetscViewerASCIIPrintf(viewer,"%7D   ", vertexCount++);
395:             for(int d = 0; d < embedDim; d++) {
396:               if (d > 0) {
397:                 PetscViewerASCIIPrintf(viewer, " ");
398:               }
399:               PetscViewerASCIIPrintf(viewer, "% 12.5E", remoteCoords[v*embedDim+d]);
400:             }
401:             PetscViewerASCIIPrintf(viewer, "\n");
402:           }
403:         }
404:       } else {
405:         ALE::Obj<ALE::Mesh::bundle_type>                           globalOrder = coordinates->getGlobalOrder();
406:         ALE::Obj<ALE::Mesh::bundle_type::order_type::coneSequence> cone        = globalOrder->getPatch(patch);
407:         const int *offsets = coordinates->getGlobalOffsets();
408:         int        embedDim = coordinates->getFiberDimension(patch, *mesh->getTopology()->depthStratum(0)->begin());
409:         int        numLocalVertices = (offsets[mesh->commRank()+1] - offsets[mesh->commRank()])/embedDim;
410:         double    *localCoords;
411:         int        k = 0;

413:         PetscMalloc(numLocalVertices*embedDim * sizeof(double), &localCoords);
414:         for(ALE::Mesh::bundle_type::order_type::coneSequence::iterator p_iter = cone->begin(); p_iter != cone->end(); ++p_iter) {
415:           int dim = globalOrder->getFiberDimension(patch, *p_iter);

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

420:             for(int i = offset; i < offset+dim; ++i) {
421:               localCoords[k++] = array[i];
422:             }
423:           }
424:         }
425:         if (k != numLocalVertices*embedDim) {
426:           SETERRQ2(PETSC_ERR_PLIB, "Invalid number of coordinates to send %d should be %d", k, numLocalVertices*embedDim);
427:         }
428:         MPI_Send(&numLocalVertices, 1, MPI_INT, 0, 1, mesh->comm());
429:         MPI_Send(localCoords, numLocalVertices*embedDim, MPI_DOUBLE, 0, 1, mesh->comm());
430:         PetscFree(localCoords);
431:       }
432: #endif
433:       return(0);
434:     };
437:     PetscErrorCode Viewer::writeElements(const ALE::Obj<ALE::Mesh>& mesh, PetscViewer viewer) {
438:       ALE::Obj<ALE::Mesh::topology_type> topology = mesh->getTopology();
439: #if 0
440:       ALE::Obj<ALE::Mesh::sieve_type::traits::heightSequence> elements = topology->heightStratum(0);
441:       ALE::Obj<ALE::Mesh::bundle_type> elementBundle = mesh->getBundle(topology->depth());
442:       ALE::Obj<ALE::Mesh::bundle_type> vertexBundle = mesh->getBundle(0);
443:       ALE::Obj<ALE::Mesh::bundle_type> globalVertex = vertexBundle->getGlobalOrder();
444:       ALE::Obj<ALE::Mesh::bundle_type> globalElement = elementBundle->getGlobalOrder();
445:       ALE::Mesh::bundle_type::patch_type patch;
446:       std::string    orderName("element");
447:       int            dim  = mesh->getDimension();
448:       int            corners = topology->nCone(*elements->begin(), topology->depth())->size();
449:       int            numElements;

453:       if (corners != dim+1) {
454:         SETERRQ(PETSC_ERR_SUP, "PCICE only supports simplicies");
455:       }
456:       if (!globalVertex) {
457:         globalVertex = vertexBundle;
458:       }
459:       if (elementBundle->getGlobalOffsets()) {
460:         numElements = elementBundle->getGlobalOffsets()[mesh->commSize()];
461:       } else {
462:         numElements = mesh->getTopology()->heightStratum(0)->size();
463:       }
464:       if (mesh->commRank() == 0) {
465:         int elementCount = 1;

467:         PetscViewerASCIIPrintf(viewer, "%d\n", numElements);
468:         for(ALE::Mesh::sieve_type::traits::heightSequence::iterator e_itor = elements->begin(); e_itor != elements->end(); ++e_itor) {
469:           ALE::Obj<ALE::Mesh::bundle_type::order_type::coneSequence> cone = vertexBundle->getPatch(orderName, *e_itor);

471:           PetscViewerASCIIPrintf(viewer, "%7d", elementCount++);
472:           for(ALE::Mesh::bundle_type::order_type::coneSequence::iterator c_itor = cone->begin(); c_itor != cone->end(); ++c_itor) {
473:             PetscViewerASCIIPrintf(viewer, " %7d", globalVertex->getIndex(patch, *c_itor).prefix);
474:           }
475:           PetscViewerASCIIPrintf(viewer, "\n");
476:         }
477:         for(int p = 1; p < mesh->commSize(); p++) {
478:           int        numLocalElements;
479:           int       *remoteVertices;
480:           MPI_Status status;

482:           MPI_Recv(&numLocalElements, 1, MPI_INT, p, 1, mesh->comm(), &status);
483:           PetscMalloc(numLocalElements*corners * sizeof(int), &remoteVertices);
484:           MPI_Recv(remoteVertices, numLocalElements*corners, MPI_INT, p, 1, mesh->comm(), &status);
485:           for(int e = 0; e < numLocalElements; e++) {
486:             PetscViewerASCIIPrintf(viewer, "%7d", elementCount++);
487:             for(int c = 0; c < corners; c++) {
488:               PetscViewerASCIIPrintf(viewer, " %7d", remoteVertices[e*corners+c]);
489:             }
490:             PetscViewerASCIIPrintf(viewer, "\n");
491:           }
492:           PetscFree(remoteVertices);
493:         }
494:       } else {
495:         const int *offsets = elementBundle->getGlobalOffsets();
496:         int        numLocalElements = offsets[mesh->commRank()+1] - offsets[mesh->commRank()];
497:         int       *localVertices;
498:         int        k = 0;

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

504:           if (globalElement->getFiberDimension(patch, *e_itor) > 0) {
505:             for(ALE::Mesh::bundle_type::order_type::coneSequence::iterator c_itor = cone->begin(); c_itor != cone->end(); ++c_itor) {
506:               localVertices[k++] = globalVertex->getIndex(patch, *c_itor).prefix;
507:             }
508:           }
509:         }
510:         if (k != numLocalElements*corners) {
511:           SETERRQ2(PETSC_ERR_PLIB, "Invalid number of vertices to send %d should be %d", k, numLocalElements*corners);
512:         }
513:         MPI_Send(&numLocalElements, 1, MPI_INT, 0, 1, mesh->comm());
514:         MPI_Send(localVertices, numLocalElements*corners, MPI_INT, 0, 1, mesh->comm());
515:         PetscFree(localVertices);
516:       }
517: #endif
518:       return(0);
519:     };
522:     PetscErrorCode Viewer::writeVerticesLocal(const Obj<Mesh>& mesh, PetscViewer viewer) {
523:       const Mesh::real_section_type::patch_type            patch       = 0;
524:       Obj<Mesh::real_section_type>                         coordinates = mesh->getRealSection("coordinates");
525:       const Obj<Mesh::topology_type>&                 topology    = mesh->getTopology();
526:       const Obj<Mesh::topology_type::label_sequence>& vertices    = topology->depthStratum(patch, 0);
527:       const Obj<Mesh::numbering_type>&                vNumbering  = mesh->getFactory()->getLocalNumbering(topology, patch, 0);
528:       int            embedDim = coordinates->getFiberDimension(patch, *vertices->begin());

532:       PetscViewerASCIIPrintf(viewer, "%D\n", vertices->size());
533:       for(Mesh::topology_type::label_sequence::iterator v_iter = vertices->begin(); v_iter != vertices->end(); ++v_iter) {
534:         const Mesh::real_section_type::value_type *array = coordinates->restrict(patch, *v_iter);

536:         PetscViewerASCIIPrintf(viewer, "%7D   ", vNumbering->getIndex(*v_iter)+1);
537:         for(int d = 0; d < embedDim; d++) {
538:           if (d > 0) {
539:             PetscViewerASCIIPrintf(viewer, " ");
540:           }
541:           PetscViewerASCIIPrintf(viewer, "% 12.5E", array[d]);
542:         }
543:         PetscViewerASCIIPrintf(viewer, "\n");
544:       }
545:       return(0);
546:     };
549:     PetscErrorCode Viewer::writeRestart(const Obj<Mesh>& mesh, PetscViewer viewer) {
550:       const Mesh::real_section_type::patch_type patch = 0;
551:       const Obj<Mesh::real_section_type>&   velocity    = mesh->getRealSection("VELN");
552:       const Obj<Mesh::real_section_type>&   pressure    = mesh->getRealSection("PN");
553:       const Obj<Mesh::real_section_type>&   temperature = mesh->getRealSection("TN");
554:       const Obj<Mesh::topology_type>&  topology    = mesh->getTopology();
555:       const Obj<Mesh::numbering_type>& cNumbering  = mesh->getFactory()->getNumbering(topology, patch, topology->depth());
556:       const Obj<Mesh::numbering_type>& vNumbering  = mesh->getFactory()->getNumbering(topology, patch, 0);
557:       const int                        numCells    = cNumbering->getGlobalSize();

561:       int          blen[2];
562:       MPI_Aint     indices[2];
563:       MPI_Datatype oldtypes[2], newtype;
564:       blen[0] = 1; indices[0] = 0;           oldtypes[0] = MPI_INT;
565:       blen[1] = 4; indices[1] = sizeof(int); oldtypes[1] = MPI_DOUBLE;
566:       MPI_Type_struct(2, blen, indices, oldtypes, &newtype);
567:       MPI_Type_commit(&newtype);

569:       if (mesh->commRank() == 0) {
570:         const Obj<Mesh::topology_type::label_sequence>& vertices = mesh->getTopology()->depthStratum(patch, 0);

572:         for(Mesh::topology_type::label_sequence::iterator v_iter = vertices->begin(); v_iter != vertices->end(); ++v_iter) {
573:           if (vNumbering->isLocal(*v_iter)) {
574:             const ALE::Mesh::real_section_type::value_type *veln = velocity->restrictPoint(patch, *v_iter);
575:             const ALE::Mesh::real_section_type::value_type *pn   = pressure->restrictPoint(patch, *v_iter);
576:             const ALE::Mesh::real_section_type::value_type *tn   = temperature->restrictPoint(patch, *v_iter);

578:             PetscViewerASCIIPrintf(viewer, "%6d% 16.8E% 16.8E% 16.8E% 16.8E\n", *v_iter-numCells+1, veln[0], veln[1], pn[0], tn[0]);
579:           }
580:         }
581:         for(int p = 1; p < mesh->commSize(); p++) {
582:           RestartType *remoteValues;
583:           int          numLocalElements;
584:           MPI_Status   status;

586:           MPI_Recv(&numLocalElements, 1, MPI_INT, p, 1, mesh->comm(), &status);
587:           PetscMalloc(numLocalElements * sizeof(RestartType), &remoteValues);
588:           MPI_Recv(remoteValues, numLocalElements, newtype, p, 1, mesh->comm(), &status);
589:           for(int e = 0; e < numLocalElements; e++) {
590:             PetscViewerASCIIPrintf(viewer, "%6d% 16.8E% 16.8E% 16.8E% 16.8E\n", remoteValues[e].vertex-numCells+1, remoteValues[e].veln_x, remoteValues[e].veln_y, remoteValues[e].pn, remoteValues[e].tn);
591:           }
592:         }
593:       } else {
594:         const Obj<Mesh::topology_type::label_sequence>& vertices = mesh->getTopology()->depthStratum(patch, 0);
595:         RestartType *localValues;
596:         int numLocalElements = vNumbering->getLocalSize();
597:         int k = 0;

599:         PetscMalloc(numLocalElements * sizeof(RestartType), &localValues);
600:         for(Mesh::topology_type::label_sequence::iterator v_iter = vertices->begin(); v_iter != vertices->end(); ++v_iter) {
601:           if (vNumbering->isLocal(*v_iter)) {
602:             const ALE::Mesh::real_section_type::value_type *veln = velocity->restrictPoint(patch, *v_iter);
603:             const ALE::Mesh::real_section_type::value_type *pn   = pressure->restrictPoint(patch, *v_iter);
604:             const ALE::Mesh::real_section_type::value_type *tn   = temperature->restrictPoint(patch, *v_iter);

606:             localValues[k].vertex = *v_iter;
607:             localValues[k].veln_x = veln[0];
608:             localValues[k].veln_y = veln[1];
609:             localValues[k].pn     = pn[0];
610:             localValues[k].tn     = tn[0];
611:             k++;
612:           }
613:         }
614:         if (k != numLocalElements) {
615:           SETERRQ2(PETSC_ERR_PLIB, "Invalid number of values to send for field, %d should be %d", k, numLocalElements);
616:         }
617:         MPI_Send(&numLocalElements, 1, MPI_INT, 0, 1, mesh->comm());
618:         MPI_Send(localValues, numLocalElements, newtype, 0, 1, mesh->comm());
619:         PetscFree(localValues);
620:       }
621:       MPI_Type_free(&newtype);
622:       return(0);
623:     };
624:   };
625: };