Actual source code: dm.c
petsc-3.8.4 2018-03-24
1: #include <petsc/private/dmimpl.h>
2: #include <petsc/private/dmlabelimpl.h>
3: #include <petsc/private/petscdsimpl.h>
4: #include <petscdmplex.h>
5: #include <petscsf.h>
6: #include <petscds.h>
8: PetscClassId DM_CLASSID;
9: PetscLogEvent DM_Convert, DM_GlobalToLocal, DM_LocalToGlobal, DM_LocalToLocal, DM_LocatePoints, DM_Coarsen, DM_Refine, DM_CreateInterpolation, DM_CreateRestriction;
11: const char *const DMBoundaryTypes[] = {"NONE","GHOSTED","MIRROR","PERIODIC","TWIST","DM_BOUNDARY_",0};
13: /*@
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: *dm = v;
77: return(0);
78: }
80: /*@
81: DMClone - Creates a DM object with the same topology as the original.
83: Collective on MPI_Comm
85: Input Parameter:
86: . dm - The original DM object
88: Output Parameter:
89: . newdm - The new DM object
91: Level: beginner
93: .keywords: DM, topology, create
94: @*/
95: PetscErrorCode DMClone(DM dm, DM *newdm)
96: {
97: PetscSF sf;
98: Vec coords;
99: void *ctx;
100: PetscInt dim, cdim;
106: DMCreate(PetscObjectComm((PetscObject) dm), newdm);
107: PetscFree((*newdm)->labels);
108: dm->labels->refct++;
109: (*newdm)->labels = dm->labels;
110: (*newdm)->depthLabel = dm->depthLabel;
111: DMGetDimension(dm, &dim);
112: DMSetDimension(*newdm, dim);
113: if (dm->ops->clone) {
114: (*dm->ops->clone)(dm, newdm);
115: }
116: (*newdm)->setupcalled = dm->setupcalled;
117: DMGetPointSF(dm, &sf);
118: DMSetPointSF(*newdm, sf);
119: DMGetApplicationContext(dm, &ctx);
120: DMSetApplicationContext(*newdm, ctx);
121: if (dm->coordinateDM) {
122: DM ncdm;
123: PetscSection cs;
124: PetscInt pEnd = -1, pEndMax = -1;
126: DMGetDefaultSection(dm->coordinateDM, &cs);
127: if (cs) {PetscSectionGetChart(cs, NULL, &pEnd);}
128: MPI_Allreduce(&pEnd,&pEndMax,1,MPIU_INT,MPI_MAX,PetscObjectComm((PetscObject)dm));
129: if (pEndMax >= 0) {
130: DMClone(dm->coordinateDM, &ncdm);
131: DMSetDefaultSection(ncdm, cs);
132: DMSetCoordinateDM(*newdm, ncdm);
133: DMDestroy(&ncdm);
134: }
135: }
136: DMGetCoordinateDim(dm, &cdim);
137: DMSetCoordinateDim(*newdm, cdim);
138: DMGetCoordinatesLocal(dm, &coords);
139: if (coords) {
140: DMSetCoordinatesLocal(*newdm, coords);
141: } else {
142: DMGetCoordinates(dm, &coords);
143: if (coords) {DMSetCoordinates(*newdm, coords);}
144: }
145: {
146: PetscBool isper;
147: const PetscReal *maxCell, *L;
148: const DMBoundaryType *bd;
149: DMGetPeriodicity(dm, &isper, &maxCell, &L, &bd);
150: DMSetPeriodicity(*newdm, isper, maxCell, L, bd);
151: }
152: return(0);
153: }
155: /*@C
156: DMSetVecType - Sets the type of vector created with DMCreateLocalVector() and DMCreateGlobalVector()
158: Logically Collective on DM
160: Input Parameter:
161: + da - initial distributed array
162: . ctype - the vector type, currently either VECSTANDARD or VECCUSP
164: Options Database:
165: . -dm_vec_type ctype
167: Level: intermediate
169: .seealso: DMCreate(), DMDestroy(), DM, DMDAInterpolationType, VecType, DMGetVecType()
170: @*/
171: PetscErrorCode DMSetVecType(DM da,VecType ctype)
172: {
177: PetscFree(da->vectype);
178: PetscStrallocpy(ctype,(char**)&da->vectype);
179: return(0);
180: }
182: /*@C
183: DMGetVecType - Gets the type of vector created with DMCreateLocalVector() and DMCreateGlobalVector()
185: Logically Collective on DM
187: Input Parameter:
188: . da - initial distributed array
190: Output Parameter:
191: . ctype - the vector type
193: Level: intermediate
195: .seealso: DMCreate(), DMDestroy(), DM, DMDAInterpolationType, VecType
196: @*/
197: PetscErrorCode DMGetVecType(DM da,VecType *ctype)
198: {
201: *ctype = da->vectype;
202: return(0);
203: }
205: /*@
206: VecGetDM - Gets the DM defining the data layout of the vector
208: Not collective
210: Input Parameter:
211: . v - The Vec
213: Output Parameter:
214: . dm - The DM
216: Level: intermediate
218: .seealso: VecSetDM(), DMGetLocalVector(), DMGetGlobalVector(), DMSetVecType()
219: @*/
220: PetscErrorCode VecGetDM(Vec v, DM *dm)
221: {
227: PetscObjectQuery((PetscObject) v, "__PETSc_dm", (PetscObject*) dm);
228: return(0);
229: }
231: /*@
232: VecSetDM - Sets the DM defining the data layout of the vector.
234: Not collective
236: Input Parameters:
237: + v - The Vec
238: - dm - The DM
240: 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.
242: Level: intermediate
244: .seealso: VecGetDM(), DMGetLocalVector(), DMGetGlobalVector(), DMSetVecType()
245: @*/
246: PetscErrorCode VecSetDM(Vec v, DM dm)
247: {
253: PetscObjectCompose((PetscObject) v, "__PETSc_dm", (PetscObject) dm);
254: return(0);
255: }
257: /*@C
258: DMSetMatType - Sets the type of matrix created with DMCreateMatrix()
260: Logically Collective on DM
262: Input Parameter:
263: + dm - the DM context
264: - ctype - the matrix type
266: Options Database:
267: . -dm_mat_type ctype
269: Level: intermediate
271: .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMCreateMatrix(), DMSetMatrixPreallocateOnly(), MatType, DMGetMatType()
272: @*/
273: PetscErrorCode DMSetMatType(DM dm,MatType ctype)
274: {
279: PetscFree(dm->mattype);
280: PetscStrallocpy(ctype,(char**)&dm->mattype);
281: return(0);
282: }
284: /*@C
285: DMGetMatType - Gets the type of matrix created with DMCreateMatrix()
287: Logically Collective on DM
289: Input Parameter:
290: . dm - the DM context
292: Output Parameter:
293: . ctype - the matrix type
295: Options Database:
296: . -dm_mat_type ctype
298: Level: intermediate
300: .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMCreateMatrix(), DMSetMatrixPreallocateOnly(), MatType, DMSetMatType()
301: @*/
302: PetscErrorCode DMGetMatType(DM dm,MatType *ctype)
303: {
306: *ctype = dm->mattype;
307: return(0);
308: }
310: /*@
311: MatGetDM - Gets the DM defining the data layout of the matrix
313: Not collective
315: Input Parameter:
316: . A - The Mat
318: Output Parameter:
319: . dm - The DM
321: Level: intermediate
323: .seealso: MatSetDM(), DMCreateMatrix(), DMSetMatType()
324: @*/
325: PetscErrorCode MatGetDM(Mat A, DM *dm)
326: {
332: PetscObjectQuery((PetscObject) A, "__PETSc_dm", (PetscObject*) dm);
333: return(0);
334: }
336: /*@
337: MatSetDM - Sets the DM defining the data layout of the matrix
339: Not collective
341: Input Parameters:
342: + A - The Mat
343: - dm - The DM
345: Level: intermediate
347: .seealso: MatGetDM(), DMCreateMatrix(), DMSetMatType()
348: @*/
349: PetscErrorCode MatSetDM(Mat A, DM dm)
350: {
356: PetscObjectCompose((PetscObject) A, "__PETSc_dm", (PetscObject) dm);
357: return(0);
358: }
360: /*@C
361: DMSetOptionsPrefix - Sets the prefix used for searching for all
362: DM options in the database.
364: Logically Collective on DM
366: Input Parameter:
367: + da - the DM context
368: - prefix - the prefix to prepend to all option names
370: Notes:
371: A hyphen (-) must NOT be given at the beginning of the prefix name.
372: The first character of all runtime options is AUTOMATICALLY the hyphen.
374: Level: advanced
376: .keywords: DM, set, options, prefix, database
378: .seealso: DMSetFromOptions()
379: @*/
380: PetscErrorCode DMSetOptionsPrefix(DM dm,const char prefix[])
381: {
386: PetscObjectSetOptionsPrefix((PetscObject)dm,prefix);
387: if (dm->sf) {
388: PetscObjectSetOptionsPrefix((PetscObject)dm->sf,prefix);
389: }
390: if (dm->defaultSF) {
391: PetscObjectSetOptionsPrefix((PetscObject)dm->defaultSF,prefix);
392: }
393: return(0);
394: }
396: /*@C
397: DMAppendOptionsPrefix - Appends to the prefix used for searching for all
398: DM options in the database.
400: Logically Collective on DM
402: Input Parameters:
403: + dm - the DM context
404: - prefix - the prefix string to prepend to all DM option requests
406: Notes:
407: A hyphen (-) must NOT be given at the beginning of the prefix name.
408: The first character of all runtime options is AUTOMATICALLY the hyphen.
410: Level: advanced
412: .keywords: DM, append, options, prefix, database
414: .seealso: DMSetOptionsPrefix(), DMGetOptionsPrefix()
415: @*/
416: PetscErrorCode DMAppendOptionsPrefix(DM dm,const char prefix[])
417: {
422: PetscObjectAppendOptionsPrefix((PetscObject)dm,prefix);
423: return(0);
424: }
426: /*@C
427: DMGetOptionsPrefix - Gets the prefix used for searching for all
428: DM options in the database.
430: Not Collective
432: Input Parameters:
433: . dm - the DM context
435: Output Parameters:
436: . prefix - pointer to the prefix string used is returned
438: Notes: On the fortran side, the user should pass in a string 'prefix' of
439: sufficient length to hold the prefix.
441: Level: advanced
443: .keywords: DM, set, options, prefix, database
445: .seealso: DMSetOptionsPrefix(), DMAppendOptionsPrefix()
446: @*/
447: PetscErrorCode DMGetOptionsPrefix(DM dm,const char *prefix[])
448: {
453: PetscObjectGetOptionsPrefix((PetscObject)dm,prefix);
454: return(0);
455: }
457: static PetscErrorCode DMCountNonCyclicReferences(DM dm, PetscBool recurseCoarse, PetscBool recurseFine, PetscInt *ncrefct)
458: {
459: PetscInt i, refct = ((PetscObject) dm)->refct;
460: DMNamedVecLink nlink;
464: *ncrefct = 0;
465: /* count all the circular references of DM and its contained Vecs */
466: for (i=0; i<DM_MAX_WORK_VECTORS; i++) {
467: if (dm->localin[i]) refct--;
468: if (dm->globalin[i]) refct--;
469: }
470: for (nlink=dm->namedglobal; nlink; nlink=nlink->next) refct--;
471: for (nlink=dm->namedlocal; nlink; nlink=nlink->next) refct--;
472: if (dm->x) {
473: DM obj;
474: VecGetDM(dm->x, &obj);
475: if (obj == dm) refct--;
476: }
477: if (dm->coarseMesh && dm->coarseMesh->fineMesh == dm) {
478: refct--;
479: if (recurseCoarse) {
480: PetscInt coarseCount;
482: DMCountNonCyclicReferences(dm->coarseMesh, PETSC_TRUE, PETSC_FALSE,&coarseCount);
483: refct += coarseCount;
484: }
485: }
486: if (dm->fineMesh && dm->fineMesh->coarseMesh == dm) {
487: refct--;
488: if (recurseFine) {
489: PetscInt fineCount;
491: DMCountNonCyclicReferences(dm->fineMesh, PETSC_FALSE, PETSC_TRUE,&fineCount);
492: refct += fineCount;
493: }
494: }
495: *ncrefct = refct;
496: return(0);
497: }
499: PetscErrorCode DMDestroyLabelLinkList(DM dm)
500: {
504: if (!--(dm->labels->refct)) {
505: DMLabelLink next = dm->labels->next;
507: /* destroy the labels */
508: while (next) {
509: DMLabelLink tmp = next->next;
511: DMLabelDestroy(&next->label);
512: PetscFree(next);
513: next = tmp;
514: }
515: PetscFree(dm->labels);
516: }
517: return(0);
518: }
520: /*@
521: DMDestroy - Destroys a vector packer or DM.
523: Collective on DM
525: Input Parameter:
526: . dm - the DM object to destroy
528: Level: developer
530: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix()
532: @*/
533: PetscErrorCode DMDestroy(DM *dm)
534: {
535: PetscInt i, cnt;
536: DMNamedVecLink nlink,nnext;
540: if (!*dm) return(0);
543: /* count all non-cyclic references in the doubly-linked list of coarse<->fine meshes */
544: DMCountNonCyclicReferences(*dm,PETSC_TRUE,PETSC_TRUE,&cnt);
545: --((PetscObject)(*dm))->refct;
546: if (--cnt > 0) {*dm = 0; return(0);}
547: /*
548: Need this test because the dm references the vectors that
549: reference the dm, so destroying the dm calls destroy on the
550: vectors that cause another destroy on the dm
551: */
552: if (((PetscObject)(*dm))->refct < 0) return(0);
553: ((PetscObject) (*dm))->refct = 0;
554: for (i=0; i<DM_MAX_WORK_VECTORS; i++) {
555: if ((*dm)->localout[i]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Destroying a DM that has a local vector obtained with DMGetLocalVector()");
556: VecDestroy(&(*dm)->localin[i]);
557: }
558: nnext=(*dm)->namedglobal;
559: (*dm)->namedglobal = NULL;
560: for (nlink=nnext; nlink; nlink=nnext) { /* Destroy the named vectors */
561: nnext = nlink->next;
562: if (nlink->status != DMVEC_STATUS_IN) SETERRQ1(((PetscObject)*dm)->comm,PETSC_ERR_ARG_WRONGSTATE,"DM still has Vec named '%s' checked out",nlink->name);
563: PetscFree(nlink->name);
564: VecDestroy(&nlink->X);
565: PetscFree(nlink);
566: }
567: nnext=(*dm)->namedlocal;
568: (*dm)->namedlocal = NULL;
569: for (nlink=nnext; nlink; nlink=nnext) { /* Destroy the named local vectors */
570: nnext = nlink->next;
571: if (nlink->status != DMVEC_STATUS_IN) SETERRQ1(((PetscObject)*dm)->comm,PETSC_ERR_ARG_WRONGSTATE,"DM still has Vec named '%s' checked out",nlink->name);
572: PetscFree(nlink->name);
573: VecDestroy(&nlink->X);
574: PetscFree(nlink);
575: }
577: /* Destroy the list of hooks */
578: {
579: DMCoarsenHookLink link,next;
580: for (link=(*dm)->coarsenhook; link; link=next) {
581: next = link->next;
582: PetscFree(link);
583: }
584: (*dm)->coarsenhook = NULL;
585: }
586: {
587: DMRefineHookLink link,next;
588: for (link=(*dm)->refinehook; link; link=next) {
589: next = link->next;
590: PetscFree(link);
591: }
592: (*dm)->refinehook = NULL;
593: }
594: {
595: DMSubDomainHookLink link,next;
596: for (link=(*dm)->subdomainhook; link; link=next) {
597: next = link->next;
598: PetscFree(link);
599: }
600: (*dm)->subdomainhook = NULL;
601: }
602: {
603: DMGlobalToLocalHookLink link,next;
604: for (link=(*dm)->gtolhook; link; link=next) {
605: next = link->next;
606: PetscFree(link);
607: }
608: (*dm)->gtolhook = NULL;
609: }
610: {
611: DMLocalToGlobalHookLink link,next;
612: for (link=(*dm)->ltoghook; link; link=next) {
613: next = link->next;
614: PetscFree(link);
615: }
616: (*dm)->ltoghook = NULL;
617: }
618: /* Destroy the work arrays */
619: {
620: DMWorkLink link,next;
621: if ((*dm)->workout) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Work array still checked out");
622: for (link=(*dm)->workin; link; link=next) {
623: next = link->next;
624: PetscFree(link->mem);
625: PetscFree(link);
626: }
627: (*dm)->workin = NULL;
628: }
629: if (!--((*dm)->labels->refct)) {
630: DMLabelLink next = (*dm)->labels->next;
632: /* destroy the labels */
633: while (next) {
634: DMLabelLink tmp = next->next;
636: DMLabelDestroy(&next->label);
637: PetscFree(next);
638: next = tmp;
639: }
640: PetscFree((*dm)->labels);
641: }
642: {
643: DMBoundary next = (*dm)->boundary;
644: while (next) {
645: DMBoundary b = next;
647: next = b->next;
648: PetscFree(b);
649: }
650: }
652: PetscObjectDestroy(&(*dm)->dmksp);
653: PetscObjectDestroy(&(*dm)->dmsnes);
654: PetscObjectDestroy(&(*dm)->dmts);
656: if ((*dm)->ctx && (*dm)->ctxdestroy) {
657: (*(*dm)->ctxdestroy)(&(*dm)->ctx);
658: }
659: VecDestroy(&(*dm)->x);
660: MatFDColoringDestroy(&(*dm)->fd);
661: DMClearGlobalVectors(*dm);
662: ISLocalToGlobalMappingDestroy(&(*dm)->ltogmap);
663: PetscFree((*dm)->vectype);
664: PetscFree((*dm)->mattype);
666: PetscSectionDestroy(&(*dm)->defaultSection);
667: PetscSectionDestroy(&(*dm)->defaultGlobalSection);
668: PetscLayoutDestroy(&(*dm)->map);
669: PetscSectionDestroy(&(*dm)->defaultConstraintSection);
670: MatDestroy(&(*dm)->defaultConstraintMat);
671: PetscSFDestroy(&(*dm)->sf);
672: PetscSFDestroy(&(*dm)->defaultSF);
673: PetscSFDestroy(&(*dm)->sfNatural);
675: if ((*dm)->coarseMesh && (*dm)->coarseMesh->fineMesh == *dm) {
676: DMSetFineDM((*dm)->coarseMesh,NULL);
677: }
678: DMDestroy(&(*dm)->coarseMesh);
679: if ((*dm)->fineMesh && (*dm)->fineMesh->coarseMesh == *dm) {
680: DMSetCoarseDM((*dm)->fineMesh,NULL);
681: }
682: DMDestroy(&(*dm)->fineMesh);
683: DMDestroy(&(*dm)->coordinateDM);
684: VecDestroy(&(*dm)->coordinates);
685: VecDestroy(&(*dm)->coordinatesLocal);
686: PetscFree3((*dm)->L,(*dm)->maxCell,(*dm)->bdtype);
688: PetscDSDestroy(&(*dm)->prob);
689: DMDestroy(&(*dm)->dmBC);
690: /* if memory was published with SAWs then destroy it */
691: PetscObjectSAWsViewOff((PetscObject)*dm);
693: (*(*dm)->ops->destroy)(*dm);
694: /* We do not destroy (*dm)->data here so that we can reference count backend objects */
695: PetscHeaderDestroy(dm);
696: return(0);
697: }
699: /*@
700: DMSetUp - sets up the data structures inside a DM object
702: Collective on DM
704: Input Parameter:
705: . dm - the DM object to setup
707: Level: developer
709: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix()
711: @*/
712: PetscErrorCode DMSetUp(DM dm)
713: {
718: if (dm->setupcalled) return(0);
719: if (dm->ops->setup) {
720: (*dm->ops->setup)(dm);
721: }
722: dm->setupcalled = PETSC_TRUE;
723: return(0);
724: }
726: /*@
727: DMSetFromOptions - sets parameters in a DM from the options database
729: Collective on DM
731: Input Parameter:
732: . dm - the DM object to set options for
734: Options Database:
735: + -dm_preallocate_only - Only preallocate the matrix for DMCreateMatrix(), but do not fill it with zeros
736: . -dm_vec_type <type> - type of vector to create inside DM
737: . -dm_mat_type <type> - type of matrix to create inside DM
738: - -dm_coloring_type - <global or ghosted>
740: Level: developer
742: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix()
744: @*/
745: PetscErrorCode DMSetFromOptions(DM dm)
746: {
747: char typeName[256];
748: PetscBool flg;
753: if (dm->prob) {
754: PetscDSSetFromOptions(dm->prob);
755: }
756: if (dm->sf) {
757: PetscSFSetFromOptions(dm->sf);
758: }
759: if (dm->defaultSF) {
760: PetscSFSetFromOptions(dm->defaultSF);
761: }
762: PetscObjectOptionsBegin((PetscObject)dm);
763: PetscOptionsBool("-dm_preallocate_only","only preallocate matrix, but do not set column indices","DMSetMatrixPreallocateOnly",dm->prealloc_only,&dm->prealloc_only,NULL);
764: PetscOptionsFList("-dm_vec_type","Vector type used for created vectors","DMSetVecType",VecList,dm->vectype,typeName,256,&flg);
765: if (flg) {
766: DMSetVecType(dm,typeName);
767: }
768: PetscOptionsFList("-dm_mat_type","Matrix type used for created matrices","DMSetMatType",MatList,dm->mattype ? dm->mattype : typeName,typeName,sizeof(typeName),&flg);
769: if (flg) {
770: DMSetMatType(dm,typeName);
771: }
772: PetscOptionsEnum("-dm_is_coloring_type","Global or local coloring of Jacobian","ISColoringType",ISColoringTypes,(PetscEnum)dm->coloringtype,(PetscEnum*)&dm->coloringtype,NULL);
773: if (dm->ops->setfromoptions) {
774: (*dm->ops->setfromoptions)(PetscOptionsObject,dm);
775: }
776: /* process any options handlers added with PetscObjectAddOptionsHandler() */
777: PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject) dm);
778: PetscOptionsEnd();
779: return(0);
780: }
782: /*@C
783: DMView - Views a DM
785: Collective on DM
787: Input Parameter:
788: + dm - the DM object to view
789: - v - the viewer
791: Level: beginner
793: .seealso DMDestroy(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix()
795: @*/
796: PetscErrorCode DMView(DM dm,PetscViewer v)
797: {
798: PetscErrorCode ierr;
799: PetscBool isbinary;
800: PetscMPIInt size;
801: PetscViewerFormat format;
805: if (!v) {
806: PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)dm),&v);
807: }
808: PetscViewerGetFormat(v,&format);
809: MPI_Comm_size(PetscObjectComm((PetscObject)dm),&size);
810: if (size == 1 && format == PETSC_VIEWER_LOAD_BALANCE) return(0);
811: PetscObjectPrintClassNamePrefixType((PetscObject)dm,v);
812: PetscObjectTypeCompare((PetscObject)v,PETSCVIEWERBINARY,&isbinary);
813: if (isbinary) {
814: PetscInt classid = DM_FILE_CLASSID;
815: char type[256];
817: PetscViewerBinaryWrite(v,&classid,1,PETSC_INT,PETSC_FALSE);
818: PetscStrncpy(type,((PetscObject)dm)->type_name,256);
819: PetscViewerBinaryWrite(v,type,256,PETSC_CHAR,PETSC_FALSE);
820: }
821: if (dm->ops->view) {
822: (*dm->ops->view)(dm,v);
823: }
824: return(0);
825: }
827: /*@
828: DMCreateGlobalVector - Creates a global vector from a DM object
830: Collective on DM
832: Input Parameter:
833: . dm - the DM object
835: Output Parameter:
836: . vec - the global vector
838: Level: beginner
840: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix()
842: @*/
843: PetscErrorCode DMCreateGlobalVector(DM dm,Vec *vec)
844: {
849: (*dm->ops->createglobalvector)(dm,vec);
850: return(0);
851: }
853: /*@
854: DMCreateLocalVector - Creates a local vector from a DM object
856: Not Collective
858: Input Parameter:
859: . dm - the DM object
861: Output Parameter:
862: . vec - the local vector
864: Level: beginner
866: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix()
868: @*/
869: PetscErrorCode DMCreateLocalVector(DM dm,Vec *vec)
870: {
875: (*dm->ops->createlocalvector)(dm,vec);
876: return(0);
877: }
879: /*@
880: DMGetLocalToGlobalMapping - Accesses the local-to-global mapping in a DM.
882: Collective on DM
884: Input Parameter:
885: . dm - the DM that provides the mapping
887: Output Parameter:
888: . ltog - the mapping
890: Level: intermediate
892: Notes:
893: This mapping can then be used by VecSetLocalToGlobalMapping() or
894: MatSetLocalToGlobalMapping().
896: .seealso: DMCreateLocalVector()
897: @*/
898: PetscErrorCode DMGetLocalToGlobalMapping(DM dm,ISLocalToGlobalMapping *ltog)
899: {
900: PetscInt bs = -1, bsLocal, bsMin, bsMax;
906: if (!dm->ltogmap) {
907: PetscSection section, sectionGlobal;
909: DMGetDefaultSection(dm, §ion);
910: if (section) {
911: const PetscInt *cdofs;
912: PetscInt *ltog;
913: PetscInt pStart, pEnd, n, p, k, l;
915: DMGetDefaultGlobalSection(dm, §ionGlobal);
916: PetscSectionGetChart(section, &pStart, &pEnd);
917: PetscSectionGetStorageSize(section, &n);
918: PetscMalloc1(n, <og); /* We want the local+overlap size */
919: for (p = pStart, l = 0; p < pEnd; ++p) {
920: PetscInt bdof, cdof, dof, off, c, cind = 0;
922: /* Should probably use constrained dofs */
923: PetscSectionGetDof(section, p, &dof);
924: PetscSectionGetConstraintDof(section, p, &cdof);
925: PetscSectionGetConstraintIndices(section, p, &cdofs);
926: PetscSectionGetOffset(sectionGlobal, p, &off);
927: /* If you have dofs, and constraints, and they are unequal, we set the blocksize to 1 */
928: bdof = cdof && (dof-cdof) ? 1 : dof;
929: if (dof) {
930: if (bs < 0) {bs = bdof;}
931: else if (bs != bdof) {bs = 1;}
932: }
933: for (c = 0; c < dof; ++c, ++l) {
934: if ((cind < cdof) && (c == cdofs[cind])) ltog[l] = off < 0 ? off-c : off+c;
935: else ltog[l] = (off < 0 ? -(off+1) : off) + c;
936: }
937: }
938: /* Must have same blocksize on all procs (some might have no points) */
939: bsLocal = bs;
940: MPIU_Allreduce(&bsLocal, &bsMax, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)dm));
941: bsLocal = bs < 0 ? bsMax : bs;
942: MPIU_Allreduce(&bsLocal, &bsMin, 1, MPIU_INT, MPI_MIN, PetscObjectComm((PetscObject)dm));
943: if (bsMin != bsMax) {bs = 1;}
944: else {bs = bsMax;}
945: bs = bs < 0 ? 1 : bs;
946: /* Must reduce indices by blocksize */
947: if (bs > 1) {
948: for (l = 0, k = 0; l < n; l += bs, ++k) ltog[k] = ltog[l]/bs;
949: n /= bs;
950: }
951: ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)dm), bs, n, ltog, PETSC_OWN_POINTER, &dm->ltogmap);
952: PetscLogObjectParent((PetscObject)dm, (PetscObject)dm->ltogmap);
953: } else {
954: if (!dm->ops->getlocaltoglobalmapping) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM can not create LocalToGlobalMapping");
955: (*dm->ops->getlocaltoglobalmapping)(dm);
956: }
957: }
958: *ltog = dm->ltogmap;
959: return(0);
960: }
962: /*@
963: DMGetBlockSize - Gets the inherent block size associated with a DM
965: Not Collective
967: Input Parameter:
968: . dm - the DM with block structure
970: Output Parameter:
971: . bs - the block size, 1 implies no exploitable block structure
973: Level: intermediate
975: .seealso: ISCreateBlock(), VecSetBlockSize(), MatSetBlockSize(), DMGetLocalToGlobalMapping()
976: @*/
977: PetscErrorCode DMGetBlockSize(DM dm,PetscInt *bs)
978: {
982: if (dm->bs < 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"DM does not have enough information to provide a block size yet");
983: *bs = dm->bs;
984: return(0);
985: }
987: /*@
988: DMCreateInterpolation - Gets interpolation matrix between two DM objects
990: Collective on DM
992: Input Parameter:
993: + dm1 - the DM object
994: - dm2 - the second, finer DM object
996: Output Parameter:
997: + mat - the interpolation
998: - vec - the scaling (optional)
1000: Level: developer
1002: 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
1003: DMCoarsen(). The coordinates set into the DMDA are completely ignored in computing the interpolation.
1005: For DMDA objects you can use this interpolation (more precisely the interpolation from the DMGetCoordinateDM()) to interpolate the mesh coordinate vectors
1006: EXCEPT in the periodic case where it does not make sense since the coordinate vectors are not periodic.
1009: .seealso DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateColoring(), DMCreateMatrix(), DMRefine(), DMCoarsen(), DMCreateRestriction()
1011: @*/
1012: PetscErrorCode DMCreateInterpolation(DM dm1,DM dm2,Mat *mat,Vec *vec)
1013: {
1019: PetscLogEventBegin(DM_CreateInterpolation,dm1,dm2,0,0);
1020: (*dm1->ops->createinterpolation)(dm1,dm2,mat,vec);
1021: PetscLogEventEnd(DM_CreateInterpolation,dm1,dm2,0,0);
1022: return(0);
1023: }
1025: /*@
1026: DMCreateRestriction - Gets restriction matrix between two DM objects
1028: Collective on DM
1030: Input Parameter:
1031: + dm1 - the DM object
1032: - dm2 - the second, finer DM object
1034: Output Parameter:
1035: . mat - the restriction
1038: Level: developer
1040: 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
1041: DMCoarsen(). The coordinates set into the DMDA are completely ignored in computing the interpolation.
1044: .seealso DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateColoring(), DMCreateMatrix(), DMRefine(), DMCoarsen(), DMCreateInterpolation()
1046: @*/
1047: PetscErrorCode DMCreateRestriction(DM dm1,DM dm2,Mat *mat)
1048: {
1054: PetscLogEventBegin(DM_CreateRestriction,dm1,dm2,0,0);
1055: if (!dm1->ops->createrestriction) SETERRQ(PetscObjectComm((PetscObject)dm1),PETSC_ERR_SUP,"DMCreateRestriction not implemented for this type");
1056: (*dm1->ops->createrestriction)(dm1,dm2,mat);
1057: PetscLogEventEnd(DM_CreateRestriction,dm1,dm2,0,0);
1058: return(0);
1059: }
1061: /*@
1062: DMCreateInjection - Gets injection matrix between two DM objects
1064: Collective on DM
1066: Input Parameter:
1067: + dm1 - the DM object
1068: - dm2 - the second, finer DM object
1070: Output Parameter:
1071: . mat - the injection
1073: Level: developer
1075: Notes: For DMDA objects this only works for "uniform refinement", that is the refined mesh was obtained DMRefine() or the coarse mesh was obtained by
1076: DMCoarsen(). The coordinates set into the DMDA are completely ignored in computing the injection.
1078: .seealso DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateColoring(), DMCreateMatrix(), DMCreateInterpolation()
1080: @*/
1081: PetscErrorCode DMCreateInjection(DM dm1,DM dm2,Mat *mat)
1082: {
1088: if (!dm1->ops->getinjection) SETERRQ(PetscObjectComm((PetscObject)dm1),PETSC_ERR_SUP,"DMCreateInjection not implemented for this type");
1089: (*dm1->ops->getinjection)(dm1,dm2,mat);
1090: return(0);
1091: }
1093: /*@
1094: DMCreateColoring - Gets coloring for a DM
1096: Collective on DM
1098: Input Parameter:
1099: + dm - the DM object
1100: - ctype - IS_COLORING_LOCAL or IS_COLORING_GLOBAL
1102: Output Parameter:
1103: . coloring - the coloring
1105: Level: developer
1107: .seealso DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateMatrix(), DMSetMatType()
1109: @*/
1110: PetscErrorCode DMCreateColoring(DM dm,ISColoringType ctype,ISColoring *coloring)
1111: {
1116: if (!dm->ops->getcoloring) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"No coloring for this type of DM yet");
1117: (*dm->ops->getcoloring)(dm,ctype,coloring);
1118: return(0);
1119: }
1121: /*@
1122: DMCreateMatrix - Gets empty Jacobian for a DM
1124: Collective on DM
1126: Input Parameter:
1127: . dm - the DM object
1129: Output Parameter:
1130: . mat - the empty Jacobian
1132: Level: beginner
1134: Notes: This properly preallocates the number of nonzeros in the sparse matrix so you
1135: do not need to do it yourself.
1137: By default it also sets the nonzero structure and puts in the zero entries. To prevent setting
1138: the nonzero pattern call DMSetMatrixPreallocateOnly()
1140: For structured grid problems, when you call MatView() on this matrix it is displayed using the global natural ordering, NOT in the ordering used
1141: internally by PETSc.
1143: For structured grid problems, in general it is easiest to use MatSetValuesStencil() or MatSetValuesLocal() to put values into the matrix because MatSetValues() requires
1144: the indices for the global numbering for DMDAs which is complicated.
1146: .seealso DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMSetMatType()
1148: @*/
1149: PetscErrorCode DMCreateMatrix(DM dm,Mat *mat)
1150: {
1155: MatInitializePackage();
1158: (*dm->ops->creatematrix)(dm,mat);
1159: return(0);
1160: }
1162: /*@
1163: DMSetMatrixPreallocateOnly - When DMCreateMatrix() is called the matrix will be properly
1164: preallocated but the nonzero structure and zero values will not be set.
1166: Logically Collective on DM
1168: Input Parameter:
1169: + dm - the DM
1170: - only - PETSC_TRUE if only want preallocation
1172: Level: developer
1173: .seealso DMCreateMatrix(), DMSetMatrixStructureOnly()
1174: @*/
1175: PetscErrorCode DMSetMatrixPreallocateOnly(DM dm, PetscBool only)
1176: {
1179: dm->prealloc_only = only;
1180: return(0);
1181: }
1183: /*@
1184: DMSetMatrixStructureOnly - When DMCreateMatrix() is called, the matrix structure will be created
1185: but the array for values will not be allocated.
1187: Logically Collective on DM
1189: Input Parameter:
1190: + dm - the DM
1191: - only - PETSC_TRUE if only want matrix stucture
1193: Level: developer
1194: .seealso DMCreateMatrix(), DMSetMatrixPreallocateOnly()
1195: @*/
1196: PetscErrorCode DMSetMatrixStructureOnly(DM dm, PetscBool only)
1197: {
1200: dm->structure_only = only;
1201: return(0);
1202: }
1204: /*@C
1205: DMGetWorkArray - Gets a work array guaranteed to be at least the input size, restore with DMRestoreWorkArray()
1207: Not Collective
1209: Input Parameters:
1210: + dm - the DM object
1211: . count - The minium size
1212: - dtype - data type (PETSC_REAL, PETSC_SCALAR, PETSC_INT)
1214: Output Parameter:
1215: . array - the work array
1217: Level: developer
1219: .seealso DMDestroy(), DMCreate()
1220: @*/
1221: PetscErrorCode DMGetWorkArray(DM dm,PetscInt count,PetscDataType dtype,void *mem)
1222: {
1224: DMWorkLink link;
1225: size_t dsize;
1230: if (dm->workin) {
1231: link = dm->workin;
1232: dm->workin = dm->workin->next;
1233: } else {
1234: PetscNewLog(dm,&link);
1235: }
1236: PetscDataTypeGetSize(dtype,&dsize);
1237: if (dsize*count > link->bytes) {
1238: PetscFree(link->mem);
1239: PetscMalloc(dsize*count,&link->mem);
1240: link->bytes = dsize*count;
1241: }
1242: link->next = dm->workout;
1243: dm->workout = link;
1244: *(void**)mem = link->mem;
1245: return(0);
1246: }
1248: /*@C
1249: DMRestoreWorkArray - Restores a work array guaranteed to be at least the input size, restore with DMRestoreWorkArray()
1251: Not Collective
1253: Input Parameters:
1254: + dm - the DM object
1255: . count - The minium size
1256: - dtype - data type (PETSC_REAL, PETSC_SCALAR, PETSC_INT)
1258: Output Parameter:
1259: . array - the work array
1261: Level: developer
1263: .seealso DMDestroy(), DMCreate()
1264: @*/
1265: PetscErrorCode DMRestoreWorkArray(DM dm,PetscInt count,PetscDataType dtype,void *mem)
1266: {
1267: DMWorkLink *p,link;
1272: for (p=&dm->workout; (link=*p); p=&link->next) {
1273: if (link->mem == *(void**)mem) {
1274: *p = link->next;
1275: link->next = dm->workin;
1276: dm->workin = link;
1277: *(void**)mem = NULL;
1278: return(0);
1279: }
1280: }
1281: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Array was not checked out");
1282: }
1284: PetscErrorCode DMSetNullSpaceConstructor(DM dm, PetscInt field, PetscErrorCode (*nullsp)(DM dm, PetscInt field, MatNullSpace *nullSpace))
1285: {
1288: if (field >= 10) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle %d >= 10 fields", field);
1289: dm->nullspaceConstructors[field] = nullsp;
1290: return(0);
1291: }
1293: /*@C
1294: DMCreateFieldIS - Creates a set of IS objects with the global indices of dofs for each field
1296: Not collective
1298: Input Parameter:
1299: . dm - the DM object
1301: Output Parameters:
1302: + numFields - The number of fields (or NULL if not requested)
1303: . fieldNames - The name for each field (or NULL if not requested)
1304: - fields - The global indices for each field (or NULL if not requested)
1306: Level: intermediate
1308: Notes:
1309: The user is responsible for freeing all requested arrays. In particular, every entry of names should be freed with
1310: PetscFree(), every entry of fields should be destroyed with ISDestroy(), and both arrays should be freed with
1311: PetscFree().
1313: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix()
1314: @*/
1315: PetscErrorCode DMCreateFieldIS(DM dm, PetscInt *numFields, char ***fieldNames, IS **fields)
1316: {
1317: PetscSection section, sectionGlobal;
1322: if (numFields) {
1324: *numFields = 0;
1325: }
1326: if (fieldNames) {
1328: *fieldNames = NULL;
1329: }
1330: if (fields) {
1332: *fields = NULL;
1333: }
1334: DMGetDefaultSection(dm, §ion);
1335: if (section) {
1336: PetscInt *fieldSizes, **fieldIndices;
1337: PetscInt nF, f, pStart, pEnd, p;
1339: DMGetDefaultGlobalSection(dm, §ionGlobal);
1340: PetscSectionGetNumFields(section, &nF);
1341: PetscMalloc2(nF,&fieldSizes,nF,&fieldIndices);
1342: PetscSectionGetChart(sectionGlobal, &pStart, &pEnd);
1343: for (f = 0; f < nF; ++f) {
1344: fieldSizes[f] = 0;
1345: }
1346: for (p = pStart; p < pEnd; ++p) {
1347: PetscInt gdof;
1349: PetscSectionGetDof(sectionGlobal, p, &gdof);
1350: if (gdof > 0) {
1351: for (f = 0; f < nF; ++f) {
1352: PetscInt fdof, fcdof;
1354: PetscSectionGetFieldDof(section, p, f, &fdof);
1355: PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);
1356: fieldSizes[f] += fdof-fcdof;
1357: }
1358: }
1359: }
1360: for (f = 0; f < nF; ++f) {
1361: PetscMalloc1(fieldSizes[f], &fieldIndices[f]);
1362: fieldSizes[f] = 0;
1363: }
1364: for (p = pStart; p < pEnd; ++p) {
1365: PetscInt gdof, goff;
1367: PetscSectionGetDof(sectionGlobal, p, &gdof);
1368: if (gdof > 0) {
1369: PetscSectionGetOffset(sectionGlobal, p, &goff);
1370: for (f = 0; f < nF; ++f) {
1371: PetscInt fdof, fcdof, fc;
1373: PetscSectionGetFieldDof(section, p, f, &fdof);
1374: PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);
1375: for (fc = 0; fc < fdof-fcdof; ++fc, ++fieldSizes[f]) {
1376: fieldIndices[f][fieldSizes[f]] = goff++;
1377: }
1378: }
1379: }
1380: }
1381: if (numFields) *numFields = nF;
1382: if (fieldNames) {
1383: PetscMalloc1(nF, fieldNames);
1384: for (f = 0; f < nF; ++f) {
1385: const char *fieldName;
1387: PetscSectionGetFieldName(section, f, &fieldName);
1388: PetscStrallocpy(fieldName, (char**) &(*fieldNames)[f]);
1389: }
1390: }
1391: if (fields) {
1392: PetscMalloc1(nF, fields);
1393: for (f = 0; f < nF; ++f) {
1394: ISCreateGeneral(PetscObjectComm((PetscObject)dm), fieldSizes[f], fieldIndices[f], PETSC_OWN_POINTER, &(*fields)[f]);
1395: }
1396: }
1397: PetscFree2(fieldSizes,fieldIndices);
1398: } else if (dm->ops->createfieldis) {
1399: (*dm->ops->createfieldis)(dm, numFields, fieldNames, fields);
1400: }
1401: return(0);
1402: }
1405: /*@C
1406: DMCreateFieldDecomposition - Returns a list of IS objects defining a decomposition of a problem into subproblems
1407: corresponding to different fields: each IS contains the global indices of the dofs of the
1408: corresponding field. The optional list of DMs define the DM for each subproblem.
1409: Generalizes DMCreateFieldIS().
1411: Not collective
1413: Input Parameter:
1414: . dm - the DM object
1416: Output Parameters:
1417: + len - The number of subproblems in the field decomposition (or NULL if not requested)
1418: . namelist - The name for each field (or NULL if not requested)
1419: . islist - The global indices for each field (or NULL if not requested)
1420: - dmlist - The DMs for each field subproblem (or NULL, if not requested; if NULL is returned, no DMs are defined)
1422: Level: intermediate
1424: Notes:
1425: The user is responsible for freeing all requested arrays. In particular, every entry of names should be freed with
1426: PetscFree(), every entry of is should be destroyed with ISDestroy(), every entry of dm should be destroyed with DMDestroy(),
1427: and all of the arrays should be freed with PetscFree().
1429: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateFieldIS()
1430: @*/
1431: PetscErrorCode DMCreateFieldDecomposition(DM dm, PetscInt *len, char ***namelist, IS **islist, DM **dmlist)
1432: {
1437: if (len) {
1439: *len = 0;
1440: }
1441: if (namelist) {
1443: *namelist = 0;
1444: }
1445: if (islist) {
1447: *islist = 0;
1448: }
1449: if (dmlist) {
1451: *dmlist = 0;
1452: }
1453: /*
1454: Is it a good idea to apply the following check across all impls?
1455: Perhaps some impls can have a well-defined decomposition before DMSetUp?
1456: This, however, follows the general principle that accessors are not well-behaved until the object is set up.
1457: */
1458: if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE, "Decomposition defined only after DMSetUp");
1459: if (!dm->ops->createfielddecomposition) {
1460: PetscSection section;
1461: PetscInt numFields, f;
1463: DMGetDefaultSection(dm, §ion);
1464: if (section) {PetscSectionGetNumFields(section, &numFields);}
1465: if (section && numFields && dm->ops->createsubdm) {
1466: if (len) *len = numFields;
1467: if (namelist) {PetscMalloc1(numFields,namelist);}
1468: if (islist) {PetscMalloc1(numFields,islist);}
1469: if (dmlist) {PetscMalloc1(numFields,dmlist);}
1470: for (f = 0; f < numFields; ++f) {
1471: const char *fieldName;
1473: DMCreateSubDM(dm, 1, &f, islist ? &(*islist)[f] : NULL, dmlist ? &(*dmlist)[f] : NULL);
1474: if (namelist) {
1475: PetscSectionGetFieldName(section, f, &fieldName);
1476: PetscStrallocpy(fieldName, (char**) &(*namelist)[f]);
1477: }
1478: }
1479: } else {
1480: DMCreateFieldIS(dm, len, namelist, islist);
1481: /* By default there are no DMs associated with subproblems. */
1482: if (dmlist) *dmlist = NULL;
1483: }
1484: } else {
1485: (*dm->ops->createfielddecomposition)(dm,len,namelist,islist,dmlist);
1486: }
1487: return(0);
1488: }
1490: /*@
1491: DMCreateSubDM - Returns an IS and DM encapsulating a subproblem defined by the fields passed in.
1492: The fields are defined by DMCreateFieldIS().
1494: Not collective
1496: Input Parameters:
1497: + dm - the DM object
1498: . numFields - the number of fields in this subproblem
1499: - fields - the fields in the subproblem
1501: Output Parameters:
1502: + is - the global indices for the subproblem
1503: - dm - the DM for the subproblem
1505: Level: intermediate
1507: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateFieldIS()
1508: @*/
1509: PetscErrorCode DMCreateSubDM(DM dm, PetscInt numFields, PetscInt fields[], IS *is, DM *subdm)
1510: {
1518: if (dm->ops->createsubdm) {
1519: (*dm->ops->createsubdm)(dm, numFields, fields, is, subdm);
1520: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "This type has no DMCreateSubDM implementation defined");
1521: return(0);
1522: }
1525: /*@C
1526: DMCreateDomainDecomposition - Returns lists of IS objects defining a decomposition of a problem into subproblems
1527: corresponding to restrictions to pairs nested subdomains: each IS contains the global
1528: indices of the dofs of the corresponding subdomains. The inner subdomains conceptually
1529: define a nonoverlapping covering, while outer subdomains can overlap.
1530: The optional list of DMs define the DM for each subproblem.
1532: Not collective
1534: Input Parameter:
1535: . dm - the DM object
1537: Output Parameters:
1538: + len - The number of subproblems in the domain decomposition (or NULL if not requested)
1539: . namelist - The name for each subdomain (or NULL if not requested)
1540: . innerislist - The global indices for each inner subdomain (or NULL, if not requested)
1541: . outerislist - The global indices for each outer subdomain (or NULL, if not requested)
1542: - dmlist - The DMs for each subdomain subproblem (or NULL, if not requested; if NULL is returned, no DMs are defined)
1544: Level: intermediate
1546: Notes:
1547: The user is responsible for freeing all requested arrays. In particular, every entry of names should be freed with
1548: PetscFree(), every entry of is should be destroyed with ISDestroy(), every entry of dm should be destroyed with DMDestroy(),
1549: and all of the arrays should be freed with PetscFree().
1551: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateDomainDecompositionDM(), DMCreateFieldDecomposition()
1552: @*/
1553: PetscErrorCode DMCreateDomainDecomposition(DM dm, PetscInt *len, char ***namelist, IS **innerislist, IS **outerislist, DM **dmlist)
1554: {
1555: PetscErrorCode ierr;
1556: DMSubDomainHookLink link;
1557: PetscInt i,l;
1566: /*
1567: Is it a good idea to apply the following check across all impls?
1568: Perhaps some impls can have a well-defined decomposition before DMSetUp?
1569: This, however, follows the general principle that accessors are not well-behaved until the object is set up.
1570: */
1571: if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE, "Decomposition defined only after DMSetUp");
1572: if (dm->ops->createdomaindecomposition) {
1573: (*dm->ops->createdomaindecomposition)(dm,&l,namelist,innerislist,outerislist,dmlist);
1574: /* copy subdomain hooks and context over to the subdomain DMs */
1575: if (dmlist && *dmlist) {
1576: for (i = 0; i < l; i++) {
1577: for (link=dm->subdomainhook; link; link=link->next) {
1578: if (link->ddhook) {(*link->ddhook)(dm,(*dmlist)[i],link->ctx);}
1579: }
1580: if (dm->ctx) (*dmlist)[i]->ctx = dm->ctx;
1581: }
1582: }
1583: if (len) *len = l;
1584: }
1585: return(0);
1586: }
1589: /*@C
1590: DMCreateDomainDecompositionScatters - Returns scatters to the subdomain vectors from the global vector
1592: Not collective
1594: Input Parameters:
1595: + dm - the DM object
1596: . n - the number of subdomain scatters
1597: - subdms - the local subdomains
1599: Output Parameters:
1600: + n - the number of scatters returned
1601: . iscat - scatter from global vector to nonoverlapping global vector entries on subdomain
1602: . oscat - scatter from global vector to overlapping global vector entries on subdomain
1603: - gscat - scatter from global vector to local vector on subdomain (fills in ghosts)
1605: Notes: This is an alternative to the iis and ois arguments in DMCreateDomainDecomposition that allow for the solution
1606: of general nonlinear problems with overlapping subdomain methods. While merely having index sets that enable subsets
1607: of the residual equations to be created is fine for linear problems, nonlinear problems require local assembly of
1608: solution and residual data.
1610: Level: developer
1612: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateFieldIS()
1613: @*/
1614: PetscErrorCode DMCreateDomainDecompositionScatters(DM dm,PetscInt n,DM *subdms,VecScatter **iscat,VecScatter **oscat,VecScatter **gscat)
1615: {
1621: if (dm->ops->createddscatters) {
1622: (*dm->ops->createddscatters)(dm,n,subdms,iscat,oscat,gscat);
1623: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "This type has no DMCreateDomainDecompositionScatter implementation defined");
1624: return(0);
1625: }
1627: /*@
1628: DMRefine - Refines a DM object
1630: Collective on DM
1632: Input Parameter:
1633: + dm - the DM object
1634: - comm - the communicator to contain the new DM object (or MPI_COMM_NULL)
1636: Output Parameter:
1637: . dmf - the refined DM, or NULL
1639: Note: If no refinement was done, the return value is NULL
1641: Level: developer
1643: .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation()
1644: @*/
1645: PetscErrorCode DMRefine(DM dm,MPI_Comm comm,DM *dmf)
1646: {
1647: PetscErrorCode ierr;
1648: DMRefineHookLink link;
1652: PetscLogEventBegin(DM_Refine,dm,0,0,0);
1653: if (!dm->ops->refine) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"This DM cannot refine");
1654: (*dm->ops->refine)(dm,comm,dmf);
1655: if (*dmf) {
1656: (*dmf)->ops->creatematrix = dm->ops->creatematrix;
1658: PetscObjectCopyFortranFunctionPointers((PetscObject)dm,(PetscObject)*dmf);
1660: (*dmf)->ctx = dm->ctx;
1661: (*dmf)->leveldown = dm->leveldown;
1662: (*dmf)->levelup = dm->levelup + 1;
1664: DMSetMatType(*dmf,dm->mattype);
1665: for (link=dm->refinehook; link; link=link->next) {
1666: if (link->refinehook) {
1667: (*link->refinehook)(dm,*dmf,link->ctx);
1668: }
1669: }
1670: }
1671: PetscLogEventEnd(DM_Refine,dm,0,0,0);
1672: return(0);
1673: }
1675: /*@C
1676: DMRefineHookAdd - adds a callback to be run when interpolating a nonlinear problem to a finer grid
1678: Logically Collective
1680: Input Arguments:
1681: + coarse - nonlinear solver context on which to run a hook when restricting to a coarser level
1682: . refinehook - function to run when setting up a coarser level
1683: . interphook - function to run to update data on finer levels (once per SNESSolve())
1684: - ctx - [optional] user-defined context for provide data for the hooks (may be NULL)
1686: Calling sequence of refinehook:
1687: $ refinehook(DM coarse,DM fine,void *ctx);
1689: + coarse - coarse level DM
1690: . fine - fine level DM to interpolate problem to
1691: - ctx - optional user-defined function context
1693: Calling sequence for interphook:
1694: $ interphook(DM coarse,Mat interp,DM fine,void *ctx)
1696: + coarse - coarse level DM
1697: . interp - matrix interpolating a coarse-level solution to the finer grid
1698: . fine - fine level DM to update
1699: - ctx - optional user-defined function context
1701: Level: advanced
1703: Notes:
1704: This function is only needed if auxiliary data needs to be passed to fine grids while grid sequencing
1706: If this function is called multiple times, the hooks will be run in the order they are added.
1708: This function is currently not available from Fortran.
1710: .seealso: DMCoarsenHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
1711: @*/
1712: PetscErrorCode DMRefineHookAdd(DM coarse,PetscErrorCode (*refinehook)(DM,DM,void*),PetscErrorCode (*interphook)(DM,Mat,DM,void*),void *ctx)
1713: {
1714: PetscErrorCode ierr;
1715: DMRefineHookLink link,*p;
1719: for (p=&coarse->refinehook; *p; p=&(*p)->next) { /* Scan to the end of the current list of hooks */
1720: if ((*p)->refinehook == refinehook && (*p)->interphook == interphook && (*p)->ctx == ctx) return(0);
1721: }
1722: PetscNew(&link);
1723: link->refinehook = refinehook;
1724: link->interphook = interphook;
1725: link->ctx = ctx;
1726: link->next = NULL;
1727: *p = link;
1728: return(0);
1729: }
1731: /*@C
1732: DMRefineHookRemove - remove a callback from the list of hooks to be run when interpolating a nonlinear problem to a finer grid
1734: Logically Collective
1736: Input Arguments:
1737: + coarse - nonlinear solver context on which to run a hook when restricting to a coarser level
1738: . refinehook - function to run when setting up a coarser level
1739: . interphook - function to run to update data on finer levels (once per SNESSolve())
1740: - ctx - [optional] user-defined context for provide data for the hooks (may be NULL)
1742: Level: advanced
1744: Notes:
1745: This function does nothing if the hook is not in the list.
1747: This function is currently not available from Fortran.
1749: .seealso: DMCoarsenHookRemove(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
1750: @*/
1751: PetscErrorCode DMRefineHookRemove(DM coarse,PetscErrorCode (*refinehook)(DM,DM,void*),PetscErrorCode (*interphook)(DM,Mat,DM,void*),void *ctx)
1752: {
1753: PetscErrorCode ierr;
1754: DMRefineHookLink link,*p;
1758: for (p=&coarse->refinehook; *p; p=&(*p)->next) { /* Search the list of current hooks */
1759: if ((*p)->refinehook == refinehook && (*p)->interphook == interphook && (*p)->ctx == ctx) {
1760: link = *p;
1761: *p = link->next;
1762: PetscFree(link);
1763: break;
1764: }
1765: }
1766: return(0);
1767: }
1769: /*@
1770: DMInterpolate - interpolates user-defined problem data to a finer DM by running hooks registered by DMRefineHookAdd()
1772: Collective if any hooks are
1774: Input Arguments:
1775: + coarse - coarser DM to use as a base
1776: . interp - interpolation matrix, apply using MatInterpolate()
1777: - fine - finer DM to update
1779: Level: developer
1781: .seealso: DMRefineHookAdd(), MatInterpolate()
1782: @*/
1783: PetscErrorCode DMInterpolate(DM coarse,Mat interp,DM fine)
1784: {
1785: PetscErrorCode ierr;
1786: DMRefineHookLink link;
1789: for (link=fine->refinehook; link; link=link->next) {
1790: if (link->interphook) {
1791: (*link->interphook)(coarse,interp,fine,link->ctx);
1792: }
1793: }
1794: return(0);
1795: }
1797: /*@
1798: DMGetRefineLevel - Get's the number of refinements that have generated this DM.
1800: Not Collective
1802: Input Parameter:
1803: . dm - the DM object
1805: Output Parameter:
1806: . level - number of refinements
1808: Level: developer
1810: .seealso DMCoarsen(), DMGetCoarsenLevel(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation()
1812: @*/
1813: PetscErrorCode DMGetRefineLevel(DM dm,PetscInt *level)
1814: {
1817: *level = dm->levelup;
1818: return(0);
1819: }
1821: /*@
1822: DMSetRefineLevel - Set's the number of refinements that have generated this DM.
1824: Not Collective
1826: Input Parameter:
1827: + dm - the DM object
1828: - level - number of refinements
1830: Level: advanced
1832: Notes: This value is used by PCMG to determine how many multigrid levels to use
1834: .seealso DMCoarsen(), DMGetCoarsenLevel(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation()
1836: @*/
1837: PetscErrorCode DMSetRefineLevel(DM dm,PetscInt level)
1838: {
1841: dm->levelup = level;
1842: return(0);
1843: }
1845: /*@C
1846: DMGlobalToLocalHookAdd - adds a callback to be run when global to local is called
1848: Logically Collective
1850: Input Arguments:
1851: + dm - the DM
1852: . beginhook - function to run at the beginning of DMGlobalToLocalBegin()
1853: . endhook - function to run after DMGlobalToLocalEnd() has completed
1854: - ctx - [optional] user-defined context for provide data for the hooks (may be NULL)
1856: Calling sequence for beginhook:
1857: $ beginhook(DM fine,VecScatter out,VecScatter in,DM coarse,void *ctx)
1859: + dm - global DM
1860: . g - global vector
1861: . mode - mode
1862: . l - local vector
1863: - ctx - optional user-defined function context
1866: Calling sequence for endhook:
1867: $ endhook(DM fine,VecScatter out,VecScatter in,DM coarse,void *ctx)
1869: + global - global DM
1870: - ctx - optional user-defined function context
1872: Level: advanced
1874: .seealso: DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
1875: @*/
1876: PetscErrorCode DMGlobalToLocalHookAdd(DM dm,PetscErrorCode (*beginhook)(DM,Vec,InsertMode,Vec,void*),PetscErrorCode (*endhook)(DM,Vec,InsertMode,Vec,void*),void *ctx)
1877: {
1878: PetscErrorCode ierr;
1879: DMGlobalToLocalHookLink link,*p;
1883: for (p=&dm->gtolhook; *p; p=&(*p)->next) {} /* Scan to the end of the current list of hooks */
1884: PetscNew(&link);
1885: link->beginhook = beginhook;
1886: link->endhook = endhook;
1887: link->ctx = ctx;
1888: link->next = NULL;
1889: *p = link;
1890: return(0);
1891: }
1893: static PetscErrorCode DMGlobalToLocalHook_Constraints(DM dm, Vec g, InsertMode mode, Vec l, void *ctx)
1894: {
1895: Mat cMat;
1896: Vec cVec;
1897: PetscSection section, cSec;
1898: PetscInt pStart, pEnd, p, dof;
1903: DMGetDefaultConstraints(dm,&cSec,&cMat);
1904: if (cMat && (mode == INSERT_VALUES || mode == INSERT_ALL_VALUES || mode == INSERT_BC_VALUES)) {
1905: PetscInt nRows;
1907: MatGetSize(cMat,&nRows,NULL);
1908: if (nRows <= 0) return(0);
1909: DMGetDefaultSection(dm,§ion);
1910: MatCreateVecs(cMat,NULL,&cVec);
1911: MatMult(cMat,l,cVec);
1912: PetscSectionGetChart(cSec,&pStart,&pEnd);
1913: for (p = pStart; p < pEnd; p++) {
1914: PetscSectionGetDof(cSec,p,&dof);
1915: if (dof) {
1916: PetscScalar *vals;
1917: VecGetValuesSection(cVec,cSec,p,&vals);
1918: VecSetValuesSection(l,section,p,vals,INSERT_ALL_VALUES);
1919: }
1920: }
1921: VecDestroy(&cVec);
1922: }
1923: return(0);
1924: }
1926: /*@
1927: DMGlobalToLocalBegin - Begins updating local vectors from global vector
1929: Neighbor-wise Collective on DM
1931: Input Parameters:
1932: + dm - the DM object
1933: . g - the global vector
1934: . mode - INSERT_VALUES or ADD_VALUES
1935: - l - the local vector
1938: Level: beginner
1940: .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMGlobalToLocalEnd(), DMLocalToGlobalBegin()
1942: @*/
1943: PetscErrorCode DMGlobalToLocalBegin(DM dm,Vec g,InsertMode mode,Vec l)
1944: {
1945: PetscSF sf;
1946: PetscErrorCode ierr;
1947: DMGlobalToLocalHookLink link;
1951: for (link=dm->gtolhook; link; link=link->next) {
1952: if (link->beginhook) {
1953: (*link->beginhook)(dm,g,mode,l,link->ctx);
1954: }
1955: }
1956: DMGetDefaultSF(dm, &sf);
1957: if (sf) {
1958: const PetscScalar *gArray;
1959: PetscScalar *lArray;
1961: if (mode == ADD_VALUES) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode);
1962: VecGetArray(l, &lArray);
1963: VecGetArrayRead(g, &gArray);
1964: PetscSFBcastBegin(sf, MPIU_SCALAR, gArray, lArray);
1965: VecRestoreArray(l, &lArray);
1966: VecRestoreArrayRead(g, &gArray);
1967: } else {
1968: (*dm->ops->globaltolocalbegin)(dm,g,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),l);
1969: }
1970: return(0);
1971: }
1973: /*@
1974: DMGlobalToLocalEnd - Ends updating local vectors from global vector
1976: Neighbor-wise Collective on DM
1978: Input Parameters:
1979: + dm - the DM object
1980: . g - the global vector
1981: . mode - INSERT_VALUES or ADD_VALUES
1982: - l - the local vector
1985: Level: beginner
1987: .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMGlobalToLocalEnd(), DMLocalToGlobalBegin()
1989: @*/
1990: PetscErrorCode DMGlobalToLocalEnd(DM dm,Vec g,InsertMode mode,Vec l)
1991: {
1992: PetscSF sf;
1993: PetscErrorCode ierr;
1994: const PetscScalar *gArray;
1995: PetscScalar *lArray;
1996: DMGlobalToLocalHookLink link;
2000: DMGetDefaultSF(dm, &sf);
2001: if (sf) {
2002: if (mode == ADD_VALUES) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode);
2004: VecGetArray(l, &lArray);
2005: VecGetArrayRead(g, &gArray);
2006: PetscSFBcastEnd(sf, MPIU_SCALAR, gArray, lArray);
2007: VecRestoreArray(l, &lArray);
2008: VecRestoreArrayRead(g, &gArray);
2009: } else {
2010: (*dm->ops->globaltolocalend)(dm,g,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),l);
2011: }
2012: DMGlobalToLocalHook_Constraints(dm,g,mode,l,NULL);
2013: for (link=dm->gtolhook; link; link=link->next) {
2014: if (link->endhook) {(*link->endhook)(dm,g,mode,l,link->ctx);}
2015: }
2016: return(0);
2017: }
2019: /*@C
2020: DMLocalToGlobalHookAdd - adds a callback to be run when a local to global is called
2022: Logically Collective
2024: Input Arguments:
2025: + dm - the DM
2026: . beginhook - function to run at the beginning of DMLocalToGlobalBegin()
2027: . endhook - function to run after DMLocalToGlobalEnd() has completed
2028: - ctx - [optional] user-defined context for provide data for the hooks (may be NULL)
2030: Calling sequence for beginhook:
2031: $ beginhook(DM fine,Vec l,InsertMode mode,Vec g,void *ctx)
2033: + dm - global DM
2034: . l - local vector
2035: . mode - mode
2036: . g - global vector
2037: - ctx - optional user-defined function context
2040: Calling sequence for endhook:
2041: $ endhook(DM fine,Vec l,InsertMode mode,Vec g,void *ctx)
2043: + global - global DM
2044: . l - local vector
2045: . mode - mode
2046: . g - global vector
2047: - ctx - optional user-defined function context
2049: Level: advanced
2051: .seealso: DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2052: @*/
2053: PetscErrorCode DMLocalToGlobalHookAdd(DM dm,PetscErrorCode (*beginhook)(DM,Vec,InsertMode,Vec,void*),PetscErrorCode (*endhook)(DM,Vec,InsertMode,Vec,void*),void *ctx)
2054: {
2055: PetscErrorCode ierr;
2056: DMLocalToGlobalHookLink link,*p;
2060: for (p=&dm->ltoghook; *p; p=&(*p)->next) {} /* Scan to the end of the current list of hooks */
2061: PetscNew(&link);
2062: link->beginhook = beginhook;
2063: link->endhook = endhook;
2064: link->ctx = ctx;
2065: link->next = NULL;
2066: *p = link;
2067: return(0);
2068: }
2070: static PetscErrorCode DMLocalToGlobalHook_Constraints(DM dm, Vec l, InsertMode mode, Vec g, void *ctx)
2071: {
2072: Mat cMat;
2073: Vec cVec;
2074: PetscSection section, cSec;
2075: PetscInt pStart, pEnd, p, dof;
2080: DMGetDefaultConstraints(dm,&cSec,&cMat);
2081: if (cMat && (mode == ADD_VALUES || mode == ADD_ALL_VALUES || mode == ADD_BC_VALUES)) {
2082: PetscInt nRows;
2084: MatGetSize(cMat,&nRows,NULL);
2085: if (nRows <= 0) return(0);
2086: DMGetDefaultSection(dm,§ion);
2087: MatCreateVecs(cMat,NULL,&cVec);
2088: PetscSectionGetChart(cSec,&pStart,&pEnd);
2089: for (p = pStart; p < pEnd; p++) {
2090: PetscSectionGetDof(cSec,p,&dof);
2091: if (dof) {
2092: PetscInt d;
2093: PetscScalar *vals;
2094: VecGetValuesSection(l,section,p,&vals);
2095: VecSetValuesSection(cVec,cSec,p,vals,mode);
2096: /* for this to be the true transpose, we have to zero the values that
2097: * we just extracted */
2098: for (d = 0; d < dof; d++) {
2099: vals[d] = 0.;
2100: }
2101: }
2102: }
2103: MatMultTransposeAdd(cMat,cVec,l,l);
2104: VecDestroy(&cVec);
2105: }
2106: return(0);
2107: }
2109: /*@
2110: DMLocalToGlobalBegin - updates global vectors from local vectors
2112: Neighbor-wise Collective on DM
2114: Input Parameters:
2115: + dm - the DM object
2116: . l - the local vector
2117: . 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.
2118: - g - the global vector
2120: Notes: In the ADD_VALUES case you normally would zero the receiving vector before beginning this operation.
2121: INSERT_VALUES is not supported for DMDA, in that case simply compute the values directly into a global vector instead of a local one.
2123: Level: beginner
2125: .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMGlobalToLocalEnd(), DMGlobalToLocalBegin()
2127: @*/
2128: PetscErrorCode DMLocalToGlobalBegin(DM dm,Vec l,InsertMode mode,Vec g)
2129: {
2130: PetscSF sf;
2131: PetscSection s, gs;
2132: DMLocalToGlobalHookLink link;
2133: const PetscScalar *lArray;
2134: PetscScalar *gArray;
2135: PetscBool isInsert;
2136: PetscErrorCode ierr;
2140: for (link=dm->ltoghook; link; link=link->next) {
2141: if (link->beginhook) {
2142: (*link->beginhook)(dm,l,mode,g,link->ctx);
2143: }
2144: }
2145: DMLocalToGlobalHook_Constraints(dm,l,mode,g,NULL);
2146: DMGetDefaultSF(dm, &sf);
2147: DMGetDefaultSection(dm, &s);
2148: switch (mode) {
2149: case INSERT_VALUES:
2150: case INSERT_ALL_VALUES:
2151: case INSERT_BC_VALUES:
2152: isInsert = PETSC_TRUE; break;
2153: case ADD_VALUES:
2154: case ADD_ALL_VALUES:
2155: case ADD_BC_VALUES:
2156: isInsert = PETSC_FALSE; break;
2157: default:
2158: SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode);
2159: }
2160: if (sf && !isInsert) {
2161: VecGetArrayRead(l, &lArray);
2162: VecGetArray(g, &gArray);
2163: PetscSFReduceBegin(sf, MPIU_SCALAR, lArray, gArray, MPIU_SUM);
2164: VecRestoreArrayRead(l, &lArray);
2165: VecRestoreArray(g, &gArray);
2166: } else if (s && isInsert) {
2167: PetscInt gStart, pStart, pEnd, p;
2169: DMGetDefaultGlobalSection(dm, &gs);
2170: PetscSectionGetChart(s, &pStart, &pEnd);
2171: VecGetOwnershipRange(g, &gStart, NULL);
2172: VecGetArrayRead(l, &lArray);
2173: VecGetArray(g, &gArray);
2174: for (p = pStart; p < pEnd; ++p) {
2175: PetscInt dof, gdof, cdof, gcdof, off, goff, d, e;
2177: PetscSectionGetDof(s, p, &dof);
2178: PetscSectionGetDof(gs, p, &gdof);
2179: PetscSectionGetConstraintDof(s, p, &cdof);
2180: PetscSectionGetConstraintDof(gs, p, &gcdof);
2181: PetscSectionGetOffset(s, p, &off);
2182: PetscSectionGetOffset(gs, p, &goff);
2183: /* Ignore off-process data and points with no global data */
2184: if (!gdof || goff < 0) continue;
2185: 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);
2186: /* If no constraints are enforced in the global vector */
2187: if (!gcdof) {
2188: for (d = 0; d < dof; ++d) gArray[goff-gStart+d] = lArray[off+d];
2189: /* If constraints are enforced in the global vector */
2190: } else if (cdof == gcdof) {
2191: const PetscInt *cdofs;
2192: PetscInt cind = 0;
2194: PetscSectionGetConstraintIndices(s, p, &cdofs);
2195: for (d = 0, e = 0; d < dof; ++d) {
2196: if ((cind < cdof) && (d == cdofs[cind])) {++cind; continue;}
2197: gArray[goff-gStart+e++] = lArray[off+d];
2198: }
2199: } 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);
2200: }
2201: VecRestoreArrayRead(l, &lArray);
2202: VecRestoreArray(g, &gArray);
2203: } else {
2204: (*dm->ops->localtoglobalbegin)(dm,l,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),g);
2205: }
2206: return(0);
2207: }
2209: /*@
2210: DMLocalToGlobalEnd - updates global vectors from local vectors
2212: Neighbor-wise Collective on DM
2214: Input Parameters:
2215: + dm - the DM object
2216: . l - the local vector
2217: . mode - INSERT_VALUES or ADD_VALUES
2218: - g - the global vector
2221: Level: beginner
2223: .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMGlobalToLocalEnd(), DMGlobalToLocalEnd()
2225: @*/
2226: PetscErrorCode DMLocalToGlobalEnd(DM dm,Vec l,InsertMode mode,Vec g)
2227: {
2228: PetscSF sf;
2229: PetscSection s;
2230: DMLocalToGlobalHookLink link;
2231: PetscBool isInsert;
2232: PetscErrorCode ierr;
2236: DMGetDefaultSF(dm, &sf);
2237: DMGetDefaultSection(dm, &s);
2238: switch (mode) {
2239: case INSERT_VALUES:
2240: case INSERT_ALL_VALUES:
2241: isInsert = PETSC_TRUE; break;
2242: case ADD_VALUES:
2243: case ADD_ALL_VALUES:
2244: isInsert = PETSC_FALSE; break;
2245: default:
2246: SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode);
2247: }
2248: if (sf && !isInsert) {
2249: const PetscScalar *lArray;
2250: PetscScalar *gArray;
2252: VecGetArrayRead(l, &lArray);
2253: VecGetArray(g, &gArray);
2254: PetscSFReduceEnd(sf, MPIU_SCALAR, lArray, gArray, MPIU_SUM);
2255: VecRestoreArrayRead(l, &lArray);
2256: VecRestoreArray(g, &gArray);
2257: } else if (s && isInsert) {
2258: } else {
2259: (*dm->ops->localtoglobalend)(dm,l,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),g);
2260: }
2261: for (link=dm->ltoghook; link; link=link->next) {
2262: if (link->endhook) {(*link->endhook)(dm,g,mode,l,link->ctx);}
2263: }
2264: return(0);
2265: }
2267: /*@
2268: DMLocalToLocalBegin - Maps from a local vector (including ghost points
2269: that contain irrelevant values) to another local vector where the ghost
2270: points in the second are set correctly. Must be followed by DMLocalToLocalEnd().
2272: Neighbor-wise Collective on DM and Vec
2274: Input Parameters:
2275: + dm - the DM object
2276: . g - the original local vector
2277: - mode - one of INSERT_VALUES or ADD_VALUES
2279: Output Parameter:
2280: . l - the local vector with correct ghost values
2282: Level: intermediate
2284: Notes:
2285: The local vectors used here need not be the same as those
2286: obtained from DMCreateLocalVector(), BUT they
2287: must have the same parallel data layout; they could, for example, be
2288: obtained with VecDuplicate() from the DM originating vectors.
2290: .keywords: DM, local-to-local, begin
2291: .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateLocalVector(), DMCreateGlobalVector(), DMCreateInterpolation(), DMLocalToLocalEnd(), DMGlobalToLocalEnd(), DMLocalToGlobalBegin()
2293: @*/
2294: PetscErrorCode DMLocalToLocalBegin(DM dm,Vec g,InsertMode mode,Vec l)
2295: {
2296: PetscErrorCode ierr;
2300: if (!dm->ops->localtolocalbegin) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"This DM does not support local to local maps");
2301: (*dm->ops->localtolocalbegin)(dm,g,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),l);
2302: return(0);
2303: }
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: if (!dm->ops->localtolocalend) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"This DM does not support local to local maps");
2339: (*dm->ops->localtolocalend)(dm,g,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),l);
2340: return(0);
2341: }
2344: /*@
2345: DMCoarsen - Coarsens a DM object
2347: Collective on DM
2349: Input Parameter:
2350: + dm - the DM object
2351: - comm - the communicator to contain the new DM object (or MPI_COMM_NULL)
2353: Output Parameter:
2354: . dmc - the coarsened DM
2356: Level: developer
2358: .seealso DMRefine(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation()
2360: @*/
2361: PetscErrorCode DMCoarsen(DM dm, MPI_Comm comm, DM *dmc)
2362: {
2363: PetscErrorCode ierr;
2364: DMCoarsenHookLink link;
2368: if (!dm->ops->coarsen) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"This DM cannot coarsen");
2369: PetscLogEventBegin(DM_Coarsen,dm,0,0,0);
2370: (*dm->ops->coarsen)(dm, comm, dmc);
2371: if (!(*dmc)) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "NULL coarse mesh produced");
2372: DMSetCoarseDM(dm,*dmc);
2373: (*dmc)->ops->creatematrix = dm->ops->creatematrix;
2374: PetscObjectCopyFortranFunctionPointers((PetscObject)dm,(PetscObject)*dmc);
2375: (*dmc)->ctx = dm->ctx;
2376: (*dmc)->levelup = dm->levelup;
2377: (*dmc)->leveldown = dm->leveldown + 1;
2378: DMSetMatType(*dmc,dm->mattype);
2379: for (link=dm->coarsenhook; link; link=link->next) {
2380: if (link->coarsenhook) {(*link->coarsenhook)(dm,*dmc,link->ctx);}
2381: }
2382: PetscLogEventEnd(DM_Coarsen,dm,0,0,0);
2383: return(0);
2384: }
2386: /*@C
2387: DMCoarsenHookAdd - adds a callback to be run when restricting a nonlinear problem to the coarse grid
2389: Logically Collective
2391: Input Arguments:
2392: + fine - nonlinear solver context on which to run a hook when restricting to a coarser level
2393: . coarsenhook - function to run when setting up a coarser level
2394: . restricthook - function to run to update data on coarser levels (once per SNESSolve())
2395: - ctx - [optional] user-defined context for provide data for the hooks (may be NULL)
2397: Calling sequence of coarsenhook:
2398: $ coarsenhook(DM fine,DM coarse,void *ctx);
2400: + fine - fine level DM
2401: . coarse - coarse level DM to restrict problem to
2402: - ctx - optional user-defined function context
2404: Calling sequence for restricthook:
2405: $ restricthook(DM fine,Mat mrestrict,Vec rscale,Mat inject,DM coarse,void *ctx)
2407: + fine - fine level DM
2408: . mrestrict - matrix restricting a fine-level solution to the coarse grid
2409: . rscale - scaling vector for restriction
2410: . inject - matrix restricting by injection
2411: . coarse - coarse level DM to update
2412: - ctx - optional user-defined function context
2414: Level: advanced
2416: Notes:
2417: This function is only needed if auxiliary data needs to be set up on coarse grids.
2419: If this function is called multiple times, the hooks will be run in the order they are added.
2421: In order to compose with nonlinear preconditioning without duplicating storage, the hook should be implemented to
2422: extract the finest level information from its context (instead of from the SNES).
2424: This function is currently not available from Fortran.
2426: .seealso: DMCoarsenHookRemove(), DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2427: @*/
2428: PetscErrorCode DMCoarsenHookAdd(DM fine,PetscErrorCode (*coarsenhook)(DM,DM,void*),PetscErrorCode (*restricthook)(DM,Mat,Vec,Mat,DM,void*),void *ctx)
2429: {
2430: PetscErrorCode ierr;
2431: DMCoarsenHookLink link,*p;
2435: for (p=&fine->coarsenhook; *p; p=&(*p)->next) { /* Scan to the end of the current list of hooks */
2436: if ((*p)->coarsenhook == coarsenhook && (*p)->restricthook == restricthook && (*p)->ctx == ctx) return(0);
2437: }
2438: PetscNew(&link);
2439: link->coarsenhook = coarsenhook;
2440: link->restricthook = restricthook;
2441: link->ctx = ctx;
2442: link->next = NULL;
2443: *p = link;
2444: return(0);
2445: }
2447: /*@C
2448: DMCoarsenHookRemove - remove a callback from the list of hooks to be run when restricting a nonlinear problem to the coarse grid
2450: Logically Collective
2452: Input Arguments:
2453: + fine - nonlinear solver context on which to run a hook when restricting to a coarser level
2454: . coarsenhook - function to run when setting up a coarser level
2455: . restricthook - function to run to update data on coarser levels (once per SNESSolve())
2456: - ctx - [optional] user-defined context for provide data for the hooks (may be NULL)
2458: Level: advanced
2460: Notes:
2461: This function does nothing if the hook is not in the list.
2463: This function is currently not available from Fortran.
2465: .seealso: DMCoarsenHookAdd(), DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2466: @*/
2467: PetscErrorCode DMCoarsenHookRemove(DM fine,PetscErrorCode (*coarsenhook)(DM,DM,void*),PetscErrorCode (*restricthook)(DM,Mat,Vec,Mat,DM,void*),void *ctx)
2468: {
2469: PetscErrorCode ierr;
2470: DMCoarsenHookLink link,*p;
2474: for (p=&fine->coarsenhook; *p; p=&(*p)->next) { /* Search the list of current hooks */
2475: if ((*p)->coarsenhook == coarsenhook && (*p)->restricthook == restricthook && (*p)->ctx == ctx) {
2476: link = *p;
2477: *p = link->next;
2478: PetscFree(link);
2479: break;
2480: }
2481: }
2482: return(0);
2483: }
2486: /*@
2487: DMRestrict - restricts user-defined problem data to a coarser DM by running hooks registered by DMCoarsenHookAdd()
2489: Collective if any hooks are
2491: Input Arguments:
2492: + fine - finer DM to use as a base
2493: . restrct - restriction matrix, apply using MatRestrict()
2494: . rscale - scaling vector for restriction
2495: . inject - injection matrix, also use MatRestrict()
2496: - coarse - coarser DM to update
2498: Level: developer
2500: .seealso: DMCoarsenHookAdd(), MatRestrict()
2501: @*/
2502: PetscErrorCode DMRestrict(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse)
2503: {
2504: PetscErrorCode ierr;
2505: DMCoarsenHookLink link;
2508: for (link=fine->coarsenhook; link; link=link->next) {
2509: if (link->restricthook) {
2510: (*link->restricthook)(fine,restrct,rscale,inject,coarse,link->ctx);
2511: }
2512: }
2513: return(0);
2514: }
2516: /*@C
2517: DMSubDomainHookAdd - adds a callback to be run when restricting a problem to the coarse grid
2519: Logically Collective
2521: Input Arguments:
2522: + global - global DM
2523: . ddhook - function to run to pass data to the decomposition DM upon its creation
2524: . restricthook - function to run to update data on block solve (at the beginning of the block solve)
2525: - ctx - [optional] user-defined context for provide data for the hooks (may be NULL)
2528: Calling sequence for ddhook:
2529: $ ddhook(DM global,DM block,void *ctx)
2531: + global - global DM
2532: . block - block DM
2533: - ctx - optional user-defined function context
2535: Calling sequence for restricthook:
2536: $ restricthook(DM global,VecScatter out,VecScatter in,DM block,void *ctx)
2538: + global - global DM
2539: . out - scatter to the outer (with ghost and overlap points) block vector
2540: . in - scatter to block vector values only owned locally
2541: . block - block DM
2542: - ctx - optional user-defined function context
2544: Level: advanced
2546: Notes:
2547: This function is only needed if auxiliary data needs to be set up on subdomain DMs.
2549: If this function is called multiple times, the hooks will be run in the order they are added.
2551: In order to compose with nonlinear preconditioning without duplicating storage, the hook should be implemented to
2552: extract the global information from its context (instead of from the SNES).
2554: This function is currently not available from Fortran.
2556: .seealso: DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2557: @*/
2558: PetscErrorCode DMSubDomainHookAdd(DM global,PetscErrorCode (*ddhook)(DM,DM,void*),PetscErrorCode (*restricthook)(DM,VecScatter,VecScatter,DM,void*),void *ctx)
2559: {
2560: PetscErrorCode ierr;
2561: DMSubDomainHookLink link,*p;
2565: for (p=&global->subdomainhook; *p; p=&(*p)->next) { /* Scan to the end of the current list of hooks */
2566: if ((*p)->ddhook == ddhook && (*p)->restricthook == restricthook && (*p)->ctx == ctx) return(0);
2567: }
2568: PetscNew(&link);
2569: link->restricthook = restricthook;
2570: link->ddhook = ddhook;
2571: link->ctx = ctx;
2572: link->next = NULL;
2573: *p = link;
2574: return(0);
2575: }
2577: /*@C
2578: DMSubDomainHookRemove - remove a callback from the list to be run when restricting a problem to the coarse grid
2580: Logically Collective
2582: Input Arguments:
2583: + global - global DM
2584: . ddhook - function to run to pass data to the decomposition DM upon its creation
2585: . restricthook - function to run to update data on block solve (at the beginning of the block solve)
2586: - ctx - [optional] user-defined context for provide data for the hooks (may be NULL)
2588: Level: advanced
2590: Notes:
2592: This function is currently not available from Fortran.
2594: .seealso: DMSubDomainHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2595: @*/
2596: PetscErrorCode DMSubDomainHookRemove(DM global,PetscErrorCode (*ddhook)(DM,DM,void*),PetscErrorCode (*restricthook)(DM,VecScatter,VecScatter,DM,void*),void *ctx)
2597: {
2598: PetscErrorCode ierr;
2599: DMSubDomainHookLink link,*p;
2603: for (p=&global->subdomainhook; *p; p=&(*p)->next) { /* Search the list of current hooks */
2604: if ((*p)->ddhook == ddhook && (*p)->restricthook == restricthook && (*p)->ctx == ctx) {
2605: link = *p;
2606: *p = link->next;
2607: PetscFree(link);
2608: break;
2609: }
2610: }
2611: return(0);
2612: }
2614: /*@
2615: DMSubDomainRestrict - restricts user-defined problem data to a block DM by running hooks registered by DMSubDomainHookAdd()
2617: Collective if any hooks are
2619: Input Arguments:
2620: + fine - finer DM to use as a base
2621: . oscatter - scatter from domain global vector filling subdomain global vector with overlap
2622: . gscatter - scatter from domain global vector filling subdomain local vector with ghosts
2623: - coarse - coarer DM to update
2625: Level: developer
2627: .seealso: DMCoarsenHookAdd(), MatRestrict()
2628: @*/
2629: PetscErrorCode DMSubDomainRestrict(DM global,VecScatter oscatter,VecScatter gscatter,DM subdm)
2630: {
2631: PetscErrorCode ierr;
2632: DMSubDomainHookLink link;
2635: for (link=global->subdomainhook; link; link=link->next) {
2636: if (link->restricthook) {
2637: (*link->restricthook)(global,oscatter,gscatter,subdm,link->ctx);
2638: }
2639: }
2640: return(0);
2641: }
2643: /*@
2644: DMGetCoarsenLevel - Get's the number of coarsenings that have generated this DM.
2646: Not Collective
2648: Input Parameter:
2649: . dm - the DM object
2651: Output Parameter:
2652: . level - number of coarsenings
2654: Level: developer
2656: .seealso DMCoarsen(), DMGetRefineLevel(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation()
2658: @*/
2659: PetscErrorCode DMGetCoarsenLevel(DM dm,PetscInt *level)
2660: {
2663: *level = dm->leveldown;
2664: return(0);
2665: }
2669: /*@C
2670: DMRefineHierarchy - Refines a DM object, all levels at once
2672: Collective on DM
2674: Input Parameter:
2675: + dm - the DM object
2676: - nlevels - the number of levels of refinement
2678: Output Parameter:
2679: . dmf - the refined DM hierarchy
2681: Level: developer
2683: .seealso DMCoarsenHierarchy(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation()
2685: @*/
2686: PetscErrorCode DMRefineHierarchy(DM dm,PetscInt nlevels,DM dmf[])
2687: {
2692: if (nlevels < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"nlevels cannot be negative");
2693: if (nlevels == 0) return(0);
2694: if (dm->ops->refinehierarchy) {
2695: (*dm->ops->refinehierarchy)(dm,nlevels,dmf);
2696: } else if (dm->ops->refine) {
2697: PetscInt i;
2699: DMRefine(dm,PetscObjectComm((PetscObject)dm),&dmf[0]);
2700: for (i=1; i<nlevels; i++) {
2701: DMRefine(dmf[i-1],PetscObjectComm((PetscObject)dm),&dmf[i]);
2702: }
2703: } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"No RefineHierarchy for this DM yet");
2704: return(0);
2705: }
2707: /*@C
2708: DMCoarsenHierarchy - Coarsens a DM object, all levels at once
2710: Collective on DM
2712: Input Parameter:
2713: + dm - the DM object
2714: - nlevels - the number of levels of coarsening
2716: Output Parameter:
2717: . dmc - the coarsened DM hierarchy
2719: Level: developer
2721: .seealso DMRefineHierarchy(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation()
2723: @*/
2724: PetscErrorCode DMCoarsenHierarchy(DM dm, PetscInt nlevels, DM dmc[])
2725: {
2730: if (nlevels < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"nlevels cannot be negative");
2731: if (nlevels == 0) return(0);
2733: if (dm->ops->coarsenhierarchy) {
2734: (*dm->ops->coarsenhierarchy)(dm, nlevels, dmc);
2735: } else if (dm->ops->coarsen) {
2736: PetscInt i;
2738: DMCoarsen(dm,PetscObjectComm((PetscObject)dm),&dmc[0]);
2739: for (i=1; i<nlevels; i++) {
2740: DMCoarsen(dmc[i-1],PetscObjectComm((PetscObject)dm),&dmc[i]);
2741: }
2742: } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"No CoarsenHierarchy for this DM yet");
2743: return(0);
2744: }
2746: /*@
2747: DMCreateAggregates - Gets the aggregates that map between
2748: grids associated with two DMs.
2750: Collective on DM
2752: Input Parameters:
2753: + dmc - the coarse grid DM
2754: - dmf - the fine grid DM
2756: Output Parameters:
2757: . rest - the restriction matrix (transpose of the projection matrix)
2759: Level: intermediate
2761: .keywords: interpolation, restriction, multigrid
2763: .seealso: DMRefine(), DMCreateInjection(), DMCreateInterpolation()
2764: @*/
2765: PetscErrorCode DMCreateAggregates(DM dmc, DM dmf, Mat *rest)
2766: {
2772: (*dmc->ops->getaggregates)(dmc, dmf, rest);
2773: return(0);
2774: }
2776: /*@C
2777: DMSetApplicationContextDestroy - Sets a user function that will be called to destroy the application context when the DM is destroyed
2779: Not Collective
2781: Input Parameters:
2782: + dm - the DM object
2783: - destroy - the destroy function
2785: Level: intermediate
2787: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext()
2789: @*/
2790: PetscErrorCode DMSetApplicationContextDestroy(DM dm,PetscErrorCode (*destroy)(void**))
2791: {
2794: dm->ctxdestroy = destroy;
2795: return(0);
2796: }
2798: /*@
2799: DMSetApplicationContext - Set a user context into a DM object
2801: Not Collective
2803: Input Parameters:
2804: + dm - the DM object
2805: - ctx - the user context
2807: Level: intermediate
2809: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext()
2811: @*/
2812: PetscErrorCode DMSetApplicationContext(DM dm,void *ctx)
2813: {
2816: dm->ctx = ctx;
2817: return(0);
2818: }
2820: /*@
2821: DMGetApplicationContext - Gets a user context from a DM object
2823: Not Collective
2825: Input Parameter:
2826: . dm - the DM object
2828: Output Parameter:
2829: . ctx - the user context
2831: Level: intermediate
2833: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext()
2835: @*/
2836: PetscErrorCode DMGetApplicationContext(DM dm,void *ctx)
2837: {
2840: *(void**)ctx = dm->ctx;
2841: return(0);
2842: }
2844: /*@C
2845: DMSetVariableBounds - sets a function to compute the lower and upper bound vectors for SNESVI.
2847: Logically Collective on DM
2849: Input Parameter:
2850: + dm - the DM object
2851: - f - the function that computes variable bounds used by SNESVI (use NULL to cancel a previous function that was set)
2853: Level: intermediate
2855: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext(),
2856: DMSetJacobian()
2858: @*/
2859: PetscErrorCode DMSetVariableBounds(DM dm,PetscErrorCode (*f)(DM,Vec,Vec))
2860: {
2862: dm->ops->computevariablebounds = f;
2863: return(0);
2864: }
2866: /*@
2867: DMHasVariableBounds - does the DM object have a variable bounds function?
2869: Not Collective
2871: Input Parameter:
2872: . dm - the DM object to destroy
2874: Output Parameter:
2875: . flg - PETSC_TRUE if the variable bounds function exists
2877: Level: developer
2879: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext()
2881: @*/
2882: PetscErrorCode DMHasVariableBounds(DM dm,PetscBool *flg)
2883: {
2885: *flg = (dm->ops->computevariablebounds) ? PETSC_TRUE : PETSC_FALSE;
2886: return(0);
2887: }
2889: /*@C
2890: DMComputeVariableBounds - compute variable bounds used by SNESVI.
2892: Logically Collective on DM
2894: Input Parameters:
2895: . dm - the DM object
2897: Output parameters:
2898: + xl - lower bound
2899: - xu - upper bound
2901: Level: advanced
2903: Notes: This is generally not called by users. It calls the function provided by the user with DMSetVariableBounds()
2905: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext()
2907: @*/
2908: PetscErrorCode DMComputeVariableBounds(DM dm, Vec xl, Vec xu)
2909: {
2915: if (dm->ops->computevariablebounds) {
2916: (*dm->ops->computevariablebounds)(dm, xl,xu);
2917: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "This DM is incapable of computing variable bounds.");
2918: return(0);
2919: }
2921: /*@
2922: DMHasColoring - does the DM object have a method of providing a coloring?
2924: Not Collective
2926: Input Parameter:
2927: . dm - the DM object
2929: Output Parameter:
2930: . flg - PETSC_TRUE if the DM has facilities for DMCreateColoring().
2932: Level: developer
2934: .seealso DMHasFunction(), DMCreateColoring()
2936: @*/
2937: PetscErrorCode DMHasColoring(DM dm,PetscBool *flg)
2938: {
2940: *flg = (dm->ops->getcoloring) ? PETSC_TRUE : PETSC_FALSE;
2941: return(0);
2942: }
2944: /*@
2945: DMHasCreateRestriction - does the DM object have a method of providing a restriction?
2947: Not Collective
2949: Input Parameter:
2950: . dm - the DM object
2952: Output Parameter:
2953: . flg - PETSC_TRUE if the DM has facilities for DMCreateRestriction().
2955: Level: developer
2957: .seealso DMHasFunction(), DMCreateRestriction()
2959: @*/
2960: PetscErrorCode DMHasCreateRestriction(DM dm,PetscBool *flg)
2961: {
2963: *flg = (dm->ops->createrestriction) ? PETSC_TRUE : PETSC_FALSE;
2964: return(0);
2965: }
2967: /*@C
2968: DMSetVec - set the vector at which to compute residual, Jacobian and VI bounds, if the problem is nonlinear.
2970: Collective on DM
2972: Input Parameter:
2973: + dm - the DM object
2974: - x - location to compute residual and Jacobian, if NULL is passed to those routines; will be NULL for linear problems.
2976: Level: developer
2978: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext()
2980: @*/
2981: PetscErrorCode DMSetVec(DM dm,Vec x)
2982: {
2986: if (x) {
2987: if (!dm->x) {
2988: DMCreateGlobalVector(dm,&dm->x);
2989: }
2990: VecCopy(x,dm->x);
2991: } else if (dm->x) {
2992: VecDestroy(&dm->x);
2993: }
2994: return(0);
2995: }
2997: PetscFunctionList DMList = NULL;
2998: PetscBool DMRegisterAllCalled = PETSC_FALSE;
3000: /*@C
3001: DMSetType - Builds a DM, for a particular DM implementation.
3003: Collective on DM
3005: Input Parameters:
3006: + dm - The DM object
3007: - method - The name of the DM type
3009: Options Database Key:
3010: . -dm_type <type> - Sets the DM type; use -help for a list of available types
3012: Notes:
3013: See "petsc/include/petscdm.h" for available DM types (for instance, DM1D, DM2D, or DM3D).
3015: Level: intermediate
3017: .keywords: DM, set, type
3018: .seealso: DMGetType(), DMCreate()
3019: @*/
3020: PetscErrorCode DMSetType(DM dm, DMType method)
3021: {
3022: PetscErrorCode (*r)(DM);
3023: PetscBool match;
3028: PetscObjectTypeCompare((PetscObject) dm, method, &match);
3029: if (match) return(0);
3031: DMRegisterAll();
3032: PetscFunctionListFind(DMList,method,&r);
3033: if (!r) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown DM type: %s", method);
3035: if (dm->ops->destroy) {
3036: (*dm->ops->destroy)(dm);
3037: dm->ops->destroy = NULL;
3038: }
3039: (*r)(dm);
3040: PetscObjectChangeTypeName((PetscObject)dm,method);
3041: return(0);
3042: }
3044: /*@C
3045: DMGetType - Gets the DM type name (as a string) from the DM.
3047: Not Collective
3049: Input Parameter:
3050: . dm - The DM
3052: Output Parameter:
3053: . type - The DM type name
3055: Level: intermediate
3057: .keywords: DM, get, type, name
3058: .seealso: DMSetType(), DMCreate()
3059: @*/
3060: PetscErrorCode DMGetType(DM dm, DMType *type)
3061: {
3067: DMRegisterAll();
3068: *type = ((PetscObject)dm)->type_name;
3069: return(0);
3070: }
3072: /*@C
3073: DMConvert - Converts a DM to another DM, either of the same or different type.
3075: Collective on DM
3077: Input Parameters:
3078: + dm - the DM
3079: - newtype - new DM type (use "same" for the same type)
3081: Output Parameter:
3082: . M - pointer to new DM
3084: Notes:
3085: Cannot be used to convert a sequential DM to parallel or parallel to sequential,
3086: the MPI communicator of the generated DM is always the same as the communicator
3087: of the input DM.
3089: Level: intermediate
3091: .seealso: DMCreate()
3092: @*/
3093: PetscErrorCode DMConvert(DM dm, DMType newtype, DM *M)
3094: {
3095: DM B;
3096: char convname[256];
3097: PetscBool sametype/*, issame */;
3104: PetscObjectTypeCompare((PetscObject) dm, newtype, &sametype);
3105: /* PetscStrcmp(newtype, "same", &issame); */
3106: if (sametype) {
3107: *M = dm;
3108: PetscObjectReference((PetscObject) dm);
3109: return(0);
3110: } else {
3111: PetscErrorCode (*conv)(DM, DMType, DM*) = NULL;
3113: /*
3114: Order of precedence:
3115: 1) See if a specialized converter is known to the current DM.
3116: 2) See if a specialized converter is known to the desired DM class.
3117: 3) See if a good general converter is registered for the desired class
3118: 4) See if a good general converter is known for the current matrix.
3119: 5) Use a really basic converter.
3120: */
3122: /* 1) See if a specialized converter is known to the current DM and the desired class */
3123: PetscStrcpy(convname,"DMConvert_");
3124: PetscStrcat(convname,((PetscObject) dm)->type_name);
3125: PetscStrcat(convname,"_");
3126: PetscStrcat(convname,newtype);
3127: PetscStrcat(convname,"_C");
3128: PetscObjectQueryFunction((PetscObject)dm,convname,&conv);
3129: if (conv) goto foundconv;
3131: /* 2) See if a specialized converter is known to the desired DM class. */
3132: DMCreate(PetscObjectComm((PetscObject)dm), &B);
3133: DMSetType(B, newtype);
3134: PetscStrcpy(convname,"DMConvert_");
3135: PetscStrcat(convname,((PetscObject) dm)->type_name);
3136: PetscStrcat(convname,"_");
3137: PetscStrcat(convname,newtype);
3138: PetscStrcat(convname,"_C");
3139: PetscObjectQueryFunction((PetscObject)B,convname,&conv);
3140: if (conv) {
3141: DMDestroy(&B);
3142: goto foundconv;
3143: }
3145: #if 0
3146: /* 3) See if a good general converter is registered for the desired class */
3147: conv = B->ops->convertfrom;
3148: DMDestroy(&B);
3149: if (conv) goto foundconv;
3151: /* 4) See if a good general converter is known for the current matrix */
3152: if (dm->ops->convert) {
3153: conv = dm->ops->convert;
3154: }
3155: if (conv) goto foundconv;
3156: #endif
3158: /* 5) Use a really basic converter. */
3159: SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "No conversion possible between DM types %s and %s", ((PetscObject) dm)->type_name, newtype);
3161: foundconv:
3162: PetscLogEventBegin(DM_Convert,dm,0,0,0);
3163: (*conv)(dm,newtype,M);
3164: /* Things that are independent of DM type: We should consult DMClone() here */
3165: {
3166: PetscBool isper;
3167: const PetscReal *maxCell, *L;
3168: const DMBoundaryType *bd;
3169: DMGetPeriodicity(dm, &isper, &maxCell, &L, &bd);
3170: DMSetPeriodicity(*M, isper, maxCell, L, bd);
3171: }
3172: PetscLogEventEnd(DM_Convert,dm,0,0,0);
3173: }
3174: PetscObjectStateIncrease((PetscObject) *M);
3175: return(0);
3176: }
3178: /*--------------------------------------------------------------------------------------------------------------------*/
3180: /*@C
3181: DMRegister - Adds a new DM component implementation
3183: Not Collective
3185: Input Parameters:
3186: + name - The name of a new user-defined creation routine
3187: - create_func - The creation routine itself
3189: Notes:
3190: DMRegister() may be called multiple times to add several user-defined DMs
3193: Sample usage:
3194: .vb
3195: DMRegister("my_da", MyDMCreate);
3196: .ve
3198: Then, your DM type can be chosen with the procedural interface via
3199: .vb
3200: DMCreate(MPI_Comm, DM *);
3201: DMSetType(DM,"my_da");
3202: .ve
3203: or at runtime via the option
3204: .vb
3205: -da_type my_da
3206: .ve
3208: Level: advanced
3210: .keywords: DM, register
3211: .seealso: DMRegisterAll(), DMRegisterDestroy()
3213: @*/
3214: PetscErrorCode DMRegister(const char sname[],PetscErrorCode (*function)(DM))
3215: {
3219: PetscFunctionListAdd(&DMList,sname,function);
3220: return(0);
3221: }
3223: /*@C
3224: DMLoad - Loads a DM that has been stored in binary with DMView().
3226: Collective on PetscViewer
3228: Input Parameters:
3229: + newdm - the newly loaded DM, this needs to have been created with DMCreate() or
3230: some related function before a call to DMLoad().
3231: - viewer - binary file viewer, obtained from PetscViewerBinaryOpen() or
3232: HDF5 file viewer, obtained from PetscViewerHDF5Open()
3234: Level: intermediate
3236: Notes:
3237: The type is determined by the data in the file, any type set into the DM before this call is ignored.
3239: Notes for advanced users:
3240: Most users should not need to know the details of the binary storage
3241: format, since DMLoad() and DMView() completely hide these details.
3242: But for anyone who's interested, the standard binary matrix storage
3243: format is
3244: .vb
3245: has not yet been determined
3246: .ve
3248: .seealso: PetscViewerBinaryOpen(), DMView(), MatLoad(), VecLoad()
3249: @*/
3250: PetscErrorCode DMLoad(DM newdm, PetscViewer viewer)
3251: {
3252: PetscBool isbinary, ishdf5;
3258: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
3259: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);
3260: if (isbinary) {
3261: PetscInt classid;
3262: char type[256];
3264: PetscViewerBinaryRead(viewer,&classid,1,NULL,PETSC_INT);
3265: if (classid != DM_FILE_CLASSID) SETERRQ1(PetscObjectComm((PetscObject)newdm),PETSC_ERR_ARG_WRONG,"Not DM next in file, classid found %d",(int)classid);
3266: PetscViewerBinaryRead(viewer,type,256,NULL,PETSC_CHAR);
3267: DMSetType(newdm, type);
3268: if (newdm->ops->load) {(*newdm->ops->load)(newdm,viewer);}
3269: } else if (ishdf5) {
3270: if (newdm->ops->load) {(*newdm->ops->load)(newdm,viewer);}
3271: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid viewer; open viewer with PetscViewerBinaryOpen() or PetscViewerHDF5Open()");
3272: return(0);
3273: }
3275: /******************************** FEM Support **********************************/
3277: PetscErrorCode DMPrintCellVector(PetscInt c, const char name[], PetscInt len, const PetscScalar x[])
3278: {
3279: PetscInt f;
3283: PetscPrintf(PETSC_COMM_SELF, "Cell %D Element %s\n", c, name);
3284: for (f = 0; f < len; ++f) {
3285: PetscPrintf(PETSC_COMM_SELF, " | %g |\n", (double)PetscRealPart(x[f]));
3286: }
3287: return(0);
3288: }
3290: PetscErrorCode DMPrintCellMatrix(PetscInt c, const char name[], PetscInt rows, PetscInt cols, const PetscScalar A[])
3291: {
3292: PetscInt f, g;
3296: PetscPrintf(PETSC_COMM_SELF, "Cell %D Element %s\n", c, name);
3297: for (f = 0; f < rows; ++f) {
3298: PetscPrintf(PETSC_COMM_SELF, " |");
3299: for (g = 0; g < cols; ++g) {
3300: PetscPrintf(PETSC_COMM_SELF, " % 9.5g", PetscRealPart(A[f*cols+g]));
3301: }
3302: PetscPrintf(PETSC_COMM_SELF, " |\n");
3303: }
3304: return(0);
3305: }
3307: PetscErrorCode DMPrintLocalVec(DM dm, const char name[], PetscReal tol, Vec X)
3308: {
3309: PetscInt localSize, bs;
3310: PetscMPIInt size;
3311: Vec x, xglob;
3312: const PetscScalar *xarray;
3313: PetscErrorCode ierr;
3316: MPI_Comm_size(PetscObjectComm((PetscObject) dm),&size);
3317: VecDuplicate(X, &x);
3318: VecCopy(X, x);
3319: VecChop(x, tol);
3320: PetscPrintf(PetscObjectComm((PetscObject) dm),"%s:\n",name);
3321: if (size > 1) {
3322: VecGetLocalSize(x,&localSize);
3323: VecGetArrayRead(x,&xarray);
3324: VecGetBlockSize(x,&bs);
3325: VecCreateMPIWithArray(PetscObjectComm((PetscObject) dm),bs,localSize,PETSC_DETERMINE,xarray,&xglob);
3326: } else {
3327: xglob = x;
3328: }
3329: VecView(xglob,PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject) dm)));
3330: if (size > 1) {
3331: VecDestroy(&xglob);
3332: VecRestoreArrayRead(x,&xarray);
3333: }
3334: VecDestroy(&x);
3335: return(0);
3336: }
3338: /*@
3339: DMGetDefaultSection - Get the PetscSection encoding the local data layout for the DM.
3341: Input Parameter:
3342: . dm - The DM
3344: Output Parameter:
3345: . section - The PetscSection
3347: Level: intermediate
3349: Note: This gets a borrowed reference, so the user should not destroy this PetscSection.
3351: .seealso: DMSetDefaultSection(), DMGetDefaultGlobalSection()
3352: @*/
3353: PetscErrorCode DMGetDefaultSection(DM dm, PetscSection *section)
3354: {
3360: if (!dm->defaultSection && dm->ops->createdefaultsection) {
3361: (*dm->ops->createdefaultsection)(dm);
3362: if (dm->defaultSection) {PetscObjectViewFromOptions((PetscObject) dm->defaultSection, NULL, "-dm_petscsection_view");}
3363: }
3364: *section = dm->defaultSection;
3365: return(0);
3366: }
3368: /*@
3369: DMSetDefaultSection - Set the PetscSection encoding the local data layout for the DM.
3371: Input Parameters:
3372: + dm - The DM
3373: - section - The PetscSection
3375: Level: intermediate
3377: Note: Any existing Section will be destroyed
3379: .seealso: DMSetDefaultSection(), DMGetDefaultGlobalSection()
3380: @*/
3381: PetscErrorCode DMSetDefaultSection(DM dm, PetscSection section)
3382: {
3383: PetscInt numFields = 0;
3384: PetscInt f;
3389: if (section) {
3391: PetscObjectReference((PetscObject)section);
3392: }
3393: PetscSectionDestroy(&dm->defaultSection);
3394: dm->defaultSection = section;
3395: if (section) {PetscSectionGetNumFields(dm->defaultSection, &numFields);}
3396: if (numFields) {
3397: DMSetNumFields(dm, numFields);
3398: for (f = 0; f < numFields; ++f) {
3399: PetscObject disc;
3400: const char *name;
3402: PetscSectionGetFieldName(dm->defaultSection, f, &name);
3403: DMGetField(dm, f, &disc);
3404: PetscObjectSetName(disc, name);
3405: }
3406: }
3407: /* The global section will be rebuilt in the next call to DMGetDefaultGlobalSection(). */
3408: PetscSectionDestroy(&dm->defaultGlobalSection);
3409: return(0);
3410: }
3412: /*@
3413: DMGetDefaultConstraints - Get the PetscSection and Mat the specify the local constraint interpolation. See DMSetDefaultConstraints() for a description of the purpose of constraint interpolation.
3415: not collective
3417: Input Parameter:
3418: . dm - The DM
3420: Output Parameter:
3421: + 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.
3422: - 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.
3424: Level: advanced
3426: Note: This gets borrowed references, so the user should not destroy the PetscSection or the Mat.
3428: .seealso: DMSetDefaultConstraints()
3429: @*/
3430: PetscErrorCode DMGetDefaultConstraints(DM dm, PetscSection *section, Mat *mat)
3431: {
3436: if (!dm->defaultConstraintSection && !dm->defaultConstraintMat && dm->ops->createdefaultconstraints) {(*dm->ops->createdefaultconstraints)(dm);}
3437: if (section) {*section = dm->defaultConstraintSection;}
3438: if (mat) {*mat = dm->defaultConstraintMat;}
3439: return(0);
3440: }
3442: /*@
3443: DMSetDefaultConstraints - Set the PetscSection and Mat the specify the local constraint interpolation.
3445: 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().
3447: 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.
3449: collective on dm
3451: Input Parameters:
3452: + dm - The DM
3453: + 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).
3454: - 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).
3456: Level: advanced
3458: Note: This increments the references of the PetscSection and the Mat, so they user can destroy them
3460: .seealso: DMGetDefaultConstraints()
3461: @*/
3462: PetscErrorCode DMSetDefaultConstraints(DM dm, PetscSection section, Mat mat)
3463: {
3464: PetscMPIInt result;
3469: if (section) {
3471: MPI_Comm_compare(PETSC_COMM_SELF,PetscObjectComm((PetscObject)section),&result);
3472: if (result != MPI_CONGRUENT && result != MPI_IDENT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMECOMM,"constraint section must have local communicator");
3473: }
3474: if (mat) {
3476: MPI_Comm_compare(PETSC_COMM_SELF,PetscObjectComm((PetscObject)mat),&result);
3477: if (result != MPI_CONGRUENT && result != MPI_IDENT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMECOMM,"constraint matrix must have local communicator");
3478: }
3479: PetscObjectReference((PetscObject)section);
3480: PetscSectionDestroy(&dm->defaultConstraintSection);
3481: dm->defaultConstraintSection = section;
3482: PetscObjectReference((PetscObject)mat);
3483: MatDestroy(&dm->defaultConstraintMat);
3484: dm->defaultConstraintMat = mat;
3485: return(0);
3486: }
3488: #ifdef PETSC_USE_DEBUG
3489: /*
3490: DMDefaultSectionCheckConsistency - Check the consistentcy of the global and local sections.
3492: Input Parameters:
3493: + dm - The DM
3494: . localSection - PetscSection describing the local data layout
3495: - globalSection - PetscSection describing the global data layout
3497: Level: intermediate
3499: .seealso: DMGetDefaultSF(), DMSetDefaultSF()
3500: */
3501: static PetscErrorCode DMDefaultSectionCheckConsistency_Internal(DM dm, PetscSection localSection, PetscSection globalSection)
3502: {
3503: MPI_Comm comm;
3504: PetscLayout layout;
3505: const PetscInt *ranges;
3506: PetscInt pStart, pEnd, p, nroots;
3507: PetscMPIInt size, rank;
3508: PetscBool valid = PETSC_TRUE, gvalid;
3509: PetscErrorCode ierr;
3512: PetscObjectGetComm((PetscObject)dm,&comm);
3514: MPI_Comm_size(comm, &size);
3515: MPI_Comm_rank(comm, &rank);
3516: PetscSectionGetChart(globalSection, &pStart, &pEnd);
3517: PetscSectionGetConstrainedStorageSize(globalSection, &nroots);
3518: PetscLayoutCreate(comm, &layout);
3519: PetscLayoutSetBlockSize(layout, 1);
3520: PetscLayoutSetLocalSize(layout, nroots);
3521: PetscLayoutSetUp(layout);
3522: PetscLayoutGetRanges(layout, &ranges);
3523: for (p = pStart; p < pEnd; ++p) {
3524: PetscInt dof, cdof, off, gdof, gcdof, goff, gsize, d;
3526: PetscSectionGetDof(localSection, p, &dof);
3527: PetscSectionGetOffset(localSection, p, &off);
3528: PetscSectionGetConstraintDof(localSection, p, &cdof);
3529: PetscSectionGetDof(globalSection, p, &gdof);
3530: PetscSectionGetConstraintDof(globalSection, p, &gcdof);
3531: PetscSectionGetOffset(globalSection, p, &goff);
3532: if (!gdof) continue; /* Censored point */
3533: 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;}
3534: 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;}
3535: if (gdof < 0) {
3536: gsize = gdof < 0 ? -(gdof+1)-gcdof : gdof-gcdof;
3537: for (d = 0; d < gsize; ++d) {
3538: PetscInt offset = -(goff+1) + d, r;
3540: PetscFindInt(offset,size+1,ranges,&r);
3541: if (r < 0) r = -(r+2);
3542: 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;}
3543: }
3544: }
3545: }
3546: PetscLayoutDestroy(&layout);
3547: PetscSynchronizedFlush(comm, NULL);
3548: MPIU_Allreduce(&valid, &gvalid, 1, MPIU_BOOL, MPI_LAND, comm);
3549: if (!gvalid) {
3550: DMView(dm, NULL);
3551: SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Inconsistent local and global sections");
3552: }
3553: return(0);
3554: }
3555: #endif
3557: /*@
3558: DMGetDefaultGlobalSection - Get the PetscSection encoding the global data layout for the DM.
3560: Collective on DM
3562: Input Parameter:
3563: . dm - The DM
3565: Output Parameter:
3566: . section - The PetscSection
3568: Level: intermediate
3570: Note: This gets a borrowed reference, so the user should not destroy this PetscSection.
3572: .seealso: DMSetDefaultSection(), DMGetDefaultSection()
3573: @*/
3574: PetscErrorCode DMGetDefaultGlobalSection(DM dm, PetscSection *section)
3575: {
3581: if (!dm->defaultGlobalSection) {
3582: PetscSection s;
3584: DMGetDefaultSection(dm, &s);
3585: if (!s) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONGSTATE, "DM must have a default PetscSection in order to create a global PetscSection");
3586: if (!dm->sf) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DM must have a default PetscSF in order to create a global PetscSection");
3587: PetscSectionCreateGlobalSection(s, dm->sf, PETSC_FALSE, PETSC_FALSE, &dm->defaultGlobalSection);
3588: PetscLayoutDestroy(&dm->map);
3589: PetscSectionGetValueLayout(PetscObjectComm((PetscObject)dm), dm->defaultGlobalSection, &dm->map);
3590: PetscSectionViewFromOptions(dm->defaultGlobalSection, NULL, "-global_section_view");
3591: }
3592: *section = dm->defaultGlobalSection;
3593: return(0);
3594: }
3596: /*@
3597: DMSetDefaultGlobalSection - Set the PetscSection encoding the global data layout for the DM.
3599: Input Parameters:
3600: + dm - The DM
3601: - section - The PetscSection, or NULL
3603: Level: intermediate
3605: Note: Any existing Section will be destroyed
3607: .seealso: DMGetDefaultGlobalSection(), DMSetDefaultSection()
3608: @*/
3609: PetscErrorCode DMSetDefaultGlobalSection(DM dm, PetscSection section)
3610: {
3616: PetscObjectReference((PetscObject)section);
3617: PetscSectionDestroy(&dm->defaultGlobalSection);
3618: dm->defaultGlobalSection = section;
3619: #ifdef PETSC_USE_DEBUG
3620: if (section) {DMDefaultSectionCheckConsistency_Internal(dm, dm->defaultSection, section);}
3621: #endif
3622: return(0);
3623: }
3625: /*@
3626: DMGetDefaultSF - Get the PetscSF encoding the parallel dof overlap for the DM. If it has not been set,
3627: it is created from the default PetscSection layouts in the DM.
3629: Input Parameter:
3630: . dm - The DM
3632: Output Parameter:
3633: . sf - The PetscSF
3635: Level: intermediate
3637: Note: This gets a borrowed reference, so the user should not destroy this PetscSF.
3639: .seealso: DMSetDefaultSF(), DMCreateDefaultSF()
3640: @*/
3641: PetscErrorCode DMGetDefaultSF(DM dm, PetscSF *sf)
3642: {
3643: PetscInt nroots;
3649: PetscSFGetGraph(dm->defaultSF, &nroots, NULL, NULL, NULL);
3650: if (nroots < 0) {
3651: PetscSection section, gSection;
3653: DMGetDefaultSection(dm, §ion);
3654: if (section) {
3655: DMGetDefaultGlobalSection(dm, &gSection);
3656: DMCreateDefaultSF(dm, section, gSection);
3657: } else {
3658: *sf = NULL;
3659: return(0);
3660: }
3661: }
3662: *sf = dm->defaultSF;
3663: return(0);
3664: }
3666: /*@
3667: DMSetDefaultSF - Set the PetscSF encoding the parallel dof overlap for the DM
3669: Input Parameters:
3670: + dm - The DM
3671: - sf - The PetscSF
3673: Level: intermediate
3675: Note: Any previous SF is destroyed
3677: .seealso: DMGetDefaultSF(), DMCreateDefaultSF()
3678: @*/
3679: PetscErrorCode DMSetDefaultSF(DM dm, PetscSF sf)
3680: {
3686: PetscSFDestroy(&dm->defaultSF);
3687: dm->defaultSF = sf;
3688: return(0);
3689: }
3691: /*@C
3692: DMCreateDefaultSF - Create the PetscSF encoding the parallel dof overlap for the DM based upon the PetscSections
3693: describing the data layout.
3695: Input Parameters:
3696: + dm - The DM
3697: . localSection - PetscSection describing the local data layout
3698: - globalSection - PetscSection describing the global data layout
3700: Level: intermediate
3702: .seealso: DMGetDefaultSF(), DMSetDefaultSF()
3703: @*/
3704: PetscErrorCode DMCreateDefaultSF(DM dm, PetscSection localSection, PetscSection globalSection)
3705: {
3706: MPI_Comm comm;
3707: PetscLayout layout;
3708: const PetscInt *ranges;
3709: PetscInt *local;
3710: PetscSFNode *remote;
3711: PetscInt pStart, pEnd, p, nroots, nleaves = 0, l;
3712: PetscMPIInt size, rank;
3716: PetscObjectGetComm((PetscObject)dm,&comm);
3718: MPI_Comm_size(comm, &size);
3719: MPI_Comm_rank(comm, &rank);
3720: PetscSectionGetChart(globalSection, &pStart, &pEnd);
3721: PetscSectionGetConstrainedStorageSize(globalSection, &nroots);
3722: PetscLayoutCreate(comm, &layout);
3723: PetscLayoutSetBlockSize(layout, 1);
3724: PetscLayoutSetLocalSize(layout, nroots);
3725: PetscLayoutSetUp(layout);
3726: PetscLayoutGetRanges(layout, &ranges);
3727: for (p = pStart; p < pEnd; ++p) {
3728: PetscInt gdof, gcdof;
3730: PetscSectionGetDof(globalSection, p, &gdof);
3731: PetscSectionGetConstraintDof(globalSection, p, &gcdof);
3732: 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));
3733: nleaves += gdof < 0 ? -(gdof+1)-gcdof : gdof-gcdof;
3734: }
3735: PetscMalloc1(nleaves, &local);
3736: PetscMalloc1(nleaves, &remote);
3737: for (p = pStart, l = 0; p < pEnd; ++p) {
3738: const PetscInt *cind;
3739: PetscInt dof, cdof, off, gdof, gcdof, goff, gsize, d, c;
3741: PetscSectionGetDof(localSection, p, &dof);
3742: PetscSectionGetOffset(localSection, p, &off);
3743: PetscSectionGetConstraintDof(localSection, p, &cdof);
3744: PetscSectionGetConstraintIndices(localSection, p, &cind);
3745: PetscSectionGetDof(globalSection, p, &gdof);
3746: PetscSectionGetConstraintDof(globalSection, p, &gcdof);
3747: PetscSectionGetOffset(globalSection, p, &goff);
3748: if (!gdof) continue; /* Censored point */
3749: gsize = gdof < 0 ? -(gdof+1)-gcdof : gdof-gcdof;
3750: if (gsize != dof-cdof) {
3751: 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);
3752: cdof = 0; /* Ignore constraints */
3753: }
3754: for (d = 0, c = 0; d < dof; ++d) {
3755: if ((c < cdof) && (cind[c] == d)) {++c; continue;}
3756: local[l+d-c] = off+d;
3757: }
3758: if (gdof < 0) {
3759: for (d = 0; d < gsize; ++d, ++l) {
3760: PetscInt offset = -(goff+1) + d, r;
3762: PetscFindInt(offset,size+1,ranges,&r);
3763: if (r < 0) r = -(r+2);
3764: 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);
3765: remote[l].rank = r;
3766: remote[l].index = offset - ranges[r];
3767: }
3768: } else {
3769: for (d = 0; d < gsize; ++d, ++l) {
3770: remote[l].rank = rank;
3771: remote[l].index = goff+d - ranges[rank];
3772: }
3773: }
3774: }
3775: if (l != nleaves) SETERRQ2(comm, PETSC_ERR_PLIB, "Iteration error, l %d != nleaves %d", l, nleaves);
3776: PetscLayoutDestroy(&layout);
3777: PetscSFSetGraph(dm->defaultSF, nroots, nleaves, local, PETSC_OWN_POINTER, remote, PETSC_OWN_POINTER);
3778: return(0);
3779: }
3781: /*@
3782: DMGetPointSF - Get the PetscSF encoding the parallel section point overlap for the DM.
3784: Input Parameter:
3785: . dm - The DM
3787: Output Parameter:
3788: . sf - The PetscSF
3790: Level: intermediate
3792: Note: This gets a borrowed reference, so the user should not destroy this PetscSF.
3794: .seealso: DMSetPointSF(), DMGetDefaultSF(), DMSetDefaultSF(), DMCreateDefaultSF()
3795: @*/
3796: PetscErrorCode DMGetPointSF(DM dm, PetscSF *sf)
3797: {
3801: *sf = dm->sf;
3802: return(0);
3803: }
3805: /*@
3806: DMSetPointSF - Set the PetscSF encoding the parallel section point overlap for the DM.
3808: Input Parameters:
3809: + dm - The DM
3810: - sf - The PetscSF
3812: Level: intermediate
3814: .seealso: DMGetPointSF(), DMGetDefaultSF(), DMSetDefaultSF(), DMCreateDefaultSF()
3815: @*/
3816: PetscErrorCode DMSetPointSF(DM dm, PetscSF sf)
3817: {
3823: PetscSFDestroy(&dm->sf);
3824: PetscObjectReference((PetscObject) sf);
3825: dm->sf = sf;
3826: return(0);
3827: }
3829: /*@
3830: DMGetDS - Get the PetscDS
3832: Input Parameter:
3833: . dm - The DM
3835: Output Parameter:
3836: . prob - The PetscDS
3838: Level: developer
3840: .seealso: DMSetDS()
3841: @*/
3842: PetscErrorCode DMGetDS(DM dm, PetscDS *prob)
3843: {
3847: *prob = dm->prob;
3848: return(0);
3849: }
3851: /*@
3852: DMSetDS - Set the PetscDS
3854: Input Parameters:
3855: + dm - The DM
3856: - prob - The PetscDS
3858: Level: developer
3860: .seealso: DMGetDS()
3861: @*/
3862: PetscErrorCode DMSetDS(DM dm, PetscDS prob)
3863: {
3869: PetscObjectReference((PetscObject) prob);
3870: PetscDSDestroy(&dm->prob);
3871: dm->prob = prob;
3872: return(0);
3873: }
3875: PetscErrorCode DMGetNumFields(DM dm, PetscInt *numFields)
3876: {
3881: PetscDSGetNumFields(dm->prob, numFields);
3882: return(0);
3883: }
3885: PetscErrorCode DMSetNumFields(DM dm, PetscInt numFields)
3886: {
3887: PetscInt Nf, f;
3892: PetscDSGetNumFields(dm->prob, &Nf);
3893: for (f = Nf; f < numFields; ++f) {
3894: PetscContainer obj;
3896: PetscContainerCreate(PetscObjectComm((PetscObject) dm), &obj);
3897: PetscDSSetDiscretization(dm->prob, f, (PetscObject) obj);
3898: PetscContainerDestroy(&obj);
3899: }
3900: return(0);
3901: }
3903: /*@
3904: DMGetField - Return the discretization object for a given DM field
3906: Not collective
3908: Input Parameters:
3909: + dm - The DM
3910: - f - The field number
3912: Output Parameter:
3913: . field - The discretization object
3915: Level: developer
3917: .seealso: DMSetField()
3918: @*/
3919: PetscErrorCode DMGetField(DM dm, PetscInt f, PetscObject *field)
3920: {
3925: PetscDSGetDiscretization(dm->prob, f, field);
3926: return(0);
3927: }
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: }
3953: PetscErrorCode DMRestrictHook_Coordinates(DM dm,DM dmc,void *ctx)
3954: {
3955: DM dm_coord,dmc_coord;
3957: Vec coords,ccoords;
3958: Mat inject;
3960: DMGetCoordinateDM(dm,&dm_coord);
3961: DMGetCoordinateDM(dmc,&dmc_coord);
3962: DMGetCoordinates(dm,&coords);
3963: DMGetCoordinates(dmc,&ccoords);
3964: if (coords && !ccoords) {
3965: DMCreateGlobalVector(dmc_coord,&ccoords);
3966: PetscObjectSetName((PetscObject)ccoords,"coordinates");
3967: DMCreateInjection(dmc_coord,dm_coord,&inject);
3968: MatRestrict(inject,coords,ccoords);
3969: MatDestroy(&inject);
3970: DMSetCoordinates(dmc,ccoords);
3971: VecDestroy(&ccoords);
3972: }
3973: return(0);
3974: }
3976: static PetscErrorCode DMSubDomainHook_Coordinates(DM dm,DM subdm,void *ctx)
3977: {
3978: DM dm_coord,subdm_coord;
3980: Vec coords,ccoords,clcoords;
3981: VecScatter *scat_i,*scat_g;
3983: DMGetCoordinateDM(dm,&dm_coord);
3984: DMGetCoordinateDM(subdm,&subdm_coord);
3985: DMGetCoordinates(dm,&coords);
3986: DMGetCoordinates(subdm,&ccoords);
3987: if (coords && !ccoords) {
3988: DMCreateGlobalVector(subdm_coord,&ccoords);
3989: PetscObjectSetName((PetscObject)ccoords,"coordinates");
3990: DMCreateLocalVector(subdm_coord,&clcoords);
3991: PetscObjectSetName((PetscObject)clcoords,"coordinates");
3992: DMCreateDomainDecompositionScatters(dm_coord,1,&subdm_coord,NULL,&scat_i,&scat_g);
3993: VecScatterBegin(scat_i[0],coords,ccoords,INSERT_VALUES,SCATTER_FORWARD);
3994: VecScatterBegin(scat_g[0],coords,clcoords,INSERT_VALUES,SCATTER_FORWARD);
3995: VecScatterEnd(scat_i[0],coords,ccoords,INSERT_VALUES,SCATTER_FORWARD);
3996: VecScatterEnd(scat_g[0],coords,clcoords,INSERT_VALUES,SCATTER_FORWARD);
3997: DMSetCoordinates(subdm,ccoords);
3998: DMSetCoordinatesLocal(subdm,clcoords);
3999: VecScatterDestroy(&scat_i[0]);
4000: VecScatterDestroy(&scat_g[0]);
4001: VecDestroy(&ccoords);
4002: VecDestroy(&clcoords);
4003: PetscFree(scat_i);
4004: PetscFree(scat_g);
4005: }
4006: return(0);
4007: }
4009: /*@
4010: DMGetDimension - Return the topological dimension of the DM
4012: Not collective
4014: Input Parameter:
4015: . dm - The DM
4017: Output Parameter:
4018: . dim - The topological dimension
4020: Level: beginner
4022: .seealso: DMSetDimension(), DMCreate()
4023: @*/
4024: PetscErrorCode DMGetDimension(DM dm, PetscInt *dim)
4025: {
4029: *dim = dm->dim;
4030: return(0);
4031: }
4033: /*@
4034: DMSetDimension - Set the topological dimension of the DM
4036: Collective on dm
4038: Input Parameters:
4039: + dm - The DM
4040: - dim - The topological dimension
4042: Level: beginner
4044: .seealso: DMGetDimension(), DMCreate()
4045: @*/
4046: PetscErrorCode DMSetDimension(DM dm, PetscInt dim)
4047: {
4051: dm->dim = dim;
4052: return(0);
4053: }
4055: /*@
4056: DMGetDimPoints - Get the half-open interval for all points of a given dimension
4058: Collective on DM
4060: Input Parameters:
4061: + dm - the DM
4062: - dim - the dimension
4064: Output Parameters:
4065: + pStart - The first point of the given dimension
4066: . pEnd - The first point following points of the given dimension
4068: Note:
4069: The points are vertices in the Hasse diagram encoding the topology. This is explained in
4070: http://arxiv.org/abs/0908.4427. If not points exist of this dimension in the storage scheme,
4071: then the interval is empty.
4073: Level: intermediate
4075: .keywords: point, Hasse Diagram, dimension
4076: .seealso: DMPLEX, DMPlexGetDepthStratum(), DMPlexGetHeightStratum()
4077: @*/
4078: PetscErrorCode DMGetDimPoints(DM dm, PetscInt dim, PetscInt *pStart, PetscInt *pEnd)
4079: {
4080: PetscInt d;
4085: DMGetDimension(dm, &d);
4086: if ((dim < 0) || (dim > d)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension %d 1", dim, d);
4087: (*dm->ops->getdimpoints)(dm, dim, pStart, pEnd);
4088: return(0);
4089: }
4091: /*@
4092: DMSetCoordinates - Sets into the DM a global vector that holds the coordinates
4094: Collective on DM
4096: Input Parameters:
4097: + dm - the DM
4098: - c - coordinate vector
4100: Note:
4101: The coordinates do include those for ghost points, which are in the local vector
4103: Level: intermediate
4105: .keywords: distributed array, get, corners, nodes, local indices, coordinates
4106: .seealso: DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLoca(), DMGetCoordinateDM()
4107: @*/
4108: PetscErrorCode DMSetCoordinates(DM dm, Vec c)
4109: {
4115: PetscObjectReference((PetscObject) c);
4116: VecDestroy(&dm->coordinates);
4117: dm->coordinates = c;
4118: VecDestroy(&dm->coordinatesLocal);
4119: DMCoarsenHookAdd(dm,DMRestrictHook_Coordinates,NULL,NULL);
4120: DMSubDomainHookAdd(dm,DMSubDomainHook_Coordinates,NULL,NULL);
4121: return(0);
4122: }
4124: /*@
4125: DMSetCoordinatesLocal - Sets into the DM a local vector that holds the coordinates
4127: Collective on DM
4129: Input Parameters:
4130: + dm - the DM
4131: - c - coordinate vector
4133: Note:
4134: The coordinates of ghost points can be set using DMSetCoordinates()
4135: followed by DMGetCoordinatesLocal(). This is intended to enable the
4136: setting of ghost coordinates outside of the domain.
4138: Level: intermediate
4140: .keywords: distributed array, get, corners, nodes, local indices, coordinates
4141: .seealso: DMGetCoordinatesLocal(), DMSetCoordinates(), DMGetCoordinates(), DMGetCoordinateDM()
4142: @*/
4143: PetscErrorCode DMSetCoordinatesLocal(DM dm, Vec c)
4144: {
4150: PetscObjectReference((PetscObject) c);
4151: VecDestroy(&dm->coordinatesLocal);
4153: dm->coordinatesLocal = c;
4155: VecDestroy(&dm->coordinates);
4156: return(0);
4157: }
4159: /*@
4160: DMGetCoordinates - Gets a global vector with the coordinates associated with the DM.
4162: Not Collective
4164: Input Parameter:
4165: . dm - the DM
4167: Output Parameter:
4168: . c - global coordinate vector
4170: Note:
4171: This is a borrowed reference, so the user should NOT destroy this vector
4173: Each process has only the local coordinates (does NOT have the ghost coordinates).
4175: For DMDA, in two and three dimensions coordinates are interlaced (x_0,y_0,x_1,y_1,...)
4176: and (x_0,y_0,z_0,x_1,y_1,z_1...)
4178: Level: intermediate
4180: .keywords: distributed array, get, corners, nodes, local indices, coordinates
4181: .seealso: DMSetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM()
4182: @*/
4183: PetscErrorCode DMGetCoordinates(DM dm, Vec *c)
4184: {
4190: if (!dm->coordinates && dm->coordinatesLocal) {
4191: DM cdm = NULL;
4193: DMGetCoordinateDM(dm, &cdm);
4194: DMCreateGlobalVector(cdm, &dm->coordinates);
4195: PetscObjectSetName((PetscObject) dm->coordinates, "coordinates");
4196: DMLocalToGlobalBegin(cdm, dm->coordinatesLocal, INSERT_VALUES, dm->coordinates);
4197: DMLocalToGlobalEnd(cdm, dm->coordinatesLocal, INSERT_VALUES, dm->coordinates);
4198: }
4199: *c = dm->coordinates;
4200: return(0);
4201: }
4203: /*@
4204: DMGetCoordinatesLocal - Gets a local vector with the coordinates associated with the DM.
4206: Collective on DM
4208: Input Parameter:
4209: . dm - the DM
4211: Output Parameter:
4212: . c - coordinate vector
4214: Note:
4215: This is a borrowed reference, so the user should NOT destroy this vector
4217: Each process has the local and ghost coordinates
4219: For DMDA, in two and three dimensions coordinates are interlaced (x_0,y_0,x_1,y_1,...)
4220: and (x_0,y_0,z_0,x_1,y_1,z_1...)
4222: Level: intermediate
4224: .keywords: distributed array, get, corners, nodes, local indices, coordinates
4225: .seealso: DMSetCoordinatesLocal(), DMGetCoordinates(), DMSetCoordinates(), DMGetCoordinateDM()
4226: @*/
4227: PetscErrorCode DMGetCoordinatesLocal(DM dm, Vec *c)
4228: {
4234: if (!dm->coordinatesLocal && dm->coordinates) {
4235: DM cdm = NULL;
4237: DMGetCoordinateDM(dm, &cdm);
4238: DMCreateLocalVector(cdm, &dm->coordinatesLocal);
4239: PetscObjectSetName((PetscObject) dm->coordinatesLocal, "coordinates");
4240: DMGlobalToLocalBegin(cdm, dm->coordinates, INSERT_VALUES, dm->coordinatesLocal);
4241: DMGlobalToLocalEnd(cdm, dm->coordinates, INSERT_VALUES, dm->coordinatesLocal);
4242: }
4243: *c = dm->coordinatesLocal;
4244: return(0);
4245: }
4247: /*@
4248: DMGetCoordinateDM - Gets the DM that prescribes coordinate layout and scatters between global and local coordinates
4250: Collective on DM
4252: Input Parameter:
4253: . dm - the DM
4255: Output Parameter:
4256: . cdm - coordinate DM
4258: Level: intermediate
4260: .keywords: distributed array, get, corners, nodes, local indices, coordinates
4261: .seealso: DMSetCoordinateDM(), DMSetCoordinates(), DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLocal()
4262: @*/
4263: PetscErrorCode DMGetCoordinateDM(DM dm, DM *cdm)
4264: {
4270: if (!dm->coordinateDM) {
4271: if (!dm->ops->createcoordinatedm) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unable to create coordinates for this DM");
4272: (*dm->ops->createcoordinatedm)(dm, &dm->coordinateDM);
4273: }
4274: *cdm = dm->coordinateDM;
4275: return(0);
4276: }
4278: /*@
4279: DMSetCoordinateDM - Sets the DM that prescribes coordinate layout and scatters between global and local coordinates
4281: Logically Collective on DM
4283: Input Parameters:
4284: + dm - the DM
4285: - cdm - coordinate DM
4287: Level: intermediate
4289: .keywords: distributed array, get, corners, nodes, local indices, coordinates
4290: .seealso: DMGetCoordinateDM(), DMSetCoordinates(), DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLocal()
4291: @*/
4292: PetscErrorCode DMSetCoordinateDM(DM dm, DM cdm)
4293: {
4299: PetscObjectReference((PetscObject)cdm);
4300: DMDestroy(&dm->coordinateDM);
4301: dm->coordinateDM = cdm;
4302: return(0);
4303: }
4305: /*@
4306: DMGetCoordinateDim - Retrieve the dimension of embedding space for coordinate values.
4308: Not Collective
4310: Input Parameter:
4311: . dm - The DM object
4313: Output Parameter:
4314: . dim - The embedding dimension
4316: Level: intermediate
4318: .keywords: mesh, coordinates
4319: .seealso: DMSetCoordinateDim(), DMGetCoordinateSection(), DMGetCoordinateDM(), DMGetDefaultSection(), DMSetDefaultSection()
4320: @*/
4321: PetscErrorCode DMGetCoordinateDim(DM dm, PetscInt *dim)
4322: {
4326: if (dm->dimEmbed == PETSC_DEFAULT) {
4327: dm->dimEmbed = dm->dim;
4328: }
4329: *dim = dm->dimEmbed;
4330: return(0);
4331: }
4333: /*@
4334: DMSetCoordinateDim - Set the dimension of the embedding space for coordinate values.
4336: Not Collective
4338: Input Parameters:
4339: + dm - The DM object
4340: - dim - The embedding dimension
4342: Level: intermediate
4344: .keywords: mesh, coordinates
4345: .seealso: DMGetCoordinateDim(), DMSetCoordinateSection(), DMGetCoordinateSection(), DMGetDefaultSection(), DMSetDefaultSection()
4346: @*/
4347: PetscErrorCode DMSetCoordinateDim(DM dm, PetscInt dim)
4348: {
4351: dm->dimEmbed = dim;
4352: return(0);
4353: }
4355: /*@
4356: DMGetCoordinateSection - Retrieve the layout of coordinate values over the mesh.
4358: Not Collective
4360: Input Parameter:
4361: . dm - The DM object
4363: Output Parameter:
4364: . section - The PetscSection object
4366: Level: intermediate
4368: .keywords: mesh, coordinates
4369: .seealso: DMGetCoordinateDM(), DMGetDefaultSection(), DMSetDefaultSection()
4370: @*/
4371: PetscErrorCode DMGetCoordinateSection(DM dm, PetscSection *section)
4372: {
4373: DM cdm;
4379: DMGetCoordinateDM(dm, &cdm);
4380: DMGetDefaultSection(cdm, section);
4381: return(0);
4382: }
4384: /*@
4385: DMSetCoordinateSection - Set the layout of coordinate values over the mesh.
4387: Not Collective
4389: Input Parameters:
4390: + dm - The DM object
4391: . dim - The embedding dimension, or PETSC_DETERMINE
4392: - section - The PetscSection object
4394: Level: intermediate
4396: .keywords: mesh, coordinates
4397: .seealso: DMGetCoordinateSection(), DMGetDefaultSection(), DMSetDefaultSection()
4398: @*/
4399: PetscErrorCode DMSetCoordinateSection(DM dm, PetscInt dim, PetscSection section)
4400: {
4401: DM cdm;
4407: DMGetCoordinateDM(dm, &cdm);
4408: DMSetDefaultSection(cdm, section);
4409: if (dim == PETSC_DETERMINE) {
4410: PetscInt d = PETSC_DEFAULT;
4411: PetscInt pStart, pEnd, vStart, vEnd, v, dd;
4413: PetscSectionGetChart(section, &pStart, &pEnd);
4414: DMGetDimPoints(dm, 0, &vStart, &vEnd);
4415: pStart = PetscMax(vStart, pStart);
4416: pEnd = PetscMin(vEnd, pEnd);
4417: for (v = pStart; v < pEnd; ++v) {
4418: PetscSectionGetDof(section, v, &dd);
4419: if (dd) {d = dd; break;}
4420: }
4421: if (d < 0) d = PETSC_DEFAULT;
4422: DMSetCoordinateDim(dm, d);
4423: }
4424: return(0);
4425: }
4427: /*@C
4428: DMGetPeriodicity - Get the description of mesh periodicity
4430: Input Parameters:
4431: . dm - The DM object
4433: Output Parameters:
4434: + per - Whether the DM is periodic or not
4435: . maxCell - Over distances greater than this, we can assume a point has crossed over to another sheet, when trying to localize cell coordinates
4436: . L - If we assume the mesh is a torus, this is the length of each coordinate
4437: - bd - This describes the type of periodicity in each topological dimension
4439: Level: developer
4441: .seealso: DMGetPeriodicity()
4442: @*/
4443: PetscErrorCode DMGetPeriodicity(DM dm, PetscBool *per, const PetscReal **maxCell, const PetscReal **L, const DMBoundaryType **bd)
4444: {
4447: if (per) *per = dm->periodic;
4448: if (L) *L = dm->L;
4449: if (maxCell) *maxCell = dm->maxCell;
4450: if (bd) *bd = dm->bdtype;
4451: return(0);
4452: }
4454: /*@C
4455: DMSetPeriodicity - Set the description of mesh periodicity
4457: Input Parameters:
4458: + dm - The DM object
4459: . per - Whether the DM is periodic or not. If maxCell is not provided, coordinates need to be localized
4460: . maxCell - Over distances greater than this, we can assume a point has crossed over to another sheet, when trying to localize cell coordinates
4461: . L - If we assume the mesh is a torus, this is the length of each coordinate
4462: - bd - This describes the type of periodicity in each topological dimension
4464: Level: developer
4466: .seealso: DMGetPeriodicity()
4467: @*/
4468: PetscErrorCode DMSetPeriodicity(DM dm, PetscBool per, const PetscReal maxCell[], const PetscReal L[], const DMBoundaryType bd[])
4469: {
4470: PetscInt dim, d;
4476: if (maxCell) {
4480: }
4481: PetscFree3(dm->L,dm->maxCell,dm->bdtype);
4482: DMGetDimension(dm, &dim);
4483: if (maxCell) {
4484: PetscMalloc3(dim,&dm->L,dim,&dm->maxCell,dim,&dm->bdtype);
4485: for (d = 0; d < dim; ++d) {dm->L[d] = L[d]; dm->maxCell[d] = maxCell[d]; dm->bdtype[d] = bd[d];}
4486: dm->periodic = PETSC_TRUE;
4487: } else {
4488: dm->periodic = per;
4489: }
4490: return(0);
4491: }
4493: /*@
4494: 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.
4496: Input Parameters:
4497: + dm - The DM
4498: . in - The input coordinate point (dim numbers)
4499: - endpoint - Include the endpoint L_i
4501: Output Parameter:
4502: . out - The localized coordinate point
4504: Level: developer
4506: .seealso: DMLocalizeCoordinates(), DMLocalizeAddCoordinate()
4507: @*/
4508: PetscErrorCode DMLocalizeCoordinate(DM dm, const PetscScalar in[], PetscBool endpoint, PetscScalar out[])
4509: {
4510: PetscInt dim, d;
4514: DMGetCoordinateDim(dm, &dim);
4515: if (!dm->maxCell) {
4516: for (d = 0; d < dim; ++d) out[d] = in[d];
4517: } else {
4518: if (endpoint) {
4519: for (d = 0; d < dim; ++d) {
4520: if ((PetscAbsReal(PetscRealPart(in[d])/dm->L[d] - PetscFloorReal(PetscRealPart(in[d])/dm->L[d])) < PETSC_SMALL) && (PetscRealPart(in[d])/dm->L[d] > PETSC_SMALL)) {
4521: out[d] = in[d] - dm->L[d]*(PetscFloorReal(PetscRealPart(in[d])/dm->L[d]) - 1);
4522: } else {
4523: out[d] = in[d] - dm->L[d]*PetscFloorReal(PetscRealPart(in[d])/dm->L[d]);
4524: }
4525: }
4526: } else {
4527: for (d = 0; d < dim; ++d) {
4528: out[d] = in[d] - dm->L[d]*PetscFloorReal(PetscRealPart(in[d])/dm->L[d]);
4529: }
4530: }
4531: }
4532: return(0);
4533: }
4535: /*
4536: 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.
4538: Input Parameters:
4539: + dm - The DM
4540: . dim - The spatial dimension
4541: . anchor - The anchor point, the input point can be no more than maxCell away from it
4542: - in - The input coordinate point (dim numbers)
4544: Output Parameter:
4545: . out - The localized coordinate point
4547: Level: developer
4549: 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
4551: .seealso: DMLocalizeCoordinates(), DMLocalizeAddCoordinate()
4552: */
4553: PetscErrorCode DMLocalizeCoordinate_Internal(DM dm, PetscInt dim, const PetscScalar anchor[], const PetscScalar in[], PetscScalar out[])
4554: {
4555: PetscInt d;
4558: if (!dm->maxCell) {
4559: for (d = 0; d < dim; ++d) out[d] = in[d];
4560: } else {
4561: for (d = 0; d < dim; ++d) {
4562: if ((dm->bdtype[d] != DM_BOUNDARY_NONE) && (PetscAbsScalar(anchor[d] - in[d]) > dm->maxCell[d])) {
4563: out[d] = PetscRealPart(anchor[d]) > PetscRealPart(in[d]) ? dm->L[d] + in[d] : in[d] - dm->L[d];
4564: } else {
4565: out[d] = in[d];
4566: }
4567: }
4568: }
4569: return(0);
4570: }
4571: PetscErrorCode DMLocalizeCoordinateReal_Internal(DM dm, PetscInt dim, const PetscReal anchor[], const PetscReal in[], PetscReal out[])
4572: {
4573: PetscInt d;
4576: if (!dm->maxCell) {
4577: for (d = 0; d < dim; ++d) out[d] = in[d];
4578: } else {
4579: for (d = 0; d < dim; ++d) {
4580: if ((dm->bdtype[d] != DM_BOUNDARY_NONE) && (PetscAbsReal(anchor[d] - in[d]) > dm->maxCell[d])) {
4581: out[d] = anchor[d] > in[d] ? dm->L[d] + in[d] : in[d] - dm->L[d];
4582: } else {
4583: out[d] = in[d];
4584: }
4585: }
4586: }
4587: return(0);
4588: }
4590: /*
4591: 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.
4593: Input Parameters:
4594: + dm - The DM
4595: . dim - The spatial dimension
4596: . anchor - The anchor point, the input point can be no more than maxCell away from it
4597: . in - The input coordinate delta (dim numbers)
4598: - out - The input coordinate point (dim numbers)
4600: Output Parameter:
4601: . out - The localized coordinate in + out
4603: Level: developer
4605: 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
4607: .seealso: DMLocalizeCoordinates(), DMLocalizeCoordinate()
4608: */
4609: PetscErrorCode DMLocalizeAddCoordinate_Internal(DM dm, PetscInt dim, const PetscScalar anchor[], const PetscScalar in[], PetscScalar out[])
4610: {
4611: PetscInt d;
4614: if (!dm->maxCell) {
4615: for (d = 0; d < dim; ++d) out[d] += in[d];
4616: } else {
4617: for (d = 0; d < dim; ++d) {
4618: if ((dm->bdtype[d] != DM_BOUNDARY_NONE) && (PetscAbsScalar(anchor[d] - in[d]) > dm->maxCell[d])) {
4619: out[d] += PetscRealPart(anchor[d]) > PetscRealPart(in[d]) ? dm->L[d] + in[d] : in[d] - dm->L[d];
4620: } else {
4621: out[d] += in[d];
4622: }
4623: }
4624: }
4625: return(0);
4626: }
4628: /*@
4629: DMGetCoordinatesLocalized - Check if the DM coordinates have been localized for cells
4631: Input Parameter:
4632: . dm - The DM
4634: Output Parameter:
4635: areLocalized - True if localized
4637: Level: developer
4639: .seealso: DMLocalizeCoordinates()
4640: @*/
4641: PetscErrorCode DMGetCoordinatesLocalized(DM dm,PetscBool *areLocalized)
4642: {
4643: DM cdm;
4644: PetscSection coordSection;
4645: PetscInt cStart, cEnd, c, sStart, sEnd, dof;
4646: PetscBool alreadyLocalized, alreadyLocalizedGlobal;
4651: if (!dm->periodic) {
4652: *areLocalized = PETSC_FALSE;
4653: return(0);
4654: }
4655: /* We need some generic way of refering to cells/vertices */
4656: DMGetCoordinateDM(dm, &cdm);
4657: {
4658: PetscBool isplex;
4660: PetscObjectTypeCompare((PetscObject) cdm, DMPLEX, &isplex);
4661: if (isplex) {
4662: DMPlexGetHeightStratum(cdm, 0, &cStart, &cEnd);
4663: } else SETERRQ(PetscObjectComm((PetscObject) cdm), PETSC_ERR_ARG_WRONG, "Coordinate localization requires a DMPLEX coordinate DM");
4664: }
4665: DMGetCoordinateSection(dm, &coordSection);
4666: PetscSectionGetChart(coordSection,&sStart,&sEnd);
4667: alreadyLocalized = alreadyLocalizedGlobal = PETSC_FALSE;
4668: for (c = cStart; c < cEnd; ++c) {
4669: if (c < sStart || c >= sEnd) {
4670: alreadyLocalized = PETSC_FALSE;
4671: break;
4672: }
4673: PetscSectionGetDof(coordSection, c, &dof);
4674: if (dof) {
4675: alreadyLocalized = PETSC_TRUE;
4676: break;
4677: }
4678: }
4679: MPI_Allreduce(&alreadyLocalized,&alreadyLocalizedGlobal,1,MPIU_BOOL,MPI_LOR,PetscObjectComm((PetscObject)dm));
4680: *areLocalized = alreadyLocalizedGlobal;
4681: return(0);
4682: }
4685: /*@
4686: DMLocalizeCoordinates - If a mesh is periodic, create local coordinates for cells having periodic faces
4688: Input Parameter:
4689: . dm - The DM
4691: Level: developer
4693: .seealso: DMLocalizeCoordinate(), DMLocalizeAddCoordinate()
4694: @*/
4695: PetscErrorCode DMLocalizeCoordinates(DM dm)
4696: {
4697: DM cdm;
4698: PetscSection coordSection, cSection;
4699: Vec coordinates, cVec;
4700: PetscScalar *coords, *coords2, *anchor, *localized;
4701: PetscInt Nc, vStart, vEnd, v, sStart, sEnd, newStart = PETSC_MAX_INT, newEnd = PETSC_MIN_INT, dof, d, off, off2, bs, coordSize;
4702: PetscBool alreadyLocalized, alreadyLocalizedGlobal;
4703: PetscInt maxHeight = 0, h;
4704: PetscInt *pStart = NULL, *pEnd = NULL;
4709: if (!dm->periodic) return(0);
4710: /* We need some generic way of refering to cells/vertices */
4711: DMGetCoordinateDM(dm, &cdm);
4712: {
4713: PetscBool isplex;
4715: PetscObjectTypeCompare((PetscObject) cdm, DMPLEX, &isplex);
4716: if (isplex) {
4717: DMPlexGetDepthStratum(cdm, 0, &vStart, &vEnd);
4718: DMPlexGetMaxProjectionHeight(cdm,&maxHeight);
4719: DMGetWorkArray(dm,2*(maxHeight + 1),PETSC_INT,&pStart);
4720: pEnd = &pStart[maxHeight + 1];
4721: newStart = vStart;
4722: newEnd = vEnd;
4723: for (h = 0; h <= maxHeight; h++) {
4724: DMPlexGetHeightStratum(cdm, h, &pStart[h], &pEnd[h]);
4725: newStart = PetscMin(newStart,pStart[h]);
4726: newEnd = PetscMax(newEnd,pEnd[h]);
4727: }
4728: } else SETERRQ(PetscObjectComm((PetscObject) cdm), PETSC_ERR_ARG_WRONG, "Coordinate localization requires a DMPLEX coordinate DM");
4729: }
4730: DMGetCoordinatesLocal(dm, &coordinates);
4731: DMGetCoordinateSection(dm, &coordSection);
4732: VecGetBlockSize(coordinates, &bs);
4733: PetscSectionGetChart(coordSection,&sStart,&sEnd);
4735: PetscSectionCreate(PetscObjectComm((PetscObject) dm), &cSection);
4736: PetscSectionSetNumFields(cSection, 1);
4737: PetscSectionGetFieldComponents(coordSection, 0, &Nc);
4738: PetscSectionSetFieldComponents(cSection, 0, Nc);
4739: PetscSectionSetChart(cSection, newStart, newEnd);
4741: DMGetWorkArray(dm, 2 * bs, PETSC_SCALAR, &anchor);
4742: localized = &anchor[bs];
4743: alreadyLocalized = alreadyLocalizedGlobal = PETSC_TRUE;
4744: for (h = 0; h <= maxHeight; h++) {
4745: PetscInt cStart = pStart[h], cEnd = pEnd[h], c;
4747: for (c = cStart; c < cEnd; ++c) {
4748: PetscScalar *cellCoords = NULL;
4749: PetscInt b;
4751: if (c < sStart || c >= sEnd) alreadyLocalized = PETSC_FALSE;
4752: DMPlexVecGetClosure(cdm, coordSection, coordinates, c, &dof, &cellCoords);
4753: for (b = 0; b < bs; ++b) anchor[b] = cellCoords[b];
4754: for (d = 0; d < dof/bs; ++d) {
4755: DMLocalizeCoordinate_Internal(dm, bs, anchor, &cellCoords[d*bs], localized);
4756: for (b = 0; b < bs; b++) {
4757: if (cellCoords[d*bs + b] != localized[b]) break;
4758: }
4759: if (b < bs) break;
4760: }
4761: if (d < dof/bs) {
4762: if (c >= sStart && c < sEnd) {
4763: PetscInt cdof;
4765: PetscSectionGetDof(coordSection, c, &cdof);
4766: if (cdof != dof) alreadyLocalized = PETSC_FALSE;
4767: }
4768: PetscSectionSetDof(cSection, c, dof);
4769: PetscSectionSetFieldDof(cSection, c, 0, dof);
4770: }
4771: DMPlexVecRestoreClosure(cdm, coordSection, coordinates, c, &dof, &cellCoords);
4772: }
4773: }
4774: MPI_Allreduce(&alreadyLocalized,&alreadyLocalizedGlobal,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)dm));
4775: if (alreadyLocalizedGlobal) {
4776: DMRestoreWorkArray(dm, 2 * bs, PETSC_SCALAR, &anchor);
4777: PetscSectionDestroy(&cSection);
4778: DMRestoreWorkArray(dm,2*(maxHeight + 1),PETSC_INT,&pStart);
4779: return(0);
4780: }
4781: for (v = vStart; v < vEnd; ++v) {
4782: PetscSectionGetDof(coordSection, v, &dof);
4783: PetscSectionSetDof(cSection, v, dof);
4784: PetscSectionSetFieldDof(cSection, v, 0, dof);
4785: }
4786: PetscSectionSetUp(cSection);
4787: PetscSectionGetStorageSize(cSection, &coordSize);
4788: VecCreate(PETSC_COMM_SELF, &cVec);
4789: PetscObjectSetName((PetscObject)cVec,"coordinates");
4790: VecSetBlockSize(cVec, bs);
4791: VecSetSizes(cVec, coordSize, PETSC_DETERMINE);
4792: VecSetType(cVec, VECSTANDARD);
4793: VecGetArrayRead(coordinates, (const PetscScalar**)&coords);
4794: VecGetArray(cVec, &coords2);
4795: for (v = vStart; v < vEnd; ++v) {
4796: PetscSectionGetDof(coordSection, v, &dof);
4797: PetscSectionGetOffset(coordSection, v, &off);
4798: PetscSectionGetOffset(cSection, v, &off2);
4799: for (d = 0; d < dof; ++d) coords2[off2+d] = coords[off+d];
4800: }
4801: for (h = 0; h <= maxHeight; h++) {
4802: PetscInt cStart = pStart[h], cEnd = pEnd[h], c;
4804: for (c = cStart; c < cEnd; ++c) {
4805: PetscScalar *cellCoords = NULL;
4806: PetscInt b, cdof;
4808: PetscSectionGetDof(cSection,c,&cdof);
4809: if (!cdof) continue;
4810: DMPlexVecGetClosure(cdm, coordSection, coordinates, c, &dof, &cellCoords);
4811: PetscSectionGetOffset(cSection, c, &off2);
4812: for (b = 0; b < bs; ++b) anchor[b] = cellCoords[b];
4813: for (d = 0; d < dof/bs; ++d) {DMLocalizeCoordinate_Internal(dm, bs, anchor, &cellCoords[d*bs], &coords2[off2+d*bs]);}
4814: DMPlexVecRestoreClosure(cdm, coordSection, coordinates, c, &dof, &cellCoords);
4815: }
4816: }
4817: DMRestoreWorkArray(dm, 2 * bs, PETSC_SCALAR, &anchor);
4818: DMRestoreWorkArray(dm,2*(maxHeight + 1),PETSC_INT,&pStart);
4819: VecRestoreArrayRead(coordinates, (const PetscScalar**)&coords);
4820: VecRestoreArray(cVec, &coords2);
4821: DMSetCoordinateSection(dm, PETSC_DETERMINE, cSection);
4822: DMSetCoordinatesLocal(dm, cVec);
4823: VecDestroy(&cVec);
4824: PetscSectionDestroy(&cSection);
4825: return(0);
4826: }
4828: /*@
4829: DMLocatePoints - Locate the points in v in the mesh and return a PetscSF of the containing cells
4831: Collective on Vec v (see explanation below)
4833: Input Parameters:
4834: + dm - The DM
4835: . v - The Vec of points
4836: . ltype - The type of point location, e.g. DM_POINTLOCATION_NONE or DM_POINTLOCATION_NEAREST
4837: - cells - Points to either NULL, or a PetscSF with guesses for which cells contain each point.
4839: Output Parameter:
4840: + v - The Vec of points, which now contains the nearest mesh points to the given points if DM_POINTLOCATION_NEAREST is used
4841: - cells - The PetscSF containing the ranks and local indices of the containing points.
4844: Level: developer
4846: Notes:
4847: To do a search of the local cells of the mesh, v should have PETSC_COMM_SELF as its communicator.
4848: To do a search of all the cells in the distributed mesh, v should have the same communicator as dm.
4850: If *cellSF is NULL on input, a PetscSF will be created.
4851: If *cellSF is not NULL on input, it should point to an existing PetscSF, whose graph will be used as initial guesses.
4853: An array that maps each point to its containing cell can be obtained with
4855: $ const PetscSFNode *cells;
4856: $ PetscInt nFound;
4857: $ const PetscSFNode *found;
4858: $
4859: $ PetscSFGetGraph(cells,NULL,&nFound,&found,&cells);
4861: Where cells[i].rank is the rank of the cell containing point found[i] (or i if found == NULL), and cells[i].index is
4862: the index of the cell in its rank's local numbering.
4864: .keywords: point location, mesh
4865: .seealso: DMSetCoordinates(), DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLocal(), DMPointLocationType
4866: @*/
4867: PetscErrorCode DMLocatePoints(DM dm, Vec v, DMPointLocationType ltype, PetscSF *cellSF)
4868: {
4875: if (*cellSF) {
4876: PetscMPIInt result;
4879: MPI_Comm_compare(PetscObjectComm((PetscObject)v),PetscObjectComm((PetscObject)*cellSF),&result);
4880: if (result != MPI_IDENT && result != MPI_CONGRUENT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"cellSF must have a communicator congruent to v's");
4881: } else {
4882: PetscSFCreate(PetscObjectComm((PetscObject)v),cellSF);
4883: }
4884: PetscLogEventBegin(DM_LocatePoints,dm,0,0,0);
4885: if (dm->ops->locatepoints) {
4886: (*dm->ops->locatepoints)(dm,v,ltype,*cellSF);
4887: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Point location not available for this DM");
4888: PetscLogEventEnd(DM_LocatePoints,dm,0,0,0);
4889: return(0);
4890: }
4892: /*@
4893: DMGetOutputDM - Retrieve the DM associated with the layout for output
4895: Input Parameter:
4896: . dm - The original DM
4898: Output Parameter:
4899: . odm - The DM which provides the layout for output
4901: Level: intermediate
4903: .seealso: VecView(), DMGetDefaultGlobalSection()
4904: @*/
4905: PetscErrorCode DMGetOutputDM(DM dm, DM *odm)
4906: {
4907: PetscSection section;
4908: PetscBool hasConstraints, ghasConstraints;
4914: DMGetDefaultSection(dm, §ion);
4915: PetscSectionHasConstraints(section, &hasConstraints);
4916: MPI_Allreduce(&hasConstraints, &ghasConstraints, 1, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject) dm));
4917: if (!ghasConstraints) {
4918: *odm = dm;
4919: return(0);
4920: }
4921: if (!dm->dmBC) {
4922: PetscDS ds;
4923: PetscSection newSection, gsection;
4924: PetscSF sf;
4926: DMClone(dm, &dm->dmBC);
4927: DMGetDS(dm, &ds);
4928: DMSetDS(dm->dmBC, ds);
4929: PetscSectionClone(section, &newSection);
4930: DMSetDefaultSection(dm->dmBC, newSection);
4931: PetscSectionDestroy(&newSection);
4932: DMGetPointSF(dm->dmBC, &sf);
4933: PetscSectionCreateGlobalSection(section, sf, PETSC_TRUE, PETSC_FALSE, &gsection);
4934: DMSetDefaultGlobalSection(dm->dmBC, gsection);
4935: PetscSectionDestroy(&gsection);
4936: }
4937: *odm = dm->dmBC;
4938: return(0);
4939: }
4941: /*@
4942: DMGetOutputSequenceNumber - Retrieve the sequence number/value for output
4944: Input Parameter:
4945: . dm - The original DM
4947: Output Parameters:
4948: + num - The output sequence number
4949: - val - The output sequence value
4951: Level: intermediate
4953: Note: This is intended for output that should appear in sequence, for instance
4954: a set of timesteps in an HDF5 file, or a set of realizations of a stochastic system.
4956: .seealso: VecView()
4957: @*/
4958: PetscErrorCode DMGetOutputSequenceNumber(DM dm, PetscInt *num, PetscReal *val)
4959: {
4964: return(0);
4965: }
4967: /*@
4968: DMSetOutputSequenceNumber - Set the sequence number/value for output
4970: Input Parameters:
4971: + dm - The original DM
4972: . num - The output sequence number
4973: - val - The output sequence value
4975: Level: intermediate
4977: Note: This is intended for output that should appear in sequence, for instance
4978: a set of timesteps in an HDF5 file, or a set of realizations of a stochastic system.
4980: .seealso: VecView()
4981: @*/
4982: PetscErrorCode DMSetOutputSequenceNumber(DM dm, PetscInt num, PetscReal val)
4983: {
4986: dm->outputSequenceNum = num;
4987: dm->outputSequenceVal = val;
4988: return(0);
4989: }
4991: /*@C
4992: DMOutputSequenceLoad - Retrieve the sequence value from a Viewer
4994: Input Parameters:
4995: + dm - The original DM
4996: . name - The sequence name
4997: - num - The output sequence number
4999: Output Parameter:
5000: . val - The output sequence value
5002: Level: intermediate
5004: Note: This is intended for output that should appear in sequence, for instance
5005: a set of timesteps in an HDF5 file, or a set of realizations of a stochastic system.
5007: .seealso: DMGetOutputSequenceNumber(), DMSetOutputSequenceNumber(), VecView()
5008: @*/
5009: PetscErrorCode DMOutputSequenceLoad(DM dm, PetscViewer viewer, const char *name, PetscInt num, PetscReal *val)
5010: {
5011: PetscBool ishdf5;
5018: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERHDF5, &ishdf5);
5019: if (ishdf5) {
5020: #if defined(PETSC_HAVE_HDF5)
5021: PetscScalar value;
5023: DMSequenceLoad_HDF5_Internal(dm, name, num, &value, viewer);
5024: *val = PetscRealPart(value);
5025: #endif
5026: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid viewer; open viewer with PetscViewerHDF5Open()");
5027: return(0);
5028: }
5030: /*@
5031: DMGetUseNatural - Get the flag for creating a mapping to the natural order on distribution
5033: Not collective
5035: Input Parameter:
5036: . dm - The DM
5038: Output Parameter:
5039: . useNatural - The flag to build the mapping to a natural order during distribution
5041: Level: beginner
5043: .seealso: DMSetUseNatural(), DMCreate()
5044: @*/
5045: PetscErrorCode DMGetUseNatural(DM dm, PetscBool *useNatural)
5046: {
5050: *useNatural = dm->useNatural;
5051: return(0);
5052: }
5054: /*@
5055: DMSetUseNatural - Set the flag for creating a mapping to the natural order on distribution
5057: Collective on dm
5059: Input Parameters:
5060: + dm - The DM
5061: - useNatural - The flag to build the mapping to a natural order during distribution
5063: Level: beginner
5065: .seealso: DMGetUseNatural(), DMCreate()
5066: @*/
5067: PetscErrorCode DMSetUseNatural(DM dm, PetscBool useNatural)
5068: {
5072: dm->useNatural = useNatural;
5073: return(0);
5074: }
5077: /*@C
5078: DMCreateLabel - Create a label of the given name if it does not already exist
5080: Not Collective
5082: Input Parameters:
5083: + dm - The DM object
5084: - name - The label name
5086: Level: intermediate
5088: .keywords: mesh
5089: .seealso: DMLabelCreate(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5090: @*/
5091: PetscErrorCode DMCreateLabel(DM dm, const char name[])
5092: {
5093: DMLabelLink next = dm->labels->next;
5094: PetscBool flg = PETSC_FALSE;
5100: while (next) {
5101: PetscStrcmp(name, next->label->name, &flg);
5102: if (flg) break;
5103: next = next->next;
5104: }
5105: if (!flg) {
5106: DMLabelLink tmpLabel;
5108: PetscCalloc1(1, &tmpLabel);
5109: DMLabelCreate(name, &tmpLabel->label);
5110: tmpLabel->output = PETSC_TRUE;
5111: tmpLabel->next = dm->labels->next;
5112: dm->labels->next = tmpLabel;
5113: }
5114: return(0);
5115: }
5117: /*@C
5118: DMGetLabelValue - Get the value in a Sieve Label for the given point, with 0 as the default
5120: Not Collective
5122: Input Parameters:
5123: + dm - The DM object
5124: . name - The label name
5125: - point - The mesh point
5127: Output Parameter:
5128: . value - The label value for this point, or -1 if the point is not in the label
5130: Level: beginner
5132: .keywords: mesh
5133: .seealso: DMLabelGetValue(), DMSetLabelValue(), DMGetStratumIS()
5134: @*/
5135: PetscErrorCode DMGetLabelValue(DM dm, const char name[], PetscInt point, PetscInt *value)
5136: {
5137: DMLabel label;
5143: DMGetLabel(dm, name, &label);
5144: if (!label) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No label named %s was found", name);
5145: DMLabelGetValue(label, point, value);
5146: return(0);
5147: }
5149: /*@C
5150: DMSetLabelValue - Add a point to a Sieve Label with given value
5152: Not Collective
5154: Input Parameters:
5155: + dm - The DM object
5156: . name - The label name
5157: . point - The mesh point
5158: - value - The label value for this point
5160: Output Parameter:
5162: Level: beginner
5164: .keywords: mesh
5165: .seealso: DMLabelSetValue(), DMGetStratumIS(), DMClearLabelValue()
5166: @*/
5167: PetscErrorCode DMSetLabelValue(DM dm, const char name[], PetscInt point, PetscInt value)
5168: {
5169: DMLabel label;
5175: DMGetLabel(dm, name, &label);
5176: if (!label) {
5177: DMCreateLabel(dm, name);
5178: DMGetLabel(dm, name, &label);
5179: }
5180: DMLabelSetValue(label, point, value);
5181: return(0);
5182: }
5184: /*@C
5185: DMClearLabelValue - Remove a point from a Sieve Label with given value
5187: Not Collective
5189: Input Parameters:
5190: + dm - The DM object
5191: . name - The label name
5192: . point - The mesh point
5193: - value - The label value for this point
5195: Output Parameter:
5197: Level: beginner
5199: .keywords: mesh
5200: .seealso: DMLabelClearValue(), DMSetLabelValue(), DMGetStratumIS()
5201: @*/
5202: PetscErrorCode DMClearLabelValue(DM dm, const char name[], PetscInt point, PetscInt value)
5203: {
5204: DMLabel label;
5210: DMGetLabel(dm, name, &label);
5211: if (!label) return(0);
5212: DMLabelClearValue(label, point, value);
5213: return(0);
5214: }
5216: /*@C
5217: DMGetLabelSize - Get the number of different integer ids in a Label
5219: Not Collective
5221: Input Parameters:
5222: + dm - The DM object
5223: - name - The label name
5225: Output Parameter:
5226: . size - The number of different integer ids, or 0 if the label does not exist
5228: Level: beginner
5230: .keywords: mesh
5231: .seealso: DMLabeGetNumValues(), DMSetLabelValue()
5232: @*/
5233: PetscErrorCode DMGetLabelSize(DM dm, const char name[], PetscInt *size)
5234: {
5235: DMLabel label;
5242: DMGetLabel(dm, name, &label);
5243: *size = 0;
5244: if (!label) return(0);
5245: DMLabelGetNumValues(label, size);
5246: return(0);
5247: }
5249: /*@C
5250: DMGetLabelIdIS - Get the integer ids in a label
5252: Not Collective
5254: Input Parameters:
5255: + mesh - The DM object
5256: - name - The label name
5258: Output Parameter:
5259: . ids - The integer ids, or NULL if the label does not exist
5261: Level: beginner
5263: .keywords: mesh
5264: .seealso: DMLabelGetValueIS(), DMGetLabelSize()
5265: @*/
5266: PetscErrorCode DMGetLabelIdIS(DM dm, const char name[], IS *ids)
5267: {
5268: DMLabel label;
5275: DMGetLabel(dm, name, &label);
5276: *ids = NULL;
5277: if (!label) return(0);
5278: DMLabelGetValueIS(label, ids);
5279: return(0);
5280: }
5282: /*@C
5283: DMGetStratumSize - Get the number of points in a label stratum
5285: Not Collective
5287: Input Parameters:
5288: + dm - The DM object
5289: . name - The label name
5290: - value - The stratum value
5292: Output Parameter:
5293: . size - The stratum size
5295: Level: beginner
5297: .keywords: mesh
5298: .seealso: DMLabelGetStratumSize(), DMGetLabelSize(), DMGetLabelIds()
5299: @*/
5300: PetscErrorCode DMGetStratumSize(DM dm, const char name[], PetscInt value, PetscInt *size)
5301: {
5302: DMLabel label;
5309: DMGetLabel(dm, name, &label);
5310: *size = 0;
5311: if (!label) return(0);
5312: DMLabelGetStratumSize(label, value, size);
5313: return(0);
5314: }
5316: /*@C
5317: DMGetStratumIS - Get the points in a label stratum
5319: Not Collective
5321: Input Parameters:
5322: + dm - The DM object
5323: . name - The label name
5324: - value - The stratum value
5326: Output Parameter:
5327: . points - The stratum points, or NULL if the label does not exist or does not have that value
5329: Level: beginner
5331: .keywords: mesh
5332: .seealso: DMLabelGetStratumIS(), DMGetStratumSize()
5333: @*/
5334: PetscErrorCode DMGetStratumIS(DM dm, const char name[], PetscInt value, IS *points)
5335: {
5336: DMLabel label;
5343: DMGetLabel(dm, name, &label);
5344: *points = NULL;
5345: if (!label) return(0);
5346: DMLabelGetStratumIS(label, value, points);
5347: return(0);
5348: }
5350: /*@C
5351: DMGetStratumIS - Set the points in a label stratum
5353: Not Collective
5355: Input Parameters:
5356: + dm - The DM object
5357: . name - The label name
5358: . value - The stratum value
5359: - points - The stratum points
5361: Level: beginner
5363: .keywords: mesh
5364: .seealso: DMLabelSetStratumIS(), DMGetStratumSize()
5365: @*/
5366: PetscErrorCode DMSetStratumIS(DM dm, const char name[], PetscInt value, IS points)
5367: {
5368: DMLabel label;
5375: DMGetLabel(dm, name, &label);
5376: if (!label) return(0);
5377: DMLabelSetStratumIS(label, value, points);
5378: return(0);
5379: }
5381: /*@C
5382: DMClearLabelStratum - Remove all points from a stratum from a Sieve Label
5384: Not Collective
5386: Input Parameters:
5387: + dm - The DM object
5388: . name - The label name
5389: - value - The label value for this point
5391: Output Parameter:
5393: Level: beginner
5395: .keywords: mesh
5396: .seealso: DMLabelClearStratum(), DMSetLabelValue(), DMGetStratumIS(), DMClearLabelValue()
5397: @*/
5398: PetscErrorCode DMClearLabelStratum(DM dm, const char name[], PetscInt value)
5399: {
5400: DMLabel label;
5406: DMGetLabel(dm, name, &label);
5407: if (!label) return(0);
5408: DMLabelClearStratum(label, value);
5409: return(0);
5410: }
5412: /*@
5413: DMGetNumLabels - Return the number of labels defined by the mesh
5415: Not Collective
5417: Input Parameter:
5418: . dm - The DM object
5420: Output Parameter:
5421: . numLabels - the number of Labels
5423: Level: intermediate
5425: .keywords: mesh
5426: .seealso: DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5427: @*/
5428: PetscErrorCode DMGetNumLabels(DM dm, PetscInt *numLabels)
5429: {
5430: DMLabelLink next = dm->labels->next;
5431: PetscInt n = 0;
5436: while (next) {++n; next = next->next;}
5437: *numLabels = n;
5438: return(0);
5439: }
5441: /*@C
5442: DMGetLabelName - Return the name of nth label
5444: Not Collective
5446: Input Parameters:
5447: + dm - The DM object
5448: - n - the label number
5450: Output Parameter:
5451: . name - the label name
5453: Level: intermediate
5455: .keywords: mesh
5456: .seealso: DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5457: @*/
5458: PetscErrorCode DMGetLabelName(DM dm, PetscInt n, const char **name)
5459: {
5460: DMLabelLink next = dm->labels->next;
5461: PetscInt l = 0;
5466: while (next) {
5467: if (l == n) {
5468: *name = next->label->name;
5469: return(0);
5470: }
5471: ++l;
5472: next = next->next;
5473: }
5474: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label %D does not exist in this DM", n);
5475: }
5477: /*@C
5478: DMHasLabel - Determine whether the mesh has a label of a given name
5480: Not Collective
5482: Input Parameters:
5483: + dm - The DM object
5484: - name - The label name
5486: Output Parameter:
5487: . hasLabel - PETSC_TRUE if the label is present
5489: Level: intermediate
5491: .keywords: mesh
5492: .seealso: DMCreateLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5493: @*/
5494: PetscErrorCode DMHasLabel(DM dm, const char name[], PetscBool *hasLabel)
5495: {
5496: DMLabelLink next = dm->labels->next;
5503: *hasLabel = PETSC_FALSE;
5504: while (next) {
5505: PetscStrcmp(name, next->label->name, hasLabel);
5506: if (*hasLabel) break;
5507: next = next->next;
5508: }
5509: return(0);
5510: }
5512: /*@C
5513: DMGetLabel - Return the label of a given name, or NULL
5515: Not Collective
5517: Input Parameters:
5518: + dm - The DM object
5519: - name - The label name
5521: Output Parameter:
5522: . label - The DMLabel, or NULL if the label is absent
5524: Level: intermediate
5526: .keywords: mesh
5527: .seealso: DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5528: @*/
5529: PetscErrorCode DMGetLabel(DM dm, const char name[], DMLabel *label)
5530: {
5531: DMLabelLink next = dm->labels->next;
5532: PetscBool hasLabel;
5539: *label = NULL;
5540: while (next) {
5541: PetscStrcmp(name, next->label->name, &hasLabel);
5542: if (hasLabel) {
5543: *label = next->label;
5544: break;
5545: }
5546: next = next->next;
5547: }
5548: return(0);
5549: }
5551: /*@C
5552: DMGetLabelByNum - Return the nth label
5554: Not Collective
5556: Input Parameters:
5557: + dm - The DM object
5558: - n - the label number
5560: Output Parameter:
5561: . label - the label
5563: Level: intermediate
5565: .keywords: mesh
5566: .seealso: DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5567: @*/
5568: PetscErrorCode DMGetLabelByNum(DM dm, PetscInt n, DMLabel *label)
5569: {
5570: DMLabelLink next = dm->labels->next;
5571: PetscInt l = 0;
5576: while (next) {
5577: if (l == n) {
5578: *label = next->label;
5579: return(0);
5580: }
5581: ++l;
5582: next = next->next;
5583: }
5584: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label %D does not exist in this DM", n);
5585: }
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: }
5619: /*@C
5620: DMRemoveLabel - Remove the label from this mesh
5622: Not Collective
5624: Input Parameters:
5625: + dm - The DM object
5626: - name - The label name
5628: Output Parameter:
5629: . label - The DMLabel, or NULL if the label is absent
5631: Level: developer
5633: .keywords: mesh
5634: .seealso: DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5635: @*/
5636: PetscErrorCode DMRemoveLabel(DM dm, const char name[], DMLabel *label)
5637: {
5638: DMLabelLink next = dm->labels->next;
5639: DMLabelLink last = NULL;
5640: PetscBool hasLabel;
5645: DMHasLabel(dm, name, &hasLabel);
5646: *label = NULL;
5647: if (!hasLabel) return(0);
5648: while (next) {
5649: PetscStrcmp(name, next->label->name, &hasLabel);
5650: if (hasLabel) {
5651: if (last) last->next = next->next;
5652: else dm->labels->next = next->next;
5653: next->next = NULL;
5654: *label = next->label;
5655: PetscStrcmp(name, "depth", &hasLabel);
5656: if (hasLabel) {
5657: dm->depthLabel = NULL;
5658: }
5659: PetscFree(next);
5660: break;
5661: }
5662: last = next;
5663: next = next->next;
5664: }
5665: return(0);
5666: }
5668: /*@C
5669: DMGetLabelOutput - Get the output flag for a given label
5671: Not Collective
5673: Input Parameters:
5674: + dm - The DM object
5675: - name - The label name
5677: Output Parameter:
5678: . output - The flag for output
5680: Level: developer
5682: .keywords: mesh
5683: .seealso: DMSetLabelOutput(), DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5684: @*/
5685: PetscErrorCode DMGetLabelOutput(DM dm, const char name[], PetscBool *output)
5686: {
5687: DMLabelLink next = dm->labels->next;
5694: while (next) {
5695: PetscBool flg;
5697: PetscStrcmp(name, next->label->name, &flg);
5698: if (flg) {*output = next->output; return(0);}
5699: next = next->next;
5700: }
5701: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No label named %s was present in this dm", name);
5702: }
5704: /*@C
5705: DMSetLabelOutput - Set the output flag for a given label
5707: Not Collective
5709: Input Parameters:
5710: + dm - The DM object
5711: . name - The label name
5712: - output - The flag for output
5714: Level: developer
5716: .keywords: mesh
5717: .seealso: DMGetLabelOutput(), DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5718: @*/
5719: PetscErrorCode DMSetLabelOutput(DM dm, const char name[], PetscBool output)
5720: {
5721: DMLabelLink next = dm->labels->next;
5727: while (next) {
5728: PetscBool flg;
5730: PetscStrcmp(name, next->label->name, &flg);
5731: if (flg) {next->output = output; return(0);}
5732: next = next->next;
5733: }
5734: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No label named %s was present in this dm", name);
5735: }
5738: /*@
5739: DMCopyLabels - Copy labels from one mesh to another with a superset of the points
5741: Collective on DM
5743: Input Parameter:
5744: . dmA - The DM object with initial labels
5746: Output Parameter:
5747: . dmB - The DM object with copied labels
5749: Level: intermediate
5751: Note: This is typically used when interpolating or otherwise adding to a mesh
5753: .keywords: mesh
5754: .seealso: DMCopyCoordinates(), DMGetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM(), DMGetCoordinateSection()
5755: @*/
5756: PetscErrorCode DMCopyLabels(DM dmA, DM dmB)
5757: {
5758: PetscInt numLabels, l;
5762: if (dmA == dmB) return(0);
5763: DMGetNumLabels(dmA, &numLabels);
5764: for (l = 0; l < numLabels; ++l) {
5765: DMLabel label, labelNew;
5766: const char *name;
5767: PetscBool flg;
5769: DMGetLabelName(dmA, l, &name);
5770: PetscStrcmp(name, "depth", &flg);
5771: if (flg) continue;
5772: DMGetLabel(dmA, name, &label);
5773: DMLabelDuplicate(label, &labelNew);
5774: DMAddLabel(dmB, labelNew);
5775: }
5776: return(0);
5777: }
5779: /*@
5780: DMGetCoarseDM - Get the coarse mesh from which this was obtained by refinement
5782: Input Parameter:
5783: . dm - The DM object
5785: Output Parameter:
5786: . cdm - The coarse DM
5788: Level: intermediate
5790: .seealso: DMSetCoarseDM()
5791: @*/
5792: PetscErrorCode DMGetCoarseDM(DM dm, DM *cdm)
5793: {
5797: *cdm = dm->coarseMesh;
5798: return(0);
5799: }
5801: /*@
5802: DMSetCoarseDM - Set the coarse mesh from which this was obtained by refinement
5804: Input Parameters:
5805: + dm - The DM object
5806: - cdm - The coarse DM
5808: Level: intermediate
5810: .seealso: DMGetCoarseDM()
5811: @*/
5812: PetscErrorCode DMSetCoarseDM(DM dm, DM cdm)
5813: {
5819: PetscObjectReference((PetscObject)cdm);
5820: DMDestroy(&dm->coarseMesh);
5821: dm->coarseMesh = cdm;
5822: return(0);
5823: }
5825: /*@
5826: DMGetFineDM - Get the fine mesh from which this was obtained by refinement
5828: Input Parameter:
5829: . dm - The DM object
5831: Output Parameter:
5832: . fdm - The fine DM
5834: Level: intermediate
5836: .seealso: DMSetFineDM()
5837: @*/
5838: PetscErrorCode DMGetFineDM(DM dm, DM *fdm)
5839: {
5843: *fdm = dm->fineMesh;
5844: return(0);
5845: }
5847: /*@
5848: DMSetFineDM - Set the fine mesh from which this was obtained by refinement
5850: Input Parameters:
5851: + dm - The DM object
5852: - fdm - The fine DM
5854: Level: intermediate
5856: .seealso: DMGetFineDM()
5857: @*/
5858: PetscErrorCode DMSetFineDM(DM dm, DM fdm)
5859: {
5865: PetscObjectReference((PetscObject)fdm);
5866: DMDestroy(&dm->fineMesh);
5867: dm->fineMesh = fdm;
5868: return(0);
5869: }
5871: /*=== DMBoundary code ===*/
5873: PetscErrorCode DMCopyBoundary(DM dm, DM dmNew)
5874: {
5878: PetscDSCopyBoundary(dm->prob,dmNew->prob);
5879: return(0);
5880: }
5882: /*@C
5883: DMAddBoundary - Add a boundary condition to the model
5885: Input Parameters:
5886: + dm - The DM, with a PetscDS that matches the problem being constrained
5887: . type - The type of condition, e.g. DM_BC_ESSENTIAL_ANALYTIC/DM_BC_ESSENTIAL_FIELD (Dirichlet), or DM_BC_NATURAL (Neumann)
5888: . name - The BC name
5889: . labelname - The label defining constrained points
5890: . field - The field to constrain
5891: . numcomps - The number of constrained field components
5892: . comps - An array of constrained component numbers
5893: . bcFunc - A pointwise function giving boundary values
5894: . numids - The number of DMLabel ids for constrained points
5895: . ids - An array of ids for constrained points
5896: - ctx - An optional user context for bcFunc
5898: Options Database Keys:
5899: + -bc_<boundary name> <num> - Overrides the boundary ids
5900: - -bc_<boundary name>_comp <num> - Overrides the boundary components
5902: Level: developer
5904: .seealso: DMGetBoundary()
5905: @*/
5906: PetscErrorCode DMAddBoundary(DM dm, DMBoundaryConditionType type, const char name[], const char labelname[], PetscInt field, PetscInt numcomps, const PetscInt *comps, void (*bcFunc)(void), PetscInt numids, const PetscInt *ids, void *ctx)
5907: {
5912: PetscDSAddBoundary(dm->prob,type,name,labelname,field,numcomps,comps,bcFunc,numids,ids,ctx);
5913: return(0);
5914: }
5916: /*@
5917: DMGetNumBoundary - Get the number of registered BC
5919: Input Parameters:
5920: . dm - The mesh object
5922: Output Parameters:
5923: . numBd - The number of BC
5925: Level: intermediate
5927: .seealso: DMAddBoundary(), DMGetBoundary()
5928: @*/
5929: PetscErrorCode DMGetNumBoundary(DM dm, PetscInt *numBd)
5930: {
5935: PetscDSGetNumBoundary(dm->prob,numBd);
5936: return(0);
5937: }
5939: /*@C
5940: DMGetBoundary - Get a model boundary condition
5942: Input Parameters:
5943: + dm - The mesh object
5944: - bd - The BC number
5946: Output Parameters:
5947: + type - The type of condition, e.g. DM_BC_ESSENTIAL_ANALYTIC/DM_BC_ESSENTIAL_FIELD (Dirichlet), or DM_BC_NATURAL (Neumann)
5948: . name - The BC name
5949: . labelname - The label defining constrained points
5950: . field - The field to constrain
5951: . numcomps - The number of constrained field components
5952: . comps - An array of constrained component numbers
5953: . bcFunc - A pointwise function giving boundary values
5954: . numids - The number of DMLabel ids for constrained points
5955: . ids - An array of ids for constrained points
5956: - ctx - An optional user context for bcFunc
5958: Options Database Keys:
5959: + -bc_<boundary name> <num> - Overrides the boundary ids
5960: - -bc_<boundary name>_comp <num> - Overrides the boundary components
5962: Level: developer
5964: .seealso: DMAddBoundary()
5965: @*/
5966: PetscErrorCode DMGetBoundary(DM dm, PetscInt bd, DMBoundaryConditionType *type, const char **name, const char **labelname, PetscInt *field, PetscInt *numcomps, const PetscInt **comps, void (**func)(void), PetscInt *numids, const PetscInt **ids, void **ctx)
5967: {
5972: PetscDSGetBoundary(dm->prob,bd,type,name,labelname,field,numcomps,comps,func,numids,ids,ctx);
5973: return(0);
5974: }
5976: static PetscErrorCode DMPopulateBoundary(DM dm)
5977: {
5978: DMBoundary *lastnext;
5979: DSBoundary dsbound;
5983: dsbound = dm->prob->boundary;
5984: if (dm->boundary) {
5985: DMBoundary next = dm->boundary;
5987: /* quick check to see if the PetscDS has changed */
5988: if (next->dsboundary == dsbound) return(0);
5989: /* the PetscDS has changed: tear down and rebuild */
5990: while (next) {
5991: DMBoundary b = next;
5993: next = b->next;
5994: PetscFree(b);
5995: }
5996: dm->boundary = NULL;
5997: }
5999: lastnext = &(dm->boundary);
6000: while (dsbound) {
6001: DMBoundary dmbound;
6003: PetscNew(&dmbound);
6004: dmbound->dsboundary = dsbound;
6005: DMGetLabel(dm, dsbound->labelname, &(dmbound->label));
6006: if (!dmbound->label) PetscInfo2(dm, "DSBoundary %s wants label %s, which is not in this dm.\n",dsbound->name,dsbound->labelname);
6007: /* push on the back instead of the front so that it is in the same order as in the PetscDS */
6008: *lastnext = dmbound;
6009: lastnext = &(dmbound->next);
6010: dsbound = dsbound->next;
6011: }
6012: return(0);
6013: }
6015: PetscErrorCode DMIsBoundaryPoint(DM dm, PetscInt point, PetscBool *isBd)
6016: {
6017: DMBoundary b;
6023: *isBd = PETSC_FALSE;
6024: DMPopulateBoundary(dm);
6025: b = dm->boundary;
6026: while (b && !(*isBd)) {
6027: DMLabel label = b->label;
6028: DSBoundary dsb = b->dsboundary;
6030: if (label) {
6031: PetscInt i;
6033: for (i = 0; i < dsb->numids && !(*isBd); ++i) {
6034: DMLabelStratumHasPoint(label, dsb->ids[i], point, isBd);
6035: }
6036: }
6037: b = b->next;
6038: }
6039: return(0);
6040: }
6042: /*@C
6043: DMProjectFunction - This projects the given function into the function space provided.
6045: Input Parameters:
6046: + dm - The DM
6047: . time - The time
6048: . funcs - The coordinate functions to evaluate, one per field
6049: . ctxs - Optional array of contexts to pass to each coordinate function. ctxs itself may be null.
6050: - mode - The insertion mode for values
6052: Output Parameter:
6053: . X - vector
6055: Calling sequence of func:
6056: $ func(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar u[], void *ctx);
6058: + dim - The spatial dimension
6059: . x - The coordinates
6060: . Nf - The number of fields
6061: . u - The output field values
6062: - ctx - optional user-defined function context
6064: Level: developer
6066: .seealso: DMComputeL2Diff()
6067: @*/
6068: PetscErrorCode DMProjectFunction(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec X)
6069: {
6070: Vec localX;
6075: DMGetLocalVector(dm, &localX);
6076: DMProjectFunctionLocal(dm, time, funcs, ctxs, mode, localX);
6077: DMLocalToGlobalBegin(dm, localX, mode, X);
6078: DMLocalToGlobalEnd(dm, localX, mode, X);
6079: DMRestoreLocalVector(dm, &localX);
6080: return(0);
6081: }
6083: PetscErrorCode DMProjectFunctionLocal(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX)
6084: {
6090: if (!dm->ops->projectfunctionlocal) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implemnt DMProjectFunctionLocal",((PetscObject)dm)->type_name);
6091: (dm->ops->projectfunctionlocal) (dm, time, funcs, ctxs, mode, localX);
6092: return(0);
6093: }
6095: PetscErrorCode DMProjectFunctionLabelLocal(DM dm, PetscReal time, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscInt Nc, const PetscInt comps[], PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX)
6096: {
6102: if (!dm->ops->projectfunctionlabellocal) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implemnt DMProjectFunctionLabelLocal",((PetscObject)dm)->type_name);
6103: (dm->ops->projectfunctionlabellocal) (dm, time, label, numIds, ids, Nc, comps, funcs, ctxs, mode, localX);
6104: return(0);
6105: }
6107: PetscErrorCode DMProjectFieldLocal(DM dm, PetscReal time, Vec localU,
6108: void (**funcs)(PetscInt, PetscInt, PetscInt,
6109: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
6110: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
6111: PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]),
6112: InsertMode mode, Vec localX)
6113: {
6120: if (!dm->ops->projectfieldlocal) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implemnt DMProjectFieldLocal",((PetscObject)dm)->type_name);
6121: (dm->ops->projectfieldlocal) (dm, time, localU, funcs, mode, localX);
6122: return(0);
6123: }
6125: PetscErrorCode DMProjectFieldLabelLocal(DM dm, PetscReal time, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscInt Nc, const PetscInt comps[], Vec localU,
6126: void (**funcs)(PetscInt, PetscInt, PetscInt,
6127: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
6128: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
6129: PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]),
6130: InsertMode mode, Vec localX)
6131: {
6138: if (!dm->ops->projectfieldlabellocal) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implemnt DMProjectFieldLocal",((PetscObject)dm)->type_name);
6139: (dm->ops->projectfieldlabellocal)(dm, time, label, numIds, ids, Nc, comps, localU, funcs, mode, localX);
6140: return(0);
6141: }
6143: /*@C
6144: DMComputeL2Diff - This function computes the L_2 difference between a function u and an FEM interpolant solution u_h.
6146: Input Parameters:
6147: + dm - The DM
6148: . time - The time
6149: . funcs - The functions to evaluate for each field component
6150: . ctxs - Optional array of contexts to pass to each function, or NULL.
6151: - X - The coefficient vector u_h
6153: Output Parameter:
6154: . diff - The diff ||u - u_h||_2
6156: Level: developer
6158: .seealso: DMProjectFunction(), DMComputeL2FieldDiff(), DMComputeL2GradientDiff()
6159: @*/
6160: PetscErrorCode DMComputeL2Diff(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff)
6161: {
6167: if (!dm->ops->computel2diff) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implemnt DMComputeL2Diff",((PetscObject)dm)->type_name);
6168: (dm->ops->computel2diff)(dm,time,funcs,ctxs,X,diff);
6169: return(0);
6170: }
6172: /*@C
6173: DMComputeL2GradientDiff - This function computes the L_2 difference between the gradient of a function u and an FEM interpolant solution grad u_h.
6175: Input Parameters:
6176: + dm - The DM
6177: , time - The time
6178: . funcs - The gradient functions to evaluate for each field component
6179: . ctxs - Optional array of contexts to pass to each function, or NULL.
6180: . X - The coefficient vector u_h
6181: - n - The vector to project along
6183: Output Parameter:
6184: . diff - The diff ||(grad u - grad u_h) . n||_2
6186: Level: developer
6188: .seealso: DMProjectFunction(), DMComputeL2Diff()
6189: @*/
6190: 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)
6191: {
6197: if (!dm->ops->computel2gradientdiff) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMComputeL2GradientDiff",((PetscObject)dm)->type_name);
6198: (dm->ops->computel2gradientdiff)(dm,time,funcs,ctxs,X,n,diff);
6199: return(0);
6200: }
6202: /*@C
6203: DMComputeL2FieldDiff - This function computes the L_2 difference between a function u and an FEM interpolant solution u_h, separated into field components.
6205: Input Parameters:
6206: + dm - The DM
6207: . time - The time
6208: . funcs - The functions to evaluate for each field component
6209: . ctxs - Optional array of contexts to pass to each function, or NULL.
6210: - X - The coefficient vector u_h
6212: Output Parameter:
6213: . diff - The array of differences, ||u^f - u^f_h||_2
6215: Level: developer
6217: .seealso: DMProjectFunction(), DMComputeL2FieldDiff(), DMComputeL2GradientDiff()
6218: @*/
6219: PetscErrorCode DMComputeL2FieldDiff(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal diff[])
6220: {
6226: if (!dm->ops->computel2fielddiff) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implemnt DMComputeL2FieldDiff",((PetscObject)dm)->type_name);
6227: (dm->ops->computel2fielddiff)(dm,time,funcs,ctxs,X,diff);
6228: return(0);
6229: }
6231: /*@C
6232: DMAdaptLabel - Adapt a dm based on a label with values interpreted as coarsening and refining flags. Specific implementations of DM maybe have
6233: specialized flags, but all implementations should accept flag values DM_ADAPT_DETERMINE, DM_ADAPT_KEEP, DM_ADAPT_REFINE, and DM_ADAPT_COARSEN.
6235: Collective on dm
6237: Input parameters:
6238: + dm - the pre-adaptation DM object
6239: - label - label with the flags
6241: Output parameters:
6242: . dmAdapt - the adapted DM object: may be NULL if an adapted DM could not be produced.
6244: Level: intermediate
6246: .seealso: DMAdaptMetric(), DMCoarsen(), DMRefine()
6247: @*/
6248: PetscErrorCode DMAdaptLabel(DM dm, DMLabel label, DM *dmAdapt)
6249: {
6256: *dmAdapt = NULL;
6257: if (!dm->ops->adaptlabel) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMAdaptLabel",((PetscObject)dm)->type_name);
6258: (dm->ops->adaptlabel)(dm, label, dmAdapt);
6259: return(0);
6260: }
6262: /*@C
6263: DMAdaptMetric - Generates a mesh adapted to the specified metric field using the pragmatic library.
6265: Input Parameters:
6266: + dm - The DM object
6267: . metric - The metric to which the mesh is adapted, defined vertex-wise.
6268: - bdLabel - Label for boundary tags, which will be preserved in the output mesh. bdLabel should be NULL if there is no such label, and should be different from "_boundary_".
6270: Output Parameter:
6271: . dmAdapt - Pointer to the DM object containing the adapted mesh
6273: Note: The label in the adapted mesh will be registered under the name of the input DMLabel object
6275: Level: advanced
6277: .seealso: DMAdaptLabel(), DMCoarsen(), DMRefine()
6278: @*/
6279: PetscErrorCode DMAdaptMetric(DM dm, Vec metric, DMLabel bdLabel, DM *dmAdapt)
6280: {
6288: *dmAdapt = NULL;
6289: if (!dm->ops->adaptmetric) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMAdaptMetric",((PetscObject)dm)->type_name);
6290: (dm->ops->adaptmetric)(dm, metric, bdLabel, dmAdapt);
6291: return(0);
6292: }
6294: /*@C
6295: DMGetNeighbors - Gets an array containing the MPI rank of all the processes neighbors
6297: Not Collective
6299: Input Parameter:
6300: . dm - The DM
6302: Output Parameter:
6303: . nranks - the number of neighbours
6304: . ranks - the neighbors ranks
6306: Notes:
6307: Do not free the array, it is freed when the DM is destroyed.
6309: Level: beginner
6311: .seealso: DMDAGetNeighbors(), PetscSFGetRanks()
6312: @*/
6313: PetscErrorCode DMGetNeighbors(DM dm,PetscInt *nranks,const PetscMPIInt *ranks[])
6314: {
6319: if (!dm->ops->getneighbors) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implemnt DMGetNeighbors",((PetscObject)dm)->type_name);
6320: (dm->ops->getneighbors)(dm,nranks,ranks);
6321: return(0);
6322: }
6324: #include <petsc/private/matimpl.h> /* Needed because of coloring->ctype below */
6326: /*
6327: Converts the input vector to a ghosted vector and then calls the standard coloring code.
6328: This has be a different function because it requires DM which is not defined in the Mat library
6329: */
6330: PetscErrorCode MatFDColoringApply_AIJDM(Mat J,MatFDColoring coloring,Vec x1,void *sctx)
6331: {
6335: if (coloring->ctype == IS_COLORING_LOCAL) {
6336: Vec x1local;
6337: DM dm;
6338: MatGetDM(J,&dm);
6339: if (!dm) SETERRQ(PetscObjectComm((PetscObject)J),PETSC_ERR_ARG_INCOMP,"IS_COLORING_LOCAL requires a DM");
6340: DMGetLocalVector(dm,&x1local);
6341: DMGlobalToLocalBegin(dm,x1,INSERT_VALUES,x1local);
6342: DMGlobalToLocalEnd(dm,x1,INSERT_VALUES,x1local);
6343: x1 = x1local;
6344: }
6345: MatFDColoringApply_AIJ(J,coloring,x1,sctx);
6346: if (coloring->ctype == IS_COLORING_LOCAL) {
6347: DM dm;
6348: MatGetDM(J,&dm);
6349: DMRestoreLocalVector(dm,&x1);
6350: }
6351: return(0);
6352: }
6354: /*@
6355: MatFDColoringUseDM - allows a MatFDColoring object to use the DM associated with the matrix to use a IS_COLORING_LOCAL coloring
6357: Input Parameter:
6358: . coloring - the MatFDColoring object
6360: Developer Notes: this routine exists because the PETSc Mat library does not know about the DM objects
6362: Level: advanced
6364: .seealso: MatFDColoring, MatFDColoringCreate(), ISColoringType
6365: @*/
6366: PetscErrorCode MatFDColoringUseDM(Mat coloring,MatFDColoring fdcoloring)
6367: {
6369: coloring->ops->fdcoloringapply = MatFDColoringApply_AIJDM;
6370: return(0);
6371: }