Actual source code: dm.c
petsc-3.9.4 2018-09-11
1: #include <petsc/private/dmimpl.h>
2: #include <petsc/private/dmlabelimpl.h>
3: #include <petsc/private/petscdsimpl.h>
4: #include <petscdmplex.h>
5: #include <petscsf.h>
6: #include <petscds.h>
8: PetscClassId DM_CLASSID;
9: PetscLogEvent DM_Convert, DM_GlobalToLocal, DM_LocalToGlobal, DM_LocalToLocal, DM_LocatePoints, DM_Coarsen, DM_Refine, DM_CreateInterpolation, DM_CreateRestriction;
11: const char *const DMBoundaryTypes[] = {"NONE","GHOSTED","MIRROR","PERIODIC","TWIST","DM_BOUNDARY_",0};
13: static PetscErrorCode DMHasCreateInjection_Default(DM dm, PetscBool *flg)
14: {
18: *flg = PETSC_FALSE;
19: return(0);
20: }
22: /*@
23: DMCreate - Creates an empty DM object. The type can then be set with DMSetType().
25: If you never call DMSetType() it will generate an
26: error when you try to use the vector.
28: Collective on MPI_Comm
30: Input Parameter:
31: . comm - The communicator for the DM object
33: Output Parameter:
34: . dm - The DM object
36: Level: beginner
38: .seealso: DMSetType(), DMDA, DMSLICED, DMCOMPOSITE, DMPLEX, DMMOAB, DMNETWORK
39: @*/
40: PetscErrorCode DMCreate(MPI_Comm comm,DM *dm)
41: {
42: DM v;
47: *dm = NULL;
48: PetscSysInitializePackage();
49: VecInitializePackage();
50: MatInitializePackage();
51: DMInitializePackage();
53: PetscHeaderCreate(v, DM_CLASSID, "DM", "Distribution Manager", "DM", comm, DMDestroy, DMView);
55: v->ltogmap = NULL;
56: v->bs = 1;
57: v->coloringtype = IS_COLORING_GLOBAL;
58: PetscSFCreate(comm, &v->sf);
59: PetscSFCreate(comm, &v->defaultSF);
60: v->labels = NULL;
61: v->depthLabel = NULL;
62: v->defaultSection = NULL;
63: v->defaultGlobalSection = NULL;
64: v->defaultConstraintSection = NULL;
65: v->defaultConstraintMat = NULL;
66: v->L = NULL;
67: v->maxCell = NULL;
68: v->bdtype = NULL;
69: v->dimEmbed = PETSC_DEFAULT;
70: {
71: PetscInt i;
72: for (i = 0; i < 10; ++i) {
73: v->nullspaceConstructors[i] = NULL;
74: }
75: }
76: PetscDSCreate(comm, &v->prob);
77: v->dmBC = NULL;
78: v->coarseMesh = NULL;
79: v->outputSequenceNum = -1;
80: v->outputSequenceVal = 0.0;
81: DMSetVecType(v,VECSTANDARD);
82: DMSetMatType(v,MATAIJ);
83: PetscNew(&(v->labels));
84: v->labels->refct = 1;
86: v->ops->hascreateinjection = DMHasCreateInjection_Default;
88: *dm = v;
89: return(0);
90: }
92: /*@
93: DMClone - Creates a DM object with the same topology as the original.
95: Collective on MPI_Comm
97: Input Parameter:
98: . dm - The original DM object
100: Output Parameter:
101: . newdm - The new DM object
103: Level: beginner
105: .keywords: DM, topology, create
106: @*/
107: PetscErrorCode DMClone(DM dm, DM *newdm)
108: {
109: PetscSF sf;
110: Vec coords;
111: void *ctx;
112: PetscInt dim, cdim;
118: DMCreate(PetscObjectComm((PetscObject) dm), newdm);
119: PetscFree((*newdm)->labels);
120: dm->labels->refct++;
121: (*newdm)->labels = dm->labels;
122: (*newdm)->depthLabel = dm->depthLabel;
123: DMGetDimension(dm, &dim);
124: DMSetDimension(*newdm, dim);
125: if (dm->ops->clone) {
126: (*dm->ops->clone)(dm, newdm);
127: }
128: (*newdm)->setupcalled = dm->setupcalled;
129: DMGetPointSF(dm, &sf);
130: DMSetPointSF(*newdm, sf);
131: DMGetApplicationContext(dm, &ctx);
132: DMSetApplicationContext(*newdm, ctx);
133: if (dm->coordinateDM) {
134: DM ncdm;
135: PetscSection cs;
136: PetscInt pEnd = -1, pEndMax = -1;
138: DMGetDefaultSection(dm->coordinateDM, &cs);
139: if (cs) {PetscSectionGetChart(cs, NULL, &pEnd);}
140: MPI_Allreduce(&pEnd,&pEndMax,1,MPIU_INT,MPI_MAX,PetscObjectComm((PetscObject)dm));
141: if (pEndMax >= 0) {
142: DMClone(dm->coordinateDM, &ncdm);
143: DMSetDefaultSection(ncdm, cs);
144: DMSetCoordinateDM(*newdm, ncdm);
145: DMDestroy(&ncdm);
146: }
147: }
148: DMGetCoordinateDim(dm, &cdim);
149: DMSetCoordinateDim(*newdm, cdim);
150: DMGetCoordinatesLocal(dm, &coords);
151: if (coords) {
152: DMSetCoordinatesLocal(*newdm, coords);
153: } else {
154: DMGetCoordinates(dm, &coords);
155: if (coords) {DMSetCoordinates(*newdm, coords);}
156: }
157: {
158: PetscBool isper;
159: const PetscReal *maxCell, *L;
160: const DMBoundaryType *bd;
161: DMGetPeriodicity(dm, &isper, &maxCell, &L, &bd);
162: DMSetPeriodicity(*newdm, isper, maxCell, L, bd);
163: }
164: return(0);
165: }
167: /*@C
168: DMSetVecType - Sets the type of vector created with DMCreateLocalVector() and DMCreateGlobalVector()
170: Logically Collective on DM
172: Input Parameter:
173: + da - initial distributed array
174: . ctype - the vector type, currently either VECSTANDARD, VECCUDA, or VECVIENNACL
176: Options Database:
177: . -dm_vec_type ctype
179: Level: intermediate
181: .seealso: DMCreate(), DMDestroy(), DM, DMDAInterpolationType, VecType, DMGetVecType()
182: @*/
183: PetscErrorCode DMSetVecType(DM da,VecType ctype)
184: {
189: PetscFree(da->vectype);
190: PetscStrallocpy(ctype,(char**)&da->vectype);
191: return(0);
192: }
194: /*@C
195: DMGetVecType - Gets the type of vector created with DMCreateLocalVector() and DMCreateGlobalVector()
197: Logically Collective on DM
199: Input Parameter:
200: . da - initial distributed array
202: Output Parameter:
203: . ctype - the vector type
205: Level: intermediate
207: .seealso: DMCreate(), DMDestroy(), DM, DMDAInterpolationType, VecType
208: @*/
209: PetscErrorCode DMGetVecType(DM da,VecType *ctype)
210: {
213: *ctype = da->vectype;
214: return(0);
215: }
217: /*@
218: VecGetDM - Gets the DM defining the data layout of the vector
220: Not collective
222: Input Parameter:
223: . v - The Vec
225: Output Parameter:
226: . dm - The DM
228: Level: intermediate
230: .seealso: VecSetDM(), DMGetLocalVector(), DMGetGlobalVector(), DMSetVecType()
231: @*/
232: PetscErrorCode VecGetDM(Vec v, DM *dm)
233: {
239: PetscObjectQuery((PetscObject) v, "__PETSc_dm", (PetscObject*) dm);
240: return(0);
241: }
243: /*@
244: VecSetDM - Sets the DM defining the data layout of the vector.
246: Not collective
248: Input Parameters:
249: + v - The Vec
250: - dm - The DM
252: 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.
254: Level: intermediate
256: .seealso: VecGetDM(), DMGetLocalVector(), DMGetGlobalVector(), DMSetVecType()
257: @*/
258: PetscErrorCode VecSetDM(Vec v, DM dm)
259: {
265: PetscObjectCompose((PetscObject) v, "__PETSc_dm", (PetscObject) dm);
266: return(0);
267: }
269: /*@C
270: DMSetISColoringType - Sets the type of coloring, global or local, that is created by the DM
272: Logically Collective on DM
274: Input Parameters:
275: + dm - the DM context
276: - ctype - the matrix type
278: Options Database:
279: . -dm_is_coloring_type - global or local
281: Level: intermediate
283: .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMCreateMatrix(), DMSetMatrixPreallocateOnly(), MatType, DMGetMatType(),
284: DMGetISColoringType()
285: @*/
286: PetscErrorCode DMSetISColoringType(DM dm,ISColoringType ctype)
287: {
290: dm->coloringtype = ctype;
291: return(0);
292: }
294: /*@C
295: DMGetISColoringType - Gets the type of coloring, global or local, that is created by the DM
297: Logically Collective on DM
299: Input Parameter:
300: . dm - the DM context
302: Output Parameter:
303: . ctype - the matrix type
305: Options Database:
306: . -dm_is_coloring_type - global or local
308: Level: intermediate
310: .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMCreateMatrix(), DMSetMatrixPreallocateOnly(), MatType, DMGetMatType(),
311: DMGetISColoringType()
312: @*/
313: PetscErrorCode DMGetISColoringType(DM dm,ISColoringType *ctype)
314: {
317: *ctype = dm->coloringtype;
318: return(0);
319: }
321: /*@C
322: DMSetMatType - Sets the type of matrix created with DMCreateMatrix()
324: Logically Collective on DM
326: Input Parameters:
327: + dm - the DM context
328: - ctype - the matrix type
330: Options Database:
331: . -dm_mat_type ctype
333: Level: intermediate
335: .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMCreateMatrix(), DMSetMatrixPreallocateOnly(), MatType, DMGetMatType()
336: @*/
337: PetscErrorCode DMSetMatType(DM dm,MatType ctype)
338: {
343: PetscFree(dm->mattype);
344: PetscStrallocpy(ctype,(char**)&dm->mattype);
345: return(0);
346: }
348: /*@C
349: DMGetMatType - Gets the type of matrix created with DMCreateMatrix()
351: Logically Collective on DM
353: Input Parameter:
354: . dm - the DM context
356: Output Parameter:
357: . ctype - the matrix type
359: Options Database:
360: . -dm_mat_type ctype
362: Level: intermediate
364: .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMCreateMatrix(), DMSetMatrixPreallocateOnly(), MatType, DMSetMatType()
365: @*/
366: PetscErrorCode DMGetMatType(DM dm,MatType *ctype)
367: {
370: *ctype = dm->mattype;
371: return(0);
372: }
374: /*@
375: MatGetDM - Gets the DM defining the data layout of the matrix
377: Not collective
379: Input Parameter:
380: . A - The Mat
382: Output Parameter:
383: . dm - The DM
385: Level: intermediate
387: Developer Note: Since the Mat class doesn't know about the DM class the DM object is associated with
388: the Mat through a PetscObjectCompose() operation
390: .seealso: MatSetDM(), DMCreateMatrix(), DMSetMatType()
391: @*/
392: PetscErrorCode MatGetDM(Mat A, DM *dm)
393: {
399: PetscObjectQuery((PetscObject) A, "__PETSc_dm", (PetscObject*) dm);
400: return(0);
401: }
403: /*@
404: MatSetDM - Sets the DM defining the data layout of the matrix
406: Not collective
408: Input Parameters:
409: + A - The Mat
410: - dm - The DM
412: Level: intermediate
414: Developer Note: Since the Mat class doesn't know about the DM class the DM object is associated with
415: the Mat through a PetscObjectCompose() operation
418: .seealso: MatGetDM(), DMCreateMatrix(), DMSetMatType()
419: @*/
420: PetscErrorCode MatSetDM(Mat A, DM dm)
421: {
427: PetscObjectCompose((PetscObject) A, "__PETSc_dm", (PetscObject) dm);
428: return(0);
429: }
431: /*@C
432: DMSetOptionsPrefix - Sets the prefix used for searching for all
433: DM options in the database.
435: Logically Collective on DM
437: Input Parameter:
438: + da - the DM context
439: - prefix - the prefix to prepend to all option names
441: Notes:
442: A hyphen (-) must NOT be given at the beginning of the prefix name.
443: The first character of all runtime options is AUTOMATICALLY the hyphen.
445: Level: advanced
447: .keywords: DM, set, options, prefix, database
449: .seealso: DMSetFromOptions()
450: @*/
451: PetscErrorCode DMSetOptionsPrefix(DM dm,const char prefix[])
452: {
457: PetscObjectSetOptionsPrefix((PetscObject)dm,prefix);
458: if (dm->sf) {
459: PetscObjectSetOptionsPrefix((PetscObject)dm->sf,prefix);
460: }
461: if (dm->defaultSF) {
462: PetscObjectSetOptionsPrefix((PetscObject)dm->defaultSF,prefix);
463: }
464: return(0);
465: }
467: /*@C
468: DMAppendOptionsPrefix - Appends to the prefix used for searching for all
469: DM options in the database.
471: Logically Collective on DM
473: Input Parameters:
474: + dm - the DM context
475: - prefix - the prefix string to prepend to all DM option requests
477: Notes:
478: A hyphen (-) must NOT be given at the beginning of the prefix name.
479: The first character of all runtime options is AUTOMATICALLY the hyphen.
481: Level: advanced
483: .keywords: DM, append, options, prefix, database
485: .seealso: DMSetOptionsPrefix(), DMGetOptionsPrefix()
486: @*/
487: PetscErrorCode DMAppendOptionsPrefix(DM dm,const char prefix[])
488: {
493: PetscObjectAppendOptionsPrefix((PetscObject)dm,prefix);
494: return(0);
495: }
497: /*@C
498: DMGetOptionsPrefix - Gets the prefix used for searching for all
499: DM options in the database.
501: Not Collective
503: Input Parameters:
504: . dm - the DM context
506: Output Parameters:
507: . prefix - pointer to the prefix string used is returned
509: Notes: On the fortran side, the user should pass in a string 'prefix' of
510: sufficient length to hold the prefix.
512: Level: advanced
514: .keywords: DM, set, options, prefix, database
516: .seealso: DMSetOptionsPrefix(), DMAppendOptionsPrefix()
517: @*/
518: PetscErrorCode DMGetOptionsPrefix(DM dm,const char *prefix[])
519: {
524: PetscObjectGetOptionsPrefix((PetscObject)dm,prefix);
525: return(0);
526: }
528: static PetscErrorCode DMCountNonCyclicReferences(DM dm, PetscBool recurseCoarse, PetscBool recurseFine, PetscInt *ncrefct)
529: {
530: PetscInt i, refct = ((PetscObject) dm)->refct;
531: DMNamedVecLink nlink;
535: *ncrefct = 0;
536: /* count all the circular references of DM and its contained Vecs */
537: for (i=0; i<DM_MAX_WORK_VECTORS; i++) {
538: if (dm->localin[i]) refct--;
539: if (dm->globalin[i]) refct--;
540: }
541: for (nlink=dm->namedglobal; nlink; nlink=nlink->next) refct--;
542: for (nlink=dm->namedlocal; nlink; nlink=nlink->next) refct--;
543: if (dm->x) {
544: DM obj;
545: VecGetDM(dm->x, &obj);
546: if (obj == dm) refct--;
547: }
548: if (dm->coarseMesh && dm->coarseMesh->fineMesh == dm) {
549: refct--;
550: if (recurseCoarse) {
551: PetscInt coarseCount;
553: DMCountNonCyclicReferences(dm->coarseMesh, PETSC_TRUE, PETSC_FALSE,&coarseCount);
554: refct += coarseCount;
555: }
556: }
557: if (dm->fineMesh && dm->fineMesh->coarseMesh == dm) {
558: refct--;
559: if (recurseFine) {
560: PetscInt fineCount;
562: DMCountNonCyclicReferences(dm->fineMesh, PETSC_FALSE, PETSC_TRUE,&fineCount);
563: refct += fineCount;
564: }
565: }
566: *ncrefct = refct;
567: return(0);
568: }
570: PetscErrorCode DMDestroyLabelLinkList(DM dm)
571: {
575: if (!--(dm->labels->refct)) {
576: DMLabelLink next = dm->labels->next;
578: /* destroy the labels */
579: while (next) {
580: DMLabelLink tmp = next->next;
582: DMLabelDestroy(&next->label);
583: PetscFree(next);
584: next = tmp;
585: }
586: PetscFree(dm->labels);
587: }
588: return(0);
589: }
591: /*@
592: DMDestroy - Destroys a vector packer or DM.
594: Collective on DM
596: Input Parameter:
597: . dm - the DM object to destroy
599: Level: developer
601: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix()
603: @*/
604: PetscErrorCode DMDestroy(DM *dm)
605: {
606: PetscInt i, cnt;
607: DMNamedVecLink nlink,nnext;
611: if (!*dm) return(0);
614: /* count all non-cyclic references in the doubly-linked list of coarse<->fine meshes */
615: DMCountNonCyclicReferences(*dm,PETSC_TRUE,PETSC_TRUE,&cnt);
616: --((PetscObject)(*dm))->refct;
617: if (--cnt > 0) {*dm = 0; return(0);}
618: /*
619: Need this test because the dm references the vectors that
620: reference the dm, so destroying the dm calls destroy on the
621: vectors that cause another destroy on the dm
622: */
623: if (((PetscObject)(*dm))->refct < 0) return(0);
624: ((PetscObject) (*dm))->refct = 0;
625: for (i=0; i<DM_MAX_WORK_VECTORS; i++) {
626: if ((*dm)->localout[i]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Destroying a DM that has a local vector obtained with DMGetLocalVector()");
627: VecDestroy(&(*dm)->localin[i]);
628: }
629: nnext=(*dm)->namedglobal;
630: (*dm)->namedglobal = NULL;
631: for (nlink=nnext; nlink; nlink=nnext) { /* Destroy the named vectors */
632: nnext = nlink->next;
633: if (nlink->status != DMVEC_STATUS_IN) SETERRQ1(((PetscObject)*dm)->comm,PETSC_ERR_ARG_WRONGSTATE,"DM still has Vec named '%s' checked out",nlink->name);
634: PetscFree(nlink->name);
635: VecDestroy(&nlink->X);
636: PetscFree(nlink);
637: }
638: nnext=(*dm)->namedlocal;
639: (*dm)->namedlocal = NULL;
640: for (nlink=nnext; nlink; nlink=nnext) { /* Destroy the named local vectors */
641: nnext = nlink->next;
642: if (nlink->status != DMVEC_STATUS_IN) SETERRQ1(((PetscObject)*dm)->comm,PETSC_ERR_ARG_WRONGSTATE,"DM still has Vec named '%s' checked out",nlink->name);
643: PetscFree(nlink->name);
644: VecDestroy(&nlink->X);
645: PetscFree(nlink);
646: }
648: /* Destroy the list of hooks */
649: {
650: DMCoarsenHookLink link,next;
651: for (link=(*dm)->coarsenhook; link; link=next) {
652: next = link->next;
653: PetscFree(link);
654: }
655: (*dm)->coarsenhook = NULL;
656: }
657: {
658: DMRefineHookLink link,next;
659: for (link=(*dm)->refinehook; link; link=next) {
660: next = link->next;
661: PetscFree(link);
662: }
663: (*dm)->refinehook = NULL;
664: }
665: {
666: DMSubDomainHookLink link,next;
667: for (link=(*dm)->subdomainhook; link; link=next) {
668: next = link->next;
669: PetscFree(link);
670: }
671: (*dm)->subdomainhook = NULL;
672: }
673: {
674: DMGlobalToLocalHookLink link,next;
675: for (link=(*dm)->gtolhook; link; link=next) {
676: next = link->next;
677: PetscFree(link);
678: }
679: (*dm)->gtolhook = NULL;
680: }
681: {
682: DMLocalToGlobalHookLink link,next;
683: for (link=(*dm)->ltoghook; link; link=next) {
684: next = link->next;
685: PetscFree(link);
686: }
687: (*dm)->ltoghook = NULL;
688: }
689: /* Destroy the work arrays */
690: {
691: DMWorkLink link,next;
692: if ((*dm)->workout) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Work array still checked out");
693: for (link=(*dm)->workin; link; link=next) {
694: next = link->next;
695: PetscFree(link->mem);
696: PetscFree(link);
697: }
698: (*dm)->workin = NULL;
699: }
700: if (!--((*dm)->labels->refct)) {
701: DMLabelLink next = (*dm)->labels->next;
703: /* destroy the labels */
704: while (next) {
705: DMLabelLink tmp = next->next;
707: DMLabelDestroy(&next->label);
708: PetscFree(next);
709: next = tmp;
710: }
711: PetscFree((*dm)->labels);
712: }
713: {
714: DMBoundary next = (*dm)->boundary;
715: while (next) {
716: DMBoundary b = next;
718: next = b->next;
719: PetscFree(b);
720: }
721: }
723: PetscObjectDestroy(&(*dm)->dmksp);
724: PetscObjectDestroy(&(*dm)->dmsnes);
725: PetscObjectDestroy(&(*dm)->dmts);
727: if ((*dm)->ctx && (*dm)->ctxdestroy) {
728: (*(*dm)->ctxdestroy)(&(*dm)->ctx);
729: }
730: VecDestroy(&(*dm)->x);
731: MatFDColoringDestroy(&(*dm)->fd);
732: DMClearGlobalVectors(*dm);
733: ISLocalToGlobalMappingDestroy(&(*dm)->ltogmap);
734: PetscFree((*dm)->vectype);
735: PetscFree((*dm)->mattype);
737: PetscSectionDestroy(&(*dm)->defaultSection);
738: PetscSectionDestroy(&(*dm)->defaultGlobalSection);
739: PetscLayoutDestroy(&(*dm)->map);
740: PetscSectionDestroy(&(*dm)->defaultConstraintSection);
741: MatDestroy(&(*dm)->defaultConstraintMat);
742: PetscSFDestroy(&(*dm)->sf);
743: PetscSFDestroy(&(*dm)->defaultSF);
744: if ((*dm)->useNatural) {
745: if ((*dm)->sfNatural) {
746: PetscSFDestroy(&(*dm)->sfNatural);
747: }
748: PetscObjectDereference((PetscObject) (*dm)->sfMigration);
749: }
750: if ((*dm)->coarseMesh && (*dm)->coarseMesh->fineMesh == *dm) {
751: DMSetFineDM((*dm)->coarseMesh,NULL);
752: }
753: DMDestroy(&(*dm)->coarseMesh);
754: if ((*dm)->fineMesh && (*dm)->fineMesh->coarseMesh == *dm) {
755: DMSetCoarseDM((*dm)->fineMesh,NULL);
756: }
757: DMDestroy(&(*dm)->fineMesh);
758: DMDestroy(&(*dm)->coordinateDM);
759: VecDestroy(&(*dm)->coordinates);
760: VecDestroy(&(*dm)->coordinatesLocal);
761: PetscFree3((*dm)->L,(*dm)->maxCell,(*dm)->bdtype);
763: PetscDSDestroy(&(*dm)->prob);
764: DMDestroy(&(*dm)->dmBC);
765: /* if memory was published with SAWs then destroy it */
766: PetscObjectSAWsViewOff((PetscObject)*dm);
768: (*(*dm)->ops->destroy)(*dm);
769: /* We do not destroy (*dm)->data here so that we can reference count backend objects */
770: PetscHeaderDestroy(dm);
771: return(0);
772: }
774: /*@
775: DMSetUp - sets up the data structures inside a DM object
777: Collective on DM
779: Input Parameter:
780: . dm - the DM object to setup
782: Level: developer
784: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix()
786: @*/
787: PetscErrorCode DMSetUp(DM dm)
788: {
793: if (dm->setupcalled) return(0);
794: if (dm->ops->setup) {
795: (*dm->ops->setup)(dm);
796: }
797: dm->setupcalled = PETSC_TRUE;
798: return(0);
799: }
801: /*@
802: DMSetFromOptions - sets parameters in a DM from the options database
804: Collective on DM
806: Input Parameter:
807: . dm - the DM object to set options for
809: Options Database:
810: + -dm_preallocate_only - Only preallocate the matrix for DMCreateMatrix(), but do not fill it with zeros
811: . -dm_vec_type <type> - type of vector to create inside DM
812: . -dm_mat_type <type> - type of matrix to create inside DM
813: - -dm_is_coloring_type - <global or local>
815: Level: developer
817: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix()
819: @*/
820: PetscErrorCode DMSetFromOptions(DM dm)
821: {
822: char typeName[256];
823: PetscBool flg;
828: if (dm->prob) {
829: PetscDSSetFromOptions(dm->prob);
830: }
831: if (dm->sf) {
832: PetscSFSetFromOptions(dm->sf);
833: }
834: if (dm->defaultSF) {
835: PetscSFSetFromOptions(dm->defaultSF);
836: }
837: PetscObjectOptionsBegin((PetscObject)dm);
838: PetscOptionsBool("-dm_preallocate_only","only preallocate matrix, but do not set column indices","DMSetMatrixPreallocateOnly",dm->prealloc_only,&dm->prealloc_only,NULL);
839: PetscOptionsFList("-dm_vec_type","Vector type used for created vectors","DMSetVecType",VecList,dm->vectype,typeName,256,&flg);
840: if (flg) {
841: DMSetVecType(dm,typeName);
842: }
843: PetscOptionsFList("-dm_mat_type","Matrix type used for created matrices","DMSetMatType",MatList,dm->mattype ? dm->mattype : typeName,typeName,sizeof(typeName),&flg);
844: if (flg) {
845: DMSetMatType(dm,typeName);
846: }
847: PetscOptionsEnum("-dm_is_coloring_type","Global or local coloring of Jacobian","DMSetISColoringType",ISColoringTypes,(PetscEnum)dm->coloringtype,(PetscEnum*)&dm->coloringtype,NULL);
848: if (dm->ops->setfromoptions) {
849: (*dm->ops->setfromoptions)(PetscOptionsObject,dm);
850: }
851: /* process any options handlers added with PetscObjectAddOptionsHandler() */
852: PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject) dm);
853: PetscOptionsEnd();
854: return(0);
855: }
857: /*@C
858: DMView - Views a DM
860: Collective on DM
862: Input Parameter:
863: + dm - the DM object to view
864: - v - the viewer
866: Level: beginner
868: .seealso DMDestroy(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix()
870: @*/
871: PetscErrorCode DMView(DM dm,PetscViewer v)
872: {
873: PetscErrorCode ierr;
874: PetscBool isbinary;
875: PetscMPIInt size;
876: PetscViewerFormat format;
880: if (!v) {
881: PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)dm),&v);
882: }
883: PetscViewerGetFormat(v,&format);
884: MPI_Comm_size(PetscObjectComm((PetscObject)dm),&size);
885: if (size == 1 && format == PETSC_VIEWER_LOAD_BALANCE) return(0);
886: PetscObjectPrintClassNamePrefixType((PetscObject)dm,v);
887: PetscObjectTypeCompare((PetscObject)v,PETSCVIEWERBINARY,&isbinary);
888: if (isbinary) {
889: PetscInt classid = DM_FILE_CLASSID;
890: char type[256];
892: PetscViewerBinaryWrite(v,&classid,1,PETSC_INT,PETSC_FALSE);
893: PetscStrncpy(type,((PetscObject)dm)->type_name,256);
894: PetscViewerBinaryWrite(v,type,256,PETSC_CHAR,PETSC_FALSE);
895: }
896: if (dm->ops->view) {
897: (*dm->ops->view)(dm,v);
898: }
899: return(0);
900: }
902: /*@
903: DMCreateGlobalVector - Creates a global vector from a DM object
905: Collective on DM
907: Input Parameter:
908: . dm - the DM object
910: Output Parameter:
911: . vec - the global vector
913: Level: beginner
915: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix()
917: @*/
918: PetscErrorCode DMCreateGlobalVector(DM dm,Vec *vec)
919: {
924: (*dm->ops->createglobalvector)(dm,vec);
925: return(0);
926: }
928: /*@
929: DMCreateLocalVector - Creates a local vector from a DM object
931: Not Collective
933: Input Parameter:
934: . dm - the DM object
936: Output Parameter:
937: . vec - the local vector
939: Level: beginner
941: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix()
943: @*/
944: PetscErrorCode DMCreateLocalVector(DM dm,Vec *vec)
945: {
950: (*dm->ops->createlocalvector)(dm,vec);
951: return(0);
952: }
954: /*@
955: DMGetLocalToGlobalMapping - Accesses the local-to-global mapping in a DM.
957: Collective on DM
959: Input Parameter:
960: . dm - the DM that provides the mapping
962: Output Parameter:
963: . ltog - the mapping
965: Level: intermediate
967: Notes:
968: This mapping can then be used by VecSetLocalToGlobalMapping() or
969: MatSetLocalToGlobalMapping().
971: .seealso: DMCreateLocalVector()
972: @*/
973: PetscErrorCode DMGetLocalToGlobalMapping(DM dm,ISLocalToGlobalMapping *ltog)
974: {
975: PetscInt bs = -1, bsLocal[2], bsMinMax[2];
981: if (!dm->ltogmap) {
982: PetscSection section, sectionGlobal;
984: DMGetDefaultSection(dm, §ion);
985: if (section) {
986: const PetscInt *cdofs;
987: PetscInt *ltog;
988: PetscInt pStart, pEnd, n, p, k, l;
990: DMGetDefaultGlobalSection(dm, §ionGlobal);
991: PetscSectionGetChart(section, &pStart, &pEnd);
992: PetscSectionGetStorageSize(section, &n);
993: PetscMalloc1(n, <og); /* We want the local+overlap size */
994: for (p = pStart, l = 0; p < pEnd; ++p) {
995: PetscInt bdof, cdof, dof, off, c, cind = 0;
997: /* Should probably use constrained dofs */
998: PetscSectionGetDof(section, p, &dof);
999: PetscSectionGetConstraintDof(section, p, &cdof);
1000: PetscSectionGetConstraintIndices(section, p, &cdofs);
1001: PetscSectionGetOffset(sectionGlobal, p, &off);
1002: /* If you have dofs, and constraints, and they are unequal, we set the blocksize to 1 */
1003: bdof = cdof && (dof-cdof) ? 1 : dof;
1004: if (dof) {
1005: if (bs < 0) {bs = bdof;}
1006: else if (bs != bdof) {bs = 1;}
1007: }
1008: for (c = 0; c < dof; ++c, ++l) {
1009: if ((cind < cdof) && (c == cdofs[cind])) ltog[l] = off < 0 ? off-c : off+c;
1010: else ltog[l] = (off < 0 ? -(off+1) : off) + c;
1011: }
1012: }
1013: /* Must have same blocksize on all procs (some might have no points) */
1014: bsLocal[0] = bs < 0 ? PETSC_MAX_INT : bs; bsLocal[1] = bs;
1015: PetscGlobalMinMaxInt(PetscObjectComm((PetscObject) dm), bsLocal, bsMinMax);
1016: if (bsMinMax[0] != bsMinMax[1]) {bs = 1;}
1017: else {bs = bsMinMax[0];}
1018: bs = bs < 0 ? 1 : bs;
1019: /* Must reduce indices by blocksize */
1020: if (bs > 1) {
1021: for (l = 0, k = 0; l < n; l += bs, ++k) ltog[k] = ltog[l]/bs;
1022: n /= bs;
1023: }
1024: ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)dm), bs, n, ltog, PETSC_OWN_POINTER, &dm->ltogmap);
1025: PetscLogObjectParent((PetscObject)dm, (PetscObject)dm->ltogmap);
1026: } else {
1027: if (!dm->ops->getlocaltoglobalmapping) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM can not create LocalToGlobalMapping");
1028: (*dm->ops->getlocaltoglobalmapping)(dm);
1029: }
1030: }
1031: *ltog = dm->ltogmap;
1032: return(0);
1033: }
1035: /*@
1036: DMGetBlockSize - Gets the inherent block size associated with a DM
1038: Not Collective
1040: Input Parameter:
1041: . dm - the DM with block structure
1043: Output Parameter:
1044: . bs - the block size, 1 implies no exploitable block structure
1046: Level: intermediate
1048: .seealso: ISCreateBlock(), VecSetBlockSize(), MatSetBlockSize(), DMGetLocalToGlobalMapping()
1049: @*/
1050: PetscErrorCode DMGetBlockSize(DM dm,PetscInt *bs)
1051: {
1055: if (dm->bs < 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"DM does not have enough information to provide a block size yet");
1056: *bs = dm->bs;
1057: return(0);
1058: }
1060: /*@
1061: DMCreateInterpolation - Gets interpolation matrix between two DM objects
1063: Collective on DM
1065: Input Parameter:
1066: + dm1 - the DM object
1067: - dm2 - the second, finer DM object
1069: Output Parameter:
1070: + mat - the interpolation
1071: - vec - the scaling (optional)
1073: Level: developer
1075: Notes: For DMDA objects this only works for "uniform refinement", that is the refined mesh was obtained DMRefine() or the coarse mesh was obtained by
1076: DMCoarsen(). The coordinates set into the DMDA are completely ignored in computing the interpolation.
1078: For DMDA objects you can use this interpolation (more precisely the interpolation from the DMGetCoordinateDM()) to interpolate the mesh coordinate vectors
1079: EXCEPT in the periodic case where it does not make sense since the coordinate vectors are not periodic.
1082: .seealso DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateColoring(), DMCreateMatrix(), DMRefine(), DMCoarsen(), DMCreateRestriction()
1084: @*/
1085: PetscErrorCode DMCreateInterpolation(DM dm1,DM dm2,Mat *mat,Vec *vec)
1086: {
1092: PetscLogEventBegin(DM_CreateInterpolation,dm1,dm2,0,0);
1093: (*dm1->ops->createinterpolation)(dm1,dm2,mat,vec);
1094: PetscLogEventEnd(DM_CreateInterpolation,dm1,dm2,0,0);
1095: return(0);
1096: }
1098: /*@
1099: DMCreateRestriction - Gets restriction matrix between two DM objects
1101: Collective on DM
1103: Input Parameter:
1104: + dm1 - the DM object
1105: - dm2 - the second, finer DM object
1107: Output Parameter:
1108: . mat - the restriction
1111: Level: developer
1113: Notes: For DMDA objects this only works for "uniform refinement", that is the refined mesh was obtained DMRefine() or the coarse mesh was obtained by
1114: DMCoarsen(). The coordinates set into the DMDA are completely ignored in computing the interpolation.
1117: .seealso DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateColoring(), DMCreateMatrix(), DMRefine(), DMCoarsen(), DMCreateInterpolation()
1119: @*/
1120: PetscErrorCode DMCreateRestriction(DM dm1,DM dm2,Mat *mat)
1121: {
1127: PetscLogEventBegin(DM_CreateRestriction,dm1,dm2,0,0);
1128: if (!dm1->ops->createrestriction) SETERRQ(PetscObjectComm((PetscObject)dm1),PETSC_ERR_SUP,"DMCreateRestriction not implemented for this type");
1129: (*dm1->ops->createrestriction)(dm1,dm2,mat);
1130: PetscLogEventEnd(DM_CreateRestriction,dm1,dm2,0,0);
1131: return(0);
1132: }
1134: /*@
1135: DMCreateInjection - Gets injection matrix between two DM objects
1137: Collective on DM
1139: Input Parameter:
1140: + dm1 - the DM object
1141: - dm2 - the second, finer DM object
1143: Output Parameter:
1144: . mat - the injection
1146: Level: developer
1148: Notes: For DMDA objects this only works for "uniform refinement", that is the refined mesh was obtained DMRefine() or the coarse mesh was obtained by
1149: DMCoarsen(). The coordinates set into the DMDA are completely ignored in computing the injection.
1151: .seealso DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateColoring(), DMCreateMatrix(), DMCreateInterpolation()
1153: @*/
1154: PetscErrorCode DMCreateInjection(DM dm1,DM dm2,Mat *mat)
1155: {
1161: if (!dm1->ops->getinjection) SETERRQ(PetscObjectComm((PetscObject)dm1),PETSC_ERR_SUP,"DMCreateInjection not implemented for this type");
1162: (*dm1->ops->getinjection)(dm1,dm2,mat);
1163: return(0);
1164: }
1166: /*@
1167: DMCreateMassMatrix - Gets mass matrix between two DM objects, M_ij = \int \phi_i \psi_j
1169: Collective on DM
1171: Input Parameter:
1172: + dm1 - the DM object
1173: - dm2 - the second, finer DM object
1175: Output Parameter:
1176: . mat - the interpolation
1178: Level: developer
1180: .seealso DMCreateMatrix(), DMRefine(), DMCoarsen(), DMCreateRestriction(), DMCreateInterpolation(), DMCreateInjection()
1181: @*/
1182: PetscErrorCode DMCreateMassMatrix(DM dm1, DM dm2, Mat *mat)
1183: {
1189: (*dm1->ops->createmassmatrix)(dm1, dm2, mat);
1190: return(0);
1191: }
1193: /*@
1194: DMCreateColoring - Gets coloring for a DM
1196: Collective on DM
1198: Input Parameter:
1199: + dm - the DM object
1200: - ctype - IS_COLORING_LOCAL or IS_COLORING_GLOBAL
1202: Output Parameter:
1203: . coloring - the coloring
1205: Level: developer
1207: .seealso DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateMatrix(), DMSetMatType()
1209: @*/
1210: PetscErrorCode DMCreateColoring(DM dm,ISColoringType ctype,ISColoring *coloring)
1211: {
1216: if (!dm->ops->getcoloring) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"No coloring for this type of DM yet");
1217: (*dm->ops->getcoloring)(dm,ctype,coloring);
1218: return(0);
1219: }
1221: /*@
1222: DMCreateMatrix - Gets empty Jacobian for a DM
1224: Collective on DM
1226: Input Parameter:
1227: . dm - the DM object
1229: Output Parameter:
1230: . mat - the empty Jacobian
1232: Level: beginner
1234: Notes: This properly preallocates the number of nonzeros in the sparse matrix so you
1235: do not need to do it yourself.
1237: By default it also sets the nonzero structure and puts in the zero entries. To prevent setting
1238: the nonzero pattern call DMSetMatrixPreallocateOnly()
1240: For structured grid problems, when you call MatView() on this matrix it is displayed using the global natural ordering, NOT in the ordering used
1241: internally by PETSc.
1243: For structured grid problems, in general it is easiest to use MatSetValuesStencil() or MatSetValuesLocal() to put values into the matrix because MatSetValues() requires
1244: the indices for the global numbering for DMDAs which is complicated.
1246: .seealso DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMSetMatType()
1248: @*/
1249: PetscErrorCode DMCreateMatrix(DM dm,Mat *mat)
1250: {
1255: MatInitializePackage();
1258: (*dm->ops->creatematrix)(dm,mat);
1259: return(0);
1260: }
1262: /*@
1263: DMSetMatrixPreallocateOnly - When DMCreateMatrix() is called the matrix will be properly
1264: preallocated but the nonzero structure and zero values will not be set.
1266: Logically Collective on DM
1268: Input Parameter:
1269: + dm - the DM
1270: - only - PETSC_TRUE if only want preallocation
1272: Level: developer
1273: .seealso DMCreateMatrix(), DMSetMatrixStructureOnly()
1274: @*/
1275: PetscErrorCode DMSetMatrixPreallocateOnly(DM dm, PetscBool only)
1276: {
1279: dm->prealloc_only = only;
1280: return(0);
1281: }
1283: /*@
1284: DMSetMatrixStructureOnly - When DMCreateMatrix() is called, the matrix structure will be created
1285: but the array for values will not be allocated.
1287: Logically Collective on DM
1289: Input Parameter:
1290: + dm - the DM
1291: - only - PETSC_TRUE if only want matrix stucture
1293: Level: developer
1294: .seealso DMCreateMatrix(), DMSetMatrixPreallocateOnly()
1295: @*/
1296: PetscErrorCode DMSetMatrixStructureOnly(DM dm, PetscBool only)
1297: {
1300: dm->structure_only = only;
1301: return(0);
1302: }
1304: /*@C
1305: DMGetWorkArray - Gets a work array guaranteed to be at least the input size, restore with DMRestoreWorkArray()
1307: Not Collective
1309: Input Parameters:
1310: + dm - the DM object
1311: . count - The minium size
1312: - dtype - MPI data type, often MPIU_REAL, MPIU_SCALAR, MPIU_INT)
1314: Output Parameter:
1315: . array - the work array
1317: Level: developer
1319: .seealso DMDestroy(), DMCreate()
1320: @*/
1321: PetscErrorCode DMGetWorkArray(DM dm,PetscInt count,MPI_Datatype dtype,void *mem)
1322: {
1324: DMWorkLink link;
1325: PetscMPIInt dsize;
1330: if (dm->workin) {
1331: link = dm->workin;
1332: dm->workin = dm->workin->next;
1333: } else {
1334: PetscNewLog(dm,&link);
1335: }
1336: MPI_Type_size(dtype,&dsize);
1337: if (((size_t)dsize*count) > link->bytes) {
1338: PetscFree(link->mem);
1339: PetscMalloc(dsize*count,&link->mem);
1340: link->bytes = dsize*count;
1341: }
1342: link->next = dm->workout;
1343: dm->workout = link;
1344: *(void**)mem = link->mem;
1345: return(0);
1346: }
1348: /*@C
1349: DMRestoreWorkArray - Restores a work array guaranteed to be at least the input size, restore with DMRestoreWorkArray()
1351: Not Collective
1353: Input Parameters:
1354: + dm - the DM object
1355: . count - The minium size
1356: - dtype - MPI data type, often MPIU_REAL, MPIU_SCALAR, MPIU_INT
1358: Output Parameter:
1359: . array - the work array
1361: Level: developer
1363: Developer Notes: count and dtype are ignored, they are only needed for DMGetWorkArray()
1364: .seealso DMDestroy(), DMCreate()
1365: @*/
1366: PetscErrorCode DMRestoreWorkArray(DM dm,PetscInt count,MPI_Datatype dtype,void *mem)
1367: {
1368: DMWorkLink *p,link;
1373: for (p=&dm->workout; (link=*p); p=&link->next) {
1374: if (link->mem == *(void**)mem) {
1375: *p = link->next;
1376: link->next = dm->workin;
1377: dm->workin = link;
1378: *(void**)mem = NULL;
1379: return(0);
1380: }
1381: }
1382: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Array was not checked out");
1383: }
1385: PetscErrorCode DMSetNullSpaceConstructor(DM dm, PetscInt field, PetscErrorCode (*nullsp)(DM dm, PetscInt field, MatNullSpace *nullSpace))
1386: {
1389: if (field >= 10) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle %d >= 10 fields", field);
1390: dm->nullspaceConstructors[field] = nullsp;
1391: return(0);
1392: }
1394: /*@C
1395: DMCreateFieldIS - Creates a set of IS objects with the global indices of dofs for each field
1397: Not collective
1399: Input Parameter:
1400: . dm - the DM object
1402: Output Parameters:
1403: + numFields - The number of fields (or NULL if not requested)
1404: . fieldNames - The name for each field (or NULL if not requested)
1405: - fields - The global indices for each field (or NULL if not requested)
1407: Level: intermediate
1409: Notes:
1410: The user is responsible for freeing all requested arrays. In particular, every entry of names should be freed with
1411: PetscFree(), every entry of fields should be destroyed with ISDestroy(), and both arrays should be freed with
1412: PetscFree().
1414: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix()
1415: @*/
1416: PetscErrorCode DMCreateFieldIS(DM dm, PetscInt *numFields, char ***fieldNames, IS **fields)
1417: {
1418: PetscSection section, sectionGlobal;
1423: if (numFields) {
1425: *numFields = 0;
1426: }
1427: if (fieldNames) {
1429: *fieldNames = NULL;
1430: }
1431: if (fields) {
1433: *fields = NULL;
1434: }
1435: DMGetDefaultSection(dm, §ion);
1436: if (section) {
1437: PetscInt *fieldSizes, **fieldIndices;
1438: PetscInt nF, f, pStart, pEnd, p;
1440: DMGetDefaultGlobalSection(dm, §ionGlobal);
1441: PetscSectionGetNumFields(section, &nF);
1442: PetscMalloc2(nF,&fieldSizes,nF,&fieldIndices);
1443: PetscSectionGetChart(sectionGlobal, &pStart, &pEnd);
1444: for (f = 0; f < nF; ++f) {
1445: fieldSizes[f] = 0;
1446: }
1447: for (p = pStart; p < pEnd; ++p) {
1448: PetscInt gdof;
1450: PetscSectionGetDof(sectionGlobal, p, &gdof);
1451: if (gdof > 0) {
1452: for (f = 0; f < nF; ++f) {
1453: PetscInt fdof, fcdof;
1455: PetscSectionGetFieldDof(section, p, f, &fdof);
1456: PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);
1457: fieldSizes[f] += fdof-fcdof;
1458: }
1459: }
1460: }
1461: for (f = 0; f < nF; ++f) {
1462: PetscMalloc1(fieldSizes[f], &fieldIndices[f]);
1463: fieldSizes[f] = 0;
1464: }
1465: for (p = pStart; p < pEnd; ++p) {
1466: PetscInt gdof, goff;
1468: PetscSectionGetDof(sectionGlobal, p, &gdof);
1469: if (gdof > 0) {
1470: PetscSectionGetOffset(sectionGlobal, p, &goff);
1471: for (f = 0; f < nF; ++f) {
1472: PetscInt fdof, fcdof, fc;
1474: PetscSectionGetFieldDof(section, p, f, &fdof);
1475: PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);
1476: for (fc = 0; fc < fdof-fcdof; ++fc, ++fieldSizes[f]) {
1477: fieldIndices[f][fieldSizes[f]] = goff++;
1478: }
1479: }
1480: }
1481: }
1482: if (numFields) *numFields = nF;
1483: if (fieldNames) {
1484: PetscMalloc1(nF, fieldNames);
1485: for (f = 0; f < nF; ++f) {
1486: const char *fieldName;
1488: PetscSectionGetFieldName(section, f, &fieldName);
1489: PetscStrallocpy(fieldName, (char**) &(*fieldNames)[f]);
1490: }
1491: }
1492: if (fields) {
1493: PetscMalloc1(nF, fields);
1494: for (f = 0; f < nF; ++f) {
1495: ISCreateGeneral(PetscObjectComm((PetscObject)dm), fieldSizes[f], fieldIndices[f], PETSC_OWN_POINTER, &(*fields)[f]);
1496: }
1497: }
1498: PetscFree2(fieldSizes,fieldIndices);
1499: } else if (dm->ops->createfieldis) {
1500: (*dm->ops->createfieldis)(dm, numFields, fieldNames, fields);
1501: }
1502: return(0);
1503: }
1506: /*@C
1507: DMCreateFieldDecomposition - Returns a list of IS objects defining a decomposition of a problem into subproblems
1508: corresponding to different fields: each IS contains the global indices of the dofs of the
1509: corresponding field. The optional list of DMs define the DM for each subproblem.
1510: Generalizes DMCreateFieldIS().
1512: Not collective
1514: Input Parameter:
1515: . dm - the DM object
1517: Output Parameters:
1518: + len - The number of subproblems in the field decomposition (or NULL if not requested)
1519: . namelist - The name for each field (or NULL if not requested)
1520: . islist - The global indices for each field (or NULL if not requested)
1521: - dmlist - The DMs for each field subproblem (or NULL, if not requested; if NULL is returned, no DMs are defined)
1523: Level: intermediate
1525: Notes:
1526: The user is responsible for freeing all requested arrays. In particular, every entry of names should be freed with
1527: PetscFree(), every entry of is should be destroyed with ISDestroy(), every entry of dm should be destroyed with DMDestroy(),
1528: and all of the arrays should be freed with PetscFree().
1530: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateFieldIS()
1531: @*/
1532: PetscErrorCode DMCreateFieldDecomposition(DM dm, PetscInt *len, char ***namelist, IS **islist, DM **dmlist)
1533: {
1538: if (len) {
1540: *len = 0;
1541: }
1542: if (namelist) {
1544: *namelist = 0;
1545: }
1546: if (islist) {
1548: *islist = 0;
1549: }
1550: if (dmlist) {
1552: *dmlist = 0;
1553: }
1554: /*
1555: Is it a good idea to apply the following check across all impls?
1556: Perhaps some impls can have a well-defined decomposition before DMSetUp?
1557: This, however, follows the general principle that accessors are not well-behaved until the object is set up.
1558: */
1559: if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE, "Decomposition defined only after DMSetUp");
1560: if (!dm->ops->createfielddecomposition) {
1561: PetscSection section;
1562: PetscInt numFields, f;
1564: DMGetDefaultSection(dm, §ion);
1565: if (section) {PetscSectionGetNumFields(section, &numFields);}
1566: if (section && numFields && dm->ops->createsubdm) {
1567: if (len) *len = numFields;
1568: if (namelist) {PetscMalloc1(numFields,namelist);}
1569: if (islist) {PetscMalloc1(numFields,islist);}
1570: if (dmlist) {PetscMalloc1(numFields,dmlist);}
1571: for (f = 0; f < numFields; ++f) {
1572: const char *fieldName;
1574: DMCreateSubDM(dm, 1, &f, islist ? &(*islist)[f] : NULL, dmlist ? &(*dmlist)[f] : NULL);
1575: if (namelist) {
1576: PetscSectionGetFieldName(section, f, &fieldName);
1577: PetscStrallocpy(fieldName, (char**) &(*namelist)[f]);
1578: }
1579: }
1580: } else {
1581: DMCreateFieldIS(dm, len, namelist, islist);
1582: /* By default there are no DMs associated with subproblems. */
1583: if (dmlist) *dmlist = NULL;
1584: }
1585: } else {
1586: (*dm->ops->createfielddecomposition)(dm,len,namelist,islist,dmlist);
1587: }
1588: return(0);
1589: }
1591: /*@
1592: DMCreateSubDM - Returns an IS and DM encapsulating a subproblem defined by the fields passed in.
1593: The fields are defined by DMCreateFieldIS().
1595: Not collective
1597: Input Parameters:
1598: + dm - The DM object
1599: . numFields - The number of fields in this subproblem
1600: - fields - The field numbers of the selected fields
1602: Output Parameters:
1603: + is - The global indices for the subproblem
1604: - subdm - The DM for the subproblem
1606: Level: intermediate
1608: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateFieldIS()
1609: @*/
1610: PetscErrorCode DMCreateSubDM(DM dm, PetscInt numFields, const PetscInt fields[], IS *is, DM *subdm)
1611: {
1619: if (dm->ops->createsubdm) {
1620: (*dm->ops->createsubdm)(dm, numFields, fields, is, subdm);
1621: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "This type has no DMCreateSubDM implementation defined");
1622: return(0);
1623: }
1625: /*@C
1626: DMCreateSuperDM - Returns an arrays of ISes and DM encapsulating a superproblem defined by the DMs passed in.
1628: Not collective
1630: Input Parameter:
1631: + dms - The DM objects
1632: - len - The number of DMs
1634: Output Parameters:
1635: + is - The global indices for the subproblem
1636: - superdm - The DM for the superproblem
1638: Level: intermediate
1640: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateFieldIS()
1641: @*/
1642: PetscErrorCode DMCreateSuperDM(DM dms[], PetscInt len, IS **is, DM *superdm)
1643: {
1644: PetscInt i;
1652: if (len < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of DMs must be nonnegative: %D", len);
1653: if (len) {
1654: if (dms[0]->ops->createsuperdm) {
1655: (*dms[0]->ops->createsuperdm)(dms, len, is, superdm);
1656: } else SETERRQ(PetscObjectComm((PetscObject) dms[0]), PETSC_ERR_SUP, "This type has no DMCreateSuperDM implementation defined");
1657: }
1658: return(0);
1659: }
1662: /*@C
1663: DMCreateDomainDecomposition - Returns lists of IS objects defining a decomposition of a problem into subproblems
1664: corresponding to restrictions to pairs nested subdomains: each IS contains the global
1665: indices of the dofs of the corresponding subdomains. The inner subdomains conceptually
1666: define a nonoverlapping covering, while outer subdomains can overlap.
1667: The optional list of DMs define the DM for each subproblem.
1669: Not collective
1671: Input Parameter:
1672: . dm - the DM object
1674: Output Parameters:
1675: + len - The number of subproblems in the domain decomposition (or NULL if not requested)
1676: . namelist - The name for each subdomain (or NULL if not requested)
1677: . innerislist - The global indices for each inner subdomain (or NULL, if not requested)
1678: . outerislist - The global indices for each outer subdomain (or NULL, if not requested)
1679: - dmlist - The DMs for each subdomain subproblem (or NULL, if not requested; if NULL is returned, no DMs are defined)
1681: Level: intermediate
1683: Notes:
1684: The user is responsible for freeing all requested arrays. In particular, every entry of names should be freed with
1685: PetscFree(), every entry of is should be destroyed with ISDestroy(), every entry of dm should be destroyed with DMDestroy(),
1686: and all of the arrays should be freed with PetscFree().
1688: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateDomainDecompositionDM(), DMCreateFieldDecomposition()
1689: @*/
1690: PetscErrorCode DMCreateDomainDecomposition(DM dm, PetscInt *len, char ***namelist, IS **innerislist, IS **outerislist, DM **dmlist)
1691: {
1692: PetscErrorCode ierr;
1693: DMSubDomainHookLink link;
1694: PetscInt i,l;
1703: /*
1704: Is it a good idea to apply the following check across all impls?
1705: Perhaps some impls can have a well-defined decomposition before DMSetUp?
1706: This, however, follows the general principle that accessors are not well-behaved until the object is set up.
1707: */
1708: if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE, "Decomposition defined only after DMSetUp");
1709: if (dm->ops->createdomaindecomposition) {
1710: (*dm->ops->createdomaindecomposition)(dm,&l,namelist,innerislist,outerislist,dmlist);
1711: /* copy subdomain hooks and context over to the subdomain DMs */
1712: if (dmlist && *dmlist) {
1713: for (i = 0; i < l; i++) {
1714: for (link=dm->subdomainhook; link; link=link->next) {
1715: if (link->ddhook) {(*link->ddhook)(dm,(*dmlist)[i],link->ctx);}
1716: }
1717: if (dm->ctx) (*dmlist)[i]->ctx = dm->ctx;
1718: }
1719: }
1720: if (len) *len = l;
1721: }
1722: return(0);
1723: }
1726: /*@C
1727: DMCreateDomainDecompositionScatters - Returns scatters to the subdomain vectors from the global vector
1729: Not collective
1731: Input Parameters:
1732: + dm - the DM object
1733: . n - the number of subdomain scatters
1734: - subdms - the local subdomains
1736: Output Parameters:
1737: + n - the number of scatters returned
1738: . iscat - scatter from global vector to nonoverlapping global vector entries on subdomain
1739: . oscat - scatter from global vector to overlapping global vector entries on subdomain
1740: - gscat - scatter from global vector to local vector on subdomain (fills in ghosts)
1742: Notes: This is an alternative to the iis and ois arguments in DMCreateDomainDecomposition that allow for the solution
1743: of general nonlinear problems with overlapping subdomain methods. While merely having index sets that enable subsets
1744: of the residual equations to be created is fine for linear problems, nonlinear problems require local assembly of
1745: solution and residual data.
1747: Level: developer
1749: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateFieldIS()
1750: @*/
1751: PetscErrorCode DMCreateDomainDecompositionScatters(DM dm,PetscInt n,DM *subdms,VecScatter **iscat,VecScatter **oscat,VecScatter **gscat)
1752: {
1758: if (dm->ops->createddscatters) {
1759: (*dm->ops->createddscatters)(dm,n,subdms,iscat,oscat,gscat);
1760: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "This type has no DMCreateDomainDecompositionScatter implementation defined");
1761: return(0);
1762: }
1764: /*@
1765: DMRefine - Refines a DM object
1767: Collective on DM
1769: Input Parameter:
1770: + dm - the DM object
1771: - comm - the communicator to contain the new DM object (or MPI_COMM_NULL)
1773: Output Parameter:
1774: . dmf - the refined DM, or NULL
1776: Note: If no refinement was done, the return value is NULL
1778: Level: developer
1780: .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation()
1781: @*/
1782: PetscErrorCode DMRefine(DM dm,MPI_Comm comm,DM *dmf)
1783: {
1784: PetscErrorCode ierr;
1785: DMRefineHookLink link;
1789: PetscLogEventBegin(DM_Refine,dm,0,0,0);
1790: if (!dm->ops->refine) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"This DM cannot refine");
1791: (*dm->ops->refine)(dm,comm,dmf);
1792: if (*dmf) {
1793: (*dmf)->ops->creatematrix = dm->ops->creatematrix;
1795: PetscObjectCopyFortranFunctionPointers((PetscObject)dm,(PetscObject)*dmf);
1797: (*dmf)->ctx = dm->ctx;
1798: (*dmf)->leveldown = dm->leveldown;
1799: (*dmf)->levelup = dm->levelup + 1;
1801: DMSetMatType(*dmf,dm->mattype);
1802: for (link=dm->refinehook; link; link=link->next) {
1803: if (link->refinehook) {
1804: (*link->refinehook)(dm,*dmf,link->ctx);
1805: }
1806: }
1807: }
1808: PetscLogEventEnd(DM_Refine,dm,0,0,0);
1809: return(0);
1810: }
1812: /*@C
1813: DMRefineHookAdd - adds a callback to be run when interpolating a nonlinear problem to a finer grid
1815: Logically Collective
1817: Input Arguments:
1818: + coarse - nonlinear solver context on which to run a hook when restricting to a coarser level
1819: . refinehook - function to run when setting up a coarser level
1820: . interphook - function to run to update data on finer levels (once per SNESSolve())
1821: - ctx - [optional] user-defined context for provide data for the hooks (may be NULL)
1823: Calling sequence of refinehook:
1824: $ refinehook(DM coarse,DM fine,void *ctx);
1826: + coarse - coarse level DM
1827: . fine - fine level DM to interpolate problem to
1828: - ctx - optional user-defined function context
1830: Calling sequence for interphook:
1831: $ interphook(DM coarse,Mat interp,DM fine,void *ctx)
1833: + coarse - coarse level DM
1834: . interp - matrix interpolating a coarse-level solution to the finer grid
1835: . fine - fine level DM to update
1836: - ctx - optional user-defined function context
1838: Level: advanced
1840: Notes:
1841: This function is only needed if auxiliary data needs to be passed to fine grids while grid sequencing
1843: If this function is called multiple times, the hooks will be run in the order they are added.
1845: This function is currently not available from Fortran.
1847: .seealso: DMCoarsenHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
1848: @*/
1849: PetscErrorCode DMRefineHookAdd(DM coarse,PetscErrorCode (*refinehook)(DM,DM,void*),PetscErrorCode (*interphook)(DM,Mat,DM,void*),void *ctx)
1850: {
1851: PetscErrorCode ierr;
1852: DMRefineHookLink link,*p;
1856: for (p=&coarse->refinehook; *p; p=&(*p)->next) { /* Scan to the end of the current list of hooks */
1857: if ((*p)->refinehook == refinehook && (*p)->interphook == interphook && (*p)->ctx == ctx) return(0);
1858: }
1859: PetscNew(&link);
1860: link->refinehook = refinehook;
1861: link->interphook = interphook;
1862: link->ctx = ctx;
1863: link->next = NULL;
1864: *p = link;
1865: return(0);
1866: }
1868: /*@C
1869: DMRefineHookRemove - remove a callback from the list of hooks to be run when interpolating a nonlinear problem to a finer grid
1871: Logically Collective
1873: Input Arguments:
1874: + coarse - nonlinear solver context on which to run a hook when restricting to a coarser level
1875: . refinehook - function to run when setting up a coarser level
1876: . interphook - function to run to update data on finer levels (once per SNESSolve())
1877: - ctx - [optional] user-defined context for provide data for the hooks (may be NULL)
1879: Level: advanced
1881: Notes:
1882: This function does nothing if the hook is not in the list.
1884: This function is currently not available from Fortran.
1886: .seealso: DMCoarsenHookRemove(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
1887: @*/
1888: PetscErrorCode DMRefineHookRemove(DM coarse,PetscErrorCode (*refinehook)(DM,DM,void*),PetscErrorCode (*interphook)(DM,Mat,DM,void*),void *ctx)
1889: {
1890: PetscErrorCode ierr;
1891: DMRefineHookLink link,*p;
1895: for (p=&coarse->refinehook; *p; p=&(*p)->next) { /* Search the list of current hooks */
1896: if ((*p)->refinehook == refinehook && (*p)->interphook == interphook && (*p)->ctx == ctx) {
1897: link = *p;
1898: *p = link->next;
1899: PetscFree(link);
1900: break;
1901: }
1902: }
1903: return(0);
1904: }
1906: /*@
1907: DMInterpolate - interpolates user-defined problem data to a finer DM by running hooks registered by DMRefineHookAdd()
1909: Collective if any hooks are
1911: Input Arguments:
1912: + coarse - coarser DM to use as a base
1913: . interp - interpolation matrix, apply using MatInterpolate()
1914: - fine - finer DM to update
1916: Level: developer
1918: .seealso: DMRefineHookAdd(), MatInterpolate()
1919: @*/
1920: PetscErrorCode DMInterpolate(DM coarse,Mat interp,DM fine)
1921: {
1922: PetscErrorCode ierr;
1923: DMRefineHookLink link;
1926: for (link=fine->refinehook; link; link=link->next) {
1927: if (link->interphook) {
1928: (*link->interphook)(coarse,interp,fine,link->ctx);
1929: }
1930: }
1931: return(0);
1932: }
1934: /*@
1935: DMGetRefineLevel - Get's the number of refinements that have generated this DM.
1937: Not Collective
1939: Input Parameter:
1940: . dm - the DM object
1942: Output Parameter:
1943: . level - number of refinements
1945: Level: developer
1947: .seealso DMCoarsen(), DMGetCoarsenLevel(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation()
1949: @*/
1950: PetscErrorCode DMGetRefineLevel(DM dm,PetscInt *level)
1951: {
1954: *level = dm->levelup;
1955: return(0);
1956: }
1958: /*@
1959: DMSetRefineLevel - Set's the number of refinements that have generated this DM.
1961: Not Collective
1963: Input Parameter:
1964: + dm - the DM object
1965: - level - number of refinements
1967: Level: advanced
1969: Notes: This value is used by PCMG to determine how many multigrid levels to use
1971: .seealso DMCoarsen(), DMGetCoarsenLevel(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation()
1973: @*/
1974: PetscErrorCode DMSetRefineLevel(DM dm,PetscInt level)
1975: {
1978: dm->levelup = level;
1979: return(0);
1980: }
1982: /*@C
1983: DMGlobalToLocalHookAdd - adds a callback to be run when global to local is called
1985: Logically Collective
1987: Input Arguments:
1988: + dm - the DM
1989: . beginhook - function to run at the beginning of DMGlobalToLocalBegin()
1990: . endhook - function to run after DMGlobalToLocalEnd() has completed
1991: - ctx - [optional] user-defined context for provide data for the hooks (may be NULL)
1993: Calling sequence for beginhook:
1994: $ beginhook(DM fine,VecScatter out,VecScatter in,DM coarse,void *ctx)
1996: + dm - global DM
1997: . g - global vector
1998: . mode - mode
1999: . l - local vector
2000: - ctx - optional user-defined function context
2003: Calling sequence for endhook:
2004: $ endhook(DM fine,VecScatter out,VecScatter in,DM coarse,void *ctx)
2006: + global - global DM
2007: - ctx - optional user-defined function context
2009: Level: advanced
2011: .seealso: DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2012: @*/
2013: PetscErrorCode DMGlobalToLocalHookAdd(DM dm,PetscErrorCode (*beginhook)(DM,Vec,InsertMode,Vec,void*),PetscErrorCode (*endhook)(DM,Vec,InsertMode,Vec,void*),void *ctx)
2014: {
2015: PetscErrorCode ierr;
2016: DMGlobalToLocalHookLink link,*p;
2020: for (p=&dm->gtolhook; *p; p=&(*p)->next) {} /* Scan to the end of the current list of hooks */
2021: PetscNew(&link);
2022: link->beginhook = beginhook;
2023: link->endhook = endhook;
2024: link->ctx = ctx;
2025: link->next = NULL;
2026: *p = link;
2027: return(0);
2028: }
2030: static PetscErrorCode DMGlobalToLocalHook_Constraints(DM dm, Vec g, InsertMode mode, Vec l, void *ctx)
2031: {
2032: Mat cMat;
2033: Vec cVec;
2034: PetscSection section, cSec;
2035: PetscInt pStart, pEnd, p, dof;
2040: DMGetDefaultConstraints(dm,&cSec,&cMat);
2041: if (cMat && (mode == INSERT_VALUES || mode == INSERT_ALL_VALUES || mode == INSERT_BC_VALUES)) {
2042: PetscInt nRows;
2044: MatGetSize(cMat,&nRows,NULL);
2045: if (nRows <= 0) return(0);
2046: DMGetDefaultSection(dm,§ion);
2047: MatCreateVecs(cMat,NULL,&cVec);
2048: MatMult(cMat,l,cVec);
2049: PetscSectionGetChart(cSec,&pStart,&pEnd);
2050: for (p = pStart; p < pEnd; p++) {
2051: PetscSectionGetDof(cSec,p,&dof);
2052: if (dof) {
2053: PetscScalar *vals;
2054: VecGetValuesSection(cVec,cSec,p,&vals);
2055: VecSetValuesSection(l,section,p,vals,INSERT_ALL_VALUES);
2056: }
2057: }
2058: VecDestroy(&cVec);
2059: }
2060: return(0);
2061: }
2063: /*@
2064: DMGlobalToLocalBegin - Begins updating local vectors from global vector
2066: Neighbor-wise Collective on DM
2068: Input Parameters:
2069: + dm - the DM object
2070: . g - the global vector
2071: . mode - INSERT_VALUES or ADD_VALUES
2072: - l - the local vector
2075: Level: beginner
2077: .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMGlobalToLocalEnd(), DMLocalToGlobalBegin()
2079: @*/
2080: PetscErrorCode DMGlobalToLocalBegin(DM dm,Vec g,InsertMode mode,Vec l)
2081: {
2082: PetscSF sf;
2083: PetscErrorCode ierr;
2084: DMGlobalToLocalHookLink link;
2088: for (link=dm->gtolhook; link; link=link->next) {
2089: if (link->beginhook) {
2090: (*link->beginhook)(dm,g,mode,l,link->ctx);
2091: }
2092: }
2093: DMGetDefaultSF(dm, &sf);
2094: if (sf) {
2095: const PetscScalar *gArray;
2096: PetscScalar *lArray;
2098: if (mode == ADD_VALUES) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode);
2099: VecGetArray(l, &lArray);
2100: VecGetArrayRead(g, &gArray);
2101: PetscSFBcastBegin(sf, MPIU_SCALAR, gArray, lArray);
2102: VecRestoreArray(l, &lArray);
2103: VecRestoreArrayRead(g, &gArray);
2104: } else {
2105: (*dm->ops->globaltolocalbegin)(dm,g,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),l);
2106: }
2107: return(0);
2108: }
2110: /*@
2111: DMGlobalToLocalEnd - Ends updating local vectors from global vector
2113: Neighbor-wise Collective on DM
2115: Input Parameters:
2116: + dm - the DM object
2117: . g - the global vector
2118: . mode - INSERT_VALUES or ADD_VALUES
2119: - l - the local vector
2122: Level: beginner
2124: .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMGlobalToLocalEnd(), DMLocalToGlobalBegin()
2126: @*/
2127: PetscErrorCode DMGlobalToLocalEnd(DM dm,Vec g,InsertMode mode,Vec l)
2128: {
2129: PetscSF sf;
2130: PetscErrorCode ierr;
2131: const PetscScalar *gArray;
2132: PetscScalar *lArray;
2133: DMGlobalToLocalHookLink link;
2137: DMGetDefaultSF(dm, &sf);
2138: if (sf) {
2139: if (mode == ADD_VALUES) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode);
2141: VecGetArray(l, &lArray);
2142: VecGetArrayRead(g, &gArray);
2143: PetscSFBcastEnd(sf, MPIU_SCALAR, gArray, lArray);
2144: VecRestoreArray(l, &lArray);
2145: VecRestoreArrayRead(g, &gArray);
2146: } else {
2147: (*dm->ops->globaltolocalend)(dm,g,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),l);
2148: }
2149: DMGlobalToLocalHook_Constraints(dm,g,mode,l,NULL);
2150: for (link=dm->gtolhook; link; link=link->next) {
2151: if (link->endhook) {(*link->endhook)(dm,g,mode,l,link->ctx);}
2152: }
2153: return(0);
2154: }
2156: /*@C
2157: DMLocalToGlobalHookAdd - adds a callback to be run when a local to global is called
2159: Logically Collective
2161: Input Arguments:
2162: + dm - the DM
2163: . beginhook - function to run at the beginning of DMLocalToGlobalBegin()
2164: . endhook - function to run after DMLocalToGlobalEnd() has completed
2165: - ctx - [optional] user-defined context for provide data for the hooks (may be NULL)
2167: Calling sequence for beginhook:
2168: $ beginhook(DM fine,Vec l,InsertMode mode,Vec g,void *ctx)
2170: + dm - global DM
2171: . l - local vector
2172: . mode - mode
2173: . g - global vector
2174: - ctx - optional user-defined function context
2177: Calling sequence for endhook:
2178: $ endhook(DM fine,Vec l,InsertMode mode,Vec g,void *ctx)
2180: + global - global DM
2181: . l - local vector
2182: . mode - mode
2183: . g - global vector
2184: - ctx - optional user-defined function context
2186: Level: advanced
2188: .seealso: DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2189: @*/
2190: PetscErrorCode DMLocalToGlobalHookAdd(DM dm,PetscErrorCode (*beginhook)(DM,Vec,InsertMode,Vec,void*),PetscErrorCode (*endhook)(DM,Vec,InsertMode,Vec,void*),void *ctx)
2191: {
2192: PetscErrorCode ierr;
2193: DMLocalToGlobalHookLink link,*p;
2197: for (p=&dm->ltoghook; *p; p=&(*p)->next) {} /* Scan to the end of the current list of hooks */
2198: PetscNew(&link);
2199: link->beginhook = beginhook;
2200: link->endhook = endhook;
2201: link->ctx = ctx;
2202: link->next = NULL;
2203: *p = link;
2204: return(0);
2205: }
2207: static PetscErrorCode DMLocalToGlobalHook_Constraints(DM dm, Vec l, InsertMode mode, Vec g, void *ctx)
2208: {
2209: Mat cMat;
2210: Vec cVec;
2211: PetscSection section, cSec;
2212: PetscInt pStart, pEnd, p, dof;
2217: DMGetDefaultConstraints(dm,&cSec,&cMat);
2218: if (cMat && (mode == ADD_VALUES || mode == ADD_ALL_VALUES || mode == ADD_BC_VALUES)) {
2219: PetscInt nRows;
2221: MatGetSize(cMat,&nRows,NULL);
2222: if (nRows <= 0) return(0);
2223: DMGetDefaultSection(dm,§ion);
2224: MatCreateVecs(cMat,NULL,&cVec);
2225: PetscSectionGetChart(cSec,&pStart,&pEnd);
2226: for (p = pStart; p < pEnd; p++) {
2227: PetscSectionGetDof(cSec,p,&dof);
2228: if (dof) {
2229: PetscInt d;
2230: PetscScalar *vals;
2231: VecGetValuesSection(l,section,p,&vals);
2232: VecSetValuesSection(cVec,cSec,p,vals,mode);
2233: /* for this to be the true transpose, we have to zero the values that
2234: * we just extracted */
2235: for (d = 0; d < dof; d++) {
2236: vals[d] = 0.;
2237: }
2238: }
2239: }
2240: MatMultTransposeAdd(cMat,cVec,l,l);
2241: VecDestroy(&cVec);
2242: }
2243: return(0);
2244: }
2246: /*@
2247: DMLocalToGlobalBegin - updates global vectors from local vectors
2249: Neighbor-wise Collective on DM
2251: Input Parameters:
2252: + dm - the DM object
2253: . l - the local vector
2254: . 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.
2255: - g - the global vector
2257: Notes: In the ADD_VALUES case you normally would zero the receiving vector before beginning this operation.
2258: INSERT_VALUES is not supported for DMDA, in that case simply compute the values directly into a global vector instead of a local one.
2260: Level: beginner
2262: .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMGlobalToLocalEnd(), DMGlobalToLocalBegin()
2264: @*/
2265: PetscErrorCode DMLocalToGlobalBegin(DM dm,Vec l,InsertMode mode,Vec g)
2266: {
2267: PetscSF sf;
2268: PetscSection s, gs;
2269: DMLocalToGlobalHookLink link;
2270: const PetscScalar *lArray;
2271: PetscScalar *gArray;
2272: PetscBool isInsert;
2273: PetscErrorCode ierr;
2277: for (link=dm->ltoghook; link; link=link->next) {
2278: if (link->beginhook) {
2279: (*link->beginhook)(dm,l,mode,g,link->ctx);
2280: }
2281: }
2282: DMLocalToGlobalHook_Constraints(dm,l,mode,g,NULL);
2283: DMGetDefaultSF(dm, &sf);
2284: DMGetDefaultSection(dm, &s);
2285: switch (mode) {
2286: case INSERT_VALUES:
2287: case INSERT_ALL_VALUES:
2288: case INSERT_BC_VALUES:
2289: isInsert = PETSC_TRUE; break;
2290: case ADD_VALUES:
2291: case ADD_ALL_VALUES:
2292: case ADD_BC_VALUES:
2293: isInsert = PETSC_FALSE; break;
2294: default:
2295: SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode);
2296: }
2297: if (sf && !isInsert) {
2298: VecGetArrayRead(l, &lArray);
2299: VecGetArray(g, &gArray);
2300: PetscSFReduceBegin(sf, MPIU_SCALAR, lArray, gArray, MPIU_SUM);
2301: VecRestoreArrayRead(l, &lArray);
2302: VecRestoreArray(g, &gArray);
2303: } else if (s && isInsert) {
2304: PetscInt gStart, pStart, pEnd, p;
2306: DMGetDefaultGlobalSection(dm, &gs);
2307: PetscSectionGetChart(s, &pStart, &pEnd);
2308: VecGetOwnershipRange(g, &gStart, NULL);
2309: VecGetArrayRead(l, &lArray);
2310: VecGetArray(g, &gArray);
2311: for (p = pStart; p < pEnd; ++p) {
2312: PetscInt dof, gdof, cdof, gcdof, off, goff, d, e;
2314: PetscSectionGetDof(s, p, &dof);
2315: PetscSectionGetDof(gs, p, &gdof);
2316: PetscSectionGetConstraintDof(s, p, &cdof);
2317: PetscSectionGetConstraintDof(gs, p, &gcdof);
2318: PetscSectionGetOffset(s, p, &off);
2319: PetscSectionGetOffset(gs, p, &goff);
2320: /* Ignore off-process data and points with no global data */
2321: if (!gdof || goff < 0) continue;
2322: 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);
2323: /* If no constraints are enforced in the global vector */
2324: if (!gcdof) {
2325: for (d = 0; d < dof; ++d) gArray[goff-gStart+d] = lArray[off+d];
2326: /* If constraints are enforced in the global vector */
2327: } else if (cdof == gcdof) {
2328: const PetscInt *cdofs;
2329: PetscInt cind = 0;
2331: PetscSectionGetConstraintIndices(s, p, &cdofs);
2332: for (d = 0, e = 0; d < dof; ++d) {
2333: if ((cind < cdof) && (d == cdofs[cind])) {++cind; continue;}
2334: gArray[goff-gStart+e++] = lArray[off+d];
2335: }
2336: } 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);
2337: }
2338: VecRestoreArrayRead(l, &lArray);
2339: VecRestoreArray(g, &gArray);
2340: } else {
2341: (*dm->ops->localtoglobalbegin)(dm,l,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),g);
2342: }
2343: return(0);
2344: }
2346: /*@
2347: DMLocalToGlobalEnd - updates global vectors from local vectors
2349: Neighbor-wise Collective on DM
2351: Input Parameters:
2352: + dm - the DM object
2353: . l - the local vector
2354: . mode - INSERT_VALUES or ADD_VALUES
2355: - g - the global vector
2358: Level: beginner
2360: .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMGlobalToLocalEnd(), DMGlobalToLocalEnd()
2362: @*/
2363: PetscErrorCode DMLocalToGlobalEnd(DM dm,Vec l,InsertMode mode,Vec g)
2364: {
2365: PetscSF sf;
2366: PetscSection s;
2367: DMLocalToGlobalHookLink link;
2368: PetscBool isInsert;
2369: PetscErrorCode ierr;
2373: DMGetDefaultSF(dm, &sf);
2374: DMGetDefaultSection(dm, &s);
2375: switch (mode) {
2376: case INSERT_VALUES:
2377: case INSERT_ALL_VALUES:
2378: isInsert = PETSC_TRUE; break;
2379: case ADD_VALUES:
2380: case ADD_ALL_VALUES:
2381: isInsert = PETSC_FALSE; break;
2382: default:
2383: SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode);
2384: }
2385: if (sf && !isInsert) {
2386: const PetscScalar *lArray;
2387: PetscScalar *gArray;
2389: VecGetArrayRead(l, &lArray);
2390: VecGetArray(g, &gArray);
2391: PetscSFReduceEnd(sf, MPIU_SCALAR, lArray, gArray, MPIU_SUM);
2392: VecRestoreArrayRead(l, &lArray);
2393: VecRestoreArray(g, &gArray);
2394: } else if (s && isInsert) {
2395: } else {
2396: (*dm->ops->localtoglobalend)(dm,l,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),g);
2397: }
2398: for (link=dm->ltoghook; link; link=link->next) {
2399: if (link->endhook) {(*link->endhook)(dm,g,mode,l,link->ctx);}
2400: }
2401: return(0);
2402: }
2404: /*@
2405: DMLocalToLocalBegin - Maps from a local vector (including ghost points
2406: that contain irrelevant values) to another local vector where the ghost
2407: points in the second are set correctly. Must be followed by DMLocalToLocalEnd().
2409: Neighbor-wise Collective on DM and Vec
2411: Input Parameters:
2412: + dm - the DM object
2413: . g - the original local vector
2414: - mode - one of INSERT_VALUES or ADD_VALUES
2416: Output Parameter:
2417: . l - the local vector with correct ghost values
2419: Level: intermediate
2421: Notes:
2422: The local vectors used here need not be the same as those
2423: obtained from DMCreateLocalVector(), BUT they
2424: must have the same parallel data layout; they could, for example, be
2425: obtained with VecDuplicate() from the DM originating vectors.
2427: .keywords: DM, local-to-local, begin
2428: .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateLocalVector(), DMCreateGlobalVector(), DMCreateInterpolation(), DMLocalToLocalEnd(), DMGlobalToLocalEnd(), DMLocalToGlobalBegin()
2430: @*/
2431: PetscErrorCode DMLocalToLocalBegin(DM dm,Vec g,InsertMode mode,Vec l)
2432: {
2433: PetscErrorCode ierr;
2437: if (!dm->ops->localtolocalbegin) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"This DM does not support local to local maps");
2438: (*dm->ops->localtolocalbegin)(dm,g,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),l);
2439: return(0);
2440: }
2442: /*@
2443: DMLocalToLocalEnd - Maps from a local vector (including ghost points
2444: that contain irrelevant values) to another local vector where the ghost
2445: points in the second are set correctly. Must be preceded by DMLocalToLocalBegin().
2447: Neighbor-wise Collective on DM and Vec
2449: Input Parameters:
2450: + da - the DM object
2451: . g - the original local vector
2452: - mode - one of INSERT_VALUES or ADD_VALUES
2454: Output Parameter:
2455: . l - the local vector with correct ghost values
2457: Level: intermediate
2459: Notes:
2460: The local vectors used here need not be the same as those
2461: obtained from DMCreateLocalVector(), BUT they
2462: must have the same parallel data layout; they could, for example, be
2463: obtained with VecDuplicate() from the DM originating vectors.
2465: .keywords: DM, local-to-local, end
2466: .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateLocalVector(), DMCreateGlobalVector(), DMCreateInterpolation(), DMLocalToLocalBegin(), DMGlobalToLocalEnd(), DMLocalToGlobalBegin()
2468: @*/
2469: PetscErrorCode DMLocalToLocalEnd(DM dm,Vec g,InsertMode mode,Vec l)
2470: {
2471: PetscErrorCode ierr;
2475: if (!dm->ops->localtolocalend) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"This DM does not support local to local maps");
2476: (*dm->ops->localtolocalend)(dm,g,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),l);
2477: return(0);
2478: }
2481: /*@
2482: DMCoarsen - Coarsens a DM object
2484: Collective on DM
2486: Input Parameter:
2487: + dm - the DM object
2488: - comm - the communicator to contain the new DM object (or MPI_COMM_NULL)
2490: Output Parameter:
2491: . dmc - the coarsened DM
2493: Level: developer
2495: .seealso DMRefine(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation()
2497: @*/
2498: PetscErrorCode DMCoarsen(DM dm, MPI_Comm comm, DM *dmc)
2499: {
2500: PetscErrorCode ierr;
2501: DMCoarsenHookLink link;
2505: if (!dm->ops->coarsen) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"This DM cannot coarsen");
2506: PetscLogEventBegin(DM_Coarsen,dm,0,0,0);
2507: (*dm->ops->coarsen)(dm, comm, dmc);
2508: if (!(*dmc)) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "NULL coarse mesh produced");
2509: DMSetCoarseDM(dm,*dmc);
2510: (*dmc)->ops->creatematrix = dm->ops->creatematrix;
2511: PetscObjectCopyFortranFunctionPointers((PetscObject)dm,(PetscObject)*dmc);
2512: (*dmc)->ctx = dm->ctx;
2513: (*dmc)->levelup = dm->levelup;
2514: (*dmc)->leveldown = dm->leveldown + 1;
2515: DMSetMatType(*dmc,dm->mattype);
2516: for (link=dm->coarsenhook; link; link=link->next) {
2517: if (link->coarsenhook) {(*link->coarsenhook)(dm,*dmc,link->ctx);}
2518: }
2519: PetscLogEventEnd(DM_Coarsen,dm,0,0,0);
2520: return(0);
2521: }
2523: /*@C
2524: DMCoarsenHookAdd - adds a callback to be run when restricting a nonlinear problem to the coarse grid
2526: Logically Collective
2528: Input Arguments:
2529: + fine - nonlinear solver context on which to run a hook when restricting to a coarser level
2530: . coarsenhook - function to run when setting up a coarser level
2531: . restricthook - function to run to update data on coarser levels (once per SNESSolve())
2532: - ctx - [optional] user-defined context for provide data for the hooks (may be NULL)
2534: Calling sequence of coarsenhook:
2535: $ coarsenhook(DM fine,DM coarse,void *ctx);
2537: + fine - fine level DM
2538: . coarse - coarse level DM to restrict problem to
2539: - ctx - optional user-defined function context
2541: Calling sequence for restricthook:
2542: $ restricthook(DM fine,Mat mrestrict,Vec rscale,Mat inject,DM coarse,void *ctx)
2544: + fine - fine level DM
2545: . mrestrict - matrix restricting a fine-level solution to the coarse grid
2546: . rscale - scaling vector for restriction
2547: . inject - matrix restricting by injection
2548: . coarse - coarse level DM to update
2549: - ctx - optional user-defined function context
2551: Level: advanced
2553: Notes:
2554: This function is only needed if auxiliary data needs to be set up on coarse grids.
2556: If this function is called multiple times, the hooks will be run in the order they are added.
2558: In order to compose with nonlinear preconditioning without duplicating storage, the hook should be implemented to
2559: extract the finest level information from its context (instead of from the SNES).
2561: This function is currently not available from Fortran.
2563: .seealso: DMCoarsenHookRemove(), DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2564: @*/
2565: PetscErrorCode DMCoarsenHookAdd(DM fine,PetscErrorCode (*coarsenhook)(DM,DM,void*),PetscErrorCode (*restricthook)(DM,Mat,Vec,Mat,DM,void*),void *ctx)
2566: {
2567: PetscErrorCode ierr;
2568: DMCoarsenHookLink link,*p;
2572: for (p=&fine->coarsenhook; *p; p=&(*p)->next) { /* Scan to the end of the current list of hooks */
2573: if ((*p)->coarsenhook == coarsenhook && (*p)->restricthook == restricthook && (*p)->ctx == ctx) return(0);
2574: }
2575: PetscNew(&link);
2576: link->coarsenhook = coarsenhook;
2577: link->restricthook = restricthook;
2578: link->ctx = ctx;
2579: link->next = NULL;
2580: *p = link;
2581: return(0);
2582: }
2584: /*@C
2585: DMCoarsenHookRemove - remove a callback from the list of hooks to be run when restricting a nonlinear problem to the coarse grid
2587: Logically Collective
2589: Input Arguments:
2590: + fine - nonlinear solver context on which to run a hook when restricting to a coarser level
2591: . coarsenhook - function to run when setting up a coarser level
2592: . restricthook - function to run to update data on coarser levels (once per SNESSolve())
2593: - ctx - [optional] user-defined context for provide data for the hooks (may be NULL)
2595: Level: advanced
2597: Notes:
2598: This function does nothing if the hook is not in the list.
2600: This function is currently not available from Fortran.
2602: .seealso: DMCoarsenHookAdd(), DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2603: @*/
2604: PetscErrorCode DMCoarsenHookRemove(DM fine,PetscErrorCode (*coarsenhook)(DM,DM,void*),PetscErrorCode (*restricthook)(DM,Mat,Vec,Mat,DM,void*),void *ctx)
2605: {
2606: PetscErrorCode ierr;
2607: DMCoarsenHookLink link,*p;
2611: for (p=&fine->coarsenhook; *p; p=&(*p)->next) { /* Search the list of current hooks */
2612: if ((*p)->coarsenhook == coarsenhook && (*p)->restricthook == restricthook && (*p)->ctx == ctx) {
2613: link = *p;
2614: *p = link->next;
2615: PetscFree(link);
2616: break;
2617: }
2618: }
2619: return(0);
2620: }
2623: /*@
2624: DMRestrict - restricts user-defined problem data to a coarser DM by running hooks registered by DMCoarsenHookAdd()
2626: Collective if any hooks are
2628: Input Arguments:
2629: + fine - finer DM to use as a base
2630: . restrct - restriction matrix, apply using MatRestrict()
2631: . rscale - scaling vector for restriction
2632: . inject - injection matrix, also use MatRestrict()
2633: - coarse - coarser DM to update
2635: Level: developer
2637: .seealso: DMCoarsenHookAdd(), MatRestrict()
2638: @*/
2639: PetscErrorCode DMRestrict(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse)
2640: {
2641: PetscErrorCode ierr;
2642: DMCoarsenHookLink link;
2645: for (link=fine->coarsenhook; link; link=link->next) {
2646: if (link->restricthook) {
2647: (*link->restricthook)(fine,restrct,rscale,inject,coarse,link->ctx);
2648: }
2649: }
2650: return(0);
2651: }
2653: /*@C
2654: DMSubDomainHookAdd - adds a callback to be run when restricting a problem to the coarse grid
2656: Logically Collective
2658: Input Arguments:
2659: + global - global DM
2660: . ddhook - function to run to pass data to the decomposition DM upon its creation
2661: . restricthook - function to run to update data on block solve (at the beginning of the block solve)
2662: - ctx - [optional] user-defined context for provide data for the hooks (may be NULL)
2665: Calling sequence for ddhook:
2666: $ ddhook(DM global,DM block,void *ctx)
2668: + global - global DM
2669: . block - block DM
2670: - ctx - optional user-defined function context
2672: Calling sequence for restricthook:
2673: $ restricthook(DM global,VecScatter out,VecScatter in,DM block,void *ctx)
2675: + global - global DM
2676: . out - scatter to the outer (with ghost and overlap points) block vector
2677: . in - scatter to block vector values only owned locally
2678: . block - block DM
2679: - ctx - optional user-defined function context
2681: Level: advanced
2683: Notes:
2684: This function is only needed if auxiliary data needs to be set up on subdomain DMs.
2686: If this function is called multiple times, the hooks will be run in the order they are added.
2688: In order to compose with nonlinear preconditioning without duplicating storage, the hook should be implemented to
2689: extract the global information from its context (instead of from the SNES).
2691: This function is currently not available from Fortran.
2693: .seealso: DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2694: @*/
2695: PetscErrorCode DMSubDomainHookAdd(DM global,PetscErrorCode (*ddhook)(DM,DM,void*),PetscErrorCode (*restricthook)(DM,VecScatter,VecScatter,DM,void*),void *ctx)
2696: {
2697: PetscErrorCode ierr;
2698: DMSubDomainHookLink link,*p;
2702: for (p=&global->subdomainhook; *p; p=&(*p)->next) { /* Scan to the end of the current list of hooks */
2703: if ((*p)->ddhook == ddhook && (*p)->restricthook == restricthook && (*p)->ctx == ctx) return(0);
2704: }
2705: PetscNew(&link);
2706: link->restricthook = restricthook;
2707: link->ddhook = ddhook;
2708: link->ctx = ctx;
2709: link->next = NULL;
2710: *p = link;
2711: return(0);
2712: }
2714: /*@C
2715: DMSubDomainHookRemove - remove a callback from the list to be run when restricting a problem to the coarse grid
2717: Logically Collective
2719: Input Arguments:
2720: + global - global DM
2721: . ddhook - function to run to pass data to the decomposition DM upon its creation
2722: . restricthook - function to run to update data on block solve (at the beginning of the block solve)
2723: - ctx - [optional] user-defined context for provide data for the hooks (may be NULL)
2725: Level: advanced
2727: Notes:
2729: This function is currently not available from Fortran.
2731: .seealso: DMSubDomainHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2732: @*/
2733: PetscErrorCode DMSubDomainHookRemove(DM global,PetscErrorCode (*ddhook)(DM,DM,void*),PetscErrorCode (*restricthook)(DM,VecScatter,VecScatter,DM,void*),void *ctx)
2734: {
2735: PetscErrorCode ierr;
2736: DMSubDomainHookLink link,*p;
2740: for (p=&global->subdomainhook; *p; p=&(*p)->next) { /* Search the list of current hooks */
2741: if ((*p)->ddhook == ddhook && (*p)->restricthook == restricthook && (*p)->ctx == ctx) {
2742: link = *p;
2743: *p = link->next;
2744: PetscFree(link);
2745: break;
2746: }
2747: }
2748: return(0);
2749: }
2751: /*@
2752: DMSubDomainRestrict - restricts user-defined problem data to a block DM by running hooks registered by DMSubDomainHookAdd()
2754: Collective if any hooks are
2756: Input Arguments:
2757: + fine - finer DM to use as a base
2758: . oscatter - scatter from domain global vector filling subdomain global vector with overlap
2759: . gscatter - scatter from domain global vector filling subdomain local vector with ghosts
2760: - coarse - coarer DM to update
2762: Level: developer
2764: .seealso: DMCoarsenHookAdd(), MatRestrict()
2765: @*/
2766: PetscErrorCode DMSubDomainRestrict(DM global,VecScatter oscatter,VecScatter gscatter,DM subdm)
2767: {
2768: PetscErrorCode ierr;
2769: DMSubDomainHookLink link;
2772: for (link=global->subdomainhook; link; link=link->next) {
2773: if (link->restricthook) {
2774: (*link->restricthook)(global,oscatter,gscatter,subdm,link->ctx);
2775: }
2776: }
2777: return(0);
2778: }
2780: /*@
2781: DMGetCoarsenLevel - Get's the number of coarsenings that have generated this DM.
2783: Not Collective
2785: Input Parameter:
2786: . dm - the DM object
2788: Output Parameter:
2789: . level - number of coarsenings
2791: Level: developer
2793: .seealso DMCoarsen(), DMGetRefineLevel(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation()
2795: @*/
2796: PetscErrorCode DMGetCoarsenLevel(DM dm,PetscInt *level)
2797: {
2800: *level = dm->leveldown;
2801: return(0);
2802: }
2806: /*@C
2807: DMRefineHierarchy - Refines a DM object, all levels at once
2809: Collective on DM
2811: Input Parameter:
2812: + dm - the DM object
2813: - nlevels - the number of levels of refinement
2815: Output Parameter:
2816: . dmf - the refined DM hierarchy
2818: Level: developer
2820: .seealso DMCoarsenHierarchy(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation()
2822: @*/
2823: PetscErrorCode DMRefineHierarchy(DM dm,PetscInt nlevels,DM dmf[])
2824: {
2829: if (nlevels < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"nlevels cannot be negative");
2830: if (nlevels == 0) return(0);
2831: if (dm->ops->refinehierarchy) {
2832: (*dm->ops->refinehierarchy)(dm,nlevels,dmf);
2833: } else if (dm->ops->refine) {
2834: PetscInt i;
2836: DMRefine(dm,PetscObjectComm((PetscObject)dm),&dmf[0]);
2837: for (i=1; i<nlevels; i++) {
2838: DMRefine(dmf[i-1],PetscObjectComm((PetscObject)dm),&dmf[i]);
2839: }
2840: } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"No RefineHierarchy for this DM yet");
2841: return(0);
2842: }
2844: /*@C
2845: DMCoarsenHierarchy - Coarsens a DM object, all levels at once
2847: Collective on DM
2849: Input Parameter:
2850: + dm - the DM object
2851: - nlevels - the number of levels of coarsening
2853: Output Parameter:
2854: . dmc - the coarsened DM hierarchy
2856: Level: developer
2858: .seealso DMRefineHierarchy(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation()
2860: @*/
2861: PetscErrorCode DMCoarsenHierarchy(DM dm, PetscInt nlevels, DM dmc[])
2862: {
2867: if (nlevels < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"nlevels cannot be negative");
2868: if (nlevels == 0) return(0);
2870: if (dm->ops->coarsenhierarchy) {
2871: (*dm->ops->coarsenhierarchy)(dm, nlevels, dmc);
2872: } else if (dm->ops->coarsen) {
2873: PetscInt i;
2875: DMCoarsen(dm,PetscObjectComm((PetscObject)dm),&dmc[0]);
2876: for (i=1; i<nlevels; i++) {
2877: DMCoarsen(dmc[i-1],PetscObjectComm((PetscObject)dm),&dmc[i]);
2878: }
2879: } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"No CoarsenHierarchy for this DM yet");
2880: return(0);
2881: }
2883: /*@
2884: DMCreateAggregates - Gets the aggregates that map between
2885: grids associated with two DMs.
2887: Collective on DM
2889: Input Parameters:
2890: + dmc - the coarse grid DM
2891: - dmf - the fine grid DM
2893: Output Parameters:
2894: . rest - the restriction matrix (transpose of the projection matrix)
2896: Level: intermediate
2898: .keywords: interpolation, restriction, multigrid
2900: .seealso: DMRefine(), DMCreateInjection(), DMCreateInterpolation()
2901: @*/
2902: PetscErrorCode DMCreateAggregates(DM dmc, DM dmf, Mat *rest)
2903: {
2909: (*dmc->ops->getaggregates)(dmc, dmf, rest);
2910: return(0);
2911: }
2913: /*@C
2914: DMSetApplicationContextDestroy - Sets a user function that will be called to destroy the application context when the DM is destroyed
2916: Not Collective
2918: Input Parameters:
2919: + dm - the DM object
2920: - destroy - the destroy function
2922: Level: intermediate
2924: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext()
2926: @*/
2927: PetscErrorCode DMSetApplicationContextDestroy(DM dm,PetscErrorCode (*destroy)(void**))
2928: {
2931: dm->ctxdestroy = destroy;
2932: return(0);
2933: }
2935: /*@
2936: DMSetApplicationContext - Set a user context into a DM object
2938: Not Collective
2940: Input Parameters:
2941: + dm - the DM object
2942: - ctx - the user context
2944: Level: intermediate
2946: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext()
2948: @*/
2949: PetscErrorCode DMSetApplicationContext(DM dm,void *ctx)
2950: {
2953: dm->ctx = ctx;
2954: return(0);
2955: }
2957: /*@
2958: DMGetApplicationContext - Gets a user context from a DM object
2960: Not Collective
2962: Input Parameter:
2963: . dm - the DM object
2965: Output Parameter:
2966: . ctx - the user context
2968: Level: intermediate
2970: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext()
2972: @*/
2973: PetscErrorCode DMGetApplicationContext(DM dm,void *ctx)
2974: {
2977: *(void**)ctx = dm->ctx;
2978: return(0);
2979: }
2981: /*@C
2982: DMSetVariableBounds - sets a function to compute the lower and upper bound vectors for SNESVI.
2984: Logically Collective on DM
2986: Input Parameter:
2987: + dm - the DM object
2988: - f - the function that computes variable bounds used by SNESVI (use NULL to cancel a previous function that was set)
2990: Level: intermediate
2992: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext(),
2993: DMSetJacobian()
2995: @*/
2996: PetscErrorCode DMSetVariableBounds(DM dm,PetscErrorCode (*f)(DM,Vec,Vec))
2997: {
2999: dm->ops->computevariablebounds = f;
3000: return(0);
3001: }
3003: /*@
3004: DMHasVariableBounds - does the DM object have a variable bounds function?
3006: Not Collective
3008: Input Parameter:
3009: . dm - the DM object to destroy
3011: Output Parameter:
3012: . flg - PETSC_TRUE if the variable bounds function exists
3014: Level: developer
3016: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext()
3018: @*/
3019: PetscErrorCode DMHasVariableBounds(DM dm,PetscBool *flg)
3020: {
3022: *flg = (dm->ops->computevariablebounds) ? PETSC_TRUE : PETSC_FALSE;
3023: return(0);
3024: }
3026: /*@C
3027: DMComputeVariableBounds - compute variable bounds used by SNESVI.
3029: Logically Collective on DM
3031: Input Parameters:
3032: . dm - the DM object
3034: Output parameters:
3035: + xl - lower bound
3036: - xu - upper bound
3038: Level: advanced
3040: Notes: This is generally not called by users. It calls the function provided by the user with DMSetVariableBounds()
3042: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext()
3044: @*/
3045: PetscErrorCode DMComputeVariableBounds(DM dm, Vec xl, Vec xu)
3046: {
3052: if (dm->ops->computevariablebounds) {
3053: (*dm->ops->computevariablebounds)(dm, xl,xu);
3054: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "This DM is incapable of computing variable bounds.");
3055: return(0);
3056: }
3058: /*@
3059: DMHasColoring - does the DM object have a method of providing a coloring?
3061: Not Collective
3063: Input Parameter:
3064: . dm - the DM object
3066: Output Parameter:
3067: . flg - PETSC_TRUE if the DM has facilities for DMCreateColoring().
3069: Level: developer
3071: .seealso DMHasFunction(), DMCreateColoring()
3073: @*/
3074: PetscErrorCode DMHasColoring(DM dm,PetscBool *flg)
3075: {
3077: *flg = (dm->ops->getcoloring) ? PETSC_TRUE : PETSC_FALSE;
3078: return(0);
3079: }
3081: /*@
3082: DMHasCreateRestriction - does the DM object have a method of providing a restriction?
3084: Not Collective
3086: Input Parameter:
3087: . dm - the DM object
3089: Output Parameter:
3090: . flg - PETSC_TRUE if the DM has facilities for DMCreateRestriction().
3092: Level: developer
3094: .seealso DMHasFunction(), DMCreateRestriction()
3096: @*/
3097: PetscErrorCode DMHasCreateRestriction(DM dm,PetscBool *flg)
3098: {
3100: *flg = (dm->ops->createrestriction) ? PETSC_TRUE : PETSC_FALSE;
3101: return(0);
3102: }
3105: /*@
3106: DMHasCreateInjection - does the DM object have a method of providing an injection?
3108: Not Collective
3110: Input Parameter:
3111: . dm - the DM object
3113: Output Parameter:
3114: . flg - PETSC_TRUE if the DM has facilities for DMCreateInjection().
3116: Level: developer
3118: .seealso DMHasFunction(), DMCreateInjection()
3120: @*/
3121: PetscErrorCode DMHasCreateInjection(DM dm,PetscBool *flg)
3122: {
3125: if (!dm->ops->hascreateinjection) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DMHasCreateInjection not implemented for this type");
3126: (*dm->ops->hascreateinjection)(dm,flg);
3127: return(0);
3128: }
3130: /*@C
3131: DMSetVec - set the vector at which to compute residual, Jacobian and VI bounds, if the problem is nonlinear.
3133: Collective on DM
3135: Input Parameter:
3136: + dm - the DM object
3137: - x - location to compute residual and Jacobian, if NULL is passed to those routines; will be NULL for linear problems.
3139: Level: developer
3141: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext()
3143: @*/
3144: PetscErrorCode DMSetVec(DM dm,Vec x)
3145: {
3149: if (x) {
3150: if (!dm->x) {
3151: DMCreateGlobalVector(dm,&dm->x);
3152: }
3153: VecCopy(x,dm->x);
3154: } else if (dm->x) {
3155: VecDestroy(&dm->x);
3156: }
3157: return(0);
3158: }
3160: PetscFunctionList DMList = NULL;
3161: PetscBool DMRegisterAllCalled = PETSC_FALSE;
3163: /*@C
3164: DMSetType - Builds a DM, for a particular DM implementation.
3166: Collective on DM
3168: Input Parameters:
3169: + dm - The DM object
3170: - method - The name of the DM type
3172: Options Database Key:
3173: . -dm_type <type> - Sets the DM type; use -help for a list of available types
3175: Notes:
3176: See "petsc/include/petscdm.h" for available DM types (for instance, DM1D, DM2D, or DM3D).
3178: Level: intermediate
3180: .keywords: DM, set, type
3181: .seealso: DMGetType(), DMCreate()
3182: @*/
3183: PetscErrorCode DMSetType(DM dm, DMType method)
3184: {
3185: PetscErrorCode (*r)(DM);
3186: PetscBool match;
3191: PetscObjectTypeCompare((PetscObject) dm, method, &match);
3192: if (match) return(0);
3194: DMRegisterAll();
3195: PetscFunctionListFind(DMList,method,&r);
3196: if (!r) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown DM type: %s", method);
3198: if (dm->ops->destroy) {
3199: (*dm->ops->destroy)(dm);
3200: dm->ops->destroy = NULL;
3201: }
3202: (*r)(dm);
3203: PetscObjectChangeTypeName((PetscObject)dm,method);
3204: return(0);
3205: }
3207: /*@C
3208: DMGetType - Gets the DM type name (as a string) from the DM.
3210: Not Collective
3212: Input Parameter:
3213: . dm - The DM
3215: Output Parameter:
3216: . type - The DM type name
3218: Level: intermediate
3220: .keywords: DM, get, type, name
3221: .seealso: DMSetType(), DMCreate()
3222: @*/
3223: PetscErrorCode DMGetType(DM dm, DMType *type)
3224: {
3230: DMRegisterAll();
3231: *type = ((PetscObject)dm)->type_name;
3232: return(0);
3233: }
3235: /*@C
3236: DMConvert - Converts a DM to another DM, either of the same or different type.
3238: Collective on DM
3240: Input Parameters:
3241: + dm - the DM
3242: - newtype - new DM type (use "same" for the same type)
3244: Output Parameter:
3245: . M - pointer to new DM
3247: Notes:
3248: Cannot be used to convert a sequential DM to parallel or parallel to sequential,
3249: the MPI communicator of the generated DM is always the same as the communicator
3250: of the input DM.
3252: Level: intermediate
3254: .seealso: DMCreate()
3255: @*/
3256: PetscErrorCode DMConvert(DM dm, DMType newtype, DM *M)
3257: {
3258: DM B;
3259: char convname[256];
3260: PetscBool sametype/*, issame */;
3267: PetscObjectTypeCompare((PetscObject) dm, newtype, &sametype);
3268: /* PetscStrcmp(newtype, "same", &issame); */
3269: if (sametype) {
3270: *M = dm;
3271: PetscObjectReference((PetscObject) dm);
3272: return(0);
3273: } else {
3274: PetscErrorCode (*conv)(DM, DMType, DM*) = NULL;
3276: /*
3277: Order of precedence:
3278: 1) See if a specialized converter is known to the current DM.
3279: 2) See if a specialized converter is known to the desired DM class.
3280: 3) See if a good general converter is registered for the desired class
3281: 4) See if a good general converter is known for the current matrix.
3282: 5) Use a really basic converter.
3283: */
3285: /* 1) See if a specialized converter is known to the current DM and the desired class */
3286: PetscStrncpy(convname,"DMConvert_",sizeof(convname));
3287: PetscStrlcat(convname,((PetscObject) dm)->type_name,sizeof(convname));
3288: PetscStrlcat(convname,"_",sizeof(convname));
3289: PetscStrlcat(convname,newtype,sizeof(convname));
3290: PetscStrlcat(convname,"_C",sizeof(convname));
3291: PetscObjectQueryFunction((PetscObject)dm,convname,&conv);
3292: if (conv) goto foundconv;
3294: /* 2) See if a specialized converter is known to the desired DM class. */
3295: DMCreate(PetscObjectComm((PetscObject)dm), &B);
3296: DMSetType(B, newtype);
3297: PetscStrncpy(convname,"DMConvert_",sizeof(convname));
3298: PetscStrlcat(convname,((PetscObject) dm)->type_name,sizeof(convname));
3299: PetscStrlcat(convname,"_",sizeof(convname));
3300: PetscStrlcat(convname,newtype,sizeof(convname));
3301: PetscStrlcat(convname,"_C",sizeof(convname));
3302: PetscObjectQueryFunction((PetscObject)B,convname,&conv);
3303: if (conv) {
3304: DMDestroy(&B);
3305: goto foundconv;
3306: }
3308: #if 0
3309: /* 3) See if a good general converter is registered for the desired class */
3310: conv = B->ops->convertfrom;
3311: DMDestroy(&B);
3312: if (conv) goto foundconv;
3314: /* 4) See if a good general converter is known for the current matrix */
3315: if (dm->ops->convert) {
3316: conv = dm->ops->convert;
3317: }
3318: if (conv) goto foundconv;
3319: #endif
3321: /* 5) Use a really basic converter. */
3322: SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "No conversion possible between DM types %s and %s", ((PetscObject) dm)->type_name, newtype);
3324: foundconv:
3325: PetscLogEventBegin(DM_Convert,dm,0,0,0);
3326: (*conv)(dm,newtype,M);
3327: /* Things that are independent of DM type: We should consult DMClone() here */
3328: {
3329: PetscBool isper;
3330: const PetscReal *maxCell, *L;
3331: const DMBoundaryType *bd;
3332: DMGetPeriodicity(dm, &isper, &maxCell, &L, &bd);
3333: DMSetPeriodicity(*M, isper, maxCell, L, bd);
3334: }
3335: PetscLogEventEnd(DM_Convert,dm,0,0,0);
3336: }
3337: PetscObjectStateIncrease((PetscObject) *M);
3338: return(0);
3339: }
3341: /*--------------------------------------------------------------------------------------------------------------------*/
3343: /*@C
3344: DMRegister - Adds a new DM component implementation
3346: Not Collective
3348: Input Parameters:
3349: + name - The name of a new user-defined creation routine
3350: - create_func - The creation routine itself
3352: Notes:
3353: DMRegister() may be called multiple times to add several user-defined DMs
3356: Sample usage:
3357: .vb
3358: DMRegister("my_da", MyDMCreate);
3359: .ve
3361: Then, your DM type can be chosen with the procedural interface via
3362: .vb
3363: DMCreate(MPI_Comm, DM *);
3364: DMSetType(DM,"my_da");
3365: .ve
3366: or at runtime via the option
3367: .vb
3368: -da_type my_da
3369: .ve
3371: Level: advanced
3373: .keywords: DM, register
3374: .seealso: DMRegisterAll(), DMRegisterDestroy()
3376: @*/
3377: PetscErrorCode DMRegister(const char sname[],PetscErrorCode (*function)(DM))
3378: {
3382: PetscFunctionListAdd(&DMList,sname,function);
3383: return(0);
3384: }
3386: /*@C
3387: DMLoad - Loads a DM that has been stored in binary with DMView().
3389: Collective on PetscViewer
3391: Input Parameters:
3392: + newdm - the newly loaded DM, this needs to have been created with DMCreate() or
3393: some related function before a call to DMLoad().
3394: - viewer - binary file viewer, obtained from PetscViewerBinaryOpen() or
3395: HDF5 file viewer, obtained from PetscViewerHDF5Open()
3397: Level: intermediate
3399: Notes:
3400: The type is determined by the data in the file, any type set into the DM before this call is ignored.
3402: Notes for advanced users:
3403: Most users should not need to know the details of the binary storage
3404: format, since DMLoad() and DMView() completely hide these details.
3405: But for anyone who's interested, the standard binary matrix storage
3406: format is
3407: .vb
3408: has not yet been determined
3409: .ve
3411: .seealso: PetscViewerBinaryOpen(), DMView(), MatLoad(), VecLoad()
3412: @*/
3413: PetscErrorCode DMLoad(DM newdm, PetscViewer viewer)
3414: {
3415: PetscBool isbinary, ishdf5;
3421: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
3422: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);
3423: if (isbinary) {
3424: PetscInt classid;
3425: char type[256];
3427: PetscViewerBinaryRead(viewer,&classid,1,NULL,PETSC_INT);
3428: if (classid != DM_FILE_CLASSID) SETERRQ1(PetscObjectComm((PetscObject)newdm),PETSC_ERR_ARG_WRONG,"Not DM next in file, classid found %d",(int)classid);
3429: PetscViewerBinaryRead(viewer,type,256,NULL,PETSC_CHAR);
3430: DMSetType(newdm, type);
3431: if (newdm->ops->load) {(*newdm->ops->load)(newdm,viewer);}
3432: } else if (ishdf5) {
3433: if (newdm->ops->load) {(*newdm->ops->load)(newdm,viewer);}
3434: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid viewer; open viewer with PetscViewerBinaryOpen() or PetscViewerHDF5Open()");
3435: return(0);
3436: }
3438: /******************************** FEM Support **********************************/
3440: PetscErrorCode DMPrintCellVector(PetscInt c, const char name[], PetscInt len, const PetscScalar x[])
3441: {
3442: PetscInt f;
3446: PetscPrintf(PETSC_COMM_SELF, "Cell %D Element %s\n", c, name);
3447: for (f = 0; f < len; ++f) {
3448: PetscPrintf(PETSC_COMM_SELF, " | %g |\n", (double)PetscRealPart(x[f]));
3449: }
3450: return(0);
3451: }
3453: PetscErrorCode DMPrintCellMatrix(PetscInt c, const char name[], PetscInt rows, PetscInt cols, const PetscScalar A[])
3454: {
3455: PetscInt f, g;
3459: PetscPrintf(PETSC_COMM_SELF, "Cell %D Element %s\n", c, name);
3460: for (f = 0; f < rows; ++f) {
3461: PetscPrintf(PETSC_COMM_SELF, " |");
3462: for (g = 0; g < cols; ++g) {
3463: PetscPrintf(PETSC_COMM_SELF, " % 9.5g", PetscRealPart(A[f*cols+g]));
3464: }
3465: PetscPrintf(PETSC_COMM_SELF, " |\n");
3466: }
3467: return(0);
3468: }
3470: PetscErrorCode DMPrintLocalVec(DM dm, const char name[], PetscReal tol, Vec X)
3471: {
3472: PetscInt localSize, bs;
3473: PetscMPIInt size;
3474: Vec x, xglob;
3475: const PetscScalar *xarray;
3476: PetscErrorCode ierr;
3479: MPI_Comm_size(PetscObjectComm((PetscObject) dm),&size);
3480: VecDuplicate(X, &x);
3481: VecCopy(X, x);
3482: VecChop(x, tol);
3483: PetscPrintf(PetscObjectComm((PetscObject) dm),"%s:\n",name);
3484: if (size > 1) {
3485: VecGetLocalSize(x,&localSize);
3486: VecGetArrayRead(x,&xarray);
3487: VecGetBlockSize(x,&bs);
3488: VecCreateMPIWithArray(PetscObjectComm((PetscObject) dm),bs,localSize,PETSC_DETERMINE,xarray,&xglob);
3489: } else {
3490: xglob = x;
3491: }
3492: VecView(xglob,PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject) dm)));
3493: if (size > 1) {
3494: VecDestroy(&xglob);
3495: VecRestoreArrayRead(x,&xarray);
3496: }
3497: VecDestroy(&x);
3498: return(0);
3499: }
3501: /*@
3502: DMGetDefaultSection - Get the PetscSection encoding the local data layout for the DM.
3504: Input Parameter:
3505: . dm - The DM
3507: Output Parameter:
3508: . section - The PetscSection
3510: Level: intermediate
3512: Note: This gets a borrowed reference, so the user should not destroy this PetscSection.
3514: .seealso: DMSetDefaultSection(), DMGetDefaultGlobalSection()
3515: @*/
3516: PetscErrorCode DMGetDefaultSection(DM dm, PetscSection *section)
3517: {
3523: if (!dm->defaultSection && dm->ops->createdefaultsection) {
3524: (*dm->ops->createdefaultsection)(dm);
3525: if (dm->defaultSection) {PetscObjectViewFromOptions((PetscObject) dm->defaultSection, NULL, "-dm_petscsection_view");}
3526: }
3527: *section = dm->defaultSection;
3528: return(0);
3529: }
3531: /*@
3532: DMSetDefaultSection - Set the PetscSection encoding the local data layout for the DM.
3534: Input Parameters:
3535: + dm - The DM
3536: - section - The PetscSection
3538: Level: intermediate
3540: Note: Any existing Section will be destroyed
3542: .seealso: DMSetDefaultSection(), DMGetDefaultGlobalSection()
3543: @*/
3544: PetscErrorCode DMSetDefaultSection(DM dm, PetscSection section)
3545: {
3546: PetscInt numFields = 0;
3547: PetscInt f;
3552: if (section) {
3554: PetscObjectReference((PetscObject)section);
3555: }
3556: PetscSectionDestroy(&dm->defaultSection);
3557: dm->defaultSection = section;
3558: if (section) {PetscSectionGetNumFields(dm->defaultSection, &numFields);}
3559: if (numFields) {
3560: DMSetNumFields(dm, numFields);
3561: for (f = 0; f < numFields; ++f) {
3562: PetscObject disc;
3563: const char *name;
3565: PetscSectionGetFieldName(dm->defaultSection, f, &name);
3566: DMGetField(dm, f, &disc);
3567: PetscObjectSetName(disc, name);
3568: }
3569: }
3570: /* The global section will be rebuilt in the next call to DMGetDefaultGlobalSection(). */
3571: PetscSectionDestroy(&dm->defaultGlobalSection);
3572: return(0);
3573: }
3575: /*@
3576: DMGetDefaultConstraints - Get the PetscSection and Mat the specify the local constraint interpolation. See DMSetDefaultConstraints() for a description of the purpose of constraint interpolation.
3578: not collective
3580: Input Parameter:
3581: . dm - The DM
3583: Output Parameter:
3584: + 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.
3585: - 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.
3587: Level: advanced
3589: Note: This gets borrowed references, so the user should not destroy the PetscSection or the Mat.
3591: .seealso: DMSetDefaultConstraints()
3592: @*/
3593: PetscErrorCode DMGetDefaultConstraints(DM dm, PetscSection *section, Mat *mat)
3594: {
3599: if (!dm->defaultConstraintSection && !dm->defaultConstraintMat && dm->ops->createdefaultconstraints) {(*dm->ops->createdefaultconstraints)(dm);}
3600: if (section) {*section = dm->defaultConstraintSection;}
3601: if (mat) {*mat = dm->defaultConstraintMat;}
3602: return(0);
3603: }
3605: /*@
3606: DMSetDefaultConstraints - Set the PetscSection and Mat the specify the local constraint interpolation.
3608: 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().
3610: 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.
3612: collective on dm
3614: Input Parameters:
3615: + dm - The DM
3616: + 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).
3617: - 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).
3619: Level: advanced
3621: Note: This increments the references of the PetscSection and the Mat, so they user can destroy them
3623: .seealso: DMGetDefaultConstraints()
3624: @*/
3625: PetscErrorCode DMSetDefaultConstraints(DM dm, PetscSection section, Mat mat)
3626: {
3627: PetscMPIInt result;
3632: if (section) {
3634: MPI_Comm_compare(PETSC_COMM_SELF,PetscObjectComm((PetscObject)section),&result);
3635: if (result != MPI_CONGRUENT && result != MPI_IDENT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMECOMM,"constraint section must have local communicator");
3636: }
3637: if (mat) {
3639: MPI_Comm_compare(PETSC_COMM_SELF,PetscObjectComm((PetscObject)mat),&result);
3640: if (result != MPI_CONGRUENT && result != MPI_IDENT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMECOMM,"constraint matrix must have local communicator");
3641: }
3642: PetscObjectReference((PetscObject)section);
3643: PetscSectionDestroy(&dm->defaultConstraintSection);
3644: dm->defaultConstraintSection = section;
3645: PetscObjectReference((PetscObject)mat);
3646: MatDestroy(&dm->defaultConstraintMat);
3647: dm->defaultConstraintMat = mat;
3648: return(0);
3649: }
3651: #ifdef PETSC_USE_DEBUG
3652: /*
3653: DMDefaultSectionCheckConsistency - Check the consistentcy of the global and local sections.
3655: Input Parameters:
3656: + dm - The DM
3657: . localSection - PetscSection describing the local data layout
3658: - globalSection - PetscSection describing the global data layout
3660: Level: intermediate
3662: .seealso: DMGetDefaultSF(), DMSetDefaultSF()
3663: */
3664: static PetscErrorCode DMDefaultSectionCheckConsistency_Internal(DM dm, PetscSection localSection, PetscSection globalSection)
3665: {
3666: MPI_Comm comm;
3667: PetscLayout layout;
3668: const PetscInt *ranges;
3669: PetscInt pStart, pEnd, p, nroots;
3670: PetscMPIInt size, rank;
3671: PetscBool valid = PETSC_TRUE, gvalid;
3672: PetscErrorCode ierr;
3675: PetscObjectGetComm((PetscObject)dm,&comm);
3677: MPI_Comm_size(comm, &size);
3678: MPI_Comm_rank(comm, &rank);
3679: PetscSectionGetChart(globalSection, &pStart, &pEnd);
3680: PetscSectionGetConstrainedStorageSize(globalSection, &nroots);
3681: PetscLayoutCreate(comm, &layout);
3682: PetscLayoutSetBlockSize(layout, 1);
3683: PetscLayoutSetLocalSize(layout, nroots);
3684: PetscLayoutSetUp(layout);
3685: PetscLayoutGetRanges(layout, &ranges);
3686: for (p = pStart; p < pEnd; ++p) {
3687: PetscInt dof, cdof, off, gdof, gcdof, goff, gsize, d;
3689: PetscSectionGetDof(localSection, p, &dof);
3690: PetscSectionGetOffset(localSection, p, &off);
3691: PetscSectionGetConstraintDof(localSection, p, &cdof);
3692: PetscSectionGetDof(globalSection, p, &gdof);
3693: PetscSectionGetConstraintDof(globalSection, p, &gcdof);
3694: PetscSectionGetOffset(globalSection, p, &goff);
3695: if (!gdof) continue; /* Censored point */
3696: 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;}
3697: 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;}
3698: if (gdof < 0) {
3699: gsize = gdof < 0 ? -(gdof+1)-gcdof : gdof-gcdof;
3700: for (d = 0; d < gsize; ++d) {
3701: PetscInt offset = -(goff+1) + d, r;
3703: PetscFindInt(offset,size+1,ranges,&r);
3704: if (r < 0) r = -(r+2);
3705: 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;}
3706: }
3707: }
3708: }
3709: PetscLayoutDestroy(&layout);
3710: PetscSynchronizedFlush(comm, NULL);
3711: MPIU_Allreduce(&valid, &gvalid, 1, MPIU_BOOL, MPI_LAND, comm);
3712: if (!gvalid) {
3713: DMView(dm, NULL);
3714: SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Inconsistent local and global sections");
3715: }
3716: return(0);
3717: }
3718: #endif
3720: /*@
3721: DMGetDefaultGlobalSection - Get the PetscSection encoding the global data layout for the DM.
3723: Collective on DM
3725: Input Parameter:
3726: . dm - The DM
3728: Output Parameter:
3729: . section - The PetscSection
3731: Level: intermediate
3733: Note: This gets a borrowed reference, so the user should not destroy this PetscSection.
3735: .seealso: DMSetDefaultSection(), DMGetDefaultSection()
3736: @*/
3737: PetscErrorCode DMGetDefaultGlobalSection(DM dm, PetscSection *section)
3738: {
3744: if (!dm->defaultGlobalSection) {
3745: PetscSection s;
3747: DMGetDefaultSection(dm, &s);
3748: if (!s) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONGSTATE, "DM must have a default PetscSection in order to create a global PetscSection");
3749: if (!dm->sf) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DM must have a default PetscSF in order to create a global PetscSection");
3750: PetscSectionCreateGlobalSection(s, dm->sf, PETSC_FALSE, PETSC_FALSE, &dm->defaultGlobalSection);
3751: PetscLayoutDestroy(&dm->map);
3752: PetscSectionGetValueLayout(PetscObjectComm((PetscObject)dm), dm->defaultGlobalSection, &dm->map);
3753: PetscSectionViewFromOptions(dm->defaultGlobalSection, NULL, "-global_section_view");
3754: }
3755: *section = dm->defaultGlobalSection;
3756: return(0);
3757: }
3759: /*@
3760: DMSetDefaultGlobalSection - Set the PetscSection encoding the global data layout for the DM.
3762: Input Parameters:
3763: + dm - The DM
3764: - section - The PetscSection, or NULL
3766: Level: intermediate
3768: Note: Any existing Section will be destroyed
3770: .seealso: DMGetDefaultGlobalSection(), DMSetDefaultSection()
3771: @*/
3772: PetscErrorCode DMSetDefaultGlobalSection(DM dm, PetscSection section)
3773: {
3779: PetscObjectReference((PetscObject)section);
3780: PetscSectionDestroy(&dm->defaultGlobalSection);
3781: dm->defaultGlobalSection = section;
3782: #ifdef PETSC_USE_DEBUG
3783: if (section) {DMDefaultSectionCheckConsistency_Internal(dm, dm->defaultSection, section);}
3784: #endif
3785: return(0);
3786: }
3788: /*@
3789: DMGetDefaultSF - Get the PetscSF encoding the parallel dof overlap for the DM. If it has not been set,
3790: it is created from the default PetscSection layouts in the DM.
3792: Input Parameter:
3793: . dm - The DM
3795: Output Parameter:
3796: . sf - The PetscSF
3798: Level: intermediate
3800: Note: This gets a borrowed reference, so the user should not destroy this PetscSF.
3802: .seealso: DMSetDefaultSF(), DMCreateDefaultSF()
3803: @*/
3804: PetscErrorCode DMGetDefaultSF(DM dm, PetscSF *sf)
3805: {
3806: PetscInt nroots;
3812: PetscSFGetGraph(dm->defaultSF, &nroots, NULL, NULL, NULL);
3813: if (nroots < 0) {
3814: PetscSection section, gSection;
3816: DMGetDefaultSection(dm, §ion);
3817: if (section) {
3818: DMGetDefaultGlobalSection(dm, &gSection);
3819: DMCreateDefaultSF(dm, section, gSection);
3820: } else {
3821: *sf = NULL;
3822: return(0);
3823: }
3824: }
3825: *sf = dm->defaultSF;
3826: return(0);
3827: }
3829: /*@
3830: DMSetDefaultSF - Set the PetscSF encoding the parallel dof overlap for the DM
3832: Input Parameters:
3833: + dm - The DM
3834: - sf - The PetscSF
3836: Level: intermediate
3838: Note: Any previous SF is destroyed
3840: .seealso: DMGetDefaultSF(), DMCreateDefaultSF()
3841: @*/
3842: PetscErrorCode DMSetDefaultSF(DM dm, PetscSF sf)
3843: {
3849: PetscSFDestroy(&dm->defaultSF);
3850: dm->defaultSF = sf;
3851: return(0);
3852: }
3854: /*@C
3855: DMCreateDefaultSF - Create the PetscSF encoding the parallel dof overlap for the DM based upon the PetscSections
3856: describing the data layout.
3858: Input Parameters:
3859: + dm - The DM
3860: . localSection - PetscSection describing the local data layout
3861: - globalSection - PetscSection describing the global data layout
3863: Level: intermediate
3865: .seealso: DMGetDefaultSF(), DMSetDefaultSF()
3866: @*/
3867: PetscErrorCode DMCreateDefaultSF(DM dm, PetscSection localSection, PetscSection globalSection)
3868: {
3869: MPI_Comm comm;
3870: PetscLayout layout;
3871: const PetscInt *ranges;
3872: PetscInt *local;
3873: PetscSFNode *remote;
3874: PetscInt pStart, pEnd, p, nroots, nleaves = 0, l;
3875: PetscMPIInt size, rank;
3879: PetscObjectGetComm((PetscObject)dm,&comm);
3881: MPI_Comm_size(comm, &size);
3882: MPI_Comm_rank(comm, &rank);
3883: PetscSectionGetChart(globalSection, &pStart, &pEnd);
3884: PetscSectionGetConstrainedStorageSize(globalSection, &nroots);
3885: PetscLayoutCreate(comm, &layout);
3886: PetscLayoutSetBlockSize(layout, 1);
3887: PetscLayoutSetLocalSize(layout, nroots);
3888: PetscLayoutSetUp(layout);
3889: PetscLayoutGetRanges(layout, &ranges);
3890: for (p = pStart; p < pEnd; ++p) {
3891: PetscInt gdof, gcdof;
3893: PetscSectionGetDof(globalSection, p, &gdof);
3894: PetscSectionGetConstraintDof(globalSection, p, &gcdof);
3895: 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));
3896: nleaves += gdof < 0 ? -(gdof+1)-gcdof : gdof-gcdof;
3897: }
3898: PetscMalloc1(nleaves, &local);
3899: PetscMalloc1(nleaves, &remote);
3900: for (p = pStart, l = 0; p < pEnd; ++p) {
3901: const PetscInt *cind;
3902: PetscInt dof, cdof, off, gdof, gcdof, goff, gsize, d, c;
3904: PetscSectionGetDof(localSection, p, &dof);
3905: PetscSectionGetOffset(localSection, p, &off);
3906: PetscSectionGetConstraintDof(localSection, p, &cdof);
3907: PetscSectionGetConstraintIndices(localSection, p, &cind);
3908: PetscSectionGetDof(globalSection, p, &gdof);
3909: PetscSectionGetConstraintDof(globalSection, p, &gcdof);
3910: PetscSectionGetOffset(globalSection, p, &goff);
3911: if (!gdof) continue; /* Censored point */
3912: gsize = gdof < 0 ? -(gdof+1)-gcdof : gdof-gcdof;
3913: if (gsize != dof-cdof) {
3914: 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);
3915: cdof = 0; /* Ignore constraints */
3916: }
3917: for (d = 0, c = 0; d < dof; ++d) {
3918: if ((c < cdof) && (cind[c] == d)) {++c; continue;}
3919: local[l+d-c] = off+d;
3920: }
3921: if (gdof < 0) {
3922: for (d = 0; d < gsize; ++d, ++l) {
3923: PetscInt offset = -(goff+1) + d, r;
3925: PetscFindInt(offset,size+1,ranges,&r);
3926: if (r < 0) r = -(r+2);
3927: 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);
3928: remote[l].rank = r;
3929: remote[l].index = offset - ranges[r];
3930: }
3931: } else {
3932: for (d = 0; d < gsize; ++d, ++l) {
3933: remote[l].rank = rank;
3934: remote[l].index = goff+d - ranges[rank];
3935: }
3936: }
3937: }
3938: if (l != nleaves) SETERRQ2(comm, PETSC_ERR_PLIB, "Iteration error, l %d != nleaves %d", l, nleaves);
3939: PetscLayoutDestroy(&layout);
3940: PetscSFSetGraph(dm->defaultSF, nroots, nleaves, local, PETSC_OWN_POINTER, remote, PETSC_OWN_POINTER);
3941: return(0);
3942: }
3944: /*@
3945: DMGetPointSF - Get the PetscSF encoding the parallel section point overlap for the DM.
3947: Input Parameter:
3948: . dm - The DM
3950: Output Parameter:
3951: . sf - The PetscSF
3953: Level: intermediate
3955: Note: This gets a borrowed reference, so the user should not destroy this PetscSF.
3957: .seealso: DMSetPointSF(), DMGetDefaultSF(), DMSetDefaultSF(), DMCreateDefaultSF()
3958: @*/
3959: PetscErrorCode DMGetPointSF(DM dm, PetscSF *sf)
3960: {
3964: *sf = dm->sf;
3965: return(0);
3966: }
3968: /*@
3969: DMSetPointSF - Set the PetscSF encoding the parallel section point overlap for the DM.
3971: Input Parameters:
3972: + dm - The DM
3973: - sf - The PetscSF
3975: Level: intermediate
3977: .seealso: DMGetPointSF(), DMGetDefaultSF(), DMSetDefaultSF(), DMCreateDefaultSF()
3978: @*/
3979: PetscErrorCode DMSetPointSF(DM dm, PetscSF sf)
3980: {
3986: PetscSFDestroy(&dm->sf);
3987: PetscObjectReference((PetscObject) sf);
3988: dm->sf = sf;
3989: return(0);
3990: }
3992: /*@
3993: DMGetDS - Get the PetscDS
3995: Input Parameter:
3996: . dm - The DM
3998: Output Parameter:
3999: . prob - The PetscDS
4001: Level: developer
4003: .seealso: DMSetDS()
4004: @*/
4005: PetscErrorCode DMGetDS(DM dm, PetscDS *prob)
4006: {
4010: *prob = dm->prob;
4011: return(0);
4012: }
4014: /*@
4015: DMSetDS - Set the PetscDS
4017: Input Parameters:
4018: + dm - The DM
4019: - prob - The PetscDS
4021: Level: developer
4023: .seealso: DMGetDS()
4024: @*/
4025: PetscErrorCode DMSetDS(DM dm, PetscDS prob)
4026: {
4027: PetscInt dimEmbed;
4033: PetscObjectReference((PetscObject) prob);
4034: PetscDSDestroy(&dm->prob);
4035: dm->prob = prob;
4036: DMGetCoordinateDim(dm, &dimEmbed);
4037: PetscDSSetCoordinateDimension(prob, dimEmbed);
4038: return(0);
4039: }
4041: PetscErrorCode DMGetNumFields(DM dm, PetscInt *numFields)
4042: {
4047: PetscDSGetNumFields(dm->prob, numFields);
4048: return(0);
4049: }
4051: PetscErrorCode DMSetNumFields(DM dm, PetscInt numFields)
4052: {
4053: PetscInt Nf, f;
4058: PetscDSGetNumFields(dm->prob, &Nf);
4059: for (f = Nf; f < numFields; ++f) {
4060: PetscContainer obj;
4062: PetscContainerCreate(PetscObjectComm((PetscObject) dm), &obj);
4063: PetscDSSetDiscretization(dm->prob, f, (PetscObject) obj);
4064: PetscContainerDestroy(&obj);
4065: }
4066: return(0);
4067: }
4069: /*@
4070: DMGetField - Return the discretization object for a given DM field
4072: Not collective
4074: Input Parameters:
4075: + dm - The DM
4076: - f - The field number
4078: Output Parameter:
4079: . field - The discretization object
4081: Level: developer
4083: .seealso: DMSetField()
4084: @*/
4085: PetscErrorCode DMGetField(DM dm, PetscInt f, PetscObject *field)
4086: {
4091: PetscDSGetDiscretization(dm->prob, f, field);
4092: return(0);
4093: }
4095: /*@
4096: DMSetField - Set the discretization object for a given DM field
4098: Logically collective on DM
4100: Input Parameters:
4101: + dm - The DM
4102: . f - The field number
4103: - field - The discretization object
4105: Level: developer
4107: .seealso: DMGetField()
4108: @*/
4109: PetscErrorCode DMSetField(DM dm, PetscInt f, PetscObject field)
4110: {
4115: PetscDSSetDiscretization(dm->prob, f, field);
4116: return(0);
4117: }
4119: PetscErrorCode DMRestrictHook_Coordinates(DM dm,DM dmc,void *ctx)
4120: {
4121: DM dm_coord,dmc_coord;
4123: Vec coords,ccoords;
4124: Mat inject;
4126: DMGetCoordinateDM(dm,&dm_coord);
4127: DMGetCoordinateDM(dmc,&dmc_coord);
4128: DMGetCoordinates(dm,&coords);
4129: DMGetCoordinates(dmc,&ccoords);
4130: if (coords && !ccoords) {
4131: DMCreateGlobalVector(dmc_coord,&ccoords);
4132: PetscObjectSetName((PetscObject)ccoords,"coordinates");
4133: DMCreateInjection(dmc_coord,dm_coord,&inject);
4134: MatRestrict(inject,coords,ccoords);
4135: MatDestroy(&inject);
4136: DMSetCoordinates(dmc,ccoords);
4137: VecDestroy(&ccoords);
4138: }
4139: return(0);
4140: }
4142: static PetscErrorCode DMSubDomainHook_Coordinates(DM dm,DM subdm,void *ctx)
4143: {
4144: DM dm_coord,subdm_coord;
4146: Vec coords,ccoords,clcoords;
4147: VecScatter *scat_i,*scat_g;
4149: DMGetCoordinateDM(dm,&dm_coord);
4150: DMGetCoordinateDM(subdm,&subdm_coord);
4151: DMGetCoordinates(dm,&coords);
4152: DMGetCoordinates(subdm,&ccoords);
4153: if (coords && !ccoords) {
4154: DMCreateGlobalVector(subdm_coord,&ccoords);
4155: PetscObjectSetName((PetscObject)ccoords,"coordinates");
4156: DMCreateLocalVector(subdm_coord,&clcoords);
4157: PetscObjectSetName((PetscObject)clcoords,"coordinates");
4158: DMCreateDomainDecompositionScatters(dm_coord,1,&subdm_coord,NULL,&scat_i,&scat_g);
4159: VecScatterBegin(scat_i[0],coords,ccoords,INSERT_VALUES,SCATTER_FORWARD);
4160: VecScatterBegin(scat_g[0],coords,clcoords,INSERT_VALUES,SCATTER_FORWARD);
4161: VecScatterEnd(scat_i[0],coords,ccoords,INSERT_VALUES,SCATTER_FORWARD);
4162: VecScatterEnd(scat_g[0],coords,clcoords,INSERT_VALUES,SCATTER_FORWARD);
4163: DMSetCoordinates(subdm,ccoords);
4164: DMSetCoordinatesLocal(subdm,clcoords);
4165: VecScatterDestroy(&scat_i[0]);
4166: VecScatterDestroy(&scat_g[0]);
4167: VecDestroy(&ccoords);
4168: VecDestroy(&clcoords);
4169: PetscFree(scat_i);
4170: PetscFree(scat_g);
4171: }
4172: return(0);
4173: }
4175: /*@
4176: DMGetDimension - Return the topological dimension of the DM
4178: Not collective
4180: Input Parameter:
4181: . dm - The DM
4183: Output Parameter:
4184: . dim - The topological dimension
4186: Level: beginner
4188: .seealso: DMSetDimension(), DMCreate()
4189: @*/
4190: PetscErrorCode DMGetDimension(DM dm, PetscInt *dim)
4191: {
4195: *dim = dm->dim;
4196: return(0);
4197: }
4199: /*@
4200: DMSetDimension - Set the topological dimension of the DM
4202: Collective on dm
4204: Input Parameters:
4205: + dm - The DM
4206: - dim - The topological dimension
4208: Level: beginner
4210: .seealso: DMGetDimension(), DMCreate()
4211: @*/
4212: PetscErrorCode DMSetDimension(DM dm, PetscInt dim)
4213: {
4219: dm->dim = dim;
4220: if (dm->prob->dimEmbed < 0) {PetscDSSetCoordinateDimension(dm->prob, dm->dim);}
4221: return(0);
4222: }
4224: /*@
4225: DMGetDimPoints - Get the half-open interval for all points of a given dimension
4227: Collective on DM
4229: Input Parameters:
4230: + dm - the DM
4231: - dim - the dimension
4233: Output Parameters:
4234: + pStart - The first point of the given dimension
4235: . pEnd - The first point following points of the given dimension
4237: Note:
4238: The points are vertices in the Hasse diagram encoding the topology. This is explained in
4239: http://arxiv.org/abs/0908.4427. If not points exist of this dimension in the storage scheme,
4240: then the interval is empty.
4242: Level: intermediate
4244: .keywords: point, Hasse Diagram, dimension
4245: .seealso: DMPLEX, DMPlexGetDepthStratum(), DMPlexGetHeightStratum()
4246: @*/
4247: PetscErrorCode DMGetDimPoints(DM dm, PetscInt dim, PetscInt *pStart, PetscInt *pEnd)
4248: {
4249: PetscInt d;
4254: DMGetDimension(dm, &d);
4255: if ((dim < 0) || (dim > d)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension %d 1", dim, d);
4256: (*dm->ops->getdimpoints)(dm, dim, pStart, pEnd);
4257: return(0);
4258: }
4260: /*@
4261: DMSetCoordinates - Sets into the DM a global vector that holds the coordinates
4263: Collective on DM
4265: Input Parameters:
4266: + dm - the DM
4267: - c - coordinate vector
4269: Note:
4270: The coordinates do include those for ghost points, which are in the local vector
4272: Level: intermediate
4274: .keywords: distributed array, get, corners, nodes, local indices, coordinates
4275: .seealso: DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLoca(), DMGetCoordinateDM()
4276: @*/
4277: PetscErrorCode DMSetCoordinates(DM dm, Vec c)
4278: {
4284: PetscObjectReference((PetscObject) c);
4285: VecDestroy(&dm->coordinates);
4286: dm->coordinates = c;
4287: VecDestroy(&dm->coordinatesLocal);
4288: DMCoarsenHookAdd(dm,DMRestrictHook_Coordinates,NULL,NULL);
4289: DMSubDomainHookAdd(dm,DMSubDomainHook_Coordinates,NULL,NULL);
4290: return(0);
4291: }
4293: /*@
4294: DMSetCoordinatesLocal - Sets into the DM a local vector that holds the coordinates
4296: Collective on DM
4298: Input Parameters:
4299: + dm - the DM
4300: - c - coordinate vector
4302: Note:
4303: The coordinates of ghost points can be set using DMSetCoordinates()
4304: followed by DMGetCoordinatesLocal(). This is intended to enable the
4305: setting of ghost coordinates outside of the domain.
4307: Level: intermediate
4309: .keywords: distributed array, get, corners, nodes, local indices, coordinates
4310: .seealso: DMGetCoordinatesLocal(), DMSetCoordinates(), DMGetCoordinates(), DMGetCoordinateDM()
4311: @*/
4312: PetscErrorCode DMSetCoordinatesLocal(DM dm, Vec c)
4313: {
4319: PetscObjectReference((PetscObject) c);
4320: VecDestroy(&dm->coordinatesLocal);
4322: dm->coordinatesLocal = c;
4324: VecDestroy(&dm->coordinates);
4325: return(0);
4326: }
4328: /*@
4329: DMGetCoordinates - Gets a global vector with the coordinates associated with the DM.
4331: Not Collective
4333: Input Parameter:
4334: . dm - the DM
4336: Output Parameter:
4337: . c - global coordinate vector
4339: Note:
4340: This is a borrowed reference, so the user should NOT destroy this vector
4342: Each process has only the local coordinates (does NOT have the ghost coordinates).
4344: For DMDA, in two and three dimensions coordinates are interlaced (x_0,y_0,x_1,y_1,...)
4345: and (x_0,y_0,z_0,x_1,y_1,z_1...)
4347: Level: intermediate
4349: .keywords: distributed array, get, corners, nodes, local indices, coordinates
4350: .seealso: DMSetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM()
4351: @*/
4352: PetscErrorCode DMGetCoordinates(DM dm, Vec *c)
4353: {
4359: if (!dm->coordinates && dm->coordinatesLocal) {
4360: DM cdm = NULL;
4362: DMGetCoordinateDM(dm, &cdm);
4363: DMCreateGlobalVector(cdm, &dm->coordinates);
4364: PetscObjectSetName((PetscObject) dm->coordinates, "coordinates");
4365: DMLocalToGlobalBegin(cdm, dm->coordinatesLocal, INSERT_VALUES, dm->coordinates);
4366: DMLocalToGlobalEnd(cdm, dm->coordinatesLocal, INSERT_VALUES, dm->coordinates);
4367: }
4368: *c = dm->coordinates;
4369: return(0);
4370: }
4372: /*@
4373: DMGetCoordinatesLocal - Gets a local vector with the coordinates associated with the DM.
4375: Collective on DM
4377: Input Parameter:
4378: . dm - the DM
4380: Output Parameter:
4381: . c - coordinate vector
4383: Note:
4384: This is a borrowed reference, so the user should NOT destroy this vector
4386: Each process has the local and ghost coordinates
4388: For DMDA, in two and three dimensions coordinates are interlaced (x_0,y_0,x_1,y_1,...)
4389: and (x_0,y_0,z_0,x_1,y_1,z_1...)
4391: Level: intermediate
4393: .keywords: distributed array, get, corners, nodes, local indices, coordinates
4394: .seealso: DMSetCoordinatesLocal(), DMGetCoordinates(), DMSetCoordinates(), DMGetCoordinateDM()
4395: @*/
4396: PetscErrorCode DMGetCoordinatesLocal(DM dm, Vec *c)
4397: {
4403: if (!dm->coordinatesLocal && dm->coordinates) {
4404: DM cdm = NULL;
4406: DMGetCoordinateDM(dm, &cdm);
4407: DMCreateLocalVector(cdm, &dm->coordinatesLocal);
4408: PetscObjectSetName((PetscObject) dm->coordinatesLocal, "coordinates");
4409: DMGlobalToLocalBegin(cdm, dm->coordinates, INSERT_VALUES, dm->coordinatesLocal);
4410: DMGlobalToLocalEnd(cdm, dm->coordinates, INSERT_VALUES, dm->coordinatesLocal);
4411: }
4412: *c = dm->coordinatesLocal;
4413: return(0);
4414: }
4416: /*@
4417: DMGetCoordinateDM - Gets the DM that prescribes coordinate layout and scatters between global and local coordinates
4419: Collective on DM
4421: Input Parameter:
4422: . dm - the DM
4424: Output Parameter:
4425: . cdm - coordinate DM
4427: Level: intermediate
4429: .keywords: distributed array, get, corners, nodes, local indices, coordinates
4430: .seealso: DMSetCoordinateDM(), DMSetCoordinates(), DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLocal()
4431: @*/
4432: PetscErrorCode DMGetCoordinateDM(DM dm, DM *cdm)
4433: {
4439: if (!dm->coordinateDM) {
4440: if (!dm->ops->createcoordinatedm) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unable to create coordinates for this DM");
4441: (*dm->ops->createcoordinatedm)(dm, &dm->coordinateDM);
4442: }
4443: *cdm = dm->coordinateDM;
4444: return(0);
4445: }
4447: /*@
4448: DMSetCoordinateDM - Sets the DM that prescribes coordinate layout and scatters between global and local coordinates
4450: Logically Collective on DM
4452: Input Parameters:
4453: + dm - the DM
4454: - cdm - coordinate DM
4456: Level: intermediate
4458: .keywords: distributed array, get, corners, nodes, local indices, coordinates
4459: .seealso: DMGetCoordinateDM(), DMSetCoordinates(), DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLocal()
4460: @*/
4461: PetscErrorCode DMSetCoordinateDM(DM dm, DM cdm)
4462: {
4468: PetscObjectReference((PetscObject)cdm);
4469: DMDestroy(&dm->coordinateDM);
4470: dm->coordinateDM = cdm;
4471: return(0);
4472: }
4474: /*@
4475: DMGetCoordinateDim - Retrieve the dimension of embedding space for coordinate values.
4477: Not Collective
4479: Input Parameter:
4480: . dm - The DM object
4482: Output Parameter:
4483: . dim - The embedding dimension
4485: Level: intermediate
4487: .keywords: mesh, coordinates
4488: .seealso: DMSetCoordinateDim(), DMGetCoordinateSection(), DMGetCoordinateDM(), DMGetDefaultSection(), DMSetDefaultSection()
4489: @*/
4490: PetscErrorCode DMGetCoordinateDim(DM dm, PetscInt *dim)
4491: {
4495: if (dm->dimEmbed == PETSC_DEFAULT) {
4496: dm->dimEmbed = dm->dim;
4497: }
4498: *dim = dm->dimEmbed;
4499: return(0);
4500: }
4502: /*@
4503: DMSetCoordinateDim - Set the dimension of the embedding space for coordinate values.
4505: Not Collective
4507: Input Parameters:
4508: + dm - The DM object
4509: - dim - The embedding dimension
4511: Level: intermediate
4513: .keywords: mesh, coordinates
4514: .seealso: DMGetCoordinateDim(), DMSetCoordinateSection(), DMGetCoordinateSection(), DMGetDefaultSection(), DMSetDefaultSection()
4515: @*/
4516: PetscErrorCode DMSetCoordinateDim(DM dm, PetscInt dim)
4517: {
4522: dm->dimEmbed = dim;
4523: PetscDSSetCoordinateDimension(dm->prob, dm->dimEmbed);
4524: return(0);
4525: }
4527: /*@
4528: DMGetCoordinateSection - Retrieve the layout of coordinate values over the mesh.
4530: Not Collective
4532: Input Parameter:
4533: . dm - The DM object
4535: Output Parameter:
4536: . section - The PetscSection object
4538: Level: intermediate
4540: .keywords: mesh, coordinates
4541: .seealso: DMGetCoordinateDM(), DMGetDefaultSection(), DMSetDefaultSection()
4542: @*/
4543: PetscErrorCode DMGetCoordinateSection(DM dm, PetscSection *section)
4544: {
4545: DM cdm;
4551: DMGetCoordinateDM(dm, &cdm);
4552: DMGetDefaultSection(cdm, section);
4553: return(0);
4554: }
4556: /*@
4557: DMSetCoordinateSection - Set the layout of coordinate values over the mesh.
4559: Not Collective
4561: Input Parameters:
4562: + dm - The DM object
4563: . dim - The embedding dimension, or PETSC_DETERMINE
4564: - section - The PetscSection object
4566: Level: intermediate
4568: .keywords: mesh, coordinates
4569: .seealso: DMGetCoordinateSection(), DMGetDefaultSection(), DMSetDefaultSection()
4570: @*/
4571: PetscErrorCode DMSetCoordinateSection(DM dm, PetscInt dim, PetscSection section)
4572: {
4573: DM cdm;
4579: DMGetCoordinateDM(dm, &cdm);
4580: DMSetDefaultSection(cdm, section);
4581: if (dim == PETSC_DETERMINE) {
4582: PetscInt d = PETSC_DEFAULT;
4583: PetscInt pStart, pEnd, vStart, vEnd, v, dd;
4585: PetscSectionGetChart(section, &pStart, &pEnd);
4586: DMGetDimPoints(dm, 0, &vStart, &vEnd);
4587: pStart = PetscMax(vStart, pStart);
4588: pEnd = PetscMin(vEnd, pEnd);
4589: for (v = pStart; v < pEnd; ++v) {
4590: PetscSectionGetDof(section, v, &dd);
4591: if (dd) {d = dd; break;}
4592: }
4593: if (d < 0) d = PETSC_DEFAULT;
4594: DMSetCoordinateDim(dm, d);
4595: }
4596: return(0);
4597: }
4599: /*@C
4600: DMGetPeriodicity - Get the description of mesh periodicity
4602: Input Parameters:
4603: . dm - The DM object
4605: Output Parameters:
4606: + per - Whether the DM is periodic or not
4607: . maxCell - Over distances greater than this, we can assume a point has crossed over to another sheet, when trying to localize cell coordinates
4608: . L - If we assume the mesh is a torus, this is the length of each coordinate
4609: - bd - This describes the type of periodicity in each topological dimension
4611: Level: developer
4613: .seealso: DMGetPeriodicity()
4614: @*/
4615: PetscErrorCode DMGetPeriodicity(DM dm, PetscBool *per, const PetscReal **maxCell, const PetscReal **L, const DMBoundaryType **bd)
4616: {
4619: if (per) *per = dm->periodic;
4620: if (L) *L = dm->L;
4621: if (maxCell) *maxCell = dm->maxCell;
4622: if (bd) *bd = dm->bdtype;
4623: return(0);
4624: }
4626: /*@C
4627: DMSetPeriodicity - Set the description of mesh periodicity
4629: Input Parameters:
4630: + dm - The DM object
4631: . per - Whether the DM is periodic or not. If maxCell is not provided, coordinates need to be localized
4632: . maxCell - Over distances greater than this, we can assume a point has crossed over to another sheet, when trying to localize cell coordinates
4633: . L - If we assume the mesh is a torus, this is the length of each coordinate
4634: - bd - This describes the type of periodicity in each topological dimension
4636: Level: developer
4638: .seealso: DMGetPeriodicity()
4639: @*/
4640: PetscErrorCode DMSetPeriodicity(DM dm, PetscBool per, const PetscReal maxCell[], const PetscReal L[], const DMBoundaryType bd[])
4641: {
4642: PetscInt dim, d;
4648: if (maxCell) {
4652: }
4653: PetscFree3(dm->L,dm->maxCell,dm->bdtype);
4654: DMGetDimension(dm, &dim);
4655: if (maxCell) {
4656: PetscMalloc3(dim,&dm->L,dim,&dm->maxCell,dim,&dm->bdtype);
4657: for (d = 0; d < dim; ++d) {dm->L[d] = L[d]; dm->maxCell[d] = maxCell[d]; dm->bdtype[d] = bd[d];}
4658: dm->periodic = PETSC_TRUE;
4659: } else {
4660: dm->periodic = per;
4661: }
4662: return(0);
4663: }
4665: /*@
4666: 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.
4668: Input Parameters:
4669: + dm - The DM
4670: . in - The input coordinate point (dim numbers)
4671: - endpoint - Include the endpoint L_i
4673: Output Parameter:
4674: . out - The localized coordinate point
4676: Level: developer
4678: .seealso: DMLocalizeCoordinates(), DMLocalizeAddCoordinate()
4679: @*/
4680: PetscErrorCode DMLocalizeCoordinate(DM dm, const PetscScalar in[], PetscBool endpoint, PetscScalar out[])
4681: {
4682: PetscInt dim, d;
4686: DMGetCoordinateDim(dm, &dim);
4687: if (!dm->maxCell) {
4688: for (d = 0; d < dim; ++d) out[d] = in[d];
4689: } else {
4690: if (endpoint) {
4691: for (d = 0; d < dim; ++d) {
4692: 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)) {
4693: out[d] = in[d] - dm->L[d]*(PetscFloorReal(PetscRealPart(in[d])/dm->L[d]) - 1);
4694: } else {
4695: out[d] = in[d] - dm->L[d]*PetscFloorReal(PetscRealPart(in[d])/dm->L[d]);
4696: }
4697: }
4698: } else {
4699: for (d = 0; d < dim; ++d) {
4700: out[d] = in[d] - dm->L[d]*PetscFloorReal(PetscRealPart(in[d])/dm->L[d]);
4701: }
4702: }
4703: }
4704: return(0);
4705: }
4707: /*
4708: 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.
4710: Input Parameters:
4711: + dm - The DM
4712: . dim - The spatial dimension
4713: . anchor - The anchor point, the input point can be no more than maxCell away from it
4714: - in - The input coordinate point (dim numbers)
4716: Output Parameter:
4717: . out - The localized coordinate point
4719: Level: developer
4721: 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
4723: .seealso: DMLocalizeCoordinates(), DMLocalizeAddCoordinate()
4724: */
4725: PetscErrorCode DMLocalizeCoordinate_Internal(DM dm, PetscInt dim, const PetscScalar anchor[], const PetscScalar in[], PetscScalar out[])
4726: {
4727: PetscInt d;
4730: if (!dm->maxCell) {
4731: for (d = 0; d < dim; ++d) out[d] = in[d];
4732: } else {
4733: for (d = 0; d < dim; ++d) {
4734: if ((dm->bdtype[d] != DM_BOUNDARY_NONE) && (PetscAbsScalar(anchor[d] - in[d]) > dm->maxCell[d])) {
4735: out[d] = PetscRealPart(anchor[d]) > PetscRealPart(in[d]) ? dm->L[d] + in[d] : in[d] - dm->L[d];
4736: } else {
4737: out[d] = in[d];
4738: }
4739: }
4740: }
4741: return(0);
4742: }
4743: PetscErrorCode DMLocalizeCoordinateReal_Internal(DM dm, PetscInt dim, const PetscReal anchor[], const PetscReal in[], PetscReal out[])
4744: {
4745: PetscInt d;
4748: if (!dm->maxCell) {
4749: for (d = 0; d < dim; ++d) out[d] = in[d];
4750: } else {
4751: for (d = 0; d < dim; ++d) {
4752: if ((dm->bdtype[d] != DM_BOUNDARY_NONE) && (PetscAbsReal(anchor[d] - in[d]) > dm->maxCell[d])) {
4753: out[d] = anchor[d] > in[d] ? dm->L[d] + in[d] : in[d] - dm->L[d];
4754: } else {
4755: out[d] = in[d];
4756: }
4757: }
4758: }
4759: return(0);
4760: }
4762: /*
4763: 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.
4765: Input Parameters:
4766: + dm - The DM
4767: . dim - The spatial dimension
4768: . anchor - The anchor point, the input point can be no more than maxCell away from it
4769: . in - The input coordinate delta (dim numbers)
4770: - out - The input coordinate point (dim numbers)
4772: Output Parameter:
4773: . out - The localized coordinate in + out
4775: Level: developer
4777: 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
4779: .seealso: DMLocalizeCoordinates(), DMLocalizeCoordinate()
4780: */
4781: PetscErrorCode DMLocalizeAddCoordinate_Internal(DM dm, PetscInt dim, const PetscScalar anchor[], const PetscScalar in[], PetscScalar out[])
4782: {
4783: PetscInt d;
4786: if (!dm->maxCell) {
4787: for (d = 0; d < dim; ++d) out[d] += in[d];
4788: } else {
4789: for (d = 0; d < dim; ++d) {
4790: if ((dm->bdtype[d] != DM_BOUNDARY_NONE) && (PetscAbsScalar(anchor[d] - in[d]) > dm->maxCell[d])) {
4791: out[d] += PetscRealPart(anchor[d]) > PetscRealPart(in[d]) ? dm->L[d] + in[d] : in[d] - dm->L[d];
4792: } else {
4793: out[d] += in[d];
4794: }
4795: }
4796: }
4797: return(0);
4798: }
4800: /*@
4801: DMGetCoordinatesLocalized - Check if the DM coordinates have been localized for cells
4803: Input Parameter:
4804: . dm - The DM
4806: Output Parameter:
4807: areLocalized - True if localized
4809: Level: developer
4811: .seealso: DMLocalizeCoordinates()
4812: @*/
4813: PetscErrorCode DMGetCoordinatesLocalized(DM dm,PetscBool *areLocalized)
4814: {
4815: DM cdm;
4816: PetscSection coordSection;
4817: PetscInt cStart, cEnd, sStart, sEnd, c, dof;
4818: PetscBool isPlex, alreadyLocalized;
4825: *areLocalized = PETSC_FALSE;
4826: if (!dm->periodic) return(0); /* This is a hideous optimization hack! */
4828: /* We need some generic way of refering to cells/vertices */
4829: DMGetCoordinateDM(dm, &cdm);
4830: DMGetCoordinateSection(dm, &coordSection);
4831: PetscObjectTypeCompare((PetscObject) cdm, DMPLEX, &isPlex);
4832: if (!isPlex) SETERRQ(PetscObjectComm((PetscObject) cdm), PETSC_ERR_ARG_WRONG, "Coordinate localization requires a DMPLEX coordinate DM");
4834: DMPlexGetHeightStratum(cdm, 0, &cStart, &cEnd);
4835: PetscSectionGetChart(coordSection, &sStart, &sEnd);
4836: alreadyLocalized = PETSC_FALSE;
4837: for (c = cStart; c < cEnd; ++c) {
4838: if (c < sStart || c >= sEnd) continue;
4839: PetscSectionGetDof(coordSection, c, &dof);
4840: if (dof) { alreadyLocalized = PETSC_TRUE; break; }
4841: }
4842: MPIU_Allreduce(&alreadyLocalized,areLocalized,1,MPIU_BOOL,MPI_LOR,PetscObjectComm((PetscObject)dm));
4843: return(0);
4844: }
4847: /*@
4848: DMLocalizeCoordinates - If a mesh is periodic, create local coordinates for cells having periodic faces
4850: Input Parameter:
4851: . dm - The DM
4853: Level: developer
4855: .seealso: DMLocalizeCoordinate(), DMLocalizeAddCoordinate()
4856: @*/
4857: PetscErrorCode DMLocalizeCoordinates(DM dm)
4858: {
4859: DM cdm;
4860: PetscSection coordSection, cSection;
4861: Vec coordinates, cVec;
4862: PetscScalar *coords, *coords2, *anchor, *localized;
4863: PetscInt Nc, vStart, vEnd, v, sStart, sEnd, newStart = PETSC_MAX_INT, newEnd = PETSC_MIN_INT, dof, d, off, off2, bs, coordSize;
4864: PetscBool alreadyLocalized, alreadyLocalizedGlobal;
4865: PetscInt maxHeight = 0, h;
4866: PetscInt *pStart = NULL, *pEnd = NULL;
4871: if (!dm->periodic) return(0);
4872: DMGetCoordinatesLocalized(dm, &alreadyLocalized);
4873: if (alreadyLocalized) return(0);
4875: /* We need some generic way of refering to cells/vertices */
4876: DMGetCoordinateDM(dm, &cdm);
4877: {
4878: PetscBool isplex;
4880: PetscObjectTypeCompare((PetscObject) cdm, DMPLEX, &isplex);
4881: if (isplex) {
4882: DMPlexGetDepthStratum(cdm, 0, &vStart, &vEnd);
4883: DMPlexGetMaxProjectionHeight(cdm,&maxHeight);
4884: DMGetWorkArray(dm,2*(maxHeight + 1),MPIU_INT,&pStart);
4885: pEnd = &pStart[maxHeight + 1];
4886: newStart = vStart;
4887: newEnd = vEnd;
4888: for (h = 0; h <= maxHeight; h++) {
4889: DMPlexGetHeightStratum(cdm, h, &pStart[h], &pEnd[h]);
4890: newStart = PetscMin(newStart,pStart[h]);
4891: newEnd = PetscMax(newEnd,pEnd[h]);
4892: }
4893: } else SETERRQ(PetscObjectComm((PetscObject) cdm), PETSC_ERR_ARG_WRONG, "Coordinate localization requires a DMPLEX coordinate DM");
4894: }
4895: DMGetCoordinatesLocal(dm, &coordinates);
4896: if (!coordinates) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Missing local coordinates vector");
4897: DMGetCoordinateSection(dm, &coordSection);
4898: VecGetBlockSize(coordinates, &bs);
4899: PetscSectionGetChart(coordSection,&sStart,&sEnd);
4901: PetscSectionCreate(PetscObjectComm((PetscObject) dm), &cSection);
4902: PetscSectionSetNumFields(cSection, 1);
4903: PetscSectionGetFieldComponents(coordSection, 0, &Nc);
4904: PetscSectionSetFieldComponents(cSection, 0, Nc);
4905: PetscSectionSetChart(cSection, newStart, newEnd);
4907: DMGetWorkArray(dm, 2 * bs, MPIU_SCALAR, &anchor);
4908: localized = &anchor[bs];
4909: alreadyLocalized = alreadyLocalizedGlobal = PETSC_TRUE;
4910: for (h = 0; h <= maxHeight; h++) {
4911: PetscInt cStart = pStart[h], cEnd = pEnd[h], c;
4913: for (c = cStart; c < cEnd; ++c) {
4914: PetscScalar *cellCoords = NULL;
4915: PetscInt b;
4917: if (c < sStart || c >= sEnd) alreadyLocalized = PETSC_FALSE;
4918: DMPlexVecGetClosure(cdm, coordSection, coordinates, c, &dof, &cellCoords);
4919: for (b = 0; b < bs; ++b) anchor[b] = cellCoords[b];
4920: for (d = 0; d < dof/bs; ++d) {
4921: DMLocalizeCoordinate_Internal(dm, bs, anchor, &cellCoords[d*bs], localized);
4922: for (b = 0; b < bs; b++) {
4923: if (cellCoords[d*bs + b] != localized[b]) break;
4924: }
4925: if (b < bs) break;
4926: }
4927: if (d < dof/bs) {
4928: if (c >= sStart && c < sEnd) {
4929: PetscInt cdof;
4931: PetscSectionGetDof(coordSection, c, &cdof);
4932: if (cdof != dof) alreadyLocalized = PETSC_FALSE;
4933: }
4934: PetscSectionSetDof(cSection, c, dof);
4935: PetscSectionSetFieldDof(cSection, c, 0, dof);
4936: }
4937: DMPlexVecRestoreClosure(cdm, coordSection, coordinates, c, &dof, &cellCoords);
4938: }
4939: }
4940: MPI_Allreduce(&alreadyLocalized,&alreadyLocalizedGlobal,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)dm));
4941: if (alreadyLocalizedGlobal) {
4942: DMRestoreWorkArray(dm, 2 * bs, MPIU_SCALAR, &anchor);
4943: PetscSectionDestroy(&cSection);
4944: DMRestoreWorkArray(dm,2*(maxHeight + 1),MPIU_INT,&pStart);
4945: return(0);
4946: }
4947: for (v = vStart; v < vEnd; ++v) {
4948: PetscSectionGetDof(coordSection, v, &dof);
4949: PetscSectionSetDof(cSection, v, dof);
4950: PetscSectionSetFieldDof(cSection, v, 0, dof);
4951: }
4952: PetscSectionSetUp(cSection);
4953: PetscSectionGetStorageSize(cSection, &coordSize);
4954: VecCreate(PETSC_COMM_SELF, &cVec);
4955: PetscObjectSetName((PetscObject)cVec,"coordinates");
4956: VecSetBlockSize(cVec, bs);
4957: VecSetSizes(cVec, coordSize, PETSC_DETERMINE);
4958: VecSetType(cVec, VECSTANDARD);
4959: VecGetArrayRead(coordinates, (const PetscScalar**)&coords);
4960: VecGetArray(cVec, &coords2);
4961: for (v = vStart; v < vEnd; ++v) {
4962: PetscSectionGetDof(coordSection, v, &dof);
4963: PetscSectionGetOffset(coordSection, v, &off);
4964: PetscSectionGetOffset(cSection, v, &off2);
4965: for (d = 0; d < dof; ++d) coords2[off2+d] = coords[off+d];
4966: }
4967: for (h = 0; h <= maxHeight; h++) {
4968: PetscInt cStart = pStart[h], cEnd = pEnd[h], c;
4970: for (c = cStart; c < cEnd; ++c) {
4971: PetscScalar *cellCoords = NULL;
4972: PetscInt b, cdof;
4974: PetscSectionGetDof(cSection,c,&cdof);
4975: if (!cdof) continue;
4976: DMPlexVecGetClosure(cdm, coordSection, coordinates, c, &dof, &cellCoords);
4977: PetscSectionGetOffset(cSection, c, &off2);
4978: for (b = 0; b < bs; ++b) anchor[b] = cellCoords[b];
4979: for (d = 0; d < dof/bs; ++d) {DMLocalizeCoordinate_Internal(dm, bs, anchor, &cellCoords[d*bs], &coords2[off2+d*bs]);}
4980: DMPlexVecRestoreClosure(cdm, coordSection, coordinates, c, &dof, &cellCoords);
4981: }
4982: }
4983: DMRestoreWorkArray(dm, 2 * bs, MPIU_SCALAR, &anchor);
4984: DMRestoreWorkArray(dm,2*(maxHeight + 1),MPIU_INT,&pStart);
4985: VecRestoreArrayRead(coordinates, (const PetscScalar**)&coords);
4986: VecRestoreArray(cVec, &coords2);
4987: DMSetCoordinateSection(dm, PETSC_DETERMINE, cSection);
4988: DMSetCoordinatesLocal(dm, cVec);
4989: VecDestroy(&cVec);
4990: PetscSectionDestroy(&cSection);
4991: return(0);
4992: }
4994: /*@
4995: DMLocatePoints - Locate the points in v in the mesh and return a PetscSF of the containing cells
4997: Collective on Vec v (see explanation below)
4999: Input Parameters:
5000: + dm - The DM
5001: . v - The Vec of points
5002: . ltype - The type of point location, e.g. DM_POINTLOCATION_NONE or DM_POINTLOCATION_NEAREST
5003: - cells - Points to either NULL, or a PetscSF with guesses for which cells contain each point.
5005: Output Parameter:
5006: + v - The Vec of points, which now contains the nearest mesh points to the given points if DM_POINTLOCATION_NEAREST is used
5007: - cells - The PetscSF containing the ranks and local indices of the containing points.
5010: Level: developer
5012: Notes:
5013: To do a search of the local cells of the mesh, v should have PETSC_COMM_SELF as its communicator.
5014: To do a search of all the cells in the distributed mesh, v should have the same communicator as dm.
5016: If *cellSF is NULL on input, a PetscSF will be created.
5017: If *cellSF is not NULL on input, it should point to an existing PetscSF, whose graph will be used as initial guesses.
5019: An array that maps each point to its containing cell can be obtained with
5021: $ const PetscSFNode *cells;
5022: $ PetscInt nFound;
5023: $ const PetscSFNode *found;
5024: $
5025: $ PetscSFGetGraph(cells,NULL,&nFound,&found,&cells);
5027: Where cells[i].rank is the rank of the cell containing point found[i] (or i if found == NULL), and cells[i].index is
5028: the index of the cell in its rank's local numbering.
5030: .keywords: point location, mesh
5031: .seealso: DMSetCoordinates(), DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLocal(), DMPointLocationType
5032: @*/
5033: PetscErrorCode DMLocatePoints(DM dm, Vec v, DMPointLocationType ltype, PetscSF *cellSF)
5034: {
5041: if (*cellSF) {
5042: PetscMPIInt result;
5045: MPI_Comm_compare(PetscObjectComm((PetscObject)v),PetscObjectComm((PetscObject)*cellSF),&result);
5046: if (result != MPI_IDENT && result != MPI_CONGRUENT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"cellSF must have a communicator congruent to v's");
5047: } else {
5048: PetscSFCreate(PetscObjectComm((PetscObject)v),cellSF);
5049: }
5050: PetscLogEventBegin(DM_LocatePoints,dm,0,0,0);
5051: if (dm->ops->locatepoints) {
5052: (*dm->ops->locatepoints)(dm,v,ltype,*cellSF);
5053: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Point location not available for this DM");
5054: PetscLogEventEnd(DM_LocatePoints,dm,0,0,0);
5055: return(0);
5056: }
5058: /*@
5059: DMGetOutputDM - Retrieve the DM associated with the layout for output
5061: Input Parameter:
5062: . dm - The original DM
5064: Output Parameter:
5065: . odm - The DM which provides the layout for output
5067: Level: intermediate
5069: .seealso: VecView(), DMGetDefaultGlobalSection()
5070: @*/
5071: PetscErrorCode DMGetOutputDM(DM dm, DM *odm)
5072: {
5073: PetscSection section;
5074: PetscBool hasConstraints, ghasConstraints;
5080: DMGetDefaultSection(dm, §ion);
5081: PetscSectionHasConstraints(section, &hasConstraints);
5082: MPI_Allreduce(&hasConstraints, &ghasConstraints, 1, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject) dm));
5083: if (!ghasConstraints) {
5084: *odm = dm;
5085: return(0);
5086: }
5087: if (!dm->dmBC) {
5088: PetscDS ds;
5089: PetscSection newSection, gsection;
5090: PetscSF sf;
5092: DMClone(dm, &dm->dmBC);
5093: DMGetDS(dm, &ds);
5094: DMSetDS(dm->dmBC, ds);
5095: PetscSectionClone(section, &newSection);
5096: DMSetDefaultSection(dm->dmBC, newSection);
5097: PetscSectionDestroy(&newSection);
5098: DMGetPointSF(dm->dmBC, &sf);
5099: PetscSectionCreateGlobalSection(section, sf, PETSC_TRUE, PETSC_FALSE, &gsection);
5100: DMSetDefaultGlobalSection(dm->dmBC, gsection);
5101: PetscSectionDestroy(&gsection);
5102: }
5103: *odm = dm->dmBC;
5104: return(0);
5105: }
5107: /*@
5108: DMGetOutputSequenceNumber - Retrieve the sequence number/value for output
5110: Input Parameter:
5111: . dm - The original DM
5113: Output Parameters:
5114: + num - The output sequence number
5115: - val - The output sequence value
5117: Level: intermediate
5119: Note: This is intended for output that should appear in sequence, for instance
5120: a set of timesteps in an HDF5 file, or a set of realizations of a stochastic system.
5122: .seealso: VecView()
5123: @*/
5124: PetscErrorCode DMGetOutputSequenceNumber(DM dm, PetscInt *num, PetscReal *val)
5125: {
5130: return(0);
5131: }
5133: /*@
5134: DMSetOutputSequenceNumber - Set the sequence number/value for output
5136: Input Parameters:
5137: + dm - The original DM
5138: . num - The output sequence number
5139: - val - The output sequence value
5141: Level: intermediate
5143: Note: This is intended for output that should appear in sequence, for instance
5144: a set of timesteps in an HDF5 file, or a set of realizations of a stochastic system.
5146: .seealso: VecView()
5147: @*/
5148: PetscErrorCode DMSetOutputSequenceNumber(DM dm, PetscInt num, PetscReal val)
5149: {
5152: dm->outputSequenceNum = num;
5153: dm->outputSequenceVal = val;
5154: return(0);
5155: }
5157: /*@C
5158: DMOutputSequenceLoad - Retrieve the sequence value from a Viewer
5160: Input Parameters:
5161: + dm - The original DM
5162: . name - The sequence name
5163: - num - The output sequence number
5165: Output Parameter:
5166: . val - The output sequence value
5168: Level: intermediate
5170: Note: This is intended for output that should appear in sequence, for instance
5171: a set of timesteps in an HDF5 file, or a set of realizations of a stochastic system.
5173: .seealso: DMGetOutputSequenceNumber(), DMSetOutputSequenceNumber(), VecView()
5174: @*/
5175: PetscErrorCode DMOutputSequenceLoad(DM dm, PetscViewer viewer, const char *name, PetscInt num, PetscReal *val)
5176: {
5177: PetscBool ishdf5;
5184: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERHDF5, &ishdf5);
5185: if (ishdf5) {
5186: #if defined(PETSC_HAVE_HDF5)
5187: PetscScalar value;
5189: DMSequenceLoad_HDF5_Internal(dm, name, num, &value, viewer);
5190: *val = PetscRealPart(value);
5191: #endif
5192: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid viewer; open viewer with PetscViewerHDF5Open()");
5193: return(0);
5194: }
5196: /*@
5197: DMGetUseNatural - Get the flag for creating a mapping to the natural order on distribution
5199: Not collective
5201: Input Parameter:
5202: . dm - The DM
5204: Output Parameter:
5205: . useNatural - The flag to build the mapping to a natural order during distribution
5207: Level: beginner
5209: .seealso: DMSetUseNatural(), DMCreate()
5210: @*/
5211: PetscErrorCode DMGetUseNatural(DM dm, PetscBool *useNatural)
5212: {
5216: *useNatural = dm->useNatural;
5217: return(0);
5218: }
5220: /*@
5221: DMSetUseNatural - Set the flag for creating a mapping to the natural order on distribution
5223: Collective on dm
5225: Input Parameters:
5226: + dm - The DM
5227: - useNatural - The flag to build the mapping to a natural order during distribution
5229: Level: beginner
5231: .seealso: DMGetUseNatural(), DMCreate()
5232: @*/
5233: PetscErrorCode DMSetUseNatural(DM dm, PetscBool useNatural)
5234: {
5238: dm->useNatural = useNatural;
5239: return(0);
5240: }
5243: /*@C
5244: DMCreateLabel - Create a label of the given name if it does not already exist
5246: Not Collective
5248: Input Parameters:
5249: + dm - The DM object
5250: - name - The label name
5252: Level: intermediate
5254: .keywords: mesh
5255: .seealso: DMLabelCreate(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5256: @*/
5257: PetscErrorCode DMCreateLabel(DM dm, const char name[])
5258: {
5259: DMLabelLink next = dm->labels->next;
5260: PetscBool flg = PETSC_FALSE;
5266: while (next) {
5267: PetscStrcmp(name, next->label->name, &flg);
5268: if (flg) break;
5269: next = next->next;
5270: }
5271: if (!flg) {
5272: DMLabelLink tmpLabel;
5274: PetscCalloc1(1, &tmpLabel);
5275: DMLabelCreate(name, &tmpLabel->label);
5276: tmpLabel->output = PETSC_TRUE;
5277: tmpLabel->next = dm->labels->next;
5278: dm->labels->next = tmpLabel;
5279: }
5280: return(0);
5281: }
5283: /*@C
5284: DMGetLabelValue - Get the value in a Sieve Label for the given point, with 0 as the default
5286: Not Collective
5288: Input Parameters:
5289: + dm - The DM object
5290: . name - The label name
5291: - point - The mesh point
5293: Output Parameter:
5294: . value - The label value for this point, or -1 if the point is not in the label
5296: Level: beginner
5298: .keywords: mesh
5299: .seealso: DMLabelGetValue(), DMSetLabelValue(), DMGetStratumIS()
5300: @*/
5301: PetscErrorCode DMGetLabelValue(DM dm, const char name[], PetscInt point, PetscInt *value)
5302: {
5303: DMLabel label;
5309: DMGetLabel(dm, name, &label);
5310: if (!label) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No label named %s was found", name);
5311: DMLabelGetValue(label, point, value);
5312: return(0);
5313: }
5315: /*@C
5316: DMSetLabelValue - Add a point to a Sieve Label with given value
5318: Not Collective
5320: Input Parameters:
5321: + dm - The DM object
5322: . name - The label name
5323: . point - The mesh point
5324: - value - The label value for this point
5326: Output Parameter:
5328: Level: beginner
5330: .keywords: mesh
5331: .seealso: DMLabelSetValue(), DMGetStratumIS(), DMClearLabelValue()
5332: @*/
5333: PetscErrorCode DMSetLabelValue(DM dm, const char name[], PetscInt point, PetscInt value)
5334: {
5335: DMLabel label;
5341: DMGetLabel(dm, name, &label);
5342: if (!label) {
5343: DMCreateLabel(dm, name);
5344: DMGetLabel(dm, name, &label);
5345: }
5346: DMLabelSetValue(label, point, value);
5347: return(0);
5348: }
5350: /*@C
5351: DMClearLabelValue - Remove a point from a Sieve Label with given value
5353: Not Collective
5355: Input Parameters:
5356: + dm - The DM object
5357: . name - The label name
5358: . point - The mesh point
5359: - value - The label value for this point
5361: Output Parameter:
5363: Level: beginner
5365: .keywords: mesh
5366: .seealso: DMLabelClearValue(), DMSetLabelValue(), DMGetStratumIS()
5367: @*/
5368: PetscErrorCode DMClearLabelValue(DM dm, const char name[], PetscInt point, PetscInt value)
5369: {
5370: DMLabel label;
5376: DMGetLabel(dm, name, &label);
5377: if (!label) return(0);
5378: DMLabelClearValue(label, point, value);
5379: return(0);
5380: }
5382: /*@C
5383: DMGetLabelSize - Get the number of different integer ids in a Label
5385: Not Collective
5387: Input Parameters:
5388: + dm - The DM object
5389: - name - The label name
5391: Output Parameter:
5392: . size - The number of different integer ids, or 0 if the label does not exist
5394: Level: beginner
5396: .keywords: mesh
5397: .seealso: DMLabeGetNumValues(), DMSetLabelValue()
5398: @*/
5399: PetscErrorCode DMGetLabelSize(DM dm, const char name[], PetscInt *size)
5400: {
5401: DMLabel label;
5408: DMGetLabel(dm, name, &label);
5409: *size = 0;
5410: if (!label) return(0);
5411: DMLabelGetNumValues(label, size);
5412: return(0);
5413: }
5415: /*@C
5416: DMGetLabelIdIS - Get the integer ids in a label
5418: Not Collective
5420: Input Parameters:
5421: + mesh - The DM object
5422: - name - The label name
5424: Output Parameter:
5425: . ids - The integer ids, or NULL if the label does not exist
5427: Level: beginner
5429: .keywords: mesh
5430: .seealso: DMLabelGetValueIS(), DMGetLabelSize()
5431: @*/
5432: PetscErrorCode DMGetLabelIdIS(DM dm, const char name[], IS *ids)
5433: {
5434: DMLabel label;
5441: DMGetLabel(dm, name, &label);
5442: *ids = NULL;
5443: if (label) {
5444: DMLabelGetValueIS(label, ids);
5445: } else {
5446: /* returning an empty IS */
5447: ISCreateGeneral(PETSC_COMM_SELF,0,NULL,PETSC_USE_POINTER,ids);
5448: }
5449: return(0);
5450: }
5452: /*@C
5453: DMGetStratumSize - Get the number of points in a label stratum
5455: Not Collective
5457: Input Parameters:
5458: + dm - The DM object
5459: . name - The label name
5460: - value - The stratum value
5462: Output Parameter:
5463: . size - The stratum size
5465: Level: beginner
5467: .keywords: mesh
5468: .seealso: DMLabelGetStratumSize(), DMGetLabelSize(), DMGetLabelIds()
5469: @*/
5470: PetscErrorCode DMGetStratumSize(DM dm, const char name[], PetscInt value, PetscInt *size)
5471: {
5472: DMLabel label;
5479: DMGetLabel(dm, name, &label);
5480: *size = 0;
5481: if (!label) return(0);
5482: DMLabelGetStratumSize(label, value, size);
5483: return(0);
5484: }
5486: /*@C
5487: DMGetStratumIS - Get the points in a label stratum
5489: Not Collective
5491: Input Parameters:
5492: + dm - The DM object
5493: . name - The label name
5494: - value - The stratum value
5496: Output Parameter:
5497: . points - The stratum points, or NULL if the label does not exist or does not have that value
5499: Level: beginner
5501: .keywords: mesh
5502: .seealso: DMLabelGetStratumIS(), DMGetStratumSize()
5503: @*/
5504: PetscErrorCode DMGetStratumIS(DM dm, const char name[], PetscInt value, IS *points)
5505: {
5506: DMLabel label;
5513: DMGetLabel(dm, name, &label);
5514: *points = NULL;
5515: if (!label) return(0);
5516: DMLabelGetStratumIS(label, value, points);
5517: return(0);
5518: }
5520: /*@C
5521: DMGetStratumIS - Set the points in a label stratum
5523: Not Collective
5525: Input Parameters:
5526: + dm - The DM object
5527: . name - The label name
5528: . value - The stratum value
5529: - points - The stratum points
5531: Level: beginner
5533: .keywords: mesh
5534: .seealso: DMLabelSetStratumIS(), DMGetStratumSize()
5535: @*/
5536: PetscErrorCode DMSetStratumIS(DM dm, const char name[], PetscInt value, IS points)
5537: {
5538: DMLabel label;
5545: DMGetLabel(dm, name, &label);
5546: if (!label) return(0);
5547: DMLabelSetStratumIS(label, value, points);
5548: return(0);
5549: }
5551: /*@C
5552: DMClearLabelStratum - Remove all points from a stratum from a Sieve Label
5554: Not Collective
5556: Input Parameters:
5557: + dm - The DM object
5558: . name - The label name
5559: - value - The label value for this point
5561: Output Parameter:
5563: Level: beginner
5565: .keywords: mesh
5566: .seealso: DMLabelClearStratum(), DMSetLabelValue(), DMGetStratumIS(), DMClearLabelValue()
5567: @*/
5568: PetscErrorCode DMClearLabelStratum(DM dm, const char name[], PetscInt value)
5569: {
5570: DMLabel label;
5576: DMGetLabel(dm, name, &label);
5577: if (!label) return(0);
5578: DMLabelClearStratum(label, value);
5579: return(0);
5580: }
5582: /*@
5583: DMGetNumLabels - Return the number of labels defined by the mesh
5585: Not Collective
5587: Input Parameter:
5588: . dm - The DM object
5590: Output Parameter:
5591: . numLabels - the number of Labels
5593: Level: intermediate
5595: .keywords: mesh
5596: .seealso: DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5597: @*/
5598: PetscErrorCode DMGetNumLabels(DM dm, PetscInt *numLabels)
5599: {
5600: DMLabelLink next = dm->labels->next;
5601: PetscInt n = 0;
5606: while (next) {++n; next = next->next;}
5607: *numLabels = n;
5608: return(0);
5609: }
5611: /*@C
5612: DMGetLabelName - Return the name of nth label
5614: Not Collective
5616: Input Parameters:
5617: + dm - The DM object
5618: - n - the label number
5620: Output Parameter:
5621: . name - the label name
5623: Level: intermediate
5625: .keywords: mesh
5626: .seealso: DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5627: @*/
5628: PetscErrorCode DMGetLabelName(DM dm, PetscInt n, const char **name)
5629: {
5630: DMLabelLink next = dm->labels->next;
5631: PetscInt l = 0;
5636: while (next) {
5637: if (l == n) {
5638: *name = next->label->name;
5639: return(0);
5640: }
5641: ++l;
5642: next = next->next;
5643: }
5644: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label %D does not exist in this DM", n);
5645: }
5647: /*@C
5648: DMHasLabel - Determine whether the mesh has a label of a given name
5650: Not Collective
5652: Input Parameters:
5653: + dm - The DM object
5654: - name - The label name
5656: Output Parameter:
5657: . hasLabel - PETSC_TRUE if the label is present
5659: Level: intermediate
5661: .keywords: mesh
5662: .seealso: DMCreateLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5663: @*/
5664: PetscErrorCode DMHasLabel(DM dm, const char name[], PetscBool *hasLabel)
5665: {
5666: DMLabelLink next = dm->labels->next;
5673: *hasLabel = PETSC_FALSE;
5674: while (next) {
5675: PetscStrcmp(name, next->label->name, hasLabel);
5676: if (*hasLabel) break;
5677: next = next->next;
5678: }
5679: return(0);
5680: }
5682: /*@C
5683: DMGetLabel - Return the label of a given name, or NULL
5685: Not Collective
5687: Input Parameters:
5688: + dm - The DM object
5689: - name - The label name
5691: Output Parameter:
5692: . label - The DMLabel, or NULL if the label is absent
5694: Level: intermediate
5696: .keywords: mesh
5697: .seealso: DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5698: @*/
5699: PetscErrorCode DMGetLabel(DM dm, const char name[], DMLabel *label)
5700: {
5701: DMLabelLink next = dm->labels->next;
5702: PetscBool hasLabel;
5709: *label = NULL;
5710: while (next) {
5711: PetscStrcmp(name, next->label->name, &hasLabel);
5712: if (hasLabel) {
5713: *label = next->label;
5714: break;
5715: }
5716: next = next->next;
5717: }
5718: return(0);
5719: }
5721: /*@C
5722: DMGetLabelByNum - Return the nth label
5724: Not Collective
5726: Input Parameters:
5727: + dm - The DM object
5728: - n - the label number
5730: Output Parameter:
5731: . label - the label
5733: Level: intermediate
5735: .keywords: mesh
5736: .seealso: DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5737: @*/
5738: PetscErrorCode DMGetLabelByNum(DM dm, PetscInt n, DMLabel *label)
5739: {
5740: DMLabelLink next = dm->labels->next;
5741: PetscInt l = 0;
5746: while (next) {
5747: if (l == n) {
5748: *label = next->label;
5749: return(0);
5750: }
5751: ++l;
5752: next = next->next;
5753: }
5754: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label %D does not exist in this DM", n);
5755: }
5757: /*@C
5758: DMAddLabel - Add the label to this mesh
5760: Not Collective
5762: Input Parameters:
5763: + dm - The DM object
5764: - label - The DMLabel
5766: Level: developer
5768: .keywords: mesh
5769: .seealso: DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5770: @*/
5771: PetscErrorCode DMAddLabel(DM dm, DMLabel label)
5772: {
5773: DMLabelLink tmpLabel;
5774: PetscBool hasLabel;
5779: DMHasLabel(dm, label->name, &hasLabel);
5780: if (hasLabel) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label %s already exists in this DM", label->name);
5781: PetscCalloc1(1, &tmpLabel);
5782: tmpLabel->label = label;
5783: tmpLabel->output = PETSC_TRUE;
5784: tmpLabel->next = dm->labels->next;
5785: dm->labels->next = tmpLabel;
5786: return(0);
5787: }
5789: /*@C
5790: DMRemoveLabel - Remove the label from this mesh
5792: Not Collective
5794: Input Parameters:
5795: + dm - The DM object
5796: - name - The label name
5798: Output Parameter:
5799: . label - The DMLabel, or NULL if the label is absent
5801: Level: developer
5803: .keywords: mesh
5804: .seealso: DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5805: @*/
5806: PetscErrorCode DMRemoveLabel(DM dm, const char name[], DMLabel *label)
5807: {
5808: DMLabelLink next = dm->labels->next;
5809: DMLabelLink last = NULL;
5810: PetscBool hasLabel;
5815: DMHasLabel(dm, name, &hasLabel);
5816: *label = NULL;
5817: if (!hasLabel) return(0);
5818: while (next) {
5819: PetscStrcmp(name, next->label->name, &hasLabel);
5820: if (hasLabel) {
5821: if (last) last->next = next->next;
5822: else dm->labels->next = next->next;
5823: next->next = NULL;
5824: *label = next->label;
5825: PetscStrcmp(name, "depth", &hasLabel);
5826: if (hasLabel) {
5827: dm->depthLabel = NULL;
5828: }
5829: PetscFree(next);
5830: break;
5831: }
5832: last = next;
5833: next = next->next;
5834: }
5835: return(0);
5836: }
5838: /*@C
5839: DMGetLabelOutput - Get the output flag for a given label
5841: Not Collective
5843: Input Parameters:
5844: + dm - The DM object
5845: - name - The label name
5847: Output Parameter:
5848: . output - The flag for output
5850: Level: developer
5852: .keywords: mesh
5853: .seealso: DMSetLabelOutput(), DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5854: @*/
5855: PetscErrorCode DMGetLabelOutput(DM dm, const char name[], PetscBool *output)
5856: {
5857: DMLabelLink next = dm->labels->next;
5864: while (next) {
5865: PetscBool flg;
5867: PetscStrcmp(name, next->label->name, &flg);
5868: if (flg) {*output = next->output; return(0);}
5869: next = next->next;
5870: }
5871: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No label named %s was present in this dm", name);
5872: }
5874: /*@C
5875: DMSetLabelOutput - Set the output flag for a given label
5877: Not Collective
5879: Input Parameters:
5880: + dm - The DM object
5881: . name - The label name
5882: - output - The flag for output
5884: Level: developer
5886: .keywords: mesh
5887: .seealso: DMGetLabelOutput(), DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5888: @*/
5889: PetscErrorCode DMSetLabelOutput(DM dm, const char name[], PetscBool output)
5890: {
5891: DMLabelLink next = dm->labels->next;
5897: while (next) {
5898: PetscBool flg;
5900: PetscStrcmp(name, next->label->name, &flg);
5901: if (flg) {next->output = output; return(0);}
5902: next = next->next;
5903: }
5904: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No label named %s was present in this dm", name);
5905: }
5908: /*@
5909: DMCopyLabels - Copy labels from one mesh to another with a superset of the points
5911: Collective on DM
5913: Input Parameter:
5914: . dmA - The DM object with initial labels
5916: Output Parameter:
5917: . dmB - The DM object with copied labels
5919: Level: intermediate
5921: Note: This is typically used when interpolating or otherwise adding to a mesh
5923: .keywords: mesh
5924: .seealso: DMCopyCoordinates(), DMGetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM(), DMGetCoordinateSection()
5925: @*/
5926: PetscErrorCode DMCopyLabels(DM dmA, DM dmB)
5927: {
5928: PetscInt numLabels, l;
5932: if (dmA == dmB) return(0);
5933: DMGetNumLabels(dmA, &numLabels);
5934: for (l = 0; l < numLabels; ++l) {
5935: DMLabel label, labelNew;
5936: const char *name;
5937: PetscBool flg;
5939: DMGetLabelName(dmA, l, &name);
5940: PetscStrcmp(name, "depth", &flg);
5941: if (flg) continue;
5942: DMGetLabel(dmA, name, &label);
5943: DMLabelDuplicate(label, &labelNew);
5944: DMAddLabel(dmB, labelNew);
5945: }
5946: return(0);
5947: }
5949: /*@
5950: DMGetCoarseDM - Get the coarse mesh from which this was obtained by refinement
5952: Input Parameter:
5953: . dm - The DM object
5955: Output Parameter:
5956: . cdm - The coarse DM
5958: Level: intermediate
5960: .seealso: DMSetCoarseDM()
5961: @*/
5962: PetscErrorCode DMGetCoarseDM(DM dm, DM *cdm)
5963: {
5967: *cdm = dm->coarseMesh;
5968: return(0);
5969: }
5971: /*@
5972: DMSetCoarseDM - Set the coarse mesh from which this was obtained by refinement
5974: Input Parameters:
5975: + dm - The DM object
5976: - cdm - The coarse DM
5978: Level: intermediate
5980: .seealso: DMGetCoarseDM()
5981: @*/
5982: PetscErrorCode DMSetCoarseDM(DM dm, DM cdm)
5983: {
5989: PetscObjectReference((PetscObject)cdm);
5990: DMDestroy(&dm->coarseMesh);
5991: dm->coarseMesh = cdm;
5992: return(0);
5993: }
5995: /*@
5996: DMGetFineDM - Get the fine mesh from which this was obtained by refinement
5998: Input Parameter:
5999: . dm - The DM object
6001: Output Parameter:
6002: . fdm - The fine DM
6004: Level: intermediate
6006: .seealso: DMSetFineDM()
6007: @*/
6008: PetscErrorCode DMGetFineDM(DM dm, DM *fdm)
6009: {
6013: *fdm = dm->fineMesh;
6014: return(0);
6015: }
6017: /*@
6018: DMSetFineDM - Set the fine mesh from which this was obtained by refinement
6020: Input Parameters:
6021: + dm - The DM object
6022: - fdm - The fine DM
6024: Level: intermediate
6026: .seealso: DMGetFineDM()
6027: @*/
6028: PetscErrorCode DMSetFineDM(DM dm, DM fdm)
6029: {
6035: PetscObjectReference((PetscObject)fdm);
6036: DMDestroy(&dm->fineMesh);
6037: dm->fineMesh = fdm;
6038: return(0);
6039: }
6041: /*=== DMBoundary code ===*/
6043: PetscErrorCode DMCopyBoundary(DM dm, DM dmNew)
6044: {
6048: PetscDSCopyBoundary(dm->prob,dmNew->prob);
6049: return(0);
6050: }
6052: /*@C
6053: DMAddBoundary - Add a boundary condition to the model
6055: Input Parameters:
6056: + dm - The DM, with a PetscDS that matches the problem being constrained
6057: . type - The type of condition, e.g. DM_BC_ESSENTIAL_ANALYTIC/DM_BC_ESSENTIAL_FIELD (Dirichlet), or DM_BC_NATURAL (Neumann)
6058: . name - The BC name
6059: . labelname - The label defining constrained points
6060: . field - The field to constrain
6061: . numcomps - The number of constrained field components
6062: . comps - An array of constrained component numbers
6063: . bcFunc - A pointwise function giving boundary values
6064: . numids - The number of DMLabel ids for constrained points
6065: . ids - An array of ids for constrained points
6066: - ctx - An optional user context for bcFunc
6068: Options Database Keys:
6069: + -bc_<boundary name> <num> - Overrides the boundary ids
6070: - -bc_<boundary name>_comp <num> - Overrides the boundary components
6072: Level: developer
6074: .seealso: DMGetBoundary()
6075: @*/
6076: 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)
6077: {
6082: PetscDSAddBoundary(dm->prob,type,name,labelname,field,numcomps,comps,bcFunc,numids,ids,ctx);
6083: return(0);
6084: }
6086: /*@
6087: DMGetNumBoundary - Get the number of registered BC
6089: Input Parameters:
6090: . dm - The mesh object
6092: Output Parameters:
6093: . numBd - The number of BC
6095: Level: intermediate
6097: .seealso: DMAddBoundary(), DMGetBoundary()
6098: @*/
6099: PetscErrorCode DMGetNumBoundary(DM dm, PetscInt *numBd)
6100: {
6105: PetscDSGetNumBoundary(dm->prob,numBd);
6106: return(0);
6107: }
6109: /*@C
6110: DMGetBoundary - Get a model boundary condition
6112: Input Parameters:
6113: + dm - The mesh object
6114: - bd - The BC number
6116: Output Parameters:
6117: + type - The type of condition, e.g. DM_BC_ESSENTIAL_ANALYTIC/DM_BC_ESSENTIAL_FIELD (Dirichlet), or DM_BC_NATURAL (Neumann)
6118: . name - The BC name
6119: . labelname - The label defining constrained points
6120: . field - The field to constrain
6121: . numcomps - The number of constrained field components
6122: . comps - An array of constrained component numbers
6123: . bcFunc - A pointwise function giving boundary values
6124: . numids - The number of DMLabel ids for constrained points
6125: . ids - An array of ids for constrained points
6126: - ctx - An optional user context for bcFunc
6128: Options Database Keys:
6129: + -bc_<boundary name> <num> - Overrides the boundary ids
6130: - -bc_<boundary name>_comp <num> - Overrides the boundary components
6132: Level: developer
6134: .seealso: DMAddBoundary()
6135: @*/
6136: 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)
6137: {
6142: PetscDSGetBoundary(dm->prob,bd,type,name,labelname,field,numcomps,comps,func,numids,ids,ctx);
6143: return(0);
6144: }
6146: static PetscErrorCode DMPopulateBoundary(DM dm)
6147: {
6148: DMBoundary *lastnext;
6149: DSBoundary dsbound;
6153: dsbound = dm->prob->boundary;
6154: if (dm->boundary) {
6155: DMBoundary next = dm->boundary;
6157: /* quick check to see if the PetscDS has changed */
6158: if (next->dsboundary == dsbound) return(0);
6159: /* the PetscDS has changed: tear down and rebuild */
6160: while (next) {
6161: DMBoundary b = next;
6163: next = b->next;
6164: PetscFree(b);
6165: }
6166: dm->boundary = NULL;
6167: }
6169: lastnext = &(dm->boundary);
6170: while (dsbound) {
6171: DMBoundary dmbound;
6173: PetscNew(&dmbound);
6174: dmbound->dsboundary = dsbound;
6175: DMGetLabel(dm, dsbound->labelname, &(dmbound->label));
6176: if (!dmbound->label) PetscInfo2(dm, "DSBoundary %s wants label %s, which is not in this dm.\n",dsbound->name,dsbound->labelname);
6177: /* push on the back instead of the front so that it is in the same order as in the PetscDS */
6178: *lastnext = dmbound;
6179: lastnext = &(dmbound->next);
6180: dsbound = dsbound->next;
6181: }
6182: return(0);
6183: }
6185: PetscErrorCode DMIsBoundaryPoint(DM dm, PetscInt point, PetscBool *isBd)
6186: {
6187: DMBoundary b;
6193: *isBd = PETSC_FALSE;
6194: DMPopulateBoundary(dm);
6195: b = dm->boundary;
6196: while (b && !(*isBd)) {
6197: DMLabel label = b->label;
6198: DSBoundary dsb = b->dsboundary;
6200: if (label) {
6201: PetscInt i;
6203: for (i = 0; i < dsb->numids && !(*isBd); ++i) {
6204: DMLabelStratumHasPoint(label, dsb->ids[i], point, isBd);
6205: }
6206: }
6207: b = b->next;
6208: }
6209: return(0);
6210: }
6212: /*@C
6213: DMProjectFunction - This projects the given function into the function space provided.
6215: Input Parameters:
6216: + dm - The DM
6217: . time - The time
6218: . funcs - The coordinate functions to evaluate, one per field
6219: . ctxs - Optional array of contexts to pass to each coordinate function. ctxs itself may be null.
6220: - mode - The insertion mode for values
6222: Output Parameter:
6223: . X - vector
6225: Calling sequence of func:
6226: $ func(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar u[], void *ctx);
6228: + dim - The spatial dimension
6229: . x - The coordinates
6230: . Nf - The number of fields
6231: . u - The output field values
6232: - ctx - optional user-defined function context
6234: Level: developer
6236: .seealso: DMComputeL2Diff()
6237: @*/
6238: PetscErrorCode DMProjectFunction(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec X)
6239: {
6240: Vec localX;
6245: DMGetLocalVector(dm, &localX);
6246: DMProjectFunctionLocal(dm, time, funcs, ctxs, mode, localX);
6247: DMLocalToGlobalBegin(dm, localX, mode, X);
6248: DMLocalToGlobalEnd(dm, localX, mode, X);
6249: DMRestoreLocalVector(dm, &localX);
6250: return(0);
6251: }
6253: PetscErrorCode DMProjectFunctionLocal(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX)
6254: {
6260: if (!dm->ops->projectfunctionlocal) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMProjectFunctionLocal",((PetscObject)dm)->type_name);
6261: (dm->ops->projectfunctionlocal) (dm, time, funcs, ctxs, mode, localX);
6262: return(0);
6263: }
6265: 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)
6266: {
6267: Vec localX;
6272: DMGetLocalVector(dm, &localX);
6273: DMProjectFunctionLabelLocal(dm, time, label, numIds, ids, Nc, comps, funcs, ctxs, mode, localX);
6274: DMLocalToGlobalBegin(dm, localX, mode, X);
6275: DMLocalToGlobalEnd(dm, localX, mode, X);
6276: DMRestoreLocalVector(dm, &localX);
6277: return(0);
6278: }
6280: 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)
6281: {
6287: if (!dm->ops->projectfunctionlabellocal) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMProjectFunctionLabelLocal",((PetscObject)dm)->type_name);
6288: (dm->ops->projectfunctionlabellocal) (dm, time, label, numIds, ids, Nc, comps, funcs, ctxs, mode, localX);
6289: return(0);
6290: }
6292: PetscErrorCode DMProjectFieldLocal(DM dm, PetscReal time, Vec localU,
6293: void (**funcs)(PetscInt, PetscInt, PetscInt,
6294: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
6295: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
6296: PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]),
6297: InsertMode mode, Vec localX)
6298: {
6305: if (!dm->ops->projectfieldlocal) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMProjectFieldLocal",((PetscObject)dm)->type_name);
6306: (dm->ops->projectfieldlocal) (dm, time, localU, funcs, mode, localX);
6307: return(0);
6308: }
6310: PetscErrorCode DMProjectFieldLabelLocal(DM dm, PetscReal time, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscInt Nc, const PetscInt comps[], Vec localU,
6311: void (**funcs)(PetscInt, PetscInt, PetscInt,
6312: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
6313: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
6314: PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]),
6315: InsertMode mode, Vec localX)
6316: {
6323: if (!dm->ops->projectfieldlabellocal) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMProjectFieldLocal",((PetscObject)dm)->type_name);
6324: (dm->ops->projectfieldlabellocal)(dm, time, label, numIds, ids, Nc, comps, localU, funcs, mode, localX);
6325: return(0);
6326: }
6328: /*@C
6329: DMComputeL2Diff - This function computes the L_2 difference between a function u and an FEM interpolant solution u_h.
6331: Input Parameters:
6332: + dm - The DM
6333: . time - The time
6334: . funcs - The functions to evaluate for each field component
6335: . ctxs - Optional array of contexts to pass to each function, or NULL.
6336: - X - The coefficient vector u_h, a global vector
6338: Output Parameter:
6339: . diff - The diff ||u - u_h||_2
6341: Level: developer
6343: .seealso: DMProjectFunction(), DMComputeL2FieldDiff(), DMComputeL2GradientDiff()
6344: @*/
6345: PetscErrorCode DMComputeL2Diff(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff)
6346: {
6352: if (!dm->ops->computel2diff) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMComputeL2Diff",((PetscObject)dm)->type_name);
6353: (dm->ops->computel2diff)(dm,time,funcs,ctxs,X,diff);
6354: return(0);
6355: }
6357: /*@C
6358: DMComputeL2GradientDiff - This function computes the L_2 difference between the gradient of a function u and an FEM interpolant solution grad u_h.
6360: Input Parameters:
6361: + dm - The DM
6362: , time - The time
6363: . funcs - The gradient functions to evaluate for each field component
6364: . ctxs - Optional array of contexts to pass to each function, or NULL.
6365: . X - The coefficient vector u_h, a global vector
6366: - n - The vector to project along
6368: Output Parameter:
6369: . diff - The diff ||(grad u - grad u_h) . n||_2
6371: Level: developer
6373: .seealso: DMProjectFunction(), DMComputeL2Diff()
6374: @*/
6375: 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)
6376: {
6382: if (!dm->ops->computel2gradientdiff) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMComputeL2GradientDiff",((PetscObject)dm)->type_name);
6383: (dm->ops->computel2gradientdiff)(dm,time,funcs,ctxs,X,n,diff);
6384: return(0);
6385: }
6387: /*@C
6388: DMComputeL2FieldDiff - This function computes the L_2 difference between a function u and an FEM interpolant solution u_h, separated into field components.
6390: Input Parameters:
6391: + dm - The DM
6392: . time - The time
6393: . funcs - The functions to evaluate for each field component
6394: . ctxs - Optional array of contexts to pass to each function, or NULL.
6395: - X - The coefficient vector u_h, a global vector
6397: Output Parameter:
6398: . diff - The array of differences, ||u^f - u^f_h||_2
6400: Level: developer
6402: .seealso: DMProjectFunction(), DMComputeL2FieldDiff(), DMComputeL2GradientDiff()
6403: @*/
6404: PetscErrorCode DMComputeL2FieldDiff(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal diff[])
6405: {
6411: if (!dm->ops->computel2fielddiff) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMComputeL2FieldDiff",((PetscObject)dm)->type_name);
6412: (dm->ops->computel2fielddiff)(dm,time,funcs,ctxs,X,diff);
6413: return(0);
6414: }
6416: /*@C
6417: DMAdaptLabel - Adapt a dm based on a label with values interpreted as coarsening and refining flags. Specific implementations of DM maybe have
6418: specialized flags, but all implementations should accept flag values DM_ADAPT_DETERMINE, DM_ADAPT_KEEP, DM_ADAPT_REFINE, and DM_ADAPT_COARSEN.
6420: Collective on dm
6422: Input parameters:
6423: + dm - the pre-adaptation DM object
6424: - label - label with the flags
6426: Output parameters:
6427: . dmAdapt - the adapted DM object: may be NULL if an adapted DM could not be produced.
6429: Level: intermediate
6431: .seealso: DMAdaptMetric(), DMCoarsen(), DMRefine()
6432: @*/
6433: PetscErrorCode DMAdaptLabel(DM dm, DMLabel label, DM *dmAdapt)
6434: {
6441: *dmAdapt = NULL;
6442: if (!dm->ops->adaptlabel) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMAdaptLabel",((PetscObject)dm)->type_name);
6443: (dm->ops->adaptlabel)(dm, label, dmAdapt);
6444: return(0);
6445: }
6447: /*@C
6448: DMAdaptMetric - Generates a mesh adapted to the specified metric field using the pragmatic library.
6450: Input Parameters:
6451: + dm - The DM object
6452: . metric - The metric to which the mesh is adapted, defined vertex-wise.
6453: - 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_".
6455: Output Parameter:
6456: . dmAdapt - Pointer to the DM object containing the adapted mesh
6458: Note: The label in the adapted mesh will be registered under the name of the input DMLabel object
6460: Level: advanced
6462: .seealso: DMAdaptLabel(), DMCoarsen(), DMRefine()
6463: @*/
6464: PetscErrorCode DMAdaptMetric(DM dm, Vec metric, DMLabel bdLabel, DM *dmAdapt)
6465: {
6473: *dmAdapt = NULL;
6474: if (!dm->ops->adaptmetric) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMAdaptMetric",((PetscObject)dm)->type_name);
6475: (dm->ops->adaptmetric)(dm, metric, bdLabel, dmAdapt);
6476: return(0);
6477: }
6479: /*@C
6480: DMGetNeighbors - Gets an array containing the MPI rank of all the processes neighbors
6482: Not Collective
6484: Input Parameter:
6485: . dm - The DM
6487: Output Parameter:
6488: . nranks - the number of neighbours
6489: . ranks - the neighbors ranks
6491: Notes:
6492: Do not free the array, it is freed when the DM is destroyed.
6494: Level: beginner
6496: .seealso: DMDAGetNeighbors(), PetscSFGetRanks()
6497: @*/
6498: PetscErrorCode DMGetNeighbors(DM dm,PetscInt *nranks,const PetscMPIInt *ranks[])
6499: {
6504: if (!dm->ops->getneighbors) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMGetNeighbors",((PetscObject)dm)->type_name);
6505: (dm->ops->getneighbors)(dm,nranks,ranks);
6506: return(0);
6507: }
6509: #include <petsc/private/matimpl.h> /* Needed because of coloring->ctype below */
6511: /*
6512: Converts the input vector to a ghosted vector and then calls the standard coloring code.
6513: This has be a different function because it requires DM which is not defined in the Mat library
6514: */
6515: PetscErrorCode MatFDColoringApply_AIJDM(Mat J,MatFDColoring coloring,Vec x1,void *sctx)
6516: {
6520: if (coloring->ctype == IS_COLORING_LOCAL) {
6521: Vec x1local;
6522: DM dm;
6523: MatGetDM(J,&dm);
6524: if (!dm) SETERRQ(PetscObjectComm((PetscObject)J),PETSC_ERR_ARG_INCOMP,"IS_COLORING_LOCAL requires a DM");
6525: DMGetLocalVector(dm,&x1local);
6526: DMGlobalToLocalBegin(dm,x1,INSERT_VALUES,x1local);
6527: DMGlobalToLocalEnd(dm,x1,INSERT_VALUES,x1local);
6528: x1 = x1local;
6529: }
6530: MatFDColoringApply_AIJ(J,coloring,x1,sctx);
6531: if (coloring->ctype == IS_COLORING_LOCAL) {
6532: DM dm;
6533: MatGetDM(J,&dm);
6534: DMRestoreLocalVector(dm,&x1);
6535: }
6536: return(0);
6537: }
6539: /*@
6540: MatFDColoringUseDM - allows a MatFDColoring object to use the DM associated with the matrix to use a IS_COLORING_LOCAL coloring
6542: Input Parameter:
6543: . coloring - the MatFDColoring object
6545: Developer Notes: this routine exists because the PETSc Mat library does not know about the DM objects
6547: Level: advanced
6549: .seealso: MatFDColoring, MatFDColoringCreate(), ISColoringType
6550: @*/
6551: PetscErrorCode MatFDColoringUseDM(Mat coloring,MatFDColoring fdcoloring)
6552: {
6554: coloring->ops->fdcoloringapply = MatFDColoringApply_AIJDM;
6555: return(0);
6556: }