Actual source code: dm.c
petsc-3.11.4 2019-09-28
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;
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: static PetscErrorCode DMHasCreateInjection_Default(DM dm, PetscBool *flg)
17: {
21: *flg = PETSC_FALSE;
22: return(0);
23: }
25: /*@
26: DMCreate - Creates an empty DM object. The type can then be set with DMSetType().
28: If you never call DMSetType() it will generate an
29: error when you try to use the vector.
31: Collective on MPI_Comm
33: Input Parameter:
34: . comm - The communicator for the DM object
36: Output Parameter:
37: . dm - The DM object
39: Level: beginner
41: .seealso: DMSetType(), DMDA, DMSLICED, DMCOMPOSITE, DMPLEX, DMMOAB, DMNETWORK
42: @*/
43: PetscErrorCode DMCreate(MPI_Comm comm,DM *dm)
44: {
45: DM v;
46: PetscDS ds;
51: *dm = NULL;
52: DMInitializePackage();
54: PetscHeaderCreate(v, DM_CLASSID, "DM", "Distribution Manager", "DM", comm, DMDestroy, DMView);
56: v->setupcalled = PETSC_FALSE;
57: v->setfromoptionscalled = PETSC_FALSE;
58: v->ltogmap = NULL;
59: v->bs = 1;
60: v->coloringtype = IS_COLORING_GLOBAL;
61: PetscSFCreate(comm, &v->sf);
62: PetscSFCreate(comm, &v->defaultSF);
63: v->labels = NULL;
64: v->adjacency[0] = PETSC_FALSE;
65: v->adjacency[1] = PETSC_TRUE;
66: v->depthLabel = NULL;
67: v->defaultSection = NULL;
68: v->defaultGlobalSection = NULL;
69: v->defaultConstraintSection = NULL;
70: v->defaultConstraintMat = NULL;
71: v->L = NULL;
72: v->maxCell = NULL;
73: v->bdtype = NULL;
74: v->dimEmbed = PETSC_DEFAULT;
75: v->dim = PETSC_DETERMINE;
76: {
77: PetscInt i;
78: for (i = 0; i < 10; ++i) {
79: v->nullspaceConstructors[i] = NULL;
80: v->nearnullspaceConstructors[i] = NULL;
81: }
82: }
83: PetscDSCreate(PetscObjectComm((PetscObject) v), &ds);
84: DMSetRegionDS(v, NULL, ds);
85: PetscDSDestroy(&ds);
86: v->dmBC = NULL;
87: v->coarseMesh = NULL;
88: v->outputSequenceNum = -1;
89: v->outputSequenceVal = 0.0;
90: DMSetVecType(v,VECSTANDARD);
91: DMSetMatType(v,MATAIJ);
92: PetscNew(&(v->labels));
93: v->labels->refct = 1;
95: v->ops->hascreateinjection = DMHasCreateInjection_Default;
97: *dm = v;
98: return(0);
99: }
101: /*@
102: DMClone - Creates a DM object with the same topology as the original.
104: Collective on MPI_Comm
106: Input Parameter:
107: . dm - The original DM object
109: Output Parameter:
110: . newdm - The new DM object
112: Level: beginner
114: .keywords: DM, topology, create
115: @*/
116: PetscErrorCode DMClone(DM dm, DM *newdm)
117: {
118: PetscSF sf;
119: Vec coords;
120: void *ctx;
121: PetscInt dim, cdim;
127: DMCreate(PetscObjectComm((PetscObject) dm), newdm);
128: PetscFree((*newdm)->labels);
129: dm->labels->refct++;
130: (*newdm)->labels = dm->labels;
131: (*newdm)->depthLabel = dm->depthLabel;
132: (*newdm)->leveldown = dm->leveldown;
133: (*newdm)->levelup = dm->levelup;
134: DMGetDimension(dm, &dim);
135: DMSetDimension(*newdm, dim);
136: if (dm->ops->clone) {
137: (*dm->ops->clone)(dm, newdm);
138: }
139: (*newdm)->setupcalled = dm->setupcalled;
140: DMGetPointSF(dm, &sf);
141: DMSetPointSF(*newdm, sf);
142: DMGetApplicationContext(dm, &ctx);
143: DMSetApplicationContext(*newdm, ctx);
144: if (dm->coordinateDM) {
145: DM ncdm;
146: PetscSection cs;
147: PetscInt pEnd = -1, pEndMax = -1;
149: DMGetSection(dm->coordinateDM, &cs);
150: if (cs) {PetscSectionGetChart(cs, NULL, &pEnd);}
151: MPI_Allreduce(&pEnd,&pEndMax,1,MPIU_INT,MPI_MAX,PetscObjectComm((PetscObject)dm));
152: if (pEndMax >= 0) {
153: DMClone(dm->coordinateDM, &ncdm);
154: DMCopyDisc(dm->coordinateDM, ncdm);
155: DMSetSection(ncdm, cs);
156: DMSetCoordinateDM(*newdm, ncdm);
157: DMDestroy(&ncdm);
158: }
159: }
160: DMGetCoordinateDim(dm, &cdim);
161: DMSetCoordinateDim(*newdm, cdim);
162: DMGetCoordinatesLocal(dm, &coords);
163: if (coords) {
164: DMSetCoordinatesLocal(*newdm, coords);
165: } else {
166: DMGetCoordinates(dm, &coords);
167: if (coords) {DMSetCoordinates(*newdm, coords);}
168: }
169: {
170: PetscBool isper;
171: const PetscReal *maxCell, *L;
172: const DMBoundaryType *bd;
173: DMGetPeriodicity(dm, &isper, &maxCell, &L, &bd);
174: DMSetPeriodicity(*newdm, isper, maxCell, L, bd);
175: }
176: {
177: PetscBool useCone, useClosure;
179: DMGetAdjacency(dm, PETSC_DEFAULT, &useCone, &useClosure);
180: DMSetAdjacency(*newdm, PETSC_DEFAULT, useCone, useClosure);
181: }
182: return(0);
183: }
185: /*@C
186: DMSetVecType - Sets the type of vector created with DMCreateLocalVector() and DMCreateGlobalVector()
188: Logically Collective on DM
190: Input Parameter:
191: + da - initial distributed array
192: . ctype - the vector type, currently either VECSTANDARD, VECCUDA, or VECVIENNACL
194: Options Database:
195: . -dm_vec_type ctype
197: Level: intermediate
199: .seealso: DMCreate(), DMDestroy(), DM, DMDAInterpolationType, VecType, DMGetVecType(), DMSetMatType(), DMGetMatType()
200: @*/
201: PetscErrorCode DMSetVecType(DM da,VecType ctype)
202: {
207: PetscFree(da->vectype);
208: PetscStrallocpy(ctype,(char**)&da->vectype);
209: return(0);
210: }
212: /*@C
213: DMGetVecType - Gets the type of vector created with DMCreateLocalVector() and DMCreateGlobalVector()
215: Logically Collective on DM
217: Input Parameter:
218: . da - initial distributed array
220: Output Parameter:
221: . ctype - the vector type
223: Level: intermediate
225: .seealso: DMCreate(), DMDestroy(), DM, DMDAInterpolationType, VecType, DMSetMatType(), DMGetMatType(), DMSetVecType()
226: @*/
227: PetscErrorCode DMGetVecType(DM da,VecType *ctype)
228: {
231: *ctype = da->vectype;
232: return(0);
233: }
235: /*@
236: VecGetDM - Gets the DM defining the data layout of the vector
238: Not collective
240: Input Parameter:
241: . v - The Vec
243: Output Parameter:
244: . dm - The DM
246: Level: intermediate
248: .seealso: VecSetDM(), DMGetLocalVector(), DMGetGlobalVector(), DMSetVecType()
249: @*/
250: PetscErrorCode VecGetDM(Vec v, DM *dm)
251: {
257: PetscObjectQuery((PetscObject) v, "__PETSc_dm", (PetscObject*) dm);
258: return(0);
259: }
261: /*@
262: VecSetDM - Sets the DM defining the data layout of the vector.
264: Not collective
266: Input Parameters:
267: + v - The Vec
268: - dm - The DM
270: 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.
272: Level: intermediate
274: .seealso: VecGetDM(), DMGetLocalVector(), DMGetGlobalVector(), DMSetVecType()
275: @*/
276: PetscErrorCode VecSetDM(Vec v, DM dm)
277: {
283: PetscObjectCompose((PetscObject) v, "__PETSc_dm", (PetscObject) dm);
284: return(0);
285: }
287: /*@C
288: DMSetISColoringType - Sets the type of coloring, global or local, that is created by the DM
290: Logically Collective on DM
292: Input Parameters:
293: + dm - the DM context
294: - ctype - the matrix type
296: Options Database:
297: . -dm_is_coloring_type - global or local
299: Level: intermediate
301: .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMCreateMatrix(), DMSetMatrixPreallocateOnly(), MatType, DMGetMatType(),
302: DMGetISColoringType()
303: @*/
304: PetscErrorCode DMSetISColoringType(DM dm,ISColoringType ctype)
305: {
308: dm->coloringtype = ctype;
309: return(0);
310: }
312: /*@C
313: DMGetISColoringType - Gets the type of coloring, global or local, that is created by the DM
315: Logically Collective on DM
317: Input Parameter:
318: . dm - the DM context
320: Output Parameter:
321: . ctype - the matrix type
323: Options Database:
324: . -dm_is_coloring_type - global or local
326: Level: intermediate
328: .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMCreateMatrix(), DMSetMatrixPreallocateOnly(), MatType, DMGetMatType(),
329: DMGetISColoringType()
330: @*/
331: PetscErrorCode DMGetISColoringType(DM dm,ISColoringType *ctype)
332: {
335: *ctype = dm->coloringtype;
336: return(0);
337: }
339: /*@C
340: DMSetMatType - Sets the type of matrix created with DMCreateMatrix()
342: Logically Collective on DM
344: Input Parameters:
345: + dm - the DM context
346: - ctype - the matrix type
348: Options Database:
349: . -dm_mat_type ctype
351: Level: intermediate
353: .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMCreateMatrix(), DMSetMatrixPreallocateOnly(), MatType, DMGetMatType(), DMSetMatType(), DMGetMatType()
354: @*/
355: PetscErrorCode DMSetMatType(DM dm,MatType ctype)
356: {
361: PetscFree(dm->mattype);
362: PetscStrallocpy(ctype,(char**)&dm->mattype);
363: return(0);
364: }
366: /*@C
367: DMGetMatType - Gets the type of matrix created with DMCreateMatrix()
369: Logically Collective on DM
371: Input Parameter:
372: . dm - the DM context
374: Output Parameter:
375: . ctype - the matrix type
377: Options Database:
378: . -dm_mat_type ctype
380: Level: intermediate
382: .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMCreateMatrix(), DMSetMatrixPreallocateOnly(), MatType, DMSetMatType(), DMSetMatType(), DMGetMatType()
383: @*/
384: PetscErrorCode DMGetMatType(DM dm,MatType *ctype)
385: {
388: *ctype = dm->mattype;
389: return(0);
390: }
392: /*@
393: MatGetDM - Gets the DM defining the data layout of the matrix
395: Not collective
397: Input Parameter:
398: . A - The Mat
400: Output Parameter:
401: . dm - The DM
403: Level: intermediate
405: Developer Note: Since the Mat class doesn't know about the DM class the DM object is associated with
406: the Mat through a PetscObjectCompose() operation
408: .seealso: MatSetDM(), DMCreateMatrix(), DMSetMatType()
409: @*/
410: PetscErrorCode MatGetDM(Mat A, DM *dm)
411: {
417: PetscObjectQuery((PetscObject) A, "__PETSc_dm", (PetscObject*) dm);
418: return(0);
419: }
421: /*@
422: MatSetDM - Sets the DM defining the data layout of the matrix
424: Not collective
426: Input Parameters:
427: + A - The Mat
428: - dm - The DM
430: Level: intermediate
432: Developer Note: Since the Mat class doesn't know about the DM class the DM object is associated with
433: the Mat through a PetscObjectCompose() operation
436: .seealso: MatGetDM(), DMCreateMatrix(), DMSetMatType()
437: @*/
438: PetscErrorCode MatSetDM(Mat A, DM dm)
439: {
445: PetscObjectCompose((PetscObject) A, "__PETSc_dm", (PetscObject) dm);
446: return(0);
447: }
449: /*@C
450: DMSetOptionsPrefix - Sets the prefix used for searching for all
451: DM options in the database.
453: Logically Collective on DM
455: Input Parameter:
456: + da - the DM context
457: - prefix - the prefix to prepend to all option names
459: Notes:
460: A hyphen (-) must NOT be given at the beginning of the prefix name.
461: The first character of all runtime options is AUTOMATICALLY the hyphen.
463: Level: advanced
465: .keywords: DM, set, options, prefix, database
467: .seealso: DMSetFromOptions()
468: @*/
469: PetscErrorCode DMSetOptionsPrefix(DM dm,const char prefix[])
470: {
475: PetscObjectSetOptionsPrefix((PetscObject)dm,prefix);
476: if (dm->sf) {
477: PetscObjectSetOptionsPrefix((PetscObject)dm->sf,prefix);
478: }
479: if (dm->defaultSF) {
480: PetscObjectSetOptionsPrefix((PetscObject)dm->defaultSF,prefix);
481: }
482: return(0);
483: }
485: /*@C
486: DMAppendOptionsPrefix - Appends to the prefix used for searching for all
487: DM options in the database.
489: Logically Collective on DM
491: Input Parameters:
492: + dm - the DM context
493: - prefix - the prefix string to prepend to all DM option requests
495: Notes:
496: A hyphen (-) must NOT be given at the beginning of the prefix name.
497: The first character of all runtime options is AUTOMATICALLY the hyphen.
499: Level: advanced
501: .keywords: DM, append, options, prefix, database
503: .seealso: DMSetOptionsPrefix(), DMGetOptionsPrefix()
504: @*/
505: PetscErrorCode DMAppendOptionsPrefix(DM dm,const char prefix[])
506: {
511: PetscObjectAppendOptionsPrefix((PetscObject)dm,prefix);
512: return(0);
513: }
515: /*@C
516: DMGetOptionsPrefix - Gets the prefix used for searching for all
517: DM options in the database.
519: Not Collective
521: Input Parameters:
522: . dm - the DM context
524: Output Parameters:
525: . prefix - pointer to the prefix string used is returned
527: Notes:
528: On the fortran side, the user should pass in a string 'prefix' of
529: sufficient length to hold the prefix.
531: Level: advanced
533: .keywords: DM, set, options, prefix, database
535: .seealso: DMSetOptionsPrefix(), DMAppendOptionsPrefix()
536: @*/
537: PetscErrorCode DMGetOptionsPrefix(DM dm,const char *prefix[])
538: {
543: PetscObjectGetOptionsPrefix((PetscObject)dm,prefix);
544: return(0);
545: }
547: static PetscErrorCode DMCountNonCyclicReferences(DM dm, PetscBool recurseCoarse, PetscBool recurseFine, PetscInt *ncrefct)
548: {
549: PetscInt i, refct = ((PetscObject) dm)->refct;
550: DMNamedVecLink nlink;
554: *ncrefct = 0;
555: /* count all the circular references of DM and its contained Vecs */
556: for (i=0; i<DM_MAX_WORK_VECTORS; i++) {
557: if (dm->localin[i]) refct--;
558: if (dm->globalin[i]) refct--;
559: }
560: for (nlink=dm->namedglobal; nlink; nlink=nlink->next) refct--;
561: for (nlink=dm->namedlocal; nlink; nlink=nlink->next) refct--;
562: if (dm->x) {
563: DM obj;
564: VecGetDM(dm->x, &obj);
565: if (obj == dm) refct--;
566: }
567: if (dm->coarseMesh && dm->coarseMesh->fineMesh == dm) {
568: refct--;
569: if (recurseCoarse) {
570: PetscInt coarseCount;
572: DMCountNonCyclicReferences(dm->coarseMesh, PETSC_TRUE, PETSC_FALSE,&coarseCount);
573: refct += coarseCount;
574: }
575: }
576: if (dm->fineMesh && dm->fineMesh->coarseMesh == dm) {
577: refct--;
578: if (recurseFine) {
579: PetscInt fineCount;
581: DMCountNonCyclicReferences(dm->fineMesh, PETSC_FALSE, PETSC_TRUE,&fineCount);
582: refct += fineCount;
583: }
584: }
585: *ncrefct = refct;
586: return(0);
587: }
589: PetscErrorCode DMDestroyLabelLinkList(DM dm)
590: {
594: if (!--(dm->labels->refct)) {
595: DMLabelLink next = dm->labels->next;
597: /* destroy the labels */
598: while (next) {
599: DMLabelLink tmp = next->next;
601: DMLabelDestroy(&next->label);
602: PetscFree(next);
603: next = tmp;
604: }
605: PetscFree(dm->labels);
606: }
607: return(0);
608: }
610: /*@
611: DMDestroy - Destroys a vector packer or DM.
613: Collective on DM
615: Input Parameter:
616: . dm - the DM object to destroy
618: Level: developer
620: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix()
622: @*/
623: PetscErrorCode DMDestroy(DM *dm)
624: {
625: PetscInt i, cnt;
626: DMNamedVecLink nlink,nnext;
630: if (!*dm) return(0);
633: /* count all non-cyclic references in the doubly-linked list of coarse<->fine meshes */
634: DMCountNonCyclicReferences(*dm,PETSC_TRUE,PETSC_TRUE,&cnt);
635: --((PetscObject)(*dm))->refct;
636: if (--cnt > 0) {*dm = 0; return(0);}
637: /*
638: Need this test because the dm references the vectors that
639: reference the dm, so destroying the dm calls destroy on the
640: vectors that cause another destroy on the dm
641: */
642: if (((PetscObject)(*dm))->refct < 0) return(0);
643: ((PetscObject) (*dm))->refct = 0;
644: for (i=0; i<DM_MAX_WORK_VECTORS; i++) {
645: if ((*dm)->localout[i]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Destroying a DM that has a local vector obtained with DMGetLocalVector()");
646: VecDestroy(&(*dm)->localin[i]);
647: }
648: nnext=(*dm)->namedglobal;
649: (*dm)->namedglobal = NULL;
650: for (nlink=nnext; nlink; nlink=nnext) { /* Destroy the named vectors */
651: nnext = nlink->next;
652: if (nlink->status != DMVEC_STATUS_IN) SETERRQ1(((PetscObject)*dm)->comm,PETSC_ERR_ARG_WRONGSTATE,"DM still has Vec named '%s' checked out",nlink->name);
653: PetscFree(nlink->name);
654: VecDestroy(&nlink->X);
655: PetscFree(nlink);
656: }
657: nnext=(*dm)->namedlocal;
658: (*dm)->namedlocal = NULL;
659: for (nlink=nnext; nlink; nlink=nnext) { /* Destroy the named local vectors */
660: nnext = nlink->next;
661: if (nlink->status != DMVEC_STATUS_IN) SETERRQ1(((PetscObject)*dm)->comm,PETSC_ERR_ARG_WRONGSTATE,"DM still has Vec named '%s' checked out",nlink->name);
662: PetscFree(nlink->name);
663: VecDestroy(&nlink->X);
664: PetscFree(nlink);
665: }
667: /* Destroy the list of hooks */
668: {
669: DMCoarsenHookLink link,next;
670: for (link=(*dm)->coarsenhook; link; link=next) {
671: next = link->next;
672: PetscFree(link);
673: }
674: (*dm)->coarsenhook = NULL;
675: }
676: {
677: DMRefineHookLink link,next;
678: for (link=(*dm)->refinehook; link; link=next) {
679: next = link->next;
680: PetscFree(link);
681: }
682: (*dm)->refinehook = NULL;
683: }
684: {
685: DMSubDomainHookLink link,next;
686: for (link=(*dm)->subdomainhook; link; link=next) {
687: next = link->next;
688: PetscFree(link);
689: }
690: (*dm)->subdomainhook = NULL;
691: }
692: {
693: DMGlobalToLocalHookLink link,next;
694: for (link=(*dm)->gtolhook; link; link=next) {
695: next = link->next;
696: PetscFree(link);
697: }
698: (*dm)->gtolhook = NULL;
699: }
700: {
701: DMLocalToGlobalHookLink link,next;
702: for (link=(*dm)->ltoghook; link; link=next) {
703: next = link->next;
704: PetscFree(link);
705: }
706: (*dm)->ltoghook = NULL;
707: }
708: /* Destroy the work arrays */
709: {
710: DMWorkLink link,next;
711: if ((*dm)->workout) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Work array still checked out");
712: for (link=(*dm)->workin; link; link=next) {
713: next = link->next;
714: PetscFree(link->mem);
715: PetscFree(link);
716: }
717: (*dm)->workin = NULL;
718: }
719: if (!--((*dm)->labels->refct)) {
720: DMLabelLink next = (*dm)->labels->next;
722: /* destroy the labels */
723: while (next) {
724: DMLabelLink tmp = next->next;
726: DMLabelDestroy(&next->label);
727: PetscFree(next);
728: next = tmp;
729: }
730: PetscFree((*dm)->labels);
731: }
732: DMClearFields(*dm);
733: {
734: DMBoundary next = (*dm)->boundary;
735: while (next) {
736: DMBoundary b = next;
738: next = b->next;
739: PetscFree(b);
740: }
741: }
743: PetscObjectDestroy(&(*dm)->dmksp);
744: PetscObjectDestroy(&(*dm)->dmsnes);
745: PetscObjectDestroy(&(*dm)->dmts);
747: if ((*dm)->ctx && (*dm)->ctxdestroy) {
748: (*(*dm)->ctxdestroy)(&(*dm)->ctx);
749: }
750: VecDestroy(&(*dm)->x);
751: MatFDColoringDestroy(&(*dm)->fd);
752: DMClearGlobalVectors(*dm);
753: ISLocalToGlobalMappingDestroy(&(*dm)->ltogmap);
754: PetscFree((*dm)->vectype);
755: PetscFree((*dm)->mattype);
757: PetscSectionDestroy(&(*dm)->defaultSection);
758: PetscSectionDestroy(&(*dm)->defaultGlobalSection);
759: PetscLayoutDestroy(&(*dm)->map);
760: PetscSectionDestroy(&(*dm)->defaultConstraintSection);
761: MatDestroy(&(*dm)->defaultConstraintMat);
762: PetscSFDestroy(&(*dm)->sf);
763: PetscSFDestroy(&(*dm)->defaultSF);
764: if ((*dm)->useNatural) {
765: if ((*dm)->sfNatural) {
766: PetscSFDestroy(&(*dm)->sfNatural);
767: }
768: PetscObjectDereference((PetscObject) (*dm)->sfMigration);
769: }
770: if ((*dm)->coarseMesh && (*dm)->coarseMesh->fineMesh == *dm) {
771: DMSetFineDM((*dm)->coarseMesh,NULL);
772: }
773: DMDestroy(&(*dm)->coarseMesh);
774: if ((*dm)->fineMesh && (*dm)->fineMesh->coarseMesh == *dm) {
775: DMSetCoarseDM((*dm)->fineMesh,NULL);
776: }
777: DMDestroy(&(*dm)->fineMesh);
778: DMFieldDestroy(&(*dm)->coordinateField);
779: DMDestroy(&(*dm)->coordinateDM);
780: VecDestroy(&(*dm)->coordinates);
781: VecDestroy(&(*dm)->coordinatesLocal);
782: PetscFree3((*dm)->L,(*dm)->maxCell,(*dm)->bdtype);
784: DMClearDS(*dm);
785: DMDestroy(&(*dm)->dmBC);
786: /* if memory was published with SAWs then destroy it */
787: PetscObjectSAWsViewOff((PetscObject)*dm);
789: if ((*dm)->ops->destroy) {
790: (*(*dm)->ops->destroy)(*dm);
791: }
792: /* We do not destroy (*dm)->data here so that we can reference count backend objects */
793: PetscHeaderDestroy(dm);
794: return(0);
795: }
797: /*@
798: DMSetUp - sets up the data structures inside a DM object
800: Collective on DM
802: Input Parameter:
803: . dm - the DM object to setup
805: Level: developer
807: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix()
809: @*/
810: PetscErrorCode DMSetUp(DM dm)
811: {
816: if (dm->setupcalled) return(0);
817: if (dm->ops->setup) {
818: (*dm->ops->setup)(dm);
819: }
820: dm->setupcalled = PETSC_TRUE;
821: return(0);
822: }
824: /*@
825: DMSetFromOptions - sets parameters in a DM from the options database
827: Collective on DM
829: Input Parameter:
830: . dm - the DM object to set options for
832: Options Database:
833: + -dm_preallocate_only - Only preallocate the matrix for DMCreateMatrix(), but do not fill it with zeros
834: . -dm_vec_type <type> - type of vector to create inside DM
835: . -dm_mat_type <type> - type of matrix to create inside DM
836: - -dm_is_coloring_type - <global or local>
838: Level: developer
840: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix()
842: @*/
843: PetscErrorCode DMSetFromOptions(DM dm)
844: {
845: char typeName[256];
846: PetscBool flg;
851: dm->setfromoptionscalled = PETSC_TRUE;
852: if (dm->sf) {PetscSFSetFromOptions(dm->sf);}
853: if (dm->defaultSF) {PetscSFSetFromOptions(dm->defaultSF);}
854: PetscObjectOptionsBegin((PetscObject)dm);
855: PetscOptionsBool("-dm_preallocate_only","only preallocate matrix, but do not set column indices","DMSetMatrixPreallocateOnly",dm->prealloc_only,&dm->prealloc_only,NULL);
856: PetscOptionsFList("-dm_vec_type","Vector type used for created vectors","DMSetVecType",VecList,dm->vectype,typeName,256,&flg);
857: if (flg) {
858: DMSetVecType(dm,typeName);
859: }
860: PetscOptionsFList("-dm_mat_type","Matrix type used for created matrices","DMSetMatType",MatList,dm->mattype ? dm->mattype : typeName,typeName,sizeof(typeName),&flg);
861: if (flg) {
862: DMSetMatType(dm,typeName);
863: }
864: PetscOptionsEnum("-dm_is_coloring_type","Global or local coloring of Jacobian","DMSetISColoringType",ISColoringTypes,(PetscEnum)dm->coloringtype,(PetscEnum*)&dm->coloringtype,NULL);
865: if (dm->ops->setfromoptions) {
866: (*dm->ops->setfromoptions)(PetscOptionsObject,dm);
867: }
868: /* process any options handlers added with PetscObjectAddOptionsHandler() */
869: PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject) dm);
870: PetscOptionsEnd();
871: return(0);
872: }
874: /*@C
875: DMView - Views a DM
877: Collective on DM
879: Input Parameter:
880: + dm - the DM object to view
881: - v - the viewer
883: Level: beginner
885: .seealso DMDestroy(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix()
887: @*/
888: PetscErrorCode DMView(DM dm,PetscViewer v)
889: {
890: PetscErrorCode ierr;
891: PetscBool isbinary;
892: PetscMPIInt size;
893: PetscViewerFormat format;
897: if (!v) {
898: PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)dm),&v);
899: }
900: PetscViewerCheckWritable(v);
901: PetscViewerGetFormat(v,&format);
902: MPI_Comm_size(PetscObjectComm((PetscObject)dm),&size);
903: if (size == 1 && format == PETSC_VIEWER_LOAD_BALANCE) return(0);
904: PetscObjectPrintClassNamePrefixType((PetscObject)dm,v);
905: PetscObjectTypeCompare((PetscObject)v,PETSCVIEWERBINARY,&isbinary);
906: if (isbinary) {
907: PetscInt classid = DM_FILE_CLASSID;
908: char type[256];
910: PetscViewerBinaryWrite(v,&classid,1,PETSC_INT,PETSC_FALSE);
911: PetscStrncpy(type,((PetscObject)dm)->type_name,256);
912: PetscViewerBinaryWrite(v,type,256,PETSC_CHAR,PETSC_FALSE);
913: }
914: if (dm->ops->view) {
915: (*dm->ops->view)(dm,v);
916: }
917: return(0);
918: }
920: /*@
921: DMCreateGlobalVector - Creates a global vector from a DM object
923: Collective on DM
925: Input Parameter:
926: . dm - the DM object
928: Output Parameter:
929: . vec - the global vector
931: Level: beginner
933: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix()
935: @*/
936: PetscErrorCode DMCreateGlobalVector(DM dm,Vec *vec)
937: {
942: (*dm->ops->createglobalvector)(dm,vec);
943: return(0);
944: }
946: /*@
947: DMCreateLocalVector - Creates a local vector from a DM object
949: Not Collective
951: Input Parameter:
952: . dm - the DM object
954: Output Parameter:
955: . vec - the local vector
957: Level: beginner
959: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix()
961: @*/
962: PetscErrorCode DMCreateLocalVector(DM dm,Vec *vec)
963: {
968: (*dm->ops->createlocalvector)(dm,vec);
969: return(0);
970: }
972: /*@
973: DMGetLocalToGlobalMapping - Accesses the local-to-global mapping in a DM.
975: Collective on DM
977: Input Parameter:
978: . dm - the DM that provides the mapping
980: Output Parameter:
981: . ltog - the mapping
983: Level: intermediate
985: Notes:
986: This mapping can then be used by VecSetLocalToGlobalMapping() or
987: MatSetLocalToGlobalMapping().
989: .seealso: DMCreateLocalVector()
990: @*/
991: PetscErrorCode DMGetLocalToGlobalMapping(DM dm,ISLocalToGlobalMapping *ltog)
992: {
993: PetscInt bs = -1, bsLocal[2], bsMinMax[2];
999: if (!dm->ltogmap) {
1000: PetscSection section, sectionGlobal;
1002: DMGetSection(dm, §ion);
1003: if (section) {
1004: const PetscInt *cdofs;
1005: PetscInt *ltog;
1006: PetscInt pStart, pEnd, n, p, k, l;
1008: DMGetGlobalSection(dm, §ionGlobal);
1009: PetscSectionGetChart(section, &pStart, &pEnd);
1010: PetscSectionGetStorageSize(section, &n);
1011: PetscMalloc1(n, <og); /* We want the local+overlap size */
1012: for (p = pStart, l = 0; p < pEnd; ++p) {
1013: PetscInt bdof, cdof, dof, off, c, cind = 0;
1015: /* Should probably use constrained dofs */
1016: PetscSectionGetDof(section, p, &dof);
1017: PetscSectionGetConstraintDof(section, p, &cdof);
1018: PetscSectionGetConstraintIndices(section, p, &cdofs);
1019: PetscSectionGetOffset(sectionGlobal, p, &off);
1020: /* If you have dofs, and constraints, and they are unequal, we set the blocksize to 1 */
1021: bdof = cdof && (dof-cdof) ? 1 : dof;
1022: if (dof) {
1023: if (bs < 0) {bs = bdof;}
1024: else if (bs != bdof) {bs = 1;}
1025: }
1026: for (c = 0; c < dof; ++c, ++l) {
1027: if ((cind < cdof) && (c == cdofs[cind])) ltog[l] = off < 0 ? off-c : off+c;
1028: else ltog[l] = (off < 0 ? -(off+1) : off) + c;
1029: }
1030: }
1031: /* Must have same blocksize on all procs (some might have no points) */
1032: bsLocal[0] = bs < 0 ? PETSC_MAX_INT : bs; bsLocal[1] = bs;
1033: PetscGlobalMinMaxInt(PetscObjectComm((PetscObject) dm), bsLocal, bsMinMax);
1034: if (bsMinMax[0] != bsMinMax[1]) {bs = 1;}
1035: else {bs = bsMinMax[0];}
1036: bs = bs < 0 ? 1 : bs;
1037: /* Must reduce indices by blocksize */
1038: if (bs > 1) {
1039: for (l = 0, k = 0; l < n; l += bs, ++k) ltog[k] = ltog[l]/bs;
1040: n /= bs;
1041: }
1042: ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)dm), bs, n, ltog, PETSC_OWN_POINTER, &dm->ltogmap);
1043: PetscLogObjectParent((PetscObject)dm, (PetscObject)dm->ltogmap);
1044: } else {
1045: if (!dm->ops->getlocaltoglobalmapping) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM can not create LocalToGlobalMapping");
1046: (*dm->ops->getlocaltoglobalmapping)(dm);
1047: }
1048: }
1049: *ltog = dm->ltogmap;
1050: return(0);
1051: }
1053: /*@
1054: DMGetBlockSize - Gets the inherent block size associated with a DM
1056: Not Collective
1058: Input Parameter:
1059: . dm - the DM with block structure
1061: Output Parameter:
1062: . bs - the block size, 1 implies no exploitable block structure
1064: Level: intermediate
1066: .seealso: ISCreateBlock(), VecSetBlockSize(), MatSetBlockSize(), DMGetLocalToGlobalMapping()
1067: @*/
1068: PetscErrorCode DMGetBlockSize(DM dm,PetscInt *bs)
1069: {
1073: if (dm->bs < 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"DM does not have enough information to provide a block size yet");
1074: *bs = dm->bs;
1075: return(0);
1076: }
1078: /*@
1079: DMCreateInterpolation - Gets interpolation matrix between two DM objects
1081: Collective on DM
1083: Input Parameter:
1084: + dm1 - the DM object
1085: - dm2 - the second, finer DM object
1087: Output Parameter:
1088: + mat - the interpolation
1089: - vec - the scaling (optional)
1091: Level: developer
1093: Notes:
1094: For DMDA objects this only works for "uniform refinement", that is the refined mesh was obtained DMRefine() or the coarse mesh was obtained by
1095: DMCoarsen(). The coordinates set into the DMDA are completely ignored in computing the interpolation.
1097: For DMDA objects you can use this interpolation (more precisely the interpolation from the DMGetCoordinateDM()) to interpolate the mesh coordinate vectors
1098: EXCEPT in the periodic case where it does not make sense since the coordinate vectors are not periodic.
1101: .seealso DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateColoring(), DMCreateMatrix(), DMRefine(), DMCoarsen(), DMCreateRestriction()
1103: @*/
1104: PetscErrorCode DMCreateInterpolation(DM dm1,DM dm2,Mat *mat,Vec *vec)
1105: {
1112: if (!dm1->ops->createinterpolation) SETERRQ1(PetscObjectComm((PetscObject)dm1),PETSC_ERR_SUP,"DMCreateInterpolation not implemented for type %s",((PetscObject)dm1)->type_name);
1113: PetscLogEventBegin(DM_CreateInterpolation,dm1,dm2,0,0);
1114: (*dm1->ops->createinterpolation)(dm1,dm2,mat,vec);
1115: PetscLogEventEnd(DM_CreateInterpolation,dm1,dm2,0,0);
1116: return(0);
1117: }
1119: /*@
1120: DMCreateRestriction - Gets restriction matrix between two DM objects
1122: Collective on DM
1124: Input Parameter:
1125: + dm1 - the DM object
1126: - dm2 - the second, finer DM object
1128: Output Parameter:
1129: . mat - the restriction
1132: Level: developer
1134: Notes:
1135: For DMDA objects this only works for "uniform refinement", that is the refined mesh was obtained DMRefine() or the coarse mesh was obtained by
1136: DMCoarsen(). The coordinates set into the DMDA are completely ignored in computing the interpolation.
1139: .seealso DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateColoring(), DMCreateMatrix(), DMRefine(), DMCoarsen(), DMCreateInterpolation()
1141: @*/
1142: PetscErrorCode DMCreateRestriction(DM dm1,DM dm2,Mat *mat)
1143: {
1149: PetscLogEventBegin(DM_CreateRestriction,dm1,dm2,0,0);
1150: if (!dm1->ops->createrestriction) SETERRQ(PetscObjectComm((PetscObject)dm1),PETSC_ERR_SUP,"DMCreateRestriction not implemented for this type");
1151: (*dm1->ops->createrestriction)(dm1,dm2,mat);
1152: PetscLogEventEnd(DM_CreateRestriction,dm1,dm2,0,0);
1153: return(0);
1154: }
1156: /*@
1157: DMCreateInjection - Gets injection matrix between two DM objects
1159: Collective on DM
1161: Input Parameter:
1162: + dm1 - the DM object
1163: - dm2 - the second, finer DM object
1165: Output Parameter:
1166: . mat - the injection
1168: Level: developer
1170: Notes:
1171: For DMDA objects this only works for "uniform refinement", that is the refined mesh was obtained DMRefine() or the coarse mesh was obtained by
1172: DMCoarsen(). The coordinates set into the DMDA are completely ignored in computing the injection.
1174: .seealso DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateColoring(), DMCreateMatrix(), DMCreateInterpolation()
1176: @*/
1177: PetscErrorCode DMCreateInjection(DM dm1,DM dm2,Mat *mat)
1178: {
1184: if (!dm1->ops->getinjection) SETERRQ(PetscObjectComm((PetscObject)dm1),PETSC_ERR_SUP,"DMCreateInjection not implemented for this type");
1185: (*dm1->ops->getinjection)(dm1,dm2,mat);
1186: return(0);
1187: }
1189: /*@
1190: DMCreateMassMatrix - Gets mass matrix between two DM objects, M_ij = \int \phi_i \psi_j
1192: Collective on DM
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: {
1212: (*dm1->ops->createmassmatrix)(dm1, dm2, mat);
1213: return(0);
1214: }
1216: /*@
1217: DMCreateColoring - Gets coloring for a DM
1219: Collective on DM
1221: Input Parameter:
1222: + dm - the DM object
1223: - ctype - IS_COLORING_LOCAL or IS_COLORING_GLOBAL
1225: Output Parameter:
1226: . coloring - the coloring
1228: Level: developer
1230: .seealso DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateMatrix(), DMSetMatType()
1232: @*/
1233: PetscErrorCode DMCreateColoring(DM dm,ISColoringType ctype,ISColoring *coloring)
1234: {
1239: if (!dm->ops->getcoloring) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"No coloring for this type of DM yet");
1240: (*dm->ops->getcoloring)(dm,ctype,coloring);
1241: return(0);
1242: }
1244: /*@
1245: DMCreateMatrix - Gets empty Jacobian for a DM
1247: Collective on DM
1249: Input Parameter:
1250: . dm - the DM object
1252: Output Parameter:
1253: . mat - the empty Jacobian
1255: Level: beginner
1257: Notes:
1258: This properly preallocates the number of nonzeros in the sparse matrix so you
1259: do not need to do it yourself.
1261: By default it also sets the nonzero structure and puts in the zero entries. To prevent setting
1262: the nonzero pattern call DMSetMatrixPreallocateOnly()
1264: For structured grid problems, when you call MatView() on this matrix it is displayed using the global natural ordering, NOT in the ordering used
1265: internally by PETSc.
1267: For structured grid problems, in general it is easiest to use MatSetValuesStencil() or MatSetValuesLocal() to put values into the matrix because MatSetValues() requires
1268: the indices for the global numbering for DMDAs which is complicated.
1270: .seealso DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMSetMatType()
1272: @*/
1273: PetscErrorCode DMCreateMatrix(DM dm,Mat *mat)
1274: {
1279: MatInitializePackage();
1282: (*dm->ops->creatematrix)(dm,mat);
1283: /* Handle nullspace and near nullspace */
1284: if (dm->Nf) {
1285: MatNullSpace nullSpace;
1286: PetscInt Nf;
1288: DMGetNumFields(dm, &Nf);
1289: if (Nf == 1) {
1290: if (dm->nullspaceConstructors[0]) {
1291: (*dm->nullspaceConstructors[0])(dm, 0, &nullSpace);
1292: MatSetNullSpace(*mat, nullSpace);
1293: MatNullSpaceDestroy(&nullSpace);
1294: }
1295: if (dm->nearnullspaceConstructors[0]) {
1296: (*dm->nearnullspaceConstructors[0])(dm, 0, &nullSpace);
1297: MatSetNearNullSpace(*mat, nullSpace);
1298: MatNullSpaceDestroy(&nullSpace);
1299: }
1300: }
1301: }
1302: return(0);
1303: }
1305: /*@
1306: DMSetMatrixPreallocateOnly - When DMCreateMatrix() is called the matrix will be properly
1307: preallocated but the nonzero structure and zero values will not be set.
1309: Logically Collective on DM
1311: Input Parameter:
1312: + dm - the DM
1313: - only - PETSC_TRUE if only want preallocation
1315: Level: developer
1316: .seealso DMCreateMatrix(), DMSetMatrixStructureOnly()
1317: @*/
1318: PetscErrorCode DMSetMatrixPreallocateOnly(DM dm, PetscBool only)
1319: {
1322: dm->prealloc_only = only;
1323: return(0);
1324: }
1326: /*@
1327: DMSetMatrixStructureOnly - When DMCreateMatrix() is called, the matrix structure will be created
1328: but the array for values will not be allocated.
1330: Logically Collective on DM
1332: Input Parameter:
1333: + dm - the DM
1334: - only - PETSC_TRUE if only want matrix stucture
1336: Level: developer
1337: .seealso DMCreateMatrix(), DMSetMatrixPreallocateOnly()
1338: @*/
1339: PetscErrorCode DMSetMatrixStructureOnly(DM dm, PetscBool only)
1340: {
1343: dm->structure_only = only;
1344: return(0);
1345: }
1347: /*@C
1348: DMGetWorkArray - Gets a work array guaranteed to be at least the input size, restore with DMRestoreWorkArray()
1350: Not Collective
1352: Input Parameters:
1353: + dm - the DM object
1354: . count - The minium size
1355: - dtype - MPI data type, often MPIU_REAL, MPIU_SCALAR, MPIU_INT)
1357: Output Parameter:
1358: . array - the work array
1360: Level: developer
1362: .seealso DMDestroy(), DMCreate()
1363: @*/
1364: PetscErrorCode DMGetWorkArray(DM dm,PetscInt count,MPI_Datatype dtype,void *mem)
1365: {
1367: DMWorkLink link;
1368: PetscMPIInt dsize;
1373: if (dm->workin) {
1374: link = dm->workin;
1375: dm->workin = dm->workin->next;
1376: } else {
1377: PetscNewLog(dm,&link);
1378: }
1379: MPI_Type_size(dtype,&dsize);
1380: if (((size_t)dsize*count) > link->bytes) {
1381: PetscFree(link->mem);
1382: PetscMalloc(dsize*count,&link->mem);
1383: link->bytes = dsize*count;
1384: }
1385: link->next = dm->workout;
1386: dm->workout = link;
1387: *(void**)mem = link->mem;
1388: return(0);
1389: }
1391: /*@C
1392: DMRestoreWorkArray - Restores a work array guaranteed to be at least the input size, restore with DMRestoreWorkArray()
1394: Not Collective
1396: Input Parameters:
1397: + dm - the DM object
1398: . count - The minium size
1399: - dtype - MPI data type, often MPIU_REAL, MPIU_SCALAR, MPIU_INT
1401: Output Parameter:
1402: . array - the work array
1404: Level: developer
1406: Developer Notes:
1407: count and dtype are ignored, they are only needed for DMGetWorkArray()
1408: .seealso DMDestroy(), DMCreate()
1409: @*/
1410: PetscErrorCode DMRestoreWorkArray(DM dm,PetscInt count,MPI_Datatype dtype,void *mem)
1411: {
1412: DMWorkLink *p,link;
1417: for (p=&dm->workout; (link=*p); p=&link->next) {
1418: if (link->mem == *(void**)mem) {
1419: *p = link->next;
1420: link->next = dm->workin;
1421: dm->workin = link;
1422: *(void**)mem = NULL;
1423: return(0);
1424: }
1425: }
1426: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Array was not checked out");
1427: }
1429: PetscErrorCode DMSetNullSpaceConstructor(DM dm, PetscInt field, PetscErrorCode (*nullsp)(DM dm, PetscInt field, MatNullSpace *nullSpace))
1430: {
1433: if (field >= 10) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle %d >= 10 fields", field);
1434: dm->nullspaceConstructors[field] = nullsp;
1435: return(0);
1436: }
1438: PetscErrorCode DMGetNullSpaceConstructor(DM dm, PetscInt field, PetscErrorCode (**nullsp)(DM dm, PetscInt field, MatNullSpace *nullSpace))
1439: {
1443: if (field >= 10) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle %d >= 10 fields", field);
1444: *nullsp = dm->nullspaceConstructors[field];
1445: return(0);
1446: }
1448: PetscErrorCode DMSetNearNullSpaceConstructor(DM dm, PetscInt field, PetscErrorCode (*nullsp)(DM dm, PetscInt field, MatNullSpace *nullSpace))
1449: {
1452: if (field >= 10) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle %d >= 10 fields", field);
1453: dm->nearnullspaceConstructors[field] = nullsp;
1454: return(0);
1455: }
1457: PetscErrorCode DMGetNearNullSpaceConstructor(DM dm, PetscInt field, PetscErrorCode (**nullsp)(DM dm, PetscInt field, MatNullSpace *nullSpace))
1458: {
1462: if (field >= 10) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle %d >= 10 fields", field);
1463: *nullsp = dm->nearnullspaceConstructors[field];
1464: return(0);
1465: }
1467: /*@C
1468: DMCreateFieldIS - Creates a set of IS objects with the global indices of dofs for each field
1470: Not collective
1472: Input Parameter:
1473: . dm - the DM object
1475: Output Parameters:
1476: + numFields - The number of fields (or NULL if not requested)
1477: . fieldNames - The name for each field (or NULL if not requested)
1478: - fields - The global indices for each field (or NULL if not requested)
1480: Level: intermediate
1482: Notes:
1483: The user is responsible for freeing all requested arrays. In particular, every entry of names should be freed with
1484: PetscFree(), every entry of fields should be destroyed with ISDestroy(), and both arrays should be freed with
1485: PetscFree().
1487: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix()
1488: @*/
1489: PetscErrorCode DMCreateFieldIS(DM dm, PetscInt *numFields, char ***fieldNames, IS **fields)
1490: {
1491: PetscSection section, sectionGlobal;
1496: if (numFields) {
1498: *numFields = 0;
1499: }
1500: if (fieldNames) {
1502: *fieldNames = NULL;
1503: }
1504: if (fields) {
1506: *fields = NULL;
1507: }
1508: DMGetSection(dm, §ion);
1509: if (section) {
1510: PetscInt *fieldSizes, *fieldNc, **fieldIndices;
1511: PetscInt nF, f, pStart, pEnd, p;
1513: DMGetGlobalSection(dm, §ionGlobal);
1514: PetscSectionGetNumFields(section, &nF);
1515: PetscMalloc3(nF,&fieldSizes,nF,&fieldNc,nF,&fieldIndices);
1516: PetscSectionGetChart(sectionGlobal, &pStart, &pEnd);
1517: for (f = 0; f < nF; ++f) {
1518: fieldSizes[f] = 0;
1519: PetscSectionGetFieldComponents(section, f, &fieldNc[f]);
1520: }
1521: for (p = pStart; p < pEnd; ++p) {
1522: PetscInt gdof;
1524: PetscSectionGetDof(sectionGlobal, p, &gdof);
1525: if (gdof > 0) {
1526: for (f = 0; f < nF; ++f) {
1527: PetscInt fdof, fcdof, fpdof;
1529: PetscSectionGetFieldDof(section, p, f, &fdof);
1530: PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);
1531: fpdof = fdof-fcdof;
1532: if (fpdof && fpdof != fieldNc[f]) {
1533: /* Layout does not admit a pointwise block size */
1534: fieldNc[f] = 1;
1535: }
1536: fieldSizes[f] += fpdof;
1537: }
1538: }
1539: }
1540: for (f = 0; f < nF; ++f) {
1541: PetscMalloc1(fieldSizes[f], &fieldIndices[f]);
1542: fieldSizes[f] = 0;
1543: }
1544: for (p = pStart; p < pEnd; ++p) {
1545: PetscInt gdof, goff;
1547: PetscSectionGetDof(sectionGlobal, p, &gdof);
1548: if (gdof > 0) {
1549: PetscSectionGetOffset(sectionGlobal, p, &goff);
1550: for (f = 0; f < nF; ++f) {
1551: PetscInt fdof, fcdof, fc;
1553: PetscSectionGetFieldDof(section, p, f, &fdof);
1554: PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);
1555: for (fc = 0; fc < fdof-fcdof; ++fc, ++fieldSizes[f]) {
1556: fieldIndices[f][fieldSizes[f]] = goff++;
1557: }
1558: }
1559: }
1560: }
1561: if (numFields) *numFields = nF;
1562: if (fieldNames) {
1563: PetscMalloc1(nF, fieldNames);
1564: for (f = 0; f < nF; ++f) {
1565: const char *fieldName;
1567: PetscSectionGetFieldName(section, f, &fieldName);
1568: PetscStrallocpy(fieldName, (char**) &(*fieldNames)[f]);
1569: }
1570: }
1571: if (fields) {
1572: PetscMalloc1(nF, fields);
1573: for (f = 0; f < nF; ++f) {
1574: PetscInt bs, in[2], out[2];
1576: ISCreateGeneral(PetscObjectComm((PetscObject)dm), fieldSizes[f], fieldIndices[f], PETSC_OWN_POINTER, &(*fields)[f]);
1577: in[0] = -fieldNc[f];
1578: in[1] = fieldNc[f];
1579: MPIU_Allreduce(in, out, 2, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)dm));
1580: bs = (-out[0] == out[1]) ? out[1] : 1;
1581: ISSetBlockSize((*fields)[f], bs);
1582: }
1583: }
1584: PetscFree3(fieldSizes,fieldNc,fieldIndices);
1585: } else if (dm->ops->createfieldis) {
1586: (*dm->ops->createfieldis)(dm, numFields, fieldNames, fields);
1587: }
1588: return(0);
1589: }
1592: /*@C
1593: DMCreateFieldDecomposition - Returns a list of IS objects defining a decomposition of a problem into subproblems
1594: corresponding to different fields: each IS contains the global indices of the dofs of the
1595: corresponding field. The optional list of DMs define the DM for each subproblem.
1596: Generalizes DMCreateFieldIS().
1598: Not collective
1600: Input Parameter:
1601: . dm - the DM object
1603: Output Parameters:
1604: + len - The number of subproblems in the field decomposition (or NULL if not requested)
1605: . namelist - The name for each field (or NULL if not requested)
1606: . islist - The global indices for each field (or NULL if not requested)
1607: - dmlist - The DMs for each field subproblem (or NULL, if not requested; if NULL is returned, no DMs are defined)
1609: Level: intermediate
1611: Notes:
1612: The user is responsible for freeing all requested arrays. In particular, every entry of names should be freed with
1613: PetscFree(), every entry of is should be destroyed with ISDestroy(), every entry of dm should be destroyed with DMDestroy(),
1614: and all of the arrays should be freed with PetscFree().
1616: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateFieldIS()
1617: @*/
1618: PetscErrorCode DMCreateFieldDecomposition(DM dm, PetscInt *len, char ***namelist, IS **islist, DM **dmlist)
1619: {
1624: if (len) {
1626: *len = 0;
1627: }
1628: if (namelist) {
1630: *namelist = 0;
1631: }
1632: if (islist) {
1634: *islist = 0;
1635: }
1636: if (dmlist) {
1638: *dmlist = 0;
1639: }
1640: /*
1641: Is it a good idea to apply the following check across all impls?
1642: Perhaps some impls can have a well-defined decomposition before DMSetUp?
1643: This, however, follows the general principle that accessors are not well-behaved until the object is set up.
1644: */
1645: if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE, "Decomposition defined only after DMSetUp");
1646: if (!dm->ops->createfielddecomposition) {
1647: PetscSection section;
1648: PetscInt numFields, f;
1650: DMGetSection(dm, §ion);
1651: if (section) {PetscSectionGetNumFields(section, &numFields);}
1652: if (section && numFields && dm->ops->createsubdm) {
1653: if (len) *len = numFields;
1654: if (namelist) {PetscMalloc1(numFields,namelist);}
1655: if (islist) {PetscMalloc1(numFields,islist);}
1656: if (dmlist) {PetscMalloc1(numFields,dmlist);}
1657: for (f = 0; f < numFields; ++f) {
1658: const char *fieldName;
1660: DMCreateSubDM(dm, 1, &f, islist ? &(*islist)[f] : NULL, dmlist ? &(*dmlist)[f] : NULL);
1661: if (namelist) {
1662: PetscSectionGetFieldName(section, f, &fieldName);
1663: PetscStrallocpy(fieldName, (char**) &(*namelist)[f]);
1664: }
1665: }
1666: } else {
1667: DMCreateFieldIS(dm, len, namelist, islist);
1668: /* By default there are no DMs associated with subproblems. */
1669: if (dmlist) *dmlist = NULL;
1670: }
1671: } else {
1672: (*dm->ops->createfielddecomposition)(dm,len,namelist,islist,dmlist);
1673: }
1674: return(0);
1675: }
1677: /*@
1678: DMCreateSubDM - Returns an IS and DM encapsulating a subproblem defined by the fields passed in.
1679: The fields are defined by DMCreateFieldIS().
1681: Not collective
1683: Input Parameters:
1684: + dm - The DM object
1685: . numFields - The number of fields in this subproblem
1686: - fields - The field numbers of the selected fields
1688: Output Parameters:
1689: + is - The global indices for the subproblem
1690: - subdm - The DM for the subproblem
1692: Note: You need to call DMPlexSetMigrationSF() on the original DM if you want the Global-To-Natural map to be automatically constructed
1694: Level: intermediate
1696: .seealso DMPlexSetMigrationSF(), DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateFieldIS()
1697: @*/
1698: PetscErrorCode DMCreateSubDM(DM dm, PetscInt numFields, const PetscInt fields[], IS *is, DM *subdm)
1699: {
1707: if (dm->ops->createsubdm) {
1708: (*dm->ops->createsubdm)(dm, numFields, fields, is, subdm);
1709: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "This type has no DMCreateSubDM implementation defined");
1710: return(0);
1711: }
1713: /*@C
1714: DMCreateSuperDM - Returns an arrays of ISes and DM encapsulating a superproblem defined by the DMs passed in.
1716: Not collective
1718: Input Parameter:
1719: + dms - The DM objects
1720: - len - The number of DMs
1722: Output Parameters:
1723: + is - The global indices for the subproblem, or NULL
1724: - superdm - The DM for the superproblem
1726: Note: You need to call DMPlexSetMigrationSF() on the original DM if you want the Global-To-Natural map to be automatically constructed
1728: Level: intermediate
1730: .seealso DMPlexSetMigrationSF(), DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateFieldIS()
1731: @*/
1732: PetscErrorCode DMCreateSuperDM(DM dms[], PetscInt len, IS **is, DM *superdm)
1733: {
1734: PetscInt i;
1742: if (len < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of DMs must be nonnegative: %D", len);
1743: if (len) {
1744: if (dms[0]->ops->createsuperdm) {(*dms[0]->ops->createsuperdm)(dms, len, is, superdm);}
1745: else SETERRQ(PetscObjectComm((PetscObject) dms[0]), PETSC_ERR_SUP, "This type has no DMCreateSuperDM implementation defined");
1746: }
1747: return(0);
1748: }
1751: /*@C
1752: DMCreateDomainDecomposition - Returns lists of IS objects defining a decomposition of a problem into subproblems
1753: corresponding to restrictions to pairs nested subdomains: each IS contains the global
1754: indices of the dofs of the corresponding subdomains. The inner subdomains conceptually
1755: define a nonoverlapping covering, while outer subdomains can overlap.
1756: The optional list of DMs define the DM for each subproblem.
1758: Not collective
1760: Input Parameter:
1761: . dm - the DM object
1763: Output Parameters:
1764: + len - The number of subproblems in the domain decomposition (or NULL if not requested)
1765: . namelist - The name for each subdomain (or NULL if not requested)
1766: . innerislist - The global indices for each inner subdomain (or NULL, if not requested)
1767: . outerislist - The global indices for each outer subdomain (or NULL, if not requested)
1768: - dmlist - The DMs for each subdomain subproblem (or NULL, if not requested; if NULL is returned, no DMs are defined)
1770: Level: intermediate
1772: Notes:
1773: The user is responsible for freeing all requested arrays. In particular, every entry of names should be freed with
1774: PetscFree(), every entry of is should be destroyed with ISDestroy(), every entry of dm should be destroyed with DMDestroy(),
1775: and all of the arrays should be freed with PetscFree().
1777: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateDomainDecompositionDM(), DMCreateFieldDecomposition()
1778: @*/
1779: PetscErrorCode DMCreateDomainDecomposition(DM dm, PetscInt *len, char ***namelist, IS **innerislist, IS **outerislist, DM **dmlist)
1780: {
1781: PetscErrorCode ierr;
1782: DMSubDomainHookLink link;
1783: PetscInt i,l;
1792: /*
1793: Is it a good idea to apply the following check across all impls?
1794: Perhaps some impls can have a well-defined decomposition before DMSetUp?
1795: This, however, follows the general principle that accessors are not well-behaved until the object is set up.
1796: */
1797: if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE, "Decomposition defined only after DMSetUp");
1798: if (dm->ops->createdomaindecomposition) {
1799: (*dm->ops->createdomaindecomposition)(dm,&l,namelist,innerislist,outerislist,dmlist);
1800: /* copy subdomain hooks and context over to the subdomain DMs */
1801: if (dmlist && *dmlist) {
1802: for (i = 0; i < l; i++) {
1803: for (link=dm->subdomainhook; link; link=link->next) {
1804: if (link->ddhook) {(*link->ddhook)(dm,(*dmlist)[i],link->ctx);}
1805: }
1806: if (dm->ctx) (*dmlist)[i]->ctx = dm->ctx;
1807: }
1808: }
1809: if (len) *len = l;
1810: }
1811: return(0);
1812: }
1815: /*@C
1816: DMCreateDomainDecompositionScatters - Returns scatters to the subdomain vectors from the global vector
1818: Not collective
1820: Input Parameters:
1821: + dm - the DM object
1822: . n - the number of subdomain scatters
1823: - subdms - the local subdomains
1825: Output Parameters:
1826: + n - the number of scatters returned
1827: . iscat - scatter from global vector to nonoverlapping global vector entries on subdomain
1828: . oscat - scatter from global vector to overlapping global vector entries on subdomain
1829: - gscat - scatter from global vector to local vector on subdomain (fills in ghosts)
1831: Notes:
1832: This is an alternative to the iis and ois arguments in DMCreateDomainDecomposition that allow for the solution
1833: of general nonlinear problems with overlapping subdomain methods. While merely having index sets that enable subsets
1834: of the residual equations to be created is fine for linear problems, nonlinear problems require local assembly of
1835: solution and residual data.
1837: Level: developer
1839: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateFieldIS()
1840: @*/
1841: PetscErrorCode DMCreateDomainDecompositionScatters(DM dm,PetscInt n,DM *subdms,VecScatter **iscat,VecScatter **oscat,VecScatter **gscat)
1842: {
1848: if (dm->ops->createddscatters) {
1849: (*dm->ops->createddscatters)(dm,n,subdms,iscat,oscat,gscat);
1850: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "This type has no DMCreateDomainDecompositionScatter implementation defined");
1851: return(0);
1852: }
1854: /*@
1855: DMRefine - Refines a DM object
1857: Collective on DM
1859: Input Parameter:
1860: + dm - the DM object
1861: - comm - the communicator to contain the new DM object (or MPI_COMM_NULL)
1863: Output Parameter:
1864: . dmf - the refined DM, or NULL
1866: Note: If no refinement was done, the return value is NULL
1868: Level: developer
1870: .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation()
1871: @*/
1872: PetscErrorCode DMRefine(DM dm,MPI_Comm comm,DM *dmf)
1873: {
1874: PetscErrorCode ierr;
1875: DMRefineHookLink link;
1879: PetscLogEventBegin(DM_Refine,dm,0,0,0);
1880: if (!dm->ops->refine) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"This DM cannot refine");
1881: (*dm->ops->refine)(dm,comm,dmf);
1882: if (*dmf) {
1883: (*dmf)->ops->creatematrix = dm->ops->creatematrix;
1885: PetscObjectCopyFortranFunctionPointers((PetscObject)dm,(PetscObject)*dmf);
1887: (*dmf)->ctx = dm->ctx;
1888: (*dmf)->leveldown = dm->leveldown;
1889: (*dmf)->levelup = dm->levelup + 1;
1891: DMSetMatType(*dmf,dm->mattype);
1892: for (link=dm->refinehook; link; link=link->next) {
1893: if (link->refinehook) {
1894: (*link->refinehook)(dm,*dmf,link->ctx);
1895: }
1896: }
1897: }
1898: PetscLogEventEnd(DM_Refine,dm,0,0,0);
1899: return(0);
1900: }
1902: /*@C
1903: DMRefineHookAdd - adds a callback to be run when interpolating a nonlinear problem to a finer grid
1905: Logically Collective
1907: Input Arguments:
1908: + coarse - nonlinear solver context on which to run a hook when restricting to a coarser level
1909: . refinehook - function to run when setting up a coarser level
1910: . interphook - function to run to update data on finer levels (once per SNESSolve())
1911: - ctx - [optional] user-defined context for provide data for the hooks (may be NULL)
1913: Calling sequence of refinehook:
1914: $ refinehook(DM coarse,DM fine,void *ctx);
1916: + coarse - coarse level DM
1917: . fine - fine level DM to interpolate problem to
1918: - ctx - optional user-defined function context
1920: Calling sequence for interphook:
1921: $ interphook(DM coarse,Mat interp,DM fine,void *ctx)
1923: + coarse - coarse level DM
1924: . interp - matrix interpolating a coarse-level solution to the finer grid
1925: . fine - fine level DM to update
1926: - ctx - optional user-defined function context
1928: Level: advanced
1930: Notes:
1931: This function is only needed if auxiliary data needs to be passed to fine grids while grid sequencing
1933: If this function is called multiple times, the hooks will be run in the order they are added.
1935: This function is currently not available from Fortran.
1937: .seealso: DMCoarsenHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
1938: @*/
1939: PetscErrorCode DMRefineHookAdd(DM coarse,PetscErrorCode (*refinehook)(DM,DM,void*),PetscErrorCode (*interphook)(DM,Mat,DM,void*),void *ctx)
1940: {
1941: PetscErrorCode ierr;
1942: DMRefineHookLink link,*p;
1946: for (p=&coarse->refinehook; *p; p=&(*p)->next) { /* Scan to the end of the current list of hooks */
1947: if ((*p)->refinehook == refinehook && (*p)->interphook == interphook && (*p)->ctx == ctx) return(0);
1948: }
1949: PetscNew(&link);
1950: link->refinehook = refinehook;
1951: link->interphook = interphook;
1952: link->ctx = ctx;
1953: link->next = NULL;
1954: *p = link;
1955: return(0);
1956: }
1958: /*@C
1959: DMRefineHookRemove - remove a callback from the list of hooks to be run when interpolating a nonlinear problem to a finer grid
1961: Logically Collective
1963: Input Arguments:
1964: + coarse - nonlinear solver context on which to run a hook when restricting to a coarser level
1965: . refinehook - function to run when setting up a coarser level
1966: . interphook - function to run to update data on finer levels (once per SNESSolve())
1967: - ctx - [optional] user-defined context for provide data for the hooks (may be NULL)
1969: Level: advanced
1971: Notes:
1972: This function does nothing if the hook is not in the list.
1974: This function is currently not available from Fortran.
1976: .seealso: DMCoarsenHookRemove(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
1977: @*/
1978: PetscErrorCode DMRefineHookRemove(DM coarse,PetscErrorCode (*refinehook)(DM,DM,void*),PetscErrorCode (*interphook)(DM,Mat,DM,void*),void *ctx)
1979: {
1980: PetscErrorCode ierr;
1981: DMRefineHookLink link,*p;
1985: for (p=&coarse->refinehook; *p; p=&(*p)->next) { /* Search the list of current hooks */
1986: if ((*p)->refinehook == refinehook && (*p)->interphook == interphook && (*p)->ctx == ctx) {
1987: link = *p;
1988: *p = link->next;
1989: PetscFree(link);
1990: break;
1991: }
1992: }
1993: return(0);
1994: }
1996: /*@
1997: DMInterpolate - interpolates user-defined problem data to a finer DM by running hooks registered by DMRefineHookAdd()
1999: Collective if any hooks are
2001: Input Arguments:
2002: + coarse - coarser DM to use as a base
2003: . interp - interpolation matrix, apply using MatInterpolate()
2004: - fine - finer DM to update
2006: Level: developer
2008: .seealso: DMRefineHookAdd(), MatInterpolate()
2009: @*/
2010: PetscErrorCode DMInterpolate(DM coarse,Mat interp,DM fine)
2011: {
2012: PetscErrorCode ierr;
2013: DMRefineHookLink link;
2016: for (link=fine->refinehook; link; link=link->next) {
2017: if (link->interphook) {
2018: (*link->interphook)(coarse,interp,fine,link->ctx);
2019: }
2020: }
2021: return(0);
2022: }
2024: /*@
2025: DMGetRefineLevel - Get's the number of refinements that have generated this DM.
2027: Not Collective
2029: Input Parameter:
2030: . dm - the DM object
2032: Output Parameter:
2033: . level - number of refinements
2035: Level: developer
2037: .seealso DMCoarsen(), DMGetCoarsenLevel(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation()
2039: @*/
2040: PetscErrorCode DMGetRefineLevel(DM dm,PetscInt *level)
2041: {
2044: *level = dm->levelup;
2045: return(0);
2046: }
2048: /*@
2049: DMSetRefineLevel - Set's the number of refinements that have generated this DM.
2051: Not Collective
2053: Input Parameter:
2054: + dm - the DM object
2055: - level - number of refinements
2057: Level: advanced
2059: Notes:
2060: This value is used by PCMG to determine how many multigrid levels to use
2062: .seealso DMCoarsen(), DMGetCoarsenLevel(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation()
2064: @*/
2065: PetscErrorCode DMSetRefineLevel(DM dm,PetscInt level)
2066: {
2069: dm->levelup = level;
2070: return(0);
2071: }
2073: /*@C
2074: DMGlobalToLocalHookAdd - adds a callback to be run when global to local is called
2076: Logically Collective
2078: Input Arguments:
2079: + dm - the DM
2080: . beginhook - function to run at the beginning of DMGlobalToLocalBegin()
2081: . endhook - function to run after DMGlobalToLocalEnd() has completed
2082: - ctx - [optional] user-defined context for provide data for the hooks (may be NULL)
2084: Calling sequence for beginhook:
2085: $ beginhook(DM fine,VecScatter out,VecScatter in,DM coarse,void *ctx)
2087: + dm - global DM
2088: . g - global vector
2089: . mode - mode
2090: . l - local vector
2091: - ctx - optional user-defined function context
2094: Calling sequence for endhook:
2095: $ endhook(DM fine,VecScatter out,VecScatter in,DM coarse,void *ctx)
2097: + global - global DM
2098: - ctx - optional user-defined function context
2100: Level: advanced
2102: .seealso: DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2103: @*/
2104: PetscErrorCode DMGlobalToLocalHookAdd(DM dm,PetscErrorCode (*beginhook)(DM,Vec,InsertMode,Vec,void*),PetscErrorCode (*endhook)(DM,Vec,InsertMode,Vec,void*),void *ctx)
2105: {
2106: PetscErrorCode ierr;
2107: DMGlobalToLocalHookLink link,*p;
2111: for (p=&dm->gtolhook; *p; p=&(*p)->next) {} /* Scan to the end of the current list of hooks */
2112: PetscNew(&link);
2113: link->beginhook = beginhook;
2114: link->endhook = endhook;
2115: link->ctx = ctx;
2116: link->next = NULL;
2117: *p = link;
2118: return(0);
2119: }
2121: static PetscErrorCode DMGlobalToLocalHook_Constraints(DM dm, Vec g, InsertMode mode, Vec l, void *ctx)
2122: {
2123: Mat cMat;
2124: Vec cVec;
2125: PetscSection section, cSec;
2126: PetscInt pStart, pEnd, p, dof;
2131: DMGetDefaultConstraints(dm,&cSec,&cMat);
2132: if (cMat && (mode == INSERT_VALUES || mode == INSERT_ALL_VALUES || mode == INSERT_BC_VALUES)) {
2133: PetscInt nRows;
2135: MatGetSize(cMat,&nRows,NULL);
2136: if (nRows <= 0) return(0);
2137: DMGetSection(dm,§ion);
2138: MatCreateVecs(cMat,NULL,&cVec);
2139: MatMult(cMat,l,cVec);
2140: PetscSectionGetChart(cSec,&pStart,&pEnd);
2141: for (p = pStart; p < pEnd; p++) {
2142: PetscSectionGetDof(cSec,p,&dof);
2143: if (dof) {
2144: PetscScalar *vals;
2145: VecGetValuesSection(cVec,cSec,p,&vals);
2146: VecSetValuesSection(l,section,p,vals,INSERT_ALL_VALUES);
2147: }
2148: }
2149: VecDestroy(&cVec);
2150: }
2151: return(0);
2152: }
2154: /*@
2155: DMGlobalToLocal - update local vectors from global vector
2157: Neighbor-wise Collective on DM
2159: Input Parameters:
2160: + dm - the DM object
2161: . g - the global vector
2162: . mode - INSERT_VALUES or ADD_VALUES
2163: - l - the local vector
2165: Notes:
2166: The communication involved in this update can be overlapped with computation by using
2167: DMGlobalToLocalBegin() and DMGlobalToLocalEnd().
2169: Level: beginner
2171: .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMGlobalToLocalEnd(), DMLocalToGlobalBegin(), DMLocalToGlobal(), DMLocalToGlobalBegin(), DMLocalToGlobalEnd()
2173: @*/
2174: PetscErrorCode DMGlobalToLocal(DM dm,Vec g,InsertMode mode,Vec l)
2175: {
2179: DMGlobalToLocalBegin(dm,g,mode,l);
2180: DMGlobalToLocalEnd(dm,g,mode,l);
2181: return(0);
2182: }
2184: /*@
2185: DMGlobalToLocalBegin - Begins updating local vectors from global vector
2187: Neighbor-wise Collective on DM
2189: Input Parameters:
2190: + dm - the DM object
2191: . g - the global vector
2192: . mode - INSERT_VALUES or ADD_VALUES
2193: - l - the local vector
2195: Level: intermediate
2197: .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMGlobalToLocal(), DMGlobalToLocalEnd(), DMLocalToGlobalBegin(), DMLocalToGlobal(), DMLocalToGlobalBegin(), DMLocalToGlobalEnd()
2199: @*/
2200: PetscErrorCode DMGlobalToLocalBegin(DM dm,Vec g,InsertMode mode,Vec l)
2201: {
2202: PetscSF sf;
2203: PetscErrorCode ierr;
2204: DMGlobalToLocalHookLink link;
2208: for (link=dm->gtolhook; link; link=link->next) {
2209: if (link->beginhook) {
2210: (*link->beginhook)(dm,g,mode,l,link->ctx);
2211: }
2212: }
2213: DMGetDefaultSF(dm, &sf);
2214: if (sf) {
2215: const PetscScalar *gArray;
2216: PetscScalar *lArray;
2218: if (mode == ADD_VALUES) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode);
2219: VecGetArray(l, &lArray);
2220: VecGetArrayRead(g, &gArray);
2221: PetscSFBcastBegin(sf, MPIU_SCALAR, gArray, lArray);
2222: VecRestoreArray(l, &lArray);
2223: VecRestoreArrayRead(g, &gArray);
2224: } else {
2225: if (!dm->ops->globaltolocalbegin) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Missing DMGlobalToLocalBegin() for type %s",((PetscObject)dm)->type_name);
2226: (*dm->ops->globaltolocalbegin)(dm,g,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),l);
2227: }
2228: return(0);
2229: }
2231: /*@
2232: DMGlobalToLocalEnd - Ends updating local vectors from global vector
2234: Neighbor-wise Collective on DM
2236: Input Parameters:
2237: + dm - the DM object
2238: . g - the global vector
2239: . mode - INSERT_VALUES or ADD_VALUES
2240: - l - the local vector
2242: Level: intermediate
2244: .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMGlobalToLocal(), DMLocalToGlobalBegin(), DMLocalToGlobal(), DMLocalToGlobalBegin(), DMLocalToGlobalEnd()
2246: @*/
2247: PetscErrorCode DMGlobalToLocalEnd(DM dm,Vec g,InsertMode mode,Vec l)
2248: {
2249: PetscSF sf;
2250: PetscErrorCode ierr;
2251: const PetscScalar *gArray;
2252: PetscScalar *lArray;
2253: DMGlobalToLocalHookLink link;
2257: DMGetDefaultSF(dm, &sf);
2258: if (sf) {
2259: if (mode == ADD_VALUES) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode);
2261: VecGetArray(l, &lArray);
2262: VecGetArrayRead(g, &gArray);
2263: PetscSFBcastEnd(sf, MPIU_SCALAR, gArray, lArray);
2264: VecRestoreArray(l, &lArray);
2265: VecRestoreArrayRead(g, &gArray);
2266: } else {
2267: if (!dm->ops->globaltolocalend) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Missing DMGlobalToLocalEnd() for type %s",((PetscObject)dm)->type_name);
2268: (*dm->ops->globaltolocalend)(dm,g,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),l);
2269: }
2270: DMGlobalToLocalHook_Constraints(dm,g,mode,l,NULL);
2271: for (link=dm->gtolhook; link; link=link->next) {
2272: if (link->endhook) {(*link->endhook)(dm,g,mode,l,link->ctx);}
2273: }
2274: return(0);
2275: }
2277: /*@C
2278: DMLocalToGlobalHookAdd - adds a callback to be run when a local to global is called
2280: Logically Collective
2282: Input Arguments:
2283: + dm - the DM
2284: . beginhook - function to run at the beginning of DMLocalToGlobalBegin()
2285: . endhook - function to run after DMLocalToGlobalEnd() has completed
2286: - ctx - [optional] user-defined context for provide data for the hooks (may be NULL)
2288: Calling sequence for beginhook:
2289: $ beginhook(DM fine,Vec l,InsertMode mode,Vec g,void *ctx)
2291: + dm - global DM
2292: . l - local vector
2293: . mode - mode
2294: . g - global vector
2295: - ctx - optional user-defined function context
2298: Calling sequence for endhook:
2299: $ endhook(DM fine,Vec l,InsertMode mode,Vec g,void *ctx)
2301: + global - global DM
2302: . l - local vector
2303: . mode - mode
2304: . g - global vector
2305: - ctx - optional user-defined function context
2307: Level: advanced
2309: .seealso: DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2310: @*/
2311: PetscErrorCode DMLocalToGlobalHookAdd(DM dm,PetscErrorCode (*beginhook)(DM,Vec,InsertMode,Vec,void*),PetscErrorCode (*endhook)(DM,Vec,InsertMode,Vec,void*),void *ctx)
2312: {
2313: PetscErrorCode ierr;
2314: DMLocalToGlobalHookLink link,*p;
2318: for (p=&dm->ltoghook; *p; p=&(*p)->next) {} /* Scan to the end of the current list of hooks */
2319: PetscNew(&link);
2320: link->beginhook = beginhook;
2321: link->endhook = endhook;
2322: link->ctx = ctx;
2323: link->next = NULL;
2324: *p = link;
2325: return(0);
2326: }
2328: static PetscErrorCode DMLocalToGlobalHook_Constraints(DM dm, Vec l, InsertMode mode, Vec g, void *ctx)
2329: {
2330: Mat cMat;
2331: Vec cVec;
2332: PetscSection section, cSec;
2333: PetscInt pStart, pEnd, p, dof;
2338: DMGetDefaultConstraints(dm,&cSec,&cMat);
2339: if (cMat && (mode == ADD_VALUES || mode == ADD_ALL_VALUES || mode == ADD_BC_VALUES)) {
2340: PetscInt nRows;
2342: MatGetSize(cMat,&nRows,NULL);
2343: if (nRows <= 0) return(0);
2344: DMGetSection(dm,§ion);
2345: MatCreateVecs(cMat,NULL,&cVec);
2346: PetscSectionGetChart(cSec,&pStart,&pEnd);
2347: for (p = pStart; p < pEnd; p++) {
2348: PetscSectionGetDof(cSec,p,&dof);
2349: if (dof) {
2350: PetscInt d;
2351: PetscScalar *vals;
2352: VecGetValuesSection(l,section,p,&vals);
2353: VecSetValuesSection(cVec,cSec,p,vals,mode);
2354: /* for this to be the true transpose, we have to zero the values that
2355: * we just extracted */
2356: for (d = 0; d < dof; d++) {
2357: vals[d] = 0.;
2358: }
2359: }
2360: }
2361: MatMultTransposeAdd(cMat,cVec,l,l);
2362: VecDestroy(&cVec);
2363: }
2364: return(0);
2365: }
2366: /*@
2367: DMLocalToGlobal - updates global vectors from local vectors
2369: Neighbor-wise Collective on DM
2371: Input Parameters:
2372: + dm - the DM object
2373: . l - the local vector
2374: . 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.
2375: - g - the global vector
2377: Notes:
2378: The communication involved in this update can be overlapped with computation by using
2379: DMLocalToGlobalBegin() and DMLocalToGlobalEnd().
2381: In the ADD_VALUES case you normally would zero the receiving vector before beginning this operation.
2382: INSERT_VALUES is not supported for DMDA; in that case simply compute the values directly into a global vector instead of a local one.
2384: Level: beginner
2386: .seealso DMLocalToGlobalBegin(), DMLocalToGlobalEnd(), DMCoarsen(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMGlobalToLocal(), DMGlobalToLocalEnd(), DMGlobalToLocalBegin()
2388: @*/
2389: PetscErrorCode DMLocalToGlobal(DM dm,Vec l,InsertMode mode,Vec g)
2390: {
2394: DMLocalToGlobalBegin(dm,l,mode,g);
2395: DMLocalToGlobalEnd(dm,l,mode,g);
2396: return(0);
2397: }
2399: /*@
2400: DMLocalToGlobalBegin - begins updating global vectors from local vectors
2402: Neighbor-wise Collective on DM
2404: Input Parameters:
2405: + dm - the DM object
2406: . l - the local vector
2407: . 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.
2408: - g - the global vector
2410: Notes:
2411: In the ADD_VALUES case you normally would zero the receiving vector before beginning this operation.
2412: INSERT_VALUES is not supported for DMDA, in that case simply compute the values directly into a global vector instead of a local one.
2414: Level: intermediate
2416: .seealso DMLocalToGlobal(), DMLocalToGlobalEnd(), DMCoarsen(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMGlobalToLocal(), DMGlobalToLocalEnd(), DMGlobalToLocalBegin()
2418: @*/
2419: PetscErrorCode DMLocalToGlobalBegin(DM dm,Vec l,InsertMode mode,Vec g)
2420: {
2421: PetscSF sf;
2422: PetscSection s, gs;
2423: DMLocalToGlobalHookLink link;
2424: const PetscScalar *lArray;
2425: PetscScalar *gArray;
2426: PetscBool isInsert;
2427: PetscErrorCode ierr;
2431: for (link=dm->ltoghook; link; link=link->next) {
2432: if (link->beginhook) {
2433: (*link->beginhook)(dm,l,mode,g,link->ctx);
2434: }
2435: }
2436: DMLocalToGlobalHook_Constraints(dm,l,mode,g,NULL);
2437: DMGetDefaultSF(dm, &sf);
2438: DMGetSection(dm, &s);
2439: switch (mode) {
2440: case INSERT_VALUES:
2441: case INSERT_ALL_VALUES:
2442: case INSERT_BC_VALUES:
2443: isInsert = PETSC_TRUE; break;
2444: case ADD_VALUES:
2445: case ADD_ALL_VALUES:
2446: case ADD_BC_VALUES:
2447: isInsert = PETSC_FALSE; break;
2448: default:
2449: SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode);
2450: }
2451: if (sf && !isInsert) {
2452: VecGetArrayRead(l, &lArray);
2453: VecGetArray(g, &gArray);
2454: PetscSFReduceBegin(sf, MPIU_SCALAR, lArray, gArray, MPIU_SUM);
2455: VecRestoreArrayRead(l, &lArray);
2456: VecRestoreArray(g, &gArray);
2457: } else if (s && isInsert) {
2458: PetscInt gStart, pStart, pEnd, p;
2460: DMGetGlobalSection(dm, &gs);
2461: PetscSectionGetChart(s, &pStart, &pEnd);
2462: VecGetOwnershipRange(g, &gStart, NULL);
2463: VecGetArrayRead(l, &lArray);
2464: VecGetArray(g, &gArray);
2465: for (p = pStart; p < pEnd; ++p) {
2466: PetscInt dof, gdof, cdof, gcdof, off, goff, d, e;
2468: PetscSectionGetDof(s, p, &dof);
2469: PetscSectionGetDof(gs, p, &gdof);
2470: PetscSectionGetConstraintDof(s, p, &cdof);
2471: PetscSectionGetConstraintDof(gs, p, &gcdof);
2472: PetscSectionGetOffset(s, p, &off);
2473: PetscSectionGetOffset(gs, p, &goff);
2474: /* Ignore off-process data and points with no global data */
2475: if (!gdof || goff < 0) continue;
2476: 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);
2477: /* If no constraints are enforced in the global vector */
2478: if (!gcdof) {
2479: for (d = 0; d < dof; ++d) gArray[goff-gStart+d] = lArray[off+d];
2480: /* If constraints are enforced in the global vector */
2481: } else if (cdof == gcdof) {
2482: const PetscInt *cdofs;
2483: PetscInt cind = 0;
2485: PetscSectionGetConstraintIndices(s, p, &cdofs);
2486: for (d = 0, e = 0; d < dof; ++d) {
2487: if ((cind < cdof) && (d == cdofs[cind])) {++cind; continue;}
2488: gArray[goff-gStart+e++] = lArray[off+d];
2489: }
2490: } 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);
2491: }
2492: VecRestoreArrayRead(l, &lArray);
2493: VecRestoreArray(g, &gArray);
2494: } else {
2495: (*dm->ops->localtoglobalbegin)(dm,l,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),g);
2496: }
2497: return(0);
2498: }
2500: /*@
2501: DMLocalToGlobalEnd - updates global vectors from local vectors
2503: Neighbor-wise Collective on DM
2505: Input Parameters:
2506: + dm - the DM object
2507: . l - the local vector
2508: . mode - INSERT_VALUES or ADD_VALUES
2509: - g - the global vector
2511: Level: intermediate
2513: .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMGlobalToLocalEnd(), DMGlobalToLocalEnd()
2515: @*/
2516: PetscErrorCode DMLocalToGlobalEnd(DM dm,Vec l,InsertMode mode,Vec g)
2517: {
2518: PetscSF sf;
2519: PetscSection s;
2520: DMLocalToGlobalHookLink link;
2521: PetscBool isInsert;
2522: PetscErrorCode ierr;
2526: DMGetDefaultSF(dm, &sf);
2527: DMGetSection(dm, &s);
2528: switch (mode) {
2529: case INSERT_VALUES:
2530: case INSERT_ALL_VALUES:
2531: isInsert = PETSC_TRUE; break;
2532: case ADD_VALUES:
2533: case ADD_ALL_VALUES:
2534: isInsert = PETSC_FALSE; break;
2535: default:
2536: SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode);
2537: }
2538: if (sf && !isInsert) {
2539: const PetscScalar *lArray;
2540: PetscScalar *gArray;
2542: VecGetArrayRead(l, &lArray);
2543: VecGetArray(g, &gArray);
2544: PetscSFReduceEnd(sf, MPIU_SCALAR, lArray, gArray, MPIU_SUM);
2545: VecRestoreArrayRead(l, &lArray);
2546: VecRestoreArray(g, &gArray);
2547: } else if (s && isInsert) {
2548: } else {
2549: (*dm->ops->localtoglobalend)(dm,l,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),g);
2550: }
2551: for (link=dm->ltoghook; link; link=link->next) {
2552: if (link->endhook) {(*link->endhook)(dm,g,mode,l,link->ctx);}
2553: }
2554: return(0);
2555: }
2557: /*@
2558: DMLocalToLocalBegin - Maps from a local vector (including ghost points
2559: that contain irrelevant values) to another local vector where the ghost
2560: points in the second are set correctly. Must be followed by DMLocalToLocalEnd().
2562: Neighbor-wise Collective on DM and Vec
2564: Input Parameters:
2565: + dm - the DM object
2566: . g - the original local vector
2567: - mode - one of INSERT_VALUES or ADD_VALUES
2569: Output Parameter:
2570: . l - the local vector with correct ghost values
2572: Level: intermediate
2574: Notes:
2575: The local vectors used here need not be the same as those
2576: obtained from DMCreateLocalVector(), BUT they
2577: must have the same parallel data layout; they could, for example, be
2578: obtained with VecDuplicate() from the DM originating vectors.
2580: .keywords: DM, local-to-local, begin
2581: .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateLocalVector(), DMCreateGlobalVector(), DMCreateInterpolation(), DMLocalToLocalEnd(), DMGlobalToLocalEnd(), DMLocalToGlobalBegin()
2583: @*/
2584: PetscErrorCode DMLocalToLocalBegin(DM dm,Vec g,InsertMode mode,Vec l)
2585: {
2586: PetscErrorCode ierr;
2590: if (!dm->ops->localtolocalbegin) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"This DM does not support local to local maps");
2591: (*dm->ops->localtolocalbegin)(dm,g,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),l);
2592: return(0);
2593: }
2595: /*@
2596: DMLocalToLocalEnd - Maps from a local vector (including ghost points
2597: that contain irrelevant values) to another local vector where the ghost
2598: points in the second are set correctly. Must be preceded by DMLocalToLocalBegin().
2600: Neighbor-wise Collective on DM and Vec
2602: Input Parameters:
2603: + da - the DM object
2604: . g - the original local vector
2605: - mode - one of INSERT_VALUES or ADD_VALUES
2607: Output Parameter:
2608: . l - the local vector with correct ghost values
2610: Level: intermediate
2612: Notes:
2613: The local vectors used here need not be the same as those
2614: obtained from DMCreateLocalVector(), BUT they
2615: must have the same parallel data layout; they could, for example, be
2616: obtained with VecDuplicate() from the DM originating vectors.
2618: .keywords: DM, local-to-local, end
2619: .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateLocalVector(), DMCreateGlobalVector(), DMCreateInterpolation(), DMLocalToLocalBegin(), DMGlobalToLocalEnd(), DMLocalToGlobalBegin()
2621: @*/
2622: PetscErrorCode DMLocalToLocalEnd(DM dm,Vec g,InsertMode mode,Vec l)
2623: {
2624: PetscErrorCode ierr;
2628: if (!dm->ops->localtolocalend) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"This DM does not support local to local maps");
2629: (*dm->ops->localtolocalend)(dm,g,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),l);
2630: return(0);
2631: }
2634: /*@
2635: DMCoarsen - Coarsens a DM object
2637: Collective on DM
2639: Input Parameter:
2640: + dm - the DM object
2641: - comm - the communicator to contain the new DM object (or MPI_COMM_NULL)
2643: Output Parameter:
2644: . dmc - the coarsened DM
2646: Level: developer
2648: .seealso DMRefine(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation()
2650: @*/
2651: PetscErrorCode DMCoarsen(DM dm, MPI_Comm comm, DM *dmc)
2652: {
2653: PetscErrorCode ierr;
2654: DMCoarsenHookLink link;
2658: if (!dm->ops->coarsen) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"This DM cannot coarsen");
2659: PetscLogEventBegin(DM_Coarsen,dm,0,0,0);
2660: (*dm->ops->coarsen)(dm, comm, dmc);
2661: if (!(*dmc)) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "NULL coarse mesh produced");
2662: DMSetCoarseDM(dm,*dmc);
2663: (*dmc)->ops->creatematrix = dm->ops->creatematrix;
2664: PetscObjectCopyFortranFunctionPointers((PetscObject)dm,(PetscObject)*dmc);
2665: (*dmc)->ctx = dm->ctx;
2666: (*dmc)->levelup = dm->levelup;
2667: (*dmc)->leveldown = dm->leveldown + 1;
2668: DMSetMatType(*dmc,dm->mattype);
2669: for (link=dm->coarsenhook; link; link=link->next) {
2670: if (link->coarsenhook) {(*link->coarsenhook)(dm,*dmc,link->ctx);}
2671: }
2672: PetscLogEventEnd(DM_Coarsen,dm,0,0,0);
2673: return(0);
2674: }
2676: /*@C
2677: DMCoarsenHookAdd - adds a callback to be run when restricting a nonlinear problem to the coarse grid
2679: Logically Collective
2681: Input Arguments:
2682: + fine - nonlinear solver context on which to run a hook when restricting to a coarser level
2683: . coarsenhook - function to run when setting up a coarser level
2684: . restricthook - function to run to update data on coarser levels (once per SNESSolve())
2685: - ctx - [optional] user-defined context for provide data for the hooks (may be NULL)
2687: Calling sequence of coarsenhook:
2688: $ coarsenhook(DM fine,DM coarse,void *ctx);
2690: + fine - fine level DM
2691: . coarse - coarse level DM to restrict problem to
2692: - ctx - optional user-defined function context
2694: Calling sequence for restricthook:
2695: $ restricthook(DM fine,Mat mrestrict,Vec rscale,Mat inject,DM coarse,void *ctx)
2697: + fine - fine level DM
2698: . mrestrict - matrix restricting a fine-level solution to the coarse grid
2699: . rscale - scaling vector for restriction
2700: . inject - matrix restricting by injection
2701: . coarse - coarse level DM to update
2702: - ctx - optional user-defined function context
2704: Level: advanced
2706: Notes:
2707: This function is only needed if auxiliary data needs to be set up on coarse grids.
2709: If this function is called multiple times, the hooks will be run in the order they are added.
2711: In order to compose with nonlinear preconditioning without duplicating storage, the hook should be implemented to
2712: extract the finest level information from its context (instead of from the SNES).
2714: This function is currently not available from Fortran.
2716: .seealso: DMCoarsenHookRemove(), DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2717: @*/
2718: PetscErrorCode DMCoarsenHookAdd(DM fine,PetscErrorCode (*coarsenhook)(DM,DM,void*),PetscErrorCode (*restricthook)(DM,Mat,Vec,Mat,DM,void*),void *ctx)
2719: {
2720: PetscErrorCode ierr;
2721: DMCoarsenHookLink link,*p;
2725: for (p=&fine->coarsenhook; *p; p=&(*p)->next) { /* Scan to the end of the current list of hooks */
2726: if ((*p)->coarsenhook == coarsenhook && (*p)->restricthook == restricthook && (*p)->ctx == ctx) return(0);
2727: }
2728: PetscNew(&link);
2729: link->coarsenhook = coarsenhook;
2730: link->restricthook = restricthook;
2731: link->ctx = ctx;
2732: link->next = NULL;
2733: *p = link;
2734: return(0);
2735: }
2737: /*@C
2738: DMCoarsenHookRemove - remove a callback from the list of hooks to be run when restricting a nonlinear problem to the coarse grid
2740: Logically Collective
2742: Input Arguments:
2743: + fine - nonlinear solver context on which to run a hook when restricting to a coarser level
2744: . coarsenhook - function to run when setting up a coarser level
2745: . restricthook - function to run to update data on coarser levels (once per SNESSolve())
2746: - ctx - [optional] user-defined context for provide data for the hooks (may be NULL)
2748: Level: advanced
2750: Notes:
2751: This function does nothing if the hook is not in the list.
2753: This function is currently not available from Fortran.
2755: .seealso: DMCoarsenHookAdd(), DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2756: @*/
2757: PetscErrorCode DMCoarsenHookRemove(DM fine,PetscErrorCode (*coarsenhook)(DM,DM,void*),PetscErrorCode (*restricthook)(DM,Mat,Vec,Mat,DM,void*),void *ctx)
2758: {
2759: PetscErrorCode ierr;
2760: DMCoarsenHookLink link,*p;
2764: for (p=&fine->coarsenhook; *p; p=&(*p)->next) { /* Search the list of current hooks */
2765: if ((*p)->coarsenhook == coarsenhook && (*p)->restricthook == restricthook && (*p)->ctx == ctx) {
2766: link = *p;
2767: *p = link->next;
2768: PetscFree(link);
2769: break;
2770: }
2771: }
2772: return(0);
2773: }
2776: /*@
2777: DMRestrict - restricts user-defined problem data to a coarser DM by running hooks registered by DMCoarsenHookAdd()
2779: Collective if any hooks are
2781: Input Arguments:
2782: + fine - finer DM to use as a base
2783: . restrct - restriction matrix, apply using MatRestrict()
2784: . rscale - scaling vector for restriction
2785: . inject - injection matrix, also use MatRestrict()
2786: - coarse - coarser DM to update
2788: Level: developer
2790: .seealso: DMCoarsenHookAdd(), MatRestrict()
2791: @*/
2792: PetscErrorCode DMRestrict(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse)
2793: {
2794: PetscErrorCode ierr;
2795: DMCoarsenHookLink link;
2798: for (link=fine->coarsenhook; link; link=link->next) {
2799: if (link->restricthook) {
2800: (*link->restricthook)(fine,restrct,rscale,inject,coarse,link->ctx);
2801: }
2802: }
2803: return(0);
2804: }
2806: /*@C
2807: DMSubDomainHookAdd - adds a callback to be run when restricting a problem to the coarse grid
2809: Logically Collective
2811: Input Arguments:
2812: + global - global DM
2813: . ddhook - function to run to pass data to the decomposition DM upon its creation
2814: . restricthook - function to run to update data on block solve (at the beginning of the block solve)
2815: - ctx - [optional] user-defined context for provide data for the hooks (may be NULL)
2818: Calling sequence for ddhook:
2819: $ ddhook(DM global,DM block,void *ctx)
2821: + global - global DM
2822: . block - block DM
2823: - ctx - optional user-defined function context
2825: Calling sequence for restricthook:
2826: $ restricthook(DM global,VecScatter out,VecScatter in,DM block,void *ctx)
2828: + global - global DM
2829: . out - scatter to the outer (with ghost and overlap points) block vector
2830: . in - scatter to block vector values only owned locally
2831: . block - block DM
2832: - ctx - optional user-defined function context
2834: Level: advanced
2836: Notes:
2837: This function is only needed if auxiliary data needs to be set up on subdomain DMs.
2839: If this function is called multiple times, the hooks will be run in the order they are added.
2841: In order to compose with nonlinear preconditioning without duplicating storage, the hook should be implemented to
2842: extract the global information from its context (instead of from the SNES).
2844: This function is currently not available from Fortran.
2846: .seealso: DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2847: @*/
2848: PetscErrorCode DMSubDomainHookAdd(DM global,PetscErrorCode (*ddhook)(DM,DM,void*),PetscErrorCode (*restricthook)(DM,VecScatter,VecScatter,DM,void*),void *ctx)
2849: {
2850: PetscErrorCode ierr;
2851: DMSubDomainHookLink link,*p;
2855: for (p=&global->subdomainhook; *p; p=&(*p)->next) { /* Scan to the end of the current list of hooks */
2856: if ((*p)->ddhook == ddhook && (*p)->restricthook == restricthook && (*p)->ctx == ctx) return(0);
2857: }
2858: PetscNew(&link);
2859: link->restricthook = restricthook;
2860: link->ddhook = ddhook;
2861: link->ctx = ctx;
2862: link->next = NULL;
2863: *p = link;
2864: return(0);
2865: }
2867: /*@C
2868: DMSubDomainHookRemove - remove a callback from the list to be run when restricting a problem to the coarse grid
2870: Logically Collective
2872: Input Arguments:
2873: + global - global DM
2874: . ddhook - function to run to pass data to the decomposition DM upon its creation
2875: . restricthook - function to run to update data on block solve (at the beginning of the block solve)
2876: - ctx - [optional] user-defined context for provide data for the hooks (may be NULL)
2878: Level: advanced
2880: Notes:
2882: This function is currently not available from Fortran.
2884: .seealso: DMSubDomainHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2885: @*/
2886: PetscErrorCode DMSubDomainHookRemove(DM global,PetscErrorCode (*ddhook)(DM,DM,void*),PetscErrorCode (*restricthook)(DM,VecScatter,VecScatter,DM,void*),void *ctx)
2887: {
2888: PetscErrorCode ierr;
2889: DMSubDomainHookLink link,*p;
2893: for (p=&global->subdomainhook; *p; p=&(*p)->next) { /* Search the list of current hooks */
2894: if ((*p)->ddhook == ddhook && (*p)->restricthook == restricthook && (*p)->ctx == ctx) {
2895: link = *p;
2896: *p = link->next;
2897: PetscFree(link);
2898: break;
2899: }
2900: }
2901: return(0);
2902: }
2904: /*@
2905: DMSubDomainRestrict - restricts user-defined problem data to a block DM by running hooks registered by DMSubDomainHookAdd()
2907: Collective if any hooks are
2909: Input Arguments:
2910: + fine - finer DM to use as a base
2911: . oscatter - scatter from domain global vector filling subdomain global vector with overlap
2912: . gscatter - scatter from domain global vector filling subdomain local vector with ghosts
2913: - coarse - coarer DM to update
2915: Level: developer
2917: .seealso: DMCoarsenHookAdd(), MatRestrict()
2918: @*/
2919: PetscErrorCode DMSubDomainRestrict(DM global,VecScatter oscatter,VecScatter gscatter,DM subdm)
2920: {
2921: PetscErrorCode ierr;
2922: DMSubDomainHookLink link;
2925: for (link=global->subdomainhook; link; link=link->next) {
2926: if (link->restricthook) {
2927: (*link->restricthook)(global,oscatter,gscatter,subdm,link->ctx);
2928: }
2929: }
2930: return(0);
2931: }
2933: /*@
2934: DMGetCoarsenLevel - Get's the number of coarsenings that have generated this DM.
2936: Not Collective
2938: Input Parameter:
2939: . dm - the DM object
2941: Output Parameter:
2942: . level - number of coarsenings
2944: Level: developer
2946: .seealso DMCoarsen(), DMGetRefineLevel(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation()
2948: @*/
2949: PetscErrorCode DMGetCoarsenLevel(DM dm,PetscInt *level)
2950: {
2953: *level = dm->leveldown;
2954: return(0);
2955: }
2957: /*@
2958: DMSetCoarsenLevel - Sets the number of coarsenings that have generated this DM.
2960: Not Collective
2962: Input Parameters:
2963: + dm - the DM object
2964: - level - number of coarsenings
2966: Level: developer
2968: .seealso DMCoarsen(), DMGetCoarsenLevel(), DMGetRefineLevel(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation()
2969: @*/
2970: PetscErrorCode DMSetCoarsenLevel(DM dm,PetscInt level)
2971: {
2974: dm->leveldown = level;
2975: return(0);
2976: }
2980: /*@C
2981: DMRefineHierarchy - Refines a DM object, all levels at once
2983: Collective on DM
2985: Input Parameter:
2986: + dm - the DM object
2987: - nlevels - the number of levels of refinement
2989: Output Parameter:
2990: . dmf - the refined DM hierarchy
2992: Level: developer
2994: .seealso DMCoarsenHierarchy(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation()
2996: @*/
2997: PetscErrorCode DMRefineHierarchy(DM dm,PetscInt nlevels,DM dmf[])
2998: {
3003: if (nlevels < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"nlevels cannot be negative");
3004: if (nlevels == 0) return(0);
3005: if (dm->ops->refinehierarchy) {
3006: (*dm->ops->refinehierarchy)(dm,nlevels,dmf);
3007: } else if (dm->ops->refine) {
3008: PetscInt i;
3010: DMRefine(dm,PetscObjectComm((PetscObject)dm),&dmf[0]);
3011: for (i=1; i<nlevels; i++) {
3012: DMRefine(dmf[i-1],PetscObjectComm((PetscObject)dm),&dmf[i]);
3013: }
3014: } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"No RefineHierarchy for this DM yet");
3015: return(0);
3016: }
3018: /*@C
3019: DMCoarsenHierarchy - Coarsens a DM object, all levels at once
3021: Collective on DM
3023: Input Parameter:
3024: + dm - the DM object
3025: - nlevels - the number of levels of coarsening
3027: Output Parameter:
3028: . dmc - the coarsened DM hierarchy
3030: Level: developer
3032: .seealso DMRefineHierarchy(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation()
3034: @*/
3035: PetscErrorCode DMCoarsenHierarchy(DM dm, PetscInt nlevels, DM dmc[])
3036: {
3041: if (nlevels < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"nlevels cannot be negative");
3042: if (nlevels == 0) return(0);
3044: if (dm->ops->coarsenhierarchy) {
3045: (*dm->ops->coarsenhierarchy)(dm, nlevels, dmc);
3046: } else if (dm->ops->coarsen) {
3047: PetscInt i;
3049: DMCoarsen(dm,PetscObjectComm((PetscObject)dm),&dmc[0]);
3050: for (i=1; i<nlevels; i++) {
3051: DMCoarsen(dmc[i-1],PetscObjectComm((PetscObject)dm),&dmc[i]);
3052: }
3053: } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"No CoarsenHierarchy for this DM yet");
3054: return(0);
3055: }
3057: /*@
3058: DMCreateAggregates - Gets the aggregates that map between
3059: grids associated with two DMs.
3061: Collective on DM
3063: Input Parameters:
3064: + dmc - the coarse grid DM
3065: - dmf - the fine grid DM
3067: Output Parameters:
3068: . rest - the restriction matrix (transpose of the projection matrix)
3070: Level: intermediate
3072: .keywords: interpolation, restriction, multigrid
3074: .seealso: DMRefine(), DMCreateInjection(), DMCreateInterpolation()
3075: @*/
3076: PetscErrorCode DMCreateAggregates(DM dmc, DM dmf, Mat *rest)
3077: {
3083: (*dmc->ops->getaggregates)(dmc, dmf, rest);
3084: return(0);
3085: }
3087: /*@C
3088: DMSetApplicationContextDestroy - Sets a user function that will be called to destroy the application context when the DM is destroyed
3090: Not Collective
3092: Input Parameters:
3093: + dm - the DM object
3094: - destroy - the destroy function
3096: Level: intermediate
3098: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext()
3100: @*/
3101: PetscErrorCode DMSetApplicationContextDestroy(DM dm,PetscErrorCode (*destroy)(void**))
3102: {
3105: dm->ctxdestroy = destroy;
3106: return(0);
3107: }
3109: /*@
3110: DMSetApplicationContext - Set a user context into a DM object
3112: Not Collective
3114: Input Parameters:
3115: + dm - the DM object
3116: - ctx - the user context
3118: Level: intermediate
3120: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext()
3122: @*/
3123: PetscErrorCode DMSetApplicationContext(DM dm,void *ctx)
3124: {
3127: dm->ctx = ctx;
3128: return(0);
3129: }
3131: /*@
3132: DMGetApplicationContext - Gets a user context from a DM object
3134: Not Collective
3136: Input Parameter:
3137: . dm - the DM object
3139: Output Parameter:
3140: . ctx - the user context
3142: Level: intermediate
3144: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext()
3146: @*/
3147: PetscErrorCode DMGetApplicationContext(DM dm,void *ctx)
3148: {
3151: *(void**)ctx = dm->ctx;
3152: return(0);
3153: }
3155: /*@C
3156: DMSetVariableBounds - sets a function to compute the lower and upper bound vectors for SNESVI.
3158: Logically Collective on DM
3160: Input Parameter:
3161: + dm - the DM object
3162: - f - the function that computes variable bounds used by SNESVI (use NULL to cancel a previous function that was set)
3164: Level: intermediate
3166: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext(),
3167: DMSetJacobian()
3169: @*/
3170: PetscErrorCode DMSetVariableBounds(DM dm,PetscErrorCode (*f)(DM,Vec,Vec))
3171: {
3173: dm->ops->computevariablebounds = f;
3174: return(0);
3175: }
3177: /*@
3178: DMHasVariableBounds - does the DM object have a variable bounds function?
3180: Not Collective
3182: Input Parameter:
3183: . dm - the DM object to destroy
3185: Output Parameter:
3186: . flg - PETSC_TRUE if the variable bounds function exists
3188: Level: developer
3190: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext()
3192: @*/
3193: PetscErrorCode DMHasVariableBounds(DM dm,PetscBool *flg)
3194: {
3196: *flg = (dm->ops->computevariablebounds) ? PETSC_TRUE : PETSC_FALSE;
3197: return(0);
3198: }
3200: /*@C
3201: DMComputeVariableBounds - compute variable bounds used by SNESVI.
3203: Logically Collective on DM
3205: Input Parameters:
3206: . dm - the DM object
3208: Output parameters:
3209: + xl - lower bound
3210: - xu - upper bound
3212: Level: advanced
3214: Notes:
3215: This is generally not called by users. It calls the function provided by the user with DMSetVariableBounds()
3217: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext()
3219: @*/
3220: PetscErrorCode DMComputeVariableBounds(DM dm, Vec xl, Vec xu)
3221: {
3227: if (dm->ops->computevariablebounds) {
3228: (*dm->ops->computevariablebounds)(dm, xl,xu);
3229: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "This DM is incapable of computing variable bounds.");
3230: return(0);
3231: }
3233: /*@
3234: DMHasColoring - does the DM object have a method of providing a coloring?
3236: Not Collective
3238: Input Parameter:
3239: . dm - the DM object
3241: Output Parameter:
3242: . flg - PETSC_TRUE if the DM has facilities for DMCreateColoring().
3244: Level: developer
3246: .seealso DMHasFunction(), DMCreateColoring()
3248: @*/
3249: PetscErrorCode DMHasColoring(DM dm,PetscBool *flg)
3250: {
3252: *flg = (dm->ops->getcoloring) ? PETSC_TRUE : PETSC_FALSE;
3253: return(0);
3254: }
3256: /*@
3257: DMHasCreateRestriction - does the DM object have a method of providing a restriction?
3259: Not Collective
3261: Input Parameter:
3262: . dm - the DM object
3264: Output Parameter:
3265: . flg - PETSC_TRUE if the DM has facilities for DMCreateRestriction().
3267: Level: developer
3269: .seealso DMHasFunction(), DMCreateRestriction()
3271: @*/
3272: PetscErrorCode DMHasCreateRestriction(DM dm,PetscBool *flg)
3273: {
3275: *flg = (dm->ops->createrestriction) ? PETSC_TRUE : PETSC_FALSE;
3276: return(0);
3277: }
3280: /*@
3281: DMHasCreateInjection - does the DM object have a method of providing an injection?
3283: Not Collective
3285: Input Parameter:
3286: . dm - the DM object
3288: Output Parameter:
3289: . flg - PETSC_TRUE if the DM has facilities for DMCreateInjection().
3291: Level: developer
3293: .seealso DMHasFunction(), DMCreateInjection()
3295: @*/
3296: PetscErrorCode DMHasCreateInjection(DM dm,PetscBool *flg)
3297: {
3300: if (!dm->ops->hascreateinjection) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DMHasCreateInjection not implemented for this type");
3301: (*dm->ops->hascreateinjection)(dm,flg);
3302: return(0);
3303: }
3305: /*@C
3306: DMSetVec - set the vector at which to compute residual, Jacobian and VI bounds, if the problem is nonlinear.
3308: Collective on DM
3310: Input Parameter:
3311: + dm - the DM object
3312: - x - location to compute residual and Jacobian, if NULL is passed to those routines; will be NULL for linear problems.
3314: Level: developer
3316: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext()
3318: @*/
3319: PetscErrorCode DMSetVec(DM dm,Vec x)
3320: {
3324: if (x) {
3325: if (!dm->x) {
3326: DMCreateGlobalVector(dm,&dm->x);
3327: }
3328: VecCopy(x,dm->x);
3329: } else if (dm->x) {
3330: VecDestroy(&dm->x);
3331: }
3332: return(0);
3333: }
3335: PetscFunctionList DMList = NULL;
3336: PetscBool DMRegisterAllCalled = PETSC_FALSE;
3338: /*@C
3339: DMSetType - Builds a DM, for a particular DM implementation.
3341: Collective on DM
3343: Input Parameters:
3344: + dm - The DM object
3345: - method - The name of the DM type
3347: Options Database Key:
3348: . -dm_type <type> - Sets the DM type; use -help for a list of available types
3350: Notes:
3351: See "petsc/include/petscdm.h" for available DM types (for instance, DM1D, DM2D, or DM3D).
3353: Level: intermediate
3355: .keywords: DM, set, type
3356: .seealso: DMGetType(), DMCreate()
3357: @*/
3358: PetscErrorCode DMSetType(DM dm, DMType method)
3359: {
3360: PetscErrorCode (*r)(DM);
3361: PetscBool match;
3366: PetscObjectTypeCompare((PetscObject) dm, method, &match);
3367: if (match) return(0);
3369: DMRegisterAll();
3370: PetscFunctionListFind(DMList,method,&r);
3371: if (!r) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown DM type: %s", method);
3373: if (dm->ops->destroy) {
3374: (*dm->ops->destroy)(dm);
3375: dm->ops->destroy = NULL;
3376: }
3377: (*r)(dm);
3378: PetscObjectChangeTypeName((PetscObject)dm,method);
3379: return(0);
3380: }
3382: /*@C
3383: DMGetType - Gets the DM type name (as a string) from the DM.
3385: Not Collective
3387: Input Parameter:
3388: . dm - The DM
3390: Output Parameter:
3391: . type - The DM type name
3393: Level: intermediate
3395: .keywords: DM, get, type, name
3396: .seealso: DMSetType(), DMCreate()
3397: @*/
3398: PetscErrorCode DMGetType(DM dm, DMType *type)
3399: {
3405: DMRegisterAll();
3406: *type = ((PetscObject)dm)->type_name;
3407: return(0);
3408: }
3410: /*@C
3411: DMConvert - Converts a DM to another DM, either of the same or different type.
3413: Collective on DM
3415: Input Parameters:
3416: + dm - the DM
3417: - newtype - new DM type (use "same" for the same type)
3419: Output Parameter:
3420: . M - pointer to new DM
3422: Notes:
3423: Cannot be used to convert a sequential DM to parallel or parallel to sequential,
3424: the MPI communicator of the generated DM is always the same as the communicator
3425: of the input DM.
3427: Level: intermediate
3429: .seealso: DMCreate()
3430: @*/
3431: PetscErrorCode DMConvert(DM dm, DMType newtype, DM *M)
3432: {
3433: DM B;
3434: char convname[256];
3435: PetscBool sametype/*, issame */;
3442: PetscObjectTypeCompare((PetscObject) dm, newtype, &sametype);
3443: /* PetscStrcmp(newtype, "same", &issame); */
3444: if (sametype) {
3445: *M = dm;
3446: PetscObjectReference((PetscObject) dm);
3447: return(0);
3448: } else {
3449: PetscErrorCode (*conv)(DM, DMType, DM*) = NULL;
3451: /*
3452: Order of precedence:
3453: 1) See if a specialized converter is known to the current DM.
3454: 2) See if a specialized converter is known to the desired DM class.
3455: 3) See if a good general converter is registered for the desired class
3456: 4) See if a good general converter is known for the current matrix.
3457: 5) Use a really basic converter.
3458: */
3460: /* 1) See if a specialized converter is known to the current DM and the desired class */
3461: PetscStrncpy(convname,"DMConvert_",sizeof(convname));
3462: PetscStrlcat(convname,((PetscObject) dm)->type_name,sizeof(convname));
3463: PetscStrlcat(convname,"_",sizeof(convname));
3464: PetscStrlcat(convname,newtype,sizeof(convname));
3465: PetscStrlcat(convname,"_C",sizeof(convname));
3466: PetscObjectQueryFunction((PetscObject)dm,convname,&conv);
3467: if (conv) goto foundconv;
3469: /* 2) See if a specialized converter is known to the desired DM class. */
3470: DMCreate(PetscObjectComm((PetscObject)dm), &B);
3471: DMSetType(B, newtype);
3472: PetscStrncpy(convname,"DMConvert_",sizeof(convname));
3473: PetscStrlcat(convname,((PetscObject) dm)->type_name,sizeof(convname));
3474: PetscStrlcat(convname,"_",sizeof(convname));
3475: PetscStrlcat(convname,newtype,sizeof(convname));
3476: PetscStrlcat(convname,"_C",sizeof(convname));
3477: PetscObjectQueryFunction((PetscObject)B,convname,&conv);
3478: if (conv) {
3479: DMDestroy(&B);
3480: goto foundconv;
3481: }
3483: #if 0
3484: /* 3) See if a good general converter is registered for the desired class */
3485: conv = B->ops->convertfrom;
3486: DMDestroy(&B);
3487: if (conv) goto foundconv;
3489: /* 4) See if a good general converter is known for the current matrix */
3490: if (dm->ops->convert) {
3491: conv = dm->ops->convert;
3492: }
3493: if (conv) goto foundconv;
3494: #endif
3496: /* 5) Use a really basic converter. */
3497: SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "No conversion possible between DM types %s and %s", ((PetscObject) dm)->type_name, newtype);
3499: foundconv:
3500: PetscLogEventBegin(DM_Convert,dm,0,0,0);
3501: (*conv)(dm,newtype,M);
3502: /* Things that are independent of DM type: We should consult DMClone() here */
3503: {
3504: PetscBool isper;
3505: const PetscReal *maxCell, *L;
3506: const DMBoundaryType *bd;
3507: DMGetPeriodicity(dm, &isper, &maxCell, &L, &bd);
3508: DMSetPeriodicity(*M, isper, maxCell, L, bd);
3509: }
3510: PetscLogEventEnd(DM_Convert,dm,0,0,0);
3511: }
3512: PetscObjectStateIncrease((PetscObject) *M);
3513: return(0);
3514: }
3516: /*--------------------------------------------------------------------------------------------------------------------*/
3518: /*@C
3519: DMRegister - Adds a new DM component implementation
3521: Not Collective
3523: Input Parameters:
3524: + name - The name of a new user-defined creation routine
3525: - create_func - The creation routine itself
3527: Notes:
3528: DMRegister() may be called multiple times to add several user-defined DMs
3531: Sample usage:
3532: .vb
3533: DMRegister("my_da", MyDMCreate);
3534: .ve
3536: Then, your DM type can be chosen with the procedural interface via
3537: .vb
3538: DMCreate(MPI_Comm, DM *);
3539: DMSetType(DM,"my_da");
3540: .ve
3541: or at runtime via the option
3542: .vb
3543: -da_type my_da
3544: .ve
3546: Level: advanced
3548: .keywords: DM, register
3549: .seealso: DMRegisterAll(), DMRegisterDestroy()
3551: @*/
3552: PetscErrorCode DMRegister(const char sname[],PetscErrorCode (*function)(DM))
3553: {
3557: DMInitializePackage();
3558: PetscFunctionListAdd(&DMList,sname,function);
3559: return(0);
3560: }
3562: /*@C
3563: DMLoad - Loads a DM that has been stored in binary with DMView().
3565: Collective on PetscViewer
3567: Input Parameters:
3568: + newdm - the newly loaded DM, this needs to have been created with DMCreate() or
3569: some related function before a call to DMLoad().
3570: - viewer - binary file viewer, obtained from PetscViewerBinaryOpen() or
3571: HDF5 file viewer, obtained from PetscViewerHDF5Open()
3573: Level: intermediate
3575: Notes:
3576: The type is determined by the data in the file, any type set into the DM before this call is ignored.
3578: Notes for advanced users:
3579: Most users should not need to know the details of the binary storage
3580: format, since DMLoad() and DMView() completely hide these details.
3581: But for anyone who's interested, the standard binary matrix storage
3582: format is
3583: .vb
3584: has not yet been determined
3585: .ve
3587: .seealso: PetscViewerBinaryOpen(), DMView(), MatLoad(), VecLoad()
3588: @*/
3589: PetscErrorCode DMLoad(DM newdm, PetscViewer viewer)
3590: {
3591: PetscBool isbinary, ishdf5;
3597: PetscViewerCheckReadable(viewer);
3598: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
3599: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);
3600: if (isbinary) {
3601: PetscInt classid;
3602: char type[256];
3604: PetscViewerBinaryRead(viewer,&classid,1,NULL,PETSC_INT);
3605: if (classid != DM_FILE_CLASSID) SETERRQ1(PetscObjectComm((PetscObject)newdm),PETSC_ERR_ARG_WRONG,"Not DM next in file, classid found %d",(int)classid);
3606: PetscViewerBinaryRead(viewer,type,256,NULL,PETSC_CHAR);
3607: DMSetType(newdm, type);
3608: if (newdm->ops->load) {(*newdm->ops->load)(newdm,viewer);}
3609: } else if (ishdf5) {
3610: if (newdm->ops->load) {(*newdm->ops->load)(newdm,viewer);}
3611: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid viewer; open viewer with PetscViewerBinaryOpen() or PetscViewerHDF5Open()");
3612: return(0);
3613: }
3615: /******************************** FEM Support **********************************/
3617: PetscErrorCode DMPrintCellVector(PetscInt c, const char name[], PetscInt len, const PetscScalar x[])
3618: {
3619: PetscInt f;
3623: PetscPrintf(PETSC_COMM_SELF, "Cell %D Element %s\n", c, name);
3624: for (f = 0; f < len; ++f) {
3625: PetscPrintf(PETSC_COMM_SELF, " | %g |\n", (double)PetscRealPart(x[f]));
3626: }
3627: return(0);
3628: }
3630: PetscErrorCode DMPrintCellMatrix(PetscInt c, const char name[], PetscInt rows, PetscInt cols, const PetscScalar A[])
3631: {
3632: PetscInt f, g;
3636: PetscPrintf(PETSC_COMM_SELF, "Cell %D Element %s\n", c, name);
3637: for (f = 0; f < rows; ++f) {
3638: PetscPrintf(PETSC_COMM_SELF, " |");
3639: for (g = 0; g < cols; ++g) {
3640: PetscPrintf(PETSC_COMM_SELF, " % 9.5g", PetscRealPart(A[f*cols+g]));
3641: }
3642: PetscPrintf(PETSC_COMM_SELF, " |\n");
3643: }
3644: return(0);
3645: }
3647: PetscErrorCode DMPrintLocalVec(DM dm, const char name[], PetscReal tol, Vec X)
3648: {
3649: PetscInt localSize, bs;
3650: PetscMPIInt size;
3651: Vec x, xglob;
3652: const PetscScalar *xarray;
3653: PetscErrorCode ierr;
3656: MPI_Comm_size(PetscObjectComm((PetscObject) dm),&size);
3657: VecDuplicate(X, &x);
3658: VecCopy(X, x);
3659: VecChop(x, tol);
3660: PetscPrintf(PetscObjectComm((PetscObject) dm),"%s:\n",name);
3661: if (size > 1) {
3662: VecGetLocalSize(x,&localSize);
3663: VecGetArrayRead(x,&xarray);
3664: VecGetBlockSize(x,&bs);
3665: VecCreateMPIWithArray(PetscObjectComm((PetscObject) dm),bs,localSize,PETSC_DETERMINE,xarray,&xglob);
3666: } else {
3667: xglob = x;
3668: }
3669: VecView(xglob,PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject) dm)));
3670: if (size > 1) {
3671: VecDestroy(&xglob);
3672: VecRestoreArrayRead(x,&xarray);
3673: }
3674: VecDestroy(&x);
3675: return(0);
3676: }
3678: /*@
3679: DMGetSection - Get the PetscSection encoding the local data layout for the DM.
3681: Input Parameter:
3682: . dm - The DM
3684: Output Parameter:
3685: . section - The PetscSection
3687: Options Database Keys:
3688: . -dm_petscsection_view - View the Section created by the DM
3690: Level: intermediate
3692: Note: This gets a borrowed reference, so the user should not destroy this PetscSection.
3694: .seealso: DMSetSection(), DMGetGlobalSection()
3695: @*/
3696: PetscErrorCode DMGetSection(DM dm, PetscSection *section)
3697: {
3703: if (!dm->defaultSection && dm->ops->createdefaultsection) {
3704: PetscInt d;
3706: if (dm->setfromoptionscalled) for (d = 0; d < dm->Nds; ++d) {PetscDSSetFromOptions(dm->probs[d].ds);}
3707: (*dm->ops->createdefaultsection)(dm);
3708: if (dm->defaultSection) {PetscObjectViewFromOptions((PetscObject) dm->defaultSection, NULL, "-dm_petscsection_view");}
3709: }
3710: *section = dm->defaultSection;
3711: return(0);
3712: }
3714: /*@
3715: DMSetSection - Set the PetscSection encoding the local data layout for the DM.
3717: Input Parameters:
3718: + dm - The DM
3719: - section - The PetscSection
3721: Level: intermediate
3723: Note: Any existing Section will be destroyed
3725: .seealso: DMSetSection(), DMGetGlobalSection()
3726: @*/
3727: PetscErrorCode DMSetSection(DM dm, PetscSection section)
3728: {
3729: PetscInt numFields = 0;
3730: PetscInt f;
3735: if (section) {
3737: PetscObjectReference((PetscObject)section);
3738: }
3739: PetscSectionDestroy(&dm->defaultSection);
3740: dm->defaultSection = section;
3741: if (section) {PetscSectionGetNumFields(dm->defaultSection, &numFields);}
3742: if (numFields) {
3743: DMSetNumFields(dm, numFields);
3744: for (f = 0; f < numFields; ++f) {
3745: PetscObject disc;
3746: const char *name;
3748: PetscSectionGetFieldName(dm->defaultSection, f, &name);
3749: DMGetField(dm, f, NULL, &disc);
3750: PetscObjectSetName(disc, name);
3751: }
3752: }
3753: /* The global section will be rebuilt in the next call to DMGetGlobalSection(). */
3754: PetscSectionDestroy(&dm->defaultGlobalSection);
3755: return(0);
3756: }
3758: /*@
3759: DMGetDefaultConstraints - Get the PetscSection and Mat the specify the local constraint interpolation. See DMSetDefaultConstraints() for a description of the purpose of constraint interpolation.
3761: not collective
3763: Input Parameter:
3764: . dm - The DM
3766: Output Parameter:
3767: + 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.
3768: - 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.
3770: Level: advanced
3772: Note: This gets borrowed references, so the user should not destroy the PetscSection or the Mat.
3774: .seealso: DMSetDefaultConstraints()
3775: @*/
3776: PetscErrorCode DMGetDefaultConstraints(DM dm, PetscSection *section, Mat *mat)
3777: {
3782: if (!dm->defaultConstraintSection && !dm->defaultConstraintMat && dm->ops->createdefaultconstraints) {(*dm->ops->createdefaultconstraints)(dm);}
3783: if (section) {*section = dm->defaultConstraintSection;}
3784: if (mat) {*mat = dm->defaultConstraintMat;}
3785: return(0);
3786: }
3788: /*@
3789: DMSetDefaultConstraints - Set the PetscSection and Mat the specify the local constraint interpolation.
3791: 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().
3793: 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.
3795: collective on dm
3797: Input Parameters:
3798: + dm - The DM
3799: + 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).
3800: - 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).
3802: Level: advanced
3804: Note: This increments the references of the PetscSection and the Mat, so they user can destroy them
3806: .seealso: DMGetDefaultConstraints()
3807: @*/
3808: PetscErrorCode DMSetDefaultConstraints(DM dm, PetscSection section, Mat mat)
3809: {
3810: PetscMPIInt result;
3815: if (section) {
3817: MPI_Comm_compare(PETSC_COMM_SELF,PetscObjectComm((PetscObject)section),&result);
3818: if (result != MPI_CONGRUENT && result != MPI_IDENT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMECOMM,"constraint section must have local communicator");
3819: }
3820: if (mat) {
3822: MPI_Comm_compare(PETSC_COMM_SELF,PetscObjectComm((PetscObject)mat),&result);
3823: if (result != MPI_CONGRUENT && result != MPI_IDENT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMECOMM,"constraint matrix must have local communicator");
3824: }
3825: PetscObjectReference((PetscObject)section);
3826: PetscSectionDestroy(&dm->defaultConstraintSection);
3827: dm->defaultConstraintSection = section;
3828: PetscObjectReference((PetscObject)mat);
3829: MatDestroy(&dm->defaultConstraintMat);
3830: dm->defaultConstraintMat = mat;
3831: return(0);
3832: }
3834: #if defined(PETSC_USE_DEBUG)
3835: /*
3836: DMDefaultSectionCheckConsistency - Check the consistentcy of the global and local sections.
3838: Input Parameters:
3839: + dm - The DM
3840: . localSection - PetscSection describing the local data layout
3841: - globalSection - PetscSection describing the global data layout
3843: Level: intermediate
3845: .seealso: DMGetDefaultSF(), DMSetDefaultSF()
3846: */
3847: static PetscErrorCode DMDefaultSectionCheckConsistency_Internal(DM dm, PetscSection localSection, PetscSection globalSection)
3848: {
3849: MPI_Comm comm;
3850: PetscLayout layout;
3851: const PetscInt *ranges;
3852: PetscInt pStart, pEnd, p, nroots;
3853: PetscMPIInt size, rank;
3854: PetscBool valid = PETSC_TRUE, gvalid;
3855: PetscErrorCode ierr;
3858: PetscObjectGetComm((PetscObject)dm,&comm);
3860: MPI_Comm_size(comm, &size);
3861: MPI_Comm_rank(comm, &rank);
3862: PetscSectionGetChart(globalSection, &pStart, &pEnd);
3863: PetscSectionGetConstrainedStorageSize(globalSection, &nroots);
3864: PetscLayoutCreate(comm, &layout);
3865: PetscLayoutSetBlockSize(layout, 1);
3866: PetscLayoutSetLocalSize(layout, nroots);
3867: PetscLayoutSetUp(layout);
3868: PetscLayoutGetRanges(layout, &ranges);
3869: for (p = pStart; p < pEnd; ++p) {
3870: PetscInt dof, cdof, off, gdof, gcdof, goff, gsize, d;
3872: PetscSectionGetDof(localSection, p, &dof);
3873: PetscSectionGetOffset(localSection, p, &off);
3874: PetscSectionGetConstraintDof(localSection, p, &cdof);
3875: PetscSectionGetDof(globalSection, p, &gdof);
3876: PetscSectionGetConstraintDof(globalSection, p, &gcdof);
3877: PetscSectionGetOffset(globalSection, p, &goff);
3878: if (!gdof) continue; /* Censored point */
3879: 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;}
3880: 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;}
3881: if (gdof < 0) {
3882: gsize = gdof < 0 ? -(gdof+1)-gcdof : gdof-gcdof;
3883: for (d = 0; d < gsize; ++d) {
3884: PetscInt offset = -(goff+1) + d, r;
3886: PetscFindInt(offset,size+1,ranges,&r);
3887: if (r < 0) r = -(r+2);
3888: 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;}
3889: }
3890: }
3891: }
3892: PetscLayoutDestroy(&layout);
3893: PetscSynchronizedFlush(comm, NULL);
3894: MPIU_Allreduce(&valid, &gvalid, 1, MPIU_BOOL, MPI_LAND, comm);
3895: if (!gvalid) {
3896: DMView(dm, NULL);
3897: SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Inconsistent local and global sections");
3898: }
3899: return(0);
3900: }
3901: #endif
3903: /*@
3904: DMGetGlobalSection - Get the PetscSection encoding the global data layout for the DM.
3906: Collective on DM
3908: Input Parameter:
3909: . dm - The DM
3911: Output Parameter:
3912: . section - The PetscSection
3914: Level: intermediate
3916: Note: This gets a borrowed reference, so the user should not destroy this PetscSection.
3918: .seealso: DMSetSection(), DMGetSection()
3919: @*/
3920: PetscErrorCode DMGetGlobalSection(DM dm, PetscSection *section)
3921: {
3927: if (!dm->defaultGlobalSection) {
3928: PetscSection s;
3930: DMGetSection(dm, &s);
3931: if (!s) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONGSTATE, "DM must have a default PetscSection in order to create a global PetscSection");
3932: if (!dm->sf) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DM must have a point PetscSF in order to create a global PetscSection");
3933: PetscSectionCreateGlobalSection(s, dm->sf, PETSC_FALSE, PETSC_FALSE, &dm->defaultGlobalSection);
3934: PetscLayoutDestroy(&dm->map);
3935: PetscSectionGetValueLayout(PetscObjectComm((PetscObject)dm), dm->defaultGlobalSection, &dm->map);
3936: PetscSectionViewFromOptions(dm->defaultGlobalSection, NULL, "-global_section_view");
3937: }
3938: *section = dm->defaultGlobalSection;
3939: return(0);
3940: }
3942: /*@
3943: DMSetGlobalSection - Set the PetscSection encoding the global data layout for the DM.
3945: Input Parameters:
3946: + dm - The DM
3947: - section - The PetscSection, or NULL
3949: Level: intermediate
3951: Note: Any existing Section will be destroyed
3953: .seealso: DMGetGlobalSection(), DMSetSection()
3954: @*/
3955: PetscErrorCode DMSetGlobalSection(DM dm, PetscSection section)
3956: {
3962: PetscObjectReference((PetscObject)section);
3963: PetscSectionDestroy(&dm->defaultGlobalSection);
3964: dm->defaultGlobalSection = section;
3965: #if defined(PETSC_USE_DEBUG)
3966: if (section) {DMDefaultSectionCheckConsistency_Internal(dm, dm->defaultSection, section);}
3967: #endif
3968: return(0);
3969: }
3971: /*@
3972: DMGetDefaultSF - Get the PetscSF encoding the parallel dof overlap for the DM. If it has not been set,
3973: it is created from the default PetscSection layouts in the DM.
3975: Input Parameter:
3976: . dm - The DM
3978: Output Parameter:
3979: . sf - The PetscSF
3981: Level: intermediate
3983: Note: This gets a borrowed reference, so the user should not destroy this PetscSF.
3985: .seealso: DMSetDefaultSF(), DMCreateDefaultSF()
3986: @*/
3987: PetscErrorCode DMGetDefaultSF(DM dm, PetscSF *sf)
3988: {
3989: PetscInt nroots;
3995: if (!dm->defaultSF) {
3996: PetscSFCreate(PetscObjectComm((PetscObject)dm),&dm->defaultSF);
3997: }
3998: PetscSFGetGraph(dm->defaultSF, &nroots, NULL, NULL, NULL);
3999: if (nroots < 0) {
4000: PetscSection section, gSection;
4002: DMGetSection(dm, §ion);
4003: if (section) {
4004: DMGetGlobalSection(dm, &gSection);
4005: DMCreateDefaultSF(dm, section, gSection);
4006: } else {
4007: *sf = NULL;
4008: return(0);
4009: }
4010: }
4011: *sf = dm->defaultSF;
4012: return(0);
4013: }
4015: /*@
4016: DMSetDefaultSF - Set the PetscSF encoding the parallel dof overlap for the DM
4018: Input Parameters:
4019: + dm - The DM
4020: - sf - The PetscSF
4022: Level: intermediate
4024: Note: Any previous SF is destroyed
4026: .seealso: DMGetDefaultSF(), DMCreateDefaultSF()
4027: @*/
4028: PetscErrorCode DMSetDefaultSF(DM dm, PetscSF sf)
4029: {
4034: if (sf) {
4036: PetscObjectReference((PetscObject) sf);
4037: }
4038: PetscSFDestroy(&dm->defaultSF);
4039: dm->defaultSF = sf;
4040: return(0);
4041: }
4043: /*@C
4044: DMCreateDefaultSF - Create the PetscSF encoding the parallel dof overlap for the DM based upon the PetscSections
4045: describing the data layout.
4047: Input Parameters:
4048: + dm - The DM
4049: . localSection - PetscSection describing the local data layout
4050: - globalSection - PetscSection describing the global data layout
4052: Level: intermediate
4054: .seealso: DMGetDefaultSF(), DMSetDefaultSF()
4055: @*/
4056: PetscErrorCode DMCreateDefaultSF(DM dm, PetscSection localSection, PetscSection globalSection)
4057: {
4058: MPI_Comm comm;
4059: PetscLayout layout;
4060: const PetscInt *ranges;
4061: PetscInt *local;
4062: PetscSFNode *remote;
4063: PetscInt pStart, pEnd, p, nroots, nleaves = 0, l;
4064: PetscMPIInt size, rank;
4069: PetscObjectGetComm((PetscObject)dm,&comm);
4070: MPI_Comm_size(comm, &size);
4071: MPI_Comm_rank(comm, &rank);
4072: PetscSectionGetChart(globalSection, &pStart, &pEnd);
4073: PetscSectionGetConstrainedStorageSize(globalSection, &nroots);
4074: PetscLayoutCreate(comm, &layout);
4075: PetscLayoutSetBlockSize(layout, 1);
4076: PetscLayoutSetLocalSize(layout, nroots);
4077: PetscLayoutSetUp(layout);
4078: PetscLayoutGetRanges(layout, &ranges);
4079: for (p = pStart; p < pEnd; ++p) {
4080: PetscInt gdof, gcdof;
4082: PetscSectionGetDof(globalSection, p, &gdof);
4083: PetscSectionGetConstraintDof(globalSection, p, &gcdof);
4084: 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));
4085: nleaves += gdof < 0 ? -(gdof+1)-gcdof : gdof-gcdof;
4086: }
4087: PetscMalloc1(nleaves, &local);
4088: PetscMalloc1(nleaves, &remote);
4089: for (p = pStart, l = 0; p < pEnd; ++p) {
4090: const PetscInt *cind;
4091: PetscInt dof, cdof, off, gdof, gcdof, goff, gsize, d, c;
4093: PetscSectionGetDof(localSection, p, &dof);
4094: PetscSectionGetOffset(localSection, p, &off);
4095: PetscSectionGetConstraintDof(localSection, p, &cdof);
4096: PetscSectionGetConstraintIndices(localSection, p, &cind);
4097: PetscSectionGetDof(globalSection, p, &gdof);
4098: PetscSectionGetConstraintDof(globalSection, p, &gcdof);
4099: PetscSectionGetOffset(globalSection, p, &goff);
4100: if (!gdof) continue; /* Censored point */
4101: gsize = gdof < 0 ? -(gdof+1)-gcdof : gdof-gcdof;
4102: if (gsize != dof-cdof) {
4103: 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);
4104: cdof = 0; /* Ignore constraints */
4105: }
4106: for (d = 0, c = 0; d < dof; ++d) {
4107: if ((c < cdof) && (cind[c] == d)) {++c; continue;}
4108: local[l+d-c] = off+d;
4109: }
4110: if (gdof < 0) {
4111: for (d = 0; d < gsize; ++d, ++l) {
4112: PetscInt offset = -(goff+1) + d, r;
4114: PetscFindInt(offset,size+1,ranges,&r);
4115: if (r < 0) r = -(r+2);
4116: 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);
4117: remote[l].rank = r;
4118: remote[l].index = offset - ranges[r];
4119: }
4120: } else {
4121: for (d = 0; d < gsize; ++d, ++l) {
4122: remote[l].rank = rank;
4123: remote[l].index = goff+d - ranges[rank];
4124: }
4125: }
4126: }
4127: if (l != nleaves) SETERRQ2(comm, PETSC_ERR_PLIB, "Iteration error, l %d != nleaves %d", l, nleaves);
4128: PetscLayoutDestroy(&layout);
4129: PetscSFSetGraph(dm->defaultSF, nroots, nleaves, local, PETSC_OWN_POINTER, remote, PETSC_OWN_POINTER);
4130: return(0);
4131: }
4133: /*@
4134: DMGetPointSF - Get the PetscSF encoding the parallel section point overlap for the DM.
4136: Input Parameter:
4137: . dm - The DM
4139: Output Parameter:
4140: . sf - The PetscSF
4142: Level: intermediate
4144: Note: This gets a borrowed reference, so the user should not destroy this PetscSF.
4146: .seealso: DMSetPointSF(), DMGetDefaultSF(), DMSetDefaultSF(), DMCreateDefaultSF()
4147: @*/
4148: PetscErrorCode DMGetPointSF(DM dm, PetscSF *sf)
4149: {
4153: *sf = dm->sf;
4154: return(0);
4155: }
4157: /*@
4158: DMSetPointSF - Set the PetscSF encoding the parallel section point overlap for the DM.
4160: Input Parameters:
4161: + dm - The DM
4162: - sf - The PetscSF
4164: Level: intermediate
4166: .seealso: DMGetPointSF(), DMGetDefaultSF(), DMSetDefaultSF(), DMCreateDefaultSF()
4167: @*/
4168: PetscErrorCode DMSetPointSF(DM dm, PetscSF sf)
4169: {
4174: if (sf) {
4176: PetscObjectReference((PetscObject) sf);
4177: }
4178: PetscSFDestroy(&dm->sf);
4179: dm->sf = sf;
4180: return(0);
4181: }
4183: static PetscErrorCode DMSetDefaultAdjacency_Private(DM dm, PetscInt f, PetscObject disc)
4184: {
4185: PetscClassId id;
4189: PetscObjectGetClassId(disc, &id);
4190: if (id == PETSCFE_CLASSID) {
4191: DMSetAdjacency(dm, f, PETSC_FALSE, PETSC_TRUE);
4192: } else if (id == PETSCFV_CLASSID) {
4193: DMSetAdjacency(dm, f, PETSC_TRUE, PETSC_FALSE);
4194: } else {
4195: DMSetAdjacency(dm, f, PETSC_FALSE, PETSC_TRUE);
4196: }
4197: return(0);
4198: }
4200: static PetscErrorCode DMFieldEnlarge_Static(DM dm, PetscInt NfNew)
4201: {
4202: RegionField *tmpr;
4203: PetscInt Nf = dm->Nf, f;
4207: if (Nf >= NfNew) return(0);
4208: PetscMalloc1(NfNew, &tmpr);
4209: for (f = 0; f < Nf; ++f) tmpr[f] = dm->fields[f];
4210: for (f = Nf; f < NfNew; ++f) {tmpr[f].disc = NULL; tmpr[f].label = NULL;}
4211: PetscFree(dm->fields);
4212: dm->Nf = NfNew;
4213: dm->fields = tmpr;
4214: return(0);
4215: }
4217: /*@
4218: DMClearFields - Remove all fields from the DM
4220: Logically collective on DM
4222: Input Parameter:
4223: . dm - The DM
4225: Level: intermediate
4227: .seealso: DMGetNumFields(), DMSetNumFields(), DMSetField()
4228: @*/
4229: PetscErrorCode DMClearFields(DM dm)
4230: {
4231: PetscInt f;
4236: for (f = 0; f < dm->Nf; ++f) {
4237: PetscObjectDestroy(&dm->fields[f].disc);
4238: DMLabelDestroy(&dm->fields[f].label);
4239: }
4240: PetscFree(dm->fields);
4241: dm->fields = NULL;
4242: dm->Nf = 0;
4243: return(0);
4244: }
4246: /*@
4247: DMGetNumFields - Get the number of fields in the DM
4249: Not collective
4251: Input Parameter:
4252: . dm - The DM
4254: Output Parameter:
4255: . Nf - The number of fields
4257: Level: intermediate
4259: .seealso: DMSetNumFields(), DMSetField()
4260: @*/
4261: PetscErrorCode DMGetNumFields(DM dm, PetscInt *numFields)
4262: {
4266: *numFields = dm->Nf;
4267: return(0);
4268: }
4270: /*@
4271: DMSetNumFields - Set the number of fields in the DM
4273: Logically collective on DM
4275: Input Parameters:
4276: + dm - The DM
4277: - Nf - The number of fields
4279: Level: intermediate
4281: .seealso: DMGetNumFields(), DMSetField()
4282: @*/
4283: PetscErrorCode DMSetNumFields(DM dm, PetscInt numFields)
4284: {
4285: PetscInt Nf, f;
4290: DMGetNumFields(dm, &Nf);
4291: for (f = Nf; f < numFields; ++f) {
4292: PetscContainer obj;
4294: PetscContainerCreate(PetscObjectComm((PetscObject) dm), &obj);
4295: DMAddField(dm, NULL, (PetscObject) obj);
4296: PetscContainerDestroy(&obj);
4297: }
4298: return(0);
4299: }
4301: /*@
4302: DMGetField - Return the discretization object for a given DM field
4304: Not collective
4306: Input Parameters:
4307: + dm - The DM
4308: - f - The field number
4310: Output Parameters:
4311: + label - The label indicating the support of the field, or NULL for the entire mesh
4312: - field - The discretization object
4314: Level: intermediate
4316: .seealso: DMAddField(), DMSetField()
4317: @*/
4318: PetscErrorCode DMGetField(DM dm, PetscInt f, DMLabel *label, PetscObject *field)
4319: {
4323: 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);
4324: if (label) *label = dm->fields[f].label;
4325: if (field) *field = dm->fields[f].disc;
4326: return(0);
4327: }
4329: /*@
4330: DMSetField - Set the discretization object for a given DM field
4332: Logically collective on DM
4334: Input Parameters:
4335: + dm - The DM
4336: . f - The field number
4337: . label - The label indicating the support of the field, or NULL for the entire mesh
4338: - field - The discretization object
4340: Level: intermediate
4342: .seealso: DMAddField(), DMGetField()
4343: @*/
4344: PetscErrorCode DMSetField(DM dm, PetscInt f, DMLabel label, PetscObject field)
4345: {
4352: if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
4353: DMFieldEnlarge_Static(dm, f+1);
4354: DMLabelDestroy(&dm->fields[f].label);
4355: PetscObjectDestroy(&dm->fields[f].disc);
4356: dm->fields[f].label = label;
4357: dm->fields[f].disc = field;
4358: PetscObjectReference((PetscObject) label);
4359: PetscObjectReference((PetscObject) field);
4360: DMSetDefaultAdjacency_Private(dm, f, field);
4361: DMClearDS(dm);
4362: return(0);
4363: }
4365: /*@
4366: DMAddField - Add the discretization object for the given DM field
4368: Logically collective on DM
4370: Input Parameters:
4371: + dm - The DM
4372: . label - The label indicating the support of the field, or NULL for the entire mesh
4373: - field - The discretization object
4375: Level: intermediate
4377: .seealso: DMSetField(), DMGetField()
4378: @*/
4379: PetscErrorCode DMAddField(DM dm, DMLabel label, PetscObject field)
4380: {
4381: PetscInt Nf = dm->Nf;
4388: DMFieldEnlarge_Static(dm, Nf+1);
4389: dm->fields[Nf].label = label;
4390: dm->fields[Nf].disc = field;
4391: PetscObjectReference((PetscObject) label);
4392: PetscObjectReference((PetscObject) field);
4393: DMSetDefaultAdjacency_Private(dm, Nf, field);
4394: DMClearDS(dm);
4395: return(0);
4396: }
4398: /*@
4399: DMCopyFields - Copy the discretizations for the DM into another DM
4401: Collective on DM
4403: Input Parameter:
4404: . dm - The DM
4406: Output Parameter:
4407: . newdm - The DM
4409: Level: advanced
4411: .seealso: DMGetField(), DMSetField(), DMAddField(), DMCopyDS(), DMGetDS(), DMGetCellDS()
4412: @*/
4413: PetscErrorCode DMCopyFields(DM dm, DM newdm)
4414: {
4415: PetscInt Nf, f;
4419: if (dm == newdm) return(0);
4420: DMGetNumFields(dm, &Nf);
4421: DMClearFields(newdm);
4422: for (f = 0; f < Nf; ++f) {
4423: DMLabel label;
4424: PetscObject field;
4425: PetscBool useCone, useClosure;
4427: DMGetField(dm, f, &label, &field);
4428: DMSetField(newdm, f, label, field);
4429: DMGetAdjacency(dm, f, &useCone, &useClosure);
4430: DMSetAdjacency(newdm, f, useCone, useClosure);
4431: }
4432: return(0);
4433: }
4435: /*@
4436: DMGetAdjacency - Returns the flags for determining variable influence
4438: Not collective
4440: Input Parameters:
4441: + dm - The DM object
4442: - f - The field number, or PETSC_DEFAULT for the default adjacency
4444: Output Parameter:
4445: + useCone - Flag for variable influence starting with the cone operation
4446: - useClosure - Flag for variable influence using transitive closure
4448: Notes:
4449: $ FEM: Two points p and q are adjacent if q \in closure(star(p)), useCone = PETSC_FALSE, useClosure = PETSC_TRUE
4450: $ FVM: Two points p and q are adjacent if q \in support(p+cone(p)), useCone = PETSC_TRUE, useClosure = PETSC_FALSE
4451: $ FVM++: Two points p and q are adjacent if q \in star(closure(p)), useCone = PETSC_TRUE, useClosure = PETSC_TRUE
4452: Further explanation can be found in the User's Manual Section on the Influence of Variables on One Another.
4454: Level: developer
4456: .seealso: DMSetAdjacency(), DMGetField(), DMSetField()
4457: @*/
4458: PetscErrorCode DMGetAdjacency(DM dm, PetscInt f, PetscBool *useCone, PetscBool *useClosure)
4459: {
4464: if (f < 0) {
4465: if (useCone) *useCone = dm->adjacency[0];
4466: if (useClosure) *useClosure = dm->adjacency[1];
4467: } else {
4468: PetscInt Nf;
4471: DMGetNumFields(dm, &Nf);
4472: if (f >= Nf) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, Nf);
4473: if (useCone) *useCone = dm->fields[f].adjacency[0];
4474: if (useClosure) *useClosure = dm->fields[f].adjacency[1];
4475: }
4476: return(0);
4477: }
4479: /*@
4480: DMSetAdjacency - Set the flags for determining variable influence
4482: Not collective
4484: Input Parameters:
4485: + dm - The DM object
4486: . f - The field number
4487: . useCone - Flag for variable influence starting with the cone operation
4488: - useClosure - Flag for variable influence using transitive closure
4490: Notes:
4491: $ FEM: Two points p and q are adjacent if q \in closure(star(p)), useCone = PETSC_FALSE, useClosure = PETSC_TRUE
4492: $ FVM: Two points p and q are adjacent if q \in support(p+cone(p)), useCone = PETSC_TRUE, useClosure = PETSC_FALSE
4493: $ FVM++: Two points p and q are adjacent if q \in star(closure(p)), useCone = PETSC_TRUE, useClosure = PETSC_TRUE
4494: Further explanation can be found in the User's Manual Section on the Influence of Variables on One Another.
4496: Level: developer
4498: .seealso: DMGetAdjacency(), DMGetField(), DMSetField()
4499: @*/
4500: PetscErrorCode DMSetAdjacency(DM dm, PetscInt f, PetscBool useCone, PetscBool useClosure)
4501: {
4504: if (f < 0) {
4505: dm->adjacency[0] = useCone;
4506: dm->adjacency[1] = useClosure;
4507: } else {
4508: PetscInt Nf;
4511: DMGetNumFields(dm, &Nf);
4512: if (f >= Nf) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, Nf);
4513: dm->fields[f].adjacency[0] = useCone;
4514: dm->fields[f].adjacency[1] = useClosure;
4515: }
4516: return(0);
4517: }
4519: /*@
4520: DMGetBasicAdjacency - Returns the flags for determining variable influence, using either the default or field 0 if it is defined
4522: Not collective
4524: Input Parameters:
4525: . dm - The DM object
4527: Output Parameter:
4528: + useCone - Flag for variable influence starting with the cone operation
4529: - useClosure - Flag for variable influence using transitive closure
4531: Notes:
4532: $ FEM: Two points p and q are adjacent if q \in closure(star(p)), useCone = PETSC_FALSE, useClosure = PETSC_TRUE
4533: $ FVM: Two points p and q are adjacent if q \in support(p+cone(p)), useCone = PETSC_TRUE, useClosure = PETSC_FALSE
4534: $ FVM++: Two points p and q are adjacent if q \in star(closure(p)), useCone = PETSC_TRUE, useClosure = PETSC_TRUE
4536: Level: developer
4538: .seealso: DMSetBasicAdjacency(), DMGetField(), DMSetField()
4539: @*/
4540: PetscErrorCode DMGetBasicAdjacency(DM dm, PetscBool *useCone, PetscBool *useClosure)
4541: {
4542: PetscInt Nf;
4549: DMGetNumFields(dm, &Nf);
4550: if (!Nf) {
4551: DMGetAdjacency(dm, PETSC_DEFAULT, useCone, useClosure);
4552: } else {
4553: DMGetAdjacency(dm, 0, useCone, useClosure);
4554: }
4555: return(0);
4556: }
4558: /*@
4559: DMSetBasicAdjacency - Set the flags for determining variable influence, using either the default or field 0 if it is defined
4561: Not collective
4563: Input Parameters:
4564: + dm - The DM object
4565: . useCone - Flag for variable influence starting with the cone operation
4566: - useClosure - Flag for variable influence using transitive closure
4568: Notes:
4569: $ FEM: Two points p and q are adjacent if q \in closure(star(p)), useCone = PETSC_FALSE, useClosure = PETSC_TRUE
4570: $ FVM: Two points p and q are adjacent if q \in support(p+cone(p)), useCone = PETSC_TRUE, useClosure = PETSC_FALSE
4571: $ FVM++: Two points p and q are adjacent if q \in star(closure(p)), useCone = PETSC_TRUE, useClosure = PETSC_TRUE
4573: Level: developer
4575: .seealso: DMGetBasicAdjacency(), DMGetField(), DMSetField()
4576: @*/
4577: PetscErrorCode DMSetBasicAdjacency(DM dm, PetscBool useCone, PetscBool useClosure)
4578: {
4579: PetscInt Nf;
4584: DMGetNumFields(dm, &Nf);
4585: if (!Nf) {
4586: DMSetAdjacency(dm, PETSC_DEFAULT, useCone, useClosure);
4587: } else {
4588: DMSetAdjacency(dm, 0, useCone, useClosure);
4589: }
4590: return(0);
4591: }
4593: static PetscErrorCode DMDSEnlarge_Static(DM dm, PetscInt NdsNew)
4594: {
4595: DMSpace *tmpd;
4596: PetscInt Nds = dm->Nds, s;
4600: if (Nds >= NdsNew) return(0);
4601: PetscMalloc1(NdsNew, &tmpd);
4602: for (s = 0; s < Nds; ++s) tmpd[s] = dm->probs[s];
4603: for (s = Nds; s < NdsNew; ++s) {tmpd[s].ds = NULL; tmpd[s].label = NULL;}
4604: PetscFree(dm->probs);
4605: dm->Nds = NdsNew;
4606: dm->probs = tmpd;
4607: return(0);
4608: }
4610: /*@
4611: DMGetNumDS - Get the number of discrete systems in the DM
4613: Not collective
4615: Input Parameter:
4616: . dm - The DM
4618: Output Parameter:
4619: . Nds - The number of PetscDS objects
4621: Level: intermediate
4623: .seealso: DMGetDS(), DMGetCellDS()
4624: @*/
4625: PetscErrorCode DMGetNumDS(DM dm, PetscInt *Nds)
4626: {
4630: *Nds = dm->Nds;
4631: return(0);
4632: }
4634: /*@
4635: DMClearDS - Remove all discrete systems from the DM
4637: Logically collective on DM
4639: Input Parameter:
4640: . dm - The DM
4642: Level: intermediate
4644: .seealso: DMGetNumDS(), DMGetDS(), DMSetField()
4645: @*/
4646: PetscErrorCode DMClearDS(DM dm)
4647: {
4648: PetscInt s;
4653: for (s = 0; s < dm->Nds; ++s) {
4654: PetscDSDestroy(&dm->probs[s].ds);
4655: DMLabelDestroy(&dm->probs[s].label);
4656: }
4657: PetscFree(dm->probs);
4658: dm->probs = NULL;
4659: dm->Nds = 0;
4660: return(0);
4661: }
4663: /*@
4664: DMGetDS - Get the default PetscDS
4666: Not collective
4668: Input Parameter:
4669: . dm - The DM
4671: Output Parameter:
4672: . prob - The default PetscDS
4674: Level: intermediate
4676: .seealso: DMGetCellDS(), DMGetRegionDS()
4677: @*/
4678: PetscErrorCode DMGetDS(DM dm, PetscDS *prob)
4679: {
4685: if (dm->Nds <= 0) {
4686: PetscDS ds;
4688: PetscDSCreate(PetscObjectComm((PetscObject) dm), &ds);
4689: DMSetRegionDS(dm, NULL, ds);
4690: PetscDSDestroy(&ds);
4691: }
4692: *prob = dm->probs[0].ds;
4693: return(0);
4694: }
4696: /*@
4697: DMGetCellDS - Get the PetscDS defined on a given cell
4699: Not collective
4701: Input Parameters:
4702: + dm - The DM
4703: - point - Cell for the DS
4705: Output Parameter:
4706: . prob - The PetscDS defined on the given cell
4708: Level: developer
4710: .seealso: DMGetDS(), DMSetRegionDS()
4711: @*/
4712: PetscErrorCode DMGetCellDS(DM dm, PetscInt point, PetscDS *prob)
4713: {
4714: PetscDS probDef = NULL;
4715: PetscInt s;
4721: *prob = NULL;
4722: for (s = 0; s < dm->Nds; ++s) {
4723: PetscInt val;
4725: if (!dm->probs[s].label) {probDef = dm->probs[s].ds;}
4726: else {
4727: DMLabelGetValue(dm->probs[s].label, point, &val);
4728: if (val >= 0) {*prob = dm->probs[s].ds; break;}
4729: }
4730: }
4731: if (!*prob) *prob = probDef;
4732: return(0);
4733: }
4735: /*@
4736: DMGetRegionDS - Get the PetscDS for a given mesh region, defined by a DMLabel
4738: Not collective
4740: Input Parameters:
4741: + dm - The DM
4742: - label - The DMLabel defining the mesh region, or NULL for the entire mesh
4744: Output Parameter:
4745: . prob - The PetscDS defined on the given region
4747: Note: If the label is missing, this function returns an error
4749: Level: advanced
4751: .seealso: DMGetRegionNumDS(), DMSetRegionDS(), DMGetDS(), DMGetCellDS()
4752: @*/
4753: PetscErrorCode DMGetRegionDS(DM dm, DMLabel label, PetscDS *ds)
4754: {
4755: PetscInt Nds = dm->Nds, s;
4761: *ds = NULL;
4762: for (s = 0; s < Nds; ++s) {
4763: if (dm->probs[s].label == label) {*ds = dm->probs[s].ds; return(0);}
4764: }
4765: return(0);
4766: }
4768: /*@
4769: DMGetRegionNumDS - Get the PetscDS for a given mesh region, defined by the region number
4771: Not collective
4773: Input Parameters:
4774: + dm - The DM
4775: - num - The region number, in [0, Nds)
4777: Output Parameters:
4778: + label - The region label
4779: - prob - The PetscDS defined on the given region
4781: Level: advanced
4783: .seealso: DMGetRegionDS(), DMSetRegionDS(), DMGetDS(), DMGetCellDS()
4784: @*/
4785: PetscErrorCode DMGetRegionNumDS(DM dm, PetscInt num, DMLabel *label, PetscDS *ds)
4786: {
4787: PetscInt Nds;
4792: DMGetNumDS(dm, &Nds);
4793: if ((num < 0) || (num >= Nds)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Region number %D is not in [0, %D)", num, Nds);
4794: if (label) {
4796: *label = dm->probs[num].label;
4797: }
4798: if (ds) {
4800: *ds = dm->probs[num].ds;
4801: }
4802: return(0);
4803: }
4805: /*@
4806: DMSetRegionDS - Set the PetscDS for a given mesh region, defined by a DMLabel
4808: Collective on DM
4810: Input Parameters:
4811: + dm - The DM
4812: . label - The DMLabel defining the mesh region, or NULL for the entire mesh
4813: - prob - The PetscDS defined on the given cell
4815: Note: If the label has a DS defined, it will be replaced. Otherwise, it will be added to the DM.
4817: Level: advanced
4819: .seealso: DMGetRegionDS(), DMGetDS(), DMGetCellDS()
4820: @*/
4821: PetscErrorCode DMSetRegionDS(DM dm, DMLabel label, PetscDS ds)
4822: {
4823: PetscInt Nds = dm->Nds, s;
4830: for (s = 0; s < Nds; ++s) {
4831: if (dm->probs[s].label == label) {
4832: dm->probs[s].ds = ds;
4833: return(0);
4834: }
4835: }
4836: DMDSEnlarge_Static(dm, Nds+1);
4837: PetscObjectReference((PetscObject) label);
4838: PetscObjectReference((PetscObject) ds);
4839: if (!label) {
4840: /* Put the NULL label at the front, so it is returned as the default */
4841: for (s = Nds-1; s >=0; --s) dm->probs[s+1] = dm->probs[s];
4842: Nds = 0;
4843: }
4844: dm->probs[Nds].label = label;
4845: dm->probs[Nds].ds = ds;
4846: return(0);
4847: }
4849: /*@
4850: DMCreateDS - Create the discrete systems for the DM based upon the fields added to the DM
4852: Collective on DM
4854: Input Parameter:
4855: . dm - The DM
4857: Note: If the label has a DS defined, it will be replaced. Otherwise, it will be added to the DM.
4859: Level: intermediate
4861: .seealso: DMSetField, DMAddField(), DMGetDS(), DMGetCellDS(), DMGetRegionDS(), DMSetRegionDS()
4862: @*/
4863: PetscErrorCode DMCreateDS(DM dm)
4864: {
4865: MPI_Comm comm;
4866: PetscDS prob, probh = NULL;
4867: PetscInt dimEmbed, f, s, field = 0, fieldh = 0;
4868: PetscBool doSetup = PETSC_TRUE;
4873: if (!dm->fields) return(0);
4874: /* Can only handle two label cases right now:
4875: 1) NULL
4876: 2) Hybrid cells
4877: */
4878: PetscObjectGetComm((PetscObject) dm, &comm);
4879: DMGetCoordinateDim(dm, &dimEmbed);
4880: /* Create default DS */
4881: DMGetRegionDS(dm, NULL, &prob);
4882: if (!prob) {
4883: PetscDSCreate(comm, &prob);
4884: DMSetRegionDS(dm, NULL, prob);
4885: PetscDSDestroy(&prob);
4886: DMGetRegionDS(dm, NULL, &prob);
4887: }
4888: PetscDSSetCoordinateDimension(prob, dimEmbed);
4889: /* Optionally create hybrid DS */
4890: for (f = 0; f < dm->Nf; ++f) {
4891: DMLabel label = dm->fields[f].label;
4892: PetscInt lStart, lEnd;
4894: if (label) {
4895: DM plex;
4896: PetscInt depth, pMax[4];
4898: DMConvert(dm, DMPLEX, &plex);
4899: DMPlexGetDepth(plex, &depth);
4900: DMPlexGetHybridBounds(plex, depth >= 0 ? &pMax[depth] : NULL, depth>1 ? &pMax[depth-1] : NULL, depth>2 ? &pMax[1] : NULL, &pMax[0]);
4901: DMDestroy(&plex);
4903: DMLabelGetBounds(label, &lStart, &lEnd);
4904: if (lStart < pMax[depth]) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only support labels over hybrid cells right now");
4905: PetscDSCreate(comm, &probh);
4906: DMSetRegionDS(dm, label, probh);
4907: PetscDSSetHybrid(probh, PETSC_TRUE);
4908: PetscDSSetCoordinateDimension(probh, dimEmbed);
4909: break;
4910: }
4911: }
4912: /* Set fields in DSes */
4913: for (f = 0; f < dm->Nf; ++f) {
4914: DMLabel label = dm->fields[f].label;
4915: PetscObject disc = dm->fields[f].disc;
4917: if (!label) {
4918: PetscDSSetDiscretization(prob, field++, disc);
4919: if (probh) {
4920: PetscFE subfe;
4922: PetscFEGetHeightSubspace((PetscFE) disc, 1, &subfe);
4923: PetscDSSetDiscretization(probh, fieldh++, (PetscObject) subfe);
4924: }
4925: } else {
4926: PetscDSSetDiscretization(probh, fieldh++, disc);
4927: }
4928: /* We allow people to have placeholder fields and construct the Section by hand */
4929: {
4930: PetscClassId id;
4932: PetscObjectGetClassId(disc, &id);
4933: if ((id != PETSCFE_CLASSID) && (id != PETSCFV_CLASSID)) doSetup = PETSC_FALSE;
4934: }
4935: }
4936: PetscDSDestroy(&probh);
4937: /* Setup DSes */
4938: if (doSetup) {
4939: for (s = 0; s < dm->Nds; ++s) {PetscDSSetUp(dm->probs[s].ds);}
4940: }
4941: return(0);
4942: }
4944: /*@
4945: DMCopyDS - Copy the discrete systems for the DM into another DM
4947: Collective on DM
4949: Input Parameter:
4950: . dm - The DM
4952: Output Parameter:
4953: . newdm - The DM
4955: Level: advanced
4957: .seealso: DMCopyFields(), DMAddField(), DMGetDS(), DMGetCellDS(), DMGetRegionDS(), DMSetRegionDS()
4958: @*/
4959: PetscErrorCode DMCopyDS(DM dm, DM newdm)
4960: {
4961: PetscInt Nds, s;
4965: if (dm == newdm) return(0);
4966: DMGetNumDS(dm, &Nds);
4967: DMClearDS(newdm);
4968: for (s = 0; s < Nds; ++s) {
4969: DMLabel label;
4970: PetscDS ds;
4972: DMGetRegionNumDS(dm, s, &label, &ds);
4973: DMSetRegionDS(newdm, label, ds);
4974: }
4975: return(0);
4976: }
4978: /*@
4979: DMCopyDisc - Copy the fields and discrete systems for the DM into another DM
4981: Collective on DM
4983: Input Parameter:
4984: . dm - The DM
4986: Output Parameter:
4987: . newdm - The DM
4989: Level: advanced
4991: .seealso: DMCopyFields(), DMCopyDS()
4992: @*/
4993: PetscErrorCode DMCopyDisc(DM dm, DM newdm)
4994: {
4998: if (dm == newdm) return(0);
4999: DMCopyFields(dm, newdm);
5000: DMCopyDS(dm, newdm);
5001: return(0);
5002: }
5004: PetscErrorCode DMRestrictHook_Coordinates(DM dm,DM dmc,void *ctx)
5005: {
5006: DM dm_coord,dmc_coord;
5008: Vec coords,ccoords;
5009: Mat inject;
5011: DMGetCoordinateDM(dm,&dm_coord);
5012: DMGetCoordinateDM(dmc,&dmc_coord);
5013: DMGetCoordinates(dm,&coords);
5014: DMGetCoordinates(dmc,&ccoords);
5015: if (coords && !ccoords) {
5016: DMCreateGlobalVector(dmc_coord,&ccoords);
5017: PetscObjectSetName((PetscObject)ccoords,"coordinates");
5018: DMCreateInjection(dmc_coord,dm_coord,&inject);
5019: MatRestrict(inject,coords,ccoords);
5020: MatDestroy(&inject);
5021: DMSetCoordinates(dmc,ccoords);
5022: VecDestroy(&ccoords);
5023: }
5024: return(0);
5025: }
5027: static PetscErrorCode DMSubDomainHook_Coordinates(DM dm,DM subdm,void *ctx)
5028: {
5029: DM dm_coord,subdm_coord;
5031: Vec coords,ccoords,clcoords;
5032: VecScatter *scat_i,*scat_g;
5034: DMGetCoordinateDM(dm,&dm_coord);
5035: DMGetCoordinateDM(subdm,&subdm_coord);
5036: DMGetCoordinates(dm,&coords);
5037: DMGetCoordinates(subdm,&ccoords);
5038: if (coords && !ccoords) {
5039: DMCreateGlobalVector(subdm_coord,&ccoords);
5040: PetscObjectSetName((PetscObject)ccoords,"coordinates");
5041: DMCreateLocalVector(subdm_coord,&clcoords);
5042: PetscObjectSetName((PetscObject)clcoords,"coordinates");
5043: DMCreateDomainDecompositionScatters(dm_coord,1,&subdm_coord,NULL,&scat_i,&scat_g);
5044: VecScatterBegin(scat_i[0],coords,ccoords,INSERT_VALUES,SCATTER_FORWARD);
5045: VecScatterBegin(scat_g[0],coords,clcoords,INSERT_VALUES,SCATTER_FORWARD);
5046: VecScatterEnd(scat_i[0],coords,ccoords,INSERT_VALUES,SCATTER_FORWARD);
5047: VecScatterEnd(scat_g[0],coords,clcoords,INSERT_VALUES,SCATTER_FORWARD);
5048: DMSetCoordinates(subdm,ccoords);
5049: DMSetCoordinatesLocal(subdm,clcoords);
5050: VecScatterDestroy(&scat_i[0]);
5051: VecScatterDestroy(&scat_g[0]);
5052: VecDestroy(&ccoords);
5053: VecDestroy(&clcoords);
5054: PetscFree(scat_i);
5055: PetscFree(scat_g);
5056: }
5057: return(0);
5058: }
5060: /*@
5061: DMGetDimension - Return the topological dimension of the DM
5063: Not collective
5065: Input Parameter:
5066: . dm - The DM
5068: Output Parameter:
5069: . dim - The topological dimension
5071: Level: beginner
5073: .seealso: DMSetDimension(), DMCreate()
5074: @*/
5075: PetscErrorCode DMGetDimension(DM dm, PetscInt *dim)
5076: {
5080: *dim = dm->dim;
5081: return(0);
5082: }
5084: /*@
5085: DMSetDimension - Set the topological dimension of the DM
5087: Collective on dm
5089: Input Parameters:
5090: + dm - The DM
5091: - dim - The topological dimension
5093: Level: beginner
5095: .seealso: DMGetDimension(), DMCreate()
5096: @*/
5097: PetscErrorCode DMSetDimension(DM dm, PetscInt dim)
5098: {
5099: PetscDS ds;
5105: dm->dim = dim;
5106: DMGetDS(dm, &ds);
5107: if (ds->dimEmbed < 0) {PetscDSSetCoordinateDimension(ds, dm->dim);}
5108: return(0);
5109: }
5111: /*@
5112: DMGetDimPoints - Get the half-open interval for all points of a given dimension
5114: Collective on DM
5116: Input Parameters:
5117: + dm - the DM
5118: - dim - the dimension
5120: Output Parameters:
5121: + pStart - The first point of the given dimension
5122: . pEnd - The first point following points of the given dimension
5124: Note:
5125: The points are vertices in the Hasse diagram encoding the topology. This is explained in
5126: http://arxiv.org/abs/0908.4427. If no points exist of this dimension in the storage scheme,
5127: then the interval is empty.
5129: Level: intermediate
5131: .keywords: point, Hasse Diagram, dimension
5132: .seealso: DMPLEX, DMPlexGetDepthStratum(), DMPlexGetHeightStratum()
5133: @*/
5134: PetscErrorCode DMGetDimPoints(DM dm, PetscInt dim, PetscInt *pStart, PetscInt *pEnd)
5135: {
5136: PetscInt d;
5141: DMGetDimension(dm, &d);
5142: if ((dim < 0) || (dim > d)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension %d 1", dim, d);
5143: (*dm->ops->getdimpoints)(dm, dim, pStart, pEnd);
5144: return(0);
5145: }
5147: /*@
5148: DMSetCoordinates - Sets into the DM a global vector that holds the coordinates
5150: Collective on DM
5152: Input Parameters:
5153: + dm - the DM
5154: - c - coordinate vector
5156: Notes:
5157: The coordinates do include those for ghost points, which are in the local vector.
5159: The vector c should be destroyed by the caller.
5161: Level: intermediate
5163: .keywords: distributed array, get, corners, nodes, local indices, coordinates
5164: .seealso: DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM()
5165: @*/
5166: PetscErrorCode DMSetCoordinates(DM dm, Vec c)
5167: {
5173: PetscObjectReference((PetscObject) c);
5174: VecDestroy(&dm->coordinates);
5175: dm->coordinates = c;
5176: VecDestroy(&dm->coordinatesLocal);
5177: DMCoarsenHookAdd(dm,DMRestrictHook_Coordinates,NULL,NULL);
5178: DMSubDomainHookAdd(dm,DMSubDomainHook_Coordinates,NULL,NULL);
5179: return(0);
5180: }
5182: /*@
5183: DMSetCoordinatesLocal - Sets into the DM a local vector that holds the coordinates
5185: Not collective
5187: Input Parameters:
5188: + dm - the DM
5189: - c - coordinate vector
5191: Notes:
5192: The coordinates of ghost points can be set using DMSetCoordinates()
5193: followed by DMGetCoordinatesLocal(). This is intended to enable the
5194: setting of ghost coordinates outside of the domain.
5196: The vector c should be destroyed by the caller.
5198: Level: intermediate
5200: .keywords: distributed array, get, corners, nodes, local indices, coordinates
5201: .seealso: DMGetCoordinatesLocal(), DMSetCoordinates(), DMGetCoordinates(), DMGetCoordinateDM()
5202: @*/
5203: PetscErrorCode DMSetCoordinatesLocal(DM dm, Vec c)
5204: {
5210: PetscObjectReference((PetscObject) c);
5211: VecDestroy(&dm->coordinatesLocal);
5213: dm->coordinatesLocal = c;
5215: VecDestroy(&dm->coordinates);
5216: return(0);
5217: }
5219: /*@
5220: DMGetCoordinates - Gets a global vector with the coordinates associated with the DM.
5222: Collective on DM
5224: Input Parameter:
5225: . dm - the DM
5227: Output Parameter:
5228: . c - global coordinate vector
5230: Note:
5231: This is a borrowed reference, so the user should NOT destroy this vector
5233: Each process has only the local coordinates (does NOT have the ghost coordinates).
5235: For DMDA, in two and three dimensions coordinates are interlaced (x_0,y_0,x_1,y_1,...)
5236: and (x_0,y_0,z_0,x_1,y_1,z_1...)
5238: Level: intermediate
5240: .keywords: distributed array, get, corners, nodes, local indices, coordinates
5241: .seealso: DMSetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM()
5242: @*/
5243: PetscErrorCode DMGetCoordinates(DM dm, Vec *c)
5244: {
5250: if (!dm->coordinates && dm->coordinatesLocal) {
5251: DM cdm = NULL;
5252: PetscBool localized;
5254: DMGetCoordinateDM(dm, &cdm);
5255: DMCreateGlobalVector(cdm, &dm->coordinates);
5256: DMGetCoordinatesLocalized(dm, &localized);
5257: /* Block size is not correctly set by CreateGlobalVector() if coordinates are localized */
5258: if (localized) {
5259: PetscInt cdim;
5261: DMGetCoordinateDim(dm, &cdim);
5262: VecSetBlockSize(dm->coordinates, cdim);
5263: }
5264: PetscObjectSetName((PetscObject) dm->coordinates, "coordinates");
5265: DMLocalToGlobalBegin(cdm, dm->coordinatesLocal, INSERT_VALUES, dm->coordinates);
5266: DMLocalToGlobalEnd(cdm, dm->coordinatesLocal, INSERT_VALUES, dm->coordinates);
5267: }
5268: *c = dm->coordinates;
5269: return(0);
5270: }
5272: /*@
5273: DMGetCoordinatesLocalSetUp - Prepares a local vector of coordinates, so that DMGetCoordinatesLocalNoncollective() can be used as non-collective afterwards.
5275: Collective on DM
5277: Input Parameter:
5278: . dm - the DM
5280: Level: advanced
5282: .keywords: distributed array, get, corners, nodes, local indices, coordinates
5283: .seealso: DMGetCoordinatesLocalNoncollective()
5284: @*/
5285: PetscErrorCode DMGetCoordinatesLocalSetUp(DM dm)
5286: {
5291: if (!dm->coordinatesLocal && dm->coordinates) {
5292: DM cdm = NULL;
5293: PetscBool localized;
5295: DMGetCoordinateDM(dm, &cdm);
5296: DMCreateLocalVector(cdm, &dm->coordinatesLocal);
5297: DMGetCoordinatesLocalized(dm, &localized);
5298: /* Block size is not correctly set by CreateLocalVector() if coordinates are localized */
5299: if (localized) {
5300: PetscInt cdim;
5302: DMGetCoordinateDim(dm, &cdim);
5303: VecSetBlockSize(dm->coordinates, cdim);
5304: }
5305: PetscObjectSetName((PetscObject) dm->coordinatesLocal, "coordinates");
5306: DMGlobalToLocalBegin(cdm, dm->coordinates, INSERT_VALUES, dm->coordinatesLocal);
5307: DMGlobalToLocalEnd(cdm, dm->coordinates, INSERT_VALUES, dm->coordinatesLocal);
5308: }
5309: return(0);
5310: }
5312: /*@
5313: DMGetCoordinatesLocal - Gets a local vector with the coordinates associated with the DM.
5315: Collective on DM
5317: Input Parameter:
5318: . dm - the DM
5320: Output Parameter:
5321: . c - coordinate vector
5323: Note:
5324: This is a borrowed reference, so the user should NOT destroy this vector
5326: Each process has the local and ghost coordinates
5328: For DMDA, in two and three dimensions coordinates are interlaced (x_0,y_0,x_1,y_1,...)
5329: and (x_0,y_0,z_0,x_1,y_1,z_1...)
5331: Level: intermediate
5333: .keywords: distributed array, get, corners, nodes, local indices, coordinates
5334: .seealso: DMSetCoordinatesLocal(), DMGetCoordinates(), DMSetCoordinates(), DMGetCoordinateDM(), DMGetCoordinatesLocalNoncollective()
5335: @*/
5336: PetscErrorCode DMGetCoordinatesLocal(DM dm, Vec *c)
5337: {
5343: DMGetCoordinatesLocalSetUp(dm);
5344: *c = dm->coordinatesLocal;
5345: return(0);
5346: }
5348: /*@
5349: DMGetCoordinatesLocalNoncollective - Non-collective version of DMGetCoordinatesLocal(). Fails if global coordinates have been set and DMGetCoordinatesLocalSetUp() not called.
5351: Not collective
5353: Input Parameter:
5354: . dm - the DM
5356: Output Parameter:
5357: . c - coordinate vector
5359: Level: advanced
5361: .keywords: distributed array, get, corners, nodes, local indices, coordinates
5362: .seealso: DMGetCoordinatesLocalSetUp(), DMGetCoordinatesLocal(), DMSetCoordinatesLocal(), DMGetCoordinates(), DMSetCoordinates(), DMGetCoordinateDM()
5363: @*/
5364: PetscErrorCode DMGetCoordinatesLocalNoncollective(DM dm, Vec *c)
5365: {
5369: if (!dm->coordinatesLocal && dm->coordinates) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DMGetCoordinatesLocalSetUp() has not been called");
5370: *c = dm->coordinatesLocal;
5371: return(0);
5372: }
5374: /*@
5375: DMGetCoordinatesLocalTuple - Gets a local vector with the coordinates of specified points and section describing its layout.
5377: Not collective
5379: Input Parameter:
5380: + dm - the DM
5381: - p - the IS of points whose coordinates will be returned
5383: Output Parameter:
5384: + pCoordSection - the PetscSection describing the layout of pCoord, i.e. each point corresponds to one point in p, and DOFs correspond to coordinates
5385: - pCoord - the Vec with coordinates of points in p
5387: Note:
5388: DMGetCoordinatesLocalSetUp() must be called first. This function employs DMGetCoordinatesLocalNoncollective() so it is not collective.
5390: This creates a new vector, so the user SHOULD destroy this vector
5392: Each process has the local and ghost coordinates
5394: For DMDA, in two and three dimensions coordinates are interlaced (x_0,y_0,x_1,y_1,...)
5395: and (x_0,y_0,z_0,x_1,y_1,z_1...)
5397: Level: advanced
5399: .keywords: distributed array, get, corners, nodes, local indices, coordinates
5400: .seealso: DMSetCoordinatesLocal(), DMGetCoordinatesLocal(), DMGetCoordinatesLocalNoncollective(), DMGetCoordinatesLocalSetUp(), DMGetCoordinates(), DMSetCoordinates(), DMGetCoordinateDM()
5401: @*/
5402: PetscErrorCode DMGetCoordinatesLocalTuple(DM dm, IS p, PetscSection *pCoordSection, Vec *pCoord)
5403: {
5404: PetscSection cs, newcs;
5405: Vec coords;
5406: const PetscScalar *arr;
5407: PetscScalar *newarr=NULL;
5408: PetscInt n;
5409: PetscErrorCode ierr;
5416: if (!dm->coordinatesLocal) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DMGetCoordinatesLocalSetUp() has not been called or coordinates not set");
5417: if (!dm->coordinateDM || !dm->coordinateDM->defaultSection) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DM not supported");
5418: cs = dm->coordinateDM->defaultSection;
5419: coords = dm->coordinatesLocal;
5420: VecGetArrayRead(coords, &arr);
5421: PetscSectionExtractDofsFromArray(cs, MPIU_SCALAR, arr, p, &newcs, pCoord ? ((void**)&newarr) : NULL);
5422: VecRestoreArrayRead(coords, &arr);
5423: if (pCoord) {
5424: PetscSectionGetStorageSize(newcs, &n);
5425: /* set array in two steps to mimic PETSC_OWN_POINTER */
5426: VecCreateSeqWithArray(PetscObjectComm((PetscObject)p), 1, n, NULL, pCoord);
5427: VecReplaceArray(*pCoord, newarr);
5428: } else {
5429: PetscFree(newarr);
5430: }
5431: if (pCoordSection) {*pCoordSection = newcs;}
5432: else {PetscSectionDestroy(&newcs);}
5433: return(0);
5434: }
5436: PetscErrorCode DMGetCoordinateField(DM dm, DMField *field)
5437: {
5443: if (!dm->coordinateField) {
5444: if (dm->ops->createcoordinatefield) {
5445: (*dm->ops->createcoordinatefield)(dm,&dm->coordinateField);
5446: }
5447: }
5448: *field = dm->coordinateField;
5449: return(0);
5450: }
5452: PetscErrorCode DMSetCoordinateField(DM dm, DMField field)
5453: {
5459: PetscObjectReference((PetscObject)field);
5460: DMFieldDestroy(&dm->coordinateField);
5461: dm->coordinateField = field;
5462: return(0);
5463: }
5465: /*@
5466: DMGetCoordinateDM - Gets the DM that prescribes coordinate layout and scatters between global and local coordinates
5468: Collective on DM
5470: Input Parameter:
5471: . dm - the DM
5473: Output Parameter:
5474: . cdm - coordinate DM
5476: Level: intermediate
5478: .keywords: distributed array, get, corners, nodes, local indices, coordinates
5479: .seealso: DMSetCoordinateDM(), DMSetCoordinates(), DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLocal()
5480: @*/
5481: PetscErrorCode DMGetCoordinateDM(DM dm, DM *cdm)
5482: {
5488: if (!dm->coordinateDM) {
5489: DM cdm;
5491: if (!dm->ops->createcoordinatedm) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unable to create coordinates for this DM");
5492: (*dm->ops->createcoordinatedm)(dm, &cdm);
5493: /* Just in case the DM sets the coordinate DM when creating it (DMP4est can do this, because it may not setup
5494: * until the call to CreateCoordinateDM) */
5495: DMDestroy(&dm->coordinateDM);
5496: dm->coordinateDM = cdm;
5497: }
5498: *cdm = dm->coordinateDM;
5499: return(0);
5500: }
5502: /*@
5503: DMSetCoordinateDM - Sets the DM that prescribes coordinate layout and scatters between global and local coordinates
5505: Logically Collective on DM
5507: Input Parameters:
5508: + dm - the DM
5509: - cdm - coordinate DM
5511: Level: intermediate
5513: .keywords: distributed array, get, corners, nodes, local indices, coordinates
5514: .seealso: DMGetCoordinateDM(), DMSetCoordinates(), DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLocal()
5515: @*/
5516: PetscErrorCode DMSetCoordinateDM(DM dm, DM cdm)
5517: {
5523: PetscObjectReference((PetscObject)cdm);
5524: DMDestroy(&dm->coordinateDM);
5525: dm->coordinateDM = cdm;
5526: return(0);
5527: }
5529: /*@
5530: DMGetCoordinateDim - Retrieve the dimension of embedding space for coordinate values.
5532: Not Collective
5534: Input Parameter:
5535: . dm - The DM object
5537: Output Parameter:
5538: . dim - The embedding dimension
5540: Level: intermediate
5542: .keywords: mesh, coordinates
5543: .seealso: DMSetCoordinateDim(), DMGetCoordinateSection(), DMGetCoordinateDM(), DMGetSection(), DMSetSection()
5544: @*/
5545: PetscErrorCode DMGetCoordinateDim(DM dm, PetscInt *dim)
5546: {
5550: if (dm->dimEmbed == PETSC_DEFAULT) {
5551: dm->dimEmbed = dm->dim;
5552: }
5553: *dim = dm->dimEmbed;
5554: return(0);
5555: }
5557: /*@
5558: DMSetCoordinateDim - Set the dimension of the embedding space for coordinate values.
5560: Not Collective
5562: Input Parameters:
5563: + dm - The DM object
5564: - dim - The embedding dimension
5566: Level: intermediate
5568: .keywords: mesh, coordinates
5569: .seealso: DMGetCoordinateDim(), DMSetCoordinateSection(), DMGetCoordinateSection(), DMGetSection(), DMSetSection()
5570: @*/
5571: PetscErrorCode DMSetCoordinateDim(DM dm, PetscInt dim)
5572: {
5573: PetscDS ds;
5578: dm->dimEmbed = dim;
5579: DMGetDS(dm, &ds);
5580: PetscDSSetCoordinateDimension(ds, dim);
5581: return(0);
5582: }
5584: /*@
5585: DMGetCoordinateSection - Retrieve the layout of coordinate values over the mesh.
5587: Collective on DM
5589: Input Parameter:
5590: . dm - The DM object
5592: Output Parameter:
5593: . section - The PetscSection object
5595: Level: intermediate
5597: .keywords: mesh, coordinates
5598: .seealso: DMGetCoordinateDM(), DMGetSection(), DMSetSection()
5599: @*/
5600: PetscErrorCode DMGetCoordinateSection(DM dm, PetscSection *section)
5601: {
5602: DM cdm;
5608: DMGetCoordinateDM(dm, &cdm);
5609: DMGetSection(cdm, section);
5610: return(0);
5611: }
5613: /*@
5614: DMSetCoordinateSection - Set the layout of coordinate values over the mesh.
5616: Not Collective
5618: Input Parameters:
5619: + dm - The DM object
5620: . dim - The embedding dimension, or PETSC_DETERMINE
5621: - section - The PetscSection object
5623: Level: intermediate
5625: .keywords: mesh, coordinates
5626: .seealso: DMGetCoordinateSection(), DMGetSection(), DMSetSection()
5627: @*/
5628: PetscErrorCode DMSetCoordinateSection(DM dm, PetscInt dim, PetscSection section)
5629: {
5630: DM cdm;
5636: DMGetCoordinateDM(dm, &cdm);
5637: DMSetSection(cdm, section);
5638: if (dim == PETSC_DETERMINE) {
5639: PetscInt d = PETSC_DEFAULT;
5640: PetscInt pStart, pEnd, vStart, vEnd, v, dd;
5642: PetscSectionGetChart(section, &pStart, &pEnd);
5643: DMGetDimPoints(dm, 0, &vStart, &vEnd);
5644: pStart = PetscMax(vStart, pStart);
5645: pEnd = PetscMin(vEnd, pEnd);
5646: for (v = pStart; v < pEnd; ++v) {
5647: PetscSectionGetDof(section, v, &dd);
5648: if (dd) {d = dd; break;}
5649: }
5650: if (d >= 0) {DMSetCoordinateDim(dm, d);}
5651: }
5652: return(0);
5653: }
5655: /*@C
5656: DMGetPeriodicity - Get the description of mesh periodicity
5658: Input Parameters:
5659: . dm - The DM object
5661: Output Parameters:
5662: + per - Whether the DM is periodic or not
5663: . maxCell - Over distances greater than this, we can assume a point has crossed over to another sheet, when trying to localize cell coordinates
5664: . L - If we assume the mesh is a torus, this is the length of each coordinate
5665: - bd - This describes the type of periodicity in each topological dimension
5667: Level: developer
5669: .seealso: DMGetPeriodicity()
5670: @*/
5671: PetscErrorCode DMGetPeriodicity(DM dm, PetscBool *per, const PetscReal **maxCell, const PetscReal **L, const DMBoundaryType **bd)
5672: {
5675: if (per) *per = dm->periodic;
5676: if (L) *L = dm->L;
5677: if (maxCell) *maxCell = dm->maxCell;
5678: if (bd) *bd = dm->bdtype;
5679: return(0);
5680: }
5682: /*@C
5683: DMSetPeriodicity - Set the description of mesh periodicity
5685: Input Parameters:
5686: + dm - The DM object
5687: . per - Whether the DM is periodic or not. If maxCell is not provided, coordinates need to be localized
5688: . maxCell - Over distances greater than this, we can assume a point has crossed over to another sheet, when trying to localize cell coordinates
5689: . L - If we assume the mesh is a torus, this is the length of each coordinate
5690: - bd - This describes the type of periodicity in each topological dimension
5692: Level: developer
5694: .seealso: DMGetPeriodicity()
5695: @*/
5696: PetscErrorCode DMSetPeriodicity(DM dm, PetscBool per, const PetscReal maxCell[], const PetscReal L[], const DMBoundaryType bd[])
5697: {
5698: PetscInt dim, d;
5704: if (maxCell) {
5708: }
5709: PetscFree3(dm->L,dm->maxCell,dm->bdtype);
5710: DMGetDimension(dm, &dim);
5711: if (maxCell) {
5712: PetscMalloc3(dim,&dm->L,dim,&dm->maxCell,dim,&dm->bdtype);
5713: for (d = 0; d < dim; ++d) {dm->L[d] = L[d]; dm->maxCell[d] = maxCell[d]; dm->bdtype[d] = bd[d];}
5714: }
5715: dm->periodic = per;
5716: return(0);
5717: }
5719: /*@
5720: 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.
5722: Input Parameters:
5723: + dm - The DM
5724: . in - The input coordinate point (dim numbers)
5725: - endpoint - Include the endpoint L_i
5727: Output Parameter:
5728: . out - The localized coordinate point
5730: Level: developer
5732: .seealso: DMLocalizeCoordinates(), DMLocalizeAddCoordinate()
5733: @*/
5734: PetscErrorCode DMLocalizeCoordinate(DM dm, const PetscScalar in[], PetscBool endpoint, PetscScalar out[])
5735: {
5736: PetscInt dim, d;
5740: DMGetCoordinateDim(dm, &dim);
5741: if (!dm->maxCell) {
5742: for (d = 0; d < dim; ++d) out[d] = in[d];
5743: } else {
5744: if (endpoint) {
5745: for (d = 0; d < dim; ++d) {
5746: 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)) {
5747: out[d] = in[d] - dm->L[d]*(PetscFloorReal(PetscRealPart(in[d])/dm->L[d]) - 1);
5748: } else {
5749: out[d] = in[d] - dm->L[d]*PetscFloorReal(PetscRealPart(in[d])/dm->L[d]);
5750: }
5751: }
5752: } else {
5753: for (d = 0; d < dim; ++d) {
5754: out[d] = in[d] - dm->L[d]*PetscFloorReal(PetscRealPart(in[d])/dm->L[d]);
5755: }
5756: }
5757: }
5758: return(0);
5759: }
5761: /*
5762: 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.
5764: Input Parameters:
5765: + dm - The DM
5766: . dim - The spatial dimension
5767: . anchor - The anchor point, the input point can be no more than maxCell away from it
5768: - in - The input coordinate point (dim numbers)
5770: Output Parameter:
5771: . out - The localized coordinate point
5773: Level: developer
5775: 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
5777: .seealso: DMLocalizeCoordinates(), DMLocalizeAddCoordinate()
5778: */
5779: PetscErrorCode DMLocalizeCoordinate_Internal(DM dm, PetscInt dim, const PetscScalar anchor[], const PetscScalar in[], PetscScalar out[])
5780: {
5781: PetscInt d;
5784: if (!dm->maxCell) {
5785: for (d = 0; d < dim; ++d) out[d] = in[d];
5786: } else {
5787: for (d = 0; d < dim; ++d) {
5788: if ((dm->bdtype[d] != DM_BOUNDARY_NONE) && (PetscAbsScalar(anchor[d] - in[d]) > dm->maxCell[d])) {
5789: out[d] = PetscRealPart(anchor[d]) > PetscRealPart(in[d]) ? dm->L[d] + in[d] : in[d] - dm->L[d];
5790: } else {
5791: out[d] = in[d];
5792: }
5793: }
5794: }
5795: return(0);
5796: }
5797: PetscErrorCode DMLocalizeCoordinateReal_Internal(DM dm, PetscInt dim, const PetscReal anchor[], const PetscReal in[], PetscReal out[])
5798: {
5799: PetscInt d;
5802: if (!dm->maxCell) {
5803: for (d = 0; d < dim; ++d) out[d] = in[d];
5804: } else {
5805: for (d = 0; d < dim; ++d) {
5806: if ((dm->bdtype[d] != DM_BOUNDARY_NONE) && (PetscAbsReal(anchor[d] - in[d]) > dm->maxCell[d])) {
5807: out[d] = anchor[d] > in[d] ? dm->L[d] + in[d] : in[d] - dm->L[d];
5808: } else {
5809: out[d] = in[d];
5810: }
5811: }
5812: }
5813: return(0);
5814: }
5816: /*
5817: 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.
5819: Input Parameters:
5820: + dm - The DM
5821: . dim - The spatial dimension
5822: . anchor - The anchor point, the input point can be no more than maxCell away from it
5823: . in - The input coordinate delta (dim numbers)
5824: - out - The input coordinate point (dim numbers)
5826: Output Parameter:
5827: . out - The localized coordinate in + out
5829: Level: developer
5831: 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
5833: .seealso: DMLocalizeCoordinates(), DMLocalizeCoordinate()
5834: */
5835: PetscErrorCode DMLocalizeAddCoordinate_Internal(DM dm, PetscInt dim, const PetscScalar anchor[], const PetscScalar in[], PetscScalar out[])
5836: {
5837: PetscInt d;
5840: if (!dm->maxCell) {
5841: for (d = 0; d < dim; ++d) out[d] += in[d];
5842: } else {
5843: for (d = 0; d < dim; ++d) {
5844: if ((dm->bdtype[d] != DM_BOUNDARY_NONE) && (PetscAbsScalar(anchor[d] - in[d]) > dm->maxCell[d])) {
5845: out[d] += PetscRealPart(anchor[d]) > PetscRealPart(in[d]) ? dm->L[d] + in[d] : in[d] - dm->L[d];
5846: } else {
5847: out[d] += in[d];
5848: }
5849: }
5850: }
5851: return(0);
5852: }
5854: /*@
5855: DMGetCoordinatesLocalizedLocal - Check if the DM coordinates have been localized for cells on this process
5857: Not collective
5859: Input Parameter:
5860: . dm - The DM
5862: Output Parameter:
5863: areLocalized - True if localized
5865: Level: developer
5867: .seealso: DMLocalizeCoordinates(), DMGetCoordinatesLocalized(), DMSetPeriodicity()
5868: @*/
5869: PetscErrorCode DMGetCoordinatesLocalizedLocal(DM dm,PetscBool *areLocalized)
5870: {
5871: DM cdm;
5872: PetscSection coordSection;
5873: PetscInt cStart, cEnd, sStart, sEnd, c, dof;
5874: PetscBool isPlex, alreadyLocalized;
5880: *areLocalized = PETSC_FALSE;
5882: /* We need some generic way of refering to cells/vertices */
5883: DMGetCoordinateDM(dm, &cdm);
5884: PetscObjectTypeCompare((PetscObject) cdm, DMPLEX, &isPlex);
5885: if (!isPlex) return(0);
5887: DMGetCoordinateSection(dm, &coordSection);
5888: DMPlexGetHeightStratum(cdm, 0, &cStart, &cEnd);
5889: PetscSectionGetChart(coordSection, &sStart, &sEnd);
5890: alreadyLocalized = PETSC_FALSE;
5891: for (c = cStart; c < cEnd; ++c) {
5892: if (c < sStart || c >= sEnd) continue;
5893: PetscSectionGetDof(coordSection, c, &dof);
5894: if (dof) { alreadyLocalized = PETSC_TRUE; break; }
5895: }
5896: *areLocalized = alreadyLocalized;
5897: return(0);
5898: }
5900: /*@
5901: DMGetCoordinatesLocalized - Check if the DM coordinates have been localized for cells
5903: Collective on dm
5905: Input Parameter:
5906: . dm - The DM
5908: Output Parameter:
5909: areLocalized - True if localized
5911: Level: developer
5913: .seealso: DMLocalizeCoordinates(), DMSetPeriodicity(), DMGetCoordinatesLocalizedLocal()
5914: @*/
5915: PetscErrorCode DMGetCoordinatesLocalized(DM dm,PetscBool *areLocalized)
5916: {
5917: PetscBool localized;
5923: DMGetCoordinatesLocalizedLocal(dm,&localized);
5924: MPIU_Allreduce(&localized,areLocalized,1,MPIU_BOOL,MPI_LOR,PetscObjectComm((PetscObject)dm));
5925: return(0);
5926: }
5928: /*@
5929: DMLocalizeCoordinates - If a mesh is periodic, create local coordinates for cells having periodic faces
5931: Collective on dm
5933: Input Parameter:
5934: . dm - The DM
5936: Level: developer
5938: .seealso: DMSetPeriodicity(), DMLocalizeCoordinate(), DMLocalizeAddCoordinate()
5939: @*/
5940: PetscErrorCode DMLocalizeCoordinates(DM dm)
5941: {
5942: DM cdm;
5943: PetscSection coordSection, cSection;
5944: Vec coordinates, cVec;
5945: PetscScalar *coords, *coords2, *anchor, *localized;
5946: PetscInt Nc, vStart, vEnd, v, sStart, sEnd, newStart = PETSC_MAX_INT, newEnd = PETSC_MIN_INT, dof, d, off, off2, bs, coordSize;
5947: PetscBool alreadyLocalized, alreadyLocalizedGlobal;
5948: PetscInt maxHeight = 0, h;
5949: PetscInt *pStart = NULL, *pEnd = NULL;
5954: if (!dm->periodic) return(0);
5955: DMGetCoordinatesLocalized(dm, &alreadyLocalized);
5956: if (alreadyLocalized) return(0);
5958: /* We need some generic way of refering to cells/vertices */
5959: DMGetCoordinateDM(dm, &cdm);
5960: {
5961: PetscBool isplex;
5963: PetscObjectTypeCompare((PetscObject) cdm, DMPLEX, &isplex);
5964: if (isplex) {
5965: DMPlexGetDepthStratum(cdm, 0, &vStart, &vEnd);
5966: DMPlexGetMaxProjectionHeight(cdm,&maxHeight);
5967: DMGetWorkArray(dm,2*(maxHeight + 1),MPIU_INT,&pStart);
5968: pEnd = &pStart[maxHeight + 1];
5969: newStart = vStart;
5970: newEnd = vEnd;
5971: for (h = 0; h <= maxHeight; h++) {
5972: DMPlexGetHeightStratum(cdm, h, &pStart[h], &pEnd[h]);
5973: newStart = PetscMin(newStart,pStart[h]);
5974: newEnd = PetscMax(newEnd,pEnd[h]);
5975: }
5976: } else SETERRQ(PetscObjectComm((PetscObject) cdm), PETSC_ERR_ARG_WRONG, "Coordinate localization requires a DMPLEX coordinate DM");
5977: }
5978: DMGetCoordinatesLocal(dm, &coordinates);
5979: if (!coordinates) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Missing local coordinates vector");
5980: DMGetCoordinateSection(dm, &coordSection);
5981: VecGetBlockSize(coordinates, &bs);
5982: PetscSectionGetChart(coordSection,&sStart,&sEnd);
5984: PetscSectionCreate(PetscObjectComm((PetscObject) dm), &cSection);
5985: PetscSectionSetNumFields(cSection, 1);
5986: PetscSectionGetFieldComponents(coordSection, 0, &Nc);
5987: PetscSectionSetFieldComponents(cSection, 0, Nc);
5988: PetscSectionSetChart(cSection, newStart, newEnd);
5990: DMGetWorkArray(dm, 2 * bs, MPIU_SCALAR, &anchor);
5991: localized = &anchor[bs];
5992: alreadyLocalized = alreadyLocalizedGlobal = PETSC_TRUE;
5993: for (h = 0; h <= maxHeight; h++) {
5994: PetscInt cStart = pStart[h], cEnd = pEnd[h], c;
5996: for (c = cStart; c < cEnd; ++c) {
5997: PetscScalar *cellCoords = NULL;
5998: PetscInt b;
6000: if (c < sStart || c >= sEnd) alreadyLocalized = PETSC_FALSE;
6001: DMPlexVecGetClosure(cdm, coordSection, coordinates, c, &dof, &cellCoords);
6002: for (b = 0; b < bs; ++b) anchor[b] = cellCoords[b];
6003: for (d = 0; d < dof/bs; ++d) {
6004: DMLocalizeCoordinate_Internal(dm, bs, anchor, &cellCoords[d*bs], localized);
6005: for (b = 0; b < bs; b++) {
6006: if (cellCoords[d*bs + b] != localized[b]) break;
6007: }
6008: if (b < bs) break;
6009: }
6010: if (d < dof/bs) {
6011: if (c >= sStart && c < sEnd) {
6012: PetscInt cdof;
6014: PetscSectionGetDof(coordSection, c, &cdof);
6015: if (cdof != dof) alreadyLocalized = PETSC_FALSE;
6016: }
6017: PetscSectionSetDof(cSection, c, dof);
6018: PetscSectionSetFieldDof(cSection, c, 0, dof);
6019: }
6020: DMPlexVecRestoreClosure(cdm, coordSection, coordinates, c, &dof, &cellCoords);
6021: }
6022: }
6023: MPI_Allreduce(&alreadyLocalized,&alreadyLocalizedGlobal,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)dm));
6024: if (alreadyLocalizedGlobal) {
6025: DMRestoreWorkArray(dm, 2 * bs, MPIU_SCALAR, &anchor);
6026: PetscSectionDestroy(&cSection);
6027: DMRestoreWorkArray(dm,2*(maxHeight + 1),MPIU_INT,&pStart);
6028: return(0);
6029: }
6030: for (v = vStart; v < vEnd; ++v) {
6031: PetscSectionGetDof(coordSection, v, &dof);
6032: PetscSectionSetDof(cSection, v, dof);
6033: PetscSectionSetFieldDof(cSection, v, 0, dof);
6034: }
6035: PetscSectionSetUp(cSection);
6036: PetscSectionGetStorageSize(cSection, &coordSize);
6037: VecCreate(PETSC_COMM_SELF, &cVec);
6038: PetscObjectSetName((PetscObject)cVec,"coordinates");
6039: VecSetBlockSize(cVec, bs);
6040: VecSetSizes(cVec, coordSize, PETSC_DETERMINE);
6041: VecSetType(cVec, VECSTANDARD);
6042: VecGetArrayRead(coordinates, (const PetscScalar**)&coords);
6043: VecGetArray(cVec, &coords2);
6044: for (v = vStart; v < vEnd; ++v) {
6045: PetscSectionGetDof(coordSection, v, &dof);
6046: PetscSectionGetOffset(coordSection, v, &off);
6047: PetscSectionGetOffset(cSection, v, &off2);
6048: for (d = 0; d < dof; ++d) coords2[off2+d] = coords[off+d];
6049: }
6050: for (h = 0; h <= maxHeight; h++) {
6051: PetscInt cStart = pStart[h], cEnd = pEnd[h], c;
6053: for (c = cStart; c < cEnd; ++c) {
6054: PetscScalar *cellCoords = NULL;
6055: PetscInt b, cdof;
6057: PetscSectionGetDof(cSection,c,&cdof);
6058: if (!cdof) continue;
6059: DMPlexVecGetClosure(cdm, coordSection, coordinates, c, &dof, &cellCoords);
6060: PetscSectionGetOffset(cSection, c, &off2);
6061: for (b = 0; b < bs; ++b) anchor[b] = cellCoords[b];
6062: for (d = 0; d < dof/bs; ++d) {DMLocalizeCoordinate_Internal(dm, bs, anchor, &cellCoords[d*bs], &coords2[off2+d*bs]);}
6063: DMPlexVecRestoreClosure(cdm, coordSection, coordinates, c, &dof, &cellCoords);
6064: }
6065: }
6066: DMRestoreWorkArray(dm, 2 * bs, MPIU_SCALAR, &anchor);
6067: DMRestoreWorkArray(dm,2*(maxHeight + 1),MPIU_INT,&pStart);
6068: VecRestoreArrayRead(coordinates, (const PetscScalar**)&coords);
6069: VecRestoreArray(cVec, &coords2);
6070: DMSetCoordinateSection(dm, PETSC_DETERMINE, cSection);
6071: DMSetCoordinatesLocal(dm, cVec);
6072: VecDestroy(&cVec);
6073: PetscSectionDestroy(&cSection);
6074: return(0);
6075: }
6077: /*@
6078: DMLocatePoints - Locate the points in v in the mesh and return a PetscSF of the containing cells
6080: Collective on Vec v (see explanation below)
6082: Input Parameters:
6083: + dm - The DM
6084: . v - The Vec of points
6085: . ltype - The type of point location, e.g. DM_POINTLOCATION_NONE or DM_POINTLOCATION_NEAREST
6086: - cells - Points to either NULL, or a PetscSF with guesses for which cells contain each point.
6088: Output Parameter:
6089: + v - The Vec of points, which now contains the nearest mesh points to the given points if DM_POINTLOCATION_NEAREST is used
6090: - cells - The PetscSF containing the ranks and local indices of the containing points.
6093: Level: developer
6095: Notes:
6096: To do a search of the local cells of the mesh, v should have PETSC_COMM_SELF as its communicator.
6097: To do a search of all the cells in the distributed mesh, v should have the same communicator as dm.
6099: If *cellSF is NULL on input, a PetscSF will be created.
6100: If *cellSF is not NULL on input, it should point to an existing PetscSF, whose graph will be used as initial guesses.
6102: An array that maps each point to its containing cell can be obtained with
6104: $ const PetscSFNode *cells;
6105: $ PetscInt nFound;
6106: $ const PetscInt *found;
6107: $
6108: $ PetscSFGetGraph(cellSF,NULL,&nFound,&found,&cells);
6110: Where cells[i].rank is the rank of the cell containing point found[i] (or i if found == NULL), and cells[i].index is
6111: the index of the cell in its rank's local numbering.
6113: .keywords: point location, mesh
6114: .seealso: DMSetCoordinates(), DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLocal(), DMPointLocationType
6115: @*/
6116: PetscErrorCode DMLocatePoints(DM dm, Vec v, DMPointLocationType ltype, PetscSF *cellSF)
6117: {
6124: if (*cellSF) {
6125: PetscMPIInt result;
6128: MPI_Comm_compare(PetscObjectComm((PetscObject)v),PetscObjectComm((PetscObject)*cellSF),&result);
6129: if (result != MPI_IDENT && result != MPI_CONGRUENT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"cellSF must have a communicator congruent to v's");
6130: } else {
6131: PetscSFCreate(PetscObjectComm((PetscObject)v),cellSF);
6132: }
6133: PetscLogEventBegin(DM_LocatePoints,dm,0,0,0);
6134: if (dm->ops->locatepoints) {
6135: (*dm->ops->locatepoints)(dm,v,ltype,*cellSF);
6136: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Point location not available for this DM");
6137: PetscLogEventEnd(DM_LocatePoints,dm,0,0,0);
6138: return(0);
6139: }
6141: /*@
6142: DMGetOutputDM - Retrieve the DM associated with the layout for output
6144: Collective on dm
6146: Input Parameter:
6147: . dm - The original DM
6149: Output Parameter:
6150: . odm - The DM which provides the layout for output
6152: Level: intermediate
6154: .seealso: VecView(), DMGetGlobalSection()
6155: @*/
6156: PetscErrorCode DMGetOutputDM(DM dm, DM *odm)
6157: {
6158: PetscSection section;
6159: PetscBool hasConstraints, ghasConstraints;
6165: DMGetSection(dm, §ion);
6166: PetscSectionHasConstraints(section, &hasConstraints);
6167: MPI_Allreduce(&hasConstraints, &ghasConstraints, 1, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject) dm));
6168: if (!ghasConstraints) {
6169: *odm = dm;
6170: return(0);
6171: }
6172: if (!dm->dmBC) {
6173: PetscSection newSection, gsection;
6174: PetscSF sf;
6176: DMClone(dm, &dm->dmBC);
6177: DMCopyDisc(dm, dm->dmBC);
6178: PetscSectionClone(section, &newSection);
6179: DMSetSection(dm->dmBC, newSection);
6180: PetscSectionDestroy(&newSection);
6181: DMGetPointSF(dm->dmBC, &sf);
6182: PetscSectionCreateGlobalSection(section, sf, PETSC_TRUE, PETSC_FALSE, &gsection);
6183: DMSetGlobalSection(dm->dmBC, gsection);
6184: PetscSectionDestroy(&gsection);
6185: }
6186: *odm = dm->dmBC;
6187: return(0);
6188: }
6190: /*@
6191: DMGetOutputSequenceNumber - Retrieve the sequence number/value for output
6193: Input Parameter:
6194: . dm - The original DM
6196: Output Parameters:
6197: + num - The output sequence number
6198: - val - The output sequence value
6200: Level: intermediate
6202: Note: This is intended for output that should appear in sequence, for instance
6203: a set of timesteps in an HDF5 file, or a set of realizations of a stochastic system.
6205: .seealso: VecView()
6206: @*/
6207: PetscErrorCode DMGetOutputSequenceNumber(DM dm, PetscInt *num, PetscReal *val)
6208: {
6213: return(0);
6214: }
6216: /*@
6217: DMSetOutputSequenceNumber - Set the sequence number/value for output
6219: Input Parameters:
6220: + dm - The original DM
6221: . num - The output sequence number
6222: - val - The output sequence value
6224: Level: intermediate
6226: Note: This is intended for output that should appear in sequence, for instance
6227: a set of timesteps in an HDF5 file, or a set of realizations of a stochastic system.
6229: .seealso: VecView()
6230: @*/
6231: PetscErrorCode DMSetOutputSequenceNumber(DM dm, PetscInt num, PetscReal val)
6232: {
6235: dm->outputSequenceNum = num;
6236: dm->outputSequenceVal = val;
6237: return(0);
6238: }
6240: /*@C
6241: DMOutputSequenceLoad - Retrieve the sequence value from a Viewer
6243: Input Parameters:
6244: + dm - The original DM
6245: . name - The sequence name
6246: - num - The output sequence number
6248: Output Parameter:
6249: . val - The output sequence value
6251: Level: intermediate
6253: Note: This is intended for output that should appear in sequence, for instance
6254: a set of timesteps in an HDF5 file, or a set of realizations of a stochastic system.
6256: .seealso: DMGetOutputSequenceNumber(), DMSetOutputSequenceNumber(), VecView()
6257: @*/
6258: PetscErrorCode DMOutputSequenceLoad(DM dm, PetscViewer viewer, const char *name, PetscInt num, PetscReal *val)
6259: {
6260: PetscBool ishdf5;
6267: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERHDF5, &ishdf5);
6268: if (ishdf5) {
6269: #if defined(PETSC_HAVE_HDF5)
6270: PetscScalar value;
6272: DMSequenceLoad_HDF5_Internal(dm, name, num, &value, viewer);
6273: *val = PetscRealPart(value);
6274: #endif
6275: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid viewer; open viewer with PetscViewerHDF5Open()");
6276: return(0);
6277: }
6279: /*@
6280: DMGetUseNatural - Get the flag for creating a mapping to the natural order on distribution
6282: Not collective
6284: Input Parameter:
6285: . dm - The DM
6287: Output Parameter:
6288: . useNatural - The flag to build the mapping to a natural order during distribution
6290: Level: beginner
6292: .seealso: DMSetUseNatural(), DMCreate()
6293: @*/
6294: PetscErrorCode DMGetUseNatural(DM dm, PetscBool *useNatural)
6295: {
6299: *useNatural = dm->useNatural;
6300: return(0);
6301: }
6303: /*@
6304: DMSetUseNatural - Set the flag for creating a mapping to the natural order after distribution
6306: Collective on dm
6308: Input Parameters:
6309: + dm - The DM
6310: - useNatural - The flag to build the mapping to a natural order during distribution
6312: Note: This also causes the map to be build after DMCreateSubDM() and DMCreateSuperDM()
6314: Level: beginner
6316: .seealso: DMGetUseNatural(), DMCreate(), DMPlexDistribute(), DMCreateSubDM(), DMCreateSuperDM()
6317: @*/
6318: PetscErrorCode DMSetUseNatural(DM dm, PetscBool useNatural)
6319: {
6323: dm->useNatural = useNatural;
6324: return(0);
6325: }
6328: /*@C
6329: DMCreateLabel - Create a label of the given name if it does not already exist
6331: Not Collective
6333: Input Parameters:
6334: + dm - The DM object
6335: - name - The label name
6337: Level: intermediate
6339: .keywords: mesh
6340: .seealso: DMLabelCreate(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
6341: @*/
6342: PetscErrorCode DMCreateLabel(DM dm, const char name[])
6343: {
6344: DMLabelLink next = dm->labels->next;
6345: PetscBool flg = PETSC_FALSE;
6346: const char *lname;
6352: while (next) {
6353: PetscObjectGetName((PetscObject) next->label, &lname);
6354: PetscStrcmp(name, lname, &flg);
6355: if (flg) break;
6356: next = next->next;
6357: }
6358: if (!flg) {
6359: DMLabelLink tmpLabel;
6361: PetscCalloc1(1, &tmpLabel);
6362: DMLabelCreate(PETSC_COMM_SELF, name, &tmpLabel->label);
6363: tmpLabel->output = PETSC_TRUE;
6364: tmpLabel->next = dm->labels->next;
6365: dm->labels->next = tmpLabel;
6366: }
6367: return(0);
6368: }
6370: /*@C
6371: DMGetLabelValue - Get the value in a Sieve Label for the given point, with 0 as the default
6373: Not Collective
6375: Input Parameters:
6376: + dm - The DM object
6377: . name - The label name
6378: - point - The mesh point
6380: Output Parameter:
6381: . value - The label value for this point, or -1 if the point is not in the label
6383: Level: beginner
6385: .keywords: mesh
6386: .seealso: DMLabelGetValue(), DMSetLabelValue(), DMGetStratumIS()
6387: @*/
6388: PetscErrorCode DMGetLabelValue(DM dm, const char name[], PetscInt point, PetscInt *value)
6389: {
6390: DMLabel label;
6396: DMGetLabel(dm, name, &label);
6397: if (!label) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No label named %s was found", name);
6398: DMLabelGetValue(label, point, value);
6399: return(0);
6400: }
6402: /*@C
6403: DMSetLabelValue - Add a point to a Sieve Label with given value
6405: Not Collective
6407: Input Parameters:
6408: + dm - The DM object
6409: . name - The label name
6410: . point - The mesh point
6411: - value - The label value for this point
6413: Output Parameter:
6415: Level: beginner
6417: .keywords: mesh
6418: .seealso: DMLabelSetValue(), DMGetStratumIS(), DMClearLabelValue()
6419: @*/
6420: PetscErrorCode DMSetLabelValue(DM dm, const char name[], PetscInt point, PetscInt value)
6421: {
6422: DMLabel label;
6428: DMGetLabel(dm, name, &label);
6429: if (!label) {
6430: DMCreateLabel(dm, name);
6431: DMGetLabel(dm, name, &label);
6432: }
6433: DMLabelSetValue(label, point, value);
6434: return(0);
6435: }
6437: /*@C
6438: DMClearLabelValue - Remove a point from a Sieve Label with given value
6440: Not Collective
6442: Input Parameters:
6443: + dm - The DM object
6444: . name - The label name
6445: . point - The mesh point
6446: - value - The label value for this point
6448: Output Parameter:
6450: Level: beginner
6452: .keywords: mesh
6453: .seealso: DMLabelClearValue(), DMSetLabelValue(), DMGetStratumIS()
6454: @*/
6455: PetscErrorCode DMClearLabelValue(DM dm, const char name[], PetscInt point, PetscInt value)
6456: {
6457: DMLabel label;
6463: DMGetLabel(dm, name, &label);
6464: if (!label) return(0);
6465: DMLabelClearValue(label, point, value);
6466: return(0);
6467: }
6469: /*@C
6470: DMGetLabelSize - Get the number of different integer ids in a Label
6472: Not Collective
6474: Input Parameters:
6475: + dm - The DM object
6476: - name - The label name
6478: Output Parameter:
6479: . size - The number of different integer ids, or 0 if the label does not exist
6481: Level: beginner
6483: .keywords: mesh
6484: .seealso: DMLabelGetNumValues(), DMSetLabelValue()
6485: @*/
6486: PetscErrorCode DMGetLabelSize(DM dm, const char name[], PetscInt *size)
6487: {
6488: DMLabel label;
6495: DMGetLabel(dm, name, &label);
6496: *size = 0;
6497: if (!label) return(0);
6498: DMLabelGetNumValues(label, size);
6499: return(0);
6500: }
6502: /*@C
6503: DMGetLabelIdIS - Get the integer ids in a label
6505: Not Collective
6507: Input Parameters:
6508: + mesh - The DM object
6509: - name - The label name
6511: Output Parameter:
6512: . ids - The integer ids, or NULL if the label does not exist
6514: Level: beginner
6516: .keywords: mesh
6517: .seealso: DMLabelGetValueIS(), DMGetLabelSize()
6518: @*/
6519: PetscErrorCode DMGetLabelIdIS(DM dm, const char name[], IS *ids)
6520: {
6521: DMLabel label;
6528: DMGetLabel(dm, name, &label);
6529: *ids = NULL;
6530: if (label) {
6531: DMLabelGetValueIS(label, ids);
6532: } else {
6533: /* returning an empty IS */
6534: ISCreateGeneral(PETSC_COMM_SELF,0,NULL,PETSC_USE_POINTER,ids);
6535: }
6536: return(0);
6537: }
6539: /*@C
6540: DMGetStratumSize - Get the number of points in a label stratum
6542: Not Collective
6544: Input Parameters:
6545: + dm - The DM object
6546: . name - The label name
6547: - value - The stratum value
6549: Output Parameter:
6550: . size - The stratum size
6552: Level: beginner
6554: .keywords: mesh
6555: .seealso: DMLabelGetStratumSize(), DMGetLabelSize(), DMGetLabelIds()
6556: @*/
6557: PetscErrorCode DMGetStratumSize(DM dm, const char name[], PetscInt value, PetscInt *size)
6558: {
6559: DMLabel label;
6566: DMGetLabel(dm, name, &label);
6567: *size = 0;
6568: if (!label) return(0);
6569: DMLabelGetStratumSize(label, value, size);
6570: return(0);
6571: }
6573: /*@C
6574: DMGetStratumIS - Get the points in a label stratum
6576: Not Collective
6578: Input Parameters:
6579: + dm - The DM object
6580: . name - The label name
6581: - value - The stratum value
6583: Output Parameter:
6584: . points - The stratum points, or NULL if the label does not exist or does not have that value
6586: Level: beginner
6588: .keywords: mesh
6589: .seealso: DMLabelGetStratumIS(), DMGetStratumSize()
6590: @*/
6591: PetscErrorCode DMGetStratumIS(DM dm, const char name[], PetscInt value, IS *points)
6592: {
6593: DMLabel label;
6600: DMGetLabel(dm, name, &label);
6601: *points = NULL;
6602: if (!label) return(0);
6603: DMLabelGetStratumIS(label, value, points);
6604: return(0);
6605: }
6607: /*@C
6608: DMSetStratumIS - Set the points in a label stratum
6610: Not Collective
6612: Input Parameters:
6613: + dm - The DM object
6614: . name - The label name
6615: . value - The stratum value
6616: - points - The stratum points
6618: Level: beginner
6620: .keywords: mesh
6621: .seealso: DMLabelSetStratumIS(), DMGetStratumSize()
6622: @*/
6623: PetscErrorCode DMSetStratumIS(DM dm, const char name[], PetscInt value, IS points)
6624: {
6625: DMLabel label;
6632: DMGetLabel(dm, name, &label);
6633: if (!label) return(0);
6634: DMLabelSetStratumIS(label, value, points);
6635: return(0);
6636: }
6638: /*@C
6639: DMClearLabelStratum - Remove all points from a stratum from a Sieve Label
6641: Not Collective
6643: Input Parameters:
6644: + dm - The DM object
6645: . name - The label name
6646: - value - The label value for this point
6648: Output Parameter:
6650: Level: beginner
6652: .keywords: mesh
6653: .seealso: DMLabelClearStratum(), DMSetLabelValue(), DMGetStratumIS(), DMClearLabelValue()
6654: @*/
6655: PetscErrorCode DMClearLabelStratum(DM dm, const char name[], PetscInt value)
6656: {
6657: DMLabel label;
6663: DMGetLabel(dm, name, &label);
6664: if (!label) return(0);
6665: DMLabelClearStratum(label, value);
6666: return(0);
6667: }
6669: /*@
6670: DMGetNumLabels - Return the number of labels defined by the mesh
6672: Not Collective
6674: Input Parameter:
6675: . dm - The DM object
6677: Output Parameter:
6678: . numLabels - the number of Labels
6680: Level: intermediate
6682: .keywords: mesh
6683: .seealso: DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
6684: @*/
6685: PetscErrorCode DMGetNumLabels(DM dm, PetscInt *numLabels)
6686: {
6687: DMLabelLink next = dm->labels->next;
6688: PetscInt n = 0;
6693: while (next) {++n; next = next->next;}
6694: *numLabels = n;
6695: return(0);
6696: }
6698: /*@C
6699: DMGetLabelName - Return the name of nth label
6701: Not Collective
6703: Input Parameters:
6704: + dm - The DM object
6705: - n - the label number
6707: Output Parameter:
6708: . name - the label name
6710: Level: intermediate
6712: .keywords: mesh
6713: .seealso: DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
6714: @*/
6715: PetscErrorCode DMGetLabelName(DM dm, PetscInt n, const char **name)
6716: {
6717: DMLabelLink next = dm->labels->next;
6718: PetscInt l = 0;
6724: while (next) {
6725: if (l == n) {
6726: PetscObjectGetName((PetscObject) next->label, name);
6727: return(0);
6728: }
6729: ++l;
6730: next = next->next;
6731: }
6732: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label %D does not exist in this DM", n);
6733: }
6735: /*@C
6736: DMHasLabel - Determine whether the mesh has a label of a given name
6738: Not Collective
6740: Input Parameters:
6741: + dm - The DM object
6742: - name - The label name
6744: Output Parameter:
6745: . hasLabel - PETSC_TRUE if the label is present
6747: Level: intermediate
6749: .keywords: mesh
6750: .seealso: DMCreateLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
6751: @*/
6752: PetscErrorCode DMHasLabel(DM dm, const char name[], PetscBool *hasLabel)
6753: {
6754: DMLabelLink next = dm->labels->next;
6755: const char *lname;
6762: *hasLabel = PETSC_FALSE;
6763: while (next) {
6764: PetscObjectGetName((PetscObject) next->label, &lname);
6765: PetscStrcmp(name, lname, hasLabel);
6766: if (*hasLabel) break;
6767: next = next->next;
6768: }
6769: return(0);
6770: }
6772: /*@C
6773: DMGetLabel - Return the label of a given name, or NULL
6775: Not Collective
6777: Input Parameters:
6778: + dm - The DM object
6779: - name - The label name
6781: Output Parameter:
6782: . label - The DMLabel, or NULL if the label is absent
6784: Level: intermediate
6786: .keywords: mesh
6787: .seealso: DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
6788: @*/
6789: PetscErrorCode DMGetLabel(DM dm, const char name[], DMLabel *label)
6790: {
6791: DMLabelLink next = dm->labels->next;
6792: PetscBool hasLabel;
6793: const char *lname;
6800: *label = NULL;
6801: while (next) {
6802: PetscObjectGetName((PetscObject) next->label, &lname);
6803: PetscStrcmp(name, lname, &hasLabel);
6804: if (hasLabel) {
6805: *label = next->label;
6806: break;
6807: }
6808: next = next->next;
6809: }
6810: return(0);
6811: }
6813: /*@C
6814: DMGetLabelByNum - Return the nth label
6816: Not Collective
6818: Input Parameters:
6819: + dm - The DM object
6820: - n - the label number
6822: Output Parameter:
6823: . label - the label
6825: Level: intermediate
6827: .keywords: mesh
6828: .seealso: DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
6829: @*/
6830: PetscErrorCode DMGetLabelByNum(DM dm, PetscInt n, DMLabel *label)
6831: {
6832: DMLabelLink next = dm->labels->next;
6833: PetscInt l = 0;
6838: while (next) {
6839: if (l == n) {
6840: *label = next->label;
6841: return(0);
6842: }
6843: ++l;
6844: next = next->next;
6845: }
6846: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label %D does not exist in this DM", n);
6847: }
6849: /*@C
6850: DMAddLabel - Add the label to this mesh
6852: Not Collective
6854: Input Parameters:
6855: + dm - The DM object
6856: - label - The DMLabel
6858: Level: developer
6860: .keywords: mesh
6861: .seealso: DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
6862: @*/
6863: PetscErrorCode DMAddLabel(DM dm, DMLabel label)
6864: {
6865: DMLabelLink tmpLabel;
6866: PetscBool hasLabel;
6867: const char *lname;
6872: PetscObjectGetName((PetscObject) label, &lname);
6873: DMHasLabel(dm, lname, &hasLabel);
6874: if (hasLabel) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label %s already exists in this DM", lname);
6875: PetscCalloc1(1, &tmpLabel);
6876: tmpLabel->label = label;
6877: tmpLabel->output = PETSC_TRUE;
6878: tmpLabel->next = dm->labels->next;
6879: dm->labels->next = tmpLabel;
6880: return(0);
6881: }
6883: /*@C
6884: DMRemoveLabel - Remove the label from this mesh
6886: Not Collective
6888: Input Parameters:
6889: + dm - The DM object
6890: - name - The label name
6892: Output Parameter:
6893: . label - The DMLabel, or NULL if the label is absent
6895: Level: developer
6897: .keywords: mesh
6898: .seealso: DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
6899: @*/
6900: PetscErrorCode DMRemoveLabel(DM dm, const char name[], DMLabel *label)
6901: {
6902: DMLabelLink next = dm->labels->next;
6903: DMLabelLink last = NULL;
6904: PetscBool hasLabel;
6905: const char *lname;
6910: DMHasLabel(dm, name, &hasLabel);
6911: *label = NULL;
6912: if (!hasLabel) return(0);
6913: while (next) {
6914: PetscObjectGetName((PetscObject) next->label, &lname);
6915: PetscStrcmp(name, lname, &hasLabel);
6916: if (hasLabel) {
6917: if (last) last->next = next->next;
6918: else dm->labels->next = next->next;
6919: next->next = NULL;
6920: *label = next->label;
6921: PetscStrcmp(name, "depth", &hasLabel);
6922: if (hasLabel) {
6923: dm->depthLabel = NULL;
6924: }
6925: PetscFree(next);
6926: break;
6927: }
6928: last = next;
6929: next = next->next;
6930: }
6931: return(0);
6932: }
6934: /*@C
6935: DMGetLabelOutput - Get the output flag for a given label
6937: Not Collective
6939: Input Parameters:
6940: + dm - The DM object
6941: - name - The label name
6943: Output Parameter:
6944: . output - The flag for output
6946: Level: developer
6948: .keywords: mesh
6949: .seealso: DMSetLabelOutput(), DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
6950: @*/
6951: PetscErrorCode DMGetLabelOutput(DM dm, const char name[], PetscBool *output)
6952: {
6953: DMLabelLink next = dm->labels->next;
6954: const char *lname;
6961: while (next) {
6962: PetscBool flg;
6964: PetscObjectGetName((PetscObject) next->label, &lname);
6965: PetscStrcmp(name, lname, &flg);
6966: if (flg) {*output = next->output; return(0);}
6967: next = next->next;
6968: }
6969: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No label named %s was present in this dm", name);
6970: }
6972: /*@C
6973: DMSetLabelOutput - Set the output flag for a given label
6975: Not Collective
6977: Input Parameters:
6978: + dm - The DM object
6979: . name - The label name
6980: - output - The flag for output
6982: Level: developer
6984: .keywords: mesh
6985: .seealso: DMGetLabelOutput(), DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
6986: @*/
6987: PetscErrorCode DMSetLabelOutput(DM dm, const char name[], PetscBool output)
6988: {
6989: DMLabelLink next = dm->labels->next;
6990: const char *lname;
6996: while (next) {
6997: PetscBool flg;
6999: PetscObjectGetName((PetscObject) next->label, &lname);
7000: PetscStrcmp(name, lname, &flg);
7001: if (flg) {next->output = output; return(0);}
7002: next = next->next;
7003: }
7004: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No label named %s was present in this dm", name);
7005: }
7008: /*@
7009: DMCopyLabels - Copy labels from one mesh to another with a superset of the points
7011: Collective on DM
7013: Input Parameter:
7014: . dmA - The DM object with initial labels
7016: Output Parameter:
7017: . dmB - The DM object with copied labels
7019: Level: intermediate
7021: Note: This is typically used when interpolating or otherwise adding to a mesh
7023: .keywords: mesh
7024: .seealso: DMGetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM(), DMGetCoordinateSection()
7025: @*/
7026: PetscErrorCode DMCopyLabels(DM dmA, DM dmB)
7027: {
7028: PetscInt numLabels, l;
7032: if (dmA == dmB) return(0);
7033: DMGetNumLabels(dmA, &numLabels);
7034: for (l = 0; l < numLabels; ++l) {
7035: DMLabel label, labelNew;
7036: const char *name;
7037: PetscBool flg;
7039: DMGetLabelName(dmA, l, &name);
7040: PetscStrcmp(name, "depth", &flg);
7041: if (flg) continue;
7042: PetscStrcmp(name, "dim", &flg);
7043: if (flg) continue;
7044: DMGetLabel(dmA, name, &label);
7045: DMLabelDuplicate(label, &labelNew);
7046: DMAddLabel(dmB, labelNew);
7047: }
7048: return(0);
7049: }
7051: /*@
7052: DMGetCoarseDM - Get the coarse mesh from which this was obtained by refinement
7054: Input Parameter:
7055: . dm - The DM object
7057: Output Parameter:
7058: . cdm - The coarse DM
7060: Level: intermediate
7062: .seealso: DMSetCoarseDM()
7063: @*/
7064: PetscErrorCode DMGetCoarseDM(DM dm, DM *cdm)
7065: {
7069: *cdm = dm->coarseMesh;
7070: return(0);
7071: }
7073: /*@
7074: DMSetCoarseDM - Set the coarse mesh from which this was obtained by refinement
7076: Input Parameters:
7077: + dm - The DM object
7078: - cdm - The coarse DM
7080: Level: intermediate
7082: .seealso: DMGetCoarseDM()
7083: @*/
7084: PetscErrorCode DMSetCoarseDM(DM dm, DM cdm)
7085: {
7091: PetscObjectReference((PetscObject)cdm);
7092: DMDestroy(&dm->coarseMesh);
7093: dm->coarseMesh = cdm;
7094: return(0);
7095: }
7097: /*@
7098: DMGetFineDM - Get the fine mesh from which this was obtained by refinement
7100: Input Parameter:
7101: . dm - The DM object
7103: Output Parameter:
7104: . fdm - The fine DM
7106: Level: intermediate
7108: .seealso: DMSetFineDM()
7109: @*/
7110: PetscErrorCode DMGetFineDM(DM dm, DM *fdm)
7111: {
7115: *fdm = dm->fineMesh;
7116: return(0);
7117: }
7119: /*@
7120: DMSetFineDM - Set the fine mesh from which this was obtained by refinement
7122: Input Parameters:
7123: + dm - The DM object
7124: - fdm - The fine DM
7126: Level: intermediate
7128: .seealso: DMGetFineDM()
7129: @*/
7130: PetscErrorCode DMSetFineDM(DM dm, DM fdm)
7131: {
7137: PetscObjectReference((PetscObject)fdm);
7138: DMDestroy(&dm->fineMesh);
7139: dm->fineMesh = fdm;
7140: return(0);
7141: }
7143: /*=== DMBoundary code ===*/
7145: PetscErrorCode DMCopyBoundary(DM dm, DM dmNew)
7146: {
7147: PetscInt d;
7151: for (d = 0; d < dm->Nds; ++d) {
7152: PetscDSCopyBoundary(dm->probs[d].ds, dmNew->probs[d].ds);
7153: }
7154: return(0);
7155: }
7157: /*@C
7158: DMAddBoundary - Add a boundary condition to the model
7160: Input Parameters:
7161: + dm - The DM, with a PetscDS that matches the problem being constrained
7162: . type - The type of condition, e.g. DM_BC_ESSENTIAL_ANALYTIC/DM_BC_ESSENTIAL_FIELD (Dirichlet), or DM_BC_NATURAL (Neumann)
7163: . name - The BC name
7164: . labelname - The label defining constrained points
7165: . field - The field to constrain
7166: . numcomps - The number of constrained field components (0 will constrain all fields)
7167: . comps - An array of constrained component numbers
7168: . bcFunc - A pointwise function giving boundary values
7169: . numids - The number of DMLabel ids for constrained points
7170: . ids - An array of ids for constrained points
7171: - ctx - An optional user context for bcFunc
7173: Options Database Keys:
7174: + -bc_<boundary name> <num> - Overrides the boundary ids
7175: - -bc_<boundary name>_comp <num> - Overrides the boundary components
7177: Level: developer
7179: .seealso: DMGetBoundary()
7180: @*/
7181: 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)
7182: {
7183: PetscDS ds;
7188: DMGetDS(dm, &ds);
7189: PetscDSAddBoundary(ds, type,name, labelname, field, numcomps, comps, bcFunc, numids, ids, ctx);
7190: return(0);
7191: }
7193: /*@
7194: DMGetNumBoundary - Get the number of registered BC
7196: Input Parameters:
7197: . dm - The mesh object
7199: Output Parameters:
7200: . numBd - The number of BC
7202: Level: intermediate
7204: .seealso: DMAddBoundary(), DMGetBoundary()
7205: @*/
7206: PetscErrorCode DMGetNumBoundary(DM dm, PetscInt *numBd)
7207: {
7208: PetscDS ds;
7213: DMGetDS(dm, &ds);
7214: PetscDSGetNumBoundary(ds, numBd);
7215: return(0);
7216: }
7218: /*@C
7219: DMGetBoundary - Get a model boundary condition
7221: Input Parameters:
7222: + dm - The mesh object
7223: - bd - The BC number
7225: Output Parameters:
7226: + type - The type of condition, e.g. DM_BC_ESSENTIAL_ANALYTIC/DM_BC_ESSENTIAL_FIELD (Dirichlet), or DM_BC_NATURAL (Neumann)
7227: . name - The BC name
7228: . labelname - The label defining constrained points
7229: . field - The field to constrain
7230: . numcomps - The number of constrained field components
7231: . comps - An array of constrained component numbers
7232: . bcFunc - A pointwise function giving boundary values
7233: . numids - The number of DMLabel ids for constrained points
7234: . ids - An array of ids for constrained points
7235: - ctx - An optional user context for bcFunc
7237: Options Database Keys:
7238: + -bc_<boundary name> <num> - Overrides the boundary ids
7239: - -bc_<boundary name>_comp <num> - Overrides the boundary components
7241: Level: developer
7243: .seealso: DMAddBoundary()
7244: @*/
7245: 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)
7246: {
7247: PetscDS ds;
7252: DMGetDS(dm, &ds);
7253: PetscDSGetBoundary(ds, bd, type, name, labelname, field, numcomps, comps, func, numids, ids, ctx);
7254: return(0);
7255: }
7257: static PetscErrorCode DMPopulateBoundary(DM dm)
7258: {
7259: PetscDS ds;
7260: DMBoundary *lastnext;
7261: DSBoundary dsbound;
7265: DMGetDS(dm, &ds);
7266: dsbound = ds->boundary;
7267: if (dm->boundary) {
7268: DMBoundary next = dm->boundary;
7270: /* quick check to see if the PetscDS has changed */
7271: if (next->dsboundary == dsbound) return(0);
7272: /* the PetscDS has changed: tear down and rebuild */
7273: while (next) {
7274: DMBoundary b = next;
7276: next = b->next;
7277: PetscFree(b);
7278: }
7279: dm->boundary = NULL;
7280: }
7282: lastnext = &(dm->boundary);
7283: while (dsbound) {
7284: DMBoundary dmbound;
7286: PetscNew(&dmbound);
7287: dmbound->dsboundary = dsbound;
7288: DMGetLabel(dm, dsbound->labelname, &(dmbound->label));
7289: if (!dmbound->label) {PetscInfo2(dm, "DSBoundary %s wants label %s, which is not in this dm.\n",dsbound->name,dsbound->labelname);}
7290: /* push on the back instead of the front so that it is in the same order as in the PetscDS */
7291: *lastnext = dmbound;
7292: lastnext = &(dmbound->next);
7293: dsbound = dsbound->next;
7294: }
7295: return(0);
7296: }
7298: PetscErrorCode DMIsBoundaryPoint(DM dm, PetscInt point, PetscBool *isBd)
7299: {
7300: DMBoundary b;
7306: *isBd = PETSC_FALSE;
7307: DMPopulateBoundary(dm);
7308: b = dm->boundary;
7309: while (b && !(*isBd)) {
7310: DMLabel label = b->label;
7311: DSBoundary dsb = b->dsboundary;
7313: if (label) {
7314: PetscInt i;
7316: for (i = 0; i < dsb->numids && !(*isBd); ++i) {
7317: DMLabelStratumHasPoint(label, dsb->ids[i], point, isBd);
7318: }
7319: }
7320: b = b->next;
7321: }
7322: return(0);
7323: }
7325: /*@C
7326: DMProjectFunction - This projects the given function into the function space provided.
7328: Input Parameters:
7329: + dm - The DM
7330: . time - The time
7331: . funcs - The coordinate functions to evaluate, one per field
7332: . ctxs - Optional array of contexts to pass to each coordinate function. ctxs itself may be null.
7333: - mode - The insertion mode for values
7335: Output Parameter:
7336: . X - vector
7338: Calling sequence of func:
7339: $ func(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar u[], void *ctx);
7341: + dim - The spatial dimension
7342: . x - The coordinates
7343: . Nf - The number of fields
7344: . u - The output field values
7345: - ctx - optional user-defined function context
7347: Level: developer
7349: .seealso: DMComputeL2Diff()
7350: @*/
7351: PetscErrorCode DMProjectFunction(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec X)
7352: {
7353: Vec localX;
7358: DMGetLocalVector(dm, &localX);
7359: DMProjectFunctionLocal(dm, time, funcs, ctxs, mode, localX);
7360: DMLocalToGlobalBegin(dm, localX, mode, X);
7361: DMLocalToGlobalEnd(dm, localX, mode, X);
7362: DMRestoreLocalVector(dm, &localX);
7363: return(0);
7364: }
7366: PetscErrorCode DMProjectFunctionLocal(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX)
7367: {
7373: if (!dm->ops->projectfunctionlocal) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMProjectFunctionLocal",((PetscObject)dm)->type_name);
7374: (dm->ops->projectfunctionlocal) (dm, time, funcs, ctxs, mode, localX);
7375: return(0);
7376: }
7378: 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)
7379: {
7380: Vec localX;
7385: DMGetLocalVector(dm, &localX);
7386: DMProjectFunctionLabelLocal(dm, time, label, numIds, ids, Nc, comps, funcs, ctxs, mode, localX);
7387: DMLocalToGlobalBegin(dm, localX, mode, X);
7388: DMLocalToGlobalEnd(dm, localX, mode, X);
7389: DMRestoreLocalVector(dm, &localX);
7390: return(0);
7391: }
7393: 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)
7394: {
7400: if (!dm->ops->projectfunctionlabellocal) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMProjectFunctionLabelLocal",((PetscObject)dm)->type_name);
7401: (dm->ops->projectfunctionlabellocal) (dm, time, label, numIds, ids, Nc, comps, funcs, ctxs, mode, localX);
7402: return(0);
7403: }
7405: PetscErrorCode DMProjectFieldLocal(DM dm, PetscReal time, Vec localU,
7406: void (**funcs)(PetscInt, PetscInt, PetscInt,
7407: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
7408: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
7409: PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]),
7410: InsertMode mode, Vec localX)
7411: {
7418: if (!dm->ops->projectfieldlocal) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMProjectFieldLocal",((PetscObject)dm)->type_name);
7419: (dm->ops->projectfieldlocal) (dm, time, localU, funcs, mode, localX);
7420: return(0);
7421: }
7423: PetscErrorCode DMProjectFieldLabelLocal(DM dm, PetscReal time, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscInt Nc, const PetscInt comps[], Vec localU,
7424: void (**funcs)(PetscInt, PetscInt, PetscInt,
7425: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
7426: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
7427: PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]),
7428: InsertMode mode, Vec localX)
7429: {
7436: if (!dm->ops->projectfieldlabellocal) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMProjectFieldLocal",((PetscObject)dm)->type_name);
7437: (dm->ops->projectfieldlabellocal)(dm, time, label, numIds, ids, Nc, comps, localU, funcs, mode, localX);
7438: return(0);
7439: }
7441: /*@C
7442: DMComputeL2Diff - This function computes the L_2 difference between a function u and an FEM interpolant solution u_h.
7444: Input Parameters:
7445: + dm - The DM
7446: . time - The time
7447: . funcs - The functions to evaluate for each field component
7448: . ctxs - Optional array of contexts to pass to each function, or NULL.
7449: - X - The coefficient vector u_h, a global vector
7451: Output Parameter:
7452: . diff - The diff ||u - u_h||_2
7454: Level: developer
7456: .seealso: DMProjectFunction(), DMComputeL2FieldDiff(), DMComputeL2GradientDiff()
7457: @*/
7458: PetscErrorCode DMComputeL2Diff(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff)
7459: {
7465: if (!dm->ops->computel2diff) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMComputeL2Diff",((PetscObject)dm)->type_name);
7466: (dm->ops->computel2diff)(dm,time,funcs,ctxs,X,diff);
7467: return(0);
7468: }
7470: /*@C
7471: DMComputeL2GradientDiff - This function computes the L_2 difference between the gradient of a function u and an FEM interpolant solution grad u_h.
7473: Input Parameters:
7474: + dm - The DM
7475: , time - The time
7476: . funcs - The gradient functions to evaluate for each field component
7477: . ctxs - Optional array of contexts to pass to each function, or NULL.
7478: . X - The coefficient vector u_h, a global vector
7479: - n - The vector to project along
7481: Output Parameter:
7482: . diff - The diff ||(grad u - grad u_h) . n||_2
7484: Level: developer
7486: .seealso: DMProjectFunction(), DMComputeL2Diff()
7487: @*/
7488: 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)
7489: {
7495: if (!dm->ops->computel2gradientdiff) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMComputeL2GradientDiff",((PetscObject)dm)->type_name);
7496: (dm->ops->computel2gradientdiff)(dm,time,funcs,ctxs,X,n,diff);
7497: return(0);
7498: }
7500: /*@C
7501: DMComputeL2FieldDiff - This function computes the L_2 difference between a function u and an FEM interpolant solution u_h, separated into field components.
7503: Input Parameters:
7504: + dm - The DM
7505: . time - The time
7506: . funcs - The functions to evaluate for each field component
7507: . ctxs - Optional array of contexts to pass to each function, or NULL.
7508: - X - The coefficient vector u_h, a global vector
7510: Output Parameter:
7511: . diff - The array of differences, ||u^f - u^f_h||_2
7513: Level: developer
7515: .seealso: DMProjectFunction(), DMComputeL2FieldDiff(), DMComputeL2GradientDiff()
7516: @*/
7517: PetscErrorCode DMComputeL2FieldDiff(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal diff[])
7518: {
7524: if (!dm->ops->computel2fielddiff) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMComputeL2FieldDiff",((PetscObject)dm)->type_name);
7525: (dm->ops->computel2fielddiff)(dm,time,funcs,ctxs,X,diff);
7526: return(0);
7527: }
7529: /*@C
7530: DMAdaptLabel - Adapt a dm based on a label with values interpreted as coarsening and refining flags. Specific implementations of DM maybe have
7531: specialized flags, but all implementations should accept flag values DM_ADAPT_DETERMINE, DM_ADAPT_KEEP, DM_ADAPT_REFINE, and DM_ADAPT_COARSEN.
7533: Collective on dm
7535: Input parameters:
7536: + dm - the pre-adaptation DM object
7537: - label - label with the flags
7539: Output parameters:
7540: . dmAdapt - the adapted DM object: may be NULL if an adapted DM could not be produced.
7542: Level: intermediate
7544: .seealso: DMAdaptMetric(), DMCoarsen(), DMRefine()
7545: @*/
7546: PetscErrorCode DMAdaptLabel(DM dm, DMLabel label, DM *dmAdapt)
7547: {
7554: *dmAdapt = NULL;
7555: if (!dm->ops->adaptlabel) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMAdaptLabel",((PetscObject)dm)->type_name);
7556: (dm->ops->adaptlabel)(dm, label, dmAdapt);
7557: return(0);
7558: }
7560: /*@C
7561: DMAdaptMetric - Generates a mesh adapted to the specified metric field using the pragmatic library.
7563: Input Parameters:
7564: + dm - The DM object
7565: . metric - The metric to which the mesh is adapted, defined vertex-wise.
7566: - 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_".
7568: Output Parameter:
7569: . dmAdapt - Pointer to the DM object containing the adapted mesh
7571: Note: The label in the adapted mesh will be registered under the name of the input DMLabel object
7573: Level: advanced
7575: .seealso: DMAdaptLabel(), DMCoarsen(), DMRefine()
7576: @*/
7577: PetscErrorCode DMAdaptMetric(DM dm, Vec metric, DMLabel bdLabel, DM *dmAdapt)
7578: {
7586: *dmAdapt = NULL;
7587: if (!dm->ops->adaptmetric) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMAdaptMetric",((PetscObject)dm)->type_name);
7588: (dm->ops->adaptmetric)(dm, metric, bdLabel, dmAdapt);
7589: return(0);
7590: }
7592: /*@C
7593: DMGetNeighbors - Gets an array containing the MPI rank of all the processes neighbors
7595: Not Collective
7597: Input Parameter:
7598: . dm - The DM
7600: Output Parameter:
7601: . nranks - the number of neighbours
7602: . ranks - the neighbors ranks
7604: Notes:
7605: Do not free the array, it is freed when the DM is destroyed.
7607: Level: beginner
7609: .seealso: DMDAGetNeighbors(), PetscSFGetRanks()
7610: @*/
7611: PetscErrorCode DMGetNeighbors(DM dm,PetscInt *nranks,const PetscMPIInt *ranks[])
7612: {
7617: if (!dm->ops->getneighbors) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMGetNeighbors",((PetscObject)dm)->type_name);
7618: (dm->ops->getneighbors)(dm,nranks,ranks);
7619: return(0);
7620: }
7622: #include <petsc/private/matimpl.h> /* Needed because of coloring->ctype below */
7624: /*
7625: Converts the input vector to a ghosted vector and then calls the standard coloring code.
7626: This has be a different function because it requires DM which is not defined in the Mat library
7627: */
7628: PetscErrorCode MatFDColoringApply_AIJDM(Mat J,MatFDColoring coloring,Vec x1,void *sctx)
7629: {
7633: if (coloring->ctype == IS_COLORING_LOCAL) {
7634: Vec x1local;
7635: DM dm;
7636: MatGetDM(J,&dm);
7637: if (!dm) SETERRQ(PetscObjectComm((PetscObject)J),PETSC_ERR_ARG_INCOMP,"IS_COLORING_LOCAL requires a DM");
7638: DMGetLocalVector(dm,&x1local);
7639: DMGlobalToLocalBegin(dm,x1,INSERT_VALUES,x1local);
7640: DMGlobalToLocalEnd(dm,x1,INSERT_VALUES,x1local);
7641: x1 = x1local;
7642: }
7643: MatFDColoringApply_AIJ(J,coloring,x1,sctx);
7644: if (coloring->ctype == IS_COLORING_LOCAL) {
7645: DM dm;
7646: MatGetDM(J,&dm);
7647: DMRestoreLocalVector(dm,&x1);
7648: }
7649: return(0);
7650: }
7652: /*@
7653: MatFDColoringUseDM - allows a MatFDColoring object to use the DM associated with the matrix to use a IS_COLORING_LOCAL coloring
7655: Input Parameter:
7656: . coloring - the MatFDColoring object
7658: Developer Notes:
7659: this routine exists because the PETSc Mat library does not know about the DM objects
7661: Level: advanced
7663: .seealso: MatFDColoring, MatFDColoringCreate(), ISColoringType
7664: @*/
7665: PetscErrorCode MatFDColoringUseDM(Mat coloring,MatFDColoring fdcoloring)
7666: {
7668: coloring->ops->fdcoloringapply = MatFDColoringApply_AIJDM;
7669: return(0);
7670: }
7672: /*@
7673: DMGetCompatibility - determine if two DMs are compatible
7675: Collective
7677: Input Parameters:
7678: + dm - the first DM
7679: - dm2 - the second DM
7681: Output Parameters:
7682: + compatible - whether or not the two DMs are compatible
7683: - set - whether or not the compatible value was set
7685: Notes:
7686: Two DMs are deemed compatible if they represent the same parallel decomposition
7687: of the same topology. This implies that the section (field data) on one
7688: "makes sense" with respect to the topology and parallel decomposition of the other.
7689: Loosely speaking, compatible DMs represent the same domain and parallel
7690: decomposition, but hold different data.
7692: Typically, one would confirm compatibility if intending to simultaneously iterate
7693: over a pair of vectors obtained from different DMs.
7695: For example, two DMDA objects are compatible if they have the same local
7696: and global sizes and the same stencil width. They can have different numbers
7697: of degrees of freedom per node. Thus, one could use the node numbering from
7698: either DM in bounds for a loop over vectors derived from either DM.
7700: Consider the operation of summing data living on a 2-dof DMDA to data living
7701: on a 1-dof DMDA, which should be compatible, as in the following snippet.
7702: .vb
7703: ...
7704: DMGetCompatibility(da1,da2,&compatible,&set);
7705: if (set && compatible) {
7706: DMDAVecGetArrayDOF(da1,vec1,&arr1);
7707: DMDAVecGetArrayDOF(da2,vec2,&arr2);
7708: DMDAGetCorners(da1,&x,&y,NULL,&m,&n,NULL);
7709: for (j=y; j<y+n; ++j) {
7710: for (i=x; i<x+m, ++i) {
7711: arr1[j][i][0] = arr2[j][i][0] + arr2[j][i][1];
7712: }
7713: }
7714: DMDAVecRestoreArrayDOF(da1,vec1,&arr1);
7715: DMDAVecRestoreArrayDOF(da2,vec2,&arr2);
7716: } else {
7717: SETERRQ(PetscObjectComm((PetscObject)da1,PETSC_ERR_ARG_INCOMP,"DMDA objects incompatible");
7718: }
7719: ...
7720: .ve
7722: Checking compatibility might be expensive for a given implementation of DM,
7723: or might be impossible to unambiguously confirm or deny. For this reason,
7724: this function may decline to determine compatibility, and hence users should
7725: always check the "set" output parameter.
7727: A DM is always compatible with itself.
7729: In the current implementation, DMs which live on "unequal" communicators
7730: (MPI_UNEQUAL in the terminology of MPI_Comm_compare()) are always deemed
7731: incompatible.
7733: This function is labeled "Collective," as information about all subdomains
7734: is required on each rank. However, in DM implementations which store all this
7735: information locally, this function may be merely "Logically Collective".
7737: Developer Notes:
7738: Compatibility is assumed to be a symmetric concept; DM A is compatible with DM B
7739: iff B is compatible with A. Thus, this function checks the implementations
7740: of both dm and dm2 (if they are of different types), attempting to determine
7741: compatibility. It is left to DM implementers to ensure that symmetry is
7742: preserved. The simplest way to do this is, when implementing type-specific
7743: logic for this function, is to check for existing logic in the implementation
7744: of other DM types and let *set = PETSC_FALSE if found.
7746: Level: advanced
7748: .seealso: DM, DMDACreateCompatibleDMDA(), DMStagCreateCompatibleDMStag()
7749: @*/
7751: PetscErrorCode DMGetCompatibility(DM dm,DM dm2,PetscBool *compatible,PetscBool *set)
7752: {
7754: PetscMPIInt compareResult;
7755: DMType type,type2;
7756: PetscBool sameType;
7762: /* Declare a DM compatible with itself */
7763: if (dm == dm2) {
7764: *set = PETSC_TRUE;
7765: *compatible = PETSC_TRUE;
7766: return(0);
7767: }
7769: /* Declare a DM incompatible with a DM that lives on an "unequal"
7770: communicator. Note that this does not preclude compatibility with
7771: DMs living on "congruent" or "similar" communicators, but this must be
7772: determined by the implementation-specific logic */
7773: MPI_Comm_compare(PetscObjectComm((PetscObject)dm),PetscObjectComm((PetscObject)dm2),&compareResult);
7774: if (compareResult == MPI_UNEQUAL) {
7775: *set = PETSC_TRUE;
7776: *compatible = PETSC_FALSE;
7777: return(0);
7778: }
7780: /* Pass to the implementation-specific routine, if one exists. */
7781: if (dm->ops->getcompatibility) {
7782: (*dm->ops->getcompatibility)(dm,dm2,compatible,set);
7783: if (*set) {
7784: return(0);
7785: }
7786: }
7788: /* If dm and dm2 are of different types, then attempt to check compatibility
7789: with an implementation of this function from dm2 */
7790: DMGetType(dm,&type);
7791: DMGetType(dm2,&type2);
7792: PetscStrcmp(type,type2,&sameType);
7793: if (!sameType && dm2->ops->getcompatibility) {
7794: (*dm2->ops->getcompatibility)(dm2,dm,compatible,set); /* Note argument order */
7795: } else {
7796: *set = PETSC_FALSE;
7797: }
7798: return(0);
7799: }