Actual source code: mesh.c

  1: 
 2:  #include src/dm/mesh/meshimpl.h
  3: #include <Distribution.hh>
  4: #include <Generator.hh>
 5:  #include src/dm/mesh/meshvtk.h
 6:  #include src/dm/mesh/meshpcice.h
 7:  #include src/dm/mesh/meshpylith.h

  9: /* Logging support */
 10: PetscCookie  MESH_COOKIE = 0;
 11: PetscEvent  Mesh_View = 0, Mesh_GetGlobalScatter = 0, Mesh_restrictVector = 0, Mesh_assembleVector = 0,
 12:             Mesh_assembleVectorComplete = 0, Mesh_assembleMatrix = 0, Mesh_updateOperator = 0;

 14: EXTERN PetscErrorCode MeshView_DM(Mesh, PetscViewer);
 15: EXTERN PetscErrorCode MeshRefine_DM(Mesh, MPI_Comm, Mesh *);
 16: EXTERN PetscErrorCode MeshCoarsenHierarchy_DM(Mesh, int, Mesh **);
 17: EXTERN PetscErrorCode MeshGetInterpolation_DM(DM, DM, Mat *, Vec *);

 21: PetscErrorCode MeshFinalize()
 22: {
 24:   ALE::Mesh::NumberingFactory::singleton(0, true);
 25:   return(0);
 26: }

 30: PetscErrorCode MeshView_Sieve_Ascii(const ALE::Obj<ALE::Mesh>& mesh, PetscViewer viewer)
 31: {
 32:   PetscViewerFormat format;
 33:   PetscErrorCode    ierr;

 36:   PetscViewerGetFormat(viewer, &format);
 37:   if (format == PETSC_VIEWER_ASCII_VTK) {
 38:     VTKViewer::writeHeader(viewer);
 39:     VTKViewer::writeVertices(mesh, 0, viewer);
 40:     VTKViewer::writeElements(mesh, 0, viewer);
 41:   } else if (format == PETSC_VIEWER_ASCII_PYLITH) {
 42:     char *filename;
 43:     char  connectFilename[2048];
 44:     char  coordFilename[2048];

 46:     PetscViewerFileGetName(viewer, &filename);
 47:     PetscViewerFileSetMode(viewer, FILE_MODE_WRITE);
 48:     PetscStrcpy(connectFilename, filename);
 49:     PetscStrcat(connectFilename, ".connect");
 50:     PetscViewerFileSetName(viewer, connectFilename);
 51:     ALE::PyLith::Viewer::writeElements(mesh, mesh->getIntSection("material"), viewer);
 52:     PetscStrcpy(coordFilename, filename);
 53:     PetscStrcat(coordFilename, ".coord");
 54:     PetscViewerFileSetName(viewer, coordFilename);
 55:     ALE::PyLith::Viewer::writeVertices(mesh, viewer);
 56:     PetscViewerFileSetMode(viewer, FILE_MODE_READ);
 57:     PetscExceptionTry1(PetscViewerFileSetName(viewer, filename), PETSC_ERR_FILE_OPEN);
 58:     if (PetscExceptionValue(ierr)) {
 59:       /* this means that a caller above me has also tryed this exception so I don't handle it here, pass it up */
 60:     } else if (PetscExceptionCaught(ierr, PETSC_ERR_FILE_OPEN)) {
 61:       0;
 62:     }
 63: 
 64:   } else if (format == PETSC_VIEWER_ASCII_PYLITH_LOCAL) {
 65:     PetscViewer connectViewer, coordViewer, splitViewer, tractionViewer;
 66:     char       *filename;
 67:     char        localFilename[2048];
 68:     int         rank = mesh->commRank();

 70:     PetscViewerFileGetName(viewer, &filename);

 72:     sprintf(localFilename, "%s.%d.connect", filename, rank);
 73:     PetscViewerCreate(PETSC_COMM_SELF, &connectViewer);
 74:     PetscViewerSetType(connectViewer, PETSC_VIEWER_ASCII);
 75:     PetscViewerSetFormat(connectViewer, PETSC_VIEWER_ASCII_PYLITH);
 76:     PetscViewerFileSetName(connectViewer, localFilename);
 77:     ALE::PyLith::Viewer::writeElementsLocal(mesh, mesh->getIntSection("material"), connectViewer);
 78:     PetscViewerDestroy(connectViewer);

 80:     sprintf(localFilename, "%s.%d.coord", filename, rank);
 81:     PetscViewerCreate(PETSC_COMM_SELF, &coordViewer);
 82:     PetscViewerSetType(coordViewer, PETSC_VIEWER_ASCII);
 83:     PetscViewerSetFormat(coordViewer, PETSC_VIEWER_ASCII_PYLITH);
 84:     PetscViewerFileSetName(coordViewer, localFilename);
 85:     ALE::PyLith::Viewer::writeVerticesLocal(mesh, coordViewer);
 86:     PetscViewerDestroy(coordViewer);

 88:     if (mesh->hasPairSection("split")) {
 89:       sprintf(localFilename, "%s.%d.split", filename, rank);
 90:       PetscViewerCreate(PETSC_COMM_SELF, &splitViewer);
 91:       PetscViewerSetType(splitViewer, PETSC_VIEWER_ASCII);
 92:       PetscViewerSetFormat(splitViewer, PETSC_VIEWER_ASCII_PYLITH);
 93:       PetscViewerFileSetName(splitViewer, localFilename);
 94:       ALE::PyLith::Viewer::writeSplitLocal(mesh, mesh->getPairSection("split"), splitViewer);
 95:       PetscViewerDestroy(splitViewer);
 96:     }

 98:     if (mesh->hasRealSection("traction")) {
 99:       sprintf(localFilename, "%s.%d.traction", filename, rank);
100:       PetscViewerCreate(PETSC_COMM_SELF, &tractionViewer);
101:       PetscViewerSetType(tractionViewer, PETSC_VIEWER_ASCII);
102:       PetscViewerSetFormat(tractionViewer, PETSC_VIEWER_ASCII_PYLITH);
103:       PetscViewerFileSetName(tractionViewer, localFilename);
104:       ALE::PyLith::Viewer::writeTractionsLocal(mesh, mesh->getRealSection("traction"), tractionViewer);
105:       PetscViewerDestroy(tractionViewer);
106:     }
107:   } else if (format == PETSC_VIEWER_ASCII_PCICE) {
108:     char      *filename;
109:     char       coordFilename[2048];
110:     PetscTruth isConnect;
111:     size_t     len;

113:     PetscViewerFileGetName(viewer, &filename);
114:     PetscStrlen(filename, &len);
115:     PetscStrcmp(&(filename[len-5]), ".lcon", &isConnect);
116:     if (!isConnect) {
117:       SETERRQ1(PETSC_ERR_ARG_WRONG, "Invalid element connectivity filename: %s", filename);
118:     }
119:     ALE::PCICE::Viewer::writeElements(mesh, viewer);
120:     PetscStrncpy(coordFilename, filename, len-5);
121:     coordFilename[len-5] = '\0';
122:     PetscStrcat(coordFilename, ".nodes");
123:     PetscViewerFileSetName(viewer, coordFilename);
124:     ALE::PCICE::Viewer::writeVertices(mesh, viewer);
125:   } else {
126:     int dim = mesh->getDimension();

128:     PetscViewerASCIIPrintf(viewer, "Mesh in %d dimensions:\n", dim);
129:     for(int d = 0; d <= dim; d++) {
130:       // FIX: Need to globalize
131:       PetscViewerASCIIPrintf(viewer, "  %d %d-cells\n", mesh->getTopology()->depthStratum(0, d)->size(), d);
132:     }
133:   }
134:   PetscViewerFlush(viewer);
135:   return(0);
136: }

140: PetscErrorCode MeshView_Sieve(const ALE::Obj<ALE::Mesh>& mesh, PetscViewer viewer)
141: {
142:   PetscTruth     iascii, isbinary, isdraw;

146:   PetscTypeCompare((PetscObject) viewer, PETSC_VIEWER_ASCII, &iascii);
147:   PetscTypeCompare((PetscObject) viewer, PETSC_VIEWER_BINARY, &isbinary);
148:   PetscTypeCompare((PetscObject) viewer, PETSC_VIEWER_DRAW, &isdraw);

150:   if (iascii){
151:     MeshView_Sieve_Ascii(mesh, viewer);
152:   } else if (isbinary) {
153:     SETERRQ(PETSC_ERR_SUP, "Binary viewer not implemented for Mesh");
154:   } else if (isdraw){
155:     SETERRQ(PETSC_ERR_SUP, "Draw viewer not implemented for Mesh");
156:   } else {
157:     SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported by this mesh object", ((PetscObject)viewer)->type_name);
158:   }
159:   return(0);
160: }

162: PetscErrorCode  MeshView_DM(Mesh mesh, PetscViewer viewer)
163: {

167:   MeshView_Sieve(mesh->m, viewer);
168:   return(0);
169: }

173: /*@C
174:    MeshView - Views a Mesh object. 

176:    Collective on Mesh

178:    Input Parameters:
179: +  mesh - the mesh
180: -  viewer - an optional visualization context

182:    Notes:
183:    The available visualization contexts include
184: +     PETSC_VIEWER_STDOUT_SELF - standard output (default)
185: -     PETSC_VIEWER_STDOUT_WORLD - synchronized standard
186:          output where only the first processor opens
187:          the file.  All other processors send their 
188:          data to the first processor to print. 

190:    You can change the format the mesh is printed using the 
191:    option PetscViewerSetFormat().

193:    The user can open alternative visualization contexts with
194: +    PetscViewerASCIIOpen() - Outputs mesh to a specified file
195: .    PetscViewerBinaryOpen() - Outputs mesh in binary to a
196:          specified file; corresponding input uses MeshLoad()
197: .    PetscViewerDrawOpen() - Outputs mesh to an X window display

199:    The user can call PetscViewerSetFormat() to specify the output
200:    format of ASCII printed objects (when using PETSC_VIEWER_STDOUT_SELF,
201:    PETSC_VIEWER_STDOUT_WORLD and PetscViewerASCIIOpen).  Available formats include
202: +    PETSC_VIEWER_ASCII_DEFAULT - default, prints mesh information
203: -    PETSC_VIEWER_ASCII_VTK - outputs a VTK file describing the mesh

205:    Level: beginner

207:    Concepts: mesh^printing
208:    Concepts: mesh^saving to disk

210: .seealso: PetscViewerASCIIOpen(), PetscViewerDrawOpen(), PetscViewerBinaryOpen(),
211:           MeshLoad(), PetscViewerCreate()
212: @*/
213: PetscErrorCode  MeshView(Mesh mesh, PetscViewer viewer)
214: {

220:   if (!viewer) {
221:     PetscViewerASCIIGetStdout(mesh->comm,&viewer);
222:   }

227:   (*mesh->ops->view)(mesh, viewer);
229:   return(0);
230: }

234: /*@C
235:     MeshLoad - Create a mesh topology from the saved data in a viewer.

237:     Collective on Viewer

239:     Input Parameter:
240: .   viewer - The viewer containing the data

242:     Output Parameters:
243: .   mesh - the mesh object

245:     Level: advanced

247: .seealso MeshView()

249: @*/
250: PetscErrorCode  MeshLoad(PetscViewer viewer, Mesh *mesh)
251: {
252:   SETERRQ(PETSC_ERR_SUP, "");
253: }

257: /*@C
258:     MeshGetMesh - Gets the internal mesh object

260:     Not collective

262:     Input Parameter:
263: .    mesh - the mesh object

265:     Output Parameter:
266: .    m - the internal mesh object
267:  
268:     Level: advanced

270: .seealso MeshCreate(), MeshSetMesh()

272: @*/
273: PetscErrorCode  MeshGetMesh(Mesh mesh, ALE::Obj<ALE::Mesh>& m)
274: {
277:   m = mesh->m;
278:   return(0);
279: }

283: /*@C
284:     MeshSetMesh - Sets the internal mesh object

286:     Not collective

288:     Input Parameters:
289: +    mesh - the mesh object
290: -    m - the internal mesh object
291:  
292:     Level: advanced

294: .seealso MeshCreate(), MeshGetMesh()

296: @*/
297: PetscErrorCode  MeshSetMesh(Mesh mesh, const ALE::Obj<ALE::Mesh>& m)
298: {
301:   mesh->m = m;
302:   return(0);
303: }

307: template<typename Atlas>
308: PetscErrorCode  MeshCreateMatrix(Mesh mesh, const Obj<Atlas>& atlas, MatType mtype, Mat *J)
309: {
310:   Obj<ALE::Mesh> m;
312:   MeshGetMesh(mesh, m);
313:   const ALE::Obj<ALE::Mesh::order_type>& order = m->getFactory()->getGlobalOrder(m->getTopology(), 0, "default", atlas);
314:   int localSize  = order->getLocalSize();
315:   int globalSize = order->getGlobalSize();

318:   MatCreate(m->comm(), J);
319:   MatSetSizes(*J, localSize, localSize, globalSize, globalSize);
320:   MatSetType(*J, mtype);
321:   MatSetFromOptions(*J);
322:   //MatSetBlockSize(*J, 1);
323:   PetscObjectCompose((PetscObject) *J, "mesh", (PetscObject) mesh);

325:   preallocateOperator(m->getTopology(), atlas, order, *J);
326:   return(0);
327: }

331: PetscErrorCode  MeshGetVertexMatrix(Mesh mesh, MatType mtype, Mat *J)
332: {
333:   ALE::Obj<ALE::Mesh> m;
334:   PetscErrorCode      ierr;

337:   MeshGetMesh(mesh, m);
338:   ALE::Obj<ALE::Mesh::real_section_type> s = new ALE::Mesh::real_section_type(m->getTopology());
339:   s->setFiberDimensionByDepth(0, 0, 1);
340:   s->allocate();
341:   MeshCreateMatrix(mesh, s->getAtlas(), mtype, J);
342:   return(0);
343: }

347: /*@C
348:     MeshGetMatrix - Creates a matrix with the correct parallel layout required for 
349:       computing the Jacobian on a function defined using the informatin in Mesh.

351:     Collective on Mesh

353:     Input Parameter:
354: +   mesh - the mesh object
355: -   mtype - Supported types are MATSEQAIJ, MATMPIAIJ, MATSEQBAIJ, MATMPIBAIJ, MATSEQSBAIJ, MATMPISBAIJ,
356:             or any type which inherits from one of these (such as MATAIJ, MATLUSOL, etc.).

358:     Output Parameters:
359: .   J  - matrix with the correct nonzero preallocation
360:         (obviously without the correct Jacobian values)

362:     Level: advanced

364:     Notes: This properly preallocates the number of nonzeros in the sparse matrix so you 
365:        do not need to do it yourself.

367: .seealso ISColoringView(), ISColoringGetIS(), MatFDColoringCreate(), DASetBlockFills()

369: @*/
370: PetscErrorCode  MeshGetMatrix(Mesh mesh, MatType mtype, Mat *J)
371: {
372:   ALE::Obj<ALE::Mesh> m;
373:   PetscTruth          flag;
374:   PetscErrorCode      ierr;

377:   MeshHasSectionReal(mesh, "default", &flag);
378:   if (!flag) SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "Must set default section");
379:   MeshGetMesh(mesh, m);
380:   MeshCreateMatrix(mesh, m->getRealSection("default")->getAtlas(), mtype, J);
381:   return(0);
382: }

386: /*@C
387:     MeshCreate - Creates a DM object, used to manage data for an unstructured problem
388:     described by a Sieve.

390:     Collective on MPI_Comm

392:     Input Parameter:
393: .   comm - the processors that will share the global vector

395:     Output Parameters:
396: .   mesh - the mesh object

398:     Level: advanced

400: .seealso MeshDestroy(), MeshCreateGlobalVector(), MeshGetGlobalIndices()

402: @*/
403: PetscErrorCode  MeshCreate(MPI_Comm comm,Mesh *mesh)
404: {
406:   Mesh         p;

410:   *mesh = PETSC_NULL;
411: #ifndef PETSC_USE_DYNAMIC_LIBRARIES
412:   DMInitializePackage(PETSC_NULL);
413: #endif

415:   PetscHeaderCreate(p,_p_Mesh,struct _MeshOps,MESH_COOKIE,0,"Mesh",comm,MeshDestroy,0);
416:   p->ops->view               = MeshView_DM;
417:   p->ops->createglobalvector = MeshCreateGlobalVector;
418:   p->ops->getcoloring        = PETSC_NULL;
419:   p->ops->getmatrix          = MeshGetMatrix;
420:   p->ops->getinterpolation   = MeshGetInterpolation_DM;
421:   p->ops->getinjection       = PETSC_NULL;
422:   p->ops->refine             = MeshRefine_DM;
423:   p->ops->coarsen            = PETSC_NULL;
424:   p->ops->refinehierarchy    = PETSC_NULL;
425:   p->ops->coarsenhierarchy   = MeshCoarsenHierarchy_DM;

427:   PetscObjectChangeTypeName((PetscObject) p, "sieve");

429:   p->m             = PETSC_NULL;
430:   p->globalScatter = PETSC_NULL;
431:   p->lf            = PETSC_NULL;
432:   p->lj            = PETSC_NULL;
433:   *mesh = p;
434:   return(0);
435: }

439: /*@C
440:     MeshDestroy - Destroys a mesh.

442:     Collective on Mesh

444:     Input Parameter:
445: .   mesh - the mesh object

447:     Level: advanced

449: .seealso MeshCreate(), MeshCreateGlobalVector(), MeshGetGlobalIndices()

451: @*/
452: PetscErrorCode  MeshDestroy(Mesh mesh)
453: {
454:   PetscErrorCode     ierr;

457:   if (--mesh->refct > 0) return(0);
458:   if (mesh->globalScatter) {VecScatterDestroy(mesh->globalScatter);}
459:   mesh->m = PETSC_NULL;
460:   PetscHeaderDestroy(mesh);
461:   return(0);
462: }

466: inline void ExpandInterval(const ALE::Point& interval, int indices[], int& indx)
467: {
468:   const int end = interval.prefix + interval.index;
469:   for(int i = interval.index; i < end; i++) {
470:     indices[indx++] = i;
471:   }
472: }

476: inline void ExpandInterval_New(ALE::Point interval, PetscInt indices[], PetscInt *indx)
477: {
478:   for(int i = 0; i < interval.prefix; i++) {
479:     indices[(*indx)++] = interval.index + i;
480:   }
481:   for(int i = 0; i < -interval.prefix; i++) {
482:     indices[(*indx)++] = -1;
483:   }
484: }

488: PetscErrorCode ExpandIntervals(ALE::Obj<ALE::Mesh::real_section_type::IndexArray> intervals, PetscInt *indices)
489: {
490:   int k = 0;

493:   for(ALE::Mesh::real_section_type::IndexArray::iterator i_itor = intervals->begin(); i_itor != intervals->end(); i_itor++) {
494:     ExpandInterval_New(*i_itor, indices, &k);
495:   }
496:   return(0);
497: }

501: /*@C
502:     MeshCreateGlobalVector - Creates a vector of the correct size to be gathered into 
503:         by the mesh.

505:     Collective on Mesh

507:     Input Parameter:
508: .    mesh - the mesh object

510:     Output Parameters:
511: .   gvec - the global vector

513:     Level: advanced

515:     Notes: Once this has been created you cannot add additional arrays or vectors to be packed.

517: .seealso MeshDestroy(), MeshCreate(), MeshGetGlobalIndices()

519: @*/
520: PetscErrorCode  MeshCreateGlobalVector(Mesh mesh, Vec *gvec)
521: {
522:   ALE::Obj<ALE::Mesh> m;
523:   PetscTruth     flag;

527:   MeshHasSectionReal(mesh, "default", &flag);
528:   if (!flag) SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "Must set default section");
529:   MeshGetMesh(mesh, m);
530:   const ALE::Obj<ALE::Mesh::order_type>& order = m->getFactory()->getGlobalOrder(m->getTopology(), 0, "default", m->getRealSection("default")->getAtlas());

532:   VecCreate(m->comm(), gvec);
533:   VecSetSizes(*gvec, order->getLocalSize(), order->getGlobalSize());
534:   VecSetFromOptions(*gvec);
535:   return(0);
536: }

540: /*@C
541:   MeshCreateLocalVector - Creates a vector with the local piece of the Section

543:   Collective on Mesh

545:   Input Parameters:
546: + mesh - the Mesh
547: - section - the Section  

549:   Output Parameter:
550: . localVec - the local vector

552:   Level: advanced

554: .seealso MeshDestroy(), MeshCreate()
555: @*/
556: PetscErrorCode  MeshCreateLocalVector(Mesh mesh, SectionReal section, Vec *localVec)
557: {
558:   ALE::Obj<ALE::Mesh::real_section_type> s;

562:   SectionRealGetSection(section, s);
563:   VecCreateSeqWithArray(PETSC_COMM_SELF, s->size(0), s->restrict(0), localVec);
564:   return(0);
565: }

569: /*@C
570:     MeshGetGlobalIndices - Gets the global indices for all the local entries

572:     Collective on Mesh

574:     Input Parameter:
575: .    mesh - the mesh object

577:     Output Parameters:
578: .    idx - the individual indices for each packed vector/array
579:  
580:     Level: advanced

582:     Notes:
583:        The idx parameters should be freed by the calling routine with PetscFree()

585: .seealso MeshDestroy(), MeshCreateGlobalVector(), MeshCreate()

587: @*/
588: PetscErrorCode  MeshGetGlobalIndices(Mesh mesh,PetscInt *idx[])
589: {
590:   SETERRQ(PETSC_ERR_SUP, "");
591: }

595: PetscErrorCode  MeshCreateGlobalScatter(Mesh mesh, SectionReal section, VecScatter *scatter)
596: {
597:   typedef ALE::Mesh::real_section_type::index_type index_type;
598:   ALE::Obj<ALE::Mesh> m;
599:   ALE::Obj<ALE::Mesh::real_section_type> s;
600:   const char         *name;
601:   PetscErrorCode      ierr;

605:   MeshGetMesh(mesh, m);
606:   SectionRealGetSection(section, s);
607:   PetscObjectGetName((PetscObject) section, &name);
608:   const ALE::Obj<ALE::Mesh::topology_type>&                   topology = m->getTopology();
609:   const ALE::Obj<ALE::Mesh::real_section_type::atlas_type>&   atlas    = s->getAtlas();
610:   const ALE::Mesh::real_section_type::patch_type              patch    = 0;
611:   const ALE::Mesh::real_section_type::atlas_type::chart_type& chart    = atlas->getPatch(patch);
612:   const ALE::Obj<ALE::Mesh::order_type>& globalOrder = m->getFactory()->getGlobalOrder(topology, patch, name, atlas);
613:   int *localIndices, *globalIndices;
614:   int  localSize = s->size(patch);
615:   int  localIndx = 0, globalIndx = 0;
616:   Vec  globalVec, localVec;
617:   IS   localIS, globalIS;

619:   VecCreate(m->comm(), &globalVec);
620:   VecSetSizes(globalVec, globalOrder->getLocalSize(), PETSC_DETERMINE);
621:   VecSetFromOptions(globalVec);
622:   // Loop over all local points
623:   PetscMalloc(localSize*sizeof(int), &localIndices);
624:   PetscMalloc(localSize*sizeof(int), &globalIndices);
625:   for(ALE::Mesh::real_section_type::atlas_type::chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
626:     const ALE::Mesh::real_section_type::index_type& idx = atlas->restrictPoint(patch, *p_iter)[0];

628:     // Map local indices to global indices
629:     ExpandInterval(idx, localIndices, localIndx);
630:     ExpandInterval(index_type(idx.prefix, globalOrder->getIndex(*p_iter)), globalIndices, globalIndx);
631:   }
632:   if (localIndx  != localSize) SETERRQ2(PETSC_ERR_ARG_SIZ, "Invalid number of local indices %d, should be %d", localIndx, localSize);
633:   if (globalIndx != localSize) SETERRQ2(PETSC_ERR_ARG_SIZ, "Invalid number of global indices %d, should be %d", globalIndx, localSize);
634:   ISCreateGeneral(PETSC_COMM_SELF, localSize, localIndices,  &localIS);
635:   ISCreateGeneral(PETSC_COMM_SELF, localSize, globalIndices, &globalIS);
636:   PetscFree(localIndices);
637:   PetscFree(globalIndices);
638:   VecCreateSeqWithArray(PETSC_COMM_SELF, localSize, s->restrict(patch), &localVec);
639:   VecScatterCreate(localVec, localIS, globalVec, globalIS, scatter);
640:   ISDestroy(globalIS);
641:   ISDestroy(localIS);
642:   VecDestroy(localVec);
643:   VecDestroy(globalVec);
645:   return(0);
646: }

650: PetscErrorCode  MeshGetGlobalScatter(Mesh mesh, VecScatter *scatter)
651: {

657:   if (!mesh->globalScatter) {
658:     SectionReal section;

660:     MeshGetSectionReal(mesh, "default", &section);
661:     MeshCreateGlobalScatter(mesh, section, &mesh->globalScatter);
662:     SectionRealDestroy(section);
663:   }
664:   *scatter = mesh->globalScatter;
665:   return(0);
666: }

670: PetscErrorCode  MeshSetLocalFunction(Mesh mesh, PetscErrorCode (*lf)(Mesh, SectionReal, SectionReal, void *))
671: {
674:   mesh->lf = lf;
675:   return(0);
676: }

680: PetscErrorCode  MeshSetLocalJacobian(Mesh mesh, PetscErrorCode (*lj)(Mesh, SectionReal, Mat, void *))
681: {
684:   mesh->lj = lj;
685:   return(0);
686: }

690: PetscErrorCode  MeshFormFunction(Mesh mesh, SectionReal X, SectionReal F, void *ctx)
691: {

698:   if (mesh->lf) {
699:     (*mesh->lf)(mesh, X, F, ctx);
700:   }
701:   return(0);
702: }

706: PetscErrorCode  MeshFormJacobian(Mesh mesh, SectionReal X, Mat J, void *ctx)
707: {

714:   if (mesh->lj) {
715:     (*mesh->lj)(mesh, X, J, ctx);
716:   }
717:   return(0);
718: }

722: // Here we assume:
723: //  - Assumes 3D and tetrahedron
724: //  - The section takes values on vertices and is P1
725: //  - Points have the same dimension as the mesh
726: //  - All values have the same dimension
727: PetscErrorCode  MeshInterpolatePoints(Mesh mesh, SectionReal section, int numPoints, double *points, double **values)
728: {
729:   Obj<ALE::Mesh> m;
730:   Obj<ALE::Mesh::real_section_type> s;
731:   double        *v0, *J, *invJ, detJ;

735:   MeshGetMesh(mesh, m);
736:   SectionRealGetSection(section, s);
737:   const Obj<ALE::Mesh::real_section_type>& coordinates = m->getRealSection("coordinates");
738:   int embedDim = coordinates->getFiberDimension(0, *m->getTopology()->depthStratum(0, 0)->begin());
739:   int dim      = s->getFiberDimension(0, *m->getTopology()->depthStratum(0, 0)->begin());

741:   PetscMalloc3(embedDim,double,&v0,embedDim*embedDim,double,&J,embedDim*embedDim,double,&invJ);
742:   PetscMalloc(numPoints*dim * sizeof(double), &values);
743:   for(int p = 0; p < numPoints; p++) {
744:     double *point = &points[p*embedDim];
745:     ALE::Mesh::point_type e = m->locatePoint(0, point);
746:     const ALE::Mesh::real_section_type::value_type *coeff = s->restrict(0, e);

748:     m->computeElementGeometry(coordinates, e, v0, J, invJ, detJ);
749:     double xi   = (invJ[0*embedDim+0]*(point[0] - v0[0]) + invJ[0*embedDim+1]*(point[1] - v0[1]) + invJ[0*embedDim+2]*(point[2] - v0[2]))*0.5;
750:     double eta  = (invJ[1*embedDim+0]*(point[0] - v0[0]) + invJ[1*embedDim+1]*(point[1] - v0[1]) + invJ[1*embedDim+2]*(point[2] - v0[2]))*0.5;
751:     double zeta = (invJ[2*embedDim+0]*(point[0] - v0[0]) + invJ[2*embedDim+1]*(point[1] - v0[1]) + invJ[2*embedDim+2]*(point[2] - v0[2]))*0.5;

753:     for(int d = 0; d < dim; d++) {
754:       (*values)[p*dim+d] = coeff[0*dim+d]*(1 - xi - eta - zeta) + coeff[1*dim+d]*xi + coeff[2*dim+d]*eta + coeff[3*dim+d]*zeta;
755:     }
756:   }
757:   PetscFree3(v0, J, invJ);
758:   return(0);
759: }

761: EXTERN PetscErrorCode assembleFullField(VecScatter, Vec, Vec, InsertMode);

765: /*@C
766:   restrictVector - Insert values from a global vector into a local ghosted vector

768:   Collective on g

770:   Input Parameters:
771: + g - The global vector
772: . l - The local vector
773: - mode - either ADD_VALUES or INSERT_VALUES, where
774:    ADD_VALUES adds values to any existing entries, and
775:    INSERT_VALUES replaces existing entries with new values

777:    Level: beginner

779: .seealso: MatSetOption()
780: @*/
781: PetscErrorCode restrictVector(Vec g, Vec l, InsertMode mode)
782: {
783:   VecScatter     injection;

788:   PetscObjectQuery((PetscObject) g, "injection", (PetscObject *) &injection);
789:   if (injection) {
790:     VecScatterBegin(g, l, mode, SCATTER_REVERSE, injection);
791:     VecScatterEnd(g, l, mode, SCATTER_REVERSE, injection);
792:   } else {
793:     if (mode == INSERT_VALUES) {
794:       VecCopy(g, l);
795:     } else {
796:       VecAXPY(l, 1.0, g);
797:     }
798:   }
800:   return(0);
801: }

805: /*@C
806:   assembleVectorComplete - Insert values from a local ghosted vector into a global vector

808:   Collective on g

810:   Input Parameters:
811: + g - The global vector
812: . l - The local vector
813: - mode - either ADD_VALUES or INSERT_VALUES, where
814:    ADD_VALUES adds values to any existing entries, and
815:    INSERT_VALUES replaces existing entries with new values

817:    Level: beginner

819: .seealso: MatSetOption()
820: @*/
821: PetscErrorCode assembleVectorComplete(Vec g, Vec l, InsertMode mode)
822: {
823:   VecScatter     injection;

828:   PetscObjectQuery((PetscObject) g, "injection", (PetscObject *) &injection);
829:   if (injection) {
830:     VecScatterBegin(l, g, mode, SCATTER_FORWARD, injection);
831:     VecScatterEnd(l, g, mode, SCATTER_FORWARD, injection);
832:   } else {
833:     if (mode == INSERT_VALUES) {
834:       VecCopy(l, g);
835:     } else {
836:       VecAXPY(g, 1.0, l);
837:     }
838:   }
840:   return(0);
841: }

845: /*@C
846:   assembleVector - Insert values into a vector

848:   Collective on A

850:   Input Parameters:
851: + b - the vector
852: . e - The element number
853: . v - The values
854: - mode - either ADD_VALUES or INSERT_VALUES, where
855:    ADD_VALUES adds values to any existing entries, and
856:    INSERT_VALUES replaces existing entries with new values

858:    Level: beginner

860: .seealso: VecSetOption()
861: @*/
862: PetscErrorCode assembleVector(Vec b, PetscInt e, PetscScalar v[], InsertMode mode)
863: {
864:   Mesh                     mesh;
865:   ALE::Obj<ALE::Mesh> m;
866:   ALE::Mesh::real_section_type::patch_type patch;
867:   PetscInt                 firstElement;
868:   PetscErrorCode           ierr;

872:   PetscObjectQuery((PetscObject) b, "mesh", (PetscObject *) &mesh);
873:   MeshGetMesh(mesh, m);
874:   //firstElement = elementBundle->getLocalSizes()[bundle->getCommRank()];
875:   firstElement = 0;
876:   // Must relate b to field
877:   if (mode == INSERT_VALUES) {
878:     m->getRealSection(std::string("x"))->update(patch, ALE::Mesh::point_type(e + firstElement), v);
879:   } else {
880:     m->getRealSection(std::string("x"))->updateAdd(patch, ALE::Mesh::point_type(e + firstElement), v);
881:   }
883:   return(0);
884: }

888: PetscErrorCode updateOperator(Mat A, const ALE::Obj<ALE::Mesh::real_section_type>& section, const ALE::Obj<ALE::Mesh::order_type>& globalOrder, const ALE::Mesh::point_type& e, PetscScalar array[], InsertMode mode)
889: {
890:   ALE::Mesh::real_section_type::patch_type patch = 0;
891:   static PetscInt  indicesSize = 0;
892:   static PetscInt *indices = NULL;
893:   PetscInt         numIndices = 0;
894:   PetscErrorCode   ierr;

897:   const ALE::Obj<ALE::Mesh::real_section_type::IndexArray> intervals = section->getIndices(patch, e, globalOrder);

900:   if (section->debug()) {printf("[%d]mat for element %d\n", section->commRank(), e);}
901:   for(ALE::Mesh::real_section_type::IndexArray::iterator i_iter = intervals->begin(); i_iter != intervals->end(); ++i_iter) {
902:     numIndices += std::abs(i_iter->prefix);
903:     if (section->debug()) {
904:       printf("[%d]mat interval (%d, %d)\n", section->commRank(), i_iter->prefix, i_iter->index);
905:     }
906:   }
907:   if (indicesSize && (indicesSize != numIndices)) {
908:     PetscFree(indices);
909:     indices = NULL;
910:   }
911:   if (!indices) {
912:     indicesSize = numIndices;
913:     PetscMalloc(indicesSize * sizeof(PetscInt), &indices);
914:   }
915:   ExpandIntervals(intervals, indices);
916:   if (section->debug()) {
917:     for(int i = 0; i < numIndices; i++) {
918:       printf("[%d]mat indices[%d] = %d\n", section->commRank(), i, indices[i]);
919:     }
920:     for(int i = 0; i < numIndices; i++) {
921:       printf("[%d]", section->commRank());
922:       for(int j = 0; j < numIndices; j++) {
923:         printf(" %g", array[i*numIndices+j]);
924:       }
925:       printf("\n");
926:     }
927:   }
928:   MatSetValues(A, numIndices, indices, numIndices, indices, array, mode);
929:   if (ierr) {
930:     printf("[%d]ERROR in updateOperator: point %d\n", section->commRank(), e);
931:     for(int i = 0; i < numIndices; i++) {
932:       printf("[%d]mat indices[%d] = %d\n", section->commRank(), i, indices[i]);
933:     }
934: 
935:   }
937:   return(0);
938: }

942: PetscErrorCode updateOperatorGeneral(Mat A, const ALE::Obj<ALE::Mesh::real_section_type>& rowSection, const ALE::Obj<ALE::Mesh::order_type>& rowGlobalOrder, const ALE::Mesh::point_type& rowE, const ALE::Obj<ALE::Mesh::real_section_type>& colSection, const ALE::Obj<ALE::Mesh::order_type>& colGlobalOrder, const ALE::Mesh::point_type& colE, PetscScalar array[], InsertMode mode)
943: {
944:   ALE::Mesh::real_section_type::patch_type patch = 0;
945:   static PetscInt  rowIndicesSize = 0;
946:   static PetscInt *rowIndices = NULL;
947:   PetscInt         numRowIndices = 0;
948:   static PetscInt  colIndicesSize = 0;
949:   static PetscInt *colIndices = NULL;
950:   PetscInt         numColIndices = 0;
951:   PetscErrorCode   ierr;

954:   const ALE::Obj<ALE::Mesh::real_section_type::IndexArray> rowIntervals = rowSection->getIndices(patch, rowE, rowGlobalOrder);
955:   const ALE::Obj<ALE::Mesh::real_section_type::IndexArray> colIntervals = colSection->getIndices(patch, colE, colGlobalOrder);

958:   if (rowSection->debug()) {printf("[%d]mat for elements %d %d\n", rowSection->commRank(), rowE, colE);}
959:   for(ALE::Mesh::real_section_type::IndexArray::iterator i_iter = rowIntervals->begin(); i_iter != rowIntervals->end(); ++i_iter) {
960:     numRowIndices += std::abs(i_iter->prefix);
961:     if (rowSection->debug()) {
962:       printf("[%d]mat row interval (%d, %d)\n", rowSection->commRank(), i_iter->prefix, i_iter->index);
963:     }
964:   }
965:   if (rowIndicesSize && (rowIndicesSize != numRowIndices)) {
966:     PetscFree(rowIndices);
967:     rowIndices = NULL;
968:   }
969:   if (!rowIndices) {
970:     rowIndicesSize = numRowIndices;
971:     PetscMalloc(rowIndicesSize * sizeof(PetscInt), &rowIndices);
972:   }
973:   ExpandIntervals(rowIntervals, rowIndices);
974:   if (rowSection->debug()) {
975:     for(int i = 0; i < numRowIndices; i++) {
976:       printf("[%d]mat row indices[%d] = %d\n", rowSection->commRank(), i, rowIndices[i]);
977:     }
978:   }
979:   for(ALE::Mesh::real_section_type::IndexArray::iterator i_iter = colIntervals->begin(); i_iter != colIntervals->end(); ++i_iter) {
980:     numColIndices += std::abs(i_iter->prefix);
981:     if (colSection->debug()) {
982:       printf("[%d]mat col interval (%d, %d)\n", colSection->commRank(), i_iter->prefix, i_iter->index);
983:     }
984:   }
985:   if (colIndicesSize && (colIndicesSize != numColIndices)) {
986:     PetscFree(colIndices);
987:     colIndices = NULL;
988:   }
989:   if (!colIndices) {
990:     colIndicesSize = numColIndices;
991:     PetscMalloc(colIndicesSize * sizeof(PetscInt), &colIndices);
992:   }
993:   ExpandIntervals(colIntervals, colIndices);
994:   if (colSection->debug()) {
995:     for(int i = 0; i < numColIndices; i++) {
996:       printf("[%d]mat col indices[%d] = %d\n", colSection->commRank(), i, colIndices[i]);
997:     }
998:     for(int i = 0; i < numRowIndices; i++) {
999:       printf("[%d]", rowSection->commRank());
1000:       for(int j = 0; j < numColIndices; j++) {
1001:         printf(" %g", array[i*numColIndices+j]);
1002:       }
1003:       printf("\n");
1004:     }
1005:   }
1006:   MatSetValues(A, numRowIndices, rowIndices, numColIndices, colIndices, array, mode);
1007:   if (ierr) {
1008:     printf("[%d]ERROR in updateOperator: points %d %d\n", colSection->commRank(), rowE, colE);
1009:     for(int i = 0; i < numRowIndices; i++) {
1010:       printf("[%d]mat row indices[%d] = %d\n", rowSection->commRank(), i, rowIndices[i]);
1011:     }
1012:     for(int i = 0; i < numColIndices; i++) {
1013:       printf("[%d]mat col indices[%d] = %d\n", colSection->commRank(), i, colIndices[i]);
1014:     }
1015: 
1016:   }
1018:   return(0);
1019: }

1023: /*@C
1024:   assembleMatrix - Insert values into a matrix

1026:   Collective on A

1028:   Input Parameters:
1029: + A - the matrix
1030: . e - The element number
1031: . v - The values
1032: - mode - either ADD_VALUES or INSERT_VALUES, where
1033:    ADD_VALUES adds values to any existing entries, and
1034:    INSERT_VALUES replaces existing entries with new values

1036:    Level: beginner

1038: .seealso: MatSetOption()
1039: @*/
1040: PetscErrorCode assembleMatrix(Mat A, PetscInt e, PetscScalar v[], InsertMode mode)
1041: {
1042:   Mesh           mesh;
1043:   Obj<ALE::Mesh> m;

1048:   PetscObjectQuery((PetscObject) A, "mesh", (PetscObject *) &mesh);
1049:   MeshGetMesh(mesh, m);
1050:   try {
1051:     const ALE::Obj<ALE::Mesh::topology_type>&     topology    = m->getTopology();
1052:     const ALE::Obj<ALE::Mesh::numbering_type>&    cNumbering  = m->getFactory()->getLocalNumbering(topology, 0, topology->depth());
1053:     const ALE::Obj<ALE::Mesh::real_section_type>& section     = m->getRealSection("default");
1054:     const ALE::Obj<ALE::Mesh::order_type>&        globalOrder = m->getFactory()->getGlobalOrder(topology, 0, "default", section->getAtlas());

1056:     if (section->debug()) {
1057:       std::cout << "Assembling matrix for element number " << e << " --> point " << cNumbering->getPoint(e) << std::endl;
1058:     }
1059:     updateOperator(A, section, globalOrder, cNumbering->getPoint(e), v, mode);
1060:   } catch (ALE::Exception e) {
1061:     std::cout << e.msg() << std::endl;
1062:   }
1064:   return(0);
1065: }

1069: template<typename Atlas>
1070: PetscErrorCode preallocateOperator(const ALE::Obj<ALE::Mesh::topology_type>& topology, const ALE::Obj<Atlas>& atlas, const ALE::Obj<ALE::Mesh::order_type>& globalOrder, Mat A)
1071: {
1072:   typedef ALE::New::NumberingFactory<ALE::Mesh::topology_type> NumberingFactory;
1073:   const ALE::Obj<ALE::Mesh::sieve_type>     adjGraph    = new ALE::Mesh::sieve_type(topology->comm(), topology->debug());
1074:   const ALE::Obj<ALE::Mesh::topology_type>  adjTopology = new ALE::Mesh::topology_type(topology->comm(), topology->debug());
1075:   const ALE::Mesh::real_section_type::patch_type patch  = 0;
1076:   const ALE::Obj<ALE::Mesh::sieve_type>&    sieve       = topology->getPatch(patch);
1077:   PetscInt       numLocalRows, firstRow;
1078:   PetscInt      *dnz, *onz;

1082:   adjTopology->setPatch(patch, adjGraph);
1083:   numLocalRows = globalOrder->getLocalSize();
1084:   firstRow     = globalOrder->getGlobalOffsets()[topology->commRank()];
1085:   PetscMalloc2(numLocalRows, PetscInt, &dnz, numLocalRows, PetscInt, &onz);
1086:   /* Create local adjacency graph */
1087:   /*   In general, we need to get FIAT info that attaches dual basis vectors to sieve points */
1088:   const ALE::Mesh::real_section_type::atlas_type::chart_type& chart = atlas->getPatch(patch);

1090:   for(ALE::Mesh::real_section_type::atlas_type::chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
1091:     const ALE::Mesh::real_section_type::atlas_type::point_type& point = *c_iter;

1093:     adjGraph->addCone(sieve->cone(sieve->support(point)), point);
1094:   }
1095:   /* Distribute adjacency graph */
1096:   topology->constructOverlap(patch);
1097:   const Obj<ALE::Mesh::send_overlap_type>& vertexSendOverlap = topology->getSendOverlap();
1098:   const Obj<ALE::Mesh::recv_overlap_type>& vertexRecvOverlap = topology->getRecvOverlap();
1099:   const Obj<ALE::Mesh::send_overlap_type>  nbrSendOverlap    = new ALE::Mesh::send_overlap_type(topology->comm(), topology->debug());
1100:   const Obj<ALE::Mesh::recv_overlap_type>  nbrRecvOverlap    = new ALE::Mesh::recv_overlap_type(topology->comm(), topology->debug());
1101:   const Obj<ALE::Mesh::send_section_type>  sendSection       = new ALE::Mesh::send_section_type(topology->comm(), topology->debug());
1102:   const Obj<ALE::Mesh::recv_section_type>  recvSection       = new ALE::Mesh::recv_section_type(topology->comm(), sendSection->getTag(), topology->debug());

1104:   ALE::New::Distribution<ALE::Mesh::topology_type>::coneCompletion(vertexSendOverlap, vertexRecvOverlap, adjTopology, sendSection, recvSection);
1105:   /* Distribute indices for new points */
1106:   ALE::New::Distribution<ALE::Mesh::topology_type>::updateOverlap(sendSection, recvSection, nbrSendOverlap, nbrRecvOverlap);
1107:   NumberingFactory::singleton(topology->debug())->completeOrder(globalOrder, nbrSendOverlap, nbrRecvOverlap, patch, true);
1108:   /* Read out adjacency graph */
1109:   const ALE::Obj<ALE::Mesh::sieve_type> graph = adjTopology->getPatch(patch);

1111:   PetscMemzero(dnz, numLocalRows * sizeof(PetscInt));
1112:   PetscMemzero(onz, numLocalRows * sizeof(PetscInt));
1113:   for(ALE::Mesh::real_section_type::atlas_type::chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
1114:     const ALE::Mesh::real_section_type::atlas_type::point_type& point = *c_iter;

1116:     if (globalOrder->isLocal(point)) {
1117:       const Obj<ALE::Mesh::sieve_type::traits::coneSequence>& adj   = graph->cone(point);
1118:       const ALE::Mesh::order_type::value_type&                rIdx  = globalOrder->restrictPoint(patch, point)[0];
1119:       const int                                               row   = rIdx.prefix;
1120:       const int                                               rSize = rIdx.index;

1122:       for(ALE::Mesh::sieve_type::traits::coneSequence::iterator v_iter = adj->begin(); v_iter != adj->end(); ++v_iter) {
1123:         const ALE::Mesh::real_section_type::atlas_type::point_type& neighbor = *v_iter;
1124:         const ALE::Mesh::order_type::value_type& cIdx     = globalOrder->restrictPoint(patch, neighbor)[0];
1125:         const int&                               cSize    = cIdx.index;

1127:         if (cSize > 0) {
1128:           if (globalOrder->isLocal(neighbor)) {
1129:             for(int r = 0; r < rSize; ++r) {dnz[row - firstRow + r] += cSize;}
1130:           } else {
1131:             for(int r = 0; r < rSize; ++r) {onz[row - firstRow + r] += cSize;}
1132:           }
1133:         }
1134:       }
1135:     }
1136:   }
1137:   if (topology->debug()) {
1138:     int rank = topology->commRank();
1139:     for(int r = 0; r < numLocalRows; r++) {
1140:       std::cout << "["<<rank<<"]: dnz["<<r<<"]: " << dnz[r] << " onz["<<r<<"]: " << onz[r] << std::endl;
1141:     }
1142:   }
1143:   MatSeqAIJSetPreallocation(A, 0, dnz);
1144:   MatMPIAIJSetPreallocation(A, 0, dnz, 0, onz);
1145:   PetscFree2(dnz, onz);
1146:   MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR);
1147:   return(0);
1148: }

1152: PetscErrorCode preallocateMatrix(const ALE::Obj<ALE::Mesh::topology_type>& topology, const ALE::Obj<ALE::Mesh::real_section_type::atlas_type>& atlas, const ALE::Obj<ALE::Mesh::order_type>& globalOrder, Mat A)
1153: {
1154:   return preallocateOperator(topology, atlas, globalOrder, A);
1155: }

1157: /******************************** C Wrappers **********************************/

1161: PetscErrorCode WriteVTKHeader(PetscViewer viewer)
1162: {
1163:   return VTKViewer::writeHeader(viewer);
1164: }

1168: PetscErrorCode WriteVTKVertices(Mesh mesh, PetscViewer viewer)
1169: {
1170:   ALE::Obj<ALE::Mesh> m;

1173:   MeshGetMesh(mesh, m);
1174:   return VTKViewer::writeVertices(m, 0, viewer);
1175: }

1179: PetscErrorCode WriteVTKElements(Mesh mesh, PetscViewer viewer)
1180: {
1181:   ALE::Obj<ALE::Mesh> m;

1184:   MeshGetMesh(mesh, m);
1185:   return VTKViewer::writeElements(m, 0, viewer);
1186: }

1190: PetscErrorCode WritePCICEVertices(Mesh mesh, PetscViewer viewer)
1191: {
1192:   ALE::Obj<ALE::Mesh> m;

1195:   MeshGetMesh(mesh, m);
1196:   return ALE::PCICE::Viewer::writeVertices(m, viewer);
1197: }

1201: PetscErrorCode WritePCICEElements(Mesh mesh, PetscViewer viewer)
1202: {
1203:   ALE::Obj<ALE::Mesh> m;

1206:   MeshGetMesh(mesh, m);
1207:   return ALE::PCICE::Viewer::writeElements(m, viewer);
1208: }

1212: PetscErrorCode WritePCICERestart(Mesh mesh, PetscViewer viewer)
1213: {
1214:   ALE::Obj<ALE::Mesh> m;

1217:   MeshGetMesh(mesh, m);
1218:   return ALE::PCICE::Viewer::writeRestart(m, viewer);
1219: }

1223: /*@C
1224:   MeshCreatePCICE - Create a Mesh from PCICE files.

1226:   Not Collective

1228:   Input Parameters:
1229: + dim - The topological mesh dimension
1230: . coordFilename - The file containing vertex coordinates
1231: . adjFilename - The file containing the vertices for each element
1232: . bcFilename - The file containing the boundary topology and conditions
1233: . numBdFaces - The number of boundary faces (or edges)
1234: - numBdVertices - The number of boundary vertices

1236:   Output Parameter:
1237: . mesh - The Mesh object

1239:   Level: beginner

1241: .keywords: mesh, PCICE
1242: .seealso: MeshCreate()
1243: @*/
1244: PetscErrorCode MeshCreatePCICE(MPI_Comm comm, const int dim, const char coordFilename[], const char adjFilename[], const char bcFilename[], const int numBdFaces, const int numBdVertices, Mesh *mesh)
1245: {
1246:   ALE::Obj<ALE::Mesh> m;
1247:   PetscInt            debug = 0;
1248:   PetscTruth          flag;
1249:   PetscErrorCode      ierr;

1252:   MeshCreate(comm, mesh);
1253:   PetscOptionsGetInt(PETSC_NULL, "-debug", &debug, &flag);
1254:   try {
1255:     m  = ALE::PCICE::Builder::readMesh(comm, dim, std::string(coordFilename), std::string(adjFilename), false, false, debug);
1256:   } catch(ALE::Exception e) {
1257:     SETERRQ(PETSC_ERR_FILE_OPEN, e.message());
1258:   }
1259:   if (bcFilename) {
1260:     ALE::PCICE::Builder::readBoundary(m, std::string(bcFilename), numBdFaces, numBdVertices);
1261:   }
1262:   MeshSetMesh(*mesh, m);
1263:   return(0);
1264: }

1268: /*@C
1269:   MeshCreatePyLith - Create a Mesh from PyLith files.

1271:   Not Collective

1273:   Input Parameters:
1274: + dim - The topological mesh dimension
1275: . baseFilename - The basename for mesh files
1276: . zeroBase - Use 0 to start numbering
1277: - interpolate - The flag for mesh interpolation

1279:   Output Parameter:
1280: . mesh - The Mesh object

1282:   Level: beginner

1284: .keywords: mesh, PCICE
1285: .seealso: MeshCreate()
1286: @*/
1287: PetscErrorCode MeshCreatePyLith(MPI_Comm comm, const int dim, const char baseFilename[], PetscTruth zeroBase, PetscTruth interpolate, Mesh *mesh)
1288: {
1289:   ALE::Obj<ALE::Mesh> m;
1290:   PetscInt            debug = 0;
1291:   PetscTruth          flag;
1292:   PetscErrorCode      ierr;

1295:   MeshCreate(comm, mesh);
1296:   PetscOptionsGetInt(PETSC_NULL, "-debug", &debug, &flag);
1297:   try {
1298:     m  = ALE::PyLith::Builder::readMesh(comm, dim, std::string(baseFilename), zeroBase, interpolate, debug);
1299:   } catch(ALE::Exception e) {
1300:     SETERRQ(PETSC_ERR_FILE_OPEN, e.message());
1301:   }
1302:   MeshSetMesh(*mesh, m);
1303:   return(0);
1304: }

1308: /*@C
1309:   MeshGetCoordinates - Creates an array holding the coordinates.

1311:   Not Collective

1313:   Input Parameter:
1314: + mesh - The Mesh object
1315: - columnMajor - Flag for column major order

1317:   Output Parameter:
1318: + numVertices - The number of vertices
1319: . dim - The embedding dimension
1320: - coords - The array holding local coordinates

1322:   Level: intermediate

1324: .keywords: mesh, coordinates
1325: .seealso: MeshCreate()
1326: @*/
1327: PetscErrorCode MeshGetCoordinates(Mesh mesh, PetscTruth columnMajor, PetscInt *numVertices, PetscInt *dim, PetscReal *coords[])
1328: {
1329:   ALE::Obj<ALE::Mesh> m;
1330:   PetscErrorCode      ierr;

1333:   MeshGetMesh(mesh, m);
1334:   ALE::PCICE::Builder::outputVerticesLocal(m, numVertices, dim, coords, columnMajor);
1335:   return(0);
1336: }

1340: /*@C
1341:   MeshGetElements - Creates an array holding the vertices on each element.

1343:   Not Collective

1345:   Input Parameters:
1346: + mesh - The Mesh object
1347: - columnMajor - Flag for column major order

1349:   Output Parameters:
1350: + numElements - The number of elements
1351: . numCorners - The number of vertices per element
1352: - vertices - The array holding vertices on each local element

1354:   Level: intermediate

1356: .keywords: mesh, elements
1357: .seealso: MeshCreate()
1358: @*/
1359: PetscErrorCode MeshGetElements(Mesh mesh, PetscTruth columnMajor, PetscInt *numElements, PetscInt *numCorners, PetscInt *vertices[])
1360: {
1361:   ALE::Obj<ALE::Mesh> m;
1362:   PetscErrorCode      ierr;

1365:   MeshGetMesh(mesh, m);
1366:   ALE::PCICE::Builder::outputElementsLocal(m, numElements, numCorners, vertices, columnMajor);
1367:   return(0);
1368: }

1372: /*@C
1373:   MeshDistribute - Distributes the mesh and any associated sections.

1375:   Not Collective

1377:   Input Parameter:
1378: . serialMesh - The original Mesh object

1380:   Output Parameter:
1381: . parallelMesh - The distributed Mesh object

1383:   Level: intermediate

1385: .keywords: mesh, elements
1386: .seealso: MeshCreate()
1387: @*/
1388: PetscErrorCode MeshDistribute(Mesh serialMesh, Mesh *parallelMesh)
1389: {
1390:   ALE::Obj<ALE::Mesh> oldMesh;
1391:   PetscErrorCode      ierr;

1394:   MeshGetMesh(serialMesh, oldMesh);
1395:   MeshCreate(oldMesh->comm(), parallelMesh);
1396:   ALE::Obj<ALE::Mesh> newMesh = ALE::New::Distribution<ALE::Mesh::topology_type>::distributeMesh(oldMesh);
1397:   MeshSetMesh(*parallelMesh, newMesh);
1398:   return(0);
1399: }

1403: /*@C
1404:   MeshGenerate - Generates a mesh.

1406:   Not Collective

1408:   Input Parameters:
1409: + boundary - The Mesh boundary object
1410: - interpolate - Flag to create intermediate mesh elements

1412:   Output Parameter:
1413: . mesh - The Mesh object

1415:   Level: intermediate

1417: .keywords: mesh, elements
1418: .seealso: MeshCreate(), MeshRefine()
1419: @*/
1420: PetscErrorCode MeshGenerate(Mesh boundary, PetscTruth interpolate, Mesh *mesh)
1421: {
1422:   ALE::Obj<ALE::Mesh> mB;
1423:   PetscErrorCode      ierr;

1426:   MeshGetMesh(boundary, mB);
1427:   MeshCreate(mB->comm(), mesh);
1428:   ALE::Obj<ALE::Mesh> m = ALE::Generator::generateMesh(mB, interpolate);
1429:   MeshSetMesh(*mesh, m);
1430:   return(0);
1431: }

1435: /*@C
1436:   MeshRefine - Refines the mesh.

1438:   Not Collective

1440:   Input Parameters:
1441: + mesh - The original Mesh object
1442: . refinementLimit - The maximum size of any cell
1443: - interpolate - Flag to create intermediate mesh elements

1445:   Output Parameter:
1446: . refinedMesh - The refined Mesh object

1448:   Level: intermediate

1450: .keywords: mesh, elements
1451: .seealso: MeshCreate(), MeshGenerate()
1452: @*/
1453: PetscErrorCode MeshRefine(Mesh mesh, double refinementLimit, PetscTruth interpolate, Mesh *refinedMesh)
1454: {
1455:   ALE::Obj<ALE::Mesh> oldMesh;
1456:   PetscErrorCode      ierr;

1459:   MeshGetMesh(mesh, oldMesh);
1460:   MeshCreate(oldMesh->comm(), refinedMesh);
1461:   ALE::Obj<ALE::Mesh> newMesh = ALE::Generator::refineMesh(oldMesh, refinementLimit, interpolate);
1462:   MeshSetMesh(*refinedMesh, newMesh);
1463:   return(0);
1464: }

1468: PetscErrorCode MeshRefine_DM(Mesh mesh, MPI_Comm comm, Mesh *refinedMesh)
1469: {
1470:   ALE::Obj<ALE::Mesh> oldMesh;
1471:   double              refinementLimit;
1472:   PetscErrorCode      ierr;

1475:   MeshGetMesh(mesh, oldMesh);
1476:   MeshCreate(comm, refinedMesh);
1477:   refinementLimit = oldMesh->getMaxVolume()/2.0;
1478:   ALE::Obj<ALE::Mesh> newMesh = ALE::Generator::refineMesh(oldMesh, refinementLimit, false);
1479:   MeshSetMesh(*refinedMesh, newMesh);
1480:   const ALE::Obj<ALE::Mesh::real_section_type>& s = newMesh->getRealSection("default");

1482:   newMesh->setDiscretization(oldMesh->getDiscretization());
1483:   newMesh->setBoundaryCondition(oldMesh->getBoundaryCondition());
1484:   newMesh->setupField(s);
1485:   return(0);
1486: }

1488: #include "src/dm/mesh/examples/tutorials/Coarsener.h"
1489: #include "src/dm/mesh/examples/tutorials/fast_coarsen.h"

1493: /*@C
1494:   MeshCoarsenHierarchy - Coarsens the mesh into a hierarchy.

1496:   Not Collective

1498:   Input Parameters:
1499: + mesh - The original Mesh object
1500: . numLevels - The number of 
1501: . coarseningFactor - The expansion factor for coarse meshes
1502: - interpolate - Flag to create intermediate mesh elements

1504:   Output Parameter:
1505: . coarseHierarchy - The coarse Mesh objects

1507:   Level: intermediate

1509: .keywords: mesh, elements
1510: .seealso: MeshCreate(), MeshGenerate()
1511: @*/
1512: PetscErrorCode MeshCoarsenHierarchy(Mesh mesh, int numLevels, double coarseningFactor, PetscTruth interpolate, Mesh **coarseHierarchy)
1513: {
1514:   ALE::Obj<ALE::Mesh> oldMesh;
1515:   PetscErrorCode      ierr;

1518:   if (numLevels < 1) {
1519:     *coarseHierarchy = PETSC_NULL;
1520:     return(0);
1521:   }
1522:   MeshGetMesh(mesh, oldMesh);
1523:   if (oldMesh->getDimension() != 2) SETERRQ(PETSC_ERR_SUP, "Coarsening only works in two dimensions right now");
1524:   ALE::Coarsener::IdentifyBoundary(oldMesh, 2);
1525:   ALE::Coarsener::make_coarsest_boundary(oldMesh, 2, numLevels+1);
1526:   ALE::Coarsener::CreateSpacingFunction(oldMesh, 2);
1527:   ALE::Coarsener::CreateCoarsenedHierarchyNew(oldMesh, 2, numLevels, coarseningFactor);
1528:   PetscMalloc(numLevels * sizeof(Mesh),coarseHierarchy);
1529:   for(int l = 0; l < numLevels; l++) {
1530:     ALE::Obj<ALE::Mesh> newMesh = new ALE::Mesh(oldMesh->comm(), oldMesh->debug());
1531:     const ALE::Obj<ALE::Mesh::real_section_type>& s = newMesh->getRealSection("default");

1533:     MeshCreate(oldMesh->comm(), &(*coarseHierarchy)[l]);
1534:     newMesh->getTopology()->setPatch(0, oldMesh->getTopology()->getPatch(l+1));
1535:     newMesh->setDiscretization(oldMesh->getDiscretization());
1536:     newMesh->setBoundaryCondition(oldMesh->getBoundaryCondition());
1537:     newMesh->setupField(s);
1538:     MeshSetMesh((*coarseHierarchy)[l], newMesh);
1539:   }
1540:   return(0);
1541: }

1543: PetscErrorCode MeshCoarsenHierarchy_DM(Mesh mesh, int numLevels, Mesh **coarseHierarchy)
1544: {

1548:   MeshCoarsenHierarchy(mesh, numLevels, 1.41, PETSC_FALSE, coarseHierarchy);
1549:   return(0);
1550: }

1554: /*
1555:   This method only handle P_1 discretizations at present.
1556: */
1557: PetscErrorCode MeshGetInterpolation_DM(DM dmFine, DM dmCoarse, Mat *interpolation, Vec *scaling)
1558: {
1559:   ALE::Obj<ALE::Mesh> coarse;
1560:   ALE::Obj<ALE::Mesh> fine;
1561:   Mesh coarseMesh = (Mesh) dmCoarse;
1562:   Mesh fineMesh   = (Mesh) dmFine;
1563:   const ALE::Mesh::patch_type patch = 0;
1564:   Mat  P;

1568:   MeshGetMesh(fineMesh,   fine);
1569:   MeshGetMesh(coarseMesh, coarse);
1570:   const ALE::Obj<ALE::Mesh::real_section_type>&             coarseCoordinates = coarse->getRealSection("coordinates");
1571:   const ALE::Obj<ALE::Mesh::real_section_type>&             fineCoordinates   = fine->getRealSection("coordinates");
1572:   const ALE::Obj<ALE::Mesh::topology_type::label_sequence>& vertices          = fine->getTopology()->depthStratum(patch, 0);
1573:   const ALE::Obj<ALE::Mesh::real_section_type>&             sCoarse           = coarse->getRealSection("default");
1574:   const ALE::Obj<ALE::Mesh::real_section_type>&             sFine             = fine->getRealSection("default");
1575:   const ALE::Obj<ALE::Mesh::real_section_type::atlas_type>& coarseAtlas       = sCoarse->getAtlas();
1576:   const ALE::Obj<ALE::Mesh::real_section_type::atlas_type>& fineAtlas         = sFine->getAtlas();
1577:   const ALE::Obj<ALE::Mesh::order_type>& coarseOrder = coarse->getFactory()->getGlobalOrder(coarse->getTopology(), patch, "default", coarseAtlas);
1578:   const ALE::Obj<ALE::Mesh::order_type>& fineOrder   = fine->getFactory()->getGlobalOrder(fine->getTopology(), patch, "default", fineAtlas);
1579:   const int dim = coarse->getDimension();
1580:   double *v0, *J, *invJ, detJ, *refCoords, *values;

1582:   MatCreate(fine->comm(), &P);
1583:   MatSetSizes(P, sFine->size(patch), sCoarse->size(patch), PETSC_DETERMINE, PETSC_DETERMINE);
1584:   MatSetFromOptions(P);
1585:   PetscMalloc5(dim,double,&v0,dim*dim,double,&J,dim*dim,double,&invJ,dim,double,&refCoords,dim+1,double,&values);
1586:   for(ALE::Mesh::topology_type::label_sequence::iterator v_iter = vertices->begin(); v_iter != vertices->end(); ++v_iter) {
1587:     const ALE::Mesh::real_section_type::value_type *coords     = fineCoordinates->restrict(patch, *v_iter);
1588:     const ALE::Mesh::point_type                     coarseCell = coarse->locatePoint(patch, coords);

1590:     coarse->computeElementGeometry(coarseCoordinates, coarseCell, v0, J, invJ, detJ);
1591:     for(int d = 0; d < dim; ++d) {
1592:       refCoords[d] = 0.0;
1593:       for(int e = 0; e < dim; ++e) {
1594:         refCoords[d] += invJ[d*dim+e]*(coords[e] - v0[e]);
1595:       }
1596:       refCoords[d] -= 1.0;
1597:     }
1598:     values[0] = 1.0/3.0 - (refCoords[0] + refCoords[1])/3.0;
1599:     values[1] = 0.5*(refCoords[0] + 1.0);
1600:     values[2] = 0.5*(refCoords[1] + 1.0);
1601:     updateOperatorGeneral(P, sFine, fineOrder, *v_iter, sCoarse, coarseOrder, coarseCell, values, INSERT_VALUES);
1602:   }
1603:   PetscFree5(v0,J,invJ,refCoords,values);
1604:   MatAssemblyBegin(P, MAT_FINAL_ASSEMBLY);
1605:   MatAssemblyEnd(P, MAT_FINAL_ASSEMBLY);
1606:   *interpolation = P;
1607:   return(0);
1608: }

1612: /*@C
1613:   MeshHasSectionReal - Determines whether this mesh has a SectionReal with the given name.

1615:   Not Collective

1617:   Input Parameters:
1618: + mesh - The Mesh object
1619: - name - The section name

1621:   Output Parameter:
1622: . flag - True if the SectionReal is present in the Mesh

1624:   Level: intermediate

1626: .keywords: mesh, elements
1627: .seealso: MeshCreate()
1628: @*/
1629: PetscErrorCode MeshHasSectionReal(Mesh mesh, const char name[], PetscTruth *flag)
1630: {
1631:   ALE::Obj<ALE::Mesh> m;
1632:   PetscErrorCode      ierr;

1635:   MeshGetMesh(mesh, m);
1636:   *flag = (PetscTruth) m->hasRealSection(std::string(name));
1637:   return(0);
1638: }

1642: /*@C
1643:   MeshGetSectionReal - Returns a SectionReal of the given name from the Mesh.

1645:   Collective on Mesh

1647:   Input Parameters:
1648: + mesh - The Mesh object
1649: - name - The section name

1651:   Output Parameter:
1652: . section - The SectionReal

1654:   Note: The section is a new object, and must be destroyed by the user

1656:   Level: intermediate

1658: .keywords: mesh, elements
1659: .seealso: MeshCreate()
1660: @*/
1661: PetscErrorCode MeshGetSectionReal(Mesh mesh, const char name[], SectionReal *section)
1662: {
1663:   ALE::Obj<ALE::Mesh> m;
1664:   PetscErrorCode      ierr;

1667:   MeshGetMesh(mesh, m);
1668:   SectionRealCreate(m->comm(), section);
1669:   PetscObjectSetName((PetscObject) *section, name);
1670:   SectionRealSetSection(*section, m->getRealSection(std::string(name)));
1671:   return(0);
1672: }

1676: /*@C
1677:   MeshSetSectionReal - Puts a SectionReal of the given name into the Mesh.

1679:   Collective on Mesh

1681:   Input Parameters:
1682: + mesh - The Mesh object
1683: - section - The SectionReal

1685:   Note: This takes the section name from the PETSc object

1687:   Level: intermediate

1689: .keywords: mesh, elements
1690: .seealso: MeshCreate()
1691: @*/
1692: PetscErrorCode MeshSetSectionReal(Mesh mesh, SectionReal section)
1693: {
1694:   ALE::Obj<ALE::Mesh> m;
1695:   ALE::Obj<ALE::Mesh::real_section_type> s;
1696:   const char         *name;
1697:   PetscErrorCode      ierr;

1700:   MeshGetMesh(mesh, m);
1701:   PetscObjectGetName((PetscObject) section, &name);
1702:   SectionRealGetSection(section, s);
1703:   m->setRealSection(std::string(name), s);
1704:   return(0);
1705: }

1709: /*@C
1710:   MeshHasSectionInt - Determines whether this mesh has a SectionInt with the given name.

1712:   Not Collective

1714:   Input Parameters:
1715: + mesh - The Mesh object
1716: - name - The section name

1718:   Output Parameter:
1719: . flag - True if the SectionInt is present in the Mesh

1721:   Level: intermediate

1723: .keywords: mesh, elements
1724: .seealso: MeshCreate()
1725: @*/
1726: PetscErrorCode MeshHasSectionInt(Mesh mesh, const char name[], PetscTruth *flag)
1727: {
1728:   ALE::Obj<ALE::Mesh> m;
1729:   PetscErrorCode      ierr;

1732:   MeshGetMesh(mesh, m);
1733:   *flag = (PetscTruth) m->hasIntSection(std::string(name));
1734:   return(0);
1735: }

1739: /*@C
1740:   MeshGetSectionInt - Returns a SectionInt of the given name from the Mesh.

1742:   Collective on Mesh

1744:   Input Parameters:
1745: + mesh - The Mesh object
1746: - name - The section name

1748:   Output Parameter:
1749: . section - The SectionInt

1751:   Note: The section is a new object, and must be destroyed by the user

1753:   Level: intermediate

1755: .keywords: mesh, elements
1756: .seealso: MeshCreate()
1757: @*/
1758: PetscErrorCode MeshGetSectionInt(Mesh mesh, const char name[], SectionInt *section)
1759: {
1760:   ALE::Obj<ALE::Mesh> m;
1761:   PetscErrorCode      ierr;

1764:   MeshGetMesh(mesh, m);
1765:   SectionIntCreate(m->comm(), section);
1766:   PetscObjectSetName((PetscObject) *section, name);
1767:   SectionIntSetSection(*section, m->getIntSection(std::string(name)));
1768:   return(0);
1769: }

1773: /*@C
1774:   MeshSetSectionInt - Puts a SectionInt of the given name into the Mesh.

1776:   Collective on Mesh

1778:   Input Parameters:
1779: + mesh - The Mesh object
1780: - section - The SectionInt

1782:   Note: This takes the section name from the PETSc object

1784:   Level: intermediate

1786: .keywords: mesh, elements
1787: .seealso: MeshCreate()
1788: @*/
1789: PetscErrorCode MeshSetSectionInt(Mesh mesh, SectionInt section)
1790: {
1791:   ALE::Obj<ALE::Mesh> m;
1792:   ALE::Obj<ALE::Mesh::int_section_type> s;
1793:   const char         *name;
1794:   PetscErrorCode      ierr;

1797:   MeshGetMesh(mesh, m);
1798:   PetscObjectGetName((PetscObject) section, &name);
1799:   SectionIntGetSection(section, s);
1800:   m->setIntSection(std::string(name), s);
1801:   return(0);
1802: }

1806: /*@C
1807:   MeshHasSectionPair - Determines whether this mesh has a SectionPair with the given name.

1809:   Not Collective

1811:   Input Parameters:
1812: + mesh - The Mesh object
1813: - name - The section name

1815:   Output Parameter:
1816: . flag - True if the SectionPair is present in the Mesh

1818:   Level: intermediate

1820: .keywords: mesh, elements
1821: .seealso: MeshCreate()
1822: @*/
1823: PetscErrorCode MeshHasSectionPair(Mesh mesh, const char name[], PetscTruth *flag)
1824: {
1825:   ALE::Obj<ALE::Mesh> m;
1826:   PetscErrorCode      ierr;

1829:   MeshGetMesh(mesh, m);
1830:   *flag = (PetscTruth) m->hasPairSection(std::string(name));
1831:   return(0);
1832: }

1836: /*@C
1837:   MeshGetSectionPair - Returns a SectionPair of the given name from the Mesh.

1839:   Collective on Mesh

1841:   Input Parameters:
1842: + mesh - The Mesh object
1843: - name - The section name

1845:   Output Parameter:
1846: . section - The SectionPair

1848:   Note: The section is a new object, and must be destroyed by the user

1850:   Level: intermediate

1852: .keywords: mesh, elements
1853: .seealso: MeshCreate()
1854: @*/
1855: PetscErrorCode MeshGetSectionPair(Mesh mesh, const char name[], SectionPair *section)
1856: {
1857:   ALE::Obj<ALE::Mesh> m;
1858:   PetscErrorCode      ierr;

1861:   MeshGetMesh(mesh, m);
1862:   SectionPairCreate(m->comm(), section);
1863:   PetscObjectSetName((PetscObject) *section, name);
1864:   SectionPairSetSection(*section, m->getPairSection(std::string(name)));
1865:   return(0);
1866: }

1870: /*@C
1871:   MeshSetSectionPair - Puts a SectionPair of the given name into the Mesh.

1873:   Collective on Mesh

1875:   Input Parameters:
1876: + mesh - The Mesh object
1877: - section - The SectionPair

1879:   Note: This takes the section name from the PETSc object

1881:   Level: intermediate

1883: .keywords: mesh, elements
1884: .seealso: MeshCreate()
1885: @*/
1886: PetscErrorCode MeshSetSectionPair(Mesh mesh, SectionPair section)
1887: {
1888:   ALE::Obj<ALE::Mesh> m;
1889:   ALE::Obj<ALE::Mesh::pair_section_type> s;
1890:   const char         *name;
1891:   PetscErrorCode      ierr;

1894:   MeshGetMesh(mesh, m);
1895:   PetscObjectGetName((PetscObject) section, &name);
1896:   SectionPairGetSection(section, s);
1897:   m->setPairSection(std::string(name), s);
1898:   return(0);
1899: }

1903: /*@C
1904:   SectionGetArray - Returns the array underlying the Section.

1906:   Not Collective

1908:   Input Parameters:
1909: + mesh - The Mesh object
1910: - name - The section name

1912:   Output Parameters:
1913: + numElements - The number of mesh element with values
1914: . fiberDim - The number of values per element
1915: - array - The array

1917:   Level: intermediate

1919: .keywords: mesh, elements
1920: .seealso: MeshCreate()
1921: @*/
1922: PetscErrorCode SectionGetArray(Mesh mesh, const char name[], PetscInt *numElements, PetscInt *fiberDim, PetscScalar *array[])
1923: {
1924:   ALE::Obj<ALE::Mesh> m;
1925:   PetscErrorCode      ierr;

1928:   MeshGetMesh(mesh, m);
1929:   const Obj<ALE::Mesh::real_section_type>&        section = m->getRealSection(std::string(name));
1930:   const ALE::Mesh::real_section_type::patch_type  patch   = 0;
1931:   if (!section->hasPatch(patch)) {
1932:     *numElements = 0;
1933:     *fiberDim    = 0;
1934:     *array       = NULL;
1935:     return(0);
1936:   }
1937:   const ALE::Mesh::real_section_type::chart_type& chart   = section->getPatch(patch);
1938: /*   const int                                  depth   = m->getTopology()->depth(patch, *chart.begin()); */
1939: /*   *numElements = m->getTopology()->depthStratum(patch, depth)->size(); */
1940: /*   *fiberDim    = section->getFiberDimension(patch, *chart.begin()); */
1941: /*   *array       = (PetscScalar *) section->restrict(patch); */
1942:   int fiberDimMin = section->getFiberDimension(patch, *chart.begin());
1943:   int numElem     = 0;

1945:   for(ALE::Mesh::real_section_type::chart_type::iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
1946:     const int fiberDim = section->getFiberDimension(patch, *c_iter);

1948:     if (fiberDim < fiberDimMin) fiberDimMin = fiberDim;
1949:   }
1950:   for(ALE::Mesh::real_section_type::chart_type::iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
1951:     const int fiberDim = section->getFiberDimension(patch, *c_iter);

1953:     numElem += fiberDim/fiberDimMin;
1954:   }
1955:   *numElements = numElem;
1956:   *fiberDim    = fiberDimMin;
1957:   *array       = (PetscScalar *) section->restrict(patch);
1958:   return(0);
1959: }

1963: /*@C
1964:   BCSectionGetArray - Returns the array underlying the BCSection.

1966:   Not Collective

1968:   Input Parameters:
1969: + mesh - The Mesh object
1970: - name - The section name

1972:   Output Parameters:
1973: + numElements - The number of mesh element with values
1974: . fiberDim - The number of values per element
1975: - array - The array

1977:   Level: intermediate

1979: .keywords: mesh, elements
1980: .seealso: MeshCreate()
1981: @*/
1982: PetscErrorCode BCSectionGetArray(Mesh mesh, const char name[], PetscInt *numElements, PetscInt *fiberDim, PetscInt *array[])
1983: {
1984:   ALE::Obj<ALE::Mesh> m;
1985:   PetscErrorCode      ierr;

1988:   MeshGetMesh(mesh, m);
1989:   const Obj<ALE::Mesh::int_section_type>&       section = m->getIntSection(std::string(name));
1990:   const ALE::Mesh::int_section_type::patch_type patch   = 0;
1991:   if (!section->hasPatch(patch)) {
1992:     *numElements = 0;
1993:     *fiberDim    = 0;
1994:     *array       = NULL;
1995:     return(0);
1996:   }
1997:   const ALE::Mesh::int_section_type::chart_type& chart = section->getPatch(patch);
1998:   int fiberDimMin = section->getFiberDimension(patch, *chart.begin());
1999:   int numElem     = 0;

2001:   for(ALE::Mesh::int_section_type::chart_type::iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
2002:     const int fiberDim = section->getFiberDimension(patch, *c_iter);

2004:     if (fiberDim < fiberDimMin) fiberDimMin = fiberDim;
2005:   }
2006:   for(ALE::Mesh::int_section_type::chart_type::iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
2007:     const int fiberDim = section->getFiberDimension(patch, *c_iter);

2009:     numElem += fiberDim/fiberDimMin;
2010:   }
2011:   *numElements = numElem;
2012:   *fiberDim    = fiberDimMin;
2013:   *array       = (PetscInt *) section->restrict(patch);
2014:   return(0);
2015: }

2019: PetscErrorCode BCFUNCGetArray(Mesh mesh, PetscInt *numElements, PetscInt *fiberDim, PetscScalar *array[])
2020: {
2021:   ALE::Obj<ALE::Mesh> m;
2022:   PetscErrorCode      ierr;

2025:   MeshGetMesh(mesh, m);
2026:   ALE::Mesh::bc_values_type& bcValues = m->getBCValues();
2027:   *numElements = bcValues.size();
2028:   *fiberDim    = 4;
2029:   *array       = new PetscScalar[(*numElements)*(*fiberDim)];
2030:   for(int bcf = 1; bcf <= (int) bcValues.size(); ++bcf) {
2031:     (*array)[(bcf-1)*4+0] = bcValues[bcf].rho;
2032:     (*array)[(bcf-1)*4+1] = bcValues[bcf].u;
2033:     (*array)[(bcf-1)*4+2] = bcValues[bcf].v;
2034:     (*array)[(bcf-1)*4+3] = bcValues[bcf].p;
2035:   }
2036:   return(0);
2037: }

2041: PetscErrorCode WritePyLithVertices(Mesh mesh, PetscViewer viewer)
2042: {
2043:   ALE::Obj<ALE::Mesh> m;

2046:   MeshGetMesh(mesh, m);
2047:   return ALE::PyLith::Viewer::writeVertices(m, viewer);
2048: }

2052: PetscErrorCode WritePyLithElements(Mesh mesh, SectionInt material, PetscViewer viewer)
2053: {
2054:   ALE::Obj<ALE::Mesh> m;
2055:   ALE::Obj<ALE::Mesh::int_section_type> s;

2058:   MeshGetMesh(mesh, m);
2059:   SectionIntGetSection(material, s);
2060:   return ALE::PyLith::Viewer::writeElements(m, s, viewer);
2061: }

2065: PetscErrorCode WritePyLithVerticesLocal(Mesh mesh, PetscViewer viewer)
2066: {
2067:   ALE::Obj<ALE::Mesh> m;

2070:   MeshGetMesh(mesh, m);
2071:   return ALE::PyLith::Viewer::writeVerticesLocal(m, viewer);
2072: }

2076: PetscErrorCode WritePyLithElementsLocal(Mesh mesh, SectionInt material, PetscViewer viewer)
2077: {
2078:   ALE::Obj<ALE::Mesh> m;
2079:   ALE::Obj<ALE::Mesh::int_section_type> s;

2082:   MeshGetMesh(mesh, m);
2083:   SectionIntGetSection(material, s);
2084:   return ALE::PyLith::Viewer::writeElementsLocal(m, s, viewer);
2085: }

2089: PetscErrorCode WritePyLithTractionsLocal(Mesh mesh, PetscViewer viewer)
2090: {
2091:   ALE::Obj<ALE::Mesh> m;

2094:   MeshGetMesh(mesh, m);
2095:   return ALE::PyLith::Viewer::writeTractionsLocal(m, m->getRealSection("tractions"), viewer);
2096: }