Actual source code: mesh.c
1:
2: #include private/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: PetscTruth MeshRegisterAllCalled = PETSC_FALSE;
15: PetscFList MeshList;
17: EXTERN PetscErrorCode MeshView_Mesh(Mesh, PetscViewer);
18: EXTERN PetscErrorCode MeshRefine_Mesh(Mesh, MPI_Comm, Mesh *);
19: EXTERN PetscErrorCode MeshCoarsenHierarchy_Mesh(Mesh, int, Mesh **);
20: EXTERN PetscErrorCode MeshGetInterpolation_Mesh(Mesh, Mesh, Mat *, Vec *);
22: EXTERN PetscErrorCode updateOperatorCompat(Mat, const ALE::Obj<ALECompat::Mesh::real_section_type>&, const ALE::Obj<ALECompat::Mesh::order_type>&, const ALECompat::Mesh::point_type&, PetscScalar[], InsertMode);
27: /*
28: Private routine to delete internal tag storage when a communicator is freed.
30: This is called by MPI, not by users.
34: we do not use PetscFree() since it is unsafe after PetscFinalize()
35: */
36: PetscMPIInt Mesh_DelTag(MPI_Comm comm,PetscMPIInt keyval,void* attr_val,void* extra_state)
37: {
38: free(attr_val);
39: return(MPI_SUCCESS);
40: }
45: PetscErrorCode MeshFinalize()
46: {
48: ALE::Mesh::NumberingFactory::singleton(0, 0, true);
49: return(0);
50: }
54: PetscErrorCode MeshView_Sieve_Ascii(const ALE::Obj<ALE::Mesh>& mesh, PetscViewer viewer)
55: {
56: PetscViewerFormat format;
57: PetscErrorCode ierr;
60: PetscViewerGetFormat(viewer, &format);
61: if (format == PETSC_VIEWER_ASCII_VTK) {
62: VTKViewer::writeHeader(viewer);
63: VTKViewer::writeVertices(mesh, viewer);
64: VTKViewer::writeElements(mesh, viewer);
65: } else if (format == PETSC_VIEWER_ASCII_PYLITH) {
66: char *filename;
67: char connectFilename[2048];
68: char coordFilename[2048];
70: PetscViewerFileGetName(viewer, &filename);
71: PetscViewerFileSetMode(viewer, FILE_MODE_WRITE);
72: PetscStrcpy(connectFilename, filename);
73: PetscStrcat(connectFilename, ".connect");
74: PetscViewerFileSetName(viewer, connectFilename);
75: ALE::PyLith::Viewer::writeElements(mesh, mesh->getIntSection("material"), viewer);
76: PetscStrcpy(coordFilename, filename);
77: PetscStrcat(coordFilename, ".coord");
78: PetscViewerFileSetName(viewer, coordFilename);
79: ALE::PyLith::Viewer::writeVertices(mesh, viewer);
80: PetscViewerFileSetMode(viewer, FILE_MODE_READ);
81: PetscExceptionTry1(PetscViewerFileSetName(viewer, filename), PETSC_ERR_FILE_OPEN);
82: if (PetscExceptionValue(ierr)) {
83: /* this means that a caller above me has also tryed this exception so I don't handle it here, pass it up */
84: } else if (PetscExceptionCaught(ierr, PETSC_ERR_FILE_OPEN)) {
85: 0;
86: }
87:
88: } else if (format == PETSC_VIEWER_ASCII_PYLITH_LOCAL) {
89: PetscViewer connectViewer, coordViewer;
90: char *filename;
91: char localFilename[2048];
92: int rank = mesh->commRank();
94: PetscViewerFileGetName(viewer, &filename);
96: sprintf(localFilename, "%s.%d.connect", filename, rank);
97: PetscViewerCreate(PETSC_COMM_SELF, &connectViewer);
98: PetscViewerSetType(connectViewer, PETSC_VIEWER_ASCII);
99: PetscViewerSetFormat(connectViewer, PETSC_VIEWER_ASCII_PYLITH);
100: PetscViewerFileSetName(connectViewer, localFilename);
101: ALE::PyLith::Viewer::writeElementsLocal(mesh, mesh->getIntSection("material"), connectViewer);
102: PetscViewerDestroy(connectViewer);
104: sprintf(localFilename, "%s.%d.coord", filename, rank);
105: PetscViewerCreate(PETSC_COMM_SELF, &coordViewer);
106: PetscViewerSetType(coordViewer, PETSC_VIEWER_ASCII);
107: PetscViewerSetFormat(coordViewer, PETSC_VIEWER_ASCII_PYLITH);
108: PetscViewerFileSetName(coordViewer, localFilename);
109: ALE::PyLith::Viewer::writeVerticesLocal(mesh, coordViewer);
110: PetscViewerDestroy(coordViewer);
111: } else if (format == PETSC_VIEWER_ASCII_PCICE) {
112: char *filename;
113: char coordFilename[2048];
114: PetscTruth isConnect;
115: size_t len;
117: PetscViewerFileGetName(viewer, &filename);
118: PetscStrlen(filename, &len);
119: PetscStrcmp(&(filename[len-5]), ".lcon", &isConnect);
120: if (!isConnect) {
121: SETERRQ1(PETSC_ERR_ARG_WRONG, "Invalid element connectivity filename: %s", filename);
122: }
123: ALE::PCICE::Viewer::writeElements(mesh, viewer);
124: PetscStrncpy(coordFilename, filename, len-5);
125: coordFilename[len-5] = '\0';
126: PetscStrcat(coordFilename, ".nodes");
127: PetscViewerFileSetName(viewer, coordFilename);
128: ALE::PCICE::Viewer::writeVertices(mesh, viewer);
129: } else {
130: int dim = mesh->getDimension();
132: PetscViewerASCIIPrintf(viewer, "Mesh in %d dimensions:\n", dim);
133: for(int d = 0; d <= dim; d++) {
134: // FIX: Need to globalize
135: PetscViewerASCIIPrintf(viewer, " %d %d-cells\n", mesh->depthStratum(d)->size(), d);
136: }
137: }
138: PetscViewerFlush(viewer);
139: return(0);
140: }
144: PetscErrorCode MeshCompatView_Sieve_Ascii(const ALE::Obj<ALECompat::Mesh>& mesh, PetscViewer viewer)
145: {
146: PetscViewerFormat format;
147: PetscErrorCode ierr;
150: PetscViewerGetFormat(viewer, &format);
151: if (format == PETSC_VIEWER_ASCII_PYLITH) {
152: char *filename;
153: char connectFilename[2048];
154: char coordFilename[2048];
156: PetscViewerFileGetName(viewer, &filename);
157: PetscViewerFileSetMode(viewer, FILE_MODE_WRITE);
158: PetscStrcpy(connectFilename, filename);
159: PetscStrcat(connectFilename, ".connect");
160: PetscViewerFileSetName(viewer, connectFilename);
161: ALECompat::PyLith::Viewer::writeElements(mesh, mesh->getIntSection("material"), viewer);
162: PetscStrcpy(coordFilename, filename);
163: PetscStrcat(coordFilename, ".coord");
164: PetscViewerFileSetName(viewer, coordFilename);
165: ALECompat::PyLith::Viewer::writeVertices(mesh, viewer);
166: PetscViewerFileSetMode(viewer, FILE_MODE_READ);
167: PetscExceptionTry1(PetscViewerFileSetName(viewer, filename), PETSC_ERR_FILE_OPEN);
168: if (PetscExceptionValue(ierr)) {
169: /* this means that a caller above me has also tryed this exception so I don't handle it here, pass it up */
170: } else if (PetscExceptionCaught(ierr, PETSC_ERR_FILE_OPEN)) {
171: 0;
172: }
173:
174: } else if (format == PETSC_VIEWER_ASCII_PYLITH_LOCAL) {
175: PetscViewer connectViewer, coordViewer;
176: char *filename;
177: char localFilename[2048];
178: int rank = mesh->commRank();
180: PetscViewerFileGetName(viewer, &filename);
182: sprintf(localFilename, "%s.%d.connect", filename, rank);
183: PetscViewerCreate(PETSC_COMM_SELF, &connectViewer);
184: PetscViewerSetType(connectViewer, PETSC_VIEWER_ASCII);
185: PetscViewerSetFormat(connectViewer, PETSC_VIEWER_ASCII_PYLITH);
186: PetscViewerFileSetName(connectViewer, localFilename);
187: ALECompat::PyLith::Viewer::writeElementsLocal(mesh, mesh->getIntSection("material"), connectViewer);
188: PetscViewerDestroy(connectViewer);
190: sprintf(localFilename, "%s.%d.coord", filename, rank);
191: PetscViewerCreate(PETSC_COMM_SELF, &coordViewer);
192: PetscViewerSetType(coordViewer, PETSC_VIEWER_ASCII);
193: PetscViewerSetFormat(coordViewer, PETSC_VIEWER_ASCII_PYLITH);
194: PetscViewerFileSetName(coordViewer, localFilename);
195: ALECompat::PyLith::Viewer::writeVerticesLocal(mesh, coordViewer);
196: PetscViewerDestroy(coordViewer);
198: if (mesh->hasPairSection("split")) {
199: PetscViewer splitViewer;
201: sprintf(localFilename, "%s.%d.split", filename, rank);
202: PetscViewerCreate(PETSC_COMM_SELF, &splitViewer);
203: PetscViewerSetType(splitViewer, PETSC_VIEWER_ASCII);
204: PetscViewerSetFormat(splitViewer, PETSC_VIEWER_ASCII_PYLITH);
205: PetscViewerFileSetName(splitViewer, localFilename);
206: ALECompat::PyLith::Viewer::writeSplitLocal(mesh, mesh->getPairSection("split"), splitViewer);
207: PetscViewerDestroy(splitViewer);
208: }
210: if (mesh->hasRealSection("traction")) {
211: PetscViewer tractionViewer;
213: sprintf(localFilename, "%s.%d.traction", filename, rank);
214: PetscViewerCreate(PETSC_COMM_SELF, &tractionViewer);
215: PetscViewerSetType(tractionViewer, PETSC_VIEWER_ASCII);
216: PetscViewerSetFormat(tractionViewer, PETSC_VIEWER_ASCII_PYLITH);
217: PetscViewerFileSetName(tractionViewer, localFilename);
218: ALECompat::PyLith::Viewer::writeTractionsLocal(mesh, mesh->getRealSection("traction"), tractionViewer);
219: PetscViewerDestroy(tractionViewer);
220: }
221: } else {
222: int dim = mesh->getDimension();
224: PetscViewerASCIIPrintf(viewer, "Mesh in %d dimensions:\n", dim);
225: for(int d = 0; d <= dim; d++) {
226: // FIX: Need to globalize
227: PetscViewerASCIIPrintf(viewer, " %d %d-cells\n", mesh->getTopology()->depthStratum(0, d)->size(), d);
228: }
229: }
230: PetscViewerFlush(viewer);
231: return(0);
232: }
236: PetscErrorCode MeshView_Sieve(const ALE::Obj<ALE::Mesh>& mesh, PetscViewer viewer)
237: {
238: PetscTruth iascii, isbinary, isdraw;
242: PetscTypeCompare((PetscObject) viewer, PETSC_VIEWER_ASCII, &iascii);
243: PetscTypeCompare((PetscObject) viewer, PETSC_VIEWER_BINARY, &isbinary);
244: PetscTypeCompare((PetscObject) viewer, PETSC_VIEWER_DRAW, &isdraw);
246: if (iascii){
247: MeshView_Sieve_Ascii(mesh, viewer);
248: } else if (isbinary) {
249: SETERRQ(PETSC_ERR_SUP, "Binary viewer not implemented for Mesh");
250: } else if (isdraw){
251: SETERRQ(PETSC_ERR_SUP, "Draw viewer not implemented for Mesh");
252: } else {
253: SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported by this mesh object", ((PetscObject)viewer)->type_name);
254: }
255: return(0);
256: }
260: PetscErrorCode MeshView_Mesh(Mesh mesh, PetscViewer viewer)
261: {
265: if (!mesh->mcompat.isNull()) {
266: MeshCompatView_Sieve_Ascii(mesh->mcompat, viewer);
267: } else {
268: MeshView_Sieve(mesh->m, viewer);
269: }
270: return(0);
271: }
275: /*@C
276: MeshView - Views a Mesh object.
278: Collective on Mesh
280: Input Parameters:
281: + mesh - the mesh
282: - viewer - an optional visualization context
284: Notes:
285: The available visualization contexts include
286: + PETSC_VIEWER_STDOUT_SELF - standard output (default)
287: - PETSC_VIEWER_STDOUT_WORLD - synchronized standard
288: output where only the first processor opens
289: the file. All other processors send their
290: data to the first processor to print.
292: You can change the format the mesh is printed using the
293: option PetscViewerSetFormat().
295: The user can open alternative visualization contexts with
296: + PetscViewerASCIIOpen() - Outputs mesh to a specified file
297: . PetscViewerBinaryOpen() - Outputs mesh in binary to a
298: specified file; corresponding input uses MeshLoad()
299: . PetscViewerDrawOpen() - Outputs mesh to an X window display
301: The user can call PetscViewerSetFormat() to specify the output
302: format of ASCII printed objects (when using PETSC_VIEWER_STDOUT_SELF,
303: PETSC_VIEWER_STDOUT_WORLD and PetscViewerASCIIOpen). Available formats include
304: + PETSC_VIEWER_ASCII_DEFAULT - default, prints mesh information
305: - PETSC_VIEWER_ASCII_VTK - outputs a VTK file describing the mesh
307: Level: beginner
309: Concepts: mesh^printing
310: Concepts: mesh^saving to disk
312: .seealso: PetscViewerASCIIOpen(), PetscViewerDrawOpen(), PetscViewerBinaryOpen(),
313: MeshLoad(), PetscViewerCreate()
314: @*/
315: PetscErrorCode MeshView(Mesh mesh, PetscViewer viewer)
316: {
322: if (!viewer) {
323: PetscViewerASCIIGetStdout(mesh->comm,&viewer);
324: }
329: (*mesh->ops->view)(mesh, viewer);
331: return(0);
332: }
336: /*@C
337: MeshLoad - Create a mesh topology from the saved data in a viewer.
339: Collective on Viewer
341: Input Parameter:
342: . viewer - The viewer containing the data
344: Output Parameters:
345: . mesh - the mesh object
347: Level: advanced
349: .seealso MeshView()
351: @*/
352: PetscErrorCode MeshLoad(PetscViewer viewer, Mesh *mesh)
353: {
354: SETERRQ(PETSC_ERR_SUP, "");
355: }
359: /*@C
360: MeshGetMesh - Gets the internal mesh object
362: Not collective
364: Input Parameter:
365: . mesh - the mesh object
367: Output Parameter:
368: . m - the internal mesh object
369:
370: Level: advanced
372: .seealso MeshCreate(), MeshSetMesh()
374: @*/
375: PetscErrorCode MeshGetMesh(Mesh mesh, ALE::Obj<ALE::Mesh>& m)
376: {
379: m = mesh->m;
380: return(0);
381: }
385: /*@C
386: MeshSetMesh - Sets the internal mesh object
388: Not collective
390: Input Parameters:
391: + mesh - the mesh object
392: - m - the internal mesh object
393:
394: Level: advanced
396: .seealso MeshCreate(), MeshGetMesh()
398: @*/
399: PetscErrorCode MeshSetMesh(Mesh mesh, const ALE::Obj<ALE::Mesh>& m)
400: {
403: mesh->m = m;
404: return(0);
405: }
409: template<typename Section>
410: PetscErrorCode MeshCreateMatrix(const Obj<ALE::Mesh>& mesh, const Obj<Section>& section, MatType mtype, Mat *J)
411: {
412: const ALE::Obj<typename ALE::Mesh::order_type>& order = mesh->getFactory()->getGlobalOrder(mesh, "default", section);
413: int localSize = order->getLocalSize();
414: int globalSize = order->getGlobalSize();
415: PetscTruth isShell, isBlock, isSeqBlock, isMPIBlock;
419: MatCreate(mesh->comm(), J);
420: MatSetSizes(*J, localSize, localSize, globalSize, globalSize);
421: MatSetType(*J, mtype);
422: MatSetFromOptions(*J);
423: PetscStrcmp(mtype, MATSHELL, &isShell);
424: PetscStrcmp(mtype, MATBAIJ, &isBlock);
425: PetscStrcmp(mtype, MATSEQBAIJ, &isSeqBlock);
426: PetscStrcmp(mtype, MATMPIBAIJ, &isMPIBlock);
427: if (!isShell) {
428: int bs = 1;
430: if (isBlock || isSeqBlock || isMPIBlock) {
431: bs = section->getFiberDimension(*section->getChart().begin());
432: }
433: preallocateOperator(mesh, bs, section->getAtlas(), order, *J);
434: }
435: return(0);
436: }
440: PetscErrorCode MeshCreateMatrix(Mesh mesh, SectionReal section, MatType mtype, Mat *J)
441: {
442: ALE::Obj<ALE::Mesh> m;
443: ALE::Obj<ALE::Mesh::real_section_type> s;
447: MeshGetMesh(mesh, m);
448: SectionRealGetSection(section, s);
449: MeshCreateMatrix(m, s, mtype, J);
450: PetscObjectCompose((PetscObject) *J, "mesh", (PetscObject) mesh);
451: return(0);
452: }
456: PetscErrorCode MeshGetVertexMatrix(Mesh mesh, MatType mtype, Mat *J)
457: {
458: SectionReal section;
462: MeshGetVertexSectionReal(mesh, 1, §ion);
463: MeshCreateMatrix(mesh, section, mtype, J);
464: SectionRealDestroy(section);
465: return(0);
466: }
470: /*@C
471: MeshGetMatrix - Creates a matrix with the correct parallel layout required for
472: computing the Jacobian on a function defined using the informatin in Mesh.
474: Collective on Mesh
476: Input Parameter:
477: + mesh - the mesh object
478: - mtype - Supported types are MATSEQAIJ, MATMPIAIJ, MATSEQBAIJ, MATMPIBAIJ, MATSEQSBAIJ, MATMPISBAIJ,
479: or any type which inherits from one of these (such as MATAIJ, MATLUSOL, etc.).
481: Output Parameters:
482: . J - matrix with the correct nonzero preallocation
483: (obviously without the correct Jacobian values)
485: Level: advanced
487: Notes: This properly preallocates the number of nonzeros in the sparse matrix so you
488: do not need to do it yourself.
490: .seealso ISColoringView(), ISColoringGetIS(), MatFDColoringCreate(), DASetBlockFills()
492: @*/
493: PetscErrorCode MeshGetMatrix(Mesh mesh, MatType mtype, Mat *J)
494: {
495: ALE::Obj<ALE::Mesh> m;
496: PetscTruth flag;
497: PetscErrorCode ierr;
500: MeshHasSectionReal(mesh, "default", &flag);
501: if (!flag) SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "Must set default section");
502: MeshGetMesh(mesh, m);
503: MeshCreateMatrix(m, m->getRealSection("default"), mtype, J);
504: PetscObjectCompose((PetscObject) *J, "mesh", (PetscObject) mesh);
505: return(0);
506: }
510: /*@C
511: MeshCreate - Creates a DM object, used to manage data for an unstructured problem
512: described by a Sieve.
514: Collective on MPI_Comm
516: Input Parameter:
517: . comm - the processors that will share the global vector
519: Output Parameters:
520: . mesh - the mesh object
522: Level: advanced
524: .seealso MeshDestroy(), MeshCreateGlobalVector(), MeshGetGlobalIndices()
526: @*/
527: PetscErrorCode MeshCreate(MPI_Comm comm,Mesh *mesh)
528: {
530: Mesh p;
534: *mesh = PETSC_NULL;
535: #ifndef PETSC_USE_DYNAMIC_LIBRARIES
536: DMInitializePackage(PETSC_NULL);
537: #endif
539: PetscHeaderCreate(p,_p_Mesh,struct _MeshOps,MESH_COOKIE,0,"Mesh",comm,MeshDestroy,MeshView);
540: p->ops->view = MeshView_Mesh;
541: p->ops->destroy = PETSC_NULL;
542: p->ops->createglobalvector = MeshCreateGlobalVector;
543: p->ops->getcoloring = PETSC_NULL;
544: p->ops->getmatrix = MeshGetMatrix;
545: p->ops->getinterpolation = MeshGetInterpolation_Mesh;
546: p->ops->getinjection = PETSC_NULL;
547: p->ops->refine = MeshRefine_Mesh;
548: p->ops->coarsen = PETSC_NULL;
549: p->ops->refinehierarchy = PETSC_NULL;
550: p->ops->coarsenhierarchy = MeshCoarsenHierarchy_Mesh;
552: PetscObjectChangeTypeName((PetscObject) p, "sieve");
554: p->m = PETSC_NULL;
555: p->globalScatter = PETSC_NULL;
556: p->lf = PETSC_NULL;
557: p->lj = PETSC_NULL;
558: p->mcompat = PETSC_NULL;
559: p->data = PETSC_NULL;
560: *mesh = p;
561: return(0);
562: }
566: /*@
567: MeshDestroy - Destroys a mesh.
569: Collective on Mesh
571: Input Parameter:
572: . mesh - the mesh object
574: Level: advanced
576: .seealso MeshCreate(), MeshCreateGlobalVector(), MeshGetGlobalIndices()
577: @*/
578: PetscErrorCode MeshDestroy(Mesh mesh)
579: {
580: PetscErrorCode ierr;
583: if (--mesh->refct > 0) return(0);
584: if (mesh->globalScatter) {VecScatterDestroy(mesh->globalScatter);}
585: mesh->m = PETSC_NULL;
586: PetscHeaderDestroy(mesh);
587: return(0);
588: }
592: /*@C
593: MeshSetType - Sets the Mesh type
595: Collective on Mesh
597: Input Parameters:
598: + mesh - the Mesh context
599: - type - the type
601: Options Database Key:
602: . -mesh_type <method> - Sets the type; use -help for a list
603: of available types (for instance, cartesian or sieve)
605: Notes:
606: See "petsc/include/petscmesh.h" for available types (for instance,
607: MESHCARTESIAN or MESHSIEVE).
609: Level: intermediate
611: .keywords: Mesh, set, typr
612: .seealso: MeshGetType(), MeshType
613: @*/
614: PetscErrorCode MeshSetType(Mesh mesh, MeshType type)
615: {
616: PetscErrorCode ierr,(*r)(Mesh);
617: PetscTruth match;
623: PetscTypeCompare((PetscObject)mesh,type,&match);
624: if (match) return(0);
626: PetscFListFind(MeshList,mesh->comm,type,(void (**)(void)) &r);
627: if (!r) SETERRQ1(PETSC_ERR_ARG_UNKNOWN_TYPE,"Unable to find requested Mesh type %s",type);
628: /* Destroy the previous private Mesh context */
629: if (mesh->ops->destroy) { (*mesh->ops->destroy)(mesh); }
630: /* Reinitialize function pointers in MeshOps structure */
631: PetscMemzero(mesh->ops, sizeof(struct _MeshOps));
632: /* Call the MeshCreate_XXX routine for this particular mesh */
633: (*r)(mesh);
634: PetscObjectChangeTypeName((PetscObject) mesh, type);
635: return(0);
636: }
640: /*@C
641: MeshGetType - Gets the Mesh type as a string from the Mesh object.
643: Not Collective
645: Input Parameter:
646: . mesh - Mesh context
648: Output Parameter:
649: . name - name of Mesh type
651: Level: intermediate
653: .keywords: Mesh, get, type
654: .seealso: MeshSetType()
655: @*/
656: PetscErrorCode MeshGetType(Mesh mesh,MeshType *type)
657: {
661: *type = mesh->type_name;
662: return(0);
663: }
667: /*@C
668: MeshRegister - See MeshRegisterDynamic()
670: Level: advanced
671: @*/
672: PetscErrorCode MeshRegister(const char sname[],const char path[],const char name[],PetscErrorCode (*function)(Mesh))
673: {
675: char fullname[PETSC_MAX_PATH_LEN];
678: PetscFListConcat(path,name,fullname);
679: PetscFListAdd(&MeshList,sname,fullname,(void (*)(void))function);
680: return(0);
681: }
684: EXTERN PetscErrorCode MeshCreate_Cartesian(Mesh);
689: /*@C
690: MeshRegisterAll - Registers all of the Mesh types in the Mesh package.
692: Not Collective
694: Level: advanced
696: .keywords: Mesh, register, all
697: .seealso: MeshRegisterDestroy()
698: @*/
699: PetscErrorCode MeshRegisterAll(const char path[])
700: {
704: MeshRegisterAllCalled = PETSC_TRUE;
706: MeshRegisterDynamic(MESHCARTESIAN, path, "MeshCreate_Cartesian", MeshCreate_Cartesian);
707: return(0);
708: }
712: /*@
713: MeshRegisterDestroy - Frees the list of Mesh types that were
714: registered by MeshRegister().
716: Not Collective
718: Level: advanced
720: .keywords: Mesh, register, destroy
721: .seealso: MeshRegister(), MeshRegisterAll()
722: @*/
723: PetscErrorCode MeshRegisterDestroy(void)
724: {
728: PetscFListDestroy(&MeshList);
729: MeshRegisterAllCalled = PETSC_FALSE;
730: return(0);
731: }
735: /*@C
736: MeshCreateGlobalVector - Creates a vector of the correct size to be gathered into
737: by the mesh.
739: Collective on Mesh
741: Input Parameter:
742: . mesh - the mesh object
744: Output Parameters:
745: . gvec - the global vector
747: Level: advanced
749: Notes: Once this has been created you cannot add additional arrays or vectors to be packed.
751: .seealso MeshDestroy(), MeshCreate(), MeshGetGlobalIndices()
753: @*/
754: PetscErrorCode MeshCreateGlobalVector(Mesh mesh, Vec *gvec)
755: {
756: ALE::Obj<ALE::Mesh> m;
757: PetscTruth flag;
761: MeshHasSectionReal(mesh, "default", &flag);
762: if (!flag) SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "Must set default section");
763: MeshGetMesh(mesh, m);
764: const ALE::Obj<ALE::Mesh::order_type>& order = m->getFactory()->getGlobalOrder(m, "default", m->getRealSection("default"));
766: VecCreate(m->comm(), gvec);
767: VecSetSizes(*gvec, order->getLocalSize(), order->getGlobalSize());
768: VecSetFromOptions(*gvec);
769: return(0);
770: }
774: /*@C
775: MeshGetGlobalIndices - Gets the global indices for all the local entries
777: Collective on Mesh
779: Input Parameter:
780: . mesh - the mesh object
782: Output Parameters:
783: . idx - the individual indices for each packed vector/array
784:
785: Level: advanced
787: Notes:
788: The idx parameters should be freed by the calling routine with PetscFree()
790: .seealso MeshDestroy(), MeshCreateGlobalVector(), MeshCreate()
792: @*/
793: PetscErrorCode MeshGetGlobalIndices(Mesh mesh,PetscInt *idx[])
794: {
795: SETERRQ(PETSC_ERR_SUP, "");
796: }
800: PetscErrorCode MeshCreateGlobalScatter(Mesh mesh, SectionReal section, VecScatter *scatter)
801: {
802: ALE::Obj<ALE::Mesh> m;
803: ALE::Obj<ALE::Mesh::real_section_type> s;
807: MeshGetMesh(mesh, m);
808: SectionRealGetSection(section, s);
809: MeshCreateGlobalScatter(m, s, scatter);
810: return(0);
811: }
815: template<typename Section>
816: PetscErrorCode MeshCreateGlobalScatter(const ALE::Obj<ALE::Mesh>& m, const ALE::Obj<Section>& s, VecScatter *scatter)
817: {
818: typedef ALE::Mesh::real_section_type::index_type index_type;
823: const ALE::Mesh::real_section_type::chart_type& chart = s->getChart();
824: const ALE::Obj<ALE::Mesh::order_type>& globalOrder = m->getFactory()->getGlobalOrder(m, s->getName(), s);
825: int *localIndices, *globalIndices;
826: int localSize = s->size();
827: int localIndx = 0, globalIndx = 0;
828: Vec globalVec, localVec;
829: IS localIS, globalIS;
831: VecCreate(m->comm(), &globalVec);
832: VecSetSizes(globalVec, globalOrder->getLocalSize(), PETSC_DETERMINE);
833: VecSetFromOptions(globalVec);
834: // Loop over all local points
835: PetscMalloc(localSize*sizeof(int), &localIndices);
836: PetscMalloc(localSize*sizeof(int), &globalIndices);
837: for(ALE::Mesh::real_section_type::chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
838: // Map local indices to global indices
839: s->getIndices(*p_iter, localIndices, &localIndx, 0, true);
840: s->getIndices(*p_iter, globalOrder, globalIndices, &globalIndx, 0, true);
841: }
842: if (localIndx != localSize) SETERRQ2(PETSC_ERR_ARG_SIZ, "Invalid number of local indices %d, should be %d", localIndx, localSize);
843: if (globalIndx != localSize) SETERRQ2(PETSC_ERR_ARG_SIZ, "Invalid number of global indices %d, should be %d", globalIndx, localSize);
844: if (m->debug()) {
845: globalOrder->view("Global Order");
846: for(int i = 0; i < localSize; ++i) {
847: printf("[%d] localIndex[%d]: %d globalIndex[%d]: %d\n", m->commRank(), i, localIndices[i], i, globalIndices[i]);
848: }
849: }
850: ISCreateGeneral(PETSC_COMM_SELF, localSize, localIndices, &localIS);
851: ISCreateGeneral(PETSC_COMM_SELF, localSize, globalIndices, &globalIS);
852: PetscFree(localIndices);
853: PetscFree(globalIndices);
854: VecCreateSeqWithArray(PETSC_COMM_SELF, s->sizeWithBC(), s->restrict(), &localVec);
855: VecScatterCreate(localVec, localIS, globalVec, globalIS, scatter);
856: ISDestroy(globalIS);
857: ISDestroy(localIS);
858: VecDestroy(localVec);
859: VecDestroy(globalVec);
861: return(0);
862: }
866: PetscErrorCode MeshGetGlobalScatter(Mesh mesh, VecScatter *scatter)
867: {
873: if (!mesh->globalScatter) {
874: SectionReal section;
876: MeshGetSectionReal(mesh, "default", §ion);
877: MeshCreateGlobalScatter(mesh, section, &mesh->globalScatter);
878: SectionRealDestroy(section);
879: }
880: *scatter = mesh->globalScatter;
881: return(0);
882: }
886: PetscErrorCode MeshSetLocalFunction(Mesh mesh, PetscErrorCode (*lf)(Mesh, SectionReal, SectionReal, void *))
887: {
890: mesh->lf = lf;
891: return(0);
892: }
896: PetscErrorCode MeshSetLocalJacobian(Mesh mesh, PetscErrorCode (*lj)(Mesh, SectionReal, Mat, void *))
897: {
900: mesh->lj = lj;
901: return(0);
902: }
906: PetscErrorCode MeshFormFunction(Mesh mesh, SectionReal X, SectionReal F, void *ctx)
907: {
914: if (mesh->lf) {
915: (*mesh->lf)(mesh, X, F, ctx);
916: }
917: return(0);
918: }
922: PetscErrorCode MeshFormJacobian(Mesh mesh, SectionReal X, Mat J, void *ctx)
923: {
930: if (mesh->lj) {
931: (*mesh->lj)(mesh, X, J, ctx);
932: }
933: return(0);
934: }
938: // Here we assume:
939: // - Assumes 3D and tetrahedron
940: // - The section takes values on vertices and is P1
941: // - Points have the same dimension as the mesh
942: // - All values have the same dimension
943: PetscErrorCode MeshInterpolatePoints(Mesh mesh, SectionReal section, int numPoints, double *points, double **values)
944: {
945: Obj<ALE::Mesh> m;
946: Obj<ALE::Mesh::real_section_type> s;
947: double *v0, *J, *invJ, detJ;
951: MeshGetMesh(mesh, m);
952: SectionRealGetSection(section, s);
953: const Obj<ALE::Mesh::real_section_type>& coordinates = m->getRealSection("coordinates");
954: int embedDim = coordinates->getFiberDimension(*m->depthStratum(0)->begin());
955: int dim = s->getFiberDimension(*m->depthStratum(0)->begin());
957: PetscMalloc3(embedDim,double,&v0,embedDim*embedDim,double,&J,embedDim*embedDim,double,&invJ);
958: PetscMalloc(numPoints*dim * sizeof(double), &values);
959: for(int p = 0; p < numPoints; p++) {
960: double *point = &points[p*embedDim];
961: ALE::Mesh::point_type e = m->locatePoint(point);
962: const ALE::Mesh::real_section_type::value_type *coeff = s->restrictPoint(e);
964: m->computeElementGeometry(coordinates, e, v0, J, invJ, detJ);
965: 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;
966: 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;
967: 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;
969: for(int d = 0; d < dim; d++) {
970: (*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;
971: }
972: }
973: PetscFree3(v0, J, invJ);
974: return(0);
975: }
979: /*@
980: MeshGetMaximumDegree - Return the maximum degree of any mesh vertex
982: Collective on mesh
984: Input Parameter:
985: . mesh - The Mesh
987: Output Parameter:
988: . maxDegree - The maximum number of edges at any vertex
990: Level: beginner
992: .seealso: MeshCreate()
993: @*/
994: PetscErrorCode MeshGetMaximumDegree(Mesh mesh, PetscInt *maxDegree)
995: {
996: Obj<ALE::Mesh> m;
1000: MeshGetMesh(mesh, m);
1001: const ALE::Obj<ALE::Mesh::label_sequence>& vertices = m->depthStratum(0);
1002: const ALE::Obj<ALE::Mesh::sieve_type>& sieve = m->getSieve();
1003: PetscInt maxDeg = -1;
1005: for(ALE::Mesh::label_sequence::iterator v_iter = vertices->begin(); v_iter != vertices->end(); ++v_iter) {
1006: maxDeg = PetscMax(maxDeg, (PetscInt) sieve->support(*v_iter)->size());
1007: }
1008: *maxDegree = maxDeg;
1009: return(0);
1010: }
1012: EXTERN PetscErrorCode assembleFullField(VecScatter, Vec, Vec, InsertMode);
1016: /*@C
1017: restrictVector - Insert values from a global vector into a local ghosted vector
1019: Collective on g
1021: Input Parameters:
1022: + g - The global vector
1023: . l - The local vector
1024: - mode - either ADD_VALUES or INSERT_VALUES, where
1025: ADD_VALUES adds values to any existing entries, and
1026: INSERT_VALUES replaces existing entries with new values
1028: Level: beginner
1030: .seealso: MatSetOption()
1031: @*/
1032: PetscErrorCode restrictVector(Vec g, Vec l, InsertMode mode)
1033: {
1034: VecScatter injection;
1039: PetscObjectQuery((PetscObject) g, "injection", (PetscObject *) &injection);
1040: if (injection) {
1041: VecScatterBegin(injection, g, l, mode, SCATTER_REVERSE);
1042: VecScatterEnd(injection, g, l, mode, SCATTER_REVERSE);
1043: } else {
1044: if (mode == INSERT_VALUES) {
1045: VecCopy(g, l);
1046: } else {
1047: VecAXPY(l, 1.0, g);
1048: }
1049: }
1051: return(0);
1052: }
1056: /*@C
1057: assembleVectorComplete - Insert values from a local ghosted vector into a global vector
1059: Collective on g
1061: Input Parameters:
1062: + g - The global vector
1063: . l - The local vector
1064: - mode - either ADD_VALUES or INSERT_VALUES, where
1065: ADD_VALUES adds values to any existing entries, and
1066: INSERT_VALUES replaces existing entries with new values
1068: Level: beginner
1070: .seealso: MatSetOption()
1071: @*/
1072: PetscErrorCode assembleVectorComplete(Vec g, Vec l, InsertMode mode)
1073: {
1074: VecScatter injection;
1079: PetscObjectQuery((PetscObject) g, "injection", (PetscObject *) &injection);
1080: if (injection) {
1081: VecScatterBegin(injection, l, g, mode, SCATTER_FORWARD);
1082: VecScatterEnd(injection, l, g, mode, SCATTER_FORWARD);
1083: } else {
1084: if (mode == INSERT_VALUES) {
1085: VecCopy(l, g);
1086: } else {
1087: VecAXPY(g, 1.0, l);
1088: }
1089: }
1091: return(0);
1092: }
1096: /*@C
1097: assembleVector - Insert values into a vector
1099: Collective on A
1101: Input Parameters:
1102: + b - the vector
1103: . e - The element number
1104: . v - The values
1105: - mode - either ADD_VALUES or INSERT_VALUES, where
1106: ADD_VALUES adds values to any existing entries, and
1107: INSERT_VALUES replaces existing entries with new values
1109: Level: beginner
1111: .seealso: VecSetOption()
1112: @*/
1113: PetscErrorCode assembleVector(Vec b, PetscInt e, PetscScalar v[], InsertMode mode)
1114: {
1115: Mesh mesh;
1116: ALE::Obj<ALE::Mesh> m;
1117: PetscInt firstElement;
1118: PetscErrorCode ierr;
1122: PetscObjectQuery((PetscObject) b, "mesh", (PetscObject *) &mesh);
1123: MeshGetMesh(mesh, m);
1124: //firstElement = elementBundle->getLocalSizes()[bundle->getCommRank()];
1125: firstElement = 0;
1126: // Must relate b to field
1127: if (mode == INSERT_VALUES) {
1128: m->update(m->getRealSection(std::string("x")), ALE::Mesh::point_type(e + firstElement), v);
1129: } else {
1130: m->updateAdd(m->getRealSection(std::string("x")), ALE::Mesh::point_type(e + firstElement), v);
1131: }
1133: return(0);
1134: }
1138: PetscErrorCode updateOperator(Mat A, const ALE::Obj<ALE::Mesh>& m, 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)
1139: {
1143: const ALE::Mesh::indices_type indicesBlock = m->getIndices(section, e, globalOrder);
1144: const PetscInt *indices = indicesBlock.first;
1145: const int& numIndices = indicesBlock.second;
1148: if (section->debug()) {
1149: printf("[%d]mat for element %d\n", section->commRank(), e);
1150: for(int i = 0; i < numIndices; i++) {
1151: printf("[%d]mat indices[%d] = %d\n", section->commRank(), i, indices[i]);
1152: }
1153: for(int i = 0; i < numIndices; i++) {
1154: printf("[%d]", section->commRank());
1155: for(int j = 0; j < numIndices; j++) {
1156: printf(" %g", array[i*numIndices+j]);
1157: }
1158: printf("\n");
1159: }
1160: }
1161: MatSetValues(A, numIndices, indices, numIndices, indices, array, mode);
1162: if (ierr) {
1163: printf("[%d]ERROR in updateOperator: point %d\n", section->commRank(), e);
1164: for(int i = 0; i < numIndices; i++) {
1165: printf("[%d]mat indices[%d] = %d\n", section->commRank(), i, indices[i]);
1166: }
1167:
1168: }
1170: return(0);
1171: }
1175: PetscErrorCode updateOperatorGeneral(Mat A, const ALE::Obj<ALE::Mesh>& rowM, 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>& colM, 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)
1176: {
1180: const ALE::Mesh::indices_type rowIndicesBlock = rowM->getIndices(rowSection, rowE, rowGlobalOrder);
1181: const PetscInt *rowIndices = rowIndicesBlock.first;
1182: const int& numRowIndices = rowIndicesBlock.second;
1183: const ALE::Mesh::indices_type colIndicesBlock = colM->getIndices(colSection, colE, colGlobalOrder);
1184: const PetscInt *colIndices = colIndicesBlock.first;
1185: const int& numColIndices = colIndicesBlock.second;
1188: if (rowSection->debug()) {
1189: printf("[%d]mat for elements %d %d\n", rowSection->commRank(), rowE, colE);
1190: for(int i = 0; i < numRowIndices; i++) {
1191: printf("[%d]mat row indices[%d] = %d\n", rowSection->commRank(), i, rowIndices[i]);
1192: }
1193: }
1194: if (colSection->debug()) {
1195: for(int i = 0; i < numColIndices; i++) {
1196: printf("[%d]mat col indices[%d] = %d\n", colSection->commRank(), i, colIndices[i]);
1197: }
1198: for(int i = 0; i < numRowIndices; i++) {
1199: printf("[%d]", rowSection->commRank());
1200: for(int j = 0; j < numColIndices; j++) {
1201: printf(" %g", array[i*numColIndices+j]);
1202: }
1203: printf("\n");
1204: }
1205: }
1206: MatSetValues(A, numRowIndices, rowIndices, numColIndices, colIndices, array, mode);
1207: if (ierr) {
1208: printf("[%d]ERROR in updateOperator: points %d %d\n", colSection->commRank(), rowE, colE);
1209: for(int i = 0; i < numRowIndices; i++) {
1210: printf("[%d]mat row indices[%d] = %d\n", rowSection->commRank(), i, rowIndices[i]);
1211: }
1212: for(int i = 0; i < numColIndices; i++) {
1213: printf("[%d]mat col indices[%d] = %d\n", colSection->commRank(), i, colIndices[i]);
1214: }
1215:
1216: }
1218: return(0);
1219: }
1223: /*@C
1224: assembleMatrix - Insert values into a matrix
1226: Collective on A
1228: Input Parameters:
1229: + A - the matrix
1230: . e - The element number
1231: . v - The values
1232: - mode - either ADD_VALUES or INSERT_VALUES, where
1233: ADD_VALUES adds values to any existing entries, and
1234: INSERT_VALUES replaces existing entries with new values
1236: Level: beginner
1238: .seealso: MatSetOption()
1239: @*/
1240: PetscErrorCode assembleMatrix(Mat A, PetscInt e, PetscScalar v[], InsertMode mode)
1241: {
1242: Mesh mesh;
1247: PetscObjectQuery((PetscObject) A, "mesh", (PetscObject *) &mesh);
1248: try {
1249: if (!mesh->mcompat.isNull()) {
1250: Obj<ALECompat::Mesh> m;
1252: MeshCompatGetMesh(mesh, m);
1253: const ALE::Obj<ALECompat::Mesh::topology_type>& topology = m->getTopology();
1254: const ALE::Obj<ALECompat::Mesh::numbering_type>& cNumbering = m->getFactory()->getLocalNumbering(topology, 0, topology->depth());
1255: const ALE::Obj<ALECompat::Mesh::real_section_type>& s = m->getRealSection("default");
1256: const ALE::Obj<ALECompat::Mesh::order_type>& globalOrder = m->getFactory()->getGlobalOrder(topology, 0, "default", s->getAtlas());
1258: if (m->debug()) {
1259: std::cout << "Assembling matrix for element number " << e << " --> point " << cNumbering->getPoint(e) << std::endl;
1260: }
1261: updateOperatorCompat(A, s, globalOrder, cNumbering->getPoint(e), v, mode);
1262: } else {
1263: Obj<ALE::Mesh> m;
1265: MeshGetMesh(mesh, m);
1266: const ALE::Obj<ALE::Mesh::numbering_type>& cNumbering = m->getFactory()->getLocalNumbering(m, m->depth());
1267: const ALE::Obj<ALE::Mesh::real_section_type>& s = m->getRealSection("default");
1268: const ALE::Obj<ALE::Mesh::order_type>& globalOrder = m->getFactory()->getGlobalOrder(m, "default", s);
1270: if (m->debug()) {
1271: std::cout << "Assembling matrix for element number " << e << " --> point " << cNumbering->getPoint(e) << std::endl;
1272: }
1273: updateOperator(A, m, s, globalOrder, cNumbering->getPoint(e), v, mode);
1274: }
1275: } catch (ALE::Exception e) {
1276: std::cout << e.msg() << std::endl;
1277: }
1279: return(0);
1280: }
1284: template<typename Atlas>
1285: PetscErrorCode preallocateOperator(const ALE::Obj<ALE::Mesh>& mesh, const int bs, const ALE::Obj<Atlas>& atlas, const ALE::Obj<ALE::Mesh::order_type>& globalOrder, Mat A)
1286: {
1287: typedef ALE::SieveAlg<ALE::Mesh> sieve_alg_type;
1288: MPI_Comm comm = mesh->comm();
1289: const ALE::Obj<ALE::Mesh> adjBundle = new ALE::Mesh(comm, mesh->debug());
1290: const ALE::Obj<ALE::Mesh::sieve_type> adjGraph = new ALE::Mesh::sieve_type(comm, mesh->debug());
1291: PetscInt numLocalRows, firstRow;
1292: PetscInt *dnz, *onz;
1296: adjBundle->setSieve(adjGraph);
1297: numLocalRows = globalOrder->getLocalSize();
1298: firstRow = globalOrder->getGlobalOffsets()[mesh->commRank()];
1299: PetscMalloc2(numLocalRows, PetscInt, &dnz, numLocalRows, PetscInt, &onz);
1300: /* Create local adjacency graph */
1301: /* In general, we need to get FIAT info that attaches dual basis vectors to sieve points */
1302: const typename Atlas::chart_type& chart = atlas->getChart();
1304: for(typename Atlas::chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
1305: const Obj<typename sieve_alg_type::supportArray>& star = sieve_alg_type::star(mesh, *c_iter);
1307: for(typename sieve_alg_type::supportArray::const_iterator s_iter = star->begin(); s_iter != star->end(); ++s_iter) {
1308: const Obj<typename sieve_alg_type::coneArray>& closure = sieve_alg_type::closure(mesh, *s_iter);
1310: for(typename sieve_alg_type::coneArray::const_iterator cl_iter = closure->begin(); cl_iter != closure->end(); ++cl_iter) {
1311: adjGraph->addCone(*cl_iter, *c_iter);
1312: }
1313: }
1314: }
1315: /* Distribute adjacency graph */
1316: adjBundle->constructOverlap();
1317: typedef typename ALE::Mesh::sieve_type::point_type point_type;
1318: typedef typename ALE::Mesh::send_overlap_type send_overlap_type;
1319: typedef typename ALE::Mesh::recv_overlap_type recv_overlap_type;
1320: typedef typename ALE::Field<send_overlap_type, int, ALE::Section<point_type, point_type> > send_section_type;
1321: typedef typename ALE::Field<recv_overlap_type, int, ALE::Section<point_type, point_type> > recv_section_type;
1322: const Obj<send_overlap_type>& vertexSendOverlap = mesh->getSendOverlap();
1323: const Obj<recv_overlap_type>& vertexRecvOverlap = mesh->getRecvOverlap();
1324: const Obj<send_overlap_type> nbrSendOverlap = new send_overlap_type(comm, mesh->debug());
1325: const Obj<recv_overlap_type> nbrRecvOverlap = new recv_overlap_type(comm, mesh->debug());
1326: const Obj<send_section_type> sendSection = new send_section_type(comm, mesh->debug());
1327: const Obj<recv_section_type> recvSection = new recv_section_type(comm, sendSection->getTag(), mesh->debug());
1329: ALE::Distribution<ALE::Mesh>::coneCompletion(vertexSendOverlap, vertexRecvOverlap, adjBundle, sendSection, recvSection);
1330: /* Distribute indices for new points */
1331: ALE::Distribution<ALE::Mesh>::updateOverlap(sendSection, recvSection, nbrSendOverlap, nbrRecvOverlap);
1332: mesh->getFactory()->completeOrder(globalOrder, nbrSendOverlap, nbrRecvOverlap, true);
1333: /* Read out adjacency graph */
1334: const ALE::Obj<ALE::Mesh::sieve_type> graph = adjBundle->getSieve();
1336: PetscMemzero(dnz, numLocalRows/bs * sizeof(PetscInt));
1337: PetscMemzero(onz, numLocalRows/bs * sizeof(PetscInt));
1338: for(typename Atlas::chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
1339: const typename Atlas::point_type& point = *c_iter;
1341: if (globalOrder->isLocal(point)) {
1342: const ALE::Obj<ALE::Mesh::sieve_type::traits::coneSequence>& adj = graph->cone(point);
1343: const ALE::Mesh::order_type::value_type& rIdx = globalOrder->restrictPoint(point)[0];
1344: const int row = rIdx.prefix;
1345: const int rSize = rIdx.index/bs;
1347: if (rSize == 0) continue;
1348: for(ALE::Mesh::sieve_type::traits::coneSequence::iterator v_iter = adj->begin(); v_iter != adj->end(); ++v_iter) {
1349: const ALE::Mesh::point_type& neighbor = *v_iter;
1350: const ALE::Mesh::order_type::value_type& cIdx = globalOrder->restrictPoint(neighbor)[0];
1351: const int& cSize = cIdx.index/bs;
1353: if (cSize > 0) {
1354: if (globalOrder->isLocal(neighbor)) {
1355: for(int r = 0; r < rSize; ++r) {dnz[(row - firstRow)/bs + r] += cSize;}
1356: } else {
1357: for(int r = 0; r < rSize; ++r) {onz[(row - firstRow)/bs + r] += cSize;}
1358: }
1359: }
1360: }
1361: }
1362: }
1363: if (mesh->debug()) {
1364: int rank = mesh->commRank();
1365: for(int r = 0; r < numLocalRows/bs; r++) {
1366: std::cout << "["<<rank<<"]: dnz["<<r<<"]: " << dnz[r] << " onz["<<r<<"]: " << onz[r] << std::endl;
1367: }
1368: }
1369: MatSeqAIJSetPreallocation(A, 0, dnz);
1370: MatMPIAIJSetPreallocation(A, 0, dnz, 0, onz);
1371: MatSeqBAIJSetPreallocation(A, bs, 0, dnz);
1372: MatMPIBAIJSetPreallocation(A, bs, 0, dnz, 0, onz);
1373: PetscFree2(dnz, onz);
1374: MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR);
1375: return(0);
1376: }
1380: PetscErrorCode preallocateMatrix(const ALE::Obj<ALE::Mesh>& mesh, const int bs, const ALE::Obj<ALE::Mesh::real_section_type::atlas_type>& atlas, const ALE::Obj<ALE::Mesh::order_type>& globalOrder, Mat A)
1381: {
1382: return preallocateOperator(mesh, bs, atlas, globalOrder, A);
1383: }
1385: /******************************** C Wrappers **********************************/
1389: PetscErrorCode WriteVTKHeader(PetscViewer viewer)
1390: {
1391: return VTKViewer::writeHeader(viewer);
1392: }
1396: PetscErrorCode WriteVTKVertices(Mesh mesh, PetscViewer viewer)
1397: {
1398: ALE::Obj<ALE::Mesh> m;
1401: MeshGetMesh(mesh, m);
1402: return VTKViewer::writeVertices(m, viewer);
1403: }
1407: PetscErrorCode WriteVTKElements(Mesh mesh, PetscViewer viewer)
1408: {
1409: ALE::Obj<ALE::Mesh> m;
1412: MeshGetMesh(mesh, m);
1413: return VTKViewer::writeElements(m, viewer);
1414: }
1418: PetscErrorCode WritePCICEVertices(Mesh mesh, PetscViewer viewer)
1419: {
1420: ALE::Obj<ALE::Mesh> m;
1423: MeshGetMesh(mesh, m);
1424: return ALE