Actual source code: dm.c
petsc-3.12.5 2020-03-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;
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};
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->localSection = NULL;
59: v->globalSection = NULL;
60: v->defaultConstraintSection = NULL;
61: v->defaultConstraintMat = NULL;
62: v->L = NULL;
63: v->maxCell = NULL;
64: v->bdtype = NULL;
65: v->dimEmbed = PETSC_DEFAULT;
66: v->dim = PETSC_DETERMINE;
67: {
68: PetscInt i;
69: for (i = 0; i < 10; ++i) {
70: v->nullspaceConstructors[i] = NULL;
71: v->nearnullspaceConstructors[i] = NULL;
72: }
73: }
74: PetscDSCreate(PetscObjectComm((PetscObject) v), &ds);
75: DMSetRegionDS(v, NULL, NULL, ds);
76: PetscDSDestroy(&ds);
77: v->dmBC = NULL;
78: v->coarseMesh = NULL;
79: v->outputSequenceNum = -1;
80: v->outputSequenceVal = 0.0;
81: DMSetVecType(v,VECSTANDARD);
82: DMSetMatType(v,MATAIJ);
83: PetscNew(&(v->labels));
84: v->labels->refct = 1;
86: *dm = v;
87: return(0);
88: }
90: /*@
91: DMClone - Creates a DM object with the same topology as the original.
93: Collective
95: Input Parameter:
96: . dm - The original DM object
98: Output Parameter:
99: . newdm - The new DM object
101: Level: beginner
103: Notes: For some DM 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
105: share the PetscSection of the original DM
107: .seealso: DMDestry(), DMCreate(), DMSetType(), DMSetLocalSection(), DMSetGlobalSection()
109: @*/
110: PetscErrorCode DMClone(DM dm, DM *newdm)
111: {
112: PetscSF sf;
113: Vec coords;
114: void *ctx;
115: PetscInt dim, cdim;
121: DMCreate(PetscObjectComm((PetscObject) dm), newdm);
122: PetscFree((*newdm)->labels);
123: dm->labels->refct++;
124: (*newdm)->labels = dm->labels;
125: (*newdm)->depthLabel = dm->depthLabel;
126: (*newdm)->leveldown = dm->leveldown;
127: (*newdm)->levelup = dm->levelup;
128: DMGetDimension(dm, &dim);
129: DMSetDimension(*newdm, dim);
130: if (dm->ops->clone) {
131: (*dm->ops->clone)(dm, newdm);
132: }
133: (*newdm)->setupcalled = dm->setupcalled;
134: DMGetPointSF(dm, &sf);
135: DMSetPointSF(*newdm, sf);
136: DMGetApplicationContext(dm, &ctx);
137: DMSetApplicationContext(*newdm, ctx);
138: if (dm->coordinateDM) {
139: DM ncdm;
140: PetscSection cs;
141: PetscInt pEnd = -1, pEndMax = -1;
143: DMGetLocalSection(dm->coordinateDM, &cs);
144: if (cs) {PetscSectionGetChart(cs, NULL, &pEnd);}
145: MPI_Allreduce(&pEnd,&pEndMax,1,MPIU_INT,MPI_MAX,PetscObjectComm((PetscObject)dm));
146: if (pEndMax >= 0) {
147: DMClone(dm->coordinateDM, &ncdm);
148: DMCopyDisc(dm->coordinateDM, ncdm);
149: DMSetLocalSection(ncdm, cs);
150: DMSetCoordinateDM(*newdm, ncdm);
151: DMDestroy(&ncdm);
152: }
153: }
154: DMGetCoordinateDim(dm, &cdim);
155: DMSetCoordinateDim(*newdm, cdim);
156: DMGetCoordinatesLocal(dm, &coords);
157: if (coords) {
158: DMSetCoordinatesLocal(*newdm, coords);
159: } else {
160: DMGetCoordinates(dm, &coords);
161: if (coords) {DMSetCoordinates(*newdm, coords);}
162: }
163: {
164: PetscBool isper;
165: const PetscReal *maxCell, *L;
166: const DMBoundaryType *bd;
167: DMGetPeriodicity(dm, &isper, &maxCell, &L, &bd);
168: DMSetPeriodicity(*newdm, isper, maxCell, L, bd);
169: }
170: {
171: PetscBool useCone, useClosure;
173: DMGetAdjacency(dm, PETSC_DEFAULT, &useCone, &useClosure);
174: DMSetAdjacency(*newdm, PETSC_DEFAULT, useCone, useClosure);
175: }
176: return(0);
177: }
179: /*@C
180: DMSetVecType - Sets the type of vector created with DMCreateLocalVector() and DMCreateGlobalVector()
182: Logically Collective on da
184: Input Parameter:
185: + da - initial distributed array
186: . ctype - the vector type, currently either VECSTANDARD, VECCUDA, or VECVIENNACL
188: Options Database:
189: . -dm_vec_type ctype
191: Level: intermediate
193: .seealso: DMCreate(), DMDestroy(), DM, DMDAInterpolationType, VecType, DMGetVecType(), DMSetMatType(), DMGetMatType()
194: @*/
195: PetscErrorCode DMSetVecType(DM da,VecType ctype)
196: {
201: PetscFree(da->vectype);
202: PetscStrallocpy(ctype,(char**)&da->vectype);
203: return(0);
204: }
206: /*@C
207: DMGetVecType - Gets the type of vector created with DMCreateLocalVector() and DMCreateGlobalVector()
209: Logically Collective on da
211: Input Parameter:
212: . da - initial distributed array
214: Output Parameter:
215: . ctype - the vector type
217: Level: intermediate
219: .seealso: DMCreate(), DMDestroy(), DM, DMDAInterpolationType, VecType, DMSetMatType(), DMGetMatType(), DMSetVecType()
220: @*/
221: PetscErrorCode DMGetVecType(DM da,VecType *ctype)
222: {
225: *ctype = da->vectype;
226: return(0);
227: }
229: /*@
230: VecGetDM - Gets the DM defining the data layout of the vector
232: Not collective
234: Input Parameter:
235: . v - The Vec
237: Output Parameter:
238: . dm - The DM
240: Level: intermediate
242: .seealso: VecSetDM(), DMGetLocalVector(), DMGetGlobalVector(), DMSetVecType()
243: @*/
244: PetscErrorCode VecGetDM(Vec v, DM *dm)
245: {
251: PetscObjectQuery((PetscObject) v, "__PETSc_dm", (PetscObject*) dm);
252: return(0);
253: }
255: /*@
256: VecSetDM - Sets the DM defining the data layout of the vector.
258: Not collective
260: Input Parameters:
261: + v - The Vec
262: - dm - The DM
264: 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.
266: Level: intermediate
268: .seealso: VecGetDM(), DMGetLocalVector(), DMGetGlobalVector(), DMSetVecType()
269: @*/
270: PetscErrorCode VecSetDM(Vec v, DM dm)
271: {
277: PetscObjectCompose((PetscObject) v, "__PETSc_dm", (PetscObject) dm);
278: return(0);
279: }
281: /*@C
282: DMSetISColoringType - Sets the type of coloring, global or local, that is created by the DM
284: Logically Collective on dm
286: Input Parameters:
287: + dm - the DM context
288: - ctype - the matrix type
290: Options Database:
291: . -dm_is_coloring_type - global or local
293: Level: intermediate
295: .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMCreateMatrix(), DMSetMatrixPreallocateOnly(), MatType, DMGetMatType(),
296: DMGetISColoringType()
297: @*/
298: PetscErrorCode DMSetISColoringType(DM dm,ISColoringType ctype)
299: {
302: dm->coloringtype = ctype;
303: return(0);
304: }
306: /*@C
307: DMGetISColoringType - Gets the type of coloring, global or local, that is created by the DM
309: Logically Collective on dm
311: Input Parameter:
312: . dm - the DM context
314: Output Parameter:
315: . ctype - the matrix type
317: Options Database:
318: . -dm_is_coloring_type - global or local
320: Level: intermediate
322: .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMCreateMatrix(), DMSetMatrixPreallocateOnly(), MatType, DMGetMatType(),
323: DMGetISColoringType()
324: @*/
325: PetscErrorCode DMGetISColoringType(DM dm,ISColoringType *ctype)
326: {
329: *ctype = dm->coloringtype;
330: return(0);
331: }
333: /*@C
334: DMSetMatType - Sets the type of matrix created with DMCreateMatrix()
336: Logically Collective on dm
338: Input Parameters:
339: + dm - the DM context
340: - ctype - the matrix type
342: Options Database:
343: . -dm_mat_type ctype
345: Level: intermediate
347: .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMCreateMatrix(), DMSetMatrixPreallocateOnly(), MatType, DMGetMatType(), DMSetMatType(), DMGetMatType()
348: @*/
349: PetscErrorCode DMSetMatType(DM dm,MatType ctype)
350: {
355: PetscFree(dm->mattype);
356: PetscStrallocpy(ctype,(char**)&dm->mattype);
357: return(0);
358: }
360: /*@C
361: DMGetMatType - Gets the type of matrix created with DMCreateMatrix()
363: Logically Collective on dm
365: Input Parameter:
366: . dm - the DM context
368: Output Parameter:
369: . ctype - the matrix type
371: Options Database:
372: . -dm_mat_type ctype
374: Level: intermediate
376: .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMCreateMatrix(), DMSetMatrixPreallocateOnly(), MatType, DMSetMatType(), DMSetMatType(), DMGetMatType()
377: @*/
378: PetscErrorCode DMGetMatType(DM dm,MatType *ctype)
379: {
382: *ctype = dm->mattype;
383: return(0);
384: }
386: /*@
387: MatGetDM - Gets the DM defining the data layout of the matrix
389: Not collective
391: Input Parameter:
392: . A - The Mat
394: Output Parameter:
395: . dm - The DM
397: Level: intermediate
399: Developer Note: Since the Mat class doesn't know about the DM class the DM object is associated with
400: the Mat through a PetscObjectCompose() operation
402: .seealso: MatSetDM(), DMCreateMatrix(), DMSetMatType()
403: @*/
404: PetscErrorCode MatGetDM(Mat A, DM *dm)
405: {
411: PetscObjectQuery((PetscObject) A, "__PETSc_dm", (PetscObject*) dm);
412: return(0);
413: }
415: /*@
416: MatSetDM - Sets the DM defining the data layout of the matrix
418: Not collective
420: Input Parameters:
421: + A - The Mat
422: - dm - The DM
424: Level: intermediate
426: Developer Note: Since the Mat class doesn't know about the DM class the DM object is associated with
427: the Mat through a PetscObjectCompose() operation
430: .seealso: MatGetDM(), DMCreateMatrix(), DMSetMatType()
431: @*/
432: PetscErrorCode MatSetDM(Mat A, DM dm)
433: {
439: PetscObjectCompose((PetscObject) A, "__PETSc_dm", (PetscObject) dm);
440: return(0);
441: }
443: /*@C
444: DMSetOptionsPrefix - Sets the prefix used for searching for all
445: DM options in the database.
447: Logically Collective on dm
449: Input Parameter:
450: + da - the DM context
451: - prefix - the prefix to prepend to all option names
453: Notes:
454: A hyphen (-) must NOT be given at the beginning of the prefix name.
455: The first character of all runtime options is AUTOMATICALLY the hyphen.
457: Level: advanced
459: .seealso: DMSetFromOptions()
460: @*/
461: PetscErrorCode DMSetOptionsPrefix(DM dm,const char prefix[])
462: {
467: PetscObjectSetOptionsPrefix((PetscObject)dm,prefix);
468: if (dm->sf) {
469: PetscObjectSetOptionsPrefix((PetscObject)dm->sf,prefix);
470: }
471: if (dm->sectionSF) {
472: PetscObjectSetOptionsPrefix((PetscObject)dm->sectionSF,prefix);
473: }
474: return(0);
475: }
477: /*@C
478: DMAppendOptionsPrefix - Appends to the prefix used for searching for all
479: DM options in the database.
481: Logically Collective on dm
483: Input Parameters:
484: + dm - the DM context
485: - prefix - the prefix string to prepend to all DM option requests
487: Notes:
488: A hyphen (-) must NOT be given at the beginning of the prefix name.
489: The first character of all runtime options is AUTOMATICALLY the hyphen.
491: Level: advanced
493: .seealso: DMSetOptionsPrefix(), DMGetOptionsPrefix()
494: @*/
495: PetscErrorCode DMAppendOptionsPrefix(DM dm,const char prefix[])
496: {
501: PetscObjectAppendOptionsPrefix((PetscObject)dm,prefix);
502: return(0);
503: }
505: /*@C
506: DMGetOptionsPrefix - Gets the prefix used for searching for all
507: DM options in the database.
509: Not Collective
511: Input Parameters:
512: . dm - the DM context
514: Output Parameters:
515: . prefix - pointer to the prefix string used is returned
517: Notes:
518: On the fortran side, the user should pass in a string 'prefix' of
519: sufficient length to hold the prefix.
521: Level: advanced
523: .seealso: DMSetOptionsPrefix(), DMAppendOptionsPrefix()
524: @*/
525: PetscErrorCode DMGetOptionsPrefix(DM dm,const char *prefix[])
526: {
531: PetscObjectGetOptionsPrefix((PetscObject)dm,prefix);
532: return(0);
533: }
535: static PetscErrorCode DMCountNonCyclicReferences(DM dm, PetscBool recurseCoarse, PetscBool recurseFine, PetscInt *ncrefct)
536: {
537: PetscInt i, refct = ((PetscObject) dm)->refct;
538: DMNamedVecLink nlink;
542: *ncrefct = 0;
543: /* count all the circular references of DM and its contained Vecs */
544: for (i=0; i<DM_MAX_WORK_VECTORS; i++) {
545: if (dm->localin[i]) refct--;
546: if (dm->globalin[i]) refct--;
547: }
548: for (nlink=dm->namedglobal; nlink; nlink=nlink->next) refct--;
549: for (nlink=dm->namedlocal; nlink; nlink=nlink->next) refct--;
550: if (dm->x) {
551: DM obj;
552: VecGetDM(dm->x, &obj);
553: if (obj == dm) refct--;
554: }
555: if (dm->coarseMesh && dm->coarseMesh->fineMesh == dm) {
556: refct--;
557: if (recurseCoarse) {
558: PetscInt coarseCount;
560: DMCountNonCyclicReferences(dm->coarseMesh, PETSC_TRUE, PETSC_FALSE,&coarseCount);
561: refct += coarseCount;
562: }
563: }
564: if (dm->fineMesh && dm->fineMesh->coarseMesh == dm) {
565: refct--;
566: if (recurseFine) {
567: PetscInt fineCount;
569: DMCountNonCyclicReferences(dm->fineMesh, PETSC_FALSE, PETSC_TRUE,&fineCount);
570: refct += fineCount;
571: }
572: }
573: *ncrefct = refct;
574: return(0);
575: }
577: PetscErrorCode DMDestroyLabelLinkList_Internal(DM dm)
578: {
582: if (!--(dm->labels->refct)) {
583: DMLabelLink next = dm->labels->next;
585: /* destroy the labels */
586: while (next) {
587: DMLabelLink tmp = next->next;
589: DMLabelDestroy(&next->label);
590: PetscFree(next);
591: next = tmp;
592: }
593: PetscFree(dm->labels);
594: }
595: return(0);
596: }
598: /*@
599: DMDestroy - Destroys a vector packer or DM.
601: Collective on dm
603: Input Parameter:
604: . dm - the DM object to destroy
606: Level: developer
608: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix()
610: @*/
611: PetscErrorCode DMDestroy(DM *dm)
612: {
613: PetscInt i, cnt;
614: DMNamedVecLink nlink,nnext;
618: if (!*dm) return(0);
621: /* count all non-cyclic references in the doubly-linked list of coarse<->fine meshes */
622: DMCountNonCyclicReferences(*dm,PETSC_TRUE,PETSC_TRUE,&cnt);
623: --((PetscObject)(*dm))->refct;
624: if (--cnt > 0) {*dm = 0; return(0);}
625: /*
626: Need this test because the dm references the vectors that
627: reference the dm, so destroying the dm calls destroy on the
628: vectors that cause another destroy on the dm
629: */
630: if (((PetscObject)(*dm))->refct < 0) return(0);
631: ((PetscObject) (*dm))->refct = 0;
632: for (i=0; i<DM_MAX_WORK_VECTORS; i++) {
633: if ((*dm)->localout[i]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Destroying a DM that has a local vector obtained with DMGetLocalVector()");
634: VecDestroy(&(*dm)->localin[i]);
635: }
636: nnext=(*dm)->namedglobal;
637: (*dm)->namedglobal = NULL;
638: for (nlink=nnext; nlink; nlink=nnext) { /* Destroy the named vectors */
639: nnext = nlink->next;
640: if (nlink->status != DMVEC_STATUS_IN) SETERRQ1(((PetscObject)*dm)->comm,PETSC_ERR_ARG_WRONGSTATE,"DM still has Vec named '%s' checked out",nlink->name);
641: PetscFree(nlink->name);
642: VecDestroy(&nlink->X);
643: PetscFree(nlink);
644: }
645: nnext=(*dm)->namedlocal;
646: (*dm)->namedlocal = NULL;
647: for (nlink=nnext; nlink; nlink=nnext) { /* Destroy the named local vectors */
648: nnext = nlink->next;
649: if (nlink->status != DMVEC_STATUS_IN) SETERRQ1(((PetscObject)*dm)->comm,PETSC_ERR_ARG_WRONGSTATE,"DM still has Vec named '%s' checked out",nlink->name);
650: PetscFree(nlink->name);
651: VecDestroy(&nlink->X);
652: PetscFree(nlink);
653: }
655: /* Destroy the list of hooks */
656: {
657: DMCoarsenHookLink link,next;
658: for (link=(*dm)->coarsenhook; link; link=next) {
659: next = link->next;
660: PetscFree(link);
661: }
662: (*dm)->coarsenhook = NULL;
663: }
664: {
665: DMRefineHookLink link,next;
666: for (link=(*dm)->refinehook; link; link=next) {
667: next = link->next;
668: PetscFree(link);
669: }
670: (*dm)->refinehook = NULL;
671: }
672: {
673: DMSubDomainHookLink link,next;
674: for (link=(*dm)->subdomainhook; link; link=next) {
675: next = link->next;
676: PetscFree(link);
677: }
678: (*dm)->subdomainhook = NULL;
679: }
680: {
681: DMGlobalToLocalHookLink link,next;
682: for (link=(*dm)->gtolhook; link; link=next) {
683: next = link->next;
684: PetscFree(link);
685: }
686: (*dm)->gtolhook = NULL;
687: }
688: {
689: DMLocalToGlobalHookLink link,next;
690: for (link=(*dm)->ltoghook; link; link=next) {
691: next = link->next;
692: PetscFree(link);
693: }
694: (*dm)->ltoghook = NULL;
695: }
696: /* Destroy the work arrays */
697: {
698: DMWorkLink link,next;
699: if ((*dm)->workout) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Work array still checked out");
700: for (link=(*dm)->workin; link; link=next) {
701: next = link->next;
702: PetscFree(link->mem);
703: PetscFree(link);
704: }
705: (*dm)->workin = NULL;
706: }
707: /* destroy the labels */
708: DMDestroyLabelLinkList_Internal(*dm);
709: /* destroy the fields */
710: DMClearFields(*dm);
711: /* destroy the boundaries */
712: {
713: DMBoundary next = (*dm)->boundary;
714: while (next) {
715: DMBoundary b = next;
717: next = b->next;
718: PetscFree(b);
719: }
720: }
722: PetscObjectDestroy(&(*dm)->dmksp);
723: PetscObjectDestroy(&(*dm)->dmsnes);
724: PetscObjectDestroy(&(*dm)->dmts);
726: if ((*dm)->ctx && (*dm)->ctxdestroy) {
727: (*(*dm)->ctxdestroy)(&(*dm)->ctx);
728: }
729: VecDestroy(&(*dm)->x);
730: MatFDColoringDestroy(&(*dm)->fd);
731: DMClearGlobalVectors(*dm);
732: ISLocalToGlobalMappingDestroy(&(*dm)->ltogmap);
733: PetscFree((*dm)->vectype);
734: PetscFree((*dm)->mattype);
736: PetscSectionDestroy(&(*dm)->localSection);
737: PetscSectionDestroy(&(*dm)->globalSection);
738: PetscLayoutDestroy(&(*dm)->map);
739: PetscSectionDestroy(&(*dm)->defaultConstraintSection);
740: MatDestroy(&(*dm)->defaultConstraintMat);
741: PetscSFDestroy(&(*dm)->sf);
742: PetscSFDestroy(&(*dm)->sectionSF);
743: if ((*dm)->useNatural) {
744: if ((*dm)->sfNatural) {
745: PetscSFDestroy(&(*dm)->sfNatural);
746: }
747: PetscObjectDereference((PetscObject) (*dm)->sfMigration);
748: }
749: if ((*dm)->coarseMesh && (*dm)->coarseMesh->fineMesh == *dm) {
750: DMSetFineDM((*dm)->coarseMesh,NULL);
751: }
752: DMDestroy(&(*dm)->coarseMesh);
753: if ((*dm)->fineMesh && (*dm)->fineMesh->coarseMesh == *dm) {
754: DMSetCoarseDM((*dm)->fineMesh,NULL);
755: }
756: DMDestroy(&(*dm)->fineMesh);
757: DMFieldDestroy(&(*dm)->coordinateField);
758: DMDestroy(&(*dm)->coordinateDM);
759: VecDestroy(&(*dm)->coordinates);
760: VecDestroy(&(*dm)->coordinatesLocal);
761: PetscFree3((*dm)->L,(*dm)->maxCell,(*dm)->bdtype);
762: if ((*dm)->transformDestroy) {(*(*dm)->transformDestroy)(*dm, (*dm)->transformCtx);}
763: DMDestroy(&(*dm)->transformDM);
764: VecDestroy(&(*dm)->transform);
766: DMClearDS(*dm);
767: DMDestroy(&(*dm)->dmBC);
768: /* if memory was published with SAWs then destroy it */
769: PetscObjectSAWsViewOff((PetscObject)*dm);
771: if ((*dm)->ops->destroy) {
772: (*(*dm)->ops->destroy)(*dm);
773: }
774: /* We do not destroy (*dm)->data here so that we can reference count backend objects */
775: PetscHeaderDestroy(dm);
776: return(0);
777: }
779: /*@
780: DMSetUp - sets up the data structures inside a DM object
782: Collective on dm
784: Input Parameter:
785: . dm - the DM object to setup
787: Level: developer
789: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix()
791: @*/
792: PetscErrorCode DMSetUp(DM dm)
793: {
798: if (dm->setupcalled) return(0);
799: if (dm->ops->setup) {
800: (*dm->ops->setup)(dm);
801: }
802: dm->setupcalled = PETSC_TRUE;
803: return(0);
804: }
806: /*@
807: DMSetFromOptions - sets parameters in a DM from the options database
809: Collective on dm
811: Input Parameter:
812: . dm - the DM object to set options for
814: Options Database:
815: + -dm_preallocate_only - Only preallocate the matrix for DMCreateMatrix(), but do not fill it with zeros
816: . -dm_vec_type <type> - type of vector to create inside DM
817: . -dm_mat_type <type> - type of matrix to create inside DM
818: - -dm_is_coloring_type - <global or local>
820: Level: developer
822: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix()
824: @*/
825: PetscErrorCode DMSetFromOptions(DM dm)
826: {
827: char typeName[256];
828: PetscBool flg;
833: dm->setfromoptionscalled = PETSC_TRUE;
834: if (dm->sf) {PetscSFSetFromOptions(dm->sf);}
835: if (dm->sectionSF) {PetscSFSetFromOptions(dm->sectionSF);}
836: PetscObjectOptionsBegin((PetscObject)dm);
837: PetscOptionsBool("-dm_preallocate_only","only preallocate matrix, but do not set column indices","DMSetMatrixPreallocateOnly",dm->prealloc_only,&dm->prealloc_only,NULL);
838: PetscOptionsFList("-dm_vec_type","Vector type used for created vectors","DMSetVecType",VecList,dm->vectype,typeName,256,&flg);
839: if (flg) {
840: DMSetVecType(dm,typeName);
841: }
842: PetscOptionsFList("-dm_mat_type","Matrix type used for created matrices","DMSetMatType",MatList,dm->mattype ? dm->mattype : typeName,typeName,sizeof(typeName),&flg);
843: if (flg) {
844: DMSetMatType(dm,typeName);
845: }
846: PetscOptionsEnum("-dm_is_coloring_type","Global or local coloring of Jacobian","DMSetISColoringType",ISColoringTypes,(PetscEnum)dm->coloringtype,(PetscEnum*)&dm->coloringtype,NULL);
847: if (dm->ops->setfromoptions) {
848: (*dm->ops->setfromoptions)(PetscOptionsObject,dm);
849: }
850: /* process any options handlers added with PetscObjectAddOptionsHandler() */
851: PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject) dm);
852: PetscOptionsEnd();
853: return(0);
854: }
856: /*@C
857: DMView - Views a DM
859: Collective on dm
861: Input Parameter:
862: + dm - the DM object to view
863: - v - the viewer
865: Level: beginner
867: .seealso DMDestroy(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix()
869: @*/
870: PetscErrorCode DMView(DM dm,PetscViewer v)
871: {
872: PetscErrorCode ierr;
873: PetscBool isbinary;
874: PetscMPIInt size;
875: PetscViewerFormat format;
879: if (!v) {
880: PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)dm),&v);
881: }
883: /* Ideally, we would like to have this test on.
884: However, it currently breaks socket viz via GLVis.
885: During DMView(parallel_mesh,glvis_viewer), each
886: process opens a sequential ASCII socket to visualize
887: the local mesh, and PetscObjectView(dm,local_socket)
888: is internally called inside VecView_GLVis, incurring
889: in an error here */
891: PetscViewerCheckWritable(v);
893: PetscViewerGetFormat(v,&format);
894: MPI_Comm_size(PetscObjectComm((PetscObject)dm),&size);
895: if (size == 1 && format == PETSC_VIEWER_LOAD_BALANCE) return(0);
896: PetscObjectPrintClassNamePrefixType((PetscObject)dm,v);
897: PetscObjectTypeCompare((PetscObject)v,PETSCVIEWERBINARY,&isbinary);
898: if (isbinary) {
899: PetscInt classid = DM_FILE_CLASSID;
900: char type[256];
902: PetscViewerBinaryWrite(v,&classid,1,PETSC_INT,PETSC_FALSE);
903: PetscStrncpy(type,((PetscObject)dm)->type_name,256);
904: PetscViewerBinaryWrite(v,type,256,PETSC_CHAR,PETSC_FALSE);
905: }
906: if (dm->ops->view) {
907: (*dm->ops->view)(dm,v);
908: }
909: return(0);
910: }
912: /*@
913: DMCreateGlobalVector - Creates a global vector from a DM object
915: Collective on dm
917: Input Parameter:
918: . dm - the DM object
920: Output Parameter:
921: . vec - the global vector
923: Level: beginner
925: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix()
927: @*/
928: PetscErrorCode DMCreateGlobalVector(DM dm,Vec *vec)
929: {
935: if (!dm->ops->createglobalvector) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMCreateGlobalVector",((PetscObject)dm)->type_name);
936: (*dm->ops->createglobalvector)(dm,vec);
937: return(0);
938: }
940: /*@
941: DMCreateLocalVector - Creates a local vector from a DM object
943: Not Collective
945: Input Parameter:
946: . dm - the DM object
948: Output Parameter:
949: . vec - the local vector
951: Level: beginner
953: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix()
955: @*/
956: PetscErrorCode DMCreateLocalVector(DM dm,Vec *vec)
957: {
963: if (!dm->ops->createlocalvector) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMCreateLocalVector",((PetscObject)dm)->type_name);
964: (*dm->ops->createlocalvector)(dm,vec);
965: return(0);
966: }
968: /*@
969: DMGetLocalToGlobalMapping - Accesses the local-to-global mapping in a DM.
971: Collective on dm
973: Input Parameter:
974: . dm - the DM that provides the mapping
976: Output Parameter:
977: . ltog - the mapping
979: Level: intermediate
981: Notes:
982: This mapping can then be used by VecSetLocalToGlobalMapping() or
983: MatSetLocalToGlobalMapping().
985: .seealso: DMCreateLocalVector()
986: @*/
987: PetscErrorCode DMGetLocalToGlobalMapping(DM dm,ISLocalToGlobalMapping *ltog)
988: {
989: PetscInt bs = -1, bsLocal[2], bsMinMax[2];
995: if (!dm->ltogmap) {
996: PetscSection section, sectionGlobal;
998: DMGetLocalSection(dm, §ion);
999: if (section) {
1000: const PetscInt *cdofs;
1001: PetscInt *ltog;
1002: PetscInt pStart, pEnd, n, p, k, l;
1004: DMGetGlobalSection(dm, §ionGlobal);
1005: PetscSectionGetChart(section, &pStart, &pEnd);
1006: PetscSectionGetStorageSize(section, &n);
1007: PetscMalloc1(n, <og); /* We want the local+overlap size */
1008: for (p = pStart, l = 0; p < pEnd; ++p) {
1009: PetscInt bdof, cdof, dof, off, c, cind = 0;
1011: /* Should probably use constrained dofs */
1012: PetscSectionGetDof(section, p, &dof);
1013: PetscSectionGetConstraintDof(section, p, &cdof);
1014: PetscSectionGetConstraintIndices(section, p, &cdofs);
1015: PetscSectionGetOffset(sectionGlobal, p, &off);
1016: /* If you have dofs, and constraints, and they are unequal, we set the blocksize to 1 */
1017: bdof = cdof && (dof-cdof) ? 1 : dof;
1018: if (dof) {
1019: if (bs < 0) {bs = bdof;}
1020: else if (bs != bdof) {bs = 1;}
1021: }
1022: for (c = 0; c < dof; ++c, ++l) {
1023: if ((cind < cdof) && (c == cdofs[cind])) ltog[l] = off < 0 ? off-c : off+c;
1024: else ltog[l] = (off < 0 ? -(off+1) : off) + c;
1025: }
1026: }
1027: /* Must have same blocksize on all procs (some might have no points) */
1028: bsLocal[0] = bs < 0 ? PETSC_MAX_INT : bs; bsLocal[1] = bs;
1029: PetscGlobalMinMaxInt(PetscObjectComm((PetscObject) dm), bsLocal, bsMinMax);
1030: if (bsMinMax[0] != bsMinMax[1]) {bs = 1;}
1031: else {bs = bsMinMax[0];}
1032: bs = bs < 0 ? 1 : bs;
1033: /* Must reduce indices by blocksize */
1034: if (bs > 1) {
1035: for (l = 0, k = 0; l < n; l += bs, ++k) ltog[k] = ltog[l]/bs;
1036: n /= bs;
1037: }
1038: ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)dm), bs, n, ltog, PETSC_OWN_POINTER, &dm->ltogmap);
1039: PetscLogObjectParent((PetscObject)dm, (PetscObject)dm->ltogmap);
1040: } else {
1041: if (!dm->ops->getlocaltoglobalmapping) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMGetLocalToGlobalMapping",((PetscObject)dm)->type_name);
1042: (*dm->ops->getlocaltoglobalmapping)(dm);
1043: }
1044: }
1045: *ltog = dm->ltogmap;
1046: return(0);
1047: }
1049: /*@
1050: DMGetBlockSize - Gets the inherent block size associated with a DM
1052: Not Collective
1054: Input Parameter:
1055: . dm - the DM with block structure
1057: Output Parameter:
1058: . bs - the block size, 1 implies no exploitable block structure
1060: Level: intermediate
1062: .seealso: ISCreateBlock(), VecSetBlockSize(), MatSetBlockSize(), DMGetLocalToGlobalMapping()
1063: @*/
1064: PetscErrorCode DMGetBlockSize(DM dm,PetscInt *bs)
1065: {
1069: if (dm->bs < 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"DM does not have enough information to provide a block size yet");
1070: *bs = dm->bs;
1071: return(0);
1072: }
1074: /*@
1075: DMCreateInterpolation - Gets interpolation matrix between two DM objects
1077: Collective on dm1
1079: Input Parameter:
1080: + dm1 - the DM object
1081: - dm2 - the second, finer DM object
1083: Output Parameter:
1084: + mat - the interpolation
1085: - vec - the scaling (optional)
1087: Level: developer
1089: Notes:
1090: For DMDA objects this only works for "uniform refinement", that is the refined mesh was obtained DMRefine() or the coarse mesh was obtained by
1091: DMCoarsen(). The coordinates set into the DMDA are completely ignored in computing the interpolation.
1093: For DMDA objects you can use this interpolation (more precisely the interpolation from the DMGetCoordinateDM()) to interpolate the mesh coordinate vectors
1094: EXCEPT in the periodic case where it does not make sense since the coordinate vectors are not periodic.
1097: .seealso DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateColoring(), DMCreateMatrix(), DMRefine(), DMCoarsen(), DMCreateRestriction()
1099: @*/
1100: PetscErrorCode DMCreateInterpolation(DM dm1,DM dm2,Mat *mat,Vec *vec)
1101: {
1108: if (!dm1->ops->createinterpolation) SETERRQ1(PetscObjectComm((PetscObject)dm1),PETSC_ERR_SUP,"DM type %s does not implement DMCreateInterpolation",((PetscObject)dm1)->type_name);
1109: PetscLogEventBegin(DM_CreateInterpolation,dm1,dm2,0,0);
1110: (*dm1->ops->createinterpolation)(dm1,dm2,mat,vec);
1111: PetscLogEventEnd(DM_CreateInterpolation,dm1,dm2,0,0);
1112: return(0);
1113: }
1115: /*@
1116: DMCreateRestriction - Gets restriction matrix between two DM objects
1118: Collective on dm1
1120: Input Parameter:
1121: + dm1 - the DM object
1122: - dm2 - the second, finer DM object
1124: Output Parameter:
1125: . mat - the restriction
1128: Level: developer
1130: Notes:
1131: For DMDA objects this only works for "uniform refinement", that is the refined mesh was obtained DMRefine() or the coarse mesh was obtained by
1132: DMCoarsen(). The coordinates set into the DMDA are completely ignored in computing the interpolation.
1135: .seealso DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateColoring(), DMCreateMatrix(), DMRefine(), DMCoarsen(), DMCreateInterpolation()
1137: @*/
1138: PetscErrorCode DMCreateRestriction(DM dm1,DM dm2,Mat *mat)
1139: {
1146: if (!dm1->ops->createrestriction) SETERRQ1(PetscObjectComm((PetscObject)dm1),PETSC_ERR_SUP,"DM type %s does not implement DMCreateRestriction",((PetscObject)dm1)->type_name);
1147: PetscLogEventBegin(DM_CreateRestriction,dm1,dm2,0,0);
1148: (*dm1->ops->createrestriction)(dm1,dm2,mat);
1149: PetscLogEventEnd(DM_CreateRestriction,dm1,dm2,0,0);
1150: return(0);
1151: }
1153: /*@
1154: DMCreateInjection - Gets injection matrix between two DM objects
1156: Collective on dm1
1158: Input Parameter:
1159: + dm1 - the DM object
1160: - dm2 - the second, finer DM object
1162: Output Parameter:
1163: . mat - the injection
1165: Level: developer
1167: Notes:
1168: For DMDA objects this only works for "uniform refinement", that is the refined mesh was obtained DMRefine() or the coarse mesh was obtained by
1169: DMCoarsen(). The coordinates set into the DMDA are completely ignored in computing the injection.
1171: .seealso DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateColoring(), DMCreateMatrix(), DMCreateInterpolation()
1173: @*/
1174: PetscErrorCode DMCreateInjection(DM dm1,DM dm2,Mat *mat)
1175: {
1182: if (!dm1->ops->createinjection) SETERRQ1(PetscObjectComm((PetscObject)dm1),PETSC_ERR_SUP,"DM type %s does not implement DMCreateInjection",((PetscObject)dm1)->type_name);
1183: PetscLogEventBegin(DM_CreateInjection,dm1,dm2,0,0);
1184: (*dm1->ops->createinjection)(dm1,dm2,mat);
1185: PetscLogEventEnd(DM_CreateInjection,dm1,dm2,0,0);
1186: return(0);
1187: }
1189: /*@
1190: DMCreateMassMatrix - Gets mass matrix between two DM objects, M_ij = \int \phi_i \psi_j
1192: Collective on dm1
1194: Input Parameter:
1195: + dm1 - the DM object
1196: - dm2 - the second, finer DM object
1198: Output Parameter:
1199: . mat - the interpolation
1201: Level: developer
1203: .seealso DMCreateMatrix(), DMRefine(), DMCoarsen(), DMCreateRestriction(), DMCreateInterpolation(), DMCreateInjection()
1204: @*/
1205: PetscErrorCode DMCreateMassMatrix(DM dm1, DM dm2, Mat *mat)
1206: {
1213: if (!dm1->ops->createmassmatrix) SETERRQ1(PetscObjectComm((PetscObject)dm1),PETSC_ERR_SUP,"DM type %s does not implement DMCreateMassMatrix",((PetscObject)dm1)->type_name);
1214: (*dm1->ops->createmassmatrix)(dm1, dm2, mat);
1215: return(0);
1216: }
1218: /*@
1219: DMCreateColoring - Gets coloring for a DM
1221: Collective on dm
1223: Input Parameter:
1224: + dm - the DM object
1225: - ctype - IS_COLORING_LOCAL or IS_COLORING_GLOBAL
1227: Output Parameter:
1228: . coloring - the coloring
1230: Level: developer
1232: .seealso DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateMatrix(), DMSetMatType()
1234: @*/
1235: PetscErrorCode DMCreateColoring(DM dm,ISColoringType ctype,ISColoring *coloring)
1236: {
1242: if (!dm->ops->getcoloring) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMCreateColoring",((PetscObject)dm)->type_name);
1243: (*dm->ops->getcoloring)(dm,ctype,coloring);
1244: return(0);
1245: }
1247: /*@
1248: DMCreateMatrix - Gets empty Jacobian for a DM
1250: Collective on dm
1252: Input Parameter:
1253: . dm - the DM object
1255: Output Parameter:
1256: . mat - the empty Jacobian
1258: Level: beginner
1260: Notes:
1261: This properly preallocates the number of nonzeros in the sparse matrix so you
1262: do not need to do it yourself.
1264: By default it also sets the nonzero structure and puts in the zero entries. To prevent setting
1265: the nonzero pattern call DMSetMatrixPreallocateOnly()
1267: For structured grid problems, when you call MatView() on this matrix it is displayed using the global natural ordering, NOT in the ordering used
1268: internally by PETSc.
1270: For structured grid problems, in general it is easiest to use MatSetValuesStencil() or MatSetValuesLocal() to put values into the matrix because MatSetValues() requires
1271: the indices for the global numbering for DMDAs which is complicated.
1273: .seealso DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMSetMatType()
1275: @*/
1276: PetscErrorCode DMCreateMatrix(DM dm,Mat *mat)
1277: {
1283: if (!dm->ops->creatematrix) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMCreateMatrix",((PetscObject)dm)->type_name);
1284: MatInitializePackage();
1285: PetscLogEventBegin(DM_CreateMatrix,0,0,0,0);
1286: (*dm->ops->creatematrix)(dm,mat);
1287: /* Handle nullspace and near nullspace */
1288: if (dm->Nf) {
1289: MatNullSpace nullSpace;
1290: PetscInt Nf;
1292: DMGetNumFields(dm, &Nf);
1293: if (Nf == 1) {
1294: if (dm->nullspaceConstructors[0]) {
1295: (*dm->nullspaceConstructors[0])(dm, 0, &nullSpace);
1296: MatSetNullSpace(*mat, nullSpace);
1297: MatNullSpaceDestroy(&nullSpace);
1298: }
1299: if (dm->nearnullspaceConstructors[0]) {
1300: (*dm->nearnullspaceConstructors[0])(dm, 0, &nullSpace);
1301: MatSetNearNullSpace(*mat, nullSpace);
1302: MatNullSpaceDestroy(&nullSpace);
1303: }
1304: }
1305: }
1306: PetscLogEventEnd(DM_CreateMatrix,0,0,0,0);
1307: return(0);
1308: }
1310: /*@
1311: DMSetMatrixPreallocateOnly - When DMCreateMatrix() is called the matrix will be properly
1312: preallocated but the nonzero structure and zero values will not be set.
1314: Logically Collective on dm
1316: Input Parameter:
1317: + dm - the DM
1318: - only - PETSC_TRUE if only want preallocation
1320: Level: developer
1321: .seealso DMCreateMatrix(), DMSetMatrixStructureOnly()
1322: @*/
1323: PetscErrorCode DMSetMatrixPreallocateOnly(DM dm, PetscBool only)
1324: {
1327: dm->prealloc_only = only;
1328: return(0);
1329: }
1331: /*@
1332: DMSetMatrixStructureOnly - When DMCreateMatrix() is called, the matrix structure will be created
1333: but the array for values will not be allocated.
1335: Logically Collective on dm
1337: Input Parameter:
1338: + dm - the DM
1339: - only - PETSC_TRUE if only want matrix stucture
1341: Level: developer
1342: .seealso DMCreateMatrix(), DMSetMatrixPreallocateOnly()
1343: @*/
1344: PetscErrorCode DMSetMatrixStructureOnly(DM dm, PetscBool only)
1345: {
1348: dm->structure_only = only;
1349: return(0);
1350: }
1352: /*@C
1353: DMGetWorkArray - Gets a work array guaranteed to be at least the input size, restore with DMRestoreWorkArray()
1355: Not Collective
1357: Input Parameters:
1358: + dm - the DM object
1359: . count - The minium size
1360: - dtype - MPI data type, often MPIU_REAL, MPIU_SCALAR, MPIU_INT)
1362: Output Parameter:
1363: . array - the work array
1365: Level: developer
1367: .seealso DMDestroy(), DMCreate()
1368: @*/
1369: PetscErrorCode DMGetWorkArray(DM dm,PetscInt count,MPI_Datatype dtype,void *mem)
1370: {
1372: DMWorkLink link;
1373: PetscMPIInt dsize;
1378: if (dm->workin) {
1379: link = dm->workin;
1380: dm->workin = dm->workin->next;
1381: } else {
1382: PetscNewLog(dm,&link);
1383: }
1384: MPI_Type_size(dtype,&dsize);
1385: if (((size_t)dsize*count) > link->bytes) {
1386: PetscFree(link->mem);
1387: PetscMalloc(dsize*count,&link->mem);
1388: link->bytes = dsize*count;
1389: }
1390: link->next = dm->workout;
1391: dm->workout = link;
1392: *(void**)mem = link->mem;
1393: return(0);
1394: }
1396: /*@C
1397: DMRestoreWorkArray - Restores a work array guaranteed to be at least the input size, restore with DMRestoreWorkArray()
1399: Not Collective
1401: Input Parameters:
1402: + dm - the DM object
1403: . count - The minium size
1404: - dtype - MPI data type, often MPIU_REAL, MPIU_SCALAR, MPIU_INT
1406: Output Parameter:
1407: . array - the work array
1409: Level: developer
1411: Developer Notes:
1412: count and dtype are ignored, they are only needed for DMGetWorkArray()
1413: .seealso DMDestroy(), DMCreate()
1414: @*/
1415: PetscErrorCode DMRestoreWorkArray(DM dm,PetscInt count,MPI_Datatype dtype,void *mem)
1416: {
1417: DMWorkLink *p,link;
1422: for (p=&dm->workout; (link=*p); p=&link->next) {
1423: if (link->mem == *(void**)mem) {
1424: *p = link->next;
1425: link->next = dm->workin;
1426: dm->workin = link;
1427: *(void**)mem = NULL;
1428: return(0);
1429: }
1430: }
1431: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Array was not checked out");
1432: }
1434: PetscErrorCode DMSetNullSpaceConstructor(DM dm, PetscInt field, PetscErrorCode (*nullsp)(DM dm, PetscInt field, MatNullSpace *nullSpace))
1435: {
1438: if (field >= 10) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle %d >= 10 fields", field);
1439: dm->nullspaceConstructors[field] = nullsp;
1440: return(0);
1441: }
1443: PetscErrorCode DMGetNullSpaceConstructor(DM dm, PetscInt field, PetscErrorCode (**nullsp)(DM dm, PetscInt field, MatNullSpace *nullSpace))
1444: {
1448: if (field >= 10) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle %d >= 10 fields", field);
1449: *nullsp = dm->nullspaceConstructors[field];
1450: return(0);
1451: }
1453: PetscErrorCode DMSetNearNullSpaceConstructor(DM dm, PetscInt field, PetscErrorCode (*nullsp)(DM dm, PetscInt field, MatNullSpace *nullSpace))
1454: {
1457: if (field >= 10) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle %d >= 10 fields", field);
1458: dm->nearnullspaceConstructors[field] = nullsp;
1459: return(0);
1460: }
1462: PetscErrorCode DMGetNearNullSpaceConstructor(DM dm, PetscInt field, PetscErrorCode (**nullsp)(DM dm, PetscInt field, MatNullSpace *nullSpace))
1463: {
1467: if (field >= 10) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle %d >= 10 fields", field);
1468: *nullsp = dm->nearnullspaceConstructors[field];
1469: return(0);
1470: }
1472: /*@C
1473: DMCreateFieldIS - Creates a set of IS objects with the global indices of dofs for each field
1475: Not collective
1477: Input Parameter:
1478: . dm - the DM object
1480: Output Parameters:
1481: + numFields - The number of fields (or NULL if not requested)
1482: . fieldNames - The name for each field (or NULL if not requested)
1483: - fields - The global indices for each field (or NULL if not requested)
1485: Level: intermediate
1487: Notes:
1488: The user is responsible for freeing all requested arrays. In particular, every entry of names should be freed with
1489: PetscFree(), every entry of fields should be destroyed with ISDestroy(), and both arrays should be freed with
1490: PetscFree().
1492: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix()
1493: @*/
1494: PetscErrorCode DMCreateFieldIS(DM dm, PetscInt *numFields, char ***fieldNames, IS **fields)
1495: {
1496: PetscSection section, sectionGlobal;
1501: if (numFields) {
1503: *numFields = 0;
1504: }
1505: if (fieldNames) {
1507: *fieldNames = NULL;
1508: }
1509: if (fields) {
1511: *fields = NULL;
1512: }
1513: DMGetLocalSection(dm, §ion);
1514: if (section) {
1515: PetscInt *fieldSizes, *fieldNc, **fieldIndices;
1516: PetscInt nF, f, pStart, pEnd, p;
1518: DMGetGlobalSection(dm, §ionGlobal);
1519: PetscSectionGetNumFields(section, &nF);
1520: PetscMalloc3(nF,&fieldSizes,nF,&fieldNc,nF,&fieldIndices);
1521: PetscSectionGetChart(sectionGlobal, &pStart, &pEnd);
1522: for (f = 0; f < nF; ++f) {
1523: fieldSizes[f] = 0;
1524: PetscSectionGetFieldComponents(section, f, &fieldNc[f]);
1525: }
1526: for (p = pStart; p < pEnd; ++p) {
1527: PetscInt gdof;
1529: PetscSectionGetDof(sectionGlobal, p, &gdof);
1530: if (gdof > 0) {
1531: for (f = 0; f < nF; ++f) {
1532: PetscInt fdof, fcdof, fpdof;
1534: PetscSectionGetFieldDof(section, p, f, &fdof);
1535: PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);
1536: fpdof = fdof-fcdof;
1537: if (fpdof && fpdof != fieldNc[f]) {
1538: /* Layout does not admit a pointwise block size */
1539: fieldNc[f] = 1;
1540: }
1541: fieldSizes[f] += fpdof;
1542: }
1543: }
1544: }
1545: for (f = 0; f < nF; ++f) {
1546: PetscMalloc1(fieldSizes[f], &fieldIndices[f]);
1547: fieldSizes[f] = 0;
1548: }
1549: for (p = pStart; p < pEnd; ++p) {
1550: PetscInt gdof, goff;
1552: PetscSectionGetDof(sectionGlobal, p, &gdof);
1553: if (gdof > 0) {
1554: PetscSectionGetOffset(sectionGlobal, p, &goff);
1555: for (f = 0; f < nF; ++f) {
1556: PetscInt fdof, fcdof, fc;
1558: PetscSectionGetFieldDof(section, p, f, &fdof);
1559: PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);
1560: for (fc = 0; fc < fdof-fcdof; ++fc, ++fieldSizes[f]) {
1561: fieldIndices[f][fieldSizes[f]] = goff++;
1562: }
1563: }
1564: }
1565: }
1566: if (numFields) *numFields = nF;
1567: if (fieldNames) {
1568: PetscMalloc1(nF, fieldNames);
1569: for (f = 0; f < nF; ++f) {
1570: const char *fieldName;
1572: PetscSectionGetFieldName(section, f, &fieldName);
1573: PetscStrallocpy(fieldName, (char**) &(*fieldNames)[f]);
1574: }
1575: }
1576: if (fields) {
1577: PetscMalloc1(nF, fields);
1578: for (f = 0; f < nF; ++f) {
1579: PetscInt bs, in[2], out[2];
1581: ISCreateGeneral(PetscObjectComm((PetscObject)dm), fieldSizes[f], fieldIndices[f], PETSC_OWN_POINTER, &(*fields)[f]);
1582: in[0] = -fieldNc[f];
1583: in[1] = fieldNc[f];
1584: MPIU_Allreduce(in, out, 2, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)dm));
1585: bs = (-out[0] == out[1]) ? out[1] : 1;
1586: ISSetBlockSize((*fields)[f], bs);
1587: }
1588: }
1589: PetscFree3(fieldSizes,fieldNc,fieldIndices);
1590: } else if (dm->ops->createfieldis) {
1591: (*dm->ops->createfieldis)(dm, numFields, fieldNames, fields);
1592: }
1593: return(0);
1594: }
1597: /*@C
1598: DMCreateFieldDecomposition - Returns a list of IS objects defining a decomposition of a problem into subproblems
1599: corresponding to different fields: each IS contains the global indices of the dofs of the
1600: corresponding field. The optional list of DMs define the DM for each subproblem.
1601: Generalizes DMCreateFieldIS().
1603: Not collective
1605: Input Parameter:
1606: . dm - the DM object
1608: Output Parameters:
1609: + len - The number of subproblems in the field decomposition (or NULL if not requested)
1610: . namelist - The name for each field (or NULL if not requested)
1611: . islist - The global indices for each field (or NULL if not requested)
1612: - dmlist - The DMs for each field subproblem (or NULL, if not requested; if NULL is returned, no DMs are defined)
1614: Level: intermediate
1616: Notes:
1617: The user is responsible for freeing all requested arrays. In particular, every entry of names should be freed with
1618: PetscFree(), every entry of is should be destroyed with ISDestroy(), every entry of dm should be destroyed with DMDestroy(),
1619: and all of the arrays should be freed with PetscFree().
1621: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateFieldIS()
1622: @*/
1623: PetscErrorCode DMCreateFieldDecomposition(DM dm, PetscInt *len, char ***namelist, IS **islist, DM **dmlist)
1624: {
1629: if (len) {
1631: *len = 0;
1632: }
1633: if (namelist) {
1635: *namelist = 0;
1636: }
1637: if (islist) {
1639: *islist = 0;
1640: }
1641: if (dmlist) {
1643: *dmlist = 0;
1644: }
1645: /*
1646: Is it a good idea to apply the following check across all impls?
1647: Perhaps some impls can have a well-defined decomposition before DMSetUp?
1648: This, however, follows the general principle that accessors are not well-behaved until the object is set up.
1649: */
1650: if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE, "Decomposition defined only after DMSetUp");
1651: if (!dm->ops->createfielddecomposition) {
1652: PetscSection section;
1653: PetscInt numFields, f;
1655: DMGetLocalSection(dm, §ion);
1656: if (section) {PetscSectionGetNumFields(section, &numFields);}
1657: if (section && numFields && dm->ops->createsubdm) {
1658: if (len) *len = numFields;
1659: if (namelist) {PetscMalloc1(numFields,namelist);}
1660: if (islist) {PetscMalloc1(numFields,islist);}
1661: if (dmlist) {PetscMalloc1(numFields,dmlist);}
1662: for (f = 0; f < numFields; ++f) {
1663: const char *fieldName;
1665: DMCreateSubDM(dm, 1, &f, islist ? &(*islist)[f] : NULL, dmlist ? &(*dmlist)[f] : NULL);
1666: if (namelist) {
1667: PetscSectionGetFieldName(section, f, &fieldName);
1668: PetscStrallocpy(fieldName, (char**) &(*namelist)[f]);
1669: }
1670: }
1671: } else {
1672: DMCreateFieldIS(dm, len, namelist, islist);
1673: /* By default there are no DMs associated with subproblems. */
1674: if (dmlist) *dmlist = NULL;
1675: }
1676: } else {
1677: (*dm->ops->createfielddecomposition)(dm,len,namelist,islist,dmlist);
1678: }
1679: return(0);
1680: }
1682: /*@
1683: DMCreateSubDM - Returns an IS and DM encapsulating a subproblem defined by the fields passed in.
1684: The fields are defined by DMCreateFieldIS().
1686: Not collective
1688: Input Parameters:
1689: + dm - The DM object
1690: . numFields - The number of fields in this subproblem
1691: - fields - The field numbers of the selected fields
1693: Output Parameters:
1694: + is - The global indices for the subproblem
1695: - subdm - The DM for the subproblem
1697: Note: You need to call DMPlexSetMigrationSF() on the original DM if you want the Global-To-Natural map to be automatically constructed
1699: Level: intermediate
1701: .seealso DMPlexSetMigrationSF(), DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateFieldIS()
1702: @*/
1703: PetscErrorCode DMCreateSubDM(DM dm, PetscInt numFields, const PetscInt fields[], IS *is, DM *subdm)
1704: {
1712: if (!dm->ops->createsubdm) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMCreateSubDM",((PetscObject)dm)->type_name);
1713: (*dm->ops->createsubdm)(dm, numFields, fields, is, subdm);
1714: return(0);
1715: }
1717: /*@C
1718: DMCreateSuperDM - Returns an arrays of ISes and DM encapsulating a superproblem defined by the DMs passed in.
1720: Not collective
1722: Input Parameter:
1723: + dms - The DM objects
1724: - len - The number of DMs
1726: Output Parameters:
1727: + is - The global indices for the subproblem, or NULL
1728: - superdm - The DM for the superproblem
1730: Note: You need to call DMPlexSetMigrationSF() on the original DM if you want the Global-To-Natural map to be automatically constructed
1732: Level: intermediate
1734: .seealso DMPlexSetMigrationSF(), DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateFieldIS()
1735: @*/
1736: PetscErrorCode DMCreateSuperDM(DM dms[], PetscInt len, IS **is, DM *superdm)
1737: {
1738: PetscInt i;
1746: if (len < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of DMs must be nonnegative: %D", len);
1747: if (len) {
1748: DM dm = dms[0];
1749: if (!dm->ops->createsuperdm) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMCreateSuperDM",((PetscObject)dm)->type_name);
1750: (*dm->ops->createsuperdm)(dms, len, is, superdm);
1751: }
1752: return(0);
1753: }
1756: /*@C
1757: DMCreateDomainDecomposition - Returns lists of IS objects defining a decomposition of a problem into subproblems
1758: corresponding to restrictions to pairs nested subdomains: each IS contains the global
1759: indices of the dofs of the corresponding subdomains. The inner subdomains conceptually
1760: define a nonoverlapping covering, while outer subdomains can overlap.
1761: The optional list of DMs define the DM for each subproblem.
1763: Not collective
1765: Input Parameter:
1766: . dm - the DM object
1768: Output Parameters:
1769: + len - The number of subproblems in the domain decomposition (or NULL if not requested)
1770: . namelist - The name for each subdomain (or NULL if not requested)
1771: . innerislist - The global indices for each inner subdomain (or NULL, if not requested)
1772: . outerislist - The global indices for each outer subdomain (or NULL, if not requested)
1773: - dmlist - The DMs for each subdomain subproblem (or NULL, if not requested; if NULL is returned, no DMs are defined)
1775: Level: intermediate
1777: Notes:
1778: The user is responsible for freeing all requested arrays. In particular, every entry of names should be freed with
1779: PetscFree(), every entry of is should be destroyed with ISDestroy(), every entry of dm should be destroyed with DMDestroy(),
1780: and all of the arrays should be freed with PetscFree().
1782: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateFieldDecomposition()
1783: @*/
1784: PetscErrorCode DMCreateDomainDecomposition(DM dm, PetscInt *len, char ***namelist, IS **innerislist, IS **outerislist, DM **dmlist)
1785: {
1786: PetscErrorCode ierr;
1787: DMSubDomainHookLink link;
1788: PetscInt i,l;
1797: /*
1798: Is it a good idea to apply the following check across all impls?
1799: Perhaps some impls can have a well-defined decomposition before DMSetUp?
1800: This, however, follows the general principle that accessors are not well-behaved until the object is set up.
1801: */
1802: if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE, "Decomposition defined only after DMSetUp");
1803: if (dm->ops->createdomaindecomposition) {
1804: (*dm->ops->createdomaindecomposition)(dm,&l,namelist,innerislist,outerislist,dmlist);
1805: /* copy subdomain hooks and context over to the subdomain DMs */
1806: if (dmlist && *dmlist) {
1807: for (i = 0; i < l; i++) {
1808: for (link=dm->subdomainhook; link; link=link->next) {
1809: if (link->ddhook) {(*link->ddhook)(dm,(*dmlist)[i],link->ctx);}
1810: }
1811: if (dm->ctx) (*dmlist)[i]->ctx = dm->ctx;
1812: }
1813: }
1814: if (len) *len = l;
1815: }
1816: return(0);
1817: }
1820: /*@C
1821: DMCreateDomainDecompositionScatters - Returns scatters to the subdomain vectors from the global vector
1823: Not collective
1825: Input Parameters:
1826: + dm - the DM object
1827: . n - the number of subdomain scatters
1828: - subdms - the local subdomains
1830: Output Parameters:
1831: + n - the number of scatters returned
1832: . iscat - scatter from global vector to nonoverlapping global vector entries on subdomain
1833: . oscat - scatter from global vector to overlapping global vector entries on subdomain
1834: - gscat - scatter from global vector to local vector on subdomain (fills in ghosts)
1836: Notes:
1837: This is an alternative to the iis and ois arguments in DMCreateDomainDecomposition that allow for the solution
1838: of general nonlinear problems with overlapping subdomain methods. While merely having index sets that enable subsets
1839: of the residual equations to be created is fine for linear problems, nonlinear problems require local assembly of
1840: solution and residual data.
1842: Level: developer
1844: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateFieldIS()
1845: @*/
1846: PetscErrorCode DMCreateDomainDecompositionScatters(DM dm,PetscInt n,DM *subdms,VecScatter **iscat,VecScatter **oscat,VecScatter **gscat)
1847: {
1853: if (!dm->ops->createddscatters) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMCreateDomainDecompositionScatters",((PetscObject)dm)->type_name);
1854: (*dm->ops->createddscatters)(dm,n,subdms,iscat,oscat,gscat);
1855: return(0);
1856: }
1858: /*@
1859: DMRefine - Refines a DM object
1861: Collective on dm
1863: Input Parameter:
1864: + dm - the DM object
1865: - comm - the communicator to contain the new DM object (or MPI_COMM_NULL)
1867: Output Parameter:
1868: . dmf - the refined DM, or NULL
1870: Note: If no refinement was done, the return value is NULL
1872: Level: developer
1874: .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation()
1875: @*/
1876: PetscErrorCode DMRefine(DM dm,MPI_Comm comm,DM *dmf)
1877: {
1878: PetscErrorCode ierr;
1879: DMRefineHookLink link;
1883: if (!dm->ops->refine) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMRefine",((PetscObject)dm)->type_name);
1884: PetscLogEventBegin(DM_Refine,dm,0,0,0);
1885: (*dm->ops->refine)(dm,comm,dmf);
1886: if (*dmf) {
1887: (*dmf)->ops->creatematrix = dm->ops->creatematrix;
1889: PetscObjectCopyFortranFunctionPointers((PetscObject)dm,(PetscObject)*dmf);
1891: (*dmf)->ctx = dm->ctx;
1892: (*dmf)->leveldown = dm->leveldown;
1893: (*dmf)->levelup = dm->levelup + 1;
1895: DMSetMatType(*dmf,dm->mattype);
1896: for (link=dm->refinehook; link; link=link->next) {
1897: if (link->refinehook) {
1898: (*link->refinehook)(dm,*dmf,link->ctx);
1899: }
1900: }
1901: }
1902: PetscLogEventEnd(DM_Refine,dm,0,0,0);
1903: return(0);
1904: }
1906: /*@C
1907: DMRefineHookAdd - adds a callback to be run when interpolating a nonlinear problem to a finer grid
1909: Logically Collective
1911: Input Arguments:
1912: + coarse - nonlinear solver context on which to run a hook when restricting to a coarser level
1913: . refinehook - function to run when setting up a coarser level
1914: . interphook - function to run to update data on finer levels (once per SNESSolve())
1915: - ctx - [optional] user-defined context for provide data for the hooks (may be NULL)
1917: Calling sequence of refinehook:
1918: $ refinehook(DM coarse,DM fine,void *ctx);
1920: + coarse - coarse level DM
1921: . fine - fine level DM to interpolate problem to
1922: - ctx - optional user-defined function context
1924: Calling sequence for interphook:
1925: $ interphook(DM coarse,Mat interp,DM fine,void *ctx)
1927: + coarse - coarse level DM
1928: . interp - matrix interpolating a coarse-level solution to the finer grid
1929: . fine - fine level DM to update
1930: - ctx - optional user-defined function context
1932: Level: advanced
1934: Notes:
1935: This function is only needed if auxiliary data needs to be passed to fine grids while grid sequencing
1937: If this function is called multiple times, the hooks will be run in the order they are added.
1939: This function is currently not available from Fortran.
1941: .seealso: DMCoarsenHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
1942: @*/
1943: PetscErrorCode DMRefineHookAdd(DM coarse,PetscErrorCode (*refinehook)(DM,DM,void*),PetscErrorCode (*interphook)(DM,Mat,DM,void*),void *ctx)
1944: {
1945: PetscErrorCode ierr;
1946: DMRefineHookLink link,*p;
1950: for (p=&coarse->refinehook; *p; p=&(*p)->next) { /* Scan to the end of the current list of hooks */
1951: if ((*p)->refinehook == refinehook && (*p)->interphook == interphook && (*p)->ctx == ctx) return(0);
1952: }
1953: PetscNew(&link);
1954: link->refinehook = refinehook;
1955: link->interphook = interphook;
1956: link->ctx = ctx;
1957: link->next = NULL;
1958: *p = link;
1959: return(0);
1960: }
1962: /*@C
1963: DMRefineHookRemove - remove a callback from the list of hooks to be run when interpolating a nonlinear problem to a finer grid
1965: Logically Collective
1967: Input Arguments:
1968: + coarse - nonlinear solver context on which to run a hook when restricting to a coarser level
1969: . refinehook - function to run when setting up a coarser level
1970: . interphook - function to run to update data on finer levels (once per SNESSolve())
1971: - ctx - [optional] user-defined context for provide data for the hooks (may be NULL)
1973: Level: advanced
1975: Notes:
1976: This function does nothing if the hook is not in the list.
1978: This function is currently not available from Fortran.
1980: .seealso: DMCoarsenHookRemove(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
1981: @*/
1982: PetscErrorCode DMRefineHookRemove(DM coarse,PetscErrorCode (*refinehook)(DM,DM,void*),PetscErrorCode (*interphook)(DM,Mat,DM,void*),void *ctx)
1983: {
1984: PetscErrorCode ierr;
1985: DMRefineHookLink link,*p;
1989: for (p=&coarse->refinehook; *p; p=&(*p)->next) { /* Search the list of current hooks */
1990: if ((*p)->refinehook == refinehook && (*p)->interphook == interphook && (*p)->ctx == ctx) {
1991: link = *p;
1992: *p = link->next;
1993: PetscFree(link);
1994: break;
1995: }
1996: }
1997: return(0);
1998: }
2000: /*@
2001: DMInterpolate - interpolates user-defined problem data to a finer DM by running hooks registered by DMRefineHookAdd()
2003: Collective if any hooks are
2005: Input Arguments:
2006: + coarse - coarser DM to use as a base
2007: . interp - interpolation matrix, apply using MatInterpolate()
2008: - fine - finer DM to update
2010: Level: developer
2012: .seealso: DMRefineHookAdd(), MatInterpolate()
2013: @*/
2014: PetscErrorCode DMInterpolate(DM coarse,Mat interp,DM fine)
2015: {
2016: PetscErrorCode ierr;
2017: DMRefineHookLink link;
2020: for (link=fine->refinehook; link; link=link->next) {
2021: if (link->interphook) {
2022: (*link->interphook)(coarse,interp,fine,link->ctx);
2023: }
2024: }
2025: return(0);
2026: }
2028: /*@
2029: DMGetRefineLevel - Get's the number of refinements that have generated this DM.
2031: Not Collective
2033: Input Parameter:
2034: . dm - the DM object
2036: Output Parameter:
2037: . level - number of refinements
2039: Level: developer
2041: .seealso DMCoarsen(), DMGetCoarsenLevel(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation()
2043: @*/
2044: PetscErrorCode DMGetRefineLevel(DM dm,PetscInt *level)
2045: {
2048: *level = dm->levelup;
2049: return(0);
2050: }
2052: /*@
2053: DMSetRefineLevel - Set's the number of refinements that have generated this DM.
2055: Not Collective
2057: Input Parameter:
2058: + dm - the DM object
2059: - level - number of refinements
2061: Level: advanced
2063: Notes:
2064: This value is used by PCMG to determine how many multigrid levels to use
2066: .seealso DMCoarsen(), DMGetCoarsenLevel(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation()
2068: @*/
2069: PetscErrorCode DMSetRefineLevel(DM dm,PetscInt level)
2070: {
2073: dm->levelup = level;
2074: return(0);
2075: }
2077: PetscErrorCode DMGetBasisTransformDM_Internal(DM dm, DM *tdm)
2078: {
2082: *tdm = dm->transformDM;
2083: return(0);
2084: }
2086: PetscErrorCode DMGetBasisTransformVec_Internal(DM dm, Vec *tv)
2087: {
2091: *tv = dm->transform;
2092: return(0);
2093: }
2095: /*@
2096: DMHasBasisTransform - Whether we employ a basis transformation from functions in global vectors to functions in local vectors
2098: Input Parameter:
2099: . dm - The DM
2101: Output Parameter:
2102: . flg - PETSC_TRUE if a basis transformation should be done
2104: Level: developer
2106: .seealso: DMPlexGlobalToLocalBasis(), DMPlexLocalToGlobalBasis()()
2107: @*/
2108: PetscErrorCode DMHasBasisTransform(DM dm, PetscBool *flg)
2109: {
2110: Vec tv;
2116: DMGetBasisTransformVec_Internal(dm, &tv);
2117: *flg = tv ? PETSC_TRUE : PETSC_FALSE;
2118: return(0);
2119: }
2121: PetscErrorCode DMConstructBasisTransform_Internal(DM dm)
2122: {
2123: PetscSection s, ts;
2124: PetscScalar *ta;
2125: PetscInt cdim, pStart, pEnd, p, Nf, f, Nc, dof;
2129: DMGetCoordinateDim(dm, &cdim);
2130: DMGetLocalSection(dm, &s);
2131: PetscSectionGetChart(s, &pStart, &pEnd);
2132: PetscSectionGetNumFields(s, &Nf);
2133: DMClone(dm, &dm->transformDM);
2134: DMGetLocalSection(dm->transformDM, &ts);
2135: PetscSectionSetNumFields(ts, Nf);
2136: PetscSectionSetChart(ts, pStart, pEnd);
2137: for (f = 0; f < Nf; ++f) {
2138: PetscSectionGetFieldComponents(s, f, &Nc);
2139: /* We could start to label fields by their transformation properties */
2140: if (Nc != cdim) continue;
2141: for (p = pStart; p < pEnd; ++p) {
2142: PetscSectionGetFieldDof(s, p, f, &dof);
2143: if (!dof) continue;
2144: PetscSectionSetFieldDof(ts, p, f, PetscSqr(cdim));
2145: PetscSectionAddDof(ts, p, PetscSqr(cdim));
2146: }
2147: }
2148: PetscSectionSetUp(ts);
2149: DMCreateLocalVector(dm->transformDM, &dm->transform);
2150: VecGetArray(dm->transform, &ta);
2151: for (p = pStart; p < pEnd; ++p) {
2152: for (f = 0; f < Nf; ++f) {
2153: PetscSectionGetFieldDof(ts, p, f, &dof);
2154: if (dof) {
2155: PetscReal x[3] = {0.0, 0.0, 0.0};
2156: PetscScalar *tva;
2157: const PetscScalar *A;
2159: /* TODO Get quadrature point for this dual basis vector for coordinate */
2160: (*dm->transformGetMatrix)(dm, x, PETSC_TRUE, &A, dm->transformCtx);
2161: DMPlexPointLocalFieldRef(dm->transformDM, p, f, ta, (void *) &tva);
2162: PetscArraycpy(tva, A, PetscSqr(cdim));
2163: }
2164: }
2165: }
2166: VecRestoreArray(dm->transform, &ta);
2167: return(0);
2168: }
2170: PetscErrorCode DMCopyTransform(DM dm, DM newdm)
2171: {
2177: newdm->transformCtx = dm->transformCtx;
2178: newdm->transformSetUp = dm->transformSetUp;
2179: newdm->transformDestroy = NULL;
2180: newdm->transformGetMatrix = dm->transformGetMatrix;
2181: if (newdm->transformSetUp) {DMConstructBasisTransform_Internal(newdm);}
2182: return(0);
2183: }
2185: /*@C
2186: DMGlobalToLocalHookAdd - adds a callback to be run when global to local is called
2188: Logically Collective
2190: Input Arguments:
2191: + dm - the DM
2192: . beginhook - function to run at the beginning of DMGlobalToLocalBegin()
2193: . endhook - function to run after DMGlobalToLocalEnd() has completed
2194: - ctx - [optional] user-defined context for provide data for the hooks (may be NULL)
2196: Calling sequence for beginhook:
2197: $ beginhook(DM fine,VecScatter out,VecScatter in,DM coarse,void *ctx)
2199: + dm - global DM
2200: . g - global vector
2201: . mode - mode
2202: . l - local vector
2203: - ctx - optional user-defined function context
2206: Calling sequence for endhook:
2207: $ endhook(DM fine,VecScatter out,VecScatter in,DM coarse,void *ctx)
2209: + global - global DM
2210: - ctx - optional user-defined function context
2212: Level: advanced
2214: .seealso: DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2215: @*/
2216: PetscErrorCode DMGlobalToLocalHookAdd(DM dm,PetscErrorCode (*beginhook)(DM,Vec,InsertMode,Vec,void*),PetscErrorCode (*endhook)(DM,Vec,InsertMode,Vec,void*),void *ctx)
2217: {
2218: PetscErrorCode ierr;
2219: DMGlobalToLocalHookLink link,*p;
2223: for (p=&dm->gtolhook; *p; p=&(*p)->next) {} /* Scan to the end of the current list of hooks */
2224: PetscNew(&link);
2225: link->beginhook = beginhook;
2226: link->endhook = endhook;
2227: link->ctx = ctx;
2228: link->next = NULL;
2229: *p = link;
2230: return(0);
2231: }
2233: static PetscErrorCode DMGlobalToLocalHook_Constraints(DM dm, Vec g, InsertMode mode, Vec l, void *ctx)
2234: {
2235: Mat cMat;
2236: Vec cVec;
2237: PetscSection section, cSec;
2238: PetscInt pStart, pEnd, p, dof;
2243: DMGetDefaultConstraints(dm,&cSec,&cMat);
2244: if (cMat && (mode == INSERT_VALUES || mode == INSERT_ALL_VALUES || mode == INSERT_BC_VALUES)) {
2245: PetscInt nRows;
2247: MatGetSize(cMat,&nRows,NULL);
2248: if (nRows <= 0) return(0);
2249: DMGetLocalSection(dm,§ion);
2250: MatCreateVecs(cMat,NULL,&cVec);
2251: MatMult(cMat,l,cVec);
2252: PetscSectionGetChart(cSec,&pStart,&pEnd);
2253: for (p = pStart; p < pEnd; p++) {
2254: PetscSectionGetDof(cSec,p,&dof);
2255: if (dof) {
2256: PetscScalar *vals;
2257: VecGetValuesSection(cVec,cSec,p,&vals);
2258: VecSetValuesSection(l,section,p,vals,INSERT_ALL_VALUES);
2259: }
2260: }
2261: VecDestroy(&cVec);
2262: }
2263: return(0);
2264: }
2266: /*@
2267: DMGlobalToLocal - update local vectors from global vector
2269: Neighbor-wise Collective on dm
2271: Input Parameters:
2272: + dm - the DM object
2273: . g - the global vector
2274: . mode - INSERT_VALUES or ADD_VALUES
2275: - l - the local vector
2277: Notes:
2278: The communication involved in this update can be overlapped with computation by using
2279: DMGlobalToLocalBegin() and DMGlobalToLocalEnd().
2281: Level: beginner
2283: .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMGlobalToLocalEnd(), DMLocalToGlobalBegin(), DMLocalToGlobal(), DMLocalToGlobalBegin(), DMLocalToGlobalEnd()
2285: @*/
2286: PetscErrorCode DMGlobalToLocal(DM dm,Vec g,InsertMode mode,Vec l)
2287: {
2291: DMGlobalToLocalBegin(dm,g,mode,l);
2292: DMGlobalToLocalEnd(dm,g,mode,l);
2293: return(0);
2294: }
2296: /*@
2297: DMGlobalToLocalBegin - Begins updating local vectors from global vector
2299: Neighbor-wise Collective on dm
2301: Input Parameters:
2302: + dm - the DM object
2303: . g - the global vector
2304: . mode - INSERT_VALUES or ADD_VALUES
2305: - l - the local vector
2307: Level: intermediate
2309: .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMGlobalToLocal(), DMGlobalToLocalEnd(), DMLocalToGlobalBegin(), DMLocalToGlobal(), DMLocalToGlobalBegin(), DMLocalToGlobalEnd()
2311: @*/
2312: PetscErrorCode DMGlobalToLocalBegin(DM dm,Vec g,InsertMode mode,Vec l)
2313: {
2314: PetscSF sf;
2315: PetscErrorCode ierr;
2316: DMGlobalToLocalHookLink link;
2320: for (link=dm->gtolhook; link; link=link->next) {
2321: if (link->beginhook) {
2322: (*link->beginhook)(dm,g,mode,l,link->ctx);
2323: }
2324: }
2325: DMGetSectionSF(dm, &sf);
2326: if (sf) {
2327: const PetscScalar *gArray;
2328: PetscScalar *lArray;
2330: if (mode == ADD_VALUES) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode);
2331: VecGetArray(l, &lArray);
2332: VecGetArrayRead(g, &gArray);
2333: PetscSFBcastBegin(sf, MPIU_SCALAR, gArray, lArray);
2334: VecRestoreArray(l, &lArray);
2335: VecRestoreArrayRead(g, &gArray);
2336: } else {
2337: if (!dm->ops->globaltolocalbegin) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Missing DMGlobalToLocalBegin() for type %s",((PetscObject)dm)->type_name);
2338: (*dm->ops->globaltolocalbegin)(dm,g,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),l);
2339: }
2340: return(0);
2341: }
2343: /*@
2344: DMGlobalToLocalEnd - Ends updating local vectors from global vector
2346: Neighbor-wise Collective on dm
2348: Input Parameters:
2349: + dm - the DM object
2350: . g - the global vector
2351: . mode - INSERT_VALUES or ADD_VALUES
2352: - l - the local vector
2354: Level: intermediate
2356: .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMGlobalToLocal(), DMLocalToGlobalBegin(), DMLocalToGlobal(), DMLocalToGlobalBegin(), DMLocalToGlobalEnd()
2358: @*/
2359: PetscErrorCode DMGlobalToLocalEnd(DM dm,Vec g,InsertMode mode,Vec l)
2360: {
2361: PetscSF sf;
2362: PetscErrorCode ierr;
2363: const PetscScalar *gArray;
2364: PetscScalar *lArray;
2365: PetscBool transform;
2366: DMGlobalToLocalHookLink link;
2370: DMGetSectionSF(dm, &sf);
2371: DMHasBasisTransform(dm, &transform);
2372: if (sf) {
2373: if (mode == ADD_VALUES) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode);
2375: VecGetArray(l, &lArray);
2376: VecGetArrayRead(g, &gArray);
2377: PetscSFBcastEnd(sf, MPIU_SCALAR, gArray, lArray);
2378: VecRestoreArray(l, &lArray);
2379: VecRestoreArrayRead(g, &gArray);
2380: if (transform) {DMPlexGlobalToLocalBasis(dm, l);}
2381: } else {
2382: if (!dm->ops->globaltolocalend) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Missing DMGlobalToLocalEnd() for type %s",((PetscObject)dm)->type_name);
2383: (*dm->ops->globaltolocalend)(dm,g,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),l);
2384: }
2385: DMGlobalToLocalHook_Constraints(dm,g,mode,l,NULL);
2386: for (link=dm->gtolhook; link; link=link->next) {
2387: if (link->endhook) {(*link->endhook)(dm,g,mode,l,link->ctx);}
2388: }
2389: return(0);
2390: }
2392: /*@C
2393: DMLocalToGlobalHookAdd - adds a callback to be run when a local to global is called
2395: Logically Collective
2397: Input Arguments:
2398: + dm - the DM
2399: . beginhook - function to run at the beginning of DMLocalToGlobalBegin()
2400: . endhook - function to run after DMLocalToGlobalEnd() has completed
2401: - ctx - [optional] user-defined context for provide data for the hooks (may be NULL)
2403: Calling sequence for beginhook:
2404: $ beginhook(DM fine,Vec l,InsertMode mode,Vec g,void *ctx)
2406: + dm - global DM
2407: . l - local vector
2408: . mode - mode
2409: . g - global vector
2410: - ctx - optional user-defined function context
2413: Calling sequence for endhook:
2414: $ endhook(DM fine,Vec l,InsertMode mode,Vec g,void *ctx)
2416: + global - global DM
2417: . l - local vector
2418: . mode - mode
2419: . g - global vector
2420: - ctx - optional user-defined function context
2422: Level: advanced
2424: .seealso: DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2425: @*/
2426: PetscErrorCode DMLocalToGlobalHookAdd(DM dm,PetscErrorCode (*beginhook)(DM,Vec,InsertMode,Vec,void*),PetscErrorCode (*endhook)(DM,Vec,InsertMode,Vec,void*),void *ctx)
2427: {
2428: PetscErrorCode ierr;
2429: DMLocalToGlobalHookLink link,*p;
2433: for (p=&dm->ltoghook; *p; p=&(*p)->next) {} /* Scan to the end of the current list of hooks */
2434: PetscNew(&link);
2435: link->beginhook = beginhook;
2436: link->endhook = endhook;
2437: link->ctx = ctx;
2438: link->next = NULL;
2439: *p = link;
2440: return(0);
2441: }
2443: static PetscErrorCode DMLocalToGlobalHook_Constraints(DM dm, Vec l, InsertMode mode, Vec g, void *ctx)
2444: {
2445: Mat cMat;
2446: Vec cVec;
2447: PetscSection section, cSec;
2448: PetscInt pStart, pEnd, p, dof;
2453: DMGetDefaultConstraints(dm,&cSec,&cMat);
2454: if (cMat && (mode == ADD_VALUES || mode == ADD_ALL_VALUES || mode == ADD_BC_VALUES)) {
2455: PetscInt nRows;
2457: MatGetSize(cMat,&nRows,NULL);
2458: if (nRows <= 0) return(0);
2459: DMGetLocalSection(dm,§ion);
2460: MatCreateVecs(cMat,NULL,&cVec);
2461: PetscSectionGetChart(cSec,&pStart,&pEnd);
2462: for (p = pStart; p < pEnd; p++) {
2463: PetscSectionGetDof(cSec,p,&dof);
2464: if (dof) {
2465: PetscInt d;
2466: PetscScalar *vals;
2467: VecGetValuesSection(l,section,p,&vals);
2468: VecSetValuesSection(cVec,cSec,p,vals,mode);
2469: /* for this to be the true transpose, we have to zero the values that
2470: * we just extracted */
2471: for (d = 0; d < dof; d++) {
2472: vals[d] = 0.;
2473: }
2474: }
2475: }
2476: MatMultTransposeAdd(cMat,cVec,l,l);
2477: VecDestroy(&cVec);
2478: }
2479: return(0);
2480: }
2481: /*@
2482: DMLocalToGlobal - updates global vectors from local vectors
2484: Neighbor-wise Collective on dm
2486: Input Parameters:
2487: + dm - the DM object
2488: . l - the local vector
2489: . 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.
2490: - g - the global vector
2492: Notes:
2493: The communication involved in this update can be overlapped with computation by using
2494: DMLocalToGlobalBegin() and DMLocalToGlobalEnd().
2496: In the ADD_VALUES case you normally would zero the receiving vector before beginning this operation.
2497: INSERT_VALUES is not supported for DMDA; in that case simply compute the values directly into a global vector instead of a local one.
2499: Level: beginner
2501: .seealso DMLocalToGlobalBegin(), DMLocalToGlobalEnd(), DMCoarsen(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMGlobalToLocal(), DMGlobalToLocalEnd(), DMGlobalToLocalBegin()
2503: @*/
2504: PetscErrorCode DMLocalToGlobal(DM dm,Vec l,InsertMode mode,Vec g)
2505: {
2509: DMLocalToGlobalBegin(dm,l,mode,g);
2510: DMLocalToGlobalEnd(dm,l,mode,g);
2511: return(0);
2512: }
2514: /*@
2515: DMLocalToGlobalBegin - begins updating global vectors from local vectors
2517: Neighbor-wise Collective on dm
2519: Input Parameters:
2520: + dm - the DM object
2521: . l - the local vector
2522: . 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.
2523: - g - the global vector
2525: Notes:
2526: In the ADD_VALUES case you normally would zero the receiving vector before beginning this operation.
2527: INSERT_VALUES is not supported for DMDA, in that case simply compute the values directly into a global vector instead of a local one.
2529: Level: intermediate
2531: .seealso DMLocalToGlobal(), DMLocalToGlobalEnd(), DMCoarsen(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMGlobalToLocal(), DMGlobalToLocalEnd(), DMGlobalToLocalBegin()
2533: @*/
2534: PetscErrorCode DMLocalToGlobalBegin(DM dm,Vec l,InsertMode mode,Vec g)
2535: {
2536: PetscSF sf;
2537: PetscSection s, gs;
2538: DMLocalToGlobalHookLink link;
2539: Vec tmpl;
2540: const PetscScalar *lArray;
2541: PetscScalar *gArray;
2542: PetscBool isInsert, transform;
2543: PetscErrorCode ierr;
2547: for (link=dm->ltoghook; link; link=link->next) {
2548: if (link->beginhook) {
2549: (*link->beginhook)(dm,l,mode,g,link->ctx);
2550: }
2551: }
2552: DMLocalToGlobalHook_Constraints(dm,l,mode,g,NULL);
2553: DMGetSectionSF(dm, &sf);
2554: DMGetLocalSection(dm, &s);
2555: switch (mode) {
2556: case INSERT_VALUES:
2557: case INSERT_ALL_VALUES:
2558: case INSERT_BC_VALUES:
2559: isInsert = PETSC_TRUE; break;
2560: case ADD_VALUES:
2561: case ADD_ALL_VALUES:
2562: case ADD_BC_VALUES:
2563: isInsert = PETSC_FALSE; break;
2564: default:
2565: SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode);
2566: }
2567: if ((sf && !isInsert) || (s && isInsert)) {
2568: DMHasBasisTransform(dm, &transform);
2569: if (transform) {
2570: DMGetNamedLocalVector(dm, "__petsc_dm_transform_local_copy", &tmpl);
2571: VecCopy(l, tmpl);
2572: DMPlexLocalToGlobalBasis(dm, tmpl);
2573: VecGetArrayRead(tmpl, &lArray);
2574: } else {
2575: VecGetArrayRead(l, &lArray);
2576: }
2577: VecGetArray(g, &gArray);
2578: if (sf && !isInsert) {
2579: PetscSFReduceBegin(sf, MPIU_SCALAR, lArray, gArray, MPIU_SUM);
2580: } else if (s && isInsert) {
2581: PetscInt gStart, pStart, pEnd, p;
2583: DMGetGlobalSection(dm, &gs);
2584: PetscSectionGetChart(s, &pStart, &pEnd);
2585: VecGetOwnershipRange(g, &gStart, NULL);
2586: for (p = pStart; p < pEnd; ++p) {
2587: PetscInt dof, gdof, cdof, gcdof, off, goff, d, e;
2589: PetscSectionGetDof(s, p, &dof);
2590: PetscSectionGetDof(gs, p, &gdof);
2591: PetscSectionGetConstraintDof(s, p, &cdof);
2592: PetscSectionGetConstraintDof(gs, p, &gcdof);
2593: PetscSectionGetOffset(s, p, &off);
2594: PetscSectionGetOffset(gs, p, &goff);
2595: /* Ignore off-process data and points with no global data */
2596: if (!gdof || goff < 0) continue;
2597: 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);
2598: /* If no constraints are enforced in the global vector */
2599: if (!gcdof) {
2600: for (d = 0; d < dof; ++d) gArray[goff-gStart+d] = lArray[off+d];
2601: /* If constraints are enforced in the global vector */
2602: } else if (cdof == gcdof) {
2603: const PetscInt *cdofs;
2604: PetscInt cind = 0;
2606: PetscSectionGetConstraintIndices(s, p, &cdofs);
2607: for (d = 0, e = 0; d < dof; ++d) {
2608: if ((cind < cdof) && (d == cdofs[cind])) {++cind; continue;}
2609: gArray[goff-gStart+e++] = lArray[off+d];
2610: }
2611: } 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);
2612: }
2613: }
2614: VecRestoreArray(g, &gArray);
2615: if (transform) {
2616: VecRestoreArrayRead(tmpl, &lArray);
2617: DMRestoreNamedLocalVector(dm, "__petsc_dm_transform_local_copy", &tmpl);
2618: } else {
2619: VecRestoreArrayRead(l, &lArray);
2620: }
2621: } else {
2622: if (!dm->ops->localtoglobalbegin) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Missing DMLocalToGlobalBegin() for type %s",((PetscObject)dm)->type_name);
2623: (*dm->ops->localtoglobalbegin)(dm,l,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),g);
2624: }
2625: return(0);
2626: }
2628: /*@
2629: DMLocalToGlobalEnd - updates global vectors from local vectors
2631: Neighbor-wise Collective on dm
2633: Input Parameters:
2634: + dm - the DM object
2635: . l - the local vector
2636: . mode - INSERT_VALUES or ADD_VALUES
2637: - g - the global vector
2639: Level: intermediate
2641: .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMGlobalToLocalEnd(), DMGlobalToLocalEnd()
2643: @*/
2644: PetscErrorCode DMLocalToGlobalEnd(DM dm,Vec l,InsertMode mode,Vec g)
2645: {
2646: PetscSF sf;
2647: PetscSection s;
2648: DMLocalToGlobalHookLink link;
2649: PetscBool isInsert, transform;
2650: PetscErrorCode ierr;
2654: DMGetSectionSF(dm, &sf);
2655: DMGetLocalSection(dm, &s);
2656: switch (mode) {
2657: case INSERT_VALUES:
2658: case INSERT_ALL_VALUES:
2659: isInsert = PETSC_TRUE; break;
2660: case ADD_VALUES:
2661: case ADD_ALL_VALUES:
2662: isInsert = PETSC_FALSE; break;
2663: default:
2664: SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode);
2665: }
2666: if (sf && !isInsert) {
2667: const PetscScalar *lArray;
2668: PetscScalar *gArray;
2669: Vec tmpl;
2671: DMHasBasisTransform(dm, &transform);
2672: if (transform) {
2673: DMGetNamedLocalVector(dm, "__petsc_dm_transform_local_copy", &tmpl);
2674: VecGetArrayRead(tmpl, &lArray);
2675: } else {
2676: VecGetArrayRead(l, &lArray);
2677: }
2678: VecGetArray(g, &gArray);
2679: PetscSFReduceEnd(sf, MPIU_SCALAR, lArray, gArray, MPIU_SUM);
2680: if (transform) {
2681: VecRestoreArrayRead(tmpl, &lArray);
2682: DMRestoreNamedLocalVector(dm, "__petsc_dm_transform_local_copy", &tmpl);
2683: } else {
2684: VecRestoreArrayRead(l, &lArray);
2685: }
2686: VecRestoreArray(g, &gArray);
2687: } else if (s && isInsert) {
2688: } else {
2689: if (!dm->ops->localtoglobalend) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Missing DMLocalToGlobalEnd() for type %s",((PetscObject)dm)->type_name);
2690: (*dm->ops->localtoglobalend)(dm,l,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),g);
2691: }
2692: for (link=dm->ltoghook; link; link=link->next) {
2693: if (link->endhook) {(*link->endhook)(dm,g,mode,l,link->ctx);}
2694: }
2695: return(0);
2696: }
2698: /*@
2699: DMLocalToLocalBegin - Maps from a local vector (including ghost points
2700: that contain irrelevant values) to another local vector where the ghost
2701: points in the second are set correctly. Must be followed by DMLocalToLocalEnd().
2703: Neighbor-wise Collective on dm
2705: Input Parameters:
2706: + dm - the DM object
2707: . g - the original local vector
2708: - mode - one of INSERT_VALUES or ADD_VALUES
2710: Output Parameter:
2711: . l - the local vector with correct ghost values
2713: Level: intermediate
2715: Notes:
2716: The local vectors used here need not be the same as those
2717: obtained from DMCreateLocalVector(), BUT they
2718: must have the same parallel data layout; they could, for example, be
2719: obtained with VecDuplicate() from the DM originating vectors.
2721: .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateLocalVector(), DMCreateGlobalVector(), DMCreateInterpolation(), DMLocalToLocalEnd(), DMGlobalToLocalEnd(), DMLocalToGlobalBegin()
2723: @*/
2724: PetscErrorCode DMLocalToLocalBegin(DM dm,Vec g,InsertMode mode,Vec l)
2725: {
2726: PetscErrorCode ierr;
2730: if (!dm->ops->localtolocalbegin) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"This DM does not support local to local maps");
2731: (*dm->ops->localtolocalbegin)(dm,g,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),l);
2732: return(0);
2733: }
2735: /*@
2736: DMLocalToLocalEnd - Maps from a local vector (including ghost points
2737: that contain irrelevant values) to another local vector where the ghost
2738: points in the second are set correctly. Must be preceded by DMLocalToLocalBegin().
2740: Neighbor-wise Collective on dm
2742: Input Parameters:
2743: + da - the DM object
2744: . g - the original local vector
2745: - mode - one of INSERT_VALUES or ADD_VALUES
2747: Output Parameter:
2748: . l - the local vector with correct ghost values
2750: Level: intermediate
2752: Notes:
2753: The local vectors used here need not be the same as those
2754: obtained from DMCreateLocalVector(), BUT they
2755: must have the same parallel data layout; they could, for example, be
2756: obtained with VecDuplicate() from the DM originating vectors.
2758: .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateLocalVector(), DMCreateGlobalVector(), DMCreateInterpolation(), DMLocalToLocalBegin(), DMGlobalToLocalEnd(), DMLocalToGlobalBegin()
2760: @*/
2761: PetscErrorCode DMLocalToLocalEnd(DM dm,Vec g,InsertMode mode,Vec l)
2762: {
2763: PetscErrorCode ierr;
2767: if (!dm->ops->localtolocalend) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"This DM does not support local to local maps");
2768: (*dm->ops->localtolocalend)(dm,g,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),l);
2769: return(0);
2770: }
2773: /*@
2774: DMCoarsen - Coarsens a DM object
2776: Collective on dm
2778: Input Parameter:
2779: + dm - the DM object
2780: - comm - the communicator to contain the new DM object (or MPI_COMM_NULL)
2782: Output Parameter:
2783: . dmc - the coarsened DM
2785: Level: developer
2787: .seealso DMRefine(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation()
2789: @*/
2790: PetscErrorCode DMCoarsen(DM dm, MPI_Comm comm, DM *dmc)
2791: {
2792: PetscErrorCode ierr;
2793: DMCoarsenHookLink link;
2797: if (!dm->ops->coarsen) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMCoarsen",((PetscObject)dm)->type_name);
2798: PetscLogEventBegin(DM_Coarsen,dm,0,0,0);
2799: (*dm->ops->coarsen)(dm, comm, dmc);
2800: if (*dmc) {
2801: DMSetCoarseDM(dm,*dmc);
2802: (*dmc)->ops->creatematrix = dm->ops->creatematrix;
2803: PetscObjectCopyFortranFunctionPointers((PetscObject)dm,(PetscObject)*dmc);
2804: (*dmc)->ctx = dm->ctx;
2805: (*dmc)->levelup = dm->levelup;
2806: (*dmc)->leveldown = dm->leveldown + 1;
2807: DMSetMatType(*dmc,dm->mattype);
2808: for (link=dm->coarsenhook; link; link=link->next) {
2809: if (link->coarsenhook) {(*link->coarsenhook)(dm,*dmc,link->ctx);}
2810: }
2811: }
2812: PetscLogEventEnd(DM_Coarsen,dm,0,0,0);
2813: if (!(*dmc)) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "NULL coarse mesh produced");
2814: return(0);
2815: }
2817: /*@C
2818: DMCoarsenHookAdd - adds a callback to be run when restricting a nonlinear problem to the coarse grid
2820: Logically Collective
2822: Input Arguments:
2823: + fine - nonlinear solver context on which to run a hook when restricting to a coarser level
2824: . coarsenhook - function to run when setting up a coarser level
2825: . restricthook - function to run to update data on coarser levels (once per SNESSolve())
2826: - ctx - [optional] user-defined context for provide data for the hooks (may be NULL)
2828: Calling sequence of coarsenhook:
2829: $ coarsenhook(DM fine,DM coarse,void *ctx);
2831: + fine - fine level DM
2832: . coarse - coarse level DM to restrict problem to
2833: - ctx - optional user-defined function context
2835: Calling sequence for restricthook:
2836: $ restricthook(DM fine,Mat mrestrict,Vec rscale,Mat inject,DM coarse,void *ctx)
2838: + fine - fine level DM
2839: . mrestrict - matrix restricting a fine-level solution to the coarse grid
2840: . rscale - scaling vector for restriction
2841: . inject - matrix restricting by injection
2842: . coarse - coarse level DM to update
2843: - ctx - optional user-defined function context
2845: Level: advanced
2847: Notes:
2848: This function is only needed if auxiliary data needs to be set up on coarse grids.
2850: If this function is called multiple times, the hooks will be run in the order they are added.
2852: In order to compose with nonlinear preconditioning without duplicating storage, the hook should be implemented to
2853: extract the finest level information from its context (instead of from the SNES).
2855: This function is currently not available from Fortran.
2857: .seealso: DMCoarsenHookRemove(), DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2858: @*/
2859: PetscErrorCode DMCoarsenHookAdd(DM fine,PetscErrorCode (*coarsenhook)(DM,DM,void*),PetscErrorCode (*restricthook)(DM,Mat,Vec,Mat,DM,void*),void *ctx)
2860: {
2861: PetscErrorCode ierr;
2862: DMCoarsenHookLink link,*p;
2866: for (p=&fine->coarsenhook; *p; p=&(*p)->next) { /* Scan to the end of the current list of hooks */
2867: if ((*p)->coarsenhook == coarsenhook && (*p)->restricthook == restricthook && (*p)->ctx == ctx) return(0);
2868: }
2869: PetscNew(&link);
2870: link->coarsenhook = coarsenhook;
2871: link->restricthook = restricthook;
2872: link->ctx = ctx;
2873: link->next = NULL;
2874: *p = link;
2875: return(0);
2876: }
2878: /*@C
2879: DMCoarsenHookRemove - remove a callback from the list of hooks to be run when restricting a nonlinear problem to the coarse grid
2881: Logically Collective
2883: Input Arguments:
2884: + fine - nonlinear solver context on which to run a hook when restricting to a coarser level
2885: . coarsenhook - function to run when setting up a coarser level
2886: . restricthook - function to run to update data on coarser levels (once per SNESSolve())
2887: - ctx - [optional] user-defined context for provide data for the hooks (may be NULL)
2889: Level: advanced
2891: Notes:
2892: This function does nothing if the hook is not in the list.
2894: This function is currently not available from Fortran.
2896: .seealso: DMCoarsenHookAdd(), DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2897: @*/
2898: PetscErrorCode DMCoarsenHookRemove(DM fine,PetscErrorCode (*coarsenhook)(DM,DM,void*),PetscErrorCode (*restricthook)(DM,Mat,Vec,Mat,DM,void*),void *ctx)
2899: {
2900: PetscErrorCode ierr;
2901: DMCoarsenHookLink link,*p;
2905: for (p=&fine->coarsenhook; *p; p=&(*p)->next) { /* Search the list of current hooks */
2906: if ((*p)->coarsenhook == coarsenhook && (*p)->restricthook == restricthook && (*p)->ctx == ctx) {
2907: link = *p;
2908: *p = link->next;
2909: PetscFree(link);
2910: break;
2911: }
2912: }
2913: return(0);
2914: }
2917: /*@
2918: DMRestrict - restricts user-defined problem data to a coarser DM by running hooks registered by DMCoarsenHookAdd()
2920: Collective if any hooks are
2922: Input Arguments:
2923: + fine - finer DM to use as a base
2924: . restrct - restriction matrix, apply using MatRestrict()
2925: . rscale - scaling vector for restriction
2926: . inject - injection matrix, also use MatRestrict()
2927: - coarse - coarser DM to update
2929: Level: developer
2931: .seealso: DMCoarsenHookAdd(), MatRestrict()
2932: @*/
2933: PetscErrorCode DMRestrict(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse)
2934: {
2935: PetscErrorCode ierr;
2936: DMCoarsenHookLink link;
2939: for (link=fine->coarsenhook; link; link=link->next) {
2940: if (link->restricthook) {
2941: (*link->restricthook)(fine,restrct,rscale,inject,coarse,link->ctx);
2942: }
2943: }
2944: return(0);
2945: }
2947: /*@C
2948: DMSubDomainHookAdd - adds a callback to be run when restricting a problem to the coarse grid
2950: Logically Collective on global
2952: Input Arguments:
2953: + global - global DM
2954: . ddhook - function to run to pass data to the decomposition DM upon its creation
2955: . restricthook - function to run to update data on block solve (at the beginning of the block solve)
2956: - ctx - [optional] user-defined context for provide data for the hooks (may be NULL)
2959: Calling sequence for ddhook:
2960: $ ddhook(DM global,DM block,void *ctx)
2962: + global - global DM
2963: . block - block DM
2964: - ctx - optional user-defined function context
2966: Calling sequence for restricthook:
2967: $ restricthook(DM global,VecScatter out,VecScatter in,DM block,void *ctx)
2969: + global - global DM
2970: . out - scatter to the outer (with ghost and overlap points) block vector
2971: . in - scatter to block vector values only owned locally
2972: . block - block DM
2973: - ctx - optional user-defined function context
2975: Level: advanced
2977: Notes:
2978: This function is only needed if auxiliary data needs to be set up on subdomain DMs.
2980: If this function is called multiple times, the hooks will be run in the order they are added.
2982: In order to compose with nonlinear preconditioning without duplicating storage, the hook should be implemented to
2983: extract the global information from its context (instead of from the SNES).
2985: This function is currently not available from Fortran.
2987: .seealso: DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2988: @*/
2989: PetscErrorCode DMSubDomainHookAdd(DM global,PetscErrorCode (*ddhook)(DM,DM,void*),PetscErrorCode (*restricthook)(DM,VecScatter,VecScatter,DM,void*),void *ctx)
2990: {
2991: PetscErrorCode ierr;
2992: DMSubDomainHookLink link,*p;
2996: for (p=&global->subdomainhook; *p; p=&(*p)->next) { /* Scan to the end of the current list of hooks */
2997: if ((*p)->ddhook == ddhook && (*p)->restricthook == restricthook && (*p)->ctx == ctx) return(0);
2998: }
2999: PetscNew(&link);
3000: link->restricthook = restricthook;
3001: link->ddhook = ddhook;
3002: link->ctx = ctx;
3003: link->next = NULL;
3004: *p = link;
3005: return(0);
3006: }
3008: /*@C
3009: DMSubDomainHookRemove - remove a callback from the list to be run when restricting a problem to the coarse grid
3011: Logically Collective
3013: Input Arguments:
3014: + global - global DM
3015: . ddhook - function to run to pass data to the decomposition DM upon its creation
3016: . restricthook - function to run to update data on block solve (at the beginning of the block solve)
3017: - ctx - [optional] user-defined context for provide data for the hooks (may be NULL)
3019: Level: advanced
3021: Notes:
3023: This function is currently not available from Fortran.
3025: .seealso: DMSubDomainHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
3026: @*/
3027: PetscErrorCode DMSubDomainHookRemove(DM global,PetscErrorCode (*ddhook)(DM,DM,void*),PetscErrorCode (*restricthook)(DM,VecScatter,VecScatter,DM,void*),void *ctx)
3028: {
3029: PetscErrorCode ierr;
3030: DMSubDomainHookLink link,*p;
3034: for (p=&global->subdomainhook; *p; p=&(*p)->next) { /* Search the list of current hooks */
3035: if ((*p)->ddhook == ddhook && (*p)->restricthook == restricthook && (*p)->ctx == ctx) {
3036: link = *p;
3037: *p = link->next;
3038: PetscFree(link);
3039: break;
3040: }
3041: }
3042: return(0);
3043: }
3045: /*@
3046: DMSubDomainRestrict - restricts user-defined problem data to a block DM by running hooks registered by DMSubDomainHookAdd()
3048: Collective if any hooks are
3050: Input Arguments:
3051: + fine - finer DM to use as a base
3052: . oscatter - scatter from domain global vector filling subdomain global vector with overlap
3053: . gscatter - scatter from domain global vector filling subdomain local vector with ghosts
3054: - coarse - coarer DM to update
3056: Level: developer
3058: .seealso: DMCoarsenHookAdd(), MatRestrict()
3059: @*/
3060: PetscErrorCode DMSubDomainRestrict(DM global,VecScatter oscatter,VecScatter gscatter,DM subdm)
3061: {
3062: PetscErrorCode ierr;
3063: DMSubDomainHookLink link;
3066: for (link=global->subdomainhook; link; link=link->next) {
3067: if (link->restricthook) {
3068: (*link->restricthook)(global,oscatter,gscatter,subdm,link->ctx);
3069: }
3070: }
3071: return(0);
3072: }
3074: /*@
3075: DMGetCoarsenLevel - Get's the number of coarsenings that have generated this DM.
3077: Not Collective
3079: Input Parameter:
3080: . dm - the DM object
3082: Output Parameter:
3083: . level - number of coarsenings
3085: Level: developer
3087: .seealso DMCoarsen(), DMGetRefineLevel(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation()
3089: @*/
3090: PetscErrorCode DMGetCoarsenLevel(DM dm,PetscInt *level)
3091: {
3095: *level = dm->leveldown;
3096: return(0);
3097: }
3099: /*@
3100: DMSetCoarsenLevel - Sets the number of coarsenings that have generated this DM.
3102: Not Collective
3104: Input Parameters:
3105: + dm - the DM object
3106: - level - number of coarsenings
3108: Level: developer
3110: .seealso DMCoarsen(), DMGetCoarsenLevel(), DMGetRefineLevel(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation()
3111: @*/
3112: PetscErrorCode DMSetCoarsenLevel(DM dm,PetscInt level)
3113: {
3116: dm->leveldown = level;
3117: return(0);
3118: }
3122: /*@C
3123: DMRefineHierarchy - Refines a DM object, all levels at once
3125: Collective on dm
3127: Input Parameter:
3128: + dm - the DM object
3129: - nlevels - the number of levels of refinement
3131: Output Parameter:
3132: . dmf - the refined DM hierarchy
3134: Level: developer
3136: .seealso DMCoarsenHierarchy(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation()
3138: @*/
3139: PetscErrorCode DMRefineHierarchy(DM dm,PetscInt nlevels,DM dmf[])
3140: {
3145: if (nlevels < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"nlevels cannot be negative");
3146: if (nlevels == 0) return(0);
3148: if (dm->ops->refinehierarchy) {
3149: (*dm->ops->refinehierarchy)(dm,nlevels,dmf);
3150: } else if (dm->ops->refine) {
3151: PetscInt i;
3153: DMRefine(dm,PetscObjectComm((PetscObject)dm),&dmf[0]);
3154: for (i=1; i<nlevels; i++) {
3155: DMRefine(dmf[i-1],PetscObjectComm((PetscObject)dm),&dmf[i]);
3156: }
3157: } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"No RefineHierarchy for this DM yet");
3158: return(0);
3159: }
3161: /*@C
3162: DMCoarsenHierarchy - Coarsens a DM object, all levels at once
3164: Collective on dm
3166: Input Parameter:
3167: + dm - the DM object
3168: - nlevels - the number of levels of coarsening
3170: Output Parameter:
3171: . dmc - the coarsened DM hierarchy
3173: Level: developer
3175: .seealso DMRefineHierarchy(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation()
3177: @*/
3178: PetscErrorCode DMCoarsenHierarchy(DM dm, PetscInt nlevels, DM dmc[])
3179: {
3184: if (nlevels < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"nlevels cannot be negative");
3185: if (nlevels == 0) return(0);
3187: if (dm->ops->coarsenhierarchy) {
3188: (*dm->ops->coarsenhierarchy)(dm, nlevels, dmc);
3189: } else if (dm->ops->coarsen) {
3190: PetscInt i;
3192: DMCoarsen(dm,PetscObjectComm((PetscObject)dm),&dmc[0]);
3193: for (i=1; i<nlevels; i++) {
3194: DMCoarsen(dmc[i-1],PetscObjectComm((PetscObject)dm),&dmc[i]);
3195: }
3196: } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"No CoarsenHierarchy for this DM yet");
3197: return(0);
3198: }
3200: /*@C
3201: DMSetApplicationContextDestroy - Sets a user function that will be called to destroy the application context when the DM is destroyed
3203: Not Collective
3205: Input Parameters:
3206: + dm - the DM object
3207: - destroy - the destroy function
3209: Level: intermediate
3211: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext()
3213: @*/
3214: PetscErrorCode DMSetApplicationContextDestroy(DM dm,PetscErrorCode (*destroy)(void**))
3215: {
3218: dm->ctxdestroy = destroy;
3219: return(0);
3220: }
3222: /*@
3223: DMSetApplicationContext - Set a user context into a DM object
3225: Not Collective
3227: Input Parameters:
3228: + dm - the DM object
3229: - ctx - the user context
3231: Level: intermediate
3233: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext()
3235: @*/
3236: PetscErrorCode DMSetApplicationContext(DM dm,void *ctx)
3237: {
3240: dm->ctx = ctx;
3241: return(0);
3242: }
3244: /*@
3245: DMGetApplicationContext - Gets a user context from a DM object
3247: Not Collective
3249: Input Parameter:
3250: . dm - the DM object
3252: Output Parameter:
3253: . ctx - the user context
3255: Level: intermediate
3257: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext()
3259: @*/
3260: PetscErrorCode DMGetApplicationContext(DM dm,void *ctx)
3261: {
3264: *(void**)ctx = dm->ctx;
3265: return(0);
3266: }
3268: /*@C
3269: DMSetVariableBounds - sets a function to compute the lower and upper bound vectors for SNESVI.
3271: Logically Collective on dm
3273: Input Parameter:
3274: + dm - the DM object
3275: - f - the function that computes variable bounds used by SNESVI (use NULL to cancel a previous function that was set)
3277: Level: intermediate
3279: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext(),
3280: DMSetJacobian()
3282: @*/
3283: PetscErrorCode DMSetVariableBounds(DM dm,PetscErrorCode (*f)(DM,Vec,Vec))
3284: {
3287: dm->ops->computevariablebounds = f;
3288: return(0);
3289: }
3291: /*@
3292: DMHasVariableBounds - does the DM object have a variable bounds function?
3294: Not Collective
3296: Input Parameter:
3297: . dm - the DM object to destroy
3299: Output Parameter:
3300: . flg - PETSC_TRUE if the variable bounds function exists
3302: Level: developer
3304: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext()
3306: @*/
3307: PetscErrorCode DMHasVariableBounds(DM dm,PetscBool *flg)
3308: {
3312: *flg = (dm->ops->computevariablebounds) ? PETSC_TRUE : PETSC_FALSE;
3313: return(0);
3314: }
3316: /*@C
3317: DMComputeVariableBounds - compute variable bounds used by SNESVI.
3319: Logically Collective on dm
3321: Input Parameters:
3322: . dm - the DM object
3324: Output parameters:
3325: + xl - lower bound
3326: - xu - upper bound
3328: Level: advanced
3330: Notes:
3331: This is generally not called by users. It calls the function provided by the user with DMSetVariableBounds()
3333: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext()
3335: @*/
3336: PetscErrorCode DMComputeVariableBounds(DM dm, Vec xl, Vec xu)
3337: {
3344: if (!dm->ops->computevariablebounds) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMComputeVariableBounds",((PetscObject)dm)->type_name);
3345: (*dm->ops->computevariablebounds)(dm, xl,xu);
3346: return(0);
3347: }
3349: /*@
3350: DMHasColoring - does the DM object have a method of providing a coloring?
3352: Not Collective
3354: Input Parameter:
3355: . dm - the DM object
3357: Output Parameter:
3358: . flg - PETSC_TRUE if the DM has facilities for DMCreateColoring().
3360: Level: developer
3362: .seealso DMCreateColoring()
3364: @*/
3365: PetscErrorCode DMHasColoring(DM dm,PetscBool *flg)
3366: {
3370: *flg = (dm->ops->getcoloring) ? PETSC_TRUE : PETSC_FALSE;
3371: return(0);
3372: }
3374: /*@
3375: DMHasCreateRestriction - does the DM object have a method of providing a restriction?
3377: Not Collective
3379: Input Parameter:
3380: . dm - the DM object
3382: Output Parameter:
3383: . flg - PETSC_TRUE if the DM has facilities for DMCreateRestriction().
3385: Level: developer
3387: .seealso DMCreateRestriction()
3389: @*/
3390: PetscErrorCode DMHasCreateRestriction(DM dm,PetscBool *flg)
3391: {
3395: *flg = (dm->ops->createrestriction) ? PETSC_TRUE : PETSC_FALSE;
3396: return(0);
3397: }
3400: /*@
3401: DMHasCreateInjection - does the DM object have a method of providing an injection?
3403: Not Collective
3405: Input Parameter:
3406: . dm - the DM object
3408: Output Parameter:
3409: . flg - PETSC_TRUE if the DM has facilities for DMCreateInjection().
3411: Level: developer
3413: .seealso DMCreateInjection()
3415: @*/
3416: PetscErrorCode DMHasCreateInjection(DM dm,PetscBool *flg)
3417: {
3423: if (dm->ops->hascreateinjection) {
3424: (*dm->ops->hascreateinjection)(dm,flg);
3425: } else {
3426: *flg = (dm->ops->createinjection) ? PETSC_TRUE : PETSC_FALSE;
3427: }
3428: return(0);
3429: }
3432: /*@C
3433: DMSetVec - set the vector at which to compute residual, Jacobian and VI bounds, if the problem is nonlinear.
3435: Collective on dm
3437: Input Parameter:
3438: + dm - the DM object
3439: - x - location to compute residual and Jacobian, if NULL is passed to those routines; will be NULL for linear problems.
3441: Level: developer
3443: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext()
3445: @*/
3446: PetscErrorCode DMSetVec(DM dm,Vec x)
3447: {
3452: if (x) {
3453: if (!dm->x) {
3454: DMCreateGlobalVector(dm,&dm->x);
3455: }
3456: VecCopy(x,dm->x);
3457: } else if (dm->x) {
3458: VecDestroy(&dm->x);
3459: }
3460: return(0);
3461: }
3463: PetscFunctionList DMList = NULL;
3464: PetscBool DMRegisterAllCalled = PETSC_FALSE;
3466: /*@C
3467: DMSetType - Builds a DM, for a particular DM implementation.
3469: Collective on dm
3471: Input Parameters:
3472: + dm - The DM object
3473: - method - The name of the DM type
3475: Options Database Key:
3476: . -dm_type <type> - Sets the DM type; use -help for a list of available types
3478: Notes:
3479: See "petsc/include/petscdm.h" for available DM types (for instance, DM1D, DM2D, or DM3D).
3481: Level: intermediate
3483: .seealso: DMGetType(), DMCreate()
3484: @*/
3485: PetscErrorCode DMSetType(DM dm, DMType method)
3486: {
3487: PetscErrorCode (*r)(DM);
3488: PetscBool match;
3493: PetscObjectTypeCompare((PetscObject) dm, method, &match);
3494: if (match) return(0);
3496: DMRegisterAll();
3497: PetscFunctionListFind(DMList,method,&r);
3498: if (!r) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown DM type: %s", method);
3500: if (dm->ops->destroy) {
3501: (*dm->ops->destroy)(dm);
3502: }
3503: PetscMemzero(dm->ops,sizeof(*dm->ops));
3504: PetscObjectChangeTypeName((PetscObject)dm,method);
3505: (*r)(dm);
3506: return(0);
3507: }
3509: /*@C
3510: DMGetType - Gets the DM type name (as a string) from the DM.
3512: Not Collective
3514: Input Parameter:
3515: . dm - The DM
3517: Output Parameter:
3518: . type - The DM type name
3520: Level: intermediate
3522: .seealso: DMSetType(), DMCreate()
3523: @*/
3524: PetscErrorCode DMGetType(DM dm, DMType *type)
3525: {
3531: DMRegisterAll();
3532: *type = ((PetscObject)dm)->type_name;
3533: return(0);
3534: }
3536: /*@C
3537: DMConvert - Converts a DM to another DM, either of the same or different type.
3539: Collective on dm
3541: Input Parameters:
3542: + dm - the DM
3543: - newtype - new DM type (use "same" for the same type)
3545: Output Parameter:
3546: . M - pointer to new DM
3548: Notes:
3549: Cannot be used to convert a sequential DM to parallel or parallel to sequential,
3550: the MPI communicator of the generated DM is always the same as the communicator
3551: of the input DM.
3553: Level: intermediate
3555: .seealso: DMCreate()
3556: @*/
3557: PetscErrorCode DMConvert(DM dm, DMType newtype, DM *M)
3558: {
3559: DM B;
3560: char convname[256];
3561: PetscBool sametype/*, issame */;
3568: PetscObjectTypeCompare((PetscObject) dm, newtype, &sametype);
3569: /* PetscStrcmp(newtype, "same", &issame); */
3570: if (sametype) {
3571: *M = dm;
3572: PetscObjectReference((PetscObject) dm);
3573: return(0);
3574: } else {
3575: PetscErrorCode (*conv)(DM, DMType, DM*) = NULL;
3577: /*
3578: Order of precedence:
3579: 1) See if a specialized converter is known to the current DM.
3580: 2) See if a specialized converter is known to the desired DM class.
3581: 3) See if a good general converter is registered for the desired class
3582: 4) See if a good general converter is known for the current matrix.
3583: 5) Use a really basic converter.
3584: */
3586: /* 1) See if a specialized converter is known to the current DM and the desired class */
3587: PetscStrncpy(convname,"DMConvert_",sizeof(convname));
3588: PetscStrlcat(convname,((PetscObject) dm)->type_name,sizeof(convname));
3589: PetscStrlcat(convname,"_",sizeof(convname));
3590: PetscStrlcat(convname,newtype,sizeof(convname));
3591: PetscStrlcat(convname,"_C",sizeof(convname));
3592: PetscObjectQueryFunction((PetscObject)dm,convname,&conv);
3593: if (conv) goto foundconv;
3595: /* 2) See if a specialized converter is known to the desired DM class. */
3596: DMCreate(PetscObjectComm((PetscObject)dm), &B);
3597: DMSetType(B, newtype);
3598: PetscStrncpy(convname,"DMConvert_",sizeof(convname));
3599: PetscStrlcat(convname,((PetscObject) dm)->type_name,sizeof(convname));
3600: PetscStrlcat(convname,"_",sizeof(convname));
3601: PetscStrlcat(convname,newtype,sizeof(convname));
3602: PetscStrlcat(convname,"_C",sizeof(convname));
3603: PetscObjectQueryFunction((PetscObject)B,convname,&conv);
3604: if (conv) {
3605: DMDestroy(&B);
3606: goto foundconv;
3607: }
3609: #if 0
3610: /* 3) See if a good general converter is registered for the desired class */
3611: conv = B->ops->convertfrom;
3612: DMDestroy(&B);
3613: if (conv) goto foundconv;
3615: /* 4) See if a good general converter is known for the current matrix */
3616: if (dm->ops->convert) {
3617: conv = dm->ops->convert;
3618: }
3619: if (conv) goto foundconv;
3620: #endif
3622: /* 5) Use a really basic converter. */
3623: SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "No conversion possible between DM types %s and %s", ((PetscObject) dm)->type_name, newtype);
3625: foundconv:
3626: PetscLogEventBegin(DM_Convert,dm,0,0,0);
3627: (*conv)(dm,newtype,M);
3628: /* Things that are independent of DM type: We should consult DMClone() here */
3629: {
3630: PetscBool isper;
3631: const PetscReal *maxCell, *L;
3632: const DMBoundaryType *bd;
3633: DMGetPeriodicity(dm, &isper, &maxCell, &L, &bd);
3634: DMSetPeriodicity(*M, isper, maxCell, L, bd);
3635: }
3636: PetscLogEventEnd(DM_Convert,dm,0,0,0);
3637: }
3638: PetscObjectStateIncrease((PetscObject) *M);
3639: return(0);
3640: }
3642: /*--------------------------------------------------------------------------------------------------------------------*/
3644: /*@C
3645: DMRegister - Adds a new DM component implementation
3647: Not Collective
3649: Input Parameters:
3650: + name - The name of a new user-defined creation routine
3651: - create_func - The creation routine itself
3653: Notes:
3654: DMRegister() may be called multiple times to add several user-defined DMs
3657: Sample usage:
3658: .vb
3659: DMRegister("my_da", MyDMCreate);
3660: .ve
3662: Then, your DM type can be chosen with the procedural interface via
3663: .vb
3664: DMCreate(MPI_Comm, DM *);
3665: DMSetType(DM,"my_da");
3666: .ve
3667: or at runtime via the option
3668: .vb
3669: -da_type my_da
3670: .ve
3672: Level: advanced
3674: .seealso: DMRegisterAll(), DMRegisterDestroy()
3676: @*/
3677: PetscErrorCode DMRegister(const char sname[],PetscErrorCode (*function)(DM))
3678: {
3682: DMInitializePackage();
3683: PetscFunctionListAdd(&DMList,sname,function);
3684: return(0);
3685: }
3687: /*@C
3688: DMLoad - Loads a DM that has been stored in binary with DMView().
3690: Collective on viewer
3692: Input Parameters:
3693: + newdm - the newly loaded DM, this needs to have been created with DMCreate() or
3694: some related function before a call to DMLoad().
3695: - viewer - binary file viewer, obtained from PetscViewerBinaryOpen() or
3696: HDF5 file viewer, obtained from PetscViewerHDF5Open()
3698: Level: intermediate
3700: Notes:
3701: The type is determined by the data in the file, any type set into the DM before this call is ignored.
3703: Notes for advanced users:
3704: Most users should not need to know the details of the binary storage
3705: format, since DMLoad() and DMView() completely hide these details.
3706: But for anyone who's interested, the standard binary matrix storage
3707: format is
3708: .vb
3709: has not yet been determined
3710: .ve
3712: .seealso: PetscViewerBinaryOpen(), DMView(), MatLoad(), VecLoad()
3713: @*/
3714: PetscErrorCode DMLoad(DM newdm, PetscViewer viewer)
3715: {
3716: PetscBool isbinary, ishdf5;
3722: PetscViewerCheckReadable(viewer);
3723: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
3724: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);
3725: if (isbinary) {
3726: PetscInt classid;
3727: char type[256];
3729: PetscViewerBinaryRead(viewer,&classid,1,NULL,PETSC_INT);
3730: if (classid != DM_FILE_CLASSID) SETERRQ1(PetscObjectComm((PetscObject)newdm),PETSC_ERR_ARG_WRONG,"Not DM next in file, classid found %d",(int)classid);
3731: PetscViewerBinaryRead(viewer,type,256,NULL,PETSC_CHAR);
3732: DMSetType(newdm, type);
3733: if (newdm->ops->load) {(*newdm->ops->load)(newdm,viewer);}
3734: } else if (ishdf5) {
3735: if (newdm->ops->load) {(*newdm->ops->load)(newdm,viewer);}
3736: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid viewer; open viewer with PetscViewerBinaryOpen() or PetscViewerHDF5Open()");
3737: return(0);
3738: }
3740: /*@
3741: DMGetLocalBoundingBox - Returns the bounding box for the piece of the DM on this process.
3743: Not collective
3745: Input Parameter:
3746: . dm - the DM
3748: Output Parameters:
3749: + lmin - local minimum coordinates (length coord dim, optional)
3750: - lmax - local maximim coordinates (length coord dim, optional)
3752: Level: beginner
3754: Note: If the DM is a DMDA and has no coordinates, the index bounds are returned instead.
3757: .seealso: DMGetCoordinates(), DMGetCoordinatesLocal(), DMGetBoundingBox()
3758: @*/
3759: PetscErrorCode DMGetLocalBoundingBox(DM dm, PetscReal lmin[], PetscReal lmax[])
3760: {
3761: Vec coords = NULL;
3762: PetscReal min[3] = {PETSC_MAX_REAL, PETSC_MAX_REAL, PETSC_MAX_REAL};
3763: PetscReal max[3] = {PETSC_MIN_REAL, PETSC_MIN_REAL, PETSC_MIN_REAL};
3764: const PetscScalar *local_coords;
3765: PetscInt N, Ni;
3766: PetscInt cdim, i, j;
3767: PetscErrorCode ierr;
3771: DMGetCoordinateDim(dm, &cdim);
3772: DMGetCoordinates(dm, &coords);
3773: if (coords) {
3774: VecGetArrayRead(coords, &local_coords);
3775: VecGetLocalSize(coords, &N);
3776: Ni = N/cdim;
3777: for (i = 0; i < Ni; ++i) {
3778: for (j = 0; j < 3; ++j) {
3779: min[j] = j < cdim ? PetscMin(min[j], PetscRealPart(local_coords[i*cdim+j])) : 0;
3780: max[j] = j < cdim ? PetscMax(max[j], PetscRealPart(local_coords[i*cdim+j])) : 0;
3781: }
3782: }
3783: VecRestoreArrayRead(coords, &local_coords);
3784: } else {
3785: PetscBool isda;
3787: PetscObjectTypeCompare((PetscObject) dm, DMDA, &isda);
3788: if (isda) {DMGetLocalBoundingIndices_DMDA(dm, min, max);}
3789: }
3790: if (lmin) {PetscArraycpy(lmin, min, cdim);}
3791: if (lmax) {PetscArraycpy(lmax, max, cdim);}
3792: return(0);
3793: }
3795: /*@
3796: DMGetBoundingBox - Returns the global bounding box for the DM.
3798: Collective
3800: Input Parameter:
3801: . dm - the DM
3803: Output Parameters:
3804: + gmin - global minimum coordinates (length coord dim, optional)
3805: - gmax - global maximim coordinates (length coord dim, optional)
3807: Level: beginner
3809: .seealso: DMGetLocalBoundingBox(), DMGetCoordinates(), DMGetCoordinatesLocal()
3810: @*/
3811: PetscErrorCode DMGetBoundingBox(DM dm, PetscReal gmin[], PetscReal gmax[])
3812: {
3813: PetscReal lmin[3], lmax[3];
3814: PetscInt cdim;
3815: PetscMPIInt count;
3820: DMGetCoordinateDim(dm, &cdim);
3821: PetscMPIIntCast(cdim, &count);
3822: DMGetLocalBoundingBox(dm, lmin, lmax);
3823: if (gmin) {MPIU_Allreduce(lmin, gmin, count, MPIU_REAL, MPIU_MIN, PetscObjectComm((PetscObject) dm));}
3824: if (gmax) {MPIU_Allreduce(lmax, gmax, count, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject) dm));}
3825: return(0);
3826: }
3828: /******************************** FEM Support **********************************/
3830: PetscErrorCode DMPrintCellVector(PetscInt c, const char name[], PetscInt len, const PetscScalar x[])
3831: {
3832: PetscInt f;
3836: PetscPrintf(PETSC_COMM_SELF, "Cell %D Element %s\n", c, name);
3837: for (f = 0; f < len; ++f) {
3838: PetscPrintf(PETSC_COMM_SELF, " | %g |\n", (double)PetscRealPart(x[f]));
3839: }
3840: return(0);
3841: }
3843: PetscErrorCode DMPrintCellMatrix(PetscInt c, const char name[], PetscInt rows, PetscInt cols, const PetscScalar A[])
3844: {
3845: PetscInt f, g;
3849: PetscPrintf(PETSC_COMM_SELF, "Cell %D Element %s\n", c, name);
3850: for (f = 0; f < rows; ++f) {
3851: PetscPrintf(PETSC_COMM_SELF, " |");
3852: for (g = 0; g < cols; ++g) {
3853: PetscPrintf(PETSC_COMM_SELF, " % 9.5g", PetscRealPart(A[f*cols+g]));
3854: }
3855: PetscPrintf(PETSC_COMM_SELF, " |\n");
3856: }
3857: return(0);
3858: }
3860: PetscErrorCode DMPrintLocalVec(DM dm, const char name[], PetscReal tol, Vec X)
3861: {
3862: PetscInt localSize, bs;
3863: PetscMPIInt size;
3864: Vec x, xglob;
3865: const PetscScalar *xarray;
3866: PetscErrorCode ierr;
3869: MPI_Comm_size(PetscObjectComm((PetscObject) dm),&size);
3870: VecDuplicate(X, &x);
3871: VecCopy(X, x);
3872: VecChop(x, tol);
3873: PetscPrintf(PetscObjectComm((PetscObject) dm),"%s:\n",name);
3874: if (size > 1) {
3875: VecGetLocalSize(x,&localSize);
3876: VecGetArrayRead(x,&xarray);
3877: VecGetBlockSize(x,&bs);
3878: VecCreateMPIWithArray(PetscObjectComm((PetscObject) dm),bs,localSize,PETSC_DETERMINE,xarray,&xglob);
3879: } else {
3880: xglob = x;
3881: }
3882: VecView(xglob,PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject) dm)));
3883: if (size > 1) {
3884: VecDestroy(&xglob);
3885: VecRestoreArrayRead(x,&xarray);
3886: }
3887: VecDestroy(&x);
3888: return(0);
3889: }
3891: /*@
3892: DMGetSection - Get the PetscSection encoding the local data layout for the DM. This is equivalent to DMGetLocalSection(). Deprecated in v3.12
3894: Input Parameter:
3895: . dm - The DM
3897: Output Parameter:
3898: . section - The PetscSection
3900: Options Database Keys:
3901: . -dm_petscsection_view - View the Section created by the DM
3903: Level: advanced
3905: Notes:
3906: Use DMGetLocalSection() in new code.
3908: This gets a borrowed reference, so the user should not destroy this PetscSection.
3910: .seealso: DMGetLocalSection(), DMSetLocalSection(), DMGetGlobalSection()
3911: @*/
3912: PetscErrorCode DMGetSection(DM dm, PetscSection *section)
3913: {
3917: DMGetLocalSection(dm,section);
3918: return(0);
3919: }
3921: /*@
3922: DMGetLocalSection - Get the PetscSection encoding the local data layout for the DM.
3924: Input Parameter:
3925: . dm - The DM
3927: Output Parameter:
3928: . section - The PetscSection
3930: Options Database Keys:
3931: . -dm_petscsection_view - View the Section created by the DM
3933: Level: intermediate
3935: Note: This gets a borrowed reference, so the user should not destroy this PetscSection.
3937: .seealso: DMSetLocalSection(), DMGetGlobalSection()
3938: @*/
3939: PetscErrorCode DMGetLocalSection(DM dm, PetscSection *section)
3940: {
3946: if (!dm->localSection && dm->ops->createlocalsection) {
3947: PetscInt d;
3949: if (dm->setfromoptionscalled) for (d = 0; d < dm->Nds; ++d) {PetscDSSetFromOptions(dm->probs[d].ds);}
3950: (*dm->ops->createlocalsection)(dm);
3951: if (dm->localSection) {PetscObjectViewFromOptions((PetscObject) dm->localSection, NULL, "-dm_petscsection_view");}
3952: }
3953: *section = dm->localSection;
3954: return(0);
3955: }
3957: /*@
3958: DMSetSection - Set the PetscSection encoding the local data layout for the DM. This is equivalent to DMSetLocalSection(). Deprecated in v3.12
3960: Input Parameters:
3961: + dm - The DM
3962: - section - The PetscSection
3964: Level: advanced
3966: Notes:
3967: Use DMSetLocalSection() in new code.
3969: Any existing Section will be destroyed
3971: .seealso: DMSetLocalSection(), DMGetLocalSection(), DMSetGlobalSection()
3972: @*/
3973: PetscErrorCode DMSetSection(DM dm, PetscSection section)
3974: {
3978: DMSetLocalSection(dm,section);
3979: return(0);
3980: }
3982: /*@
3983: DMSetLocalSection - Set the PetscSection encoding the local data layout for the DM.
3985: Input Parameters:
3986: + dm - The DM
3987: - section - The PetscSection
3989: Level: intermediate
3991: Note: Any existing Section will be destroyed
3993: .seealso: DMGetLocalSection(), DMSetGlobalSection()
3994: @*/
3995: PetscErrorCode DMSetLocalSection(DM dm, PetscSection section)
3996: {
3997: PetscInt numFields = 0;
3998: PetscInt f;
4004: PetscObjectReference((PetscObject)section);
4005: PetscSectionDestroy(&dm->localSection);
4006: dm->localSection = section;
4007: if (section) {PetscSectionGetNumFields(dm->localSection, &numFields);}
4008: if (numFields) {
4009: DMSetNumFields(dm, numFields);
4010: for (f = 0; f < numFields; ++f) {
4011: PetscObject disc;
4012: const char *name;
4014: PetscSectionGetFieldName(dm->localSection, f, &name);
4015: DMGetField(dm, f, NULL, &disc);
4016: PetscObjectSetName(disc, name);
4017: }
4018: }
4019: /* The global section will be rebuilt in the next call to DMGetGlobalSection(). */
4020: PetscSectionDestroy(&dm->globalSection);
4021: return(0);
4022: }
4024: /*@
4025: DMGetDefaultConstraints - Get the PetscSection and Mat that specify the local constraint interpolation. See DMSetDefaultConstraints() for a description of the purpose of constraint interpolation.
4027: not collective
4029: Input Parameter:
4030: . dm - The DM
4032: Output Parameter:
4033: + 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.
4034: - 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.
4036: Level: advanced
4038: Note: This gets borrowed references, so the user should not destroy the PetscSection or the Mat.
4040: .seealso: DMSetDefaultConstraints()
4041: @*/
4042: PetscErrorCode DMGetDefaultConstraints(DM dm, PetscSection *section, Mat *mat)
4043: {
4048: if (!dm->defaultConstraintSection && !dm->defaultConstraintMat && dm->ops->createdefaultconstraints) {(*dm->ops->createdefaultconstraints)(dm);}
4049: if (section) {*section = dm->defaultConstraintSection;}
4050: if (mat) {*mat = dm->defaultConstraintMat;}
4051: return(0);
4052: }
4054: /*@
4055: DMSetDefaultConstraints - Set the PetscSection and Mat that specify the local constraint interpolation.
4057: 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().
4059: 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.
4061: collective on dm
4063: Input Parameters:
4064: + dm - The DM
4065: + 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).
4066: - 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).
4068: Level: advanced
4070: Note: This increments the references of the PetscSection and the Mat, so they user can destroy them
4072: .seealso: DMGetDefaultConstraints()
4073: @*/
4074: PetscErrorCode DMSetDefaultConstraints(DM dm, PetscSection section, Mat mat)
4075: {
4076: PetscMPIInt result;
4081: if (section) {
4083: MPI_Comm_compare(PETSC_COMM_SELF,PetscObjectComm((PetscObject)section),&result);
4084: if (result != MPI_CONGRUENT && result != MPI_IDENT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMECOMM,"constraint section must have local communicator");
4085: }
4086: if (mat) {
4088: MPI_Comm_compare(PETSC_COMM_SELF,PetscObjectComm((PetscObject)mat),&result);
4089: if (result != MPI_CONGRUENT && result != MPI_IDENT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMECOMM,"constraint matrix must have local communicator");
4090: }
4091: PetscObjectReference((PetscObject)section);
4092: PetscSectionDestroy(&dm->defaultConstraintSection);
4093: dm->defaultConstraintSection = section;
4094: PetscObjectReference((PetscObject)mat);
4095: MatDestroy(&dm->defaultConstraintMat);
4096: dm->defaultConstraintMat = mat;
4097: return(0);
4098: }
4100: #if defined(PETSC_USE_DEBUG)
4101: /*
4102: DMDefaultSectionCheckConsistency - Check the consistentcy of the global and local sections.
4104: Input Parameters:
4105: + dm - The DM
4106: . localSection - PetscSection describing the local data layout
4107: - globalSection - PetscSection describing the global data layout
4109: Level: intermediate
4111: .seealso: DMGetSectionSF(), DMSetSectionSF()
4112: */
4113: static PetscErrorCode DMDefaultSectionCheckConsistency_Internal(DM dm, PetscSection localSection, PetscSection globalSection)
4114: {
4115: MPI_Comm comm;
4116: PetscLayout layout;
4117: const PetscInt *ranges;
4118: PetscInt pStart, pEnd, p, nroots;
4119: PetscMPIInt size, rank;
4120: PetscBool valid = PETSC_TRUE, gvalid;
4121: PetscErrorCode ierr;
4124: PetscObjectGetComm((PetscObject)dm,&comm);
4126: MPI_Comm_size(comm, &size);
4127: MPI_Comm_rank(comm, &rank);
4128: PetscSectionGetChart(globalSection, &pStart, &pEnd);
4129: PetscSectionGetConstrainedStorageSize(globalSection, &nroots);
4130: PetscLayoutCreate(comm, &layout);
4131: PetscLayoutSetBlockSize(layout, 1);
4132: PetscLayoutSetLocalSize(layout, nroots);
4133: PetscLayoutSetUp(layout);
4134: PetscLayoutGetRanges(layout, &ranges);
4135: for (p = pStart; p < pEnd; ++p) {
4136: PetscInt dof, cdof, off, gdof, gcdof, goff, gsize, d;
4138: PetscSectionGetDof(localSection, p, &dof);
4139: PetscSectionGetOffset(localSection, p, &off);
4140: PetscSectionGetConstraintDof(localSection, p, &cdof);
4141: PetscSectionGetDof(globalSection, p, &gdof);
4142: PetscSectionGetConstraintDof(globalSection, p, &gcdof);
4143: PetscSectionGetOffset(globalSection, p, &goff);
4144: if (!gdof) continue; /* Censored point */
4145: 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;}
4146: 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;}
4147: if (gdof < 0) {
4148: gsize = gdof < 0 ? -(gdof+1)-gcdof : gdof-gcdof;
4149: for (d = 0; d < gsize; ++d) {
4150: PetscInt offset = -(goff+1) + d, r;
4152: PetscFindInt(offset,size+1,ranges,&r);
4153: if (r < 0) r = -(r+2);
4154: 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;}
4155: }
4156: }
4157: }
4158: PetscLayoutDestroy(&layout);
4159: PetscSynchronizedFlush(comm, NULL);
4160: MPIU_Allreduce(&valid, &gvalid, 1, MPIU_BOOL, MPI_LAND, comm);
4161: if (!gvalid) {
4162: DMView(dm, NULL);
4163: SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Inconsistent local and global sections");
4164: }
4165: return(0);
4166: }
4167: #endif
4169: /*@
4170: DMGetGlobalSection - Get the PetscSection encoding the global data layout for the DM.
4172: Collective on dm
4174: Input Parameter:
4175: . dm - The DM
4177: Output Parameter:
4178: . section - The PetscSection
4180: Level: intermediate
4182: Note: This gets a borrowed reference, so the user should not destroy this PetscSection.
4184: .seealso: DMSetLocalSection(), DMGetLocalSection()
4185: @*/
4186: PetscErrorCode DMGetGlobalSection(DM dm, PetscSection *section)
4187: {
4193: if (!dm->globalSection) {
4194: PetscSection s;
4196: DMGetLocalSection(dm, &s);
4197: if (!s) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONGSTATE, "DM must have a default PetscSection in order to create a global PetscSection");
4198: if (!dm->sf) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DM must have a point PetscSF in order to create a global PetscSection");
4199: PetscSectionCreateGlobalSection(s, dm->sf, PETSC_FALSE, PETSC_FALSE, &dm->globalSection);
4200: PetscLayoutDestroy(&dm->map);
4201: PetscSectionGetValueLayout(PetscObjectComm((PetscObject)dm), dm->globalSection, &dm->map);
4202: PetscSectionViewFromOptions(dm->globalSection, NULL, "-global_section_view");
4203: }
4204: *section = dm->globalSection;
4205: return(0);
4206: }
4208: /*@
4209: DMSetGlobalSection - Set the PetscSection encoding the global data layout for the DM.
4211: Input Parameters:
4212: + dm - The DM
4213: - section - The PetscSection, or NULL
4215: Level: intermediate
4217: Note: Any existing Section will be destroyed
4219: .seealso: DMGetGlobalSection(), DMSetLocalSection()
4220: @*/
4221: PetscErrorCode DMSetGlobalSection(DM dm, PetscSection section)
4222: {
4228: PetscObjectReference((PetscObject)section);
4229: PetscSectionDestroy(&dm->globalSection);
4230: dm->globalSection = section;
4231: #if defined(PETSC_USE_DEBUG)
4232: if (section) {DMDefaultSectionCheckConsistency_Internal(dm, dm->localSection, section);}
4233: #endif
4234: return(0);
4235: }
4237: /*@
4238: DMGetSectionSF - Get the PetscSF encoding the parallel dof overlap for the DM. If it has not been set,
4239: it is created from the default PetscSection layouts in the DM.
4241: Input Parameter:
4242: . dm - The DM
4244: Output Parameter:
4245: . sf - The PetscSF
4247: Level: intermediate
4249: Note: This gets a borrowed reference, so the user should not destroy this PetscSF.
4251: .seealso: DMSetSectionSF(), DMCreateSectionSF()
4252: @*/
4253: PetscErrorCode DMGetSectionSF(DM dm, PetscSF *sf)
4254: {
4255: PetscInt nroots;
4261: if (!dm->sectionSF) {
4262: PetscSFCreate(PetscObjectComm((PetscObject)dm),&dm->sectionSF);
4263: }
4264: PetscSFGetGraph(dm->sectionSF, &nroots, NULL, NULL, NULL);
4265: if (nroots < 0) {
4266: PetscSection section, gSection;
4268: DMGetLocalSection(dm, §ion);
4269: if (section) {
4270: DMGetGlobalSection(dm, &gSection);
4271: DMCreateSectionSF(dm, section, gSection);
4272: } else {
4273: *sf = NULL;
4274: return(0);
4275: }
4276: }
4277: *sf = dm->sectionSF;
4278: return(0);
4279: }
4281: /*@
4282: DMSetSectionSF - Set the PetscSF encoding the parallel dof overlap for the DM
4284: Input Parameters:
4285: + dm - The DM
4286: - sf - The PetscSF
4288: Level: intermediate
4290: Note: Any previous SF is destroyed
4292: .seealso: DMGetSectionSF(), DMCreateSectionSF()
4293: @*/
4294: PetscErrorCode DMSetSectionSF(DM dm, PetscSF sf)
4295: {
4301: PetscObjectReference((PetscObject) sf);
4302: PetscSFDestroy(&dm->sectionSF);
4303: dm->sectionSF = sf;
4304: return(0);
4305: }
4307: /*@C
4308: DMCreateSectionSF - Create the PetscSF encoding the parallel dof overlap for the DM based upon the PetscSections
4309: describing the data layout.
4311: Input Parameters:
4312: + dm - The DM
4313: . localSection - PetscSection describing the local data layout
4314: - globalSection - PetscSection describing the global data layout
4316: Notes: One usually uses DMGetSectionSF() to obtain the PetscSF
4318: Level: developer
4320: Developer Note: Since this routine has for arguments the two sections from the DM and puts the resulting PetscSF
4321: directly into the DM, perhaps this function should not take the local and global sections as
4322: input and should just obtain them from the DM?
4324: .seealso: DMGetSectionSF(), DMSetSectionSF(), DMGetLocalSection(), DMGetGlobalSection()
4325: @*/
4326: PetscErrorCode DMCreateSectionSF(DM dm, PetscSection localSection, PetscSection globalSection)
4327: {
4328: MPI_Comm comm;
4329: PetscLayout layout;
4330: const PetscInt *ranges;
4331: PetscInt *local;
4332: PetscSFNode *remote;
4333: PetscInt pStart, pEnd, p, nroots, nleaves = 0, l;
4334: PetscMPIInt size, rank;
4339: PetscObjectGetComm((PetscObject)dm,&comm);
4340: MPI_Comm_size(comm, &size);
4341: MPI_Comm_rank(comm, &rank);
4342: PetscSectionGetChart(globalSection, &pStart, &pEnd);
4343: PetscSectionGetConstrainedStorageSize(globalSection, &nroots);
4344: PetscLayoutCreate(comm, &layout);
4345: PetscLayoutSetBlockSize(layout, 1);
4346: PetscLayoutSetLocalSize(layout, nroots);
4347: PetscLayoutSetUp(layout);
4348: PetscLayoutGetRanges(layout, &ranges);
4349: for (p = pStart; p < pEnd; ++p) {
4350: PetscInt gdof, gcdof;
4352: PetscSectionGetDof(globalSection, p, &gdof);
4353: PetscSectionGetConstraintDof(globalSection, p, &gcdof);
4354: 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));
4355: nleaves += gdof < 0 ? -(gdof+1)-gcdof : gdof-gcdof;
4356: }
4357: PetscMalloc1(nleaves, &local);
4358: PetscMalloc1(nleaves, &remote);
4359: for (p = pStart, l = 0; p < pEnd; ++p) {
4360: const PetscInt *cind;
4361: PetscInt dof, cdof, off, gdof, gcdof, goff, gsize, d, c;
4363: PetscSectionGetDof(localSection, p, &dof);
4364: PetscSectionGetOffset(localSection, p, &off);
4365: PetscSectionGetConstraintDof(localSection, p, &cdof);
4366: PetscSectionGetConstraintIndices(localSection, p, &cind);
4367: PetscSectionGetDof(globalSection, p, &gdof);
4368: PetscSectionGetConstraintDof(globalSection, p, &gcdof);
4369: PetscSectionGetOffset(globalSection, p, &goff);
4370: if (!gdof) continue; /* Censored point */
4371: gsize = gdof < 0 ? -(gdof+1)-gcdof : gdof-gcdof;
4372: if (gsize != dof-cdof) {
4373: 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);
4374: cdof = 0; /* Ignore constraints */
4375: }
4376: for (d = 0, c = 0; d < dof; ++d) {
4377: if ((c < cdof) && (cind[c] == d)) {++c; continue;}
4378: local[l+d-c] = off+d;
4379: }
4380: if (gdof < 0) {
4381: for (d = 0; d < gsize; ++d, ++l) {
4382: PetscInt offset = -(goff+1) + d, r;
4384: PetscFindInt(offset,size+1,ranges,&r);
4385: if (r < 0) r = -(r+2);
4386: 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);
4387: remote[l].rank = r;
4388: remote[l].index = offset - ranges[r];
4389: }
4390: } else {
4391: for (d = 0; d < gsize; ++d, ++l) {
4392: remote[l].rank = rank;
4393: remote[l].index = goff+d - ranges[rank];
4394: }
4395: }
4396: }
4397: if (l != nleaves) SETERRQ2(comm, PETSC_ERR_PLIB, "Iteration error, l %d != nleaves %d", l, nleaves);
4398: PetscLayoutDestroy(&layout);
4399: PetscSFSetGraph(dm->sectionSF, nroots, nleaves, local, PETSC_OWN_POINTER, remote, PETSC_OWN_POINTER);
4400: return(0);
4401: }
4403: /*@
4404: DMGetPointSF - Get the PetscSF encoding the parallel section point overlap for the DM.
4406: Input Parameter:
4407: . dm - The DM
4409: Output Parameter:
4410: . sf - The PetscSF
4412: Level: intermediate
4414: Note: This gets a borrowed reference, so the user should not destroy this PetscSF.
4416: .seealso: DMSetPointSF(), DMGetSectionSF(), DMSetSectionSF(), DMCreateSectionSF()
4417: @*/
4418: PetscErrorCode DMGetPointSF(DM dm, PetscSF *sf)
4419: {
4423: *sf = dm->sf;
4424: return(0);
4425: }
4427: /*@
4428: DMSetPointSF - Set the PetscSF encoding the parallel section point overlap for the DM.
4430: Input Parameters:
4431: + dm - The DM
4432: - sf - The PetscSF
4434: Level: intermediate
4436: .seealso: DMGetPointSF(), DMGetSectionSF(), DMSetSectionSF(), DMCreateSectionSF()
4437: @*/
4438: PetscErrorCode DMSetPointSF(DM dm, PetscSF sf)
4439: {
4445: PetscObjectReference((PetscObject) sf);
4446: PetscSFDestroy(&dm->sf);
4447: dm->sf = sf;
4448: return(0);
4449: }
4451: static PetscErrorCode DMSetDefaultAdjacency_Private(DM dm, PetscInt f, PetscObject disc)
4452: {
4453: PetscClassId id;
4457: PetscObjectGetClassId(disc, &id);
4458: if (id == PETSCFE_CLASSID) {
4459: DMSetAdjacency(dm, f, PETSC_FALSE, PETSC_TRUE);
4460: } else if (id == PETSCFV_CLASSID) {
4461: DMSetAdjacency(dm, f, PETSC_TRUE, PETSC_FALSE);
4462: } else {
4463: DMSetAdjacency(dm, f, PETSC_FALSE, PETSC_TRUE);
4464: }
4465: return(0);
4466: }
4468: static PetscErrorCode DMFieldEnlarge_Static(DM dm, PetscInt NfNew)
4469: {
4470: RegionField *tmpr;
4471: PetscInt Nf = dm->Nf, f;
4475: if (Nf >= NfNew) return(0);
4476: PetscMalloc1(NfNew, &tmpr);
4477: for (f = 0; f < Nf; ++f) tmpr[f] = dm->fields[f];
4478: for (f = Nf; f < NfNew; ++f) {tmpr[f].disc = NULL; tmpr[f].label = NULL;}
4479: PetscFree(dm->fields);
4480: dm->Nf = NfNew;
4481: dm->fields = tmpr;
4482: return(0);
4483: }
4485: /*@
4486: DMClearFields - Remove all fields from the DM
4488: Logically collective on dm
4490: Input Parameter:
4491: . dm - The DM
4493: Level: intermediate
4495: .seealso: DMGetNumFields(), DMSetNumFields(), DMSetField()
4496: @*/
4497: PetscErrorCode DMClearFields(DM dm)
4498: {
4499: PetscInt f;
4504: for (f = 0; f < dm->Nf; ++f) {
4505: PetscObjectDestroy(&dm->fields[f].disc);
4506: DMLabelDestroy(&dm->fields[f].label);
4507: }
4508: PetscFree(dm->fields);
4509: dm->fields = NULL;
4510: dm->Nf = 0;
4511: return(0);
4512: }
4514: /*@
4515: DMGetNumFields - Get the number of fields in the DM
4517: Not collective
4519: Input Parameter:
4520: . dm - The DM
4522: Output Parameter:
4523: . Nf - The number of fields
4525: Level: intermediate
4527: .seealso: DMSetNumFields(), DMSetField()
4528: @*/
4529: PetscErrorCode DMGetNumFields(DM dm, PetscInt *numFields)
4530: {
4534: *numFields = dm->Nf;
4535: return(0);
4536: }
4538: /*@
4539: DMSetNumFields - Set the number of fields in the DM
4541: Logically collective on dm
4543: Input Parameters:
4544: + dm - The DM
4545: - Nf - The number of fields
4547: Level: intermediate
4549: .seealso: DMGetNumFields(), DMSetField()
4550: @*/
4551: PetscErrorCode DMSetNumFields(DM dm, PetscInt numFields)
4552: {
4553: PetscInt Nf, f;
4558: DMGetNumFields(dm, &Nf);
4559: for (f = Nf; f < numFields; ++f) {
4560: PetscContainer obj;
4562: PetscContainerCreate(PetscObjectComm((PetscObject) dm), &obj);
4563: DMAddField(dm, NULL, (PetscObject) obj);
4564: PetscContainerDestroy(&obj);
4565: }
4566: return(0);
4567: }
4569: /*@
4570: DMGetField - Return the discretization object for a given DM field
4572: Not collective
4574: Input Parameters:
4575: + dm - The DM
4576: - f - The field number
4578: Output Parameters:
4579: + label - The label indicating the support of the field, or NULL for the entire mesh
4580: - field - The discretization object
4582: Level: intermediate
4584: .seealso: DMAddField(), DMSetField()
4585: @*/
4586: PetscErrorCode DMGetField(DM dm, PetscInt f, DMLabel *label, PetscObject *field)
4587: {
4591: 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);
4592: if (label) *label = dm->fields[f].label;
4593: if (field) *field = dm->fields[f].disc;
4594: return(0);
4595: }
4597: /*@
4598: DMSetField - Set the discretization object for a given DM field
4600: Logically collective on dm
4602: Input Parameters:
4603: + dm - The DM
4604: . f - The field number
4605: . label - The label indicating the support of the field, or NULL for the entire mesh
4606: - field - The discretization object
4608: Level: intermediate
4610: .seealso: DMAddField(), DMGetField()
4611: @*/
4612: PetscErrorCode DMSetField(DM dm, PetscInt f, DMLabel label, PetscObject field)
4613: {
4620: if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
4621: DMFieldEnlarge_Static(dm, f+1);
4622: DMLabelDestroy(&dm->fields[f].label);
4623: PetscObjectDestroy(&dm->fields[f].disc);
4624: dm->fields[f].label = label;
4625: dm->fields[f].disc = field;
4626: PetscObjectReference((PetscObject) label);
4627: PetscObjectReference((PetscObject) field);
4628: DMSetDefaultAdjacency_Private(dm, f, field);
4629: DMClearDS(dm);
4630: return(0);
4631: }
4633: /*@
4634: DMAddField - Add the discretization object for the given DM field
4636: Logically collective on dm
4638: Input Parameters:
4639: + dm - The DM
4640: . label - The label indicating the support of the field, or NULL for the entire mesh
4641: - field - The discretization object
4643: Level: intermediate
4645: .seealso: DMSetField(), DMGetField()
4646: @*/
4647: PetscErrorCode DMAddField(DM dm, DMLabel label, PetscObject field)
4648: {
4649: PetscInt Nf = dm->Nf;
4656: DMFieldEnlarge_Static(dm, Nf+1);
4657: dm->fields[Nf].label = label;
4658: dm->fields[Nf].disc = field;
4659: PetscObjectReference((PetscObject) label);
4660: PetscObjectReference((PetscObject) field);
4661: DMSetDefaultAdjacency_Private(dm, Nf, field);
4662: DMClearDS(dm);
4663: return(0);
4664: }
4666: /*@
4667: DMCopyFields - Copy the discretizations for the DM into another DM
4669: Collective on dm
4671: Input Parameter:
4672: . dm - The DM
4674: Output Parameter:
4675: . newdm - The DM
4677: Level: advanced
4679: .seealso: DMGetField(), DMSetField(), DMAddField(), DMCopyDS(), DMGetDS(), DMGetCellDS()
4680: @*/
4681: PetscErrorCode DMCopyFields(DM dm, DM newdm)
4682: {
4683: PetscInt Nf, f;
4687: if (dm == newdm) return(0);
4688: DMGetNumFields(dm, &Nf);
4689: DMClearFields(newdm);
4690: for (f = 0; f < Nf; ++f) {
4691: DMLabel label;
4692: PetscObject field;
4693: PetscBool useCone, useClosure;
4695: DMGetField(dm, f, &label, &field);
4696: DMSetField(newdm, f, label, field);
4697: DMGetAdjacency(dm, f, &useCone, &useClosure);
4698: DMSetAdjacency(newdm, f, useCone, useClosure);
4699: }
4700: return(0);
4701: }
4703: /*@
4704: DMGetAdjacency - Returns the flags for determining variable influence
4706: Not collective
4708: Input Parameters:
4709: + dm - The DM object
4710: - f - The field number, or PETSC_DEFAULT for the default adjacency
4712: Output Parameter:
4713: + useCone - Flag for variable influence starting with the cone operation
4714: - useClosure - Flag for variable influence using transitive closure
4716: Notes:
4717: $ FEM: Two points p and q are adjacent if q \in closure(star(p)), useCone = PETSC_FALSE, useClosure = PETSC_TRUE
4718: $ FVM: Two points p and q are adjacent if q \in support(p+cone(p)), useCone = PETSC_TRUE, useClosure = PETSC_FALSE
4719: $ FVM++: Two points p and q are adjacent if q \in star(closure(p)), useCone = PETSC_TRUE, useClosure = PETSC_TRUE
4720: Further explanation can be found in the User's Manual Section on the Influence of Variables on One Another.
4722: Level: developer
4724: .seealso: DMSetAdjacency(), DMGetField(), DMSetField()
4725: @*/
4726: PetscErrorCode DMGetAdjacency(DM dm, PetscInt f, PetscBool *useCone, PetscBool *useClosure)
4727: {
4732: if (f < 0) {
4733: if (useCone) *useCone = dm->adjacency[0];
4734: if (useClosure) *useClosure = dm->adjacency[1];
4735: } else {
4736: PetscInt Nf;
4739: DMGetNumFields(dm, &Nf);
4740: if (f >= Nf) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, Nf);
4741: if (useCone) *useCone = dm->fields[f].adjacency[0];
4742: if (useClosure) *useClosure = dm->fields[f].adjacency[1];
4743: }
4744: return(0);
4745: }
4747: /*@
4748: DMSetAdjacency - Set the flags for determining variable influence
4750: Not collective
4752: Input Parameters:
4753: + dm - The DM object
4754: . f - The field number
4755: . useCone - Flag for variable influence starting with the cone operation
4756: - useClosure - Flag for variable influence using transitive closure
4758: Notes:
4759: $ FEM: Two points p and q are adjacent if q \in closure(star(p)), useCone = PETSC_FALSE, useClosure = PETSC_TRUE
4760: $ FVM: Two points p and q are adjacent if q \in support(p+cone(p)), useCone = PETSC_TRUE, useClosure = PETSC_FALSE
4761: $ FVM++: Two points p and q are adjacent if q \in star(closure(p)), useCone = PETSC_TRUE, useClosure = PETSC_TRUE
4762: Further explanation can be found in the User's Manual Section on the Influence of Variables on One Another.
4764: Level: developer
4766: .seealso: DMGetAdjacency(), DMGetField(), DMSetField()
4767: @*/
4768: PetscErrorCode DMSetAdjacency(DM dm, PetscInt f, PetscBool useCone, PetscBool useClosure)
4769: {
4772: if (f < 0) {
4773: dm->adjacency[0] = useCone;
4774: dm->adjacency[1] = useClosure;
4775: } else {
4776: PetscInt Nf;
4779: DMGetNumFields(dm, &Nf);
4780: if (f >= Nf) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, Nf);
4781: dm->fields[f].adjacency[0] = useCone;
4782: dm->fields[f].adjacency[1] = useClosure;
4783: }
4784: return(0);
4785: }
4787: /*@
4788: DMGetBasicAdjacency - Returns the flags for determining variable influence, using either the default or field 0 if it is defined
4790: Not collective
4792: Input Parameters:
4793: . dm - The DM object
4795: Output Parameter:
4796: + useCone - Flag for variable influence starting with the cone operation
4797: - useClosure - Flag for variable influence using transitive closure
4799: Notes:
4800: $ FEM: Two points p and q are adjacent if q \in closure(star(p)), useCone = PETSC_FALSE, useClosure = PETSC_TRUE
4801: $ FVM: Two points p and q are adjacent if q \in support(p+cone(p)), useCone = PETSC_TRUE, useClosure = PETSC_FALSE
4802: $ FVM++: Two points p and q are adjacent if q \in star(closure(p)), useCone = PETSC_TRUE, useClosure = PETSC_TRUE
4804: Level: developer
4806: .seealso: DMSetBasicAdjacency(), DMGetField(), DMSetField()
4807: @*/
4808: PetscErrorCode DMGetBasicAdjacency(DM dm, PetscBool *useCone, PetscBool *useClosure)
4809: {
4810: PetscInt Nf;
4817: DMGetNumFields(dm, &Nf);
4818: if (!Nf) {
4819: DMGetAdjacency(dm, PETSC_DEFAULT, useCone, useClosure);
4820: } else {
4821: DMGetAdjacency(dm, 0, useCone, useClosure);
4822: }
4823: return(0);
4824: }
4826: /*@
4827: DMSetBasicAdjacency - Set the flags for determining variable influence, using either the default or field 0 if it is defined
4829: Not collective
4831: Input Parameters:
4832: + dm - The DM object
4833: . useCone - Flag for variable influence starting with the cone operation
4834: - useClosure - Flag for variable influence using transitive closure
4836: Notes:
4837: $ FEM: Two points p and q are adjacent if q \in closure(star(p)), useCone = PETSC_FALSE, useClosure = PETSC_TRUE
4838: $ FVM: Two points p and q are adjacent if q \in support(p+cone(p)), useCone = PETSC_TRUE, useClosure = PETSC_FALSE
4839: $ FVM++: Two points p and q are adjacent if q \in star(closure(p)), useCone = PETSC_TRUE, useClosure = PETSC_TRUE
4841: Level: developer
4843: .seealso: DMGetBasicAdjacency(), DMGetField(), DMSetField()
4844: @*/
4845: PetscErrorCode DMSetBasicAdjacency(DM dm, PetscBool useCone, PetscBool useClosure)
4846: {
4847: PetscInt Nf;
4852: DMGetNumFields(dm, &Nf);
4853: if (!Nf) {
4854: DMSetAdjacency(dm, PETSC_DEFAULT, useCone, useClosure);
4855: } else {
4856: DMSetAdjacency(dm, 0, useCone, useClosure);
4857: }
4858: return(0);
4859: }
4861: static PetscErrorCode DMDSEnlarge_Static(DM dm, PetscInt NdsNew)
4862: {
4863: DMSpace *tmpd;
4864: PetscInt Nds = dm->Nds, s;
4868: if (Nds >= NdsNew) return(0);
4869: PetscMalloc1(NdsNew, &tmpd);
4870: for (s = 0; s < Nds; ++s) tmpd[s] = dm->probs[s];
4871: for (s = Nds; s < NdsNew; ++s) {tmpd[s].ds = NULL; tmpd[s].label = NULL; tmpd[s].fields = NULL;}
4872: PetscFree(dm->probs);
4873: dm->Nds = NdsNew;
4874: dm->probs = tmpd;
4875: return(0);
4876: }
4878: /*@
4879: DMGetNumDS - Get the number of discrete systems in the DM
4881: Not collective
4883: Input Parameter:
4884: . dm - The DM
4886: Output Parameter:
4887: . Nds - The number of PetscDS objects
4889: Level: intermediate
4891: .seealso: DMGetDS(), DMGetCellDS()
4892: @*/
4893: PetscErrorCode DMGetNumDS(DM dm, PetscInt *Nds)
4894: {
4898: *Nds = dm->Nds;
4899: return(0);
4900: }
4902: /*@
4903: DMClearDS - Remove all discrete systems from the DM
4905: Logically collective on dm
4907: Input Parameter:
4908: . dm - The DM
4910: Level: intermediate
4912: .seealso: DMGetNumDS(), DMGetDS(), DMSetField()
4913: @*/
4914: PetscErrorCode DMClearDS(DM dm)
4915: {
4916: PetscInt s;
4921: for (s = 0; s < dm->Nds; ++s) {
4922: PetscDSDestroy(&dm->probs[s].ds);
4923: DMLabelDestroy(&dm->probs[s].label);
4924: ISDestroy(&dm->probs[s].fields);
4925: }
4926: PetscFree(dm->probs);
4927: dm->probs = NULL;
4928: dm->Nds = 0;
4929: return(0);
4930: }
4932: /*@
4933: DMGetDS - Get the default PetscDS
4935: Not collective
4937: Input Parameter:
4938: . dm - The DM
4940: Output Parameter:
4941: . prob - The default PetscDS
4943: Level: intermediate
4945: .seealso: DMGetCellDS(), DMGetRegionDS()
4946: @*/
4947: PetscErrorCode DMGetDS(DM dm, PetscDS *prob)
4948: {
4954: if (dm->Nds <= 0) {
4955: PetscDS ds;
4957: PetscDSCreate(PetscObjectComm((PetscObject) dm), &ds);
4958: DMSetRegionDS(dm, NULL, NULL, ds);
4959: PetscDSDestroy(&ds);
4960: }
4961: *prob = dm->probs[0].ds;
4962: return(0);
4963: }
4965: /*@
4966: DMGetCellDS - Get the PetscDS defined on a given cell
4968: Not collective
4970: Input Parameters:
4971: + dm - The DM
4972: - point - Cell for the DS
4974: Output Parameter:
4975: . prob - The PetscDS defined on the given cell
4977: Level: developer
4979: .seealso: DMGetDS(), DMSetRegionDS()
4980: @*/
4981: PetscErrorCode DMGetCellDS(DM dm, PetscInt point, PetscDS *prob)
4982: {
4983: PetscDS probDef = NULL;
4984: PetscInt s;
4990: *prob = NULL;
4991: for (s = 0; s < dm->Nds; ++s) {
4992: PetscInt val;
4994: if (!dm->probs[s].label) {probDef = dm->probs[s].ds;}
4995: else {
4996: DMLabelGetValue(dm->probs[s].label, point, &val);
4997: if (val >= 0) {*prob = dm->probs[s].ds; break;}
4998: }
4999: }
5000: if (!*prob) *prob = probDef;
5001: return(0);
5002: }
5004: /*@
5005: DMGetRegionDS - Get the PetscDS for a given mesh region, defined by a DMLabel
5007: Not collective
5009: Input Parameters:
5010: + dm - The DM
5011: - label - The DMLabel defining the mesh region, or NULL for the entire mesh
5013: Output Parameters:
5014: + fields - The IS containing the DM field numbers for the fields in this DS, or NULL
5015: - prob - The PetscDS defined on the given region, or NULL
5017: Note: If the label is missing, this function returns an error
5019: Level: advanced
5021: .seealso: DMGetRegionNumDS(), DMSetRegionDS(), DMGetDS(), DMGetCellDS()
5022: @*/
5023: PetscErrorCode DMGetRegionDS(DM dm, DMLabel label, IS *fields, PetscDS *ds)
5024: {
5025: PetscInt Nds = dm->Nds, s;
5032: for (s = 0; s < Nds; ++s) {
5033: if (dm->probs[s].label == label) {
5034: if (fields) *fields = dm->probs[s].fields;
5035: if (ds) *ds = dm->probs[s].ds;
5036: return(0);
5037: }
5038: }
5039: return(0);
5040: }
5042: /*@
5043: DMGetRegionNumDS - Get the PetscDS for a given mesh region, defined by the region number
5045: Not collective
5047: Input Parameters:
5048: + dm - The DM
5049: - num - The region number, in [0, Nds)
5051: Output Parameters:
5052: + label - The region label, or NULL
5053: . fields - The IS containing the DM field numbers for the fields in this DS, or NULL
5054: - prob - The PetscDS defined on the given region, or NULL
5056: Level: advanced
5058: .seealso: DMGetRegionDS(), DMSetRegionDS(), DMGetDS(), DMGetCellDS()
5059: @*/
5060: PetscErrorCode DMGetRegionNumDS(DM dm, PetscInt num, DMLabel *label, IS *fields, PetscDS *ds)
5061: {
5062: PetscInt Nds;
5067: DMGetNumDS(dm, &Nds);
5068: if ((num < 0) || (num >= Nds)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Region number %D is not in [0, %D)", num, Nds);
5069: if (label) {
5071: *label = dm->probs[num].label;
5072: }
5073: if (fields) {
5075: *fields = dm->probs[num].fields;
5076: }
5077: if (ds) {
5079: *ds = dm->probs[num].ds;
5080: }
5081: return(0);
5082: }
5084: /*@
5085: DMSetRegionDS - Set the PetscDS for a given mesh region, defined by a DMLabel
5087: Collective on dm
5089: Input Parameters:
5090: + dm - The DM
5091: . label - The DMLabel defining the mesh region, or NULL for the entire mesh
5092: . fields - The IS containing the DM field numbers for the fields in this DS, or NULL for all fields
5093: - prob - The PetscDS defined on the given cell
5095: Note: If the label has a DS defined, it will be replaced. Otherwise, it will be added to the DM. If DS is replaced,
5096: the fields argument is ignored.
5098: Level: advanced
5100: .seealso: DMGetRegionDS(), DMGetDS(), DMGetCellDS()
5101: @*/
5102: PetscErrorCode DMSetRegionDS(DM dm, DMLabel label, IS fields, PetscDS ds)
5103: {
5104: PetscInt Nds = dm->Nds, s;
5111: for (s = 0; s < Nds; ++s) {
5112: if (dm->probs[s].label == label) {
5113: PetscDSDestroy(&dm->probs[s].ds);
5114: dm->probs[s].ds = ds;
5115: return(0);
5116: }
5117: }
5118: DMDSEnlarge_Static(dm, Nds+1);
5119: PetscObjectReference((PetscObject) label);
5120: PetscObjectReference((PetscObject) fields);
5121: PetscObjectReference((PetscObject) ds);
5122: if (!label) {
5123: /* Put the NULL label at the front, so it is returned as the default */
5124: for (s = Nds-1; s >=0; --s) dm->probs[s+1] = dm->probs[s];
5125: Nds = 0;
5126: }
5127: dm->probs[Nds].label = label;
5128: dm->probs[Nds].fields = fields;
5129: dm->probs[Nds].ds = ds;
5130: return(0);
5131: }
5133: /*@
5134: DMCreateDS - Create the discrete systems for the DM based upon the fields added to the DM
5136: Collective on dm
5138: Input Parameter:
5139: . dm - The DM
5141: Note: If the label has a DS defined, it will be replaced. Otherwise, it will be added to the DM.
5143: Level: intermediate
5145: .seealso: DMSetField, DMAddField(), DMGetDS(), DMGetCellDS(), DMGetRegionDS(), DMSetRegionDS()
5146: @*/
5147: PetscErrorCode DMCreateDS(DM dm)
5148: {
5149: MPI_Comm comm;
5150: PetscDS prob, probh = NULL;
5151: PetscInt dimEmbed, Nf = dm->Nf, f, s, field = 0, fieldh = 0;
5152: PetscBool doSetup = PETSC_TRUE;
5157: if (!dm->fields) return(0);
5158: /* Can only handle two label cases right now:
5159: 1) NULL
5160: 2) Hybrid cells
5161: */
5162: PetscObjectGetComm((PetscObject) dm, &comm);
5163: DMGetCoordinateDim(dm, &dimEmbed);
5164: /* Create default DS */
5165: DMGetRegionDS(dm, NULL, NULL, &prob);
5166: if (!prob) {
5167: IS fields;
5168: PetscInt *fld, nf;
5170: for (f = 0, nf = 0; f < Nf; ++f) if (!dm->fields[f].label) ++nf;
5171: PetscMalloc1(nf, &fld);
5172: for (f = 0, nf = 0; f < Nf; ++f) if (!dm->fields[f].label) fld[nf++] = f;
5173: ISCreate(PETSC_COMM_SELF, &fields);
5174: PetscObjectSetOptionsPrefix((PetscObject) fields, "dm_fields_");
5175: ISSetType(fields, ISGENERAL);
5176: ISGeneralSetIndices(fields, nf, fld, PETSC_OWN_POINTER);
5178: PetscDSCreate(comm, &prob);
5179: DMSetRegionDS(dm, NULL, fields, prob);
5180: PetscDSDestroy(&prob);
5181: ISDestroy(&fields);
5182: DMGetRegionDS(dm, NULL, NULL, &prob);
5183: }
5184: PetscDSSetCoordinateDimension(prob, dimEmbed);
5185: /* Optionally create hybrid DS */
5186: for (f = 0; f < Nf; ++f) {
5187: DMLabel label = dm->fields[f].label;
5188: PetscInt lStart, lEnd;
5190: if (label) {
5191: DM plex;
5192: IS fields;
5193: PetscInt *fld;
5194: PetscInt depth, pMax[4];
5196: DMConvert(dm, DMPLEX, &plex);
5197: DMPlexGetDepth(plex, &depth);
5198: DMPlexGetHybridBounds(plex, depth >= 0 ? &pMax[depth] : NULL, depth>1 ? &pMax[depth-1] : NULL, depth>2 ? &pMax[1] : NULL, &pMax[0]);
5199: DMDestroy(&plex);
5201: DMLabelGetBounds(label, &lStart, &lEnd);
5202: if (lStart < pMax[depth]) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only support labels over hybrid cells right now");
5203: PetscDSCreate(comm, &probh);
5204: PetscMalloc1(1, &fld);
5205: fld[0] = f;
5206: ISCreate(PETSC_COMM_SELF, &fields);
5207: PetscObjectSetOptionsPrefix((PetscObject) fields, "dm_fields_");
5208: ISSetType(fields, ISGENERAL);
5209: ISGeneralSetIndices(fields, 1, fld, PETSC_OWN_POINTER);
5210: DMSetRegionDS(dm, label, fields, probh);
5211: ISDestroy(&fields);
5212: PetscDSSetHybrid(probh, PETSC_TRUE);
5213: PetscDSSetCoordinateDimension(probh, dimEmbed);
5214: break;
5215: }
5216: }
5217: /* Set fields in DSes */
5218: for (f = 0; f < Nf; ++f) {
5219: DMLabel label = dm->fields[f].label;
5220: PetscObject disc = dm->fields[f].disc;
5222: if (!label) {
5223: PetscDSSetDiscretization(prob, field++, disc);
5224: if (probh) {
5225: PetscFE subfe;
5227: PetscFEGetHeightSubspace((PetscFE) disc, 1, &subfe);
5228: PetscDSSetDiscretization(probh, fieldh++, (PetscObject) subfe);
5229: }
5230: } else {
5231: PetscDSSetDiscretization(probh, fieldh++, disc);
5232: }
5233: /* We allow people to have placeholder fields and construct the Section by hand */
5234: {
5235: PetscClassId id;
5237: PetscObjectGetClassId(disc, &id);
5238: if ((id != PETSCFE_CLASSID) && (id != PETSCFV_CLASSID)) doSetup = PETSC_FALSE;
5239: }
5240: }
5241: PetscDSDestroy(&probh);
5242: /* Setup DSes */
5243: if (doSetup) {
5244: for (s = 0; s < dm->Nds; ++s) {PetscDSSetUp(dm->probs[s].ds);}
5245: }
5246: return(0);
5247: }
5249: /*@
5250: DMCopyDS - Copy the discrete systems for the DM into another DM
5252: Collective on dm
5254: Input Parameter:
5255: . dm - The DM
5257: Output Parameter:
5258: . newdm - The DM
5260: Level: advanced
5262: .seealso: DMCopyFields(), DMAddField(), DMGetDS(), DMGetCellDS(), DMGetRegionDS(), DMSetRegionDS()
5263: @*/
5264: PetscErrorCode DMCopyDS(DM dm, DM newdm)
5265: {
5266: PetscInt Nds, s;
5270: if (dm == newdm) return(0);
5271: DMGetNumDS(dm, &Nds);
5272: DMClearDS(newdm);
5273: for (s = 0; s < Nds; ++s) {
5274: DMLabel label;
5275: IS fields;
5276: PetscDS ds;
5278: DMGetRegionNumDS(dm, s, &label, &fields, &ds);
5279: DMSetRegionDS(newdm, label, fields, ds);
5280: }
5281: return(0);
5282: }
5284: /*@
5285: DMCopyDisc - Copy the fields and discrete systems for the DM into another DM
5287: Collective on dm
5289: Input Parameter:
5290: . dm - The DM
5292: Output Parameter:
5293: . newdm - The DM
5295: Level: advanced
5297: .seealso: DMCopyFields(), DMCopyDS()
5298: @*/
5299: PetscErrorCode DMCopyDisc(DM dm, DM newdm)
5300: {
5304: if (dm == newdm) return(0);
5305: DMCopyFields(dm, newdm);
5306: DMCopyDS(dm, newdm);
5307: return(0);
5308: }
5310: PetscErrorCode DMRestrictHook_Coordinates(DM dm,DM dmc,void *ctx)
5311: {
5312: DM dm_coord,dmc_coord;
5314: Vec coords,ccoords;
5315: Mat inject;
5317: DMGetCoordinateDM(dm,&dm_coord);
5318: DMGetCoordinateDM(dmc,&dmc_coord);
5319: DMGetCoordinates(dm,&coords);
5320: DMGetCoordinates(dmc,&ccoords);
5321: if (coords && !ccoords) {
5322: DMCreateGlobalVector(dmc_coord,&ccoords);
5323: PetscObjectSetName((PetscObject)ccoords,"coordinates");
5324: DMCreateInjection(dmc_coord,dm_coord,&inject);
5325: MatRestrict(inject,coords,ccoords);
5326: MatDestroy(&inject);
5327: DMSetCoordinates(dmc,ccoords);
5328: VecDestroy(&ccoords);
5329: }
5330: return(0);
5331: }
5333: static PetscErrorCode DMSubDomainHook_Coordinates(DM dm,DM subdm,void *ctx)
5334: {
5335: DM dm_coord,subdm_coord;
5337: Vec coords,ccoords,clcoords;
5338: VecScatter *scat_i,*scat_g;
5340: DMGetCoordinateDM(dm,&dm_coord);
5341: DMGetCoordinateDM(subdm,&subdm_coord);
5342: DMGetCoordinates(dm,&coords);
5343: DMGetCoordinates(subdm,&ccoords);
5344: if (coords && !ccoords) {
5345: DMCreateGlobalVector(subdm_coord,&ccoords);
5346: PetscObjectSetName((PetscObject)ccoords,"coordinates");
5347: DMCreateLocalVector(subdm_coord,&clcoords);
5348: PetscObjectSetName((PetscObject)clcoords,"coordinates");
5349: DMCreateDomainDecompositionScatters(dm_coord,1,&subdm_coord,NULL,&scat_i,&scat_g);
5350: VecScatterBegin(scat_i[0],coords,ccoords,INSERT_VALUES,SCATTER_FORWARD);
5351: VecScatterEnd(scat_i[0],coords,ccoords,INSERT_VALUES,SCATTER_FORWARD);
5352: VecScatterBegin(scat_g[0],coords,clcoords,INSERT_VALUES,SCATTER_FORWARD);
5353: VecScatterEnd(scat_g[0],coords,clcoords,INSERT_VALUES,SCATTER_FORWARD);
5354: DMSetCoordinates(subdm,ccoords);
5355: DMSetCoordinatesLocal(subdm,clcoords);
5356: VecScatterDestroy(&scat_i[0]);
5357: VecScatterDestroy(&scat_g[0]);
5358: VecDestroy(&ccoords);
5359: VecDestroy(&clcoords);
5360: PetscFree(scat_i);
5361: PetscFree(scat_g);
5362: }
5363: return(0);
5364: }
5366: /*@
5367: DMGetDimension - Return the topological dimension of the DM
5369: Not collective
5371: Input Parameter:
5372: . dm - The DM
5374: Output Parameter:
5375: . dim - The topological dimension
5377: Level: beginner
5379: .seealso: DMSetDimension(), DMCreate()
5380: @*/
5381: PetscErrorCode DMGetDimension(DM dm, PetscInt *dim)
5382: {
5386: *dim = dm->dim;
5387: return(0);
5388: }
5390: /*@
5391: DMSetDimension - Set the topological dimension of the DM
5393: Collective on dm
5395: Input Parameters:
5396: + dm - The DM
5397: - dim - The topological dimension
5399: Level: beginner
5401: .seealso: DMGetDimension(), DMCreate()
5402: @*/
5403: PetscErrorCode DMSetDimension(DM dm, PetscInt dim)
5404: {
5405: PetscDS ds;
5411: dm->dim = dim;
5412: DMGetDS(dm, &ds);
5413: if (ds->dimEmbed < 0) {PetscDSSetCoordinateDimension(ds, dm->dim);}
5414: return(0);
5415: }
5417: /*@
5418: DMGetDimPoints - Get the half-open interval for all points of a given dimension
5420: Collective on dm
5422: Input Parameters:
5423: + dm - the DM
5424: - dim - the dimension
5426: Output Parameters:
5427: + pStart - The first point of the given dimension
5428: - pEnd - The first point following points of the given dimension
5430: Note:
5431: The points are vertices in the Hasse diagram encoding the topology. This is explained in
5432: https://arxiv.org/abs/0908.4427. If no points exist of this dimension in the storage scheme,
5433: then the interval is empty.
5435: Level: intermediate
5437: .seealso: DMPLEX, DMPlexGetDepthStratum(), DMPlexGetHeightStratum()
5438: @*/
5439: PetscErrorCode DMGetDimPoints(DM dm, PetscInt dim, PetscInt *pStart, PetscInt *pEnd)
5440: {
5441: PetscInt d;
5446: DMGetDimension(dm, &d);
5447: if ((dim < 0) || (dim > d)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension %d 1", dim, d);
5448: if (!dm->ops->getdimpoints) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "DM type %s does not implement DMGetDimPoints",((PetscObject)dm)->type_name);
5449: (*dm->ops->getdimpoints)(dm, dim, pStart, pEnd);
5450: return(0);
5451: }
5453: /*@
5454: DMSetCoordinates - Sets into the DM a global vector that holds the coordinates
5456: Collective on dm
5458: Input Parameters:
5459: + dm - the DM
5460: - c - coordinate vector
5462: Notes:
5463: The coordinates do include those for ghost points, which are in the local vector.
5465: The vector c should be destroyed by the caller.
5467: Level: intermediate
5469: .seealso: DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM()
5470: @*/
5471: PetscErrorCode DMSetCoordinates(DM dm, Vec c)
5472: {
5478: PetscObjectReference((PetscObject) c);
5479: VecDestroy(&dm->coordinates);
5480: dm->coordinates = c;
5481: VecDestroy(&dm->coordinatesLocal);
5482: DMCoarsenHookAdd(dm,DMRestrictHook_Coordinates,NULL,NULL);
5483: DMSubDomainHookAdd(dm,DMSubDomainHook_Coordinates,NULL,NULL);
5484: return(0);
5485: }
5487: /*@
5488: DMSetCoordinatesLocal - Sets into the DM a local vector that holds the coordinates
5490: Not collective
5492: Input Parameters:
5493: + dm - the DM
5494: - c - coordinate vector
5496: Notes:
5497: The coordinates of ghost points can be set using DMSetCoordinates()
5498: followed by DMGetCoordinatesLocal(). This is intended to enable the
5499: setting of ghost coordinates outside of the domain.
5501: The vector c should be destroyed by the caller.
5503: Level: intermediate
5505: .seealso: DMGetCoordinatesLocal(), DMSetCoordinates(), DMGetCoordinates(), DMGetCoordinateDM()
5506: @*/
5507: PetscErrorCode DMSetCoordinatesLocal(DM dm, Vec c)
5508: {
5514: PetscObjectReference((PetscObject) c);
5515: VecDestroy(&dm->coordinatesLocal);
5517: dm->coordinatesLocal = c;
5519: VecDestroy(&dm->coordinates);
5520: return(0);
5521: }
5523: /*@
5524: DMGetCoordinates - Gets a global vector with the coordinates associated with the DM.
5526: Collective on dm
5528: Input Parameter:
5529: . dm - the DM
5531: Output Parameter:
5532: . c - global coordinate vector
5534: Note:
5535: This is a borrowed reference, so the user should NOT destroy this vector
5537: Each process has only the local coordinates (does NOT have the ghost coordinates).
5539: For DMDA, in two and three dimensions coordinates are interlaced (x_0,y_0,x_1,y_1,...)
5540: and (x_0,y_0,z_0,x_1,y_1,z_1...)
5542: Level: intermediate
5544: .seealso: DMSetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM()
5545: @*/
5546: PetscErrorCode DMGetCoordinates(DM dm, Vec *c)
5547: {
5553: if (!dm->coordinates && dm->coordinatesLocal) {
5554: DM cdm = NULL;
5555: PetscBool localized;
5557: DMGetCoordinateDM(dm, &cdm);
5558: DMCreateGlobalVector(cdm, &dm->coordinates);
5559: DMGetCoordinatesLocalized(dm, &localized);
5560: /* Block size is not correctly set by CreateGlobalVector() if coordinates are localized */
5561: if (localized) {
5562: PetscInt cdim;
5564: DMGetCoordinateDim(dm, &cdim);
5565: VecSetBlockSize(dm->coordinates, cdim);
5566: }
5567: PetscObjectSetName((PetscObject) dm->coordinates, "coordinates");
5568: DMLocalToGlobalBegin(cdm, dm->coordinatesLocal, INSERT_VALUES, dm->coordinates);
5569: DMLocalToGlobalEnd(cdm, dm->coordinatesLocal, INSERT_VALUES, dm->coordinates);
5570: }
5571: *c = dm->coordinates;
5572: return(0);
5573: }
5575: /*@
5576: DMGetCoordinatesLocalSetUp - Prepares a local vector of coordinates, so that DMGetCoordinatesLocalNoncollective() can be used as non-collective afterwards.
5578: Collective on dm
5580: Input Parameter:
5581: . dm - the DM
5583: Level: advanced
5585: .seealso: DMGetCoordinatesLocalNoncollective()
5586: @*/
5587: PetscErrorCode DMGetCoordinatesLocalSetUp(DM dm)
5588: {
5593: if (!dm->coordinatesLocal && dm->coordinates) {
5594: DM cdm = NULL;
5595: PetscBool localized;
5597: DMGetCoordinateDM(dm, &cdm);
5598: DMCreateLocalVector(cdm, &dm->coordinatesLocal);
5599: DMGetCoordinatesLocalized(dm, &localized);
5600: /* Block size is not correctly set by CreateLocalVector() if coordinates are localized */
5601: if (localized) {
5602: PetscInt cdim;
5604: DMGetCoordinateDim(dm, &cdim);
5605: VecSetBlockSize(dm->coordinates, cdim);
5606: }
5607: PetscObjectSetName((PetscObject) dm->coordinatesLocal, "coordinates");
5608: DMGlobalToLocalBegin(cdm, dm->coordinates, INSERT_VALUES, dm->coordinatesLocal);
5609: DMGlobalToLocalEnd(cdm, dm->coordinates, INSERT_VALUES, dm->coordinatesLocal);
5610: }
5611: return(0);
5612: }
5614: /*@
5615: DMGetCoordinatesLocal - Gets a local vector with the coordinates associated with the DM.
5617: Collective on dm
5619: Input Parameter:
5620: . dm - the DM
5622: Output Parameter:
5623: . c - coordinate vector
5625: Note:
5626: This is a borrowed reference, so the user should NOT destroy this vector
5628: Each process has the local and ghost coordinates
5630: For DMDA, in two and three dimensions coordinates are interlaced (x_0,y_0,x_1,y_1,...)
5631: and (x_0,y_0,z_0,x_1,y_1,z_1...)
5633: Level: intermediate
5635: .seealso: DMSetCoordinatesLocal(), DMGetCoordinates(), DMSetCoordinates(), DMGetCoordinateDM(), DMGetCoordinatesLocalNoncollective()
5636: @*/
5637: PetscErrorCode DMGetCoordinatesLocal(DM dm, Vec *c)
5638: {
5644: DMGetCoordinatesLocalSetUp(dm);
5645: *c = dm->coordinatesLocal;
5646: return(0);
5647: }
5649: /*@
5650: DMGetCoordinatesLocalNoncollective - Non-collective version of DMGetCoordinatesLocal(). Fails if global coordinates have been set and DMGetCoordinatesLocalSetUp() not called.
5652: Not collective
5654: Input Parameter:
5655: . dm - the DM
5657: Output Parameter:
5658: . c - coordinate vector
5660: Level: advanced
5662: .seealso: DMGetCoordinatesLocalSetUp(), DMGetCoordinatesLocal(), DMSetCoordinatesLocal(), DMGetCoordinates(), DMSetCoordinates(), DMGetCoordinateDM()
5663: @*/
5664: PetscErrorCode DMGetCoordinatesLocalNoncollective(DM dm, Vec *c)
5665: {
5669: if (!dm->coordinatesLocal && dm->coordinates) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DMGetCoordinatesLocalSetUp() has not been called");
5670: *c = dm->coordinatesLocal;
5671: return(0);
5672: }
5674: /*@
5675: DMGetCoordinatesLocalTuple - Gets a local vector with the coordinates of specified points and section describing its layout.
5677: Not collective
5679: Input Parameter:
5680: + dm - the DM
5681: - p - the IS of points whose coordinates will be returned
5683: Output Parameter:
5684: + pCoordSection - the PetscSection describing the layout of pCoord, i.e. each point corresponds to one point in p, and DOFs correspond to coordinates
5685: - pCoord - the Vec with coordinates of points in p
5687: Note:
5688: DMGetCoordinatesLocalSetUp() must be called first. This function employs DMGetCoordinatesLocalNoncollective() so it is not collective.
5690: This creates a new vector, so the user SHOULD destroy this vector
5692: Each process has the local and ghost coordinates
5694: For DMDA, in two and three dimensions coordinates are interlaced (x_0,y_0,x_1,y_1,...)
5695: and (x_0,y_0,z_0,x_1,y_1,z_1...)
5697: Level: advanced
5699: .seealso: DMSetCoordinatesLocal(), DMGetCoordinatesLocal(), DMGetCoordinatesLocalNoncollective(), DMGetCoordinatesLocalSetUp(), DMGetCoordinates(), DMSetCoordinates(), DMGetCoordinateDM()
5700: @*/
5701: PetscErrorCode DMGetCoordinatesLocalTuple(DM dm, IS p, PetscSection *pCoordSection, Vec *pCoord)
5702: {
5703: PetscSection cs, newcs;
5704: Vec coords;
5705: const PetscScalar *arr;
5706: PetscScalar *newarr=NULL;
5707: PetscInt n;
5708: PetscErrorCode ierr;
5715: if (!dm->coordinatesLocal) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DMGetCoordinatesLocalSetUp() has not been called or coordinates not set");
5716: if (!dm->coordinateDM || !dm->coordinateDM->localSection) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DM not supported");
5717: cs = dm->coordinateDM->localSection;
5718: coords = dm->coordinatesLocal;
5719: VecGetArrayRead(coords, &arr);
5720: PetscSectionExtractDofsFromArray(cs, MPIU_SCALAR, arr, p, &newcs, pCoord ? ((void**)&newarr) : NULL);
5721: VecRestoreArrayRead(coords, &arr);
5722: if (pCoord) {
5723: PetscSectionGetStorageSize(newcs, &n);
5724: /* set array in two steps to mimic PETSC_OWN_POINTER */
5725: VecCreateSeqWithArray(PetscObjectComm((PetscObject)p), 1, n, NULL, pCoord);
5726: VecReplaceArray(*pCoord, newarr);
5727: } else {
5728: PetscFree(newarr);
5729: }
5730: if (pCoordSection) {*pCoordSection = newcs;}
5731: else {PetscSectionDestroy(&newcs);}
5732: return(0);
5733: }
5735: PetscErrorCode DMGetCoordinateField(DM dm, DMField *field)
5736: {
5742: if (!dm->coordinateField) {
5743: if (dm->ops->createcoordinatefield) {
5744: (*dm->ops->createcoordinatefield)(dm,&dm->coordinateField);
5745: }
5746: }
5747: *field = dm->coordinateField;
5748: return(0);
5749: }
5751: PetscErrorCode DMSetCoordinateField(DM dm, DMField field)
5752: {
5758: PetscObjectReference((PetscObject)field);
5759: DMFieldDestroy(&dm->coordinateField);
5760: dm->coordinateField = field;
5761: return(0);
5762: }
5764: /*@
5765: DMGetCoordinateDM - Gets the DM that prescribes coordinate layout and scatters between global and local coordinates
5767: Collective on dm
5769: Input Parameter:
5770: . dm - the DM
5772: Output Parameter:
5773: . cdm - coordinate DM
5775: Level: intermediate
5777: .seealso: DMSetCoordinateDM(), DMSetCoordinates(), DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLocal()
5778: @*/
5779: PetscErrorCode DMGetCoordinateDM(DM dm, DM *cdm)
5780: {
5786: if (!dm->coordinateDM) {
5787: DM cdm;
5789: if (!dm->ops->createcoordinatedm) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unable to create coordinates for this DM");
5790: (*dm->ops->createcoordinatedm)(dm, &cdm);
5791: /* Just in case the DM sets the coordinate DM when creating it (DMP4est can do this, because it may not setup
5792: * until the call to CreateCoordinateDM) */
5793: DMDestroy(&dm->coordinateDM);
5794: dm->coordinateDM = cdm;
5795: }
5796: *cdm = dm->coordinateDM;
5797: return(0);
5798: }
5800: /*@
5801: DMSetCoordinateDM - Sets the DM that prescribes coordinate layout and scatters between global and local coordinates
5803: Logically Collective on dm
5805: Input Parameters:
5806: + dm - the DM
5807: - cdm - coordinate DM
5809: Level: intermediate
5811: .seealso: DMGetCoordinateDM(), DMSetCoordinates(), DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLocal()
5812: @*/
5813: PetscErrorCode DMSetCoordinateDM(DM dm, DM cdm)
5814: {
5820: PetscObjectReference((PetscObject)cdm);
5821: DMDestroy(&dm->coordinateDM);
5822: dm->coordinateDM = cdm;
5823: return(0);
5824: }
5826: /*@
5827: DMGetCoordinateDim - Retrieve the dimension of embedding space for coordinate values.
5829: Not Collective
5831: Input Parameter:
5832: . dm - The DM object
5834: Output Parameter:
5835: . dim - The embedding dimension
5837: Level: intermediate
5839: .seealso: DMSetCoordinateDim(), DMGetCoordinateSection(), DMGetCoordinateDM(), DMGetLocalSection(), DMSetLocalSection()
5840: @*/
5841: PetscErrorCode DMGetCoordinateDim(DM dm, PetscInt *dim)
5842: {
5846: if (dm->dimEmbed == PETSC_DEFAULT) {
5847: dm->dimEmbed = dm->dim;
5848: }
5849: *dim = dm->dimEmbed;
5850: return(0);
5851: }
5853: /*@
5854: DMSetCoordinateDim - Set the dimension of the embedding space for coordinate values.
5856: Not Collective
5858: Input Parameters:
5859: + dm - The DM object
5860: - dim - The embedding dimension
5862: Level: intermediate
5864: .seealso: DMGetCoordinateDim(), DMSetCoordinateSection(), DMGetCoordinateSection(), DMGetLocalSection(), DMSetLocalSection()
5865: @*/
5866: PetscErrorCode DMSetCoordinateDim(DM dm, PetscInt dim)
5867: {
5868: PetscDS ds;
5873: dm->dimEmbed = dim;
5874: DMGetDS(dm, &ds);
5875: PetscDSSetCoordinateDimension(ds, dim);
5876: return(0);
5877: }
5879: /*@
5880: DMGetCoordinateSection - Retrieve the layout of coordinate values over the mesh.
5882: Collective on dm
5884: Input Parameter:
5885: . dm - The DM object
5887: Output Parameter:
5888: . section - The PetscSection object
5890: Level: intermediate
5892: .seealso: DMGetCoordinateDM(), DMGetLocalSection(), DMSetLocalSection()
5893: @*/
5894: PetscErrorCode DMGetCoordinateSection(DM dm, PetscSection *section)
5895: {
5896: DM cdm;
5902: DMGetCoordinateDM(dm, &cdm);
5903: DMGetLocalSection(cdm, section);
5904: return(0);
5905: }
5907: /*@
5908: DMSetCoordinateSection - Set the layout of coordinate values over the mesh.
5910: Not Collective
5912: Input Parameters:
5913: + dm - The DM object
5914: . dim - The embedding dimension, or PETSC_DETERMINE
5915: - section - The PetscSection object
5917: Level: intermediate
5919: .seealso: DMGetCoordinateSection(), DMGetLocalSection(), DMSetLocalSection()
5920: @*/
5921: PetscErrorCode DMSetCoordinateSection(DM dm, PetscInt dim, PetscSection section)
5922: {
5923: DM cdm;
5929: DMGetCoordinateDM(dm, &cdm);
5930: DMSetLocalSection(cdm, section);
5931: if (dim == PETSC_DETERMINE) {
5932: PetscInt d = PETSC_DEFAULT;
5933: PetscInt pStart, pEnd, vStart, vEnd, v, dd;
5935: PetscSectionGetChart(section, &pStart, &pEnd);
5936: DMGetDimPoints(dm, 0, &vStart, &vEnd);
5937: pStart = PetscMax(vStart, pStart);
5938: pEnd = PetscMin(vEnd, pEnd);
5939: for (v = pStart; v < pEnd; ++v) {
5940: PetscSectionGetDof(section, v, &dd);
5941: if (dd) {d = dd; break;}
5942: }
5943: if (d >= 0) {DMSetCoordinateDim(dm, d);}
5944: }
5945: return(0);
5946: }
5948: /*@C
5949: DMGetPeriodicity - Get the description of mesh periodicity
5951: Input Parameters:
5952: . dm - The DM object
5954: Output Parameters:
5955: + per - Whether the DM is periodic or not
5956: . maxCell - Over distances greater than this, we can assume a point has crossed over to another sheet, when trying to localize cell coordinates
5957: . L - If we assume the mesh is a torus, this is the length of each coordinate
5958: - bd - This describes the type of periodicity in each topological dimension
5960: Level: developer
5962: .seealso: DMGetPeriodicity()
5963: @*/
5964: PetscErrorCode DMGetPeriodicity(DM dm, PetscBool *per, const PetscReal **maxCell, const PetscReal **L, const DMBoundaryType **bd)
5965: {
5968: if (per) *per = dm->periodic;
5969: if (L) *L = dm->L;
5970: if (maxCell) *maxCell = dm->maxCell;
5971: if (bd) *bd = dm->bdtype;
5972: return(0);
5973: }
5975: /*@C
5976: DMSetPeriodicity - Set the description of mesh periodicity
5978: Input Parameters:
5979: + dm - The DM object
5980: . per - Whether the DM is periodic or not. If maxCell is not provided, coordinates need to be localized
5981: . maxCell - Over distances greater than this, we can assume a point has crossed over to another sheet, when trying to localize cell coordinates
5982: . L - If we assume the mesh is a torus, this is the length of each coordinate
5983: - bd - This describes the type of periodicity in each topological dimension
5985: Level: developer
5987: .seealso: DMGetPeriodicity()
5988: @*/
5989: PetscErrorCode DMSetPeriodicity(DM dm, PetscBool per, const PetscReal maxCell[], const PetscReal L[], const DMBoundaryType bd[])
5990: {
5991: PetscInt dim, d;
5997: if (maxCell) {
6001: }
6002: PetscFree3(dm->L,dm->maxCell,dm->bdtype);
6003: DMGetDimension(dm, &dim);
6004: if (maxCell) {
6005: PetscMalloc3(dim,&dm->L,dim,&dm->maxCell,dim,&dm->bdtype);
6006: for (d = 0; d < dim; ++d) {dm->L[d] = L[d]; dm->maxCell[d] = maxCell[d]; dm->bdtype[d] = bd[d];}
6007: }
6008: dm->periodic = per;
6009: return(0);
6010: }
6012: /*@
6013: 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.
6015: Input Parameters:
6016: + dm - The DM
6017: . in - The input coordinate point (dim numbers)
6018: - endpoint - Include the endpoint L_i
6020: Output Parameter:
6021: . out - The localized coordinate point
6023: Level: developer
6025: .seealso: DMLocalizeCoordinates(), DMLocalizeAddCoordinate()
6026: @*/
6027: PetscErrorCode DMLocalizeCoordinate(DM dm, const PetscScalar in[], PetscBool endpoint, PetscScalar out[])
6028: {
6029: PetscInt dim, d;
6033: DMGetCoordinateDim(dm, &dim);
6034: if (!dm->maxCell) {
6035: for (d = 0; d < dim; ++d) out[d] = in[d];
6036: } else {
6037: if (endpoint) {
6038: for (d = 0; d < dim; ++d) {
6039: 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)) {
6040: out[d] = in[d] - dm->L[d]*(PetscFloorReal(PetscRealPart(in[d])/dm->L[d]) - 1);
6041: } else {
6042: out[d] = in[d] - dm->L[d]*PetscFloorReal(PetscRealPart(in[d])/dm->L[d]);
6043: }
6044: }
6045: } else {
6046: for (d = 0; d < dim; ++d) {
6047: out[d] = in[d] - dm->L[d]*PetscFloorReal(PetscRealPart(in[d])/dm->L[d]);
6048: }
6049: }
6050: }
6051: return(0);
6052: }
6054: /*
6055: 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.
6057: Input Parameters:
6058: + dm - The DM
6059: . dim - The spatial dimension
6060: . anchor - The anchor point, the input point can be no more than maxCell away from it
6061: - in - The input coordinate point (dim numbers)
6063: Output Parameter:
6064: . out - The localized coordinate point
6066: Level: developer
6068: 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
6070: .seealso: DMLocalizeCoordinates(), DMLocalizeAddCoordinate()
6071: */
6072: PetscErrorCode DMLocalizeCoordinate_Internal(DM dm, PetscInt dim, const PetscScalar anchor[], const PetscScalar in[], PetscScalar out[])
6073: {
6074: PetscInt d;
6077: if (!dm->maxCell) {
6078: for (d = 0; d < dim; ++d) out[d] = in[d];
6079: } else {
6080: for (d = 0; d < dim; ++d) {
6081: if ((dm->bdtype[d] != DM_BOUNDARY_NONE) && (PetscAbsScalar(anchor[d] - in[d]) > dm->maxCell[d])) {
6082: out[d] = PetscRealPart(anchor[d]) > PetscRealPart(in[d]) ? dm->L[d] + in[d] : in[d] - dm->L[d];
6083: } else {
6084: out[d] = in[d];
6085: }
6086: }
6087: }
6088: return(0);
6089: }
6090: PetscErrorCode DMLocalizeCoordinateReal_Internal(DM dm, PetscInt dim, const PetscReal anchor[], const PetscReal in[], PetscReal out[])
6091: {
6092: PetscInt d;
6095: if (!dm->maxCell) {
6096: for (d = 0; d < dim; ++d) out[d] = in[d];
6097: } else {
6098: for (d = 0; d < dim; ++d) {
6099: if ((dm->bdtype[d] != DM_BOUNDARY_NONE) && (PetscAbsReal(anchor[d] - in[d]) > dm->maxCell[d])) {
6100: out[d] = anchor[d] > in[d] ? dm->L[d] + in[d] : in[d] - dm->L[d];
6101: } else {
6102: out[d] = in[d];
6103: }
6104: }
6105: }
6106: return(0);
6107: }
6109: /*
6110: 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.
6112: Input Parameters:
6113: + dm - The DM
6114: . dim - The spatial dimension
6115: . anchor - The anchor point, the input point can be no more than maxCell away from it
6116: . in - The input coordinate delta (dim numbers)
6117: - out - The input coordinate point (dim numbers)
6119: Output Parameter:
6120: . out - The localized coordinate in + out
6122: Level: developer
6124: 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
6126: .seealso: DMLocalizeCoordinates(), DMLocalizeCoordinate()
6127: */
6128: PetscErrorCode DMLocalizeAddCoordinate_Internal(DM dm, PetscInt dim, const PetscScalar anchor[], const PetscScalar in[], PetscScalar out[])
6129: {
6130: PetscInt d;
6133: if (!dm->maxCell) {
6134: for (d = 0; d < dim; ++d) out[d] += in[d];
6135: } else {
6136: for (d = 0; d < dim; ++d) {
6137: if ((dm->bdtype[d] != DM_BOUNDARY_NONE) && (PetscAbsScalar(anchor[d] - in[d]) > dm->maxCell[d])) {
6138: out[d] += PetscRealPart(anchor[d]) > PetscRealPart(in[d]) ? dm->L[d] + in[d] : in[d] - dm->L[d];
6139: } else {
6140: out[d] += in[d];
6141: }
6142: }
6143: }
6144: return(0);
6145: }
6147: /*@
6148: DMGetCoordinatesLocalizedLocal - Check if the DM coordinates have been localized for cells on this process
6150: Not collective
6152: Input Parameter:
6153: . dm - The DM
6155: Output Parameter:
6156: areLocalized - True if localized
6158: Level: developer
6160: .seealso: DMLocalizeCoordinates(), DMGetCoordinatesLocalized(), DMSetPeriodicity()
6161: @*/
6162: PetscErrorCode DMGetCoordinatesLocalizedLocal(DM dm,PetscBool *areLocalized)
6163: {
6164: DM cdm;
6165: PetscSection coordSection;
6166: PetscInt cStart, cEnd, sStart, sEnd, c, dof;
6167: PetscBool isPlex, alreadyLocalized;
6173: *areLocalized = PETSC_FALSE;
6175: /* We need some generic way of refering to cells/vertices */
6176: DMGetCoordinateDM(dm, &cdm);
6177: PetscObjectTypeCompare((PetscObject) cdm, DMPLEX, &isPlex);
6178: if (!isPlex) return(0);
6180: DMGetCoordinateSection(dm, &coordSection);
6181: DMPlexGetHeightStratum(cdm, 0, &cStart, &cEnd);
6182: PetscSectionGetChart(coordSection, &sStart, &sEnd);
6183: alreadyLocalized = PETSC_FALSE;
6184: for (c = cStart; c < cEnd; ++c) {
6185: if (c < sStart || c >= sEnd) continue;
6186: PetscSectionGetDof(coordSection, c, &dof);
6187: if (dof) { alreadyLocalized = PETSC_TRUE; break; }
6188: }
6189: *areLocalized = alreadyLocalized;
6190: return(0);
6191: }
6193: /*@
6194: DMGetCoordinatesLocalized - Check if the DM coordinates have been localized for cells
6196: Collective on dm
6198: Input Parameter:
6199: . dm - The DM
6201: Output Parameter:
6202: areLocalized - True if localized
6204: Level: developer
6206: .seealso: DMLocalizeCoordinates(), DMSetPeriodicity(), DMGetCoordinatesLocalizedLocal()
6207: @*/
6208: PetscErrorCode DMGetCoordinatesLocalized(DM dm,PetscBool *areLocalized)
6209: {
6210: PetscBool localized;
6216: DMGetCoordinatesLocalizedLocal(dm,&localized);
6217: MPIU_Allreduce(&localized,areLocalized,1,MPIU_BOOL,MPI_LOR,PetscObjectComm((PetscObject)dm));
6218: return(0);
6219: }
6221: /*@
6222: DMLocalizeCoordinates - If a mesh is periodic, create local coordinates for cells having periodic faces
6224: Collective on dm
6226: Input Parameter:
6227: . dm - The DM
6229: Level: developer
6231: .seealso: DMSetPeriodicity(), DMLocalizeCoordinate(), DMLocalizeAddCoordinate()
6232: @*/
6233: PetscErrorCode DMLocalizeCoordinates(DM dm)
6234: {
6235: DM cdm;
6236: PetscSection coordSection, cSection;
6237: Vec coordinates, cVec;
6238: PetscScalar *coords, *coords2, *anchor, *localized;
6239: PetscInt Nc, vStart, vEnd, v, sStart, sEnd, newStart = PETSC_MAX_INT, newEnd = PETSC_MIN_INT, dof, d, off, off2, bs, coordSize;
6240: PetscBool alreadyLocalized, alreadyLocalizedGlobal;
6241: PetscInt maxHeight = 0, h;
6242: PetscInt *pStart = NULL, *pEnd = NULL;
6247: if (!dm->periodic) return(0);
6248: DMGetCoordinatesLocalized(dm, &alreadyLocalized);
6249: if (alreadyLocalized) return(0);
6251: /* We need some generic way of refering to cells/vertices */
6252: DMGetCoordinateDM(dm, &cdm);
6253: {
6254: PetscBool isplex;
6256: PetscObjectTypeCompare((PetscObject) cdm, DMPLEX, &isplex);
6257: if (isplex) {
6258: DMPlexGetDepthStratum(cdm, 0, &vStart, &vEnd);
6259: DMPlexGetMaxProjectionHeight(cdm,&maxHeight);
6260: DMGetWorkArray(dm,2*(maxHeight + 1),MPIU_INT,&pStart);
6261: pEnd = &pStart[maxHeight + 1];
6262: newStart = vStart;
6263: newEnd = vEnd;
6264: for (h = 0; h <= maxHeight; h++) {
6265: DMPlexGetHeightStratum(cdm, h, &pStart[h], &pEnd[h]);
6266: newStart = PetscMin(newStart,pStart[h]);
6267: newEnd = PetscMax(newEnd,pEnd[h]);
6268: }
6269: } else SETERRQ(PetscObjectComm((PetscObject) cdm), PETSC_ERR_ARG_WRONG, "Coordinate localization requires a DMPLEX coordinate DM");
6270: }
6271: DMGetCoordinatesLocal(dm, &coordinates);
6272: if (!coordinates) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Missing local coordinates vector");
6273: DMGetCoordinateSection(dm, &coordSection);
6274: VecGetBlockSize(coordinates, &bs);
6275: PetscSectionGetChart(coordSection,&sStart,&sEnd);
6277: PetscSectionCreate(PetscObjectComm((PetscObject) dm), &cSection);
6278: PetscSectionSetNumFields(cSection, 1);
6279: PetscSectionGetFieldComponents(coordSection, 0, &Nc);
6280: PetscSectionSetFieldComponents(cSection, 0, Nc);
6281: PetscSectionSetChart(cSection, newStart, newEnd);
6283: DMGetWorkArray(dm, 2 * bs, MPIU_SCALAR, &anchor);
6284: localized = &anchor[bs];
6285: alreadyLocalized = alreadyLocalizedGlobal = PETSC_TRUE;
6286: for (h = 0; h <= maxHeight; h++) {
6287: PetscInt cStart = pStart[h], cEnd = pEnd[h], c;
6289: for (c = cStart; c < cEnd; ++c) {
6290: PetscScalar *cellCoords = NULL;
6291: PetscInt b;
6293: if (c < sStart || c >= sEnd) alreadyLocalized = PETSC_FALSE;
6294: DMPlexVecGetClosure(cdm, coordSection, coordinates, c, &dof, &cellCoords);
6295: for (b = 0; b < bs; ++b) anchor[b] = cellCoords[b];
6296: for (d = 0; d < dof/bs; ++d) {
6297: DMLocalizeCoordinate_Internal(dm, bs, anchor, &cellCoords[d*bs], localized);
6298: for (b = 0; b < bs; b++) {
6299: if (cellCoords[d*bs + b] != localized[b]) break;
6300: }
6301: if (b < bs) break;
6302: }
6303: if (d < dof/bs) {
6304: if (c >= sStart && c < sEnd) {
6305: PetscInt cdof;
6307: PetscSectionGetDof(coordSection, c, &cdof);
6308: if (cdof != dof) alreadyLocalized = PETSC_FALSE;
6309: }
6310: PetscSectionSetDof(cSection, c, dof);
6311: PetscSectionSetFieldDof(cSection, c, 0, dof);
6312: }
6313: DMPlexVecRestoreClosure(cdm, coordSection, coordinates, c, &dof, &cellCoords);
6314: }
6315: }
6316: MPI_Allreduce(&alreadyLocalized,&alreadyLocalizedGlobal,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)dm));
6317: if (alreadyLocalizedGlobal) {
6318: DMRestoreWorkArray(dm, 2 * bs, MPIU_SCALAR, &anchor);
6319: PetscSectionDestroy(&cSection);
6320: DMRestoreWorkArray(dm,2*(maxHeight + 1),MPIU_INT,&pStart);
6321: return(0);
6322: }
6323: for (v = vStart; v < vEnd; ++v) {
6324: PetscSectionGetDof(coordSection, v, &dof);
6325: PetscSectionSetDof(cSection, v, dof);
6326: PetscSectionSetFieldDof(cSection, v, 0, dof);
6327: }
6328: PetscSectionSetUp(cSection);
6329: PetscSectionGetStorageSize(cSection, &coordSize);
6330: VecCreate(PETSC_COMM_SELF, &cVec);
6331: PetscObjectSetName((PetscObject)cVec,"coordinates");
6332: VecSetBlockSize(cVec, bs);
6333: VecSetSizes(cVec, coordSize, PETSC_DETERMINE);
6334: VecSetType(cVec, VECSTANDARD);
6335: VecGetArrayRead(coordinates, (const PetscScalar**)&coords);
6336: VecGetArray(cVec, &coords2);
6337: for (v = vStart; v < vEnd; ++v) {
6338: PetscSectionGetDof(coordSection, v, &dof);
6339: PetscSectionGetOffset(coordSection, v, &off);
6340: PetscSectionGetOffset(cSection, v, &off2);
6341: for (d = 0; d < dof; ++d) coords2[off2+d] = coords[off+d];
6342: }
6343: for (h = 0; h <= maxHeight; h++) {
6344: PetscInt cStart = pStart[h], cEnd = pEnd[h], c;
6346: for (c = cStart; c < cEnd; ++c) {
6347: PetscScalar *cellCoords = NULL;
6348: PetscInt b, cdof;
6350: PetscSectionGetDof(cSection,c,&cdof);
6351: if (!cdof) continue;
6352: DMPlexVecGetClosure(cdm, coordSection, coordinates, c, &dof, &cellCoords);
6353: PetscSectionGetOffset(cSection, c, &off2);
6354: for (b = 0; b < bs; ++b) anchor[b] = cellCoords[b];
6355: for (d = 0; d < dof/bs; ++d) {DMLocalizeCoordinate_Internal(dm, bs, anchor, &cellCoords[d*bs], &coords2[off2+d*bs]);}
6356: DMPlexVecRestoreClosure(cdm, coordSection, coordinates, c, &dof, &cellCoords);
6357: }
6358: }
6359: DMRestoreWorkArray(dm, 2 * bs, MPIU_SCALAR, &anchor);
6360: DMRestoreWorkArray(dm,2*(maxHeight + 1),MPIU_INT,&pStart);
6361: VecRestoreArrayRead(coordinates, (const PetscScalar**)&coords);
6362: VecRestoreArray(cVec, &coords2);
6363: DMSetCoordinateSection(dm, PETSC_DETERMINE, cSection);
6364: DMSetCoordinatesLocal(dm, cVec);
6365: VecDestroy(&cVec);
6366: PetscSectionDestroy(&cSection);
6367: return(0);
6368: }
6370: /*@
6371: DMLocatePoints - Locate the points in v in the mesh and return a PetscSF of the containing cells
6373: Collective on v (see explanation below)
6375: Input Parameters:
6376: + dm - The DM
6377: . v - The Vec of points
6378: . ltype - The type of point location, e.g. DM_POINTLOCATION_NONE or DM_POINTLOCATION_NEAREST
6379: - cells - Points to either NULL, or a PetscSF with guesses for which cells contain each point.
6381: Output Parameter:
6382: + v - The Vec of points, which now contains the nearest mesh points to the given points if DM_POINTLOCATION_NEAREST is used
6383: - cells - The PetscSF containing the ranks and local indices of the containing points.
6386: Level: developer
6388: Notes:
6389: To do a search of the local cells of the mesh, v should have PETSC_COMM_SELF as its communicator.
6390: To do a search of all the cells in the distributed mesh, v should have the same communicator as dm.
6392: If *cellSF is NULL on input, a PetscSF will be created.
6393: If *cellSF is not NULL on input, it should point to an existing PetscSF, whose graph will be used as initial guesses.
6395: An array that maps each point to its containing cell can be obtained with
6397: $ const PetscSFNode *cells;
6398: $ PetscInt nFound;
6399: $ const PetscInt *found;
6400: $
6401: $ PetscSFGetGraph(cellSF,NULL,&nFound,&found,&cells);
6403: Where cells[i].rank is the rank of the cell containing point found[i] (or i if found == NULL), and cells[i].index is
6404: the index of the cell in its rank's local numbering.
6406: .seealso: DMSetCoordinates(), DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLocal(), DMPointLocationType
6407: @*/
6408: PetscErrorCode DMLocatePoints(DM dm, Vec v, DMPointLocationType ltype, PetscSF *cellSF)
6409: {
6416: if (*cellSF) {
6417: PetscMPIInt result;
6420: MPI_Comm_compare(PetscObjectComm((PetscObject)v),PetscObjectComm((PetscObject)*cellSF),&result);
6421: if (result != MPI_IDENT && result != MPI_CONGRUENT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"cellSF must have a communicator congruent to v's");
6422: } else {
6423: PetscSFCreate(PetscObjectComm((PetscObject)v),cellSF);
6424: }
6425: if (!dm->ops->locatepoints) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Point location not available for this DM");
6426: PetscLogEventBegin(DM_LocatePoints,dm,0,0,0);
6427: (*dm->ops->locatepoints)(dm,v,ltype,*cellSF);
6428: PetscLogEventEnd(DM_LocatePoints,dm,0,0,0);
6429: return(0);
6430: }
6432: /*@
6433: DMGetOutputDM - Retrieve the DM associated with the layout for output
6435: Collective on dm
6437: Input Parameter:
6438: . dm - The original DM
6440: Output Parameter:
6441: . odm - The DM which provides the layout for output
6443: Level: intermediate
6445: .seealso: VecView(), DMGetGlobalSection()
6446: @*/
6447: PetscErrorCode DMGetOutputDM(DM dm, DM *odm)
6448: {
6449: PetscSection section;
6450: PetscBool hasConstraints, ghasConstraints;
6456: DMGetLocalSection(dm, §ion);
6457: PetscSectionHasConstraints(section, &hasConstraints);
6458: MPI_Allreduce(&hasConstraints, &ghasConstraints, 1, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject) dm));
6459: if (!ghasConstraints) {
6460: *odm = dm;
6461: return(0);
6462: }
6463: if (!dm->dmBC) {
6464: PetscSection newSection, gsection;
6465: PetscSF sf;
6467: DMClone(dm, &dm->dmBC);
6468: DMCopyDisc(dm, dm->dmBC);
6469: PetscSectionClone(section, &newSection);
6470: DMSetLocalSection(dm->dmBC, newSection);
6471: PetscSectionDestroy(&newSection);
6472: DMGetPointSF(dm->dmBC, &sf);
6473: PetscSectionCreateGlobalSection(section, sf, PETSC_TRUE, PETSC_FALSE, &gsection);
6474: DMSetGlobalSection(dm->dmBC, gsection);
6475: PetscSectionDestroy(&gsection);
6476: }
6477: *odm = dm->dmBC;
6478: return(0);
6479: }
6481: /*@
6482: DMGetOutputSequenceNumber - Retrieve the sequence number/value for output
6484: Input Parameter:
6485: . dm - The original DM
6487: Output Parameters:
6488: + num - The output sequence number
6489: - val - The output sequence value
6491: Level: intermediate
6493: Note: This is intended for output that should appear in sequence, for instance
6494: a set of timesteps in an HDF5 file, or a set of realizations of a stochastic system.
6496: .seealso: VecView()
6497: @*/
6498: PetscErrorCode DMGetOutputSequenceNumber(DM dm, PetscInt *num, PetscReal *val)
6499: {
6504: return(0);
6505: }
6507: /*@
6508: DMSetOutputSequenceNumber - Set the sequence number/value for output
6510: Input Parameters:
6511: + dm - The original DM
6512: . num - The output sequence number
6513: - val - The output sequence value
6515: Level: intermediate
6517: Note: This is intended for output that should appear in sequence, for instance
6518: a set of timesteps in an HDF5 file, or a set of realizations of a stochastic system.
6520: .seealso: VecView()
6521: @*/
6522: PetscErrorCode DMSetOutputSequenceNumber(DM dm, PetscInt num, PetscReal val)
6523: {
6526: dm->outputSequenceNum = num;
6527: dm->outputSequenceVal = val;
6528: return(0);
6529: }
6531: /*@C
6532: DMOutputSequenceLoad - Retrieve the sequence value from a Viewer
6534: Input Parameters:
6535: + dm - The original DM
6536: . name - The sequence name
6537: - num - The output sequence number
6539: Output Parameter:
6540: . val - The output sequence value
6542: Level: intermediate
6544: Note: This is intended for output that should appear in sequence, for instance
6545: a set of timesteps in an HDF5 file, or a set of realizations of a stochastic system.
6547: .seealso: DMGetOutputSequenceNumber(), DMSetOutputSequenceNumber(), VecView()
6548: @*/
6549: PetscErrorCode DMOutputSequenceLoad(DM dm, PetscViewer viewer, const char *name, PetscInt num, PetscReal *val)
6550: {
6551: PetscBool ishdf5;
6558: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERHDF5, &ishdf5);
6559: if (ishdf5) {
6560: #if defined(PETSC_HAVE_HDF5)
6561: PetscScalar value;
6563: DMSequenceLoad_HDF5_Internal(dm, name, num, &value, viewer);
6564: *val = PetscRealPart(value);
6565: #endif
6566: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid viewer; open viewer with PetscViewerHDF5Open()");
6567: return(0);
6568: }
6570: /*@
6571: DMGetUseNatural - Get the flag for creating a mapping to the natural order on distribution
6573: Not collective
6575: Input Parameter:
6576: . dm - The DM
6578: Output Parameter:
6579: . useNatural - The flag to build the mapping to a natural order during distribution
6581: Level: beginner
6583: .seealso: DMSetUseNatural(), DMCreate()
6584: @*/
6585: PetscErrorCode DMGetUseNatural(DM dm, PetscBool *useNatural)
6586: {
6590: *useNatural = dm->useNatural;
6591: return(0);
6592: }
6594: /*@
6595: DMSetUseNatural - Set the flag for creating a mapping to the natural order after distribution
6597: Collective on dm
6599: Input Parameters:
6600: + dm - The DM
6601: - useNatural - The flag to build the mapping to a natural order during distribution
6603: Note: This also causes the map to be build after DMCreateSubDM() and DMCreateSuperDM()
6605: Level: beginner
6607: .seealso: DMGetUseNatural(), DMCreate(), DMPlexDistribute(), DMCreateSubDM(), DMCreateSuperDM()
6608: @*/
6609: PetscErrorCode DMSetUseNatural(DM dm, PetscBool useNatural)
6610: {
6614: dm->useNatural = useNatural;
6615: return(0);
6616: }
6619: /*@C
6620: DMCreateLabel - Create a label of the given name if it does not already exist
6622: Not Collective
6624: Input Parameters:
6625: + dm - The DM object
6626: - name - The label name
6628: Level: intermediate
6630: .seealso: DMLabelCreate(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
6631: @*/
6632: PetscErrorCode DMCreateLabel(DM dm, const char name[])
6633: {
6634: DMLabelLink next = dm->labels->next;
6635: PetscBool flg = PETSC_FALSE;
6636: const char *lname;
6642: while (next) {
6643: PetscObjectGetName((PetscObject) next->label, &lname);
6644: PetscStrcmp(name, lname, &flg);
6645: if (flg) break;
6646: next = next->next;
6647: }
6648: if (!flg) {
6649: DMLabelLink tmpLabel;
6651: PetscCalloc1(1, &tmpLabel);
6652: DMLabelCreate(PETSC_COMM_SELF, name, &tmpLabel->label);
6653: tmpLabel->output = PETSC_TRUE;
6654: tmpLabel->next = dm->labels->next;
6655: dm->labels->next = tmpLabel;
6656: }
6657: return(0);
6658: }
6660: /*@C
6661: DMGetLabelValue - Get the value in a Sieve Label for the given point, with 0 as the default
6663: Not Collective
6665: Input Parameters:
6666: + dm - The DM object
6667: . name - The label name
6668: - point - The mesh point
6670: Output Parameter:
6671: . value - The label value for this point, or -1 if the point is not in the label
6673: Level: beginner
6675: .seealso: DMLabelGetValue(), DMSetLabelValue(), DMGetStratumIS()
6676: @*/
6677: PetscErrorCode DMGetLabelValue(DM dm, const char name[], PetscInt point, PetscInt *value)
6678: {
6679: DMLabel label;
6685: DMGetLabel(dm, name, &label);
6686: if (!label) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No label named %s was found", name);
6687: DMLabelGetValue(label, point, value);
6688: return(0);
6689: }
6691: /*@C
6692: DMSetLabelValue - Add a point to a Sieve Label with given value
6694: Not Collective
6696: Input Parameters:
6697: + dm - The DM object
6698: . name - The label name
6699: . point - The mesh point
6700: - value - The label value for this point
6702: Output Parameter:
6704: Level: beginner
6706: .seealso: DMLabelSetValue(), DMGetStratumIS(), DMClearLabelValue()
6707: @*/
6708: PetscErrorCode DMSetLabelValue(DM dm, const char name[], PetscInt point, PetscInt value)
6709: {
6710: DMLabel label;
6716: DMGetLabel(dm, name, &label);
6717: if (!label) {
6718: DMCreateLabel(dm, name);
6719: DMGetLabel(dm, name, &label);
6720: }
6721: DMLabelSetValue(label, point, value);
6722: return(0);
6723: }
6725: /*@C
6726: DMClearLabelValue - Remove a point from a Sieve Label with given value
6728: Not Collective
6730: Input Parameters:
6731: + dm - The DM object
6732: . name - The label name
6733: . point - The mesh point
6734: - value - The label value for this point
6736: Output Parameter:
6738: Level: beginner
6740: .seealso: DMLabelClearValue(), DMSetLabelValue(), DMGetStratumIS()
6741: @*/
6742: PetscErrorCode DMClearLabelValue(DM dm, const char name[], PetscInt point, PetscInt value)
6743: {
6744: DMLabel label;
6750: DMGetLabel(dm, name, &label);
6751: if (!label) return(0);
6752: DMLabelClearValue(label, point, value);
6753: return(0);
6754: }
6756: /*@C
6757: DMGetLabelSize - Get the number of different integer ids in a Label
6759: Not Collective
6761: Input Parameters:
6762: + dm - The DM object
6763: - name - The label name
6765: Output Parameter:
6766: . size - The number of different integer ids, or 0 if the label does not exist
6768: Level: beginner
6770: .seealso: DMLabelGetNumValues(), DMSetLabelValue()
6771: @*/
6772: PetscErrorCode DMGetLabelSize(DM dm, const char name[], PetscInt *size)
6773: {
6774: DMLabel label;
6781: DMGetLabel(dm, name, &label);
6782: *size = 0;
6783: if (!label) return(0);
6784: DMLabelGetNumValues(label, size);
6785: return(0);
6786: }
6788: /*@C
6789: DMGetLabelIdIS - Get the integer ids in a label
6791: Not Collective
6793: Input Parameters:
6794: + mesh - The DM object
6795: - name - The label name
6797: Output Parameter:
6798: . ids - The integer ids, or NULL if the label does not exist
6800: Level: beginner
6802: .seealso: DMLabelGetValueIS(), DMGetLabelSize()
6803: @*/
6804: PetscErrorCode DMGetLabelIdIS(DM dm, const char name[], IS *ids)
6805: {
6806: DMLabel label;
6813: DMGetLabel(dm, name, &label);
6814: *ids = NULL;
6815: if (label) {
6816: DMLabelGetValueIS(label, ids);
6817: } else {
6818: /* returning an empty IS */
6819: ISCreateGeneral(PETSC_COMM_SELF,0,NULL,PETSC_USE_POINTER,ids);
6820: }
6821: return(0);
6822: }
6824: /*@C
6825: DMGetStratumSize - Get the number of points in a label stratum
6827: Not Collective
6829: Input Parameters:
6830: + dm - The DM object
6831: . name - The label name
6832: - value - The stratum value
6834: Output Parameter:
6835: . size - The stratum size
6837: Level: beginner
6839: .seealso: DMLabelGetStratumSize(), DMGetLabelSize(), DMGetLabelIds()
6840: @*/
6841: PetscErrorCode DMGetStratumSize(DM dm, const char name[], PetscInt value, PetscInt *size)
6842: {
6843: DMLabel label;
6850: DMGetLabel(dm, name, &label);
6851: *size = 0;
6852: if (!label) return(0);
6853: DMLabelGetStratumSize(label, value, size);
6854: return(0);
6855: }
6857: /*@C
6858: DMGetStratumIS - Get the points in a label stratum
6860: Not Collective
6862: Input Parameters:
6863: + dm - The DM object
6864: . name - The label name
6865: - value - The stratum value
6867: Output Parameter:
6868: . points - The stratum points, or NULL if the label does not exist or does not have that value
6870: Level: beginner
6872: .seealso: DMLabelGetStratumIS(), DMGetStratumSize()
6873: @*/
6874: PetscErrorCode DMGetStratumIS(DM dm, const char name[], PetscInt value, IS *points)
6875: {
6876: DMLabel label;
6883: DMGetLabel(dm, name, &label);
6884: *points = NULL;
6885: if (!label) return(0);
6886: DMLabelGetStratumIS(label, value, points);
6887: return(0);
6888: }
6890: /*@C
6891: DMSetStratumIS - Set the points in a label stratum
6893: Not Collective
6895: Input Parameters:
6896: + dm - The DM object
6897: . name - The label name
6898: . value - The stratum value
6899: - points - The stratum points
6901: Level: beginner
6903: .seealso: DMLabelSetStratumIS(), DMGetStratumSize()
6904: @*/
6905: PetscErrorCode DMSetStratumIS(DM dm, const char name[], PetscInt value, IS points)
6906: {
6907: DMLabel label;
6914: DMGetLabel(dm, name, &label);
6915: if (!label) return(0);
6916: DMLabelSetStratumIS(label, value, points);
6917: return(0);
6918: }
6920: /*@C
6921: DMClearLabelStratum - Remove all points from a stratum from a Sieve Label
6923: Not Collective
6925: Input Parameters:
6926: + dm - The DM object
6927: . name - The label name
6928: - value - The label value for this point
6930: Output Parameter:
6932: Level: beginner
6934: .seealso: DMLabelClearStratum(), DMSetLabelValue(), DMGetStratumIS(), DMClearLabelValue()
6935: @*/
6936: PetscErrorCode DMClearLabelStratum(DM dm, const char name[], PetscInt value)
6937: {
6938: DMLabel label;
6944: DMGetLabel(dm, name, &label);
6945: if (!label) return(0);
6946: DMLabelClearStratum(label, value);
6947: return(0);
6948: }
6950: /*@
6951: DMGetNumLabels - Return the number of labels defined by the mesh
6953: Not Collective
6955: Input Parameter:
6956: . dm - The DM object
6958: Output Parameter:
6959: . numLabels - the number of Labels
6961: Level: intermediate
6963: .seealso: DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
6964: @*/
6965: PetscErrorCode DMGetNumLabels(DM dm, PetscInt *numLabels)
6966: {
6967: DMLabelLink next = dm->labels->next;
6968: PetscInt n = 0;
6973: while (next) {++n; next = next->next;}
6974: *numLabels = n;
6975: return(0);
6976: }
6978: /*@C
6979: DMGetLabelName - Return the name of nth label
6981: Not Collective
6983: Input Parameters:
6984: + dm - The DM object
6985: - n - the label number
6987: Output Parameter:
6988: . name - the label name
6990: Level: intermediate
6992: .seealso: DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
6993: @*/
6994: PetscErrorCode DMGetLabelName(DM dm, PetscInt n, const char **name)
6995: {
6996: DMLabelLink next = dm->labels->next;
6997: PetscInt l = 0;
7003: while (next) {
7004: if (l == n) {
7005: PetscObjectGetName((PetscObject) next->label, name);
7006: return(0);
7007: }
7008: ++l;
7009: next = next->next;
7010: }
7011: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label %D does not exist in this DM", n);
7012: }
7014: /*@C
7015: DMHasLabel - Determine whether the mesh has a label of a given name
7017: Not Collective
7019: Input Parameters:
7020: + dm - The DM object
7021: - name - The label name
7023: Output Parameter:
7024: . hasLabel - PETSC_TRUE if the label is present
7026: Level: intermediate
7028: .seealso: DMCreateLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
7029: @*/
7030: PetscErrorCode DMHasLabel(DM dm, const char name[], PetscBool *hasLabel)
7031: {
7032: DMLabelLink next = dm->labels->next;
7033: const char *lname;
7040: *hasLabel = PETSC_FALSE;
7041: while (next) {
7042: PetscObjectGetName((PetscObject) next->label, &lname);
7043: PetscStrcmp(name, lname, hasLabel);
7044: if (*hasLabel) break;
7045: next = next->next;
7046: }
7047: return(0);
7048: }
7050: /*@C
7051: DMGetLabel - Return the label of a given name, or NULL
7053: Not Collective
7055: Input Parameters:
7056: + dm - The DM object
7057: - name - The label name
7059: Output Parameter:
7060: . label - The DMLabel, or NULL if the label is absent
7062: Level: intermediate
7064: .seealso: DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
7065: @*/
7066: PetscErrorCode DMGetLabel(DM dm, const char name[], DMLabel *label)
7067: {
7068: DMLabelLink next = dm->labels->next;
7069: PetscBool hasLabel;
7070: const char *lname;
7077: *label = NULL;
7078: while (next) {
7079: PetscObjectGetName((PetscObject) next->label, &lname);
7080: PetscStrcmp(name, lname, &hasLabel);
7081: if (hasLabel) {
7082: *label = next->label;
7083: break;
7084: }
7085: next = next->next;
7086: }
7087: return(0);
7088: }
7090: /*@C
7091: DMGetLabelByNum - Return the nth label
7093: Not Collective
7095: Input Parameters:
7096: + dm - The DM object
7097: - n - the label number
7099: Output Parameter:
7100: . label - the label
7102: Level: intermediate
7104: .seealso: DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
7105: @*/
7106: PetscErrorCode DMGetLabelByNum(DM dm, PetscInt n, DMLabel *label)
7107: {
7108: DMLabelLink next = dm->labels->next;
7109: PetscInt l = 0;
7114: while (next) {
7115: if (l == n) {
7116: *label = next->label;
7117: return(0);
7118: }
7119: ++l;
7120: next = next->next;
7121: }
7122: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label %D does not exist in this DM", n);
7123: }
7125: /*@C
7126: DMAddLabel - Add the label to this mesh
7128: Not Collective
7130: Input Parameters:
7131: + dm - The DM object
7132: - label - The DMLabel
7134: Level: developer
7136: .seealso: DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
7137: @*/
7138: PetscErrorCode DMAddLabel(DM dm, DMLabel label)
7139: {
7140: DMLabelLink tmpLabel;
7141: PetscBool hasLabel;
7142: const char *lname;
7147: PetscObjectGetName((PetscObject) label, &lname);
7148: DMHasLabel(dm, lname, &hasLabel);
7149: if (hasLabel) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label %s already exists in this DM", lname);
7150: PetscCalloc1(1, &tmpLabel);
7151: tmpLabel->label = label;
7152: tmpLabel->output = PETSC_TRUE;
7153: tmpLabel->next = dm->labels->next;
7154: dm->labels->next = tmpLabel;
7155: PetscObjectReference((PetscObject)label);
7156: return(0);
7157: }
7159: /*@C
7160: DMRemoveLabel - Remove the label given by name from this mesh
7162: Not Collective
7164: Input Parameters:
7165: + dm - The DM object
7166: - name - The label name
7168: Output Parameter:
7169: . label - The DMLabel, or NULL if the label is absent
7171: Level: developer
7173: Notes:
7174: DMRemoveLabel(dm,name,NULL) removes the label from dm and calls
7175: DMLabelDestroy() on the label.
7177: DMRemoveLabel(dm,name,&label) removes the label from dm, but it DOES NOT
7178: call DMLabelDestroy(). Instead, the label is returned and the user is
7179: responsible of calling DMLabelDestroy() at some point.
7181: .seealso: DMCreateLabel(), DMHasLabel(), DMGetLabel(), DMGetLabelValue(), DMSetLabelValue(), DMLabelDestroy(), DMRemoveLabelBySelf()
7182: @*/
7183: PetscErrorCode DMRemoveLabel(DM dm, const char name[], DMLabel *label)
7184: {
7185: DMLabelLink link, *pnext;
7186: PetscBool hasLabel;
7187: const char *lname;
7193: if (label) {
7195: *label = NULL;
7196: }
7197: for (pnext=&dm->labels->next; (link=*pnext); pnext=&link->next) {
7198: PetscObjectGetName((PetscObject) link->label, &lname);
7199: PetscStrcmp(name, lname, &hasLabel);
7200: if (hasLabel) {
7201: *pnext = link->next; /* Remove from list */
7202: PetscStrcmp(name, "depth", &hasLabel);
7203: if (hasLabel) dm->depthLabel = NULL;
7204: if (label) *label = link->label;
7205: else {DMLabelDestroy(&link->label);}
7206: PetscFree(link);
7207: break;
7208: }
7209: }
7210: return(0);
7211: }
7213: /*@
7214: DMRemoveLabelBySelf - Remove the label from this mesh
7216: Not Collective
7218: Input Parameters:
7219: + dm - The DM object
7220: . label - (Optional) The DMLabel to be removed from the DM
7221: - failNotFound - Should it fail if the label is not found in the DM?
7223: Level: developer
7225: Notes:
7226: Only exactly the same instance is removed if found, name match is ignored.
7227: If the DM has an exclusive reference to the label, it gets destroyed and
7228: *label nullified.
7230: .seealso: DMCreateLabel(), DMHasLabel(), DMGetLabel() DMGetLabelValue(), DMSetLabelValue(), DMLabelDestroy(), DMRemoveLabel()
7231: @*/
7232: PetscErrorCode DMRemoveLabelBySelf(DM dm, DMLabel *label, PetscBool failNotFound)
7233: {
7234: DMLabelLink link, *pnext;
7235: PetscBool hasLabel = PETSC_FALSE;
7241: if (!*label && !failNotFound) return(0);
7244: for (pnext=&dm->labels->next; (link=*pnext); pnext=&link->next) {
7245: if (*label == link->label) {
7246: hasLabel = PETSC_TRUE;
7247: *pnext = link->next; /* Remove from list */
7248: if (*label == dm->depthLabel) dm->depthLabel = NULL;
7249: if (((PetscObject) link->label)->refct < 2) *label = NULL; /* nullify if exclusive reference */
7250: DMLabelDestroy(&link->label);
7251: PetscFree(link);
7252: break;
7253: }
7254: }
7255: if (!hasLabel && failNotFound) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Given label not found in DM");
7256: return(0);
7257: }
7259: /*@C
7260: DMGetLabelOutput - Get the output flag for a given label
7262: Not Collective
7264: Input Parameters:
7265: + dm - The DM object
7266: - name - The label name
7268: Output Parameter:
7269: . output - The flag for output
7271: Level: developer
7273: .seealso: DMSetLabelOutput(), DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
7274: @*/
7275: PetscErrorCode DMGetLabelOutput(DM dm, const char name[], PetscBool *output)
7276: {
7277: DMLabelLink next = dm->labels->next;
7278: const char *lname;
7285: while (next) {
7286: PetscBool flg;
7288: PetscObjectGetName((PetscObject) next->label, &lname);
7289: PetscStrcmp(name, lname, &flg);
7290: if (flg) {*output = next->output; return(0);}
7291: next = next->next;
7292: }
7293: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No label named %s was present in this dm", name);
7294: }
7296: /*@C
7297: DMSetLabelOutput - Set the output flag for a given label
7299: Not Collective
7301: Input Parameters:
7302: + dm - The DM object
7303: . name - The label name
7304: - output - The flag for output
7306: Level: developer
7308: .seealso: DMGetLabelOutput(), DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
7309: @*/
7310: PetscErrorCode DMSetLabelOutput(DM dm, const char name[], PetscBool output)
7311: {
7312: DMLabelLink next = dm->labels->next;
7313: const char *lname;
7319: while (next) {
7320: PetscBool flg;
7322: PetscObjectGetName((PetscObject) next->label, &lname);
7323: PetscStrcmp(name, lname, &flg);
7324: if (flg) {next->output = output; return(0);}
7325: next = next->next;
7326: }
7327: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No label named %s was present in this dm", name);
7328: }
7331: /*@
7332: DMCopyLabels - Copy labels from one mesh to another with a superset of the points
7334: Collective on dmA
7336: Input Parameter:
7337: . dmA - The DM object with initial labels
7339: Output Parameter:
7340: . dmB - The DM object with copied labels
7342: Level: intermediate
7344: Note: This is typically used when interpolating or otherwise adding to a mesh
7346: .seealso: DMGetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM(), DMGetCoordinateSection()
7347: @*/
7348: PetscErrorCode DMCopyLabels(DM dmA, DM dmB)
7349: {
7350: PetscInt numLabels, l;
7354: if (dmA == dmB) return(0);
7355: DMGetNumLabels(dmA, &numLabels);
7356: for (l = 0; l < numLabels; ++l) {
7357: DMLabel label, labelNew;
7358: const char *name;
7359: PetscBool flg;
7361: DMGetLabelName(dmA, l, &name);
7362: PetscStrcmp(name, "depth", &flg);
7363: if (flg) continue;
7364: PetscStrcmp(name, "dim", &flg);
7365: if (flg) continue;
7366: DMGetLabel(dmA, name, &label);
7367: DMLabelDuplicate(label, &labelNew);
7368: DMAddLabel(dmB, labelNew);
7369: DMLabelDestroy(&labelNew);
7370: }
7371: return(0);
7372: }
7374: /*@
7375: DMGetCoarseDM - Get the coarse mesh from which this was obtained by refinement
7377: Input Parameter:
7378: . dm - The DM object
7380: Output Parameter:
7381: . cdm - The coarse DM
7383: Level: intermediate
7385: .seealso: DMSetCoarseDM()
7386: @*/
7387: PetscErrorCode DMGetCoarseDM(DM dm, DM *cdm)
7388: {
7392: *cdm = dm->coarseMesh;
7393: return(0);
7394: }
7396: /*@
7397: DMSetCoarseDM - Set the coarse mesh from which this was obtained by refinement
7399: Input Parameters:
7400: + dm - The DM object
7401: - cdm - The coarse DM
7403: Level: intermediate
7405: .seealso: DMGetCoarseDM()
7406: @*/
7407: PetscErrorCode DMSetCoarseDM(DM dm, DM cdm)
7408: {
7414: PetscObjectReference((PetscObject)cdm);
7415: DMDestroy(&dm->coarseMesh);
7416: dm->coarseMesh = cdm;
7417: return(0);
7418: }
7420: /*@
7421: DMGetFineDM - Get the fine mesh from which this was obtained by refinement
7423: Input Parameter:
7424: . dm - The DM object
7426: Output Parameter:
7427: . fdm - The fine DM
7429: Level: intermediate
7431: .seealso: DMSetFineDM()
7432: @*/
7433: PetscErrorCode DMGetFineDM(DM dm, DM *fdm)
7434: {
7438: *fdm = dm->fineMesh;
7439: return(0);
7440: }
7442: /*@
7443: DMSetFineDM - Set the fine mesh from which this was obtained by refinement
7445: Input Parameters:
7446: + dm - The DM object
7447: - fdm - The fine DM
7449: Level: intermediate
7451: .seealso: DMGetFineDM()
7452: @*/
7453: PetscErrorCode DMSetFineDM(DM dm, DM fdm)
7454: {
7460: PetscObjectReference((PetscObject)fdm);
7461: DMDestroy(&dm->fineMesh);
7462: dm->fineMesh = fdm;
7463: return(0);
7464: }
7466: /*=== DMBoundary code ===*/
7468: PetscErrorCode DMCopyBoundary(DM dm, DM dmNew)
7469: {
7470: PetscInt d;
7474: for (d = 0; d < dm->Nds; ++d) {
7475: PetscDSCopyBoundary(dm->probs[d].ds, dmNew->probs[d].ds);
7476: }
7477: return(0);
7478: }
7480: /*@C
7481: DMAddBoundary - Add a boundary condition to the model
7483: Input Parameters:
7484: + dm - The DM, with a PetscDS that matches the problem being constrained
7485: . type - The type of condition, e.g. DM_BC_ESSENTIAL_ANALYTIC/DM_BC_ESSENTIAL_FIELD (Dirichlet), or DM_BC_NATURAL (Neumann)
7486: . name - The BC name
7487: . labelname - The label defining constrained points
7488: . field - The field to constrain
7489: . numcomps - The number of constrained field components (0 will constrain all fields)
7490: . comps - An array of constrained component numbers
7491: . bcFunc - A pointwise function giving boundary values
7492: . numids - The number of DMLabel ids for constrained points
7493: . ids - An array of ids for constrained points
7494: - ctx - An optional user context for bcFunc
7496: Options Database Keys:
7497: + -bc_<boundary name> <num> - Overrides the boundary ids
7498: - -bc_<boundary name>_comp <num> - Overrides the boundary components
7500: Level: developer
7502: .seealso: DMGetBoundary()
7503: @*/
7504: 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)
7505: {
7506: PetscDS ds;
7511: DMGetDS(dm, &ds);
7512: PetscDSAddBoundary(ds, type,name, labelname, field, numcomps, comps, bcFunc, numids, ids, ctx);
7513: return(0);
7514: }
7516: /*@
7517: DMGetNumBoundary - Get the number of registered BC
7519: Input Parameters:
7520: . dm - The mesh object
7522: Output Parameters:
7523: . numBd - The number of BC
7525: Level: intermediate
7527: .seealso: DMAddBoundary(), DMGetBoundary()
7528: @*/
7529: PetscErrorCode DMGetNumBoundary(DM dm, PetscInt *numBd)
7530: {
7531: PetscDS ds;
7536: DMGetDS(dm, &ds);
7537: PetscDSGetNumBoundary(ds, numBd);
7538: return(0);
7539: }
7541: /*@C
7542: DMGetBoundary - Get a model boundary condition
7544: Input Parameters:
7545: + dm - The mesh object
7546: - bd - The BC number
7548: Output Parameters:
7549: + type - The type of condition, e.g. DM_BC_ESSENTIAL_ANALYTIC/DM_BC_ESSENTIAL_FIELD (Dirichlet), or DM_BC_NATURAL (Neumann)
7550: . name - The BC name
7551: . labelname - The label defining constrained points
7552: . field - The field to constrain
7553: . numcomps - The number of constrained field components
7554: . comps - An array of constrained component numbers
7555: . bcFunc - A pointwise function giving boundary values
7556: . numids - The number of DMLabel ids for constrained points
7557: . ids - An array of ids for constrained points
7558: - ctx - An optional user context for bcFunc
7560: Options Database Keys:
7561: + -bc_<boundary name> <num> - Overrides the boundary ids
7562: - -bc_<boundary name>_comp <num> - Overrides the boundary components
7564: Level: developer
7566: .seealso: DMAddBoundary()
7567: @*/
7568: 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)
7569: {
7570: PetscDS ds;
7575: DMGetDS(dm, &ds);
7576: PetscDSGetBoundary(ds, bd, type, name, labelname, field, numcomps, comps, func, numids, ids, ctx);
7577: return(0);
7578: }
7580: static PetscErrorCode DMPopulateBoundary(DM dm)
7581: {
7582: PetscDS ds;
7583: DMBoundary *lastnext;
7584: DSBoundary dsbound;
7588: DMGetDS(dm, &ds);
7589: dsbound = ds->boundary;
7590: if (dm->boundary) {
7591: DMBoundary next = dm->boundary;
7593: /* quick check to see if the PetscDS has changed */
7594: if (next->dsboundary == dsbound) return(0);
7595: /* the PetscDS has changed: tear down and rebuild */
7596: while (next) {
7597: DMBoundary b = next;
7599: next = b->next;
7600: PetscFree(b);
7601: }
7602: dm->boundary = NULL;
7603: }
7605: lastnext = &(dm->boundary);
7606: while (dsbound) {
7607: DMBoundary dmbound;
7609: PetscNew(&dmbound);
7610: dmbound->dsboundary = dsbound;
7611: DMGetLabel(dm, dsbound->labelname, &(dmbound->label));
7612: if (!dmbound->label) {PetscInfo2(dm, "DSBoundary %s wants label %s, which is not in this dm.\n",dsbound->name,dsbound->labelname);}
7613: /* push on the back instead of the front so that it is in the same order as in the PetscDS */
7614: *lastnext = dmbound;
7615: lastnext = &(dmbound->next);
7616: dsbound = dsbound->next;
7617: }
7618: return(0);
7619: }
7621: PetscErrorCode DMIsBoundaryPoint(DM dm, PetscInt point, PetscBool *isBd)
7622: {
7623: DMBoundary b;
7629: *isBd = PETSC_FALSE;
7630: DMPopulateBoundary(dm);
7631: b = dm->boundary;
7632: while (b && !(*isBd)) {
7633: DMLabel label = b->label;
7634: DSBoundary dsb = b->dsboundary;
7636: if (label) {
7637: PetscInt i;
7639: for (i = 0; i < dsb->numids && !(*isBd); ++i) {
7640: DMLabelStratumHasPoint(label, dsb->ids[i], point, isBd);
7641: }
7642: }
7643: b = b->next;
7644: }
7645: return(0);
7646: }
7648: /*@C
7649: DMProjectFunction - This projects the given function into the function space provided.
7651: Input Parameters:
7652: + dm - The DM
7653: . time - The time
7654: . funcs - The coordinate functions to evaluate, one per field
7655: . ctxs - Optional array of contexts to pass to each coordinate function. ctxs itself may be null.
7656: - mode - The insertion mode for values
7658: Output Parameter:
7659: . X - vector
7661: Calling sequence of func:
7662: $ func(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar u[], void *ctx);
7664: + dim - The spatial dimension
7665: . x - The coordinates
7666: . Nf - The number of fields
7667: . u - The output field values
7668: - ctx - optional user-defined function context
7670: Level: developer
7672: .seealso: DMComputeL2Diff()
7673: @*/
7674: PetscErrorCode DMProjectFunction(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec X)
7675: {
7676: Vec localX;
7681: DMGetLocalVector(dm, &localX);
7682: DMProjectFunctionLocal(dm, time, funcs, ctxs, mode, localX);
7683: DMLocalToGlobalBegin(dm, localX, mode, X);
7684: DMLocalToGlobalEnd(dm, localX, mode, X);
7685: DMRestoreLocalVector(dm, &localX);
7686: return(0);
7687: }
7689: PetscErrorCode DMProjectFunctionLocal(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX)
7690: {
7696: if (!dm->ops->projectfunctionlocal) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMProjectFunctionLocal",((PetscObject)dm)->type_name);
7697: (dm->ops->projectfunctionlocal) (dm, time, funcs, ctxs, mode, localX);
7698: return(0);
7699: }
7701: 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)
7702: {
7703: Vec localX;
7708: DMGetLocalVector(dm, &localX);
7709: DMProjectFunctionLabelLocal(dm, time, label, numIds, ids, Nc, comps, funcs, ctxs, mode, localX);
7710: DMLocalToGlobalBegin(dm, localX, mode, X);
7711: DMLocalToGlobalEnd(dm, localX, mode, X);
7712: DMRestoreLocalVector(dm, &localX);
7713: return(0);
7714: }
7716: 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)
7717: {
7723: if (!dm->ops->projectfunctionlabellocal) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMProjectFunctionLabelLocal",((PetscObject)dm)->type_name);
7724: (dm->ops->projectfunctionlabellocal) (dm, time, label, numIds, ids, Nc, comps, funcs, ctxs, mode, localX);
7725: return(0);
7726: }
7728: PetscErrorCode DMProjectFieldLocal(DM dm, PetscReal time, Vec localU,
7729: void (**funcs)(PetscInt, PetscInt, PetscInt,
7730: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
7731: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
7732: PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]),
7733: InsertMode mode, Vec localX)
7734: {
7741: if (!dm->ops->projectfieldlocal) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMProjectFieldLocal",((PetscObject)dm)->type_name);
7742: (dm->ops->projectfieldlocal) (dm, time, localU, funcs, mode, localX);
7743: return(0);
7744: }
7746: PetscErrorCode DMProjectFieldLabelLocal(DM dm, PetscReal time, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscInt Nc, const PetscInt comps[], Vec localU,
7747: void (**funcs)(PetscInt, PetscInt, PetscInt,
7748: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
7749: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
7750: PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]),
7751: InsertMode mode, Vec localX)
7752: {
7759: if (!dm->ops->projectfieldlabellocal) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMProjectFieldLocal",((PetscObject)dm)->type_name);
7760: (dm->ops->projectfieldlabellocal)(dm, time, label, numIds, ids, Nc, comps, localU, funcs, mode, localX);
7761: return(0);
7762: }
7764: /*@C
7765: DMComputeL2Diff - This function computes the L_2 difference between a function u and an FEM interpolant solution u_h.
7767: Input Parameters:
7768: + dm - The DM
7769: . time - The time
7770: . funcs - The functions to evaluate for each field component
7771: . ctxs - Optional array of contexts to pass to each function, or NULL.
7772: - X - The coefficient vector u_h, a global vector
7774: Output Parameter:
7775: . diff - The diff ||u - u_h||_2
7777: Level: developer
7779: .seealso: DMProjectFunction(), DMComputeL2FieldDiff(), DMComputeL2GradientDiff()
7780: @*/
7781: PetscErrorCode DMComputeL2Diff(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff)
7782: {
7788: if (!dm->ops->computel2diff) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMComputeL2Diff",((PetscObject)dm)->type_name);
7789: (dm->ops->computel2diff)(dm,time,funcs,ctxs,X,diff);
7790: return(0);
7791: }
7793: /*@C
7794: DMComputeL2GradientDiff - This function computes the L_2 difference between the gradient of a function u and an FEM interpolant solution grad u_h.
7796: Collective on dm
7798: Input Parameters:
7799: + dm - The DM
7800: , time - The time
7801: . funcs - The gradient functions to evaluate for each field component
7802: . ctxs - Optional array of contexts to pass to each function, or NULL.
7803: . X - The coefficient vector u_h, a global vector
7804: - n - The vector to project along
7806: Output Parameter:
7807: . diff - The diff ||(grad u - grad u_h) . n||_2
7809: Level: developer
7811: .seealso: DMProjectFunction(), DMComputeL2Diff()
7812: @*/
7813: 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)
7814: {
7820: if (!dm->ops->computel2gradientdiff) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMComputeL2GradientDiff",((PetscObject)dm)->type_name);
7821: (dm->ops->computel2gradientdiff)(dm,time,funcs,ctxs,X,n,diff);
7822: return(0);
7823: }
7825: /*@C
7826: DMComputeL2FieldDiff - This function computes the L_2 difference between a function u and an FEM interpolant solution u_h, separated into field components.
7828: Collective on dm
7830: Input Parameters:
7831: + dm - The DM
7832: . time - The time
7833: . funcs - The functions to evaluate for each field component
7834: . ctxs - Optional array of contexts to pass to each function, or NULL.
7835: - X - The coefficient vector u_h, a global vector
7837: Output Parameter:
7838: . diff - The array of differences, ||u^f - u^f_h||_2
7840: Level: developer
7842: .seealso: DMProjectFunction(), DMComputeL2FieldDiff(), DMComputeL2GradientDiff()
7843: @*/
7844: PetscErrorCode DMComputeL2FieldDiff(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal diff[])
7845: {
7851: if (!dm->ops->computel2fielddiff) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMComputeL2FieldDiff",((PetscObject)dm)->type_name);
7852: (dm->ops->computel2fielddiff)(dm,time,funcs,ctxs,X,diff);
7853: return(0);
7854: }
7856: /*@C
7857: DMAdaptLabel - Adapt a dm based on a label with values interpreted as coarsening and refining flags. Specific implementations of DM maybe have
7858: specialized flags, but all implementations should accept flag values DM_ADAPT_DETERMINE, DM_ADAPT_KEEP, DM_ADAPT_REFINE, and DM_ADAPT_COARSEN.
7860: Collective on dm
7862: Input parameters:
7863: + dm - the pre-adaptation DM object
7864: - label - label with the flags
7866: Output parameters:
7867: . dmAdapt - the adapted DM object: may be NULL if an adapted DM could not be produced.
7869: Level: intermediate
7871: .seealso: DMAdaptMetric(), DMCoarsen(), DMRefine()
7872: @*/
7873: PetscErrorCode DMAdaptLabel(DM dm, DMLabel label, DM *dmAdapt)
7874: {
7881: *dmAdapt = NULL;
7882: if (!dm->ops->adaptlabel) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMAdaptLabel",((PetscObject)dm)->type_name);
7883: (dm->ops->adaptlabel)(dm, label, dmAdapt);
7884: return(0);
7885: }
7887: /*@C
7888: DMAdaptMetric - Generates a mesh adapted to the specified metric field using the pragmatic library.
7890: Input Parameters:
7891: + dm - The DM object
7892: . metric - The metric to which the mesh is adapted, defined vertex-wise.
7893: - 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_".
7895: Output Parameter:
7896: . dmAdapt - Pointer to the DM object containing the adapted mesh
7898: Note: The label in the adapted mesh will be registered under the name of the input DMLabel object
7900: Level: advanced
7902: .seealso: DMAdaptLabel(), DMCoarsen(), DMRefine()
7903: @*/
7904: PetscErrorCode DMAdaptMetric(DM dm, Vec metric, DMLabel bdLabel, DM *dmAdapt)
7905: {
7913: *dmAdapt = NULL;
7914: if (!dm->ops->adaptmetric) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMAdaptMetric",((PetscObject)dm)->type_name);
7915: (dm->ops->adaptmetric)(dm, metric, bdLabel, dmAdapt);
7916: return(0);
7917: }
7919: /*@C
7920: DMGetNeighbors - Gets an array containing the MPI rank of all the processes neighbors
7922: Not Collective
7924: Input Parameter:
7925: . dm - The DM
7927: Output Parameter:
7928: . nranks - the number of neighbours
7929: . ranks - the neighbors ranks
7931: Notes:
7932: Do not free the array, it is freed when the DM is destroyed.
7934: Level: beginner
7936: .seealso: DMDAGetNeighbors(), PetscSFGetRootRanks()
7937: @*/
7938: PetscErrorCode DMGetNeighbors(DM dm,PetscInt *nranks,const PetscMPIInt *ranks[])
7939: {
7944: if (!dm->ops->getneighbors) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMGetNeighbors",((PetscObject)dm)->type_name);
7945: (dm->ops->getneighbors)(dm,nranks,ranks);
7946: return(0);
7947: }
7949: #include <petsc/private/matimpl.h> /* Needed because of coloring->ctype below */
7951: /*
7952: Converts the input vector to a ghosted vector and then calls the standard coloring code.
7953: This has be a different function because it requires DM which is not defined in the Mat library
7954: */
7955: PetscErrorCode MatFDColoringApply_AIJDM(Mat J,MatFDColoring coloring,Vec x1,void *sctx)
7956: {
7960: if (coloring->ctype == IS_COLORING_LOCAL) {
7961: Vec x1local;
7962: DM dm;
7963: MatGetDM(J,&dm);
7964: if (!dm) SETERRQ(PetscObjectComm((PetscObject)J),PETSC_ERR_ARG_INCOMP,"IS_COLORING_LOCAL requires a DM");
7965: DMGetLocalVector(dm,&x1local);
7966: DMGlobalToLocalBegin(dm,x1,INSERT_VALUES,x1local);
7967: DMGlobalToLocalEnd(dm,x1,INSERT_VALUES,x1local);
7968: x1 = x1local;
7969: }
7970: MatFDColoringApply_AIJ(J,coloring,x1,sctx);
7971: if (coloring->ctype == IS_COLORING_LOCAL) {
7972: DM dm;
7973: MatGetDM(J,&dm);
7974: DMRestoreLocalVector(dm,&x1);
7975: }
7976: return(0);
7977: }
7979: /*@
7980: MatFDColoringUseDM - allows a MatFDColoring object to use the DM associated with the matrix to use a IS_COLORING_LOCAL coloring
7982: Input Parameter:
7983: . coloring - the MatFDColoring object
7985: Developer Notes:
7986: this routine exists because the PETSc Mat library does not know about the DM objects
7988: Level: advanced
7990: .seealso: MatFDColoring, MatFDColoringCreate(), ISColoringType
7991: @*/
7992: PetscErrorCode MatFDColoringUseDM(Mat coloring,MatFDColoring fdcoloring)
7993: {
7995: coloring->ops->fdcoloringapply = MatFDColoringApply_AIJDM;
7996: return(0);
7997: }
7999: /*@
8000: DMGetCompatibility - determine if two DMs are compatible
8002: Collective
8004: Input Parameters:
8005: + dm - the first DM
8006: - dm2 - the second DM
8008: Output Parameters:
8009: + compatible - whether or not the two DMs are compatible
8010: - set - whether or not the compatible value was set
8012: Notes:
8013: Two DMs are deemed compatible if they represent the same parallel decomposition
8014: of the same topology. This implies that the section (field data) on one
8015: "makes sense" with respect to the topology and parallel decomposition of the other.
8016: Loosely speaking, compatible DMs represent the same domain and parallel
8017: decomposition, but hold different data.
8019: Typically, one would confirm compatibility if intending to simultaneously iterate
8020: over a pair of vectors obtained from different DMs.
8022: For example, two DMDA objects are compatible if they have the same local
8023: and global sizes and the same stencil width. They can have different numbers
8024: of degrees of freedom per node. Thus, one could use the node numbering from
8025: either DM in bounds for a loop over vectors derived from either DM.
8027: Consider the operation of summing data living on a 2-dof DMDA to data living
8028: on a 1-dof DMDA, which should be compatible, as in the following snippet.
8029: .vb
8030: ...
8031: DMGetCompatibility(da1,da2,&compatible,&set);
8032: if (set && compatible) {
8033: DMDAVecGetArrayDOF(da1,vec1,&arr1);
8034: DMDAVecGetArrayDOF(da2,vec2,&arr2);
8035: DMDAGetCorners(da1,&x,&y,NULL,&m,&n,NULL);
8036: for (j=y; j<y+n; ++j) {
8037: for (i=x; i<x+m, ++i) {
8038: arr1[j][i][0] = arr2[j][i][0] + arr2[j][i][1];
8039: }
8040: }
8041: DMDAVecRestoreArrayDOF(da1,vec1,&arr1);
8042: DMDAVecRestoreArrayDOF(da2,vec2,&arr2);
8043: } else {
8044: SETERRQ(PetscObjectComm((PetscObject)da1,PETSC_ERR_ARG_INCOMP,"DMDA objects incompatible");
8045: }
8046: ...
8047: .ve
8049: Checking compatibility might be expensive for a given implementation of DM,
8050: or might be impossible to unambiguously confirm or deny. For this reason,
8051: this function may decline to determine compatibility, and hence users should
8052: always check the "set" output parameter.
8054: A DM is always compatible with itself.
8056: In the current implementation, DMs which live on "unequal" communicators
8057: (MPI_UNEQUAL in the terminology of MPI_Comm_compare()) are always deemed
8058: incompatible.
8060: This function is labeled "Collective," as information about all subdomains
8061: is required on each rank. However, in DM implementations which store all this
8062: information locally, this function may be merely "Logically Collective".
8064: Developer Notes:
8065: Compatibility is assumed to be a symmetric concept; DM A is compatible with DM B
8066: iff B is compatible with A. Thus, this function checks the implementations
8067: of both dm and dm2 (if they are of different types), attempting to determine
8068: compatibility. It is left to DM implementers to ensure that symmetry is
8069: preserved. The simplest way to do this is, when implementing type-specific
8070: logic for this function, is to check for existing logic in the implementation
8071: of other DM types and let *set = PETSC_FALSE if found.
8073: Level: advanced
8075: .seealso: DM, DMDACreateCompatibleDMDA(), DMStagCreateCompatibleDMStag()
8076: @*/
8078: PetscErrorCode DMGetCompatibility(DM dm,DM dm2,PetscBool *compatible,PetscBool *set)
8079: {
8081: PetscMPIInt compareResult;
8082: DMType type,type2;
8083: PetscBool sameType;
8089: /* Declare a DM compatible with itself */
8090: if (dm == dm2) {
8091: *set = PETSC_TRUE;
8092: *compatible = PETSC_TRUE;
8093: return(0);
8094: }
8096: /* Declare a DM incompatible with a DM that lives on an "unequal"
8097: communicator. Note that this does not preclude compatibility with
8098: DMs living on "congruent" or "similar" communicators, but this must be
8099: determined by the implementation-specific logic */
8100: MPI_Comm_compare(PetscObjectComm((PetscObject)dm),PetscObjectComm((PetscObject)dm2),&compareResult);
8101: if (compareResult == MPI_UNEQUAL) {
8102: *set = PETSC_TRUE;
8103: *compatible = PETSC_FALSE;
8104: return(0);
8105: }
8107: /* Pass to the implementation-specific routine, if one exists. */
8108: if (dm->ops->getcompatibility) {
8109: (*dm->ops->getcompatibility)(dm,dm2,compatible,set);
8110: if (*set) return(0);
8111: }
8113: /* If dm and dm2 are of different types, then attempt to check compatibility
8114: with an implementation of this function from dm2 */
8115: DMGetType(dm,&type);
8116: DMGetType(dm2,&type2);
8117: PetscStrcmp(type,type2,&sameType);
8118: if (!sameType && dm2->ops->getcompatibility) {
8119: (*dm2->ops->getcompatibility)(dm2,dm,compatible,set); /* Note argument order */
8120: } else {
8121: *set = PETSC_FALSE;
8122: }
8123: return(0);
8124: }