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", §ion);
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: }