Actual source code: dm.c
petsc-3.7.5 2017-01-01
1: #include <petsc/private/dmimpl.h> /*I "petscdm.h" I*/
2: #include <petsc/private/dmlabelimpl.h> /*I "petscdmlabel.h" I*/
3: #include <petscsf.h>
4: #include <petscds.h>
6: PetscClassId DM_CLASSID;
7: PetscLogEvent DM_Convert, DM_GlobalToLocal, DM_LocalToGlobal, DM_LocalToLocal, DM_LocatePoints, DM_Coarsen, DM_CreateInterpolation, DM_CreateRestriction;
9: const char *const DMBoundaryTypes[] = {"NONE","GHOSTED","MIRROR","PERIODIC","TWIST","DM_BOUNDARY_",0};
13: /*@
14: DMCreate - Creates an empty DM object. The type can then be set with DMSetType().
16: If you never call DMSetType() it will generate an
17: error when you try to use the vector.
19: Collective on MPI_Comm
21: Input Parameter:
22: . comm - The communicator for the DM object
24: Output Parameter:
25: . dm - The DM object
27: Level: beginner
29: .seealso: DMSetType(), DMDA, DMSLICED, DMCOMPOSITE, DMPLEX, DMMOAB, DMNETWORK
30: @*/
31: PetscErrorCode DMCreate(MPI_Comm comm,DM *dm)
32: {
33: DM v;
38: *dm = NULL;
39: PetscSysInitializePackage();
40: VecInitializePackage();
41: MatInitializePackage();
42: DMInitializePackage();
44: PetscHeaderCreate(v, DM_CLASSID, "DM", "Distribution Manager", "DM", comm, DMDestroy, DMView);
46: v->ltogmap = NULL;
47: v->bs = 1;
48: v->coloringtype = IS_COLORING_GLOBAL;
49: PetscSFCreate(comm, &v->sf);
50: PetscSFCreate(comm, &v->defaultSF);
51: v->labels = NULL;
52: v->depthLabel = NULL;
53: v->defaultSection = NULL;
54: v->defaultGlobalSection = NULL;
55: v->defaultConstraintSection = NULL;
56: v->defaultConstraintMat = NULL;
57: v->L = NULL;
58: v->maxCell = NULL;
59: v->bdtype = NULL;
60: v->dimEmbed = PETSC_DEFAULT;
61: {
62: PetscInt i;
63: for (i = 0; i < 10; ++i) {
64: v->nullspaceConstructors[i] = NULL;
65: }
66: }
67: PetscDSCreate(comm, &v->prob);
68: v->dmBC = NULL;
69: v->coarseMesh = NULL;
70: v->outputSequenceNum = -1;
71: v->outputSequenceVal = 0.0;
72: DMSetVecType(v,VECSTANDARD);
73: DMSetMatType(v,MATAIJ);
74: PetscNew(&(v->labels));
75: v->labels->refct = 1;
76: PetscNew(&(v->boundary));
77: v->boundary->refct = 1;
78: *dm = v;
79: return(0);
80: }
84: /*@
85: DMClone - Creates a DM object with the same topology as the original.
87: Collective on MPI_Comm
89: Input Parameter:
90: . dm - The original DM object
92: Output Parameter:
93: . newdm - The new DM object
95: Level: beginner
97: .keywords: DM, topology, create
98: @*/
99: PetscErrorCode DMClone(DM dm, DM *newdm)
100: {
101: PetscSF sf;
102: Vec coords;
103: void *ctx;
104: PetscInt dim;
110: DMCreate(PetscObjectComm((PetscObject) dm), newdm);
111: PetscFree((*newdm)->labels);
112: dm->labels->refct++;
113: (*newdm)->labels = dm->labels;
114: (*newdm)->depthLabel = dm->depthLabel;
115: DMGetDimension(dm, &dim);
116: DMSetDimension(*newdm, dim);
117: if (dm->ops->clone) {
118: (*dm->ops->clone)(dm, newdm);
119: }
120: (*newdm)->setupcalled = dm->setupcalled;
121: DMGetPointSF(dm, &sf);
122: DMSetPointSF(*newdm, sf);
123: DMGetApplicationContext(dm, &ctx);
124: DMSetApplicationContext(*newdm, ctx);
125: if (dm->coordinateDM) {
126: DM ncdm;
127: PetscSection cs;
128: PetscInt pEnd = -1, pEndMax = -1;
130: DMGetDefaultSection(dm->coordinateDM, &cs);
131: if (cs) {PetscSectionGetChart(cs, NULL, &pEnd);}
132: MPI_Allreduce(&pEnd,&pEndMax,1,MPIU_INT,MPI_MAX,PetscObjectComm((PetscObject)dm));
133: if (pEndMax >= 0) {
134: DMClone(dm->coordinateDM, &ncdm);
135: DMSetDefaultSection(ncdm, cs);
136: DMSetCoordinateDM(*newdm, ncdm);
137: DMDestroy(&ncdm);
138: }
139: }
140: DMGetCoordinatesLocal(dm, &coords);
141: if (coords) {
142: DMSetCoordinatesLocal(*newdm, coords);
143: } else {
144: DMGetCoordinates(dm, &coords);
145: if (coords) {DMSetCoordinates(*newdm, coords);}
146: }
147: if (dm->maxCell) {
148: const PetscReal *maxCell, *L;
149: const DMBoundaryType *bd;
150: DMGetPeriodicity(dm, &maxCell, &L, &bd);
151: DMSetPeriodicity(*newdm, maxCell, L, bd);
152: }
153: return(0);
154: }
158: /*@C
159: DMSetVecType - Sets the type of vector created with DMCreateLocalVector() and DMCreateGlobalVector()
161: Logically Collective on DM
163: Input Parameter:
164: + da - initial distributed array
165: . ctype - the vector type, currently either VECSTANDARD or VECCUSP
167: Options Database:
168: . -dm_vec_type ctype
170: Level: intermediate
172: .seealso: DMCreate(), DMDestroy(), DM, DMDAInterpolationType, VecType, DMGetVecType()
173: @*/
174: PetscErrorCode DMSetVecType(DM da,VecType ctype)
175: {
180: PetscFree(da->vectype);
181: PetscStrallocpy(ctype,(char**)&da->vectype);
182: return(0);
183: }
187: /*@C
188: DMGetVecType - Gets the type of vector created with DMCreateLocalVector() and DMCreateGlobalVector()
190: Logically Collective on DM
192: Input Parameter:
193: . da - initial distributed array
195: Output Parameter:
196: . ctype - the vector type
198: Level: intermediate
200: .seealso: DMCreate(), DMDestroy(), DM, DMDAInterpolationType, VecType
201: @*/
202: PetscErrorCode DMGetVecType(DM da,VecType *ctype)
203: {
206: *ctype = da->vectype;
207: return(0);
208: }
212: /*@
213: VecGetDM - Gets the DM defining the data layout of the vector
215: Not collective
217: Input Parameter:
218: . v - The Vec
220: Output Parameter:
221: . dm - The DM
223: Level: intermediate
225: .seealso: VecSetDM(), DMGetLocalVector(), DMGetGlobalVector(), DMSetVecType()
226: @*/
227: PetscErrorCode VecGetDM(Vec v, DM *dm)
228: {
234: PetscObjectQuery((PetscObject) v, "__PETSc_dm", (PetscObject*) dm);
235: return(0);
236: }
240: /*@
241: VecSetDM - Sets the DM defining the data layout of the vector.
243: Not collective
245: Input Parameters:
246: + v - The Vec
247: - dm - The DM
249: 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.
251: Level: intermediate
253: .seealso: VecGetDM(), DMGetLocalVector(), DMGetGlobalVector(), DMSetVecType()
254: @*/
255: PetscErrorCode VecSetDM(Vec v, DM dm)
256: {
262: PetscObjectCompose((PetscObject) v, "__PETSc_dm", (PetscObject) dm);
263: return(0);
264: }
268: /*@C
269: DMSetMatType - Sets the type of matrix created with DMCreateMatrix()
271: Logically Collective on DM
273: Input Parameter:
274: + dm - the DM context
275: . ctype - the matrix type
277: Options Database:
278: . -dm_mat_type ctype
280: Level: intermediate
282: .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMCreateMatrix(), DMSetMatrixPreallocateOnly(), MatType, DMGetMatType()
283: @*/
284: PetscErrorCode DMSetMatType(DM dm,MatType ctype)
285: {
290: PetscFree(dm->mattype);
291: PetscStrallocpy(ctype,(char**)&dm->mattype);
292: return(0);
293: }
297: /*@C
298: DMGetMatType - Gets the type of matrix created with DMCreateMatrix()
300: Logically Collective on DM
302: Input Parameter:
303: . dm - the DM context
305: Output Parameter:
306: . ctype - the matrix type
308: Options Database:
309: . -dm_mat_type ctype
311: Level: intermediate
313: .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMCreateMatrix(), DMSetMatrixPreallocateOnly(), MatType, DMSetMatType()
314: @*/
315: PetscErrorCode DMGetMatType(DM dm,MatType *ctype)
316: {
319: *ctype = dm->mattype;
320: return(0);
321: }
325: /*@
326: MatGetDM - Gets the DM defining the data layout of the matrix
328: Not collective
330: Input Parameter:
331: . A - The Mat
333: Output Parameter:
334: . dm - The DM
336: Level: intermediate
338: .seealso: MatSetDM(), DMCreateMatrix(), DMSetMatType()
339: @*/
340: PetscErrorCode MatGetDM(Mat A, DM *dm)
341: {
347: PetscObjectQuery((PetscObject) A, "__PETSc_dm", (PetscObject*) dm);
348: return(0);
349: }
353: /*@
354: MatSetDM - Sets the DM defining the data layout of the matrix
356: Not collective
358: Input Parameters:
359: + A - The Mat
360: - dm - The DM
362: Level: intermediate
364: .seealso: MatGetDM(), DMCreateMatrix(), DMSetMatType()
365: @*/
366: PetscErrorCode MatSetDM(Mat A, DM dm)
367: {
373: PetscObjectCompose((PetscObject) A, "__PETSc_dm", (PetscObject) dm);
374: return(0);
375: }
379: /*@C
380: DMSetOptionsPrefix - Sets the prefix used for searching for all
381: DM options in the database.
383: Logically Collective on DM
385: Input Parameter:
386: + da - the DM context
387: - prefix - the prefix to prepend to all option names
389: Notes:
390: A hyphen (-) must NOT be given at the beginning of the prefix name.
391: The first character of all runtime options is AUTOMATICALLY the hyphen.
393: Level: advanced
395: .keywords: DM, set, options, prefix, database
397: .seealso: DMSetFromOptions()
398: @*/
399: PetscErrorCode DMSetOptionsPrefix(DM dm,const char prefix[])
400: {
405: PetscObjectSetOptionsPrefix((PetscObject)dm,prefix);
406: if (dm->sf) {
407: PetscObjectSetOptionsPrefix((PetscObject)dm->sf,prefix);
408: }
409: if (dm->defaultSF) {
410: PetscObjectSetOptionsPrefix((PetscObject)dm->defaultSF,prefix);
411: }
412: return(0);
413: }
417: /*@C
418: DMAppendOptionsPrefix - Appends to the prefix used for searching for all
419: DM options in the database.
421: Logically Collective on DM
423: Input Parameters:
424: + dm - the DM context
425: - prefix - the prefix string to prepend to all DM option requests
427: Notes:
428: A hyphen (-) must NOT be given at the beginning of the prefix name.
429: The first character of all runtime options is AUTOMATICALLY the hyphen.
431: Level: advanced
433: .keywords: DM, append, options, prefix, database
435: .seealso: DMSetOptionsPrefix(), DMGetOptionsPrefix()
436: @*/
437: PetscErrorCode DMAppendOptionsPrefix(DM dm,const char prefix[])
438: {
443: PetscObjectAppendOptionsPrefix((PetscObject)dm,prefix);
444: return(0);
445: }
449: /*@C
450: DMGetOptionsPrefix - Gets the prefix used for searching for all
451: DM options in the database.
453: Not Collective
455: Input Parameters:
456: . dm - the DM context
458: Output Parameters:
459: . prefix - pointer to the prefix string used is returned
461: Notes: On the fortran side, the user should pass in a string 'prefix' of
462: sufficient length to hold the prefix.
464: Level: advanced
466: .keywords: DM, set, options, prefix, database
468: .seealso: DMSetOptionsPrefix(), DMAppendOptionsPrefix()
469: @*/
470: PetscErrorCode DMGetOptionsPrefix(DM dm,const char *prefix[])
471: {
476: PetscObjectGetOptionsPrefix((PetscObject)dm,prefix);
477: return(0);
478: }
482: static PetscErrorCode DMCountNonCyclicReferences(DM dm, PetscBool recurseCoarse, PetscBool recurseFine, PetscInt *ncrefct)
483: {
484: PetscInt i, refct = ((PetscObject) dm)->refct;
485: DMNamedVecLink nlink;
489: /* count all the circular references of DM and its contained Vecs */
490: for (i=0; i<DM_MAX_WORK_VECTORS; i++) {
491: if (dm->localin[i]) refct--;
492: if (dm->globalin[i]) refct--;
493: }
494: for (nlink=dm->namedglobal; nlink; nlink=nlink->next) refct--;
495: for (nlink=dm->namedlocal; nlink; nlink=nlink->next) refct--;
496: if (dm->x) {
497: DM obj;
498: VecGetDM(dm->x, &obj);
499: if (obj == dm) refct--;
500: }
501: if (dm->coarseMesh && dm->coarseMesh->fineMesh == dm) {
502: refct--;
503: if (recurseCoarse) {
504: PetscInt coarseCount;
506: DMCountNonCyclicReferences(dm->coarseMesh, PETSC_TRUE, PETSC_FALSE,&coarseCount);
507: refct += coarseCount;
508: }
509: }
510: if (dm->fineMesh && dm->fineMesh->coarseMesh == dm) {
511: refct--;
512: if (recurseFine) {
513: PetscInt fineCount;
515: DMCountNonCyclicReferences(dm->fineMesh, PETSC_FALSE, PETSC_TRUE,&fineCount);
516: refct += fineCount;
517: }
518: }
519: *ncrefct = refct;
520: return(0);
521: }
523: PetscErrorCode DMBoundaryDestroy(DMBoundaryLinkList *boundary);
527: PetscErrorCode DMDestroyLabelLinkList(DM dm)
528: {
532: if (!--(dm->labels->refct)) {
533: DMLabelLink next = dm->labels->next;
535: /* destroy the labels */
536: while (next) {
537: DMLabelLink tmp = next->next;
539: DMLabelDestroy(&next->label);
540: PetscFree(next);
541: next = tmp;
542: }
543: PetscFree(dm->labels);
544: }
545: return(0);
546: }
550: /*@
551: DMDestroy - Destroys a vector packer or DM.
553: Collective on DM
555: Input Parameter:
556: . dm - the DM object to destroy
558: Level: developer
560: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix()
562: @*/
563: PetscErrorCode DMDestroy(DM *dm)
564: {
565: PetscInt i, cnt;
566: DMNamedVecLink nlink,nnext;
570: if (!*dm) return(0);
573: /* count all non-cyclic references in the doubly-linked list of coarse<->fine meshes */
574: DMCountNonCyclicReferences(*dm,PETSC_TRUE,PETSC_TRUE,&cnt);
575: --((PetscObject)(*dm))->refct;
576: if (--cnt > 0) {*dm = 0; return(0);}
577: /*
578: Need this test because the dm references the vectors that
579: reference the dm, so destroying the dm calls destroy on the
580: vectors that cause another destroy on the dm
581: */
582: if (((PetscObject)(*dm))->refct < 0) return(0);
583: ((PetscObject) (*dm))->refct = 0;
584: for (i=0; i<DM_MAX_WORK_VECTORS; i++) {
585: if ((*dm)->localout[i]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Destroying a DM that has a local vector obtained with DMGetLocalVector()");
586: VecDestroy(&(*dm)->localin[i]);
587: }
588: nnext=(*dm)->namedglobal;
589: (*dm)->namedglobal = NULL;
590: for (nlink=nnext; nlink; nlink=nnext) { /* Destroy the named vectors */
591: nnext = nlink->next;
592: if (nlink->status != DMVEC_STATUS_IN) SETERRQ1(((PetscObject)*dm)->comm,PETSC_ERR_ARG_WRONGSTATE,"DM still has Vec named '%s' checked out",nlink->name);
593: PetscFree(nlink->name);
594: VecDestroy(&nlink->X);
595: PetscFree(nlink);
596: }
597: nnext=(*dm)->namedlocal;
598: (*dm)->namedlocal = NULL;
599: for (nlink=nnext; nlink; nlink=nnext) { /* Destroy the named local vectors */
600: nnext = nlink->next;
601: if (nlink->status != DMVEC_STATUS_IN) SETERRQ1(((PetscObject)*dm)->comm,PETSC_ERR_ARG_WRONGSTATE,"DM still has Vec named '%s' checked out",nlink->name);
602: PetscFree(nlink->name);
603: VecDestroy(&nlink->X);
604: PetscFree(nlink);
605: }
607: /* Destroy the list of hooks */
608: {
609: DMCoarsenHookLink link,next;
610: for (link=(*dm)->coarsenhook; link; link=next) {
611: next = link->next;
612: PetscFree(link);
613: }
614: (*dm)->coarsenhook = NULL;
615: }
616: {
617: DMRefineHookLink link,next;
618: for (link=(*dm)->refinehook; link; link=next) {
619: next = link->next;
620: PetscFree(link);
621: }
622: (*dm)->refinehook = NULL;
623: }
624: {
625: DMSubDomainHookLink link,next;
626: for (link=(*dm)->subdomainhook; link; link=next) {
627: next = link->next;
628: PetscFree(link);
629: }
630: (*dm)->subdomainhook = NULL;
631: }
632: {
633: DMGlobalToLocalHookLink link,next;
634: for (link=(*dm)->gtolhook; link; link=next) {
635: next = link->next;
636: PetscFree(link);
637: }
638: (*dm)->gtolhook = NULL;
639: }
640: {
641: DMLocalToGlobalHookLink link,next;
642: for (link=(*dm)->ltoghook; link; link=next) {
643: next = link->next;
644: PetscFree(link);
645: }
646: (*dm)->ltoghook = NULL;
647: }
648: /* Destroy the work arrays */
649: {
650: DMWorkLink link,next;
651: if ((*dm)->workout) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Work array still checked out");
652: for (link=(*dm)->workin; link; link=next) {
653: next = link->next;
654: PetscFree(link->mem);
655: PetscFree(link);
656: }
657: (*dm)->workin = NULL;
658: }
659: if (!--((*dm)->labels->refct)) {
660: DMLabelLink next = (*dm)->labels->next;
662: /* destroy the labels */
663: while (next) {
664: DMLabelLink tmp = next->next;
666: DMLabelDestroy(&next->label);
667: PetscFree(next);
668: next = tmp;
669: }
670: PetscFree((*dm)->labels);
671: }
672: DMBoundaryDestroy(&(*dm)->boundary);
674: PetscObjectDestroy(&(*dm)->dmksp);
675: PetscObjectDestroy(&(*dm)->dmsnes);
676: PetscObjectDestroy(&(*dm)->dmts);
678: if ((*dm)->ctx && (*dm)->ctxdestroy) {
679: (*(*dm)->ctxdestroy)(&(*dm)->ctx);
680: }
681: VecDestroy(&(*dm)->x);
682: MatFDColoringDestroy(&(*dm)->fd);
683: DMClearGlobalVectors(*dm);
684: ISLocalToGlobalMappingDestroy(&(*dm)->ltogmap);
685: PetscFree((*dm)->vectype);
686: PetscFree((*dm)->mattype);
688: PetscSectionDestroy(&(*dm)->defaultSection);
689: PetscSectionDestroy(&(*dm)->defaultGlobalSection);
690: PetscLayoutDestroy(&(*dm)->map);
691: PetscSectionDestroy(&(*dm)->defaultConstraintSection);
692: MatDestroy(&(*dm)->defaultConstraintMat);
693: PetscSFDestroy(&(*dm)->sf);
694: PetscSFDestroy(&(*dm)->defaultSF);
695: PetscSFDestroy(&(*dm)->sfNatural);
697: if ((*dm)->coarseMesh && (*dm)->coarseMesh->fineMesh == *dm) {
698: DMSetFineDM((*dm)->coarseMesh,NULL);
699: }
700: DMDestroy(&(*dm)->coarseMesh);
701: if ((*dm)->fineMesh && (*dm)->fineMesh->coarseMesh == *dm) {
702: DMSetCoarseDM((*dm)->fineMesh,NULL);
703: }
704: DMDestroy(&(*dm)->fineMesh);
705: DMDestroy(&(*dm)->coordinateDM);
706: VecDestroy(&(*dm)->coordinates);
707: VecDestroy(&(*dm)->coordinatesLocal);
708: PetscFree3((*dm)->L,(*dm)->maxCell,(*dm)->bdtype);
710: PetscDSDestroy(&(*dm)->prob);
711: DMDestroy(&(*dm)->dmBC);
712: /* if memory was published with SAWs then destroy it */
713: PetscObjectSAWsViewOff((PetscObject)*dm);
715: (*(*dm)->ops->destroy)(*dm);
716: /* We do not destroy (*dm)->data here so that we can reference count backend objects */
717: PetscHeaderDestroy(dm);
718: return(0);
719: }
723: /*@
724: DMSetUp - sets up the data structures inside a DM object
726: Collective on DM
728: Input Parameter:
729: . dm - the DM object to setup
731: Level: developer
733: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix()
735: @*/
736: PetscErrorCode DMSetUp(DM dm)
737: {
742: if (dm->setupcalled) return(0);
743: if (dm->ops->setup) {
744: (*dm->ops->setup)(dm);
745: }
746: dm->setupcalled = PETSC_TRUE;
747: return(0);
748: }
752: /*@
753: DMSetFromOptions - sets parameters in a DM from the options database
755: Collective on DM
757: Input Parameter:
758: . dm - the DM object to set options for
760: Options Database:
761: + -dm_preallocate_only - Only preallocate the matrix for DMCreateMatrix(), but do not fill it with zeros
762: . -dm_vec_type <type> - type of vector to create inside DM
763: . -dm_mat_type <type> - type of matrix to create inside DM
764: - -dm_coloring_type - <global or ghosted>
766: Level: developer
768: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix()
770: @*/
771: PetscErrorCode DMSetFromOptions(DM dm)
772: {
773: char typeName[256];
774: PetscBool flg;
779: if (dm->sf) {
780: PetscSFSetFromOptions(dm->sf);
781: }
782: if (dm->defaultSF) {
783: PetscSFSetFromOptions(dm->defaultSF);
784: }
785: PetscObjectOptionsBegin((PetscObject)dm);
786: PetscOptionsBool("-dm_preallocate_only","only preallocate matrix, but do not set column indices","DMSetMatrixPreallocateOnly",dm->prealloc_only,&dm->prealloc_only,NULL);
787: PetscOptionsFList("-dm_vec_type","Vector type used for created vectors","DMSetVecType",VecList,dm->vectype,typeName,256,&flg);
788: if (flg) {
789: DMSetVecType(dm,typeName);
790: }
791: PetscOptionsFList("-dm_mat_type","Matrix type used for created matrices","DMSetMatType",MatList,dm->mattype ? dm->mattype : typeName,typeName,sizeof(typeName),&flg);
792: if (flg) {
793: DMSetMatType(dm,typeName);
794: }
795: PetscOptionsEnum("-dm_is_coloring_type","Global or local coloring of Jacobian","ISColoringType",ISColoringTypes,(PetscEnum)dm->coloringtype,(PetscEnum*)&dm->coloringtype,NULL);
796: if (dm->ops->setfromoptions) {
797: (*dm->ops->setfromoptions)(PetscOptionsObject,dm);
798: }
799: /* process any options handlers added with PetscObjectAddOptionsHandler() */
800: PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject) dm);
801: PetscOptionsEnd();
802: return(0);
803: }
807: /*@C
808: DMView - Views a DM
810: Collective on DM
812: Input Parameter:
813: + dm - the DM object to view
814: - v - the viewer
816: Level: beginner
818: .seealso DMDestroy(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix()
820: @*/
821: PetscErrorCode DMView(DM dm,PetscViewer v)
822: {
824: PetscBool isbinary;
828: if (!v) {
829: PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)dm),&v);
830: }
831: PetscObjectPrintClassNamePrefixType((PetscObject)dm,v);
832: PetscObjectTypeCompare((PetscObject)v,PETSCVIEWERBINARY,&isbinary);
833: if (isbinary) {
834: PetscInt classid = DM_FILE_CLASSID;
835: char type[256];
837: PetscViewerBinaryWrite(v,&classid,1,PETSC_INT,PETSC_FALSE);
838: PetscStrncpy(type,((PetscObject)dm)->type_name,256);
839: PetscViewerBinaryWrite(v,type,256,PETSC_CHAR,PETSC_FALSE);
840: }
841: if (dm->ops->view) {
842: (*dm->ops->view)(dm,v);
843: }
844: return(0);
845: }
849: /*@
850: DMCreateGlobalVector - Creates a global vector from a DM object
852: Collective on DM
854: Input Parameter:
855: . dm - the DM object
857: Output Parameter:
858: . vec - the global vector
860: Level: beginner
862: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix()
864: @*/
865: PetscErrorCode DMCreateGlobalVector(DM dm,Vec *vec)
866: {
871: (*dm->ops->createglobalvector)(dm,vec);
872: return(0);
873: }
877: /*@
878: DMCreateLocalVector - Creates a local vector from a DM object
880: Not Collective
882: Input Parameter:
883: . dm - the DM object
885: Output Parameter:
886: . vec - the local vector
888: Level: beginner
890: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix()
892: @*/
893: PetscErrorCode DMCreateLocalVector(DM dm,Vec *vec)
894: {
899: (*dm->ops->createlocalvector)(dm,vec);
900: return(0);
901: }
905: /*@
906: DMGetLocalToGlobalMapping - Accesses the local-to-global mapping in a DM.
908: Collective on DM
910: Input Parameter:
911: . dm - the DM that provides the mapping
913: Output Parameter:
914: . ltog - the mapping
916: Level: intermediate
918: Notes:
919: This mapping can then be used by VecSetLocalToGlobalMapping() or
920: MatSetLocalToGlobalMapping().
922: .seealso: DMCreateLocalVector()
923: @*/
924: PetscErrorCode DMGetLocalToGlobalMapping(DM dm,ISLocalToGlobalMapping *ltog)
925: {
926: PetscInt bs = -1;
932: if (!dm->ltogmap) {
933: PetscSection section, sectionGlobal;
935: DMGetDefaultSection(dm, §ion);
936: if (section) {
937: PetscInt *ltog;
938: PetscInt pStart, pEnd, size, p, l;
940: DMGetDefaultGlobalSection(dm, §ionGlobal);
941: PetscSectionGetChart(section, &pStart, &pEnd);
942: PetscSectionGetStorageSize(section, &size);
943: PetscMalloc1(size, <og); /* We want the local+overlap size */
944: for (p = pStart, l = 0; p < pEnd; ++p) {
945: PetscInt bdof, cdof, dof, off, c;
947: /* Should probably use constrained dofs */
948: PetscSectionGetDof(section, p, &dof);
949: PetscSectionGetOffset(sectionGlobal, p, &off);
950: PetscSectionGetConstraintDof(sectionGlobal, p, &cdof);
951: bdof = cdof ? 1 : dof;
952: if (bs < 0) {bs = bdof;}
953: else if (bs != bdof) {bs = 1;}
954: for (c = 0; c < dof; ++c, ++l) {
955: ltog[l] = off+c;
956: }
957: }
958: ISLocalToGlobalMappingCreate(PETSC_COMM_SELF, bs < 0 ? 1 : bs, size, ltog, PETSC_OWN_POINTER, &dm->ltogmap);
959: PetscLogObjectParent((PetscObject)dm, (PetscObject)dm->ltogmap);
960: } else {
961: if (!dm->ops->getlocaltoglobalmapping) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM can not create LocalToGlobalMapping");
962: (*dm->ops->getlocaltoglobalmapping)(dm);
963: }
964: }
965: *ltog = dm->ltogmap;
966: return(0);
967: }
971: /*@
972: DMGetBlockSize - Gets the inherent block size associated with a DM
974: Not Collective
976: Input Parameter:
977: . dm - the DM with block structure
979: Output Parameter:
980: . bs - the block size, 1 implies no exploitable block structure
982: Level: intermediate
984: .seealso: ISCreateBlock(), VecSetBlockSize(), MatSetBlockSize(), DMGetLocalToGlobalMapping()
985: @*/
986: PetscErrorCode DMGetBlockSize(DM dm,PetscInt *bs)
987: {
991: if (dm->bs < 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"DM does not have enough information to provide a block size yet");
992: *bs = dm->bs;
993: return(0);
994: }
998: /*@
999: DMCreateInterpolation - Gets interpolation matrix between two DM objects
1001: Collective on DM
1003: Input Parameter:
1004: + dm1 - the DM object
1005: - dm2 - the second, finer DM object
1007: Output Parameter:
1008: + mat - the interpolation
1009: - vec - the scaling (optional)
1011: Level: developer
1013: 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
1014: DMCoarsen(). The coordinates set into the DMDA are completely ignored in computing the interpolation.
1016: For DMDA objects you can use this interpolation (more precisely the interpolation from the DMGetCoordinateDM()) to interpolate the mesh coordinate vectors
1017: EXCEPT in the periodic case where it does not make sense since the coordinate vectors are not periodic.
1020: .seealso DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateColoring(), DMCreateMatrix(), DMRefine(), DMCoarsen(), DMCreateRestriction()
1022: @*/
1023: PetscErrorCode DMCreateInterpolation(DM dm1,DM dm2,Mat *mat,Vec *vec)
1024: {
1030: PetscLogEventBegin(DM_CreateInterpolation,dm1,dm2,0,0);
1031: (*dm1->ops->createinterpolation)(dm1,dm2,mat,vec);
1032: PetscLogEventEnd(DM_CreateInterpolation,dm1,dm2,0,0);
1033: return(0);
1034: }
1038: /*@
1039: DMCreateRestriction - Gets restriction matrix between two DM objects
1041: Collective on DM
1043: Input Parameter:
1044: + dm1 - the DM object
1045: - dm2 - the second, finer DM object
1047: Output Parameter:
1048: . mat - the restriction
1051: Level: developer
1053: 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
1054: DMCoarsen(). The coordinates set into the DMDA are completely ignored in computing the interpolation.
1056:
1057: .seealso DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateColoring(), DMCreateMatrix(), DMRefine(), DMCoarsen(), DMCreateInterpolation()
1059: @*/
1060: PetscErrorCode DMCreateRestriction(DM dm1,DM dm2,Mat *mat)
1061: {
1067: PetscLogEventBegin(DM_CreateRestriction,dm1,dm2,0,0);
1068: (*dm1->ops->createrestriction)(dm1,dm2,mat);
1069: PetscLogEventEnd(DM_CreateRestriction,dm1,dm2,0,0);
1070: return(0);
1071: }
1075: /*@
1076: DMCreateInjection - Gets injection matrix between two DM objects
1078: Collective on DM
1080: Input Parameter:
1081: + dm1 - the DM object
1082: - dm2 - the second, finer DM object
1084: Output Parameter:
1085: . mat - the injection
1087: Level: developer
1089: 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
1090: DMCoarsen(). The coordinates set into the DMDA are completely ignored in computing the injection.
1092: .seealso DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateColoring(), DMCreateMatrix(), DMCreateInterpolation()
1094: @*/
1095: PetscErrorCode DMCreateInjection(DM dm1,DM dm2,Mat *mat)
1096: {
1102: if (!*dm1->ops->getinjection) SETERRQ(PetscObjectComm((PetscObject)dm1),PETSC_ERR_SUP,"DMCreateInjection not implemented for this type");
1103: (*dm1->ops->getinjection)(dm1,dm2,mat);
1104: return(0);
1105: }
1109: /*@
1110: DMCreateColoring - Gets coloring for a DM
1112: Collective on DM
1114: Input Parameter:
1115: + dm - the DM object
1116: - ctype - IS_COLORING_GHOSTED or IS_COLORING_GLOBAL
1118: Output Parameter:
1119: . coloring - the coloring
1121: Level: developer
1123: .seealso DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateMatrix(), DMSetMatType()
1125: @*/
1126: PetscErrorCode DMCreateColoring(DM dm,ISColoringType ctype,ISColoring *coloring)
1127: {
1132: if (!dm->ops->getcoloring) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"No coloring for this type of DM yet");
1133: (*dm->ops->getcoloring)(dm,ctype,coloring);
1134: return(0);
1135: }
1139: /*@
1140: DMCreateMatrix - Gets empty Jacobian for a DM
1142: Collective on DM
1144: Input Parameter:
1145: . dm - the DM object
1147: Output Parameter:
1148: . mat - the empty Jacobian
1150: Level: beginner
1152: Notes: This properly preallocates the number of nonzeros in the sparse matrix so you
1153: do not need to do it yourself.
1155: By default it also sets the nonzero structure and puts in the zero entries. To prevent setting
1156: the nonzero pattern call DMDASetMatPreallocateOnly()
1158: For structured grid problems, when you call MatView() on this matrix it is displayed using the global natural ordering, NOT in the ordering used
1159: internally by PETSc.
1161: For structured grid problems, in general it is easiest to use MatSetValuesStencil() or MatSetValuesLocal() to put values into the matrix because MatSetValues() requires
1162: the indices for the global numbering for DMDAs which is complicated.
1164: .seealso DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMSetMatType()
1166: @*/
1167: PetscErrorCode DMCreateMatrix(DM dm,Mat *mat)
1168: {
1173: MatInitializePackage();
1176: (*dm->ops->creatematrix)(dm,mat);
1177: return(0);
1178: }
1182: /*@
1183: DMSetMatrixPreallocateOnly - When DMCreateMatrix() is called the matrix will be properly
1184: preallocated but the nonzero structure and zero values will not be set.
1186: Logically Collective on DM
1188: Input Parameter:
1189: + dm - the DM
1190: - only - PETSC_TRUE if only want preallocation
1192: Level: developer
1193: .seealso DMCreateMatrix()
1194: @*/
1195: PetscErrorCode DMSetMatrixPreallocateOnly(DM dm, PetscBool only)
1196: {
1199: dm->prealloc_only = only;
1200: return(0);
1201: }
1205: /*@C
1206: DMGetWorkArray - Gets a work array guaranteed to be at least the input size, restore with DMRestoreWorkArray()
1208: Not Collective
1210: Input Parameters:
1211: + dm - the DM object
1212: . count - The minium size
1213: - dtype - data type (PETSC_REAL, PETSC_SCALAR, PETSC_INT)
1215: Output Parameter:
1216: . array - the work array
1218: Level: developer
1220: .seealso DMDestroy(), DMCreate()
1221: @*/
1222: PetscErrorCode DMGetWorkArray(DM dm,PetscInt count,PetscDataType dtype,void *mem)
1223: {
1225: DMWorkLink link;
1226: size_t dsize;
1231: if (dm->workin) {
1232: link = dm->workin;
1233: dm->workin = dm->workin->next;
1234: } else {
1235: PetscNewLog(dm,&link);
1236: }
1237: PetscDataTypeGetSize(dtype,&dsize);
1238: if (dsize*count > link->bytes) {
1239: PetscFree(link->mem);
1240: PetscMalloc(dsize*count,&link->mem);
1241: link->bytes = dsize*count;
1242: }
1243: link->next = dm->workout;
1244: dm->workout = link;
1245: *(void**)mem = link->mem;
1246: return(0);
1247: }
1251: /*@C
1252: DMRestoreWorkArray - Restores a work array guaranteed to be at least the input size, restore with DMRestoreWorkArray()
1254: Not Collective
1256: Input Parameters:
1257: + dm - the DM object
1258: . count - The minium size
1259: - dtype - data type (PETSC_REAL, PETSC_SCALAR, PETSC_INT)
1261: Output Parameter:
1262: . array - the work array
1264: Level: developer
1266: .seealso DMDestroy(), DMCreate()
1267: @*/
1268: PetscErrorCode DMRestoreWorkArray(DM dm,PetscInt count,PetscDataType dtype,void *mem)
1269: {
1270: DMWorkLink *p,link;
1275: for (p=&dm->workout; (link=*p); p=&link->next) {
1276: if (link->mem == *(void**)mem) {
1277: *p = link->next;
1278: link->next = dm->workin;
1279: dm->workin = link;
1280: *(void**)mem = NULL;
1281: return(0);
1282: }
1283: }
1284: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Array was not checked out");
1285: }
1289: PetscErrorCode DMSetNullSpaceConstructor(DM dm, PetscInt field, PetscErrorCode (*nullsp)(DM dm, PetscInt field, MatNullSpace *nullSpace))
1290: {
1293: if (field >= 10) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle %d >= 10 fields", field);
1294: dm->nullspaceConstructors[field] = nullsp;
1295: return(0);
1296: }
1300: /*@C
1301: DMCreateFieldIS - Creates a set of IS objects with the global indices of dofs for each field
1303: Not collective
1305: Input Parameter:
1306: . dm - the DM object
1308: Output Parameters:
1309: + numFields - The number of fields (or NULL if not requested)
1310: . fieldNames - The name for each field (or NULL if not requested)
1311: - fields - The global indices for each field (or NULL if not requested)
1313: Level: intermediate
1315: Notes:
1316: The user is responsible for freeing all requested arrays. In particular, every entry of names should be freed with
1317: PetscFree(), every entry of fields should be destroyed with ISDestroy(), and both arrays should be freed with
1318: PetscFree().
1320: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix()
1321: @*/
1322: PetscErrorCode DMCreateFieldIS(DM dm, PetscInt *numFields, char ***fieldNames, IS **fields)
1323: {
1324: PetscSection section, sectionGlobal;
1329: if (numFields) {
1331: *numFields = 0;
1332: }
1333: if (fieldNames) {
1335: *fieldNames = NULL;
1336: }
1337: if (fields) {
1339: *fields = NULL;
1340: }
1341: DMGetDefaultSection(dm, §ion);
1342: if (section) {
1343: PetscInt *fieldSizes, **fieldIndices;
1344: PetscInt nF, f, pStart, pEnd, p;
1346: DMGetDefaultGlobalSection(dm, §ionGlobal);
1347: PetscSectionGetNumFields(section, &nF);
1348: PetscMalloc2(nF,&fieldSizes,nF,&fieldIndices);
1349: PetscSectionGetChart(sectionGlobal, &pStart, &pEnd);
1350: for (f = 0; f < nF; ++f) {
1351: fieldSizes[f] = 0;
1352: }
1353: for (p = pStart; p < pEnd; ++p) {
1354: PetscInt gdof;
1356: PetscSectionGetDof(sectionGlobal, p, &gdof);
1357: if (gdof > 0) {
1358: for (f = 0; f < nF; ++f) {
1359: PetscInt fdof, fcdof;
1361: PetscSectionGetFieldDof(section, p, f, &fdof);
1362: PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);
1363: fieldSizes[f] += fdof-fcdof;
1364: }
1365: }
1366: }
1367: for (f = 0; f < nF; ++f) {
1368: PetscMalloc1(fieldSizes[f], &fieldIndices[f]);
1369: fieldSizes[f] = 0;
1370: }
1371: for (p = pStart; p < pEnd; ++p) {
1372: PetscInt gdof, goff;
1374: PetscSectionGetDof(sectionGlobal, p, &gdof);
1375: if (gdof > 0) {
1376: PetscSectionGetOffset(sectionGlobal, p, &goff);
1377: for (f = 0; f < nF; ++f) {
1378: PetscInt fdof, fcdof, fc;
1380: PetscSectionGetFieldDof(section, p, f, &fdof);
1381: PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);
1382: for (fc = 0; fc < fdof-fcdof; ++fc, ++fieldSizes[f]) {
1383: fieldIndices[f][fieldSizes[f]] = goff++;
1384: }
1385: }
1386: }
1387: }
1388: if (numFields) *numFields = nF;
1389: if (fieldNames) {
1390: PetscMalloc1(nF, fieldNames);
1391: for (f = 0; f < nF; ++f) {
1392: const char *fieldName;
1394: PetscSectionGetFieldName(section, f, &fieldName);
1395: PetscStrallocpy(fieldName, (char**) &(*fieldNames)[f]);
1396: }
1397: }
1398: if (fields) {
1399: PetscMalloc1(nF, fields);
1400: for (f = 0; f < nF; ++f) {
1401: ISCreateGeneral(PetscObjectComm((PetscObject)dm), fieldSizes[f], fieldIndices[f], PETSC_OWN_POINTER, &(*fields)[f]);
1402: }
1403: }
1404: PetscFree2(fieldSizes,fieldIndices);
1405: } else if (dm->ops->createfieldis) {
1406: (*dm->ops->createfieldis)(dm, numFields, fieldNames, fields);
1407: }
1408: return(0);
1409: }
1414: /*@C
1415: DMCreateFieldDecomposition - Returns a list of IS objects defining a decomposition of a problem into subproblems
1416: corresponding to different fields: each IS contains the global indices of the dofs of the
1417: corresponding field. The optional list of DMs define the DM for each subproblem.
1418: Generalizes DMCreateFieldIS().
1420: Not collective
1422: Input Parameter:
1423: . dm - the DM object
1425: Output Parameters:
1426: + len - The number of subproblems in the field decomposition (or NULL if not requested)
1427: . namelist - The name for each field (or NULL if not requested)
1428: . islist - The global indices for each field (or NULL if not requested)
1429: - dmlist - The DMs for each field subproblem (or NULL, if not requested; if NULL is returned, no DMs are defined)
1431: Level: intermediate
1433: Notes:
1434: The user is responsible for freeing all requested arrays. In particular, every entry of names should be freed with
1435: PetscFree(), every entry of is should be destroyed with ISDestroy(), every entry of dm should be destroyed with DMDestroy(),
1436: and all of the arrays should be freed with PetscFree().
1438: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateFieldIS()
1439: @*/
1440: PetscErrorCode DMCreateFieldDecomposition(DM dm, PetscInt *len, char ***namelist, IS **islist, DM **dmlist)
1441: {
1446: if (len) {
1448: *len = 0;
1449: }
1450: if (namelist) {
1452: *namelist = 0;
1453: }
1454: if (islist) {
1456: *islist = 0;
1457: }
1458: if (dmlist) {
1460: *dmlist = 0;
1461: }
1462: /*
1463: Is it a good idea to apply the following check across all impls?
1464: Perhaps some impls can have a well-defined decomposition before DMSetUp?
1465: This, however, follows the general principle that accessors are not well-behaved until the object is set up.
1466: */
1467: if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE, "Decomposition defined only after DMSetUp");
1468: if (!dm->ops->createfielddecomposition) {
1469: PetscSection section;
1470: PetscInt numFields, f;
1472: DMGetDefaultSection(dm, §ion);
1473: if (section) {PetscSectionGetNumFields(section, &numFields);}
1474: if (section && numFields && dm->ops->createsubdm) {
1475: if (len) *len = numFields;
1476: if (namelist) {PetscMalloc1(numFields,namelist);}
1477: if (islist) {PetscMalloc1(numFields,islist);}
1478: if (dmlist) {PetscMalloc1(numFields,dmlist);}
1479: for (f = 0; f < numFields; ++f) {
1480: const char *fieldName;
1482: DMCreateSubDM(dm, 1, &f, islist ? &(*islist)[f] : NULL, dmlist ? &(*dmlist)[f] : NULL);
1483: if (namelist) {
1484: PetscSectionGetFieldName(section, f, &fieldName);
1485: PetscStrallocpy(fieldName, (char**) &(*namelist)[f]);
1486: }
1487: }
1488: } else {
1489: DMCreateFieldIS(dm, len, namelist, islist);
1490: /* By default there are no DMs associated with subproblems. */
1491: if (dmlist) *dmlist = NULL;
1492: }
1493: } else {
1494: (*dm->ops->createfielddecomposition)(dm,len,namelist,islist,dmlist);
1495: }
1496: return(0);
1497: }
1501: /*@
1502: DMCreateSubDM - Returns an IS and DM encapsulating a subproblem defined by the fields passed in.
1503: The fields are defined by DMCreateFieldIS().
1505: Not collective
1507: Input Parameters:
1508: + dm - the DM object
1509: . numFields - number of fields in this subproblem
1510: - len - The number of subproblems in the decomposition (or NULL if not requested)
1512: Output Parameters:
1513: . is - The global indices for the subproblem
1514: - dm - The DM for the subproblem
1516: Level: intermediate
1518: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateFieldIS()
1519: @*/
1520: PetscErrorCode DMCreateSubDM(DM dm, PetscInt numFields, PetscInt fields[], IS *is, DM *subdm)
1521: {
1529: if (dm->ops->createsubdm) {
1530: (*dm->ops->createsubdm)(dm, numFields, fields, is, subdm);
1531: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "This type has no DMCreateSubDM implementation defined");
1532: return(0);
1533: }
1538: /*@C
1539: DMCreateDomainDecomposition - Returns lists of IS objects defining a decomposition of a problem into subproblems
1540: corresponding to restrictions to pairs nested subdomains: each IS contains the global
1541: indices of the dofs of the corresponding subdomains. The inner subdomains conceptually
1542: define a nonoverlapping covering, while outer subdomains can overlap.
1543: The optional list of DMs define the DM for each subproblem.
1545: Not collective
1547: Input Parameter:
1548: . dm - the DM object
1550: Output Parameters:
1551: + len - The number of subproblems in the domain decomposition (or NULL if not requested)
1552: . namelist - The name for each subdomain (or NULL if not requested)
1553: . innerislist - The global indices for each inner subdomain (or NULL, if not requested)
1554: . outerislist - The global indices for each outer subdomain (or NULL, if not requested)
1555: - dmlist - The DMs for each subdomain subproblem (or NULL, if not requested; if NULL is returned, no DMs are defined)
1557: Level: intermediate
1559: Notes:
1560: The user is responsible for freeing all requested arrays. In particular, every entry of names should be freed with
1561: PetscFree(), every entry of is should be destroyed with ISDestroy(), every entry of dm should be destroyed with DMDestroy(),
1562: and all of the arrays should be freed with PetscFree().
1564: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateDomainDecompositionDM(), DMCreateFieldDecomposition()
1565: @*/
1566: PetscErrorCode DMCreateDomainDecomposition(DM dm, PetscInt *len, char ***namelist, IS **innerislist, IS **outerislist, DM **dmlist)
1567: {
1568: PetscErrorCode ierr;
1569: DMSubDomainHookLink link;
1570: PetscInt i,l;
1579: /*
1580: Is it a good idea to apply the following check across all impls?
1581: Perhaps some impls can have a well-defined decomposition before DMSetUp?
1582: This, however, follows the general principle that accessors are not well-behaved until the object is set up.
1583: */
1584: if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE, "Decomposition defined only after DMSetUp");
1585: if (dm->ops->createdomaindecomposition) {
1586: (*dm->ops->createdomaindecomposition)(dm,&l,namelist,innerislist,outerislist,dmlist);
1587: /* copy subdomain hooks and context over to the subdomain DMs */
1588: if (dmlist && *dmlist) {
1589: for (i = 0; i < l; i++) {
1590: for (link=dm->subdomainhook; link; link=link->next) {
1591: if (link->ddhook) {(*link->ddhook)(dm,(*dmlist)[i],link->ctx);}
1592: }
1593: (*dmlist)[i]->ctx = dm->ctx;
1594: }
1595: }
1596: if (len) *len = l;
1597: }
1598: return(0);
1599: }
1604: /*@C
1605: DMCreateDomainDecompositionScatters - Returns scatters to the subdomain vectors from the global vector
1607: Not collective
1609: Input Parameters:
1610: + dm - the DM object
1611: . n - the number of subdomain scatters
1612: - subdms - the local subdomains
1614: Output Parameters:
1615: + n - the number of scatters returned
1616: . iscat - scatter from global vector to nonoverlapping global vector entries on subdomain
1617: . oscat - scatter from global vector to overlapping global vector entries on subdomain
1618: - gscat - scatter from global vector to local vector on subdomain (fills in ghosts)
1620: Notes: This is an alternative to the iis and ois arguments in DMCreateDomainDecomposition that allow for the solution
1621: of general nonlinear problems with overlapping subdomain methods. While merely having index sets that enable subsets
1622: of the residual equations to be created is fine for linear problems, nonlinear problems require local assembly of
1623: solution and residual data.
1625: Level: developer
1627: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateFieldIS()
1628: @*/
1629: PetscErrorCode DMCreateDomainDecompositionScatters(DM dm,PetscInt n,DM *subdms,VecScatter **iscat,VecScatter **oscat,VecScatter **gscat)
1630: {
1636: if (dm->ops->createddscatters) {
1637: (*dm->ops->createddscatters)(dm,n,subdms,iscat,oscat,gscat);
1638: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "This type has no DMCreateDomainDecompositionScatter implementation defined");
1639: return(0);
1640: }
1644: /*@
1645: DMRefine - Refines a DM object
1647: Collective on DM
1649: Input Parameter:
1650: + dm - the DM object
1651: - comm - the communicator to contain the new DM object (or MPI_COMM_NULL)
1653: Output Parameter:
1654: . dmf - the refined DM, or NULL
1656: Note: If no refinement was done, the return value is NULL
1658: Level: developer
1660: .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation()
1661: @*/
1662: PetscErrorCode DMRefine(DM dm,MPI_Comm comm,DM *dmf)
1663: {
1664: PetscErrorCode ierr;
1665: DMRefineHookLink link;
1669: if (!dm->ops->refine) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"This DM cannot refine");
1670: (*dm->ops->refine)(dm,comm,dmf);
1671: if (*dmf) {
1672: (*dmf)->ops->creatematrix = dm->ops->creatematrix;
1674: PetscObjectCopyFortranFunctionPointers((PetscObject)dm,(PetscObject)*dmf);
1676: (*dmf)->ctx = dm->ctx;
1677: (*dmf)->leveldown = dm->leveldown;
1678: (*dmf)->levelup = dm->levelup + 1;
1680: DMSetMatType(*dmf,dm->mattype);
1681: for (link=dm->refinehook; link; link=link->next) {
1682: if (link->refinehook) {
1683: (*link->refinehook)(dm,*dmf,link->ctx);
1684: }
1685: }
1686: }
1687: return(0);
1688: }
1692: /*@C
1693: DMRefineHookAdd - adds a callback to be run when interpolating a nonlinear problem to a finer grid
1695: Logically Collective
1697: Input Arguments:
1698: + coarse - nonlinear solver context on which to run a hook when restricting to a coarser level
1699: . refinehook - function to run when setting up a coarser level
1700: . interphook - function to run to update data on finer levels (once per SNESSolve())
1701: - ctx - [optional] user-defined context for provide data for the hooks (may be NULL)
1703: Calling sequence of refinehook:
1704: $ refinehook(DM coarse,DM fine,void *ctx);
1706: + coarse - coarse level DM
1707: . fine - fine level DM to interpolate problem to
1708: - ctx - optional user-defined function context
1710: Calling sequence for interphook:
1711: $ interphook(DM coarse,Mat interp,DM fine,void *ctx)
1713: + coarse - coarse level DM
1714: . interp - matrix interpolating a coarse-level solution to the finer grid
1715: . fine - fine level DM to update
1716: - ctx - optional user-defined function context
1718: Level: advanced
1720: Notes:
1721: This function is only needed if auxiliary data needs to be passed to fine grids while grid sequencing
1723: If this function is called multiple times, the hooks will be run in the order they are added.
1725: This function is currently not available from Fortran.
1727: .seealso: DMCoarsenHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
1728: @*/
1729: PetscErrorCode DMRefineHookAdd(DM coarse,PetscErrorCode (*refinehook)(DM,DM,void*),PetscErrorCode (*interphook)(DM,Mat,DM,void*),void *ctx)
1730: {
1731: PetscErrorCode ierr;
1732: DMRefineHookLink link,*p;
1736: for (p=&coarse->refinehook; *p; p=&(*p)->next) {} /* Scan to the end of the current list of hooks */
1737: PetscMalloc(sizeof(struct _DMRefineHookLink),&link);
1738: link->refinehook = refinehook;
1739: link->interphook = interphook;
1740: link->ctx = ctx;
1741: link->next = NULL;
1742: *p = link;
1743: return(0);
1744: }
1748: /*@
1749: DMInterpolate - interpolates user-defined problem data to a finer DM by running hooks registered by DMRefineHookAdd()
1751: Collective if any hooks are
1753: Input Arguments:
1754: + coarse - coarser DM to use as a base
1755: . restrct - interpolation matrix, apply using MatInterpolate()
1756: - fine - finer DM to update
1758: Level: developer
1760: .seealso: DMRefineHookAdd(), MatInterpolate()
1761: @*/
1762: PetscErrorCode DMInterpolate(DM coarse,Mat interp,DM fine)
1763: {
1764: PetscErrorCode ierr;
1765: DMRefineHookLink link;
1768: for (link=fine->refinehook; link; link=link->next) {
1769: if (link->interphook) {
1770: (*link->interphook)(coarse,interp,fine,link->ctx);
1771: }
1772: }
1773: return(0);
1774: }
1778: /*@
1779: DMGetRefineLevel - Get's the number of refinements that have generated this DM.
1781: Not Collective
1783: Input Parameter:
1784: . dm - the DM object
1786: Output Parameter:
1787: . level - number of refinements
1789: Level: developer
1791: .seealso DMCoarsen(), DMGetCoarsenLevel(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation()
1793: @*/
1794: PetscErrorCode DMGetRefineLevel(DM dm,PetscInt *level)
1795: {
1798: *level = dm->levelup;
1799: return(0);
1800: }
1804: /*@
1805: DMSetRefineLevel - Set's the number of refinements that have generated this DM.
1807: Not Collective
1809: Input Parameter:
1810: + dm - the DM object
1811: - level - number of refinements
1813: Level: advanced
1815: Notes: This value is used by PCMG to determine how many multigrid levels to use
1817: .seealso DMCoarsen(), DMGetCoarsenLevel(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation()
1819: @*/
1820: PetscErrorCode DMSetRefineLevel(DM dm,PetscInt level)
1821: {
1824: dm->levelup = level;
1825: return(0);
1826: }
1830: /*@C
1831: DMGlobalToLocalHookAdd - adds a callback to be run when global to local is called
1833: Logically Collective
1835: Input Arguments:
1836: + dm - the DM
1837: . beginhook - function to run at the beginning of DMGlobalToLocalBegin()
1838: . endhook - function to run after DMGlobalToLocalEnd() has completed
1839: - ctx - [optional] user-defined context for provide data for the hooks (may be NULL)
1841: Calling sequence for beginhook:
1842: $ beginhook(DM fine,VecScatter out,VecScatter in,DM coarse,void *ctx)
1844: + dm - global DM
1845: . g - global vector
1846: . mode - mode
1847: . l - local vector
1848: - ctx - optional user-defined function context
1851: Calling sequence for endhook:
1852: $ endhook(DM fine,VecScatter out,VecScatter in,DM coarse,void *ctx)
1854: + global - global DM
1855: - ctx - optional user-defined function context
1857: Level: advanced
1859: .seealso: DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
1860: @*/
1861: PetscErrorCode DMGlobalToLocalHookAdd(DM dm,PetscErrorCode (*beginhook)(DM,Vec,InsertMode,Vec,void*),PetscErrorCode (*endhook)(DM,Vec,InsertMode,Vec,void*),void *ctx)
1862: {
1863: PetscErrorCode ierr;
1864: DMGlobalToLocalHookLink link,*p;
1868: for (p=&dm->gtolhook; *p; p=&(*p)->next) {} /* Scan to the end of the current list of hooks */
1869: PetscMalloc(sizeof(struct _DMGlobalToLocalHookLink),&link);
1870: link->beginhook = beginhook;
1871: link->endhook = endhook;
1872: link->ctx = ctx;
1873: link->next = NULL;
1874: *p = link;
1875: return(0);
1876: }
1880: static PetscErrorCode DMGlobalToLocalHook_Constraints(DM dm, Vec g, InsertMode mode, Vec l, void *ctx)
1881: {
1882: Mat cMat;
1883: Vec cVec;
1884: PetscSection section, cSec;
1885: PetscInt pStart, pEnd, p, dof;
1890: DMGetDefaultConstraints(dm,&cSec,&cMat);
1891: if (cMat && (mode == INSERT_VALUES || mode == INSERT_ALL_VALUES || mode == INSERT_BC_VALUES)) {
1892: PetscInt nRows;
1894: MatGetSize(cMat,&nRows,NULL);
1895: if (nRows <= 0) return(0);
1896: DMGetDefaultSection(dm,§ion);
1897: MatCreateVecs(cMat,NULL,&cVec);
1898: MatMult(cMat,l,cVec);
1899: PetscSectionGetChart(cSec,&pStart,&pEnd);
1900: for (p = pStart; p < pEnd; p++) {
1901: PetscSectionGetDof(cSec,p,&dof);
1902: if (dof) {
1903: PetscScalar *vals;
1904: VecGetValuesSection(cVec,cSec,p,&vals);
1905: VecSetValuesSection(l,section,p,vals,INSERT_ALL_VALUES);
1906: }
1907: }
1908: VecDestroy(&cVec);
1909: }
1910: return(0);
1911: }
1915: /*@
1916: DMGlobalToLocalBegin - Begins updating local vectors from global vector
1918: Neighbor-wise Collective on DM
1920: Input Parameters:
1921: + dm - the DM object
1922: . g - the global vector
1923: . mode - INSERT_VALUES or ADD_VALUES
1924: - l - the local vector
1927: Level: beginner
1929: .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMGlobalToLocalEnd(), DMLocalToGlobalBegin()
1931: @*/
1932: PetscErrorCode DMGlobalToLocalBegin(DM dm,Vec g,InsertMode mode,Vec l)
1933: {
1934: PetscSF sf;
1935: PetscErrorCode ierr;
1936: DMGlobalToLocalHookLink link;
1940: for (link=dm->gtolhook; link; link=link->next) {
1941: if (link->beginhook) {
1942: (*link->beginhook)(dm,g,mode,l,link->ctx);
1943: }
1944: }
1945: DMGetDefaultSF(dm, &sf);
1946: if (sf) {
1947: const PetscScalar *gArray;
1948: PetscScalar *lArray;
1950: if (mode == ADD_VALUES) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode);
1951: VecGetArray(l, &lArray);
1952: VecGetArrayRead(g, &gArray);
1953: PetscSFBcastBegin(sf, MPIU_SCALAR, gArray, lArray);
1954: VecRestoreArray(l, &lArray);
1955: VecRestoreArrayRead(g, &gArray);
1956: } else {
1957: (*dm->ops->globaltolocalbegin)(dm,g,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),l);
1958: }
1959: return(0);
1960: }
1964: /*@
1965: DMGlobalToLocalEnd - Ends updating local vectors from global vector
1967: Neighbor-wise Collective on DM
1969: Input Parameters:
1970: + dm - the DM object
1971: . g - the global vector
1972: . mode - INSERT_VALUES or ADD_VALUES
1973: - l - the local vector
1976: Level: beginner
1978: .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMGlobalToLocalEnd(), DMLocalToGlobalBegin()
1980: @*/
1981: PetscErrorCode DMGlobalToLocalEnd(DM dm,Vec g,InsertMode mode,Vec l)
1982: {
1983: PetscSF sf;
1984: PetscErrorCode ierr;
1985: const PetscScalar *gArray;
1986: PetscScalar *lArray;
1987: DMGlobalToLocalHookLink link;
1991: DMGetDefaultSF(dm, &sf);
1992: if (sf) {
1993: if (mode == ADD_VALUES) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode);
1995: VecGetArray(l, &lArray);
1996: VecGetArrayRead(g, &gArray);
1997: PetscSFBcastEnd(sf, MPIU_SCALAR, gArray, lArray);
1998: VecRestoreArray(l, &lArray);
1999: VecRestoreArrayRead(g, &gArray);
2000: } else {
2001: (*dm->ops->globaltolocalend)(dm,g,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),l);
2002: }
2003: DMGlobalToLocalHook_Constraints(dm,g,mode,l,NULL);
2004: for (link=dm->gtolhook; link; link=link->next) {
2005: if (link->endhook) {(*link->endhook)(dm,g,mode,l,link->ctx);}
2006: }
2007: return(0);
2008: }
2012: /*@C
2013: DMLocalToGlobalHookAdd - adds a callback to be run when a local to global is called
2015: Logically Collective
2017: Input Arguments:
2018: + dm - the DM
2019: . beginhook - function to run at the beginning of DMLocalToGlobalBegin()
2020: . endhook - function to run after DMLocalToGlobalEnd() has completed
2021: - ctx - [optional] user-defined context for provide data for the hooks (may be NULL)
2023: Calling sequence for beginhook:
2024: $ beginhook(DM fine,Vec l,InsertMode mode,Vec g,void *ctx)
2026: + dm - global DM
2027: . l - local vector
2028: . mode - mode
2029: . g - global vector
2030: - ctx - optional user-defined function context
2033: Calling sequence for endhook:
2034: $ endhook(DM fine,Vec l,InsertMode mode,Vec g,void *ctx)
2036: + global - global DM
2037: . l - local vector
2038: . mode - mode
2039: . g - global vector
2040: - ctx - optional user-defined function context
2042: Level: advanced
2044: .seealso: DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2045: @*/
2046: PetscErrorCode DMLocalToGlobalHookAdd(DM dm,PetscErrorCode (*beginhook)(DM,Vec,InsertMode,Vec,void*),PetscErrorCode (*endhook)(DM,Vec,InsertMode,Vec,void*),void *ctx)
2047: {
2048: PetscErrorCode ierr;
2049: DMLocalToGlobalHookLink link,*p;
2053: for (p=&dm->ltoghook; *p; p=&(*p)->next) {} /* Scan to the end of the current list of hooks */
2054: PetscMalloc(sizeof(struct _DMLocalToGlobalHookLink),&link);
2055: link->beginhook = beginhook;
2056: link->endhook = endhook;
2057: link->ctx = ctx;
2058: link->next = NULL;
2059: *p = link;
2060: return(0);
2061: }
2065: static PetscErrorCode DMLocalToGlobalHook_Constraints(DM dm, Vec l, InsertMode mode, Vec g, void *ctx)
2066: {
2067: Mat cMat;
2068: Vec cVec;
2069: PetscSection section, cSec;
2070: PetscInt pStart, pEnd, p, dof;
2075: DMGetDefaultConstraints(dm,&cSec,&cMat);
2076: if (cMat && (mode == ADD_VALUES || mode == ADD_ALL_VALUES || mode == ADD_BC_VALUES)) {
2077: PetscInt nRows;
2079: MatGetSize(cMat,&nRows,NULL);
2080: if (nRows <= 0) return(0);
2081: DMGetDefaultSection(dm,§ion);
2082: MatCreateVecs(cMat,NULL,&cVec);
2083: PetscSectionGetChart(cSec,&pStart,&pEnd);
2084: for (p = pStart; p < pEnd; p++) {
2085: PetscSectionGetDof(cSec,p,&dof);
2086: if (dof) {
2087: PetscInt d;
2088: PetscScalar *vals;
2089: VecGetValuesSection(l,section,p,&vals);
2090: VecSetValuesSection(cVec,cSec,p,vals,mode);
2091: /* for this to be the true transpose, we have to zero the values that
2092: * we just extracted */
2093: for (d = 0; d < dof; d++) {
2094: vals[d] = 0.;
2095: }
2096: }
2097: }
2098: MatMultTransposeAdd(cMat,cVec,l,l);
2099: VecDestroy(&cVec);
2100: }
2101: return(0);
2102: }
2106: /*@
2107: DMLocalToGlobalBegin - updates global vectors from local vectors
2109: Neighbor-wise Collective on DM
2111: Input Parameters:
2112: + dm - the DM object
2113: . l - the local vector
2114: . 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.
2115: - g - the global vector
2117: Notes: In the ADD_VALUES case you normally would zero the receiving vector before beginning this operation.
2118: INSERT_VALUES is not supported for DMDA, in that case simply compute the values directly into a global vector instead of a local one.
2120: Level: beginner
2122: .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMGlobalToLocalEnd(), DMGlobalToLocalBegin()
2124: @*/
2125: PetscErrorCode DMLocalToGlobalBegin(DM dm,Vec l,InsertMode mode,Vec g)
2126: {
2127: PetscSF sf;
2128: PetscSection s, gs;
2129: DMLocalToGlobalHookLink link;
2130: const PetscScalar *lArray;
2131: PetscScalar *gArray;
2132: PetscBool isInsert;
2133: PetscErrorCode ierr;
2137: for (link=dm->ltoghook; link; link=link->next) {
2138: if (link->beginhook) {
2139: (*link->beginhook)(dm,l,mode,g,link->ctx);
2140: }
2141: }
2142: DMLocalToGlobalHook_Constraints(dm,l,mode,g,NULL);
2143: DMGetDefaultSF(dm, &sf);
2144: DMGetDefaultSection(dm, &s);
2145: switch (mode) {
2146: case INSERT_VALUES:
2147: case INSERT_ALL_VALUES:
2148: isInsert = PETSC_TRUE; break;
2149: case ADD_VALUES:
2150: case ADD_ALL_VALUES:
2151: isInsert = PETSC_FALSE; break;
2152: default:
2153: SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode);
2154: }
2155: if (sf && !isInsert) {
2156: VecGetArrayRead(l, &lArray);
2157: VecGetArray(g, &gArray);
2158: PetscSFReduceBegin(sf, MPIU_SCALAR, lArray, gArray, MPIU_SUM);
2159: VecRestoreArrayRead(l, &lArray);
2160: VecRestoreArray(g, &gArray);
2161: } else if (s && isInsert) {
2162: PetscInt gStart, pStart, pEnd, p;
2164: DMGetDefaultGlobalSection(dm, &gs);
2165: PetscSectionGetChart(s, &pStart, &pEnd);
2166: VecGetOwnershipRange(g, &gStart, NULL);
2167: VecGetArrayRead(l, &lArray);
2168: VecGetArray(g, &gArray);
2169: for (p = pStart; p < pEnd; ++p) {
2170: PetscInt dof, gdof, cdof, gcdof, off, goff, d, e;
2172: PetscSectionGetDof(s, p, &dof);
2173: PetscSectionGetDof(gs, p, &gdof);
2174: PetscSectionGetConstraintDof(s, p, &cdof);
2175: PetscSectionGetConstraintDof(gs, p, &gcdof);
2176: PetscSectionGetOffset(s, p, &off);
2177: PetscSectionGetOffset(gs, p, &goff);
2178: /* Ignore off-process data and points with no global data */
2179: if (!gdof || goff < 0) continue;
2180: 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);
2181: /* If no constraints are enforced in the global vector */
2182: if (!gcdof) {
2183: for (d = 0; d < dof; ++d) gArray[goff-gStart+d] = lArray[off+d];
2184: /* If constraints are enforced in the global vector */
2185: } else if (cdof == gcdof) {
2186: const PetscInt *cdofs;
2187: PetscInt cind = 0;
2189: PetscSectionGetConstraintIndices(s, p, &cdofs);
2190: for (d = 0, e = 0; d < dof; ++d) {
2191: if ((cind < cdof) && (d == cdofs[cind])) {++cind; continue;}
2192: gArray[goff-gStart+e++] = lArray[off+d];
2193: }
2194: } 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);
2195: }
2196: VecRestoreArrayRead(l, &lArray);
2197: VecRestoreArray(g, &gArray);
2198: } else {
2199: (*dm->ops->localtoglobalbegin)(dm,l,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),g);
2200: }
2201: return(0);
2202: }
2206: /*@
2207: DMLocalToGlobalEnd - updates global vectors from local vectors
2209: Neighbor-wise Collective on DM
2211: Input Parameters:
2212: + dm - the DM object
2213: . l - the local vector
2214: . mode - INSERT_VALUES or ADD_VALUES
2215: - g - the global vector
2218: Level: beginner
2220: .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMGlobalToLocalEnd(), DMGlobalToLocalEnd()
2222: @*/
2223: PetscErrorCode DMLocalToGlobalEnd(DM dm,Vec l,InsertMode mode,Vec g)
2224: {
2225: PetscSF sf;
2226: PetscSection s;
2227: DMLocalToGlobalHookLink link;
2228: PetscBool isInsert;
2229: PetscErrorCode ierr;
2233: DMGetDefaultSF(dm, &sf);
2234: DMGetDefaultSection(dm, &s);
2235: switch (mode) {
2236: case INSERT_VALUES:
2237: case INSERT_ALL_VALUES:
2238: isInsert = PETSC_TRUE; break;
2239: case ADD_VALUES:
2240: case ADD_ALL_VALUES:
2241: isInsert = PETSC_FALSE; break;
2242: default:
2243: SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode);
2244: }
2245: if (sf && !isInsert) {
2246: const PetscScalar *lArray;
2247: PetscScalar *gArray;
2249: VecGetArrayRead(l, &lArray);
2250: VecGetArray(g, &gArray);
2251: PetscSFReduceEnd(sf, MPIU_SCALAR, lArray, gArray, MPIU_SUM);
2252: VecRestoreArrayRead(l, &lArray);
2253: VecRestoreArray(g, &gArray);
2254: } else if (s && isInsert) {
2255: } else {
2256: (*dm->ops->localtoglobalend)(dm,l,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),g);
2257: }
2258: for (link=dm->ltoghook; link; link=link->next) {
2259: if (link->endhook) {(*link->endhook)(dm,g,mode,l,link->ctx);}
2260: }
2261: return(0);
2262: }
2266: /*@
2267: DMLocalToLocalBegin - Maps from a local vector (including ghost points
2268: that contain irrelevant values) to another local vector where the ghost
2269: points in the second are set correctly. Must be followed by DMLocalToLocalEnd().
2271: Neighbor-wise Collective on DM and Vec
2273: Input Parameters:
2274: + dm - the DM object
2275: . g - the original local vector
2276: - mode - one of INSERT_VALUES or ADD_VALUES
2278: Output Parameter:
2279: . l - the local vector with correct ghost values
2281: Level: intermediate
2283: Notes:
2284: The local vectors used here need not be the same as those
2285: obtained from DMCreateLocalVector(), BUT they
2286: must have the same parallel data layout; they could, for example, be
2287: obtained with VecDuplicate() from the DM originating vectors.
2289: .keywords: DM, local-to-local, begin
2290: .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateLocalVector(), DMCreateGlobalVector(), DMCreateInterpolation(), DMLocalToLocalEnd(), DMGlobalToLocalEnd(), DMLocalToGlobalBegin()
2292: @*/
2293: PetscErrorCode DMLocalToLocalBegin(DM dm,Vec g,InsertMode mode,Vec l)
2294: {
2295: PetscErrorCode ierr;
2299: (*dm->ops->localtolocalbegin)(dm,g,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),l);
2300: return(0);
2301: }
2305: /*@
2306: DMLocalToLocalEnd - Maps from a local vector (including ghost points
2307: that contain irrelevant values) to another local vector where the ghost
2308: points in the second are set correctly. Must be preceded by DMLocalToLocalBegin().
2310: Neighbor-wise Collective on DM and Vec
2312: Input Parameters:
2313: + da - the DM object
2314: . g - the original local vector
2315: - mode - one of INSERT_VALUES or ADD_VALUES
2317: Output Parameter:
2318: . l - the local vector with correct ghost values
2320: Level: intermediate
2322: Notes:
2323: The local vectors used here need not be the same as those
2324: obtained from DMCreateLocalVector(), BUT they
2325: must have the same parallel data layout; they could, for example, be
2326: obtained with VecDuplicate() from the DM originating vectors.
2328: .keywords: DM, local-to-local, end
2329: .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateLocalVector(), DMCreateGlobalVector(), DMCreateInterpolation(), DMLocalToLocalBegin(), DMGlobalToLocalEnd(), DMLocalToGlobalBegin()
2331: @*/
2332: PetscErrorCode DMLocalToLocalEnd(DM dm,Vec g,InsertMode mode,Vec l)
2333: {
2334: PetscErrorCode ierr;
2338: (*dm->ops->localtolocalend)(dm,g,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),l);
2339: return(0);
2340: }
2345: /*@
2346: DMCoarsen - Coarsens a DM object
2348: Collective on DM
2350: Input Parameter:
2351: + dm - the DM object
2352: - comm - the communicator to contain the new DM object (or MPI_COMM_NULL)
2354: Output Parameter:
2355: . dmc - the coarsened DM
2357: Level: developer
2359: .seealso DMRefine(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation()
2361: @*/
2362: PetscErrorCode DMCoarsen(DM dm, MPI_Comm comm, DM *dmc)
2363: {
2364: PetscErrorCode ierr;
2365: DMCoarsenHookLink link;
2369: if (!dm->ops->coarsen) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"This DM cannot coarsen");
2370: PetscLogEventBegin(DM_Coarsen,dm,0,0,0);
2371: (*dm->ops->coarsen)(dm, comm, dmc);
2372: if (!(*dmc)) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "NULL coarse mesh produced");
2373: DMSetCoarseDM(dm,*dmc);
2374: (*dmc)->ops->creatematrix = dm->ops->creatematrix;
2375: PetscObjectCopyFortranFunctionPointers((PetscObject)dm,(PetscObject)*dmc);
2376: (*dmc)->ctx = dm->ctx;
2377: (*dmc)->levelup = dm->levelup;
2378: (*dmc)->leveldown = dm->leveldown + 1;
2379: DMSetMatType(*dmc,dm->mattype);
2380: for (link=dm->coarsenhook; link; link=link->next) {
2381: if (link->coarsenhook) {(*link->coarsenhook)(dm,*dmc,link->ctx);}
2382: }
2383: PetscLogEventEnd(DM_Coarsen,dm,0,0,0);
2384: return(0);
2385: }
2389: /*@C
2390: DMCoarsenHookAdd - adds a callback to be run when restricting a nonlinear problem to the coarse grid
2392: Logically Collective
2394: Input Arguments:
2395: + fine - nonlinear solver context on which to run a hook when restricting to a coarser level
2396: . coarsenhook - function to run when setting up a coarser level
2397: . restricthook - function to run to update data on coarser levels (once per SNESSolve())
2398: - ctx - [optional] user-defined context for provide data for the hooks (may be NULL)
2400: Calling sequence of coarsenhook:
2401: $ coarsenhook(DM fine,DM coarse,void *ctx);
2403: + fine - fine level DM
2404: . coarse - coarse level DM to restrict problem to
2405: - ctx - optional user-defined function context
2407: Calling sequence for restricthook:
2408: $ restricthook(DM fine,Mat mrestrict,Vec rscale,Mat inject,DM coarse,void *ctx)
2410: + fine - fine level DM
2411: . mrestrict - matrix restricting a fine-level solution to the coarse grid
2412: . rscale - scaling vector for restriction
2413: . inject - matrix restricting by injection
2414: . coarse - coarse level DM to update
2415: - ctx - optional user-defined function context
2417: Level: advanced
2419: Notes:
2420: This function is only needed if auxiliary data needs to be set up on coarse grids.
2422: If this function is called multiple times, the hooks will be run in the order they are added.
2424: In order to compose with nonlinear preconditioning without duplicating storage, the hook should be implemented to
2425: extract the finest level information from its context (instead of from the SNES).
2427: This function is currently not available from Fortran.
2429: .seealso: DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2430: @*/
2431: PetscErrorCode DMCoarsenHookAdd(DM fine,PetscErrorCode (*coarsenhook)(DM,DM,void*),PetscErrorCode (*restricthook)(DM,Mat,Vec,Mat,DM,void*),void *ctx)
2432: {
2433: PetscErrorCode ierr;
2434: DMCoarsenHookLink link,*p;
2438: for (p=&fine->coarsenhook; *p; p=&(*p)->next) {} /* Scan to the end of the current list of hooks */
2439: PetscMalloc(sizeof(struct _DMCoarsenHookLink),&link);
2440: link->coarsenhook = coarsenhook;
2441: link->restricthook = restricthook;
2442: link->ctx = ctx;
2443: link->next = NULL;
2444: *p = link;
2445: return(0);
2446: }
2450: /*@
2451: DMRestrict - restricts user-defined problem data to a coarser DM by running hooks registered by DMCoarsenHookAdd()
2453: Collective if any hooks are
2455: Input Arguments:
2456: + fine - finer DM to use as a base
2457: . restrct - restriction matrix, apply using MatRestrict()
2458: . inject - injection matrix, also use MatRestrict()
2459: - coarse - coarer DM to update
2461: Level: developer
2463: .seealso: DMCoarsenHookAdd(), MatRestrict()
2464: @*/
2465: PetscErrorCode DMRestrict(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse)
2466: {
2467: PetscErrorCode ierr;
2468: DMCoarsenHookLink link;
2471: for (link=fine->coarsenhook; link; link=link->next) {
2472: if (link->restricthook) {
2473: (*link->restricthook)(fine,restrct,rscale,inject,coarse,link->ctx);
2474: }
2475: }
2476: return(0);
2477: }
2481: /*@C
2482: DMSubDomainHookAdd - adds a callback to be run when restricting a problem to the coarse grid
2484: Logically Collective
2486: Input Arguments:
2487: + global - global DM
2488: . ddhook - function to run to pass data to the decomposition DM upon its creation
2489: . restricthook - function to run to update data on block solve (at the beginning of the block solve)
2490: - ctx - [optional] user-defined context for provide data for the hooks (may be NULL)
2493: Calling sequence for ddhook:
2494: $ ddhook(DM global,DM block,void *ctx)
2496: + global - global DM
2497: . block - block DM
2498: - ctx - optional user-defined function context
2500: Calling sequence for restricthook:
2501: $ restricthook(DM global,VecScatter out,VecScatter in,DM block,void *ctx)
2503: + global - global DM
2504: . out - scatter to the outer (with ghost and overlap points) block vector
2505: . in - scatter to block vector values only owned locally
2506: . block - block DM
2507: - ctx - optional user-defined function context
2509: Level: advanced
2511: Notes:
2512: This function is only needed if auxiliary data needs to be set up on subdomain DMs.
2514: If this function is called multiple times, the hooks will be run in the order they are added.
2516: In order to compose with nonlinear preconditioning without duplicating storage, the hook should be implemented to
2517: extract the global information from its context (instead of from the SNES).
2519: This function is currently not available from Fortran.
2521: .seealso: DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2522: @*/
2523: PetscErrorCode DMSubDomainHookAdd(DM global,PetscErrorCode (*ddhook)(DM,DM,void*),PetscErrorCode (*restricthook)(DM,VecScatter,VecScatter,DM,void*),void *ctx)
2524: {
2525: PetscErrorCode ierr;
2526: DMSubDomainHookLink link,*p;
2530: for (p=&global->subdomainhook; *p; p=&(*p)->next) {} /* Scan to the end of the current list of hooks */
2531: PetscMalloc(sizeof(struct _DMSubDomainHookLink),&link);
2532: link->restricthook = restricthook;
2533: link->ddhook = ddhook;
2534: link->ctx = ctx;
2535: link->next = NULL;
2536: *p = link;
2537: return(0);
2538: }
2542: /*@
2543: DMSubDomainRestrict - restricts user-defined problem data to a block DM by running hooks registered by DMSubDomainHookAdd()
2545: Collective if any hooks are
2547: Input Arguments:
2548: + fine - finer DM to use as a base
2549: . oscatter - scatter from domain global vector filling subdomain global vector with overlap
2550: . gscatter - scatter from domain global vector filling subdomain local vector with ghosts
2551: - coarse - coarer DM to update
2553: Level: developer
2555: .seealso: DMCoarsenHookAdd(), MatRestrict()
2556: @*/
2557: PetscErrorCode DMSubDomainRestrict(DM global,VecScatter oscatter,VecScatter gscatter,DM subdm)
2558: {
2559: PetscErrorCode ierr;
2560: DMSubDomainHookLink link;
2563: for (link=global->subdomainhook; link; link=link->next) {
2564: if (link->restricthook) {
2565: (*link->restricthook)(global,oscatter,gscatter,subdm,link->ctx);
2566: }
2567: }
2568: return(0);
2569: }
2573: /*@
2574: DMGetCoarsenLevel - Get's the number of coarsenings that have generated this DM.
2576: Not Collective
2578: Input Parameter:
2579: . dm - the DM object
2581: Output Parameter:
2582: . level - number of coarsenings
2584: Level: developer
2586: .seealso DMCoarsen(), DMGetRefineLevel(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation()
2588: @*/
2589: PetscErrorCode DMGetCoarsenLevel(DM dm,PetscInt *level)
2590: {
2593: *level = dm->leveldown;
2594: return(0);
2595: }
2601: /*@C
2602: DMRefineHierarchy - Refines a DM object, all levels at once
2604: Collective on DM
2606: Input Parameter:
2607: + dm - the DM object
2608: - nlevels - the number of levels of refinement
2610: Output Parameter:
2611: . dmf - the refined DM hierarchy
2613: Level: developer
2615: .seealso DMCoarsenHierarchy(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation()
2617: @*/
2618: PetscErrorCode DMRefineHierarchy(DM dm,PetscInt nlevels,DM dmf[])
2619: {
2624: if (nlevels < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"nlevels cannot be negative");
2625: if (nlevels == 0) return(0);
2626: if (dm->ops->refinehierarchy) {
2627: (*dm->ops->refinehierarchy)(dm,nlevels,dmf);
2628: } else if (dm->ops->refine) {
2629: PetscInt i;
2631: DMRefine(dm,PetscObjectComm((PetscObject)dm),&dmf[0]);
2632: for (i=1; i<nlevels; i++) {
2633: DMRefine(dmf[i-1],PetscObjectComm((PetscObject)dm),&dmf[i]);
2634: }
2635: } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"No RefineHierarchy for this DM yet");
2636: return(0);
2637: }
2641: /*@C
2642: DMCoarsenHierarchy - Coarsens a DM object, all levels at once
2644: Collective on DM
2646: Input Parameter:
2647: + dm - the DM object
2648: - nlevels - the number of levels of coarsening
2650: Output Parameter:
2651: . dmc - the coarsened DM hierarchy
2653: Level: developer
2655: .seealso DMRefineHierarchy(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation()
2657: @*/
2658: PetscErrorCode DMCoarsenHierarchy(DM dm, PetscInt nlevels, DM dmc[])
2659: {
2664: if (nlevels < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"nlevels cannot be negative");
2665: if (nlevels == 0) return(0);
2667: if (dm->ops->coarsenhierarchy) {
2668: (*dm->ops->coarsenhierarchy)(dm, nlevels, dmc);
2669: } else if (dm->ops->coarsen) {
2670: PetscInt i;
2672: DMCoarsen(dm,PetscObjectComm((PetscObject)dm),&dmc[0]);
2673: for (i=1; i<nlevels; i++) {
2674: DMCoarsen(dmc[i-1],PetscObjectComm((PetscObject)dm),&dmc[i]);
2675: }
2676: } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"No CoarsenHierarchy for this DM yet");
2677: return(0);
2678: }
2682: /*@
2683: DMCreateAggregates - Gets the aggregates that map between
2684: grids associated with two DMs.
2686: Collective on DM
2688: Input Parameters:
2689: + dmc - the coarse grid DM
2690: - dmf - the fine grid DM
2692: Output Parameters:
2693: . rest - the restriction matrix (transpose of the projection matrix)
2695: Level: intermediate
2697: .keywords: interpolation, restriction, multigrid
2699: .seealso: DMRefine(), DMCreateInjection(), DMCreateInterpolation()
2700: @*/
2701: PetscErrorCode DMCreateAggregates(DM dmc, DM dmf, Mat *rest)
2702: {
2708: (*dmc->ops->getaggregates)(dmc, dmf, rest);
2709: return(0);
2710: }
2714: /*@C
2715: DMSetApplicationContextDestroy - Sets a user function that will be called to destroy the application context when the DM is destroyed
2717: Not Collective
2719: Input Parameters:
2720: + dm - the DM object
2721: - destroy - the destroy function
2723: Level: intermediate
2725: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext()
2727: @*/
2728: PetscErrorCode DMSetApplicationContextDestroy(DM dm,PetscErrorCode (*destroy)(void**))
2729: {
2732: dm->ctxdestroy = destroy;
2733: return(0);
2734: }
2738: /*@
2739: DMSetApplicationContext - Set a user context into a DM object
2741: Not Collective
2743: Input Parameters:
2744: + dm - the DM object
2745: - ctx - the user context
2747: Level: intermediate
2749: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext()
2751: @*/
2752: PetscErrorCode DMSetApplicationContext(DM dm,void *ctx)
2753: {
2756: dm->ctx = ctx;
2757: return(0);
2758: }
2762: /*@
2763: DMGetApplicationContext - Gets a user context from a DM object
2765: Not Collective
2767: Input Parameter:
2768: . dm - the DM object
2770: Output Parameter:
2771: . ctx - the user context
2773: Level: intermediate
2775: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext()
2777: @*/
2778: PetscErrorCode DMGetApplicationContext(DM dm,void *ctx)
2779: {
2782: *(void**)ctx = dm->ctx;
2783: return(0);
2784: }
2788: /*@C
2789: DMSetVariableBounds - sets a function to compute the lower and upper bound vectors for SNESVI.
2791: Logically Collective on DM
2793: Input Parameter:
2794: + dm - the DM object
2795: - f - the function that computes variable bounds used by SNESVI (use NULL to cancel a previous function that was set)
2797: Level: intermediate
2799: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext(),
2800: DMSetJacobian()
2802: @*/
2803: PetscErrorCode DMSetVariableBounds(DM dm,PetscErrorCode (*f)(DM,Vec,Vec))
2804: {
2806: dm->ops->computevariablebounds = f;
2807: return(0);
2808: }
2812: /*@
2813: DMHasVariableBounds - does the DM object have a variable bounds function?
2815: Not Collective
2817: Input Parameter:
2818: . dm - the DM object to destroy
2820: Output Parameter:
2821: . flg - PETSC_TRUE if the variable bounds function exists
2823: Level: developer
2825: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext()
2827: @*/
2828: PetscErrorCode DMHasVariableBounds(DM dm,PetscBool *flg)
2829: {
2831: *flg = (dm->ops->computevariablebounds) ? PETSC_TRUE : PETSC_FALSE;
2832: return(0);
2833: }
2837: /*@C
2838: DMComputeVariableBounds - compute variable bounds used by SNESVI.
2840: Logically Collective on DM
2842: Input Parameters:
2843: . dm - the DM object
2845: Output parameters:
2846: + xl - lower bound
2847: - xu - upper bound
2849: Level: advanced
2851: Notes: This is generally not called by users. It calls the function provided by the user with DMSetVariableBounds()
2853: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext()
2855: @*/
2856: PetscErrorCode DMComputeVariableBounds(DM dm, Vec xl, Vec xu)
2857: {
2863: if (dm->ops->computevariablebounds) {
2864: (*dm->ops->computevariablebounds)(dm, xl,xu);
2865: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "This DM is incapable of computing variable bounds.");
2866: return(0);
2867: }
2871: /*@
2872: DMHasColoring - does the DM object have a method of providing a coloring?
2874: Not Collective
2876: Input Parameter:
2877: . dm - the DM object
2879: Output Parameter:
2880: . flg - PETSC_TRUE if the DM has facilities for DMCreateColoring().
2882: Level: developer
2884: .seealso DMHasFunction(), DMCreateColoring()
2886: @*/
2887: PetscErrorCode DMHasColoring(DM dm,PetscBool *flg)
2888: {
2890: *flg = (dm->ops->getcoloring) ? PETSC_TRUE : PETSC_FALSE;
2891: return(0);
2892: }
2896: /*@
2897: DMHasCreateRestriction - does the DM object have a method of providing a restriction?
2899: Not Collective
2901: Input Parameter:
2902: . dm - the DM object
2904: Output Parameter:
2905: . flg - PETSC_TRUE if the DM has facilities for DMCreateRestriction().
2907: Level: developer
2909: .seealso DMHasFunction(), DMCreateRestriction()
2911: @*/
2912: PetscErrorCode DMHasCreateRestriction(DM dm,PetscBool *flg)
2913: {
2915: *flg = (dm->ops->createrestriction) ? PETSC_TRUE : PETSC_FALSE;
2916: return(0);
2917: }
2919: #undef __FUNCT__
2921: /*@C
2922: DMSetVec - set the vector at which to compute residual, Jacobian and VI bounds, if the problem is nonlinear.
2924: Collective on DM
2926: Input Parameter:
2927: + dm - the DM object
2928: - x - location to compute residual and Jacobian, if NULL is passed to those routines; will be NULL for linear problems.
2930: Level: developer
2932: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext()
2934: @*/
2935: PetscErrorCode DMSetVec(DM dm,Vec x)
2936: {
2940: if (x) {
2941: if (!dm->x) {
2942: DMCreateGlobalVector(dm,&dm->x);
2943: }
2944: VecCopy(x,dm->x);
2945: } else if (dm->x) {
2946: VecDestroy(&dm->x);
2947: }
2948: return(0);
2949: }
2951: PetscFunctionList DMList = NULL;
2952: PetscBool DMRegisterAllCalled = PETSC_FALSE;
2956: /*@C
2957: DMSetType - Builds a DM, for a particular DM implementation.
2959: Collective on DM
2961: Input Parameters:
2962: + dm - The DM object
2963: - method - The name of the DM type
2965: Options Database Key:
2966: . -dm_type <type> - Sets the DM type; use -help for a list of available types
2968: Notes:
2969: See "petsc/include/petscdm.h" for available DM types (for instance, DM1D, DM2D, or DM3D).
2971: Level: intermediate
2973: .keywords: DM, set, type
2974: .seealso: DMGetType(), DMCreate()
2975: @*/
2976: PetscErrorCode DMSetType(DM dm, DMType method)
2977: {
2978: PetscErrorCode (*r)(DM);
2979: PetscBool match;
2984: PetscObjectTypeCompare((PetscObject) dm, method, &match);
2985: if (match) return(0);
2987: DMRegisterAll();
2988: PetscFunctionListFind(DMList,method,&r);
2989: if (!r) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown DM type: %s", method);
2991: if (dm->ops->destroy) {
2992: (*dm->ops->destroy)(dm);
2993: dm->ops->destroy = NULL;
2994: }
2995: (*r)(dm);
2996: PetscObjectChangeTypeName((PetscObject)dm,method);
2997: return(0);
2998: }
3002: /*@C
3003: DMGetType - Gets the DM type name (as a string) from the DM.
3005: Not Collective
3007: Input Parameter:
3008: . dm - The DM
3010: Output Parameter:
3011: . type - The DM type name
3013: Level: intermediate
3015: .keywords: DM, get, type, name
3016: .seealso: DMSetType(), DMCreate()
3017: @*/
3018: PetscErrorCode DMGetType(DM dm, DMType *type)
3019: {
3025: DMRegisterAll();
3026: *type = ((PetscObject)dm)->type_name;
3027: return(0);
3028: }
3032: /*@C
3033: DMConvert - Converts a DM to another DM, either of the same or different type.
3035: Collective on DM
3037: Input Parameters:
3038: + dm - the DM
3039: - newtype - new DM type (use "same" for the same type)
3041: Output Parameter:
3042: . M - pointer to new DM
3044: Notes:
3045: Cannot be used to convert a sequential DM to parallel or parallel to sequential,
3046: the MPI communicator of the generated DM is always the same as the communicator
3047: of the input DM.
3049: Level: intermediate
3051: .seealso: DMCreate()
3052: @*/
3053: PetscErrorCode DMConvert(DM dm, DMType newtype, DM *M)
3054: {
3055: DM B;
3056: char convname[256];
3057: PetscBool sametype/*, issame */;
3064: PetscObjectTypeCompare((PetscObject) dm, newtype, &sametype);
3065: /* PetscStrcmp(newtype, "same", &issame); */
3066: if (sametype) {
3067: *M = dm;
3068: PetscObjectReference((PetscObject) dm);
3069: return(0);
3070: } else {
3071: PetscErrorCode (*conv)(DM, DMType, DM*) = NULL;
3073: /*
3074: Order of precedence:
3075: 1) See if a specialized converter is known to the current DM.
3076: 2) See if a specialized converter is known to the desired DM class.
3077: 3) See if a good general converter is registered for the desired class
3078: 4) See if a good general converter is known for the current matrix.
3079: 5) Use a really basic converter.
3080: */
3082: /* 1) See if a specialized converter is known to the current DM and the desired class */
3083: PetscStrcpy(convname,"DMConvert_");
3084: PetscStrcat(convname,((PetscObject) dm)->type_name);
3085: PetscStrcat(convname,"_");
3086: PetscStrcat(convname,newtype);
3087: PetscStrcat(convname,"_C");
3088: PetscObjectQueryFunction((PetscObject)dm,convname,&conv);
3089: if (conv) goto foundconv;
3091: /* 2) See if a specialized converter is known to the desired DM class. */
3092: DMCreate(PetscObjectComm((PetscObject)dm), &B);
3093: DMSetType(B, newtype);
3094: PetscStrcpy(convname,"DMConvert_");
3095: PetscStrcat(convname,((PetscObject) dm)->type_name);
3096: PetscStrcat(convname,"_");
3097: PetscStrcat(convname,newtype);
3098: PetscStrcat(convname,"_C");
3099: PetscObjectQueryFunction((PetscObject)B,convname,&conv);
3100: if (conv) {
3101: DMDestroy(&B);
3102: goto foundconv;
3103: }
3105: #if 0
3106: /* 3) See if a good general converter is registered for the desired class */
3107: conv = B->ops->convertfrom;
3108: DMDestroy(&B);
3109: if (conv) goto foundconv;
3111: /* 4) See if a good general converter is known for the current matrix */
3112: if (dm->ops->convert) {
3113: conv = dm->ops->convert;
3114: }
3115: if (conv) goto foundconv;
3116: #endif
3118: /* 5) Use a really basic converter. */
3119: SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "No conversion possible between DM types %s and %s", ((PetscObject) dm)->type_name, newtype);
3121: foundconv:
3122: PetscLogEventBegin(DM_Convert,dm,0,0,0);
3123: (*conv)(dm,newtype,M);
3124: /* Things that are independent of DM type: We should consult DMClone() here */
3125: if (dm->maxCell) {
3126: const PetscReal *maxCell, *L;
3127: const DMBoundaryType *bd;
3128: DMGetPeriodicity(dm, &maxCell, &L, &bd);
3129: DMSetPeriodicity(*M, maxCell, L, bd);
3130: }
3131: PetscLogEventEnd(DM_Convert,dm,0,0,0);
3132: }
3133: PetscObjectStateIncrease((PetscObject) *M);
3134: return(0);
3135: }
3137: /*--------------------------------------------------------------------------------------------------------------------*/
3141: /*@C
3142: DMRegister - Adds a new DM component implementation
3144: Not Collective
3146: Input Parameters:
3147: + name - The name of a new user-defined creation routine
3148: - create_func - The creation routine itself
3150: Notes:
3151: DMRegister() may be called multiple times to add several user-defined DMs
3154: Sample usage:
3155: .vb
3156: DMRegister("my_da", MyDMCreate);
3157: .ve
3159: Then, your DM type can be chosen with the procedural interface via
3160: .vb
3161: DMCreate(MPI_Comm, DM *);
3162: DMSetType(DM,"my_da");
3163: .ve
3164: or at runtime via the option
3165: .vb
3166: -da_type my_da
3167: .ve
3169: Level: advanced
3171: .keywords: DM, register
3172: .seealso: DMRegisterAll(), DMRegisterDestroy()
3174: @*/
3175: PetscErrorCode DMRegister(const char sname[],PetscErrorCode (*function)(DM))
3176: {
3180: PetscFunctionListAdd(&DMList,sname,function);
3181: return(0);
3182: }
3186: /*@C
3187: DMLoad - Loads a DM that has been stored in binary with DMView().
3189: Collective on PetscViewer
3191: Input Parameters:
3192: + newdm - the newly loaded DM, this needs to have been created with DMCreate() or
3193: some related function before a call to DMLoad().
3194: - viewer - binary file viewer, obtained from PetscViewerBinaryOpen() or
3195: HDF5 file viewer, obtained from PetscViewerHDF5Open()
3197: Level: intermediate
3199: Notes:
3200: The type is determined by the data in the file, any type set into the DM before this call is ignored.
3202: Notes for advanced users:
3203: Most users should not need to know the details of the binary storage
3204: format, since DMLoad() and DMView() completely hide these details.
3205: But for anyone who's interested, the standard binary matrix storage
3206: format is
3207: .vb
3208: has not yet been determined
3209: .ve
3211: .seealso: PetscViewerBinaryOpen(), DMView(), MatLoad(), VecLoad()
3212: @*/
3213: PetscErrorCode DMLoad(DM newdm, PetscViewer viewer)
3214: {
3215: PetscBool isbinary, ishdf5;
3221: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
3222: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);
3223: if (isbinary) {
3224: PetscInt classid;
3225: char type[256];
3227: PetscViewerBinaryRead(viewer,&classid,1,NULL,PETSC_INT);
3228: if (classid != DM_FILE_CLASSID) SETERRQ1(PetscObjectComm((PetscObject)newdm),PETSC_ERR_ARG_WRONG,"Not DM next in file, classid found %d",(int)classid);
3229: PetscViewerBinaryRead(viewer,type,256,NULL,PETSC_CHAR);
3230: DMSetType(newdm, type);
3231: if (newdm->ops->load) {(*newdm->ops->load)(newdm,viewer);}
3232: } else if (ishdf5) {
3233: if (newdm->ops->load) {(*newdm->ops->load)(newdm,viewer);}
3234: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid viewer; open viewer with PetscViewerBinaryOpen() or PetscViewerHDF5Open()");
3235: return(0);
3236: }
3238: /******************************** FEM Support **********************************/
3242: PetscErrorCode DMPrintCellVector(PetscInt c, const char name[], PetscInt len, const PetscScalar x[])
3243: {
3244: PetscInt f;
3248: PetscPrintf(PETSC_COMM_SELF, "Cell %D Element %s\n", c, name);
3249: for (f = 0; f < len; ++f) {
3250: PetscPrintf(PETSC_COMM_SELF, " | %g |\n", (double)PetscRealPart(x[f]));
3251: }
3252: return(0);
3253: }
3257: PetscErrorCode DMPrintCellMatrix(PetscInt c, const char name[], PetscInt rows, PetscInt cols, const PetscScalar A[])
3258: {
3259: PetscInt f, g;
3263: PetscPrintf(PETSC_COMM_SELF, "Cell %D Element %s\n", c, name);
3264: for (f = 0; f < rows; ++f) {
3265: PetscPrintf(PETSC_COMM_SELF, " |");
3266: for (g = 0; g < cols; ++g) {
3267: PetscPrintf(PETSC_COMM_SELF, " % 9.5g", PetscRealPart(A[f*cols+g]));
3268: }
3269: PetscPrintf(PETSC_COMM_SELF, " |\n");
3270: }
3271: return(0);
3272: }
3276: PetscErrorCode DMPrintLocalVec(DM dm, const char name[], PetscReal tol, Vec X)
3277: {
3278: PetscMPIInt rank, numProcs;
3279: PetscInt p;
3283: MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
3284: MPI_Comm_size(PetscObjectComm((PetscObject) dm), &numProcs);
3285: PetscPrintf(PetscObjectComm((PetscObject) dm), "%s:\n", name);
3286: for (p = 0; p < numProcs; ++p) {
3287: if (p == rank) {
3288: Vec x;
3290: VecDuplicate(X, &x);
3291: VecCopy(X, x);
3292: VecChop(x, tol);
3293: VecView(x, PETSC_VIEWER_STDOUT_SELF);
3294: VecDestroy(&x);
3295: PetscViewerFlush(PETSC_VIEWER_STDOUT_SELF);
3296: }
3297: PetscBarrier((PetscObject) dm);
3298: }
3299: return(0);
3300: }
3304: /*@
3305: DMGetDefaultSection - Get the PetscSection encoding the local data layout for the DM.
3307: Input Parameter:
3308: . dm - The DM
3310: Output Parameter:
3311: . section - The PetscSection
3313: Level: intermediate
3315: Note: This gets a borrowed reference, so the user should not destroy this PetscSection.
3317: .seealso: DMSetDefaultSection(), DMGetDefaultGlobalSection()
3318: @*/
3319: PetscErrorCode DMGetDefaultSection(DM dm, PetscSection *section)
3320: {
3326: if (!dm->defaultSection && dm->ops->createdefaultsection) {
3327: (*dm->ops->createdefaultsection)(dm);
3328: if (dm->defaultSection) {PetscObjectViewFromOptions((PetscObject) dm->defaultSection, NULL, "-dm_petscsection_view");}
3329: }
3330: *section = dm->defaultSection;
3331: return(0);
3332: }
3336: /*@
3337: DMSetDefaultSection - Set the PetscSection encoding the local data layout for the DM.
3339: Input Parameters:
3340: + dm - The DM
3341: - section - The PetscSection
3343: Level: intermediate
3345: Note: Any existing Section will be destroyed
3347: .seealso: DMSetDefaultSection(), DMGetDefaultGlobalSection()
3348: @*/
3349: PetscErrorCode DMSetDefaultSection(DM dm, PetscSection section)
3350: {
3351: PetscInt numFields = 0;
3352: PetscInt f;
3357: if (section) {
3359: PetscObjectReference((PetscObject)section);
3360: }
3361: PetscSectionDestroy(&dm->defaultSection);
3362: dm->defaultSection = section;
3363: if (section) {PetscSectionGetNumFields(dm->defaultSection, &numFields);}
3364: if (numFields) {
3365: DMSetNumFields(dm, numFields);
3366: for (f = 0; f < numFields; ++f) {
3367: PetscObject disc;
3368: const char *name;
3370: PetscSectionGetFieldName(dm->defaultSection, f, &name);
3371: DMGetField(dm, f, &disc);
3372: PetscObjectSetName(disc, name);
3373: }
3374: }
3375: /* The global section will be rebuilt in the next call to DMGetDefaultGlobalSection(). */
3376: PetscSectionDestroy(&dm->defaultGlobalSection);
3377: return(0);
3378: }
3382: /*@
3383: DMGetDefaultConstraints - Get the PetscSection and Mat the specify the local constraint interpolation. See DMSetDefaultConstraints() for a description of the purpose of constraint interpolation.
3385: not collective
3387: Input Parameter:
3388: . dm - The DM
3390: Output Parameter:
3391: + 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.
3392: - 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.
3394: Level: advanced
3396: Note: This gets borrowed references, so the user should not destroy the PetscSection or the Mat.
3398: .seealso: DMSetDefaultConstraints()
3399: @*/
3400: PetscErrorCode DMGetDefaultConstraints(DM dm, PetscSection *section, Mat *mat)
3401: {
3406: if (!dm->defaultConstraintSection && !dm->defaultConstraintMat && dm->ops->createdefaultconstraints) {(*dm->ops->createdefaultconstraints)(dm);}
3407: if (section) {*section = dm->defaultConstraintSection;}
3408: if (mat) {*mat = dm->defaultConstraintMat;}
3409: return(0);
3410: }
3414: /*@
3415: DMSetDefaultConstraints - Set the PetscSection and Mat the specify the local constraint interpolation.
3417: 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().
3419: 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.
3421: collective on dm
3423: Input Parameters:
3424: + dm - The DM
3425: + 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).
3426: - 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).
3428: Level: advanced
3430: Note: This increments the references of the PetscSection and the Mat, so they user can destroy them
3432: .seealso: DMGetDefaultConstraints()
3433: @*/
3434: PetscErrorCode DMSetDefaultConstraints(DM dm, PetscSection section, Mat mat)
3435: {
3436: PetscMPIInt result;
3441: if (section) {
3443: MPI_Comm_compare(PETSC_COMM_SELF,PetscObjectComm((PetscObject)section),&result);
3444: if (result != MPI_CONGRUENT && result != MPI_IDENT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMECOMM,"constraint section must have local communicator");
3445: }
3446: if (mat) {
3448: MPI_Comm_compare(PETSC_COMM_SELF,PetscObjectComm((PetscObject)mat),&result);
3449: if (result != MPI_CONGRUENT && result != MPI_IDENT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMECOMM,"constraint matrix must have local communicator");
3450: }
3451: PetscObjectReference((PetscObject)section);
3452: PetscSectionDestroy(&dm->defaultConstraintSection);
3453: dm->defaultConstraintSection = section;
3454: PetscObjectReference((PetscObject)mat);
3455: MatDestroy(&dm->defaultConstraintMat);
3456: dm->defaultConstraintMat = mat;
3457: return(0);
3458: }
3460: #ifdef PETSC_USE_DEBUG
3463: /*
3464: DMDefaultSectionCheckConsistency - Check the consistentcy of the global and local sections.
3466: Input Parameters:
3467: + dm - The DM
3468: . localSection - PetscSection describing the local data layout
3469: - globalSection - PetscSection describing the global data layout
3471: Level: intermediate
3473: .seealso: DMGetDefaultSF(), DMSetDefaultSF()
3474: */
3475: static PetscErrorCode DMDefaultSectionCheckConsistency_Internal(DM dm, PetscSection localSection, PetscSection globalSection)
3476: {
3477: MPI_Comm comm;
3478: PetscLayout layout;
3479: const PetscInt *ranges;
3480: PetscInt pStart, pEnd, p, nroots;
3481: PetscMPIInt size, rank;
3482: PetscBool valid = PETSC_TRUE, gvalid;
3483: PetscErrorCode ierr;
3486: PetscObjectGetComm((PetscObject)dm,&comm);
3488: MPI_Comm_size(comm, &size);
3489: MPI_Comm_rank(comm, &rank);
3490: PetscSectionGetChart(globalSection, &pStart, &pEnd);
3491: PetscSectionGetConstrainedStorageSize(globalSection, &nroots);
3492: PetscLayoutCreate(comm, &layout);
3493: PetscLayoutSetBlockSize(layout, 1);
3494: PetscLayoutSetLocalSize(layout, nroots);
3495: PetscLayoutSetUp(layout);
3496: PetscLayoutGetRanges(layout, &ranges);
3497: for (p = pStart; p < pEnd; ++p) {
3498: PetscInt dof, cdof, off, gdof, gcdof, goff, gsize, d;
3500: PetscSectionGetDof(localSection, p, &dof);
3501: PetscSectionGetOffset(localSection, p, &off);
3502: PetscSectionGetConstraintDof(localSection, p, &cdof);
3503: PetscSectionGetDof(globalSection, p, &gdof);
3504: PetscSectionGetConstraintDof(globalSection, p, &gcdof);
3505: PetscSectionGetOffset(globalSection, p, &goff);
3506: if (!gdof) continue; /* Censored point */
3507: 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;}
3508: 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;}
3509: if (gdof < 0) {
3510: gsize = gdof < 0 ? -(gdof+1)-gcdof : gdof-gcdof;
3511: for (d = 0; d < gsize; ++d) {
3512: PetscInt offset = -(goff+1) + d, r;
3514: PetscFindInt(offset,size+1,ranges,&r);
3515: if (r < 0) r = -(r+2);
3516: 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;}
3517: }
3518: }
3519: }
3520: PetscLayoutDestroy(&layout);
3521: PetscSynchronizedFlush(comm, NULL);
3522: MPIU_Allreduce(&valid, &gvalid, 1, MPIU_BOOL, MPI_LAND, comm);
3523: if (!gvalid) {
3524: DMView(dm, NULL);
3525: SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Inconsistent local and global sections");
3526: }
3527: return(0);
3528: }
3529: #endif
3533: /*@
3534: DMGetDefaultGlobalSection - Get the PetscSection encoding the global data layout for the DM.
3536: Collective on DM
3538: Input Parameter:
3539: . dm - The DM
3541: Output Parameter:
3542: . section - The PetscSection
3544: Level: intermediate
3546: Note: This gets a borrowed reference, so the user should not destroy this PetscSection.
3548: .seealso: DMSetDefaultSection(), DMGetDefaultSection()
3549: @*/
3550: PetscErrorCode DMGetDefaultGlobalSection(DM dm, PetscSection *section)
3551: {
3557: if (!dm->defaultGlobalSection) {
3558: PetscSection s;
3560: DMGetDefaultSection(dm, &s);
3561: if (!s) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONGSTATE, "DM must have a default PetscSection in order to create a global PetscSection");
3562: if (!dm->sf) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DM must have a default PetscSF in order to create a global PetscSection");
3563: PetscSectionCreateGlobalSection(s, dm->sf, PETSC_FALSE, PETSC_FALSE, &dm->defaultGlobalSection);
3564: PetscLayoutDestroy(&dm->map);
3565: PetscSectionGetValueLayout(PetscObjectComm((PetscObject)dm), dm->defaultGlobalSection, &dm->map);
3566: PetscSectionViewFromOptions(dm->defaultGlobalSection, NULL, "-global_section_view");
3567: }
3568: *section = dm->defaultGlobalSection;
3569: return(0);
3570: }
3574: /*@
3575: DMSetDefaultGlobalSection - Set the PetscSection encoding the global data layout for the DM.
3577: Input Parameters:
3578: + dm - The DM
3579: - section - The PetscSection, or NULL
3581: Level: intermediate
3583: Note: Any existing Section will be destroyed
3585: .seealso: DMGetDefaultGlobalSection(), DMSetDefaultSection()
3586: @*/
3587: PetscErrorCode DMSetDefaultGlobalSection(DM dm, PetscSection section)
3588: {
3594: PetscObjectReference((PetscObject)section);
3595: PetscSectionDestroy(&dm->defaultGlobalSection);
3596: dm->defaultGlobalSection = section;
3597: #ifdef PETSC_USE_DEBUG
3598: if (section) {DMDefaultSectionCheckConsistency_Internal(dm, dm->defaultSection, section);}
3599: #endif
3600: return(0);
3601: }
3605: /*@
3606: DMGetDefaultSF - Get the PetscSF encoding the parallel dof overlap for the DM. If it has not been set,
3607: it is created from the default PetscSection layouts in the DM.
3609: Input Parameter:
3610: . dm - The DM
3612: Output Parameter:
3613: . sf - The PetscSF
3615: Level: intermediate
3617: Note: This gets a borrowed reference, so the user should not destroy this PetscSF.
3619: .seealso: DMSetDefaultSF(), DMCreateDefaultSF()
3620: @*/
3621: PetscErrorCode DMGetDefaultSF(DM dm, PetscSF *sf)
3622: {
3623: PetscInt nroots;
3629: PetscSFGetGraph(dm->defaultSF, &nroots, NULL, NULL, NULL);
3630: if (nroots < 0) {
3631: PetscSection section, gSection;
3633: DMGetDefaultSection(dm, §ion);
3634: if (section) {
3635: DMGetDefaultGlobalSection(dm, &gSection);
3636: DMCreateDefaultSF(dm, section, gSection);
3637: } else {
3638: *sf = NULL;
3639: return(0);
3640: }
3641: }
3642: *sf = dm->defaultSF;
3643: return(0);
3644: }
3648: /*@
3649: DMSetDefaultSF - Set the PetscSF encoding the parallel dof overlap for the DM
3651: Input Parameters:
3652: + dm - The DM
3653: - sf - The PetscSF
3655: Level: intermediate
3657: Note: Any previous SF is destroyed
3659: .seealso: DMGetDefaultSF(), DMCreateDefaultSF()
3660: @*/
3661: PetscErrorCode DMSetDefaultSF(DM dm, PetscSF sf)
3662: {
3668: PetscSFDestroy(&dm->defaultSF);
3669: dm->defaultSF = sf;
3670: return(0);
3671: }
3675: /*@C
3676: DMCreateDefaultSF - Create the PetscSF encoding the parallel dof overlap for the DM based upon the PetscSections
3677: describing the data layout.
3679: Input Parameters:
3680: + dm - The DM
3681: . localSection - PetscSection describing the local data layout
3682: - globalSection - PetscSection describing the global data layout
3684: Level: intermediate
3686: .seealso: DMGetDefaultSF(), DMSetDefaultSF()
3687: @*/
3688: PetscErrorCode DMCreateDefaultSF(DM dm, PetscSection localSection, PetscSection globalSection)
3689: {
3690: MPI_Comm comm;
3691: PetscLayout layout;
3692: const PetscInt *ranges;
3693: PetscInt *local;
3694: PetscSFNode *remote;
3695: PetscInt pStart, pEnd, p, nroots, nleaves = 0, l;
3696: PetscMPIInt size, rank;
3700: PetscObjectGetComm((PetscObject)dm,&comm);
3702: MPI_Comm_size(comm, &size);
3703: MPI_Comm_rank(comm, &rank);
3704: PetscSectionGetChart(globalSection, &pStart, &pEnd);
3705: PetscSectionGetConstrainedStorageSize(globalSection, &nroots);
3706: PetscLayoutCreate(comm, &layout);
3707: PetscLayoutSetBlockSize(layout, 1);
3708: PetscLayoutSetLocalSize(layout, nroots);
3709: PetscLayoutSetUp(layout);
3710: PetscLayoutGetRanges(layout, &ranges);
3711: for (p = pStart; p < pEnd; ++p) {
3712: PetscInt gdof, gcdof;
3714: PetscSectionGetDof(globalSection, p, &gdof);
3715: PetscSectionGetConstraintDof(globalSection, p, &gcdof);
3716: 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));
3717: nleaves += gdof < 0 ? -(gdof+1)-gcdof : gdof-gcdof;
3718: }
3719: PetscMalloc1(nleaves, &local);
3720: PetscMalloc1(nleaves, &remote);
3721: for (p = pStart, l = 0; p < pEnd; ++p) {
3722: const PetscInt *cind;
3723: PetscInt dof, cdof, off, gdof, gcdof, goff, gsize, d, c;
3725: PetscSectionGetDof(localSection, p, &dof);
3726: PetscSectionGetOffset(localSection, p, &off);
3727: PetscSectionGetConstraintDof(localSection, p, &cdof);
3728: PetscSectionGetConstraintIndices(localSection, p, &cind);
3729: PetscSectionGetDof(globalSection, p, &gdof);
3730: PetscSectionGetConstraintDof(globalSection, p, &gcdof);
3731: PetscSectionGetOffset(globalSection, p, &goff);
3732: if (!gdof) continue; /* Censored point */
3733: gsize = gdof < 0 ? -(gdof+1)-gcdof : gdof-gcdof;
3734: if (gsize != dof-cdof) {
3735: 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);
3736: cdof = 0; /* Ignore constraints */
3737: }
3738: for (d = 0, c = 0; d < dof; ++d) {
3739: if ((c < cdof) && (cind[c] == d)) {++c; continue;}
3740: local[l+d-c] = off+d;
3741: }
3742: if (gdof < 0) {
3743: for (d = 0; d < gsize; ++d, ++l) {
3744: PetscInt offset = -(goff+1) + d, r;
3746: PetscFindInt(offset,size+1,ranges,&r);
3747: if (r < 0) r = -(r+2);
3748: 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);
3749: remote[l].rank = r;
3750: remote[l].index = offset - ranges[r];
3751: }
3752: } else {
3753: for (d = 0; d < gsize; ++d, ++l) {
3754: remote[l].rank = rank;
3755: remote[l].index = goff+d - ranges[rank];
3756: }
3757: }
3758: }
3759: if (l != nleaves) SETERRQ2(comm, PETSC_ERR_PLIB, "Iteration error, l %d != nleaves %d", l, nleaves);
3760: PetscLayoutDestroy(&layout);
3761: PetscSFSetGraph(dm->defaultSF, nroots, nleaves, local, PETSC_OWN_POINTER, remote, PETSC_OWN_POINTER);
3762: return(0);
3763: }
3767: /*@
3768: DMGetPointSF - Get the PetscSF encoding the parallel section point overlap for the DM.
3770: Input Parameter:
3771: . dm - The DM
3773: Output Parameter:
3774: . sf - The PetscSF
3776: Level: intermediate
3778: Note: This gets a borrowed reference, so the user should not destroy this PetscSF.
3780: .seealso: DMSetPointSF(), DMGetDefaultSF(), DMSetDefaultSF(), DMCreateDefaultSF()
3781: @*/
3782: PetscErrorCode DMGetPointSF(DM dm, PetscSF *sf)
3783: {
3787: *sf = dm->sf;
3788: return(0);
3789: }
3793: /*@
3794: DMSetPointSF - Set the PetscSF encoding the parallel section point overlap for the DM.
3796: Input Parameters:
3797: + dm - The DM
3798: - sf - The PetscSF
3800: Level: intermediate
3802: .seealso: DMGetPointSF(), DMGetDefaultSF(), DMSetDefaultSF(), DMCreateDefaultSF()
3803: @*/
3804: PetscErrorCode DMSetPointSF(DM dm, PetscSF sf)
3805: {
3811: PetscSFDestroy(&dm->sf);
3812: PetscObjectReference((PetscObject) sf);
3813: dm->sf = sf;
3814: return(0);
3815: }
3819: /*@
3820: DMGetDS - Get the PetscDS
3822: Input Parameter:
3823: . dm - The DM
3825: Output Parameter:
3826: . prob - The PetscDS
3828: Level: developer
3830: .seealso: DMSetDS()
3831: @*/
3832: PetscErrorCode DMGetDS(DM dm, PetscDS *prob)
3833: {
3837: *prob = dm->prob;
3838: return(0);
3839: }
3843: /*@
3844: DMSetDS - Set the PetscDS
3846: Input Parameters:
3847: + dm - The DM
3848: - prob - The PetscDS
3850: Level: developer
3852: .seealso: DMGetDS()
3853: @*/
3854: PetscErrorCode DMSetDS(DM dm, PetscDS prob)
3855: {
3861: PetscDSDestroy(&dm->prob);
3862: dm->prob = prob;
3863: PetscObjectReference((PetscObject) dm->prob);
3864: return(0);
3865: }
3869: PetscErrorCode DMGetNumFields(DM dm, PetscInt *numFields)
3870: {
3875: PetscDSGetNumFields(dm->prob, numFields);
3876: return(0);
3877: }
3881: PetscErrorCode DMSetNumFields(DM dm, PetscInt numFields)
3882: {
3883: PetscInt Nf, f;
3888: PetscDSGetNumFields(dm->prob, &Nf);
3889: for (f = Nf; f < numFields; ++f) {
3890: PetscContainer obj;
3892: PetscContainerCreate(PetscObjectComm((PetscObject) dm), &obj);
3893: PetscDSSetDiscretization(dm->prob, f, (PetscObject) obj);
3894: PetscContainerDestroy(&obj);
3895: }
3896: return(0);
3897: }
3901: /*@
3902: DMGetField - Return the discretization object for a given DM field
3904: Not collective
3906: Input Parameters:
3907: + dm - The DM
3908: - f - The field number
3910: Output Parameter:
3911: . field - The discretization object
3913: Level: developer
3915: .seealso: DMSetField()
3916: @*/
3917: PetscErrorCode DMGetField(DM dm, PetscInt f, PetscObject *field)
3918: {
3923: PetscDSGetDiscretization(dm->prob, f, field);
3924: return(0);
3925: }
3929: /*@
3930: DMSetField - Set the discretization object for a given DM field
3932: Logically collective on DM
3934: Input Parameters:
3935: + dm - The DM
3936: . f - The field number
3937: - field - The discretization object
3939: Level: developer
3941: .seealso: DMGetField()
3942: @*/
3943: PetscErrorCode DMSetField(DM dm, PetscInt f, PetscObject field)
3944: {
3949: PetscDSSetDiscretization(dm->prob, f, field);
3950: return(0);
3951: }
3955: PetscErrorCode DMRestrictHook_Coordinates(DM dm,DM dmc,void *ctx)
3956: {
3957: DM dm_coord,dmc_coord;
3959: Vec coords,ccoords;
3960: Mat inject;
3962: DMGetCoordinateDM(dm,&dm_coord);
3963: DMGetCoordinateDM(dmc,&dmc_coord);
3964: DMGetCoordinates(dm,&coords);
3965: DMGetCoordinates(dmc,&ccoords);
3966: if (coords && !ccoords) {
3967: DMCreateGlobalVector(dmc_coord,&ccoords);
3968: PetscObjectSetName((PetscObject)ccoords,"coordinates");
3969: DMCreateInjection(dmc_coord,dm_coord,&inject);
3970: MatRestrict(inject,coords,ccoords);
3971: MatDestroy(&inject);
3972: DMSetCoordinates(dmc,ccoords);
3973: VecDestroy(&ccoords);
3974: }
3975: return(0);
3976: }
3980: static PetscErrorCode DMSubDomainHook_Coordinates(DM dm,DM subdm,void *ctx)
3981: {
3982: DM dm_coord,subdm_coord;
3984: Vec coords,ccoords,clcoords;
3985: VecScatter *scat_i,*scat_g;
3987: DMGetCoordinateDM(dm,&dm_coord);
3988: DMGetCoordinateDM(subdm,&subdm_coord);
3989: DMGetCoordinates(dm,&coords);
3990: DMGetCoordinates(subdm,&ccoords);
3991: if (coords && !ccoords) {
3992: DMCreateGlobalVector(subdm_coord,&ccoords);
3993: PetscObjectSetName((PetscObject)ccoords,"coordinates");
3994: DMCreateLocalVector(subdm_coord,&clcoords);
3995: PetscObjectSetName((PetscObject)clcoords,"coordinates");
3996: DMCreateDomainDecompositionScatters(dm_coord,1,&subdm_coord,NULL,&scat_i,&scat_g);
3997: VecScatterBegin(scat_i[0],coords,ccoords,INSERT_VALUES,SCATTER_FORWARD);
3998: VecScatterBegin(scat_g[0],coords,clcoords,INSERT_VALUES,SCATTER_FORWARD);
3999: VecScatterEnd(scat_i[0],coords,ccoords,INSERT_VALUES,SCATTER_FORWARD);
4000: VecScatterEnd(scat_g[0],coords,clcoords,INSERT_VALUES,SCATTER_FORWARD);
4001: DMSetCoordinates(subdm,ccoords);
4002: DMSetCoordinatesLocal(subdm,clcoords);
4003: VecScatterDestroy(&scat_i[0]);
4004: VecScatterDestroy(&scat_g[0]);
4005: VecDestroy(&ccoords);
4006: VecDestroy(&clcoords);
4007: PetscFree(scat_i);
4008: PetscFree(scat_g);
4009: }
4010: return(0);
4011: }
4015: /*@
4016: DMGetDimension - Return the topological dimension of the DM
4018: Not collective
4020: Input Parameter:
4021: . dm - The DM
4023: Output Parameter:
4024: . dim - The topological dimension
4026: Level: beginner
4028: .seealso: DMSetDimension(), DMCreate()
4029: @*/
4030: PetscErrorCode DMGetDimension(DM dm, PetscInt *dim)
4031: {
4035: *dim = dm->dim;
4036: return(0);
4037: }
4041: /*@
4042: DMSetDimension - Set the topological dimension of the DM
4044: Collective on dm
4046: Input Parameters:
4047: + dm - The DM
4048: - dim - The topological dimension
4050: Level: beginner
4052: .seealso: DMGetDimension(), DMCreate()
4053: @*/
4054: PetscErrorCode DMSetDimension(DM dm, PetscInt dim)
4055: {
4059: dm->dim = dim;
4060: return(0);
4061: }
4065: /*@
4066: DMGetDimPoints - Get the half-open interval for all points of a given dimension
4068: Collective on DM
4070: Input Parameters:
4071: + dm - the DM
4072: - dim - the dimension
4074: Output Parameters:
4075: + pStart - The first point of the given dimension
4076: . pEnd - The first point following points of the given dimension
4078: Note:
4079: The points are vertices in the Hasse diagram encoding the topology. This is explained in
4080: http://arxiv.org/abs/0908.4427. If not points exist of this dimension in the storage scheme,
4081: then the interval is empty.
4083: Level: intermediate
4085: .keywords: point, Hasse Diagram, dimension
4086: .seealso: DMPLEX, DMPlexGetDepthStratum(), DMPlexGetHeightStratum()
4087: @*/
4088: PetscErrorCode DMGetDimPoints(DM dm, PetscInt dim, PetscInt *pStart, PetscInt *pEnd)
4089: {
4090: PetscInt d;
4095: DMGetDimension(dm, &d);
4096: if ((dim < 0) || (dim > d)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension %d 1", dim, d);
4097: (*dm->ops->getdimpoints)(dm, dim, pStart, pEnd);
4098: return(0);
4099: }
4103: /*@
4104: DMSetCoordinates - Sets into the DM a global vector that holds the coordinates
4106: Collective on DM
4108: Input Parameters:
4109: + dm - the DM
4110: - c - coordinate vector
4112: Note:
4113: The coordinates do include those for ghost points, which are in the local vector
4115: Level: intermediate
4117: .keywords: distributed array, get, corners, nodes, local indices, coordinates
4118: .seealso: DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLoca(), DMGetCoordinateDM()
4119: @*/
4120: PetscErrorCode DMSetCoordinates(DM dm, Vec c)
4121: {
4127: PetscObjectReference((PetscObject) c);
4128: VecDestroy(&dm->coordinates);
4129: dm->coordinates = c;
4130: VecDestroy(&dm->coordinatesLocal);
4131: DMCoarsenHookAdd(dm,DMRestrictHook_Coordinates,NULL,NULL);
4132: DMSubDomainHookAdd(dm,DMSubDomainHook_Coordinates,NULL,NULL);
4133: return(0);
4134: }
4138: /*@
4139: DMSetCoordinatesLocal - Sets into the DM a local vector that holds the coordinates
4141: Collective on DM
4143: Input Parameters:
4144: + dm - the DM
4145: - c - coordinate vector
4147: Note:
4148: The coordinates of ghost points can be set using DMSetCoordinates()
4149: followed by DMGetCoordinatesLocal(). This is intended to enable the
4150: setting of ghost coordinates outside of the domain.
4152: Level: intermediate
4154: .keywords: distributed array, get, corners, nodes, local indices, coordinates
4155: .seealso: DMGetCoordinatesLocal(), DMSetCoordinates(), DMGetCoordinates(), DMGetCoordinateDM()
4156: @*/
4157: PetscErrorCode DMSetCoordinatesLocal(DM dm, Vec c)
4158: {
4164: PetscObjectReference((PetscObject) c);
4165: VecDestroy(&dm->coordinatesLocal);
4167: dm->coordinatesLocal = c;
4169: VecDestroy(&dm->coordinates);
4170: return(0);
4171: }
4175: /*@
4176: DMGetCoordinates - Gets a global vector with the coordinates associated with the DM.
4178: Not Collective
4180: Input Parameter:
4181: . dm - the DM
4183: Output Parameter:
4184: . c - global coordinate vector
4186: Note:
4187: This is a borrowed reference, so the user should NOT destroy this vector
4189: Each process has only the local coordinates (does NOT have the ghost coordinates).
4191: For DMDA, in two and three dimensions coordinates are interlaced (x_0,y_0,x_1,y_1,...)
4192: and (x_0,y_0,z_0,x_1,y_1,z_1...)
4194: Level: intermediate
4196: .keywords: distributed array, get, corners, nodes, local indices, coordinates
4197: .seealso: DMSetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM()
4198: @*/
4199: PetscErrorCode DMGetCoordinates(DM dm, Vec *c)
4200: {
4206: if (!dm->coordinates && dm->coordinatesLocal) {
4207: DM cdm = NULL;
4209: DMGetCoordinateDM(dm, &cdm);
4210: DMCreateGlobalVector(cdm, &dm->coordinates);
4211: PetscObjectSetName((PetscObject) dm->coordinates, "coordinates");
4212: DMLocalToGlobalBegin(cdm, dm->coordinatesLocal, INSERT_VALUES, dm->coordinates);
4213: DMLocalToGlobalEnd(cdm, dm->coordinatesLocal, INSERT_VALUES, dm->coordinates);
4214: }
4215: *c = dm->coordinates;
4216: return(0);
4217: }
4221: /*@
4222: DMGetCoordinatesLocal - Gets a local vector with the coordinates associated with the DM.
4224: Collective on DM
4226: Input Parameter:
4227: . dm - the DM
4229: Output Parameter:
4230: . c - coordinate vector
4232: Note:
4233: This is a borrowed reference, so the user should NOT destroy this vector
4235: Each process has the local and ghost coordinates
4237: For DMDA, in two and three dimensions coordinates are interlaced (x_0,y_0,x_1,y_1,...)
4238: and (x_0,y_0,z_0,x_1,y_1,z_1...)
4240: Level: intermediate
4242: .keywords: distributed array, get, corners, nodes, local indices, coordinates
4243: .seealso: DMSetCoordinatesLocal(), DMGetCoordinates(), DMSetCoordinates(), DMGetCoordinateDM()
4244: @*/
4245: PetscErrorCode DMGetCoordinatesLocal(DM dm, Vec *c)
4246: {
4252: if (!dm->coordinatesLocal && dm->coordinates) {
4253: DM cdm = NULL;
4255: DMGetCoordinateDM(dm, &cdm);
4256: DMCreateLocalVector(cdm, &dm->coordinatesLocal);
4257: PetscObjectSetName((PetscObject) dm->coordinatesLocal, "coordinates");
4258: DMGlobalToLocalBegin(cdm, dm->coordinates, INSERT_VALUES, dm->coordinatesLocal);
4259: DMGlobalToLocalEnd(cdm, dm->coordinates, INSERT_VALUES, dm->coordinatesLocal);
4260: }
4261: *c = dm->coordinatesLocal;
4262: return(0);
4263: }
4267: /*@
4268: DMGetCoordinateDM - Gets the DM that prescribes coordinate layout and scatters between global and local coordinates
4270: Collective on DM
4272: Input Parameter:
4273: . dm - the DM
4275: Output Parameter:
4276: . cdm - coordinate DM
4278: Level: intermediate
4280: .keywords: distributed array, get, corners, nodes, local indices, coordinates
4281: .seealso: DMSetCoordinateDM(), DMSetCoordinates(), DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLocal()
4282: @*/
4283: PetscErrorCode DMGetCoordinateDM(DM dm, DM *cdm)
4284: {
4290: if (!dm->coordinateDM) {
4291: if (!dm->ops->createcoordinatedm) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unable to create coordinates for this DM");
4292: (*dm->ops->createcoordinatedm)(dm, &dm->coordinateDM);
4293: }
4294: *cdm = dm->coordinateDM;
4295: return(0);
4296: }
4300: /*@
4301: DMSetCoordinateDM - Sets the DM that prescribes coordinate layout and scatters between global and local coordinates
4303: Logically Collective on DM
4305: Input Parameters:
4306: + dm - the DM
4307: - cdm - coordinate DM
4309: Level: intermediate
4311: .keywords: distributed array, get, corners, nodes, local indices, coordinates
4312: .seealso: DMGetCoordinateDM(), DMSetCoordinates(), DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLocal()
4313: @*/
4314: PetscErrorCode DMSetCoordinateDM(DM dm, DM cdm)
4315: {
4321: DMDestroy(&dm->coordinateDM);
4322: dm->coordinateDM = cdm;
4323: PetscObjectReference((PetscObject) dm->coordinateDM);
4324: return(0);
4325: }
4329: /*@
4330: DMGetCoordinateDim - Retrieve the dimension of embedding space for coordinate values.
4332: Not Collective
4334: Input Parameter:
4335: . dm - The DM object
4337: Output Parameter:
4338: . dim - The embedding dimension
4340: Level: intermediate
4342: .keywords: mesh, coordinates
4343: .seealso: DMSetCoordinateDim(), DMGetCoordinateSection(), DMGetCoordinateDM(), DMGetDefaultSection(), DMSetDefaultSection()
4344: @*/
4345: PetscErrorCode DMGetCoordinateDim(DM dm, PetscInt *dim)
4346: {
4350: if (dm->dimEmbed == PETSC_DEFAULT) {
4351: dm->dimEmbed = dm->dim;
4352: }
4353: *dim = dm->dimEmbed;
4354: return(0);
4355: }
4359: /*@
4360: DMSetCoordinateDim - Set the dimension of the embedding space for coordinate values.
4362: Not Collective
4364: Input Parameters:
4365: + dm - The DM object
4366: - dim - The embedding dimension
4368: Level: intermediate
4370: .keywords: mesh, coordinates
4371: .seealso: DMGetCoordinateDim(), DMSetCoordinateSection(), DMGetCoordinateSection(), DMGetDefaultSection(), DMSetDefaultSection()
4372: @*/
4373: PetscErrorCode DMSetCoordinateDim(DM dm, PetscInt dim)
4374: {
4377: dm->dimEmbed = dim;
4378: return(0);
4379: }
4383: /*@
4384: DMGetCoordinateSection - Retrieve the layout of coordinate values over the mesh.
4386: Not Collective
4388: Input Parameter:
4389: . dm - The DM object
4391: Output Parameter:
4392: . section - The PetscSection object
4394: Level: intermediate
4396: .keywords: mesh, coordinates
4397: .seealso: DMGetCoordinateDM(), DMGetDefaultSection(), DMSetDefaultSection()
4398: @*/
4399: PetscErrorCode DMGetCoordinateSection(DM dm, PetscSection *section)
4400: {
4401: DM cdm;
4407: DMGetCoordinateDM(dm, &cdm);
4408: DMGetDefaultSection(cdm, section);
4409: return(0);
4410: }
4414: /*@
4415: DMSetCoordinateSection - Set the layout of coordinate values over the mesh.
4417: Not Collective
4419: Input Parameters:
4420: + dm - The DM object
4421: . dim - The embedding dimension, or PETSC_DETERMINE
4422: - section - The PetscSection object
4424: Level: intermediate
4426: .keywords: mesh, coordinates
4427: .seealso: DMGetCoordinateSection(), DMGetDefaultSection(), DMSetDefaultSection()
4428: @*/
4429: PetscErrorCode DMSetCoordinateSection(DM dm, PetscInt dim, PetscSection section)
4430: {
4431: DM cdm;
4437: DMGetCoordinateDM(dm, &cdm);
4438: DMSetDefaultSection(cdm, section);
4439: if (dim == PETSC_DETERMINE) {
4440: PetscInt d = dim;
4441: PetscInt pStart, pEnd, vStart, vEnd, v, dd;
4443: PetscSectionGetChart(section, &pStart, &pEnd);
4444: DMGetDimPoints(dm, 0, &vStart, &vEnd);
4445: pStart = PetscMax(vStart, pStart);
4446: pEnd = PetscMin(vEnd, pEnd);
4447: for (v = pStart; v < pEnd; ++v) {
4448: PetscSectionGetDof(section, v, &dd);
4449: if (dd) {d = dd; break;}
4450: }
4451: DMSetCoordinateDim(dm, d);
4452: }
4453: return(0);
4454: }
4458: /*@C
4459: DMSetPeriodicity - Set the description of mesh periodicity
4461: Input Parameters:
4462: + dm - The DM object
4463: . maxCell - Over distances greater than this, we can assume a point has crossed over to another sheet, when trying to localize cell coordinates
4464: . L - If we assume the mesh is a torus, this is the length of each coordinate
4465: - bd - This describes the type of periodicity in each topological dimension
4467: Level: developer
4469: .seealso: DMGetPeriodicity()
4470: @*/
4471: PetscErrorCode DMGetPeriodicity(DM dm, const PetscReal **maxCell, const PetscReal **L, const DMBoundaryType **bd)
4472: {
4475: if (L) *L = dm->L;
4476: if (maxCell) *maxCell = dm->maxCell;
4477: if (bd) *bd = dm->bdtype;
4478: return(0);
4479: }
4483: /*@C
4484: DMSetPeriodicity - Set the description of mesh periodicity
4486: Input Parameters:
4487: + dm - The DM object
4488: . maxCell - Over distances greater than this, we can assume a point has crossed over to another sheet, when trying to localize cell coordinates
4489: . L - If we assume the mesh is a torus, this is the length of each coordinate
4490: - bd - This describes the type of periodicity in each topological dimension
4492: Level: developer
4494: .seealso: DMGetPeriodicity()
4495: @*/
4496: PetscErrorCode DMSetPeriodicity(DM dm, const PetscReal maxCell[], const PetscReal L[], const DMBoundaryType bd[])
4497: {
4498: PetscInt dim, d;
4504: PetscFree3(dm->L,dm->maxCell,dm->bdtype);
4505: DMGetDimension(dm, &dim);
4506: PetscMalloc3(dim,&dm->L,dim,&dm->maxCell,dim,&dm->bdtype);
4507: for (d = 0; d < dim; ++d) {dm->L[d] = L[d]; dm->maxCell[d] = maxCell[d]; dm->bdtype[d] = bd[d];}
4508: return(0);
4509: }
4513: /*@
4514: 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.
4516: Input Parameters:
4517: + dm - The DM
4518: - in - The input coordinate point (dim numbers)
4520: Output Parameter:
4521: . out - The localized coordinate point
4523: Level: developer
4525: .seealso: DMLocalizeCoordinates(), DMLocalizeAddCoordinate()
4526: @*/
4527: PetscErrorCode DMLocalizeCoordinate(DM dm, const PetscScalar in[], PetscScalar out[])
4528: {
4529: PetscInt dim, d;
4533: DMGetCoordinateDim(dm, &dim);
4534: if (!dm->maxCell) {
4535: for (d = 0; d < dim; ++d) out[d] = in[d];
4536: } else {
4537: for (d = 0; d < dim; ++d) {
4538: out[d] = in[d] - dm->L[d]*floor(PetscRealPart(in[d])/dm->L[d]);
4539: }
4540: }
4541: return(0);
4542: }
4546: /*
4547: 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.
4549: Input Parameters:
4550: + dm - The DM
4551: . dim - The spatial dimension
4552: . anchor - The anchor point, the input point can be no more than maxCell away from it
4553: - in - The input coordinate point (dim numbers)
4555: Output Parameter:
4556: . out - The localized coordinate point
4558: Level: developer
4560: 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
4562: .seealso: DMLocalizeCoordinates(), DMLocalizeAddCoordinate()
4563: */
4564: PetscErrorCode DMLocalizeCoordinate_Internal(DM dm, PetscInt dim, const PetscScalar anchor[], const PetscScalar in[], PetscScalar out[])
4565: {
4566: PetscInt d;
4569: if (!dm->maxCell) {
4570: for (d = 0; d < dim; ++d) out[d] = in[d];
4571: } else {
4572: for (d = 0; d < dim; ++d) {
4573: if (PetscAbsScalar(anchor[d] - in[d]) > dm->maxCell[d]) {
4574: out[d] = PetscRealPart(anchor[d]) > PetscRealPart(in[d]) ? dm->L[d] + in[d] : in[d] - dm->L[d];
4575: } else {
4576: out[d] = in[d];
4577: }
4578: }
4579: }
4580: return(0);
4581: }
4584: PetscErrorCode DMLocalizeCoordinateReal_Internal(DM dm, PetscInt dim, const PetscReal anchor[], const PetscReal in[], PetscReal out[])
4585: {
4586: PetscInt d;
4589: if (!dm->maxCell) {
4590: for (d = 0; d < dim; ++d) out[d] = in[d];
4591: } else {
4592: for (d = 0; d < dim; ++d) {
4593: if (PetscAbsReal(anchor[d] - in[d]) > dm->maxCell[d]) {
4594: out[d] = anchor[d] > in[d] ? dm->L[d] + in[d] : in[d] - dm->L[d];
4595: } else {
4596: out[d] = in[d];
4597: }
4598: }
4599: }
4600: return(0);
4601: }
4605: /*
4606: 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.
4608: Input Parameters:
4609: + dm - The DM
4610: . dim - The spatial dimension
4611: . anchor - The anchor point, the input point can be no more than maxCell away from it
4612: . in - The input coordinate delta (dim numbers)
4613: - out - The input coordinate point (dim numbers)
4615: Output Parameter:
4616: . out - The localized coordinate in + out
4618: Level: developer
4620: 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
4622: .seealso: DMLocalizeCoordinates(), DMLocalizeCoordinate()
4623: */
4624: PetscErrorCode DMLocalizeAddCoordinate_Internal(DM dm, PetscInt dim, const PetscScalar anchor[], const PetscScalar in[], PetscScalar out[])
4625: {
4626: PetscInt d;
4629: if (!dm->maxCell) {
4630: for (d = 0; d < dim; ++d) out[d] += in[d];
4631: } else {
4632: for (d = 0; d < dim; ++d) {
4633: if (PetscAbsScalar(anchor[d] - in[d]) > dm->maxCell[d]) {
4634: out[d] += PetscRealPart(anchor[d]) > PetscRealPart(in[d]) ? dm->L[d] + in[d] : in[d] - dm->L[d];
4635: } else {
4636: out[d] += in[d];
4637: }
4638: }
4639: }
4640: return(0);
4641: }
4643: PETSC_EXTERN PetscErrorCode DMPlexGetDepthStratum(DM, PetscInt, PetscInt *, PetscInt *);
4644: PETSC_EXTERN PetscErrorCode DMPlexGetHeightStratum(DM, PetscInt, PetscInt *, PetscInt *);
4645: PETSC_EXTERN PetscErrorCode DMPlexVecGetClosure(DM, PetscSection, Vec, PetscInt, PetscInt *, PetscScalar *[]);
4646: PETSC_EXTERN PetscErrorCode DMPlexVecRestoreClosure(DM, PetscSection, Vec, PetscInt, PetscInt *, PetscScalar *[]);
4650: /*@
4651: DMGetCoordinatesLocalized - Check if the DM coordinates have been localized for cells
4653: Input Parameter:
4654: . dm - The DM
4656: Output Parameter:
4657: areLocalized - True if localized
4659: Level: developer
4661: .seealso: DMLocalizeCoordinates()
4662: @*/
4663: PetscErrorCode DMGetCoordinatesLocalized(DM dm,PetscBool *areLocalized)
4664: {
4665: DM cdm;
4666: PetscSection coordSection;
4667: PetscInt cStart, cEnd, c, sStart, sEnd, dof;
4668: PetscBool alreadyLocalized, alreadyLocalizedGlobal;
4673: if (!dm->maxCell) {
4674: *areLocalized = PETSC_FALSE;
4675: return(0);
4676: }
4677: /* We need some generic way of refering to cells/vertices */
4678: DMGetCoordinateDM(dm, &cdm);
4679: {
4680: PetscBool isplex;
4682: PetscObjectTypeCompare((PetscObject) cdm, DMPLEX, &isplex);
4683: if (isplex) {
4684: DMPlexGetHeightStratum(cdm, 0, &cStart, &cEnd);
4685: } else SETERRQ(PetscObjectComm((PetscObject) cdm), PETSC_ERR_ARG_WRONG, "Coordinate localization requires a DMPLEX coordinate DM");
4686: }
4687: DMGetCoordinateSection(dm, &coordSection);
4688: PetscSectionGetChart(coordSection,&sStart,&sEnd);
4689: alreadyLocalized = alreadyLocalizedGlobal = PETSC_TRUE;
4690: for (c = cStart; c < cEnd; ++c) {
4691: if (c < sStart || c >= sEnd) {
4692: alreadyLocalized = PETSC_FALSE;
4693: break;
4694: }
4695: PetscSectionGetDof(coordSection, c, &dof);
4696: if (!dof) {
4697: alreadyLocalized = PETSC_FALSE;
4698: break;
4699: }
4700: }
4701: MPI_Allreduce(&alreadyLocalized,&alreadyLocalizedGlobal,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)dm));
4702: *areLocalized = alreadyLocalizedGlobal;
4703: return(0);
4704: }
4709: /*@
4710: DMLocalizeCoordinates - If a mesh is periodic, create local coordinates for each cell
4712: Input Parameter:
4713: . dm - The DM
4715: Level: developer
4717: .seealso: DMLocalizeCoordinate(), DMLocalizeAddCoordinate()
4718: @*/
4719: PetscErrorCode DMLocalizeCoordinates(DM dm)
4720: {
4721: DM cdm;
4722: PetscSection coordSection, cSection;
4723: Vec coordinates, cVec;
4724: PetscScalar *coords, *coords2, *anchor;
4725: PetscInt Nc, cStart, cEnd, c, vStart, vEnd, v, sStart, sEnd, dof, d, off, off2, bs, coordSize;
4726: PetscBool alreadyLocalized, alreadyLocalizedGlobal;
4731: if (!dm->maxCell) return(0);
4732: /* We need some generic way of refering to cells/vertices */
4733: DMGetCoordinateDM(dm, &cdm);
4734: {
4735: PetscBool isplex;
4737: PetscObjectTypeCompare((PetscObject) cdm, DMPLEX, &isplex);
4738: if (isplex) {
4739: DMPlexGetHeightStratum(cdm, 0, &cStart, &cEnd);
4740: DMPlexGetDepthStratum(cdm, 0, &vStart, &vEnd);
4741: } else SETERRQ(PetscObjectComm((PetscObject) cdm), PETSC_ERR_ARG_WRONG, "Coordinate localization requires a DMPLEX coordinate DM");
4742: }
4743: DMGetCoordinatesLocal(dm, &coordinates);
4744: DMGetCoordinateSection(dm, &coordSection);
4745: PetscSectionGetChart(coordSection,&sStart,&sEnd);
4746: alreadyLocalized = alreadyLocalizedGlobal = PETSC_TRUE;
4747: for (c = cStart; c < cEnd; ++c) {
4748: if (c < sStart || c >= sEnd) {
4749: alreadyLocalized = PETSC_FALSE;
4750: break;
4751: }
4752: PetscSectionGetDof(coordSection, c, &dof);
4753: if (!dof) {
4754: alreadyLocalized = PETSC_FALSE;
4755: break;
4756: }
4757: }
4758: MPI_Allreduce(&alreadyLocalized,&alreadyLocalizedGlobal,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)dm));
4759: if (alreadyLocalizedGlobal) return(0);
4760: PetscSectionCreate(PetscObjectComm((PetscObject) dm), &cSection);
4761: PetscSectionSetNumFields(cSection, 1);
4762: PetscSectionGetFieldComponents(coordSection, 0, &Nc);
4763: PetscSectionSetFieldComponents(cSection, 0, Nc);
4764: PetscSectionSetChart(cSection, cStart, vEnd);
4765: for (v = vStart; v < vEnd; ++v) {
4766: PetscSectionGetDof(coordSection, v, &dof);
4767: PetscSectionSetDof(cSection, v, dof);
4768: PetscSectionSetFieldDof(cSection, v, 0, dof);
4769: }
4770: for (c = cStart; c < cEnd; ++c) {
4771: DMPlexVecGetClosure(cdm, coordSection, coordinates, c, &dof, NULL);
4772: PetscSectionSetDof(cSection, c, dof);
4773: PetscSectionSetFieldDof(cSection, c, 0, dof);
4774: }
4775: PetscSectionSetUp(cSection);
4776: PetscSectionGetStorageSize(cSection, &coordSize);
4777: VecCreate(PETSC_COMM_SELF, &cVec);
4778: PetscObjectSetName((PetscObject)cVec,"coordinates");
4779: VecGetBlockSize(coordinates, &bs);
4780: VecSetBlockSize(cVec, bs);
4781: VecSetSizes(cVec, coordSize, PETSC_DETERMINE);
4782: VecSetType(cVec, VECSTANDARD);
4783: VecGetArrayRead(coordinates, (const PetscScalar**)&coords);
4784: VecGetArray(cVec, &coords2);
4785: for (v = vStart; v < vEnd; ++v) {
4786: PetscSectionGetDof(coordSection, v, &dof);
4787: PetscSectionGetOffset(coordSection, v, &off);
4788: PetscSectionGetOffset(cSection, v, &off2);
4789: for (d = 0; d < dof; ++d) coords2[off2+d] = coords[off+d];
4790: }
4791: DMGetWorkArray(dm, 3, PETSC_SCALAR, &anchor);
4792: for (c = cStart; c < cEnd; ++c) {
4793: PetscScalar *cellCoords = NULL;
4794: PetscInt b;
4796: DMPlexVecGetClosure(cdm, coordSection, coordinates, c, &dof, &cellCoords);
4797: PetscSectionGetOffset(cSection, c, &off2);
4798: for (b = 0; b < bs; ++b) anchor[b] = cellCoords[b];
4799: for (d = 0; d < dof/bs; ++d) {DMLocalizeCoordinate_Internal(dm, bs, anchor, &cellCoords[d*bs], &coords2[off2+d*bs]);}
4800: DMPlexVecRestoreClosure(cdm, coordSection, coordinates, c, &dof, &cellCoords);
4801: }
4802: DMRestoreWorkArray(dm, 3, PETSC_SCALAR, &anchor);
4803: VecRestoreArrayRead(coordinates, (const PetscScalar**)&coords);
4804: VecRestoreArray(cVec, &coords2);
4805: DMSetCoordinateSection(dm, PETSC_DETERMINE, cSection);
4806: DMSetCoordinatesLocal(dm, cVec);
4807: VecDestroy(&cVec);
4808: PetscSectionDestroy(&cSection);
4809: return(0);
4810: }
4814: /*@
4815: DMLocatePoints - Locate the points in v in the mesh and return a PetscSF of the containing cells
4817: Collective on Vec v (see explanation below)
4819: Input Parameters:
4820: + dm - The DM
4821: . v - The Vec of points
4822: - cells - Points to either NULL, or a PetscSF with guesses for which cells contain each point.
4824: Output Parameter:
4825: . cells - The PetscSF containing the ranks and local indices of the containing points.
4828: Level: developer
4830: To do a search of the local cells of the mesh, v should have PETSC_COMM_SELF as its communicator.
4832: To do a search of all the cells in the distributed mesh, v should have the same communicator as
4833: dm.
4835: If *cellSF is NULL on input, a PetscSF will be created.
4837: If *cellSF is not NULL on input, it should point to an existing PetscSF, whose graph will be used as initial
4838: guesses.
4840: An array that maps each point to its containing cell can be obtained with
4842: const PetscSFNode *cells;
4843: PetscInt nFound;
4844: const PetscSFNode *found;
4846: PetscSFGetGraph(cells,NULL,&nFound,&found,&cells);
4848: Where cells[i].rank is the rank of the cell containing point found[i] (or i if found == NULL), and cells[i].index is
4849: the index of the cell in its rank's local numbering.
4851: .keywords: point location, mesh
4852: .seealso: DMSetCoordinates(), DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLocal()
4853: @*/
4854: PetscErrorCode DMLocatePoints(DM dm, Vec v, PetscSF *cellSF)
4855: {
4862: if (*cellSF) {
4863: PetscMPIInt result;
4866: MPI_Comm_compare(PetscObjectComm((PetscObject)v),PetscObjectComm((PetscObject)cellSF),&result);
4867: if (result != MPI_IDENT && result != MPI_CONGRUENT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"cellSF must have a communicator congruent to v's");
4868: }
4869: else {
4870: PetscSFCreate(PetscObjectComm((PetscObject)v),cellSF);
4871: }
4872: PetscLogEventBegin(DM_LocatePoints,dm,0,0,0);
4873: if (dm->ops->locatepoints) {
4874: (*dm->ops->locatepoints)(dm,v,*cellSF);
4875: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Point location not available for this DM");
4876: PetscLogEventEnd(DM_LocatePoints,dm,0,0,0);
4877: return(0);
4878: }
4882: /*@
4883: DMGetOutputDM - Retrieve the DM associated with the layout for output
4885: Input Parameter:
4886: . dm - The original DM
4888: Output Parameter:
4889: . odm - The DM which provides the layout for output
4891: Level: intermediate
4893: .seealso: VecView(), DMGetDefaultGlobalSection()
4894: @*/
4895: PetscErrorCode DMGetOutputDM(DM dm, DM *odm)
4896: {
4897: PetscSection section;
4898: PetscBool hasConstraints;
4904: DMGetDefaultSection(dm, §ion);
4905: PetscSectionHasConstraints(section, &hasConstraints);
4906: if (!hasConstraints) {
4907: *odm = dm;
4908: return(0);
4909: }
4910: if (!dm->dmBC) {
4911: PetscDS ds;
4912: PetscSection newSection, gsection;
4913: PetscSF sf;
4915: DMClone(dm, &dm->dmBC);
4916: DMGetDS(dm, &ds);
4917: DMSetDS(dm->dmBC, ds);
4918: PetscSectionClone(section, &newSection);
4919: DMSetDefaultSection(dm->dmBC, newSection);
4920: PetscSectionDestroy(&newSection);
4921: DMGetPointSF(dm->dmBC, &sf);
4922: PetscSectionCreateGlobalSection(section, sf, PETSC_TRUE, PETSC_FALSE, &gsection);
4923: DMSetDefaultGlobalSection(dm->dmBC, gsection);
4924: PetscSectionDestroy(&gsection);
4925: }
4926: *odm = dm->dmBC;
4927: return(0);
4928: }
4932: /*@
4933: DMGetOutputSequenceNumber - Retrieve the sequence number/value for output
4935: Input Parameter:
4936: . dm - The original DM
4938: Output Parameters:
4939: + num - The output sequence number
4940: - val - The output sequence value
4942: Level: intermediate
4944: Note: This is intended for output that should appear in sequence, for instance
4945: a set of timesteps in an HDF5 file, or a set of realizations of a stochastic system.
4947: .seealso: VecView()
4948: @*/
4949: PetscErrorCode DMGetOutputSequenceNumber(DM dm, PetscInt *num, PetscReal *val)
4950: {
4955: return(0);
4956: }
4960: /*@
4961: DMSetOutputSequenceNumber - Set the sequence number/value for output
4963: Input Parameters:
4964: + dm - The original DM
4965: . num - The output sequence number
4966: - val - The output sequence value
4968: Level: intermediate
4970: Note: This is intended for output that should appear in sequence, for instance
4971: a set of timesteps in an HDF5 file, or a set of realizations of a stochastic system.
4973: .seealso: VecView()
4974: @*/
4975: PetscErrorCode DMSetOutputSequenceNumber(DM dm, PetscInt num, PetscReal val)
4976: {
4979: dm->outputSequenceNum = num;
4980: dm->outputSequenceVal = val;
4981: return(0);
4982: }
4986: /*@C
4987: DMOutputSequenceLoad - Retrieve the sequence value from a Viewer
4989: Input Parameters:
4990: + dm - The original DM
4991: . name - The sequence name
4992: - num - The output sequence number
4994: Output Parameter:
4995: . val - The output sequence value
4997: Level: intermediate
4999: Note: This is intended for output that should appear in sequence, for instance
5000: a set of timesteps in an HDF5 file, or a set of realizations of a stochastic system.
5002: .seealso: DMGetOutputSequenceNumber(), DMSetOutputSequenceNumber(), VecView()
5003: @*/
5004: PetscErrorCode DMOutputSequenceLoad(DM dm, PetscViewer viewer, const char *name, PetscInt num, PetscReal *val)
5005: {
5006: PetscBool ishdf5;
5013: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERHDF5, &ishdf5);
5014: if (ishdf5) {
5015: #if defined(PETSC_HAVE_HDF5)
5016: PetscScalar value;
5018: DMSequenceLoad_HDF5(dm, name, num, &value, viewer);
5019: *val = PetscRealPart(value);
5020: #endif
5021: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid viewer; open viewer with PetscViewerHDF5Open()");
5022: return(0);
5023: }
5027: /*@
5028: DMGetUseNatural - Get the flag for creating a mapping to the natural order on distribution
5030: Not collective
5032: Input Parameter:
5033: . dm - The DM
5035: Output Parameter:
5036: . useNatural - The flag to build the mapping to a natural order during distribution
5038: Level: beginner
5040: .seealso: DMSetUseNatural(), DMCreate()
5041: @*/
5042: PetscErrorCode DMGetUseNatural(DM dm, PetscBool *useNatural)
5043: {
5047: *useNatural = dm->useNatural;
5048: return(0);
5049: }
5053: /*@
5054: DMSetUseNatural - Set the flag for creating a mapping to the natural order on distribution
5056: Collective on dm
5058: Input Parameters:
5059: + dm - The DM
5060: - useNatural - The flag to build the mapping to a natural order during distribution
5062: Level: beginner
5064: .seealso: DMGetUseNatural(), DMCreate()
5065: @*/
5066: PetscErrorCode DMSetUseNatural(DM dm, PetscBool useNatural)
5067: {
5071: dm->useNatural = useNatural;
5072: return(0);
5073: }
5080: /*@C
5081: DMCreateLabel - Create a label of the given name if it does not already exist
5083: Not Collective
5085: Input Parameters:
5086: + dm - The DM object
5087: - name - The label name
5089: Level: intermediate
5091: .keywords: mesh
5092: .seealso: DMLabelCreate(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5093: @*/
5094: PetscErrorCode DMCreateLabel(DM dm, const char name[])
5095: {
5096: DMLabelLink next = dm->labels->next;
5097: PetscBool flg = PETSC_FALSE;
5103: while (next) {
5104: PetscStrcmp(name, next->label->name, &flg);
5105: if (flg) break;
5106: next = next->next;
5107: }
5108: if (!flg) {
5109: DMLabelLink tmpLabel;
5111: PetscCalloc1(1, &tmpLabel);
5112: DMLabelCreate(name, &tmpLabel->label);
5113: tmpLabel->output = PETSC_TRUE;
5114: tmpLabel->next = dm->labels->next;
5115: dm->labels->next = tmpLabel;
5116: }
5117: return(0);
5118: }
5122: /*@C
5123: DMGetLabelValue - Get the value in a Sieve Label for the given point, with 0 as the default
5125: Not Collective
5127: Input Parameters:
5128: + dm - The DM object
5129: . name - The label name
5130: - point - The mesh point
5132: Output Parameter:
5133: . value - The label value for this point, or -1 if the point is not in the label
5135: Level: beginner
5137: .keywords: mesh
5138: .seealso: DMLabelGetValue(), DMSetLabelValue(), DMGetStratumIS()
5139: @*/
5140: PetscErrorCode DMGetLabelValue(DM dm, const char name[], PetscInt point, PetscInt *value)
5141: {
5142: DMLabel label;
5148: DMGetLabel(dm, name, &label);
5149: if (!label) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No label named %s was found", name);
5150: DMLabelGetValue(label, point, value);
5151: return(0);
5152: }
5156: /*@C
5157: DMSetLabelValue - Add a point to a Sieve Label with given value
5159: Not Collective
5161: Input Parameters:
5162: + dm - The DM object
5163: . name - The label name
5164: . point - The mesh point
5165: - value - The label value for this point
5167: Output Parameter:
5169: Level: beginner
5171: .keywords: mesh
5172: .seealso: DMLabelSetValue(), DMGetStratumIS(), DMClearLabelValue()
5173: @*/
5174: PetscErrorCode DMSetLabelValue(DM dm, const char name[], PetscInt point, PetscInt value)
5175: {
5176: DMLabel label;
5182: DMGetLabel(dm, name, &label);
5183: if (!label) {
5184: DMCreateLabel(dm, name);
5185: DMGetLabel(dm, name, &label);
5186: }
5187: DMLabelSetValue(label, point, value);
5188: return(0);
5189: }
5193: /*@C
5194: DMClearLabelValue - Remove a point from a Sieve Label with given value
5196: Not Collective
5198: Input Parameters:
5199: + dm - The DM object
5200: . name - The label name
5201: . point - The mesh point
5202: - value - The label value for this point
5204: Output Parameter:
5206: Level: beginner
5208: .keywords: mesh
5209: .seealso: DMLabelClearValue(), DMSetLabelValue(), DMGetStratumIS()
5210: @*/
5211: PetscErrorCode DMClearLabelValue(DM dm, const char name[], PetscInt point, PetscInt value)
5212: {
5213: DMLabel label;
5219: DMGetLabel(dm, name, &label);
5220: if (!label) return(0);
5221: DMLabelClearValue(label, point, value);
5222: return(0);
5223: }
5227: /*@C
5228: DMGetLabelSize - Get the number of different integer ids in a Label
5230: Not Collective
5232: Input Parameters:
5233: + dm - The DM object
5234: - name - The label name
5236: Output Parameter:
5237: . size - The number of different integer ids, or 0 if the label does not exist
5239: Level: beginner
5241: .keywords: mesh
5242: .seealso: DMLabeGetNumValues(), DMSetLabelValue()
5243: @*/
5244: PetscErrorCode DMGetLabelSize(DM dm, const char name[], PetscInt *size)
5245: {
5246: DMLabel label;
5253: DMGetLabel(dm, name, &label);
5254: *size = 0;
5255: if (!label) return(0);
5256: DMLabelGetNumValues(label, size);
5257: return(0);
5258: }
5262: /*@C
5263: DMGetLabelIdIS - Get the integer ids in a label
5265: Not Collective
5267: Input Parameters:
5268: + mesh - The DM object
5269: - name - The label name
5271: Output Parameter:
5272: . ids - The integer ids, or NULL if the label does not exist
5274: Level: beginner
5276: .keywords: mesh
5277: .seealso: DMLabelGetValueIS(), DMGetLabelSize()
5278: @*/
5279: PetscErrorCode DMGetLabelIdIS(DM dm, const char name[], IS *ids)
5280: {
5281: DMLabel label;
5288: DMGetLabel(dm, name, &label);
5289: *ids = NULL;
5290: if (!label) return(0);
5291: DMLabelGetValueIS(label, ids);
5292: return(0);
5293: }
5297: /*@C
5298: DMGetStratumSize - Get the number of points in a label stratum
5300: Not Collective
5302: Input Parameters:
5303: + dm - The DM object
5304: . name - The label name
5305: - value - The stratum value
5307: Output Parameter:
5308: . size - The stratum size
5310: Level: beginner
5312: .keywords: mesh
5313: .seealso: DMLabelGetStratumSize(), DMGetLabelSize(), DMGetLabelIds()
5314: @*/
5315: PetscErrorCode DMGetStratumSize(DM dm, const char name[], PetscInt value, PetscInt *size)
5316: {
5317: DMLabel label;
5324: DMGetLabel(dm, name, &label);
5325: *size = 0;
5326: if (!label) return(0);
5327: DMLabelGetStratumSize(label, value, size);
5328: return(0);
5329: }
5333: /*@C
5334: DMGetStratumIS - Get the points in a label stratum
5336: Not Collective
5338: Input Parameters:
5339: + dm - The DM object
5340: . name - The label name
5341: - value - The stratum value
5343: Output Parameter:
5344: . points - The stratum points, or NULL if the label does not exist or does not have that value
5346: Level: beginner
5348: .keywords: mesh
5349: .seealso: DMLabelGetStratumIS(), DMGetStratumSize()
5350: @*/
5351: PetscErrorCode DMGetStratumIS(DM dm, const char name[], PetscInt value, IS *points)
5352: {
5353: DMLabel label;
5360: DMGetLabel(dm, name, &label);
5361: *points = NULL;
5362: if (!label) return(0);
5363: DMLabelGetStratumIS(label, value, points);
5364: return(0);
5365: }
5369: /*@C
5370: DMClearLabelStratum - Remove all points from a stratum from a Sieve Label
5372: Not Collective
5374: Input Parameters:
5375: + dm - The DM object
5376: . name - The label name
5377: - value - The label value for this point
5379: Output Parameter:
5381: Level: beginner
5383: .keywords: mesh
5384: .seealso: DMLabelClearStratum(), DMSetLabelValue(), DMGetStratumIS(), DMClearLabelValue()
5385: @*/
5386: PetscErrorCode DMClearLabelStratum(DM dm, const char name[], PetscInt value)
5387: {
5388: DMLabel label;
5394: DMGetLabel(dm, name, &label);
5395: if (!label) return(0);
5396: DMLabelClearStratum(label, value);
5397: return(0);
5398: }
5402: /*@
5403: DMGetNumLabels - Return the number of labels defined by the mesh
5405: Not Collective
5407: Input Parameter:
5408: . dm - The DM object
5410: Output Parameter:
5411: . numLabels - the number of Labels
5413: Level: intermediate
5415: .keywords: mesh
5416: .seealso: DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5417: @*/
5418: PetscErrorCode DMGetNumLabels(DM dm, PetscInt *numLabels)
5419: {
5420: DMLabelLink next = dm->labels->next;
5421: PetscInt n = 0;
5426: while (next) {++n; next = next->next;}
5427: *numLabels = n;
5428: return(0);
5429: }
5433: /*@C
5434: DMGetLabelName - Return the name of nth label
5436: Not Collective
5438: Input Parameters:
5439: + dm - The DM object
5440: - n - the label number
5442: Output Parameter:
5443: . name - the label name
5445: Level: intermediate
5447: .keywords: mesh
5448: .seealso: DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5449: @*/
5450: PetscErrorCode DMGetLabelName(DM dm, PetscInt n, const char **name)
5451: {
5452: DMLabelLink next = dm->labels->next;
5453: PetscInt l = 0;
5458: while (next) {
5459: if (l == n) {
5460: *name = next->label->name;
5461: return(0);
5462: }
5463: ++l;
5464: next = next->next;
5465: }
5466: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label %D does not exist in this DM", n);
5467: }
5471: /*@C
5472: DMHasLabel - Determine whether the mesh has a label of a given name
5474: Not Collective
5476: Input Parameters:
5477: + dm - The DM object
5478: - name - The label name
5480: Output Parameter:
5481: . hasLabel - PETSC_TRUE if the label is present
5483: Level: intermediate
5485: .keywords: mesh
5486: .seealso: DMCreateLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5487: @*/
5488: PetscErrorCode DMHasLabel(DM dm, const char name[], PetscBool *hasLabel)
5489: {
5490: DMLabelLink next = dm->labels->next;
5497: *hasLabel = PETSC_FALSE;
5498: while (next) {
5499: PetscStrcmp(name, next->label->name, hasLabel);
5500: if (*hasLabel) break;
5501: next = next->next;
5502: }
5503: return(0);
5504: }
5508: /*@C
5509: DMGetLabel - Return the label of a given name, or NULL
5511: Not Collective
5513: Input Parameters:
5514: + dm - The DM object
5515: - name - The label name
5517: Output Parameter:
5518: . label - The DMLabel, or NULL if the label is absent
5520: Level: intermediate
5522: .keywords: mesh
5523: .seealso: DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5524: @*/
5525: PetscErrorCode DMGetLabel(DM dm, const char name[], DMLabel *label)
5526: {
5527: DMLabelLink next = dm->labels->next;
5528: PetscBool hasLabel;
5535: *label = NULL;
5536: while (next) {
5537: PetscStrcmp(name, next->label->name, &hasLabel);
5538: if (hasLabel) {
5539: *label = next->label;
5540: break;
5541: }
5542: next = next->next;
5543: }
5544: return(0);
5545: }
5549: /*@C
5550: DMGetLabelByNum - Return the nth label
5552: Not Collective
5554: Input Parameters:
5555: + dm - The DM object
5556: - n - the label number
5558: Output Parameter:
5559: . label - the label
5561: Level: intermediate
5563: .keywords: mesh
5564: .seealso: DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5565: @*/
5566: PetscErrorCode DMGetLabelByNum(DM dm, PetscInt n, DMLabel *label)
5567: {
5568: DMLabelLink next = dm->labels->next;
5569: PetscInt l = 0;
5574: while (next) {
5575: if (l == n) {
5576: *label = next->label;
5577: return(0);
5578: }
5579: ++l;
5580: next = next->next;
5581: }
5582: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label %D does not exist in this DM", n);
5583: }
5587: /*@C
5588: DMAddLabel - Add the label to this mesh
5590: Not Collective
5592: Input Parameters:
5593: + dm - The DM object
5594: - label - The DMLabel
5596: Level: developer
5598: .keywords: mesh
5599: .seealso: DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5600: @*/
5601: PetscErrorCode DMAddLabel(DM dm, DMLabel label)
5602: {
5603: DMLabelLink tmpLabel;
5604: PetscBool hasLabel;
5609: DMHasLabel(dm, label->name, &hasLabel);
5610: if (hasLabel) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label %s already exists in this DM", label->name);
5611: PetscCalloc1(1, &tmpLabel);
5612: tmpLabel->label = label;
5613: tmpLabel->output = PETSC_TRUE;
5614: tmpLabel->next = dm->labels->next;
5615: dm->labels->next = tmpLabel;
5616: return(0);
5617: }
5621: /*@C
5622: DMRemoveLabel - Remove the label from this mesh
5624: Not Collective
5626: Input Parameters:
5627: + dm - The DM object
5628: - name - The label name
5630: Output Parameter:
5631: . label - The DMLabel, or NULL if the label is absent
5633: Level: developer
5635: .keywords: mesh
5636: .seealso: DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5637: @*/
5638: PetscErrorCode DMRemoveLabel(DM dm, const char name[], DMLabel *label)
5639: {
5640: DMLabelLink next = dm->labels->next;
5641: DMLabelLink last = NULL;
5642: PetscBool hasLabel;
5647: DMHasLabel(dm, name, &hasLabel);
5648: *label = NULL;
5649: if (!hasLabel) return(0);
5650: while (next) {
5651: PetscStrcmp(name, next->label->name, &hasLabel);
5652: if (hasLabel) {
5653: if (last) last->next = next->next;
5654: else dm->labels->next = next->next;
5655: next->next = NULL;
5656: *label = next->label;
5657: PetscStrcmp(name, "depth", &hasLabel);
5658: if (hasLabel) {
5659: dm->depthLabel = NULL;
5660: }
5661: PetscFree(next);
5662: break;
5663: }
5664: last = next;
5665: next = next->next;
5666: }
5667: return(0);
5668: }
5672: /*@C
5673: DMGetLabelOutput - Get the output flag for a given label
5675: Not Collective
5677: Input Parameters:
5678: + dm - The DM object
5679: - name - The label name
5681: Output Parameter:
5682: . output - The flag for output
5684: Level: developer
5686: .keywords: mesh
5687: .seealso: DMSetLabelOutput(), DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5688: @*/
5689: PetscErrorCode DMGetLabelOutput(DM dm, const char name[], PetscBool *output)
5690: {
5691: DMLabelLink next = dm->labels->next;
5698: while (next) {
5699: PetscBool flg;
5701: PetscStrcmp(name, next->label->name, &flg);
5702: if (flg) {*output = next->output; return(0);}
5703: next = next->next;
5704: }
5705: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No label named %s was present in this dm", name);
5706: }
5710: /*@C
5711: DMSetLabelOutput - Set the output flag for a given label
5713: Not Collective
5715: Input Parameters:
5716: + dm - The DM object
5717: . name - The label name
5718: - output - The flag for output
5720: Level: developer
5722: .keywords: mesh
5723: .seealso: DMGetLabelOutput(), DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5724: @*/
5725: PetscErrorCode DMSetLabelOutput(DM dm, const char name[], PetscBool output)
5726: {
5727: DMLabelLink next = dm->labels->next;
5733: while (next) {
5734: PetscBool flg;
5736: PetscStrcmp(name, next->label->name, &flg);
5737: if (flg) {next->output = output; return(0);}
5738: next = next->next;
5739: }
5740: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No label named %s was present in this dm", name);
5741: }
5746: /*@
5747: DMCopyLabels - Copy labels from one mesh to another with a superset of the points
5749: Collective on DM
5751: Input Parameter:
5752: . dmA - The DM object with initial labels
5754: Output Parameter:
5755: . dmB - The DM object with copied labels
5757: Level: intermediate
5759: Note: This is typically used when interpolating or otherwise adding to a mesh
5761: .keywords: mesh
5762: .seealso: DMCopyCoordinates(), DMGetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM(), DMGetCoordinateSection()
5763: @*/
5764: PetscErrorCode DMCopyLabels(DM dmA, DM dmB)
5765: {
5766: PetscInt numLabels, l;
5770: if (dmA == dmB) return(0);
5771: DMGetNumLabels(dmA, &numLabels);
5772: for (l = 0; l < numLabels; ++l) {
5773: DMLabel label, labelNew;
5774: const char *name;
5775: PetscBool flg;
5777: DMGetLabelName(dmA, l, &name);
5778: PetscStrcmp(name, "depth", &flg);
5779: if (flg) continue;
5780: DMGetLabel(dmA, name, &label);
5781: DMLabelDuplicate(label, &labelNew);
5782: DMAddLabel(dmB, labelNew);
5783: }
5784: return(0);
5785: }
5789: /*@
5790: DMGetCoarseDM - Get the coarse mesh from which this was obtained by refinement
5792: Input Parameter:
5793: . dm - The DM object
5795: Output Parameter:
5796: . cdm - The coarse DM
5798: Level: intermediate
5800: .seealso: DMSetCoarseDM()
5801: @*/
5802: PetscErrorCode DMGetCoarseDM(DM dm, DM *cdm)
5803: {
5807: *cdm = dm->coarseMesh;
5808: return(0);
5809: }
5813: /*@
5814: DMSetCoarseDM - Set the coarse mesh from which this was obtained by refinement
5816: Input Parameters:
5817: + dm - The DM object
5818: - cdm - The coarse DM
5820: Level: intermediate
5822: .seealso: DMGetCoarseDM()
5823: @*/
5824: PetscErrorCode DMSetCoarseDM(DM dm, DM cdm)
5825: {
5831: PetscObjectReference((PetscObject)cdm);
5832: DMDestroy(&dm->coarseMesh);
5833: dm->coarseMesh = cdm;
5834: return(0);
5835: }
5839: /*@
5840: DMGetFineDM - Get the fine mesh from which this was obtained by refinement
5842: Input Parameter:
5843: . dm - The DM object
5845: Output Parameter:
5846: . fdm - The fine DM
5848: Level: intermediate
5850: .seealso: DMSetFineDM()
5851: @*/
5852: PetscErrorCode DMGetFineDM(DM dm, DM *fdm)
5853: {
5857: *fdm = dm->fineMesh;
5858: return(0);
5859: }
5863: /*@
5864: DMSetFineDM - Set the fine mesh from which this was obtained by refinement
5866: Input Parameters:
5867: + dm - The DM object
5868: - fdm - The fine DM
5870: Level: intermediate
5872: .seealso: DMGetFineDM()
5873: @*/
5874: PetscErrorCode DMSetFineDM(DM dm, DM fdm)
5875: {
5881: PetscObjectReference((PetscObject)fdm);
5882: DMDestroy(&dm->fineMesh);
5883: dm->fineMesh = fdm;
5884: return(0);
5885: }
5887: /*=== DMBoundary code ===*/
5891: PetscErrorCode DMBoundaryDuplicate(DMBoundaryLinkList bd, DMBoundaryLinkList *boundary)
5892: {
5893: DMBoundary b = bd->next, b2, bold = NULL;
5897: PetscNew(boundary);
5898: (*boundary)->refct = 1;
5899: (*boundary)->next = NULL;
5900: for (; b; b = b->next, bold = b2) {
5901: PetscNew(&b2);
5902: PetscStrallocpy(b->name, (char **) &b2->name);
5903: PetscStrallocpy(b->labelname, (char **) &b2->labelname);
5904: PetscMalloc1(b->numids, &b2->ids);
5905: PetscMemcpy(b2->ids, b->ids, b->numids*sizeof(PetscInt));
5906: PetscMalloc1(b->numcomps, &b2->comps);
5907: PetscMemcpy(b2->comps, b->comps, b->numcomps*sizeof(PetscInt));
5908: b2->label = NULL;
5909: b2->essential = b->essential;
5910: b2->field = b->field;
5911: b2->numcomps = b->numcomps;
5912: b2->func = b->func;
5913: b2->numids = b->numids;
5914: b2->ctx = b->ctx;
5915: b2->next = NULL;
5916: if (!(*boundary)->next) (*boundary)->next = b2;
5917: if (bold) bold->next = b2;
5918: }
5919: return(0);
5920: }
5924: PetscErrorCode DMBoundaryDestroy(DMBoundaryLinkList *boundary)
5925: {
5926: DMBoundary b, next;
5930: if (!boundary) return(0);
5931: if (--((*boundary)->refct)) {
5932: *boundary = NULL;
5933: return(0);
5934: }
5935: b = (*boundary)->next;
5936: for (; b; b = next) {
5937: next = b->next;
5938: PetscFree(b->comps);
5939: PetscFree(b->ids);
5940: PetscFree(b->name);
5941: PetscFree(b->labelname);
5942: PetscFree(b);
5943: }
5944: PetscFree(*boundary);
5945: return(0);
5946: }
5950: PetscErrorCode DMCopyBoundary(DM dm, DM dmNew)
5951: {
5952: DMBoundary b;
5956: DMBoundaryDestroy(&dmNew->boundary);
5957: DMBoundaryDuplicate(dm->boundary, &dmNew->boundary);
5958: for (b = dmNew->boundary->next; b; b = b->next) {
5959: if (b->labelname) {
5960: DMGetLabel(dmNew, b->labelname, &b->label);
5961: if (!b->label) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Label %s does not exist in this DM", b->labelname);
5962: }
5963: }
5964: return(0);
5965: }
5969: /*@C
5970: DMAddBoundary - Add a boundary condition to the model
5972: Input Parameters:
5973: + dm - The mesh object
5974: . isEssential - Flag for an essential (Dirichlet) condition, as opposed to a natural (Neumann) condition
5975: . name - The BC name
5976: . labelname - The label defining constrained points
5977: . field - The field to constrain
5978: . numcomps - The number of constrained field components
5979: . comps - An array of constrained component numbers
5980: . bcFunc - A pointwise function giving boundary values
5981: . numids - The number of DMLabel ids for constrained points
5982: . ids - An array of ids for constrained points
5983: - ctx - An optional user context for bcFunc
5985: Options Database Keys:
5986: + -bc_<boundary name> <num> - Overrides the boundary ids
5987: - -bc_<boundary name>_comp <num> - Overrides the boundary components
5989: Level: developer
5991: .seealso: DMGetBoundary()
5992: @*/
5993: PetscErrorCode DMAddBoundary(DM dm, PetscBool isEssential, const char name[], const char labelname[], PetscInt field, PetscInt numcomps, const PetscInt *comps, void (*bcFunc)(), PetscInt numids, const PetscInt *ids, void *ctx)
5994: {
5995: DMBoundary b;
6000: PetscNew(&b);
6001: PetscStrallocpy(name, (char **) &b->name);
6002: PetscStrallocpy(labelname, (char **) &b->labelname);
6003: PetscMalloc1(numcomps, &b->comps);
6004: if (numcomps) {PetscMemcpy(b->comps, comps, numcomps*sizeof(PetscInt));}
6005: PetscMalloc1(numids, &b->ids);
6006: if (numids) {PetscMemcpy(b->ids, ids, numids*sizeof(PetscInt));}
6007: if (b->labelname) {
6008: DMGetLabel(dm, b->labelname, &b->label);
6009: if (!b->label) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Label %s does not exist in this DM", b->labelname);
6010: }
6011: b->essential = isEssential;
6012: b->field = field;
6013: b->numcomps = numcomps;
6014: b->func = bcFunc;
6015: b->numids = numids;
6016: b->ctx = ctx;
6017: b->next = dm->boundary->next;
6018: dm->boundary->next = b;
6019: return(0);
6020: }
6024: /*@
6025: DMGetNumBoundary - Get the number of registered BC
6027: Input Parameters:
6028: . dm - The mesh object
6030: Output Parameters:
6031: . numBd - The number of BC
6033: Level: intermediate
6035: .seealso: DMAddBoundary(), DMGetBoundary()
6036: @*/
6037: PetscErrorCode DMGetNumBoundary(DM dm, PetscInt *numBd)
6038: {
6039: DMBoundary b = dm->boundary->next;
6044: *numBd = 0;
6045: while (b) {++(*numBd); b = b->next;}
6046: return(0);
6047: }
6051: /*@C
6052: DMGetBoundary - Add a boundary condition to the model
6054: Input Parameters:
6055: + dm - The mesh object
6056: - bd - The BC number
6058: Output Parameters:
6059: + isEssential - Flag for an essential (Dirichlet) condition, as opposed to a natural (Neumann) condition
6060: . name - The BC name
6061: . labelname - The label defining constrained points
6062: . field - The field to constrain
6063: . numcomps - The number of constrained field components
6064: . comps - An array of constrained component numbers
6065: . bcFunc - A pointwise function giving boundary values
6066: . numids - The number of DMLabel ids for constrained points
6067: . ids - An array of ids for constrained points
6068: - ctx - An optional user context for bcFunc
6070: Options Database Keys:
6071: + -bc_<boundary name> <num> - Overrides the boundary ids
6072: - -bc_<boundary name>_comp <num> - Overrides the boundary components
6074: Level: developer
6076: .seealso: DMAddBoundary()
6077: @*/
6078: PetscErrorCode DMGetBoundary(DM dm, PetscInt bd, PetscBool *isEssential, const char **name, const char **labelname, PetscInt *field, PetscInt *numcomps, const PetscInt **comps, void (**func)(), PetscInt *numids, const PetscInt **ids, void **ctx)
6079: {
6080: DMBoundary b = dm->boundary->next;
6081: PetscInt n = 0;
6085: while (b) {
6086: if (n == bd) break;
6087: b = b->next;
6088: ++n;
6089: }
6090: if (!b) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Boundary %d is not in [0, %d)", bd, n);
6091: if (isEssential) {
6093: *isEssential = b->essential;
6094: }
6095: if (name) {
6097: *name = b->name;
6098: }
6099: if (labelname) {
6101: *labelname = b->labelname;
6102: }
6103: if (field) {
6105: *field = b->field;
6106: }
6107: if (numcomps) {
6109: *numcomps = b->numcomps;
6110: }
6111: if (comps) {
6113: *comps = b->comps;
6114: }
6115: if (func) {
6117: *func = b->func;
6118: }
6119: if (numids) {
6121: *numids = b->numids;
6122: }
6123: if (ids) {
6125: *ids = b->ids;
6126: }
6127: if (ctx) {
6129: *ctx = b->ctx;
6130: }
6131: return(0);
6132: }
6136: PetscErrorCode DMIsBoundaryPoint(DM dm, PetscInt point, PetscBool *isBd)
6137: {
6138: DMBoundary b = dm->boundary->next;
6144: *isBd = PETSC_FALSE;
6145: while (b && !(*isBd)) {
6146: if (b->label) {
6147: PetscInt i;
6149: for (i = 0; i < b->numids && !(*isBd); ++i) {
6150: DMLabelStratumHasPoint(b->label, b->ids[i], point, isBd);
6151: }
6152: }
6153: b = b->next;
6154: }
6155: return(0);
6156: }
6160: /*@C
6161: DMProjectFunction - This projects the given function into the function space provided.
6163: Input Parameters:
6164: + dm - The DM
6165: . time - The time
6166: . funcs - The coordinate functions to evaluate, one per field
6167: . ctxs - Optional array of contexts to pass to each coordinate function. ctxs itself may be null.
6168: - mode - The insertion mode for values
6170: Output Parameter:
6171: . X - vector
6173: Calling sequence of func:
6174: $ func(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar u[], void *ctx);
6176: + dim - The spatial dimension
6177: . x - The coordinates
6178: . Nf - The number of fields
6179: . u - The output field values
6180: - ctx - optional user-defined function context
6182: Level: developer
6184: .seealso: DMComputeL2Diff()
6185: @*/
6186: PetscErrorCode DMProjectFunction(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec X)
6187: {
6188: Vec localX;
6193: DMGetLocalVector(dm, &localX);
6194: DMProjectFunctionLocal(dm, time, funcs, ctxs, mode, localX);
6195: DMLocalToGlobalBegin(dm, localX, mode, X);
6196: DMLocalToGlobalEnd(dm, localX, mode, X);
6197: DMRestoreLocalVector(dm, &localX);
6198: return(0);
6199: }
6203: PetscErrorCode DMProjectFunctionLocal(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX)
6204: {
6210: if (!dm->ops->projectfunctionlocal) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implemnt DMProjectFunctionLocal",((PetscObject)dm)->type_name);
6211: (dm->ops->projectfunctionlocal) (dm, time, funcs, ctxs, mode, localX);
6212: return(0);
6213: }
6217: PetscErrorCode DMProjectFieldLocal(DM dm, Vec localU,
6218: void (**funcs)(PetscInt, PetscInt, PetscInt,
6219: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
6220: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
6221: PetscReal, const PetscReal[], PetscScalar[]),
6222: InsertMode mode, Vec localX)
6223: {
6230: if (!dm->ops->projectfieldlocal) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implemnt DMProjectFieldLocal",((PetscObject)dm)->type_name);
6231: (dm->ops->projectfieldlocal) (dm, localU, funcs, mode, localX);
6232: return(0);
6233: }
6237: PetscErrorCode DMProjectFunctionLabelLocal(DM dm, PetscReal time, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX)
6238: {
6244: if (!dm->ops->projectfunctionlabellocal) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implemnt DMProjectFunctionLabelLocal",((PetscObject)dm)->type_name);
6245: (dm->ops->projectfunctionlabellocal) (dm, time, label, numIds, ids, funcs, ctxs, mode, localX);
6246: return(0);
6247: }
6251: /*@C
6252: DMComputeL2Diff - This function computes the L_2 difference between a function u and an FEM interpolant solution u_h.
6254: Input Parameters:
6255: + dm - The DM
6256: . time - The time
6257: . funcs - The functions to evaluate for each field component
6258: . ctxs - Optional array of contexts to pass to each function, or NULL.
6259: - X - The coefficient vector u_h
6261: Output Parameter:
6262: . diff - The diff ||u - u_h||_2
6264: Level: developer
6266: .seealso: DMProjectFunction(), DMComputeL2FieldDiff(), DMComputeL2GradientDiff()
6267: @*/
6268: PetscErrorCode DMComputeL2Diff(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff)
6269: {
6275: if (!dm->ops->computel2diff) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implemnt DMComputeL2Diff",((PetscObject)dm)->type_name);
6276: (dm->ops->computel2diff)(dm,time,funcs,ctxs,X,diff);
6277: return(0);
6278: }
6282: /*@C
6283: DMComputeL2GradientDiff - This function computes the L_2 difference between the gradient of a function u and an FEM interpolant solution grad u_h.
6285: Input Parameters:
6286: + dm - The DM
6287: , time - The time
6288: . funcs - The gradient functions to evaluate for each field component
6289: . ctxs - Optional array of contexts to pass to each function, or NULL.
6290: . X - The coefficient vector u_h
6291: - n - The vector to project along
6293: Output Parameter:
6294: . diff - The diff ||(grad u - grad u_h) . n||_2
6296: Level: developer
6298: .seealso: DMProjectFunction(), DMComputeL2Diff()
6299: @*/
6300: 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)
6301: {
6307: if (!dm->ops->computel2gradientdiff) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMComputeL2GradientDiff",((PetscObject)dm)->type_name);
6308: (dm->ops->computel2gradientdiff)(dm,time,funcs,ctxs,X,n,diff);
6309: return(0);
6310: }
6314: /*@C
6315: DMComputeL2FieldDiff - This function computes the L_2 difference between a function u and an FEM interpolant solution u_h, separated into field components.
6317: Input Parameters:
6318: + dm - The DM
6319: . time - The time
6320: . funcs - The functions to evaluate for each field component
6321: . ctxs - Optional array of contexts to pass to each function, or NULL.
6322: - X - The coefficient vector u_h
6324: Output Parameter:
6325: . diff - The array of differences, ||u^f - u^f_h||_2
6327: Level: developer
6329: .seealso: DMProjectFunction(), DMComputeL2FieldDiff(), DMComputeL2GradientDiff()
6330: @*/
6331: PetscErrorCode DMComputeL2FieldDiff(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal diff[])
6332: {
6338: if (!dm->ops->computel2fielddiff) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implemnt DMComputeL2FieldDiff",((PetscObject)dm)->type_name);
6339: (dm->ops->computel2fielddiff)(dm,time,funcs,ctxs,X,diff);
6340: return(0);
6341: }