Actual source code: mesh.c
petsc-3.4.5 2014-06-29
1: #include <petsc-private/meshimpl.h> /*I "petscdmmesh.h" I*/
2: #include <petscdmmesh_viewers.hh>
3: #include <petscdmmesh_formats.hh>
4: #include <petscsf.h>
6: /* Logging support */
7: PetscLogEvent DMMesh_View, DMMesh_GetGlobalScatter, DMMesh_restrictVector, DMMesh_assembleVector, DMMesh_assembleVectorComplete, DMMesh_assembleMatrix, DMMesh_updateOperator;
9: ALE::MemoryLogger Petsc_MemoryLogger;
11: EXTERN_C_BEGIN
14: /*
15: Private routine to delete internal tag storage when a communicator is freed.
17: This is called by MPI, not by users.
19: Note: this is declared extern "C" because it is passed to MPI_Keyval_create
21: we do not use PetscFree() since it is unsafe after PetscFinalize()
22: */
23: PetscMPIInt DMMesh_DelTag(MPI_Comm comm,PetscMPIInt keyval,void *attr_val,void *extra_state)
24: {
25: free(attr_val);
26: return(MPI_SUCCESS);
27: }
28: EXTERN_C_END
32: PetscErrorCode DMMeshFinalize()
33: {
35: PETSC_MESH_TYPE::MeshNumberingFactory::singleton(0, 0, true);
36: return(0);
37: }
41: /*@C
42: DMMeshGetMesh - Gets the internal mesh object
44: Not collective
46: Input Parameter:
47: . mesh - the mesh object
49: Output Parameter:
50: . m - the internal mesh object
52: Level: advanced
54: .seealso DMMeshCreate(), DMMeshSetMesh()
55: @*/
56: PetscErrorCode DMMeshGetMesh(DM dm, ALE::Obj<PETSC_MESH_TYPE>& m)
57: {
58: DM_Mesh *mesh = (DM_Mesh*) dm->data;
62: if (mesh->useNewImpl) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "This method is only valid for C++ implementation meshes.");
63: m = mesh->m;
64: return(0);
65: }
69: /*@C
70: DMMeshSetMesh - Sets the internal mesh object
72: Not collective
74: Input Parameters:
75: + mesh - the mesh object
76: - m - the internal mesh object
78: Level: advanced
80: .seealso DMMeshCreate(), DMMeshGetMesh()
81: @*/
82: PetscErrorCode DMMeshSetMesh(DM dm, const ALE::Obj<PETSC_MESH_TYPE>& m)
83: {
84: DM_Mesh *mesh = (DM_Mesh*) dm->data;
89: mesh->m = m;
90: VecScatterDestroy(&mesh->globalScatter);
91: return(0);
92: }
96: PetscErrorCode DMMeshView_Sieve_Ascii(const ALE::Obj<PETSC_MESH_TYPE>& mesh, PetscViewer viewer)
97: {
98: PetscViewerFormat format;
99: PetscErrorCode ierr;
102: PetscViewerGetFormat(viewer, &format);
103: if (format == PETSC_VIEWER_ASCII_VTK) {
104: VTKViewer::writeHeader(mesh, viewer);
105: VTKViewer::writeVertices(mesh, viewer);
106: VTKViewer::writeElements(mesh, viewer);
107: const ALE::Obj<PETSC_MESH_TYPE::int_section_type>& p = mesh->getIntSection("Partition");
108: const ALE::Obj<PETSC_MESH_TYPE::label_sequence>& cells = mesh->heightStratum(0);
109: const PETSC_MESH_TYPE::label_sequence::iterator end = cells->end();
110: const int rank = mesh->commRank();
112: p->setChart(PETSC_MESH_TYPE::int_section_type::chart_type(*cells));
113: p->setFiberDimension(cells, 1);
114: p->allocatePoint();
115: for (PETSC_MESH_TYPE::label_sequence::iterator c_iter = cells->begin(); c_iter != end; ++c_iter) {
116: p->updatePoint(*c_iter, &rank);
117: }
118: PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_VTK_CELL);
119: SectionView_Sieve_Ascii(mesh, p, "Partition", viewer);
120: PetscViewerPopFormat(viewer);
121: } else if (format == PETSC_VIEWER_ASCII_PCICE) {
122: char *filename;
123: char coordFilename[2048];
124: PetscBool isConnect;
125: size_t len;
127: PetscViewerFileGetName(viewer, (const char **) &filename);
128: PetscStrlen(filename, &len);
129: PetscStrcmp(&(filename[len-5]), ".lcon", &isConnect);
130: if (!isConnect) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG, "Invalid element connectivity filename: %s", filename);
132: ALE::PCICE::Viewer::writeElements(mesh, viewer);
133: PetscStrncpy(coordFilename, filename, len-5);
135: coordFilename[len-5] = '\0';
137: PetscStrcat(coordFilename, ".nodes");
138: PetscViewerFileSetName(viewer, coordFilename);
139: ALE::PCICE::Viewer::writeVertices(mesh, viewer);
140: } else if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
141: mesh->view("Mesh");
142: } else {
143: PetscInt dim = mesh->getDimension();
144: PetscInt size = mesh->commSize();
145: PetscInt locDepth = mesh->depth(), depth;
146: PetscInt num = 0;
147: PetscInt *sizes;
149: PetscMalloc(size * sizeof(PetscInt), &sizes);
150: PetscViewerASCIIPrintf(viewer, "Mesh in %d dimensions:\n", dim);
151: MPI_Allreduce(&locDepth, &depth, 1, MPIU_INT, MPI_MAX, mesh->comm());
152: if (depth == 1) {
153: num = mesh->depthStratum(0)->size();
154: MPI_Gather(&num, 1, MPIU_INT, sizes, 1, MPIU_INT, 0, mesh->comm());
155: PetscViewerASCIIPrintf(viewer, " %d-cells:", 0);
156: for (PetscInt p = 0; p < size; ++p) {PetscViewerASCIIPrintf(viewer, " %d", sizes[p]);}
157: PetscViewerASCIIPrintf(viewer, "\n");
158: num = mesh->heightStratum(0)->size();
159: MPI_Gather(&num, 1, MPIU_INT, sizes, 1, MPIU_INT, 0, mesh->comm());
160: PetscViewerASCIIPrintf(viewer, " %d-cells:", dim);
161: for (PetscInt p = 0; p < size; ++p) {PetscViewerASCIIPrintf(viewer, " %d", sizes[p]);}
162: PetscViewerASCIIPrintf(viewer, "\n");
163: } else {
164: for (int d = 0; d <= dim; d++) {
165: num = mesh->depthStratum(d)->size();
166: MPI_Gather(&num, 1, MPIU_INT, sizes, 1, MPIU_INT, 0, mesh->comm());
167: PetscViewerASCIIPrintf(viewer, " %d-cells:", d);
168: for (PetscInt p = 0; p < size; ++p) {PetscViewerASCIIPrintf(viewer, " %d", sizes[p]);}
169: PetscViewerASCIIPrintf(viewer, "\n");
170: }
171: }
172: PetscFree(sizes);
173: }
174: PetscViewerFlush(viewer);
175: return(0);
176: }
181: PetscErrorCode DMMeshView_Sieve_Binary(const ALE::Obj<PETSC_MESH_TYPE>& mesh, PetscViewer viewer)
182: {
183: char *filename;
187: PetscViewerFileGetName(viewer, (const char**) &filename);
188: ALE::MeshSerializer::writeMesh(filename, *mesh);
189: return(0);
190: }
194: PetscErrorCode DMMeshView_Sieve(const ALE::Obj<PETSC_MESH_TYPE>& mesh, PetscViewer viewer)
195: {
196: PetscBool iascii, isbinary, isdraw;
200: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
201: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERBINARY, &isbinary);
202: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERDRAW, &isdraw);
204: if (iascii) {
205: DMMeshView_Sieve_Ascii(mesh, viewer);
206: } else if (isbinary) {
207: DMMeshView_Sieve_Binary(mesh, viewer);
208: } else if (isdraw) SETERRQ(PetscObjectComm((PetscObject)viewer),PETSC_ERR_SUP, "Draw viewer not implemented for DMMesh");
209: else SETERRQ1(PetscObjectComm((PetscObject)viewer),PETSC_ERR_SUP,"Viewer type %s not supported by this mesh object", ((PetscObject)viewer)->type_name);
210: return(0);
211: }
215: PetscErrorCode DMMeshView_Mesh_Ascii(DM dm, PetscViewer viewer)
216: {
217: DM_Mesh *mesh = (DM_Mesh*) dm->data;
218: PetscViewerFormat format;
219: PetscErrorCode ierr;
222: PetscViewerGetFormat(viewer, &format);
223: if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
224: const char *name;
225: PetscInt maxConeSize, maxSupportSize;
226: PetscInt pStart, pEnd, p;
227: PetscMPIInt rank;
229: MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
230: PetscObjectGetName((PetscObject) dm, &name);
231: DMMeshGetChart(dm, &pStart, &pEnd);
232: DMMeshGetMaxSizes(dm, &maxConeSize, &maxSupportSize);
233: PetscViewerASCIISynchronizedAllow(viewer, PETSC_TRUE);
234: PetscViewerASCIIPrintf(viewer, "Mesh '%s':\n", name);
235: PetscViewerASCIISynchronizedPrintf(viewer, "Max sizes cone: %d support: %d\n", maxConeSize, maxSupportSize);
236: PetscViewerASCIIPrintf(viewer, "orientation is missing\n", name);
237: PetscViewerASCIIPrintf(viewer, "cap --> base:\n", name);
238: for (p = pStart; p < pEnd; ++p) {
239: PetscInt dof, off, s;
241: PetscSectionGetDof(mesh->supportSection, p, &dof);
242: PetscSectionGetOffset(mesh->supportSection, p, &off);
243: for (s = off; s < off+dof; ++s) {
244: PetscViewerASCIISynchronizedPrintf(viewer, "[%d]: %d ----> %d\n", rank, p, mesh->supports[s]);
245: }
246: }
247: PetscViewerFlush(viewer);
248: PetscViewerASCIIPrintf(viewer, "base <-- cap:\n", name);
249: for (p = pStart; p < pEnd; ++p) {
250: PetscInt dof, off, c;
252: PetscSectionGetDof(mesh->coneSection, p, &dof);
253: PetscSectionGetOffset(mesh->coneSection, p, &off);
254: for (c = off; c < off+dof; ++c) {
255: PetscViewerASCIISynchronizedPrintf(viewer, "[%d]: %d <---- %d\n", rank, p, mesh->cones[c]);
256: /* PetscViewerASCIISynchronizedPrintf(viewer, "[%d]: %d <---- %d: %d\n", rank, p, mesh->cones[c], mesh->coneOrientations[c]); */
257: }
258: }
259: PetscViewerFlush(viewer);
260: PetscSectionVecView(mesh->coordSection, mesh->coordinates, viewer);
261: } else {
262: MPI_Comm comm;
263: PetscInt *sizes;
264: PetscInt depth, dim, d;
265: PetscInt pStart, pEnd, p;
266: PetscMPIInt size;
268: PetscObjectGetComm((PetscObject)dm,&comm);
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: DMMeshView_Sieve(mesh->m, viewer);
327: }
328: }
329: return(0);
330: }
334: /*@
335: DMMeshLoad - Create a mesh topology from the saved data in a viewer.
337: Collective on Viewer
339: Input Parameter:
340: . viewer - The viewer containing the data
342: Output Parameters:
343: . mesh - the mesh object
345: Level: advanced
347: .seealso DMView()
348: @*/
349: PetscErrorCode DMMeshLoad(PetscViewer viewer, DM dm)
350: {
351: DM_Mesh *mesh = (DM_Mesh*) dm->data;
352: char *filename;
356: if (!mesh->m) {
357: MPI_Comm comm;
359: PetscObjectGetComm((PetscObject) viewer, &comm);
360: ALE::Obj<PETSC_MESH_TYPE> m = new PETSC_MESH_TYPE(comm, 1);
361: DMMeshSetMesh(dm, m);
362: }
363: PetscViewerFileGetName(viewer, (const char**) &filename);
364: ALE::MeshSerializer::loadMesh(filename, *mesh->m);
365: return(0);
366: }
370: PetscErrorCode GetAdjacentDof_Private()
371: {
373: return(0);
374: }
378: PetscErrorCode DMCreateMatrix_Mesh(DM dm, MatType mtype, Mat *J)
379: {
380: DM_Mesh *mesh = (DM_Mesh*) dm->data;
381: ISLocalToGlobalMapping ltog;
382: PetscErrorCode ierr;
385: #if !defined(PETSC_USE_DYNAMIC_LIBRARIES)
386: MatInitializePackage();
387: #endif
388: if (!mtype) mtype = MATAIJ;
389: if (mesh->useNewImpl) {
390: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not supported");
391: } else {
392: ALE::Obj<PETSC_MESH_TYPE> m;
393: ALE::Obj<PETSC_MESH_TYPE::real_section_type> s;
394: SectionReal section;
395: PetscBool flag;
396: DMMeshHasSectionReal(dm, "default", &flag);
397: if (!flag) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "Must set default section");
398: DMMeshGetSectionReal(dm, "default", §ion);
399: DMMeshGetMesh(dm, m);
400: SectionRealGetSection(section, s);
401: try {
402: DMMeshCreateMatrix(m, s, mtype, J, -1, !dm->prealloc_only);
403: } catch(ALE::Exception e) {
404: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_LIB, e.message());
405: }
406: SectionRealDestroy(§ion);
407: }
408: MatSetDM(*J, dm);
409: DMGetLocalToGlobalMapping(dm, <og);
410: MatSetLocalToGlobalMapping(*J, ltog, ltog);
411: return(0);
412: }
416: /*@C
417: DMMeshCreateMatrix - Creates a matrix with the correct parallel layout required for
418: computing the Jacobian on a function defined using the information in the Section.
420: Collective on DMMesh
422: Input Parameters:
423: + mesh - the mesh object
424: . section - the section which determines data layout
425: - mtype - Supported types are MATSEQAIJ, MATMPIAIJ, MATSEQBAIJ, MATMPIBAIJ, MATSEQSBAIJ, MATMPISBAIJ,
426: or any type which inherits from one of these (such as MATAIJ, MATLUSOL, etc.).
428: Output Parameter:
429: . J - matrix with the correct nonzero preallocation
430: (obviously without the correct Jacobian values)
432: Level: advanced
434: Notes: This properly preallocates the number of nonzeros in the sparse matrix so you
435: do not need to do it yourself.
437: .seealso ISColoringView(), ISColoringGetIS(), MatFDColoringCreate(), DMDASetBlockFills()
438: @*/
439: PetscErrorCode DMMeshCreateMatrix(DM dm, SectionReal section, MatType mtype, Mat *J)
440: {
441: ALE::Obj<PETSC_MESH_TYPE> m;
442: ALE::Obj<PETSC_MESH_TYPE::real_section_type> s;
443: PetscErrorCode ierr;
446: #if !defined(PETSC_USE_DYNAMIC_LIBRARIES)
447: MatInitializePackage();
448: #endif
449: if (!mtype) mtype = MATAIJ;
450: DMMeshGetMesh(dm, m);
451: SectionRealGetSection(section, s);
452: try {
453: DMMeshCreateMatrix(m, s, mtype, J, -1, !dm->prealloc_only);
454: } catch(ALE::Exception e) {
455: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB, e.message());
456: }
457: MatSetDM(*J, dm);
458: return(0);
459: }
463: /*@
464: DMMeshCreateVector - Creates a global vector matching the input section
466: Collective on DMMesh
468: Input Parameters:
469: + mesh - the DMMesh
470: - section - the Section
472: Output Parameter:
473: . vec - the global vector
475: Level: advanced
477: Notes: The vector can safely be destroyed using VecDestroy().
478: .seealso DMMeshCreate()
479: @*/
480: PetscErrorCode DMMeshCreateVector(DM mesh, SectionReal section, Vec *vec)
481: {
482: ALE::Obj<PETSC_MESH_TYPE> m;
483: ALE::Obj<PETSC_MESH_TYPE::real_section_type> s;
484: PetscErrorCode ierr;
487: DMMeshGetMesh(mesh, m);
488: SectionRealGetSection(section, s);
489: const ALE::Obj<PETSC_MESH_TYPE::order_type>& order = m->getFactory()->getGlobalOrder(m, s->getName(), s);
491: VecCreate(m->comm(), vec);
492: VecSetSizes(*vec, order->getLocalSize(), order->getGlobalSize());
493: VecSetFromOptions(*vec);
494: return(0);
495: }
499: PetscErrorCode DMDestroy_Mesh(DM dm)
500: {
501: DM_Mesh *mesh = (DM_Mesh*) dm->data;
502: SieveLabel next = mesh->labels;
506: mesh->m = NULL;
507: PetscSectionDestroy(&mesh->defaultSection);
508: VecScatterDestroy(&mesh->globalScatter);
510: /* NEW_MESH_IMPL */
511: PetscSFDestroy(&mesh->sf);
512: PetscSectionDestroy(&mesh->coneSection);
513: PetscFree(mesh->cones);
514: PetscSectionDestroy(&mesh->supportSection);
515: PetscFree(mesh->supports);
516: PetscSectionDestroy(&mesh->coordSection);
517: VecDestroy(&mesh->coordinates);
518: PetscFree2(mesh->meetTmpA,mesh->meetTmpB);
519: PetscFree2(mesh->joinTmpA,mesh->joinTmpB);
520: PetscFree2(mesh->closureTmpA,mesh->closureTmpB);
521: while (next) {
522: SieveLabel tmp;
524: PetscFree(next->name);
525: PetscFree(next->stratumValues);
526: PetscFree(next->stratumOffsets);
527: PetscFree(next->stratumSizes);
528: PetscFree(next->points);
529: tmp = next->next;
530: PetscFree(next);
531: next = tmp;
532: }
533: /* This was originally freed in DMDestroy(), but that prevents reference counting of backend objects */
534: PetscFree(mesh);
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 DMPLEX */
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(PetscObjectComm((PetscObject)dm),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(PetscObjectComm((PetscObject)dm), gvec);
566: VecSetSizes(*gvec, localSize, globalSize);
567: VecSetFromOptions(*gvec);
568: VecSetDM(*gvec, 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(PetscObjectComm((PetscObject)dm),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: VecSetDM(*lvec, dm);
599: return(0);
600: }
604: PetscErrorCode DMGetLocalToGlobalMapping_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;
610: PetscErrorCode ierr;
611: PetscInt *ltog;
614: DMMeshHasSectionReal(dm, "default", &flag);
615: if (!flag) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE, "Must set default section");
616: DMMeshGetSectionReal(dm,"default", §ion);
617: DMMeshGetMesh(dm, m);
618: SectionRealGetSection(section, s);
619: const ALE::Obj<PETSC_MESH_TYPE::order_type>& globalOrder = m->getFactory()->getGlobalOrder(m, s->getName(), s);
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;
658: PetscErrorCode ierr;
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) *dim = mesh->dim;
833: else {
834: Obj<PETSC_MESH_TYPE> m;
835: DMMeshGetMesh(dm, m);
836: *dim = m->getDimension();
837: }
838: return(0);
839: }
843: /*@
844: DMMeshSetDimension - Set the topological mesh dimension
846: Collective on mesh
848: Input Parameters:
849: + mesh - The DMMesh
850: - dim - The topological mesh dimension
852: Level: beginner
854: .seealso: DMMeshCreate()
855: @*/
856: PetscErrorCode DMMeshSetDimension(DM dm, PetscInt dim)
857: {
858: DM_Mesh *mesh = (DM_Mesh*) dm->data;
863: if (mesh->useNewImpl) mesh->dim = dim;
864: else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Cannot reset dimension of C++ mesh");
865: return(0);
866: }
870: /*@
871: DMMeshGetChart - Return the interval for all mesh points [pStart, pEnd)
873: Not collective
875: Input Parameter:
876: . mesh - The DMMesh
878: Output Parameters:
879: + pStart - The first mesh point
880: - pEnd - The upper bound for mesh points
882: Level: beginner
884: .seealso: DMMeshCreate(), DMMeshSetChart()
885: @*/
886: PetscErrorCode DMMeshGetChart(DM dm, PetscInt *pStart, PetscInt *pEnd)
887: {
888: DM_Mesh *mesh = (DM_Mesh*) dm->data;
893: if (mesh->useNewImpl) {
894: PetscSectionGetChart(mesh->coneSection, pStart, pEnd);
895: } else {
896: Obj<PETSC_MESH_TYPE> m;
897: DMMeshGetMesh(dm, m);
898: if (pStart) *pStart = m->getSieve()->getChart().min();
899: if (pEnd) *pEnd = m->getSieve()->getChart().max();
900: }
901: return(0);
902: }
906: /*@
907: DMMeshSetChart - Set the interval for all mesh points [pStart, pEnd)
909: Not collective
911: Input Parameters:
912: + mesh - The DMMesh
913: . pStart - The first mesh point
914: - pEnd - The upper bound for mesh points
916: Output Parameters:
918: Level: beginner
920: .seealso: DMMeshCreate(), DMMeshGetChart()
921: @*/
922: PetscErrorCode DMMeshSetChart(DM dm, PetscInt pStart, PetscInt pEnd)
923: {
924: DM_Mesh *mesh = (DM_Mesh*) dm->data;
929: if (mesh->useNewImpl) {
930: PetscSectionSetChart(mesh->coneSection, pStart, pEnd);
931: } else {
932: Obj<PETSC_MESH_TYPE> m;
933: DMMeshGetMesh(dm, m);
934: m->getSieve()->setChart(PETSC_MESH_TYPE::sieve_type::chart_type(pStart, pEnd));
935: }
936: return(0);
937: }
941: /*@
942: DMMeshGetConeSize - Return the number of in-edges for this point in the Sieve DAG
944: Not collective
946: Input Parameters:
947: + mesh - The DMMesh
948: - p - The Sieve point, which must lie in the chart set with DMMeshSetChart()
950: Output Parameter:
951: . size - The cone size for point p
953: Level: beginner
955: .seealso: DMMeshCreate(), DMMeshSetConeSize(), DMMeshSetChart()
956: @*/
957: PetscErrorCode DMMeshGetConeSize(DM dm, PetscInt p, PetscInt *size)
958: {
959: DM_Mesh *mesh = (DM_Mesh*) dm->data;
965: if (mesh->useNewImpl) {
966: PetscSectionGetDof(mesh->coneSection, p, size);
967: } else {
968: Obj<PETSC_MESH_TYPE> m;
969: DMMeshGetMesh(dm, m);
970: *size = m->getSieve()->getConeSize(p);
971: }
972: return(0);
973: }
977: /*@
978: DMMeshSetConeSize - Set the number of in-edges for this point in the Sieve DAG
980: Not collective
982: Input Parameters:
983: + mesh - The DMMesh
984: . p - The Sieve point, which must lie in the chart set with DMMeshSetChart()
985: - size - The cone size for point p
987: Output Parameter:
989: Note:
990: This should be called after DMMeshSetChart().
992: Level: beginner
994: .seealso: DMMeshCreate(), DMMeshGetConeSize(), DMMeshSetChart()
995: @*/
996: PetscErrorCode DMMeshSetConeSize(DM dm, PetscInt p, PetscInt size)
997: {
998: DM_Mesh *mesh = (DM_Mesh*) dm->data;
1003: if (mesh->useNewImpl) {
1004: PetscSectionSetDof(mesh->coneSection, p, size);
1005: mesh->maxConeSize = PetscMax(mesh->maxConeSize, size);
1006: } else {
1007: Obj<PETSC_MESH_TYPE> m;
1008: DMMeshGetMesh(dm, m);
1009: m->getSieve()->setConeSize(p, size);
1010: }
1011: return(0);
1012: }
1016: /*@C
1017: DMMeshGetCone - Return the points on the in-edges for this point in the Sieve DAG
1019: Not collective
1021: Input Parameters:
1022: + mesh - The DMMesh
1023: - p - The Sieve point, which must lie in the chart set with DMMeshSetChart()
1025: Output Parameter:
1026: . cone - An array of points which are on the in-edges for point p
1028: Level: beginner
1030: Note:
1031: This routine is not available in Fortran.
1033: .seealso: DMMeshCreate(), DMMeshSetCone(), DMMeshSetChart()
1034: @*/
1035: PetscErrorCode DMMeshGetCone(DM dm, PetscInt p, const PetscInt *cone[])
1036: {
1037: DM_Mesh *mesh = (DM_Mesh*) dm->data;
1043: if (mesh->useNewImpl) {
1044: PetscInt off;
1046: PetscSectionGetOffset(mesh->coneSection, p, &off);
1047: *cone = &mesh->cones[off];
1048: } else {
1049: Obj<PETSC_MESH_TYPE> m;
1050: DMMeshGetMesh(dm, m);
1051: ALE::ISieveVisitor::PointRetriever<PETSC_MESH_TYPE::sieve_type> v(m->getSieve()->getConeSize(p));
1053: m->getSieve()->cone(p, v);
1054: if (!mesh->meetTmpA) {PetscMalloc2(m->getSieve()->getMaxConeSize(),PetscInt,&mesh->meetTmpA,m->getSieve()->getMaxConeSize(),PetscInt,&mesh->meetTmpB);}
1055: for (size_t c = 0; c < v.getSize(); ++c) {
1056: mesh->meetTmpA[c] = v.getPoints()[c];
1057: }
1058: *cone = mesh->meetTmpA;
1059: }
1060: return(0);
1061: }
1065: /*@
1066: DMMeshSetCone - Set the points on the in-edges for this point in the Sieve DAG
1068: Not collective
1070: Input Parameters:
1071: + mesh - The DMMesh
1072: . p - The Sieve point, which must lie in the chart set with DMMeshSetChart()
1073: - cone - An array of points which are on the in-edges for point p
1075: Output Parameter:
1077: Note:
1078: This should be called after all calls to DMMeshSetConeSize() and DMMeshSetUp().
1080: Level: beginner
1082: .seealso: DMMeshCreate(), DMMeshGetCone(), DMMeshSetChart(), DMMeshSetConeSize(), DMMeshSetUp()
1083: @*/
1084: PetscErrorCode DMMeshSetCone(DM dm, PetscInt p, const PetscInt cone[])
1085: {
1086: DM_Mesh *mesh = (DM_Mesh*) dm->data;
1092: if (mesh->useNewImpl) {
1093: PetscInt pStart, pEnd;
1094: PetscInt dof, off, c;
1096: PetscSectionGetChart(mesh->coneSection, &pStart, &pEnd);
1097: PetscSectionGetDof(mesh->coneSection, p, &dof);
1098: PetscSectionGetOffset(mesh->coneSection, p, &off);
1099: for (c = 0; c < dof; ++c) {
1100: if ((cone[c] < pStart) || (cone[c] >= pEnd)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cone point %d is not in the valid range [%d. %d)", cone[c], pStart, pEnd);
1101: mesh->cones[off+c] = cone[c];
1102: }
1103: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "This method does not make sense for the C++ Sieve implementation");
1104: return(0);
1105: }
1109: /*@
1110: DMMeshGetSupportSize - Return the number of out-edges for this point in the Sieve DAG
1112: Not collective
1114: Input Parameters:
1115: + mesh - The DMMesh
1116: - p - The Sieve point, which must lie in the chart set with DMMeshSetChart()
1118: Output Parameter:
1119: . size - The support size for point p
1121: Level: beginner
1123: .seealso: DMMeshCreate(), DMMeshSetConeSize(), DMMeshSetChart(), DMMeshGetConeSize()
1124: @*/
1125: PetscErrorCode DMMeshGetSupportSize(DM dm, PetscInt p, PetscInt *size)
1126: {
1127: DM_Mesh *mesh = (DM_Mesh*) dm->data;
1133: if (mesh->useNewImpl) {
1134: PetscSectionGetDof(mesh->supportSection, p, size);
1135: } else {
1136: Obj<PETSC_MESH_TYPE> m;
1137: DMMeshGetMesh(dm, m);
1138: *size = m->getSieve()->getSupportSize(p);
1139: }
1140: return(0);
1141: }
1145: /*@C
1146: DMMeshGetSupport - Return the points on the out-edges for this point in the Sieve DAG
1148: Not collective
1150: Input Parameters:
1151: + mesh - The DMMesh
1152: - p - The Sieve point, which must lie in the chart set with DMMeshSetChart()
1154: Output Parameter:
1155: . support - An array of points which are on the out-edges for point p
1157: Level: beginner
1159: .seealso: DMMeshCreate(), DMMeshSetCone(), DMMeshSetChart(), DMMeshGetCone()
1160: @*/
1161: PetscErrorCode DMMeshGetSupport(DM dm, PetscInt p, const PetscInt *support[])
1162: {
1163: DM_Mesh *mesh = (DM_Mesh*) dm->data;
1169: if (mesh->useNewImpl) {
1170: PetscInt off;
1172: PetscSectionGetOffset(mesh->supportSection, p, &off);
1173: *support = &mesh->supports[off];
1174: } else {
1175: Obj<PETSC_MESH_TYPE> m;
1176: DMMeshGetMesh(dm, m);
1177: ALE::ISieveVisitor::PointRetriever<PETSC_MESH_TYPE::sieve_type> v(m->getSieve()->getSupportSize(p));
1179: m->getSieve()->support(p, v);
1180: if (!mesh->joinTmpA) {PetscMalloc2(m->getSieve()->getMaxSupportSize(),PetscInt,&mesh->joinTmpA,m->getSieve()->getMaxSupportSize(),PetscInt,&mesh->joinTmpB);}
1181: for (size_t s = 0; s < v.getSize(); ++s) {
1182: mesh->joinTmpA[s] = v.getPoints()[s];
1183: }
1184: *support = mesh->joinTmpA;
1185: }
1186: return(0);
1187: }
1191: /*@C
1192: DMMeshGetTransitiveClosure - Return the points on the transitive closure of the in-edges or out-edges for this point in the Sieve DAG
1194: Not collective
1196: Input Parameters:
1197: + mesh - The DMMesh
1198: . p - The Sieve point, which must lie in the chart set with DMMeshSetChart()
1199: - useCone - PETSC_TRUE for in-edges, otherwise use out-edges
1201: Output Parameters:
1202: + numPoints - The number of points in the closure
1203: - points - The points
1205: Level: beginner
1207: .seealso: DMMeshCreate(), DMMeshSetCone(), DMMeshSetChart(), DMMeshGetCone()
1208: @*/
1209: PetscErrorCode DMMeshGetTransitiveClosure(DM dm, PetscInt p, PetscBool useCone, PetscInt *numPoints, const PetscInt *points[])
1210: {
1211: DM_Mesh *mesh = (DM_Mesh*) dm->data;
1216: if (!mesh->closureTmpA) {
1217: PetscInt maxSize;
1218: if (mesh->useNewImpl) {
1219: maxSize = PetscMax(mesh->maxConeSize, mesh->maxSupportSize)+1;
1220: } else {
1221: Obj<PETSC_MESH_TYPE> m;
1222: DMMeshGetMesh(dm, m);
1223: maxSize = PetscMax(m->getSieve()->getMaxConeSize(), m->getSieve()->getMaxSupportSize())+1;
1224: }
1225: PetscMalloc2(maxSize,PetscInt,&mesh->closureTmpA,maxSize,PetscInt,&mesh->closureTmpB);
1226: }
1227: if (mesh->useNewImpl) {
1228: const PetscInt *tmp;
1229: PetscInt tmpSize, t;
1230: PetscInt closureSize = 1;
1232: mesh->closureTmpA[0] = p;
1233: /* This is only 1-level */
1234: if (useCone) {
1235: DMMeshGetConeSize(dm, p, &tmpSize);
1236: DMMeshGetCone(dm, p, &tmp);
1237: } else {
1238: DMMeshGetSupportSize(dm, p, &tmpSize);
1239: DMMeshGetSupport(dm, p, &tmp);
1240: }
1241: for (t = 0; t < tmpSize; ++t) {
1242: mesh->closureTmpA[closureSize++] = tmp[t];
1243: }
1244: if (numPoints) *numPoints = closureSize;
1245: if (points) *points = mesh->closureTmpA;
1246: } else {
1247: Obj<PETSC_MESH_TYPE> m;
1248: DMMeshGetMesh(dm, m);
1249: typedef ALE::ISieveVisitor::TransitiveClosureVisitor<PETSC_MESH_TYPE::sieve_type> visitor_type;
1250: visitor_type::visitor_type nV;
1251: visitor_type cV(*m->getSieve(), nV);
1253: if (useCone) {
1254: m->getSieve()->cone(p, cV);
1255: } else {
1256: cV.setIsCone(false);
1257: m->getSieve()->support(p, cV);
1258: }
1259: int i = 0;
1261: for (std::set<PETSC_MESH_TYPE::point_type>::const_iterator p_iter = cV.getPoints().begin(); p_iter != cV.getPoints().end(); ++p_iter, ++i) {
1262: mesh->closureTmpA[i] = *p_iter;
1263: }
1264: if (numPoints) *numPoints = cV.getPoints().size();
1265: if (points) *points = mesh->closureTmpA;
1266: }
1267: return(0);
1268: }
1272: /*@
1273: DMMeshGetMaxSizes - Return the maximum number of in-edges (cone) and out-edges (support) for any point in the Sieve DAG
1275: Not collective
1277: Input Parameter:
1278: . mesh - The DMMesh
1280: Output Parameters:
1281: + maxConeSize - The maximum number of in-edges
1282: - maxSupportSize - The maximum number of out-edges
1284: Level: beginner
1286: .seealso: DMMeshCreate(), DMMeshSetConeSize(), DMMeshSetChart()
1287: @*/
1288: PetscErrorCode DMMeshGetMaxSizes(DM dm, PetscInt *maxConeSize, PetscInt *maxSupportSize)
1289: {
1290: DM_Mesh *mesh = (DM_Mesh*) dm->data;
1295: if (mesh->useNewImpl) {
1296: if (maxConeSize) *maxConeSize = mesh->maxConeSize;
1297: if (maxSupportSize) *maxSupportSize = mesh->maxSupportSize;
1298: } else {
1299: Obj<PETSC_MESH_TYPE> m;
1300: DMMeshGetMesh(dm, m);
1301: if (maxConeSize) *maxConeSize = m->getSieve()->getMaxConeSize();
1302: if (maxSupportSize) *maxSupportSize = m->getSieve()->getMaxSupportSize();
1303: }
1304: return(0);
1305: }
1309: /*@
1310: DMMeshSetUp - Allocate space for the Sieve DAG
1312: Not collective
1314: Input Parameter:
1315: . mesh - The DMMesh
1317: Output Parameter:
1319: Note:
1320: This should be called after DMMeshSetChart() and all calls to DMMeshSetConeSize()
1322: Level: beginner
1324: .seealso: DMMeshCreate(), DMMeshSetChart(), DMMeshSetConeSize()
1325: @*/
1326: PetscErrorCode DMMeshSetUp(DM dm)
1327: {
1328: DM_Mesh *mesh = (DM_Mesh*) dm->data;
1333: if (mesh->useNewImpl) {
1334: PetscInt size;
1336: PetscSectionSetUp(mesh->coneSection);
1337: PetscSectionGetStorageSize(mesh->coneSection, &size);
1338: PetscMalloc(size * sizeof(PetscInt), &mesh->cones);
1339: }
1340: return(0);
1341: }
1345: /*@
1346: DMMeshSymmetrize - Creates support (out-edge) information from cone (in-edge) inoformation
1348: Not collective
1350: Input Parameter:
1351: . mesh - The DMMesh
1353: Output Parameter:
1355: Note:
1356: This should be called after all calls to DMMeshSetCone()
1358: Level: beginner
1360: .seealso: DMMeshCreate(), DMMeshSetChart(), DMMeshSetConeSize(), DMMeshSetCone()
1361: @*/
1362: PetscErrorCode DMMeshSymmetrize(DM dm)
1363: {
1364: DM_Mesh *mesh = (DM_Mesh*) dm->data;
1369: if (mesh->useNewImpl) {
1370: PetscInt *offsets;
1371: PetscInt supportSize;
1372: PetscInt pStart, pEnd, p;
1374: /* Calculate support sizes */
1375: DMMeshGetChart(dm, &pStart, &pEnd);
1376: PetscSectionSetChart(mesh->supportSection, pStart, pEnd);
1377: for (p = pStart; p < pEnd; ++p) {
1378: PetscInt dof, off, c;
1380: PetscSectionGetDof(mesh->coneSection, p, &dof);
1381: PetscSectionGetOffset(mesh->coneSection, p, &off);
1382: for (c = off; c < off+dof; ++c) {
1383: PetscSectionAddDof(mesh->supportSection, mesh->cones[c], 1);
1384: }
1385: }
1386: for (p = pStart; p < pEnd; ++p) {
1387: PetscInt dof;
1389: PetscSectionGetDof(mesh->supportSection, p, &dof);
1391: mesh->maxSupportSize = PetscMax(mesh->maxSupportSize, dof);
1392: }
1393: PetscSectionSetUp(mesh->supportSection);
1394: /* Calculate supports */
1395: PetscSectionGetStorageSize(mesh->supportSection, &supportSize);
1396: PetscMalloc(supportSize * sizeof(PetscInt), &mesh->supports);
1397: PetscMalloc((pEnd - pStart) * sizeof(PetscInt), &offsets);
1398: PetscMemzero(offsets, (pEnd - pStart) * sizeof(PetscInt));
1399: for (p = pStart; p < pEnd; ++p) {
1400: PetscInt dof, off, c;
1402: PetscSectionGetDof(mesh->coneSection, p, &dof);
1403: PetscSectionGetOffset(mesh->coneSection, p, &off);
1404: for (c = off; c < off+dof; ++c) {
1405: const PetscInt q = mesh->cones[c];
1406: PetscInt offS;
1408: PetscSectionGetOffset(mesh->supportSection, q, &offS);
1410: mesh->supports[offS+offsets[q]] = p;
1411: ++offsets[q];
1412: }
1413: }
1414: PetscFree(offsets);
1415: } else {
1416: Obj<PETSC_MESH_TYPE> m;
1417: DMMeshGetMesh(dm, m);
1418: m->getSieve()->symmetrize();
1419: }
1420: return(0);
1421: }
1425: /*@
1426: DMMeshStratify - The Sieve DAG for most topologies is a graded poset (http://en.wikipedia.org/wiki/Graded_poset), and
1427: can be illustrated by Hasse Diagram (a http://en.wikipedia.org/wiki/Hasse_diagram). The strata group all points of the
1428: same grade, and this function calculates the strata. This grade can be seen as the height (or depth) of the point in
1429: the DAG.
1431: Not collective
1433: Input Parameter:
1434: . mesh - The DMMesh
1436: Output Parameter:
1438: Notes:
1439: The normal association for the point grade is element dimension (or co-dimension). For instance, all vertices would
1440: have depth 0, and all edges depth 1. Likewise, all cells heights would have height 0, and all faces height 1.
1442: This should be called after all calls to DMMeshSymmetrize()
1444: Level: beginner
1446: .seealso: DMMeshCreate(), DMMeshSymmetrize()
1447: @*/
1448: PetscErrorCode DMMeshStratify(DM dm)
1449: {
1450: DM_Mesh *mesh = (DM_Mesh*) dm->data;
1455: if (mesh->useNewImpl) {
1456: PetscInt pStart, pEnd, p;
1457: PetscInt numRoots = 0, numLeaves = 0;
1459: /* Calculate depth */
1460: PetscSectionGetChart(mesh->coneSection, &pStart, &pEnd);
1461: /* Initialize roots and count leaves */
1462: for (p = pStart; p < pEnd; ++p) {
1463: PetscInt coneSize, supportSize;
1465: PetscSectionGetDof(mesh->coneSection, p, &coneSize);
1466: PetscSectionGetDof(mesh->supportSection, p, &supportSize);
1467: if (!coneSize && supportSize) {
1468: ++numRoots;
1469: DMMeshSetLabelValue(dm, "depth", p, 0);
1470: } else if (!supportSize && coneSize) {
1471: ++numLeaves;
1472: }
1473: }
1474: if (numRoots + numLeaves == (pEnd - pStart)) {
1475: for (p = pStart; p < pEnd; ++p) {
1476: PetscInt coneSize, supportSize;
1478: PetscSectionGetDof(mesh->coneSection, p, &coneSize);
1479: PetscSectionGetDof(mesh->supportSection, p, &supportSize);
1480: if (!supportSize && coneSize) {
1481: DMMeshSetLabelValue(dm, "depth", p, 1);
1482: }
1483: }
1484: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Have not yet coded general stratification");
1485: } else {
1486: Obj<PETSC_MESH_TYPE> m;
1487: DMMeshGetMesh(dm, m);
1488: m->stratify();
1489: }
1490: return(0);
1491: }
1495: /*@C
1496: DMMeshGetLabelValue - Get the value in a Sieve Label for the given point, with 0 as the default
1498: Not Collective
1500: Input Parameters:
1501: + dm - The DMMesh object
1502: . name - The label name
1503: - point - The mesh point
1505: Output Parameter:
1506: . value - The label value for this point, or 0 if the point is not in the label
1508: Level: beginner
1510: .keywords: mesh
1511: .seealso: DMMeshSetLabelValue(), DMMeshGetLabelStratum()
1512: @*/
1513: PetscErrorCode DMMeshGetLabelValue(DM dm, const char name[], PetscInt point, PetscInt *value)
1514: {
1515: DM_Mesh *mesh = (DM_Mesh*) dm->data;
1521: if (mesh->useNewImpl) {
1522: SieveLabel next = mesh->labels;
1523: PetscBool flg = PETSC_FALSE;
1524: PetscInt v, p;
1526: *value = 0;
1527: while (next) {
1528: PetscStrcmp(name, next->name, &flg);
1529: if (flg) break;
1530: next = next->next;
1531: }
1532: if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "No label named %s was found", name);
1533: /* Find, or add, label value */
1534: for (v = 0; v < next->numStrata; ++v) {
1535: for (p = next->stratumOffsets[v]; p < next->stratumOffsets[v]+next->stratumSizes[v]; ++p) {
1536: if (next->points[p] == point) {
1537: *value = next->stratumValues[v];
1538: break;
1539: }
1540: }
1541: }
1542: } else {
1543: ALE::Obj<PETSC_MESH_TYPE> m;
1544: DMMeshGetMesh(dm, m);
1545: *value = m->getValue(m->getLabel(name), point);
1546: }
1547: return(0);
1548: }
1552: /*@C
1553: DMMeshSetLabelValue - Add a point to a Sieve Label with given value
1555: Not Collective
1557: Input Parameters:
1558: + dm - The DMMesh object
1559: . name - The label name
1560: . point - The mesh point
1561: - value - The label value for this point
1563: Output Parameter:
1565: Level: beginner
1567: .keywords: mesh
1568: .seealso: DMMeshGetLabelStratum()
1569: @*/
1570: PetscErrorCode DMMeshSetLabelValue(DM dm, const char name[], PetscInt point, PetscInt value)
1571: {
1572: DM_Mesh *mesh = (DM_Mesh*) dm->data;
1578: if (mesh->useNewImpl) {
1579: SieveLabel next = mesh->labels;
1580: PetscBool flg = PETSC_FALSE;
1581: PetscInt v, p;
1583: /* Find, or create, label */
1584: while (next) {
1585: PetscStrcmp(name, next->name, &flg);
1586: if (flg) break;
1587: next = next->next;
1588: }
1589: if (!flg) {
1590: SieveLabel tmpLabel = mesh->labels;
1591: PetscNew(struct Sieve_Label, &mesh->labels);
1592: mesh->labels->next = tmpLabel;
1593: next = mesh->labels;
1594: PetscStrallocpy(name, &next->name);
1595: }
1596: /* Find, or add, label value */
1597: for (v = 0; v < next->numStrata; ++v) {
1598: if (next->stratumValues[v] == value) break;
1599: }
1600: if (v >= next->numStrata) {
1601: PetscInt *tmpV, *tmpO, *tmpS;
1602: PetscMalloc3(next->numStrata+1,PetscInt,&tmpV,next->numStrata+2,PetscInt,&tmpO,next->numStrata+1,PetscInt,&tmpS);
1603: for (v = 0; v < next->numStrata; ++v) {
1604: tmpV[v] = next->stratumValues[v];
1605: tmpO[v] = next->stratumOffsets[v];
1606: tmpS[v] = next->stratumSizes[v];
1607: }
1608: tmpV[v] = value;
1609: tmpO[v] = v == 0 ? 0 : next->stratumOffsets[v];
1610: tmpS[v] = 0;
1611: tmpO[v+1] = tmpO[v];
1612: ++next->numStrata;
1613: PetscFree3(next->stratumValues,next->stratumOffsets,next->stratumSizes);
1615: next->stratumValues = tmpV;
1616: next->stratumOffsets = tmpO;
1617: next->stratumSizes = tmpS;
1618: }
1619: /* Check whether point exists */
1620: for (p = next->stratumOffsets[v]; p < next->stratumOffsets[v]+next->stratumSizes[v]; ++p) {
1621: if (next->points[p] == point) break;
1622: }
1623: /* Add point: NEED TO OPTIMIZE */
1624: if (p >= next->stratumOffsets[v]+next->stratumSizes[v]) {
1625: /* Check for reallocation */
1626: if (next->stratumSizes[v] >= next->stratumOffsets[v+1]-next->stratumOffsets[v]) {
1627: PetscInt oldSize = next->stratumOffsets[v+1]-next->stratumOffsets[v];
1628: PetscInt newSize = PetscMax(10, 2*oldSize); /* Double the size, since 2 is the optimal base for this online algorithm */
1629: PetscInt shift = newSize - oldSize;
1630: PetscInt allocSize = next->stratumOffsets[next->numStrata] + shift;
1631: PetscInt *newPoints;
1632: PetscInt w, q;
1634: PetscMalloc(allocSize * sizeof(PetscInt), &newPoints);
1635: for (q = 0; q < next->stratumOffsets[v]+next->stratumSizes[v]; ++q) {
1636: newPoints[q] = next->points[q];
1637: }
1638: for (w = v+1; w < next->numStrata; ++w) {
1639: for (q = next->stratumOffsets[w]; q < next->stratumOffsets[w]+next->stratumSizes[w]; ++q) {
1640: newPoints[q+shift] = next->points[q];
1641: }
1642: next->stratumOffsets[w] += shift;
1643: }
1644: next->stratumOffsets[next->numStrata] += shift;
1646: PetscFree(next->points);
1648: next->points = newPoints;
1649: }
1650: /* Insert point and resort */
1651: next->points[next->stratumOffsets[v]+next->stratumSizes[v]] = point;
1652: ++next->stratumSizes[v];
1653: PetscSortInt(next->stratumSizes[v], &next->points[next->stratumOffsets[v]]);
1654: }
1655: } else {
1656: ALE::Obj<PETSC_MESH_TYPE> m;
1657: DMMeshGetMesh(dm, m);
1658: m->setValue(m->getLabel(name), point, value);
1659: }
1660: return(0);
1661: }
1665: /*@C
1666: DMMeshGetLabelSize - Get the number of different integer ids in a Label
1668: Not Collective
1670: Input Parameters:
1671: + dm - The DMMesh object
1672: - name - The label name
1674: Output Parameter:
1675: . size - The label size (number of different integer ids)
1677: Level: beginner
1679: .keywords: mesh
1680: .seealso: DMMeshSetLabelValue()
1681: @*/
1682: PetscErrorCode DMMeshGetLabelSize(DM dm, const char name[], PetscInt *size)
1683: {
1684: DM_Mesh *mesh = (DM_Mesh*) dm->data;
1691: if (mesh->useNewImpl) {
1692: SieveLabel next = mesh->labels;
1693: PetscBool flg;
1695: *size = 0;
1696: while (next) {
1697: PetscStrcmp(name, next->name, &flg);
1698: if (flg) {
1699: *size = next->numStrata;
1700: break;
1701: }
1702: next = next->next;
1703: }
1704: } else {
1705: ALE::Obj<PETSC_MESH_TYPE> m;
1706: DMMeshGetMesh(dm, m);
1707: *size = m->getLabel(name)->getCapSize();
1708: }
1709: return(0);
1710: }
1714: /*@C
1715: DMMeshGetLabelIdIS - Get the integer ids in a label
1717: Not Collective
1719: Input Parameters:
1720: + mesh - The DMMesh object
1721: - name - The label name
1723: Output Parameter:
1724: . ids - The integer ids
1726: Level: beginner
1728: .keywords: mesh
1729: .seealso: DMMeshGetLabelSize()
1730: @*/
1731: PetscErrorCode DMMeshGetLabelIdIS(DM dm, const char name[], IS *ids)
1732: {
1733: DM_Mesh *mesh = (DM_Mesh*) dm->data;
1734: PetscInt *values;
1735: PetscInt size, i = 0;
1742: if (mesh->useNewImpl) {
1743: SieveLabel next = mesh->labels;
1744: PetscBool flg;
1746: while (next) {
1747: PetscStrcmp(name, next->name, &flg);
1748: if (flg) {
1749: size = next->numStrata;
1750: PetscMalloc(size * sizeof(PetscInt), &values);
1751: for (i = 0; i < next->numStrata; ++i) {
1752: values[i] = next->stratumValues[i];
1753: }
1754: break;
1755: }
1756: next = next->next;
1757: }
1758: } else {
1759: ALE::Obj<PETSC_MESH_TYPE> m;
1760: DMMeshGetMesh(dm, m);
1761: const ALE::Obj<PETSC_MESH_TYPE::label_type::capSequence>& labelIds = m->getLabel(name)->cap();
1762: const PETSC_MESH_TYPE::label_type::capSequence::const_iterator iEnd = labelIds->end();
1764: size = labelIds->size();
1765: PetscMalloc(size * sizeof(PetscInt), &values);
1766: for (PETSC_MESH_TYPE::label_type::capSequence::const_iterator i_iter = labelIds->begin(); i_iter != iEnd; ++i_iter, ++i) {
1767: values[i] = *i_iter;
1768: }
1769: }
1770: ISCreateGeneral(PetscObjectComm((PetscObject)dm), size, values, PETSC_OWN_POINTER, ids);
1771: return(0);
1772: }
1776: /*@C
1777: DMMeshGetStratumSize - Get the number of points in a label stratum
1779: Not Collective
1781: Input Parameters:
1782: + dm - The DMMesh object
1783: . name - The label name
1784: - value - The stratum value
1786: Output Parameter:
1787: . size - The stratum size
1789: Level: beginner
1791: .keywords: mesh
1792: .seealso: DMMeshGetLabelSize(), DMMeshGetLabelIds()
1793: @*/
1794: PetscErrorCode DMMeshGetStratumSize(DM dm, const char name[], PetscInt value, PetscInt *size)
1795: {
1796: DM_Mesh *mesh = (DM_Mesh*) dm->data;
1803: if (mesh->useNewImpl) {
1804: SieveLabel next = mesh->labels;
1805: PetscBool flg;
1807: *size = 0;
1808: while (next) {
1809: PetscStrcmp(name, next->name, &flg);
1810: if (flg) {
1811: PetscInt v;
1813: for (v = 0; v < next->numStrata; ++v) {
1814: if (next->stratumValues[v] == value) {
1815: *size = next->stratumSizes[v];
1816: break;
1817: }
1818: }
1819: break;
1820: }
1821: next = next->next;
1822: }
1823: } else {
1824: ALE::Obj<PETSC_MESH_TYPE> m;
1825: DMMeshGetMesh(dm, m);
1826: *size = m->getLabelStratum(name, value)->size();
1827: }
1828: return(0);
1829: }
1833: /*@C
1834: DMMeshGetStratumIS - Get the points in a label stratum
1836: Not Collective
1838: Input Parameters:
1839: + dm - The DMMesh object
1840: . name - The label name
1841: - value - The stratum value
1843: Output Parameter:
1844: . is - The stratum points
1846: Level: beginner
1848: .keywords: mesh
1849: .seealso: DMMeshGetStratumSize()
1850: @*/
1851: PetscErrorCode DMMeshGetStratumIS(DM dm, const char name[], PetscInt value, IS *is)
1852: {
1853: DM_Mesh *mesh = (DM_Mesh*) dm->data;
1860: *is = NULL;
1861: if (mesh->useNewImpl) {
1862: SieveLabel next = mesh->labels;
1863: PetscBool flg;
1865: while (next) {
1866: PetscStrcmp(name, next->name, &flg);
1867: if (flg) {
1868: PetscInt v;
1870: for (v = 0; v < next->numStrata; ++v) {
1871: if (next->stratumValues[v] == value) {
1872: ISCreateGeneral(PETSC_COMM_SELF, next->stratumSizes[v], &next->points[next->stratumOffsets[v]], PETSC_COPY_VALUES, is);
1873: break;
1874: }
1875: }
1876: break;
1877: }
1878: next = next->next;
1879: }
1880: } else {
1881: ALE::Obj<PETSC_MESH_TYPE> mesh;
1882: DMMeshGetMesh(dm, mesh);
1883: if (mesh->hasLabel(name)) {
1884: const Obj<PETSC_MESH_TYPE::label_sequence>& stratum = mesh->getLabelStratum(name, value);
1885: PetscInt *idx, i = 0;
1887: PetscMalloc(stratum->size() * sizeof(PetscInt), &idx);
1888: for (PETSC_MESH_TYPE::label_sequence::iterator e_iter = stratum->begin(); e_iter != stratum->end(); ++e_iter, ++i) {
1889: idx[i] = *e_iter;
1890: }
1891: ISCreateGeneral(PETSC_COMM_SELF, stratum->size(), idx, PETSC_OWN_POINTER, is);
1892: }
1893: }
1894: return(0);
1895: }
1899: /* This is a 1-level join */
1900: PetscErrorCode DMMeshJoinPoints(DM dm, const PetscInt points[], PetscInt *coveredPoint)
1901: {
1902: DM_Mesh *mesh = (DM_Mesh*) dm->data;
1909: if (mesh->useNewImpl) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not yet supported");
1910: else {
1911: ALE::Obj<PETSC_MESH_TYPE> m;
1912: DMMeshGetMesh(dm, m);
1913: /* const Obj<typename Mesh::sieve_type::supportSet> edge = m->getSieve()->nJoin(points[0], points[1], 1); */
1914: *coveredPoint = -1;
1915: }
1916: return(0);
1917: }
1921: /* This is a 1-level meet */
1922: PetscErrorCode DMMeshMeetPoints(DM dm, PetscInt numPoints, const PetscInt points[], PetscInt *numCoveringPoints, const PetscInt **coveringPoints)
1923: {
1924: DM_Mesh *mesh = (DM_Mesh*) dm->data;
1925: PetscInt *meet[2];
1926: PetscInt meetSize, i = 0;
1927: PetscInt dof, off, p, c, m;
1935: if (mesh->useNewImpl) {
1936: if (!mesh->meetTmpA) {PetscMalloc2(mesh->maxConeSize,PetscInt,&mesh->meetTmpA,mesh->maxConeSize,PetscInt,&mesh->meetTmpB);}
1937: meet[0] = mesh->meetTmpA; meet[1] = mesh->meetTmpB;
1938: /* Copy in cone of first point */
1939: PetscSectionGetDof(mesh->coneSection, points[0], &dof);
1940: PetscSectionGetOffset(mesh->coneSection, points[0], &off);
1941: for (meetSize = 0; meetSize < dof; ++meetSize) {
1942: meet[i][meetSize] = mesh->cones[off+meetSize];
1943: }
1944: /* Check each successive cone */
1945: for (p = 1; p < numPoints; ++p) {
1946: PetscInt newMeetSize = 0;
1948: PetscSectionGetDof(mesh->coneSection, points[p], &dof);
1949: PetscSectionGetOffset(mesh->coneSection, points[p], &off);
1950: for (c = 0; c < dof; ++c) {
1951: const PetscInt point = mesh->cones[off+c];
1953: for (m = 0; m < meetSize; ++m) {
1954: if (point == meet[i][m]) {
1955: meet[1-i][newMeetSize++] = point;
1956: break;
1957: }
1958: }
1959: }
1960: meetSize = newMeetSize;
1961: i = 1-i;
1962: }
1963: *numCoveringPoints = meetSize;
1964: *coveringPoints = meet[1-i];
1965: } else {
1966: ALE::Obj<PETSC_MESH_TYPE> m;
1967: DMMeshGetMesh(dm, m);
1968: /* const Obj<typename Mesh::sieve_type::supportSet> edge = m->getSieve()->nJoin(points[0], points[1], 1); */
1969: *numCoveringPoints = 0;
1970: *coveringPoints = NULL;
1971: }
1972: return(0);
1973: }
1977: /*@C
1978: DMMeshGetMaximumDegree - Return the maximum degree of any mesh vertex
1980: Collective on mesh
1982: Input Parameter:
1983: . mesh - The DMMesh
1985: Output Parameter:
1986: . maxDegree - The maximum number of edges at any vertex
1988: Level: beginner
1990: .seealso: DMMeshCreate()
1991: @*/
1992: PetscErrorCode DMMeshGetMaximumDegree(DM dm, PetscInt *maxDegree)
1993: {
1994: Obj<PETSC_MESH_TYPE> m;
1995: PetscErrorCode ierr;
1998: DMMeshGetMesh(dm, m);
1999: const ALE::Obj<PETSC_MESH_TYPE::label_sequence>& vertices = m->depthStratum(0);
2000: const ALE::Obj<PETSC_MESH_TYPE::sieve_type>& sieve = m->getSieve();
2001: PetscInt maxDeg = -1;
2003: for (PETSC_MESH_TYPE::label_sequence::iterator v_iter = vertices->begin(); v_iter != vertices->end(); ++v_iter) {
2004: maxDeg = PetscMax(maxDeg, (PetscInt) sieve->getSupportSize(*v_iter));
2005: }
2006: *maxDegree = maxDeg;
2007: return(0);
2008: }
2010: extern PetscErrorCode assembleFullField(VecScatter, Vec, Vec, InsertMode);
2014: /*@
2015: DMMeshRestrictVector - Insert values from a global vector into a local ghosted vector
2017: Collective on g
2019: Input Parameters:
2020: + g - The global vector
2021: . l - The local vector
2022: - mode - either ADD_VALUES or INSERT_VALUES, where
2023: ADD_VALUES adds values to any existing entries, and
2024: INSERT_VALUES replaces existing entries with new values
2026: Level: beginner
2028: .seealso: MatSetOption()
2029: @*/
2030: PetscErrorCode DMMeshRestrictVector(Vec g, Vec l, InsertMode mode)
2031: {
2032: VecScatter injection;
2036: PetscLogEventBegin(DMMesh_restrictVector,0,0,0,0);
2037: PetscObjectQuery((PetscObject) g, "injection", (PetscObject*) &injection);
2038: if (injection) {
2039: VecScatterBegin(injection, g, l, mode, SCATTER_REVERSE);
2040: VecScatterEnd(injection, g, l, mode, SCATTER_REVERSE);
2041: } else {
2042: if (mode == INSERT_VALUES) {
2043: VecCopy(g, l);
2044: } else {
2045: VecAXPY(l, 1.0, g);
2046: }
2047: }
2048: PetscLogEventEnd(DMMesh_restrictVector,0,0,0,0);
2049: return(0);
2050: }
2054: /*@
2055: DMMeshAssembleVectorComplete - Insert values from a local ghosted vector into a global vector
2057: Collective on g
2059: Input Parameters:
2060: + g - The global vector
2061: . l - The local vector
2062: - mode - either ADD_VALUES or INSERT_VALUES, where
2063: ADD_VALUES adds values to any existing entries, and
2064: INSERT_VALUES replaces existing entries with new values
2066: Level: beginner
2068: .seealso: MatSetOption()
2069: @*/
2070: PetscErrorCode DMMeshAssembleVectorComplete(Vec g, Vec l, InsertMode mode)
2071: {
2072: VecScatter injection;
2076: PetscLogEventBegin(DMMesh_assembleVectorComplete,0,0,0,0);
2077: PetscObjectQuery((PetscObject) g, "injection", (PetscObject*) &injection);
2078: if (injection) {
2079: VecScatterBegin(injection, l, g, mode, SCATTER_FORWARD);
2080: VecScatterEnd(injection, l, g, mode, SCATTER_FORWARD);
2081: } else {
2082: if (mode == INSERT_VALUES) {
2083: VecCopy(l, g);
2084: } else {
2085: VecAXPY(g, 1.0, l);
2086: }
2087: }
2088: PetscLogEventEnd(DMMesh_assembleVectorComplete,0,0,0,0);
2089: return(0);
2090: }
2094: /*@
2095: DMMeshAssembleVector - Insert values into a vector
2097: Collective on A
2099: Input Parameters:
2100: + b - the vector
2101: . e - The element number
2102: . v - The values
2103: - mode - either ADD_VALUES or INSERT_VALUES, where
2104: ADD_VALUES adds values to any existing entries, and
2105: INSERT_VALUES replaces existing entries with new values
2107: Level: beginner
2109: .seealso: VecSetOption()
2110: @*/
2111: PetscErrorCode DMMeshAssembleVector(Vec b, PetscInt e, PetscScalar v[], InsertMode mode)
2112: {
2113: DM dm;
2114: SectionReal section;
2118: VecGetDM(b, &dm);
2119: DMMeshGetSectionReal(dm, "x", §ion);
2120: DMMeshAssembleVectorDM(b, dm, section, e, v, mode);
2121: SectionRealDestroy(§ion);
2122: return(0);
2123: }
2125: PetscErrorCode DMMeshAssembleVectorDM(Vec b, DM dm, SectionReal section, PetscInt e, PetscScalar v[], InsertMode mode)
2126: {
2127: ALE::Obj<PETSC_MESH_TYPE> m;
2128: ALE::Obj<PETSC_MESH_TYPE::real_section_type> s;
2129: PetscInt firstElement;
2130: PetscErrorCode ierr;
2133: PetscLogEventBegin(DMMesh_assembleVector,0,0,0,0);
2134: DMMeshGetMesh(dm, m);
2135: SectionRealGetSection(section, s);
2136: //firstElement = elementBundle->getLocalSizes()[bundle->getCommRank()];
2137: firstElement = 0;
2138: #if defined(PETSC_USE_COMPLEX)
2139: SETERRQ(PetscObjectComm((PetscObject)mesh),PETSC_ERR_SUP, "SectionReal does not support complex update");
2140: #else
2141: if (mode == INSERT_VALUES) {
2142: m->update(s, PETSC_MESH_TYPE::point_type(e + firstElement), v);
2143: } else {
2144: m->updateAdd(s, PETSC_MESH_TYPE::point_type(e + firstElement), v);
2145: }
2146: #endif
2147: PetscLogEventEnd(DMMesh_assembleVector,0,0,0,0);
2148: return(0);
2149: }
2153: /*@C
2154: MatSetValuesTopology - Sets values in a matrix using DM Mesh points rather than indices
2156: Not Collective
2158: Input Parameters:
2159: + mat - the matrix
2160: . dmr - The row DM
2161: . nrow, rowPoints - number of rows and their local Sieve points
2162: . dmc - The column DM
2163: . ncol, colPoints - number of columns and their local Sieve points
2164: . v - a logically two-dimensional array of values
2165: - mode - either ADD_VALUES or INSERT_VALUES, where
2166: ADD_VALUES adds values to any existing entries, and
2167: INSERT_VALUES replaces existing entries with new values
2169: Level: intermediate
2171: .seealso: DMMeshCreate(), MatSetValuesStencil()
2172: @*/
2173: PetscErrorCode MatSetValuesTopology(Mat mat, DM dmr, PetscInt nrow, const PetscInt rowPoints[], DM dmc, PetscInt ncol, const PetscInt colPoints[], const PetscScalar v[], InsertMode mode)
2174: {
2175: ALE::Obj<PETSC_MESH_TYPE> mr;
2176: ALE::Obj<PETSC_MESH_TYPE> mc;
2177: PetscErrorCode ierr;
2182: if (!nrow || !ncol) return(0); /* no values to insert */
2188: DMMeshGetMesh(dmr, mr);
2189: DMMeshGetMesh(dmc, mc);
2190: typedef ALE::ISieveVisitor::IndicesVisitor<PETSC_MESH_TYPE::real_section_type,PETSC_MESH_TYPE::order_type,PetscInt> visitor_type;
2191: visitor_type rV(*mr->getRealSection("default"), *mr->getFactory()->getLocalOrder(mr, "default", mr->getRealSection("default")),(int) pow((double) mr->getSieve()->getMaxConeSize(), mr->depth())*mr->getMaxDof() *nrow, mr->depth() > 1);
2192: visitor_type cV(*mc->getRealSection("default"), *mc->getFactory()->getLocalOrder(mc, "default", mc->getRealSection("default")),(int) pow((double) mc->getSieve()->getMaxConeSize(), mc->depth())*mc->getMaxDof() *ncol, mc->depth() > 1);
2194: try {
2195: for (PetscInt r = 0; r < nrow; ++r) {
2196: ALE::ISieveTraversal<PETSC_MESH_TYPE::sieve_type>::orientedClosure(*mr->getSieve(), rowPoints[r], rV);
2197: }
2198: } catch(ALE::Exception e) {
2199: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB, e.message());
2200: }
2201: const PetscInt *rowIndices = rV.getValues();
2202: const int numRowIndices = rV.getSize();
2203: try {
2204: for (PetscInt c = 0; c < ncol; ++c) {
2205: ALE::ISieveTraversal<PETSC_MESH_TYPE::sieve_type>::orientedClosure(*mc->getSieve(), colPoints[c], cV);
2206: }
2207: } catch(ALE::Exception e) {
2208: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB, e.message());
2209: }
2210: const PetscInt *colIndices = cV.getValues();
2211: const int numColIndices = cV.getSize();
2213: MatSetValuesLocal(mat, numRowIndices, rowIndices, numColIndices, colIndices, v, mode);
2214: return(0);
2215: }
2219: 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)
2220: {
2224: typedef ALE::ISieveVisitor::IndicesVisitor<PETSC_MESH_TYPE::real_section_type,PETSC_MESH_TYPE::order_type,PetscInt> visitor_type;
2225: visitor_type iV(*section, *globalOrder, (int) pow((double) m->getSieve()->getMaxConeSize(), m->depth())*m->getMaxDof(), m->depth() > 1);
2227: updateOperator(A, *m->getSieve(), iV, e, array, mode);
2228: return(0);
2229: }
2233: 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)
2234: {
2235: typedef ALE::ISieveVisitor::IndicesVisitor<PETSC_MESH_TYPE::real_section_type,PETSC_MESH_TYPE::order_type,PetscInt> visitor_type;
2236: visitor_type iVr(*rowSection, *rowGlobalOrder, (int) pow((double) rowM->getSieve()->getMaxConeSize(), rowM->depth())*rowM->getMaxDof(), rowM->depth() > 1);
2237: visitor_type iVc(*colSection, *colGlobalOrder, (int) pow((double) colM->getSieve()->getMaxConeSize(), colM->depth())*colM->getMaxDof(), colM->depth() > 1);
2239: PetscErrorCode updateOperator(A, *rowM->getSieve(), iVr, rowE, *colM->getSieve(), iVc, colE, array, mode);
2240: return(0);
2241: }
2245: /*@
2246: DMMeshSetMaxDof - Sets the maximum number of degrees of freedom on any sieve point
2248: Logically Collective on A
2250: Input Parameters:
2251: + A - the matrix
2252: . mesh - DMMesh needed for orderings
2253: . section - A Section which describes the layout
2254: . e - The element number
2255: . v - The values
2256: - mode - either ADD_VALUES or INSERT_VALUES, where
2257: ADD_VALUES adds values to any existing entries, and
2258: INSERT_VALUES replaces existing entries with new values
2260: Notes: This is used by routines like DMMeshUpdateOperator() to bound buffer sizes
2262: Level: developer
2264: .seealso: DMMeshUpdateOperator(), DMMeshAssembleMatrixDM()
2265: @*/
2266: PetscErrorCode DMMeshSetMaxDof(DM dm, PetscInt maxDof)
2267: {
2268: Obj<PETSC_MESH_TYPE> m;
2269: PetscErrorCode ierr;
2272: DMMeshGetMesh(dm, m);
2273: m->setMaxDof(maxDof);
2274: return(0);
2275: }
2279: /*@
2280: DMMeshAssembleMatrixDM - Insert values into a matrix
2282: Collective on A
2284: Input Parameters:
2285: + A - the matrix
2286: . dm - DMMesh needed for orderings
2287: . section - A Section which describes the layout
2288: . e - The element
2289: . v - The values
2290: - mode - either ADD_VALUES or INSERT_VALUES, where
2291: ADD_VALUES adds values to any existing entries, and
2292: INSERT_VALUES replaces existing entries with new values
2294: Level: beginner
2296: .seealso: MatSetOption()
2297: @*/
2298: PetscErrorCode DMMeshAssembleMatrixDM(Mat A, DM dm, SectionReal section, PetscInt e, PetscScalar v[], InsertMode mode)
2299: {
2303: PetscLogEventBegin(DMMesh_assembleMatrix,0,0,0,0);
2304: try {
2305: Obj<PETSC_MESH_TYPE> m;
2306: Obj<PETSC_MESH_TYPE::real_section_type> s;
2308: DMMeshGetMesh(dm, m);
2309: SectionRealGetSection(section, s);
2310: const ALE::Obj<PETSC_MESH_TYPE::order_type>& globalOrder = m->getFactory()->getGlobalOrder(m, s->getName(), s);
2312: if (m->debug()) {
2313: std::cout << "Assembling matrix for element number " << e << " --> point " << e << std::endl;
2314: }
2315: DMMeshUpdateOperator(A, m, s, globalOrder, e, v, mode);
2316: } catch (ALE::Exception e) {
2317: std::cout << e.msg() << std::endl;
2318: }
2319: PetscLogEventEnd(DMMesh_assembleMatrix,0,0,0,0);
2320: return(0);
2321: }
2323: /******************************** C Wrappers **********************************/
2327: PetscErrorCode WriteVTKHeader(DM dm, PetscViewer viewer)
2328: {
2329: ALE::Obj<PETSC_MESH_TYPE> m;
2330: PetscErrorCode ierr;
2332: DMMeshGetMesh(dm, m);
2333: return VTKViewer::writeHeader(m, viewer);
2334: }
2338: PetscErrorCode WriteVTKVertices(DM dm, PetscViewer viewer)
2339: {
2340: ALE::Obj<PETSC_MESH_TYPE> m;
2341: PetscErrorCode ierr;
2343: DMMeshGetMesh(dm, m);
2344: return VTKViewer::writeVertices(m, viewer);
2345: }
2349: PetscErrorCode WriteVTKElements(DM dm, PetscViewer viewer)
2350: {
2351: ALE::Obj<PETSC_MESH_TYPE> m;
2352: PetscErrorCode ierr;
2354: DMMeshGetMesh(dm, m);
2355: return VTKViewer::writeElements(m, viewer);
2356: }
2360: /*@C
2361: DMMeshGetCoordinates - Creates an array holding the coordinates.
2363: Not Collective
2365: Input Parameter:
2366: + dm - The DMMesh object
2367: - columnMajor - Flag for column major order
2369: Output Parameter:
2370: + numVertices - The number of vertices
2371: . dim - The embedding dimension
2372: - coords - The array holding local coordinates
2374: Level: intermediate
2376: .keywords: mesh, coordinates
2377: .seealso: DMMeshCreate()
2378: @*/
2379: PetscErrorCode DMMeshGetCoordinates(DM dm, PetscBool columnMajor, PetscInt *numVertices, PetscInt *dim, PetscReal *coords[])
2380: {
2381: ALE::Obj<PETSC_MESH_TYPE> m;
2382: PetscErrorCode ierr;
2385: DMMeshGetMesh(dm, m);
2386: ALE::PCICE::Builder::outputVerticesLocal(m, numVertices, dim, coords, columnMajor);
2387: return(0);
2388: }
2392: /*@C
2393: DMMeshGetElements - Creates an array holding the vertices on each element.
2395: Not Collective
2397: Input Parameters:
2398: + dm - The DMMesh object
2399: - columnMajor - Flag for column major order
2401: Output Parameters:
2402: + numElements - The number of elements
2403: . numCorners - The number of vertices per element
2404: - vertices - The array holding vertices on each local element
2406: Level: intermediate
2408: .keywords: mesh, elements
2409: .seealso: DMMeshCreate()
2410: @*/
2411: PetscErrorCode DMMeshGetElements(DM dm, PetscBool columnMajor, PetscInt *numElements, PetscInt *numCorners, PetscInt *vertices[])
2412: {
2413: ALE::Obj<PETSC_MESH_TYPE> m;
2414: PetscErrorCode ierr;
2417: DMMeshGetMesh(dm, m);
2418: ALE::PCICE::Builder::outputElementsLocal(m, numElements, numCorners, vertices, columnMajor);
2419: return(0);
2420: }
2424: PetscErrorCode DMMeshCreateNeighborCSR(DM dm, PetscInt *numVertices, PetscInt **offsets, PetscInt **adjacency)
2425: {
2426: const PetscInt maxFaceCases = 30;
2427: PetscInt numFaceCases = 0;
2428: PetscInt numFaceVertices[maxFaceCases];
2429: PetscInt *off, *adj;
2430: PetscInt dim, depth, cStart, cEnd, c, numCells, cell;
2434: /* For parallel partitioning, I think you have to communicate supports */
2435: DMMeshGetDimension(dm, &dim);
2436: DMMeshGetLabelSize(dm, "depth", &depth);
2437: --depth;
2438: DMMeshGetHeightStratum(dm, 0, &cStart, &cEnd);
2439: if (cEnd - cStart == 0) {
2440: *numVertices = 0;
2441: *offsets = NULL;
2442: *adjacency = NULL;
2443: return(0);
2444: }
2445: numCells = cEnd - cStart;
2446: /* Setup face recognition */
2447: {
2448: 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 */
2450: for (c = cStart; c < cEnd; ++c) {
2451: PetscInt corners;
2453: DMMeshGetConeSize(dm, c, &corners);
2454: if (!cornersSeen[corners]) {
2455: if (numFaceCases >= maxFaceCases) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Exceeded maximum number of face recognition cases");
2456: cornersSeen[corners] = 1;
2457: if (corners == dim+1) {
2458: numFaceVertices[numFaceCases] = dim;
2459: PetscInfo(dm, "Recognizing simplices\n");
2460: } else if ((dim == 1) && (corners == 3)) {
2461: numFaceVertices[numFaceCases] = 3;
2462: PetscInfo(dm, "Recognizing quadratic edges\n");
2463: } else if ((dim == 2) && (corners == 4)) {
2464: numFaceVertices[numFaceCases] = 2;
2465: PetscInfo(dm, "Recognizing quads\n");
2466: } else if ((dim == 2) && (corners == 6)) {
2467: numFaceVertices[numFaceCases] = 3;
2468: PetscInfo(dm, "Recognizing tri and quad cohesive Lagrange cells\n");
2469: } else if ((dim == 2) && (corners == 9)) {
2470: numFaceVertices[numFaceCases] = 3;
2471: PetscInfo(dm, "Recognizing quadratic quads and quadratic quad cohesive Lagrange cells\n");
2472: } else if ((dim == 3) && (corners == 6)) {
2473: numFaceVertices[numFaceCases] = 4;
2474: PetscInfo(dm, "Recognizing tet cohesive cells\n");
2475: } else if ((dim == 3) && (corners == 8)) {
2476: numFaceVertices[numFaceCases] = 4;
2477: PetscInfo(dm, "Recognizing hexes\n");
2478: } else if ((dim == 3) && (corners == 9)) {
2479: numFaceVertices[numFaceCases] = 6;
2480: PetscInfo(dm, "Recognizing tet cohesive Lagrange cells\n");
2481: } else if ((dim == 3) && (corners == 10)) {
2482: numFaceVertices[numFaceCases] = 6;
2483: PetscInfo(dm, "Recognizing quadratic tets\n");
2484: } else if ((dim == 3) && (corners == 12)) {
2485: numFaceVertices[numFaceCases] = 6;
2486: PetscInfo(dm, "Recognizing hex cohesive Lagrange cells\n");
2487: } else if ((dim == 3) && (corners == 18)) {
2488: numFaceVertices[numFaceCases] = 6;
2489: PetscInfo(dm, "Recognizing quadratic tet cohesive Lagrange cells\n");
2490: } else if ((dim == 3) && (corners == 27)) {
2491: numFaceVertices[numFaceCases] = 9;
2492: PetscInfo(dm, "Recognizing quadratic hexes and quadratic hex cohesive Lagrange cells\n");
2493: } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Could not recognize number of face vertices for %d corners", corners);
2494: ++numFaceCases;
2495: }
2496: }
2497: }
2498: /* Check for optimized depth 1 construction */
2499: PetscMalloc((numCells+1) * sizeof(PetscInt), &off);
2500: PetscMemzero(off, (numCells+1) * sizeof(PetscInt));
2501: if (depth == 1) {
2502: PetscInt *neighborCells;
2503: PetscInt n;
2504: PetscInt maxConeSize, maxSupportSize;
2506: /* Temp space for point adj <= maxConeSize*maxSupportSize */
2507: DMMeshGetMaxSizes(dm, &maxConeSize, &maxSupportSize);
2508: PetscMalloc(maxConeSize*maxSupportSize * sizeof(PetscInt), &neighborCells);
2509: /* Count neighboring cells */
2510: for (cell = cStart; cell < cEnd; ++cell) {
2511: const PetscInt *cone;
2512: PetscInt numNeighbors = 0;
2513: PetscInt coneSize, c;
2515: /* Get support of the cone, and make a set of the cells */
2516: DMMeshGetConeSize(dm, cell, &coneSize);
2517: DMMeshGetCone(dm, cell, &cone);
2518: for (c = 0; c < coneSize; ++c) {
2519: const PetscInt *support;
2520: PetscInt supportSize, s;
2522: DMMeshGetSupportSize(dm, cone[c], &supportSize);
2523: DMMeshGetSupport(dm, cone[c], &support);
2524: for (s = 0; s < supportSize; ++s) {
2525: const PetscInt point = support[s];
2527: if (point == cell) continue;
2528: for (n = 0; n < numNeighbors; ++n) {
2529: if (neighborCells[n] == point) break;
2530: }
2531: if (n == numNeighbors) {
2532: neighborCells[n] = point;
2533: ++numNeighbors;
2534: }
2535: }
2536: }
2537: /* Get meet with each cell, and check with recognizer (could optimize to check each pair only once) */
2538: for (n = 0; n < numNeighbors; ++n) {
2539: PetscInt cellPair[2] = {cell, neighborCells[n]};
2540: PetscInt meetSize;
2541: const PetscInt *meet;
2543: DMMeshMeetPoints(dm, 2, cellPair, &meetSize, &meet);
2544: if (meetSize) {
2545: PetscInt f;
2547: for (f = 0; f < numFaceCases; ++f) {
2548: if (numFaceVertices[f] == meetSize) {
2549: ++off[cell-cStart+1];
2550: break;
2551: }
2552: }
2553: }
2554: }
2555: }
2556: /* Prefix sum */
2557: for (cell = 1; cell <= numCells; ++cell) off[cell] += off[cell-1];
2558: PetscMalloc(off[numCells] * sizeof(PetscInt), &adj);
2559: /* Get neighboring cells */
2560: for (cell = cStart; cell < cEnd; ++cell) {
2561: const PetscInt *cone;
2562: PetscInt numNeighbors = 0;
2563: PetscInt cellOffset = 0;
2564: PetscInt coneSize, c;
2566: /* Get support of the cone, and make a set of the cells */
2567: DMMeshGetConeSize(dm, cell, &coneSize);
2568: DMMeshGetCone(dm, cell, &cone);
2569: for (c = 0; c < coneSize; ++c) {
2570: const PetscInt *support;
2571: PetscInt supportSize, s;
2573: DMMeshGetSupportSize(dm, cone[c], &supportSize);
2574: DMMeshGetSupport(dm, cone[c], &support);
2575: for (s = 0; s < supportSize; ++s) {
2576: const PetscInt point = support[s];
2578: if (point == cell) continue;
2579: for (n = 0; n < numNeighbors; ++n) {
2580: if (neighborCells[n] == point) break;
2581: }
2582: if (n == numNeighbors) {
2583: neighborCells[n] = point;
2584: ++numNeighbors;
2585: }
2586: }
2587: }
2588: /* Get meet with each cell, and check with recognizer (could optimize to check each pair only once) */
2589: for (n = 0; n < numNeighbors; ++n) {
2590: PetscInt cellPair[2] = {cell, neighborCells[n]};
2591: PetscInt meetSize;
2592: const PetscInt *meet;
2594: DMMeshMeetPoints(dm, 2, cellPair, &meetSize, &meet);
2595: if (meetSize) {
2596: PetscInt f;
2598: for (f = 0; f < numFaceCases; ++f) {
2599: if (numFaceVertices[f] == meetSize) {
2600: adj[off[cell-cStart]+cellOffset] = neighborCells[n];
2601: ++cellOffset;
2602: break;
2603: }
2604: }
2605: }
2606: }
2607: }
2608: PetscFree(neighborCells);
2609: } else if (depth == dim) {
2610: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Neighbor graph creation not implemented for interpolated meshes");
2611: #if 0
2612: OffsetVisitor<typename Mesh::sieve_type> oV(*sieve, *overlapSieve, off);
2613: PetscInt p;
2615: for (typename Mesh::label_sequence::iterator c_iter = cells->begin(); c_iter != cEnd; ++c_iter) {
2616: sieve->cone(*c_iter, oV);
2617: }
2618: for (p = 1; p <= numCells; ++p) off[p] = off[p] + off[p-1];
2619: PetscMalloc(off[numCells] * sizeof(PetscInt), &adj);
2620: AdjVisitor<typename Mesh::sieve_type> aV(adj, zeroBase);
2621: ISieveVisitor::SupportVisitor<typename Mesh::sieve_type, AdjVisitor<typename Mesh::sieve_type> > sV(*sieve, aV);
2622: ISieveVisitor::SupportVisitor<typename Mesh::sieve_type, AdjVisitor<typename Mesh::sieve_type> > ovSV(*overlapSieve, aV);
2624: for (typename Mesh::label_sequence::iterator c_iter = cells->begin(); c_iter != cEnd; ++c_iter) {
2625: aV.setCell(*c_iter);
2626: sieve->cone(*c_iter, sV);
2627: sieve->cone(*c_iter, ovSV);
2628: }
2629: offset = aV.getOffset();
2630: #endif
2631: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Neighbor graph creation not defined for partially interpolated meshes");
2633: *numVertices = numCells;
2634: *offsets = off;
2635: *adjacency = adj;
2636: return(0);
2637: }
2639: #if defined(PETSC_HAVE_CHACO)
2640: #if defined(PETSC_HAVE_UNISTD_H)
2641: #include <unistd.h>
2642: #endif
2643: /* Chaco does not have an include file */
2644: extern "C" {
2645: extern int interface(int nvtxs, int *start, int *adjacency, int *vwgts,
2646: float *ewgts, float *x, float *y, float *z, char *outassignname,
2647: char *outfilename, short *assignment, int architecture, int ndims_tot,
2648: int mesh_dims[3], double *goal, int global_method, int local_method,
2649: int rqi_flag, int vmax, int ndims, double eigtol, long seed);
2651: extern int FREE_GRAPH;
2652: }
2656: PetscErrorCode DMMeshPartition_Chaco(DM dm, PetscInt numVertices, PetscInt start[], PetscInt adjacency[], PetscSection *partSection, IS *partition)
2657: {
2658: enum {DEFAULT_METHOD = 1, INERTIAL_METHOD = 3};
2659: MPI_Comm comm = ((PetscObject) dm)->comm;
2660: int nvtxs = numVertices; /* number of vertices in full graph */
2661: int *vwgts = NULL; /* weights for all vertices */
2662: float *ewgts = NULL; /* weights for all edges */
2663: float *x = NULL, *y = NULL, *z = NULL; /* coordinates for inertial method */
2664: char *outassignname = NULL; /* name of assignment output file */
2665: char *outfilename = NULL; /* output file name */
2666: int architecture = 1; /* 0 => hypercube, d => d-dimensional mesh */
2667: int ndims_tot = 0; /* total number of cube dimensions to divide */
2668: int mesh_dims[3]; /* dimensions of mesh of processors */
2669: double *goal = NULL; /* desired set sizes for each set */
2670: int global_method = 1; /* global partitioning algorithm */
2671: int local_method = 1; /* local partitioning algorithm */
2672: int rqi_flag = 0; /* should I use RQI/Symmlq eigensolver? */
2673: int vmax = 200; /* how many vertices to coarsen down to? */
2674: int ndims = 1; /* number of eigenvectors (2^d sets) */
2675: double eigtol = 0.001; /* tolerance on eigenvectors */
2676: long seed = 123636512; /* for random graph mutations */
2677: short int *assignment; /* Output partition */
2678: int fd_stdout, fd_pipe[2];
2679: PetscInt *points;
2680: PetscMPIInt commSize;
2681: int i, v, p;
2685: MPI_Comm_size(comm, &commSize);
2686: if (!numVertices) {
2687: PetscSectionCreate(comm, partSection);
2688: PetscSectionSetChart(*partSection, 0, commSize);
2689: PetscSectionSetUp(*partSection);
2690: ISCreateGeneral(comm, 0, NULL, PETSC_OWN_POINTER, partition);
2691: return(0);
2692: }
2693: FREE_GRAPH = 0; /* Do not let Chaco free my memory */
2694: for (i = 0; i < start[numVertices]; ++i) ++adjacency[i];
2696: if (global_method == INERTIAL_METHOD) SETERRQ(comm, PETSC_ERR_SUP, "Inertial partitioning not yet supported");
2698: mesh_dims[0] = commSize;
2699: mesh_dims[1] = 1;
2700: mesh_dims[2] = 1;
2701: PetscMalloc(nvtxs * sizeof(short int), &assignment);
2702: /* Chaco outputs to stdout. We redirect this to a buffer. */
2703: /* TODO: check error codes for UNIX calls */
2704: #if defined(PETSC_HAVE_UNISTD_H)
2705: {
2706: fd_stdout = dup(1);
2707: pipe(fd_pipe);
2708: close(1);
2709: dup2(fd_pipe[1], 1);
2710: }
2711: #endif
2712: interface(nvtxs, (int*) start, (int*) adjacency, vwgts, ewgts, x, y, z, outassignname, outfilename,
2713: assignment, architecture, ndims_tot, mesh_dims, goal, global_method, local_method, rqi_flag,
2714: vmax, ndims, eigtol, seed);
2715: #if defined(PETSC_HAVE_UNISTD_H)
2716: {
2717: char msgLog[10000];
2718: int count;
2720: fflush(stdout);
2721: count = read(fd_pipe[0], msgLog, (10000-1)*sizeof(char));
2722: if (count < 0) count = 0;
2723: msgLog[count] = 0;
2724: close(1);
2725: dup2(fd_stdout, 1);
2726: close(fd_stdout);
2727: close(fd_pipe[0]);
2728: close(fd_pipe[1]);
2729: if (ierr) SETERRQ1(comm, PETSC_ERR_LIB, "Error in Chaco library: %s", msgLog);
2730: }
2731: #endif
2732: /* Convert to PetscSection+IS */
2733: PetscSectionCreate(comm, partSection);
2734: PetscSectionSetChart(*partSection, 0, commSize);
2735: for (v = 0; v < nvtxs; ++v) {
2736: PetscSectionAddDof(*partSection, assignment[v], 1);
2737: }
2738: PetscSectionSetUp(*partSection);
2739: PetscMalloc(nvtxs * sizeof(PetscInt), &points);
2740: for (p = 0, i = 0; p < commSize; ++p) {
2741: for (v = 0; v < nvtxs; ++v) {
2742: if (assignment[v] == p) points[i++] = v;
2743: }
2744: }
2745: if (i != nvtxs) SETERRQ2(comm, PETSC_ERR_PLIB, "Number of points %d should be %d", i, nvtxs);
2746: ISCreateGeneral(comm, nvtxs, points, PETSC_OWN_POINTER, partition);
2747: if (global_method == INERTIAL_METHOD) {
2748: /* manager.destroyCellCoordinates(nvtxs, &x, &y, &z); */
2749: }
2750: PetscFree(assignment);
2751: for (i = 0; i < start[numVertices]; ++i) --adjacency[i];
2752: return(0);
2753: }
2754: #endif
2756: #if defined(PETSC_HAVE_PARMETIS)
2759: PetscErrorCode DMMeshPartition_ParMetis(DM dm, PetscInt numVertices, PetscInt start[], PetscInt adjacency[], PetscSection *partSection, IS *partition)
2760: {
2762: return(0);
2763: }
2764: #endif
2768: PetscErrorCode DMMeshCreatePartition(DM dm, PetscSection *partSection, IS *partition, PetscInt height)
2769: {
2770: PetscMPIInt size;
2774: MPI_Comm_size(((PetscObject) dm)->comm, &size);
2775: if (size == 1) {
2776: PetscInt *points;
2777: PetscInt cStart, cEnd, c;
2779: DMMeshGetHeightStratum(dm, 0, &cStart, &cEnd);
2780: PetscSectionCreate(((PetscObject) dm)->comm, partSection);
2781: PetscSectionSetChart(*partSection, 0, size);
2782: PetscSectionSetDof(*partSection, 0, cEnd-cStart);
2783: PetscSectionSetUp(*partSection);
2784: PetscMalloc((cEnd - cStart) * sizeof(PetscInt), &points);
2785: for (c = cStart; c < cEnd; ++c) points[c] = c;
2786: ISCreateGeneral(((PetscObject) dm)->comm, cEnd-cStart, points, PETSC_OWN_POINTER, partition);
2787: return(0);
2788: }
2789: if (height == 0) {
2790: PetscInt numVertices;
2791: PetscInt *start = NULL;
2792: PetscInt *adjacency = NULL;
2794: if (1) {
2795: DMMeshCreateNeighborCSR(dm, &numVertices, &start, &adjacency);
2796: #if defined(PETSC_HAVE_CHACO)
2797: DMMeshPartition_Chaco(dm, numVertices, start, adjacency, partSection, partition);
2798: #endif
2799: } else {
2800: DMMeshCreateNeighborCSR(dm, &numVertices, &start, &adjacency);
2801: #if defined(PETSC_HAVE_PARMETIS)
2802: DMMeshPartition_ParMetis(dm, numVertices, start, adjacency, partSection, partition);
2803: #endif
2804: }
2805: PetscFree(start);
2806: PetscFree(adjacency);
2807: # if 0
2808: } else if (height == 1) {
2809: /* Build the dual graph for faces and partition the hypergraph */
2810: PetscInt numEdges;
2812: buildFaceCSRV(mesh, mesh->getFactory()->getNumbering(mesh, mesh->depth()-1), &numEdges, &start, &adjacency, GraphPartitioner::zeroBase());
2813: GraphPartitioner().partition(numEdges, start, adjacency, partition, manager);
2814: destroyCSR(numEdges, start, adjacency);
2815: #endif
2816: } else SETERRQ1(((PetscObject) dm)->comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid partition height %d", height);
2817: return(0);
2818: }
2822: PetscErrorCode DMMeshCreatePartitionClosure(DM dm, PetscSection pointSection, IS pointPartition, PetscSection *section, IS *partition)
2823: {
2824: const PetscInt *partArray;
2825: PetscInt *allPoints, *partPoints = NULL;
2826: PetscInt rStart, rEnd, rank, maxPartSize = 0, newSize;
2830: PetscSectionGetChart(pointSection, &rStart, &rEnd);
2831: ISGetIndices(pointPartition, &partArray);
2832: PetscSectionCreate(((PetscObject) dm)->comm, section);
2833: PetscSectionSetChart(*section, rStart, rEnd);
2834: for (rank = rStart; rank < rEnd; ++rank) {
2835: PetscInt partSize = 0;
2836: PetscInt numPoints, offset, p;
2838: PetscSectionGetDof(pointSection, rank, &numPoints);
2839: PetscSectionGetOffset(pointSection, rank, &offset);
2840: for (p = 0; p < numPoints; ++p) {
2841: PetscInt point = partArray[offset+p], closureSize;
2842: const PetscInt *closure;
2844: /* TODO Include support for height > 0 case */
2845: DMMeshGetTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
2846: /* Merge into existing points */
2847: if (partSize+closureSize > maxPartSize) {
2848: PetscInt *tmpPoints;
2850: maxPartSize = PetscMax(partSize+closureSize, 2*maxPartSize);
2851: PetscMalloc(maxPartSize * sizeof(PetscInt), &tmpPoints);
2852: PetscMemcpy(tmpPoints, partPoints, partSize * sizeof(PetscInt));
2853: PetscFree(partPoints);
2854: partPoints = tmpPoints;
2855: }
2856: PetscMemcpy(&partPoints[partSize], closure, closureSize * sizeof(PetscInt));
2857: partSize += closureSize;
2858: PetscSortRemoveDupsInt(&partSize, partPoints);
2859: }
2860: PetscSectionSetDof(*section, rank, partSize);
2861: }
2862: PetscSectionSetUp(*section);
2863: PetscSectionGetStorageSize(*section, &newSize);
2864: PetscMalloc(newSize * sizeof(PetscInt), &allPoints);
2866: for (rank = rStart; rank < rEnd; ++rank) {
2867: PetscInt partSize = 0, newOffset;
2868: PetscInt numPoints, offset, p;
2870: PetscSectionGetDof(pointSection, rank, &numPoints);
2871: PetscSectionGetOffset(pointSection, rank, &offset);
2872: for (p = 0; p < numPoints; ++p) {
2873: PetscInt point = partArray[offset+p], closureSize;
2874: const PetscInt *closure;
2876: /* TODO Include support for height > 0 case */
2877: DMMeshGetTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
2878: /* Merge into existing points */
2879: PetscMemcpy(&partPoints[partSize], closure, closureSize * sizeof(PetscInt));
2880: partSize += closureSize;
2881: PetscSortRemoveDupsInt(&partSize, partPoints);
2882: }
2883: PetscSectionGetOffset(*section, rank, &newOffset);
2884: PetscMemcpy(&allPoints[newOffset], partPoints, partSize * sizeof(PetscInt));
2885: }
2886: ISRestoreIndices(pointPartition, &partArray);
2887: PetscFree(partPoints);
2888: ISCreateGeneral(((PetscObject) dm)->comm, newSize, allPoints, PETSC_OWN_POINTER, partition);
2889: return(0);
2890: }
2894: /*
2895: Input Parameters:
2896: . originalSection
2897: , originalVec
2899: Output Parameters:
2900: . newSection
2901: . newVec
2902: */
2903: PetscErrorCode DMMeshDistributeField(DM dm, PetscSF pointSF, PetscSection originalSection, Vec originalVec, PetscSection newSection, Vec newVec)
2904: {
2905: PetscSF fieldSF;
2906: PetscInt *remoteOffsets, fieldSize;
2907: PetscScalar *originalValues, *newValues;
2911: PetscSFDistributeSection(pointSF, originalSection, &remoteOffsets, newSection);
2913: PetscSectionGetStorageSize(newSection, &fieldSize);
2914: VecSetSizes(newVec, fieldSize, PETSC_DETERMINE);
2915: VecSetFromOptions(newVec);
2917: VecGetArray(originalVec, &originalValues);
2918: VecGetArray(newVec, &newValues);
2919: PetscSFCreateSectionSF(pointSF, originalSection, remoteOffsets, newSection, &fieldSF);
2920: PetscSFBcastBegin(fieldSF, MPIU_SCALAR, originalValues, newValues);
2921: PetscSFBcastEnd(fieldSF, MPIU_SCALAR, originalValues, newValues);
2922: PetscSFDestroy(&fieldSF);
2923: VecRestoreArray(newVec, &newValues);
2924: VecRestoreArray(originalVec, &originalValues);
2925: return(0);
2926: }
2930: /*@C
2931: DMMeshDistribute - Distributes the mesh and any associated sections.
2933: Not Collective
2935: Input Parameter:
2936: + dm - The original DMMesh object
2937: - partitioner - The partitioning package, or NULL for the default
2939: Output Parameter:
2940: . parallelMesh - The distributed DMMesh object
2942: Level: intermediate
2944: .keywords: mesh, elements
2946: .seealso: DMMeshCreate(), DMMeshDistributeByFace()
2947: @*/
2948: PetscErrorCode DMMeshDistribute(DM dm, const char partitioner[], DM *dmParallel)
2949: {
2950: DM_Mesh *mesh = (DM_Mesh*) dm->data, *pmesh;
2951: MPI_Comm comm = ((PetscObject) dm)->comm;
2952: PetscMPIInt rank, numProcs, p;
2958: MPI_Comm_rank(comm, &rank);
2959: MPI_Comm_size(comm, &numProcs);
2960: if (numProcs == 1) return(0);
2961: if (mesh->useNewImpl) {
2962: const PetscInt height = 0;
2963: PetscInt dim, numRemoteRanks;
2964: IS cellPart, part;
2965: PetscSection cellPartSection, partSection;
2966: PetscSFNode *remoteRanks;
2967: PetscSF partSF, pointSF, coneSF;
2968: ISLocalToGlobalMapping renumbering;
2969: PetscSection originalConeSection, newConeSection;
2970: PetscInt *remoteOffsets;
2971: PetscInt *cones, *newCones, newConesSize;
2972: PetscBool flg;
2974: DMMeshGetDimension(dm, &dim);
2975: /* Create cell partition - We need to rewrite to use IS, use the MatPartition stuff */
2976: DMMeshCreatePartition(dm, &cellPartSection, &cellPart, height);
2977: /* Create SF assuming a serial partition for all processes: Could check for IS length here */
2978: if (!rank) numRemoteRanks = numProcs;
2979: else numRemoteRanks = 0;
2980: PetscMalloc(numRemoteRanks * sizeof(PetscSFNode), &remoteRanks);
2981: for (p = 0; p < numRemoteRanks; ++p) {
2982: remoteRanks[p].rank = p;
2983: remoteRanks[p].index = 0;
2984: }
2985: PetscSFCreate(comm, &partSF);
2986: PetscSFSetGraph(partSF, 1, numRemoteRanks, NULL, PETSC_OWN_POINTER, remoteRanks, PETSC_OWN_POINTER);
2987: PetscOptionsHasName(((PetscObject) dm)->prefix, "-partition_view", &flg);
2988: if (flg) {
2989: PetscPrintf(comm, "Cell Partition:\n");
2990: PetscSectionView(cellPartSection, PETSC_VIEWER_STDOUT_WORLD);
2991: ISView(cellPart, NULL);
2992: PetscSFView(partSF, NULL);
2993: }
2994: /* Close the partition over the mesh */
2995: DMMeshCreatePartitionClosure(dm, cellPartSection, cellPart, &partSection, &part);
2996: ISDestroy(&cellPart);
2997: PetscSectionDestroy(&cellPartSection);
2998: /* Create new mesh */
2999: DMMeshCreate(comm, dmParallel);
3000: DMMeshSetDimension(*dmParallel, dim);
3001: PetscObjectSetName((PetscObject) *dmParallel, "Parallel Mesh");
3002: pmesh = (DM_Mesh*) (*dmParallel)->data;
3003: /* Distribute sieve points and the global point numbering (replaces creating remote bases) */
3004: PetscSFConvertPartition(partSF, partSection, part, &renumbering, &pointSF);
3005: if (flg) {
3006: PetscPrintf(comm, "Point Partition:\n");
3007: PetscSectionView(partSection, PETSC_VIEWER_STDOUT_WORLD);
3008: ISView(part, NULL);
3009: PetscSFView(pointSF, NULL);
3010: PetscPrintf(comm, "Point Renumbering after partition:\n");
3011: ISLocalToGlobalMappingView(renumbering, NULL);
3012: }
3013: /* Distribute cone section */
3014: DMMeshGetConeSection(dm, &originalConeSection);
3015: DMMeshGetConeSection(*dmParallel, &newConeSection);
3016: PetscSFDistributeSection(pointSF, originalConeSection, &remoteOffsets, newConeSection);
3017: DMMeshSetUp(*dmParallel);
3018: {
3019: PetscInt pStart, pEnd, p;
3021: PetscSectionGetChart(newConeSection, &pStart, &pEnd);
3022: for (p = pStart; p < pEnd; ++p) {
3023: PetscInt coneSize;
3024: PetscSectionGetDof(newConeSection, p, &coneSize);
3025: pmesh->maxConeSize = PetscMax(pmesh->maxConeSize, coneSize);
3026: }
3027: }
3028: /* Communicate and renumber cones */
3029: PetscSFCreateSectionSF(pointSF, originalConeSection, remoteOffsets, newConeSection, &coneSF);
3030: DMMeshGetCones(dm, &cones);
3031: DMMeshGetCones(*dmParallel, &newCones);
3032: PetscSFBcastBegin(coneSF, MPIU_INT, cones, newCones);
3033: PetscSFBcastEnd(coneSF, MPIU_INT, cones, newCones);
3034: PetscSectionGetStorageSize(newConeSection, &newConesSize);
3035: ISGlobalToLocalMappingApply(renumbering, IS_GTOLM_MASK, newConesSize, newCones, NULL, newCones);
3036: PetscOptionsHasName(((PetscObject) dm)->prefix, "-cones_view", &flg);
3037: if (flg) {
3038: PetscPrintf(comm, "Serial Cone Section:\n");
3039: PetscSectionView(originalConeSection, PETSC_VIEWER_STDOUT_WORLD);
3040: PetscPrintf(comm, "Parallel Cone Section:\n");
3041: PetscSectionView(newConeSection, PETSC_VIEWER_STDOUT_WORLD);
3042: PetscSFView(coneSF, NULL);
3043: }
3044: PetscSFDestroy(&coneSF);
3045: /* Create supports and stratify sieve */
3046: DMMeshSymmetrize(*dmParallel);
3047: DMMeshStratify(*dmParallel);
3048: /* Distribute Coordinates */
3049: {
3050: PetscSection originalCoordSection, newCoordSection;
3051: Vec originalCoordinates, newCoordinates;
3053: DMMeshGetCoordinateSection(dm, &originalCoordSection);
3054: DMMeshGetCoordinateSection(*dmParallel, &newCoordSection);
3055: DMMeshGetCoordinateVec(dm, &originalCoordinates);
3056: DMMeshGetCoordinateVec(*dmParallel, &newCoordinates);
3058: DMMeshDistributeField(dm, pointSF, originalCoordSection, originalCoordinates, newCoordSection, newCoordinates);
3059: }
3060: /* Distribute labels */
3061: {
3062: SieveLabel next = mesh->labels, newNext = NULL;
3063: PetscInt numLabels = 0, l;
3065: /* Bcast number of labels */
3066: while (next) {
3067: ++numLabels; next = next->next;
3068: }
3069: MPI_Bcast(&numLabels, 1, MPIU_INT, 0, comm);
3070: next = mesh->labels;
3071: for (l = 0; l < numLabels; ++l) {
3072: SieveLabel newLabel;
3073: const PetscInt *partArray;
3074: PetscInt *stratumSizes = NULL, *points = NULL;
3075: PetscMPIInt *sendcnts = NULL, *offsets = NULL, *displs = NULL;
3076: PetscInt nameSize, s, p;
3077: size_t len = 0;
3079: PetscNew(struct Sieve_Label, &newLabel);
3080: /* Bcast name (could filter for no points) */
3081: if (!rank) {PetscStrlen(next->name, &len);}
3082: nameSize = len;
3083: MPI_Bcast(&nameSize, 1, MPIU_INT, 0, comm);
3084: PetscMalloc(nameSize+1, &newLabel->name);
3085: if (!rank) {PetscMemcpy(newLabel->name, next->name, nameSize+1);}
3086: MPI_Bcast(newLabel->name, nameSize+1, MPI_CHAR, 0, comm);
3087: /* Bcast numStrata (could filter for no points in stratum) */
3088: if (!rank) newLabel->numStrata = next->numStrata;
3089: MPI_Bcast(&newLabel->numStrata, 1, MPIU_INT, 0, comm);
3090: PetscMalloc(newLabel->numStrata * sizeof(PetscInt), &newLabel->stratumValues);
3091: PetscMalloc(newLabel->numStrata * sizeof(PetscInt), &newLabel->stratumSizes);
3092: PetscMalloc((newLabel->numStrata+1) * sizeof(PetscInt), &newLabel->stratumOffsets);
3093: /* Bcast stratumValues (could filter for no points in stratum) */
3094: if (!rank) {PetscMemcpy(newLabel->stratumValues, next->stratumValues, next->numStrata * sizeof(PetscInt));}
3095: MPI_Bcast(newLabel->stratumValues, newLabel->numStrata, MPIU_INT, 0, comm);
3096: /* Find size on each process and Scatter */
3097: if (!rank) {
3098: ISGetIndices(part, &partArray);
3099: PetscMalloc(numProcs*next->numStrata * sizeof(PetscInt), &stratumSizes);
3100: PetscMemzero(stratumSizes, numProcs*next->numStrata * sizeof(PetscInt));
3101: for (s = 0; s < next->numStrata; ++s) {
3102: for (p = next->stratumOffsets[s]; p < next->stratumOffsets[s]+next->stratumSizes[s]; ++p) {
3103: const PetscInt point = next->points[p];
3104: PetscInt proc;
3106: for (proc = 0; proc < numProcs; ++proc) {
3107: PetscInt dof, off, pPart;
3109: PetscSectionGetDof(partSection, proc, &dof);
3110: PetscSectionGetOffset(partSection, proc, &off);
3111: for (pPart = off; pPart < off+dof; ++pPart) {
3112: if (partArray[pPart] == point) {
3113: ++stratumSizes[proc*next->numStrata+s];
3114: break;
3115: }
3116: }
3117: }
3118: }
3119: }
3120: ISRestoreIndices(part, &partArray);
3121: }
3122: MPI_Scatter(stratumSizes, newLabel->numStrata, MPI_INT, newLabel->stratumSizes, newLabel->numStrata, MPI_INT, 0, comm);
3123: /* Calculate stratumOffsets */
3124: newLabel->stratumOffsets[0] = 0;
3125: for (s = 0; s < newLabel->numStrata; ++s) {
3126: newLabel->stratumOffsets[s+1] = newLabel->stratumSizes[s] + newLabel->stratumOffsets[s];
3127: }
3128: /* Pack points and Scatter */
3129: if (!rank) {
3130: PetscMalloc3(numProcs,PetscMPIInt,&sendcnts,numProcs,PetscMPIInt,&offsets,numProcs+1,PetscMPIInt,&displs);
3131: displs[0] = 0;
3132: for (p = 0; p < numProcs; ++p) {
3133: sendcnts[p] = 0;
3134: for (s = 0; s < next->numStrata; ++s) {
3135: sendcnts[p] += stratumSizes[p*next->numStrata+s];
3136: }
3137: offsets[p] = displs[p];
3138: displs[p+1] = displs[p] + sendcnts[p];
3139: }
3140: PetscMalloc(displs[numProcs] * sizeof(PetscInt), &points);
3141: for (s = 0; s < next->numStrata; ++s) {
3142: for (p = next->stratumOffsets[s]; p < next->stratumOffsets[s]+next->stratumSizes[s]; ++p) {
3143: const PetscInt point = next->points[p];
3144: PetscInt proc;
3146: for (proc = 0; proc < numProcs; ++proc) {
3147: PetscInt dof, off, pPart;
3149: PetscSectionGetDof(partSection, proc, &dof);
3150: PetscSectionGetOffset(partSection, proc, &off);
3151: for (pPart = off; pPart < off+dof; ++pPart) {
3152: if (partArray[pPart] == point) {
3153: points[offsets[proc]++] = point;
3154: break;
3155: }
3156: }
3157: }
3158: }
3159: }
3160: }
3161: PetscMalloc(newLabel->stratumOffsets[newLabel->numStrata] * sizeof(PetscInt), &newLabel->points);
3162: MPI_Scatterv(points, sendcnts, displs, MPIU_INT, newLabel->points, newLabel->stratumOffsets[newLabel->numStrata], MPIU_INT, 0, comm);
3163: PetscFree(points);
3164: PetscFree3(sendcnts,offsets,displs);
3165: PetscFree(stratumSizes);
3166: /* Renumber points */
3167: ISGlobalToLocalMappingApply(renumbering, IS_GTOLM_MASK, newLabel->stratumOffsets[newLabel->numStrata], newLabel->points, NULL, newLabel->points);
3168: /* Sort points */
3169: for (s = 0; s < newLabel->numStrata; ++s) {
3170: PetscSortInt(newLabel->stratumSizes[s], &newLabel->points[newLabel->stratumOffsets[s]]);
3171: }
3172: /* Insert into list */
3173: if (newNext) newNext->next = newLabel;
3174: else pmesh->labels = newLabel;
3175: newNext = newLabel;
3176: if (!rank) next = next->next;
3177: }
3178: }
3179: /* Cleanup Partition */
3180: ISLocalToGlobalMappingDestroy(&renumbering);
3181: PetscSFDestroy(&partSF);
3182: PetscSectionDestroy(&partSection);
3183: ISDestroy(&part);
3184: /* Create point SF for parallel mesh */
3185: {
3186: const PetscInt *leaves;
3187: PetscSFNode *remotePoints;
3188: PetscMPIInt *rowners, *lowners;
3189: PetscInt *ghostPoints, numRoots, numLeaves, numGhostPoints = 0, p, gp;
3191: PetscSFGetGraph(pointSF, &numRoots, &numLeaves, &leaves, NULL);
3192: PetscMalloc2(numRoots*2,PetscInt,&rowners,numLeaves*2,PetscInt,&lowners);
3193: for (p = 0; p < numRoots*2; ++p) rowners[p] = 0;
3194: for (p = 0; p < numLeaves; ++p) {
3195: lowners[p*2+0] = rank;
3196: PetscMPIIntCast(leaves ? leaves[p] : p,lowners + 2*p + 1);
3197: }
3198: PetscSFReduceBegin(pointSF, MPI_2INT, lowners, rowners, MPI_MAXLOC);
3199: PetscSFReduceEnd(pointSF, MPI_2INT, lowners, rowners, MPI_MAXLOC);
3200: PetscSFBcastBegin(pointSF, MPI_2INT, rowners, lowners);
3201: PetscSFBcastEnd(pointSF, MPI_2INT, rowners, lowners);
3202: for (p = 0; p < numLeaves; ++p) {
3203: if (lowners[p*2+0] != rank) ++numGhostPoints;
3204: }
3205: PetscMalloc(numGhostPoints * sizeof(PetscInt), &ghostPoints);
3206: PetscMalloc(numGhostPoints * sizeof(PetscSFNode), &remotePoints);
3207: for (p = 0, gp = 0; p < numLeaves; ++p) {
3208: if (lowners[p*2+0] != rank) {
3209: ghostPoints[gp] = leaves ? leaves[p] : p;
3210: remotePoints[gp].rank = lowners[p*2+0];
3211: remotePoints[gp].index = lowners[p*2+1];
3212: ++gp;
3213: }
3214: }
3215: PetscFree2(rowners,lowners);
3216: PetscSFCreate(comm, &pmesh->sf);
3217: PetscSFSetGraph(pmesh->sf, PETSC_DETERMINE, numGhostPoints, ghostPoints, PETSC_OWN_POINTER, remotePoints, PETSC_OWN_POINTER);
3218: PetscSFSetFromOptions(pmesh->sf);
3219: }
3220: /* Cleanup */
3221: PetscSFDestroy(&pointSF);
3222: DMSetFromOptions(*dmParallel);
3223: } else {
3224: ALE::Obj<PETSC_MESH_TYPE> oldMesh;
3225: DMMeshGetMesh(dm, oldMesh);
3226: DMMeshCreate(comm, dmParallel);
3227: const Obj<PETSC_MESH_TYPE> newMesh = new PETSC_MESH_TYPE(comm, oldMesh->getDimension(), oldMesh->debug());
3228: const Obj<PETSC_MESH_TYPE::sieve_type> newSieve = new PETSC_MESH_TYPE::sieve_type(comm, oldMesh->debug());
3230: newMesh->setSieve(newSieve);
3231: ALE::DistributionNew<PETSC_MESH_TYPE>::distributeMeshAndSectionsV(oldMesh, newMesh);
3232: DMMeshSetMesh(*dmParallel, newMesh);
3233: }
3234: return(0);
3235: }
3239: /*@C
3240: DMMeshDistribute - Distributes the mesh and any associated sections.
3242: Not Collective
3244: Input Parameter:
3245: + serialMesh - The original DMMesh object
3246: - partitioner - The partitioning package, or NULL for the default
3248: Output Parameter:
3249: . parallelMesh - The distributed DMMesh object
3251: Level: intermediate
3253: .keywords: mesh, elements
3255: .seealso: DMMeshCreate(), DMMeshDistribute()
3256: @*/
3257: PetscErrorCode DMMeshDistributeByFace(DM serialMesh, const char partitioner[], DM *parallelMesh)
3258: {
3259: ALE::Obj<PETSC_MESH_TYPE> oldMesh;
3260: PetscErrorCode ierr;
3263: DMMeshGetMesh(serialMesh, oldMesh);
3264: DMMeshCreate(oldMesh->comm(), parallelMesh);
3265: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP, "I am being lazy, bug me.");
3266: #if 0
3267: ALE::DistributionNew<PETSC_MESH_TYPE>::distributeMeshAndSectionsV(oldMesh, newMesh, height = 1);
3268: #endif
3269: return(0);
3270: }
3272: #if defined(PETSC_HAVE_TRIANGLE)
3273: /* Already included since C++ is turned on #include <triangle.h> */
3277: PetscErrorCode TriangleInitInput(struct triangulateio *inputCtx)
3278: {
3280: inputCtx->numberofpoints = 0;
3281: inputCtx->numberofpointattributes = 0;
3282: inputCtx->pointlist = NULL;
3283: inputCtx->pointattributelist = NULL;
3284: inputCtx->pointmarkerlist = NULL;
3285: inputCtx->numberofsegments = 0;
3286: inputCtx->segmentlist = NULL;
3287: inputCtx->segmentmarkerlist = NULL;
3288: inputCtx->numberoftriangleattributes = 0;
3289: inputCtx->trianglelist = NULL;
3290: inputCtx->numberofholes = 0;
3291: inputCtx->holelist = NULL;
3292: inputCtx->numberofregions = 0;
3293: inputCtx->regionlist = NULL;
3294: return(0);
3295: }
3299: PetscErrorCode TriangleInitOutput(struct triangulateio *outputCtx)
3300: {
3302: outputCtx->numberofpoints = 0;
3303: outputCtx->pointlist = NULL;
3304: outputCtx->pointattributelist = NULL;
3305: outputCtx->pointmarkerlist = NULL;
3306: outputCtx->numberoftriangles = 0;
3307: outputCtx->trianglelist = NULL;
3308: outputCtx->triangleattributelist = NULL;
3309: outputCtx->neighborlist = NULL;
3310: outputCtx->segmentlist = NULL;
3311: outputCtx->segmentmarkerlist = NULL;
3312: outputCtx->edgelist = NULL;
3313: outputCtx->edgemarkerlist = NULL;
3314: return(0);
3315: }
3319: PetscErrorCode TriangleFiniOutput(struct triangulateio *outputCtx)
3320: {
3322: free(outputCtx->pointmarkerlist);
3323: free(outputCtx->edgelist);
3324: free(outputCtx->edgemarkerlist);
3325: free(outputCtx->trianglelist);
3326: free(outputCtx->neighborlist);
3327: return(0);
3328: }
3332: PetscErrorCode DMMeshGenerate_Triangle(DM boundary, PetscBool interpolate, DM *dm)
3333: {
3334: MPI_Comm comm = ((PetscObject) boundary)->comm;
3335: DM_Mesh *bd = (DM_Mesh*) boundary->data;
3336: PetscInt dim = 2;
3337: const PetscBool createConvexHull = PETSC_FALSE;
3338: const PetscBool constrained = PETSC_FALSE;
3339: struct triangulateio in;
3340: struct triangulateio out;
3341: PetscInt vStart, vEnd, v, eStart, eEnd, e;
3342: PetscMPIInt rank;
3343: PetscErrorCode ierr;
3346: MPI_Comm_rank(comm, &rank);
3347: TriangleInitInput(&in);
3348: TriangleInitOutput(&out);
3349: DMMeshGetDepthStratum(boundary, 0, &vStart, &vEnd);
3351: in.numberofpoints = vEnd - vStart;
3352: if (in.numberofpoints > 0) {
3353: PetscScalar *array;
3355: PetscMalloc(in.numberofpoints*dim * sizeof(double), &in.pointlist);
3356: PetscMalloc(in.numberofpoints * sizeof(int), &in.pointmarkerlist);
3357: VecGetArray(bd->coordinates, &array);
3358: for (v = vStart; v < vEnd; ++v) {
3359: const PetscInt idx = v - vStart;
3360: PetscInt off, d;
3362: PetscSectionGetOffset(bd->coordSection, v, &off);
3363: for (d = 0; d < dim; ++d) {
3364: in.pointlist[idx*dim + d] = array[off+d];
3365: }
3366: DMMeshGetLabelValue(boundary, "marker", v, &in.pointmarkerlist[idx]);
3367: }
3368: VecRestoreArray(bd->coordinates, &array);
3369: }
3370: DMMeshGetHeightStratum(boundary, 0, &eStart, &eEnd);
3372: in.numberofsegments = eEnd - eStart;
3373: if (in.numberofsegments > 0) {
3374: PetscMalloc(in.numberofsegments*2 * sizeof(int), &in.segmentlist);
3375: PetscMalloc(in.numberofsegments * sizeof(int), &in.segmentmarkerlist);
3376: for (e = eStart; e < eEnd; ++e) {
3377: const PetscInt idx = e - eStart;
3378: const PetscInt *cone;
3380: DMMeshGetCone(boundary, e, &cone);
3382: in.segmentlist[idx*2+0] = cone[0] - vStart;
3383: in.segmentlist[idx*2+1] = cone[1] - vStart;
3385: DMMeshGetLabelValue(boundary, "marker", e, &in.segmentmarkerlist[idx]);
3386: }
3387: }
3388: #if 0 /* Do not currently support holes */
3389: PetscReal *holeCoords;
3390: PetscInt h, d;
3392: DMMeshGetHoles(boundary, &in.numberofholes, &holeCords);
3393: if (in.numberofholes > 0) {
3394: PetscMalloc(in.numberofholes*dim * sizeof(double), &in.holelist);CHKERRXX(ierr);
3395: for (h = 0; h < in.numberofholes; ++h) {
3396: for (d = 0; d < dim; ++d) {
3397: in.holelist[h*dim+d] = holeCoords[h*dim+d];
3398: }
3399: }
3400: }
3401: #endif
3402: if (!rank) {
3403: char args[32];
3405: /* Take away 'Q' for verbose output */
3406: PetscStrcpy(args, "pqezQ");
3407: if (createConvexHull) {
3408: PetscStrcat(args, "c");
3409: }
3410: if (constrained) {
3411: PetscStrcpy(args, "zepDQ");
3412: }
3413: triangulate(args, &in, &out, NULL);
3414: }
3415: PetscFree(in.pointlist);
3416: PetscFree(in.pointmarkerlist);
3417: PetscFree(in.segmentlist);
3418: PetscFree(in.segmentmarkerlist);
3419: PetscFree(in.holelist);
3421: DMCreate(comm, dm);
3422: DMSetType(*dm, DMMESH);
3423: DMMeshSetDimension(*dm, dim);
3424: {
3425: DM_Mesh *mesh = (DM_Mesh*) (*dm)->data;
3426: const PetscInt numCorners = 3;
3427: const PetscInt numCells = out.numberoftriangles;
3428: const PetscInt numVertices = out.numberofpoints;
3429: int *cells = out.trianglelist;
3430: double *meshCoords = out.pointlist;
3431: PetscInt coordSize, c;
3432: PetscScalar *coords;
3434: DMMeshSetChart(*dm, 0, numCells+numVertices);
3435: for (c = 0; c < numCells; ++c) {
3436: DMMeshSetConeSize(*dm, c, numCorners);
3437: }
3438: DMMeshSetUp(*dm);
3439: for (c = 0; c < numCells; ++c) {
3440: PetscInt cone[numCorners] = {cells[c*numCorners+0]+numCells, cells[c*numCorners+1]+numCells, cells[c*numCorners+2]+numCells};
3442: DMMeshSetCone(*dm, c, cone);
3443: }
3444: DMMeshSymmetrize(*dm);
3445: DMMeshStratify(*dm);
3446: PetscSectionSetChart(mesh->coordSection, numCells, numCells + numVertices);
3447: for (v = numCells; v < numCells+numVertices; ++v) {
3448: PetscSectionSetDof(mesh->coordSection, v, dim);
3449: }
3450: PetscSectionSetUp(mesh->coordSection);
3451: PetscSectionGetStorageSize(mesh->coordSection, &coordSize);
3452: VecSetSizes(mesh->coordinates, coordSize, PETSC_DETERMINE);
3453: VecSetFromOptions(mesh->coordinates);
3454: VecGetArray(mesh->coordinates, &coords);
3455: for (v = 0; v < numVertices; ++v) {
3456: coords[v*dim+0] = meshCoords[v*dim+0];
3457: coords[v*dim+1] = meshCoords[v*dim+1];
3458: }
3459: VecRestoreArray(mesh->coordinates, &coords);
3460: for (v = 0; v < numVertices; ++v) {
3461: if (out.pointmarkerlist[v]) {
3462: DMMeshSetLabelValue(*dm, "marker", v+numCells, out.pointmarkerlist[v]);
3463: }
3464: }
3465: if (interpolate) {
3466: for (e = 0; e < out.numberofedges; e++) {
3467: if (out.edgemarkerlist[e]) {
3468: const PetscInt vertices[2] = {out.edgelist[e*2+0]+numCells, out.edgelist[e*2+1]+numCells};
3469: PetscInt edge;
3471: DMMeshJoinPoints(*dm, vertices, &edge); /* 1-level join */
3472: DMMeshSetLabelValue(*dm, "marker", edge, out.edgemarkerlist[e]);
3473: }
3474: }
3475: }
3476: }
3477: #if 0 /* Do not currently support holes */
3478: DMMeshCopyHoles(*dm, boundary);
3479: #endif
3480: TriangleFiniOutput(&out);
3481: return(0);
3482: }
3483: #endif
3487: /*@C
3488: DMMeshGenerate - Generates a mesh.
3490: Not Collective
3492: Input Parameters:
3493: + boundary - The DMMesh boundary object
3494: - interpolate - Flag to create intermediate mesh elements
3496: Output Parameter:
3497: . mesh - The DMMesh object
3499: Level: intermediate
3501: .keywords: mesh, elements
3502: .seealso: DMMeshCreate(), DMMeshRefine()
3503: @*/
3504: PetscErrorCode DMMeshGenerate(DM boundary, PetscBool interpolate, DM *mesh)
3505: {
3506: DM_Mesh *bd = (DM_Mesh*) boundary->data;
3512: if (bd->useNewImpl) {
3513: if (interpolate) SETERRQ(((PetscObject) boundary)->comm, PETSC_ERR_SUP, "Interpolation (creation of faces and edges) is not yet supported.");
3514: #if defined(PETSC_HAVE_TRIANGLE)
3515: DMMeshGenerate_Triangle(boundary, interpolate, mesh);
3516: #endif
3517: } else {
3518: ALE::Obj<PETSC_MESH_TYPE> mB;
3520: DMMeshGetMesh(boundary, mB);
3521: DMMeshCreate(mB->comm(), mesh);
3522: ALE::Obj<PETSC_MESH_TYPE> m = ALE::Generator<PETSC_MESH_TYPE>::generateMeshV(mB, interpolate, false, true);
3523: DMMeshSetMesh(*mesh, m);
3524: }
3525: return(0);
3526: }
3530: /*@C
3531: DMMeshRefine - Refines the mesh.
3533: Not Collective
3535: Input Parameters:
3536: + mesh - The original DMMesh object
3537: . refinementLimit - The maximum size of any cell
3538: - interpolate - Flag to create intermediate mesh elements
3540: Output Parameter:
3541: . refinedMesh - The refined DMMesh object
3543: Level: intermediate
3545: .keywords: mesh, elements
3546: .seealso: DMMeshCreate(), DMMeshGenerate()
3547: @*/
3548: PetscErrorCode DMMeshRefine(DM mesh, double refinementLimit, PetscBool interpolate, DM *refinedMesh)
3549: {
3550: ALE::Obj<PETSC_MESH_TYPE> oldMesh;
3551: PetscErrorCode ierr;
3554: if (refinementLimit == 0.0) return(0);
3555: DMMeshGetMesh(mesh, oldMesh);
3556: DMMeshCreate(oldMesh->comm(), refinedMesh);
3557: ALE::Obj<PETSC_MESH_TYPE> newMesh = ALE::Generator<PETSC_MESH_TYPE>::refineMeshV(oldMesh, refinementLimit, interpolate, false, true);
3558: DMMeshSetMesh(*refinedMesh, newMesh);
3559: return(0);
3560: }
3564: PetscErrorCode DMRefine_Mesh(DM dm, MPI_Comm comm, DM *dmRefined)
3565: {
3566: ALE::Obj<PETSC_MESH_TYPE> oldMesh;
3567: double refinementLimit;
3568: PetscErrorCode ierr;
3571: DMMeshGetMesh(dm, oldMesh);
3572: DMMeshCreate(comm, dmRefined);
3574: refinementLimit = oldMesh->getMaxVolume()/2.0;
3575: ALE::Obj<PETSC_MESH_TYPE> newMesh = ALE::Generator<PETSC_MESH_TYPE>::refineMeshV(oldMesh, refinementLimit, true);
3576: DMMeshSetMesh(*dmRefined, newMesh);
3577: return(0);
3578: }
3582: PetscErrorCode DMCoarsenHierarchy_Mesh(DM mesh, int numLevels, DM *coarseHierarchy)
3583: {
3584: PetscReal cfactor = 1.5;
3588: PetscOptionsReal("-mesh_coarsen_factor", "The coarsening factor", NULL, cfactor, &cfactor, NULL);
3589: SETERRQ(((PetscObject) mesh)->comm, PETSC_ERR_SUP, "Peter needs to incorporate his code.");
3590: return(0);
3591: }
3595: PetscErrorCode DMCreateInterpolation_Mesh(DM dmCoarse, DM dmFine, Mat *interpolation, Vec *scaling)
3596: {
3597: SETERRQ(((PetscObject) dmCoarse)->comm, PETSC_ERR_SUP, "Peter needs to incorporate his code.");
3598: }
3602: PetscErrorCode DMMeshMarkBoundaryCells(DM dm, const char labelName[], PetscInt marker, PetscInt newMarker)
3603: {
3604: ALE::Obj<PETSC_MESH_TYPE> mesh;
3605: PetscErrorCode ierr;
3608: DMMeshGetMesh(dm, mesh);
3609: mesh->markBoundaryCells(labelName, marker, newMarker);
3610: return(0);
3611: }
3615: PetscErrorCode DMMeshGetDepthStratum(DM dm, PetscInt stratumValue, PetscInt *start, PetscInt *end)
3616: {
3617: DM_Mesh *mesh = (DM_Mesh*) dm->data;
3622: if (mesh->useNewImpl) {
3623: SieveLabel next = mesh->labels;
3624: PetscBool flg = PETSC_FALSE;
3625: PetscInt depth;
3627: if (stratumValue < 0) {
3628: DMMeshGetChart(dm, start, end);
3629: return(0);
3630: } else {
3631: PetscInt pStart, pEnd;
3633: if (start) *start = 0;
3634: if (end) *end = 0;
3635: DMMeshGetChart(dm, &pStart, &pEnd);
3636: if (pStart == pEnd) return(0);
3637: }
3638: while (next) {
3639: PetscStrcmp("depth", next->name, &flg);
3640: if (flg) break;
3641: next = next->next;
3642: }
3643: if (!flg) SETERRQ(((PetscObject) dm)->comm, PETSC_ERR_ARG_WRONG, "No label named depth was found");
3644: /* Strata are sorted and contiguous -- In addition, depth/height is either full or 1-level */
3645: depth = stratumValue;
3646: if ((depth < 0) || (depth >= next->numStrata)) {
3647: if (start) *start = 0;
3648: if (end) *end = 0;
3649: } else {
3650: if (start) *start = next->points[next->stratumOffsets[depth]];
3651: if (end) *end = next->points[next->stratumOffsets[depth]+next->stratumSizes[depth]-1]+1;
3652: }
3653: } else {
3654: ALE::Obj<PETSC_MESH_TYPE> mesh;
3655: DMMeshGetMesh(dm, mesh);
3656: if (stratumValue < 0) {
3657: if (start) *start = mesh->getSieve()->getChart().min();
3658: if (end) *end = mesh->getSieve()->getChart().max();
3659: } else {
3660: const Obj<PETSC_MESH_TYPE::label_sequence>& stratum = mesh->depthStratum(stratumValue);
3661: if (start) *start = *stratum->begin();
3662: if (end) *end = *stratum->rbegin()+1;
3663: }
3664: }
3665: return(0);
3666: }
3670: PetscErrorCode DMMeshGetHeightStratum(DM dm, PetscInt stratumValue, PetscInt *start, PetscInt *end)
3671: {
3672: DM_Mesh *mesh = (DM_Mesh*) dm->data;
3677: if (mesh->useNewImpl) {
3678: SieveLabel next = mesh->labels;
3679: PetscBool flg = PETSC_FALSE;
3680: PetscInt depth;
3682: if (stratumValue < 0) {
3683: DMMeshGetChart(dm, start, end);
3684: } else {
3685: PetscInt pStart, pEnd;
3687: if (start) *start = 0;
3688: if (end) *end = 0;
3689: DMMeshGetChart(dm, &pStart, &pEnd);
3690: if (pStart == pEnd) return(0);
3691: }
3692: while (next) {
3693: PetscStrcmp("depth", next->name, &flg);
3694: if (flg) break;
3695: next = next->next;
3696: }
3697: if (!flg) SETERRQ(((PetscObject) dm)->comm, PETSC_ERR_ARG_WRONG, "No label named depth was found");
3698: /* Strata are sorted and contiguous -- In addition, depth/height is either full or 1-level */
3699: depth = next->stratumValues[next->numStrata-1] - stratumValue;
3700: if ((depth < 0) || (depth >= next->numStrata)) {
3701: if (start) *start = 0;
3702: if (end) *end = 0;
3703: } else {
3704: if (start) *start = next->points[next->stratumOffsets[depth]];
3705: if (end) *end = next->points[next->stratumOffsets[depth]+next->stratumSizes[depth]-1]+1;
3706: }
3707: } else {
3708: ALE::Obj<PETSC_MESH_TYPE> mesh;
3709: DMMeshGetMesh(dm, mesh);
3710: if (mesh->getLabel("height")->size()) {
3711: const Obj<PETSC_MESH_TYPE::label_sequence>& stratum = mesh->heightStratum(stratumValue);
3712: if (start) *start = *stratum->begin();
3713: if (end) *end = *stratum->rbegin()+1;
3714: } else {
3715: if (start) *start = 0;
3716: if (end) *end = 0;
3717: }
3718: }
3719: return(0);
3720: }
3724: PetscErrorCode DMMeshCreateSection(DM dm, PetscInt dim, PetscInt numFields, PetscInt numComp[], PetscInt numDof[], PetscInt numBC, PetscInt bcField[], IS bcPoints[], PetscSection *section)
3725: {
3726: PetscInt *numDofTot, *maxConstraints;
3727: PetscInt pStart = 0, pEnd = 0;
3731: PetscMalloc2(dim+1,PetscInt,&numDofTot,numFields+1,PetscInt,&maxConstraints);
3732: for (PetscInt d = 0; d <= dim; ++d) {
3733: numDofTot[d] = 0;
3734: for (PetscInt f = 0; f < numFields; ++f) numDofTot[d] += numDof[f*(dim+1)+d];
3735: }
3736: PetscSectionCreate(((PetscObject) dm)->comm, section);
3737: if (numFields > 1) {
3738: PetscSectionSetNumFields(*section, numFields);
3739: if (numComp) {
3740: for (PetscInt f = 0; f < numFields; ++f) {
3741: PetscSectionSetFieldComponents(*section, f, numComp[f]);
3742: }
3743: }
3744: } else {
3745: numFields = 0;
3746: }
3747: DMMeshGetChart(dm, &pStart, &pEnd);
3748: PetscSectionSetChart(*section, pStart, pEnd);
3749: for (PetscInt d = 0; d <= dim; ++d) {
3750: DMMeshGetDepthStratum(dm, d, &pStart, &pEnd);
3751: for (PetscInt p = pStart; p < pEnd; ++p) {
3752: for (PetscInt f = 0; f < numFields; ++f) {
3753: PetscSectionSetFieldDof(*section, p, f, numDof[f*(dim+1)+d]);
3754: }
3755: PetscSectionSetDof(*section, p, numDofTot[d]);
3756: }
3757: }
3758: if (numBC) {
3759: for (PetscInt f = 0; f <= numFields; ++f) maxConstraints[f] = 0;
3760: for (PetscInt bc = 0; bc < numBC; ++bc) {
3761: PetscInt field = 0;
3762: const PetscInt *idx;
3763: PetscInt n;
3765: if (numFields) field = bcField[bc];
3766: ISGetLocalSize(bcPoints[bc], &n);
3767: ISGetIndices(bcPoints[bc], &idx);
3768: for (PetscInt i = 0; i < n; ++i) {
3769: const PetscInt p = idx[i];
3770: PetscInt depth, numConst;
3772: DMMeshGetLabelValue(dm, "depth", p, &depth);
3773: numConst = numDof[field*(dim+1)+depth];
3774: maxConstraints[field] = PetscMax(maxConstraints[field], numConst);
3775: if (numFields) {
3776: PetscSectionSetFieldConstraintDof(*section, p, field, numConst);
3777: }
3778: PetscSectionAddConstraintDof(*section, p, numConst);
3779: }
3780: ISRestoreIndices(bcPoints[bc], &idx);
3781: }
3782: for (PetscInt f = 0; f < numFields; ++f) {
3783: maxConstraints[numFields] += maxConstraints[f];
3784: }
3785: }
3786: PetscSectionSetUp(*section);
3787: if (maxConstraints[numFields]) {
3788: PetscInt *indices;
3790: PetscMalloc(maxConstraints[numFields] * sizeof(PetscInt), &indices);
3791: PetscSectionGetChart(*section, &pStart, &pEnd);
3792: for (PetscInt p = pStart; p < pEnd; ++p) {
3793: PetscInt cDof;
3795: PetscSectionGetConstraintDof(*section, p, &cDof);
3796: if (cDof) {
3797: if (cDof > maxConstraints[numFields]) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_LIB, "Likely memory corruption, point %d cDof %d > maxConstraints %d", p, cDof, maxConstraints);
3798: if (numFields) {
3799: PetscInt numConst = 0, fOff = 0;
3801: for (PetscInt f = 0; f < numFields; ++f) {
3802: PetscInt cfDof, fDof;
3804: PetscSectionGetFieldDof(*section, p, f, &fDof);
3805: PetscSectionGetFieldConstraintDof(*section, p, f, &cfDof);
3806: for (PetscInt d = 0; d < cfDof; ++d) {
3807: indices[numConst+d] = fOff+d;
3808: }
3809: PetscSectionSetFieldConstraintIndices(*section, p, f, &indices[numConst]);
3810: numConst += cfDof;
3811: fOff += fDof;
3812: }
3813: if (cDof != numConst) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_LIB, "Total number of field constraints %d should be %d", numConst, cDof);
3814: } else {
3815: for (PetscInt d = 0; d < cDof; ++d) indices[d] = d;
3816: }
3817: PetscSectionSetConstraintIndices(*section, p, indices);
3818: }
3819: }
3820: PetscFree(indices);
3821: }
3822: PetscFree2(numDofTot,maxConstraints);
3823: {
3824: PetscBool view = PETSC_FALSE;
3826: PetscOptionsHasName(((PetscObject) dm)->prefix, "-section_view", &view);
3827: if (view) {
3828: PetscSectionView(*section, PETSC_VIEWER_STDOUT_WORLD);
3829: }
3830: }
3831: return(0);
3832: }
3836: PetscErrorCode DMMeshConvertSection(const ALE::Obj<PETSC_MESH_TYPE>& mesh, const Obj<PETSC_MESH_TYPE::real_section_type>& s, PetscSection *section)
3837: {
3838: const PetscInt pStart = s->getChart().min();
3839: const PetscInt pEnd = s->getChart().max();
3840: PetscInt numFields = s->getNumSpaces();
3844: PetscSectionCreate(mesh->comm(), section);
3845: if (numFields) {
3846: PetscSectionSetNumFields(*section, numFields);
3847: for (PetscInt f = 0; f < numFields; ++f) {
3848: PetscSectionSetFieldComponents(*section, f, s->getSpaceComponents(f));
3849: }
3850: }
3851: PetscSectionSetChart(*section, pStart, pEnd);
3852: for (PetscInt p = pStart; p < pEnd; ++p) {
3853: PetscSectionSetDof(*section, p, s->getFiberDimension(p));
3854: for (PetscInt f = 0; f < numFields; ++f) {
3855: PetscSectionSetFieldDof(*section, p, f, s->getFiberDimension(p, f));
3856: }
3857: PetscSectionSetConstraintDof(*section, p, s->getConstraintDimension(p));
3858: for (PetscInt f = 0; f < numFields; ++f) {
3859: PetscSectionSetFieldConstraintDof(*section, p, f, s->getConstraintDimension(p, f));
3860: }
3861: }
3862: PetscSectionSetUp(*section);
3863: for (PetscInt p = pStart; p < pEnd; ++p) {
3864: PetscSectionSetConstraintIndices(*section, p, (PetscInt*) s->getConstraintDof(p));
3865: for (PetscInt f = 0; f < numFields; ++f) {
3866: PetscSectionSetFieldConstraintIndices(*section, p, f, (PetscInt*) s->getConstraintDof(p, f));
3867: }
3868: }
3869: return(0);
3870: }
3874: PetscErrorCode DMMeshGetSection(DM dm, const char name[], PetscSection *section)
3875: {
3876: ALE::Obj<PETSC_MESH_TYPE> mesh;
3877: PetscErrorCode ierr;
3880: DMMeshGetMesh(dm, mesh);
3881: DMMeshConvertSection(mesh, mesh->getRealSection(name), section);
3882: return(0);
3883: }
3887: PetscErrorCode DMMeshSetSection(DM dm, const char name[], PetscSection section)
3888: {
3889: ALE::Obj<PETSC_MESH_TYPE> mesh;
3890: PetscErrorCode ierr;
3893: DMMeshGetMesh(dm, mesh);
3894: {
3895: const Obj<PETSC_MESH_TYPE::real_section_type>& s = mesh->getRealSection(name);
3896: PetscInt pStart, pEnd, numFields;
3898: PetscSectionGetChart(section, &pStart, &pEnd);
3899: s->setChart(PETSC_MESH_TYPE::real_section_type::chart_type(pStart, pEnd));
3900: PetscSectionGetNumFields(section, &numFields);
3901: for (PetscInt f = 0; f < numFields; ++f) {
3902: PetscInt comp;
3903: PetscSectionGetFieldComponents(section, f, &comp);
3904: s->addSpace(comp);
3905: }
3906: for (PetscInt p = pStart; p < pEnd; ++p) {
3907: PetscInt fDim, cDim;
3909: PetscSectionGetDof(section, p, &fDim);
3910: s->setFiberDimension(p, fDim);
3911: for (PetscInt f = 0; f < numFields; ++f) {
3912: PetscSectionGetFieldDof(section, p, f, &fDim);
3913: s->setFiberDimension(p, fDim, f);
3914: }
3915: PetscSectionGetConstraintDof(section, p, &cDim);
3916: if (cDim) {
3917: s->setConstraintDimension(p, cDim);
3918: for (PetscInt f = 0; f < numFields; ++f) {
3919: PetscSectionGetFieldConstraintDof(section, p, f, &cDim);
3920: s->setConstraintDimension(p, cDim, f);
3921: }
3922: }
3923: }
3924: s->allocatePoint();
3925: for (PetscInt p = pStart; p < pEnd; ++p) {
3926: const PetscInt *indices;
3928: PetscSectionGetConstraintIndices(section, p, &indices);
3929: s->setConstraintDof(p, indices);
3930: for (PetscInt f = 0; f < numFields; ++f) {
3931: PetscSectionGetFieldConstraintIndices(section, p, f, &indices);
3932: s->setConstraintDof(p, indices, f);
3933: }
3934: }
3935: {
3936: PetscBool isDefault;
3938: PetscStrcmp(name, "default", &isDefault);
3939: if (isDefault) {
3940: PetscInt maxDof = 0;
3942: for (PetscInt p = pStart; p < pEnd; ++p) {
3943: PetscInt fDim;
3945: PetscSectionGetDof(section, p, &fDim);
3946: maxDof = PetscMax(maxDof, fDim);
3947: }
3948: mesh->setMaxDof(maxDof);
3949: }
3950: }
3951: }
3952: return(0);
3953: }
3957: /*
3958: Note: This gets a borrowed reference, so the user should not destroy this PetscSection.
3959: */
3960: PetscErrorCode DMMeshGetDefaultSection(DM dm, PetscSection *section)
3961: {
3962: DM_Mesh *mesh = (DM_Mesh*) dm->data;
3966: if (!mesh->defaultSection && !mesh->useNewImpl) {
3967: DMMeshGetSection(dm, "default", &mesh->defaultSection);
3968: }
3969: *section = mesh->defaultSection;
3970: return(0);
3971: }
3975: /*
3976: Note: This reference will be stolen, so the user should not destroy this PetscSection.
3977: */
3978: PetscErrorCode DMMeshSetDefaultSection(DM dm, PetscSection section)
3979: {
3980: DM_Mesh *mesh = (DM_Mesh*) dm->data;
3984: mesh->defaultSection = section;
3985: if (!mesh->useNewImpl) {
3986: DMMeshSetSection(dm, "default", mesh->defaultSection);
3987: }
3988: return(0);
3989: }
3993: PetscErrorCode DMMeshGetCoordinateSection(DM dm, PetscSection *section)
3994: {
3995: DM_Mesh *mesh = (DM_Mesh*) dm->data;
4001: if (mesh->useNewImpl) *section = mesh->coordSection;
4002: else {
4003: DMMeshGetSection(dm, "coordinates", section);
4004: }
4005: return(0);
4006: }
4010: PetscErrorCode DMMeshSetCoordinateSection(DM dm, PetscSection section)
4011: {
4015: DMMeshSetSection(dm, "coordinates", section);
4016: return(0);
4017: }
4021: PetscErrorCode DMMeshGetConeSection(DM dm, PetscSection *section)
4022: {
4023: DM_Mesh *mesh = (DM_Mesh*) dm->data;
4027: if (!mesh->useNewImpl) SETERRQ(((PetscObject) dm)->comm, PETSC_ERR_ARG_WRONG, "This method is only valid for C implementation meshes.");
4028: if (section) *section = mesh->coneSection;
4029: return(0);
4030: }
4034: PetscErrorCode DMMeshGetCones(DM dm, PetscInt *cones[])
4035: {
4036: DM_Mesh *mesh = (DM_Mesh*) dm->data;
4040: if (!mesh->useNewImpl) SETERRQ(((PetscObject) dm)->comm, PETSC_ERR_ARG_WRONG, "This method is only valid for C implementation meshes.");
4041: if (cones) *cones = mesh->cones;
4042: return(0);
4043: }
4047: PetscErrorCode DMMeshCreateConeSection(DM dm, PetscSection *section)
4048: {
4049: ALE::Obj<PETSC_MESH_TYPE> mesh;
4050: PetscInt p;
4051: PetscErrorCode ierr;
4054: DMMeshGetMesh(dm, mesh);
4055: PetscSectionCreate(((PetscObject) dm)->comm, section);
4056: PetscSectionSetChart(*section, mesh->getSieve()->getChart().min(), mesh->getSieve()->getChart().max());
4057: for (p = mesh->getSieve()->getChart().min(); p < mesh->getSieve()->getChart().max(); ++p) {
4058: PetscSectionSetDof(*section, p, mesh->getSieve()->getConeSize(p));
4059: }
4060: PetscSectionSetUp(*section);
4061: return(0);
4062: }
4066: PetscErrorCode DMMeshGetCoordinateVec(DM dm, Vec *coordinates)
4067: {
4068: DM_Mesh *mesh = (DM_Mesh*) dm->data;
4074: if (mesh->useNewImpl) *coordinates = mesh->coordinates;
4075: else {
4076: ALE::Obj<PETSC_MESH_TYPE> mesh;
4077: DMMeshGetMesh(dm, mesh);
4078: const Obj<PETSC_MESH_TYPE::real_section_type>& coords = mesh->getRealSection("coordinates");
4079: VecCreateSeqWithArray(PETSC_COMM_SELF,1, coords->getStorageSize(), coords->restrictSpace(), coordinates);
4080: }
4081: return(0);
4082: }
4086: PetscErrorCode DMMeshComputeCellGeometry(DM dm, PetscInt cell, PetscReal *v0, PetscReal *J, PetscReal *invJ, PetscReal *detJ)
4087: {
4088: ALE::Obj<PETSC_MESH_TYPE> mesh;
4089: PetscErrorCode ierr;
4092: DMMeshGetMesh(dm, mesh);
4093: {
4094: ALE::Obj<PETSC_MESH_TYPE::real_section_type> coordinates = mesh->getRealSection("coordinates");
4096: mesh->computeElementGeometry(coordinates, cell, v0, J, invJ, *detJ);
4097: }
4098: return(0);
4099: }
4103: PetscErrorCode DMMeshVecGetClosure(DM dm, Vec v, PetscInt point, const PetscScalar *values[])
4104: {
4105: ALE::Obj<PETSC_MESH_TYPE> mesh;
4106: PetscErrorCode ierr;
4109: #if defined(PETSC_USE_COMPLEX)
4110: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP, "DMMesh does not support complex closure");
4111: #else
4112: DMMeshGetMesh(dm, mesh);
4113: /* Peeling back IMesh::restrictClosure() */
4114: try {
4115: typedef ALE::ISieveVisitor::RestrictVecVisitor<PetscScalar> visitor_type;
4116: PetscSection section;
4117: PetscScalar *array;
4118: PetscInt numFields;
4120: DMMeshGetDefaultSection(dm, §ion);
4121: PetscSectionGetNumFields(section, &numFields);
4122: const PetscInt size = mesh->sizeWithBC(section, point); /* OPT: This can be precomputed */
4123: DMGetWorkArray(dm, 2*size+numFields+1, PETSC_SCALAR, &array);
4124: visitor_type rV(v, section, size, array, (PetscInt*) &array[2*size], (PetscInt*) &array[size]);
4125: if (mesh->depth() == 1) {
4126: rV.visitPoint(point, 0);
4127: // Cone is guarateed to be ordered correctly
4128: mesh->getSieve()->orientedCone(point, rV);
4129: } else {
4130: ALE::ISieveVisitor::PointRetriever<PETSC_MESH_TYPE::sieve_type,visitor_type> pV((int) pow((double) mesh->getSieve()->getMaxConeSize(), mesh->depth())+1, rV, true);
4132: ALE::ISieveTraversal<PETSC_MESH_TYPE::sieve_type>::orientedClosure(*mesh->getSieve(), point, pV);
4133: }
4134: *values = rV.getValues();
4135: } catch(ALE::Exception e) {
4136: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid argument: %s", e.message());
4137: } catch(PETSc::Exception e) {
4138: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid argument: %s", e.message());
4139: }
4140: #endif
4141: return(0);
4142: }
4146: PetscErrorCode DMMeshVecRestoreClosure(DM dm, Vec v, PetscInt point, const PetscScalar *values[])
4147: {
4151: DMRestoreWorkArray(dm, 0, PETSC_SCALAR, values);
4152: return(0);
4153: }
4157: PetscErrorCode DMMeshVecSetClosure(DM dm, Vec v, PetscInt point, const PetscScalar values[], InsertMode mode)
4158: {
4159: ALE::Obj<PETSC_MESH_TYPE> mesh;
4160: PetscErrorCode ierr;
4163: #if defined(PETSC_USE_COMPLEX)
4164: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP, "DMMesh does not support complex closure");
4165: #else
4166: DMMeshGetMesh(dm, mesh);
4167: /* Peeling back IMesh::update() and IMesh::updateAdd() */
4168: try {
4169: typedef ALE::ISieveVisitor::UpdateVecVisitor<PetscScalar> visitor_type;
4170: PetscSection section;
4171: PetscInt *fieldSize;
4172: PetscInt numFields;
4174: DMMeshGetDefaultSection(dm, §ion);
4175: PetscSectionGetNumFields(section, &numFields);
4176: DMGetWorkArray(dm, numFields, PETSC_INT, &fieldSize);
4177: mesh->sizeWithBC(section, point, fieldSize); /* OPT: This can be precomputed */
4178: visitor_type uV(v, section, values, mode, numFields, fieldSize);
4179: if (mesh->depth() == 1) {
4180: uV.visitPoint(point, 0);
4181: // Cone is guarateed to be ordered correctly
4182: mesh->getSieve()->orientedCone(point, uV);
4183: } else {
4184: ALE::ISieveVisitor::PointRetriever<PETSC_MESH_TYPE::sieve_type,visitor_type> pV((int) pow((double) mesh->getSieve()->getMaxConeSize(), mesh->depth())+1, uV, true);
4186: ALE::ISieveTraversal<PETSC_MESH_TYPE::sieve_type>::orientedClosure(*mesh->getSieve(), point, pV);
4187: }
4188: DMRestoreWorkArray(dm, numFields, PETSC_INT, &fieldSize);
4189: } catch(ALE::Exception e) {
4190: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid argument: %s", e.message());
4191: }
4192: #endif
4193: return(0);
4194: }
4198: PetscErrorCode DMMeshMatSetClosure(DM dm, Mat A, PetscInt point, PetscScalar values[], InsertMode mode)
4199: {
4200: ALE::Obj<PETSC_MESH_TYPE> mesh;
4201: PetscErrorCode ierr;
4204: #if defined(PETSC_USE_COMPLEX)
4205: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP, "DMMesh does not support complex closure");
4206: #else
4207: DMMeshGetMesh(dm, mesh);
4208: /* Copying from updateOperator() */
4209: try {
4210: typedef ALE::ISieveVisitor::IndicesVisitor<PetscSection,PETSC_MESH_TYPE::order_type,PetscInt> visitor_type;
4211: ALE::Obj<PETSC_MESH_TYPE::real_section_type> s = mesh->getRealSection("default");
4212: const ALE::Obj<PETSC_MESH_TYPE::order_type>& globalOrder = mesh->getFactory()->getGlobalOrder(mesh, s->getName(), s);
4213: PetscSection section;
4214: PetscInt numFields;
4215: PetscInt *fieldSize = NULL;
4217: DMMeshGetDefaultSection(dm, §ion);
4218: PetscSectionGetNumFields(section, &numFields);
4219: if (numFields) {
4220: DMGetWorkArray(dm, numFields, PETSC_INT, &fieldSize);
4221: mesh->sizeWithBC(section, point, fieldSize); /* OPT: This can be precomputed */
4222: }
4223: visitor_type iV(section, *globalOrder, (int) pow((double) mesh->getSieve()->getMaxConeSize(), mesh->depth())*mesh->getMaxDof(), mesh->depth() > 1, fieldSize);
4225: updateOperator(A, *mesh->getSieve(), iV, point, values, mode);
4226: if (numFields) {
4227: DMRestoreWorkArray(dm, numFields, PETSC_INT, &fieldSize);
4228: }
4229: } catch(ALE::Exception e) {
4230: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid argument: %s", e.message());
4231: }
4232: #endif
4233: return(0);
4234: }
4238: /*@C
4239: DMMeshHasSectionReal - Determines whether this mesh has a SectionReal with the given name.
4241: Not Collective
4243: Input Parameters:
4244: + mesh - The DMMesh object
4245: - name - The section name
4247: Output Parameter:
4248: . flag - True if the SectionReal is present in the DMMesh
4250: Level: intermediate
4252: .keywords: mesh, elements
4253: .seealso: DMMeshCreate()
4254: @*/
4255: PetscErrorCode DMMeshHasSectionReal(DM dm, const char name[], PetscBool *flag)
4256: {
4257: ALE::Obj<PETSC_MESH_TYPE> m;
4258: PetscErrorCode ierr;
4261: DMMeshGetMesh(dm, m);
4262: *flag = (PetscBool) m->hasRealSection(std::string(name));
4263: return(0);
4264: }
4268: /*@C
4269: DMMeshGetSectionReal - Returns a SectionReal of the given name from the DMMesh.
4271: Collective on DMMesh
4273: Input Parameters:
4274: + mesh - The DMMesh object
4275: - name - The section name
4277: Output Parameter:
4278: . section - The SectionReal
4280: Note: The section is a new object, and must be destroyed by the user
4282: Level: intermediate
4284: .keywords: mesh, elements
4286: .seealso: DMMeshCreate(), SectionRealDestroy()
4287: @*/
4288: PetscErrorCode DMMeshGetSectionReal(DM dm, const char name[], SectionReal *section)
4289: {
4290: ALE::Obj<PETSC_MESH_TYPE> m;
4291: bool has;
4292: PetscErrorCode ierr;
4295: DMMeshGetMesh(dm, m);
4296: SectionRealCreate(m->comm(), section);
4297: PetscObjectSetName((PetscObject) *section, name);
4298: has = m->hasRealSection(std::string(name));
4299: SectionRealSetSection(*section, m->getRealSection(std::string(name)));
4300: SectionRealSetBundle(*section, m);
4301: if (!has) {
4302: m->getRealSection(std::string(name))->setChart(m->getSieve()->getChart());
4303: }
4304: return(0);
4305: }
4309: /*@C
4310: DMMeshSetSectionReal - Puts a SectionReal of the given name into the DMMesh.
4312: Collective on DMMesh
4314: Input Parameters:
4315: + mesh - The DMMesh object
4316: - section - The SectionReal
4318: Note: This takes the section name from the PETSc object
4320: Level: intermediate
4322: .keywords: mesh, elements
4323: .seealso: DMMeshCreate()
4324: @*/
4325: PetscErrorCode DMMeshSetSectionReal(DM dm, const char name[], SectionReal section)
4326: {
4327: ALE::Obj<PETSC_MESH_TYPE> m;
4328: ALE::Obj<PETSC_MESH_TYPE::real_section_type> s;
4329: PetscErrorCode ierr;
4332: DMMeshGetMesh(dm, m);
4333: PetscObjectGetName((PetscObject) section, &name);
4334: SectionRealGetSection(section, s);
4335: m->setRealSection(std::string(name), s);
4336: return(0);
4337: }
4341: /*@C
4342: DMMeshHasSectionInt - Determines whether this mesh has a SectionInt with the given name.
4344: Not Collective
4346: Input Parameters:
4347: + mesh - The DMMesh object
4348: - name - The section name
4350: Output Parameter:
4351: . flag - True if the SectionInt is present in the DMMesh
4353: Level: intermediate
4355: .keywords: mesh, elements
4356: .seealso: DMMeshCreate()
4357: @*/
4358: PetscErrorCode DMMeshHasSectionInt(DM dm, const char name[], PetscBool *flag)
4359: {
4360: ALE::Obj<PETSC_MESH_TYPE> m;
4361: PetscErrorCode ierr;
4364: DMMeshGetMesh(dm, m);
4365: *flag = (PetscBool) m->hasIntSection(std::string(name));
4366: return(0);
4367: }
4371: /*@C
4372: DMMeshGetSectionInt - Returns a SectionInt of the given name from the DMMesh.
4374: Collective on DMMesh
4376: Input Parameters:
4377: + mesh - The DMMesh object
4378: - name - The section name
4380: Output Parameter:
4381: . section - The SectionInt
4383: Note: The section is a new object, and must be destroyed by the user
4385: Level: intermediate
4387: .keywords: mesh, elements
4388: .seealso: DMMeshCreate()
4389: @*/
4390: PetscErrorCode DMMeshGetSectionInt(DM dm, const char name[], SectionInt *section)
4391: {
4392: ALE::Obj<PETSC_MESH_TYPE> m;
4393: bool has;
4394: PetscErrorCode ierr;
4397: DMMeshGetMesh(dm, m);
4398: SectionIntCreate(m->comm(), section);
4399: PetscObjectSetName((PetscObject) *section, name);
4400: has = m->hasIntSection(std::string(name));
4401: SectionIntSetSection(*section, m->getIntSection(std::string(name)));
4402: SectionIntSetBundle(*section, m);
4403: if (!has) m->getIntSection(std::string(name))->setChart(m->getSieve()->getChart());
4404: return(0);
4405: }
4409: /*@C
4410: DMMeshSetSectionInt - Puts a SectionInt of the given name into the DMMesh.
4412: Collective on DMMesh
4414: Input Parameters:
4415: + mesh - The DMMesh object
4416: - section - The SectionInt
4418: Note: This takes the section name from the PETSc object
4420: Level: intermediate
4422: .keywords: mesh, elements
4423: .seealso: DMMeshCreate()
4424: @*/
4425: PetscErrorCode DMMeshSetSectionInt(DM dm, SectionInt section)
4426: {
4427: ALE::Obj<PETSC_MESH_TYPE> m;
4428: ALE::Obj<PETSC_MESH_TYPE::int_section_type> s;
4429: const char *name;
4430: PetscErrorCode ierr;
4433: DMMeshGetMesh(dm, m);
4434: PetscObjectGetName((PetscObject) section, &name);
4435: SectionIntGetSection(section, s);
4436: m->setIntSection(std::string(name), s);
4437: return(0);
4438: }
4442: /*@C
4443: SectionGetArray - Returns the array underlying the Section.
4445: Not Collective
4447: Input Parameters:
4448: + mesh - The DMMesh object
4449: - name - The section name
4451: Output Parameters:
4452: + numElements - The number of mesh element with values
4453: . fiberDim - The number of values per element
4454: - array - The array
4456: Level: intermediate
4458: .keywords: mesh, elements
4459: .seealso: DMMeshCreate()
4460: @*/
4461: PetscErrorCode SectionGetArray(DM dm, const char name[], PetscInt *numElements, PetscInt *fiberDim, PetscScalar *array[])
4462: {
4463: ALE::Obj<PETSC_MESH_TYPE> m;
4464: PetscErrorCode ierr;
4467: DMMeshGetMesh(dm, m);
4468: const Obj<PETSC_MESH_TYPE::real_section_type>& section = m->getRealSection(std::string(name));
4469: if (section->size() == 0) {
4470: *numElements = 0;
4471: *fiberDim = 0;
4472: *array = NULL;
4473: return(0);
4474: }
4475: const PETSC_MESH_TYPE::real_section_type::chart_type& chart = section->getChart();
4476: /* const int depth = m->depth(*chart.begin()); */
4477: /* *numElements = m->depthStratum(depth)->size(); */
4478: /* *fiberDim = section->getFiberDimension(*chart.begin()); */
4479: /* *array = (PetscScalar*) m->restrict(section); */
4480: int fiberDimMin = section->getFiberDimension(*chart.begin());
4481: int numElem = 0;
4483: for (PETSC_MESH_TYPE::real_section_type::chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
4484: const int fiberDim = section->getFiberDimension(*c_iter);
4486: if (fiberDim < fiberDimMin) fiberDimMin = fiberDim;
4487: }
4488: for (PETSC_MESH_TYPE::real_section_type::chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
4489: const int fiberDim = section->getFiberDimension(*c_iter);
4491: numElem += fiberDim/fiberDimMin;
4492: }
4493: *numElements = numElem;
4494: *fiberDim = fiberDimMin;
4495: *array = (PetscScalar*) section->restrictSpace();
4496: return(0);
4497: }
4501: inline void ExpandInterval(const ALE::Point& interval, int indices[], int& indx)
4502: {
4503: const int end = interval.prefix + interval.index;
4504: for (int i = interval.index; i < end; i++) indices[indx++] = i;
4505: }
4509: inline void ExpandInterval_New(ALE::Point interval, PetscInt indices[], PetscInt *indx)
4510: {
4511: for (int i = 0; i < interval.prefix; i++) indices[(*indx)++] = interval.index + i;
4512: for (int i = 0; i < -interval.prefix; i++) indices[(*indx)++] = -1;
4513: }
4517: /*@
4518: DMMeshClone - Creates a DMMesh object with the same mesh as the original.
4520: Collective on MPI_Comm
4522: Input Parameter:
4523: . dm - The original DMMesh object
4525: Output Parameter:
4526: . newdm - The new DMMesh object
4528: Level: beginner
4530: .keywords: DMMesh, create
4531: @*/
4532: PetscErrorCode DMMeshClone(DM dm, DM *newdm)
4533: {
4534: ALE::Obj<PETSC_MESH_TYPE> m;
4535: void *ctx;
4536: PetscErrorCode ierr;
4541: DMCreate(((PetscObject) dm)->comm, newdm);
4542: DMSetType(*newdm, DMMESH);
4543: DMMeshGetMesh(dm, m);
4544: ALE::Obj<PETSC_MESH_TYPE> newm = new PETSC_MESH_TYPE(m->comm(), m->getDimension(), m->debug());
4545: newm->copy(m);
4546: DMMeshSetMesh(*newdm, newm);
4547: DMGetApplicationContext(dm, &ctx);
4548: DMSetApplicationContext(*newdm, ctx);
4549: return(0);
4550: }