Actual source code: dm.c
petsc-3.13.6 2020-09-29
1: #include <petsc/private/dmimpl.h>
2: #include <petsc/private/dmlabelimpl.h>
3: #include <petsc/private/petscdsimpl.h>
4: #include <petscdmplex.h>
5: #include <petscdmfield.h>
6: #include <petscsf.h>
7: #include <petscds.h>
9: PetscClassId DM_CLASSID;
10: PetscClassId DMLABEL_CLASSID;
11: PetscLogEvent DM_Convert, DM_GlobalToLocal, DM_LocalToGlobal, DM_LocalToLocal, DM_LocatePoints, DM_Coarsen, DM_Refine, DM_CreateInterpolation, DM_CreateRestriction, DM_CreateInjection, DM_CreateMatrix, DM_Load;
13: const char *const DMBoundaryTypes[] = {"NONE","GHOSTED","MIRROR","PERIODIC","TWIST","DMBoundaryType","DM_BOUNDARY_",0};
14: const char *const DMBoundaryConditionTypes[] = {"INVALID","ESSENTIAL","NATURAL","INVALID","INVALID","ESSENTIAL_FIELD","NATURAL_FIELD","INVALID","INVALID","INVALID","NATURAL_RIEMANN","DMBoundaryConditionType","DM_BC_",0};
15: const char *const DMPolytopeTypes[] = {"vertex", "segment", "tensor_segment", "triangle", "quadrilateral", "tensor_quad", "tetrahedron", "hexahedron", "triangular_prism", "tensor_triangular_prism", "tensor_quadrilateral_prism", "FV_ghost_cell", "unknown", "DMPolytopeType", "DM_POLYTOPE_", 0};
16: /*@
17: DMCreate - Creates an empty DM object. The type can then be set with DMSetType().
19: If you never call DMSetType() it will generate an
20: error when you try to use the vector.
22: Collective
24: Input Parameter:
25: . comm - The communicator for the DM object
27: Output Parameter:
28: . dm - The DM object
30: Level: beginner
32: .seealso: DMSetType(), DMDA, DMSLICED, DMCOMPOSITE, DMPLEX, DMMOAB, DMNETWORK
33: @*/
34: PetscErrorCode DMCreate(MPI_Comm comm,DM *dm)
35: {
36: DM v;
37: PetscDS ds;
42: *dm = NULL;
43: DMInitializePackage();
45: PetscHeaderCreate(v, DM_CLASSID, "DM", "Distribution Manager", "DM", comm, DMDestroy, DMView);
47: v->setupcalled = PETSC_FALSE;
48: v->setfromoptionscalled = PETSC_FALSE;
49: v->ltogmap = NULL;
50: v->bs = 1;
51: v->coloringtype = IS_COLORING_GLOBAL;
52: PetscSFCreate(comm, &v->sf);
53: PetscSFCreate(comm, &v->sectionSF);
54: v->labels = NULL;
55: v->adjacency[0] = PETSC_FALSE;
56: v->adjacency[1] = PETSC_TRUE;
57: v->depthLabel = NULL;
58: v->celltypeLabel = NULL;
59: v->localSection = NULL;
60: v->globalSection = NULL;
61: v->defaultConstraintSection = NULL;
62: v->defaultConstraintMat = NULL;
63: v->L = NULL;
64: v->maxCell = NULL;
65: v->bdtype = NULL;
66: v->dimEmbed = PETSC_DEFAULT;
67: v->dim = PETSC_DETERMINE;
68: {
69: PetscInt i;
70: for (i = 0; i < 10; ++i) {
71: v->nullspaceConstructors[i] = NULL;
72: v->nearnullspaceConstructors[i] = NULL;
73: }
74: }
75: PetscDSCreate(PetscObjectComm((PetscObject) v), &ds);
76: DMSetRegionDS(v, NULL, NULL, ds);
77: PetscDSDestroy(&ds);
78: v->dmBC = NULL;
79: v->coarseMesh = NULL;
80: v->outputSequenceNum = -1;
81: v->outputSequenceVal = 0.0;
82: DMSetVecType(v,VECSTANDARD);
83: DMSetMatType(v,MATAIJ);
85: *dm = v;
86: return(0);
87: }
89: /*@
90: DMClone - Creates a DM object with the same topology as the original.
92: Collective
94: Input Parameter:
95: . dm - The original DM object
97: Output Parameter:
98: . newdm - The new DM object
100: Level: beginner
102: Notes:
103: For some DM implementations this is a shallow clone, the result of which may share (referent counted) information with its parent. For example,
104: DMClone() applied to a DMPLEX object will result in a new DMPLEX that shares the topology with the original DMPLEX. It does not
105: share the PetscSection of the original DM.
107: The clone is considered set up iff the original is.
109: .seealso: DMDestroy(), DMCreate(), DMSetType(), DMSetLocalSection(), DMSetGlobalSection()
111: @*/
112: PetscErrorCode DMClone(DM dm, DM *newdm)
113: {
114: PetscSF sf;
115: Vec coords;
116: void *ctx;
117: PetscInt dim, cdim;
123: DMCreate(PetscObjectComm((PetscObject) dm), newdm);
124: DMCopyLabels(dm, *newdm, PETSC_COPY_VALUES, PETSC_TRUE);
125: (*newdm)->leveldown = dm->leveldown;
126: (*newdm)->levelup = dm->levelup;
127: DMGetDimension(dm, &dim);
128: DMSetDimension(*newdm, dim);
129: if (dm->ops->clone) {
130: (*dm->ops->clone)(dm, newdm);
131: }
132: (*newdm)->setupcalled = dm->setupcalled;
133: DMGetPointSF(dm, &sf);
134: DMSetPointSF(*newdm, sf);
135: DMGetApplicationContext(dm, &ctx);
136: DMSetApplicationContext(*newdm, ctx);
137: if (dm->coordinateDM) {
138: DM ncdm;
139: PetscSection cs;
140: PetscInt pEnd = -1, pEndMax = -1;
142: DMGetLocalSection(dm->coordinateDM, &cs);
143: if (cs) {PetscSectionGetChart(cs, NULL, &pEnd);}
144: MPI_Allreduce(&pEnd,&pEndMax,1,MPIU_INT,MPI_MAX,PetscObjectComm((PetscObject)dm));
145: if (pEndMax >= 0) {
146: DMClone(dm->coordinateDM, &ncdm);
147: DMCopyDisc(dm->coordinateDM, ncdm);
148: DMSetLocalSection(ncdm, cs);
149: DMSetCoordinateDM(*newdm, ncdm);
150: DMDestroy(&ncdm);
151: }
152: }
153: DMGetCoordinateDim(dm, &cdim);
154: DMSetCoordinateDim(*newdm, cdim);
155: DMGetCoordinatesLocal(dm, &coords);
156: if (coords) {
157: DMSetCoordinatesLocal(*newdm, coords);
158: } else {
159: DMGetCoordinates(dm, &coords);
160: if (coords) {DMSetCoordinates(*newdm, coords);}
161: }
162: {
163: PetscBool isper;
164: const PetscReal *maxCell, *L;
165: const DMBoundaryType *bd;
166: DMGetPeriodicity(dm, &isper, &maxCell, &L, &bd);
167: DMSetPeriodicity(*newdm, isper, maxCell, L, bd);
168: }
169: {
170: PetscBool useCone, useClosure;
172: DMGetAdjacency(dm, PETSC_DEFAULT, &useCone, &useClosure);
173: DMSetAdjacency(*newdm, PETSC_DEFAULT, useCone, useClosure);
174: }
175: return(0);
176: }
178: /*@C
179: DMSetVecType - Sets the type of vector created with DMCreateLocalVector() and DMCreateGlobalVector()
181: Logically Collective on da
183: Input Parameter:
184: + da - initial distributed array
185: . ctype - the vector type, currently either VECSTANDARD, VECCUDA, or VECVIENNACL
187: Options Database:
188: . -dm_vec_type ctype
190: Level: intermediate
192: .seealso: DMCreate(), DMDestroy(), DM, DMDAInterpolationType, VecType, DMGetVecType(), DMSetMatType(), DMGetMatType()
193: @*/
194: PetscErrorCode DMSetVecType(DM da,VecType ctype)
195: {
200: PetscFree(da->vectype);
201: PetscStrallocpy(ctype,(char**)&da->vectype);
202: return(0);
203: }
205: /*@C
206: DMGetVecType - Gets the type of vector created with DMCreateLocalVector() and DMCreateGlobalVector()
208: Logically Collective on da
210: Input Parameter:
211: . da - initial distributed array
213: Output Parameter:
214: . ctype - the vector type
216: Level: intermediate
218: .seealso: DMCreate(), DMDestroy(), DM, DMDAInterpolationType, VecType, DMSetMatType(), DMGetMatType(), DMSetVecType()
219: @*/
220: PetscErrorCode DMGetVecType(DM da,VecType *ctype)
221: {
224: *ctype = da->vectype;
225: return(0);
226: }
228: /*@
229: VecGetDM - Gets the DM defining the data layout of the vector
231: Not collective
233: Input Parameter:
234: . v - The Vec
236: Output Parameter:
237: . dm - The DM
239: Level: intermediate
241: .seealso: VecSetDM(), DMGetLocalVector(), DMGetGlobalVector(), DMSetVecType()
242: @*/
243: PetscErrorCode VecGetDM(Vec v, DM *dm)
244: {
250: PetscObjectQuery((PetscObject) v, "__PETSc_dm", (PetscObject*) dm);
251: return(0);
252: }
254: /*@
255: VecSetDM - Sets the DM defining the data layout of the vector.
257: Not collective
259: Input Parameters:
260: + v - The Vec
261: - dm - The DM
263: Note: This is NOT the same as DMCreateGlobalVector() since it does not change the view methods or perform other customization, but merely sets the DM member.
265: Level: intermediate
267: .seealso: VecGetDM(), DMGetLocalVector(), DMGetGlobalVector(), DMSetVecType()
268: @*/
269: PetscErrorCode VecSetDM(Vec v, DM dm)
270: {
276: PetscObjectCompose((PetscObject) v, "__PETSc_dm", (PetscObject) dm);
277: return(0);
278: }
280: /*@C
281: DMSetISColoringType - Sets the type of coloring, global or local, that is created by the DM
283: Logically Collective on dm
285: Input Parameters:
286: + dm - the DM context
287: - ctype - the matrix type
289: Options Database:
290: . -dm_is_coloring_type - global or local
292: Level: intermediate
294: .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMCreateMatrix(), DMSetMatrixPreallocateOnly(), MatType, DMGetMatType(),
295: DMGetISColoringType()
296: @*/
297: PetscErrorCode DMSetISColoringType(DM dm,ISColoringType ctype)
298: {
301: dm->coloringtype = ctype;
302: return(0);
303: }
305: /*@C
306: DMGetISColoringType - Gets the type of coloring, global or local, that is created by the DM
308: Logically Collective on dm
310: Input Parameter:
311: . dm - the DM context
313: Output Parameter:
314: . ctype - the matrix type
316: Options Database:
317: . -dm_is_coloring_type - global or local
319: Level: intermediate
321: .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMCreateMatrix(), DMSetMatrixPreallocateOnly(), MatType, DMGetMatType(),
322: DMGetISColoringType()
323: @*/
324: PetscErrorCode DMGetISColoringType(DM dm,ISColoringType *ctype)
325: {
328: *ctype = dm->coloringtype;
329: return(0);
330: }
332: /*@C
333: DMSetMatType - Sets the type of matrix created with DMCreateMatrix()
335: Logically Collective on dm
337: Input Parameters:
338: + dm - the DM context
339: - ctype - the matrix type
341: Options Database:
342: . -dm_mat_type ctype
344: Level: intermediate
346: .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMCreateMatrix(), DMSetMatrixPreallocateOnly(), MatType, DMGetMatType(), DMSetMatType(), DMGetMatType()
347: @*/
348: PetscErrorCode DMSetMatType(DM dm,MatType ctype)
349: {
354: PetscFree(dm->mattype);
355: PetscStrallocpy(ctype,(char**)&dm->mattype);
356: return(0);
357: }
359: /*@C
360: DMGetMatType - Gets the type of matrix created with DMCreateMatrix()
362: Logically Collective on dm
364: Input Parameter:
365: . dm - the DM context
367: Output Parameter:
368: . ctype - the matrix type
370: Options Database:
371: . -dm_mat_type ctype
373: Level: intermediate
375: .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMCreateMatrix(), DMSetMatrixPreallocateOnly(), MatType, DMSetMatType(), DMSetMatType(), DMGetMatType()
376: @*/
377: PetscErrorCode DMGetMatType(DM dm,MatType *ctype)
378: {
381: *ctype = dm->mattype;
382: return(0);
383: }
385: /*@
386: MatGetDM - Gets the DM defining the data layout of the matrix
388: Not collective
390: Input Parameter:
391: . A - The Mat
393: Output Parameter:
394: . dm - The DM
396: Level: intermediate
398: Developer Note: Since the Mat class doesn't know about the DM class the DM object is associated with
399: the Mat through a PetscObjectCompose() operation
401: .seealso: MatSetDM(), DMCreateMatrix(), DMSetMatType()
402: @*/
403: PetscErrorCode MatGetDM(Mat A, DM *dm)
404: {
410: PetscObjectQuery((PetscObject) A, "__PETSc_dm", (PetscObject*) dm);
411: return(0);
412: }
414: /*@
415: MatSetDM - Sets the DM defining the data layout of the matrix
417: Not collective
419: Input Parameters:
420: + A - The Mat
421: - dm - The DM
423: Level: intermediate
425: Developer Note: Since the Mat class doesn't know about the DM class the DM object is associated with
426: the Mat through a PetscObjectCompose() operation
429: .seealso: MatGetDM(), DMCreateMatrix(), DMSetMatType()
430: @*/
431: PetscErrorCode MatSetDM(Mat A, DM dm)
432: {
438: PetscObjectCompose((PetscObject) A, "__PETSc_dm", (PetscObject) dm);
439: return(0);
440: }
442: /*@C
443: DMSetOptionsPrefix - Sets the prefix used for searching for all
444: DM options in the database.
446: Logically Collective on dm
448: Input Parameter:
449: + da - the DM context
450: - prefix - the prefix to prepend to all option names
452: Notes:
453: A hyphen (-) must NOT be given at the beginning of the prefix name.
454: The first character of all runtime options is AUTOMATICALLY the hyphen.
456: Level: advanced
458: .seealso: DMSetFromOptions()
459: @*/
460: PetscErrorCode DMSetOptionsPrefix(DM dm,const char prefix[])
461: {
466: PetscObjectSetOptionsPrefix((PetscObject)dm,prefix);
467: if (dm->sf) {
468: PetscObjectSetOptionsPrefix((PetscObject)dm->sf,prefix);
469: }
470: if (dm->sectionSF) {
471: PetscObjectSetOptionsPrefix((PetscObject)dm->sectionSF,prefix);
472: }
473: return(0);
474: }
476: /*@C
477: DMAppendOptionsPrefix - Appends to the prefix used for searching for all
478: DM options in the database.
480: Logically Collective on dm
482: Input Parameters:
483: + dm - the DM context
484: - prefix - the prefix string to prepend to all DM option requests
486: Notes:
487: A hyphen (-) must NOT be given at the beginning of the prefix name.
488: The first character of all runtime options is AUTOMATICALLY the hyphen.
490: Level: advanced
492: .seealso: DMSetOptionsPrefix(), DMGetOptionsPrefix()
493: @*/
494: PetscErrorCode DMAppendOptionsPrefix(DM dm,const char prefix[])
495: {
500: PetscObjectAppendOptionsPrefix((PetscObject)dm,prefix);
501: return(0);
502: }
504: /*@C
505: DMGetOptionsPrefix - Gets the prefix used for searching for all
506: DM options in the database.
508: Not Collective
510: Input Parameters:
511: . dm - the DM context
513: Output Parameters:
514: . prefix - pointer to the prefix string used is returned
516: Notes:
517: On the fortran side, the user should pass in a string 'prefix' of
518: sufficient length to hold the prefix.
520: Level: advanced
522: .seealso: DMSetOptionsPrefix(), DMAppendOptionsPrefix()
523: @*/
524: PetscErrorCode DMGetOptionsPrefix(DM dm,const char *prefix[])
525: {
530: PetscObjectGetOptionsPrefix((PetscObject)dm,prefix);
531: return(0);
532: }
534: static PetscErrorCode DMCountNonCyclicReferences(DM dm, PetscBool recurseCoarse, PetscBool recurseFine, PetscInt *ncrefct)
535: {
536: PetscInt refct = ((PetscObject) dm)->refct;
540: *ncrefct = 0;
541: if (dm->coarseMesh && dm->coarseMesh->fineMesh == dm) {
542: refct--;
543: if (recurseCoarse) {
544: PetscInt coarseCount;
546: DMCountNonCyclicReferences(dm->coarseMesh, PETSC_TRUE, PETSC_FALSE,&coarseCount);
547: refct += coarseCount;
548: }
549: }
550: if (dm->fineMesh && dm->fineMesh->coarseMesh == dm) {
551: refct--;
552: if (recurseFine) {
553: PetscInt fineCount;
555: DMCountNonCyclicReferences(dm->fineMesh, PETSC_FALSE, PETSC_TRUE,&fineCount);
556: refct += fineCount;
557: }
558: }
559: *ncrefct = refct;
560: return(0);
561: }
563: PetscErrorCode DMDestroyLabelLinkList_Internal(DM dm)
564: {
565: DMLabelLink next = dm->labels;
569: /* destroy the labels */
570: while (next) {
571: DMLabelLink tmp = next->next;
573: if (next->label == dm->depthLabel) dm->depthLabel = NULL;
574: if (next->label == dm->celltypeLabel) dm->celltypeLabel = NULL;
575: DMLabelDestroy(&next->label);
576: PetscFree(next);
577: next = tmp;
578: }
579: dm->labels = NULL;
580: return(0);
581: }
583: /*@
584: DMDestroy - Destroys a vector packer or DM.
586: Collective on dm
588: Input Parameter:
589: . dm - the DM object to destroy
591: Level: developer
593: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix()
595: @*/
596: PetscErrorCode DMDestroy(DM *dm)
597: {
598: PetscInt cnt;
599: DMNamedVecLink nlink,nnext;
603: if (!*dm) return(0);
606: /* count all non-cyclic references in the doubly-linked list of coarse<->fine meshes */
607: DMCountNonCyclicReferences(*dm,PETSC_TRUE,PETSC_TRUE,&cnt);
608: --((PetscObject)(*dm))->refct;
609: if (--cnt > 0) {*dm = 0; return(0);}
610: if (((PetscObject)(*dm))->refct < 0) return(0);
611: ((PetscObject)(*dm))->refct = 0;
613: DMClearGlobalVectors(*dm);
614: DMClearLocalVectors(*dm);
616: nnext=(*dm)->namedglobal;
617: (*dm)->namedglobal = NULL;
618: for (nlink=nnext; nlink; nlink=nnext) { /* Destroy the named vectors */
619: nnext = nlink->next;
620: if (nlink->status != DMVEC_STATUS_IN) SETERRQ1(((PetscObject)*dm)->comm,PETSC_ERR_ARG_WRONGSTATE,"DM still has Vec named '%s' checked out",nlink->name);
621: PetscFree(nlink->name);
622: VecDestroy(&nlink->X);
623: PetscFree(nlink);
624: }
625: nnext=(*dm)->namedlocal;
626: (*dm)->namedlocal = NULL;
627: for (nlink=nnext; nlink; nlink=nnext) { /* Destroy the named local vectors */
628: nnext = nlink->next;
629: if (nlink->status != DMVEC_STATUS_IN) SETERRQ1(((PetscObject)*dm)->comm,PETSC_ERR_ARG_WRONGSTATE,"DM still has Vec named '%s' checked out",nlink->name);
630: PetscFree(nlink->name);
631: VecDestroy(&nlink->X);
632: PetscFree(nlink);
633: }
635: /* Destroy the list of hooks */
636: {
637: DMCoarsenHookLink link,next;
638: for (link=(*dm)->coarsenhook; link; link=next) {
639: next = link->next;
640: PetscFree(link);
641: }
642: (*dm)->coarsenhook = NULL;
643: }
644: {
645: DMRefineHookLink link,next;
646: for (link=(*dm)->refinehook; link; link=next) {
647: next = link->next;
648: PetscFree(link);
649: }
650: (*dm)->refinehook = NULL;
651: }
652: {
653: DMSubDomainHookLink link,next;
654: for (link=(*dm)->subdomainhook; link; link=next) {
655: next = link->next;
656: PetscFree(link);
657: }
658: (*dm)->subdomainhook = NULL;
659: }
660: {
661: DMGlobalToLocalHookLink link,next;
662: for (link=(*dm)->gtolhook; link; link=next) {
663: next = link->next;
664: PetscFree(link);
665: }
666: (*dm)->gtolhook = NULL;
667: }
668: {
669: DMLocalToGlobalHookLink link,next;
670: for (link=(*dm)->ltoghook; link; link=next) {
671: next = link->next;
672: PetscFree(link);
673: }
674: (*dm)->ltoghook = NULL;
675: }
676: /* Destroy the work arrays */
677: {
678: DMWorkLink link,next;
679: if ((*dm)->workout) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Work array still checked out");
680: for (link=(*dm)->workin; link; link=next) {
681: next = link->next;
682: PetscFree(link->mem);
683: PetscFree(link);
684: }
685: (*dm)->workin = NULL;
686: }
687: /* destroy the labels */
688: DMDestroyLabelLinkList_Internal(*dm);
689: /* destroy the fields */
690: DMClearFields(*dm);
691: /* destroy the boundaries */
692: {
693: DMBoundary next = (*dm)->boundary;
694: while (next) {
695: DMBoundary b = next;
697: next = b->next;
698: PetscFree(b);
699: }
700: }
702: PetscObjectDestroy(&(*dm)->dmksp);
703: PetscObjectDestroy(&(*dm)->dmsnes);
704: PetscObjectDestroy(&(*dm)->dmts);
706: if ((*dm)->ctx && (*dm)->ctxdestroy) {
707: (*(*dm)->ctxdestroy)(&(*dm)->ctx);
708: }
709: MatFDColoringDestroy(&(*dm)->fd);
710: ISLocalToGlobalMappingDestroy(&(*dm)->ltogmap);
711: PetscFree((*dm)->vectype);
712: PetscFree((*dm)->mattype);
714: PetscSectionDestroy(&(*dm)->localSection);
715: PetscSectionDestroy(&(*dm)->globalSection);
716: PetscLayoutDestroy(&(*dm)->map);
717: PetscSectionDestroy(&(*dm)->defaultConstraintSection);
718: MatDestroy(&(*dm)->defaultConstraintMat);
719: PetscSFDestroy(&(*dm)->sf);
720: PetscSFDestroy(&(*dm)->sectionSF);
721: if ((*dm)->useNatural) {
722: if ((*dm)->sfNatural) {
723: PetscSFDestroy(&(*dm)->sfNatural);
724: }
725: PetscObjectDereference((PetscObject) (*dm)->sfMigration);
726: }
727: if ((*dm)->coarseMesh && (*dm)->coarseMesh->fineMesh == *dm) {
728: DMSetFineDM((*dm)->coarseMesh,NULL);
729: }
731: DMDestroy(&(*dm)->coarseMesh);
732: if ((*dm)->fineMesh && (*dm)->fineMesh->coarseMesh == *dm) {
733: DMSetCoarseDM((*dm)->fineMesh,NULL);
734: }
735: DMDestroy(&(*dm)->fineMesh);
736: DMFieldDestroy(&(*dm)->coordinateField);
737: DMDestroy(&(*dm)->coordinateDM);
738: VecDestroy(&(*dm)->coordinates);
739: VecDestroy(&(*dm)->coordinatesLocal);
740: PetscFree((*dm)->L);
741: PetscFree((*dm)->maxCell);
742: PetscFree((*dm)->bdtype);
743: if ((*dm)->transformDestroy) {(*(*dm)->transformDestroy)(*dm, (*dm)->transformCtx);}
744: DMDestroy(&(*dm)->transformDM);
745: VecDestroy(&(*dm)->transform);
747: DMClearDS(*dm);
748: DMDestroy(&(*dm)->dmBC);
749: /* if memory was published with SAWs then destroy it */
750: PetscObjectSAWsViewOff((PetscObject)*dm);
752: if ((*dm)->ops->destroy) {
753: (*(*dm)->ops->destroy)(*dm);
754: }
755: DMMonitorCancel(*dm);
756: /* We do not destroy (*dm)->data here so that we can reference count backend objects */
757: PetscHeaderDestroy(dm);
758: return(0);
759: }
761: /*@
762: DMSetUp - sets up the data structures inside a DM object
764: Collective on dm
766: Input Parameter:
767: . dm - the DM object to setup
769: Level: developer
771: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix()
773: @*/
774: PetscErrorCode DMSetUp(DM dm)
775: {
780: if (dm->setupcalled) return(0);
781: if (dm->ops->setup) {
782: (*dm->ops->setup)(dm);
783: }
784: dm->setupcalled = PETSC_TRUE;
785: return(0);
786: }
788: /*@
789: DMSetFromOptions - sets parameters in a DM from the options database
791: Collective on dm
793: Input Parameter:
794: . dm - the DM object to set options for
796: Options Database:
797: + -dm_preallocate_only - Only preallocate the matrix for DMCreateMatrix(), but do not fill it with zeros
798: . -dm_vec_type <type> - type of vector to create inside DM
799: . -dm_mat_type <type> - type of matrix to create inside DM
800: - -dm_is_coloring_type - <global or local>
802: DMPLEX Specific Checks
803: + -dm_plex_check_symmetry - Check that the adjacency information in the mesh is symmetric - DMPlexCheckSymmetry()
804: . -dm_plex_check_skeleton - Check that each cell has the correct number of vertices (only for homogeneous simplex or tensor meshes) - DMPlexCheckSkeleton()
805: . -dm_plex_check_faces - Check that the faces of each cell give a vertex order this is consistent with what we expect from the cell type - DMPlexCheckFaces()
806: . -dm_plex_check_geometry - Check that cells have positive volume - DMPlexCheckGeometry()
807: . -dm_plex_check_pointsf - Check some necessary conditions for PointSF - DMPlexCheckPointSF()
808: . -dm_plex_check_interface_cones - Check points on inter-partition interfaces have conforming order of cone points - DMPlexCheckInterfaceCones()
809: - -dm_plex_check_all - Perform all the checks above
811: Level: intermediate
813: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(),
814: DMPlexCheckSymmetry(), DMPlexCheckSkeleton(), DMPlexCheckFaces(), DMPlexCheckGeometry(), DMPlexCheckPointSF(), DMPlexCheckInterfaceCones()
816: @*/
817: PetscErrorCode DMSetFromOptions(DM dm)
818: {
819: char typeName[256];
820: PetscBool flg;
825: dm->setfromoptionscalled = PETSC_TRUE;
826: if (dm->sf) {PetscSFSetFromOptions(dm->sf);}
827: if (dm->sectionSF) {PetscSFSetFromOptions(dm->sectionSF);}
828: PetscObjectOptionsBegin((PetscObject)dm);
829: PetscOptionsBool("-dm_preallocate_only","only preallocate matrix, but do not set column indices","DMSetMatrixPreallocateOnly",dm->prealloc_only,&dm->prealloc_only,NULL);
830: PetscOptionsFList("-dm_vec_type","Vector type used for created vectors","DMSetVecType",VecList,dm->vectype,typeName,256,&flg);
831: if (flg) {
832: DMSetVecType(dm,typeName);
833: }
834: PetscOptionsFList("-dm_mat_type","Matrix type used for created matrices","DMSetMatType",MatList,dm->mattype ? dm->mattype : typeName,typeName,sizeof(typeName),&flg);
835: if (flg) {
836: DMSetMatType(dm,typeName);
837: }
838: PetscOptionsEnum("-dm_is_coloring_type","Global or local coloring of Jacobian","DMSetISColoringType",ISColoringTypes,(PetscEnum)dm->coloringtype,(PetscEnum*)&dm->coloringtype,NULL);
839: if (dm->ops->setfromoptions) {
840: (*dm->ops->setfromoptions)(PetscOptionsObject,dm);
841: }
842: /* process any options handlers added with PetscObjectAddOptionsHandler() */
843: PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject) dm);
844: PetscOptionsEnd();
845: return(0);
846: }
848: /*@C
849: DMViewFromOptions - View from Options
851: Collective on DM
853: Input Parameters:
854: + dm - the DM object
855: . obj - Optional object
856: - name - command line option
858: Level: intermediate
859: .seealso: DM, DMView, PetscObjectViewFromOptions(), DMCreate()
860: @*/
861: PetscErrorCode DMViewFromOptions(DM dm,PetscObject obj,const char name[])
862: {
867: PetscObjectViewFromOptions((PetscObject)dm,obj,name);
868: return(0);
869: }
871: /*@C
872: DMView - Views a DM
874: Collective on dm
876: Input Parameter:
877: + dm - the DM object to view
878: - v - the viewer
880: Level: beginner
882: .seealso DMDestroy(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix()
884: @*/
885: PetscErrorCode DMView(DM dm,PetscViewer v)
886: {
887: PetscErrorCode ierr;
888: PetscBool isbinary;
889: PetscMPIInt size;
890: PetscViewerFormat format;
894: if (!v) {
895: PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)dm),&v);
896: }
898: /* Ideally, we would like to have this test on.
899: However, it currently breaks socket viz via GLVis.
900: During DMView(parallel_mesh,glvis_viewer), each
901: process opens a sequential ASCII socket to visualize
902: the local mesh, and PetscObjectView(dm,local_socket)
903: is internally called inside VecView_GLVis, incurring
904: in an error here */
906: PetscViewerCheckWritable(v);
908: PetscViewerGetFormat(v,&format);
909: MPI_Comm_size(PetscObjectComm((PetscObject)dm),&size);
910: if (size == 1 && format == PETSC_VIEWER_LOAD_BALANCE) return(0);
911: PetscObjectPrintClassNamePrefixType((PetscObject)dm,v);
912: PetscObjectTypeCompare((PetscObject)v,PETSCVIEWERBINARY,&isbinary);
913: if (isbinary) {
914: PetscInt classid = DM_FILE_CLASSID;
915: char type[256];
917: PetscViewerBinaryWrite(v,&classid,1,PETSC_INT);
918: PetscStrncpy(type,((PetscObject)dm)->type_name,256);
919: PetscViewerBinaryWrite(v,type,256,PETSC_CHAR);
920: }
921: if (dm->ops->view) {
922: (*dm->ops->view)(dm,v);
923: }
924: return(0);
925: }
927: /*@
928: DMCreateGlobalVector - Creates a global vector from a DM object
930: Collective on dm
932: Input Parameter:
933: . dm - the DM object
935: Output Parameter:
936: . vec - the global vector
938: Level: beginner
940: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix()
942: @*/
943: PetscErrorCode DMCreateGlobalVector(DM dm,Vec *vec)
944: {
950: if (!dm->ops->createglobalvector) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMCreateGlobalVector",((PetscObject)dm)->type_name);
951: (*dm->ops->createglobalvector)(dm,vec);
952: #if defined(PETSC_USE_DEBUG)
953: {
954: DM vdm;
956: VecGetDM(*vec,&vdm);
957: if (!vdm) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"DM type '%s' did not attach the DM to the vector\n",((PetscObject)dm)->type_name);
958: }
959: #endif
960: return(0);
961: }
963: /*@
964: DMCreateLocalVector - Creates a local vector from a DM object
966: Not Collective
968: Input Parameter:
969: . dm - the DM object
971: Output Parameter:
972: . vec - the local vector
974: Level: beginner
976: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix()
978: @*/
979: PetscErrorCode DMCreateLocalVector(DM dm,Vec *vec)
980: {
986: if (!dm->ops->createlocalvector) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMCreateLocalVector",((PetscObject)dm)->type_name);
987: (*dm->ops->createlocalvector)(dm,vec);
988: #if defined(PETSC_USE_DEBUG)
989: {
990: DM vdm;
992: VecGetDM(*vec,&vdm);
993: if (!vdm) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"DM type '%s' did not attach the DM to the vector\n",((PetscObject)dm)->type_name);
994: }
995: #endif
996: return(0);
997: }
999: /*@
1000: DMGetLocalToGlobalMapping - Accesses the local-to-global mapping in a DM.
1002: Collective on dm
1004: Input Parameter:
1005: . dm - the DM that provides the mapping
1007: Output Parameter:
1008: . ltog - the mapping
1010: Level: intermediate
1012: Notes:
1013: This mapping can then be used by VecSetLocalToGlobalMapping() or
1014: MatSetLocalToGlobalMapping().
1016: .seealso: DMCreateLocalVector()
1017: @*/
1018: PetscErrorCode DMGetLocalToGlobalMapping(DM dm,ISLocalToGlobalMapping *ltog)
1019: {
1020: PetscInt bs = -1, bsLocal[2], bsMinMax[2];
1026: if (!dm->ltogmap) {
1027: PetscSection section, sectionGlobal;
1029: DMGetLocalSection(dm, §ion);
1030: if (section) {
1031: const PetscInt *cdofs;
1032: PetscInt *ltog;
1033: PetscInt pStart, pEnd, n, p, k, l;
1035: DMGetGlobalSection(dm, §ionGlobal);
1036: PetscSectionGetChart(section, &pStart, &pEnd);
1037: PetscSectionGetStorageSize(section, &n);
1038: PetscMalloc1(n, <og); /* We want the local+overlap size */
1039: for (p = pStart, l = 0; p < pEnd; ++p) {
1040: PetscInt bdof, cdof, dof, off, c, cind = 0;
1042: /* Should probably use constrained dofs */
1043: PetscSectionGetDof(section, p, &dof);
1044: PetscSectionGetConstraintDof(section, p, &cdof);
1045: PetscSectionGetConstraintIndices(section, p, &cdofs);
1046: PetscSectionGetOffset(sectionGlobal, p, &off);
1047: /* If you have dofs, and constraints, and they are unequal, we set the blocksize to 1 */
1048: bdof = cdof && (dof-cdof) ? 1 : dof;
1049: if (dof) {
1050: if (bs < 0) {bs = bdof;}
1051: else if (bs != bdof) {bs = 1;}
1052: }
1053: for (c = 0; c < dof; ++c, ++l) {
1054: if ((cind < cdof) && (c == cdofs[cind])) ltog[l] = off < 0 ? off-c : off+c;
1055: else ltog[l] = (off < 0 ? -(off+1) : off) + c;
1056: }
1057: }
1058: /* Must have same blocksize on all procs (some might have no points) */
1059: bsLocal[0] = bs < 0 ? PETSC_MAX_INT : bs; bsLocal[1] = bs;
1060: PetscGlobalMinMaxInt(PetscObjectComm((PetscObject) dm), bsLocal, bsMinMax);
1061: if (bsMinMax[0] != bsMinMax[1]) {bs = 1;}
1062: else {bs = bsMinMax[0];}
1063: bs = bs < 0 ? 1 : bs;
1064: /* Must reduce indices by blocksize */
1065: if (bs > 1) {
1066: for (l = 0, k = 0; l < n; l += bs, ++k) ltog[k] = ltog[l]/bs;
1067: n /= bs;
1068: }
1069: ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)dm), bs, n, ltog, PETSC_OWN_POINTER, &dm->ltogmap);
1070: PetscLogObjectParent((PetscObject)dm, (PetscObject)dm->ltogmap);
1071: } else {
1072: if (!dm->ops->getlocaltoglobalmapping) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMGetLocalToGlobalMapping",((PetscObject)dm)->type_name);
1073: (*dm->ops->getlocaltoglobalmapping)(dm);
1074: }
1075: }
1076: *ltog = dm->ltogmap;
1077: return(0);
1078: }
1080: /*@
1081: DMGetBlockSize - Gets the inherent block size associated with a DM
1083: Not Collective
1085: Input Parameter:
1086: . dm - the DM with block structure
1088: Output Parameter:
1089: . bs - the block size, 1 implies no exploitable block structure
1091: Level: intermediate
1093: .seealso: ISCreateBlock(), VecSetBlockSize(), MatSetBlockSize(), DMGetLocalToGlobalMapping()
1094: @*/
1095: PetscErrorCode DMGetBlockSize(DM dm,PetscInt *bs)
1096: {
1100: if (dm->bs < 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"DM does not have enough information to provide a block size yet");
1101: *bs = dm->bs;
1102: return(0);
1103: }
1105: /*@
1106: DMCreateInterpolation - Gets interpolation matrix between two DM objects
1108: Collective on dmc
1110: Input Parameter:
1111: + dmc - the DM object
1112: - dmf - the second, finer DM object
1114: Output Parameter:
1115: + mat - the interpolation
1116: - vec - the scaling (optional)
1118: Level: developer
1120: Notes:
1121: For DMDA objects this only works for "uniform refinement", that is the refined mesh was obtained DMRefine() or the coarse mesh was obtained by
1122: DMCoarsen(). The coordinates set into the DMDA are completely ignored in computing the interpolation.
1124: For DMDA objects you can use this interpolation (more precisely the interpolation from the DMGetCoordinateDM()) to interpolate the mesh coordinate vectors
1125: EXCEPT in the periodic case where it does not make sense since the coordinate vectors are not periodic.
1128: .seealso DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateColoring(), DMCreateMatrix(), DMRefine(), DMCoarsen(), DMCreateRestriction(), DMCreateInterpolationScale()
1130: @*/
1131: PetscErrorCode DMCreateInterpolation(DM dmc,DM dmf,Mat *mat,Vec *vec)
1132: {
1139: if (!dmc->ops->createinterpolation) SETERRQ1(PetscObjectComm((PetscObject)dmc),PETSC_ERR_SUP,"DM type %s does not implement DMCreateInterpolation",((PetscObject)dmc)->type_name);
1140: PetscLogEventBegin(DM_CreateInterpolation,dmc,dmf,0,0);
1141: (*dmc->ops->createinterpolation)(dmc,dmf,mat,vec);
1142: PetscLogEventEnd(DM_CreateInterpolation,dmc,dmf,0,0);
1143: return(0);
1144: }
1146: /*@
1147: DMCreateInterpolationScale - Forms L = 1/(R*1) such that diag(L)*R preserves scale and is thus suitable for state (versus residual) restriction.
1149: Input Parameters:
1150: + dac - DM that defines a coarse mesh
1151: . daf - DM that defines a fine mesh
1152: - mat - the restriction (or interpolation operator) from fine to coarse
1154: Output Parameter:
1155: . scale - the scaled vector
1157: Level: developer
1159: .seealso: DMCreateInterpolation()
1161: @*/
1162: PetscErrorCode DMCreateInterpolationScale(DM dac,DM daf,Mat mat,Vec *scale)
1163: {
1165: Vec fine;
1166: PetscScalar one = 1.0;
1169: DMCreateGlobalVector(daf,&fine);
1170: DMCreateGlobalVector(dac,scale);
1171: VecSet(fine,one);
1172: MatRestrict(mat,fine,*scale);
1173: VecDestroy(&fine);
1174: VecReciprocal(*scale);
1175: return(0);
1176: }
1178: /*@
1179: DMCreateRestriction - Gets restriction matrix between two DM objects
1181: Collective on dmc
1183: Input Parameter:
1184: + dmc - the DM object
1185: - dmf - the second, finer DM object
1187: Output Parameter:
1188: . mat - the restriction
1191: Level: developer
1193: Notes:
1194: For DMDA objects this only works for "uniform refinement", that is the refined mesh was obtained DMRefine() or the coarse mesh was obtained by
1195: DMCoarsen(). The coordinates set into the DMDA are completely ignored in computing the interpolation.
1198: .seealso DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateColoring(), DMCreateMatrix(), DMRefine(), DMCoarsen(), DMCreateInterpolation()
1200: @*/
1201: PetscErrorCode DMCreateRestriction(DM dmc,DM dmf,Mat *mat)
1202: {
1209: if (!dmc->ops->createrestriction) SETERRQ1(PetscObjectComm((PetscObject)dmc),PETSC_ERR_SUP,"DM type %s does not implement DMCreateRestriction",((PetscObject)dmc)->type_name);
1210: PetscLogEventBegin(DM_CreateRestriction,dmc,dmf,0,0);
1211: (*dmc->ops->createrestriction)(dmc,dmf,mat);
1212: PetscLogEventEnd(DM_CreateRestriction,dmc,dmf,0,0);
1213: return(0);
1214: }
1216: /*@
1217: DMCreateInjection - Gets injection matrix between two DM objects
1219: Collective on dac
1221: Input Parameter:
1222: + dac - the DM object
1223: - daf - the second, finer DM object
1225: Output Parameter:
1226: . mat - the injection
1228: Level: developer
1230: Notes:
1231: For DMDA objects this only works for "uniform refinement", that is the refined mesh was obtained DMRefine() or the coarse mesh was obtained by
1232: DMCoarsen(). The coordinates set into the DMDA are completely ignored in computing the injection.
1234: .seealso DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateColoring(), DMCreateMatrix(), DMCreateInterpolation()
1236: @*/
1237: PetscErrorCode DMCreateInjection(DM dac,DM daf,Mat *mat)
1238: {
1245: if (!dac->ops->createinjection) SETERRQ1(PetscObjectComm((PetscObject)dac),PETSC_ERR_SUP,"DM type %s does not implement DMCreateInjection",((PetscObject)dac)->type_name);
1246: PetscLogEventBegin(DM_CreateInjection,dac,daf,0,0);
1247: (*dac->ops->createinjection)(dac,daf,mat);
1248: PetscLogEventEnd(DM_CreateInjection,dac,daf,0,0);
1249: return(0);
1250: }
1252: /*@
1253: DMCreateMassMatrix - Gets mass matrix between two DM objects, M_ij = \int \phi_i \psi_j
1255: Collective on dac
1257: Input Parameter:
1258: + dac - the DM object
1259: - daf - the second, finer DM object
1261: Output Parameter:
1262: . mat - the interpolation
1264: Level: developer
1266: .seealso DMCreateMatrix(), DMRefine(), DMCoarsen(), DMCreateRestriction(), DMCreateInterpolation(), DMCreateInjection()
1267: @*/
1268: PetscErrorCode DMCreateMassMatrix(DM dac, DM daf, Mat *mat)
1269: {
1276: if (!dac->ops->createmassmatrix) SETERRQ1(PetscObjectComm((PetscObject)dac),PETSC_ERR_SUP,"DM type %s does not implement DMCreateMassMatrix",((PetscObject)dac)->type_name);
1277: (*dac->ops->createmassmatrix)(dac, daf, mat);
1278: return(0);
1279: }
1281: /*@
1282: DMCreateColoring - Gets coloring for a DM
1284: Collective on dm
1286: Input Parameter:
1287: + dm - the DM object
1288: - ctype - IS_COLORING_LOCAL or IS_COLORING_GLOBAL
1290: Output Parameter:
1291: . coloring - the coloring
1293: Notes:
1294: Coloring of matrices can be computed directly from the sparse matrix nonzero structure via the MatColoring object or from the mesh from which the
1295: matrix comes from. In general using the mesh produces a more optimal coloring (fewer colors).
1297: This produces a coloring with the distance of 2, see MatSetColoringDistance() which can be used for efficiently computing Jacobians with MatFDColoringCreate()
1299: Level: developer
1301: .seealso DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateMatrix(), DMSetMatType(), MatColoring, MatFDColoringCreate()
1303: @*/
1304: PetscErrorCode DMCreateColoring(DM dm,ISColoringType ctype,ISColoring *coloring)
1305: {
1311: if (!dm->ops->getcoloring) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMCreateColoring",((PetscObject)dm)->type_name);
1312: (*dm->ops->getcoloring)(dm,ctype,coloring);
1313: return(0);
1314: }
1316: /*@
1317: DMCreateMatrix - Gets empty Jacobian for a DM
1319: Collective on dm
1321: Input Parameter:
1322: . dm - the DM object
1324: Output Parameter:
1325: . mat - the empty Jacobian
1327: Level: beginner
1329: Notes:
1330: This properly preallocates the number of nonzeros in the sparse matrix so you
1331: do not need to do it yourself.
1333: By default it also sets the nonzero structure and puts in the zero entries. To prevent setting
1334: the nonzero pattern call DMSetMatrixPreallocateOnly()
1336: For structured grid problems, when you call MatView() on this matrix it is displayed using the global natural ordering, NOT in the ordering used
1337: internally by PETSc.
1339: For structured grid problems, in general it is easiest to use MatSetValuesStencil() or MatSetValuesLocal() to put values into the matrix because MatSetValues() requires
1340: the indices for the global numbering for DMDAs which is complicated.
1342: .seealso DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMSetMatType()
1344: @*/
1345: PetscErrorCode DMCreateMatrix(DM dm,Mat *mat)
1346: {
1352: if (!dm->ops->creatematrix) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMCreateMatrix",((PetscObject)dm)->type_name);
1353: MatInitializePackage();
1354: PetscLogEventBegin(DM_CreateMatrix,0,0,0,0);
1355: (*dm->ops->creatematrix)(dm,mat);
1356: #if defined(PETSC_USE_DEBUG)
1357: {
1358: DM mdm;
1360: MatGetDM(*mat,&mdm);
1361: if (!mdm) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"DM type '%s' did not attach the DM to the matrix\n",((PetscObject)dm)->type_name);
1362: }
1363: #endif
1364: /* Handle nullspace and near nullspace */
1365: if (dm->Nf) {
1366: MatNullSpace nullSpace;
1367: PetscInt Nf;
1369: DMGetNumFields(dm, &Nf);
1370: if (Nf == 1) {
1371: if (dm->nullspaceConstructors[0]) {
1372: (*dm->nullspaceConstructors[0])(dm, 0, &nullSpace);
1373: MatSetNullSpace(*mat, nullSpace);
1374: MatNullSpaceDestroy(&nullSpace);
1375: }
1376: if (dm->nearnullspaceConstructors[0]) {
1377: (*dm->nearnullspaceConstructors[0])(dm, 0, &nullSpace);
1378: MatSetNearNullSpace(*mat, nullSpace);
1379: MatNullSpaceDestroy(&nullSpace);
1380: }
1381: }
1382: }
1383: PetscLogEventEnd(DM_CreateMatrix,0,0,0,0);
1384: return(0);
1385: }
1387: /*@
1388: DMSetMatrixPreallocateOnly - When DMCreateMatrix() is called the matrix will be properly
1389: preallocated but the nonzero structure and zero values will not be set.
1391: Logically Collective on dm
1393: Input Parameter:
1394: + dm - the DM
1395: - only - PETSC_TRUE if only want preallocation
1397: Level: developer
1398: .seealso DMCreateMatrix(), DMSetMatrixStructureOnly()
1399: @*/
1400: PetscErrorCode DMSetMatrixPreallocateOnly(DM dm, PetscBool only)
1401: {
1404: dm->prealloc_only = only;
1405: return(0);
1406: }
1408: /*@
1409: DMSetMatrixStructureOnly - When DMCreateMatrix() is called, the matrix structure will be created
1410: but the array for values will not be allocated.
1412: Logically Collective on dm
1414: Input Parameter:
1415: + dm - the DM
1416: - only - PETSC_TRUE if only want matrix stucture
1418: Level: developer
1419: .seealso DMCreateMatrix(), DMSetMatrixPreallocateOnly()
1420: @*/
1421: PetscErrorCode DMSetMatrixStructureOnly(DM dm, PetscBool only)
1422: {
1425: dm->structure_only = only;
1426: return(0);
1427: }
1429: /*@C
1430: DMGetWorkArray - Gets a work array guaranteed to be at least the input size, restore with DMRestoreWorkArray()
1432: Not Collective
1434: Input Parameters:
1435: + dm - the DM object
1436: . count - The minium size
1437: - dtype - MPI data type, often MPIU_REAL, MPIU_SCALAR, MPIU_INT)
1439: Output Parameter:
1440: . array - the work array
1442: Level: developer
1444: .seealso DMDestroy(), DMCreate()
1445: @*/
1446: PetscErrorCode DMGetWorkArray(DM dm,PetscInt count,MPI_Datatype dtype,void *mem)
1447: {
1449: DMWorkLink link;
1450: PetscMPIInt dsize;
1455: if (dm->workin) {
1456: link = dm->workin;
1457: dm->workin = dm->workin->next;
1458: } else {
1459: PetscNewLog(dm,&link);
1460: }
1461: MPI_Type_size(dtype,&dsize);
1462: if (((size_t)dsize*count) > link->bytes) {
1463: PetscFree(link->mem);
1464: PetscMalloc(dsize*count,&link->mem);
1465: link->bytes = dsize*count;
1466: }
1467: link->next = dm->workout;
1468: dm->workout = link;
1469: *(void**)mem = link->mem;
1470: return(0);
1471: }
1473: /*@C
1474: DMRestoreWorkArray - Restores a work array guaranteed to be at least the input size, restore with DMRestoreWorkArray()
1476: Not Collective
1478: Input Parameters:
1479: + dm - the DM object
1480: . count - The minium size
1481: - dtype - MPI data type, often MPIU_REAL, MPIU_SCALAR, MPIU_INT
1483: Output Parameter:
1484: . array - the work array
1486: Level: developer
1488: Developer Notes:
1489: count and dtype are ignored, they are only needed for DMGetWorkArray()
1490: .seealso DMDestroy(), DMCreate()
1491: @*/
1492: PetscErrorCode DMRestoreWorkArray(DM dm,PetscInt count,MPI_Datatype dtype,void *mem)
1493: {
1494: DMWorkLink *p,link;
1499: for (p=&dm->workout; (link=*p); p=&link->next) {
1500: if (link->mem == *(void**)mem) {
1501: *p = link->next;
1502: link->next = dm->workin;
1503: dm->workin = link;
1504: *(void**)mem = NULL;
1505: return(0);
1506: }
1507: }
1508: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Array was not checked out");
1509: }
1511: PetscErrorCode DMSetNullSpaceConstructor(DM dm, PetscInt field, PetscErrorCode (*nullsp)(DM dm, PetscInt field, MatNullSpace *nullSpace))
1512: {
1515: if (field >= 10) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle %d >= 10 fields", field);
1516: dm->nullspaceConstructors[field] = nullsp;
1517: return(0);
1518: }
1520: PetscErrorCode DMGetNullSpaceConstructor(DM dm, PetscInt field, PetscErrorCode (**nullsp)(DM dm, PetscInt field, MatNullSpace *nullSpace))
1521: {
1525: if (field >= 10) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle %d >= 10 fields", field);
1526: *nullsp = dm->nullspaceConstructors[field];
1527: return(0);
1528: }
1530: PetscErrorCode DMSetNearNullSpaceConstructor(DM dm, PetscInt field, PetscErrorCode (*nullsp)(DM dm, PetscInt field, MatNullSpace *nullSpace))
1531: {
1534: if (field >= 10) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle %d >= 10 fields", field);
1535: dm->nearnullspaceConstructors[field] = nullsp;
1536: return(0);
1537: }
1539: PetscErrorCode DMGetNearNullSpaceConstructor(DM dm, PetscInt field, PetscErrorCode (**nullsp)(DM dm, PetscInt field, MatNullSpace *nullSpace))
1540: {
1544: if (field >= 10) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle %d >= 10 fields", field);
1545: *nullsp = dm->nearnullspaceConstructors[field];
1546: return(0);
1547: }
1549: /*@C
1550: DMCreateFieldIS - Creates a set of IS objects with the global indices of dofs for each field
1552: Not collective
1554: Input Parameter:
1555: . dm - the DM object
1557: Output Parameters:
1558: + numFields - The number of fields (or NULL if not requested)
1559: . fieldNames - The name for each field (or NULL if not requested)
1560: - fields - The global indices for each field (or NULL if not requested)
1562: Level: intermediate
1564: Notes:
1565: The user is responsible for freeing all requested arrays. In particular, every entry of names should be freed with
1566: PetscFree(), every entry of fields should be destroyed with ISDestroy(), and both arrays should be freed with
1567: PetscFree().
1569: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix()
1570: @*/
1571: PetscErrorCode DMCreateFieldIS(DM dm, PetscInt *numFields, char ***fieldNames, IS **fields)
1572: {
1573: PetscSection section, sectionGlobal;
1578: if (numFields) {
1580: *numFields = 0;
1581: }
1582: if (fieldNames) {
1584: *fieldNames = NULL;
1585: }
1586: if (fields) {
1588: *fields = NULL;
1589: }
1590: DMGetLocalSection(dm, §ion);
1591: if (section) {
1592: PetscInt *fieldSizes, *fieldNc, **fieldIndices;
1593: PetscInt nF, f, pStart, pEnd, p;
1595: DMGetGlobalSection(dm, §ionGlobal);
1596: PetscSectionGetNumFields(section, &nF);
1597: PetscMalloc3(nF,&fieldSizes,nF,&fieldNc,nF,&fieldIndices);
1598: PetscSectionGetChart(sectionGlobal, &pStart, &pEnd);
1599: for (f = 0; f < nF; ++f) {
1600: fieldSizes[f] = 0;
1601: PetscSectionGetFieldComponents(section, f, &fieldNc[f]);
1602: }
1603: for (p = pStart; p < pEnd; ++p) {
1604: PetscInt gdof;
1606: PetscSectionGetDof(sectionGlobal, p, &gdof);
1607: if (gdof > 0) {
1608: for (f = 0; f < nF; ++f) {
1609: PetscInt fdof, fcdof, fpdof;
1611: PetscSectionGetFieldDof(section, p, f, &fdof);
1612: PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);
1613: fpdof = fdof-fcdof;
1614: if (fpdof && fpdof != fieldNc[f]) {
1615: /* Layout does not admit a pointwise block size */
1616: fieldNc[f] = 1;
1617: }
1618: fieldSizes[f] += fpdof;
1619: }
1620: }
1621: }
1622: for (f = 0; f < nF; ++f) {
1623: PetscMalloc1(fieldSizes[f], &fieldIndices[f]);
1624: fieldSizes[f] = 0;
1625: }
1626: for (p = pStart; p < pEnd; ++p) {
1627: PetscInt gdof, goff;
1629: PetscSectionGetDof(sectionGlobal, p, &gdof);
1630: if (gdof > 0) {
1631: PetscSectionGetOffset(sectionGlobal, p, &goff);
1632: for (f = 0; f < nF; ++f) {
1633: PetscInt fdof, fcdof, fc;
1635: PetscSectionGetFieldDof(section, p, f, &fdof);
1636: PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);
1637: for (fc = 0; fc < fdof-fcdof; ++fc, ++fieldSizes[f]) {
1638: fieldIndices[f][fieldSizes[f]] = goff++;
1639: }
1640: }
1641: }
1642: }
1643: if (numFields) *numFields = nF;
1644: if (fieldNames) {
1645: PetscMalloc1(nF, fieldNames);
1646: for (f = 0; f < nF; ++f) {
1647: const char *fieldName;
1649: PetscSectionGetFieldName(section, f, &fieldName);
1650: PetscStrallocpy(fieldName, (char**) &(*fieldNames)[f]);
1651: }
1652: }
1653: if (fields) {
1654: PetscMalloc1(nF, fields);
1655: for (f = 0; f < nF; ++f) {
1656: PetscInt bs, in[2], out[2];
1658: ISCreateGeneral(PetscObjectComm((PetscObject)dm), fieldSizes[f], fieldIndices[f], PETSC_OWN_POINTER, &(*fields)[f]);
1659: in[0] = -fieldNc[f];
1660: in[1] = fieldNc[f];
1661: MPIU_Allreduce(in, out, 2, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)dm));
1662: bs = (-out[0] == out[1]) ? out[1] : 1;
1663: ISSetBlockSize((*fields)[f], bs);
1664: }
1665: }
1666: PetscFree3(fieldSizes,fieldNc,fieldIndices);
1667: } else if (dm->ops->createfieldis) {
1668: (*dm->ops->createfieldis)(dm, numFields, fieldNames, fields);
1669: }
1670: return(0);
1671: }
1674: /*@C
1675: DMCreateFieldDecomposition - Returns a list of IS objects defining a decomposition of a problem into subproblems
1676: corresponding to different fields: each IS contains the global indices of the dofs of the
1677: corresponding field. The optional list of DMs define the DM for each subproblem.
1678: Generalizes DMCreateFieldIS().
1680: Not collective
1682: Input Parameter:
1683: . dm - the DM object
1685: Output Parameters:
1686: + len - The number of subproblems in the field decomposition (or NULL if not requested)
1687: . namelist - The name for each field (or NULL if not requested)
1688: . islist - The global indices for each field (or NULL if not requested)
1689: - dmlist - The DMs for each field subproblem (or NULL, if not requested; if NULL is returned, no DMs are defined)
1691: Level: intermediate
1693: Notes:
1694: The user is responsible for freeing all requested arrays. In particular, every entry of names should be freed with
1695: PetscFree(), every entry of is should be destroyed with ISDestroy(), every entry of dm should be destroyed with DMDestroy(),
1696: and all of the arrays should be freed with PetscFree().
1698: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateFieldIS()
1699: @*/
1700: PetscErrorCode DMCreateFieldDecomposition(DM dm, PetscInt *len, char ***namelist, IS **islist, DM **dmlist)
1701: {
1706: if (len) {
1708: *len = 0;
1709: }
1710: if (namelist) {
1712: *namelist = 0;
1713: }
1714: if (islist) {
1716: *islist = 0;
1717: }
1718: if (dmlist) {
1720: *dmlist = 0;
1721: }
1722: /*
1723: Is it a good idea to apply the following check across all impls?
1724: Perhaps some impls can have a well-defined decomposition before DMSetUp?
1725: This, however, follows the general principle that accessors are not well-behaved until the object is set up.
1726: */
1727: if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE, "Decomposition defined only after DMSetUp");
1728: if (!dm->ops->createfielddecomposition) {
1729: PetscSection section;
1730: PetscInt numFields, f;
1732: DMGetLocalSection(dm, §ion);
1733: if (section) {PetscSectionGetNumFields(section, &numFields);}
1734: if (section && numFields && dm->ops->createsubdm) {
1735: if (len) *len = numFields;
1736: if (namelist) {PetscMalloc1(numFields,namelist);}
1737: if (islist) {PetscMalloc1(numFields,islist);}
1738: if (dmlist) {PetscMalloc1(numFields,dmlist);}
1739: for (f = 0; f < numFields; ++f) {
1740: const char *fieldName;
1742: DMCreateSubDM(dm, 1, &f, islist ? &(*islist)[f] : NULL, dmlist ? &(*dmlist)[f] : NULL);
1743: if (namelist) {
1744: PetscSectionGetFieldName(section, f, &fieldName);
1745: PetscStrallocpy(fieldName, (char**) &(*namelist)[f]);
1746: }
1747: }
1748: } else {
1749: DMCreateFieldIS(dm, len, namelist, islist);
1750: /* By default there are no DMs associated with subproblems. */
1751: if (dmlist) *dmlist = NULL;
1752: }
1753: } else {
1754: (*dm->ops->createfielddecomposition)(dm,len,namelist,islist,dmlist);
1755: }
1756: return(0);
1757: }
1759: /*@
1760: DMCreateSubDM - Returns an IS and DM encapsulating a subproblem defined by the fields passed in.
1761: The fields are defined by DMCreateFieldIS().
1763: Not collective
1765: Input Parameters:
1766: + dm - The DM object
1767: . numFields - The number of fields in this subproblem
1768: - fields - The field numbers of the selected fields
1770: Output Parameters:
1771: + is - The global indices for the subproblem
1772: - subdm - The DM for the subproblem
1774: Note: You need to call DMPlexSetMigrationSF() on the original DM if you want the Global-To-Natural map to be automatically constructed
1776: Level: intermediate
1778: .seealso DMPlexSetMigrationSF(), DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateFieldIS()
1779: @*/
1780: PetscErrorCode DMCreateSubDM(DM dm, PetscInt numFields, const PetscInt fields[], IS *is, DM *subdm)
1781: {
1789: if (!dm->ops->createsubdm) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMCreateSubDM",((PetscObject)dm)->type_name);
1790: (*dm->ops->createsubdm)(dm, numFields, fields, is, subdm);
1791: return(0);
1792: }
1794: /*@C
1795: DMCreateSuperDM - Returns an arrays of ISes and DM encapsulating a superproblem defined by the DMs passed in.
1797: Not collective
1799: Input Parameter:
1800: + dms - The DM objects
1801: - len - The number of DMs
1803: Output Parameters:
1804: + is - The global indices for the subproblem, or NULL
1805: - superdm - The DM for the superproblem
1807: Note: You need to call DMPlexSetMigrationSF() on the original DM if you want the Global-To-Natural map to be automatically constructed
1809: Level: intermediate
1811: .seealso DMPlexSetMigrationSF(), DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateFieldIS()
1812: @*/
1813: PetscErrorCode DMCreateSuperDM(DM dms[], PetscInt len, IS **is, DM *superdm)
1814: {
1815: PetscInt i;
1823: if (len < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of DMs must be nonnegative: %D", len);
1824: if (len) {
1825: DM dm = dms[0];
1826: if (!dm->ops->createsuperdm) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMCreateSuperDM",((PetscObject)dm)->type_name);
1827: (*dm->ops->createsuperdm)(dms, len, is, superdm);
1828: }
1829: return(0);
1830: }
1833: /*@C
1834: DMCreateDomainDecomposition - Returns lists of IS objects defining a decomposition of a problem into subproblems
1835: corresponding to restrictions to pairs nested subdomains: each IS contains the global
1836: indices of the dofs of the corresponding subdomains. The inner subdomains conceptually
1837: define a nonoverlapping covering, while outer subdomains can overlap.
1838: The optional list of DMs define the DM for each subproblem.
1840: Not collective
1842: Input Parameter:
1843: . dm - the DM object
1845: Output Parameters:
1846: + len - The number of subproblems in the domain decomposition (or NULL if not requested)
1847: . namelist - The name for each subdomain (or NULL if not requested)
1848: . innerislist - The global indices for each inner subdomain (or NULL, if not requested)
1849: . outerislist - The global indices for each outer subdomain (or NULL, if not requested)
1850: - dmlist - The DMs for each subdomain subproblem (or NULL, if not requested; if NULL is returned, no DMs are defined)
1852: Level: intermediate
1854: Notes:
1855: The user is responsible for freeing all requested arrays. In particular, every entry of names should be freed with
1856: PetscFree(), every entry of is should be destroyed with ISDestroy(), every entry of dm should be destroyed with DMDestroy(),
1857: and all of the arrays should be freed with PetscFree().
1859: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateFieldDecomposition()
1860: @*/
1861: PetscErrorCode DMCreateDomainDecomposition(DM dm, PetscInt *len, char ***namelist, IS **innerislist, IS **outerislist, DM **dmlist)
1862: {
1863: PetscErrorCode ierr;
1864: DMSubDomainHookLink link;
1865: PetscInt i,l;
1874: /*
1875: Is it a good idea to apply the following check across all impls?
1876: Perhaps some impls can have a well-defined decomposition before DMSetUp?
1877: This, however, follows the general principle that accessors are not well-behaved until the object is set up.
1878: */
1879: if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE, "Decomposition defined only after DMSetUp");
1880: if (dm->ops->createdomaindecomposition) {
1881: (*dm->ops->createdomaindecomposition)(dm,&l,namelist,innerislist,outerislist,dmlist);
1882: /* copy subdomain hooks and context over to the subdomain DMs */
1883: if (dmlist && *dmlist) {
1884: for (i = 0; i < l; i++) {
1885: for (link=dm->subdomainhook; link; link=link->next) {
1886: if (link->ddhook) {(*link->ddhook)(dm,(*dmlist)[i],link->ctx);}
1887: }
1888: if (dm->ctx) (*dmlist)[i]->ctx = dm->ctx;
1889: }
1890: }
1891: if (len) *len = l;
1892: }
1893: return(0);
1894: }
1897: /*@C
1898: DMCreateDomainDecompositionScatters - Returns scatters to the subdomain vectors from the global vector
1900: Not collective
1902: Input Parameters:
1903: + dm - the DM object
1904: . n - the number of subdomain scatters
1905: - subdms - the local subdomains
1907: Output Parameters:
1908: + n - the number of scatters returned
1909: . iscat - scatter from global vector to nonoverlapping global vector entries on subdomain
1910: . oscat - scatter from global vector to overlapping global vector entries on subdomain
1911: - gscat - scatter from global vector to local vector on subdomain (fills in ghosts)
1913: Notes:
1914: This is an alternative to the iis and ois arguments in DMCreateDomainDecomposition that allow for the solution
1915: of general nonlinear problems with overlapping subdomain methods. While merely having index sets that enable subsets
1916: of the residual equations to be created is fine for linear problems, nonlinear problems require local assembly of
1917: solution and residual data.
1919: Level: developer
1921: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateFieldIS()
1922: @*/
1923: PetscErrorCode DMCreateDomainDecompositionScatters(DM dm,PetscInt n,DM *subdms,VecScatter **iscat,VecScatter **oscat,VecScatter **gscat)
1924: {
1930: if (!dm->ops->createddscatters) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMCreateDomainDecompositionScatters",((PetscObject)dm)->type_name);
1931: (*dm->ops->createddscatters)(dm,n,subdms,iscat,oscat,gscat);
1932: return(0);
1933: }
1935: /*@
1936: DMRefine - Refines a DM object
1938: Collective on dm
1940: Input Parameter:
1941: + dm - the DM object
1942: - comm - the communicator to contain the new DM object (or MPI_COMM_NULL)
1944: Output Parameter:
1945: . dmf - the refined DM, or NULL
1947: Options Dtabase Keys:
1948: . -dm_plex_cell_refiner <strategy> - chooses the refinement strategy, e.g. regular, tohex
1950: Note: If no refinement was done, the return value is NULL
1952: Level: developer
1954: .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation()
1955: @*/
1956: PetscErrorCode DMRefine(DM dm,MPI_Comm comm,DM *dmf)
1957: {
1958: PetscErrorCode ierr;
1959: DMRefineHookLink link;
1963: if (!dm->ops->refine) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMRefine",((PetscObject)dm)->type_name);
1964: PetscLogEventBegin(DM_Refine,dm,0,0,0);
1965: (*dm->ops->refine)(dm,comm,dmf);
1966: if (*dmf) {
1967: (*dmf)->ops->creatematrix = dm->ops->creatematrix;
1969: PetscObjectCopyFortranFunctionPointers((PetscObject)dm,(PetscObject)*dmf);
1971: (*dmf)->ctx = dm->ctx;
1972: (*dmf)->leveldown = dm->leveldown;
1973: (*dmf)->levelup = dm->levelup + 1;
1975: DMSetMatType(*dmf,dm->mattype);
1976: for (link=dm->refinehook; link; link=link->next) {
1977: if (link->refinehook) {
1978: (*link->refinehook)(dm,*dmf,link->ctx);
1979: }
1980: }
1981: }
1982: PetscLogEventEnd(DM_Refine,dm,0,0,0);
1983: return(0);
1984: }
1986: /*@C
1987: DMRefineHookAdd - adds a callback to be run when interpolating a nonlinear problem to a finer grid
1989: Logically Collective
1991: Input Arguments:
1992: + coarse - nonlinear solver context on which to run a hook when restricting to a coarser level
1993: . refinehook - function to run when setting up a coarser level
1994: . interphook - function to run to update data on finer levels (once per SNESSolve())
1995: - ctx - [optional] user-defined context for provide data for the hooks (may be NULL)
1997: Calling sequence of refinehook:
1998: $ refinehook(DM coarse,DM fine,void *ctx);
2000: + coarse - coarse level DM
2001: . fine - fine level DM to interpolate problem to
2002: - ctx - optional user-defined function context
2004: Calling sequence for interphook:
2005: $ interphook(DM coarse,Mat interp,DM fine,void *ctx)
2007: + coarse - coarse level DM
2008: . interp - matrix interpolating a coarse-level solution to the finer grid
2009: . fine - fine level DM to update
2010: - ctx - optional user-defined function context
2012: Level: advanced
2014: Notes:
2015: This function is only needed if auxiliary data needs to be passed to fine grids while grid sequencing
2017: If this function is called multiple times, the hooks will be run in the order they are added.
2019: This function is currently not available from Fortran.
2021: .seealso: DMCoarsenHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2022: @*/
2023: PetscErrorCode DMRefineHookAdd(DM coarse,PetscErrorCode (*refinehook)(DM,DM,void*),PetscErrorCode (*interphook)(DM,Mat,DM,void*),void *ctx)
2024: {
2025: PetscErrorCode ierr;
2026: DMRefineHookLink link,*p;
2030: for (p=&coarse->refinehook; *p; p=&(*p)->next) { /* Scan to the end of the current list of hooks */
2031: if ((*p)->refinehook == refinehook && (*p)->interphook == interphook && (*p)->ctx == ctx) return(0);
2032: }
2033: PetscNew(&link);
2034: link->refinehook = refinehook;
2035: link->interphook = interphook;
2036: link->ctx = ctx;
2037: link->next = NULL;
2038: *p = link;
2039: return(0);
2040: }
2042: /*@C
2043: DMRefineHookRemove - remove a callback from the list of hooks to be run when interpolating a nonlinear problem to a finer grid
2045: Logically Collective
2047: Input Arguments:
2048: + coarse - nonlinear solver context on which to run a hook when restricting to a coarser level
2049: . refinehook - function to run when setting up a coarser level
2050: . interphook - function to run to update data on finer levels (once per SNESSolve())
2051: - ctx - [optional] user-defined context for provide data for the hooks (may be NULL)
2053: Level: advanced
2055: Notes:
2056: This function does nothing if the hook is not in the list.
2058: This function is currently not available from Fortran.
2060: .seealso: DMCoarsenHookRemove(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2061: @*/
2062: PetscErrorCode DMRefineHookRemove(DM coarse,PetscErrorCode (*refinehook)(DM,DM,void*),PetscErrorCode (*interphook)(DM,Mat,DM,void*),void *ctx)
2063: {
2064: PetscErrorCode ierr;
2065: DMRefineHookLink link,*p;
2069: for (p=&coarse->refinehook; *p; p=&(*p)->next) { /* Search the list of current hooks */
2070: if ((*p)->refinehook == refinehook && (*p)->interphook == interphook && (*p)->ctx == ctx) {
2071: link = *p;
2072: *p = link->next;
2073: PetscFree(link);
2074: break;
2075: }
2076: }
2077: return(0);
2078: }
2080: /*@
2081: DMInterpolate - interpolates user-defined problem data to a finer DM by running hooks registered by DMRefineHookAdd()
2083: Collective if any hooks are
2085: Input Arguments:
2086: + coarse - coarser DM to use as a base
2087: . interp - interpolation matrix, apply using MatInterpolate()
2088: - fine - finer DM to update
2090: Level: developer
2092: .seealso: DMRefineHookAdd(), MatInterpolate()
2093: @*/
2094: PetscErrorCode DMInterpolate(DM coarse,Mat interp,DM fine)
2095: {
2096: PetscErrorCode ierr;
2097: DMRefineHookLink link;
2100: for (link=fine->refinehook; link; link=link->next) {
2101: if (link->interphook) {
2102: (*link->interphook)(coarse,interp,fine,link->ctx);
2103: }
2104: }
2105: return(0);
2106: }
2108: /*@
2109: DMGetRefineLevel - Gets the number of refinements that have generated this DM.
2111: Not Collective
2113: Input Parameter:
2114: . dm - the DM object
2116: Output Parameter:
2117: . level - number of refinements
2119: Level: developer
2121: .seealso DMCoarsen(), DMGetCoarsenLevel(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation()
2123: @*/
2124: PetscErrorCode DMGetRefineLevel(DM dm,PetscInt *level)
2125: {
2128: *level = dm->levelup;
2129: return(0);
2130: }
2132: /*@
2133: DMSetRefineLevel - Sets the number of refinements that have generated this DM.
2135: Not Collective
2137: Input Parameter:
2138: + dm - the DM object
2139: - level - number of refinements
2141: Level: advanced
2143: Notes:
2144: This value is used by PCMG to determine how many multigrid levels to use
2146: .seealso DMCoarsen(), DMGetCoarsenLevel(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation()
2148: @*/
2149: PetscErrorCode DMSetRefineLevel(DM dm,PetscInt level)
2150: {
2153: dm->levelup = level;
2154: return(0);
2155: }
2157: PetscErrorCode DMGetBasisTransformDM_Internal(DM dm, DM *tdm)
2158: {
2162: *tdm = dm->transformDM;
2163: return(0);
2164: }
2166: PetscErrorCode DMGetBasisTransformVec_Internal(DM dm, Vec *tv)
2167: {
2171: *tv = dm->transform;
2172: return(0);
2173: }
2175: /*@
2176: DMHasBasisTransform - Whether we employ a basis transformation from functions in global vectors to functions in local vectors
2178: Input Parameter:
2179: . dm - The DM
2181: Output Parameter:
2182: . flg - PETSC_TRUE if a basis transformation should be done
2184: Level: developer
2186: .seealso: DMPlexGlobalToLocalBasis(), DMPlexLocalToGlobalBasis(), DMPlexCreateBasisRotation()
2187: @*/
2188: PetscErrorCode DMHasBasisTransform(DM dm, PetscBool *flg)
2189: {
2190: Vec tv;
2196: DMGetBasisTransformVec_Internal(dm, &tv);
2197: *flg = tv ? PETSC_TRUE : PETSC_FALSE;
2198: return(0);
2199: }
2201: PetscErrorCode DMConstructBasisTransform_Internal(DM dm)
2202: {
2203: PetscSection s, ts;
2204: PetscScalar *ta;
2205: PetscInt cdim, pStart, pEnd, p, Nf, f, Nc, dof;
2209: DMGetCoordinateDim(dm, &cdim);
2210: DMGetLocalSection(dm, &s);
2211: PetscSectionGetChart(s, &pStart, &pEnd);
2212: PetscSectionGetNumFields(s, &Nf);
2213: DMClone(dm, &dm->transformDM);
2214: DMGetLocalSection(dm->transformDM, &ts);
2215: PetscSectionSetNumFields(ts, Nf);
2216: PetscSectionSetChart(ts, pStart, pEnd);
2217: for (f = 0; f < Nf; ++f) {
2218: PetscSectionGetFieldComponents(s, f, &Nc);
2219: /* We could start to label fields by their transformation properties */
2220: if (Nc != cdim) continue;
2221: for (p = pStart; p < pEnd; ++p) {
2222: PetscSectionGetFieldDof(s, p, f, &dof);
2223: if (!dof) continue;
2224: PetscSectionSetFieldDof(ts, p, f, PetscSqr(cdim));
2225: PetscSectionAddDof(ts, p, PetscSqr(cdim));
2226: }
2227: }
2228: PetscSectionSetUp(ts);
2229: DMCreateLocalVector(dm->transformDM, &dm->transform);
2230: VecGetArray(dm->transform, &ta);
2231: for (p = pStart; p < pEnd; ++p) {
2232: for (f = 0; f < Nf; ++f) {
2233: PetscSectionGetFieldDof(ts, p, f, &dof);
2234: if (dof) {
2235: PetscReal x[3] = {0.0, 0.0, 0.0};
2236: PetscScalar *tva;
2237: const PetscScalar *A;
2239: /* TODO Get quadrature point for this dual basis vector for coordinate */
2240: (*dm->transformGetMatrix)(dm, x, PETSC_TRUE, &A, dm->transformCtx);
2241: DMPlexPointLocalFieldRef(dm->transformDM, p, f, ta, (void *) &tva);
2242: PetscArraycpy(tva, A, PetscSqr(cdim));
2243: }
2244: }
2245: }
2246: VecRestoreArray(dm->transform, &ta);
2247: return(0);
2248: }
2250: PetscErrorCode DMCopyTransform(DM dm, DM newdm)
2251: {
2257: newdm->transformCtx = dm->transformCtx;
2258: newdm->transformSetUp = dm->transformSetUp;
2259: newdm->transformDestroy = NULL;
2260: newdm->transformGetMatrix = dm->transformGetMatrix;
2261: if (newdm->transformSetUp) {DMConstructBasisTransform_Internal(newdm);}
2262: return(0);
2263: }
2265: /*@C
2266: DMGlobalToLocalHookAdd - adds a callback to be run when global to local is called
2268: Logically Collective
2270: Input Arguments:
2271: + dm - the DM
2272: . beginhook - function to run at the beginning of DMGlobalToLocalBegin()
2273: . endhook - function to run after DMGlobalToLocalEnd() has completed
2274: - ctx - [optional] user-defined context for provide data for the hooks (may be NULL)
2276: Calling sequence for beginhook:
2277: $ beginhook(DM fine,VecScatter out,VecScatter in,DM coarse,void *ctx)
2279: + dm - global DM
2280: . g - global vector
2281: . mode - mode
2282: . l - local vector
2283: - ctx - optional user-defined function context
2286: Calling sequence for endhook:
2287: $ endhook(DM fine,VecScatter out,VecScatter in,DM coarse,void *ctx)
2289: + global - global DM
2290: - ctx - optional user-defined function context
2292: Level: advanced
2294: .seealso: DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2295: @*/
2296: PetscErrorCode DMGlobalToLocalHookAdd(DM dm,PetscErrorCode (*beginhook)(DM,Vec,InsertMode,Vec,void*),PetscErrorCode (*endhook)(DM,Vec,InsertMode,Vec,void*),void *ctx)
2297: {
2298: PetscErrorCode ierr;
2299: DMGlobalToLocalHookLink link,*p;
2303: for (p=&dm->gtolhook; *p; p=&(*p)->next) {} /* Scan to the end of the current list of hooks */
2304: PetscNew(&link);
2305: link->beginhook = beginhook;
2306: link->endhook = endhook;
2307: link->ctx = ctx;
2308: link->next = NULL;
2309: *p = link;
2310: return(0);
2311: }
2313: static PetscErrorCode DMGlobalToLocalHook_Constraints(DM dm, Vec g, InsertMode mode, Vec l, void *ctx)
2314: {
2315: Mat cMat;
2316: Vec cVec;
2317: PetscSection section, cSec;
2318: PetscInt pStart, pEnd, p, dof;
2323: DMGetDefaultConstraints(dm,&cSec,&cMat);
2324: if (cMat && (mode == INSERT_VALUES || mode == INSERT_ALL_VALUES || mode == INSERT_BC_VALUES)) {
2325: PetscInt nRows;
2327: MatGetSize(cMat,&nRows,NULL);
2328: if (nRows <= 0) return(0);
2329: DMGetLocalSection(dm,§ion);
2330: MatCreateVecs(cMat,NULL,&cVec);
2331: MatMult(cMat,l,cVec);
2332: PetscSectionGetChart(cSec,&pStart,&pEnd);
2333: for (p = pStart; p < pEnd; p++) {
2334: PetscSectionGetDof(cSec,p,&dof);
2335: if (dof) {
2336: PetscScalar *vals;
2337: VecGetValuesSection(cVec,cSec,p,&vals);
2338: VecSetValuesSection(l,section,p,vals,INSERT_ALL_VALUES);
2339: }
2340: }
2341: VecDestroy(&cVec);
2342: }
2343: return(0);
2344: }
2346: /*@
2347: DMGlobalToLocal - update local vectors from global vector
2349: Neighbor-wise Collective on dm
2351: Input Parameters:
2352: + dm - the DM object
2353: . g - the global vector
2354: . mode - INSERT_VALUES or ADD_VALUES
2355: - l - the local vector
2357: Notes:
2358: The communication involved in this update can be overlapped with computation by using
2359: DMGlobalToLocalBegin() and DMGlobalToLocalEnd().
2361: Level: beginner
2363: .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMGlobalToLocalEnd(), DMLocalToGlobalBegin(), DMLocalToGlobal(), DMLocalToGlobalBegin(), DMLocalToGlobalEnd()
2365: @*/
2366: PetscErrorCode DMGlobalToLocal(DM dm,Vec g,InsertMode mode,Vec l)
2367: {
2371: DMGlobalToLocalBegin(dm,g,mode,l);
2372: DMGlobalToLocalEnd(dm,g,mode,l);
2373: return(0);
2374: }
2376: /*@
2377: DMGlobalToLocalBegin - Begins updating local vectors from global vector
2379: Neighbor-wise Collective on dm
2381: Input Parameters:
2382: + dm - the DM object
2383: . g - the global vector
2384: . mode - INSERT_VALUES or ADD_VALUES
2385: - l - the local vector
2387: Level: intermediate
2389: .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMGlobalToLocal(), DMGlobalToLocalEnd(), DMLocalToGlobalBegin(), DMLocalToGlobal(), DMLocalToGlobalBegin(), DMLocalToGlobalEnd()
2391: @*/
2392: PetscErrorCode DMGlobalToLocalBegin(DM dm,Vec g,InsertMode mode,Vec l)
2393: {
2394: PetscSF sf;
2395: PetscErrorCode ierr;
2396: DMGlobalToLocalHookLink link;
2400: for (link=dm->gtolhook; link; link=link->next) {
2401: if (link->beginhook) {
2402: (*link->beginhook)(dm,g,mode,l,link->ctx);
2403: }
2404: }
2405: DMGetSectionSF(dm, &sf);
2406: if (sf) {
2407: const PetscScalar *gArray;
2408: PetscScalar *lArray;
2410: if (mode == ADD_VALUES) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode);
2411: VecGetArrayInPlace(l, &lArray);
2412: VecGetArrayReadInPlace(g, &gArray);
2413: PetscSFBcastBegin(sf, MPIU_SCALAR, gArray, lArray);
2414: VecRestoreArrayInPlace(l, &lArray);
2415: VecRestoreArrayReadInPlace(g, &gArray);
2416: } else {
2417: if (!dm->ops->globaltolocalbegin) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Missing DMGlobalToLocalBegin() for type %s",((PetscObject)dm)->type_name);
2418: (*dm->ops->globaltolocalbegin)(dm,g,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),l);
2419: }
2420: return(0);
2421: }
2423: /*@
2424: DMGlobalToLocalEnd - Ends updating local vectors from global vector
2426: Neighbor-wise Collective on dm
2428: Input Parameters:
2429: + dm - the DM object
2430: . g - the global vector
2431: . mode - INSERT_VALUES or ADD_VALUES
2432: - l - the local vector
2434: Level: intermediate
2436: .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMGlobalToLocal(), DMLocalToGlobalBegin(), DMLocalToGlobal(), DMLocalToGlobalBegin(), DMLocalToGlobalEnd()
2438: @*/
2439: PetscErrorCode DMGlobalToLocalEnd(DM dm,Vec g,InsertMode mode,Vec l)
2440: {
2441: PetscSF sf;
2442: PetscErrorCode ierr;
2443: const PetscScalar *gArray;
2444: PetscScalar *lArray;
2445: PetscBool transform;
2446: DMGlobalToLocalHookLink link;
2450: DMGetSectionSF(dm, &sf);
2451: DMHasBasisTransform(dm, &transform);
2452: if (sf) {
2453: if (mode == ADD_VALUES) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode);
2455: VecGetArrayInPlace(l, &lArray);
2456: VecGetArrayReadInPlace(g, &gArray);
2457: PetscSFBcastEnd(sf, MPIU_SCALAR, gArray, lArray);
2458: VecRestoreArrayInPlace(l, &lArray);
2459: VecRestoreArrayReadInPlace(g, &gArray);
2460: if (transform) {DMPlexGlobalToLocalBasis(dm, l);}
2461: } else {
2462: if (!dm->ops->globaltolocalend) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Missing DMGlobalToLocalEnd() for type %s",((PetscObject)dm)->type_name);
2463: (*dm->ops->globaltolocalend)(dm,g,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),l);
2464: }
2465: DMGlobalToLocalHook_Constraints(dm,g,mode,l,NULL);
2466: for (link=dm->gtolhook; link; link=link->next) {
2467: if (link->endhook) {(*link->endhook)(dm,g,mode,l,link->ctx);}
2468: }
2469: return(0);
2470: }
2472: /*@C
2473: DMLocalToGlobalHookAdd - adds a callback to be run when a local to global is called
2475: Logically Collective
2477: Input Arguments:
2478: + dm - the DM
2479: . beginhook - function to run at the beginning of DMLocalToGlobalBegin()
2480: . endhook - function to run after DMLocalToGlobalEnd() has completed
2481: - ctx - [optional] user-defined context for provide data for the hooks (may be NULL)
2483: Calling sequence for beginhook:
2484: $ beginhook(DM fine,Vec l,InsertMode mode,Vec g,void *ctx)
2486: + dm - global DM
2487: . l - local vector
2488: . mode - mode
2489: . g - global vector
2490: - ctx - optional user-defined function context
2493: Calling sequence for endhook:
2494: $ endhook(DM fine,Vec l,InsertMode mode,Vec g,void *ctx)
2496: + global - global DM
2497: . l - local vector
2498: . mode - mode
2499: . g - global vector
2500: - ctx - optional user-defined function context
2502: Level: advanced
2504: .seealso: DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2505: @*/
2506: PetscErrorCode DMLocalToGlobalHookAdd(DM dm,PetscErrorCode (*beginhook)(DM,Vec,InsertMode,Vec,void*),PetscErrorCode (*endhook)(DM,Vec,InsertMode,Vec,void*),void *ctx)
2507: {
2508: PetscErrorCode ierr;
2509: DMLocalToGlobalHookLink link,*p;
2513: for (p=&dm->ltoghook; *p; p=&(*p)->next) {} /* Scan to the end of the current list of hooks */
2514: PetscNew(&link);
2515: link->beginhook = beginhook;
2516: link->endhook = endhook;
2517: link->ctx = ctx;
2518: link->next = NULL;
2519: *p = link;
2520: return(0);
2521: }
2523: static PetscErrorCode DMLocalToGlobalHook_Constraints(DM dm, Vec l, InsertMode mode, Vec g, void *ctx)
2524: {
2525: Mat cMat;
2526: Vec cVec;
2527: PetscSection section, cSec;
2528: PetscInt pStart, pEnd, p, dof;
2533: DMGetDefaultConstraints(dm,&cSec,&cMat);
2534: if (cMat && (mode == ADD_VALUES || mode == ADD_ALL_VALUES || mode == ADD_BC_VALUES)) {
2535: PetscInt nRows;
2537: MatGetSize(cMat,&nRows,NULL);
2538: if (nRows <= 0) return(0);
2539: DMGetLocalSection(dm,§ion);
2540: MatCreateVecs(cMat,NULL,&cVec);
2541: PetscSectionGetChart(cSec,&pStart,&pEnd);
2542: for (p = pStart; p < pEnd; p++) {
2543: PetscSectionGetDof(cSec,p,&dof);
2544: if (dof) {
2545: PetscInt d;
2546: PetscScalar *vals;
2547: VecGetValuesSection(l,section,p,&vals);
2548: VecSetValuesSection(cVec,cSec,p,vals,mode);
2549: /* for this to be the true transpose, we have to zero the values that
2550: * we just extracted */
2551: for (d = 0; d < dof; d++) {
2552: vals[d] = 0.;
2553: }
2554: }
2555: }
2556: MatMultTransposeAdd(cMat,cVec,l,l);
2557: VecDestroy(&cVec);
2558: }
2559: return(0);
2560: }
2561: /*@
2562: DMLocalToGlobal - updates global vectors from local vectors
2564: Neighbor-wise Collective on dm
2566: Input Parameters:
2567: + dm - the DM object
2568: . l - the local vector
2569: . mode - if INSERT_VALUES then no parallel communication is used, if ADD_VALUES then all ghost points from the same base point accumulate into that base point.
2570: - g - the global vector
2572: Notes:
2573: The communication involved in this update can be overlapped with computation by using
2574: DMLocalToGlobalBegin() and DMLocalToGlobalEnd().
2576: In the ADD_VALUES case you normally would zero the receiving vector before beginning this operation.
2577: INSERT_VALUES is not supported for DMDA; in that case simply compute the values directly into a global vector instead of a local one.
2579: Level: beginner
2581: .seealso DMLocalToGlobalBegin(), DMLocalToGlobalEnd(), DMCoarsen(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMGlobalToLocal(), DMGlobalToLocalEnd(), DMGlobalToLocalBegin()
2583: @*/
2584: PetscErrorCode DMLocalToGlobal(DM dm,Vec l,InsertMode mode,Vec g)
2585: {
2589: DMLocalToGlobalBegin(dm,l,mode,g);
2590: DMLocalToGlobalEnd(dm,l,mode,g);
2591: return(0);
2592: }
2594: /*@
2595: DMLocalToGlobalBegin - begins updating global vectors from local vectors
2597: Neighbor-wise Collective on dm
2599: Input Parameters:
2600: + dm - the DM object
2601: . l - the local vector
2602: . mode - if INSERT_VALUES then no parallel communication is used, if ADD_VALUES then all ghost points from the same base point accumulate into that base point.
2603: - g - the global vector
2605: Notes:
2606: In the ADD_VALUES case you normally would zero the receiving vector before beginning this operation.
2607: INSERT_VALUES is not supported for DMDA, in that case simply compute the values directly into a global vector instead of a local one.
2609: Level: intermediate
2611: .seealso DMLocalToGlobal(), DMLocalToGlobalEnd(), DMCoarsen(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMGlobalToLocal(), DMGlobalToLocalEnd(), DMGlobalToLocalBegin()
2613: @*/
2614: PetscErrorCode DMLocalToGlobalBegin(DM dm,Vec l,InsertMode mode,Vec g)
2615: {
2616: PetscSF sf;
2617: PetscSection s, gs;
2618: DMLocalToGlobalHookLink link;
2619: Vec tmpl;
2620: const PetscScalar *lArray;
2621: PetscScalar *gArray;
2622: PetscBool isInsert, transform, l_inplace = PETSC_FALSE, g_inplace = PETSC_FALSE;
2623: PetscErrorCode ierr;
2627: for (link=dm->ltoghook; link; link=link->next) {
2628: if (link->beginhook) {
2629: (*link->beginhook)(dm,l,mode,g,link->ctx);
2630: }
2631: }
2632: DMLocalToGlobalHook_Constraints(dm,l,mode,g,NULL);
2633: DMGetSectionSF(dm, &sf);
2634: DMGetLocalSection(dm, &s);
2635: switch (mode) {
2636: case INSERT_VALUES:
2637: case INSERT_ALL_VALUES:
2638: case INSERT_BC_VALUES:
2639: isInsert = PETSC_TRUE; break;
2640: case ADD_VALUES:
2641: case ADD_ALL_VALUES:
2642: case ADD_BC_VALUES:
2643: isInsert = PETSC_FALSE; break;
2644: default:
2645: SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode);
2646: }
2647: if ((sf && !isInsert) || (s && isInsert)) {
2648: DMHasBasisTransform(dm, &transform);
2649: if (transform) {
2650: DMGetNamedLocalVector(dm, "__petsc_dm_transform_local_copy", &tmpl);
2651: VecCopy(l, tmpl);
2652: DMPlexLocalToGlobalBasis(dm, tmpl);
2653: VecGetArrayRead(tmpl, &lArray);
2654: } else if (isInsert) {
2655: VecGetArrayRead(l, &lArray);
2656: } else {
2657: VecGetArrayReadInPlace(l, &lArray);
2658: l_inplace = PETSC_TRUE;
2659: }
2660: if (s && isInsert) {
2661: VecGetArray(g, &gArray);
2662: } else {
2663: VecGetArrayInPlace(g, &gArray);
2664: g_inplace = PETSC_TRUE;
2665: }
2666: if (sf && !isInsert) {
2667: PetscSFReduceBegin(sf, MPIU_SCALAR, lArray, gArray, MPIU_SUM);
2668: } else if (s && isInsert) {
2669: PetscInt gStart, pStart, pEnd, p;
2671: DMGetGlobalSection(dm, &gs);
2672: PetscSectionGetChart(s, &pStart, &pEnd);
2673: VecGetOwnershipRange(g, &gStart, NULL);
2674: for (p = pStart; p < pEnd; ++p) {
2675: PetscInt dof, gdof, cdof, gcdof, off, goff, d, e;
2677: PetscSectionGetDof(s, p, &dof);
2678: PetscSectionGetDof(gs, p, &gdof);
2679: PetscSectionGetConstraintDof(s, p, &cdof);
2680: PetscSectionGetConstraintDof(gs, p, &gcdof);
2681: PetscSectionGetOffset(s, p, &off);
2682: PetscSectionGetOffset(gs, p, &goff);
2683: /* Ignore off-process data and points with no global data */
2684: if (!gdof || goff < 0) continue;
2685: if (dof != gdof) SETERRQ5(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Inconsistent sizes, p: %d dof: %d gdof: %d cdof: %d gcdof: %d", p, dof, gdof, cdof, gcdof);
2686: /* If no constraints are enforced in the global vector */
2687: if (!gcdof) {
2688: for (d = 0; d < dof; ++d) gArray[goff-gStart+d] = lArray[off+d];
2689: /* If constraints are enforced in the global vector */
2690: } else if (cdof == gcdof) {
2691: const PetscInt *cdofs;
2692: PetscInt cind = 0;
2694: PetscSectionGetConstraintIndices(s, p, &cdofs);
2695: for (d = 0, e = 0; d < dof; ++d) {
2696: if ((cind < cdof) && (d == cdofs[cind])) {++cind; continue;}
2697: gArray[goff-gStart+e++] = lArray[off+d];
2698: }
2699: } else SETERRQ5(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Inconsistent sizes, p: %d dof: %d gdof: %d cdof: %d gcdof: %d", p, dof, gdof, cdof, gcdof);
2700: }
2701: }
2702: if (g_inplace) {
2703: VecRestoreArrayInPlace(g, &gArray);
2704: } else {
2705: VecRestoreArray(g, &gArray);
2706: }
2707: if (transform) {
2708: VecRestoreArrayRead(tmpl, &lArray);
2709: DMRestoreNamedLocalVector(dm, "__petsc_dm_transform_local_copy", &tmpl);
2710: } else if (l_inplace) {
2711: VecRestoreArrayReadInPlace(l, &lArray);
2712: } else {
2713: VecRestoreArrayRead(l, &lArray);
2714: }
2715: } else {
2716: if (!dm->ops->localtoglobalbegin) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Missing DMLocalToGlobalBegin() for type %s",((PetscObject)dm)->type_name);
2717: (*dm->ops->localtoglobalbegin)(dm,l,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),g);
2718: }
2719: return(0);
2720: }
2722: /*@
2723: DMLocalToGlobalEnd - updates global vectors from local vectors
2725: Neighbor-wise Collective on dm
2727: Input Parameters:
2728: + dm - the DM object
2729: . l - the local vector
2730: . mode - INSERT_VALUES or ADD_VALUES
2731: - g - the global vector
2733: Level: intermediate
2735: .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMGlobalToLocalEnd(), DMGlobalToLocalEnd()
2737: @*/
2738: PetscErrorCode DMLocalToGlobalEnd(DM dm,Vec l,InsertMode mode,Vec g)
2739: {
2740: PetscSF sf;
2741: PetscSection s;
2742: DMLocalToGlobalHookLink link;
2743: PetscBool isInsert, transform;
2744: PetscErrorCode ierr;
2748: DMGetSectionSF(dm, &sf);
2749: DMGetLocalSection(dm, &s);
2750: switch (mode) {
2751: case INSERT_VALUES:
2752: case INSERT_ALL_VALUES:
2753: isInsert = PETSC_TRUE; break;
2754: case ADD_VALUES:
2755: case ADD_ALL_VALUES:
2756: isInsert = PETSC_FALSE; break;
2757: default:
2758: SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode);
2759: }
2760: if (sf && !isInsert) {
2761: const PetscScalar *lArray;
2762: PetscScalar *gArray;
2763: Vec tmpl;
2765: DMHasBasisTransform(dm, &transform);
2766: if (transform) {
2767: DMGetNamedLocalVector(dm, "__petsc_dm_transform_local_copy", &tmpl);
2768: VecGetArrayRead(tmpl, &lArray);
2769: } else {
2770: VecGetArrayReadInPlace(l, &lArray);
2771: }
2772: VecGetArrayInPlace(g, &gArray);
2773: PetscSFReduceEnd(sf, MPIU_SCALAR, lArray, gArray, MPIU_SUM);
2774: if (transform) {
2775: VecRestoreArrayRead(tmpl, &lArray);
2776: DMRestoreNamedLocalVector(dm, "__petsc_dm_transform_local_copy", &tmpl);
2777: } else {
2778: VecRestoreArrayReadInPlace(l, &lArray);
2779: }
2780: VecRestoreArrayInPlace(g, &gArray);
2781: } else if (s && isInsert) {
2782: } else {
2783: if (!dm->ops->localtoglobalend) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Missing DMLocalToGlobalEnd() for type %s",((PetscObject)dm)->type_name);
2784: (*dm->ops->localtoglobalend)(dm,l,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),g);
2785: }
2786: for (link=dm->ltoghook; link; link=link->next) {
2787: if (link->endhook) {(*link->endhook)(dm,g,mode,l,link->ctx);}
2788: }
2789: return(0);
2790: }
2792: /*@
2793: DMLocalToLocalBegin - Maps from a local vector (including ghost points
2794: that contain irrelevant values) to another local vector where the ghost
2795: points in the second are set correctly. Must be followed by DMLocalToLocalEnd().
2797: Neighbor-wise Collective on dm
2799: Input Parameters:
2800: + dm - the DM object
2801: . g - the original local vector
2802: - mode - one of INSERT_VALUES or ADD_VALUES
2804: Output Parameter:
2805: . l - the local vector with correct ghost values
2807: Level: intermediate
2809: Notes:
2810: The local vectors used here need not be the same as those
2811: obtained from DMCreateLocalVector(), BUT they
2812: must have the same parallel data layout; they could, for example, be
2813: obtained with VecDuplicate() from the DM originating vectors.
2815: .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateLocalVector(), DMCreateGlobalVector(), DMCreateInterpolation(), DMLocalToLocalEnd(), DMGlobalToLocalEnd(), DMLocalToGlobalBegin()
2817: @*/
2818: PetscErrorCode DMLocalToLocalBegin(DM dm,Vec g,InsertMode mode,Vec l)
2819: {
2820: PetscErrorCode ierr;
2824: if (!dm->ops->localtolocalbegin) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"This DM does not support local to local maps");
2825: (*dm->ops->localtolocalbegin)(dm,g,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),l);
2826: return(0);
2827: }
2829: /*@
2830: DMLocalToLocalEnd - Maps from a local vector (including ghost points
2831: that contain irrelevant values) to another local vector where the ghost
2832: points in the second are set correctly. Must be preceded by DMLocalToLocalBegin().
2834: Neighbor-wise Collective on dm
2836: Input Parameters:
2837: + da - the DM object
2838: . g - the original local vector
2839: - mode - one of INSERT_VALUES or ADD_VALUES
2841: Output Parameter:
2842: . l - the local vector with correct ghost values
2844: Level: intermediate
2846: Notes:
2847: The local vectors used here need not be the same as those
2848: obtained from DMCreateLocalVector(), BUT they
2849: must have the same parallel data layout; they could, for example, be
2850: obtained with VecDuplicate() from the DM originating vectors.
2852: .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateLocalVector(), DMCreateGlobalVector(), DMCreateInterpolation(), DMLocalToLocalBegin(), DMGlobalToLocalEnd(), DMLocalToGlobalBegin()
2854: @*/
2855: PetscErrorCode DMLocalToLocalEnd(DM dm,Vec g,InsertMode mode,Vec l)
2856: {
2857: PetscErrorCode ierr;
2861: if (!dm->ops->localtolocalend) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"This DM does not support local to local maps");
2862: (*dm->ops->localtolocalend)(dm,g,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),l);
2863: return(0);
2864: }
2867: /*@
2868: DMCoarsen - Coarsens a DM object
2870: Collective on dm
2872: Input Parameter:
2873: + dm - the DM object
2874: - comm - the communicator to contain the new DM object (or MPI_COMM_NULL)
2876: Output Parameter:
2877: . dmc - the coarsened DM
2879: Level: developer
2881: .seealso DMRefine(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation()
2883: @*/
2884: PetscErrorCode DMCoarsen(DM dm, MPI_Comm comm, DM *dmc)
2885: {
2886: PetscErrorCode ierr;
2887: DMCoarsenHookLink link;
2891: if (!dm->ops->coarsen) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMCoarsen",((PetscObject)dm)->type_name);
2892: PetscLogEventBegin(DM_Coarsen,dm,0,0,0);
2893: (*dm->ops->coarsen)(dm, comm, dmc);
2894: if (*dmc) {
2895: DMSetCoarseDM(dm,*dmc);
2896: (*dmc)->ops->creatematrix = dm->ops->creatematrix;
2897: PetscObjectCopyFortranFunctionPointers((PetscObject)dm,(PetscObject)*dmc);
2898: (*dmc)->ctx = dm->ctx;
2899: (*dmc)->levelup = dm->levelup;
2900: (*dmc)->leveldown = dm->leveldown + 1;
2901: DMSetMatType(*dmc,dm->mattype);
2902: for (link=dm->coarsenhook; link; link=link->next) {
2903: if (link->coarsenhook) {(*link->coarsenhook)(dm,*dmc,link->ctx);}
2904: }
2905: }
2906: PetscLogEventEnd(DM_Coarsen,dm,0,0,0);
2907: if (!(*dmc)) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "NULL coarse mesh produced");
2908: return(0);
2909: }
2911: /*@C
2912: DMCoarsenHookAdd - adds a callback to be run when restricting a nonlinear problem to the coarse grid
2914: Logically Collective
2916: Input Arguments:
2917: + fine - nonlinear solver context on which to run a hook when restricting to a coarser level
2918: . coarsenhook - function to run when setting up a coarser level
2919: . restricthook - function to run to update data on coarser levels (once per SNESSolve())
2920: - ctx - [optional] user-defined context for provide data for the hooks (may be NULL)
2922: Calling sequence of coarsenhook:
2923: $ coarsenhook(DM fine,DM coarse,void *ctx);
2925: + fine - fine level DM
2926: . coarse - coarse level DM to restrict problem to
2927: - ctx - optional user-defined function context
2929: Calling sequence for restricthook:
2930: $ restricthook(DM fine,Mat mrestrict,Vec rscale,Mat inject,DM coarse,void *ctx)
2932: + fine - fine level DM
2933: . mrestrict - matrix restricting a fine-level solution to the coarse grid
2934: . rscale - scaling vector for restriction
2935: . inject - matrix restricting by injection
2936: . coarse - coarse level DM to update
2937: - ctx - optional user-defined function context
2939: Level: advanced
2941: Notes:
2942: This function is only needed if auxiliary data needs to be set up on coarse grids.
2944: If this function is called multiple times, the hooks will be run in the order they are added.
2946: In order to compose with nonlinear preconditioning without duplicating storage, the hook should be implemented to
2947: extract the finest level information from its context (instead of from the SNES).
2949: This function is currently not available from Fortran.
2951: .seealso: DMCoarsenHookRemove(), DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2952: @*/
2953: PetscErrorCode DMCoarsenHookAdd(DM fine,PetscErrorCode (*coarsenhook)(DM,DM,void*),PetscErrorCode (*restricthook)(DM,Mat,Vec,Mat,DM,void*),void *ctx)
2954: {
2955: PetscErrorCode ierr;
2956: DMCoarsenHookLink link,*p;
2960: for (p=&fine->coarsenhook; *p; p=&(*p)->next) { /* Scan to the end of the current list of hooks */
2961: if ((*p)->coarsenhook == coarsenhook && (*p)->restricthook == restricthook && (*p)->ctx == ctx) return(0);
2962: }
2963: PetscNew(&link);
2964: link->coarsenhook = coarsenhook;
2965: link->restricthook = restricthook;
2966: link->ctx = ctx;
2967: link->next = NULL;
2968: *p = link;
2969: return(0);
2970: }
2972: /*@C
2973: DMCoarsenHookRemove - remove a callback from the list of hooks to be run when restricting a nonlinear problem to the coarse grid
2975: Logically Collective
2977: Input Arguments:
2978: + fine - nonlinear solver context on which to run a hook when restricting to a coarser level
2979: . coarsenhook - function to run when setting up a coarser level
2980: . restricthook - function to run to update data on coarser levels (once per SNESSolve())
2981: - ctx - [optional] user-defined context for provide data for the hooks (may be NULL)
2983: Level: advanced
2985: Notes:
2986: This function does nothing if the hook is not in the list.
2988: This function is currently not available from Fortran.
2990: .seealso: DMCoarsenHookAdd(), DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2991: @*/
2992: PetscErrorCode DMCoarsenHookRemove(DM fine,PetscErrorCode (*coarsenhook)(DM,DM,void*),PetscErrorCode (*restricthook)(DM,Mat,Vec,Mat,DM,void*),void *ctx)
2993: {
2994: PetscErrorCode ierr;
2995: DMCoarsenHookLink link,*p;
2999: for (p=&fine->coarsenhook; *p; p=&(*p)->next) { /* Search the list of current hooks */
3000: if ((*p)->coarsenhook == coarsenhook && (*p)->restricthook == restricthook && (*p)->ctx == ctx) {
3001: link = *p;
3002: *p = link->next;
3003: PetscFree(link);
3004: break;
3005: }
3006: }
3007: return(0);
3008: }
3011: /*@
3012: DMRestrict - restricts user-defined problem data to a coarser DM by running hooks registered by DMCoarsenHookAdd()
3014: Collective if any hooks are
3016: Input Arguments:
3017: + fine - finer DM to use as a base
3018: . restrct - restriction matrix, apply using MatRestrict()
3019: . rscale - scaling vector for restriction
3020: . inject - injection matrix, also use MatRestrict()
3021: - coarse - coarser DM to update
3023: Level: developer
3025: .seealso: DMCoarsenHookAdd(), MatRestrict()
3026: @*/
3027: PetscErrorCode DMRestrict(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse)
3028: {
3029: PetscErrorCode ierr;
3030: DMCoarsenHookLink link;
3033: for (link=fine->coarsenhook; link; link=link->next) {
3034: if (link->restricthook) {
3035: (*link->restricthook)(fine,restrct,rscale,inject,coarse,link->ctx);
3036: }
3037: }
3038: return(0);
3039: }
3041: /*@C
3042: DMSubDomainHookAdd - adds a callback to be run when restricting a problem to the coarse grid
3044: Logically Collective on global
3046: Input Arguments:
3047: + global - global DM
3048: . ddhook - function to run to pass data to the decomposition DM upon its creation
3049: . restricthook - function to run to update data on block solve (at the beginning of the block solve)
3050: - ctx - [optional] user-defined context for provide data for the hooks (may be NULL)
3053: Calling sequence for ddhook:
3054: $ ddhook(DM global,DM block,void *ctx)
3056: + global - global DM
3057: . block - block DM
3058: - ctx - optional user-defined function context
3060: Calling sequence for restricthook:
3061: $ restricthook(DM global,VecScatter out,VecScatter in,DM block,void *ctx)
3063: + global - global DM
3064: . out - scatter to the outer (with ghost and overlap points) block vector
3065: . in - scatter to block vector values only owned locally
3066: . block - block DM
3067: - ctx - optional user-defined function context
3069: Level: advanced
3071: Notes:
3072: This function is only needed if auxiliary data needs to be set up on subdomain DMs.
3074: If this function is called multiple times, the hooks will be run in the order they are added.
3076: In order to compose with nonlinear preconditioning without duplicating storage, the hook should be implemented to
3077: extract the global information from its context (instead of from the SNES).
3079: This function is currently not available from Fortran.
3081: .seealso: DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
3082: @*/
3083: PetscErrorCode DMSubDomainHookAdd(DM global,PetscErrorCode (*ddhook)(DM,DM,void*),PetscErrorCode (*restricthook)(DM,VecScatter,VecScatter,DM,void*),void *ctx)
3084: {
3085: PetscErrorCode ierr;
3086: DMSubDomainHookLink link,*p;
3090: for (p=&global->subdomainhook; *p; p=&(*p)->next) { /* Scan to the end of the current list of hooks */
3091: if ((*p)->ddhook == ddhook && (*p)->restricthook == restricthook && (*p)->ctx == ctx) return(0);
3092: }
3093: PetscNew(&link);
3094: link->restricthook = restricthook;
3095: link->ddhook = ddhook;
3096: link->ctx = ctx;
3097: link->next = NULL;
3098: *p = link;
3099: return(0);
3100: }
3102: /*@C
3103: DMSubDomainHookRemove - remove a callback from the list to be run when restricting a problem to the coarse grid
3105: Logically Collective
3107: Input Arguments:
3108: + global - global DM
3109: . ddhook - function to run to pass data to the decomposition DM upon its creation
3110: . restricthook - function to run to update data on block solve (at the beginning of the block solve)
3111: - ctx - [optional] user-defined context for provide data for the hooks (may be NULL)
3113: Level: advanced
3115: Notes:
3117: This function is currently not available from Fortran.
3119: .seealso: DMSubDomainHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
3120: @*/
3121: PetscErrorCode DMSubDomainHookRemove(DM global,PetscErrorCode (*ddhook)(DM,DM,void*),PetscErrorCode (*restricthook)(DM,VecScatter,VecScatter,DM,void*),void *ctx)
3122: {
3123: PetscErrorCode ierr;
3124: DMSubDomainHookLink link,*p;
3128: for (p=&global->subdomainhook; *p; p=&(*p)->next) { /* Search the list of current hooks */
3129: if ((*p)->ddhook == ddhook && (*p)->restricthook == restricthook && (*p)->ctx == ctx) {
3130: link = *p;
3131: *p = link->next;
3132: PetscFree(link);
3133: break;
3134: }
3135: }
3136: return(0);
3137: }
3139: /*@
3140: DMSubDomainRestrict - restricts user-defined problem data to a block DM by running hooks registered by DMSubDomainHookAdd()
3142: Collective if any hooks are
3144: Input Arguments:
3145: + fine - finer DM to use as a base
3146: . oscatter - scatter from domain global vector filling subdomain global vector with overlap
3147: . gscatter - scatter from domain global vector filling subdomain local vector with ghosts
3148: - coarse - coarer DM to update
3150: Level: developer
3152: .seealso: DMCoarsenHookAdd(), MatRestrict()
3153: @*/
3154: PetscErrorCode DMSubDomainRestrict(DM global,VecScatter oscatter,VecScatter gscatter,DM subdm)
3155: {
3156: PetscErrorCode ierr;
3157: DMSubDomainHookLink link;
3160: for (link=global->subdomainhook; link; link=link->next) {
3161: if (link->restricthook) {
3162: (*link->restricthook)(global,oscatter,gscatter,subdm,link->ctx);
3163: }
3164: }
3165: return(0);
3166: }
3168: /*@
3169: DMGetCoarsenLevel - Get's the number of coarsenings that have generated this DM.
3171: Not Collective
3173: Input Parameter:
3174: . dm - the DM object
3176: Output Parameter:
3177: . level - number of coarsenings
3179: Level: developer
3181: .seealso DMCoarsen(), DMGetRefineLevel(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation()
3183: @*/
3184: PetscErrorCode DMGetCoarsenLevel(DM dm,PetscInt *level)
3185: {
3189: *level = dm->leveldown;
3190: return(0);
3191: }
3193: /*@
3194: DMSetCoarsenLevel - Sets the number of coarsenings that have generated this DM.
3196: Not Collective
3198: Input Parameters:
3199: + dm - the DM object
3200: - level - number of coarsenings
3202: Level: developer
3204: .seealso DMCoarsen(), DMGetCoarsenLevel(), DMGetRefineLevel(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation()
3205: @*/
3206: PetscErrorCode DMSetCoarsenLevel(DM dm,PetscInt level)
3207: {
3210: dm->leveldown = level;
3211: return(0);
3212: }
3216: /*@C
3217: DMRefineHierarchy - Refines a DM object, all levels at once
3219: Collective on dm
3221: Input Parameter:
3222: + dm - the DM object
3223: - nlevels - the number of levels of refinement
3225: Output Parameter:
3226: . dmf - the refined DM hierarchy
3228: Level: developer
3230: .seealso DMCoarsenHierarchy(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation()
3232: @*/
3233: PetscErrorCode DMRefineHierarchy(DM dm,PetscInt nlevels,DM dmf[])
3234: {
3239: if (nlevels < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"nlevels cannot be negative");
3240: if (nlevels == 0) return(0);
3242: if (dm->ops->refinehierarchy) {
3243: (*dm->ops->refinehierarchy)(dm,nlevels,dmf);
3244: } else if (dm->ops->refine) {
3245: PetscInt i;
3247: DMRefine(dm,PetscObjectComm((PetscObject)dm),&dmf[0]);
3248: for (i=1; i<nlevels; i++) {
3249: DMRefine(dmf[i-1],PetscObjectComm((PetscObject)dm),&dmf[i]);
3250: }
3251: } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"No RefineHierarchy for this DM yet");
3252: return(0);
3253: }
3255: /*@C
3256: DMCoarsenHierarchy - Coarsens a DM object, all levels at once
3258: Collective on dm
3260: Input Parameter:
3261: + dm - the DM object
3262: - nlevels - the number of levels of coarsening
3264: Output Parameter:
3265: . dmc - the coarsened DM hierarchy
3267: Level: developer
3269: .seealso DMRefineHierarchy(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation()
3271: @*/
3272: PetscErrorCode DMCoarsenHierarchy(DM dm, PetscInt nlevels, DM dmc[])
3273: {
3278: if (nlevels < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"nlevels cannot be negative");
3279: if (nlevels == 0) return(0);
3281: if (dm->ops->coarsenhierarchy) {
3282: (*dm->ops->coarsenhierarchy)(dm, nlevels, dmc);
3283: } else if (dm->ops->coarsen) {
3284: PetscInt i;
3286: DMCoarsen(dm,PetscObjectComm((PetscObject)dm),&dmc[0]);
3287: for (i=1; i<nlevels; i++) {
3288: DMCoarsen(dmc[i-1],PetscObjectComm((PetscObject)dm),&dmc[i]);
3289: }
3290: } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"No CoarsenHierarchy for this DM yet");
3291: return(0);
3292: }
3294: /*@C
3295: DMSetApplicationContextDestroy - Sets a user function that will be called to destroy the Section 1.5 Writing Application Codes with PETSc context when the DM is destroyed
3297: Not Collective
3299: Input Parameters:
3300: + dm - the DM object
3301: - destroy - the destroy function
3303: Level: intermediate
3305: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext()
3307: @*/
3308: PetscErrorCode DMSetApplicationContextDestroy(DM dm,PetscErrorCode (*destroy)(void**))
3309: {
3312: dm->ctxdestroy = destroy;
3313: return(0);
3314: }
3316: /*@
3317: DMSetApplicationContext - Set a user context into a DM object
3319: Not Collective
3321: Input Parameters:
3322: + dm - the DM object
3323: - ctx - the user context
3325: Level: intermediate
3327: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext()
3329: @*/
3330: PetscErrorCode DMSetApplicationContext(DM dm,void *ctx)
3331: {
3334: dm->ctx = ctx;
3335: return(0);
3336: }
3338: /*@
3339: DMGetApplicationContext - Gets a user context from a DM object
3341: Not Collective
3343: Input Parameter:
3344: . dm - the DM object
3346: Output Parameter:
3347: . ctx - the user context
3349: Level: intermediate
3351: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext()
3353: @*/
3354: PetscErrorCode DMGetApplicationContext(DM dm,void *ctx)
3355: {
3358: *(void**)ctx = dm->ctx;
3359: return(0);
3360: }
3362: /*@C
3363: DMSetVariableBounds - sets a function to compute the lower and upper bound vectors for SNESVI.
3365: Logically Collective on dm
3367: Input Parameter:
3368: + dm - the DM object
3369: - f - the function that computes variable bounds used by SNESVI (use NULL to cancel a previous function that was set)
3371: Level: intermediate
3373: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext(),
3374: DMSetJacobian()
3376: @*/
3377: PetscErrorCode DMSetVariableBounds(DM dm,PetscErrorCode (*f)(DM,Vec,Vec))
3378: {
3381: dm->ops->computevariablebounds = f;
3382: return(0);
3383: }
3385: /*@
3386: DMHasVariableBounds - does the DM object have a variable bounds function?
3388: Not Collective
3390: Input Parameter:
3391: . dm - the DM object to destroy
3393: Output Parameter:
3394: . flg - PETSC_TRUE if the variable bounds function exists
3396: Level: developer
3398: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext()
3400: @*/
3401: PetscErrorCode DMHasVariableBounds(DM dm,PetscBool *flg)
3402: {
3406: *flg = (dm->ops->computevariablebounds) ? PETSC_TRUE : PETSC_FALSE;
3407: return(0);
3408: }
3410: /*@C
3411: DMComputeVariableBounds - compute variable bounds used by SNESVI.
3413: Logically Collective on dm
3415: Input Parameters:
3416: . dm - the DM object
3418: Output parameters:
3419: + xl - lower bound
3420: - xu - upper bound
3422: Level: advanced
3424: Notes:
3425: This is generally not called by users. It calls the function provided by the user with DMSetVariableBounds()
3427: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext()
3429: @*/
3430: PetscErrorCode DMComputeVariableBounds(DM dm, Vec xl, Vec xu)
3431: {
3438: if (!dm->ops->computevariablebounds) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMComputeVariableBounds",((PetscObject)dm)->type_name);
3439: (*dm->ops->computevariablebounds)(dm, xl,xu);
3440: return(0);
3441: }
3443: /*@
3444: DMHasColoring - does the DM object have a method of providing a coloring?
3446: Not Collective
3448: Input Parameter:
3449: . dm - the DM object
3451: Output Parameter:
3452: . flg - PETSC_TRUE if the DM has facilities for DMCreateColoring().
3454: Level: developer
3456: .seealso DMCreateColoring()
3458: @*/
3459: PetscErrorCode DMHasColoring(DM dm,PetscBool *flg)
3460: {
3464: *flg = (dm->ops->getcoloring) ? PETSC_TRUE : PETSC_FALSE;
3465: return(0);
3466: }
3468: /*@
3469: DMHasCreateRestriction - does the DM object have a method of providing a restriction?
3471: Not Collective
3473: Input Parameter:
3474: . dm - the DM object
3476: Output Parameter:
3477: . flg - PETSC_TRUE if the DM has facilities for DMCreateRestriction().
3479: Level: developer
3481: .seealso DMCreateRestriction()
3483: @*/
3484: PetscErrorCode DMHasCreateRestriction(DM dm,PetscBool *flg)
3485: {
3489: *flg = (dm->ops->createrestriction) ? PETSC_TRUE : PETSC_FALSE;
3490: return(0);
3491: }
3494: /*@
3495: DMHasCreateInjection - does the DM object have a method of providing an injection?
3497: Not Collective
3499: Input Parameter:
3500: . dm - the DM object
3502: Output Parameter:
3503: . flg - PETSC_TRUE if the DM has facilities for DMCreateInjection().
3505: Level: developer
3507: .seealso DMCreateInjection()
3509: @*/
3510: PetscErrorCode DMHasCreateInjection(DM dm,PetscBool *flg)
3511: {
3517: if (dm->ops->hascreateinjection) {
3518: (*dm->ops->hascreateinjection)(dm,flg);
3519: } else {
3520: *flg = (dm->ops->createinjection) ? PETSC_TRUE : PETSC_FALSE;
3521: }
3522: return(0);
3523: }
3525: PetscFunctionList DMList = NULL;
3526: PetscBool DMRegisterAllCalled = PETSC_FALSE;
3528: /*@C
3529: DMSetType - Builds a DM, for a particular DM implementation.
3531: Collective on dm
3533: Input Parameters:
3534: + dm - The DM object
3535: - method - The name of the DM type
3537: Options Database Key:
3538: . -dm_type <type> - Sets the DM type; use -help for a list of available types
3540: Notes:
3541: See "petsc/include/petscdm.h" for available DM types (for instance, DM1D, DM2D, or DM3D).
3543: Level: intermediate
3545: .seealso: DMGetType(), DMCreate()
3546: @*/
3547: PetscErrorCode DMSetType(DM dm, DMType method)
3548: {
3549: PetscErrorCode (*r)(DM);
3550: PetscBool match;
3555: PetscObjectTypeCompare((PetscObject) dm, method, &match);
3556: if (match) return(0);
3558: DMRegisterAll();
3559: PetscFunctionListFind(DMList,method,&r);
3560: if (!r) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown DM type: %s", method);
3562: if (dm->ops->destroy) {
3563: (*dm->ops->destroy)(dm);
3564: }
3565: PetscMemzero(dm->ops,sizeof(*dm->ops));
3566: PetscObjectChangeTypeName((PetscObject)dm,method);
3567: (*r)(dm);
3568: return(0);
3569: }
3571: /*@C
3572: DMGetType - Gets the DM type name (as a string) from the DM.
3574: Not Collective
3576: Input Parameter:
3577: . dm - The DM
3579: Output Parameter:
3580: . type - The DM type name
3582: Level: intermediate
3584: .seealso: DMSetType(), DMCreate()
3585: @*/
3586: PetscErrorCode DMGetType(DM dm, DMType *type)
3587: {
3593: DMRegisterAll();
3594: *type = ((PetscObject)dm)->type_name;
3595: return(0);
3596: }
3598: /*@C
3599: DMConvert - Converts a DM to another DM, either of the same or different type.
3601: Collective on dm
3603: Input Parameters:
3604: + dm - the DM
3605: - newtype - new DM type (use "same" for the same type)
3607: Output Parameter:
3608: . M - pointer to new DM
3610: Notes:
3611: Cannot be used to convert a sequential DM to parallel or parallel to sequential,
3612: the MPI communicator of the generated DM is always the same as the communicator
3613: of the input DM.
3615: Level: intermediate
3617: .seealso: DMCreate()
3618: @*/
3619: PetscErrorCode DMConvert(DM dm, DMType newtype, DM *M)
3620: {
3621: DM B;
3622: char convname[256];
3623: PetscBool sametype/*, issame */;
3630: PetscObjectTypeCompare((PetscObject) dm, newtype, &sametype);
3631: /* PetscStrcmp(newtype, "same", &issame); */
3632: if (sametype) {
3633: *M = dm;
3634: PetscObjectReference((PetscObject) dm);
3635: return(0);
3636: } else {
3637: PetscErrorCode (*conv)(DM, DMType, DM*) = NULL;
3639: /*
3640: Order of precedence:
3641: 1) See if a specialized converter is known to the current DM.
3642: 2) See if a specialized converter is known to the desired DM class.
3643: 3) See if a good general converter is registered for the desired class
3644: 4) See if a good general converter is known for the current matrix.
3645: 5) Use a really basic converter.
3646: */
3648: /* 1) See if a specialized converter is known to the current DM and the desired class */
3649: PetscStrncpy(convname,"DMConvert_",sizeof(convname));
3650: PetscStrlcat(convname,((PetscObject) dm)->type_name,sizeof(convname));
3651: PetscStrlcat(convname,"_",sizeof(convname));
3652: PetscStrlcat(convname,newtype,sizeof(convname));
3653: PetscStrlcat(convname,"_C",sizeof(convname));
3654: PetscObjectQueryFunction((PetscObject)dm,convname,&conv);
3655: if (conv) goto foundconv;
3657: /* 2) See if a specialized converter is known to the desired DM class. */
3658: DMCreate(PetscObjectComm((PetscObject)dm), &B);
3659: DMSetType(B, newtype);
3660: PetscStrncpy(convname,"DMConvert_",sizeof(convname));
3661: PetscStrlcat(convname,((PetscObject) dm)->type_name,sizeof(convname));
3662: PetscStrlcat(convname,"_",sizeof(convname));
3663: PetscStrlcat(convname,newtype,sizeof(convname));
3664: PetscStrlcat(convname,"_C",sizeof(convname));
3665: PetscObjectQueryFunction((PetscObject)B,convname,&conv);
3666: if (conv) {
3667: DMDestroy(&B);
3668: goto foundconv;
3669: }
3671: #if 0
3672: /* 3) See if a good general converter is registered for the desired class */
3673: conv = B->ops->convertfrom;
3674: DMDestroy(&B);
3675: if (conv) goto foundconv;
3677: /* 4) See if a good general converter is known for the current matrix */
3678: if (dm->ops->convert) {
3679: conv = dm->ops->convert;
3680: }
3681: if (conv) goto foundconv;
3682: #endif
3684: /* 5) Use a really basic converter. */
3685: SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "No conversion possible between DM types %s and %s", ((PetscObject) dm)->type_name, newtype);
3687: foundconv:
3688: PetscLogEventBegin(DM_Convert,dm,0,0,0);
3689: (*conv)(dm,newtype,M);
3690: /* Things that are independent of DM type: We should consult DMClone() here */
3691: {
3692: PetscBool isper;
3693: const PetscReal *maxCell, *L;
3694: const DMBoundaryType *bd;
3695: DMGetPeriodicity(dm, &isper, &maxCell, &L, &bd);
3696: DMSetPeriodicity(*M, isper, maxCell, L, bd);
3697: }
3698: PetscLogEventEnd(DM_Convert,dm,0,0,0);
3699: }
3700: PetscObjectStateIncrease((PetscObject) *M);
3701: return(0);
3702: }
3704: /*--------------------------------------------------------------------------------------------------------------------*/
3706: /*@C
3707: DMRegister - Adds a new DM component implementation
3709: Not Collective
3711: Input Parameters:
3712: + name - The name of a new user-defined creation routine
3713: - create_func - The creation routine itself
3715: Notes:
3716: DMRegister() may be called multiple times to add several user-defined DMs
3719: Sample usage:
3720: .vb
3721: DMRegister("my_da", MyDMCreate);
3722: .ve
3724: Then, your DM type can be chosen with the procedural interface via
3725: .vb
3726: DMCreate(MPI_Comm, DM *);
3727: DMSetType(DM,"my_da");
3728: .ve
3729: or at runtime via the option
3730: .vb
3731: -da_type my_da
3732: .ve
3734: Level: advanced
3736: .seealso: DMRegisterAll(), DMRegisterDestroy()
3738: @*/
3739: PetscErrorCode DMRegister(const char sname[],PetscErrorCode (*function)(DM))
3740: {
3744: DMInitializePackage();
3745: PetscFunctionListAdd(&DMList,sname,function);
3746: return(0);
3747: }
3749: /*@C
3750: DMLoad - Loads a DM that has been stored in binary with DMView().
3752: Collective on viewer
3754: Input Parameters:
3755: + newdm - the newly loaded DM, this needs to have been created with DMCreate() or
3756: some related function before a call to DMLoad().
3757: - viewer - binary file viewer, obtained from PetscViewerBinaryOpen() or
3758: HDF5 file viewer, obtained from PetscViewerHDF5Open()
3760: Level: intermediate
3762: Notes:
3763: The type is determined by the data in the file, any type set into the DM before this call is ignored.
3765: Notes for advanced users:
3766: Most users should not need to know the details of the binary storage
3767: format, since DMLoad() and DMView() completely hide these details.
3768: But for anyone who's interested, the standard binary matrix storage
3769: format is
3770: .vb
3771: has not yet been determined
3772: .ve
3774: .seealso: PetscViewerBinaryOpen(), DMView(), MatLoad(), VecLoad()
3775: @*/
3776: PetscErrorCode DMLoad(DM newdm, PetscViewer viewer)
3777: {
3778: PetscBool isbinary, ishdf5;
3784: PetscViewerCheckReadable(viewer);
3785: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
3786: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);
3787: PetscLogEventBegin(DM_Load,viewer,0,0,0);
3788: if (isbinary) {
3789: PetscInt classid;
3790: char type[256];
3792: PetscViewerBinaryRead(viewer,&classid,1,NULL,PETSC_INT);
3793: if (classid != DM_FILE_CLASSID) SETERRQ1(PetscObjectComm((PetscObject)newdm),PETSC_ERR_ARG_WRONG,"Not DM next in file, classid found %d",(int)classid);
3794: PetscViewerBinaryRead(viewer,type,256,NULL,PETSC_CHAR);
3795: DMSetType(newdm, type);
3796: if (newdm->ops->load) {(*newdm->ops->load)(newdm,viewer);}
3797: } else if (ishdf5) {
3798: if (newdm->ops->load) {(*newdm->ops->load)(newdm,viewer);}
3799: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid viewer; open viewer with PetscViewerBinaryOpen() or PetscViewerHDF5Open()");
3800: PetscLogEventEnd(DM_Load,viewer,0,0,0);
3801: return(0);
3802: }
3804: /*@
3805: DMGetLocalBoundingBox - Returns the bounding box for the piece of the DM on this process.
3807: Not collective
3809: Input Parameter:
3810: . dm - the DM
3812: Output Parameters:
3813: + lmin - local minimum coordinates (length coord dim, optional)
3814: - lmax - local maximim coordinates (length coord dim, optional)
3816: Level: beginner
3818: Note: If the DM is a DMDA and has no coordinates, the index bounds are returned instead.
3821: .seealso: DMGetCoordinates(), DMGetCoordinatesLocal(), DMGetBoundingBox()
3822: @*/
3823: PetscErrorCode DMGetLocalBoundingBox(DM dm, PetscReal lmin[], PetscReal lmax[])
3824: {
3825: Vec coords = NULL;
3826: PetscReal min[3] = {PETSC_MAX_REAL, PETSC_MAX_REAL, PETSC_MAX_REAL};
3827: PetscReal max[3] = {PETSC_MIN_REAL, PETSC_MIN_REAL, PETSC_MIN_REAL};
3828: const PetscScalar *local_coords;
3829: PetscInt N, Ni;
3830: PetscInt cdim, i, j;
3831: PetscErrorCode ierr;
3835: DMGetCoordinateDim(dm, &cdim);
3836: DMGetCoordinates(dm, &coords);
3837: if (coords) {
3838: VecGetArrayRead(coords, &local_coords);
3839: VecGetLocalSize(coords, &N);
3840: Ni = N/cdim;
3841: for (i = 0; i < Ni; ++i) {
3842: for (j = 0; j < 3; ++j) {
3843: min[j] = j < cdim ? PetscMin(min[j], PetscRealPart(local_coords[i*cdim+j])) : 0;
3844: max[j] = j < cdim ? PetscMax(max[j], PetscRealPart(local_coords[i*cdim+j])) : 0;
3845: }
3846: }
3847: VecRestoreArrayRead(coords, &local_coords);
3848: } else {
3849: PetscBool isda;
3851: PetscObjectTypeCompare((PetscObject) dm, DMDA, &isda);
3852: if (isda) {DMGetLocalBoundingIndices_DMDA(dm, min, max);}
3853: }
3854: if (lmin) {PetscArraycpy(lmin, min, cdim);}
3855: if (lmax) {PetscArraycpy(lmax, max, cdim);}
3856: return(0);
3857: }
3859: /*@
3860: DMGetBoundingBox - Returns the global bounding box for the DM.
3862: Collective
3864: Input Parameter:
3865: . dm - the DM
3867: Output Parameters:
3868: + gmin - global minimum coordinates (length coord dim, optional)
3869: - gmax - global maximim coordinates (length coord dim, optional)
3871: Level: beginner
3873: .seealso: DMGetLocalBoundingBox(), DMGetCoordinates(), DMGetCoordinatesLocal()
3874: @*/
3875: PetscErrorCode DMGetBoundingBox(DM dm, PetscReal gmin[], PetscReal gmax[])
3876: {
3877: PetscReal lmin[3], lmax[3];
3878: PetscInt cdim;
3879: PetscMPIInt count;
3884: DMGetCoordinateDim(dm, &cdim);
3885: PetscMPIIntCast(cdim, &count);
3886: DMGetLocalBoundingBox(dm, lmin, lmax);
3887: if (gmin) {MPIU_Allreduce(lmin, gmin, count, MPIU_REAL, MPIU_MIN, PetscObjectComm((PetscObject) dm));}
3888: if (gmax) {MPIU_Allreduce(lmax, gmax, count, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject) dm));}
3889: return(0);
3890: }
3892: /******************************** FEM Support **********************************/
3894: PetscErrorCode DMPrintCellVector(PetscInt c, const char name[], PetscInt len, const PetscScalar x[])
3895: {
3896: PetscInt f;
3900: PetscPrintf(PETSC_COMM_SELF, "Cell %D Element %s\n", c, name);
3901: for (f = 0; f < len; ++f) {
3902: PetscPrintf(PETSC_COMM_SELF, " | %g |\n", (double)PetscRealPart(x[f]));
3903: }
3904: return(0);
3905: }
3907: PetscErrorCode DMPrintCellMatrix(PetscInt c, const char name[], PetscInt rows, PetscInt cols, const PetscScalar A[])
3908: {
3909: PetscInt f, g;
3913: PetscPrintf(PETSC_COMM_SELF, "Cell %D Element %s\n", c, name);
3914: for (f = 0; f < rows; ++f) {
3915: PetscPrintf(PETSC_COMM_SELF, " |");
3916: for (g = 0; g < cols; ++g) {
3917: PetscPrintf(PETSC_COMM_SELF, " % 9.5g", PetscRealPart(A[f*cols+g]));
3918: }
3919: PetscPrintf(PETSC_COMM_SELF, " |\n");
3920: }
3921: return(0);
3922: }
3924: PetscErrorCode DMPrintLocalVec(DM dm, const char name[], PetscReal tol, Vec X)
3925: {
3926: PetscInt localSize, bs;
3927: PetscMPIInt size;
3928: Vec x, xglob;
3929: const PetscScalar *xarray;
3930: PetscErrorCode ierr;
3933: MPI_Comm_size(PetscObjectComm((PetscObject) dm),&size);
3934: VecDuplicate(X, &x);
3935: VecCopy(X, x);
3936: VecChop(x, tol);
3937: PetscPrintf(PetscObjectComm((PetscObject) dm),"%s:\n",name);
3938: if (size > 1) {
3939: VecGetLocalSize(x,&localSize);
3940: VecGetArrayRead(x,&xarray);
3941: VecGetBlockSize(x,&bs);
3942: VecCreateMPIWithArray(PetscObjectComm((PetscObject) dm),bs,localSize,PETSC_DETERMINE,xarray,&xglob);
3943: } else {
3944: xglob = x;
3945: }
3946: VecView(xglob,PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject) dm)));
3947: if (size > 1) {
3948: VecDestroy(&xglob);
3949: VecRestoreArrayRead(x,&xarray);
3950: }
3951: VecDestroy(&x);
3952: return(0);
3953: }
3955: /*@
3956: DMGetSection - Get the PetscSection encoding the local data layout for the DM. This is equivalent to DMGetLocalSection(). Deprecated in v3.12
3958: Input Parameter:
3959: . dm - The DM
3961: Output Parameter:
3962: . section - The PetscSection
3964: Options Database Keys:
3965: . -dm_petscsection_view - View the Section created by the DM
3967: Level: advanced
3969: Notes:
3970: Use DMGetLocalSection() in new code.
3972: This gets a borrowed reference, so the user should not destroy this PetscSection.
3974: .seealso: DMGetLocalSection(), DMSetLocalSection(), DMGetGlobalSection()
3975: @*/
3976: PetscErrorCode DMGetSection(DM dm, PetscSection *section)
3977: {
3981: DMGetLocalSection(dm,section);
3982: return(0);
3983: }
3985: /*@
3986: DMGetLocalSection - Get the PetscSection encoding the local data layout for the DM.
3988: Input Parameter:
3989: . dm - The DM
3991: Output Parameter:
3992: . section - The PetscSection
3994: Options Database Keys:
3995: . -dm_petscsection_view - View the Section created by the DM
3997: Level: intermediate
3999: Note: This gets a borrowed reference, so the user should not destroy this PetscSection.
4001: .seealso: DMSetLocalSection(), DMGetGlobalSection()
4002: @*/
4003: PetscErrorCode DMGetLocalSection(DM dm, PetscSection *section)
4004: {
4010: if (!dm->localSection && dm->ops->createlocalsection) {
4011: PetscInt d;
4013: if (dm->setfromoptionscalled) for (d = 0; d < dm->Nds; ++d) {PetscDSSetFromOptions(dm->probs[d].ds);}
4014: (*dm->ops->createlocalsection)(dm);
4015: if (dm->localSection) {PetscObjectViewFromOptions((PetscObject) dm->localSection, NULL, "-dm_petscsection_view");}
4016: }
4017: *section = dm->localSection;
4018: return(0);
4019: }
4021: /*@
4022: DMSetSection - Set the PetscSection encoding the local data layout for the DM. This is equivalent to DMSetLocalSection(). Deprecated in v3.12
4024: Input Parameters:
4025: + dm - The DM
4026: - section - The PetscSection
4028: Level: advanced
4030: Notes:
4031: Use DMSetLocalSection() in new code.
4033: Any existing Section will be destroyed
4035: .seealso: DMSetLocalSection(), DMGetLocalSection(), DMSetGlobalSection()
4036: @*/
4037: PetscErrorCode DMSetSection(DM dm, PetscSection section)
4038: {
4042: DMSetLocalSection(dm,section);
4043: return(0);
4044: }
4046: /*@
4047: DMSetLocalSection - Set the PetscSection encoding the local data layout for the DM.
4049: Input Parameters:
4050: + dm - The DM
4051: - section - The PetscSection
4053: Level: intermediate
4055: Note: Any existing Section will be destroyed
4057: .seealso: DMGetLocalSection(), DMSetGlobalSection()
4058: @*/
4059: PetscErrorCode DMSetLocalSection(DM dm, PetscSection section)
4060: {
4061: PetscInt numFields = 0;
4062: PetscInt f;
4068: PetscObjectReference((PetscObject)section);
4069: PetscSectionDestroy(&dm->localSection);
4070: dm->localSection = section;
4071: if (section) {PetscSectionGetNumFields(dm->localSection, &numFields);}
4072: if (numFields) {
4073: DMSetNumFields(dm, numFields);
4074: for (f = 0; f < numFields; ++f) {
4075: PetscObject disc;
4076: const char *name;
4078: PetscSectionGetFieldName(dm->localSection, f, &name);
4079: DMGetField(dm, f, NULL, &disc);
4080: PetscObjectSetName(disc, name);
4081: }
4082: }
4083: /* The global section will be rebuilt in the next call to DMGetGlobalSection(). */
4084: PetscSectionDestroy(&dm->globalSection);
4085: return(0);
4086: }
4088: /*@
4089: DMGetDefaultConstraints - Get the PetscSection and Mat that specify the local constraint interpolation. See DMSetDefaultConstraints() for a description of the purpose of constraint interpolation.
4091: not collective
4093: Input Parameter:
4094: . dm - The DM
4096: Output Parameter:
4097: + section - The PetscSection describing the range of the constraint matrix: relates rows of the constraint matrix to dofs of the default section. Returns NULL if there are no local constraints.
4098: - mat - The Mat that interpolates local constraints: its width should be the layout size of the default section. Returns NULL if there are no local constraints.
4100: Level: advanced
4102: Note: This gets borrowed references, so the user should not destroy the PetscSection or the Mat.
4104: .seealso: DMSetDefaultConstraints()
4105: @*/
4106: PetscErrorCode DMGetDefaultConstraints(DM dm, PetscSection *section, Mat *mat)
4107: {
4112: if (!dm->defaultConstraintSection && !dm->defaultConstraintMat && dm->ops->createdefaultconstraints) {(*dm->ops->createdefaultconstraints)(dm);}
4113: if (section) {*section = dm->defaultConstraintSection;}
4114: if (mat) {*mat = dm->defaultConstraintMat;}
4115: return(0);
4116: }
4118: /*@
4119: DMSetDefaultConstraints - Set the PetscSection and Mat that specify the local constraint interpolation.
4121: If a constraint matrix is specified, then it is applied during DMGlobalToLocalEnd() when mode is INSERT_VALUES, INSERT_BC_VALUES, or INSERT_ALL_VALUES. Without a constraint matrix, the local vector l returned by DMGlobalToLocalEnd() contains values that have been scattered from a global vector without modification; with a constraint matrix A, l is modified by computing c = A * l, l[s[i]] = c[i], where the scatter s is defined by the PetscSection returned by DMGetDefaultConstraintMatrix().
4123: If a constraint matrix is specified, then its adjoint is applied during DMLocalToGlobalBegin() when mode is ADD_VALUES, ADD_BC_VALUES, or ADD_ALL_VALUES. Without a constraint matrix, the local vector l is accumulated into a global vector without modification; with a constraint matrix A, l is first modified by computing c[i] = l[s[i]], l[s[i]] = 0, l = l + A'*c, which is the adjoint of the operation described above.
4125: collective on dm
4127: Input Parameters:
4128: + dm - The DM
4129: + section - The PetscSection describing the range of the constraint matrix: relates rows of the constraint matrix to dofs of the default section. Must have a local communicator (PETSC_COMM_SELF or derivative).
4130: - mat - The Mat that interpolates local constraints: its width should be the layout size of the default section: NULL indicates no constraints. Must have a local communicator (PETSC_COMM_SELF or derivative).
4132: Level: advanced
4134: Note: This increments the references of the PetscSection and the Mat, so they user can destroy them
4136: .seealso: DMGetDefaultConstraints()
4137: @*/
4138: PetscErrorCode DMSetDefaultConstraints(DM dm, PetscSection section, Mat mat)
4139: {
4140: PetscMPIInt result;
4145: if (section) {
4147: MPI_Comm_compare(PETSC_COMM_SELF,PetscObjectComm((PetscObject)section),&result);
4148: if (result != MPI_CONGRUENT && result != MPI_IDENT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMECOMM,"constraint section must have local communicator");
4149: }
4150: if (mat) {
4152: MPI_Comm_compare(PETSC_COMM_SELF,PetscObjectComm((PetscObject)mat),&result);
4153: if (result != MPI_CONGRUENT && result != MPI_IDENT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMECOMM,"constraint matrix must have local communicator");
4154: }
4155: PetscObjectReference((PetscObject)section);
4156: PetscSectionDestroy(&dm->defaultConstraintSection);
4157: dm->defaultConstraintSection = section;
4158: PetscObjectReference((PetscObject)mat);
4159: MatDestroy(&dm->defaultConstraintMat);
4160: dm->defaultConstraintMat = mat;
4161: return(0);
4162: }
4164: #if defined(PETSC_USE_DEBUG)
4165: /*
4166: DMDefaultSectionCheckConsistency - Check the consistentcy of the global and local sections.
4168: Input Parameters:
4169: + dm - The DM
4170: . localSection - PetscSection describing the local data layout
4171: - globalSection - PetscSection describing the global data layout
4173: Level: intermediate
4175: .seealso: DMGetSectionSF(), DMSetSectionSF()
4176: */
4177: static PetscErrorCode DMDefaultSectionCheckConsistency_Internal(DM dm, PetscSection localSection, PetscSection globalSection)
4178: {
4179: MPI_Comm comm;
4180: PetscLayout layout;
4181: const PetscInt *ranges;
4182: PetscInt pStart, pEnd, p, nroots;
4183: PetscMPIInt size, rank;
4184: PetscBool valid = PETSC_TRUE, gvalid;
4185: PetscErrorCode ierr;
4188: PetscObjectGetComm((PetscObject)dm,&comm);
4190: MPI_Comm_size(comm, &size);
4191: MPI_Comm_rank(comm, &rank);
4192: PetscSectionGetChart(globalSection, &pStart, &pEnd);
4193: PetscSectionGetConstrainedStorageSize(globalSection, &nroots);
4194: PetscLayoutCreate(comm, &layout);
4195: PetscLayoutSetBlockSize(layout, 1);
4196: PetscLayoutSetLocalSize(layout, nroots);
4197: PetscLayoutSetUp(layout);
4198: PetscLayoutGetRanges(layout, &ranges);
4199: for (p = pStart; p < pEnd; ++p) {
4200: PetscInt dof, cdof, off, gdof, gcdof, goff, gsize, d;
4202: PetscSectionGetDof(localSection, p, &dof);
4203: PetscSectionGetOffset(localSection, p, &off);
4204: PetscSectionGetConstraintDof(localSection, p, &cdof);
4205: PetscSectionGetDof(globalSection, p, &gdof);
4206: PetscSectionGetConstraintDof(globalSection, p, &gcdof);
4207: PetscSectionGetOffset(globalSection, p, &goff);
4208: if (!gdof) continue; /* Censored point */
4209: if ((gdof < 0 ? -(gdof+1) : gdof) != dof) {PetscSynchronizedPrintf(comm, "[%d]Global dof %d for point %d not equal to local dof %d\n", rank, gdof, p, dof); valid = PETSC_FALSE;}
4210: if (gcdof && (gcdof != cdof)) {PetscSynchronizedPrintf(comm, "[%d]Global constraints %d for point %d not equal to local constraints %d\n", rank, gcdof, p, cdof); valid = PETSC_FALSE;}
4211: if (gdof < 0) {
4212: gsize = gdof < 0 ? -(gdof+1)-gcdof : gdof-gcdof;
4213: for (d = 0; d < gsize; ++d) {
4214: PetscInt offset = -(goff+1) + d, r;
4216: PetscFindInt(offset,size+1,ranges,&r);
4217: if (r < 0) r = -(r+2);
4218: if ((r < 0) || (r >= size)) {PetscSynchronizedPrintf(comm, "[%d]Point %d mapped to invalid process %d (%d, %d)\n", rank, p, r, gdof, goff); valid = PETSC_FALSE;break;}
4219: }
4220: }
4221: }
4222: PetscLayoutDestroy(&layout);
4223: PetscSynchronizedFlush(comm, NULL);
4224: MPIU_Allreduce(&valid, &gvalid, 1, MPIU_BOOL, MPI_LAND, comm);
4225: if (!gvalid) {
4226: DMView(dm, NULL);
4227: SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Inconsistent local and global sections");
4228: }
4229: return(0);
4230: }
4231: #endif
4233: /*@
4234: DMGetGlobalSection - Get the PetscSection encoding the global data layout for the DM.
4236: Collective on dm
4238: Input Parameter:
4239: . dm - The DM
4241: Output Parameter:
4242: . section - The PetscSection
4244: Level: intermediate
4246: Note: This gets a borrowed reference, so the user should not destroy this PetscSection.
4248: .seealso: DMSetLocalSection(), DMGetLocalSection()
4249: @*/
4250: PetscErrorCode DMGetGlobalSection(DM dm, PetscSection *section)
4251: {
4257: if (!dm->globalSection) {
4258: PetscSection s;
4260: DMGetLocalSection(dm, &s);
4261: if (!s) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONGSTATE, "DM must have a default PetscSection in order to create a global PetscSection");
4262: if (!dm->sf) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DM must have a point PetscSF in order to create a global PetscSection");
4263: PetscSectionCreateGlobalSection(s, dm->sf, PETSC_FALSE, PETSC_FALSE, &dm->globalSection);
4264: PetscLayoutDestroy(&dm->map);
4265: PetscSectionGetValueLayout(PetscObjectComm((PetscObject)dm), dm->globalSection, &dm->map);
4266: PetscSectionViewFromOptions(dm->globalSection, NULL, "-global_section_view");
4267: }
4268: *section = dm->globalSection;
4269: return(0);
4270: }
4272: /*@
4273: DMSetGlobalSection - Set the PetscSection encoding the global data layout for the DM.
4275: Input Parameters:
4276: + dm - The DM
4277: - section - The PetscSection, or NULL
4279: Level: intermediate
4281: Note: Any existing Section will be destroyed
4283: .seealso: DMGetGlobalSection(), DMSetLocalSection()
4284: @*/
4285: PetscErrorCode DMSetGlobalSection(DM dm, PetscSection section)
4286: {
4292: PetscObjectReference((PetscObject)section);
4293: PetscSectionDestroy(&dm->globalSection);
4294: dm->globalSection = section;
4295: #if defined(PETSC_USE_DEBUG)
4296: if (section) {DMDefaultSectionCheckConsistency_Internal(dm, dm->localSection, section);}
4297: #endif
4298: return(0);
4299: }
4301: /*@
4302: DMGetSectionSF - Get the PetscSF encoding the parallel dof overlap for the DM. If it has not been set,
4303: it is created from the default PetscSection layouts in the DM.
4305: Input Parameter:
4306: . dm - The DM
4308: Output Parameter:
4309: . sf - The PetscSF
4311: Level: intermediate
4313: Note: This gets a borrowed reference, so the user should not destroy this PetscSF.
4315: .seealso: DMSetSectionSF(), DMCreateSectionSF()
4316: @*/
4317: PetscErrorCode DMGetSectionSF(DM dm, PetscSF *sf)
4318: {
4319: PetscInt nroots;
4325: if (!dm->sectionSF) {
4326: PetscSFCreate(PetscObjectComm((PetscObject)dm),&dm->sectionSF);
4327: }
4328: PetscSFGetGraph(dm->sectionSF, &nroots, NULL, NULL, NULL);
4329: if (nroots < 0) {
4330: PetscSection section, gSection;
4332: DMGetLocalSection(dm, §ion);
4333: if (section) {
4334: DMGetGlobalSection(dm, &gSection);
4335: DMCreateSectionSF(dm, section, gSection);
4336: } else {
4337: *sf = NULL;
4338: return(0);
4339: }
4340: }
4341: *sf = dm->sectionSF;
4342: return(0);
4343: }
4345: /*@
4346: DMSetSectionSF - Set the PetscSF encoding the parallel dof overlap for the DM
4348: Input Parameters:
4349: + dm - The DM
4350: - sf - The PetscSF
4352: Level: intermediate
4354: Note: Any previous SF is destroyed
4356: .seealso: DMGetSectionSF(), DMCreateSectionSF()
4357: @*/
4358: PetscErrorCode DMSetSectionSF(DM dm, PetscSF sf)
4359: {
4365: PetscObjectReference((PetscObject) sf);
4366: PetscSFDestroy(&dm->sectionSF);
4367: dm->sectionSF = sf;
4368: return(0);
4369: }
4371: /*@C
4372: DMCreateSectionSF - Create the PetscSF encoding the parallel dof overlap for the DM based upon the PetscSections
4373: describing the data layout.
4375: Input Parameters:
4376: + dm - The DM
4377: . localSection - PetscSection describing the local data layout
4378: - globalSection - PetscSection describing the global data layout
4380: Notes: One usually uses DMGetSectionSF() to obtain the PetscSF
4382: Level: developer
4384: Developer Note: Since this routine has for arguments the two sections from the DM and puts the resulting PetscSF
4385: directly into the DM, perhaps this function should not take the local and global sections as
4386: input and should just obtain them from the DM?
4388: .seealso: DMGetSectionSF(), DMSetSectionSF(), DMGetLocalSection(), DMGetGlobalSection()
4389: @*/
4390: PetscErrorCode DMCreateSectionSF(DM dm, PetscSection localSection, PetscSection globalSection)
4391: {
4392: MPI_Comm comm;
4393: PetscLayout layout;
4394: const PetscInt *ranges;
4395: PetscInt *local;
4396: PetscSFNode *remote;
4397: PetscInt pStart, pEnd, p, nroots, nleaves = 0, l;
4398: PetscMPIInt size, rank;
4403: PetscObjectGetComm((PetscObject)dm,&comm);
4404: MPI_Comm_size(comm, &size);
4405: MPI_Comm_rank(comm, &rank);
4406: PetscSectionGetChart(globalSection, &pStart, &pEnd);
4407: PetscSectionGetConstrainedStorageSize(globalSection, &nroots);
4408: PetscLayoutCreate(comm, &layout);
4409: PetscLayoutSetBlockSize(layout, 1);
4410: PetscLayoutSetLocalSize(layout, nroots);
4411: PetscLayoutSetUp(layout);
4412: PetscLayoutGetRanges(layout, &ranges);
4413: for (p = pStart; p < pEnd; ++p) {
4414: PetscInt gdof, gcdof;
4416: PetscSectionGetDof(globalSection, p, &gdof);
4417: PetscSectionGetConstraintDof(globalSection, p, &gcdof);
4418: if (gcdof > (gdof < 0 ? -(gdof+1) : gdof)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %d has %d constraints > %d dof", p, gcdof, (gdof < 0 ? -(gdof+1) : gdof));
4419: nleaves += gdof < 0 ? -(gdof+1)-gcdof : gdof-gcdof;
4420: }
4421: PetscMalloc1(nleaves, &local);
4422: PetscMalloc1(nleaves, &remote);
4423: for (p = pStart, l = 0; p < pEnd; ++p) {
4424: const PetscInt *cind;
4425: PetscInt dof, cdof, off, gdof, gcdof, goff, gsize, d, c;
4427: PetscSectionGetDof(localSection, p, &dof);
4428: PetscSectionGetOffset(localSection, p, &off);
4429: PetscSectionGetConstraintDof(localSection, p, &cdof);
4430: PetscSectionGetConstraintIndices(localSection, p, &cind);
4431: PetscSectionGetDof(globalSection, p, &gdof);
4432: PetscSectionGetConstraintDof(globalSection, p, &gcdof);
4433: PetscSectionGetOffset(globalSection, p, &goff);
4434: if (!gdof) continue; /* Censored point */
4435: gsize = gdof < 0 ? -(gdof+1)-gcdof : gdof-gcdof;
4436: if (gsize != dof-cdof) {
4437: if (gsize != dof) SETERRQ4(comm, PETSC_ERR_ARG_WRONG, "Global dof %d for point %d is neither the constrained size %d, nor the unconstrained %d", gsize, p, dof-cdof, dof);
4438: cdof = 0; /* Ignore constraints */
4439: }
4440: for (d = 0, c = 0; d < dof; ++d) {
4441: if ((c < cdof) && (cind[c] == d)) {++c; continue;}
4442: local[l+d-c] = off+d;
4443: }
4444: if (gdof < 0) {
4445: for (d = 0; d < gsize; ++d, ++l) {
4446: PetscInt offset = -(goff+1) + d, r;
4448: PetscFindInt(offset,size+1,ranges,&r);
4449: if (r < 0) r = -(r+2);
4450: if ((r < 0) || (r >= size)) SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %d mapped to invalid process %d (%d, %d)", p, r, gdof, goff);
4451: remote[l].rank = r;
4452: remote[l].index = offset - ranges[r];
4453: }
4454: } else {
4455: for (d = 0; d < gsize; ++d, ++l) {
4456: remote[l].rank = rank;
4457: remote[l].index = goff+d - ranges[rank];
4458: }
4459: }
4460: }
4461: if (l != nleaves) SETERRQ2(comm, PETSC_ERR_PLIB, "Iteration error, l %d != nleaves %d", l, nleaves);
4462: PetscLayoutDestroy(&layout);
4463: PetscSFSetGraph(dm->sectionSF, nroots, nleaves, local, PETSC_OWN_POINTER, remote, PETSC_OWN_POINTER);
4464: return(0);
4465: }
4467: /*@
4468: DMGetPointSF - Get the PetscSF encoding the parallel section point overlap for the DM.
4470: Input Parameter:
4471: . dm - The DM
4473: Output Parameter:
4474: . sf - The PetscSF
4476: Level: intermediate
4478: Note: This gets a borrowed reference, so the user should not destroy this PetscSF.
4480: .seealso: DMSetPointSF(), DMGetSectionSF(), DMSetSectionSF(), DMCreateSectionSF()
4481: @*/
4482: PetscErrorCode DMGetPointSF(DM dm, PetscSF *sf)
4483: {
4487: *sf = dm->sf;
4488: return(0);
4489: }
4491: /*@
4492: DMSetPointSF - Set the PetscSF encoding the parallel section point overlap for the DM.
4494: Input Parameters:
4495: + dm - The DM
4496: - sf - The PetscSF
4498: Level: intermediate
4500: .seealso: DMGetPointSF(), DMGetSectionSF(), DMSetSectionSF(), DMCreateSectionSF()
4501: @*/
4502: PetscErrorCode DMSetPointSF(DM dm, PetscSF sf)
4503: {
4509: PetscObjectReference((PetscObject) sf);
4510: PetscSFDestroy(&dm->sf);
4511: dm->sf = sf;
4512: return(0);
4513: }
4515: static PetscErrorCode DMSetDefaultAdjacency_Private(DM dm, PetscInt f, PetscObject disc)
4516: {
4517: PetscClassId id;
4521: PetscObjectGetClassId(disc, &id);
4522: if (id == PETSCFE_CLASSID) {
4523: DMSetAdjacency(dm, f, PETSC_FALSE, PETSC_TRUE);
4524: } else if (id == PETSCFV_CLASSID) {
4525: DMSetAdjacency(dm, f, PETSC_TRUE, PETSC_FALSE);
4526: } else {
4527: DMSetAdjacency(dm, f, PETSC_FALSE, PETSC_TRUE);
4528: }
4529: return(0);
4530: }
4532: static PetscErrorCode DMFieldEnlarge_Static(DM dm, PetscInt NfNew)
4533: {
4534: RegionField *tmpr;
4535: PetscInt Nf = dm->Nf, f;
4539: if (Nf >= NfNew) return(0);
4540: PetscMalloc1(NfNew, &tmpr);
4541: for (f = 0; f < Nf; ++f) tmpr[f] = dm->fields[f];
4542: for (f = Nf; f < NfNew; ++f) {tmpr[f].disc = NULL; tmpr[f].label = NULL;}
4543: PetscFree(dm->fields);
4544: dm->Nf = NfNew;
4545: dm->fields = tmpr;
4546: return(0);
4547: }
4549: /*@
4550: DMClearFields - Remove all fields from the DM
4552: Logically collective on dm
4554: Input Parameter:
4555: . dm - The DM
4557: Level: intermediate
4559: .seealso: DMGetNumFields(), DMSetNumFields(), DMSetField()
4560: @*/
4561: PetscErrorCode DMClearFields(DM dm)
4562: {
4563: PetscInt f;
4568: for (f = 0; f < dm->Nf; ++f) {
4569: PetscObjectDestroy(&dm->fields[f].disc);
4570: DMLabelDestroy(&dm->fields[f].label);
4571: }
4572: PetscFree(dm->fields);
4573: dm->fields = NULL;
4574: dm->Nf = 0;
4575: return(0);
4576: }
4578: /*@
4579: DMGetNumFields - Get the number of fields in the DM
4581: Not collective
4583: Input Parameter:
4584: . dm - The DM
4586: Output Parameter:
4587: . Nf - The number of fields
4589: Level: intermediate
4591: .seealso: DMSetNumFields(), DMSetField()
4592: @*/
4593: PetscErrorCode DMGetNumFields(DM dm, PetscInt *numFields)
4594: {
4598: *numFields = dm->Nf;
4599: return(0);
4600: }
4602: /*@
4603: DMSetNumFields - Set the number of fields in the DM
4605: Logically collective on dm
4607: Input Parameters:
4608: + dm - The DM
4609: - Nf - The number of fields
4611: Level: intermediate
4613: .seealso: DMGetNumFields(), DMSetField()
4614: @*/
4615: PetscErrorCode DMSetNumFields(DM dm, PetscInt numFields)
4616: {
4617: PetscInt Nf, f;
4622: DMGetNumFields(dm, &Nf);
4623: for (f = Nf; f < numFields; ++f) {
4624: PetscContainer obj;
4626: PetscContainerCreate(PetscObjectComm((PetscObject) dm), &obj);
4627: DMAddField(dm, NULL, (PetscObject) obj);
4628: PetscContainerDestroy(&obj);
4629: }
4630: return(0);
4631: }
4633: /*@
4634: DMGetField - Return the discretization object for a given DM field
4636: Not collective
4638: Input Parameters:
4639: + dm - The DM
4640: - f - The field number
4642: Output Parameters:
4643: + label - The label indicating the support of the field, or NULL for the entire mesh
4644: - field - The discretization object
4646: Level: intermediate
4648: .seealso: DMAddField(), DMSetField()
4649: @*/
4650: PetscErrorCode DMGetField(DM dm, PetscInt f, DMLabel *label, PetscObject *field)
4651: {
4655: if ((f < 0) || (f >= dm->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, dm->Nf);
4656: if (label) *label = dm->fields[f].label;
4657: if (field) *field = dm->fields[f].disc;
4658: return(0);
4659: }
4661: /*@
4662: DMSetField - Set the discretization object for a given DM field
4664: Logically collective on dm
4666: Input Parameters:
4667: + dm - The DM
4668: . f - The field number
4669: . label - The label indicating the support of the field, or NULL for the entire mesh
4670: - field - The discretization object
4672: Level: intermediate
4674: .seealso: DMAddField(), DMGetField()
4675: @*/
4676: PetscErrorCode DMSetField(DM dm, PetscInt f, DMLabel label, PetscObject field)
4677: {
4684: if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
4685: DMFieldEnlarge_Static(dm, f+1);
4686: DMLabelDestroy(&dm->fields[f].label);
4687: PetscObjectDestroy(&dm->fields[f].disc);
4688: dm->fields[f].label = label;
4689: dm->fields[f].disc = field;
4690: PetscObjectReference((PetscObject) label);
4691: PetscObjectReference((PetscObject) field);
4692: DMSetDefaultAdjacency_Private(dm, f, field);
4693: DMClearDS(dm);
4694: return(0);
4695: }
4697: /*@
4698: DMAddField - Add the discretization object for the given DM field
4700: Logically collective on dm
4702: Input Parameters:
4703: + dm - The DM
4704: . label - The label indicating the support of the field, or NULL for the entire mesh
4705: - field - The discretization object
4707: Level: intermediate
4709: .seealso: DMSetField(), DMGetField()
4710: @*/
4711: PetscErrorCode DMAddField(DM dm, DMLabel label, PetscObject field)
4712: {
4713: PetscInt Nf = dm->Nf;
4720: DMFieldEnlarge_Static(dm, Nf+1);
4721: dm->fields[Nf].label = label;
4722: dm->fields[Nf].disc = field;
4723: PetscObjectReference((PetscObject) label);
4724: PetscObjectReference((PetscObject) field);
4725: DMSetDefaultAdjacency_Private(dm, Nf, field);
4726: DMClearDS(dm);
4727: return(0);
4728: }
4730: /*@
4731: DMCopyFields - Copy the discretizations for the DM into another DM
4733: Collective on dm
4735: Input Parameter:
4736: . dm - The DM
4738: Output Parameter:
4739: . newdm - The DM
4741: Level: advanced
4743: .seealso: DMGetField(), DMSetField(), DMAddField(), DMCopyDS(), DMGetDS(), DMGetCellDS()
4744: @*/
4745: PetscErrorCode DMCopyFields(DM dm, DM newdm)
4746: {
4747: PetscInt Nf, f;
4751: if (dm == newdm) return(0);
4752: DMGetNumFields(dm, &Nf);
4753: DMClearFields(newdm);
4754: for (f = 0; f < Nf; ++f) {
4755: DMLabel label;
4756: PetscObject field;
4757: PetscBool useCone, useClosure;
4759: DMGetField(dm, f, &label, &field);
4760: DMSetField(newdm, f, label, field);
4761: DMGetAdjacency(dm, f, &useCone, &useClosure);
4762: DMSetAdjacency(newdm, f, useCone, useClosure);
4763: }
4764: return(0);
4765: }
4767: /*@
4768: DMGetAdjacency - Returns the flags for determining variable influence
4770: Not collective
4772: Input Parameters:
4773: + dm - The DM object
4774: - f - The field number, or PETSC_DEFAULT for the default adjacency
4776: Output Parameter:
4777: + useCone - Flag for variable influence starting with the cone operation
4778: - useClosure - Flag for variable influence using transitive closure
4780: Notes:
4781: $ FEM: Two points p and q are adjacent if q \in closure(star(p)), useCone = PETSC_FALSE, useClosure = PETSC_TRUE
4782: $ FVM: Two points p and q are adjacent if q \in support(p+cone(p)), useCone = PETSC_TRUE, useClosure = PETSC_FALSE
4783: $ FVM++: Two points p and q are adjacent if q \in star(closure(p)), useCone = PETSC_TRUE, useClosure = PETSC_TRUE
4784: Further explanation can be found in the User's Manual Section on the Influence of Variables on One Another.
4786: Level: developer
4788: .seealso: DMSetAdjacency(), DMGetField(), DMSetField()
4789: @*/
4790: PetscErrorCode DMGetAdjacency(DM dm, PetscInt f, PetscBool *useCone, PetscBool *useClosure)
4791: {
4796: if (f < 0) {
4797: if (useCone) *useCone = dm->adjacency[0];
4798: if (useClosure) *useClosure = dm->adjacency[1];
4799: } else {
4800: PetscInt Nf;
4803: DMGetNumFields(dm, &Nf);
4804: if (f >= Nf) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, Nf);
4805: if (useCone) *useCone = dm->fields[f].adjacency[0];
4806: if (useClosure) *useClosure = dm->fields[f].adjacency[1];
4807: }
4808: return(0);
4809: }
4811: /*@
4812: DMSetAdjacency - Set the flags for determining variable influence
4814: Not collective
4816: Input Parameters:
4817: + dm - The DM object
4818: . f - The field number
4819: . useCone - Flag for variable influence starting with the cone operation
4820: - useClosure - Flag for variable influence using transitive closure
4822: Notes:
4823: $ FEM: Two points p and q are adjacent if q \in closure(star(p)), useCone = PETSC_FALSE, useClosure = PETSC_TRUE
4824: $ FVM: Two points p and q are adjacent if q \in support(p+cone(p)), useCone = PETSC_TRUE, useClosure = PETSC_FALSE
4825: $ FVM++: Two points p and q are adjacent if q \in star(closure(p)), useCone = PETSC_TRUE, useClosure = PETSC_TRUE
4826: Further explanation can be found in the User's Manual Section on the Influence of Variables on One Another.
4828: Level: developer
4830: .seealso: DMGetAdjacency(), DMGetField(), DMSetField()
4831: @*/
4832: PetscErrorCode DMSetAdjacency(DM dm, PetscInt f, PetscBool useCone, PetscBool useClosure)
4833: {
4836: if (f < 0) {
4837: dm->adjacency[0] = useCone;
4838: dm->adjacency[1] = useClosure;
4839: } else {
4840: PetscInt Nf;
4843: DMGetNumFields(dm, &Nf);
4844: if (f >= Nf) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, Nf);
4845: dm->fields[f].adjacency[0] = useCone;
4846: dm->fields[f].adjacency[1] = useClosure;
4847: }
4848: return(0);
4849: }
4851: /*@
4852: DMGetBasicAdjacency - Returns the flags for determining variable influence, using either the default or field 0 if it is defined
4854: Not collective
4856: Input Parameters:
4857: . dm - The DM object
4859: Output Parameter:
4860: + useCone - Flag for variable influence starting with the cone operation
4861: - useClosure - Flag for variable influence using transitive closure
4863: Notes:
4864: $ FEM: Two points p and q are adjacent if q \in closure(star(p)), useCone = PETSC_FALSE, useClosure = PETSC_TRUE
4865: $ FVM: Two points p and q are adjacent if q \in support(p+cone(p)), useCone = PETSC_TRUE, useClosure = PETSC_FALSE
4866: $ FVM++: Two points p and q are adjacent if q \in star(closure(p)), useCone = PETSC_TRUE, useClosure = PETSC_TRUE
4868: Level: developer
4870: .seealso: DMSetBasicAdjacency(), DMGetField(), DMSetField()
4871: @*/
4872: PetscErrorCode DMGetBasicAdjacency(DM dm, PetscBool *useCone, PetscBool *useClosure)
4873: {
4874: PetscInt Nf;
4881: DMGetNumFields(dm, &Nf);
4882: if (!Nf) {
4883: DMGetAdjacency(dm, PETSC_DEFAULT, useCone, useClosure);
4884: } else {
4885: DMGetAdjacency(dm, 0, useCone, useClosure);
4886: }
4887: return(0);
4888: }
4890: /*@
4891: DMSetBasicAdjacency - Set the flags for determining variable influence, using either the default or field 0 if it is defined
4893: Not collective
4895: Input Parameters:
4896: + dm - The DM object
4897: . useCone - Flag for variable influence starting with the cone operation
4898: - useClosure - Flag for variable influence using transitive closure
4900: Notes:
4901: $ FEM: Two points p and q are adjacent if q \in closure(star(p)), useCone = PETSC_FALSE, useClosure = PETSC_TRUE
4902: $ FVM: Two points p and q are adjacent if q \in support(p+cone(p)), useCone = PETSC_TRUE, useClosure = PETSC_FALSE
4903: $ FVM++: Two points p and q are adjacent if q \in star(closure(p)), useCone = PETSC_TRUE, useClosure = PETSC_TRUE
4905: Level: developer
4907: .seealso: DMGetBasicAdjacency(), DMGetField(), DMSetField()
4908: @*/
4909: PetscErrorCode DMSetBasicAdjacency(DM dm, PetscBool useCone, PetscBool useClosure)
4910: {
4911: PetscInt Nf;
4916: DMGetNumFields(dm, &Nf);
4917: if (!Nf) {
4918: DMSetAdjacency(dm, PETSC_DEFAULT, useCone, useClosure);
4919: } else {
4920: DMSetAdjacency(dm, 0, useCone, useClosure);
4921: }
4922: return(0);
4923: }
4925: static PetscErrorCode DMDSEnlarge_Static(DM dm, PetscInt NdsNew)
4926: {
4927: DMSpace *tmpd;
4928: PetscInt Nds = dm->Nds, s;
4932: if (Nds >= NdsNew) return(0);
4933: PetscMalloc1(NdsNew, &tmpd);
4934: for (s = 0; s < Nds; ++s) tmpd[s] = dm->probs[s];
4935: for (s = Nds; s < NdsNew; ++s) {tmpd[s].ds = NULL; tmpd[s].label = NULL; tmpd[s].fields = NULL;}
4936: PetscFree(dm->probs);
4937: dm->Nds = NdsNew;
4938: dm->probs = tmpd;
4939: return(0);
4940: }
4942: /*@
4943: DMGetNumDS - Get the number of discrete systems in the DM
4945: Not collective
4947: Input Parameter:
4948: . dm - The DM
4950: Output Parameter:
4951: . Nds - The number of PetscDS objects
4953: Level: intermediate
4955: .seealso: DMGetDS(), DMGetCellDS()
4956: @*/
4957: PetscErrorCode DMGetNumDS(DM dm, PetscInt *Nds)
4958: {
4962: *Nds = dm->Nds;
4963: return(0);
4964: }
4966: /*@
4967: DMClearDS - Remove all discrete systems from the DM
4969: Logically collective on dm
4971: Input Parameter:
4972: . dm - The DM
4974: Level: intermediate
4976: .seealso: DMGetNumDS(), DMGetDS(), DMSetField()
4977: @*/
4978: PetscErrorCode DMClearDS(DM dm)
4979: {
4980: PetscInt s;
4985: for (s = 0; s < dm->Nds; ++s) {
4986: PetscDSDestroy(&dm->probs[s].ds);
4987: DMLabelDestroy(&dm->probs[s].label);
4988: ISDestroy(&dm->probs[s].fields);
4989: }
4990: PetscFree(dm->probs);
4991: dm->probs = NULL;
4992: dm->Nds = 0;
4993: return(0);
4994: }
4996: /*@
4997: DMGetDS - Get the default PetscDS
4999: Not collective
5001: Input Parameter:
5002: . dm - The DM
5004: Output Parameter:
5005: . prob - The default PetscDS
5007: Level: intermediate
5009: .seealso: DMGetCellDS(), DMGetRegionDS()
5010: @*/
5011: PetscErrorCode DMGetDS(DM dm, PetscDS *prob)
5012: {
5018: if (dm->Nds <= 0) {
5019: PetscDS ds;
5021: PetscDSCreate(PetscObjectComm((PetscObject) dm), &ds);
5022: DMSetRegionDS(dm, NULL, NULL, ds);
5023: PetscDSDestroy(&ds);
5024: }
5025: *prob = dm->probs[0].ds;
5026: return(0);
5027: }
5029: /*@
5030: DMGetCellDS - Get the PetscDS defined on a given cell
5032: Not collective
5034: Input Parameters:
5035: + dm - The DM
5036: - point - Cell for the DS
5038: Output Parameter:
5039: . prob - The PetscDS defined on the given cell
5041: Level: developer
5043: .seealso: DMGetDS(), DMSetRegionDS()
5044: @*/
5045: PetscErrorCode DMGetCellDS(DM dm, PetscInt point, PetscDS *prob)
5046: {
5047: PetscDS probDef = NULL;
5048: PetscInt s;
5054: *prob = NULL;
5055: for (s = 0; s < dm->Nds; ++s) {
5056: PetscInt val;
5058: if (!dm->probs[s].label) {probDef = dm->probs[s].ds;}
5059: else {
5060: DMLabelGetValue(dm->probs[s].label, point, &val);
5061: if (val >= 0) {*prob = dm->probs[s].ds; break;}
5062: }
5063: }
5064: if (!*prob) *prob = probDef;
5065: return(0);
5066: }
5068: /*@
5069: DMGetRegionDS - Get the PetscDS for a given mesh region, defined by a DMLabel
5071: Not collective
5073: Input Parameters:
5074: + dm - The DM
5075: - label - The DMLabel defining the mesh region, or NULL for the entire mesh
5077: Output Parameters:
5078: + fields - The IS containing the DM field numbers for the fields in this DS, or NULL
5079: - prob - The PetscDS defined on the given region, or NULL
5081: Note: If the label is missing, this function returns an error
5083: Level: advanced
5085: .seealso: DMGetRegionNumDS(), DMSetRegionDS(), DMGetDS(), DMGetCellDS()
5086: @*/
5087: PetscErrorCode DMGetRegionDS(DM dm, DMLabel label, IS *fields, PetscDS *ds)
5088: {
5089: PetscInt Nds = dm->Nds, s;
5096: for (s = 0; s < Nds; ++s) {
5097: if (dm->probs[s].label == label) {
5098: if (fields) *fields = dm->probs[s].fields;
5099: if (ds) *ds = dm->probs[s].ds;
5100: return(0);
5101: }
5102: }
5103: return(0);
5104: }
5106: /*@
5107: DMGetRegionNumDS - Get the PetscDS for a given mesh region, defined by the region number
5109: Not collective
5111: Input Parameters:
5112: + dm - The DM
5113: - num - The region number, in [0, Nds)
5115: Output Parameters:
5116: + label - The region label, or NULL
5117: . fields - The IS containing the DM field numbers for the fields in this DS, or NULL
5118: - prob - The PetscDS defined on the given region, or NULL
5120: Level: advanced
5122: .seealso: DMGetRegionDS(), DMSetRegionDS(), DMGetDS(), DMGetCellDS()
5123: @*/
5124: PetscErrorCode DMGetRegionNumDS(DM dm, PetscInt num, DMLabel *label, IS *fields, PetscDS *ds)
5125: {
5126: PetscInt Nds;
5131: DMGetNumDS(dm, &Nds);
5132: if ((num < 0) || (num >= Nds)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Region number %D is not in [0, %D)", num, Nds);
5133: if (label) {
5135: *label = dm->probs[num].label;
5136: }
5137: if (fields) {
5139: *fields = dm->probs[num].fields;
5140: }
5141: if (ds) {
5143: *ds = dm->probs[num].ds;
5144: }
5145: return(0);
5146: }
5148: /*@
5149: DMSetRegionDS - Set the PetscDS for a given mesh region, defined by a DMLabel
5151: Collective on dm
5153: Input Parameters:
5154: + dm - The DM
5155: . label - The DMLabel defining the mesh region, or NULL for the entire mesh
5156: . fields - The IS containing the DM field numbers for the fields in this DS, or NULL for all fields
5157: - prob - The PetscDS defined on the given cell
5159: Note: If the label has a DS defined, it will be replaced. Otherwise, it will be added to the DM. If DS is replaced,
5160: the fields argument is ignored.
5162: Level: advanced
5164: .seealso: DMGetRegionDS(), DMGetDS(), DMGetCellDS()
5165: @*/
5166: PetscErrorCode DMSetRegionDS(DM dm, DMLabel label, IS fields, PetscDS ds)
5167: {
5168: PetscInt Nds = dm->Nds, s;
5175: for (s = 0; s < Nds; ++s) {
5176: if (dm->probs[s].label == label) {
5177: PetscDSDestroy(&dm->probs[s].ds);
5178: dm->probs[s].ds = ds;
5179: return(0);
5180: }
5181: }
5182: DMDSEnlarge_Static(dm, Nds+1);
5183: PetscObjectReference((PetscObject) label);
5184: PetscObjectReference((PetscObject) fields);
5185: PetscObjectReference((PetscObject) ds);
5186: if (!label) {
5187: /* Put the NULL label at the front, so it is returned as the default */
5188: for (s = Nds-1; s >=0; --s) dm->probs[s+1] = dm->probs[s];
5189: Nds = 0;
5190: }
5191: dm->probs[Nds].label = label;
5192: dm->probs[Nds].fields = fields;
5193: dm->probs[Nds].ds = ds;
5194: return(0);
5195: }
5197: /*@
5198: DMCreateDS - Create the discrete systems for the DM based upon the fields added to the DM
5200: Collective on dm
5202: Input Parameter:
5203: . dm - The DM
5205: Note: If the label has a DS defined, it will be replaced. Otherwise, it will be added to the DM.
5207: Level: intermediate
5209: .seealso: DMSetField, DMAddField(), DMGetDS(), DMGetCellDS(), DMGetRegionDS(), DMSetRegionDS()
5210: @*/
5211: PetscErrorCode DMCreateDS(DM dm)
5212: {
5213: MPI_Comm comm;
5214: PetscDS prob, probh = NULL;
5215: PetscInt dimEmbed, Nf = dm->Nf, f, s, field = 0, fieldh = 0;
5216: PetscBool doSetup = PETSC_TRUE;
5221: if (!dm->fields) return(0);
5222: /* Can only handle two label cases right now:
5223: 1) NULL
5224: 2) Hybrid cells
5225: */
5226: PetscObjectGetComm((PetscObject) dm, &comm);
5227: DMGetCoordinateDim(dm, &dimEmbed);
5228: /* Create default DS */
5229: DMGetRegionDS(dm, NULL, NULL, &prob);
5230: if (!prob) {
5231: IS fields;
5232: PetscInt *fld, nf;
5234: for (f = 0, nf = 0; f < Nf; ++f) if (!dm->fields[f].label) ++nf;
5235: PetscMalloc1(nf, &fld);
5236: for (f = 0, nf = 0; f < Nf; ++f) if (!dm->fields[f].label) fld[nf++] = f;
5237: ISCreate(PETSC_COMM_SELF, &fields);
5238: PetscObjectSetOptionsPrefix((PetscObject) fields, "dm_fields_");
5239: ISSetType(fields, ISGENERAL);
5240: ISGeneralSetIndices(fields, nf, fld, PETSC_OWN_POINTER);
5242: PetscDSCreate(comm, &prob);
5243: DMSetRegionDS(dm, NULL, fields, prob);
5244: PetscDSDestroy(&prob);
5245: ISDestroy(&fields);
5246: DMGetRegionDS(dm, NULL, NULL, &prob);
5247: }
5248: PetscDSSetCoordinateDimension(prob, dimEmbed);
5249: /* Optionally create hybrid DS */
5250: for (f = 0; f < Nf; ++f) {
5251: DMLabel label = dm->fields[f].label;
5252: PetscInt lStart, lEnd;
5254: if (label) {
5255: DM plex;
5256: DMPolytopeType ct;
5257: IS fields;
5258: PetscInt *fld;
5259: PetscInt depth;
5261: DMConvert(dm, DMPLEX, &plex);
5262: DMPlexGetDepth(plex, &depth);
5263: DMDestroy(&plex);
5265: DMLabelGetBounds(label, &lStart, &lEnd);
5266: DMPlexGetCellType(dm, lStart, &ct);
5267: switch (ct) {
5268: case DM_POLYTOPE_POINT_PRISM_TENSOR:
5269: case DM_POLYTOPE_SEG_PRISM_TENSOR:
5270: case DM_POLYTOPE_TRI_PRISM_TENSOR:
5271: case DM_POLYTOPE_QUAD_PRISM_TENSOR:
5272: break;
5273: default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only support labels over tensor prism cells right now");
5274: }
5275: PetscDSCreate(comm, &probh);
5276: PetscMalloc1(1, &fld);
5277: fld[0] = f;
5278: ISCreate(PETSC_COMM_SELF, &fields);
5279: PetscObjectSetOptionsPrefix((PetscObject) fields, "dm_fields_");
5280: ISSetType(fields, ISGENERAL);
5281: ISGeneralSetIndices(fields, 1, fld, PETSC_OWN_POINTER);
5282: DMSetRegionDS(dm, label, fields, probh);
5283: ISDestroy(&fields);
5284: PetscDSSetHybrid(probh, PETSC_TRUE);
5285: PetscDSSetCoordinateDimension(probh, dimEmbed);
5286: break;
5287: }
5288: }
5289: /* Set fields in DSes */
5290: for (f = 0; f < Nf; ++f) {
5291: DMLabel label = dm->fields[f].label;
5292: PetscObject disc = dm->fields[f].disc;
5294: if (!label) {
5295: PetscDSSetDiscretization(prob, field++, disc);
5296: if (probh) {
5297: PetscFE subfe;
5299: PetscFEGetHeightSubspace((PetscFE) disc, 1, &subfe);
5300: PetscDSSetDiscretization(probh, fieldh++, (PetscObject) subfe);
5301: }
5302: } else {
5303: PetscDSSetDiscretization(probh, fieldh++, disc);
5304: }
5305: /* We allow people to have placeholder fields and construct the Section by hand */
5306: {
5307: PetscClassId id;
5309: PetscObjectGetClassId(disc, &id);
5310: if ((id != PETSCFE_CLASSID) && (id != PETSCFV_CLASSID)) doSetup = PETSC_FALSE;
5311: }
5312: }
5313: PetscDSDestroy(&probh);
5314: /* Setup DSes */
5315: if (doSetup) {
5316: for (s = 0; s < dm->Nds; ++s) {PetscDSSetUp(dm->probs[s].ds);}
5317: }
5318: return(0);
5319: }
5321: /*@
5322: DMCopyDS - Copy the discrete systems for the DM into another DM
5324: Collective on dm
5326: Input Parameter:
5327: . dm - The DM
5329: Output Parameter:
5330: . newdm - The DM
5332: Level: advanced
5334: .seealso: DMCopyFields(), DMAddField(), DMGetDS(), DMGetCellDS(), DMGetRegionDS(), DMSetRegionDS()
5335: @*/
5336: PetscErrorCode DMCopyDS(DM dm, DM newdm)
5337: {
5338: PetscInt Nds, s;
5342: if (dm == newdm) return(0);
5343: DMGetNumDS(dm, &Nds);
5344: DMClearDS(newdm);
5345: for (s = 0; s < Nds; ++s) {
5346: DMLabel label;
5347: IS fields;
5348: PetscDS ds;
5350: DMGetRegionNumDS(dm, s, &label, &fields, &ds);
5351: DMSetRegionDS(newdm, label, fields, ds);
5352: }
5353: return(0);
5354: }
5356: /*@
5357: DMCopyDisc - Copy the fields and discrete systems for the DM into another DM
5359: Collective on dm
5361: Input Parameter:
5362: . dm - The DM
5364: Output Parameter:
5365: . newdm - The DM
5367: Level: advanced
5369: .seealso: DMCopyFields(), DMCopyDS()
5370: @*/
5371: PetscErrorCode DMCopyDisc(DM dm, DM newdm)
5372: {
5376: DMCopyFields(dm, newdm);
5377: DMCopyDS(dm, newdm);
5378: return(0);
5379: }
5381: PetscErrorCode DMRestrictHook_Coordinates(DM dm,DM dmc,void *ctx)
5382: {
5383: DM dm_coord,dmc_coord;
5385: Vec coords,ccoords;
5386: Mat inject;
5388: DMGetCoordinateDM(dm,&dm_coord);
5389: DMGetCoordinateDM(dmc,&dmc_coord);
5390: DMGetCoordinates(dm,&coords);
5391: DMGetCoordinates(dmc,&ccoords);
5392: if (coords && !ccoords) {
5393: DMCreateGlobalVector(dmc_coord,&ccoords);
5394: PetscObjectSetName((PetscObject)ccoords,"coordinates");
5395: DMCreateInjection(dmc_coord,dm_coord,&inject);
5396: MatRestrict(inject,coords,ccoords);
5397: MatDestroy(&inject);
5398: DMSetCoordinates(dmc,ccoords);
5399: VecDestroy(&ccoords);
5400: }
5401: return(0);
5402: }
5404: static PetscErrorCode DMSubDomainHook_Coordinates(DM dm,DM subdm,void *ctx)
5405: {
5406: DM dm_coord,subdm_coord;
5408: Vec coords,ccoords,clcoords;
5409: VecScatter *scat_i,*scat_g;
5411: DMGetCoordinateDM(dm,&dm_coord);
5412: DMGetCoordinateDM(subdm,&subdm_coord);
5413: DMGetCoordinates(dm,&coords);
5414: DMGetCoordinates(subdm,&ccoords);
5415: if (coords && !ccoords) {
5416: DMCreateGlobalVector(subdm_coord,&ccoords);
5417: PetscObjectSetName((PetscObject)ccoords,"coordinates");
5418: DMCreateLocalVector(subdm_coord,&clcoords);
5419: PetscObjectSetName((PetscObject)clcoords,"coordinates");
5420: DMCreateDomainDecompositionScatters(dm_coord,1,&subdm_coord,NULL,&scat_i,&scat_g);
5421: VecScatterBegin(scat_i[0],coords,ccoords,INSERT_VALUES,SCATTER_FORWARD);
5422: VecScatterEnd(scat_i[0],coords,ccoords,INSERT_VALUES,SCATTER_FORWARD);
5423: VecScatterBegin(scat_g[0],coords,clcoords,INSERT_VALUES,SCATTER_FORWARD);
5424: VecScatterEnd(scat_g[0],coords,clcoords,INSERT_VALUES,SCATTER_FORWARD);
5425: DMSetCoordinates(subdm,ccoords);
5426: DMSetCoordinatesLocal(subdm,clcoords);
5427: VecScatterDestroy(&scat_i[0]);
5428: VecScatterDestroy(&scat_g[0]);
5429: VecDestroy(&ccoords);
5430: VecDestroy(&clcoords);
5431: PetscFree(scat_i);
5432: PetscFree(scat_g);
5433: }
5434: return(0);
5435: }
5437: /*@
5438: DMGetDimension - Return the topological dimension of the DM
5440: Not collective
5442: Input Parameter:
5443: . dm - The DM
5445: Output Parameter:
5446: . dim - The topological dimension
5448: Level: beginner
5450: .seealso: DMSetDimension(), DMCreate()
5451: @*/
5452: PetscErrorCode DMGetDimension(DM dm, PetscInt *dim)
5453: {
5457: *dim = dm->dim;
5458: return(0);
5459: }
5461: /*@
5462: DMSetDimension - Set the topological dimension of the DM
5464: Collective on dm
5466: Input Parameters:
5467: + dm - The DM
5468: - dim - The topological dimension
5470: Level: beginner
5472: .seealso: DMGetDimension(), DMCreate()
5473: @*/
5474: PetscErrorCode DMSetDimension(DM dm, PetscInt dim)
5475: {
5476: PetscDS ds;
5482: dm->dim = dim;
5483: DMGetDS(dm, &ds);
5484: if (ds->dimEmbed < 0) {PetscDSSetCoordinateDimension(ds, dm->dim);}
5485: return(0);
5486: }
5488: /*@
5489: DMGetDimPoints - Get the half-open interval for all points of a given dimension
5491: Collective on dm
5493: Input Parameters:
5494: + dm - the DM
5495: - dim - the dimension
5497: Output Parameters:
5498: + pStart - The first point of the given dimension
5499: - pEnd - The first point following points of the given dimension
5501: Note:
5502: The points are vertices in the Hasse diagram encoding the topology. This is explained in
5503: https://arxiv.org/abs/0908.4427. If no points exist of this dimension in the storage scheme,
5504: then the interval is empty.
5506: Level: intermediate
5508: .seealso: DMPLEX, DMPlexGetDepthStratum(), DMPlexGetHeightStratum()
5509: @*/
5510: PetscErrorCode DMGetDimPoints(DM dm, PetscInt dim, PetscInt *pStart, PetscInt *pEnd)
5511: {
5512: PetscInt d;
5517: DMGetDimension(dm, &d);
5518: if ((dim < 0) || (dim > d)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension %d 1", dim, d);
5519: if (!dm->ops->getdimpoints) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "DM type %s does not implement DMGetDimPoints",((PetscObject)dm)->type_name);
5520: (*dm->ops->getdimpoints)(dm, dim, pStart, pEnd);
5521: return(0);
5522: }
5524: /*@
5525: DMSetCoordinates - Sets into the DM a global vector that holds the coordinates
5527: Collective on dm
5529: Input Parameters:
5530: + dm - the DM
5531: - c - coordinate vector
5533: Notes:
5534: The coordinates do include those for ghost points, which are in the local vector.
5536: The vector c should be destroyed by the caller.
5538: Level: intermediate
5540: .seealso: DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM()
5541: @*/
5542: PetscErrorCode DMSetCoordinates(DM dm, Vec c)
5543: {
5549: PetscObjectReference((PetscObject) c);
5550: VecDestroy(&dm->coordinates);
5551: dm->coordinates = c;
5552: VecDestroy(&dm->coordinatesLocal);
5553: DMCoarsenHookAdd(dm,DMRestrictHook_Coordinates,NULL,NULL);
5554: DMSubDomainHookAdd(dm,DMSubDomainHook_Coordinates,NULL,NULL);
5555: return(0);
5556: }
5558: /*@
5559: DMSetCoordinatesLocal - Sets into the DM a local vector that holds the coordinates
5561: Not collective
5563: Input Parameters:
5564: + dm - the DM
5565: - c - coordinate vector
5567: Notes:
5568: The coordinates of ghost points can be set using DMSetCoordinates()
5569: followed by DMGetCoordinatesLocal(). This is intended to enable the
5570: setting of ghost coordinates outside of the domain.
5572: The vector c should be destroyed by the caller.
5574: Level: intermediate
5576: .seealso: DMGetCoordinatesLocal(), DMSetCoordinates(), DMGetCoordinates(), DMGetCoordinateDM()
5577: @*/
5578: PetscErrorCode DMSetCoordinatesLocal(DM dm, Vec c)
5579: {
5585: PetscObjectReference((PetscObject) c);
5586: VecDestroy(&dm->coordinatesLocal);
5588: dm->coordinatesLocal = c;
5590: VecDestroy(&dm->coordinates);
5591: return(0);
5592: }
5594: /*@
5595: DMGetCoordinates - Gets a global vector with the coordinates associated with the DM.
5597: Collective on dm
5599: Input Parameter:
5600: . dm - the DM
5602: Output Parameter:
5603: . c - global coordinate vector
5605: Note:
5606: This is a borrowed reference, so the user should NOT destroy this vector
5608: Each process has only the local coordinates (does NOT have the ghost coordinates).
5610: For DMDA, in two and three dimensions coordinates are interlaced (x_0,y_0,x_1,y_1,...)
5611: and (x_0,y_0,z_0,x_1,y_1,z_1...)
5613: Level: intermediate
5615: .seealso: DMSetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM()
5616: @*/
5617: PetscErrorCode DMGetCoordinates(DM dm, Vec *c)
5618: {
5624: if (!dm->coordinates && dm->coordinatesLocal) {
5625: DM cdm = NULL;
5626: PetscBool localized;
5628: DMGetCoordinateDM(dm, &cdm);
5629: DMCreateGlobalVector(cdm, &dm->coordinates);
5630: DMGetCoordinatesLocalized(dm, &localized);
5631: /* Block size is not correctly set by CreateGlobalVector() if coordinates are localized */
5632: if (localized) {
5633: PetscInt cdim;
5635: DMGetCoordinateDim(dm, &cdim);
5636: VecSetBlockSize(dm->coordinates, cdim);
5637: }
5638: PetscObjectSetName((PetscObject) dm->coordinates, "coordinates");
5639: DMLocalToGlobalBegin(cdm, dm->coordinatesLocal, INSERT_VALUES, dm->coordinates);
5640: DMLocalToGlobalEnd(cdm, dm->coordinatesLocal, INSERT_VALUES, dm->coordinates);
5641: }
5642: *c = dm->coordinates;
5643: return(0);
5644: }
5646: /*@
5647: DMGetCoordinatesLocalSetUp - Prepares a local vector of coordinates, so that DMGetCoordinatesLocalNoncollective() can be used as non-collective afterwards.
5649: Collective on dm
5651: Input Parameter:
5652: . dm - the DM
5654: Level: advanced
5656: .seealso: DMGetCoordinatesLocalNoncollective()
5657: @*/
5658: PetscErrorCode DMGetCoordinatesLocalSetUp(DM dm)
5659: {
5664: if (!dm->coordinatesLocal && dm->coordinates) {
5665: DM cdm = NULL;
5666: PetscBool localized;
5668: DMGetCoordinateDM(dm, &cdm);
5669: DMCreateLocalVector(cdm, &dm->coordinatesLocal);
5670: DMGetCoordinatesLocalized(dm, &localized);
5671: /* Block size is not correctly set by CreateLocalVector() if coordinates are localized */
5672: if (localized) {
5673: PetscInt cdim;
5675: DMGetCoordinateDim(dm, &cdim);
5676: VecSetBlockSize(dm->coordinates, cdim);
5677: }
5678: PetscObjectSetName((PetscObject) dm->coordinatesLocal, "coordinates");
5679: DMGlobalToLocalBegin(cdm, dm->coordinates, INSERT_VALUES, dm->coordinatesLocal);
5680: DMGlobalToLocalEnd(cdm, dm->coordinates, INSERT_VALUES, dm->coordinatesLocal);
5681: }
5682: return(0);
5683: }
5685: /*@
5686: DMGetCoordinatesLocal - Gets a local vector with the coordinates associated with the DM.
5688: Collective on dm
5690: Input Parameter:
5691: . dm - the DM
5693: Output Parameter:
5694: . c - coordinate vector
5696: Note:
5697: This is a borrowed reference, so the user should NOT destroy this vector
5699: Each process has the local and ghost coordinates
5701: For DMDA, in two and three dimensions coordinates are interlaced (x_0,y_0,x_1,y_1,...)
5702: and (x_0,y_0,z_0,x_1,y_1,z_1...)
5704: Level: intermediate
5706: .seealso: DMSetCoordinatesLocal(), DMGetCoordinates(), DMSetCoordinates(), DMGetCoordinateDM(), DMGetCoordinatesLocalNoncollective()
5707: @*/
5708: PetscErrorCode DMGetCoordinatesLocal(DM dm, Vec *c)
5709: {
5715: DMGetCoordinatesLocalSetUp(dm);
5716: *c = dm->coordinatesLocal;
5717: return(0);
5718: }
5720: /*@
5721: DMGetCoordinatesLocalNoncollective - Non-collective version of DMGetCoordinatesLocal(). Fails if global coordinates have been set and DMGetCoordinatesLocalSetUp() not called.
5723: Not collective
5725: Input Parameter:
5726: . dm - the DM
5728: Output Parameter:
5729: . c - coordinate vector
5731: Level: advanced
5733: .seealso: DMGetCoordinatesLocalSetUp(), DMGetCoordinatesLocal(), DMSetCoordinatesLocal(), DMGetCoordinates(), DMSetCoordinates(), DMGetCoordinateDM()
5734: @*/
5735: PetscErrorCode DMGetCoordinatesLocalNoncollective(DM dm, Vec *c)
5736: {
5740: if (!dm->coordinatesLocal && dm->coordinates) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DMGetCoordinatesLocalSetUp() has not been called");
5741: *c = dm->coordinatesLocal;
5742: return(0);
5743: }
5745: /*@
5746: DMGetCoordinatesLocalTuple - Gets a local vector with the coordinates of specified points and section describing its layout.
5748: Not collective
5750: Input Parameter:
5751: + dm - the DM
5752: - p - the IS of points whose coordinates will be returned
5754: Output Parameter:
5755: + pCoordSection - the PetscSection describing the layout of pCoord, i.e. each point corresponds to one point in p, and DOFs correspond to coordinates
5756: - pCoord - the Vec with coordinates of points in p
5758: Note:
5759: DMGetCoordinatesLocalSetUp() must be called first. This function employs DMGetCoordinatesLocalNoncollective() so it is not collective.
5761: This creates a new vector, so the user SHOULD destroy this vector
5763: Each process has the local and ghost coordinates
5765: For DMDA, in two and three dimensions coordinates are interlaced (x_0,y_0,x_1,y_1,...)
5766: and (x_0,y_0,z_0,x_1,y_1,z_1...)
5768: Level: advanced
5770: .seealso: DMSetCoordinatesLocal(), DMGetCoordinatesLocal(), DMGetCoordinatesLocalNoncollective(), DMGetCoordinatesLocalSetUp(), DMGetCoordinates(), DMSetCoordinates(), DMGetCoordinateDM()
5771: @*/
5772: PetscErrorCode DMGetCoordinatesLocalTuple(DM dm, IS p, PetscSection *pCoordSection, Vec *pCoord)
5773: {
5774: PetscSection cs, newcs;
5775: Vec coords;
5776: const PetscScalar *arr;
5777: PetscScalar *newarr=NULL;
5778: PetscInt n;
5779: PetscErrorCode ierr;
5786: if (!dm->coordinatesLocal) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DMGetCoordinatesLocalSetUp() has not been called or coordinates not set");
5787: if (!dm->coordinateDM || !dm->coordinateDM->localSection) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DM not supported");
5788: cs = dm->coordinateDM->localSection;
5789: coords = dm->coordinatesLocal;
5790: VecGetArrayRead(coords, &arr);
5791: PetscSectionExtractDofsFromArray(cs, MPIU_SCALAR, arr, p, &newcs, pCoord ? ((void**)&newarr) : NULL);
5792: VecRestoreArrayRead(coords, &arr);
5793: if (pCoord) {
5794: PetscSectionGetStorageSize(newcs, &n);
5795: /* set array in two steps to mimic PETSC_OWN_POINTER */
5796: VecCreateSeqWithArray(PetscObjectComm((PetscObject)p), 1, n, NULL, pCoord);
5797: VecReplaceArray(*pCoord, newarr);
5798: } else {
5799: PetscFree(newarr);
5800: }
5801: if (pCoordSection) {*pCoordSection = newcs;}
5802: else {PetscSectionDestroy(&newcs);}
5803: return(0);
5804: }
5806: PetscErrorCode DMGetCoordinateField(DM dm, DMField *field)
5807: {
5813: if (!dm->coordinateField) {
5814: if (dm->ops->createcoordinatefield) {
5815: (*dm->ops->createcoordinatefield)(dm,&dm->coordinateField);
5816: }
5817: }
5818: *field = dm->coordinateField;
5819: return(0);
5820: }
5822: PetscErrorCode DMSetCoordinateField(DM dm, DMField field)
5823: {
5829: PetscObjectReference((PetscObject)field);
5830: DMFieldDestroy(&dm->coordinateField);
5831: dm->coordinateField = field;
5832: return(0);
5833: }
5835: /*@
5836: DMGetCoordinateDM - Gets the DM that prescribes coordinate layout and scatters between global and local coordinates
5838: Collective on dm
5840: Input Parameter:
5841: . dm - the DM
5843: Output Parameter:
5844: . cdm - coordinate DM
5846: Level: intermediate
5848: .seealso: DMSetCoordinateDM(), DMSetCoordinates(), DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLocal()
5849: @*/
5850: PetscErrorCode DMGetCoordinateDM(DM dm, DM *cdm)
5851: {
5857: if (!dm->coordinateDM) {
5858: DM cdm;
5860: if (!dm->ops->createcoordinatedm) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unable to create coordinates for this DM");
5861: (*dm->ops->createcoordinatedm)(dm, &cdm);
5862: /* Just in case the DM sets the coordinate DM when creating it (DMP4est can do this, because it may not setup
5863: * until the call to CreateCoordinateDM) */
5864: DMDestroy(&dm->coordinateDM);
5865: dm->coordinateDM = cdm;
5866: }
5867: *cdm = dm->coordinateDM;
5868: return(0);
5869: }
5871: /*@
5872: DMSetCoordinateDM - Sets the DM that prescribes coordinate layout and scatters between global and local coordinates
5874: Logically Collective on dm
5876: Input Parameters:
5877: + dm - the DM
5878: - cdm - coordinate DM
5880: Level: intermediate
5882: .seealso: DMGetCoordinateDM(), DMSetCoordinates(), DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLocal()
5883: @*/
5884: PetscErrorCode DMSetCoordinateDM(DM dm, DM cdm)
5885: {
5891: PetscObjectReference((PetscObject)cdm);
5892: DMDestroy(&dm->coordinateDM);
5893: dm->coordinateDM = cdm;
5894: return(0);
5895: }
5897: /*@
5898: DMGetCoordinateDim - Retrieve the dimension of embedding space for coordinate values.
5900: Not Collective
5902: Input Parameter:
5903: . dm - The DM object
5905: Output Parameter:
5906: . dim - The embedding dimension
5908: Level: intermediate
5910: .seealso: DMSetCoordinateDim(), DMGetCoordinateSection(), DMGetCoordinateDM(), DMGetLocalSection(), DMSetLocalSection()
5911: @*/
5912: PetscErrorCode DMGetCoordinateDim(DM dm, PetscInt *dim)
5913: {
5917: if (dm->dimEmbed == PETSC_DEFAULT) {
5918: dm->dimEmbed = dm->dim;
5919: }
5920: *dim = dm->dimEmbed;
5921: return(0);
5922: }
5924: /*@
5925: DMSetCoordinateDim - Set the dimension of the embedding space for coordinate values.
5927: Not Collective
5929: Input Parameters:
5930: + dm - The DM object
5931: - dim - The embedding dimension
5933: Level: intermediate
5935: .seealso: DMGetCoordinateDim(), DMSetCoordinateSection(), DMGetCoordinateSection(), DMGetLocalSection(), DMSetLocalSection()
5936: @*/
5937: PetscErrorCode DMSetCoordinateDim(DM dm, PetscInt dim)
5938: {
5939: PetscDS ds;
5944: dm->dimEmbed = dim;
5945: DMGetDS(dm, &ds);
5946: PetscDSSetCoordinateDimension(ds, dim);
5947: return(0);
5948: }
5950: /*@
5951: DMGetCoordinateSection - Retrieve the layout of coordinate values over the mesh.
5953: Collective on dm
5955: Input Parameter:
5956: . dm - The DM object
5958: Output Parameter:
5959: . section - The PetscSection object
5961: Level: intermediate
5963: .seealso: DMGetCoordinateDM(), DMGetLocalSection(), DMSetLocalSection()
5964: @*/
5965: PetscErrorCode DMGetCoordinateSection(DM dm, PetscSection *section)
5966: {
5967: DM cdm;
5973: DMGetCoordinateDM(dm, &cdm);
5974: DMGetLocalSection(cdm, section);
5975: return(0);
5976: }
5978: /*@
5979: DMSetCoordinateSection - Set the layout of coordinate values over the mesh.
5981: Not Collective
5983: Input Parameters:
5984: + dm - The DM object
5985: . dim - The embedding dimension, or PETSC_DETERMINE
5986: - section - The PetscSection object
5988: Level: intermediate
5990: .seealso: DMGetCoordinateSection(), DMGetLocalSection(), DMSetLocalSection()
5991: @*/
5992: PetscErrorCode DMSetCoordinateSection(DM dm, PetscInt dim, PetscSection section)
5993: {
5994: DM cdm;
6000: DMGetCoordinateDM(dm, &cdm);
6001: DMSetLocalSection(cdm, section);
6002: if (dim == PETSC_DETERMINE) {
6003: PetscInt d = PETSC_DEFAULT;
6004: PetscInt pStart, pEnd, vStart, vEnd, v, dd;
6006: PetscSectionGetChart(section, &pStart, &pEnd);
6007: DMGetDimPoints(dm, 0, &vStart, &vEnd);
6008: pStart = PetscMax(vStart, pStart);
6009: pEnd = PetscMin(vEnd, pEnd);
6010: for (v = pStart; v < pEnd; ++v) {
6011: PetscSectionGetDof(section, v, &dd);
6012: if (dd) {d = dd; break;}
6013: }
6014: if (d >= 0) {DMSetCoordinateDim(dm, d);}
6015: }
6016: return(0);
6017: }
6019: /*@C
6020: DMGetPeriodicity - Get the description of mesh periodicity
6022: Input Parameters:
6023: . dm - The DM object
6025: Output Parameters:
6026: + per - Whether the DM is periodic or not
6027: . maxCell - Over distances greater than this, we can assume a point has crossed over to another sheet, when trying to localize cell coordinates
6028: . L - If we assume the mesh is a torus, this is the length of each coordinate
6029: - bd - This describes the type of periodicity in each topological dimension
6031: Level: developer
6033: .seealso: DMGetPeriodicity()
6034: @*/
6035: PetscErrorCode DMGetPeriodicity(DM dm, PetscBool *per, const PetscReal **maxCell, const PetscReal **L, const DMBoundaryType **bd)
6036: {
6039: if (per) *per = dm->periodic;
6040: if (L) *L = dm->L;
6041: if (maxCell) *maxCell = dm->maxCell;
6042: if (bd) *bd = dm->bdtype;
6043: return(0);
6044: }
6046: /*@C
6047: DMSetPeriodicity - Set the description of mesh periodicity
6049: Input Parameters:
6050: + dm - The DM object
6051: . per - Whether the DM is periodic or not.
6052: . maxCell - Over distances greater than this, we can assume a point has crossed over to another sheet, when trying to localize cell coordinates. Pass NULL to remove such information.
6053: . L - If we assume the mesh is a torus, this is the length of each coordinate
6054: - bd - This describes the type of periodicity in each topological dimension
6056: Notes: If per is PETSC_TRUE and maxCell is not provided, coordinates need to be already localized, or must be localized by hand by the user.
6058: Level: developer
6060: .seealso: DMGetPeriodicity()
6061: @*/
6062: PetscErrorCode DMSetPeriodicity(DM dm, PetscBool per, const PetscReal maxCell[], const PetscReal L[], const DMBoundaryType bd[])
6063: {
6064: PetscInt dim, d;
6073: DMGetDimension(dm, &dim);
6074: if (maxCell) {
6075: if (!dm->maxCell) {PetscMalloc1(dim, &dm->maxCell);}
6076: for (d = 0; d < dim; ++d) dm->maxCell[d] = maxCell[d];
6077: } else { /* remove maxCell information to disable automatic computation of localized vertices */
6078: PetscFree(dm->maxCell);
6079: }
6081: if (L) {
6082: if (!dm->L) {PetscMalloc1(dim, &dm->L);}
6083: for (d = 0; d < dim; ++d) dm->L[d] = L[d];
6084: }
6085: if (bd) {
6086: if (!dm->bdtype) {PetscMalloc1(dim, &dm->bdtype);}
6087: for (d = 0; d < dim; ++d) dm->bdtype[d] = bd[d];
6088: }
6089: dm->periodic = per;
6090: return(0);
6091: }
6093: /*@
6094: DMLocalizeCoordinate - If a mesh is periodic (a torus with lengths L_i, some of which can be infinite), project the coordinate onto [0, L_i) in each dimension.
6096: Input Parameters:
6097: + dm - The DM
6098: . in - The input coordinate point (dim numbers)
6099: - endpoint - Include the endpoint L_i
6101: Output Parameter:
6102: . out - The localized coordinate point
6104: Level: developer
6106: .seealso: DMLocalizeCoordinates(), DMLocalizeAddCoordinate()
6107: @*/
6108: PetscErrorCode DMLocalizeCoordinate(DM dm, const PetscScalar in[], PetscBool endpoint, PetscScalar out[])
6109: {
6110: PetscInt dim, d;
6114: DMGetCoordinateDim(dm, &dim);
6115: if (!dm->maxCell) {
6116: for (d = 0; d < dim; ++d) out[d] = in[d];
6117: } else {
6118: if (endpoint) {
6119: for (d = 0; d < dim; ++d) {
6120: if ((PetscAbsReal(PetscRealPart(in[d])/dm->L[d] - PetscFloorReal(PetscRealPart(in[d])/dm->L[d])) < PETSC_SMALL) && (PetscRealPart(in[d])/dm->L[d] > PETSC_SMALL)) {
6121: out[d] = in[d] - dm->L[d]*(PetscFloorReal(PetscRealPart(in[d])/dm->L[d]) - 1);
6122: } else {
6123: out[d] = in[d] - dm->L[d]*PetscFloorReal(PetscRealPart(in[d])/dm->L[d]);
6124: }
6125: }
6126: } else {
6127: for (d = 0; d < dim; ++d) {
6128: out[d] = in[d] - dm->L[d]*PetscFloorReal(PetscRealPart(in[d])/dm->L[d]);
6129: }
6130: }
6131: }
6132: return(0);
6133: }
6135: /*
6136: DMLocalizeCoordinate_Internal - If a mesh is periodic, and the input point is far from the anchor, pick the coordinate sheet of the torus which moves it closer.
6138: Input Parameters:
6139: + dm - The DM
6140: . dim - The spatial dimension
6141: . anchor - The anchor point, the input point can be no more than maxCell away from it
6142: - in - The input coordinate point (dim numbers)
6144: Output Parameter:
6145: . out - The localized coordinate point
6147: Level: developer
6149: Note: This is meant to get a set of coordinates close to each other, as in a cell. The anchor is usually the one of the vertices on a containing cell
6151: .seealso: DMLocalizeCoordinates(), DMLocalizeAddCoordinate()
6152: */
6153: PetscErrorCode DMLocalizeCoordinate_Internal(DM dm, PetscInt dim, const PetscScalar anchor[], const PetscScalar in[], PetscScalar out[])
6154: {
6155: PetscInt d;
6158: if (!dm->maxCell) {
6159: for (d = 0; d < dim; ++d) out[d] = in[d];
6160: } else {
6161: for (d = 0; d < dim; ++d) {
6162: if ((dm->bdtype[d] != DM_BOUNDARY_NONE) && (PetscAbsScalar(anchor[d] - in[d]) > dm->maxCell[d])) {
6163: out[d] = PetscRealPart(anchor[d]) > PetscRealPart(in[d]) ? dm->L[d] + in[d] : in[d] - dm->L[d];
6164: } else {
6165: out[d] = in[d];
6166: }
6167: }
6168: }
6169: return(0);
6170: }
6171: PetscErrorCode DMLocalizeCoordinateReal_Internal(DM dm, PetscInt dim, const PetscReal anchor[], const PetscReal in[], PetscReal out[])
6172: {
6173: PetscInt d;
6176: if (!dm->maxCell) {
6177: for (d = 0; d < dim; ++d) out[d] = in[d];
6178: } else {
6179: for (d = 0; d < dim; ++d) {
6180: if ((dm->bdtype[d] != DM_BOUNDARY_NONE) && (PetscAbsReal(anchor[d] - in[d]) > dm->maxCell[d])) {
6181: out[d] = anchor[d] > in[d] ? dm->L[d] + in[d] : in[d] - dm->L[d];
6182: } else {
6183: out[d] = in[d];
6184: }
6185: }
6186: }
6187: return(0);
6188: }
6190: /*
6191: DMLocalizeAddCoordinate_Internal - If a mesh is periodic, and the input point is far from the anchor, pick the coordinate sheet of the torus which moves it closer.
6193: Input Parameters:
6194: + dm - The DM
6195: . dim - The spatial dimension
6196: . anchor - The anchor point, the input point can be no more than maxCell away from it
6197: . in - The input coordinate delta (dim numbers)
6198: - out - The input coordinate point (dim numbers)
6200: Output Parameter:
6201: . out - The localized coordinate in + out
6203: Level: developer
6205: Note: This is meant to get a set of coordinates close to each other, as in a cell. The anchor is usually the one of the vertices on a containing cell
6207: .seealso: DMLocalizeCoordinates(), DMLocalizeCoordinate()
6208: */
6209: PetscErrorCode DMLocalizeAddCoordinate_Internal(DM dm, PetscInt dim, const PetscScalar anchor[], const PetscScalar in[], PetscScalar out[])
6210: {
6211: PetscInt d;
6214: if (!dm->maxCell) {
6215: for (d = 0; d < dim; ++d) out[d] += in[d];
6216: } else {
6217: for (d = 0; d < dim; ++d) {
6218: const PetscReal maxC = dm->maxCell[d];
6220: if ((dm->bdtype[d] != DM_BOUNDARY_NONE) && (PetscAbsScalar(anchor[d] - in[d]) > maxC)) {
6221: const PetscScalar newCoord = PetscRealPart(anchor[d]) > PetscRealPart(in[d]) ? dm->L[d] + in[d] : in[d] - dm->L[d];
6223: if (PetscAbsScalar(newCoord - anchor[d]) > maxC)
6224: SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "%D-Coordinate %g more than %g away from anchor %g", d, (double) PetscRealPart(in[d]), (double) maxC, (double) PetscRealPart(anchor[d]));
6225: out[d] += newCoord;
6226: } else {
6227: out[d] += in[d];
6228: }
6229: }
6230: }
6231: return(0);
6232: }
6234: /*@
6235: DMGetCoordinatesLocalizedLocal - Check if the DM coordinates have been localized for cells on this process
6237: Not collective
6239: Input Parameter:
6240: . dm - The DM
6242: Output Parameter:
6243: areLocalized - True if localized
6245: Level: developer
6247: .seealso: DMLocalizeCoordinates(), DMGetCoordinatesLocalized(), DMSetPeriodicity()
6248: @*/
6249: PetscErrorCode DMGetCoordinatesLocalizedLocal(DM dm,PetscBool *areLocalized)
6250: {
6251: DM cdm;
6252: PetscSection coordSection;
6253: PetscInt cStart, cEnd, sStart, sEnd, c, dof;
6254: PetscBool isPlex, alreadyLocalized;
6260: *areLocalized = PETSC_FALSE;
6262: /* We need some generic way of refering to cells/vertices */
6263: DMGetCoordinateDM(dm, &cdm);
6264: PetscObjectTypeCompare((PetscObject) cdm, DMPLEX, &isPlex);
6265: if (!isPlex) return(0);
6267: DMGetCoordinateSection(dm, &coordSection);
6268: DMPlexGetHeightStratum(cdm, 0, &cStart, &cEnd);
6269: PetscSectionGetChart(coordSection, &sStart, &sEnd);
6270: alreadyLocalized = PETSC_FALSE;
6271: for (c = cStart; c < cEnd; ++c) {
6272: if (c < sStart || c >= sEnd) continue;
6273: PetscSectionGetDof(coordSection, c, &dof);
6274: if (dof) { alreadyLocalized = PETSC_TRUE; break; }
6275: }
6276: *areLocalized = alreadyLocalized;
6277: return(0);
6278: }
6280: /*@
6281: DMGetCoordinatesLocalized - Check if the DM coordinates have been localized for cells
6283: Collective on dm
6285: Input Parameter:
6286: . dm - The DM
6288: Output Parameter:
6289: areLocalized - True if localized
6291: Level: developer
6293: .seealso: DMLocalizeCoordinates(), DMSetPeriodicity(), DMGetCoordinatesLocalizedLocal()
6294: @*/
6295: PetscErrorCode DMGetCoordinatesLocalized(DM dm,PetscBool *areLocalized)
6296: {
6297: PetscBool localized;
6303: DMGetCoordinatesLocalizedLocal(dm,&localized);
6304: MPIU_Allreduce(&localized,areLocalized,1,MPIU_BOOL,MPI_LOR,PetscObjectComm((PetscObject)dm));
6305: return(0);
6306: }
6308: /*@
6309: DMLocalizeCoordinates - If a mesh is periodic, create local coordinates for cells having periodic faces
6311: Collective on dm
6313: Input Parameter:
6314: . dm - The DM
6316: Level: developer
6318: .seealso: DMSetPeriodicity(), DMLocalizeCoordinate(), DMLocalizeAddCoordinate()
6319: @*/
6320: PetscErrorCode DMLocalizeCoordinates(DM dm)
6321: {
6322: DM cdm;
6323: PetscSection coordSection, cSection;
6324: Vec coordinates, cVec;
6325: PetscScalar *coords, *coords2, *anchor, *localized;
6326: PetscInt Nc, vStart, vEnd, v, sStart, sEnd, newStart = PETSC_MAX_INT, newEnd = PETSC_MIN_INT, dof, d, off, off2, bs, coordSize;
6327: PetscBool alreadyLocalized, alreadyLocalizedGlobal;
6328: PetscInt maxHeight = 0, h;
6329: PetscInt *pStart = NULL, *pEnd = NULL;
6334: if (!dm->periodic) return(0);
6335: DMGetCoordinatesLocalized(dm, &alreadyLocalized);
6336: if (alreadyLocalized) return(0);
6338: /* We need some generic way of refering to cells/vertices */
6339: DMGetCoordinateDM(dm, &cdm);
6340: {
6341: PetscBool isplex;
6343: PetscObjectTypeCompare((PetscObject) cdm, DMPLEX, &isplex);
6344: if (isplex) {
6345: DMPlexGetDepthStratum(cdm, 0, &vStart, &vEnd);
6346: DMPlexGetMaxProjectionHeight(cdm,&maxHeight);
6347: DMGetWorkArray(dm,2*(maxHeight + 1),MPIU_INT,&pStart);
6348: pEnd = &pStart[maxHeight + 1];
6349: newStart = vStart;
6350: newEnd = vEnd;
6351: for (h = 0; h <= maxHeight; h++) {
6352: DMPlexGetHeightStratum(cdm, h, &pStart[h], &pEnd[h]);
6353: newStart = PetscMin(newStart,pStart[h]);
6354: newEnd = PetscMax(newEnd,pEnd[h]);
6355: }
6356: } else SETERRQ(PetscObjectComm((PetscObject) cdm), PETSC_ERR_ARG_WRONG, "Coordinate localization requires a DMPLEX coordinate DM");
6357: }
6358: DMGetCoordinatesLocal(dm, &coordinates);
6359: if (!coordinates) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Missing local coordinates vector");
6360: DMGetCoordinateSection(dm, &coordSection);
6361: VecGetBlockSize(coordinates, &bs);
6362: PetscSectionGetChart(coordSection,&sStart,&sEnd);
6364: PetscSectionCreate(PetscObjectComm((PetscObject) dm), &cSection);
6365: PetscSectionSetNumFields(cSection, 1);
6366: PetscSectionGetFieldComponents(coordSection, 0, &Nc);
6367: PetscSectionSetFieldComponents(cSection, 0, Nc);
6368: PetscSectionSetChart(cSection, newStart, newEnd);
6370: DMGetWorkArray(dm, 2 * bs, MPIU_SCALAR, &anchor);
6371: localized = &anchor[bs];
6372: alreadyLocalized = alreadyLocalizedGlobal = PETSC_TRUE;
6373: for (h = 0; h <= maxHeight; h++) {
6374: PetscInt cStart = pStart[h], cEnd = pEnd[h], c;
6376: for (c = cStart; c < cEnd; ++c) {
6377: PetscScalar *cellCoords = NULL;
6378: PetscInt b;
6380: if (c < sStart || c >= sEnd) alreadyLocalized = PETSC_FALSE;
6381: DMPlexVecGetClosure(cdm, coordSection, coordinates, c, &dof, &cellCoords);
6382: for (b = 0; b < bs; ++b) anchor[b] = cellCoords[b];
6383: for (d = 0; d < dof/bs; ++d) {
6384: DMLocalizeCoordinate_Internal(dm, bs, anchor, &cellCoords[d*bs], localized);
6385: for (b = 0; b < bs; b++) {
6386: if (cellCoords[d*bs + b] != localized[b]) break;
6387: }
6388: if (b < bs) break;
6389: }
6390: if (d < dof/bs) {
6391: if (c >= sStart && c < sEnd) {
6392: PetscInt cdof;
6394: PetscSectionGetDof(coordSection, c, &cdof);
6395: if (cdof != dof) alreadyLocalized = PETSC_FALSE;
6396: }
6397: PetscSectionSetDof(cSection, c, dof);
6398: PetscSectionSetFieldDof(cSection, c, 0, dof);
6399: }
6400: DMPlexVecRestoreClosure(cdm, coordSection, coordinates, c, &dof, &cellCoords);
6401: }
6402: }
6403: MPI_Allreduce(&alreadyLocalized,&alreadyLocalizedGlobal,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)dm));
6404: if (alreadyLocalizedGlobal) {
6405: DMRestoreWorkArray(dm, 2 * bs, MPIU_SCALAR, &anchor);
6406: PetscSectionDestroy(&cSection);
6407: DMRestoreWorkArray(dm,2*(maxHeight + 1),MPIU_INT,&pStart);
6408: return(0);
6409: }
6410: for (v = vStart; v < vEnd; ++v) {
6411: PetscSectionGetDof(coordSection, v, &dof);
6412: PetscSectionSetDof(cSection, v, dof);
6413: PetscSectionSetFieldDof(cSection, v, 0, dof);
6414: }
6415: PetscSectionSetUp(cSection);
6416: PetscSectionGetStorageSize(cSection, &coordSize);
6417: VecCreate(PETSC_COMM_SELF, &cVec);
6418: PetscObjectSetName((PetscObject)cVec,"coordinates");
6419: VecSetBlockSize(cVec, bs);
6420: VecSetSizes(cVec, coordSize, PETSC_DETERMINE);
6421: VecSetType(cVec, VECSTANDARD);
6422: VecGetArrayRead(coordinates, (const PetscScalar**)&coords);
6423: VecGetArray(cVec, &coords2);
6424: for (v = vStart; v < vEnd; ++v) {
6425: PetscSectionGetDof(coordSection, v, &dof);
6426: PetscSectionGetOffset(coordSection, v, &off);
6427: PetscSectionGetOffset(cSection, v, &off2);
6428: for (d = 0; d < dof; ++d) coords2[off2+d] = coords[off+d];
6429: }
6430: for (h = 0; h <= maxHeight; h++) {
6431: PetscInt cStart = pStart[h], cEnd = pEnd[h], c;
6433: for (c = cStart; c < cEnd; ++c) {
6434: PetscScalar *cellCoords = NULL;
6435: PetscInt b, cdof;
6437: PetscSectionGetDof(cSection,c,&cdof);
6438: if (!cdof) continue;
6439: DMPlexVecGetClosure(cdm, coordSection, coordinates, c, &dof, &cellCoords);
6440: PetscSectionGetOffset(cSection, c, &off2);
6441: for (b = 0; b < bs; ++b) anchor[b] = cellCoords[b];
6442: for (d = 0; d < dof/bs; ++d) {DMLocalizeCoordinate_Internal(dm, bs, anchor, &cellCoords[d*bs], &coords2[off2+d*bs]);}
6443: DMPlexVecRestoreClosure(cdm, coordSection, coordinates, c, &dof, &cellCoords);
6444: }
6445: }
6446: DMRestoreWorkArray(dm, 2 * bs, MPIU_SCALAR, &anchor);
6447: DMRestoreWorkArray(dm,2*(maxHeight + 1),MPIU_INT,&pStart);
6448: VecRestoreArrayRead(coordinates, (const PetscScalar**)&coords);
6449: VecRestoreArray(cVec, &coords2);
6450: DMSetCoordinateSection(dm, PETSC_DETERMINE, cSection);
6451: DMSetCoordinatesLocal(dm, cVec);
6452: VecDestroy(&cVec);
6453: PetscSectionDestroy(&cSection);
6454: return(0);
6455: }
6457: /*@
6458: DMLocatePoints - Locate the points in v in the mesh and return a PetscSF of the containing cells
6460: Collective on v (see explanation below)
6462: Input Parameters:
6463: + dm - The DM
6464: . v - The Vec of points
6465: . ltype - The type of point location, e.g. DM_POINTLOCATION_NONE or DM_POINTLOCATION_NEAREST
6466: - cells - Points to either NULL, or a PetscSF with guesses for which cells contain each point.
6468: Output Parameter:
6469: + v - The Vec of points, which now contains the nearest mesh points to the given points if DM_POINTLOCATION_NEAREST is used
6470: - cells - The PetscSF containing the ranks and local indices of the containing points.
6473: Level: developer
6475: Notes:
6476: To do a search of the local cells of the mesh, v should have PETSC_COMM_SELF as its communicator.
6477: To do a search of all the cells in the distributed mesh, v should have the same communicator as dm.
6479: If *cellSF is NULL on input, a PetscSF will be created.
6480: If *cellSF is not NULL on input, it should point to an existing PetscSF, whose graph will be used as initial guesses.
6482: An array that maps each point to its containing cell can be obtained with
6484: $ const PetscSFNode *cells;
6485: $ PetscInt nFound;
6486: $ const PetscInt *found;
6487: $
6488: $ PetscSFGetGraph(cellSF,NULL,&nFound,&found,&cells);
6490: Where cells[i].rank is the rank of the cell containing point found[i] (or i if found == NULL), and cells[i].index is
6491: the index of the cell in its rank's local numbering.
6493: .seealso: DMSetCoordinates(), DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLocal(), DMPointLocationType
6494: @*/
6495: PetscErrorCode DMLocatePoints(DM dm, Vec v, DMPointLocationType ltype, PetscSF *cellSF)
6496: {
6503: if (*cellSF) {
6504: PetscMPIInt result;
6507: MPI_Comm_compare(PetscObjectComm((PetscObject)v),PetscObjectComm((PetscObject)*cellSF),&result);
6508: if (result != MPI_IDENT && result != MPI_CONGRUENT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"cellSF must have a communicator congruent to v's");
6509: } else {
6510: PetscSFCreate(PetscObjectComm((PetscObject)v),cellSF);
6511: }
6512: if (!dm->ops->locatepoints) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Point location not available for this DM");
6513: PetscLogEventBegin(DM_LocatePoints,dm,0,0,0);
6514: (*dm->ops->locatepoints)(dm,v,ltype,*cellSF);
6515: PetscLogEventEnd(DM_LocatePoints,dm,0,0,0);
6516: return(0);
6517: }
6519: /*@
6520: DMGetOutputDM - Retrieve the DM associated with the layout for output
6522: Collective on dm
6524: Input Parameter:
6525: . dm - The original DM
6527: Output Parameter:
6528: . odm - The DM which provides the layout for output
6530: Level: intermediate
6532: .seealso: VecView(), DMGetGlobalSection()
6533: @*/
6534: PetscErrorCode DMGetOutputDM(DM dm, DM *odm)
6535: {
6536: PetscSection section;
6537: PetscBool hasConstraints, ghasConstraints;
6543: DMGetLocalSection(dm, §ion);
6544: PetscSectionHasConstraints(section, &hasConstraints);
6545: MPI_Allreduce(&hasConstraints, &ghasConstraints, 1, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject) dm));
6546: if (!ghasConstraints) {
6547: *odm = dm;
6548: return(0);
6549: }
6550: if (!dm->dmBC) {
6551: PetscSection newSection, gsection;
6552: PetscSF sf;
6554: DMClone(dm, &dm->dmBC);
6555: DMCopyDisc(dm, dm->dmBC);
6556: PetscSectionClone(section, &newSection);
6557: DMSetLocalSection(dm->dmBC, newSection);
6558: PetscSectionDestroy(&newSection);
6559: DMGetPointSF(dm->dmBC, &sf);
6560: PetscSectionCreateGlobalSection(section, sf, PETSC_TRUE, PETSC_FALSE, &gsection);
6561: DMSetGlobalSection(dm->dmBC, gsection);
6562: PetscSectionDestroy(&gsection);
6563: }
6564: *odm = dm->dmBC;
6565: return(0);
6566: }
6568: /*@
6569: DMGetOutputSequenceNumber - Retrieve the sequence number/value for output
6571: Input Parameter:
6572: . dm - The original DM
6574: Output Parameters:
6575: + num - The output sequence number
6576: - val - The output sequence value
6578: Level: intermediate
6580: Note: This is intended for output that should appear in sequence, for instance
6581: a set of timesteps in an HDF5 file, or a set of realizations of a stochastic system.
6583: .seealso: VecView()
6584: @*/
6585: PetscErrorCode DMGetOutputSequenceNumber(DM dm, PetscInt *num, PetscReal *val)
6586: {
6591: return(0);
6592: }
6594: /*@
6595: DMSetOutputSequenceNumber - Set the sequence number/value for output
6597: Input Parameters:
6598: + dm - The original DM
6599: . num - The output sequence number
6600: - val - The output sequence value
6602: Level: intermediate
6604: Note: This is intended for output that should appear in sequence, for instance
6605: a set of timesteps in an HDF5 file, or a set of realizations of a stochastic system.
6607: .seealso: VecView()
6608: @*/
6609: PetscErrorCode DMSetOutputSequenceNumber(DM dm, PetscInt num, PetscReal val)
6610: {
6613: dm->outputSequenceNum = num;
6614: dm->outputSequenceVal = val;
6615: return(0);
6616: }
6618: /*@C
6619: DMOutputSequenceLoad - Retrieve the sequence value from a Viewer
6621: Input Parameters:
6622: + dm - The original DM
6623: . name - The sequence name
6624: - num - The output sequence number
6626: Output Parameter:
6627: . val - The output sequence value
6629: Level: intermediate
6631: Note: This is intended for output that should appear in sequence, for instance
6632: a set of timesteps in an HDF5 file, or a set of realizations of a stochastic system.
6634: .seealso: DMGetOutputSequenceNumber(), DMSetOutputSequenceNumber(), VecView()
6635: @*/
6636: PetscErrorCode DMOutputSequenceLoad(DM dm, PetscViewer viewer, const char *name, PetscInt num, PetscReal *val)
6637: {
6638: PetscBool ishdf5;
6645: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERHDF5, &ishdf5);
6646: if (ishdf5) {
6647: #if defined(PETSC_HAVE_HDF5)
6648: PetscScalar value;
6650: DMSequenceLoad_HDF5_Internal(dm, name, num, &value, viewer);
6651: *val = PetscRealPart(value);
6652: #endif
6653: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid viewer; open viewer with PetscViewerHDF5Open()");
6654: return(0);
6655: }
6657: /*@
6658: DMGetUseNatural - Get the flag for creating a mapping to the natural order on distribution
6660: Not collective
6662: Input Parameter:
6663: . dm - The DM
6665: Output Parameter:
6666: . useNatural - The flag to build the mapping to a natural order during distribution
6668: Level: beginner
6670: .seealso: DMSetUseNatural(), DMCreate()
6671: @*/
6672: PetscErrorCode DMGetUseNatural(DM dm, PetscBool *useNatural)
6673: {
6677: *useNatural = dm->useNatural;
6678: return(0);
6679: }
6681: /*@
6682: DMSetUseNatural - Set the flag for creating a mapping to the natural order after distribution
6684: Collective on dm
6686: Input Parameters:
6687: + dm - The DM
6688: - useNatural - The flag to build the mapping to a natural order during distribution
6690: Note: This also causes the map to be build after DMCreateSubDM() and DMCreateSuperDM()
6692: Level: beginner
6694: .seealso: DMGetUseNatural(), DMCreate(), DMPlexDistribute(), DMCreateSubDM(), DMCreateSuperDM()
6695: @*/
6696: PetscErrorCode DMSetUseNatural(DM dm, PetscBool useNatural)
6697: {
6701: dm->useNatural = useNatural;
6702: return(0);
6703: }
6706: /*@C
6707: DMCreateLabel - Create a label of the given name if it does not already exist
6709: Not Collective
6711: Input Parameters:
6712: + dm - The DM object
6713: - name - The label name
6715: Level: intermediate
6717: .seealso: DMLabelCreate(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
6718: @*/
6719: PetscErrorCode DMCreateLabel(DM dm, const char name[])
6720: {
6721: PetscBool flg;
6722: DMLabel label;
6728: DMHasLabel(dm, name, &flg);
6729: if (!flg) {
6730: DMLabelCreate(PETSC_COMM_SELF, name, &label);
6731: DMAddLabel(dm, label);
6732: DMLabelDestroy(&label);
6733: }
6734: return(0);
6735: }
6737: /*@C
6738: DMGetLabelValue - Get the value in a Sieve Label for the given point, with 0 as the default
6740: Not Collective
6742: Input Parameters:
6743: + dm - The DM object
6744: . name - The label name
6745: - point - The mesh point
6747: Output Parameter:
6748: . value - The label value for this point, or -1 if the point is not in the label
6750: Level: beginner
6752: .seealso: DMLabelGetValue(), DMSetLabelValue(), DMGetStratumIS()
6753: @*/
6754: PetscErrorCode DMGetLabelValue(DM dm, const char name[], PetscInt point, PetscInt *value)
6755: {
6756: DMLabel label;
6762: DMGetLabel(dm, name, &label);
6763: if (!label) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No label named %s was found", name);
6764: DMLabelGetValue(label, point, value);
6765: return(0);
6766: }
6768: /*@C
6769: DMSetLabelValue - Add a point to a Sieve Label with given value
6771: Not Collective
6773: Input Parameters:
6774: + dm - The DM object
6775: . name - The label name
6776: . point - The mesh point
6777: - value - The label value for this point
6779: Output Parameter:
6781: Level: beginner
6783: .seealso: DMLabelSetValue(), DMGetStratumIS(), DMClearLabelValue()
6784: @*/
6785: PetscErrorCode DMSetLabelValue(DM dm, const char name[], PetscInt point, PetscInt value)
6786: {
6787: DMLabel label;
6793: DMGetLabel(dm, name, &label);
6794: if (!label) {
6795: DMCreateLabel(dm, name);
6796: DMGetLabel(dm, name, &label);
6797: }
6798: DMLabelSetValue(label, point, value);
6799: return(0);
6800: }
6802: /*@C
6803: DMClearLabelValue - Remove a point from a Sieve Label with given value
6805: Not Collective
6807: Input Parameters:
6808: + dm - The DM object
6809: . name - The label name
6810: . point - The mesh point
6811: - value - The label value for this point
6813: Output Parameter:
6815: Level: beginner
6817: .seealso: DMLabelClearValue(), DMSetLabelValue(), DMGetStratumIS()
6818: @*/
6819: PetscErrorCode DMClearLabelValue(DM dm, const char name[], PetscInt point, PetscInt value)
6820: {
6821: DMLabel label;
6827: DMGetLabel(dm, name, &label);
6828: if (!label) return(0);
6829: DMLabelClearValue(label, point, value);
6830: return(0);
6831: }
6833: /*@C
6834: DMGetLabelSize - Get the number of different integer ids in a Label
6836: Not Collective
6838: Input Parameters:
6839: + dm - The DM object
6840: - name - The label name
6842: Output Parameter:
6843: . size - The number of different integer ids, or 0 if the label does not exist
6845: Level: beginner
6847: .seealso: DMLabelGetNumValues(), DMSetLabelValue()
6848: @*/
6849: PetscErrorCode DMGetLabelSize(DM dm, const char name[], PetscInt *size)
6850: {
6851: DMLabel label;
6858: DMGetLabel(dm, name, &label);
6859: *size = 0;
6860: if (!label) return(0);
6861: DMLabelGetNumValues(label, size);
6862: return(0);
6863: }
6865: /*@C
6866: DMGetLabelIdIS - Get the integer ids in a label
6868: Not Collective
6870: Input Parameters:
6871: + mesh - The DM object
6872: - name - The label name
6874: Output Parameter:
6875: . ids - The integer ids, or NULL if the label does not exist
6877: Level: beginner
6879: .seealso: DMLabelGetValueIS(), DMGetLabelSize()
6880: @*/
6881: PetscErrorCode DMGetLabelIdIS(DM dm, const char name[], IS *ids)
6882: {
6883: DMLabel label;
6890: DMGetLabel(dm, name, &label);
6891: *ids = NULL;
6892: if (label) {
6893: DMLabelGetValueIS(label, ids);
6894: } else {
6895: /* returning an empty IS */
6896: ISCreateGeneral(PETSC_COMM_SELF,0,NULL,PETSC_USE_POINTER,ids);
6897: }
6898: return(0);
6899: }
6901: /*@C
6902: DMGetStratumSize - Get the number of points in a label stratum
6904: Not Collective
6906: Input Parameters:
6907: + dm - The DM object
6908: . name - The label name
6909: - value - The stratum value
6911: Output Parameter:
6912: . size - The stratum size
6914: Level: beginner
6916: .seealso: DMLabelGetStratumSize(), DMGetLabelSize(), DMGetLabelIds()
6917: @*/
6918: PetscErrorCode DMGetStratumSize(DM dm, const char name[], PetscInt value, PetscInt *size)
6919: {
6920: DMLabel label;
6927: DMGetLabel(dm, name, &label);
6928: *size = 0;
6929: if (!label) return(0);
6930: DMLabelGetStratumSize(label, value, size);
6931: return(0);
6932: }
6934: /*@C
6935: DMGetStratumIS - Get the points in a label stratum
6937: Not Collective
6939: Input Parameters:
6940: + dm - The DM object
6941: . name - The label name
6942: - value - The stratum value
6944: Output Parameter:
6945: . points - The stratum points, or NULL if the label does not exist or does not have that value
6947: Level: beginner
6949: .seealso: DMLabelGetStratumIS(), DMGetStratumSize()
6950: @*/
6951: PetscErrorCode DMGetStratumIS(DM dm, const char name[], PetscInt value, IS *points)
6952: {
6953: DMLabel label;
6960: DMGetLabel(dm, name, &label);
6961: *points = NULL;
6962: if (!label) return(0);
6963: DMLabelGetStratumIS(label, value, points);
6964: return(0);
6965: }
6967: /*@C
6968: DMSetStratumIS - Set the points in a label stratum
6970: Not Collective
6972: Input Parameters:
6973: + dm - The DM object
6974: . name - The label name
6975: . value - The stratum value
6976: - points - The stratum points
6978: Level: beginner
6980: .seealso: DMLabelSetStratumIS(), DMGetStratumSize()
6981: @*/
6982: PetscErrorCode DMSetStratumIS(DM dm, const char name[], PetscInt value, IS points)
6983: {
6984: DMLabel label;
6991: DMGetLabel(dm, name, &label);
6992: if (!label) return(0);
6993: DMLabelSetStratumIS(label, value, points);
6994: return(0);
6995: }
6997: /*@C
6998: DMClearLabelStratum - Remove all points from a stratum from a Sieve Label
7000: Not Collective
7002: Input Parameters:
7003: + dm - The DM object
7004: . name - The label name
7005: - value - The label value for this point
7007: Output Parameter:
7009: Level: beginner
7011: .seealso: DMLabelClearStratum(), DMSetLabelValue(), DMGetStratumIS(), DMClearLabelValue()
7012: @*/
7013: PetscErrorCode DMClearLabelStratum(DM dm, const char name[], PetscInt value)
7014: {
7015: DMLabel label;
7021: DMGetLabel(dm, name, &label);
7022: if (!label) return(0);
7023: DMLabelClearStratum(label, value);
7024: return(0);
7025: }
7027: /*@
7028: DMGetNumLabels - Return the number of labels defined by the mesh
7030: Not Collective
7032: Input Parameter:
7033: . dm - The DM object
7035: Output Parameter:
7036: . numLabels - the number of Labels
7038: Level: intermediate
7040: .seealso: DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
7041: @*/
7042: PetscErrorCode DMGetNumLabels(DM dm, PetscInt *numLabels)
7043: {
7044: DMLabelLink next = dm->labels;
7045: PetscInt n = 0;
7050: while (next) {++n; next = next->next;}
7051: *numLabels = n;
7052: return(0);
7053: }
7055: /*@C
7056: DMGetLabelName - Return the name of nth label
7058: Not Collective
7060: Input Parameters:
7061: + dm - The DM object
7062: - n - the label number
7064: Output Parameter:
7065: . name - the label name
7067: Level: intermediate
7069: .seealso: DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
7070: @*/
7071: PetscErrorCode DMGetLabelName(DM dm, PetscInt n, const char **name)
7072: {
7073: DMLabelLink next = dm->labels;
7074: PetscInt l = 0;
7080: while (next) {
7081: if (l == n) {
7082: PetscObjectGetName((PetscObject) next->label, name);
7083: return(0);
7084: }
7085: ++l;
7086: next = next->next;
7087: }
7088: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label %D does not exist in this DM", n);
7089: }
7091: /*@C
7092: DMHasLabel - Determine whether the mesh has a label of a given name
7094: Not Collective
7096: Input Parameters:
7097: + dm - The DM object
7098: - name - The label name
7100: Output Parameter:
7101: . hasLabel - PETSC_TRUE if the label is present
7103: Level: intermediate
7105: .seealso: DMCreateLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
7106: @*/
7107: PetscErrorCode DMHasLabel(DM dm, const char name[], PetscBool *hasLabel)
7108: {
7109: DMLabelLink next = dm->labels;
7110: const char *lname;
7117: *hasLabel = PETSC_FALSE;
7118: while (next) {
7119: PetscObjectGetName((PetscObject) next->label, &lname);
7120: PetscStrcmp(name, lname, hasLabel);
7121: if (*hasLabel) break;
7122: next = next->next;
7123: }
7124: return(0);
7125: }
7127: /*@C
7128: DMGetLabel - Return the label of a given name, or NULL
7130: Not Collective
7132: Input Parameters:
7133: + dm - The DM object
7134: - name - The label name
7136: Output Parameter:
7137: . label - The DMLabel, or NULL if the label is absent
7139: Level: intermediate
7141: .seealso: DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
7142: @*/
7143: PetscErrorCode DMGetLabel(DM dm, const char name[], DMLabel *label)
7144: {
7145: DMLabelLink next = dm->labels;
7146: PetscBool hasLabel;
7147: const char *lname;
7154: *label = NULL;
7155: while (next) {
7156: PetscObjectGetName((PetscObject) next->label, &lname);
7157: PetscStrcmp(name, lname, &hasLabel);
7158: if (hasLabel) {
7159: *label = next->label;
7160: break;
7161: }
7162: next = next->next;
7163: }
7164: return(0);
7165: }
7167: /*@C
7168: DMGetLabelByNum - Return the nth label
7170: Not Collective
7172: Input Parameters:
7173: + dm - The DM object
7174: - n - the label number
7176: Output Parameter:
7177: . label - the label
7179: Level: intermediate
7181: .seealso: DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
7182: @*/
7183: PetscErrorCode DMGetLabelByNum(DM dm, PetscInt n, DMLabel *label)
7184: {
7185: DMLabelLink next = dm->labels;
7186: PetscInt l = 0;
7191: while (next) {
7192: if (l == n) {
7193: *label = next->label;
7194: return(0);
7195: }
7196: ++l;
7197: next = next->next;
7198: }
7199: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label %D does not exist in this DM", n);
7200: }
7202: /*@C
7203: DMAddLabel - Add the label to this mesh
7205: Not Collective
7207: Input Parameters:
7208: + dm - The DM object
7209: - label - The DMLabel
7211: Level: developer
7213: .seealso: DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
7214: @*/
7215: PetscErrorCode DMAddLabel(DM dm, DMLabel label)
7216: {
7217: DMLabelLink l, *p, tmpLabel;
7218: PetscBool hasLabel;
7219: const char *lname;
7220: PetscBool flg;
7225: PetscObjectGetName((PetscObject) label, &lname);
7226: DMHasLabel(dm, lname, &hasLabel);
7227: if (hasLabel) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label %s already exists in this DM", lname);
7228: PetscCalloc1(1, &tmpLabel);
7229: tmpLabel->label = label;
7230: tmpLabel->output = PETSC_TRUE;
7231: for (p=&dm->labels; (l=*p); p=&l->next) {}
7232: *p = tmpLabel;
7233: PetscObjectReference((PetscObject)label);
7234: PetscStrcmp(lname, "depth", &flg);
7235: if (flg) dm->depthLabel = label;
7236: PetscStrcmp(lname, "celltype", &flg);
7237: if (flg) dm->celltypeLabel = label;
7238: return(0);
7239: }
7241: /*@C
7242: DMRemoveLabel - Remove the label given by name from this mesh
7244: Not Collective
7246: Input Parameters:
7247: + dm - The DM object
7248: - name - The label name
7250: Output Parameter:
7251: . label - The DMLabel, or NULL if the label is absent
7253: Level: developer
7255: Notes:
7256: DMRemoveLabel(dm,name,NULL) removes the label from dm and calls
7257: DMLabelDestroy() on the label.
7259: DMRemoveLabel(dm,name,&label) removes the label from dm, but it DOES NOT
7260: call DMLabelDestroy(). Instead, the label is returned and the user is
7261: responsible of calling DMLabelDestroy() at some point.
7263: .seealso: DMCreateLabel(), DMHasLabel(), DMGetLabel(), DMGetLabelValue(), DMSetLabelValue(), DMLabelDestroy(), DMRemoveLabelBySelf()
7264: @*/
7265: PetscErrorCode DMRemoveLabel(DM dm, const char name[], DMLabel *label)
7266: {
7267: DMLabelLink link, *pnext;
7268: PetscBool hasLabel;
7269: const char *lname;
7275: if (label) {
7277: *label = NULL;
7278: }
7279: for (pnext=&dm->labels; (link=*pnext); pnext=&link->next) {
7280: PetscObjectGetName((PetscObject) link->label, &lname);
7281: PetscStrcmp(name, lname, &hasLabel);
7282: if (hasLabel) {
7283: *pnext = link->next; /* Remove from list */
7284: PetscStrcmp(name, "depth", &hasLabel);
7285: if (hasLabel) dm->depthLabel = NULL;
7286: PetscStrcmp(name, "celltype", &hasLabel);
7287: if (hasLabel) dm->celltypeLabel = NULL;
7288: if (label) *label = link->label;
7289: else {DMLabelDestroy(&link->label);}
7290: PetscFree(link);
7291: break;
7292: }
7293: }
7294: return(0);
7295: }
7297: /*@
7298: DMRemoveLabelBySelf - Remove the label from this mesh
7300: Not Collective
7302: Input Parameters:
7303: + dm - The DM object
7304: . label - (Optional) The DMLabel to be removed from the DM
7305: - failNotFound - Should it fail if the label is not found in the DM?
7307: Level: developer
7309: Notes:
7310: Only exactly the same instance is removed if found, name match is ignored.
7311: If the DM has an exclusive reference to the label, it gets destroyed and
7312: *label nullified.
7314: .seealso: DMCreateLabel(), DMHasLabel(), DMGetLabel() DMGetLabelValue(), DMSetLabelValue(), DMLabelDestroy(), DMRemoveLabel()
7315: @*/
7316: PetscErrorCode DMRemoveLabelBySelf(DM dm, DMLabel *label, PetscBool failNotFound)
7317: {
7318: DMLabelLink link, *pnext;
7319: PetscBool hasLabel = PETSC_FALSE;
7325: if (!*label && !failNotFound) return(0);
7328: for (pnext=&dm->labels; (link=*pnext); pnext=&link->next) {
7329: if (*label == link->label) {
7330: hasLabel = PETSC_TRUE;
7331: *pnext = link->next; /* Remove from list */
7332: if (*label == dm->depthLabel) dm->depthLabel = NULL;
7333: if (*label == dm->celltypeLabel) dm->celltypeLabel = NULL;
7334: if (((PetscObject) link->label)->refct < 2) *label = NULL; /* nullify if exclusive reference */
7335: DMLabelDestroy(&link->label);
7336: PetscFree(link);
7337: break;
7338: }
7339: }
7340: if (!hasLabel && failNotFound) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Given label not found in DM");
7341: return(0);
7342: }
7344: /*@C
7345: DMGetLabelOutput - Get the output flag for a given label
7347: Not Collective
7349: Input Parameters:
7350: + dm - The DM object
7351: - name - The label name
7353: Output Parameter:
7354: . output - The flag for output
7356: Level: developer
7358: .seealso: DMSetLabelOutput(), DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
7359: @*/
7360: PetscErrorCode DMGetLabelOutput(DM dm, const char name[], PetscBool *output)
7361: {
7362: DMLabelLink next = dm->labels;
7363: const char *lname;
7370: while (next) {
7371: PetscBool flg;
7373: PetscObjectGetName((PetscObject) next->label, &lname);
7374: PetscStrcmp(name, lname, &flg);
7375: if (flg) {*output = next->output; return(0);}
7376: next = next->next;
7377: }
7378: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No label named %s was present in this dm", name);
7379: }
7381: /*@C
7382: DMSetLabelOutput - Set the output flag for a given label
7384: Not Collective
7386: Input Parameters:
7387: + dm - The DM object
7388: . name - The label name
7389: - output - The flag for output
7391: Level: developer
7393: .seealso: DMGetLabelOutput(), DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
7394: @*/
7395: PetscErrorCode DMSetLabelOutput(DM dm, const char name[], PetscBool output)
7396: {
7397: DMLabelLink next = dm->labels;
7398: const char *lname;
7404: while (next) {
7405: PetscBool flg;
7407: PetscObjectGetName((PetscObject) next->label, &lname);
7408: PetscStrcmp(name, lname, &flg);
7409: if (flg) {next->output = output; return(0);}
7410: next = next->next;
7411: }
7412: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No label named %s was present in this dm", name);
7413: }
7415: /*@
7416: DMCopyLabels - Copy labels from one mesh to another with a superset of the points
7418: Collective on dmA
7420: Input Parameter:
7421: + dmA - The DM object with initial labels
7422: . dmB - The DM object with copied labels
7423: . mode - Copy labels by pointers (PETSC_OWN_POINTER) or duplicate them (PETSC_COPY_VALUES)
7424: - all - Copy all labels including "depth", "dim", and "celltype" (PETSC_TRUE) which are otherwise ignored (PETSC_FALSE)
7426: Level: intermediate
7428: Note: This is typically used when interpolating or otherwise adding to a mesh
7430: .seealso: DMGetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM(), DMGetCoordinateSection(), DMShareLabels()
7431: @*/
7432: PetscErrorCode DMCopyLabels(DM dmA, DM dmB, PetscCopyMode mode, PetscBool all)
7433: {
7434: DMLabel label, labelNew;
7435: const char *name;
7436: PetscBool flg;
7437: DMLabelLink link;
7445: if (mode==PETSC_USE_POINTER) SETERRQ(PetscObjectComm((PetscObject)dmA), PETSC_ERR_SUP, "PETSC_USE_POINTER not supported for objects");
7446: if (dmA == dmB) return(0);
7447: for (link=dmA->labels; link; link=link->next) {
7448: label=link->label;
7449: PetscObjectGetName((PetscObject)label, &name);
7450: if (!all) {
7451: PetscStrcmp(name, "depth", &flg);
7452: if (flg) continue;
7453: PetscStrcmp(name, "dim", &flg);
7454: if (flg) continue;
7455: PetscStrcmp(name, "celltype", &flg);
7456: if (flg) continue;
7457: }
7458: if (mode==PETSC_COPY_VALUES) {
7459: DMLabelDuplicate(label, &labelNew);
7460: } else {
7461: labelNew = label;
7462: }
7463: DMAddLabel(dmB, labelNew);
7464: if (mode==PETSC_COPY_VALUES) {DMLabelDestroy(&labelNew);}
7465: }
7466: return(0);
7467: }
7469: /*@
7470: DMGetCoarseDM - Get the coarse mesh from which this was obtained by refinement
7472: Input Parameter:
7473: . dm - The DM object
7475: Output Parameter:
7476: . cdm - The coarse DM
7478: Level: intermediate
7480: .seealso: DMSetCoarseDM()
7481: @*/
7482: PetscErrorCode DMGetCoarseDM(DM dm, DM *cdm)
7483: {
7487: *cdm = dm->coarseMesh;
7488: return(0);
7489: }
7491: /*@
7492: DMSetCoarseDM - Set the coarse mesh from which this was obtained by refinement
7494: Input Parameters:
7495: + dm - The DM object
7496: - cdm - The coarse DM
7498: Level: intermediate
7500: .seealso: DMGetCoarseDM()
7501: @*/
7502: PetscErrorCode DMSetCoarseDM(DM dm, DM cdm)
7503: {
7509: PetscObjectReference((PetscObject)cdm);
7510: DMDestroy(&dm->coarseMesh);
7511: dm->coarseMesh = cdm;
7512: return(0);
7513: }
7515: /*@
7516: DMGetFineDM - Get the fine mesh from which this was obtained by refinement
7518: Input Parameter:
7519: . dm - The DM object
7521: Output Parameter:
7522: . fdm - The fine DM
7524: Level: intermediate
7526: .seealso: DMSetFineDM()
7527: @*/
7528: PetscErrorCode DMGetFineDM(DM dm, DM *fdm)
7529: {
7533: *fdm = dm->fineMesh;
7534: return(0);
7535: }
7537: /*@
7538: DMSetFineDM - Set the fine mesh from which this was obtained by refinement
7540: Input Parameters:
7541: + dm - The DM object
7542: - fdm - The fine DM
7544: Level: intermediate
7546: .seealso: DMGetFineDM()
7547: @*/
7548: PetscErrorCode DMSetFineDM(DM dm, DM fdm)
7549: {
7555: PetscObjectReference((PetscObject)fdm);
7556: DMDestroy(&dm->fineMesh);
7557: dm->fineMesh = fdm;
7558: return(0);
7559: }
7561: /*=== DMBoundary code ===*/
7563: PetscErrorCode DMCopyBoundary(DM dm, DM dmNew)
7564: {
7565: PetscInt d;
7569: for (d = 0; d < dm->Nds; ++d) {
7570: PetscDSCopyBoundary(dm->probs[d].ds, dmNew->probs[d].ds);
7571: }
7572: return(0);
7573: }
7575: /*@C
7576: DMAddBoundary - Add a boundary condition to the model
7578: Input Parameters:
7579: + dm - The DM, with a PetscDS that matches the problem being constrained
7580: . type - The type of condition, e.g. DM_BC_ESSENTIAL_ANALYTIC/DM_BC_ESSENTIAL_FIELD (Dirichlet), or DM_BC_NATURAL (Neumann)
7581: . name - The BC name
7582: . labelname - The label defining constrained points
7583: . field - The field to constrain
7584: . numcomps - The number of constrained field components (0 will constrain all fields)
7585: . comps - An array of constrained component numbers
7586: . bcFunc - A pointwise function giving boundary values
7587: . numids - The number of DMLabel ids for constrained points
7588: . ids - An array of ids for constrained points
7589: - ctx - An optional user context for bcFunc
7591: Options Database Keys:
7592: + -bc_<boundary name> <num> - Overrides the boundary ids
7593: - -bc_<boundary name>_comp <num> - Overrides the boundary components
7595: Level: developer
7597: .seealso: DMGetBoundary()
7598: @*/
7599: PetscErrorCode DMAddBoundary(DM dm, DMBoundaryConditionType type, const char name[], const char labelname[], PetscInt field, PetscInt numcomps, const PetscInt *comps, void (*bcFunc)(void), PetscInt numids, const PetscInt *ids, void *ctx)
7600: {
7601: PetscDS ds;
7606: DMGetDS(dm, &ds);
7607: PetscDSAddBoundary(ds, type,name, labelname, field, numcomps, comps, bcFunc, numids, ids, ctx);
7608: return(0);
7609: }
7611: /*@
7612: DMGetNumBoundary - Get the number of registered BC
7614: Input Parameters:
7615: . dm - The mesh object
7617: Output Parameters:
7618: . numBd - The number of BC
7620: Level: intermediate
7622: .seealso: DMAddBoundary(), DMGetBoundary()
7623: @*/
7624: PetscErrorCode DMGetNumBoundary(DM dm, PetscInt *numBd)
7625: {
7626: PetscDS ds;
7631: DMGetDS(dm, &ds);
7632: PetscDSGetNumBoundary(ds, numBd);
7633: return(0);
7634: }
7636: /*@C
7637: DMGetBoundary - Get a model boundary condition
7639: Input Parameters:
7640: + dm - The mesh object
7641: - bd - The BC number
7643: Output Parameters:
7644: + type - The type of condition, e.g. DM_BC_ESSENTIAL_ANALYTIC/DM_BC_ESSENTIAL_FIELD (Dirichlet), or DM_BC_NATURAL (Neumann)
7645: . name - The BC name
7646: . labelname - The label defining constrained points
7647: . field - The field to constrain
7648: . numcomps - The number of constrained field components
7649: . comps - An array of constrained component numbers
7650: . bcFunc - A pointwise function giving boundary values
7651: . numids - The number of DMLabel ids for constrained points
7652: . ids - An array of ids for constrained points
7653: - ctx - An optional user context for bcFunc
7655: Options Database Keys:
7656: + -bc_<boundary name> <num> - Overrides the boundary ids
7657: - -bc_<boundary name>_comp <num> - Overrides the boundary components
7659: Level: developer
7661: .seealso: DMAddBoundary()
7662: @*/
7663: PetscErrorCode DMGetBoundary(DM dm, PetscInt bd, DMBoundaryConditionType *type, const char **name, const char **labelname, PetscInt *field, PetscInt *numcomps, const PetscInt **comps, void (**func)(void), PetscInt *numids, const PetscInt **ids, void **ctx)
7664: {
7665: PetscDS ds;
7670: DMGetDS(dm, &ds);
7671: PetscDSGetBoundary(ds, bd, type, name, labelname, field, numcomps, comps, func, numids, ids, ctx);
7672: return(0);
7673: }
7675: static PetscErrorCode DMPopulateBoundary(DM dm)
7676: {
7677: PetscDS ds;
7678: DMBoundary *lastnext;
7679: DSBoundary dsbound;
7683: DMGetDS(dm, &ds);
7684: dsbound = ds->boundary;
7685: if (dm->boundary) {
7686: DMBoundary next = dm->boundary;
7688: /* quick check to see if the PetscDS has changed */
7689: if (next->dsboundary == dsbound) return(0);
7690: /* the PetscDS has changed: tear down and rebuild */
7691: while (next) {
7692: DMBoundary b = next;
7694: next = b->next;
7695: PetscFree(b);
7696: }
7697: dm->boundary = NULL;
7698: }
7700: lastnext = &(dm->boundary);
7701: while (dsbound) {
7702: DMBoundary dmbound;
7704: PetscNew(&dmbound);
7705: dmbound->dsboundary = dsbound;
7706: DMGetLabel(dm, dsbound->labelname, &(dmbound->label));
7707: if (!dmbound->label) {PetscInfo2(dm, "DSBoundary %s wants label %s, which is not in this dm.\n",dsbound->name,dsbound->labelname);}
7708: /* push on the back instead of the front so that it is in the same order as in the PetscDS */
7709: *lastnext = dmbound;
7710: lastnext = &(dmbound->next);
7711: dsbound = dsbound->next;
7712: }
7713: return(0);
7714: }
7716: PetscErrorCode DMIsBoundaryPoint(DM dm, PetscInt point, PetscBool *isBd)
7717: {
7718: DMBoundary b;
7724: *isBd = PETSC_FALSE;
7725: DMPopulateBoundary(dm);
7726: b = dm->boundary;
7727: while (b && !(*isBd)) {
7728: DMLabel label = b->label;
7729: DSBoundary dsb = b->dsboundary;
7731: if (label) {
7732: PetscInt i;
7734: for (i = 0; i < dsb->numids && !(*isBd); ++i) {
7735: DMLabelStratumHasPoint(label, dsb->ids[i], point, isBd);
7736: }
7737: }
7738: b = b->next;
7739: }
7740: return(0);
7741: }
7743: /*@C
7744: DMProjectFunction - This projects the given function into the function space provided, putting the coefficients in a global vector.
7746: Collective on DM
7748: Input Parameters:
7749: + dm - The DM
7750: . time - The time
7751: . funcs - The coordinate functions to evaluate, one per field
7752: . ctxs - Optional array of contexts to pass to each coordinate function. ctxs itself may be null.
7753: - mode - The insertion mode for values
7755: Output Parameter:
7756: . X - vector
7758: Calling sequence of func:
7759: $ func(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar u[], void *ctx);
7761: + dim - The spatial dimension
7762: . x - The coordinates
7763: . Nf - The number of fields
7764: . u - The output field values
7765: - ctx - optional user-defined function context
7767: Level: developer
7769: .seealso: DMProjectFunctionLocal(), DMProjectFunctionLabel(), DMComputeL2Diff()
7770: @*/
7771: PetscErrorCode DMProjectFunction(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec X)
7772: {
7773: Vec localX;
7778: DMGetLocalVector(dm, &localX);
7779: DMProjectFunctionLocal(dm, time, funcs, ctxs, mode, localX);
7780: DMLocalToGlobalBegin(dm, localX, mode, X);
7781: DMLocalToGlobalEnd(dm, localX, mode, X);
7782: DMRestoreLocalVector(dm, &localX);
7783: return(0);
7784: }
7786: /*@C
7787: DMProjectFunctionLocal - This projects the given function into the function space provided, putting the coefficients in a local vector.
7789: Not collective
7791: Input Parameters:
7792: + dm - The DM
7793: . time - The time
7794: . funcs - The coordinate functions to evaluate, one per field
7795: . ctxs - Optional array of contexts to pass to each coordinate function. ctxs itself may be null.
7796: - mode - The insertion mode for values
7798: Output Parameter:
7799: . localX - vector
7801: Calling sequence of func:
7802: $ func(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar u[], void *ctx);
7804: + dim - The spatial dimension
7805: . x - The coordinates
7806: . Nf - The number of fields
7807: . u - The output field values
7808: - ctx - optional user-defined function context
7810: Level: developer
7812: .seealso: DMProjectFunction(), DMProjectFunctionLabel(), DMComputeL2Diff()
7813: @*/
7814: PetscErrorCode DMProjectFunctionLocal(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX)
7815: {
7821: if (!dm->ops->projectfunctionlocal) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMProjectFunctionLocal",((PetscObject)dm)->type_name);
7822: (dm->ops->projectfunctionlocal) (dm, time, funcs, ctxs, mode, localX);
7823: return(0);
7824: }
7826: /*@C
7827: DMProjectFunctionLabel - This projects the given function into the function space provided, putting the coefficients in a global vector, setting values only for points in the given label.
7829: Collective on DM
7831: Input Parameters:
7832: + dm - The DM
7833: . time - The time
7834: . label - The DMLabel selecting the portion of the mesh for projection
7835: . funcs - The coordinate functions to evaluate, one per field
7836: . ctxs - Optional array of contexts to pass to each coordinate function. ctxs itself may be null.
7837: - mode - The insertion mode for values
7839: Output Parameter:
7840: . X - vector
7842: Calling sequence of func:
7843: $ func(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar u[], void *ctx);
7845: + dim - The spatial dimension
7846: . x - The coordinates
7847: . Nf - The number of fields
7848: . u - The output field values
7849: - ctx - optional user-defined function context
7851: Level: developer
7853: .seealso: DMProjectFunction(), DMProjectFunctionLocal(), DMProjectFunctionLabelLocal(), DMComputeL2Diff()
7854: @*/
7855: PetscErrorCode DMProjectFunctionLabel(DM dm, PetscReal time, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscInt Nc, const PetscInt comps[], PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec X)
7856: {
7857: Vec localX;
7862: DMGetLocalVector(dm, &localX);
7863: DMProjectFunctionLabelLocal(dm, time, label, numIds, ids, Nc, comps, funcs, ctxs, mode, localX);
7864: DMLocalToGlobalBegin(dm, localX, mode, X);
7865: DMLocalToGlobalEnd(dm, localX, mode, X);
7866: DMRestoreLocalVector(dm, &localX);
7867: return(0);
7868: }
7870: /*@C
7871: DMProjectFunctionLabelLocal - This projects the given function into the function space provided, putting the coefficients in a local vector, setting values only for points in the given label.
7873: Not collective
7875: Input Parameters:
7876: + dm - The DM
7877: . time - The time
7878: . label - The DMLabel selecting the portion of the mesh for projection
7879: . funcs - The coordinate functions to evaluate, one per field
7880: . ctxs - Optional array of contexts to pass to each coordinate function. ctxs itself may be null.
7881: - mode - The insertion mode for values
7883: Output Parameter:
7884: . localX - vector
7886: Calling sequence of func:
7887: $ func(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar u[], void *ctx);
7889: + dim - The spatial dimension
7890: . x - The coordinates
7891: . Nf - The number of fields
7892: . u - The output field values
7893: - ctx - optional user-defined function context
7895: Level: developer
7897: .seealso: DMProjectFunction(), DMProjectFunctionLocal(), DMProjectFunctionLabel(), DMComputeL2Diff()
7898: @*/
7899: PetscErrorCode DMProjectFunctionLabelLocal(DM dm, PetscReal time, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscInt Nc, const PetscInt comps[], PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX)
7900: {
7906: if (!dm->ops->projectfunctionlabellocal) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMProjectFunctionLabelLocal",((PetscObject)dm)->type_name);
7907: (dm->ops->projectfunctionlabellocal) (dm, time, label, numIds, ids, Nc, comps, funcs, ctxs, mode, localX);
7908: return(0);
7909: }
7911: /*@C
7912: DMProjectFieldLocal - This projects the given function of the input fields into the function space provided, putting the coefficients in a local vector.
7914: Not collective
7916: Input Parameters:
7917: + dm - The DM
7918: . time - The time
7919: . localU - The input field vector
7920: . funcs - The functions to evaluate, one per field
7921: - mode - The insertion mode for values
7923: Output Parameter:
7924: . localX - The output vector
7926: Calling sequence of func:
7927: $ func(PetscInt dim, PetscInt Nf, PetscInt NfAux,
7928: $ const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
7929: $ const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
7930: $ PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f[]);
7932: + dim - The spatial dimension
7933: . Nf - The number of input fields
7934: . NfAux - The number of input auxiliary fields
7935: . uOff - The offset of each field in u[]
7936: . uOff_x - The offset of each field in u_x[]
7937: . u - The field values at this point in space
7938: . u_t - The field time derivative at this point in space (or NULL)
7939: . u_x - The field derivatives at this point in space
7940: . aOff - The offset of each auxiliary field in u[]
7941: . aOff_x - The offset of each auxiliary field in u_x[]
7942: . a - The auxiliary field values at this point in space
7943: . a_t - The auxiliary field time derivative at this point in space (or NULL)
7944: . a_x - The auxiliary field derivatives at this point in space
7945: . t - The current time
7946: . x - The coordinates of this point
7947: . numConstants - The number of constants
7948: . constants - The value of each constant
7949: - f - The value of the function at this point in space
7951: Note: There are three different DMs that potentially interact in this function. The output DM, dm, specifies the layout of the values calculates by funcs.
7952: The input DM, attached to U, may be different. For example, you can input the solution over the full domain, but output over a piece of the boundary, or
7953: a subdomain. You can also output a different number of fields than the input, with different discretizations. Last the auxiliary DM, attached to the
7954: auxiliary field vector, which is attached to dm, can also be different. It can have a different topology, number of fields, and discretizations.
7956: Level: intermediate
7958: .seealso: DMProjectField(), DMProjectFieldLabelLocal(), DMProjectFunction(), DMComputeL2Diff()
7959: @*/
7960: PetscErrorCode DMProjectFieldLocal(DM dm, PetscReal time, Vec localU,
7961: void (**funcs)(PetscInt, PetscInt, PetscInt,
7962: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
7963: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
7964: PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]),
7965: InsertMode mode, Vec localX)
7966: {
7973: if (!dm->ops->projectfieldlocal) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMProjectFieldLocal",((PetscObject)dm)->type_name);
7974: (dm->ops->projectfieldlocal) (dm, time, localU, funcs, mode, localX);
7975: return(0);
7976: }
7978: /*@C
7979: DMProjectFieldLabelLocal - This projects the given function of the input fields into the function space provided, putting the coefficients in a local vector, calculating only over the portion of the domain specified by the label.
7981: Not collective
7983: Input Parameters:
7984: + dm - The DM
7985: . time - The time
7986: . label - The DMLabel marking the portion of the domain to output
7987: . numIds - The number of label ids to use
7988: . ids - The label ids to use for marking
7989: . Nc - The number of components to set in the output, or PETSC_DETERMINE for all components
7990: . comps - The components to set in the output, or NULL for all components
7991: . localU - The input field vector
7992: . funcs - The functions to evaluate, one per field
7993: - mode - The insertion mode for values
7995: Output Parameter:
7996: . localX - The output vector
7998: Calling sequence of func:
7999: $ func(PetscInt dim, PetscInt Nf, PetscInt NfAux,
8000: $ const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
8001: $ const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
8002: $ PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f[]);
8004: + dim - The spatial dimension
8005: . Nf - The number of input fields
8006: . NfAux - The number of input auxiliary fields
8007: . uOff - The offset of each field in u[]
8008: . uOff_x - The offset of each field in u_x[]
8009: . u - The field values at this point in space
8010: . u_t - The field time derivative at this point in space (or NULL)
8011: . u_x - The field derivatives at this point in space
8012: . aOff - The offset of each auxiliary field in u[]
8013: . aOff_x - The offset of each auxiliary field in u_x[]
8014: . a - The auxiliary field values at this point in space
8015: . a_t - The auxiliary field time derivative at this point in space (or NULL)
8016: . a_x - The auxiliary field derivatives at this point in space
8017: . t - The current time
8018: . x - The coordinates of this point
8019: . numConstants - The number of constants
8020: . constants - The value of each constant
8021: - f - The value of the function at this point in space
8023: Note: There are three different DMs that potentially interact in this function. The output DM, dm, specifies the layout of the values calculates by funcs.
8024: The input DM, attached to U, may be different. For example, you can input the solution over the full domain, but output over a piece of the boundary, or
8025: a subdomain. You can also output a different number of fields than the input, with different discretizations. Last the auxiliary DM, attached to the
8026: auxiliary field vector, which is attached to dm, can also be different. It can have a different topology, number of fields, and discretizations.
8028: Level: intermediate
8030: .seealso: DMProjectField(), DMProjectFieldLabelLocal(), DMProjectFunction(), DMComputeL2Diff()
8031: @*/
8032: PetscErrorCode DMProjectFieldLabelLocal(DM dm, PetscReal time, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscInt Nc, const PetscInt comps[], Vec localU,
8033: void (**funcs)(PetscInt, PetscInt, PetscInt,
8034: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
8035: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
8036: PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]),
8037: InsertMode mode, Vec localX)
8038: {
8045: if (!dm->ops->projectfieldlabellocal) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMProjectFieldLocal",((PetscObject)dm)->type_name);
8046: (dm->ops->projectfieldlabellocal)(dm, time, label, numIds, ids, Nc, comps, localU, funcs, mode, localX);
8047: return(0);
8048: }
8050: /*@C
8051: DMComputeL2Diff - This function computes the L_2 difference between a function u and an FEM interpolant solution u_h.
8053: Input Parameters:
8054: + dm - The DM
8055: . time - The time
8056: . funcs - The functions to evaluate for each field component
8057: . ctxs - Optional array of contexts to pass to each function, or NULL.
8058: - X - The coefficient vector u_h, a global vector
8060: Output Parameter:
8061: . diff - The diff ||u - u_h||_2
8063: Level: developer
8065: .seealso: DMProjectFunction(), DMComputeL2FieldDiff(), DMComputeL2GradientDiff()
8066: @*/
8067: PetscErrorCode DMComputeL2Diff(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff)
8068: {
8074: if (!dm->ops->computel2diff) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMComputeL2Diff",((PetscObject)dm)->type_name);
8075: (dm->ops->computel2diff)(dm,time,funcs,ctxs,X,diff);
8076: return(0);
8077: }
8079: /*@C
8080: DMComputeL2GradientDiff - This function computes the L_2 difference between the gradient of a function u and an FEM interpolant solution grad u_h.
8082: Collective on dm
8084: Input Parameters:
8085: + dm - The DM
8086: , time - The time
8087: . funcs - The gradient functions to evaluate for each field component
8088: . ctxs - Optional array of contexts to pass to each function, or NULL.
8089: . X - The coefficient vector u_h, a global vector
8090: - n - The vector to project along
8092: Output Parameter:
8093: . diff - The diff ||(grad u - grad u_h) . n||_2
8095: Level: developer
8097: .seealso: DMProjectFunction(), DMComputeL2Diff()
8098: @*/
8099: PetscErrorCode DMComputeL2GradientDiff(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], const PetscReal[], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, const PetscReal n[], PetscReal *diff)
8100: {
8106: if (!dm->ops->computel2gradientdiff) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMComputeL2GradientDiff",((PetscObject)dm)->type_name);
8107: (dm->ops->computel2gradientdiff)(dm,time,funcs,ctxs,X,n,diff);
8108: return(0);
8109: }
8111: /*@C
8112: DMComputeL2FieldDiff - This function computes the L_2 difference between a function u and an FEM interpolant solution u_h, separated into field components.
8114: Collective on dm
8116: Input Parameters:
8117: + dm - The DM
8118: . time - The time
8119: . funcs - The functions to evaluate for each field component
8120: . ctxs - Optional array of contexts to pass to each function, or NULL.
8121: - X - The coefficient vector u_h, a global vector
8123: Output Parameter:
8124: . diff - The array of differences, ||u^f - u^f_h||_2
8126: Level: developer
8128: .seealso: DMProjectFunction(), DMComputeL2FieldDiff(), DMComputeL2GradientDiff()
8129: @*/
8130: PetscErrorCode DMComputeL2FieldDiff(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal diff[])
8131: {
8137: if (!dm->ops->computel2fielddiff) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMComputeL2FieldDiff",((PetscObject)dm)->type_name);
8138: (dm->ops->computel2fielddiff)(dm,time,funcs,ctxs,X,diff);
8139: return(0);
8140: }
8142: /*@C
8143: DMAdaptLabel - Adapt a dm based on a label with values interpreted as coarsening and refining flags. Specific implementations of DM maybe have
8144: specialized flags, but all implementations should accept flag values DM_ADAPT_DETERMINE, DM_ADAPT_KEEP, DM_ADAPT_REFINE, and DM_ADAPT_COARSEN.
8146: Collective on dm
8148: Input parameters:
8149: + dm - the pre-adaptation DM object
8150: - label - label with the flags
8152: Output parameters:
8153: . dmAdapt - the adapted DM object: may be NULL if an adapted DM could not be produced.
8155: Level: intermediate
8157: .seealso: DMAdaptMetric(), DMCoarsen(), DMRefine()
8158: @*/
8159: PetscErrorCode DMAdaptLabel(DM dm, DMLabel label, DM *dmAdapt)
8160: {
8167: *dmAdapt = NULL;
8168: if (!dm->ops->adaptlabel) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMAdaptLabel",((PetscObject)dm)->type_name);
8169: (dm->ops->adaptlabel)(dm, label, dmAdapt);
8170: return(0);
8171: }
8173: /*@C
8174: DMAdaptMetric - Generates a mesh adapted to the specified metric field using the pragmatic library.
8176: Input Parameters:
8177: + dm - The DM object
8178: . metric - The metric to which the mesh is adapted, defined vertex-wise.
8179: - bdLabel - Label for boundary tags, which will be preserved in the output mesh. bdLabel should be NULL if there is no such label, and should be different from "_boundary_".
8181: Output Parameter:
8182: . dmAdapt - Pointer to the DM object containing the adapted mesh
8184: Note: The label in the adapted mesh will be registered under the name of the input DMLabel object
8186: Level: advanced
8188: .seealso: DMAdaptLabel(), DMCoarsen(), DMRefine()
8189: @*/
8190: PetscErrorCode DMAdaptMetric(DM dm, Vec metric, DMLabel bdLabel, DM *dmAdapt)
8191: {
8199: *dmAdapt = NULL;
8200: if (!dm->ops->adaptmetric) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMAdaptMetric",((PetscObject)dm)->type_name);
8201: (dm->ops->adaptmetric)(dm, metric, bdLabel, dmAdapt);
8202: return(0);
8203: }
8205: /*@C
8206: DMGetNeighbors - Gets an array containing the MPI rank of all the processes neighbors
8208: Not Collective
8210: Input Parameter:
8211: . dm - The DM
8213: Output Parameters:
8214: + nranks - the number of neighbours
8215: - ranks - the neighbors ranks
8217: Notes:
8218: Do not free the array, it is freed when the DM is destroyed.
8220: Level: beginner
8222: .seealso: DMDAGetNeighbors(), PetscSFGetRootRanks()
8223: @*/
8224: PetscErrorCode DMGetNeighbors(DM dm,PetscInt *nranks,const PetscMPIInt *ranks[])
8225: {
8230: if (!dm->ops->getneighbors) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMGetNeighbors",((PetscObject)dm)->type_name);
8231: (dm->ops->getneighbors)(dm,nranks,ranks);
8232: return(0);
8233: }
8235: #include <petsc/private/matimpl.h>
8237: /*
8238: Converts the input vector to a ghosted vector and then calls the standard coloring code.
8239: This has be a different function because it requires DM which is not defined in the Mat library
8240: */
8241: PetscErrorCode MatFDColoringApply_AIJDM(Mat J,MatFDColoring coloring,Vec x1,void *sctx)
8242: {
8246: if (coloring->ctype == IS_COLORING_LOCAL) {
8247: Vec x1local;
8248: DM dm;
8249: MatGetDM(J,&dm);
8250: if (!dm) SETERRQ(PetscObjectComm((PetscObject)J),PETSC_ERR_ARG_INCOMP,"IS_COLORING_LOCAL requires a DM");
8251: DMGetLocalVector(dm,&x1local);
8252: DMGlobalToLocalBegin(dm,x1,INSERT_VALUES,x1local);
8253: DMGlobalToLocalEnd(dm,x1,INSERT_VALUES,x1local);
8254: x1 = x1local;
8255: }
8256: MatFDColoringApply_AIJ(J,coloring,x1,sctx);
8257: if (coloring->ctype == IS_COLORING_LOCAL) {
8258: DM dm;
8259: MatGetDM(J,&dm);
8260: DMRestoreLocalVector(dm,&x1);
8261: }
8262: return(0);
8263: }
8265: /*@
8266: MatFDColoringUseDM - allows a MatFDColoring object to use the DM associated with the matrix to use a IS_COLORING_LOCAL coloring
8268: Input Parameter:
8269: . coloring - the MatFDColoring object
8271: Developer Notes:
8272: this routine exists because the PETSc Mat library does not know about the DM objects
8274: Level: advanced
8276: .seealso: MatFDColoring, MatFDColoringCreate(), ISColoringType
8277: @*/
8278: PetscErrorCode MatFDColoringUseDM(Mat coloring,MatFDColoring fdcoloring)
8279: {
8281: coloring->ops->fdcoloringapply = MatFDColoringApply_AIJDM;
8282: return(0);
8283: }
8285: /*@
8286: DMGetCompatibility - determine if two DMs are compatible
8288: Collective
8290: Input Parameters:
8291: + dm1 - the first DM
8292: - dm2 - the second DM
8294: Output Parameters:
8295: + compatible - whether or not the two DMs are compatible
8296: - set - whether or not the compatible value was set
8298: Notes:
8299: Two DMs are deemed compatible if they represent the same parallel decomposition
8300: of the same topology. This implies that the section (field data) on one
8301: "makes sense" with respect to the topology and parallel decomposition of the other.
8302: Loosely speaking, compatible DMs represent the same domain and parallel
8303: decomposition, but hold different data.
8305: Typically, one would confirm compatibility if intending to simultaneously iterate
8306: over a pair of vectors obtained from different DMs.
8308: For example, two DMDA objects are compatible if they have the same local
8309: and global sizes and the same stencil width. They can have different numbers
8310: of degrees of freedom per node. Thus, one could use the node numbering from
8311: either DM in bounds for a loop over vectors derived from either DM.
8313: Consider the operation of summing data living on a 2-dof DMDA to data living
8314: on a 1-dof DMDA, which should be compatible, as in the following snippet.
8315: .vb
8316: ...
8317: DMGetCompatibility(da1,da2,&compatible,&set);
8318: if (set && compatible) {
8319: DMDAVecGetArrayDOF(da1,vec1,&arr1);
8320: DMDAVecGetArrayDOF(da2,vec2,&arr2);
8321: DMDAGetCorners(da1,&x,&y,NULL,&m,&n,NULL);
8322: for (j=y; j<y+n; ++j) {
8323: for (i=x; i<x+m, ++i) {
8324: arr1[j][i][0] = arr2[j][i][0] + arr2[j][i][1];
8325: }
8326: }
8327: DMDAVecRestoreArrayDOF(da1,vec1,&arr1);
8328: DMDAVecRestoreArrayDOF(da2,vec2,&arr2);
8329: } else {
8330: SETERRQ(PetscObjectComm((PetscObject)da1,PETSC_ERR_ARG_INCOMP,"DMDA objects incompatible");
8331: }
8332: ...
8333: .ve
8335: Checking compatibility might be expensive for a given implementation of DM,
8336: or might be impossible to unambiguously confirm or deny. For this reason,
8337: this function may decline to determine compatibility, and hence users should
8338: always check the "set" output parameter.
8340: A DM is always compatible with itself.
8342: In the current implementation, DMs which live on "unequal" communicators
8343: (MPI_UNEQUAL in the terminology of MPI_Comm_compare()) are always deemed
8344: incompatible.
8346: This function is labeled "Collective," as information about all subdomains
8347: is required on each rank. However, in DM implementations which store all this
8348: information locally, this function may be merely "Logically Collective".
8350: Developer Notes:
8351: Compatibility is assumed to be a symmetric concept; DM A is compatible with DM B
8352: iff B is compatible with A. Thus, this function checks the implementations
8353: of both dm and dmc (if they are of different types), attempting to determine
8354: compatibility. It is left to DM implementers to ensure that symmetry is
8355: preserved. The simplest way to do this is, when implementing type-specific
8356: logic for this function, is to check for existing logic in the implementation
8357: of other DM types and let *set = PETSC_FALSE if found.
8359: Level: advanced
8361: .seealso: DM, DMDACreateCompatibleDMDA(), DMStagCreateCompatibleDMStag()
8362: @*/
8364: PetscErrorCode DMGetCompatibility(DM dm1,DM dm2,PetscBool *compatible,PetscBool *set)
8365: {
8367: PetscMPIInt compareResult;
8368: DMType type,type2;
8369: PetscBool sameType;
8375: /* Declare a DM compatible with itself */
8376: if (dm1 == dm2) {
8377: *set = PETSC_TRUE;
8378: *compatible = PETSC_TRUE;
8379: return(0);
8380: }
8382: /* Declare a DM incompatible with a DM that lives on an "unequal"
8383: communicator. Note that this does not preclude compatibility with
8384: DMs living on "congruent" or "similar" communicators, but this must be
8385: determined by the implementation-specific logic */
8386: MPI_Comm_compare(PetscObjectComm((PetscObject)dm1),PetscObjectComm((PetscObject)dm2),&compareResult);
8387: if (compareResult == MPI_UNEQUAL) {
8388: *set = PETSC_TRUE;
8389: *compatible = PETSC_FALSE;
8390: return(0);
8391: }
8393: /* Pass to the implementation-specific routine, if one exists. */
8394: if (dm1->ops->getcompatibility) {
8395: (*dm1->ops->getcompatibility)(dm1,dm2,compatible,set);
8396: if (*set) return(0);
8397: }
8399: /* If dm1 and dm2 are of different types, then attempt to check compatibility
8400: with an implementation of this function from dm2 */
8401: DMGetType(dm1,&type);
8402: DMGetType(dm2,&type2);
8403: PetscStrcmp(type,type2,&sameType);
8404: if (!sameType && dm2->ops->getcompatibility) {
8405: (*dm2->ops->getcompatibility)(dm2,dm1,compatible,set); /* Note argument order */
8406: } else {
8407: *set = PETSC_FALSE;
8408: }
8409: return(0);
8410: }
8412: /*@C
8413: DMMonitorSet - Sets an ADDITIONAL function that is to be used after a solve to monitor discretization performance.
8415: Logically Collective on DM
8417: Input Parameters:
8418: + DM - the DM
8419: . f - the monitor function
8420: . mctx - [optional] user-defined context for private data for the monitor routine (use NULL if no context is desired)
8421: - monitordestroy - [optional] routine that frees monitor context (may be NULL)
8423: Options Database Keys:
8424: - -dm_monitor_cancel - cancels all monitors that have been hardwired into a code by calls to DMMonitorSet(), but
8425: does not cancel those set via the options database.
8427: Notes:
8428: Several different monitoring routines may be set by calling
8429: DMMonitorSet() multiple times; all will be called in the
8430: order in which they were set.
8432: Fortran Notes:
8433: Only a single monitor function can be set for each DM object
8435: Level: intermediate
8437: .seealso: DMMonitorCancel()
8438: @*/
8439: PetscErrorCode DMMonitorSet(DM dm, PetscErrorCode (*f)(DM, void *), void *mctx, PetscErrorCode (*monitordestroy)(void**))
8440: {
8441: PetscInt m;
8446: for (m = 0; m < dm->numbermonitors; ++m) {
8447: PetscBool identical;
8449: PetscMonitorCompare((PetscErrorCode (*)(void)) f, mctx, monitordestroy, (PetscErrorCode (*)(void)) dm->monitor[m], dm->monitorcontext[m], dm->monitordestroy[m], &identical);
8450: if (identical) return(0);
8451: }
8452: if (dm->numbermonitors >= MAXDMMONITORS) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Too many monitors set");
8453: dm->monitor[dm->numbermonitors] = f;
8454: dm->monitordestroy[dm->numbermonitors] = monitordestroy;
8455: dm->monitorcontext[dm->numbermonitors++] = (void *) mctx;
8456: return(0);
8457: }
8459: /*@
8460: DMMonitorCancel - Clears all the monitor functions for a DM object.
8462: Logically Collective on DM
8464: Input Parameter:
8465: . dm - the DM
8467: Options Database Key:
8468: . -dm_monitor_cancel - cancels all monitors that have been hardwired
8469: into a code by calls to DMonitorSet(), but does not cancel those
8470: set via the options database
8472: Notes:
8473: There is no way to clear one specific monitor from a DM object.
8475: Level: intermediate
8477: .seealso: DMMonitorSet()
8478: @*/
8479: PetscErrorCode DMMonitorCancel(DM dm)
8480: {
8482: PetscInt m;
8486: for (m = 0; m < dm->numbermonitors; ++m) {
8487: if (dm->monitordestroy[m]) {(*dm->monitordestroy[m])(&dm->monitorcontext[m]);}
8488: }
8489: dm->numbermonitors = 0;
8490: return(0);
8491: }
8493: /*@C
8494: DMMonitorSetFromOptions - Sets a monitor function and viewer appropriate for the type indicated by the user
8496: Collective on DM
8498: Input Parameters:
8499: + dm - DM object you wish to monitor
8500: . name - the monitor type one is seeking
8501: . help - message indicating what monitoring is done
8502: . manual - manual page for the monitor
8503: . monitor - the monitor function
8504: - monitorsetup - a function that is called once ONLY if the user selected this monitor that may set additional features of the DM or PetscViewer objects
8506: Output Parameter:
8507: . flg - Flag set if the monitor was created
8509: Level: developer
8511: .seealso: PetscOptionsGetViewer(), PetscOptionsGetReal(), PetscOptionsHasName(), PetscOptionsGetString(),
8512: PetscOptionsGetIntArray(), PetscOptionsGetRealArray(), PetscOptionsBool()
8513: PetscOptionsInt(), PetscOptionsString(), PetscOptionsReal(), PetscOptionsBool(),
8514: PetscOptionsName(), PetscOptionsBegin(), PetscOptionsEnd(), PetscOptionsHead(),
8515: PetscOptionsStringArray(),PetscOptionsRealArray(), PetscOptionsScalar(),
8516: PetscOptionsBoolGroupBegin(), PetscOptionsBoolGroup(), PetscOptionsBoolGroupEnd(),
8517: PetscOptionsFList(), PetscOptionsEList()
8518: @*/
8519: PetscErrorCode DMMonitorSetFromOptions(DM dm, const char name[], const char help[], const char manual[], PetscErrorCode (*monitor)(DM, void *), PetscErrorCode (*monitorsetup)(DM, PetscViewerAndFormat *), PetscBool *flg)
8520: {
8521: PetscViewer viewer;
8522: PetscViewerFormat format;
8523: PetscErrorCode ierr;
8527: PetscOptionsGetViewer(PetscObjectComm((PetscObject) dm), ((PetscObject) dm)->options, ((PetscObject) dm)->prefix, name, &viewer, &format, flg);
8528: if (*flg) {
8529: PetscViewerAndFormat *vf;
8531: PetscViewerAndFormatCreate(viewer, format, &vf);
8532: PetscObjectDereference((PetscObject) viewer);
8533: if (monitorsetup) {(*monitorsetup)(dm, vf);}
8534: DMMonitorSet(dm,(PetscErrorCode (*)(DM, void *)) monitor, vf, (PetscErrorCode (*)(void **)) PetscViewerAndFormatDestroy);
8535: }
8536: return(0);
8537: }
8539: /*@
8540: DMMonitor - runs the user provided monitor routines, if they exist
8542: Collective on DM
8544: Input Parameters:
8545: . dm - The DM
8547: Level: developer
8549: .seealso: DMMonitorSet()
8550: @*/
8551: PetscErrorCode DMMonitor(DM dm)
8552: {
8553: PetscInt m;
8557: if (!dm) return(0);
8559: for (m = 0; m < dm->numbermonitors; ++m) {
8560: (*dm->monitor[m])(dm, dm->monitorcontext[m]);
8561: }
8562: return(0);
8563: }