Actual source code: meshpcice.c

petsc-3.4.4 2014-03-13
  1: #include <petscdmmesh_formats.hh>   /*I      "petscmesh.h"   I*/

  5: PetscErrorCode WritePCICEVertices(DM dm, PetscViewer viewer)
  6: {
  7:   ALE::Obj<PETSC_MESH_TYPE> m;
  8:   PetscErrorCode            ierr;

 10:   DMMeshGetMesh(dm, m);
 11:   return ALE::PCICE::Viewer::writeVertices(m, viewer);
 12: }

 16: PetscErrorCode WritePCICEElements(DM dm, PetscViewer viewer)
 17: {
 18:   ALE::Obj<PETSC_MESH_TYPE> m;
 19:   PetscErrorCode            ierr;

 21:   DMMeshGetMesh(dm, m);
 22:   return ALE::PCICE::Viewer::writeElements(m, viewer);
 23: }

 27: PetscErrorCode WritePCICERestart(DM dm, PetscViewer viewer)
 28: {
 29:   ALE::Obj<PETSC_MESH_TYPE> m;
 30:   PetscErrorCode            ierr;

 32:   DMMeshGetMesh(dm, m);
 33:   return ALE::PCICE::Viewer::writeRestart(m, viewer);
 34: }

 38: /*@C
 39:   DMMeshCreatePCICE - Create a DMMesh from PCICE files.

 41:   Not Collective

 43:   Input Parameters:
 44: + dim - The topological mesh dimension
 45: . coordFilename - The file containing vertex coordinates
 46: . adjFilename - The file containing the vertices for each element
 47: . interpolate - The flag for construction of intermediate elements
 48: . bcFilename - The file containing the boundary topology and conditions
 49: . numBdFaces - The number of boundary faces (or edges)
 50: - numBdVertices - The number of boundary vertices

 52:   Output Parameter:
 53: . mesh - The DMMesh object

 55:   Level: beginner

 57: .keywords: mesh, PCICE
 58: .seealso: DMMeshCreate()
 59: @*/
 60: PetscErrorCode DMMeshCreatePCICE(MPI_Comm comm, const int dim, const char coordFilename[], const char adjFilename[], PetscBool interpolate, const char bcFilename[], DM *mesh)
 61: {
 62:   ALE::Obj<PETSC_MESH_TYPE> m;
 63:   PetscInt                  debug = 0;
 64:   PetscBool                 flag;
 65:   PetscErrorCode            ierr;

 68:   DMMeshCreate(comm, mesh);
 69:   PetscOptionsGetInt(NULL, "-debug", &debug, &flag);
 70:   try {
 71:     m = ALE::PCICE::Builder::readMesh(comm, dim, std::string(coordFilename), std::string(adjFilename), true, interpolate, debug);
 72:     if (debug) m->view("Mesh");
 73:   } catch(ALE::Exception e) {
 74:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_OPEN, e.message());
 75:   }
 76:   if (bcFilename) ALE::PCICE::Builder::readBoundary(m, std::string(bcFilename));
 77:   DMMeshSetMesh(*mesh, m);
 78:   return(0);
 79: }

 83: /*@C
 84:   PCICERenumberBoundary - Change global element names into offsets

 86:   Collective on DMMesh

 88:   Input Parameters:
 89: . mesh - the mesh

 91:   Level: advanced

 93:   .seealso: DMMeshCreate()
 94: @*/
 95: PetscErrorCode  PCICERenumberBoundary(DM dm)
 96: {
 97:   ALE::Obj<PETSC_MESH_TYPE> m;
 98:   PetscErrorCode            ierr;

101:   DMMeshGetMesh(dm, m);
102:   try {
103:     ALE::PCICE::fuseBoundary(m);
104:   } catch(ALE::Exception e) {
105:     SETERRQ(PETSC_COMM_SELF,100, e.msg().c_str());
106:   }
107:   return(0);
108: }

112: /*@C
113:   BCSectionGetArray - Returns the array underlying the BCSection.

115:   Not Collective

117:   Input Parameters:
118: + mesh - The DMMesh object
119: - name - The section name

121:   Output Parameters:
122: + numElements - The number of mesh element with values
123: . fiberDim - The number of values per element
124: - array - The array

126:   Level: intermediate

128: .keywords: mesh, elements
129: .seealso: DMMeshCreate()
130: @*/
131: PetscErrorCode BCSectionGetArray(DM dm, const char name[], PetscInt *numElements, PetscInt *fiberDim, PetscInt *array[])
132: {
133:   ALE::Obj<PETSC_MESH_TYPE> m;
134:   PetscErrorCode            ierr;

137:   DMMeshGetMesh(dm, m);
138:   const ALE::Obj<PETSC_MESH_TYPE::int_section_type>& section = m->getIntSection(std::string(name));
139:   if (!section->size()) {
140:     *numElements = 0;
141:     *fiberDim    = 0;
142:     *array       = NULL;
143:     return(0);
144:   }
145:   const PETSC_MESH_TYPE::int_section_type::chart_type& chart = section->getChart();

147:   int fiberDimMin = section->getFiberDimension(*chart.begin());
148:   int numElem     = 0;

150:   for (PETSC_MESH_TYPE::int_section_type::chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
151:     const int fiberDim = section->getFiberDimension(*c_iter);

153:     if (fiberDim < fiberDimMin) fiberDimMin = fiberDim;
154:   }
155:   for (PETSC_MESH_TYPE::int_section_type::chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
156:     const int fiberDim = section->getFiberDimension(*c_iter);

158:     numElem += fiberDim/fiberDimMin;
159:   }
160:   *numElements = numElem;
161:   *fiberDim    = fiberDimMin;
162:   *array       = (PetscInt*) section->restrictSpace();
163:   return(0);
164: }

168: /*@C
169:   BCSectionRealCreate - Creates a BCSection.

171:   Not Collective

173:   Input Parameters:
174: + mesh - The DMMesh object
175: . name - The section name
176: - fiberDim - The number of values per element

178:   Level: intermediate

180: .keywords: mesh, elements
181: .seealso: DMMeshCreate()
182: @*/
183: PetscErrorCode BCSectionRealCreate(DM dm, const char name[], PetscInt fiberDim)
184: {
185:   ALE::Obj<PETSC_MESH_TYPE> m;
186:   PetscErrorCode            ierr;

189:   DMMeshGetMesh(dm, m);
190:   const ALE::Obj<PETSC_MESH_TYPE::real_section_type>&  section = m->getRealSection(std::string(name));
191:   const ALE::Obj<PETSC_MESH_TYPE::int_section_type>&   ibc     = m->getIntSection("IBC");
192:   const PETSC_MESH_TYPE::int_section_type::chart_type& chart   = ibc->getChart();

194:   for (PETSC_MESH_TYPE::int_section_type::chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
195:     section->setFiberDimension(*p_iter, ibc->getFiberDimension(*p_iter));
196:   }
197:   m->allocate(section);
198:   return(0);
199: }

203: /*@C
204:   BCSectionRealGetArray - Returns the array underlying the BCSection.

206:   Not Collective

208:   Input Parameters:
209: + mesh - The DMMesh object
210: - name - The section name

212:   Output Parameters:
213: + numElements - The number of mesh element with values
214: . fiberDim - The number of values per element
215: - array - The array

217:   Level: intermediate

219: .keywords: mesh, elements
220: .seealso: DMMeshCreate()
221: @*/
222: PetscErrorCode BCSectionRealGetArray(DM dm, const char name[], PetscInt *numElements, PetscInt *fiberDim, PetscReal *array[])
223: {
224:   ALE::Obj<PETSC_MESH_TYPE> m;
225:   PetscErrorCode            ierr;

228:   DMMeshGetMesh(dm, m);
229:   const ALE::Obj<PETSC_MESH_TYPE::real_section_type>& section = m->getRealSection(std::string(name));
230:   if (!section->size()) {
231:     *numElements = 0;
232:     *fiberDim    = 0;
233:     *array       = NULL;
234:     return(0);
235:   }
236:   const PETSC_MESH_TYPE::real_section_type::chart_type& chart = section->getChart();

238:   int fiberDimMin = section->getFiberDimension(*chart.begin());
239:   int numElem     = 0;

241:   for (PETSC_MESH_TYPE::real_section_type::chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
242:     const int fiberDim = section->getFiberDimension(*c_iter);

244:     if (fiberDim < fiberDimMin) fiberDimMin = fiberDim;
245:   }
246:   for (PETSC_MESH_TYPE::real_section_type::chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
247:     const int fiberDim = section->getFiberDimension(*c_iter);

249:     numElem += fiberDim/fiberDimMin;
250:   }
251:   *numElements = numElem;
252:   *fiberDim    = fiberDimMin;
253:   *array       = (PetscReal*) section->restrictSpace();
254:   return(0);
255: }

259: PetscErrorCode BCFUNCGetArray(DM dm, PetscInt *numElements, PetscInt *fiberDim, PetscScalar *array[])
260: {
261:   ALE::Obj<PETSC_MESH_TYPE> m;
262:   PetscErrorCode            ierr;

265:   DMMeshGetMesh(dm, m);
266: #if 0
267:   PETSC_MESH_TYPE::bc_values_type& bcValues = m->getBCValues();
268:   *numElements = bcValues.size();
269:   *fiberDim    = 4;
270:   *array       = new PetscScalar[(*numElements)*(*fiberDim)];
271:   for (int bcf = 1; bcf <= (int) bcValues.size(); ++bcf) {
272:     (*array)[(bcf-1)*4+0] = bcValues[bcf].rho;
273:     (*array)[(bcf-1)*4+1] = bcValues[bcf].u;
274:     (*array)[(bcf-1)*4+2] = bcValues[bcf].v;
275:     (*array)[(bcf-1)*4+3] = bcValues[bcf].p;
276:   }
277: #else
278:   *numElements = 0;
279:   *fiberDim    = 0;
280:   *array       = NULL;
281: #endif
282:   return(0);
283: }

285: namespace ALE {
286: namespace PCICE {
287: //
288: // Builder methods
289: //
290: void Builder::readConnectivity(MPI_Comm comm, const std::string& filename, int& corners, const bool useZeroBase, int& numElements, int *vertices[])
291: {
292:   PetscViewer    viewer;
293:   FILE           *f;
294:   PetscInt       numCells, cellCount = 0;
295:   PetscInt       *verts;
296:   char           buf[2048];
297:   PetscInt       c;
298:   PetscInt       commRank;

301:   MPI_Comm_rank(comm, &commRank);

303:   if (commRank != 0) return;
304:   PetscViewerCreate(PETSC_COMM_SELF, &viewer);
305:   PetscViewerSetType(viewer, PETSCVIEWERASCII);
306:   PetscViewerFileSetMode(viewer, FILE_MODE_READ);
307:   PetscViewerFileSetName(viewer, filename.c_str());
308:   if (ierr) {
309:     ostringstream txt;
310:     txt << "Could not open PCICE connectivity file: " << filename;
311:     throw ALE::Exception(txt.str().c_str());
312:   }
313:   PetscViewerASCIIGetPointer(viewer, &f);
314:   if (!fgets(buf, 2048, f)) {
315:     throw ALE::Exception("Invalid connectivity file: Missing number of elements");
316:   }
317:   const char *sizes = strtok(buf, " ");
318:   numCells = atoi(sizes);
319:   sizes    = strtok(NULL, " ");
320:   if (sizes != NULL) {
321:     corners = atoi(sizes);
322:     std::cout << "Reset corners to " << corners << std::endl;
323:   }
324:   PetscMalloc(numCells*corners * sizeof(PetscInt), &verts);
325:   while (fgets(buf, 2048, f) != NULL) {
326:     const char *v = strtok(buf, " ");

328:     /* Ignore cell number */
329:     v = strtok(NULL, " ");
330:     for (c = 0; c < corners; c++) {
331:       int vertex = atoi(v);

333:       if (!useZeroBase) vertex -= 1;
334:       verts[cellCount*corners+c] = vertex;

336:       v = strtok(NULL, " ");
337:     }
338:     cellCount++;
339:   }
340:   PetscViewerDestroy(&viewer);
341:   numElements = numCells;
342:   *vertices   = verts;
343: };
344: void Builder::readCoordinates(MPI_Comm comm, const std::string& filename, const int dim, int& numVertices, PetscReal *coordinates[])
345: {
346:   PetscViewer    viewer;
347:   FILE           *f;
348:   PetscInt       numVerts, vertexCount = 0;
349:   PetscReal      *coords;
350:   char           buf[2048];
351:   PetscInt       c;
352:   PetscInt       commRank;

355:   MPI_Comm_rank(comm, &commRank);

357:   if (commRank != 0) return;
358:   PetscViewerCreate(PETSC_COMM_SELF, &viewer);
359:   PetscViewerSetType(viewer, PETSCVIEWERASCII);
360:   PetscViewerFileSetMode(viewer, FILE_MODE_READ);
361:   PetscViewerFileSetName(viewer, filename.c_str());
362:   if (ierr) {
363:     ostringstream txt;
364:     txt << "Could not open PCICE coordinate file: " << filename;
365:     throw ALE::Exception(txt.str().c_str());
366:   }
367:   PetscViewerASCIIGetPointer(viewer, &f);
368:   numVerts = atoi(fgets(buf, 2048, f));
369:   PetscMalloc(numVerts*dim * sizeof(PetscReal), &coords);
370:   while (fgets(buf, 2048, f) != NULL) {
371:     const char *x = strtok(buf, " ");

373:     /* Ignore vertex number */
374:     x = strtok(NULL, " ");
375:     for (c = 0; c < dim; c++) {
376:       coords[vertexCount*dim+c] = atof(x);

378:       x = strtok(NULL, " ");
379:     }
380:     vertexCount++;
381:   }
382:   PetscViewerDestroy(&viewer);
383:   numVertices  = numVerts;
384:   *coordinates = coords;
385: };
386: Obj<PETSC_MESH_TYPE> Builder::readMesh(MPI_Comm comm, const int dim, const std::string& basename, const bool useZeroBase = true, const bool interpolate = true, const int debug = 0)
387: {
388:   return readMesh(comm, dim, basename+".nodes", basename+".lcon", useZeroBase, interpolate, debug);
389: };
390: #if defined(PETSC_OPT_SIEVE)
391: Obj<PETSC_MESH_TYPE> 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)
392: {
393:   typedef ALE::Mesh<PetscInt,PetscScalar> FlexMesh;
394:   Obj<Mesh>                       mesh         = new Mesh(comm, dim, debug);
395:   Obj<sieve_type>                 sieve        = new sieve_type(comm, debug);
396:   const Obj<FlexMesh>             m            = new FlexMesh(comm, dim, debug);
397:   const Obj<FlexMesh::sieve_type> s            = new FlexMesh::sieve_type(comm, debug);
398:   int                             *cells       = NULL;
399:   PetscReal                       *coordinates = NULL;
400:   int                             numCells     = 0, numVertices = 0, numCorners = dim+1;
401:   PetscErrorCode                  ierr;

403:   ALE::PCICE::Builder::readConnectivity(comm, adjFilename, numCorners, useZeroBase, numCells, &cells);
404:   ALE::PCICE::Builder::readCoordinates(comm, coordFilename, dim, numVertices, &coordinates);
405:   ALE::SieveBuilder<FlexMesh>::buildTopology(s, dim, numCells, cells, numVertices, interpolate, numCorners, -1, m->getArrowSection("orientation"));
406:   m->setSieve(s);
407:   m->stratify();
408:   mesh->setSieve(sieve);
409:   std::map<Mesh::point_type,Mesh::point_type> renumbering;
410:   ALE::ISieveConverter::convertSieve(*s, *sieve, renumbering, false);
411:   mesh->stratify();
412:   ALE::ISieveConverter::convertOrientation(*s, *sieve, renumbering, m->getArrowSection("orientation").ptr());
413:   ALE::SieveBuilder<PETSC_MESH_TYPE>::buildCoordinates(mesh, dim, coordinates);
414:   if (cells) {PetscFree(cells);CHKERRXX(ierr);}
415:   if (coordinates) {PetscFree(coordinates);CHKERRXX(ierr);}
416:   return mesh;
417: };
418: void Builder::readBoundary(const Obj<Mesh>& mesh, const std::string& bcFilename)
419: {
420:   throw ALE::Exception("Not implemented for optimized sieves");
421: };
422: void Builder::outputVerticesLocal(const Obj<Mesh>& mesh, int *numVertices, int *dim, PetscReal *coordinates[], const bool columnMajor)
423: {
424:   const Obj<Mesh::real_section_type>& coordSec = mesh->getRealSection("coordinates");
425:   if (!coordSec->size()) {
426:     *numVertices = 0;
427:     *dim         = 0;
428:     *coordinates = NULL;
429:     return;
430:   }
431:   const Obj<Mesh::label_sequence>& vertices   = mesh->depthStratum(0);
432:   const Obj<Mesh::numbering_type>& vNumbering = mesh->getFactory()->getLocalNumbering(mesh, 0);
433:   int                              size       = vertices->size();
434:   int                              embedDim   = coordSec->getFiberDimension(*vertices->begin());
435:   PetscReal                        *coords;
436:   PetscErrorCode                   ierr;

438:   PetscMalloc(vertices->size()*embedDim * sizeof(double), &coords);CHKERRXX(ierr);
439:   for (Mesh::label_sequence::iterator v_iter = vertices->begin(); v_iter != vertices->end(); ++v_iter) {
440:     const Mesh::real_section_type::value_type *array = coordSec->restrictPoint(*v_iter);
441:     const int                                 row    = vNumbering->getIndex(*v_iter);

443:     if (columnMajor) {
444:       for (int d = 0; d < embedDim; d++) {
445:         coords[d*size + row] = array[d];
446:       }
447:     } else {
448:       for (int d = 0; d < embedDim; d++) {
449:         coords[row*embedDim + d] = array[d];
450:       }
451:     }
452:   }
453:   *numVertices = size;
454:   *dim         = embedDim;
455:   *coordinates = coords;
456: };
457: void Builder::outputElementsLocal(const Obj<Mesh>& mesh, int *numElements, int *numCorners, int *vertices[], const bool columnMajor)
458: {
459:   if (!mesh->heightStratum(0)->size()) {
460:     *numElements = 0;
461:     *numCorners  = 0;
462:     *vertices    = NULL;
463:     return;
464:   }
465:   const Obj<Mesh::sieve_type>&     sieve      = mesh->getSieve();
466:   const Obj<Mesh::label_sequence>& elements   = mesh->heightStratum(0);
467:   const Obj<Mesh::numbering_type>& eNumbering = mesh->getFactory()->getLocalNumbering(mesh, mesh->depth());
468:   const Obj<Mesh::numbering_type>& vNumbering = mesh->getFactory()->getLocalNumbering(mesh, 0);
469:   int                              size       = elements->size();
470:   int                              corners    = sieve->getConeSize(*elements->begin());
471:   int                              *v;
472:   PetscErrorCode                   ierr;

474:   PetscMalloc(size*corners * sizeof(int), &v);CHKERRXX(ierr);
475:   for (Mesh::label_sequence::iterator e_iter = elements->begin(); e_iter != elements->end(); ++e_iter) {
476:     const Obj<Mesh::sieve_type::coneSequence>      cone  = sieve->cone(*e_iter);
477:     Mesh::sieve_type::coneSequence::const_iterator begin = cone->begin();
478:     Mesh::sieve_type::coneSequence::const_iterator end   = cone->end();

480:     const int row = eNumbering->getIndex(*e_iter);
481:     int       c   = -1;
482:     if (columnMajor) {
483:       for (Mesh::sieve_type::coneSequence::iterator c_iter = begin; c_iter != end; ++c_iter) {
484:         v[(++c)*size + row] = vNumbering->getIndex(*c_iter)+1;
485:       }
486:     } else {
487:       for (Mesh::sieve_type::coneSequence::iterator c_iter = begin; c_iter != end; ++c_iter) {
488:         v[row*corners + ++c] = vNumbering->getIndex(*c_iter)+1;
489:       }
490:     }
491:   }
492:   *numElements = size;
493:   *numCorners  = corners;
494:   *vertices    = v;
495: };
496: PetscErrorCode Viewer::writeVertices(const ALE::Obj<Mesh>& mesh, PetscViewer viewer)
497: {
498:   throw ALE::Exception("Not implemented for optimized sieves");
499: };
500: PetscErrorCode Viewer::writeElements(const ALE::Obj<Mesh>& mesh, PetscViewer viewer)
501: {
502:   throw ALE::Exception("Not implemented for optimized sieves");
503: };
504: PetscErrorCode Viewer::writeVerticesLocal(const Obj<Mesh>& mesh, PetscViewer viewer)
505: {
506:   throw ALE::Exception("Not implemented for optimized sieves");
507: };
508: PetscErrorCode Viewer::writeRestart(const Obj<Mesh>& mesh, PetscViewer viewer)
509: {
510:   throw ALE::Exception("Not implemented for optimized sieves");
511: };
512: void fuseBoundary(const ALE::Obj<PETSC_MESH_TYPE>& mesh)
513: {
514:   throw ALE::Exception("Not implemented for optimized sieves");
515: };
516: #else
517: Obj<PETSC_MESH_TYPE> 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)
518: {
519:   Obj<Mesh>       mesh         = new DMMesh(comm, dim, debug);
520:   Obj<sieve_type> sieve        = new sieve_type(comm, debug);
521:   int             *cells       = NULL;
522:   double          *coordinates = NULL;
523:   int             numCells     = 0, numVertices = 0, numCorners = dim+1;
524:   PetscErrorCode  ierr;

526:   ALE::PCICE::Builder::readConnectivity(comm, adjFilename, numCorners, useZeroBase, numCells, &cells);
527:   ALE::PCICE::Builder::readCoordinates(comm, coordFilename, dim, numVertices, &coordinates);
528:   ALE::SieveBuilder<PETSC_MESH_TYPE>::buildTopology(sieve, dim, numCells, cells, numVertices, interpolate, numCorners, -1, mesh->getArrowSection("orientation"));
529:   mesh->setSieve(sieve);
530:   mesh->stratify();
531:   ALE::SieveBuilder<PETSC_MESH_TYPE>::buildCoordinates(mesh, dim, coordinates);
532:   if (cells) {PetscFree(cells);CHKERRXX(ierr);}
533:   if (coordinates) {PetscFree(coordinates);CHKERRXX(ierr);}
534:   return mesh;
535: };
536: // Creates boundary sections:
537: //   IBC[NBFS,2]:     ALL
538: //     BL[NBFS,1]:
539: //     BNVEC[NBFS,2]:
540: //   BCFUNC[NBCF,NV]: ALL
541: //   IBNDFS[NBN,2]:   STILL NEED 4-5
542: //     BNNV[NBN,2]
543: void Builder::readBoundary(const Obj<Mesh>& mesh, const std::string& bcFilename)
544: {
545:   PetscViewer    viewer;
546:   FILE           *f;
547:   char           buf[2048];

550:   const Obj<Mesh::int_section_type>&  ibc    = mesh->getIntSection("IBC");
551:   const Obj<Mesh::int_section_type>&  ibndfs = mesh->getIntSection("IBNDFS");
552:   const Obj<Mesh::int_section_type>&  ibcnum = mesh->getIntSection("IBCNUM");
553:   const Obj<Mesh::int_section_type>&  ibfcon = mesh->getIntSection("IBFCON");
554:   const Obj<Mesh::real_section_type>& bl     = mesh->getRealSection("BL");
555:   const Obj<Mesh::real_section_type>& bnvec  = mesh->getRealSection("BNVEC");
556:   const Obj<Mesh::real_section_type>& bnnv   = mesh->getRealSection("BNNV");
557:   if (mesh->commRank() != 0) {
558: #if 0
559:     mesh->distributeBCValues();
560: #endif
561:     return;
562:   }
563:   PetscViewerCreate(PETSC_COMM_SELF, &viewer);CHKERRXX(ierr);
564:   PetscViewerSetType(viewer, PETSCVIEWERASCII);CHKERRXX(ierr);
565:   PetscViewerFileSetMode(viewer, FILE_MODE_READ);CHKERRXX(ierr);
566:   PetscViewerFileSetName(viewer, bcFilename.c_str());CHKERRXX(ierr);
567:   PetscViewerASCIIGetPointer(viewer, &f);CHKERRXX(ierr);
568:   // Create IBC section
569:   int                          numBdFaces = atoi(strtok(fgets(buf, 2048, f), " "));
570:   int                          *tmpIBC    = new int[numBdFaces*4];
571:   std::map<int,std::set<int> > elem2Idx;
572:   std::map<int,int>            bfReorder;
573:   for (int bf = 0; bf < numBdFaces; bf++) {
574:     const char *x = strtok(fgets(buf, 2048, f), " ");

576:     // Ignore boundary face number
577:     x = strtok(NULL, " "); tmpIBC[bf*4+0] = atoi(x);
578:     x = strtok(NULL, " "); tmpIBC[bf*4+1] = atoi(x);
579:     x = strtok(NULL, " "); tmpIBC[bf*4+2] = atoi(x);
580:     x = strtok(NULL, " "); tmpIBC[bf*4+3] = atoi(x);
581:     const int elem = tmpIBC[bf*4+0]-1;

583:     ibc->addFiberDimension(elem, 4);
584:     ibcnum->addFiberDimension(elem, 1);
585:     ibfcon->addFiberDimension(elem, 2);
586:     bl->addFiberDimension(elem, 1);
587:     bnvec->addFiberDimension(elem, 2);
588:     elem2Idx[elem].insert(bf);
589:   }
590:   mesh->allocate(ibc);
591:   mesh->allocate(ibcnum);
592:   mesh->allocate(ibfcon);
593:   mesh->allocate(bl);
594:   mesh->allocate(bnvec);
595:   const DMMesh::int_section_type::chart_type& chart = ibc->getChart();

597:   int num = 1;

599:   for (Mesh::int_section_type::chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
600:     const int elem = *p_iter;
601:     int       bfNum[2];
602:     int       k = 0;

604:     for (std::set<int>::const_iterator i_iter = elem2Idx[elem].begin(); i_iter != elem2Idx[elem].end(); ++i_iter) {
605:       bfReorder[(*i_iter)+1] = num;
606:       bfNum[k++]             = num;
607:       num++;
608:     }
609:     ibcnum->updatePoint(elem, bfNum);
610:   }
611:   for (int bf = 0; bf < numBdFaces; bf++) {
612:     const int elem = tmpIBC[bf*4]-1;

614:     if (elem2Idx[elem].size() > 1) {
615:       if (*elem2Idx[elem].begin() == bf) {
616:         int values[8];
617:         int k = 0;

619:         for (std::set<int>::const_iterator i_iter = elem2Idx[elem].begin(); i_iter != elem2Idx[elem].end(); ++i_iter) {
620:           for (int v = 0; v < 4; ++v) {
621:             values[k*4+v] = tmpIBC[*i_iter*4+v];
622:           }
623:           k++;
624:         }
625:         ibc->updatePoint(elem, values);
626:       }
627:     } else {
628:       ibc->updatePoint(elem, &tmpIBC[bf*4]);
629:     }
630:   }
631:   delete [] tmpIBC;
632:   // Create BCFUNC section
633:   int numBcFunc = atoi(strtok(fgets(buf, 2048, f), " "));
634:   if (numBcFunc != 0) throw ALE::Exception("Cannot handle BCFUNCS after rewrite");
635:   for (int bc = 0; bc < numBcFunc; bc++) {
636: #if 0
637:     const char            *x = strtok(fgets(buf, 2048, f), " ");
638:     DMMesh::bc_value_type value;

640:     // Ignore function number
641:     x = strtok(NULL, " "); value.rho = atof(x);
642:     x = strtok(NULL, " "); value.u   = atof(x);
643:     x = strtok(NULL, " "); value.v   = atof(x);
644:     x = strtok(NULL, " "); value.p   = atof(x);
645:     mesh->setBCValue(bc+1, value);
646: #endif
647:   }
648: #if 0
649:   mesh->distributeBCValues();
650: #endif
651:   // Create IBNDFS section
652:   int       numBdVertices = atoi(strtok(fgets(buf, 2048, f), " "));
653:   const int numElements   = mesh->heightStratum(0)->size();
654:   int       *tmpIBNDFS    = new int[numBdVertices*3];

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

659:     // Ignore boundary node number
660:     x = strtok(NULL, " "); tmpIBNDFS[bv*3+0] = atoi(x);
661:     x = strtok(NULL, " "); tmpIBNDFS[bv*3+1] = atoi(x);
662:     x = strtok(NULL, " "); tmpIBNDFS[bv*3+2] = atoi(x);
663:     ibndfs->setFiberDimension(tmpIBNDFS[bv*3+0]-1+numElements, 6);
664:   }
665:   mesh->allocate(ibndfs);
666:   for (int bv = 0; bv < numBdVertices; bv++) {
667:     int values[5];

669:     values[0] = tmpIBNDFS[bv*3+0];
670:     // Covert to new boundary face numbers
671:     values[1] = bfReorder[tmpIBNDFS[bv*3+1]];
672:     values[2] = bfReorder[tmpIBNDFS[bv*3+2]];
673:     values[3] = 0;
674:     values[4] = 0;
675:     ibndfs->updatePoint(values[0]-1+numElements, values);
676:   }
677:   PetscViewerDestroy(&viewer);CHKERRXX(ierr);
678:   // Create BNNV[NBN,2]
679:   const int dim = mesh->getDimension();

681:   for (int bv = 0; bv < numBdVertices; bv++) {
682:     bnnv->setFiberDimension(tmpIBNDFS[bv*3+0]-1+numElements, dim);
683:   }
684:   mesh->allocate(bnnv);
685:   delete [] tmpIBNDFS;
686: };
687: void Builder::outputVerticesLocal(const Obj<Mesh>& mesh, int *numVertices, int *dim, double *coordinates[], const bool columnMajor)
688: {
689:   const Obj<Mesh::real_section_type>& coordSec = mesh->getRealSection("coordinates");
690:   if (!coordSec->size()) {
691:     *numVertices = 0;
692:     *dim         = 0;
693:     *coordinates = NULL;
694:     return;
695:   }
696:   const Obj<Mesh::label_sequence>& vertices   = mesh->depthStratum(0);
697:   const Obj<Mesh::numbering_type>& vNumbering = mesh->getFactory()->getLocalNumbering(mesh, 0);

699:   int            size     = vertices->size();
700:   int            embedDim = coordSec->getFiberDimension(*vertices->begin());
701:   double        *coords;

704:   PetscMalloc(vertices->size()*embedDim * sizeof(double), &coords);CHKERRXX(ierr);
705:   for (Mesh::label_sequence::iterator v_iter = vertices->begin(); v_iter != vertices->end(); ++v_iter) {
706:     const DMMesh::real_section_type::value_type *array = coordSec->restrictPoint(*v_iter);

708:     const int row = vNumbering->getIndex(*v_iter);

710:     if (columnMajor) {
711:       for (int d = 0; d < embedDim; d++) coords[d*size + row] = array[d];
712:     } else {
713:       for (int d = 0; d < embedDim; d++) coords[row*embedDim + d] = array[d];
714:     }
715:   }
716:   *numVertices = size;
717:   *dim         = embedDim;
718:   *coordinates = coords;
719: };
720: void Builder::outputElementsLocal(const Obj<Mesh>& mesh, int *numElements, int *numCorners, int *vertices[], const bool columnMajor)
721: {
722:   if (!mesh->heightStratum(0)->size()) {
723:     *numElements = 0;
724:     *numCorners  = 0;
725:     *vertices    = NULL;
726:     return;
727:   }
728:   const Obj<Mesh::sieve_type>&     sieve      = mesh->getSieve();
729:   const Obj<Mesh::label_sequence>& elements   = mesh->heightStratum(0);
730:   const Obj<Mesh::numbering_type>& eNumbering = mesh->getFactory()->getLocalNumbering(mesh, mesh->depth());
731:   const Obj<Mesh::numbering_type>& vNumbering = mesh->getFactory()->getLocalNumbering(mesh, 0);

733:   int            size    = elements->size();
734:   int            corners = sieve->getConeSize(*elements->begin());
735:   int            *v;

738:   PetscMalloc(elements->size()*corners * sizeof(int), &v);CHKERRXX(ierr);
739:   for (Mesh::label_sequence::iterator e_iter = elements->begin(); e_iter != elements->end(); ++e_iter) {
740:     const Obj<Mesh::sieve_type::traits::coneSequence>  cone  = sieve->cone(*e_iter);
741:     DMMesh::sieve_type::traits::coneSequence::iterator begin = cone->begin();
742:     DMMesh::sieve_type::traits::coneSequence::iterator end   = cone->end();

744:     const int row = eNumbering->getIndex(*e_iter);
745:     int       c   = -1;
746:     if (columnMajor) {
747:       for (Mesh::sieve_type::traits::coneSequence::iterator c_iter = begin; c_iter != end; ++c_iter) {
748:         v[(++c)*size + row] = vNumbering->getIndex(*c_iter)+1;
749:       }
750:     } else {
751:       for (Mesh::sieve_type::traits::coneSequence::iterator c_iter = begin; c_iter != end; ++c_iter) {
752:         v[row*corners + ++c] = vNumbering->getIndex(*c_iter)+1;
753:       }
754:     }
755:   }
756:   *numElements = size;
757:   *numCorners  = corners;
758:   *vertices    = v;
759: };

763: PetscErrorCode Viewer::writeVertices(const ALE::Obj<Mesh>& mesh, PetscViewer viewer)
764: {
765:   ALE::Obj<Mesh::real_section_type> coordinates = mesh->getRealSection("coordinates");
766: #if 0
767:   DMMesh::field_type::patch_type patch;
768:   const double                   *array = coordinates->restrict (patch);
769:   int                            numVertices;
770:   PetscErrorCode                 ierr;

773:   //FIX:
774:   if (vertexBundle->getGlobalOffsets()) {
775:     numVertices = vertexBundle->getGlobalOffsets()[mesh->commSize()];
776:   } else {
777:     numVertices = mesh->getTopology()->depthStratum(0)->size();
778:   }
779:   PetscViewerASCIIPrintf(viewer, "%D\n", numVertices);
780:   if (mesh->commRank() == 0) {
781:     int numLocalVertices = mesh->getTopology()->depthStratum(0)->size();
782:     int embedDim         = coordinates->getFiberDimension(patch, *mesh->getTopology()->depthStratum(0)->begin());
783:     int vertexCount      = 1;

785:     for (int v = 0; v < numLocalVertices; v++) {
786:       PetscViewerASCIIPrintf(viewer, "%7D   ", vertexCount++);
787:       for (int d = 0; d < embedDim; d++) {
788:         if (d > 0) {
789:           PetscViewerASCIIPrintf(viewer, " ");
790:         }
791:         PetscViewerASCIIPrintf(viewer, "% 12.5E", array[v*embedDim+d]);
792:       }
793:       PetscViewerASCIIPrintf(viewer, "\n");
794:     }
795:     for (int p = 1; p < mesh->commSize(); p++) {
796:       double     *remoteCoords;
797:       MPI_Status status;

799:       MPI_Recv(&numLocalVertices, 1, MPI_INT, p, 1, mesh->comm(), &status);
800:       PetscMalloc(numLocalVertices*embedDim * sizeof(double), &remoteCoords);
801:       MPI_Recv(remoteCoords, numLocalVertices*embedDim, MPI_DOUBLE, p, 1, mesh->comm(), &status);
802:       for (int v = 0; v < numLocalVertices; v++) {
803:         PetscViewerASCIIPrintf(viewer,"%7D   ", vertexCount++);
804:         for (int d = 0; d < embedDim; d++) {
805:           if (d > 0) {
806:             PetscViewerASCIIPrintf(viewer, " ");
807:           }
808:           PetscViewerASCIIPrintf(viewer, "% 12.5E", remoteCoords[v*embedDim+d]);
809:         }
810:         PetscViewerASCIIPrintf(viewer, "\n");
811:       }
812:     }
813:   } else {
814:     ALE::Obj<Mesh::bundle_type>                           globalOrder      = coordinates->getGlobalOrder();
815:     ALE::Obj<Mesh::bundle_type::order_type::coneSequence> cone             = globalOrder->getPatch(patch);
816:     const int                                             *offsets         = coordinates->getGlobalOffsets();
817:     int                                                   embedDim         = coordinates->getFiberDimension(patch, *mesh->getTopology()->depthStratum(0)->begin());
818:     int                                                   numLocalVertices = (offsets[mesh->commRank()+1] - offsets[mesh->commRank()])/embedDim;
819:     double                                                *localCoords;
820:     int                                                   k = 0;

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

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

829:         for (int i = offset; i < offset+dim; ++i) localCoords[k++] = array[i];
830:       }
831:     }
832:     if (k != numLocalVertices*embedDim) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB, "Invalid number of coordinates to send %d should be %d", k, numLocalVertices*embedDim);
833:     MPI_Send(&numLocalVertices, 1, MPI_INT, 0, 1, mesh->comm());
834:     MPI_Send(localCoords, numLocalVertices*embedDim, MPI_DOUBLE, 0, 1, mesh->comm());
835:     PetscFree(localCoords);
836:   }
837: #endif
838:   return(0);
839: };

843: PetscErrorCode Viewer::writeElements(const ALE::Obj<Mesh>& mesh, PetscViewer viewer)
844: {
845: #if 0
846:   ALE::Obj<Mesh::sieve_type::traits::heightSequence> elements      = topology->heightStratum(0);
847:   ALE::Obj<Mesh::bundle_type>                        elementBundle = mesh->getBundle(topology->depth());
848:   ALE::Obj<Mesh::bundle_type>                        vertexBundle  = mesh->getBundle(0);
849:   ALE::Obj<Mesh::bundle_type>                        globalVertex  = vertexBundle->getGlobalOrder();
850:   ALE::Obj<Mesh::bundle_type>                        globalElement = elementBundle->getGlobalOrder();
851:   DMMesh::bundle_type::patch_type                    patch;
852:   std::string    orderName("element");
853:   int            dim     = mesh->getDimension();
854:   int            corners = topology->nCone(*elements->begin(), topology->depth())->size();
855:   int            numElements;

859:   if (corners != dim+1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP, "PCICE only supports simplicies");
860:   if (!globalVertex) globalVertex = vertexBundle;

862:   if (elementBundle->getGlobalOffsets()) {
863:     numElements = elementBundle->getGlobalOffsets()[mesh->commSize()];
864:   } else {
865:     numElements = mesh->getTopology()->heightStratum(0)->size();
866:   }
867:   if (mesh->commRank() == 0) {
868:     int elementCount = 1;

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

874:       PetscViewerASCIIPrintf(viewer, "%7d", elementCount++);
875:       for (Mesh::bundle_type::order_type::coneSequence::iterator c_itor = cone->begin(); c_itor != cone->end(); ++c_itor) {
876:         PetscViewerASCIIPrintf(viewer, " %7d", globalVertex->getIndex(patch, *c_itor).prefix);
877:       }
878:       PetscViewerASCIIPrintf(viewer, "\n");
879:     }
880:     for (int p = 1; p < mesh->commSize(); p++) {
881:       int        numLocalElements;
882:       int        *remoteVertices;
883:       MPI_Status status;

885:       MPI_Recv(&numLocalElements, 1, MPI_INT, p, 1, mesh->comm(), &status);
886:       PetscMalloc(numLocalElements*corners * sizeof(int), &remoteVertices);
887:       MPI_Recv(remoteVertices, numLocalElements*corners, MPI_INT, p, 1, mesh->comm(), &status);
888:       for (int e = 0; e < numLocalElements; e++) {
889:         PetscViewerASCIIPrintf(viewer, "%7d", elementCount++);
890:         for (int c = 0; c < corners; c++) {
891:           PetscViewerASCIIPrintf(viewer, " %7d", remoteVertices[e*corners+c]);
892:         }
893:         PetscViewerASCIIPrintf(viewer, "\n");
894:       }
895:       PetscFree(remoteVertices);
896:     }
897:   } else {
898:     const int *offsets         = elementBundle->getGlobalOffsets();
899:     int       numLocalElements = offsets[mesh->commRank()+1] - offsets[mesh->commRank()];
900:     int       *localVertices;
901:     int       k = 0;

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

907:       if (globalElement->getFiberDimension(patch, *e_itor) > 0) {
908:         for (Mesh::bundle_type::order_type::coneSequence::iterator c_itor = cone->begin(); c_itor != cone->end(); ++c_itor) {
909:           localVertices[k++] = globalVertex->getIndex(patch, *c_itor).prefix;
910:         }
911:       }
912:     }
913:     if (k != numLocalElements*corners) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB, "Invalid number of vertices to send %d should be %d", k, numLocalElements*corners);
914:     MPI_Send(&numLocalElements, 1, MPI_INT, 0, 1, mesh->comm());
915:     MPI_Send(localVertices, numLocalElements*corners, MPI_INT, 0, 1, mesh->comm());
916:     PetscFree(localVertices);
917:   }
918: #endif
919:   return(0);
920: };

924: PetscErrorCode Viewer::writeVerticesLocal(const Obj<Mesh>& mesh, PetscViewer viewer)
925: {
926:   Obj<Mesh::real_section_type>     coordinates = mesh->getRealSection("coordinates");
927:   const Obj<Mesh::label_sequence>& vertices    = mesh->depthStratum(0);
928:   const Obj<Mesh::numbering_type>& vNumbering  = mesh->getFactory()->getLocalNumbering(mesh, 0);
929:   int                              embedDim    = coordinates->getFiberDimension(*vertices->begin());
930:   PetscErrorCode                   ierr;

933:   PetscViewerASCIIPrintf(viewer, "%D\n", vertices->size());
934:   for (Mesh::label_sequence::iterator v_iter = vertices->begin(); v_iter != vertices->end(); ++v_iter) {
935:     const DMMesh::real_section_type::value_type *array = coordinates->restrictPoint(*v_iter);

937:     PetscViewerASCIIPrintf(viewer, "%7D   ", vNumbering->getIndex(*v_iter)+1);
938:     for (int d = 0; d < embedDim; d++) {
939:       if (d > 0) PetscViewerASCIIPrintf(viewer, " ");
940:       PetscViewerASCIIPrintf(viewer, "% 12.5E", array[d]);
941:     }
942:     PetscViewerASCIIPrintf(viewer, "\n");
943:   }
944:   return(0);
945: };

949: PetscErrorCode Viewer::writeRestart(const Obj<Mesh>& mesh, PetscViewer viewer)
950: {
951:   const Obj<Mesh::real_section_type>& velocity    = mesh->getRealSection("VELN");
952:   const Obj<Mesh::real_section_type>& pressure    = mesh->getRealSection("PN");
953:   const Obj<Mesh::real_section_type>& temperature = mesh->getRealSection("TN");
954:   const Obj<Mesh::numbering_type>&    cNumbering  = mesh->getFactory()->getNumbering(mesh, mesh->depth());
955:   const Obj<Mesh::numbering_type>&    vNumbering  = mesh->getFactory()->getNumbering(mesh, 0);
956:   const int                           numCells    = cNumbering->getGlobalSize();
957:   PetscErrorCode                      ierr;

960:   int          blen[2];
961:   MPI_Aint     indices[2];
962:   MPI_Datatype oldtypes[2], newtype;
963:   blen[0] = 1; indices[0] = 0;           oldtypes[0] = MPI_INT;
964:   blen[1] = 4; indices[1] = sizeof(int); oldtypes[1] = MPI_DOUBLE;

966:   MPI_Type_create_struct(2, blen, indices, oldtypes, &newtype);
967:   MPI_Type_commit(&newtype);

969:   if (mesh->commRank() == 0) {
970:     const Obj<Mesh::label_sequence>& vertices = mesh->depthStratum(0);

972:     for (Mesh::label_sequence::iterator v_iter = vertices->begin(); v_iter != vertices->end(); ++v_iter) {
973:       if (vNumbering->isLocal(*v_iter)) {
974:         const DMMesh::real_section_type::value_type *veln = velocity->restrictPoint(*v_iter);
975:         const DMMesh::real_section_type::value_type *pn   = pressure->restrictPoint(*v_iter);
976:         const DMMesh::real_section_type::value_type *tn   = temperature->restrictPoint(*v_iter);

978:         PetscViewerASCIIPrintf(viewer, "%6d% 16.8E% 16.8E% 16.8E% 16.8E\n", *v_iter-numCells+1, veln[0], veln[1], pn[0], tn[0]);
979:       }
980:     }
981:     for (int p = 1; p < mesh->commSize(); p++) {
982:       RestartType *remoteValues;
983:       int         numLocalElements;
984:       MPI_Status  status;

986:       MPI_Recv(&numLocalElements, 1, MPI_INT, p, 1, mesh->comm(), &status);
987:       PetscMalloc(numLocalElements * sizeof(RestartType), &remoteValues);
988:       MPI_Recv(remoteValues, numLocalElements, newtype, p, 1, mesh->comm(), &status);
989:       for (int e = 0; e < numLocalElements; e++) {
990:         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);
991:       }
992:     }
993:   } else {
994:     const Obj<Mesh::label_sequence>& vertices = mesh->depthStratum(0);
995:     RestartType                      *localValues;
996:     int                              numLocalElements = vNumbering->getLocalSize();
997:     int                              k                = 0;

999:     PetscMalloc(numLocalElements * sizeof(RestartType), &localValues);
1000:     for (Mesh::label_sequence::iterator v_iter = vertices->begin(); v_iter != vertices->end(); ++v_iter) {
1001:       if (vNumbering->isLocal(*v_iter)) {
1002:         const DMMesh::real_section_type::value_type *veln = velocity->restrictPoint(*v_iter);
1003:         const DMMesh::real_section_type::value_type *pn   = pressure->restrictPoint(*v_iter);
1004:         const DMMesh::real_section_type::value_type *tn   = temperature->restrictPoint(*v_iter);

1006:         localValues[k].vertex = *v_iter;
1007:         localValues[k].veln_x = veln[0];
1008:         localValues[k].veln_y = veln[1];
1009:         localValues[k].pn     = pn[0];
1010:         localValues[k].tn     = tn[0];
1011:         k++;
1012:       }
1013:     }
1014:     if (k != numLocalElements) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB, "Invalid number of values to send for field, %d should be %d", k, numLocalElements);
1015:     MPI_Send(&numLocalElements, 1, MPI_INT, 0, 1, mesh->comm());
1016:     MPI_Send(localValues, numLocalElements, newtype, 0, 1, mesh->comm());
1017:     PetscFree(localValues);
1018:   }
1019:   MPI_Type_free(&newtype);
1020:   return(0);
1021: };

1023: //   This class reconstructs the local pieces of the boundary that distributed PCICE needs.
1024: // The boundary along with the boundary conditions is encoded in a collection of sections
1025: // over the PCICE mesh.  These sections contain a description of the boundary topology
1026: // using elements' global names.  This is unacceptable for PCICE, since it interprets
1027: // elements of the connectivity data arrays as local offsets into (some other of) these arrays.
1028: //   This subroutine performs the renumbering based on the local numbering on the distributed
1029: // mesh.  Essentially, we replace each global element id by its local number.
1030: //
1031: // Note: Any vertex or element number from PCICE is 1-based, but in Sieve we are 0-based. Thus
1032: //       we add and subtract 1 during conversion. Also, Sieve vertices are numbered after cells.
1033: void fuseBoundary(const ALE::Obj<PETSC_MESH_TYPE>& mesh)
1034: {
1035:   // Extract PCICE boundary sections
1036:   ALE::Obj<PETSC_MESH_TYPE::int_section_type> IBCsec    = mesh->getIntSection("IBC");
1037:   ALE::Obj<PETSC_MESH_TYPE::int_section_type> IBNDFSsec = mesh->getIntSection("IBNDFS");
1038:   ALE::Obj<PETSC_MESH_TYPE::int_section_type> IBCNUMsec = mesh->getIntSection("IBCNUM");

1040:   // Look at the sections, if debugging
1041:   if (mesh->debug()) {
1042:     IBCsec->view("IBCsec", mesh->comm());
1043:     IBNDFSsec->view("IBNDFSsec", mesh->comm());
1044:   }

1046:   // Extract the local numberings
1047:   ALE::Obj<PETSC_MESH_TYPE::numbering_type> vertexNum   = mesh->getFactory()->getLocalNumbering(mesh, 0);
1048:   ALE::Obj<PETSC_MESH_TYPE::numbering_type> elementNum  = mesh->getFactory()->getLocalNumbering(mesh, mesh->depth());
1049:   const int                                 numElements = mesh->getFactory()->getNumbering(mesh, mesh->depth())->getGlobalSize();
1050:   std::map<int,int>                         bfMap;
1051:   // Declare points to the extracted numbering data
1052:   const PETSC_MESH_TYPE::numbering_type::value_type *vNum, *eNum;

1054:   // Create map from serial bdFace numbers to local bdFace numbers
1055:   {
1056:     const PETSC_MESH_TYPE::int_section_type::chart_type           chart = IBCNUMsec->getChart();
1057:     PETSC_MESH_TYPE::int_section_type::chart_type::const_iterator begin = chart.begin();
1058:     PETSC_MESH_TYPE::int_section_type::chart_type::const_iterator end   = chart.end();

1060:     int num = 0;
1061:     for (PETSC_MESH_TYPE::int_section_type::chart_type::const_iterator p_iter = begin; p_iter != end; ++p_iter) {
1062:       const int fiberDim   = IBCNUMsec->getFiberDimension(*p_iter);
1063:       const int *globalNum = IBCNUMsec->restrictPoint(*p_iter);

1065:       for (int n = 0; n < fiberDim; ++n) bfMap[globalNum[n]] = ++num;
1066:     }
1067:   }
1068:   // Renumber vertex section IBC
1069:   {
1070:     const PETSC_MESH_TYPE::int_section_type::chart_type           IBCchart = IBCsec->getChart();
1071:     PETSC_MESH_TYPE::int_section_type::chart_type::const_iterator begin    = IBCchart.begin();
1072:     PETSC_MESH_TYPE::int_section_type::chart_type::const_iterator end      = IBCchart.end();
1073:     for (PETSC_MESH_TYPE::int_section_type::chart_type::const_iterator p_iter = begin; p_iter != end; ++p_iter) {
1074:       PETSC_MESH_TYPE::point_type                         p              = *p_iter;
1075:       const PETSC_MESH_TYPE::int_section_type::value_type *ibc_in        = IBCsec->restrictPoint(p);
1076:       int                                                 fiberDimension = IBCsec->getFiberDimension(p);
1077:       PETSC_MESH_TYPE::int_section_type::value_type       ibc_out[8];
1078:       // k controls the update of edge connectivity for each containing element;
1079:       // if fiberDimension is 4, only one boundary face is connected to the element, and that edge's data
1080:       // are contained in entries 0 - 3 of the section over the element p;
1081:       // if fiberDimension is 8, two boundary faces are connected to the element, so the second edge's data
1082:       // are contained in entries 4 - 7
1083:       for (int k = 0; k < fiberDimension/4; k++) {
1084:         // Extract IBC entry 1 (entry kk*4) for edge kk connected to element p.
1085:         // This is the entry that needs renumbering for renumbering (2,3 & 4 are invariant under distribution),
1086:         // see IBC's description.
1087:         // Here we assume that elementNum's domain contains all boundary elements found in IBC,
1088:         // so we can restrict to the extracted entry.
1089:         eNum           = elementNum->restrictPoint((PETSC_MESH_TYPE::numbering_type::point_type) ibc_in[k*4]-1);
1090:         ibc_out[k*4+0] = eNum[0]+1;
1091:         // Copy the other entries right over
1092:         ibc_out[k*4+1] = ibc_in[k*4+1];
1093:         ibc_out[k*4+2] = ibc_in[k*4+2];
1094:         ibc_out[k*4+3] = ibc_in[k*4+3];
1095:       }
1096:       // Update IBC
1097:       IBCsec->updatePoint(p, ibc_out);
1098:     }
1099:   }
1100:   {
1101:     // Renumber vertex section IBNDFS
1102:     const PETSC_MESH_TYPE::int_section_type::chart_type           IBNDFSchart = IBNDFSsec->getChart();
1103:     PETSC_MESH_TYPE::int_section_type::chart_type::const_iterator begin       = IBNDFSchart.begin();
1104:     PETSC_MESH_TYPE::int_section_type::chart_type::const_iterator end         = IBNDFSchart.end();
1105:     for (PETSC_MESH_TYPE::int_section_type::chart_type::const_iterator p_iter = begin; p_iter != end; ++p_iter) {
1106:       PETSC_MESH_TYPE::point_type                         p          = *p_iter;
1107:       const PETSC_MESH_TYPE::int_section_type::value_type *ibndfs_in = IBNDFSsec->restrictPoint(p);
1108:       // Here we assume the fiber dimension is 5
1109:       PETSC_MESH_TYPE::int_section_type::value_type ibndfs_out[5];
1110:       // Renumber entries 1,2 & 3 (4 & 5 are invariant under distribution), see IBNDFS's description
1111:       // Here we assume that vertexNum's domain contains all boundary verticies found in IBNDFS, so we can restrict to the first extracted entry
1112:       vNum          = vertexNum->restrictPoint((PETSC_MESH_TYPE::numbering_type::point_type)ibndfs_in[0]-1+numElements);
1113:       ibndfs_out[0] = vNum[0]+1;
1114:       // Map serial bdFace numbers to local bdFace numbers
1115:       ibndfs_out[1] = bfMap[ibndfs_in[1]];
1116:       ibndfs_out[2] = bfMap[ibndfs_in[2]];
1117:       // Copy the other entries right over
1118:       ibndfs_out[3] = ibndfs_in[3];
1119:       ibndfs_out[4] = ibndfs_in[4];
1120:       // Update IBNDFS
1121:       IBNDFSsec->updatePoint(p,ibndfs_out);
1122:     }
1123:   }
1124:   if (mesh->debug()) {
1125:     IBCsec->view("Renumbered IBCsec", mesh->comm());
1126:     IBNDFSsec->view("Renumbered IBNDFSsec", mesh->comm());
1127:   }
1128:   // It's not clear whether IBFCON needs to be renumbered (what does it mean that its entries are not "GLOBAL NODE NUMBER(s)" -- see IBFCON's descriptions
1129: };
1130: #endif
1131: };
1132: };