Actual source code: mesh.c
petsc-3.3-p7 2013-05-11
1: #include <petsc-private/meshimpl.h> /*I "petscdmmesh.h" I*/
2: #include <petscdmmesh_viewers.hh>
3: #include <petscdmmesh_formats.hh>
5: /* Logging support */
6: PetscLogEvent DMMesh_View, DMMesh_GetGlobalScatter, DMMesh_restrictVector, DMMesh_assembleVector, DMMesh_assembleVectorComplete, DMMesh_assembleMatrix, DMMesh_updateOperator;
8: ALE::MemoryLogger Petsc_MemoryLogger;
10: EXTERN_C_BEGIN
13: /*
14: Private routine to delete internal tag storage when a communicator is freed.
16: This is called by MPI, not by users.
18: Note: this is declared extern "C" because it is passed to MPI_Keyval_create
20: we do not use PetscFree() since it is unsafe after PetscFinalize()
21: */
22: PetscMPIInt DMMesh_DelTag(MPI_Comm comm,PetscMPIInt keyval,void* attr_val,void* extra_state)
23: {
24: free(attr_val);
25: return(MPI_SUCCESS);
26: }
27: EXTERN_C_END
31: PetscErrorCode DMMeshFinalize()
32: {
34: PETSC_MESH_TYPE::MeshNumberingFactory::singleton(0, 0, true);
35: return(0);
36: }
40: /*@C
41: DMMeshGetMesh - Gets the internal mesh object
43: Not collective
45: Input Parameter:
46: . mesh - the mesh object
48: Output Parameter:
49: . m - the internal mesh object
51: Level: advanced
53: .seealso DMMeshCreate(), DMMeshSetMesh()
54: @*/
55: PetscErrorCode DMMeshGetMesh(DM dm, ALE::Obj<PETSC_MESH_TYPE>& m)
56: {
57: DM_Mesh *mesh = (DM_Mesh*) dm->data;
61: if (mesh->useNewImpl) {SETERRQ(((PetscObject) dm)->comm, PETSC_ERR_ARG_WRONG, "This method is only valid for C++ implementation meshes.");}
62: m = mesh->m;
63: return(0);
64: }
68: /*@C
69: DMMeshSetMesh - Sets the internal mesh object
71: Not collective
73: Input Parameters:
74: + mesh - the mesh object
75: - m - the internal mesh object
77: Level: advanced
79: .seealso DMMeshCreate(), DMMeshGetMesh()
80: @*/
81: PetscErrorCode DMMeshSetMesh(DM dm, const ALE::Obj<PETSC_MESH_TYPE>& m)
82: {
83: DM_Mesh *mesh = (DM_Mesh *) dm->data;
88: mesh->m = m;
89: VecScatterDestroy(&mesh->globalScatter);
90: return(0);
91: }
95: PetscErrorCode DMMeshView_Sieve_Ascii(const ALE::Obj<PETSC_MESH_TYPE>& mesh, PetscViewer viewer)
96: {
97: PetscViewerFormat format;
98: PetscErrorCode ierr;
101: PetscViewerGetFormat(viewer, &format);
102: if (format == PETSC_VIEWER_ASCII_VTK) {
103: VTKViewer::writeHeader(mesh, viewer);
104: VTKViewer::writeVertices(mesh, viewer);
105: VTKViewer::writeElements(mesh, viewer);
106: const ALE::Obj<PETSC_MESH_TYPE::int_section_type>& p = mesh->getIntSection("Partition");
107: const ALE::Obj<PETSC_MESH_TYPE::label_sequence>& cells = mesh->heightStratum(0);
108: const PETSC_MESH_TYPE::label_sequence::iterator end = cells->end();
109: const int rank = mesh->commRank();
111: p->setChart(PETSC_MESH_TYPE::int_section_type::chart_type(*cells));
112: p->setFiberDimension(cells, 1);
113: p->allocatePoint();
114: for(PETSC_MESH_TYPE::label_sequence::iterator c_iter = cells->begin(); c_iter != end; ++c_iter) {
115: p->updatePoint(*c_iter, &rank);
116: }
117: PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_VTK_CELL);
118: SectionView_Sieve_Ascii(mesh, p, "Partition", viewer);
119: PetscViewerPopFormat(viewer);
120: } else if (format == PETSC_VIEWER_ASCII_PCICE) {
121: char *filename;
122: char coordFilename[2048];
123: PetscBool isConnect;
124: size_t len;
126: PetscViewerFileGetName(viewer, (const char **) &filename);
127: PetscStrlen(filename, &len);
128: PetscStrcmp(&(filename[len-5]), ".lcon", &isConnect);
129: if (!isConnect) {
130: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG, "Invalid element connectivity filename: %s", filename);
131: }
132: ALE::PCICE::Viewer::writeElements(mesh, viewer);
133: PetscStrncpy(coordFilename, filename, len-5);
134: coordFilename[len-5] = '\0';
135: PetscStrcat(coordFilename, ".nodes");
136: PetscViewerFileSetName(viewer, coordFilename);
137: ALE::PCICE::Viewer::writeVertices(mesh, viewer);
138: } else if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
139: mesh->view("Mesh");
140: } else {
141: PetscInt dim = mesh->getDimension();
142: PetscInt size = mesh->commSize();
143: PetscInt locDepth = mesh->depth(), depth;
144: PetscInt num = 0;
145: PetscInt *sizes;
147: PetscMalloc(size * sizeof(PetscInt), &sizes);
148: PetscViewerASCIIPrintf(viewer, "Mesh in %d dimensions:\n", dim);
149: MPI_Allreduce(&locDepth, &depth, 1, MPIU_INT, MPI_MAX, mesh->comm());
150: if (depth == 1) {
151: num = mesh->depthStratum(0)->size();
152: MPI_Gather(&num, 1, MPIU_INT, sizes, 1, MPIU_INT, 0, mesh->comm());
153: PetscViewerASCIIPrintf(viewer, " %d-cells:", 0);
154: for(PetscInt p = 0; p < size; ++p) {PetscViewerASCIIPrintf(viewer, " %d", sizes[p]);}
155: PetscViewerASCIIPrintf(viewer, "\n");
156: num = mesh->heightStratum(0)->size();
157: MPI_Gather(&num, 1, MPIU_INT, sizes, 1, MPIU_INT, 0, mesh->comm());
158: PetscViewerASCIIPrintf(viewer, " %d-cells:", dim);
159: for(PetscInt p = 0; p < size; ++p) {PetscViewerASCIIPrintf(viewer, " %d", sizes[p]);}
160: PetscViewerASCIIPrintf(viewer, "\n");
161: } else {
162: for(int d = 0; d <= dim; d++) {
163: num = mesh->depthStratum(d)->size();
164: MPI_Gather(&num, 1, MPIU_INT, sizes, 1, MPIU_INT, 0, mesh->comm());
165: PetscViewerASCIIPrintf(viewer, " %d-cells:", d);
166: for(PetscInt p = 0; p < size; ++p) {PetscViewerASCIIPrintf(viewer, " %d", sizes[p]);}
167: PetscViewerASCIIPrintf(viewer, "\n");
168: }
169: }
170: PetscFree(sizes);
171: }
172: PetscViewerFlush(viewer);
173: return(0);
174: }
179: PetscErrorCode DMMeshView_Sieve_Binary(const ALE::Obj<PETSC_MESH_TYPE>& mesh, PetscViewer viewer)
180: {
181: char *filename;
182: PetscErrorCode ierr;
185: PetscViewerFileGetName(viewer, (const char **) &filename);
186: ALE::MeshSerializer::writeMesh(filename, *mesh);
187: return(0);
188: }
192: PetscErrorCode DMMeshView_Sieve(const ALE::Obj<PETSC_MESH_TYPE>& mesh, PetscViewer viewer)
193: {
194: PetscBool iascii, isbinary, isdraw;
198: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
199: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERBINARY, &isbinary);
200: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERDRAW, &isdraw);
202: if (iascii){
203: DMMeshView_Sieve_Ascii(mesh, viewer);
204: } else if (isbinary) {
205: DMMeshView_Sieve_Binary(mesh, viewer);
206: } else if (isdraw){
207: SETERRQ(((PetscObject)viewer)->comm,PETSC_ERR_SUP, "Draw viewer not implemented for DMMesh");
208: } else {
209: SETERRQ1(((PetscObject)viewer)->comm,PETSC_ERR_SUP,"Viewer type %s not supported by this mesh object", ((PetscObject)viewer)->type_name);
210: }
211: return(0);
212: }
216: PetscErrorCode DMMeshView_Mesh_Ascii(DM dm, PetscViewer viewer)
217: {
218: DM_Mesh *mesh = (DM_Mesh *) dm->data;
219: PetscViewerFormat format;
220: PetscErrorCode ierr;
223: PetscViewerGetFormat(viewer, &format);
224: if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
225: const char *name;
226: PetscInt maxConeSize, maxSupportSize;
227: PetscInt pStart, pEnd, p;
228: PetscMPIInt rank;
230: MPI_Comm_rank(((PetscObject) dm)->comm, &rank);
231: PetscObjectGetName((PetscObject) dm, &name);
232: DMMeshGetChart(dm, &pStart, &pEnd);
233: DMMeshGetMaxSizes(dm, &maxConeSize, &maxSupportSize);
234: PetscViewerASCIISynchronizedAllow(viewer, PETSC_TRUE);
235: PetscViewerASCIIPrintf(viewer, "Mesh '%s':\n", name);
236: PetscViewerASCIISynchronizedPrintf(viewer, "Max sizes cone: %d support: %d\n", maxConeSize, maxSupportSize);
237: PetscViewerASCIIPrintf(viewer, "orientation is missing\n", name);
238: PetscViewerASCIIPrintf(viewer, "cap --> base:\n", name);
239: for(p = pStart; p < pEnd; ++p) {
240: PetscInt dof, off, s;
242: PetscSectionGetDof(mesh->supportSection, p, &dof);
243: PetscSectionGetOffset(mesh->supportSection, p, &off);
244: for(s = off; s < off+dof; ++s) {
245: PetscViewerASCIISynchronizedPrintf(viewer, "[%d]: %d ----> %d\n", rank, p, mesh->supports[s]);
246: }
247: }
248: PetscViewerFlush(viewer);
249: PetscViewerASCIIPrintf(viewer, "base <-- cap:\n", name);
250: for(p = pStart; p < pEnd; ++p) {
251: PetscInt dof, off, c;
253: PetscSectionGetDof(mesh->coneSection, p, &dof);
254: PetscSectionGetOffset(mesh->coneSection, p, &off);
255: for(c = off; c < off+dof; ++c) {
256: PetscViewerASCIISynchronizedPrintf(viewer, "[%d]: %d <---- %d\n", rank, p, mesh->cones[c]);
257: /* PetscViewerASCIISynchronizedPrintf(viewer, "[%d]: %d <---- %d: %d\n", rank, p, mesh->cones[c], mesh->coneOrientations[c]); */
258: }
259: }
260: PetscViewerFlush(viewer);
261: PetscSectionVecView(mesh->coordSection, mesh->coordinates, viewer);
262: } else {
263: MPI_Comm comm = ((PetscObject) dm)->comm;
264: PetscInt *sizes;
265: PetscInt depth, dim, d;
266: PetscInt pStart, pEnd, p;
267: PetscMPIInt size;
269: MPI_Comm_size(comm, &size);
270: DMMeshGetDimension(dm, &dim);
271: PetscViewerASCIIPrintf(viewer, "Mesh in %d dimensions:\n", dim);
272: MPI_Allreduce(&depth, &depth, 1, MPIU_INT, MPI_MAX, comm);
273: PetscMalloc(size * sizeof(PetscInt), &sizes);
274: DMMeshGetLabelSize(dm, "depth", &depth);
275: if (depth == 2) {
276: DMMeshGetDepthStratum(dm, 0, &pStart, &pEnd);
277: pEnd = pEnd - pStart;
278: MPI_Gather(&pEnd, 1, MPIU_INT, sizes, 1, MPIU_INT, 0, comm);
279: PetscViewerASCIIPrintf(viewer, " %d-cells:", 0);
280: for(p = 0; p < size; ++p) {PetscViewerASCIIPrintf(viewer, " %d", sizes[p]);}
281: PetscViewerASCIIPrintf(viewer, "\n");
282: DMMeshGetHeightStratum(dm, 0, &pStart, &pEnd);
283: pEnd = pEnd - pStart;
284: MPI_Gather(&pEnd, 1, MPIU_INT, sizes, 1, MPIU_INT, 0, comm);
285: PetscViewerASCIIPrintf(viewer, " %d-cells:", dim);
286: for(p = 0; p < size; ++p) {PetscViewerASCIIPrintf(viewer, " %d", sizes[p]);}
287: PetscViewerASCIIPrintf(viewer, "\n");
288: } else {
289: for(d = 0; d <= dim; d++) {
290: DMMeshGetDepthStratum(dm, d, &pStart, &pEnd);
291: pEnd = pEnd - pStart;
292: MPI_Gather(&pEnd, 1, MPIU_INT, sizes, 1, MPIU_INT, 0, comm);
293: PetscViewerASCIIPrintf(viewer, " %d-cells:", d);
294: for(p = 0; p < size; ++p) {PetscViewerASCIIPrintf(viewer, " %d", sizes[p]);}
295: PetscViewerASCIIPrintf(viewer, "\n");
296: }
297: }
298: PetscFree(sizes);
299: }
300: return(0);
301: }
305: PetscErrorCode DMView_Mesh(DM dm, PetscViewer viewer)
306: {
307: DM_Mesh *mesh = (DM_Mesh *) dm->data;
313: if (mesh->useNewImpl) {
314: PetscBool iascii, isbinary;
316: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
317: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERBINARY, &isbinary);
319: if (iascii) {
320: DMMeshView_Mesh_Ascii(dm, viewer);
321: #if 0
322: } else if (isbinary) {
323: DMMeshView_Mesh_Binary(dm, viewer);
324: #endif
325: } else {
326: SETERRQ1(((PetscObject)viewer)->comm,PETSC_ERR_SUP,"Viewer type %s not supported by this mesh object", ((PetscObject)viewer)->type_name);
327: }
328: } else {
329: DMMeshView_Sieve(mesh->m, viewer);
330: }
331: return(0);
332: }
336: /*@
337: DMMeshLoad - 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 DMView()
350: @*/
351: PetscErrorCode DMMeshLoad(PetscViewer viewer, DM dm)
352: {
353: DM_Mesh *mesh = (DM_Mesh *) dm->data;
354: char *filename;
358: if (!mesh->m) {
359: MPI_Comm comm;
361: PetscObjectGetComm((PetscObject) viewer, &comm);
362: ALE::Obj<PETSC_MESH_TYPE> m = new PETSC_MESH_TYPE(comm, 1);
363: DMMeshSetMesh(dm, m);
364: }
365: PetscViewerFileGetName(viewer, (const char **) &filename);
366: ALE::MeshSerializer::loadMesh(filename, *mesh->m);
367: return(0);
368: }
372: PetscErrorCode GetAdjacentDof_Private()
373: {
375: return(0);
376: }
380: PetscErrorCode DMCreateMatrix_Mesh(DM dm, const MatType mtype, Mat *J)
381: {
382: DM_Mesh *mesh = (DM_Mesh *) dm->data;
383: ISLocalToGlobalMapping ltog;
384: PetscErrorCode ierr;
387: #ifndef PETSC_USE_DYNAMIC_LIBRARIES
388: MatInitializePackage(PETSC_NULL);
389: #endif
390: if (!mtype) mtype = MATAIJ;
391: if (mesh->useNewImpl) {
392: SETERRQ(((PetscObject) dm)->comm, PETSC_ERR_SUP, "Not supported");
393: } else {
394: ALE::Obj<PETSC_MESH_TYPE> m;
395: ALE::Obj<PETSC_MESH_TYPE::real_section_type> s;
396: SectionReal section;
397: PetscBool flag;
398: DMMeshHasSectionReal(dm, "default", &flag);
399: if (!flag) SETERRQ(((PetscObject) dm)->comm, PETSC_ERR_ARG_WRONGSTATE, "Must set default section");
400: DMMeshGetSectionReal(dm, "default", §ion);
401: DMMeshGetMesh(dm, m);
402: SectionRealGetSection(section, s);
403: try {
404: DMMeshCreateMatrix(m, s, mtype, J, -1, !dm->prealloc_only);
405: } catch(ALE::Exception e) {
406: SETERRQ(((PetscObject) dm)->comm, PETSC_ERR_LIB, e.message());
407: }
408: SectionRealDestroy(§ion);
409: }
410: PetscObjectCompose((PetscObject) *J, "DM", (PetscObject) dm);
411: DMGetLocalToGlobalMapping(dm, <og);
412: MatSetLocalToGlobalMapping(*J, ltog, ltog);
413: return(0);
414: }
418: /*@C
419: DMMeshCreateMatrix - Creates a matrix with the correct parallel layout required for
420: computing the Jacobian on a function defined using the information in the Section.
422: Collective on DMMesh
424: Input Parameters:
425: + mesh - the mesh object
426: . section - the section which determines data layout
427: - mtype - Supported types are MATSEQAIJ, MATMPIAIJ, MATSEQBAIJ, MATMPIBAIJ, MATSEQSBAIJ, MATMPISBAIJ,
428: or any type which inherits from one of these (such as MATAIJ, MATLUSOL, etc.).
430: Output Parameter:
431: . J - matrix with the correct nonzero preallocation
432: (obviously without the correct Jacobian values)
434: Level: advanced
436: Notes: This properly preallocates the number of nonzeros in the sparse matrix so you
437: do not need to do it yourself.
439: .seealso ISColoringView(), ISColoringGetIS(), MatFDColoringCreate(), DMDASetBlockFills()
440: @*/
441: PetscErrorCode DMMeshCreateMatrix(DM dm, SectionReal section, const MatType mtype, Mat *J)
442: {
443: ALE::Obj<PETSC_MESH_TYPE> m;
444: ALE::Obj<PETSC_MESH_TYPE::real_section_type> s;
448: #ifndef PETSC_USE_DYNAMIC_LIBRARIES
449: MatInitializePackage(PETSC_NULL);
450: #endif
451: if (!mtype) mtype = MATAIJ;
452: DMMeshGetMesh(dm, m);
453: SectionRealGetSection(section, s);
454: try {
455: DMMeshCreateMatrix(m, s, mtype, J, -1, !dm->prealloc_only);
456: } catch(ALE::Exception e) {
457: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB, e.message());
458: }
459: PetscObjectCompose((PetscObject) *J, "DM", (PetscObject) dm);
460: return(0);
461: }
465: /*@
466: DMMeshCreateVector - Creates a global vector matching the input section
468: Collective on DMMesh
470: Input Parameters:
471: + mesh - the DMMesh
472: - section - the Section
474: Output Parameter:
475: . vec - the global vector
477: Level: advanced
479: Notes: The vector can safely be destroyed using VecDestroy().
480: .seealso DMMeshCreate()
481: @*/
482: PetscErrorCode DMMeshCreateVector(DM mesh, SectionReal section, Vec *vec)
483: {
484: ALE::Obj<PETSC_MESH_TYPE> m;
485: ALE::Obj<PETSC_MESH_TYPE::real_section_type> s;
489: DMMeshGetMesh(mesh, m);
490: SectionRealGetSection(section, s);
491: const ALE::Obj<PETSC_MESH_TYPE::order_type>& order = m->getFactory()->getGlobalOrder(m, s->getName(), s);
493: VecCreate(m->comm(), vec);
494: VecSetSizes(*vec, order->getLocalSize(), order->getGlobalSize());
495: VecSetFromOptions(*vec);
496: return(0);
497: }
501: PetscErrorCode DMDestroy_Mesh(DM dm)
502: {
503: DM_Mesh *mesh = (DM_Mesh *) dm->data;
504: SieveLabel next = mesh->labels;
508: mesh->m = PETSC_NULL;
509: PetscSectionDestroy(&mesh->defaultSection);
510: VecScatterDestroy(&mesh->globalScatter);
512: /* NEW_MESH_IMPL */
513: PetscSFDestroy(&mesh->sf);
514: PetscSectionDestroy(&mesh->coneSection);
515: PetscFree(mesh->cones);
516: PetscSectionDestroy(&mesh->supportSection);
517: PetscFree(mesh->supports);
518: PetscSectionDestroy(&mesh->coordSection);
519: VecDestroy(&mesh->coordinates);
520: PetscFree2(mesh->meetTmpA,mesh->meetTmpB);
521: PetscFree2(mesh->joinTmpA,mesh->joinTmpB);
522: PetscFree2(mesh->closureTmpA,mesh->closureTmpB);
523: while(next) {
524: SieveLabel tmp;
526: PetscFree(next->name);
527: PetscFree(next->stratumValues);
528: PetscFree(next->stratumOffsets);
529: PetscFree(next->stratumSizes);
530: PetscFree(next->points);
531: tmp = next->next;
532: PetscFree(next);
533: next = tmp;
534: }
535: return(0);
536: }
540: PetscErrorCode DMCreateGlobalVector_Mesh(DM dm, Vec *gvec)
541: {
542: DM_Mesh *mesh = (DM_Mesh *) dm->data;
543: PetscInt localSize, globalSize;
547: if (mesh->useNewImpl) {
548: PetscSection s;
550: /* DOES NOT WORK IN PARALLEL, CHANGE TO DMCOMPLEX */
551: DMMeshGetDefaultSection(dm, &s);
552: PetscSectionGetStorageSize(s, &localSize);
553: globalSize = PETSC_DETERMINE;
554: } else {
555: ALE::Obj<PETSC_MESH_TYPE> m;
556: PetscBool flag;
557: DMMeshGetMesh(dm, m);
558: DMMeshHasSectionReal(dm, "default", &flag);
559: if (!flag) SETERRQ(((PetscObject) dm)->comm,PETSC_ERR_ARG_WRONGSTATE, "Must set default section");
560: const ALE::Obj<PETSC_MESH_TYPE::order_type>& order = m->getFactory()->getGlobalOrder(m, "default", m->getRealSection("default"));
562: localSize = order->getLocalSize();
563: globalSize = order->getGlobalSize();
564: }
565: VecCreate(((PetscObject) dm)->comm, gvec);
566: VecSetSizes(*gvec, localSize, globalSize);
567: VecSetFromOptions(*gvec);
568: PetscObjectCompose((PetscObject) *gvec, "DM", (PetscObject) dm);
569: return(0);
570: }
574: PetscErrorCode DMCreateLocalVector_Mesh(DM dm, Vec *lvec)
575: {
576: DM_Mesh *mesh = (DM_Mesh *) dm->data;
577: PetscInt size;
581: if (mesh->useNewImpl) {
582: PetscSection s;
584: DMMeshGetDefaultSection(dm, &s);
585: PetscSectionGetStorageSize(s, &size);
586: } else {
587: ALE::Obj<PETSC_MESH_TYPE> m;
588: DMMeshGetMesh(dm, m);
589: PetscBool flag;
590: DMMeshGetMesh(dm, m);
591: DMMeshHasSectionReal(dm, "default", &flag);
592: if (!flag) SETERRQ(((PetscObject) dm)->comm,PETSC_ERR_ARG_WRONGSTATE, "Must set default section");
593: size = m->getRealSection("default")->getStorageSize();
594: }
595: VecCreate(PETSC_COMM_SELF, lvec);
596: VecSetSizes(*lvec, size, size);
597: VecSetFromOptions(*lvec);
598: PetscObjectCompose((PetscObject) *lvec, "DM", (PetscObject) dm);
599: return(0);
600: }
604: PetscErrorCode DMCreateLocalToGlobalMapping_Mesh(DM dm)
605: {
606: ALE::Obj<PETSC_MESH_TYPE> m;
607: ALE::Obj<PETSC_MESH_TYPE::real_section_type> s;
608: SectionReal section;
609: PetscBool flag;
613: DMMeshHasSectionReal(dm, "default", &flag);
614: if (!flag) SETERRQ(((PetscObject) dm)->comm,PETSC_ERR_ARG_WRONGSTATE, "Must set default section");
615: DMMeshGetSectionReal(dm ,"default", §ion);
616: DMMeshGetMesh(dm, m);
617: SectionRealGetSection(section, s);
618: const ALE::Obj<PETSC_MESH_TYPE::order_type>& globalOrder = m->getFactory()->getGlobalOrder(m, s->getName(), s);
619: PetscInt *ltog;
621: PetscMalloc(s->size() * sizeof(PetscInt), <og); // We want the local+overlap size
622: for(PetscInt p = s->getChart().min(), l = 0; p < s->getChart().max(); ++p) {
623: PetscInt g = globalOrder->getIndex(p);
625: for(PetscInt c = 0; c < s->getConstrainedFiberDimension(p); ++c, ++l) {
626: ltog[l] = g+c;
627: }
628: }
629: ISLocalToGlobalMappingCreate(PETSC_COMM_SELF, s->size(), ltog, PETSC_OWN_POINTER, &dm->ltogmap);
630: PetscLogObjectParent(dm, dm->ltogmap);
631: SectionRealDestroy(§ion);
632: return(0);
633: }
637: /*@
638: DMMeshCreateGlobalScatter - Create a VecScatter which maps from local, overlapping
639: storage in the Section to a global Vec
641: Collective on DMMesh
643: Input Parameters:
644: + mesh - the mesh object
645: - section - The Scetion which determines data layout
647: Output Parameter:
648: . scatter - the VecScatter
650: Level: advanced
652: .seealso DMDestroy(), DMMeshCreateGlobalRealVector(), DMMeshCreate()
653: @*/
654: PetscErrorCode DMMeshCreateGlobalScatter(DM dm, SectionReal section, VecScatter *scatter)
655: {
656: ALE::Obj<PETSC_MESH_TYPE> m;
657: ALE::Obj<PETSC_MESH_TYPE::real_section_type> s;
661: DMMeshGetMesh(dm, m);
662: SectionRealGetSection(section, s);
663: if (m->hasLabel("marker")) {
664: DMMeshCreateGlobalScatter(m, s, m->getLabel("marker"), scatter);
665: } else {
666: DMMeshCreateGlobalScatter(m, s, scatter);
667: }
668: return(0);
669: }
673: /*@
674: DMMeshGetGlobalScatter - Retrieve the VecScatter which maps from local, overlapping storage in the default Section to a global Vec
676: Collective on DMMesh
678: Input Parameters:
679: . mesh - the mesh object
681: Output Parameter:
682: . scatter - the VecScatter
684: Level: advanced
686: .seealso MeshDestroy(), DMMeshCreateGlobalrealVector(), DMMeshCreate()
687: @*/
688: PetscErrorCode DMMeshGetGlobalScatter(DM dm, VecScatter *scatter)
689: {
690: DM_Mesh *mesh = (DM_Mesh *) dm->data;
696: if (!mesh->globalScatter) {
697: SectionReal section;
699: DMMeshGetSectionReal(dm, "default", §ion);
700: DMMeshCreateGlobalScatter(dm, section, &mesh->globalScatter);
701: SectionRealDestroy(§ion);
702: }
703: *scatter = mesh->globalScatter;
704: return(0);
705: }
709: PetscErrorCode DMGlobalToLocalBegin_Mesh(DM dm, Vec g, InsertMode mode, Vec l)
710: {
711: VecScatter injection;
715: DMMeshGetGlobalScatter(dm, &injection);
716: VecScatterBegin(injection, g, l, mode, SCATTER_REVERSE);
717: return(0);
718: }
722: PetscErrorCode DMGlobalToLocalEnd_Mesh(DM dm, Vec g, InsertMode mode, Vec l)
723: {
724: VecScatter injection;
728: DMMeshGetGlobalScatter(dm, &injection);
729: VecScatterEnd(injection, g, l, mode, SCATTER_REVERSE);
730: return(0);
731: }
735: PetscErrorCode DMLocalToGlobalBegin_Mesh(DM dm, Vec l, InsertMode mode, Vec g)
736: {
737: VecScatter injection;
741: DMMeshGetGlobalScatter(dm, &injection);
742: VecScatterBegin(injection, l, g, mode, SCATTER_FORWARD);
743: return(0);
744: }
748: PetscErrorCode DMLocalToGlobalEnd_Mesh(DM dm, Vec l, InsertMode mode, Vec g)
749: {
750: VecScatter injection;
754: DMMeshGetGlobalScatter(dm, &injection);
755: VecScatterEnd(injection, l, g, mode, SCATTER_FORWARD);
756: return(0);
757: }
761: PetscErrorCode DMMeshGetLocalFunction(DM dm, PetscErrorCode (**lf)(DM, Vec, Vec, void *))
762: {
763: DM_Mesh *mesh = (DM_Mesh *) dm->data;
767: if (lf) *lf = mesh->lf;
768: return(0);
769: }
773: PetscErrorCode DMMeshSetLocalFunction(DM dm, PetscErrorCode (*lf)(DM, Vec, Vec, void *))
774: {
775: DM_Mesh *mesh = (DM_Mesh *) dm->data;
779: mesh->lf = lf;
780: return(0);
781: }
785: PetscErrorCode DMMeshGetLocalJacobian(DM dm, PetscErrorCode (**lj)(DM, Vec, Mat, void *))
786: {
787: DM_Mesh *mesh = (DM_Mesh *) dm->data;
791: if (lj) *lj = mesh->lj;
792: return(0);
793: }
797: PetscErrorCode DMMeshSetLocalJacobian(DM dm, PetscErrorCode (*lj)(DM, Vec, Mat, void *))
798: {
799: DM_Mesh *mesh = (DM_Mesh *) dm->data;
803: mesh->lj = lj;
804: return(0);
805: }
809: /*@
810: DMMeshGetDimension - Return the topological mesh dimension
812: Not collective
814: Input Parameter:
815: . mesh - The DMMesh
817: Output Parameter:
818: . dim - The topological mesh dimension
820: Level: beginner
822: .seealso: DMMeshCreate()
823: @*/
824: PetscErrorCode DMMeshGetDimension(DM dm, PetscInt *dim)
825: {
826: DM_Mesh *mesh = (DM_Mesh *) dm->data;
832: if (mesh->useNewImpl) {
833: *dim = mesh->dim;
834: } else {
835: Obj<PETSC_MESH_TYPE> m;
836: DMMeshGetMesh(dm, m);
837: *dim = m->getDimension();
838: }
839: return(0);
840: }
844: /*@
845: DMMeshSetDimension - Set the topological mesh dimension
847: Collective on mesh
849: Input Parameters:
850: + mesh - The DMMesh
851: - dim - The topological mesh dimension
853: Level: beginner
855: .seealso: DMMeshCreate()
856: @*/
857: PetscErrorCode DMMeshSetDimension(DM dm, PetscInt dim)
858: {
859: DM_Mesh *mesh = (DM_Mesh *) dm->data;
864: if (mesh->useNewImpl) {
865: mesh->dim = dim;
866: } else {
867: SETERRQ(((PetscObject) dm)->comm, PETSC_ERR_SUP, "Cannot reset dimension of C++ mesh");
868: }
869: return(0);
870: }
874: /*@
875: DMMeshGetChart - Return the interval for all mesh points [pStart, pEnd)
877: Not collective
879: Input Parameter:
880: . mesh - The DMMesh
882: Output Parameters:
883: + pStart - The first mesh point
884: - pEnd - The upper bound for mesh points
886: Level: beginner
888: .seealso: DMMeshCreate(), DMMeshSetChart()
889: @*/
890: PetscErrorCode DMMeshGetChart(DM dm, PetscInt *pStart, PetscInt *pEnd)
891: {
892: DM_Mesh *mesh = (DM_Mesh *) dm->data;
897: if (mesh->useNewImpl) {
898: PetscSectionGetChart(mesh->coneSection, pStart, pEnd);
899: } else {
900: Obj<PETSC_MESH_TYPE> m;
901: DMMeshGetMesh(dm, m);
902: if (pStart) {*pStart = m->getSieve()->getChart().min();}
903: if (pEnd) {*pEnd = m->getSieve()->getChart().max();}
904: }
905: return(0);
906: }
910: /*@
911: DMMeshSetChart - Set the interval for all mesh points [pStart, pEnd)
913: Not collective
915: Input Parameters:
916: + mesh - The DMMesh
917: . pStart - The first mesh point
918: - pEnd - The upper bound for mesh points
920: Output Parameters:
922: Level: beginner
924: .seealso: DMMeshCreate(), DMMeshGetChart()
925: @*/
926: PetscErrorCode DMMeshSetChart(DM dm, PetscInt pStart, PetscInt pEnd)
927: {
928: DM_Mesh *mesh = (DM_Mesh *) dm->data;
933: if (mesh->useNewImpl) {
934: PetscSectionSetChart(mesh->coneSection, pStart, pEnd);
935: } else {
936: Obj<PETSC_MESH_TYPE> m;
937: DMMeshGetMesh(dm, m);
938: m->getSieve()->setChart(PETSC_MESH_TYPE::sieve_type::chart_type(pStart, pEnd));
939: }
940: return(0);
941: }
945: /*@
946: DMMeshGetConeSize - Return the number of in-edges for this point in the Sieve DAG
948: Not collective
950: Input Parameters:
951: + mesh - The DMMesh
952: - p - The Sieve point, which must lie in the chart set with DMMeshSetChart()
954: Output Parameter:
955: . size - The cone size for point p
957: Level: beginner
959: .seealso: DMMeshCreate(), DMMeshSetConeSize(), DMMeshSetChart()
960: @*/
961: PetscErrorCode DMMeshGetConeSize(DM dm, PetscInt p, PetscInt *size)
962: {
963: DM_Mesh *mesh = (DM_Mesh *) dm->data;
969: if (mesh->useNewImpl) {
970: PetscSectionGetDof(mesh->coneSection, p, size);
971: } else {
972: Obj<PETSC_MESH_TYPE> m;
973: DMMeshGetMesh(dm, m);
974: *size = m->getSieve()->getConeSize(p);
975: }
976: return(0);
977: }
981: /*@
982: DMMeshSetConeSize - Set the number of in-edges for this point in the Sieve DAG
984: Not collective
986: Input Parameters:
987: + mesh - The DMMesh
988: . p - The Sieve point, which must lie in the chart set with DMMeshSetChart()
989: - size - The cone size for point p
991: Output Parameter:
993: Note:
994: This should be called after DMMeshSetChart().
996: Level: beginner
998: .seealso: DMMeshCreate(), DMMeshGetConeSize(), DMMeshSetChart()
999: @*/
1000: PetscErrorCode DMMeshSetConeSize(DM dm, PetscInt p, PetscInt size)
1001: {
1002: DM_Mesh *mesh = (DM_Mesh *) dm->data;
1007: if (mesh->useNewImpl) {
1008: PetscSectionSetDof(mesh->coneSection, p, size);
1009: mesh->maxConeSize = PetscMax(mesh->maxConeSize, size);
1010: } else {
1011: Obj<PETSC_MESH_TYPE> m;
1012: DMMeshGetMesh(dm, m);
1013: m->getSieve()->setConeSize(p, size);
1014: }
1015: return(0);
1016: }
1020: /*@C
1021: DMMeshGetCone - Return the points on the in-edges for this point in the Sieve DAG
1023: Not collective
1025: Input Parameters:
1026: + mesh - The DMMesh
1027: - p - The Sieve point, which must lie in the chart set with DMMeshSetChart()
1029: Output Parameter:
1030: . cone - An array of points which are on the in-edges for point p
1032: Level: beginner
1034: Note:
1035: This routine is not available in Fortran.
1037: .seealso: DMMeshCreate(), DMMeshSetCone(), DMMeshSetChart()
1038: @*/
1039: PetscErrorCode DMMeshGetCone(DM dm, PetscInt p, const PetscInt *cone[])
1040: {
1041: DM_Mesh *mesh = (DM_Mesh *) dm->data;
1047: if (mesh->useNewImpl) {
1048: PetscInt off;
1050: PetscSectionGetOffset(mesh->coneSection, p, &off);
1051: *cone = &mesh->cones[off];
1052: } else {
1053: Obj<PETSC_MESH_TYPE> m;
1054: DMMeshGetMesh(dm, m);
1055: ALE::ISieveVisitor::PointRetriever<PETSC_MESH_TYPE::sieve_type> v(m->getSieve()->getConeSize(p));
1057: m->getSieve()->cone(p, v);
1058: if (!mesh->meetTmpA) {PetscMalloc2(m->getSieve()->getMaxConeSize(),PetscInt,&mesh->meetTmpA,m->getSieve()->getMaxConeSize(),PetscInt,&mesh->meetTmpB);}
1059: for(size_t c = 0; c < v.getSize(); ++c) {
1060: mesh->meetTmpA[c] = v.getPoints()[c];
1061: }
1062: *cone = mesh->meetTmpA;
1063: }
1064: return(0);
1065: }
1069: /*@
1070: DMMeshSetCone - Set the points on the in-edges for this point in the Sieve DAG
1072: Not collective
1074: Input Parameters:
1075: + mesh - The DMMesh
1076: . p - The Sieve point, which must lie in the chart set with DMMeshSetChart()
1077: - cone - An array of points which are on the in-edges for point p
1079: Output Parameter:
1081: Note:
1082: This should be called after all calls to DMMeshSetConeSize() and DMMeshSetUp().
1084: Level: beginner
1086: .seealso: DMMeshCreate(), DMMeshGetCone(), DMMeshSetChart(), DMMeshSetConeSize(), DMMeshSetUp()
1087: @*/
1088: PetscErrorCode DMMeshSetCone(DM dm, PetscInt p, const PetscInt cone[])
1089: {
1090: DM_Mesh *mesh = (DM_Mesh *) dm->data;
1096: if (mesh->useNewImpl) {
1097: PetscInt pStart, pEnd;
1098: PetscInt dof, off, c;
1100: PetscSectionGetChart(mesh->coneSection, &pStart, &pEnd);
1101: PetscSectionGetDof(mesh->coneSection, p, &dof);
1102: PetscSectionGetOffset(mesh->coneSection, p, &off);
1103: for(c = 0; c < dof; ++c) {
1104: if ((cone[c] < pStart) || (cone[c] >= pEnd)) {SETERRQ3(((PetscObject) dm)->comm, PETSC_ERR_ARG_OUTOFRANGE, "Cone point %d is not in the valid range [%d. %d)", cone[c], pStart, pEnd);}
1105: mesh->cones[off+c] = cone[c];
1106: }
1107: } else {
1108: SETERRQ(((PetscObject) dm)->comm, PETSC_ERR_SUP, "This method does not make sense for the C++ Sieve implementation");
1109: }
1110: return(0);
1111: }
1115: /*@
1116: DMMeshGetSupportSize - Return the number of out-edges for this point in the Sieve DAG
1118: Not collective
1120: Input Parameters:
1121: + mesh - The DMMesh
1122: - p - The Sieve point, which must lie in the chart set with DMMeshSetChart()
1124: Output Parameter:
1125: . size - The support size for point p
1127: Level: beginner
1129: .seealso: DMMeshCreate(), DMMeshSetConeSize(), DMMeshSetChart(), DMMeshGetConeSize()
1130: @*/
1131: PetscErrorCode DMMeshGetSupportSize(DM dm, PetscInt p, PetscInt *size)
1132: {
1133: DM_Mesh *mesh = (DM_Mesh *) dm->data;
1139: if (mesh->useNewImpl) {
1140: PetscSectionGetDof(mesh->supportSection, p, size);
1141: } else {
1142: Obj<PETSC_MESH_TYPE> m;
1143: DMMeshGetMesh(dm, m);
1144: *size = m->getSieve()->getSupportSize(p);
1145: }
1146: return(0);
1147: }
1151: /*@C
1152: DMMeshGetSupport - Return the points on the out-edges for this point in the Sieve DAG
1154: Not collective
1156: Input Parameters:
1157: + mesh - The DMMesh
1158: - p - The Sieve point, which must lie in the chart set with DMMeshSetChart()
1160: Output Parameter:
1161: . support - An array of points which are on the out-edges for point p
1163: Level: beginner
1165: .seealso: DMMeshCreate(), DMMeshSetCone(), DMMeshSetChart(), DMMeshGetCone()
1166: @*/
1167: PetscErrorCode DMMeshGetSupport(DM dm, PetscInt p, const PetscInt *support[])
1168: {
1169: DM_Mesh *mesh = (DM_Mesh *) dm->data;
1175: if (mesh->useNewImpl) {
1176: PetscInt off;
1178: PetscSectionGetOffset(mesh->supportSection, p, &off);
1179: *support = &mesh->supports[off];
1180: } else {
1181: Obj<PETSC_MESH_TYPE> m;
1182: DMMeshGetMesh(dm, m);
1183: ALE::ISieveVisitor::PointRetriever<PETSC_MESH_TYPE::sieve_type> v(m->getSieve()->getSupportSize(p));
1185: m->getSieve()->support(p, v);
1186: if (!mesh->joinTmpA) {PetscMalloc2(m->getSieve()->getMaxSupportSize(),PetscInt,&mesh->joinTmpA,m->getSieve()->getMaxSupportSize(),PetscInt,&mesh->joinTmpB);}
1187: for(size_t s = 0; s < v.getSize(); ++s) {
1188: mesh->joinTmpA[s] = v.getPoints()[s];
1189: }
1190: *support = mesh->joinTmpA;
1191: }
1192: return(0);
1193: }
1197: /*@C
1198: DMMeshGetTransitiveClosure - Return the points on the transitive closure of the in-edges or out-edges for this point in the Sieve DAG
1200: Not collective
1202: Input Parameters:
1203: + mesh - The DMMesh
1204: . p - The Sieve point, which must lie in the chart set with DMMeshSetChart()
1205: - useCone - PETSC_TRUE for in-edges, otherwise use out-edges
1207: Output Parameters:
1208: + numPoints - The number of points in the closure
1209: - points - The points
1211: Level: beginner
1213: .seealso: DMMeshCreate(), DMMeshSetCone(), DMMeshSetChart(), DMMeshGetCone()
1214: @*/
1215: PetscErrorCode DMMeshGetTransitiveClosure(DM dm, PetscInt p, PetscBool useCone, PetscInt *numPoints, const PetscInt *points[])
1216: {
1217: DM_Mesh *mesh = (DM_Mesh *) dm->data;
1222: if (!mesh->closureTmpA) {
1223: PetscInt maxSize;
1224: if (mesh->useNewImpl) {
1225: maxSize = PetscMax(mesh->maxConeSize, mesh->maxSupportSize)+1;
1226: } else {
1227: Obj<PETSC_MESH_TYPE> m;
1228: DMMeshGetMesh(dm, m);
1229: maxSize = PetscMax(m->getSieve()->getMaxConeSize(), m->getSieve()->getMaxSupportSize())+1;
1230: }
1231: PetscMalloc2(maxSize,PetscInt,&mesh->closureTmpA,maxSize,PetscInt,&mesh->closureTmpB);
1232: }
1233: if (mesh->useNewImpl) {
1234: const PetscInt *tmp;
1235: PetscInt tmpSize, t;
1236: PetscInt closureSize = 1;
1238: mesh->closureTmpA[0] = p;
1239: /* This is only 1-level */
1240: if (useCone) {
1241: DMMeshGetConeSize(dm, p, &tmpSize);
1242: DMMeshGetCone(dm, p, &tmp);
1243: } else {
1244: DMMeshGetSupportSize(dm, p, &tmpSize);
1245: DMMeshGetSupport(dm, p, &tmp);
1246: }
1247: for(t = 0; t < tmpSize; ++t) {
1248: mesh->closureTmpA[closureSize++] = tmp[t];
1249: }
1250: if (numPoints) *numPoints = closureSize;
1251: if (points) *points = mesh->closureTmpA;
1252: } else {
1253: Obj<PETSC_MESH_TYPE> m;
1254: DMMeshGetMesh(dm, m);
1255: typedef ALE::ISieveVisitor::TransitiveClosureVisitor<PETSC_MESH_TYPE::sieve_type> visitor_type;
1256: visitor_type::visitor_type nV;
1257: visitor_type cV(*m->getSieve(), nV);
1259: if (useCone) {
1260: m->getSieve()->cone(p, cV);
1261: } else {
1262: cV.setIsCone(false);
1263: m->getSieve()->support(p, cV);
1264: }
1265: int i = 0;
1267: for(std::set<PETSC_MESH_TYPE::point_type>::const_iterator p_iter = cV.getPoints().begin(); p_iter != cV.getPoints().end(); ++p_iter, ++i) {
1268: mesh->closureTmpA[i] = *p_iter;
1269: }
1270: if (numPoints) *numPoints = cV.getPoints().size();
1271: if (points) *points = mesh->closureTmpA;
1272: }
1273: return(0);
1274: }
1278: /*@
1279: DMMeshGetMaxSizes - Return the maximum number of in-edges (cone) and out-edges (support) for any point in the Sieve DAG
1281: Not collective
1283: Input Parameter:
1284: . mesh - The DMMesh
1286: Output Parameters:
1287: + maxConeSize - The maximum number of in-edges
1288: - maxSupportSize - The maximum number of out-edges
1290: Level: beginner
1292: .seealso: DMMeshCreate(), DMMeshSetConeSize(), DMMeshSetChart()
1293: @*/
1294: PetscErrorCode DMMeshGetMaxSizes(DM dm, PetscInt *maxConeSize, PetscInt *maxSupportSize)
1295: {
1296: DM_Mesh *mesh = (DM_Mesh *) dm->data;
1301: if (mesh->useNewImpl) {
1302: if (maxConeSize) *maxConeSize = mesh->maxConeSize;
1303: if (maxSupportSize) *maxSupportSize = mesh->maxSupportSize;
1304: } else {
1305: Obj<PETSC_MESH_TYPE> m;
1306: DMMeshGetMesh(dm, m);
1307: if (maxConeSize) *maxConeSize = m->getSieve()->getMaxConeSize();
1308: if (maxSupportSize) *maxSupportSize = m->getSieve()->getMaxSupportSize();
1309: }
1310: return(0);
1311: }
1315: /*@
1316: DMMeshSetUp - Allocate space for the Sieve DAG
1318: Not collective
1320: Input Parameter:
1321: . mesh - The DMMesh
1323: Output Parameter:
1325: Note:
1326: This should be called after DMMeshSetChart() and all calls to DMMeshSetConeSize()
1328: Level: beginner
1330: .seealso: DMMeshCreate(), DMMeshSetChart(), DMMeshSetConeSize()
1331: @*/
1332: PetscErrorCode DMMeshSetUp(DM dm)
1333: {
1334: DM_Mesh *mesh = (DM_Mesh *) dm->data;
1339: if (mesh->useNewImpl) {
1340: PetscInt size;
1342: PetscSectionSetUp(mesh->coneSection);
1343: PetscSectionGetStorageSize(mesh->coneSection, &size);
1344: PetscMalloc(size * sizeof(PetscInt), &mesh->cones);
1345: }
1346: return(0);
1347: }
1351: /*@
1352: DMMeshSymmetrize - Creates support (out-edge) information from cone (in-edge) inoformation
1354: Not collective
1356: Input Parameter:
1357: . mesh - The DMMesh
1359: Output Parameter:
1361: Note:
1362: This should be called after all calls to DMMeshSetCone()
1364: Level: beginner
1366: .seealso: DMMeshCreate(), DMMeshSetChart(), DMMeshSetConeSize(), DMMeshSetCone()
1367: @*/
1368: PetscErrorCode DMMeshSymmetrize(DM dm)
1369: {
1370: DM_Mesh *mesh = (DM_Mesh *) dm->data;
1375: if (mesh->useNewImpl) {
1376: PetscInt *offsets;
1377: PetscInt supportSize;
1378: PetscInt pStart, pEnd, p;
1380: /* Calculate support sizes */
1381: DMMeshGetChart(dm, &pStart, &pEnd);
1382: PetscSectionSetChart(mesh->supportSection, pStart, pEnd);
1383: for(p = pStart; p < pEnd; ++p) {
1384: PetscInt dof, off, c;
1386: PetscSectionGetDof(mesh->coneSection, p, &dof);
1387: PetscSectionGetOffset(mesh->coneSection, p, &off);
1388: for(c = off; c < off+dof; ++c) {
1389: PetscSectionAddDof(mesh->supportSection, mesh->cones[c], 1);
1390: }
1391: }
1392: for(p = pStart; p < pEnd; ++p) {
1393: PetscInt dof;
1395: PetscSectionGetDof(mesh->supportSection, p, &dof);
1396: mesh->maxSupportSize = PetscMax(mesh->maxSupportSize, dof);
1397: }
1398: PetscSectionSetUp(mesh->supportSection);
1399: /* Calculate supports */
1400: PetscSectionGetStorageSize(mesh->supportSection, &supportSize);
1401: PetscMalloc(supportSize * sizeof(PetscInt), &mesh->supports);
1402: PetscMalloc((pEnd - pStart) * sizeof(PetscInt), &offsets);
1403: PetscMemzero(offsets, (pEnd - pStart) * sizeof(PetscInt));
1404: for(p = pStart; p < pEnd; ++p) {
1405: PetscInt dof, off, c;
1407: PetscSectionGetDof(mesh->coneSection, p, &dof);
1408: PetscSectionGetOffset(mesh->coneSection, p, &off);
1409: for(c = off; c < off+dof; ++c) {
1410: const PetscInt q = mesh->cones[c];
1411: PetscInt offS;
1413: PetscSectionGetOffset(mesh->supportSection, q, &offS);
1414: mesh->supports[offS+offsets[q]] = p;
1415: ++offsets[q];
1416: }
1417: }
1418: PetscFree(offsets);
1419: } else {
1420: Obj<PETSC_MESH_TYPE> m;
1421: DMMeshGetMesh(dm, m);
1422: m->getSieve()->symmetrize();
1423: }
1424: return(0);
1425: }
1429: /*@
1430: DMMeshStratify - The Sieve DAG for most topologies is a graded poset (http://en.wikipedia.org/wiki/Graded_poset), and
1431: can be illustrated by Hasse Diagram (a http://en.wikipedia.org/wiki/Hasse_diagram). The strata group all points of the
1432: same grade, and this function calculates the strata. This grade can be seen as the height (or depth) of the point in
1433: the DAG.
1435: Not collective
1437: Input Parameter:
1438: . mesh - The DMMesh
1440: Output Parameter:
1442: Notes:
1443: The normal association for the point grade is element dimension (or co-dimension). For instance, all vertices would
1444: have depth 0, and all edges depth 1. Likewise, all cells heights would have height 0, and all faces height 1.
1446: This should be called after all calls to DMMeshSymmetrize()
1448: Level: beginner
1450: .seealso: DMMeshCreate(), DMMeshSymmetrize()
1451: @*/
1452: PetscErrorCode DMMeshStratify(DM dm)
1453: {
1454: DM_Mesh *mesh = (DM_Mesh *) dm->data;
1459: if (mesh->useNewImpl) {
1460: PetscInt pStart, pEnd, p;
1461: PetscInt numRoots = 0, numLeaves = 0;
1463: /* Calculate depth */
1464: PetscSectionGetChart(mesh->coneSection, &pStart, &pEnd);
1465: /* Initialize roots and count leaves */
1466: for(p = pStart; p < pEnd; ++p) {
1467: PetscInt coneSize, supportSize;
1469: PetscSectionGetDof(mesh->coneSection, p, &coneSize);
1470: PetscSectionGetDof(mesh->supportSection, p, &supportSize);
1471: if (!coneSize && supportSize) {
1472: ++numRoots;
1473: DMMeshSetLabelValue(dm, "depth", p, 0);
1474: } else if (!supportSize && coneSize) {
1475: ++numLeaves;
1476: }
1477: }
1478: if (numRoots + numLeaves == (pEnd - pStart)) {
1479: for(p = pStart; p < pEnd; ++p) {
1480: PetscInt coneSize, supportSize;
1482: PetscSectionGetDof(mesh->coneSection, p, &coneSize);
1483: PetscSectionGetDof(mesh->supportSection, p, &supportSize);
1484: if (!supportSize && coneSize) {
1485: DMMeshSetLabelValue(dm, "depth", p, 1);
1486: }
1487: }
1488: } else {
1489: SETERRQ(((PetscObject) dm)->comm, PETSC_ERR_SUP, "Have not yet coded general stratification");
1490: }
1491: } else {
1492: Obj<PETSC_MESH_TYPE> m;
1493: DMMeshGetMesh(dm, m);
1494: m->stratify();
1495: }
1496: return(0);
1497: }
1501: /*@C
1502: DMMeshGetLabelValue - Get the value in a Sieve Label for the given point, with 0 as the default
1504: Not Collective
1506: Input Parameters:
1507: + dm - The DMMesh object
1508: . name - The label name
1509: - point - The mesh point
1511: Output Parameter:
1512: . value - The label value for this point, or 0 if the point is not in the label
1514: Level: beginner
1516: .keywords: mesh
1517: .seealso: DMMeshSetLabelValue(), DMMeshGetLabelStratum()
1518: @*/
1519: PetscErrorCode DMMeshGetLabelValue(DM dm, const char name[], PetscInt point, PetscInt *value)
1520: {
1521: DM_Mesh *mesh = (DM_Mesh *) dm->data;
1527: if (mesh->useNewImpl) {
1528: SieveLabel next = mesh->labels;
1529: PetscBool flg = PETSC_FALSE;
1530: PetscInt v, p;
1532: *value = 0;
1533: while(next) {
1534: PetscStrcmp(name, next->name, &flg);
1535: if (flg) break;
1536: next = next->next;
1537: }
1538: if (!flg) {SETERRQ1(((PetscObject) dm)->comm, PETSC_ERR_ARG_WRONG, "No label named %s was found", name);}
1539: /* Find, or add, label value */
1540: for(v = 0; v < next->numStrata; ++v) {
1541: for(p = next->stratumOffsets[v]; p < next->stratumOffsets[v]+next->stratumSizes[v]; ++p) {
1542: if (next->points[p] == point) {
1543: *value = next->stratumValues[v];
1544: break;
1545: }
1546: }
1547: }
1548: } else {
1549: ALE::Obj<PETSC_MESH_TYPE> m;
1550: DMMeshGetMesh(dm, m);
1551: *value = m->getValue(m->getLabel(name), point);
1552: }
1553: return(0);
1554: }
1558: /*@C
1559: DMMeshSetLabelValue - Add a point to a Sieve Label with given value
1561: Not Collective
1563: Input Parameters:
1564: + dm - The DMMesh object
1565: . name - The label name
1566: . point - The mesh point
1567: - value - The label value for this point
1569: Output Parameter:
1571: Level: beginner
1573: .keywords: mesh
1574: .seealso: DMMeshGetLabelStratum()
1575: @*/
1576: PetscErrorCode DMMeshSetLabelValue(DM dm, const char name[], PetscInt point, PetscInt value)
1577: {
1578: DM_Mesh *mesh = (DM_Mesh *) dm->data;
1584: if (mesh->useNewImpl) {
1585: SieveLabel next = mesh->labels;
1586: PetscBool flg = PETSC_FALSE;
1587: PetscInt v, p;
1589: /* Find, or create, label */
1590: while(next) {
1591: PetscStrcmp(name, next->name, &flg);
1592: if (flg) break;
1593: next = next->next;
1594: }
1595: if (!flg) {
1596: SieveLabel tmpLabel = mesh->labels;
1597: PetscNew(struct Sieve_Label, &mesh->labels);
1598: mesh->labels->next = tmpLabel;
1599: next = mesh->labels;
1600: PetscStrallocpy(name, &next->name);
1601: }
1602: /* Find, or add, label value */
1603: for(v = 0; v < next->numStrata; ++v) {
1604: if (next->stratumValues[v] == value) break;
1605: }
1606: if (v >= next->numStrata) {
1607: PetscInt *tmpV, *tmpO, *tmpS;
1608: PetscMalloc3(next->numStrata+1,PetscInt,&tmpV,next->numStrata+2,PetscInt,&tmpO,next->numStrata+1,PetscInt,&tmpS);
1609: for(v = 0; v < next->numStrata; ++v) {
1610: tmpV[v] = next->stratumValues[v];
1611: tmpO[v] = next->stratumOffsets[v];
1612: tmpS[v] = next->stratumSizes[v];
1613: }
1614: tmpV[v] = value;
1615: tmpO[v] = v == 0 ? 0 : next->stratumOffsets[v];
1616: tmpS[v] = 0;
1617: tmpO[v+1] = tmpO[v];
1618: ++next->numStrata;
1619: PetscFree3(next->stratumValues,next->stratumOffsets,next->stratumSizes);
1620: next->stratumValues = tmpV;
1621: next->stratumOffsets = tmpO;
1622: next->stratumSizes = tmpS;
1623: }
1624: /* Check whether point exists */
1625: for(p = next->stratumOffsets[v]; p < next->stratumOffsets[v]+next->stratumSizes[v]; ++p) {
1626: if (next->points[p] == point) {
1627: break;
1628: }
1629: }
1630: /* Add point: NEED TO OPTIMIZE */
1631: if (p >= next->stratumOffsets[v]+next->stratumSizes[v]) {
1632: /* Check for reallocation */
1633: if (next->stratumSizes[v] >= next->stratumOffsets[v+1]-next->stratumOffsets[v]) {
1634: PetscInt oldSize = next->stratumOffsets[v+1]-next->stratumOffsets[v];
1635: PetscInt newSize = PetscMax(10, 2*oldSize); /* Double the size, since 2 is the optimal base for this online algorithm */
1636: PetscInt shift = newSize - oldSize;
1637: PetscInt allocSize = next->stratumOffsets[next->numStrata] + shift;
1638: PetscInt *newPoints;
1639: PetscInt w, q;
1641: PetscMalloc(allocSize * sizeof(PetscInt), &newPoints);
1642: for(q = 0; q < next->stratumOffsets[v]+next->stratumSizes[v]; ++q) {
1643: newPoints[q] = next->points[q];
1644: }
1645: for(w = v+1; w < next->numStrata; ++w) {
1646: for(q = next->stratumOffsets[w]; q < next->stratumOffsets[w]+next->stratumSizes[w]; ++q) {
1647: newPoints[q+shift] = next->points[q];
1648: }
1649: next->stratumOffsets[w] += shift;
1650: }
1651: next->stratumOffsets[next->numStrata] += shift;
1652: PetscFree(next->points);
1653: next->points = newPoints;
1654: }
1655: /* Insert point and resort */
1656: next->points[next->stratumOffsets[v]+next->stratumSizes[v]] = point;
1657: ++next->stratumSizes[v];
1658: PetscSortInt(next->stratumSizes[v], &next->points[next->stratumOffsets[v]]);
1659: }
1660: } else {
1661: ALE::Obj<PETSC_MESH_TYPE> m;
1662: DMMeshGetMesh(dm, m);
1663: m->setValue(m->getLabel(name), point, value);
1664: }
1665: return(0);
1666: }
1670: /*@C
1671: DMMeshGetLabelSize - Get the number of different integer ids in a Label
1673: Not Collective
1675: Input Parameters:
1676: + dm - The DMMesh object
1677: - name - The label name
1679: Output Parameter:
1680: . size - The label size (number of different integer ids)
1682: Level: beginner
1684: .keywords: mesh
1685: .seealso: DMMeshSetLabelValue()
1686: @*/
1687: PetscErrorCode DMMeshGetLabelSize(DM dm, const char name[], PetscInt *size)
1688: {
1689: DM_Mesh *mesh = (DM_Mesh *) dm->data;
1696: if (mesh->useNewImpl) {
1697: SieveLabel next = mesh->labels;
1698: PetscBool flg;
1700: *size = 0;
1701: while(next) {
1702: PetscStrcmp(name, next->name, &flg);
1703: if (flg) {
1704: *size = next->numStrata;
1705: break;
1706: }
1707: next = next->next;
1708: }
1709: } else {
1710: ALE::Obj<PETSC_MESH_TYPE> m;
1711: DMMeshGetMesh(dm, m);
1712: *size = m->getLabel(name)->getCapSize();
1713: }
1714: return(0);
1715: }
1719: /*@C
1720: DMMeshGetLabelIdIS - Get the integer ids in a label
1722: Not Collective
1724: Input Parameters:
1725: + mesh - The DMMesh object
1726: - name - The label name
1728: Output Parameter:
1729: . ids - The integer ids
1731: Level: beginner
1733: .keywords: mesh
1734: .seealso: DMMeshGetLabelSize()
1735: @*/
1736: PetscErrorCode DMMeshGetLabelIdIS(DM dm, const char name[], IS *ids)
1737: {
1738: DM_Mesh *mesh = (DM_Mesh *) dm->data;
1739: PetscInt *values;
1740: PetscInt size, i = 0;
1747: if (mesh->useNewImpl) {
1748: SieveLabel next = mesh->labels;
1749: PetscBool flg;
1751: while(next) {
1752: PetscStrcmp(name, next->name, &flg);
1753: if (flg) {
1754: size = next->numStrata;
1755: PetscMalloc(size * sizeof(PetscInt), &values);
1756: for(i = 0; i < next->numStrata; ++i) {
1757: values[i] = next->stratumValues[i];
1758: }
1759: break;
1760: }
1761: next = next->next;
1762: }
1763: } else {
1764: ALE::Obj<PETSC_MESH_TYPE> m;
1765: DMMeshGetMesh(dm, m);
1766: const ALE::Obj<PETSC_MESH_TYPE::label_type::capSequence>& labelIds = m->getLabel(name)->cap();
1767: const PETSC_MESH_TYPE::label_type::capSequence::const_iterator iEnd = labelIds->end();
1769: size = labelIds->size();
1770: PetscMalloc(size * sizeof(PetscInt), &values);
1771: for(PETSC_MESH_TYPE::label_type::capSequence::const_iterator i_iter = labelIds->begin(); i_iter != iEnd; ++i_iter, ++i) {
1772: values[i] = *i_iter;
1773: }
1774: }
1775: ISCreateGeneral(((PetscObject) dm)->comm, size, values, PETSC_OWN_POINTER, ids);
1776: return(0);
1777: }
1781: /*@C
1782: DMMeshGetStratumSize - Get the number of points in a label stratum
1784: Not Collective
1786: Input Parameters:
1787: + dm - The DMMesh object
1788: . name - The label name
1789: - value - The stratum value
1791: Output Parameter:
1792: . size - The stratum size
1794: Level: beginner
1796: .keywords: mesh
1797: .seealso: DMMeshGetLabelSize(), DMMeshGetLabelIds()
1798: @*/
1799: PetscErrorCode DMMeshGetStratumSize(DM dm, const char name[], PetscInt value, PetscInt *size)
1800: {
1801: DM_Mesh *mesh = (DM_Mesh *) dm->data;
1808: if (mesh->useNewImpl) {
1809: SieveLabel next = mesh->labels;
1810: PetscBool flg;
1812: *size = 0;
1813: while(next) {
1814: PetscStrcmp(name, next->name, &flg);
1815: if (flg) {
1816: PetscInt v;
1818: for(v = 0; v < next->numStrata; ++v) {
1819: if (next->stratumValues[v] == value) {
1820: *size = next->stratumSizes[v];
1821: break;
1822: }
1823: }
1824: break;
1825: }
1826: next = next->next;
1827: }
1828: } else {
1829: ALE::Obj<PETSC_MESH_TYPE> m;
1830: DMMeshGetMesh(dm, m);
1831: *size = m->getLabelStratum(name, value)->size();
1832: }
1833: return(0);
1834: }
1838: /*@C
1839: DMMeshGetStratumIS - Get the points in a label stratum
1841: Not Collective
1843: Input Parameters:
1844: + dm - The DMMesh object
1845: . name - The label name
1846: - value - The stratum value
1848: Output Parameter:
1849: . is - The stratum points
1851: Level: beginner
1853: .keywords: mesh
1854: .seealso: DMMeshGetStratumSize()
1855: @*/
1856: PetscErrorCode DMMeshGetStratumIS(DM dm, const char name[], PetscInt value, IS *is) {
1857: DM_Mesh *mesh = (DM_Mesh *) dm->data;
1864: *is = PETSC_NULL;
1865: if (mesh->useNewImpl) {
1866: SieveLabel next = mesh->labels;
1867: PetscBool flg;
1869: while(next) {
1870: PetscStrcmp(name, next->name, &flg);
1871: if (flg) {
1872: PetscInt v;
1874: for(v = 0; v < next->numStrata; ++v) {
1875: if (next->stratumValues[v] == value) {
1876: ISCreateGeneral(PETSC_COMM_SELF, next->stratumSizes[v], &next->points[next->stratumOffsets[v]], PETSC_COPY_VALUES, is);
1877: break;
1878: }
1879: }
1880: break;
1881: }
1882: next = next->next;
1883: }
1884: } else {
1885: ALE::Obj<PETSC_MESH_TYPE> mesh;
1886: DMMeshGetMesh(dm, mesh);
1887: if (mesh->hasLabel(name)) {
1888: const Obj<PETSC_MESH_TYPE::label_sequence>& stratum = mesh->getLabelStratum(name, value);
1889: PetscInt *idx, i = 0;
1891: PetscMalloc(stratum->size() * sizeof(PetscInt), &idx);
1892: for(PETSC_MESH_TYPE::label_sequence::iterator e_iter = stratum->begin(); e_iter != stratum->end(); ++e_iter, ++i) {
1893: idx[i] = *e_iter;
1894: }
1895: ISCreateGeneral(PETSC_COMM_SELF, stratum->size(), idx, PETSC_OWN_POINTER, is);
1896: }
1897: }
1898: return(0);
1899: }
1903: /* This is a 1-level join */
1904: PetscErrorCode DMMeshJoinPoints(DM dm, const PetscInt points[], PetscInt *coveredPoint)
1905: {
1906: DM_Mesh *mesh = (DM_Mesh *) dm->data;
1913: if (mesh->useNewImpl) {
1914: SETERRQ(((PetscObject) dm)->comm, PETSC_ERR_SUP, "Not yet supported");
1915: } else {
1916: ALE::Obj<PETSC_MESH_TYPE> m;
1917: DMMeshGetMesh(dm, m);
1918: /* const Obj<typename Mesh::sieve_type::supportSet> edge = m->getSieve()->nJoin(points[0], points[1], 1); */
1919: *coveredPoint = -1;
1920: }
1921: return(0);
1922: }
1926: /* This is a 1-level meet */
1927: PetscErrorCode DMMeshMeetPoints(DM dm, PetscInt numPoints, const PetscInt points[], PetscInt *numCoveringPoints, const PetscInt **coveringPoints)
1928: {
1929: DM_Mesh *mesh = (DM_Mesh *) dm->data;
1930: PetscInt *meet[2];
1931: PetscInt meetSize, i = 0;
1932: PetscInt dof, off, p, c, m;
1940: if (mesh->useNewImpl) {
1941: if (!mesh->meetTmpA) {PetscMalloc2(mesh->maxConeSize,PetscInt,&mesh->meetTmpA,mesh->maxConeSize,PetscInt,&mesh->meetTmpB);}
1942: meet[0] = mesh->meetTmpA; meet[1] = mesh->meetTmpB;
1943: /* Copy in cone of first point */
1944: PetscSectionGetDof(mesh->coneSection, points[0], &dof);
1945: PetscSectionGetOffset(mesh->coneSection, points[0], &off);
1946: for(meetSize = 0; meetSize < dof; ++meetSize) {
1947: meet[i][meetSize] = mesh->cones[off+meetSize];
1948: }
1949: /* Check each successive cone */
1950: for(p = 1; p < numPoints; ++p) {
1951: PetscInt newMeetSize = 0;
1953: PetscSectionGetDof(mesh->coneSection, points[p], &dof);
1954: PetscSectionGetOffset(mesh->coneSection, points[p], &off);
1955: for(c = 0; c < dof; ++c) {
1956: const PetscInt point = mesh->cones[off+c];
1958: for(m = 0; m < meetSize; ++m) {
1959: if (point == meet[i][m]) {
1960: meet[1-i][newMeetSize++] = point;
1961: break;
1962: }
1963: }
1964: }
1965: meetSize = newMeetSize;
1966: i = 1-i;
1967: }
1968: *numCoveringPoints = meetSize;
1969: *coveringPoints = meet[1-i];
1970: } else {
1971: ALE::Obj<PETSC_MESH_TYPE> m;
1972: DMMeshGetMesh(dm, m);
1973: /* const Obj<typename Mesh::sieve_type::supportSet> edge = m->getSieve()->nJoin(points[0], points[1], 1); */
1974: *numCoveringPoints = 0;
1975: *coveringPoints = PETSC_NULL;
1976: }
1977: return(0);
1978: }
1982: /*@C
1983: DMMeshGetMaximumDegree - Return the maximum degree of any mesh vertex
1985: Collective on mesh
1987: Input Parameter:
1988: . mesh - The DMMesh
1990: Output Parameter:
1991: . maxDegree - The maximum number of edges at any vertex
1993: Level: beginner
1995: .seealso: DMMeshCreate()
1996: @*/
1997: PetscErrorCode DMMeshGetMaximumDegree(DM dm, PetscInt *maxDegree)
1998: {
1999: Obj<PETSC_MESH_TYPE> m;
2003: DMMeshGetMesh(dm, m);
2004: const ALE::Obj<PETSC_MESH_TYPE::label_sequence>& vertices = m->depthStratum(0);
2005: const ALE::Obj<PETSC_MESH_TYPE::sieve_type>& sieve = m->getSieve();
2006: PetscInt maxDeg = -1;
2008: for(PETSC_MESH_TYPE::label_sequence::iterator v_iter = vertices->begin(); v_iter != vertices->end(); ++v_iter) {
2009: maxDeg = PetscMax(maxDeg, (PetscInt) sieve->getSupportSize(*v_iter));
2010: }
2011: *maxDegree = maxDeg;
2012: return(0);
2013: }
2015: extern PetscErrorCode assembleFullField(VecScatter, Vec, Vec, InsertMode);
2019: /*@
2020: DMMeshRestrictVector - Insert values from a global vector into a local ghosted vector
2022: Collective on g
2024: Input Parameters:
2025: + g - The global vector
2026: . l - The local vector
2027: - mode - either ADD_VALUES or INSERT_VALUES, where
2028: ADD_VALUES adds values to any existing entries, and
2029: INSERT_VALUES replaces existing entries with new values
2031: Level: beginner
2033: .seealso: MatSetOption()
2034: @*/
2035: PetscErrorCode DMMeshRestrictVector(Vec g, Vec l, InsertMode mode)
2036: {
2037: VecScatter injection;
2041: PetscLogEventBegin(DMMesh_restrictVector,0,0,0,0);
2042: PetscObjectQuery((PetscObject) g, "injection", (PetscObject *) &injection);
2043: if (injection) {
2044: VecScatterBegin(injection, g, l, mode, SCATTER_REVERSE);
2045: VecScatterEnd(injection, g, l, mode, SCATTER_REVERSE);
2046: } else {
2047: if (mode == INSERT_VALUES) {
2048: VecCopy(g, l);
2049: } else {
2050: VecAXPY(l, 1.0, g);
2051: }
2052: }
2053: PetscLogEventEnd(DMMesh_restrictVector,0,0,0,0);
2054: return(0);
2055: }
2059: /*@
2060: DMMeshAssembleVectorComplete - Insert values from a local ghosted vector into a global vector
2062: Collective on g
2064: Input Parameters:
2065: + g - The global vector
2066: . l - The local vector
2067: - mode - either ADD_VALUES or INSERT_VALUES, where
2068: ADD_VALUES adds values to any existing entries, and
2069: INSERT_VALUES replaces existing entries with new values
2071: Level: beginner
2073: .seealso: MatSetOption()
2074: @*/
2075: PetscErrorCode DMMeshAssembleVectorComplete(Vec g, Vec l, InsertMode mode)
2076: {
2077: VecScatter injection;
2081: PetscLogEventBegin(DMMesh_assembleVectorComplete,0,0,0,0);
2082: PetscObjectQuery((PetscObject) g, "injection", (PetscObject *) &injection);
2083: if (injection) {
2084: VecScatterBegin(injection, l, g, mode, SCATTER_FORWARD);
2085: VecScatterEnd(injection, l, g, mode, SCATTER_FORWARD);
2086: } else {
2087: if (mode == INSERT_VALUES) {
2088: VecCopy(l, g);
2089: } else {
2090: VecAXPY(g, 1.0, l);
2091: }
2092: }
2093: PetscLogEventEnd(DMMesh_assembleVectorComplete,0,0,0,0);
2094: return(0);
2095: }
2099: /*@
2100: DMMeshAssembleVector - Insert values into a vector
2102: Collective on A
2104: Input Parameters:
2105: + b - the vector
2106: . e - The element number
2107: . v - The values
2108: - mode - either ADD_VALUES or INSERT_VALUES, where
2109: ADD_VALUES adds values to any existing entries, and
2110: INSERT_VALUES replaces existing entries with new values
2112: Level: beginner
2114: .seealso: VecSetOption()
2115: @*/
2116: PetscErrorCode DMMeshAssembleVector(Vec b, PetscInt e, PetscScalar v[], InsertMode mode)
2117: {
2118: DM dm;
2119: SectionReal section;
2123: PetscObjectQuery((PetscObject) b, "DM", (PetscObject *) &dm);
2124: DMMeshGetSectionReal(dm, "x", §ion);
2125: DMMeshAssembleVector(b, dm, section, e, v, mode);
2126: SectionRealDestroy(§ion);
2127: return(0);
2128: }
2130: PetscErrorCode DMMeshAssembleVector(Vec b, DM dm, SectionReal section, PetscInt e, PetscScalar v[], InsertMode mode)
2131: {
2132: ALE::Obj<PETSC_MESH_TYPE> m;
2133: ALE::Obj<PETSC_MESH_TYPE::real_section_type> s;
2134: PetscInt firstElement;
2135: PetscErrorCode ierr;
2138: PetscLogEventBegin(DMMesh_assembleVector,0,0,0,0);
2139: DMMeshGetMesh(dm, m);
2140: SectionRealGetSection(section, s);
2141: //firstElement = elementBundle->getLocalSizes()[bundle->getCommRank()];
2142: firstElement = 0;
2143: #ifdef PETSC_USE_COMPLEX
2144: SETERRQ(((PetscObject)mesh)->comm,PETSC_ERR_SUP, "SectionReal does not support complex update");
2145: #else
2146: if (mode == INSERT_VALUES) {
2147: m->update(s, PETSC_MESH_TYPE::point_type(e + firstElement), v);
2148: } else {
2149: m->updateAdd(s, PETSC_MESH_TYPE::point_type(e + firstElement), v);
2150: }
2151: #endif
2152: PetscLogEventEnd(DMMesh_assembleVector,0,0,0,0);
2153: return(0);
2154: }
2158: /*@C
2159: MatSetValuesTopology - Sets values in a matrix using DM Mesh points rather than indices
2161: Not Collective
2163: Input Parameters:
2164: + mat - the matrix
2165: . dmr - The row DM
2166: . nrow, rowPoints - number of rows and their local Sieve points
2167: . dmc - The column DM
2168: . ncol, colPoints - number of columns and their local Sieve points
2169: . v - a logically two-dimensional array of values
2170: - mode - either ADD_VALUES or INSERT_VALUES, where
2171: ADD_VALUES adds values to any existing entries, and
2172: INSERT_VALUES replaces existing entries with new values
2174: Level: intermediate
2176: .seealso: DMMeshCreate(), MatSetValuesStencil()
2177: @*/
2178: PetscErrorCode MatSetValuesTopology(Mat mat, DM dmr, PetscInt nrow, const PetscInt rowPoints[], DM dmc, PetscInt ncol, const PetscInt colPoints[], const PetscScalar v[], InsertMode mode)
2179: {
2180: ALE::Obj<PETSC_MESH_TYPE> mr;
2181: ALE::Obj<PETSC_MESH_TYPE> mc;
2187: if (!nrow || !ncol) return(0); /* no values to insert */
2193: DMMeshGetMesh(dmr, mr);
2194: DMMeshGetMesh(dmc, mc);
2195: typedef ALE::ISieveVisitor::IndicesVisitor<PETSC_MESH_TYPE::real_section_type,PETSC_MESH_TYPE::order_type,PetscInt> visitor_type;
2196: visitor_type rV(*mr->getRealSection("default"), *mr->getFactory()->getLocalOrder(mr, "default", mr->getRealSection("default")),
2197: (int) pow((double) mr->getSieve()->getMaxConeSize(), mr->depth())*mr->getMaxDof()*nrow, mr->depth() > 1);
2198: visitor_type cV(*mc->getRealSection("default"), *mc->getFactory()->getLocalOrder(mc, "default", mc->getRealSection("default")),
2199: (int) pow((double) mc->getSieve()->getMaxConeSize(), mc->depth())*mc->getMaxDof()*ncol, mc->depth() > 1);
2201: try {
2202: for(PetscInt r = 0; r < nrow; ++r) {
2203: ALE::ISieveTraversal<PETSC_MESH_TYPE::sieve_type>::orientedClosure(*mr->getSieve(), rowPoints[r], rV);
2204: }
2205: } catch(ALE::Exception e) {
2206: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB, e.message());
2207: }
2208: const PetscInt *rowIndices = rV.getValues();
2209: const int numRowIndices = rV.getSize();
2210: try {
2211: for(PetscInt c = 0; c < ncol; ++c) {
2212: ALE::ISieveTraversal<PETSC_MESH_TYPE::sieve_type>::orientedClosure(*mc->getSieve(), colPoints[c], cV);
2213: }
2214: } catch(ALE::Exception e) {
2215: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB, e.message());
2216: }
2217: const PetscInt *colIndices = cV.getValues();
2218: const int numColIndices = cV.getSize();
2220: MatSetValuesLocal(mat, numRowIndices, rowIndices, numColIndices, colIndices, v, mode);
2221: return(0);
2222: }
2226: PetscErrorCode DMMeshUpdateOperator(Mat A, const ALE::Obj<PETSC_MESH_TYPE>& m, const ALE::Obj<PETSC_MESH_TYPE::real_section_type>& section, const ALE::Obj<PETSC_MESH_TYPE::order_type>& globalOrder, const PETSC_MESH_TYPE::point_type& e, PetscScalar array[], InsertMode mode)
2227: {
2231: typedef ALE::ISieveVisitor::IndicesVisitor<PETSC_MESH_TYPE::real_section_type,PETSC_MESH_TYPE::order_type,PetscInt> visitor_type;
2232: visitor_type iV(*section, *globalOrder, (int) pow((double) m->getSieve()->getMaxConeSize(), m->depth())*m->getMaxDof(), m->depth() > 1);
2234: updateOperator(A, *m->getSieve(), iV, e, array, mode);
2235: return(0);
2236: }
2240: PetscErrorCode DMMeshUpdateOperatorGeneral(Mat A, const ALE::Obj<PETSC_MESH_TYPE>& rowM, const ALE::Obj<PETSC_MESH_TYPE::real_section_type>& rowSection, const ALE::Obj<PETSC_MESH_TYPE::order_type>& rowGlobalOrder, const PETSC_MESH_TYPE::point_type& rowE, const ALE::Obj<PETSC_MESH_TYPE>& colM, const ALE::Obj<PETSC_MESH_TYPE::real_section_type>& colSection, const ALE::Obj<PETSC_MESH_TYPE::order_type>& colGlobalOrder, const PETSC_MESH_TYPE::point_type& colE, PetscScalar array[], InsertMode mode)
2241: {
2242: typedef ALE::ISieveVisitor::IndicesVisitor<PETSC_MESH_TYPE::real_section_type,PETSC_MESH_TYPE::order_type,PetscInt> visitor_type;
2243: visitor_type iVr(*rowSection, *rowGlobalOrder, (int) pow((double) rowM->getSieve()->getMaxConeSize(), rowM->depth())*rowM->getMaxDof(), rowM->depth() > 1);
2244: visitor_type iVc(*colSection, *colGlobalOrder, (int) pow((double) colM->getSieve()->getMaxConeSize(), colM->depth())*colM->getMaxDof(), colM->depth() > 1);
2246: PetscErrorCode updateOperator(A, *rowM->getSieve(), iVr, rowE, *colM->getSieve(), iVc, colE, array, mode);
2247: return(0);
2248: }
2252: /*@
2253: DMMeshSetMaxDof - Sets the maximum number of degrees of freedom on any sieve point
2255: Logically Collective on A
2257: Input Parameters:
2258: + A - the matrix
2259: . mesh - DMMesh needed for orderings
2260: . section - A Section which describes the layout
2261: . e - The element number
2262: . v - The values
2263: - mode - either ADD_VALUES or INSERT_VALUES, where
2264: ADD_VALUES adds values to any existing entries, and
2265: INSERT_VALUES replaces existing entries with new values
2267: Notes: This is used by routines like DMMeshUpdateOperator() to bound buffer sizes
2269: Level: developer
2271: .seealso: DMMeshUpdateOperator(), DMMeshAssembleMatrix()
2272: @*/
2273: PetscErrorCode DMMeshSetMaxDof(DM dm, PetscInt maxDof)
2274: {
2275: Obj<PETSC_MESH_TYPE> m;
2276: PetscErrorCode ierr;
2279: DMMeshGetMesh(dm, m);
2280: m->setMaxDof(maxDof);
2281: return(0);
2282: }
2286: /*@
2287: DMMeshAssembleMatrix - Insert values into a matrix
2289: Collective on A
2291: Input Parameters:
2292: + A - the matrix
2293: . dm - DMMesh needed for orderings
2294: . section - A Section which describes the layout
2295: . e - The element
2296: . v - The values
2297: - mode - either ADD_VALUES or INSERT_VALUES, where
2298: ADD_VALUES adds values to any existing entries, and
2299: INSERT_VALUES replaces existing entries with new values
2301: Level: beginner
2303: .seealso: MatSetOption()
2304: @*/
2305: PetscErrorCode DMMeshAssembleMatrix(Mat A, DM dm, SectionReal section, PetscInt e, PetscScalar v[], InsertMode mode)
2306: {
2310: PetscLogEventBegin(DMMesh_assembleMatrix,0,0,0,0);
2311: try {
2312: Obj<PETSC_MESH_TYPE> m;
2313: Obj<PETSC_MESH_TYPE::real_section_type> s;
2315: DMMeshGetMesh(dm, m);
2316: SectionRealGetSection(section, s);
2317: const ALE::Obj<PETSC_MESH_TYPE::order_type>& globalOrder = m->getFactory()->getGlobalOrder(m, s->getName(), s);
2319: if (m->debug()) {
2320: std::cout << "Assembling matrix for element number " << e << " --> point " << e << std::endl;
2321: }
2322: DMMeshUpdateOperator(A, m, s, globalOrder, e, v, mode);
2323: } catch (ALE::Exception e) {
2324: std::cout << e.msg() << std::endl;
2325: }
2326: PetscLogEventEnd(DMMesh_assembleMatrix,0,0,0,0);
2327: return(0);
2328: }
2330: /******************************** C Wrappers **********************************/
2334: PetscErrorCode WriteVTKHeader(DM dm, PetscViewer viewer)
2335: {
2336: ALE::Obj<PETSC_MESH_TYPE> m;
2339: DMMeshGetMesh(dm, m);
2340: return VTKViewer::writeHeader(m, viewer);
2341: }
2345: PetscErrorCode WriteVTKVertices(DM dm, PetscViewer viewer)
2346: {
2347: ALE::Obj<PETSC_MESH_TYPE> m;
2350: DMMeshGetMesh(dm, m);
2351: return VTKViewer::writeVertices(m, viewer);
2352: }
2356: PetscErrorCode WriteVTKElements(DM dm, PetscViewer viewer)
2357: {
2358: ALE::Obj<PETSC_MESH_TYPE> m;
2361: DMMeshGetMesh(dm, m);
2362: return VTKViewer::writeElements(m, viewer);
2363: }
2367: /*@C
2368: DMMeshGetCoordinates - Creates an array holding the coordinates.
2370: Not Collective
2372: Input Parameter:
2373: + dm - The DMMesh object
2374: - columnMajor - Flag for column major order
2376: Output Parameter:
2377: + numVertices - The number of vertices
2378: . dim - The embedding dimension
2379: - coords - The array holding local coordinates
2381: Level: intermediate
2383: .keywords: mesh, coordinates
2384: .seealso: DMMeshCreate()
2385: @*/
2386: PetscErrorCode DMMeshGetCoordinates(DM dm, PetscBool columnMajor, PetscInt *numVertices, PetscInt *dim, PetscReal *coords[])
2387: {
2388: ALE::Obj<PETSC_MESH_TYPE> m;
2389: PetscErrorCode ierr;
2392: DMMeshGetMesh(dm, m);
2393: ALE::PCICE::Builder::outputVerticesLocal(m, numVertices, dim, coords, columnMajor);
2394: return(0);
2395: }
2399: /*@C
2400: DMMeshGetElements - Creates an array holding the vertices on each element.
2402: Not Collective
2404: Input Parameters:
2405: + dm - The DMMesh object
2406: - columnMajor - Flag for column major order
2408: Output Parameters:
2409: + numElements - The number of elements
2410: . numCorners - The number of vertices per element
2411: - vertices - The array holding vertices on each local element
2413: Level: intermediate
2415: .keywords: mesh, elements
2416: .seealso: DMMeshCreate()
2417: @*/
2418: PetscErrorCode DMMeshGetElements(DM dm, PetscBool columnMajor, PetscInt *numElements, PetscInt *numCorners, PetscInt *vertices[])
2419: {
2420: ALE::Obj<PETSC_MESH_TYPE> m;
2421: PetscErrorCode ierr;
2424: DMMeshGetMesh(dm, m);
2425: ALE::PCICE::Builder::outputElementsLocal(m, numElements, numCorners, vertices, columnMajor);
2426: return(0);
2427: }
2431: PetscErrorCode DMMeshCreateNeighborCSR(DM dm, PetscInt *numVertices, PetscInt **offsets, PetscInt **adjacency) {
2432: const PetscInt maxFaceCases = 30;
2433: PetscInt numFaceCases = 0;
2434: PetscInt numFaceVertices[maxFaceCases];
2435: PetscInt *off, *adj;
2436: PetscInt dim, depth, cStart, cEnd, c, numCells, cell;
2440: /* For parallel partitioning, I think you have to communicate supports */
2441: DMMeshGetDimension(dm, &dim);
2442: DMMeshGetLabelSize(dm, "depth", &depth);
2443: --depth;
2444: DMMeshGetHeightStratum(dm, 0, &cStart, &cEnd);
2445: if (cEnd - cStart == 0) {
2446: *numVertices = 0;
2447: *offsets = PETSC_NULL;
2448: *adjacency = PETSC_NULL;
2449: return(0);
2450: }
2451: numCells = cEnd - cStart;
2452: /* Setup face recognition */
2453: {
2454: PetscInt cornersSeen[30] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0}; /* Could use PetscBT */
2456: for(c = cStart; c < cEnd; ++c) {
2457: PetscInt corners;
2459: DMMeshGetConeSize(dm, c, &corners);
2460: if (!cornersSeen[corners]) {
2461: if (numFaceCases >= maxFaceCases) {SETERRQ(((PetscObject) dm)->comm, PETSC_ERR_PLIB, "Exceeded maximum number of face recognition cases");}
2462: cornersSeen[corners] = 1;
2463: if (corners == dim+1) {
2464: numFaceVertices[numFaceCases] = dim;
2465: PetscInfo(dm, "Recognizing simplices\n");
2466: } else if ((dim == 1) && (corners == 3)) {
2467: numFaceVertices[numFaceCases] = 3;
2468: PetscInfo(dm, "Recognizing quadratic edges\n");
2469: } else if ((dim == 2) && (corners == 4)) {
2470: numFaceVertices[numFaceCases] = 2;
2471: PetscInfo(dm, "Recognizing quads\n");
2472: } else if ((dim == 2) && (corners == 6)) {
2473: numFaceVertices[numFaceCases] = 3;
2474: PetscInfo(dm, "Recognizing tri and quad cohesive Lagrange cells\n");
2475: } else if ((dim == 2) && (corners == 9)) {
2476: numFaceVertices[numFaceCases] = 3;
2477: PetscInfo(dm, "Recognizing quadratic quads and quadratic quad cohesive Lagrange cells\n");
2478: } else if ((dim == 3) && (corners == 6)) {
2479: numFaceVertices[numFaceCases] = 4;
2480: PetscInfo(dm, "Recognizing tet cohesive cells\n");
2481: } else if ((dim == 3) && (corners == 8)) {
2482: numFaceVertices[numFaceCases] = 4;
2483: PetscInfo(dm, "Recognizing hexes\n");
2484: } else if ((dim == 3) && (corners == 9)) {
2485: numFaceVertices[numFaceCases] = 6;
2486: PetscInfo(dm, "Recognizing tet cohesive Lagrange cells\n");
2487: } else if ((dim == 3) && (corners == 10)) {
2488: numFaceVertices[numFaceCases] = 6;
2489: PetscInfo(dm, "Recognizing quadratic tets\n");
2490: } else if ((dim == 3) && (corners == 12)) {
2491: numFaceVertices[numFaceCases] = 6;
2492: PetscInfo(dm, "Recognizing hex cohesive Lagrange cells\n");
2493: } else if ((dim == 3) && (corners == 18)) {
2494: numFaceVertices[numFaceCases] = 6;
2495: PetscInfo(dm, "Recognizing quadratic tet cohesive Lagrange cells\n");
2496: } else if ((dim == 3) && (corners == 27)) {
2497: numFaceVertices[numFaceCases] = 9;
2498: PetscInfo(dm, "Recognizing quadratic hexes and quadratic hex cohesive Lagrange cells\n");
2499: } else {
2500: SETERRQ1(((PetscObject) dm)->comm, PETSC_ERR_ARG_OUTOFRANGE, "Could not recognize number of face vertices for %d corners", corners);
2501: }
2502: ++numFaceCases;
2503: }
2504: }
2505: }
2506: /* Check for optimized depth 1 construction */
2507: PetscMalloc((numCells+1) * sizeof(PetscInt), &off);
2508: PetscMemzero(off, (numCells+1) * sizeof(PetscInt));
2509: if (depth == 1) {
2510: PetscInt *neighborCells;
2511: PetscInt n;
2512: PetscInt maxConeSize, maxSupportSize;
2514: /* Temp space for point adj <= maxConeSize*maxSupportSize */
2515: DMMeshGetMaxSizes(dm, &maxConeSize, &maxSupportSize);
2516: PetscMalloc(maxConeSize*maxSupportSize * sizeof(PetscInt), &neighborCells);
2517: /* Count neighboring cells */
2518: for(cell = cStart; cell < cEnd; ++cell) {
2519: const PetscInt *cone;
2520: PetscInt numNeighbors = 0;
2521: PetscInt coneSize, c;
2523: /* Get support of the cone, and make a set of the cells */
2524: DMMeshGetConeSize(dm, cell, &coneSize);
2525: DMMeshGetCone(dm, cell, &cone);
2526: for(c = 0; c < coneSize; ++c) {
2527: const PetscInt *support;
2528: PetscInt supportSize, s;
2530: DMMeshGetSupportSize(dm, cone[c], &supportSize);
2531: DMMeshGetSupport(dm, cone[c], &support);
2532: for(s = 0; s < supportSize; ++s) {
2533: const PetscInt point = support[s];
2535: if (point == cell) continue;
2536: for(n = 0; n < numNeighbors; ++n) {
2537: if (neighborCells[n] == point) break;
2538: }
2539: if (n == numNeighbors) {
2540: neighborCells[n] = point;
2541: ++numNeighbors;
2542: }
2543: }
2544: }
2545: /* Get meet with each cell, and check with recognizer (could optimize to check each pair only once) */
2546: for(n = 0; n < numNeighbors; ++n) {
2547: PetscInt cellPair[2] = {cell, neighborCells[n]};
2548: PetscInt meetSize;
2549: const PetscInt *meet;
2551: DMMeshMeetPoints(dm, 2, cellPair, &meetSize, &meet);
2552: if (meetSize) {
2553: PetscInt f;
2555: for(f = 0; f < numFaceCases; ++f) {
2556: if (numFaceVertices[f] == meetSize) {
2557: ++off[cell-cStart+1];
2558: break;
2559: }
2560: }
2561: }
2562: }
2563: }
2564: /* Prefix sum */
2565: for(cell = 1; cell <= numCells; ++cell) {
2566: off[cell] += off[cell-1];
2567: }
2568: PetscMalloc(off[numCells] * sizeof(PetscInt), &adj);
2569: /* Get neighboring cells */
2570: for(cell = cStart; cell < cEnd; ++cell) {
2571: const PetscInt *cone;
2572: PetscInt numNeighbors = 0;
2573: PetscInt cellOffset = 0;
2574: PetscInt coneSize, c;
2576: /* Get support of the cone, and make a set of the cells */
2577: DMMeshGetConeSize(dm, cell, &coneSize);
2578: DMMeshGetCone(dm, cell, &cone);
2579: for(c = 0; c < coneSize; ++c) {
2580: const PetscInt *support;
2581: PetscInt supportSize, s;
2583: DMMeshGetSupportSize(dm, cone[c], &supportSize);
2584: DMMeshGetSupport(dm, cone[c], &support);
2585: for(s = 0; s < supportSize; ++s) {
2586: const PetscInt point = support[s];
2588: if (point == cell) continue;
2589: for(n = 0; n < numNeighbors; ++n) {
2590: if (neighborCells[n] == point) break;
2591: }
2592: if (n == numNeighbors) {
2593: neighborCells[n] = point;
2594: ++numNeighbors;
2595: }
2596: }
2597: }
2598: /* Get meet with each cell, and check with recognizer (could optimize to check each pair only once) */
2599: for(n = 0; n < numNeighbors; ++n) {
2600: PetscInt cellPair[2] = {cell, neighborCells[n]};
2601: PetscInt meetSize;
2602: const PetscInt *meet;
2604: DMMeshMeetPoints(dm, 2, cellPair, &meetSize, &meet);
2605: if (meetSize) {
2606: PetscInt f;
2608: for(f = 0; f < numFaceCases; ++f) {
2609: if (numFaceVertices[f] == meetSize) {
2610: adj[off[cell-cStart]+cellOffset] = neighborCells[n];
2611: ++cellOffset;
2612: break;
2613: }
2614: }
2615: }
2616: }
2617: }
2618: PetscFree(neighborCells);
2619: } else if (depth == dim) {
2620: SETERRQ(((PetscObject) dm)->comm, PETSC_ERR_ARG_WRONG, "Neighbor graph creation not implemented for interpolated meshes");
2621: #if 0
2622: OffsetVisitor<typename Mesh::sieve_type> oV(*sieve, *overlapSieve, off);
2623: PetscInt p;
2625: for(typename Mesh::label_sequence::iterator c_iter = cells->begin(); c_iter != cEnd; ++c_iter) {
2626: sieve->cone(*c_iter, oV);
2627: }
2628: for(p = 1; p <= numCells; ++p) {
2629: off[p] = off[p] + off[p-1];
2630: }
2631: PetscMalloc(off[numCells] * sizeof(PetscInt), &adj);
2632: AdjVisitor<typename Mesh::sieve_type> aV(adj, zeroBase);
2633: ISieveVisitor::SupportVisitor<typename Mesh::sieve_type, AdjVisitor<typename Mesh::sieve_type> > sV(*sieve, aV);
2634: ISieveVisitor::SupportVisitor<typename Mesh::sieve_type, AdjVisitor<typename Mesh::sieve_type> > ovSV(*overlapSieve, aV);
2636: for(typename Mesh::label_sequence::iterator c_iter = cells->begin(); c_iter != cEnd; ++c_iter) {
2637: aV.setCell(*c_iter);
2638: sieve->cone(*c_iter, sV);
2639: sieve->cone(*c_iter, ovSV);
2640: }
2641: offset = aV.getOffset();
2642: #endif
2643: } else {
2644: SETERRQ(((PetscObject) dm)->comm, PETSC_ERR_ARG_WRONG, "Neighbor graph creation not defined for partially interpolated meshes");
2645: }
2646: *numVertices = numCells;
2647: *offsets = off;
2648: *adjacency = adj;
2649: return(0);
2650: }
2652: #ifdef PETSC_HAVE_CHACO
2653: #ifdef PETSC_HAVE_UNISTD_H
2654: #include <unistd.h>
2655: #endif
2656: /* Chaco does not have an include file */
2657: extern "C" {
2658: extern int interface(int nvtxs, int *start, int *adjacency, int *vwgts,
2659: float *ewgts, float *x, float *y, float *z, char *outassignname,
2660: char *outfilename, short *assignment, int architecture, int ndims_tot,
2661: int mesh_dims[3], double *goal, int global_method, int local_method,
2662: int rqi_flag, int vmax, int ndims, double eigtol, long seed);
2664: extern int FREE_GRAPH;
2665: }
2669: PetscErrorCode DMMeshPartition_Chaco(DM dm, PetscInt numVertices, PetscInt start[], PetscInt adjacency[], PetscSection *partSection, IS *partition)
2670: {
2671: enum {DEFAULT_METHOD = 1, INERTIAL_METHOD = 3};
2672: MPI_Comm comm = ((PetscObject) dm)->comm;
2673: int nvtxs = numVertices; /* number of vertices in full graph */
2674: int *vwgts = NULL; /* weights for all vertices */
2675: float *ewgts = NULL; /* weights for all edges */
2676: float *x = NULL, *y = NULL, *z = NULL; /* coordinates for inertial method */
2677: char *outassignname = NULL; /* name of assignment output file */
2678: char *outfilename = NULL; /* output file name */
2679: int architecture = 1; /* 0 => hypercube, d => d-dimensional mesh */
2680: int ndims_tot = 0; /* total number of cube dimensions to divide */
2681: int mesh_dims[3]; /* dimensions of mesh of processors */
2682: double *goal = NULL; /* desired set sizes for each set */
2683: int global_method = 1; /* global partitioning algorithm */
2684: int local_method = 1; /* local partitioning algorithm */
2685: int rqi_flag = 0; /* should I use RQI/Symmlq eigensolver? */
2686: int vmax = 200; /* how many vertices to coarsen down to? */
2687: int ndims = 1; /* number of eigenvectors (2^d sets) */
2688: double eigtol = 0.001; /* tolerance on eigenvectors */
2689: long seed = 123636512; /* for random graph mutations */
2690: short int *assignment; /* Output partition */
2691: int fd_stdout, fd_pipe[2];
2692: PetscInt *points;
2693: PetscMPIInt commSize;
2694: int i, v, p;
2698: MPI_Comm_size(comm, &commSize);
2699: if (!numVertices) {
2700: PetscSectionCreate(comm, partSection);
2701: PetscSectionSetChart(*partSection, 0, commSize);
2702: PetscSectionSetUp(*partSection);
2703: ISCreateGeneral(comm, 0, PETSC_NULL, PETSC_OWN_POINTER, partition);
2704: return(0);
2705: }
2706: FREE_GRAPH = 0; /* Do not let Chaco free my memory */
2707: for(i = 0; i < start[numVertices]; ++i) {
2708: ++adjacency[i];
2709: }
2710: if (global_method == INERTIAL_METHOD) {
2711: /* manager.createCellCoordinates(nvtxs, &x, &y, &z); */
2712: SETERRQ(comm, PETSC_ERR_SUP, "Inertial partitioning not yet supported");
2713: }
2714: mesh_dims[0] = commSize;
2715: mesh_dims[1] = 1;
2716: mesh_dims[2] = 1;
2717: PetscMalloc(nvtxs * sizeof(short int), &assignment);
2718: /* Chaco outputs to stdout. We redirect this to a buffer. */
2719: /* TODO: check error codes for UNIX calls */
2720: #ifdef PETSC_HAVE_UNISTD_H
2721: {
2722: fd_stdout = dup(1);
2723: pipe(fd_pipe);
2724: close(1);
2725: dup2(fd_pipe[1], 1);
2726: }
2727: #endif
2728: interface(nvtxs, (int *) start, (int *) adjacency, vwgts, ewgts, x, y, z, outassignname, outfilename,
2729: assignment, architecture, ndims_tot, mesh_dims, goal, global_method, local_method, rqi_flag,
2730: vmax, ndims, eigtol, seed);
2731: #ifdef PETSC_HAVE_UNISTD_H
2732: {
2733: char msgLog[10000];
2734: int count;
2736: fflush(stdout);
2737: count = read(fd_pipe[0], msgLog, (10000-1)*sizeof(char));
2738: if (count < 0) count = 0;
2739: msgLog[count] = 0;
2740: close(1);
2741: dup2(fd_stdout, 1);
2742: close(fd_stdout);
2743: close(fd_pipe[0]);
2744: close(fd_pipe[1]);
2745: if (ierr) {SETERRQ1(comm, PETSC_ERR_LIB, "Error in Chaco library: %s", msgLog);}
2746: }
2747: #endif
2748: /* Convert to PetscSection+IS */
2749: PetscSectionCreate(comm, partSection);
2750: PetscSectionSetChart(*partSection, 0, commSize);
2751: for(v = 0; v < nvtxs; ++v) {
2752: PetscSectionAddDof(*partSection, assignment[v], 1);
2753: }
2754: PetscSectionSetUp(*partSection);
2755: PetscMalloc(nvtxs * sizeof(PetscInt), &points);
2756: for(p = 0, i = 0; p < commSize; ++p) {
2757: for(v = 0; v < nvtxs; ++v) {
2758: if (assignment[v] == p) points[i++] = v;
2759: }
2760: }
2761: if (i != nvtxs) {SETERRQ2(comm, PETSC_ERR_PLIB, "Number of points %d should be %d", i, nvtxs);}
2762: ISCreateGeneral(comm, nvtxs, points, PETSC_OWN_POINTER, partition);
2763: if (global_method == INERTIAL_METHOD) {
2764: /* manager.destroyCellCoordinates(nvtxs, &x, &y, &z); */
2765: }
2766: PetscFree(assignment);
2767: for(i = 0; i < start[numVertices]; ++i) {
2768: --adjacency[i];
2769: }
2770: return(0);
2771: }
2772: #endif
2774: #ifdef PETSC_HAVE_PARMETIS
2777: PetscErrorCode DMMeshPartition_ParMetis(DM dm, PetscInt numVertices, PetscInt start[], PetscInt adjacency[], PetscSection *partSection, IS *partition)
2778: {
2780: return(0);
2781: }
2782: #endif
2786: PetscErrorCode DMMeshCreatePartition(DM dm, PetscSection *partSection, IS *partition, PetscInt height) {
2787: PetscMPIInt size;
2791: MPI_Comm_size(((PetscObject) dm)->comm, &size);
2792: if (size == 1) {
2793: PetscInt *points;
2794: PetscInt cStart, cEnd, c;
2796: DMMeshGetHeightStratum(dm, 0, &cStart, &cEnd);
2797: PetscSectionCreate(((PetscObject) dm)->comm, partSection);
2798: PetscSectionSetChart(*partSection, 0, size);
2799: PetscSectionSetDof(*partSection, 0, cEnd-cStart);
2800: PetscSectionSetUp(*partSection);
2801: PetscMalloc((cEnd - cStart) * sizeof(PetscInt), &points);
2802: for(c = cStart; c < cEnd; ++c) {
2803: points[c] = c;
2804: }
2805: ISCreateGeneral(((PetscObject) dm)->comm, cEnd-cStart, points, PETSC_OWN_POINTER, partition);
2806: return(0);
2807: }
2808: if (height == 0) {
2809: PetscInt numVertices;
2810: PetscInt *start = PETSC_NULL;
2811: PetscInt *adjacency = PETSC_NULL;
2813: if (1) {
2814: DMMeshCreateNeighborCSR(dm, &numVertices, &start, &adjacency);
2815: #ifdef PETSC_HAVE_CHACO
2816: DMMeshPartition_Chaco(dm, numVertices, start, adjacency, partSection, partition);
2817: #endif
2818: } else {
2819: DMMeshCreateNeighborCSR(dm, &numVertices, &start, &adjacency);
2820: #ifdef PETSC_HAVE_PARMETIS
2821: DMMeshPartition_ParMetis(dm, numVertices, start, adjacency, partSection, partition);
2822: #endif
2823: }
2824: PetscFree(start);
2825: PetscFree(adjacency);
2826: # if 0
2827: } else if (height == 1) {
2828: /* Build the dual graph for faces and partition the hypergraph */
2829: PetscInt numEdges;
2831: buildFaceCSRV(mesh, mesh->getFactory()->getNumbering(mesh, mesh->depth()-1), &numEdges, &start, &adjacency, GraphPartitioner::zeroBase());
2832: GraphPartitioner().partition(numEdges, start, adjacency, partition, manager);
2833: destroyCSR(numEdges, start, adjacency);
2834: #endif
2835: } else {
2836: SETERRQ1(((PetscObject) dm)->comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid partition height %d", height);
2837: }
2838: return(0);
2839: }
2843: PetscErrorCode DMMeshCreatePartitionClosure(DM dm, PetscSection pointSection, IS pointPartition, PetscSection *section, IS *partition) {
2844: const PetscInt *partArray;
2845: PetscInt *allPoints, *partPoints = PETSC_NULL;
2846: PetscInt rStart, rEnd, rank, maxPartSize = 0, newSize;
2847: PetscErrorCode ierr;
2850: PetscSectionGetChart(pointSection, &rStart, &rEnd);
2851: ISGetIndices(pointPartition, &partArray);
2852: PetscSectionCreate(((PetscObject) dm)->comm, section);
2853: PetscSectionSetChart(*section, rStart, rEnd);
2854: for(rank = rStart; rank < rEnd; ++rank) {
2855: PetscInt partSize = 0;
2856: PetscInt numPoints, offset, p;
2858: PetscSectionGetDof(pointSection, rank, &numPoints);
2859: PetscSectionGetOffset(pointSection, rank, &offset);
2860: for(p = 0; p < numPoints; ++p) {
2861: PetscInt point = partArray[offset+p], closureSize;
2862: const PetscInt *closure;
2864: /* TODO Include support for height > 0 case */
2865: DMMeshGetTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
2866: /* Merge into existing points */
2867: if (partSize+closureSize > maxPartSize) {
2868: PetscInt *tmpPoints;
2870: maxPartSize = PetscMax(partSize+closureSize, 2*maxPartSize);
2871: PetscMalloc(maxPartSize * sizeof(PetscInt), &tmpPoints);
2872: PetscMemcpy(tmpPoints, partPoints, partSize * sizeof(PetscInt));
2873: PetscFree(partPoints);
2874: partPoints = tmpPoints;
2875: }
2876: PetscMemcpy(&partPoints[partSize], closure, closureSize * sizeof(PetscInt));
2877: partSize += closureSize;
2878: PetscSortRemoveDupsInt(&partSize, partPoints);
2879: }
2880: PetscSectionSetDof(*section, rank, partSize);
2881: }
2882: PetscSectionSetUp(*section);
2883: PetscSectionGetStorageSize(*section, &newSize);
2884: PetscMalloc(newSize * sizeof(PetscInt), &allPoints);
2886: for(rank = rStart; rank < rEnd; ++rank) {
2887: PetscInt partSize = 0, newOffset;
2888: PetscInt numPoints, offset, p;
2890: PetscSectionGetDof(pointSection, rank, &numPoints);
2891: PetscSectionGetOffset(pointSection, rank, &offset);
2892: for(p = 0; p < numPoints; ++p) {
2893: PetscInt point = partArray[offset+p], closureSize;
2894: const PetscInt *closure;
2896: /* TODO Include support for height > 0 case */
2897: DMMeshGetTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
2898: /* Merge into existing points */
2899: PetscMemcpy(&partPoints[partSize], closure, closureSize * sizeof(PetscInt));
2900: partSize += closureSize;
2901: PetscSortRemoveDupsInt(&partSize, partPoints);
2902: }
2903: PetscSectionGetOffset(*section, rank, &newOffset);
2904: PetscMemcpy(&allPoints[newOffset], partPoints, partSize * sizeof(PetscInt));
2905: }
2906: ISRestoreIndices(pointPartition, &partArray);
2907: PetscFree(partPoints);
2908: ISCreateGeneral(((PetscObject) dm)->comm, newSize, allPoints, PETSC_OWN_POINTER, partition);
2909: return(0);
2910: }
2914: /*
2915: Input Parameters:
2916: . originalSection
2917: , originalVec
2919: Output Parameters:
2920: . newSection
2921: . newVec
2922: */
2923: PetscErrorCode DMMeshDistributeField(DM dm, PetscSF pointSF, PetscSection originalSection, Vec originalVec, PetscSection newSection, Vec newVec)
2924: {
2925: PetscSF fieldSF;
2926: PetscInt *remoteOffsets, fieldSize;
2927: PetscScalar *originalValues, *newValues;
2928: PetscErrorCode ierr;
2931: PetscSFDistributeSection(pointSF, originalSection, &remoteOffsets, newSection);
2933: PetscSectionGetStorageSize(newSection, &fieldSize);
2934: VecSetSizes(newVec, fieldSize, PETSC_DETERMINE);
2935: VecSetFromOptions(newVec);
2937: VecGetArray(originalVec, &originalValues);
2938: VecGetArray(newVec, &newValues);
2939: PetscSFCreateSectionSF(pointSF, originalSection, remoteOffsets, newSection, &fieldSF);
2940: PetscSFBcastBegin(fieldSF, MPIU_SCALAR, originalValues, newValues);
2941: PetscSFBcastEnd(fieldSF, MPIU_SCALAR, originalValues, newValues);
2942: PetscSFDestroy(&fieldSF);
2943: VecRestoreArray(newVec, &newValues);
2944: VecRestoreArray(originalVec, &originalValues);
2945: return(0);
2946: }
2950: /*@C
2951: DMMeshDistribute - Distributes the mesh and any associated sections.
2953: Not Collective
2955: Input Parameter:
2956: + dm - The original DMMesh object
2957: - partitioner - The partitioning package, or NULL for the default
2959: Output Parameter:
2960: . parallelMesh - The distributed DMMesh object
2962: Level: intermediate
2964: .keywords: mesh, elements
2966: .seealso: DMMeshCreate(), DMMeshDistributeByFace()
2967: @*/
2968: PetscErrorCode DMMeshDistribute(DM dm, const char partitioner[], DM *dmParallel)
2969: {
2970: DM_Mesh *mesh = (DM_Mesh *) dm->data, *pmesh;
2971: MPI_Comm comm = ((PetscObject) dm)->comm;
2972: PetscMPIInt rank, numProcs, p;
2978: MPI_Comm_rank(comm, &rank);
2979: MPI_Comm_size(comm, &numProcs);
2980: if (numProcs == 1) return(0);
2981: if (mesh->useNewImpl) {
2982: const PetscInt height = 0;
2983: PetscInt dim, numRemoteRanks;
2984: IS cellPart, part;
2985: PetscSection cellPartSection, partSection;
2986: PetscSFNode *remoteRanks;
2987: PetscSF partSF, pointSF, coneSF;
2988: ISLocalToGlobalMapping renumbering;
2989: PetscSection originalConeSection, newConeSection;
2990: PetscInt *remoteOffsets;
2991: PetscInt *cones, *newCones, newConesSize;
2992: PetscBool flg;
2994: DMMeshGetDimension(dm, &dim);
2995: /* Create cell partition - We need to rewrite to use IS, use the MatPartition stuff */
2996: DMMeshCreatePartition(dm, &cellPartSection, &cellPart, height);
2997: /* Create SF assuming a serial partition for all processes: Could check for IS length here */
2998: if (!rank) {
2999: numRemoteRanks = numProcs;
3000: } else {
3001: numRemoteRanks = 0;
3002: }
3003: PetscMalloc(numRemoteRanks * sizeof(PetscSFNode), &remoteRanks);
3004: for(p = 0; p < numRemoteRanks; ++p) {
3005: remoteRanks[p].rank = p;
3006: remoteRanks[p].index = 0;
3007: }
3008: PetscSFCreate(comm, &partSF);
3009: PetscSFSetGraph(partSF, 1, numRemoteRanks, PETSC_NULL, PETSC_OWN_POINTER, remoteRanks, PETSC_OWN_POINTER);
3010: PetscOptionsHasName(((PetscObject) dm)->prefix, "-partition_view", &flg);
3011: if (flg) {
3012: PetscPrintf(comm, "Cell Partition:\n");
3013: PetscSectionView(cellPartSection, PETSC_VIEWER_STDOUT_WORLD);
3014: ISView(cellPart, PETSC_NULL);
3015: PetscSFView(partSF, PETSC_NULL);
3016: }
3017: /* Close the partition over the mesh */
3018: DMMeshCreatePartitionClosure(dm, cellPartSection, cellPart, &partSection, &part);
3019: ISDestroy(&cellPart);
3020: PetscSectionDestroy(&cellPartSection);
3021: /* Create new mesh */
3022: DMMeshCreate(comm, dmParallel);
3023: DMMeshSetDimension(*dmParallel, dim);
3024: PetscObjectSetName((PetscObject) *dmParallel, "Parallel Mesh");
3025: pmesh = (DM_Mesh *) (*dmParallel)->data;
3026: /* Distribute sieve points and the global point numbering (replaces creating remote bases) */
3027: PetscSFConvertPartition(partSF, partSection, part, &renumbering, &pointSF);
3028: if (flg) {
3029: PetscPrintf(comm, "Point Partition:\n");
3030: PetscSectionView(partSection, PETSC_VIEWER_STDOUT_WORLD);
3031: ISView(part, PETSC_NULL);
3032: PetscSFView(pointSF, PETSC_NULL);
3033: PetscPrintf(comm, "Point Renumbering after partition:\n");
3034: ISLocalToGlobalMappingView(renumbering, PETSC_NULL);
3035: }
3036: /* Distribute cone section */
3037: DMMeshGetConeSection(dm, &originalConeSection);
3038: DMMeshGetConeSection(*dmParallel, &newConeSection);
3039: PetscSFDistributeSection(pointSF, originalConeSection, &remoteOffsets, newConeSection);
3040: DMMeshSetUp(*dmParallel);
3041: {
3042: PetscInt pStart, pEnd, p;
3044: PetscSectionGetChart(newConeSection, &pStart, &pEnd);
3045: for(p = pStart; p < pEnd; ++p) {
3046: PetscInt coneSize;
3047: PetscSectionGetDof(newConeSection, p, &coneSize);
3048: pmesh->maxConeSize = PetscMax(pmesh->maxConeSize, coneSize);
3049: }
3050: }
3051: /* Communicate and renumber cones */
3052: PetscSFCreateSectionSF(pointSF, originalConeSection, remoteOffsets, newConeSection, &coneSF);
3053: DMMeshGetCones(dm, &cones);
3054: DMMeshGetCones(*dmParallel, &newCones);
3055: PetscSFBcastBegin(coneSF, MPIU_INT, cones, newCones);
3056: PetscSFBcastEnd(coneSF, MPIU_INT, cones, newCones);
3057: PetscSectionGetStorageSize(newConeSection, &newConesSize);
3058: ISGlobalToLocalMappingApply(renumbering, IS_GTOLM_MASK, newConesSize, newCones, PETSC_NULL, newCones);
3059: PetscOptionsHasName(((PetscObject) dm)->prefix, "-cones_view", &flg);
3060: if (flg) {
3061: PetscPrintf(comm, "Serial Cone Section:\n");
3062: PetscSectionView(originalConeSection, PETSC_VIEWER_STDOUT_WORLD);
3063: PetscPrintf(comm, "Parallel Cone Section:\n");
3064: PetscSectionView(newConeSection, PETSC_VIEWER_STDOUT_WORLD);
3065: PetscSFView(coneSF, PETSC_NULL);
3066: }
3067: PetscSFDestroy(&coneSF);
3068: /* Create supports and stratify sieve */
3069: DMMeshSymmetrize(*dmParallel);
3070: DMMeshStratify(*dmParallel);
3071: /* Distribute Coordinates */
3072: {
3073: PetscSection originalCoordSection, newCoordSection;
3074: Vec originalCoordinates, newCoordinates;
3076: DMMeshGetCoordinateSection(dm, &originalCoordSection);
3077: DMMeshGetCoordinateSection(*dmParallel, &newCoordSection);
3078: DMMeshGetCoordinateVec(dm, &originalCoordinates);
3079: DMMeshGetCoordinateVec(*dmParallel, &newCoordinates);
3081: DMMeshDistributeField(dm, pointSF, originalCoordSection, originalCoordinates, newCoordSection, newCoordinates);
3082: }
3083: /* Distribute labels */
3084: {
3085: SieveLabel next = mesh->labels, newNext = PETSC_NULL;
3086: PetscInt numLabels = 0, l;
3088: /* Bcast number of labels */
3089: while(next) {++numLabels; next = next->next;}
3090: MPI_Bcast(&numLabels, 1, MPIU_INT, 0, comm);
3091: next = mesh->labels;
3092: for(l = 0; l < numLabels; ++l) {
3093: SieveLabel newLabel;
3094: const PetscInt *partArray;
3095: PetscInt *stratumSizes = PETSC_NULL, *points = PETSC_NULL;
3096: PetscMPIInt *sendcnts = PETSC_NULL, *offsets = PETSC_NULL, *displs = PETSC_NULL;
3097: PetscInt nameSize, s, p;
3098: size_t len = 0;
3100: PetscNew(struct Sieve_Label, &newLabel);
3101: /* Bcast name (could filter for no points) */
3102: if (!rank) {PetscStrlen(next->name, &len);}
3103: nameSize = len;
3104: MPI_Bcast(&nameSize, 1, MPIU_INT, 0, comm);
3105: PetscMalloc(nameSize+1, &newLabel->name);
3106: if (!rank) {PetscMemcpy(newLabel->name, next->name, nameSize+1);}
3107: MPI_Bcast(newLabel->name, nameSize+1, MPI_CHAR, 0, comm);
3108: /* Bcast numStrata (could filter for no points in stratum) */
3109: if (!rank) {newLabel->numStrata = next->numStrata;}
3110: MPI_Bcast(&newLabel->numStrata, 1, MPIU_INT, 0, comm);
3111: PetscMalloc(newLabel->numStrata * sizeof(PetscInt), &newLabel->stratumValues);
3112: PetscMalloc(newLabel->numStrata * sizeof(PetscInt), &newLabel->stratumSizes);
3113: PetscMalloc((newLabel->numStrata+1) * sizeof(PetscInt), &newLabel->stratumOffsets);
3114: /* Bcast stratumValues (could filter for no points in stratum) */
3115: if (!rank) {PetscMemcpy(newLabel->stratumValues, next->stratumValues, next->numStrata * sizeof(PetscInt));}
3116: MPI_Bcast(newLabel->stratumValues, newLabel->numStrata, MPIU_INT, 0, comm);
3117: /* Find size on each process and Scatter */
3118: if (!rank) {
3119: ISGetIndices(part, &partArray);
3120: PetscMalloc(numProcs*next->numStrata * sizeof(PetscInt), &stratumSizes);
3121: PetscMemzero(stratumSizes, numProcs*next->numStrata * sizeof(PetscInt));
3122: for(s = 0; s < next->numStrata; ++s) {
3123: for(p = next->stratumOffsets[s]; p < next->stratumOffsets[s]+next->stratumSizes[s]; ++p) {
3124: const PetscInt point = next->points[p];
3125: PetscInt proc;
3127: for(proc = 0; proc < numProcs; ++proc) {
3128: PetscInt dof, off, pPart;
3130: PetscSectionGetDof(partSection, proc, &dof);
3131: PetscSectionGetOffset(partSection, proc, &off);
3132: for(pPart = off; pPart < off+dof; ++pPart) {
3133: if (partArray[pPart] == point) {
3134: ++stratumSizes[proc*next->numStrata+s];
3135: break;
3136: }
3137: }
3138: }
3139: }
3140: }
3141: ISRestoreIndices(part, &partArray);
3142: }
3143: MPI_Scatter(stratumSizes, newLabel->numStrata, MPI_INT, newLabel->stratumSizes, newLabel->numStrata, MPI_INT, 0, comm);
3144: /* Calculate stratumOffsets */
3145: newLabel->stratumOffsets[0] = 0;
3146: for(s = 0; s < newLabel->numStrata; ++s) {
3147: newLabel->stratumOffsets[s+1] = newLabel->stratumSizes[s] + newLabel->stratumOffsets[s];
3148: }
3149: /* Pack points and Scatter */
3150: if (!rank) {
3151: PetscMalloc3(numProcs,PetscMPIInt,&sendcnts,numProcs,PetscMPIInt,&offsets,numProcs+1,PetscMPIInt,&displs);
3152: displs[0] = 0;
3153: for(p = 0; p < numProcs; ++p) {
3154: sendcnts[p] = 0;
3155: for(s = 0; s < next->numStrata; ++s) {
3156: sendcnts[p] += stratumSizes[p*next->numStrata+s];
3157: }
3158: offsets[p] = displs[p];
3159: displs[p+1] = displs[p] + sendcnts[p];
3160: }
3161: PetscMalloc(displs[numProcs] * sizeof(PetscInt), &points);
3162: for(s = 0; s < next->numStrata; ++s) {
3163: for(p = next->stratumOffsets[s]; p < next->stratumOffsets[s]+next->stratumSizes[s]; ++p) {
3164: const PetscInt point = next->points[p];
3165: PetscInt proc;
3167: for(proc = 0; proc < numProcs; ++proc) {
3168: PetscInt dof, off, pPart;
3170: PetscSectionGetDof(partSection, proc, &dof);
3171: PetscSectionGetOffset(partSection, proc, &off);
3172: for(pPart = off; pPart < off+dof; ++pPart) {
3173: if (partArray[pPart] == point) {
3174: points[offsets[proc]++] = point;
3175: break;
3176: }
3177: }
3178: }
3179: }
3180: }
3181: }
3182: PetscMalloc(newLabel->stratumOffsets[newLabel->numStrata] * sizeof(PetscInt), &newLabel->points);
3183: MPI_Scatterv(points, sendcnts, displs, MPIU_INT, newLabel->points, newLabel->stratumOffsets[newLabel->numStrata], MPIU_INT, 0, comm);
3184: PetscFree(points);
3185: PetscFree3(sendcnts,offsets,displs);
3186: PetscFree(stratumSizes);
3187: /* Renumber points */
3188: ISGlobalToLocalMappingApply(renumbering, IS_GTOLM_MASK, newLabel->stratumOffsets[newLabel->numStrata], newLabel->points, PETSC_NULL, newLabel->points);
3189: /* Sort points */
3190: for(s = 0; s < newLabel->numStrata; ++s) {
3191: PetscSortInt(newLabel->stratumSizes[s], &newLabel->points[newLabel->stratumOffsets[s]]);
3192: }
3193: /* Insert into list */
3194: if (newNext) {
3195: newNext->next = newLabel;
3196: } else {
3197: pmesh->labels = newLabel;
3198: }
3199: newNext = newLabel;
3200: if (!rank) {next = next->next;}
3201: }
3202: }
3203: /* Cleanup Partition */
3204: ISLocalToGlobalMappingDestroy(&renumbering);
3205: PetscSFDestroy(&partSF);
3206: PetscSectionDestroy(&partSection);
3207: ISDestroy(&part);
3208: /* Create point SF for parallel mesh */
3209: {
3210: const PetscInt *leaves;
3211: PetscSFNode *remotePoints;
3212: PetscMPIInt *rowners, *lowners;
3213: PetscInt *ghostPoints, numRoots, numLeaves, numGhostPoints = 0, p, gp;
3215: PetscSFGetGraph(pointSF, &numRoots, &numLeaves, &leaves, PETSC_NULL);
3216: PetscMalloc2(numRoots*2,PetscInt,&rowners,numLeaves*2,PetscInt,&lowners);
3217: for(p = 0; p < numRoots*2; ++p) {
3218: rowners[p] = 0;
3219: }
3220: for(p = 0; p < numLeaves; ++p) {
3221: lowners[p*2+0] = rank;
3222: lowners[p*2+1] = PetscMPIIntCast(leaves ? leaves[p] : p);
3223: }
3224: PetscSFReduceBegin(pointSF, MPI_2INT, lowners, rowners, MPI_MAXLOC);
3225: PetscSFReduceEnd(pointSF, MPI_2INT, lowners, rowners, MPI_MAXLOC);
3226: PetscSFBcastBegin(pointSF, MPI_2INT, rowners, lowners);
3227: PetscSFBcastEnd(pointSF, MPI_2INT, rowners, lowners);
3228: for(p = 0; p < numLeaves; ++p) {
3229: if (lowners[p*2+0] != rank) ++numGhostPoints;
3230: }
3231: PetscMalloc(numGhostPoints * sizeof(PetscInt), &ghostPoints);
3232: PetscMalloc(numGhostPoints * sizeof(PetscSFNode), &remotePoints);
3233: for(p = 0, gp = 0; p < numLeaves; ++p) {
3234: if (lowners[p*2+0] != rank) {
3235: ghostPoints[gp] = leaves ? leaves[p] : p;
3236: remotePoints[gp].rank = lowners[p*2+0];
3237: remotePoints[gp].index = lowners[p*2+1];
3238: ++gp;
3239: }
3240: }
3241: PetscFree2(rowners,lowners);
3242: PetscSFCreate(comm, &pmesh->sf);
3243: PetscSFSetGraph(pmesh->sf, PETSC_DETERMINE, numGhostPoints, ghostPoints, PETSC_OWN_POINTER, remotePoints, PETSC_OWN_POINTER);
3244: PetscSFSetFromOptions(pmesh->sf);
3245: }
3246: /* Cleanup */
3247: PetscSFDestroy(&pointSF);
3248: DMSetFromOptions(*dmParallel);
3249: } else {
3250: ALE::Obj<PETSC_MESH_TYPE> oldMesh;
3251: DMMeshGetMesh(dm, oldMesh);
3252: DMMeshCreate(comm, dmParallel);
3253: const Obj<PETSC_MESH_TYPE> newMesh = new PETSC_MESH_TYPE(comm, oldMesh->getDimension(), oldMesh->debug());
3254: const Obj<PETSC_MESH_TYPE::sieve_type> newSieve = new PETSC_MESH_TYPE::sieve_type(comm, oldMesh->debug());
3256: newMesh->setSieve(newSieve);
3257: ALE::DistributionNew<PETSC_MESH_TYPE>::distributeMeshAndSectionsV(oldMesh, newMesh);
3258: DMMeshSetMesh(*dmParallel, newMesh);
3259: }
3260: return(0);
3261: }
3265: /*@C
3266: DMMeshDistribute - Distributes the mesh and any associated sections.
3268: Not Collective
3270: Input Parameter:
3271: + serialMesh - The original DMMesh object
3272: - partitioner - The partitioning package, or NULL for the default
3274: Output Parameter:
3275: . parallelMesh - The distributed DMMesh object
3277: Level: intermediate
3279: .keywords: mesh, elements
3281: .seealso: DMMeshCreate(), DMMeshDistribute()
3282: @*/
3283: PetscErrorCode DMMeshDistributeByFace(DM serialMesh, const char partitioner[], DM *parallelMesh)
3284: {
3285: ALE::Obj<PETSC_MESH_TYPE> oldMesh;
3286: PetscErrorCode ierr;
3289: DMMeshGetMesh(serialMesh, oldMesh);
3290: DMMeshCreate(oldMesh->comm(), parallelMesh);
3291: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP, "I am being lazy, bug me.");
3292: #if 0
3293: ALE::DistributionNew<PETSC_MESH_TYPE>::distributeMeshAndSectionsV(oldMesh, newMesh, height = 1);
3294: #endif
3295: return(0);
3296: }
3298: #ifdef PETSC_HAVE_TRIANGLE
3299: /* Already included since C++ is turned on #include <triangle.h> */
3303: PetscErrorCode TriangleInitInput(struct triangulateio *inputCtx) {
3305: inputCtx->numberofpoints = 0;
3306: inputCtx->numberofpointattributes = 0;
3307: inputCtx->pointlist = PETSC_NULL;
3308: inputCtx->pointattributelist = PETSC_NULL;
3309: inputCtx->pointmarkerlist = PETSC_NULL;
3310: inputCtx->numberofsegments = 0;
3311: inputCtx->segmentlist = PETSC_NULL;
3312: inputCtx->segmentmarkerlist = PETSC_NULL;
3313: inputCtx->numberoftriangleattributes = 0;
3314: inputCtx->trianglelist = PETSC_NULL;
3315: inputCtx->numberofholes = 0;
3316: inputCtx->holelist = PETSC_NULL;
3317: inputCtx->numberofregions = 0;
3318: inputCtx->regionlist = PETSC_NULL;
3319: return(0);
3320: }
3324: PetscErrorCode TriangleInitOutput(struct triangulateio *outputCtx) {
3326: outputCtx->numberofpoints = 0;
3327: outputCtx->pointlist = PETSC_NULL;
3328: outputCtx->pointattributelist = PETSC_NULL;
3329: outputCtx->pointmarkerlist = PETSC_NULL;
3330: outputCtx->numberoftriangles = 0;
3331: outputCtx->trianglelist = PETSC_NULL;
3332: outputCtx->triangleattributelist = PETSC_NULL;
3333: outputCtx->neighborlist = PETSC_NULL;
3334: outputCtx->segmentlist = PETSC_NULL;
3335: outputCtx->segmentmarkerlist = PETSC_NULL;
3336: outputCtx->edgelist = PETSC_NULL;
3337: outputCtx->edgemarkerlist = PETSC_NULL;
3338: return(0);
3339: }
3343: PetscErrorCode TriangleFiniOutput(struct triangulateio *outputCtx) {
3345: free(outputCtx->pointmarkerlist);
3346: free(outputCtx->edgelist);
3347: free(outputCtx->edgemarkerlist);
3348: free(outputCtx->trianglelist);
3349: free(outputCtx->neighborlist);
3350: return(0);
3351: }
3355: PetscErrorCode DMMeshGenerate_Triangle(DM boundary, PetscBool interpolate, DM *dm)
3356: {
3357: MPI_Comm comm = ((PetscObject) boundary)->comm;
3358: DM_Mesh *bd = (DM_Mesh *) boundary->data;
3359: PetscInt dim = 2;
3360: const PetscBool createConvexHull = PETSC_FALSE;
3361: const PetscBool constrained = PETSC_FALSE;
3362: struct triangulateio in;
3363: struct triangulateio out;
3364: PetscInt vStart, vEnd, v, eStart, eEnd, e;
3365: PetscMPIInt rank;
3366: PetscErrorCode ierr;
3369: MPI_Comm_rank(comm, &rank);
3370: TriangleInitInput(&in);
3371: TriangleInitOutput(&out);
3372: DMMeshGetDepthStratum(boundary, 0, &vStart, &vEnd);
3373: in.numberofpoints = vEnd - vStart;
3374: if (in.numberofpoints > 0) {
3375: PetscScalar *array;
3377: PetscMalloc(in.numberofpoints*dim * sizeof(double), &in.pointlist);
3378: PetscMalloc(in.numberofpoints * sizeof(int), &in.pointmarkerlist);
3379: VecGetArray(bd->coordinates, &array);
3380: for(v = vStart; v < vEnd; ++v) {
3381: const PetscInt idx = v - vStart;
3382: PetscInt off, d;
3384: PetscSectionGetOffset(bd->coordSection, v, &off);
3385: for(d = 0; d < dim; ++d) {
3386: in.pointlist[idx*dim + d] = array[off+d];
3387: }
3388: DMMeshGetLabelValue(boundary, "marker", v, &in.pointmarkerlist[idx]);
3389: }
3390: VecRestoreArray(bd->coordinates, &array);
3391: }
3392: DMMeshGetHeightStratum(boundary, 0, &eStart, &eEnd);
3393: in.numberofsegments = eEnd - eStart;
3394: if (in.numberofsegments > 0) {
3395: PetscMalloc(in.numberofsegments*2 * sizeof(int), &in.segmentlist);
3396: PetscMalloc(in.numberofsegments * sizeof(int), &in.segmentmarkerlist);
3397: for(e = eStart; e < eEnd; ++e) {
3398: const PetscInt idx = e - eStart;
3399: const PetscInt *cone;
3401: DMMeshGetCone(boundary, e, &cone);
3402: in.segmentlist[idx*2+0] = cone[0] - vStart;
3403: in.segmentlist[idx*2+1] = cone[1] - vStart;
3404: DMMeshGetLabelValue(boundary, "marker", e, &in.segmentmarkerlist[idx]);
3405: }
3406: }
3407: #if 0 /* Do not currently support holes */
3408: PetscReal *holeCoords;
3409: PetscInt h, d;
3411: DMMeshGetHoles(boundary, &in.numberofholes, &holeCords);
3412: if (in.numberofholes > 0) {
3413: PetscMalloc(in.numberofholes*dim * sizeof(double), &in.holelist);CHKERRXX(ierr);
3414: for(h = 0; h < in.numberofholes; ++h) {
3415: for(d = 0; d < dim; ++d) {
3416: in.holelist[h*dim+d] = holeCoords[h*dim+d];
3417: }
3418: }
3419: }
3420: #endif
3421: if (!rank) {
3422: char args[32];
3424: /* Take away 'Q' for verbose output */
3425: PetscStrcpy(args, "pqezQ");
3426: if (createConvexHull) {
3427: PetscStrcat(args, "c");
3428: }
3429: if (constrained) {
3430: PetscStrcpy(args, "zepDQ");
3431: }
3432: triangulate(args, &in, &out, PETSC_NULL);
3433: }
3434: PetscFree(in.pointlist);
3435: PetscFree(in.pointmarkerlist);
3436: PetscFree(in.segmentlist);
3437: PetscFree(in.segmentmarkerlist);
3438: PetscFree(in.holelist);
3440: DMCreate(comm, dm);
3441: DMSetType(*dm, DMMESH);
3442: DMMeshSetDimension(*dm, dim);
3443: {
3444: DM_Mesh *mesh = (DM_Mesh *) (*dm)->data;
3445: const PetscInt numCorners = 3;
3446: const PetscInt numCells = out.numberoftriangles;
3447: const PetscInt numVertices = out.numberofpoints;
3448: int *cells = out.trianglelist;
3449: double *meshCoords = out.pointlist;
3450: PetscInt coordSize, c;
3451: PetscScalar *coords;
3453: DMMeshSetChart(*dm, 0, numCells+numVertices);
3454: for(c = 0; c < numCells; ++c) {
3455: DMMeshSetConeSize(*dm, c, numCorners);
3456: }
3457: DMMeshSetUp(*dm);
3458: for(c = 0; c < numCells; ++c) {
3459: PetscInt cone[numCorners] = {cells[c*numCorners+0]+numCells, cells[c*numCorners+1]+numCells, cells[c*numCorners+2]+numCells};
3461: DMMeshSetCone(*dm, c, cone);
3462: }
3463: DMMeshSymmetrize(*dm);
3464: DMMeshStratify(*dm);
3465: PetscSectionSetChart(mesh->coordSection, numCells, numCells + numVertices);
3466: for(v = numCells; v < numCells+numVertices; ++v) {
3467: PetscSectionSetDof(mesh->coordSection, v, dim);
3468: }
3469: PetscSectionSetUp(mesh->coordSection);
3470: PetscSectionGetStorageSize(mesh->coordSection, &coordSize);
3471: VecSetSizes(mesh->coordinates, coordSize, PETSC_DETERMINE);
3472: VecSetFromOptions(mesh->coordinates);
3473: VecGetArray(mesh->coordinates, &coords);
3474: for(v = 0; v < numVertices; ++v) {
3475: coords[v*dim+0] = meshCoords[v*dim+0];
3476: coords[v*dim+1] = meshCoords[v*dim+1];
3477: }
3478: VecRestoreArray(mesh->coordinates, &coords);
3479: for(v = 0; v < numVertices; ++v) {
3480: if (out.pointmarkerlist[v]) {
3481: DMMeshSetLabelValue(*dm, "marker", v+numCells, out.pointmarkerlist[v]);
3482: }
3483: }
3484: if (interpolate) {
3485: for(e = 0; e < out.numberofedges; e++) {
3486: if (out.edgemarkerlist[e]) {
3487: const PetscInt vertices[2] = {out.edgelist[e*2+0]+numCells, out.edgelist[e*2+1]+numCells};
3488: PetscInt edge;
3490: DMMeshJoinPoints(*dm, vertices, &edge); /* 1-level join */
3491: DMMeshSetLabelValue(*dm, "marker", edge, out.edgemarkerlist[e]);
3492: }
3493: }
3494: }
3495: }
3496: #if 0 /* Do not currently support holes */
3497: DMMeshCopyHoles(*dm, boundary);
3498: #endif
3499: TriangleFiniOutput(&out);
3500: return(0);
3501: }
3502: #endif
3506: /*@C
3507: DMMeshGenerate - Generates a mesh.
3509: Not Collective
3511: Input Parameters:
3512: + boundary - The DMMesh boundary object
3513: - interpolate - Flag to create intermediate mesh elements
3515: Output Parameter:
3516: . mesh - The DMMesh object
3518: Level: intermediate
3520: .keywords: mesh, elements
3521: .seealso: DMMeshCreate(), DMMeshRefine()
3522: @*/
3523: PetscErrorCode DMMeshGenerate(DM boundary, PetscBool interpolate, DM *mesh)
3524: {
3525: DM_Mesh *bd = (DM_Mesh *) boundary->data;
3531: if (bd->useNewImpl) {
3532: if (interpolate) {SETERRQ(((PetscObject) boundary)->comm, PETSC_ERR_SUP, "Interpolation (creation of faces and edges) is not yet supported.");}
3533: #ifdef PETSC_HAVE_TRIANGLE
3534: DMMeshGenerate_Triangle(boundary, interpolate, mesh);
3535: #endif
3536: } else {
3537: ALE::Obj<PETSC_MESH_TYPE> mB;
3539: DMMeshGetMesh(boundary, mB);
3540: DMMeshCreate(mB->comm(), mesh);
3541: ALE::Obj<PETSC_MESH_TYPE> m = ALE::Generator<PETSC_MESH_TYPE>::generateMeshV(mB, interpolate, false, true);
3542: DMMeshSetMesh(*mesh, m);
3543: }
3544: return(0);
3545: }
3549: /*@C
3550: DMMeshRefine - Refines the mesh.
3552: Not Collective
3554: Input Parameters:
3555: + mesh - The original DMMesh object
3556: . refinementLimit - The maximum size of any cell
3557: - interpolate - Flag to create intermediate mesh elements
3559: Output Parameter:
3560: . refinedMesh - The refined DMMesh object
3562: Level: intermediate
3564: .keywords: mesh, elements
3565: .seealso: DMMeshCreate(), DMMeshGenerate()
3566: @*/
3567: PetscErrorCode DMMeshRefine(DM mesh, double refinementLimit, PetscBool interpolate, DM *refinedMesh)
3568: {
3569: ALE::Obj<PETSC_MESH_TYPE> oldMesh;
3573: if (refinementLimit == 0.0) return(0);
3574: DMMeshGetMesh(mesh, oldMesh);
3575: DMMeshCreate(oldMesh->comm(), refinedMesh);
3576: ALE::Obj<PETSC_MESH_TYPE> newMesh = ALE::Generator<PETSC_MESH_TYPE>::refineMeshV(oldMesh, refinementLimit, interpolate, false, true);
3577: DMMeshSetMesh(*refinedMesh, newMesh);
3578: return(0);
3579: }
3583: PetscErrorCode DMRefine_Mesh(DM dm, MPI_Comm comm, DM *dmRefined)
3584: {
3585: ALE::Obj<PETSC_MESH_TYPE> oldMesh;
3586: double refinementLimit;
3587: PetscErrorCode ierr;
3590: DMMeshGetMesh(dm, oldMesh);
3591: DMMeshCreate(comm, dmRefined);
3592: refinementLimit = oldMesh->getMaxVolume()/2.0;
3593: ALE::Obj<PETSC_MESH_TYPE> newMesh = ALE::Generator<PETSC_MESH_TYPE>::refineMeshV(oldMesh, refinementLimit, true);
3594: DMMeshSetMesh(*dmRefined, newMesh);
3595: return(0);
3596: }
3600: PetscErrorCode DMCoarsenHierarchy_Mesh(DM mesh, int numLevels, DM *coarseHierarchy)
3601: {
3602: PetscReal cfactor = 1.5;
3606: PetscOptionsReal("-mesh_coarsen_factor", "The coarsening factor", PETSC_NULL, cfactor, &cfactor, PETSC_NULL);
3607: SETERRQ(((PetscObject) mesh)->comm, PETSC_ERR_SUP, "Peter needs to incorporate his code.");
3608: return(0);
3609: }
3613: PetscErrorCode DMCreateInterpolation_Mesh(DM dmCoarse, DM dmFine, Mat *interpolation, Vec *scaling) {
3614: SETERRQ(((PetscObject) dmCoarse)->comm, PETSC_ERR_SUP, "Peter needs to incorporate his code.");
3615: }
3619: PetscErrorCode DMMeshMarkBoundaryCells(DM dm, const char labelName[], PetscInt marker, PetscInt newMarker) {
3620: ALE::Obj<PETSC_MESH_TYPE> mesh;
3624: DMMeshGetMesh(dm, mesh);
3625: mesh->markBoundaryCells(labelName, marker, newMarker);
3626: return(0);
3627: }
3631: PetscErrorCode DMMeshGetDepthStratum(DM dm, PetscInt stratumValue, PetscInt *start, PetscInt *end) {
3632: DM_Mesh *mesh = (DM_Mesh *) dm->data;
3637: if (mesh->useNewImpl) {
3638: SieveLabel next = mesh->labels;
3639: PetscBool flg = PETSC_FALSE;
3640: PetscInt depth;
3642: if (stratumValue < 0) {
3643: DMMeshGetChart(dm, start, end);
3644: return(0);
3645: } else {
3646: PetscInt pStart, pEnd;
3648: if (start) {*start = 0;}
3649: if (end) {*end = 0;}
3650: DMMeshGetChart(dm, &pStart, &pEnd);
3651: if (pStart == pEnd) {return(0);}
3652: }
3653: while(next) {
3654: PetscStrcmp("depth", next->name, &flg);
3655: if (flg) break;
3656: next = next->next;
3657: }
3658: if (!flg) {SETERRQ(((PetscObject) dm)->comm, PETSC_ERR_ARG_WRONG, "No label named depth was found");}
3659: /* Strata are sorted and contiguous -- In addition, depth/height is either full or 1-level */
3660: depth = stratumValue;
3661: if ((depth < 0) || (depth >= next->numStrata)) {
3662: if (start) {*start = 0;}
3663: if (end) {*end = 0;}
3664: } else {
3665: if (start) {*start = next->points[next->stratumOffsets[depth]];}
3666: if (end) {*end = next->points[next->stratumOffsets[depth]+next->stratumSizes[depth]-1]+1;}
3667: }
3668: } else {
3669: ALE::Obj<PETSC_MESH_TYPE> mesh;
3670: DMMeshGetMesh(dm, mesh);
3671: if (stratumValue < 0) {
3672: if (start) *start = mesh->getSieve()->getChart().min();
3673: if (end) *end = mesh->getSieve()->getChart().max();
3674: } else {
3675: const Obj<PETSC_MESH_TYPE::label_sequence>& stratum = mesh->depthStratum(stratumValue);
3676: if (start) *start = *stratum->begin();
3677: if (end) *end = *stratum->rbegin()+1;
3678: }
3679: }
3680: return(0);
3681: }
3685: PetscErrorCode DMMeshGetHeightStratum(DM dm, PetscInt stratumValue, PetscInt *start, PetscInt *end) {
3686: DM_Mesh *mesh = (DM_Mesh *) dm->data;
3691: if (mesh->useNewImpl) {
3692: SieveLabel next = mesh->labels;
3693: PetscBool flg = PETSC_FALSE;
3694: PetscInt depth;
3696: if (stratumValue < 0) {
3697: DMMeshGetChart(dm, start, end);
3698: } else {
3699: PetscInt pStart, pEnd;
3701: if (start) {*start = 0;}
3702: if (end) {*end = 0;}
3703: DMMeshGetChart(dm, &pStart, &pEnd);
3704: if (pStart == pEnd) {return(0);}
3705: }
3706: while(next) {
3707: PetscStrcmp("depth", next->name, &flg);
3708: if (flg) break;
3709: next = next->next;
3710: }
3711: if (!flg) {SETERRQ(((PetscObject) dm)->comm, PETSC_ERR_ARG_WRONG, "No label named depth was found");}
3712: /* Strata are sorted and contiguous -- In addition, depth/height is either full or 1-level */
3713: depth = next->stratumValues[next->numStrata-1] - stratumValue;
3714: if ((depth < 0) || (depth >= next->numStrata)) {
3715: if (start) {*start = 0;}
3716: if (end) {*end = 0;}
3717: } else {
3718: if (start) {*start = next->points[next->stratumOffsets[depth]];}
3719: if (end) {*end = next->points[next->stratumOffsets[depth]+next->stratumSizes[depth]-1]+1;}
3720: }
3721: } else {
3722: ALE::Obj<PETSC_MESH_TYPE> mesh;
3723: DMMeshGetMesh(dm, mesh);
3724: if (mesh->getLabel("height")->size()) {
3725: const Obj<PETSC_MESH_TYPE::label_sequence>& stratum = mesh->heightStratum(stratumValue);
3726: if (start) *start = *stratum->begin();
3727: if (end) *end = *stratum->rbegin()+1;
3728: } else {
3729: if (start) *start = 0;
3730: if (end) *end = 0;
3731: }
3732: }
3733: return(0);
3734: }
3738: PetscErrorCode DMMeshCreateSection(DM dm, PetscInt dim, PetscInt numFields, PetscInt numComp[], PetscInt numDof[], PetscInt numBC, PetscInt bcField[], IS bcPoints[], PetscSection *section) {
3739: PetscInt *numDofTot, *maxConstraints;
3740: PetscInt pStart = 0, pEnd = 0;
3744: PetscMalloc2(dim+1,PetscInt,&numDofTot,numFields+1,PetscInt,&maxConstraints);
3745: for(PetscInt d = 0; d <= dim; ++d) {
3746: numDofTot[d] = 0;
3747: for(PetscInt f = 0; f < numFields; ++f) {
3748: numDofTot[d] += numDof[f*(dim+1)+d];
3749: }
3750: }
3751: PetscSectionCreate(((PetscObject) dm)->comm, section);
3752: if (numFields > 1) {
3753: PetscSectionSetNumFields(*section, numFields);
3754: if (numComp) {
3755: for(PetscInt f = 0; f < numFields; ++f) {
3756: PetscSectionSetFieldComponents(*section, f, numComp[f]);
3757: }
3758: }
3759: } else {
3760: numFields = 0;
3761: }
3762: DMMeshGetChart(dm, &pStart, &pEnd);
3763: PetscSectionSetChart(*section, pStart, pEnd);
3764: for(PetscInt d = 0; d <= dim; ++d) {
3765: DMMeshGetDepthStratum(dm, d, &pStart, &pEnd);
3766: for(PetscInt p = pStart; p < pEnd; ++p) {
3767: for(PetscInt f = 0; f < numFields; ++f) {
3768: PetscSectionSetFieldDof(*section, p, f, numDof[f*(dim+1)+d]);
3769: }
3770: PetscSectionSetDof(*section, p, numDofTot[d]);
3771: }
3772: }
3773: if (numBC) {
3774: for(PetscInt f = 0; f <= numFields; ++f) {maxConstraints[f] = 0;}
3775: for(PetscInt bc = 0; bc < numBC; ++bc) {
3776: PetscInt field = 0;
3777: const PetscInt *idx;
3778: PetscInt n;
3780: if (numFields) {field = bcField[bc];}
3781: ISGetLocalSize(bcPoints[bc], &n);
3782: ISGetIndices(bcPoints[bc], &idx);
3783: for(PetscInt i = 0; i < n; ++i) {
3784: const PetscInt p = idx[i];
3785: PetscInt depth, numConst;
3787: DMMeshGetLabelValue(dm, "depth", p, &depth);
3788: numConst = numDof[field*(dim+1)+depth];
3789: maxConstraints[field] = PetscMax(maxConstraints[field], numConst);
3790: if (numFields) {
3791: PetscSectionSetFieldConstraintDof(*section, p, field, numConst);
3792: }
3793: PetscSectionAddConstraintDof(*section, p, numConst);
3794: }
3795: ISRestoreIndices(bcPoints[bc], &idx);
3796: }
3797: for(PetscInt f = 0; f < numFields; ++f) {
3798: maxConstraints[numFields] += maxConstraints[f];
3799: }
3800: }
3801: PetscSectionSetUp(*section);
3802: if (maxConstraints[numFields]) {
3803: PetscInt *indices;
3805: PetscMalloc(maxConstraints[numFields] * sizeof(PetscInt), &indices);
3806: PetscSectionGetChart(*section, &pStart, &pEnd);
3807: for(PetscInt p = pStart; p < pEnd; ++p) {
3808: PetscInt cDof;
3810: PetscSectionGetConstraintDof(*section, p, &cDof);
3811: if (cDof) {
3812: if (cDof > maxConstraints[numFields]) {SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_LIB, "Likely memory corruption, point %d cDof %d > maxConstraints %d", p, cDof, maxConstraints);}
3813: if (numFields) {
3814: PetscInt numConst = 0, fOff = 0;
3816: for(PetscInt f = 0; f < numFields; ++f) {
3817: PetscInt cfDof, fDof;
3819: PetscSectionGetFieldDof(*section, p, f, &fDof);
3820: PetscSectionGetFieldConstraintDof(*section, p, f, &cfDof);
3821: for(PetscInt d = 0; d < cfDof; ++d) {
3822: indices[numConst+d] = fOff+d;
3823: }
3824: PetscSectionSetFieldConstraintIndices(*section, p, f, &indices[numConst]);
3825: numConst += cfDof;
3826: fOff += fDof;
3827: }
3828: if (cDof != numConst) {SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_LIB, "Total number of field constraints %d should be %d", numConst, cDof);}
3829: } else {
3830: for(PetscInt d = 0; d < cDof; ++d) {
3831: indices[d] = d;
3832: }
3833: }
3834: PetscSectionSetConstraintIndices(*section, p, indices);
3835: }
3836: }
3837: PetscFree(indices);
3838: }
3839: PetscFree2(numDofTot,maxConstraints);
3840: {
3841: PetscBool view = PETSC_FALSE;
3843: PetscOptionsHasName(((PetscObject) dm)->prefix, "-section_view", &view);
3844: if (view) {
3845: PetscSectionView(*section, PETSC_VIEWER_STDOUT_WORLD);
3846: }
3847: }
3848: return(0);
3849: }
3853: PetscErrorCode DMMeshConvertSection(const ALE::Obj<PETSC_MESH_TYPE>& mesh, const Obj<PETSC_MESH_TYPE::real_section_type>& s, PetscSection *section)
3854: {
3855: const PetscInt pStart = s->getChart().min();
3856: const PetscInt pEnd = s->getChart().max();
3857: PetscInt numFields = s->getNumSpaces();
3861: PetscSectionCreate(mesh->comm(), section);
3862: if (numFields) {
3863: PetscSectionSetNumFields(*section, numFields);
3864: for(PetscInt f = 0; f < numFields; ++f) {
3865: PetscSectionSetFieldComponents(*section, f, s->getSpaceComponents(f));
3866: }
3867: }
3868: PetscSectionSetChart(*section, pStart, pEnd);
3869: for(PetscInt p = pStart; p < pEnd; ++p) {
3870: PetscSectionSetDof(*section, p, s->getFiberDimension(p));
3871: for(PetscInt f = 0; f < numFields; ++f) {
3872: PetscSectionSetFieldDof(*section, p, f, s->getFiberDimension(p, f));
3873: }
3874: PetscSectionSetConstraintDof(*section, p, s->getConstraintDimension(p));
3875: for(PetscInt f = 0; f < numFields; ++f) {
3876: PetscSectionSetFieldConstraintDof(*section, p, f, s->getConstraintDimension(p, f));
3877: }
3878: }
3879: PetscSectionSetUp(*section);
3880: for(PetscInt p = pStart; p < pEnd; ++p) {
3881: PetscSectionSetConstraintIndices(*section, p, (PetscInt *) s->getConstraintDof(p));
3882: for(PetscInt f = 0; f < numFields; ++f) {
3883: PetscSectionSetFieldConstraintIndices(*section, p, f, (PetscInt *) s->getConstraintDof(p, f));
3884: }
3885: }
3886: return(0);
3887: }
3891: PetscErrorCode DMMeshGetSection(DM dm, const char name[], PetscSection *section) {
3892: ALE::Obj<PETSC_MESH_TYPE> mesh;
3896: DMMeshGetMesh(dm, mesh);
3897: DMMeshConvertSection(mesh, mesh->getRealSection(name), section);
3898: return(0);
3899: }
3903: PetscErrorCode DMMeshSetSection(DM dm, const char name[], PetscSection section) {
3904: ALE::Obj<PETSC_MESH_TYPE> mesh;
3908: DMMeshGetMesh(dm, mesh);
3909: {
3910: const Obj<PETSC_MESH_TYPE::real_section_type>& s = mesh->getRealSection(name);
3911: PetscInt pStart, pEnd, numFields;
3913: PetscSectionGetChart(section, &pStart, &pEnd);
3914: s->setChart(PETSC_MESH_TYPE::real_section_type::chart_type(pStart, pEnd));
3915: PetscSectionGetNumFields(section, &numFields);
3916: for(PetscInt f = 0; f < numFields; ++f) {
3917: PetscInt comp;
3918: PetscSectionGetFieldComponents(section, f, &comp);
3919: s->addSpace(comp);
3920: }
3921: for(PetscInt p = pStart; p < pEnd; ++p) {
3922: PetscInt fDim, cDim;
3924: PetscSectionGetDof(section, p, &fDim);
3925: s->setFiberDimension(p, fDim);
3926: for(PetscInt f = 0; f < numFields; ++f) {
3927: PetscSectionGetFieldDof(section, p, f, &fDim);
3928: s->setFiberDimension(p, fDim, f);
3929: }
3930: PetscSectionGetConstraintDof(section, p, &cDim);
3931: if (cDim) {
3932: s->setConstraintDimension(p, cDim);
3933: for(PetscInt f = 0; f < numFields; ++f) {
3934: PetscSectionGetFieldConstraintDof(section, p, f, &cDim);
3935: s->setConstraintDimension(p, cDim, f);
3936: }
3937: }
3938: }
3939: s->allocatePoint();
3940: for(PetscInt p = pStart; p < pEnd; ++p) {
3941: PetscInt *indices;
3943: PetscSectionGetConstraintIndices(section, p, &indices);
3944: s->setConstraintDof(p, indices);
3945: for(PetscInt f = 0; f < numFields; ++f) {
3946: PetscSectionGetFieldConstraintIndices(section, p, f, &indices);
3947: s->setConstraintDof(p, indices, f);
3948: }
3949: }
3950: {
3951: PetscBool isDefault;
3953: PetscStrcmp(name, "default", &isDefault);
3954: if (isDefault) {
3955: PetscInt maxDof = 0;
3957: for(PetscInt p = pStart; p < pEnd; ++p) {
3958: PetscInt fDim;
3960: PetscSectionGetDof(section, p, &fDim);
3961: maxDof = PetscMax(maxDof, fDim);
3962: }
3963: mesh->setMaxDof(maxDof);
3964: }
3965: }
3966: }
3967: return(0);
3968: }
3972: /*
3973: Note: This gets a borrowed reference, so the user should not destroy this PetscSection.
3974: */
3975: PetscErrorCode DMMeshGetDefaultSection(DM dm, PetscSection *section) {
3976: DM_Mesh *mesh = (DM_Mesh *) dm->data;
3980: if (!mesh->defaultSection && !mesh->useNewImpl) {
3981: DMMeshGetSection(dm, "default", &mesh->defaultSection);
3982: }
3983: *section = mesh->defaultSection;
3984: return(0);
3985: }
3989: /*
3990: Note: This reference will be stolen, so the user should not destroy this PetscSection.
3991: */
3992: PetscErrorCode DMMeshSetDefaultSection(DM dm, PetscSection section) {
3993: DM_Mesh *mesh = (DM_Mesh *) dm->data;
3997: mesh->defaultSection = section;
3998: if (!mesh->useNewImpl) {
3999: DMMeshSetSection(dm, "default", mesh->defaultSection);
4000: }
4001: return(0);
4002: }
4006: PetscErrorCode DMMeshGetCoordinateSection(DM dm, PetscSection *section) {
4007: DM_Mesh *mesh = (DM_Mesh *) dm->data;
4013: if (mesh->useNewImpl) {
4014: *section = mesh->coordSection;
4015: } else {
4016: DMMeshGetSection(dm, "coordinates", section);
4017: }
4018: return(0);
4019: }
4023: PetscErrorCode DMMeshSetCoordinateSection(DM dm, PetscSection section) {
4027: DMMeshSetSection(dm, "coordinates", section);
4028: return(0);
4029: }
4033: PetscErrorCode DMMeshGetConeSection(DM dm, PetscSection *section) {
4034: DM_Mesh *mesh = (DM_Mesh *) dm->data;
4038: if (!mesh->useNewImpl) {SETERRQ(((PetscObject) dm)->comm, PETSC_ERR_ARG_WRONG, "This method is only valid for C implementation meshes.");}
4039: if (section) *section = mesh->coneSection;
4040: return(0);
4041: }
4045: PetscErrorCode DMMeshGetCones(DM dm, PetscInt *cones[]) {
4046: DM_Mesh *mesh = (DM_Mesh *) dm->data;
4050: if (!mesh->useNewImpl) {SETERRQ(((PetscObject) dm)->comm, PETSC_ERR_ARG_WRONG, "This method is only valid for C implementation meshes.");}
4051: if (cones) *cones = mesh->cones;
4052: return(0);
4053: }
4057: PetscErrorCode DMMeshCreateConeSection(DM dm, PetscSection *section) {
4058: ALE::Obj<PETSC_MESH_TYPE> mesh;
4059: PetscInt p;
4063: DMMeshGetMesh(dm, mesh);
4064: PetscSectionCreate(((PetscObject) dm)->comm, section);
4065: PetscSectionSetChart(*section, mesh->getSieve()->getChart().min(), mesh->getSieve()->getChart().max());
4066: for(p = mesh->getSieve()->getChart().min(); p < mesh->getSieve()->getChart().max(); ++p) {
4067: PetscSectionSetDof(*section, p, mesh->getSieve()->getConeSize(p));
4068: }
4069: PetscSectionSetUp(*section);
4070: return(0);
4071: }
4075: PetscErrorCode DMMeshGetCoordinateVec(DM dm, Vec *coordinates) {
4076: DM_Mesh *mesh = (DM_Mesh *) dm->data;
4082: if (mesh->useNewImpl) {
4083: *coordinates = mesh->coordinates;
4084: } else {
4085: ALE::Obj<PETSC_MESH_TYPE> mesh;
4086: DMMeshGetMesh(dm, mesh);
4087: const Obj<PETSC_MESH_TYPE::real_section_type>& coords = mesh->getRealSection("coordinates");
4088: VecCreateSeqWithArray(PETSC_COMM_SELF,1, coords->getStorageSize(), coords->restrictSpace(), coordinates);
4089: }
4090: return(0);
4091: }
4095: PetscErrorCode DMMeshComputeCellGeometry(DM dm, PetscInt cell, PetscReal *v0, PetscReal *J, PetscReal *invJ, PetscReal *detJ) {
4096: ALE::Obj<PETSC_MESH_TYPE> mesh;
4100: DMMeshGetMesh(dm, mesh);
4101: {
4102: ALE::Obj<PETSC_MESH_TYPE::real_section_type> coordinates = mesh->getRealSection("coordinates");
4104: mesh->computeElementGeometry(coordinates, cell, v0, J, invJ, *detJ);
4105: }
4106: return(0);
4107: }
4111: PetscErrorCode DMMeshVecGetClosure(DM dm, Vec v, PetscInt point, const PetscScalar *values[]) {
4112: ALE::Obj<PETSC_MESH_TYPE> mesh;
4116: #ifdef PETSC_USE_COMPLEX
4117: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP, "DMMesh does not support complex closure");
4118: #else
4119: DMMeshGetMesh(dm, mesh);
4120: /* Peeling back IMesh::restrictClosure() */
4121: try {
4122: typedef ALE::ISieveVisitor::RestrictVecVisitor<PetscScalar> visitor_type;
4123: PetscSection section;
4124: PetscScalar *array;
4125: PetscInt numFields;
4127: DMMeshGetDefaultSection(dm, §ion);
4128: PetscSectionGetNumFields(section, &numFields);
4129: const PetscInt size = mesh->sizeWithBC(section, point); /* OPT: This can be precomputed */
4130: DMGetWorkArray(dm, 2*size+numFields+1, &array);
4131: visitor_type rV(v, section, size, array, (PetscInt *) &array[2*size], (PetscInt *) &array[size]);
4132: if (mesh->depth() == 1) {
4133: rV.visitPoint(point, 0);
4134: // Cone is guarateed to be ordered correctly
4135: mesh->getSieve()->orientedCone(point, rV);
4136: } else {
4137: ALE::ISieveVisitor::PointRetriever<PETSC_MESH_TYPE::sieve_type,visitor_type> pV((int) pow((double) mesh->getSieve()->getMaxConeSize(), mesh->depth())+1, rV, true);
4139: ALE::ISieveTraversal<PETSC_MESH_TYPE::sieve_type>::orientedClosure(*mesh->getSieve(), point, pV);
4140: }
4141: *values = rV.getValues();
4142: } catch(ALE::Exception e) {
4143: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid argument: %s", e.message());
4144: } catch(PETSc::Exception e) {
4145: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid argument: %s", e.message());
4146: }
4147: #endif
4148: return(0);
4149: }
4153: PetscErrorCode DMMeshVecSetClosure(DM dm, Vec v, PetscInt point, const PetscScalar values[], InsertMode mode) {
4154: ALE::Obj<PETSC_MESH_TYPE> mesh;
4158: #ifdef PETSC_USE_COMPLEX
4159: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP, "DMMesh does not support complex closure");
4160: #else
4161: DMMeshGetMesh(dm, mesh);
4162: /* Peeling back IMesh::update() and IMesh::updateAdd() */
4163: try {
4164: typedef ALE::ISieveVisitor::UpdateVecVisitor<PetscScalar> visitor_type;
4165: PetscSection section;
4166: PetscInt *fieldSize;
4167: PetscInt numFields;
4169: DMMeshGetDefaultSection(dm, §ion);
4170: PetscSectionGetNumFields(section, &numFields);
4171: DMGetWorkArray(dm, numFields, (PetscScalar **) &fieldSize);
4172: mesh->sizeWithBC(section, point, fieldSize); /* OPT: This can be precomputed */
4173: visitor_type uV(v, section, values, mode, numFields, fieldSize);
4174: if (mesh->depth() == 1) {
4175: uV.visitPoint(point, 0);
4176: // Cone is guarateed to be ordered correctly
4177: mesh->getSieve()->orientedCone(point, uV);
4178: } else {
4179: ALE::ISieveVisitor::PointRetriever<PETSC_MESH_TYPE::sieve_type,visitor_type> pV((int) pow((double) mesh->getSieve()->getMaxConeSize(), mesh->depth())+1, uV, true);
4181: ALE::ISieveTraversal<PETSC_MESH_TYPE::sieve_type>::orientedClosure(*mesh->getSieve(), point, pV);
4182: }
4183: } catch(ALE::Exception e) {
4184: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid argument: %s", e.message());
4185: }
4186: #endif
4187: return(0);
4188: }
4192: PetscErrorCode DMMeshMatSetClosure(DM dm, Mat A, PetscInt point, PetscScalar values[], InsertMode mode) {
4193: ALE::Obj<PETSC_MESH_TYPE> mesh;
4197: #ifdef PETSC_USE_COMPLEX
4198: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP, "DMMesh does not support complex closure");
4199: #else
4200: DMMeshGetMesh(dm, mesh);
4201: /* Copying from updateOperator() */
4202: try {
4203: typedef ALE::ISieveVisitor::IndicesVisitor<PetscSection,PETSC_MESH_TYPE::order_type,PetscInt> visitor_type;
4204: ALE::Obj<PETSC_MESH_TYPE::real_section_type> s = mesh->getRealSection("default");
4205: const ALE::Obj<PETSC_MESH_TYPE::order_type>& globalOrder = mesh->getFactory()->getGlobalOrder(mesh, s->getName(), s);
4206: PetscSection section;
4207: PetscInt numFields;
4208: PetscInt *fieldSize = PETSC_NULL;
4210: DMMeshGetDefaultSection(dm, §ion);
4211: PetscSectionGetNumFields(section, &numFields);
4212: if (numFields) {
4213: DMGetWorkArray(dm, numFields, (PetscScalar **) &fieldSize);
4214: mesh->sizeWithBC(section, point, fieldSize); /* OPT: This can be precomputed */
4215: }
4216: visitor_type iV(section, *globalOrder, (int) pow((double) mesh->getSieve()->getMaxConeSize(), mesh->depth())*mesh->getMaxDof(), mesh->depth() > 1, fieldSize);
4218: updateOperator(A, *mesh->getSieve(), iV, point, values, mode);
4219: } catch(ALE::Exception e) {
4220: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid argument: %s", e.message());
4221: }
4222: #endif
4223: return(0);
4224: }
4228: /*@C
4229: DMMeshHasSectionReal - Determines whether this mesh has a SectionReal with the given name.
4231: Not Collective
4233: Input Parameters:
4234: + mesh - The DMMesh object
4235: - name - The section name
4237: Output Parameter:
4238: . flag - True if the SectionReal is present in the DMMesh
4240: Level: intermediate
4242: .keywords: mesh, elements
4243: .seealso: DMMeshCreate()
4244: @*/
4245: PetscErrorCode DMMeshHasSectionReal(DM dm, const char name[], PetscBool *flag)
4246: {
4247: ALE::Obj<PETSC_MESH_TYPE> m;
4251: DMMeshGetMesh(dm, m);
4252: *flag = (PetscBool) m->hasRealSection(std::string(name));
4253: return(0);
4254: }
4258: /*@C
4259: DMMeshGetSectionReal - Returns a SectionReal of the given name from the DMMesh.
4261: Collective on DMMesh
4263: Input Parameters:
4264: + mesh - The DMMesh object
4265: - name - The section name
4267: Output Parameter:
4268: . section - The SectionReal
4270: Note: The section is a new object, and must be destroyed by the user
4272: Level: intermediate
4274: .keywords: mesh, elements
4276: .seealso: DMMeshCreate(), SectionRealDestroy()
4277: @*/
4278: PetscErrorCode DMMeshGetSectionReal(DM dm, const char name[], SectionReal *section)
4279: {
4280: ALE::Obj<PETSC_MESH_TYPE> m;
4281: bool has;
4282: PetscErrorCode ierr;
4285: DMMeshGetMesh(dm, m);
4286: SectionRealCreate(m->comm(), section);
4287: PetscObjectSetName((PetscObject) *section, name);
4288: has = m->hasRealSection(std::string(name));
4289: SectionRealSetSection(*section, m->getRealSection(std::string(name)));
4290: SectionRealSetBundle(*section, m);
4291: if (!has) {
4292: m->getRealSection(std::string(name))->setChart(m->getSieve()->getChart());
4293: }
4294: return(0);
4295: }
4299: /*@C
4300: DMMeshSetSectionReal - Puts a SectionReal of the given name into the DMMesh.
4302: Collective on DMMesh
4304: Input Parameters:
4305: + mesh - The DMMesh object
4306: - section - The SectionReal
4308: Note: This takes the section name from the PETSc object
4310: Level: intermediate
4312: .keywords: mesh, elements
4313: .seealso: DMMeshCreate()
4314: @*/
4315: PetscErrorCode DMMeshSetSectionReal(DM dm, const char name[], SectionReal section)
4316: {
4317: ALE::Obj<PETSC_MESH_TYPE> m;
4318: ALE::Obj<PETSC_MESH_TYPE::real_section_type> s;
4322: DMMeshGetMesh(dm, m);
4323: PetscObjectGetName((PetscObject) section, &name);
4324: SectionRealGetSection(section, s);
4325: m->setRealSection(std::string(name), s);
4326: return(0);
4327: }
4331: /*@C
4332: DMMeshHasSectionInt - Determines whether this mesh has a SectionInt with the given name.
4334: Not Collective
4336: Input Parameters:
4337: + mesh - The DMMesh object
4338: - name - The section name
4340: Output Parameter:
4341: . flag - True if the SectionInt is present in the DMMesh
4343: Level: intermediate
4345: .keywords: mesh, elements
4346: .seealso: DMMeshCreate()
4347: @*/
4348: PetscErrorCode DMMeshHasSectionInt(DM dm, const char name[], PetscBool *flag)
4349: {
4350: ALE::Obj<PETSC_MESH_TYPE> m;
4354: DMMeshGetMesh(dm, m);
4355: *flag = (PetscBool) m->hasIntSection(std::string(name));
4356: return(0);
4357: }
4361: /*@C
4362: DMMeshGetSectionInt - Returns a SectionInt of the given name from the DMMesh.
4364: Collective on DMMesh
4366: Input Parameters:
4367: + mesh - The DMMesh object
4368: - name - The section name
4370: Output Parameter:
4371: . section - The SectionInt
4373: Note: The section is a new object, and must be destroyed by the user
4375: Level: intermediate
4377: .keywords: mesh, elements
4378: .seealso: DMMeshCreate()
4379: @*/
4380: PetscErrorCode DMMeshGetSectionInt(DM dm, const char name[], SectionInt *section)
4381: {
4382: ALE::Obj<PETSC_MESH_TYPE> m;
4383: bool has;
4384: PetscErrorCode ierr;
4387: DMMeshGetMesh(dm, m);
4388: SectionIntCreate(m->comm(), section);
4389: PetscObjectSetName((PetscObject) *section, name);
4390: has = m->hasIntSection(std::string(name));
4391: SectionIntSetSection(*section, m->getIntSection(std::string(name)));
4392: SectionIntSetBundle(*section, m);
4393: if (!has) {
4394: m->getIntSection(std::string(name))->setChart(m->getSieve()->getChart());
4395: }
4396: return(0);
4397: }
4401: /*@C
4402: DMMeshSetSectionInt - Puts a SectionInt of the given name into the DMMesh.
4404: Collective on DMMesh
4406: Input Parameters:
4407: + mesh - The DMMesh object
4408: - section - The SectionInt
4410: Note: This takes the section name from the PETSc object
4412: Level: intermediate
4414: .keywords: mesh, elements
4415: .seealso: DMMeshCreate()
4416: @*/
4417: PetscErrorCode DMMeshSetSectionInt(DM dm, SectionInt section)
4418: {
4419: ALE::Obj<PETSC_MESH_TYPE> m;
4420: ALE::Obj<PETSC_MESH_TYPE::int_section_type> s;
4421: const char *name;
4425: DMMeshGetMesh(dm, m);
4426: PetscObjectGetName((PetscObject) section, &name);
4427: SectionIntGetSection(section, s);
4428: m->setIntSection(std::string(name), s);
4429: return(0);
4430: }
4434: /*@C
4435: SectionGetArray - Returns the array underlying the Section.
4437: Not Collective
4439: Input Parameters:
4440: + mesh - The DMMesh object
4441: - name - The section name
4443: Output Parameters:
4444: + numElements - The number of mesh element with values
4445: . fiberDim - The number of values per element
4446: - array - The array
4448: Level: intermediate
4450: .keywords: mesh, elements
4451: .seealso: DMMeshCreate()
4452: @*/
4453: PetscErrorCode SectionGetArray(DM dm, const char name[], PetscInt *numElements, PetscInt *fiberDim, PetscScalar *array[])
4454: {
4455: ALE::Obj<PETSC_MESH_TYPE> m;
4456: PetscErrorCode ierr;
4459: DMMeshGetMesh(dm, m);
4460: const Obj<PETSC_MESH_TYPE::real_section_type>& section = m->getRealSection(std::string(name));
4461: if (section->size() == 0) {
4462: *numElements = 0;
4463: *fiberDim = 0;
4464: *array = NULL;
4465: return(0);
4466: }
4467: const PETSC_MESH_TYPE::real_section_type::chart_type& chart = section->getChart();
4468: /* const int depth = m->depth(*chart.begin()); */
4469: /* *numElements = m->depthStratum(depth)->size(); */
4470: /* *fiberDim = section->getFiberDimension(*chart.begin()); */
4471: /* *array = (PetscScalar *) m->restrict(section); */
4472: int fiberDimMin = section->getFiberDimension(*chart.begin());
4473: int numElem = 0;
4475: for(PETSC_MESH_TYPE::real_section_type::chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
4476: const int fiberDim = section->getFiberDimension(*c_iter);
4478: if (fiberDim < fiberDimMin) fiberDimMin = fiberDim;
4479: }
4480: for(PETSC_MESH_TYPE::real_section_type::chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
4481: const int fiberDim = section->getFiberDimension(*c_iter);
4483: numElem += fiberDim/fiberDimMin;
4484: }
4485: *numElements = numElem;
4486: *fiberDim = fiberDimMin;
4487: *array = (PetscScalar *) section->restrictSpace();
4488: return(0);
4489: }
4493: inline void ExpandInterval(const ALE::Point& interval, int indices[], int& indx)
4494: {
4495: const int end = interval.prefix + interval.index;
4496: for(int i = interval.index; i < end; i++) {
4497: indices[indx++] = i;
4498: }
4499: }
4503: inline void ExpandInterval_New(ALE::Point interval, PetscInt indices[], PetscInt *indx)
4504: {
4505: for(int i = 0; i < interval.prefix; i++) {
4506: indices[(*indx)++] = interval.index + i;
4507: }
4508: for(int i = 0; i < -interval.prefix; i++) {
4509: indices[(*indx)++] = -1;
4510: }
4511: }