Actual source code: dm.c
petsc-3.5.4 2015-05-23
1: #include <petsc-private/dmimpl.h> /*I "petscdm.h" I*/
2: #include <petscsf.h>
3: #include <petscds.h>
5: PetscClassId DM_CLASSID;
6: PetscLogEvent DM_Convert, DM_GlobalToLocal, DM_LocalToGlobal, DM_LocalToLocal;
8: const char *const DMBoundaryTypes[] = {"NONE","GHOSTED","MIRROR","PERIODIC","TWIST","DM_BOUNDARY_",0};
12: /*@
13: DMCreate - Creates an empty DM object. The type can then be set with DMSetType().
15: If you never call DMSetType() it will generate an
16: error when you try to use the vector.
18: Collective on MPI_Comm
20: Input Parameter:
21: . comm - The communicator for the DM object
23: Output Parameter:
24: . dm - The DM object
26: Level: beginner
28: .seealso: DMSetType(), DMDA, DMSLICED, DMCOMPOSITE
29: @*/
30: PetscErrorCode DMCreate(MPI_Comm comm,DM *dm)
31: {
32: DM v;
37: *dm = NULL;
38: PetscSysInitializePackage();
39: VecInitializePackage();
40: MatInitializePackage();
41: DMInitializePackage();
43: PetscHeaderCreate(v, _p_DM, struct _DMOps, DM_CLASSID, "DM", "Distribution Manager", "DM", comm, DMDestroy, DMView);
44: PetscMemzero(v->ops, sizeof(struct _DMOps));
47: v->ltogmap = NULL;
48: v->bs = 1;
49: v->coloringtype = IS_COLORING_GLOBAL;
50: PetscSFCreate(comm, &v->sf);
51: PetscSFCreate(comm, &v->defaultSF);
52: v->defaultSection = NULL;
53: v->defaultGlobalSection = NULL;
54: v->L = NULL;
55: v->maxCell = NULL;
56: {
57: PetscInt i;
58: for (i = 0; i < 10; ++i) {
59: v->nullspaceConstructors[i] = NULL;
60: }
61: }
62: PetscDSCreate(comm, &v->prob);
63: v->dmBC = NULL;
64: v->outputSequenceNum = -1;
65: v->outputSequenceVal = 0.0;
66: DMSetVecType(v,VECSTANDARD);
67: DMSetMatType(v,MATAIJ);
68: *dm = v;
69: return(0);
70: }
74: /*@
75: DMClone - Creates a DM object with the same topology as the original.
77: Collective on MPI_Comm
79: Input Parameter:
80: . dm - The original DM object
82: Output Parameter:
83: . newdm - The new DM object
85: Level: beginner
87: .keywords: DM, topology, create
88: @*/
89: PetscErrorCode DMClone(DM dm, DM *newdm)
90: {
91: PetscSF sf;
92: Vec coords;
93: void *ctx;
99: DMCreate(PetscObjectComm((PetscObject)dm), newdm);
100: if (dm->ops->clone) {
101: (*dm->ops->clone)(dm, newdm);
102: }
103: (*newdm)->setupcalled = PETSC_TRUE;
104: DMGetPointSF(dm, &sf);
105: DMSetPointSF(*newdm, sf);
106: DMGetApplicationContext(dm, &ctx);
107: DMSetApplicationContext(*newdm, ctx);
108: DMGetCoordinatesLocal(dm, &coords);
109: if (coords) {
110: DMSetCoordinatesLocal(*newdm, coords);
111: } else {
112: DMGetCoordinates(dm, &coords);
113: if (coords) {DMSetCoordinates(*newdm, coords);}
114: }
115: if (dm->maxCell) {
116: const PetscReal *maxCell, *L;
117: DMGetPeriodicity(dm, &maxCell, &L);
118: DMSetPeriodicity(*newdm, maxCell, L);
119: }
120: return(0);
121: }
125: /*@C
126: DMSetVecType - Sets the type of vector created with DMCreateLocalVector() and DMCreateGlobalVector()
128: Logically Collective on DMDA
130: Input Parameter:
131: + da - initial distributed array
132: . ctype - the vector type, currently either VECSTANDARD or VECCUSP
134: Options Database:
135: . -dm_vec_type ctype
137: Level: intermediate
139: .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMDestroy(), DMDA, DMDAInterpolationType, VecType, DMGetVecType()
140: @*/
141: PetscErrorCode DMSetVecType(DM da,VecType ctype)
142: {
147: PetscFree(da->vectype);
148: PetscStrallocpy(ctype,(char**)&da->vectype);
149: return(0);
150: }
154: /*@C
155: DMGetVecType - Gets the type of vector created with DMCreateLocalVector() and DMCreateGlobalVector()
157: Logically Collective on DMDA
159: Input Parameter:
160: . da - initial distributed array
162: Output Parameter:
163: . ctype - the vector type
165: Level: intermediate
167: .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMDestroy(), DMDA, DMDAInterpolationType, VecType
168: @*/
169: PetscErrorCode DMGetVecType(DM da,VecType *ctype)
170: {
173: *ctype = da->vectype;
174: return(0);
175: }
179: /*@
180: VecGetDM - Gets the DM defining the data layout of the vector
182: Not collective
184: Input Parameter:
185: . v - The Vec
187: Output Parameter:
188: . dm - The DM
190: Level: intermediate
192: .seealso: VecSetDM(), DMGetLocalVector(), DMGetGlobalVector(), DMSetVecType()
193: @*/
194: PetscErrorCode VecGetDM(Vec v, DM *dm)
195: {
201: PetscObjectQuery((PetscObject) v, "__PETSc_dm", (PetscObject*) dm);
202: return(0);
203: }
207: /*@
208: VecSetDM - Sets the DM defining the data layout of the vector
210: Not collective
212: Input Parameters:
213: + v - The Vec
214: - dm - The DM
216: Level: intermediate
218: .seealso: VecGetDM(), DMGetLocalVector(), DMGetGlobalVector(), DMSetVecType()
219: @*/
220: PetscErrorCode VecSetDM(Vec v, DM dm)
221: {
227: PetscObjectCompose((PetscObject) v, "__PETSc_dm", (PetscObject) dm);
228: return(0);
229: }
233: /*@C
234: DMSetMatType - Sets the type of matrix created with DMCreateMatrix()
236: Logically Collective on DM
238: Input Parameter:
239: + dm - the DM context
240: . ctype - the matrix type
242: Options Database:
243: . -dm_mat_type ctype
245: Level: intermediate
247: .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMCreateMatrix(), DMSetMatrixPreallocateOnly(), MatType, DMGetMatType()
248: @*/
249: PetscErrorCode DMSetMatType(DM dm,MatType ctype)
250: {
255: PetscFree(dm->mattype);
256: PetscStrallocpy(ctype,(char**)&dm->mattype);
257: return(0);
258: }
262: /*@C
263: DMGetMatType - Gets the type of matrix created with DMCreateMatrix()
265: Logically Collective on DM
267: Input Parameter:
268: . dm - the DM context
270: Output Parameter:
271: . ctype - the matrix type
273: Options Database:
274: . -dm_mat_type ctype
276: Level: intermediate
278: .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMCreateMatrix(), DMSetMatrixPreallocateOnly(), MatType, DMSetMatType()
279: @*/
280: PetscErrorCode DMGetMatType(DM dm,MatType *ctype)
281: {
284: *ctype = dm->mattype;
285: return(0);
286: }
290: /*@
291: MatGetDM - Gets the DM defining the data layout of the matrix
293: Not collective
295: Input Parameter:
296: . A - The Mat
298: Output Parameter:
299: . dm - The DM
301: Level: intermediate
303: .seealso: MatSetDM(), DMCreateMatrix(), DMSetMatType()
304: @*/
305: PetscErrorCode MatGetDM(Mat A, DM *dm)
306: {
312: PetscObjectQuery((PetscObject) A, "__PETSc_dm", (PetscObject*) dm);
313: return(0);
314: }
318: /*@
319: MatSetDM - Sets the DM defining the data layout of the matrix
321: Not collective
323: Input Parameters:
324: + A - The Mat
325: - dm - The DM
327: Level: intermediate
329: .seealso: MatGetDM(), DMCreateMatrix(), DMSetMatType()
330: @*/
331: PetscErrorCode MatSetDM(Mat A, DM dm)
332: {
338: PetscObjectCompose((PetscObject) A, "__PETSc_dm", (PetscObject) dm);
339: return(0);
340: }
344: /*@C
345: DMSetOptionsPrefix - Sets the prefix used for searching for all
346: DMDA options in the database.
348: Logically Collective on DMDA
350: Input Parameter:
351: + da - the DMDA context
352: - prefix - the prefix to prepend to all option names
354: Notes:
355: A hyphen (-) must NOT be given at the beginning of the prefix name.
356: The first character of all runtime options is AUTOMATICALLY the hyphen.
358: Level: advanced
360: .keywords: DMDA, set, options, prefix, database
362: .seealso: DMSetFromOptions()
363: @*/
364: PetscErrorCode DMSetOptionsPrefix(DM dm,const char prefix[])
365: {
370: PetscObjectSetOptionsPrefix((PetscObject)dm,prefix);
371: return(0);
372: }
376: /*@
377: DMDestroy - Destroys a vector packer or DMDA.
379: Collective on DM
381: Input Parameter:
382: . dm - the DM object to destroy
384: Level: developer
386: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix()
388: @*/
389: PetscErrorCode DMDestroy(DM *dm)
390: {
391: PetscInt i, cnt = 0, Nf = 0, f;
392: DMNamedVecLink nlink,nnext;
396: if (!*dm) return(0);
399: if ((*dm)->prob) {
400: PetscObject disc;
402: /* I think it makes sense to dump all attached things when you are destroyed, which also eliminates circular references */
403: PetscDSGetNumFields((*dm)->prob, &Nf);
404: for (f = 0; f < Nf; ++f) {
405: PetscDSGetDiscretization((*dm)->prob, f, &disc);
406: PetscObjectCompose(disc, "pmat", NULL);
407: PetscObjectCompose(disc, "nullspace", NULL);
408: PetscObjectCompose(disc, "nearnullspace", NULL);
409: }
410: }
411: /* count all the circular references of DM and its contained Vecs */
412: for (i=0; i<DM_MAX_WORK_VECTORS; i++) {
413: if ((*dm)->localin[i]) cnt++;
414: if ((*dm)->globalin[i]) cnt++;
415: }
416: for (nlink=(*dm)->namedglobal; nlink; nlink=nlink->next) cnt++;
417: for (nlink=(*dm)->namedlocal; nlink; nlink=nlink->next) cnt++;
418: if ((*dm)->x) {
419: DM obj;
420: VecGetDM((*dm)->x, &obj);
421: if (obj == *dm) cnt++;
422: }
424: if (--((PetscObject)(*dm))->refct - cnt > 0) {*dm = 0; return(0);}
425: /*
426: Need this test because the dm references the vectors that
427: reference the dm, so destroying the dm calls destroy on the
428: vectors that cause another destroy on the dm
429: */
430: if (((PetscObject)(*dm))->refct < 0) return(0);
431: ((PetscObject) (*dm))->refct = 0;
432: for (i=0; i<DM_MAX_WORK_VECTORS; i++) {
433: if ((*dm)->localout[i]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Destroying a DM that has a local vector obtained with DMGetLocalVector()");
434: VecDestroy(&(*dm)->localin[i]);
435: }
436: nnext=(*dm)->namedglobal;
437: (*dm)->namedglobal = NULL;
438: for (nlink=nnext; nlink; nlink=nnext) { /* Destroy the named vectors */
439: nnext = nlink->next;
440: if (nlink->status != DMVEC_STATUS_IN) SETERRQ1(((PetscObject)*dm)->comm,PETSC_ERR_ARG_WRONGSTATE,"DM still has Vec named '%s' checked out",nlink->name);
441: PetscFree(nlink->name);
442: VecDestroy(&nlink->X);
443: PetscFree(nlink);
444: }
445: nnext=(*dm)->namedlocal;
446: (*dm)->namedlocal = NULL;
447: for (nlink=nnext; nlink; nlink=nnext) { /* Destroy the named local vectors */
448: nnext = nlink->next;
449: if (nlink->status != DMVEC_STATUS_IN) SETERRQ1(((PetscObject)*dm)->comm,PETSC_ERR_ARG_WRONGSTATE,"DM still has Vec named '%s' checked out",nlink->name);
450: PetscFree(nlink->name);
451: VecDestroy(&nlink->X);
452: PetscFree(nlink);
453: }
455: /* Destroy the list of hooks */
456: {
457: DMCoarsenHookLink link,next;
458: for (link=(*dm)->coarsenhook; link; link=next) {
459: next = link->next;
460: PetscFree(link);
461: }
462: (*dm)->coarsenhook = NULL;
463: }
464: {
465: DMRefineHookLink link,next;
466: for (link=(*dm)->refinehook; link; link=next) {
467: next = link->next;
468: PetscFree(link);
469: }
470: (*dm)->refinehook = NULL;
471: }
472: {
473: DMSubDomainHookLink link,next;
474: for (link=(*dm)->subdomainhook; link; link=next) {
475: next = link->next;
476: PetscFree(link);
477: }
478: (*dm)->subdomainhook = NULL;
479: }
480: {
481: DMGlobalToLocalHookLink link,next;
482: for (link=(*dm)->gtolhook; link; link=next) {
483: next = link->next;
484: PetscFree(link);
485: }
486: (*dm)->gtolhook = NULL;
487: }
488: /* Destroy the work arrays */
489: {
490: DMWorkLink link,next;
491: if ((*dm)->workout) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Work array still checked out");
492: for (link=(*dm)->workin; link; link=next) {
493: next = link->next;
494: PetscFree(link->mem);
495: PetscFree(link);
496: }
497: (*dm)->workin = NULL;
498: }
500: PetscObjectDestroy(&(*dm)->dmksp);
501: PetscObjectDestroy(&(*dm)->dmsnes);
502: PetscObjectDestroy(&(*dm)->dmts);
504: if ((*dm)->ctx && (*dm)->ctxdestroy) {
505: (*(*dm)->ctxdestroy)(&(*dm)->ctx);
506: }
507: VecDestroy(&(*dm)->x);
508: MatFDColoringDestroy(&(*dm)->fd);
509: DMClearGlobalVectors(*dm);
510: ISLocalToGlobalMappingDestroy(&(*dm)->ltogmap);
511: PetscFree((*dm)->vectype);
512: PetscFree((*dm)->mattype);
514: PetscSectionDestroy(&(*dm)->defaultSection);
515: PetscSectionDestroy(&(*dm)->defaultGlobalSection);
516: PetscLayoutDestroy(&(*dm)->map);
517: PetscSFDestroy(&(*dm)->sf);
518: PetscSFDestroy(&(*dm)->defaultSF);
520: DMDestroy(&(*dm)->coordinateDM);
521: VecDestroy(&(*dm)->coordinates);
522: VecDestroy(&(*dm)->coordinatesLocal);
523: PetscFree((*dm)->maxCell);
524: PetscFree((*dm)->L);
526: PetscDSDestroy(&(*dm)->prob);
527: DMDestroy(&(*dm)->dmBC);
528: /* if memory was published with SAWs then destroy it */
529: PetscObjectSAWsViewOff((PetscObject)*dm);
531: (*(*dm)->ops->destroy)(*dm);
532: /* We do not destroy (*dm)->data here so that we can reference count backend objects */
533: PetscHeaderDestroy(dm);
534: return(0);
535: }
539: /*@
540: DMSetUp - sets up the data structures inside a DM object
542: Collective on DM
544: Input Parameter:
545: . dm - the DM object to setup
547: Level: developer
549: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix()
551: @*/
552: PetscErrorCode DMSetUp(DM dm)
553: {
558: if (dm->setupcalled) return(0);
559: if (dm->ops->setup) {
560: (*dm->ops->setup)(dm);
561: }
562: dm->setupcalled = PETSC_TRUE;
563: return(0);
564: }
568: /*@
569: DMSetFromOptions - sets parameters in a DM from the options database
571: Collective on DM
573: Input Parameter:
574: . dm - the DM object to set options for
576: Options Database:
577: + -dm_preallocate_only: Only preallocate the matrix for DMCreateMatrix(), but do not fill it with zeros
578: . -dm_vec_type <type> type of vector to create inside DM
579: . -dm_mat_type <type> type of matrix to create inside DM
580: - -dm_coloring_type <global or ghosted>
582: Level: developer
584: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix()
586: @*/
587: PetscErrorCode DMSetFromOptions(DM dm)
588: {
589: char typeName[256];
590: PetscBool flg;
595: PetscObjectOptionsBegin((PetscObject)dm);
596: PetscOptionsBool("-dm_preallocate_only","only preallocate matrix, but do not set column indices","DMSetMatrixPreallocateOnly",dm->prealloc_only,&dm->prealloc_only,NULL);
597: PetscOptionsFList("-dm_vec_type","Vector type used for created vectors","DMSetVecType",VecList,dm->vectype,typeName,256,&flg);
598: if (flg) {
599: DMSetVecType(dm,typeName);
600: }
601: PetscOptionsFList("-dm_mat_type","Matrix type used for created matrices","DMSetMatType",MatList,dm->mattype ? dm->mattype : typeName,typeName,sizeof(typeName),&flg);
602: if (flg) {
603: DMSetMatType(dm,typeName);
604: }
605: PetscOptionsEnum("-dm_is_coloring_type","Global or local coloring of Jacobian","ISColoringType",ISColoringTypes,(PetscEnum)dm->coloringtype,(PetscEnum*)&dm->coloringtype,NULL);
606: if (dm->ops->setfromoptions) {
607: (*dm->ops->setfromoptions)(dm);
608: }
609: /* process any options handlers added with PetscObjectAddOptionsHandler() */
610: PetscObjectProcessOptionsHandlers((PetscObject) dm);
611: PetscOptionsEnd();
612: return(0);
613: }
617: /*@C
618: DMView - Views a vector packer or DMDA.
620: Collective on DM
622: Input Parameter:
623: + dm - the DM object to view
624: - v - the viewer
626: Level: developer
628: .seealso DMDestroy(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix()
630: @*/
631: PetscErrorCode DMView(DM dm,PetscViewer v)
632: {
634: PetscBool isbinary;
638: if (!v) {
639: PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)dm),&v);
640: }
641: PetscObjectPrintClassNamePrefixType((PetscObject)dm,v);
642: PetscObjectTypeCompare((PetscObject)v,PETSCVIEWERBINARY,&isbinary);
643: if (isbinary) {
644: PetscInt classid = DM_FILE_CLASSID;
645: char type[256];
647: PetscViewerBinaryWrite(v,&classid,1,PETSC_INT,PETSC_FALSE);
648: PetscStrncpy(type,((PetscObject)dm)->type_name,256);
649: PetscViewerBinaryWrite(v,type,256,PETSC_CHAR,PETSC_FALSE);
650: }
651: if (dm->ops->view) {
652: (*dm->ops->view)(dm,v);
653: }
654: return(0);
655: }
659: /*@
660: DMCreateGlobalVector - Creates a global vector from a DMDA or DMComposite object
662: Collective on DM
664: Input Parameter:
665: . dm - the DM object
667: Output Parameter:
668: . vec - the global vector
670: Level: beginner
672: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix()
674: @*/
675: PetscErrorCode DMCreateGlobalVector(DM dm,Vec *vec)
676: {
681: (*dm->ops->createglobalvector)(dm,vec);
682: return(0);
683: }
687: /*@
688: DMCreateLocalVector - Creates a local vector from a DMDA or DMComposite object
690: Not Collective
692: Input Parameter:
693: . dm - the DM object
695: Output Parameter:
696: . vec - the local vector
698: Level: beginner
700: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix()
702: @*/
703: PetscErrorCode DMCreateLocalVector(DM dm,Vec *vec)
704: {
709: (*dm->ops->createlocalvector)(dm,vec);
710: return(0);
711: }
715: /*@
716: DMGetLocalToGlobalMapping - Accesses the local-to-global mapping in a DM.
718: Collective on DM
720: Input Parameter:
721: . dm - the DM that provides the mapping
723: Output Parameter:
724: . ltog - the mapping
726: Level: intermediate
728: Notes:
729: This mapping can then be used by VecSetLocalToGlobalMapping() or
730: MatSetLocalToGlobalMapping().
732: .seealso: DMCreateLocalVector()
733: @*/
734: PetscErrorCode DMGetLocalToGlobalMapping(DM dm,ISLocalToGlobalMapping *ltog)
735: {
741: if (!dm->ltogmap) {
742: PetscSection section, sectionGlobal;
744: DMGetDefaultSection(dm, §ion);
745: if (section) {
746: PetscInt *ltog;
747: PetscInt pStart, pEnd, size, p, l;
749: DMGetDefaultGlobalSection(dm, §ionGlobal);
750: PetscSectionGetChart(section, &pStart, &pEnd);
751: PetscSectionGetStorageSize(section, &size);
752: PetscMalloc1(size, <og); /* We want the local+overlap size */
753: for (p = pStart, l = 0; p < pEnd; ++p) {
754: PetscInt dof, off, c;
756: /* Should probably use constrained dofs */
757: PetscSectionGetDof(section, p, &dof);
758: PetscSectionGetOffset(sectionGlobal, p, &off);
759: for (c = 0; c < dof; ++c, ++l) {
760: ltog[l] = off+c;
761: }
762: }
763: ISLocalToGlobalMappingCreate(PETSC_COMM_SELF, 1,size, ltog, PETSC_OWN_POINTER, &dm->ltogmap);
764: PetscLogObjectParent((PetscObject)dm, (PetscObject)dm->ltogmap);
765: } else {
766: if (!dm->ops->getlocaltoglobalmapping) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM can not create LocalToGlobalMapping");
767: (*dm->ops->getlocaltoglobalmapping)(dm);
768: }
769: }
770: *ltog = dm->ltogmap;
771: return(0);
772: }
776: /*@
777: DMGetBlockSize - Gets the inherent block size associated with a DM
779: Not Collective
781: Input Parameter:
782: . dm - the DM with block structure
784: Output Parameter:
785: . bs - the block size, 1 implies no exploitable block structure
787: Level: intermediate
789: .seealso: ISCreateBlock(), VecSetBlockSize(), MatSetBlockSize(), DMGetLocalToGlobalMapping()
790: @*/
791: PetscErrorCode DMGetBlockSize(DM dm,PetscInt *bs)
792: {
796: if (dm->bs < 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"DM does not have enough information to provide a block size yet");
797: *bs = dm->bs;
798: return(0);
799: }
803: /*@
804: DMCreateInterpolation - Gets interpolation matrix between two DMDA or DMComposite objects
806: Collective on DM
808: Input Parameter:
809: + dm1 - the DM object
810: - dm2 - the second, finer DM object
812: Output Parameter:
813: + mat - the interpolation
814: - vec - the scaling (optional)
816: Level: developer
818: 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
819: DMCoarsen(). The coordinates set into the DMDA are completely ignored in computing the interpolation.
821: For DMDA objects you can use this interpolation (more precisely the interpolation from the DMGetCoordinateDM()) to interpolate the mesh coordinate vectors
822: EXCEPT in the periodic case where it does not make sense since the coordinate vectors are not periodic.
825: .seealso DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateColoring(), DMCreateMatrix(), DMRefine(), DMCoarsen()
827: @*/
828: PetscErrorCode DMCreateInterpolation(DM dm1,DM dm2,Mat *mat,Vec *vec)
829: {
835: (*dm1->ops->createinterpolation)(dm1,dm2,mat,vec);
836: return(0);
837: }
841: /*@
842: DMCreateInjection - Gets injection matrix between two DMDA or DMComposite objects
844: Collective on DM
846: Input Parameter:
847: + dm1 - the DM object
848: - dm2 - the second, finer DM object
850: Output Parameter:
851: . ctx - the injection
853: Level: developer
855: 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
856: DMCoarsen(). The coordinates set into the DMDA are completely ignored in computing the injection.
858: .seealso DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateColoring(), DMCreateMatrix(), DMCreateInterpolation()
860: @*/
861: PetscErrorCode DMCreateInjection(DM dm1,DM dm2,VecScatter *ctx)
862: {
868: (*dm1->ops->getinjection)(dm1,dm2,ctx);
869: return(0);
870: }
874: /*@
875: DMCreateColoring - Gets coloring for a DM
877: Collective on DM
879: Input Parameter:
880: + dm - the DM object
881: - ctype - IS_COLORING_GHOSTED or IS_COLORING_GLOBAL
883: Output Parameter:
884: . coloring - the coloring
886: Level: developer
888: .seealso DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateMatrix(), DMSetMatType()
890: @*/
891: PetscErrorCode DMCreateColoring(DM dm,ISColoringType ctype,ISColoring *coloring)
892: {
897: if (!dm->ops->getcoloring) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"No coloring for this type of DM yet");
898: (*dm->ops->getcoloring)(dm,ctype,coloring);
899: return(0);
900: }
904: /*@
905: DMCreateMatrix - Gets empty Jacobian for a DMDA or DMComposite
907: Collective on DM
909: Input Parameter:
910: . dm - the DM object
912: Output Parameter:
913: . mat - the empty Jacobian
915: Level: beginner
917: Notes: This properly preallocates the number of nonzeros in the sparse matrix so you
918: do not need to do it yourself.
920: By default it also sets the nonzero structure and puts in the zero entries. To prevent setting
921: the nonzero pattern call DMDASetMatPreallocateOnly()
923: For structured grid problems, when you call MatView() on this matrix it is displayed using the global natural ordering, NOT in the ordering used
924: internally by PETSc.
926: For structured grid problems, in general it is easiest to use MatSetValuesStencil() or MatSetValuesLocal() to put values into the matrix because MatSetValues() requires
927: the indices for the global numbering for DMDAs which is complicated.
929: .seealso DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMSetMatType()
931: @*/
932: PetscErrorCode DMCreateMatrix(DM dm,Mat *mat)
933: {
938: MatInitializePackage();
941: (*dm->ops->creatematrix)(dm,mat);
942: return(0);
943: }
947: /*@
948: DMSetMatrixPreallocateOnly - When DMCreateMatrix() is called the matrix will be properly
949: preallocated but the nonzero structure and zero values will not be set.
951: Logically Collective on DMDA
953: Input Parameter:
954: + dm - the DM
955: - only - PETSC_TRUE if only want preallocation
957: Level: developer
958: .seealso DMCreateMatrix()
959: @*/
960: PetscErrorCode DMSetMatrixPreallocateOnly(DM dm, PetscBool only)
961: {
964: dm->prealloc_only = only;
965: return(0);
966: }
970: /*@C
971: DMGetWorkArray - Gets a work array guaranteed to be at least the input size, restore with DMRestoreWorkArray()
973: Not Collective
975: Input Parameters:
976: + dm - the DM object
977: . count - The minium size
978: - dtype - data type (PETSC_REAL, PETSC_SCALAR, PETSC_INT)
980: Output Parameter:
981: . array - the work array
983: Level: developer
985: .seealso DMDestroy(), DMCreate()
986: @*/
987: PetscErrorCode DMGetWorkArray(DM dm,PetscInt count,PetscDataType dtype,void *mem)
988: {
990: DMWorkLink link;
991: size_t size;
996: if (dm->workin) {
997: link = dm->workin;
998: dm->workin = dm->workin->next;
999: } else {
1000: PetscNewLog(dm,&link);
1001: }
1002: PetscDataTypeGetSize(dtype,&size);
1003: if (size*count > link->bytes) {
1004: PetscFree(link->mem);
1005: PetscMalloc(size*count,&link->mem);
1006: link->bytes = size*count;
1007: }
1008: link->next = dm->workout;
1009: dm->workout = link;
1010: *(void**)mem = link->mem;
1011: return(0);
1012: }
1016: /*@C
1017: DMRestoreWorkArray - Restores a work array guaranteed to be at least the input size, restore with DMRestoreWorkArray()
1019: Not Collective
1021: Input Parameters:
1022: + dm - the DM object
1023: . count - The minium size
1024: - dtype - data type (PETSC_REAL, PETSC_SCALAR, PETSC_INT)
1026: Output Parameter:
1027: . array - the work array
1029: Level: developer
1031: .seealso DMDestroy(), DMCreate()
1032: @*/
1033: PetscErrorCode DMRestoreWorkArray(DM dm,PetscInt count,PetscDataType dtype,void *mem)
1034: {
1035: DMWorkLink *p,link;
1040: for (p=&dm->workout; (link=*p); p=&link->next) {
1041: if (link->mem == *(void**)mem) {
1042: *p = link->next;
1043: link->next = dm->workin;
1044: dm->workin = link;
1045: *(void**)mem = NULL;
1046: return(0);
1047: }
1048: }
1049: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Array was not checked out");
1050: return(0);
1051: }
1055: PetscErrorCode DMSetNullSpaceConstructor(DM dm, PetscInt field, PetscErrorCode (*nullsp)(DM dm, PetscInt field, MatNullSpace *nullSpace))
1056: {
1059: if (field >= 10) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle %d >= 10 fields", field);
1060: dm->nullspaceConstructors[field] = nullsp;
1061: return(0);
1062: }
1066: /*@C
1067: DMCreateFieldIS - Creates a set of IS objects with the global indices of dofs for each field
1069: Not collective
1071: Input Parameter:
1072: . dm - the DM object
1074: Output Parameters:
1075: + numFields - The number of fields (or NULL if not requested)
1076: . fieldNames - The name for each field (or NULL if not requested)
1077: - fields - The global indices for each field (or NULL if not requested)
1079: Level: intermediate
1081: Notes:
1082: The user is responsible for freeing all requested arrays. In particular, every entry of names should be freed with
1083: PetscFree(), every entry of fields should be destroyed with ISDestroy(), and both arrays should be freed with
1084: PetscFree().
1086: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix()
1087: @*/
1088: PetscErrorCode DMCreateFieldIS(DM dm, PetscInt *numFields, char ***fieldNames, IS **fields)
1089: {
1090: PetscSection section, sectionGlobal;
1095: if (numFields) {
1097: *numFields = 0;
1098: }
1099: if (fieldNames) {
1101: *fieldNames = NULL;
1102: }
1103: if (fields) {
1105: *fields = NULL;
1106: }
1107: DMGetDefaultSection(dm, §ion);
1108: if (section) {
1109: PetscInt *fieldSizes, **fieldIndices;
1110: PetscInt nF, f, pStart, pEnd, p;
1112: DMGetDefaultGlobalSection(dm, §ionGlobal);
1113: PetscSectionGetNumFields(section, &nF);
1114: PetscMalloc2(nF,&fieldSizes,nF,&fieldIndices);
1115: PetscSectionGetChart(sectionGlobal, &pStart, &pEnd);
1116: for (f = 0; f < nF; ++f) {
1117: fieldSizes[f] = 0;
1118: }
1119: for (p = pStart; p < pEnd; ++p) {
1120: PetscInt gdof;
1122: PetscSectionGetDof(sectionGlobal, p, &gdof);
1123: if (gdof > 0) {
1124: for (f = 0; f < nF; ++f) {
1125: PetscInt fdof, fcdof;
1127: PetscSectionGetFieldDof(section, p, f, &fdof);
1128: PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);
1129: fieldSizes[f] += fdof-fcdof;
1130: }
1131: }
1132: }
1133: for (f = 0; f < nF; ++f) {
1134: PetscMalloc1(fieldSizes[f], &fieldIndices[f]);
1135: fieldSizes[f] = 0;
1136: }
1137: for (p = pStart; p < pEnd; ++p) {
1138: PetscInt gdof, goff;
1140: PetscSectionGetDof(sectionGlobal, p, &gdof);
1141: if (gdof > 0) {
1142: PetscSectionGetOffset(sectionGlobal, p, &goff);
1143: for (f = 0; f < nF; ++f) {
1144: PetscInt fdof, fcdof, fc;
1146: PetscSectionGetFieldDof(section, p, f, &fdof);
1147: PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);
1148: for (fc = 0; fc < fdof-fcdof; ++fc, ++fieldSizes[f]) {
1149: fieldIndices[f][fieldSizes[f]] = goff++;
1150: }
1151: }
1152: }
1153: }
1154: if (numFields) *numFields = nF;
1155: if (fieldNames) {
1156: PetscMalloc1(nF, fieldNames);
1157: for (f = 0; f < nF; ++f) {
1158: const char *fieldName;
1160: PetscSectionGetFieldName(section, f, &fieldName);
1161: PetscStrallocpy(fieldName, (char**) &(*fieldNames)[f]);
1162: }
1163: }
1164: if (fields) {
1165: PetscMalloc1(nF, fields);
1166: for (f = 0; f < nF; ++f) {
1167: ISCreateGeneral(PetscObjectComm((PetscObject)dm), fieldSizes[f], fieldIndices[f], PETSC_OWN_POINTER, &(*fields)[f]);
1168: }
1169: }
1170: PetscFree2(fieldSizes,fieldIndices);
1171: } else if (dm->ops->createfieldis) {
1172: (*dm->ops->createfieldis)(dm, numFields, fieldNames, fields);
1173: }
1174: return(0);
1175: }
1180: /*@C
1181: DMCreateFieldDecomposition - Returns a list of IS objects defining a decomposition of a problem into subproblems
1182: corresponding to different fields: each IS contains the global indices of the dofs of the
1183: corresponding field. The optional list of DMs define the DM for each subproblem.
1184: Generalizes DMCreateFieldIS().
1186: Not collective
1188: Input Parameter:
1189: . dm - the DM object
1191: Output Parameters:
1192: + len - The number of subproblems in the field decomposition (or NULL if not requested)
1193: . namelist - The name for each field (or NULL if not requested)
1194: . islist - The global indices for each field (or NULL if not requested)
1195: - dmlist - The DMs for each field subproblem (or NULL, if not requested; if NULL is returned, no DMs are defined)
1197: Level: intermediate
1199: Notes:
1200: The user is responsible for freeing all requested arrays. In particular, every entry of names should be freed with
1201: PetscFree(), every entry of is should be destroyed with ISDestroy(), every entry of dm should be destroyed with DMDestroy(),
1202: and all of the arrays should be freed with PetscFree().
1204: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateFieldIS()
1205: @*/
1206: PetscErrorCode DMCreateFieldDecomposition(DM dm, PetscInt *len, char ***namelist, IS **islist, DM **dmlist)
1207: {
1212: if (len) {
1214: *len = 0;
1215: }
1216: if (namelist) {
1218: *namelist = 0;
1219: }
1220: if (islist) {
1222: *islist = 0;
1223: }
1224: if (dmlist) {
1226: *dmlist = 0;
1227: }
1228: /*
1229: Is it a good idea to apply the following check across all impls?
1230: Perhaps some impls can have a well-defined decomposition before DMSetUp?
1231: This, however, follows the general principle that accessors are not well-behaved until the object is set up.
1232: */
1233: if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE, "Decomposition defined only after DMSetUp");
1234: if (!dm->ops->createfielddecomposition) {
1235: PetscSection section;
1236: PetscInt numFields, f;
1238: DMGetDefaultSection(dm, §ion);
1239: if (section) {PetscSectionGetNumFields(section, &numFields);}
1240: if (section && numFields && dm->ops->createsubdm) {
1241: *len = numFields;
1242: if (namelist) {PetscMalloc1(numFields,namelist);}
1243: if (islist) {PetscMalloc1(numFields,islist);}
1244: if (dmlist) {PetscMalloc1(numFields,dmlist);}
1245: for (f = 0; f < numFields; ++f) {
1246: const char *fieldName;
1248: DMCreateSubDM(dm, 1, &f, islist ? &(*islist)[f] : NULL, dmlist ? &(*dmlist)[f] : NULL);
1249: if (namelist) {
1250: PetscSectionGetFieldName(section, f, &fieldName);
1251: PetscStrallocpy(fieldName, (char**) &(*namelist)[f]);
1252: }
1253: }
1254: } else {
1255: DMCreateFieldIS(dm, len, namelist, islist);
1256: /* By default there are no DMs associated with subproblems. */
1257: if (dmlist) *dmlist = NULL;
1258: }
1259: } else {
1260: (*dm->ops->createfielddecomposition)(dm,len,namelist,islist,dmlist);
1261: }
1262: return(0);
1263: }
1267: /*@C
1268: DMCreateSubDM - Returns an IS and DM encapsulating a subproblem defined by the fields passed in.
1269: The fields are defined by DMCreateFieldIS().
1271: Not collective
1273: Input Parameters:
1274: + dm - the DM object
1275: . numFields - number of fields in this subproblem
1276: - len - The number of subproblems in the decomposition (or NULL if not requested)
1278: Output Parameters:
1279: . is - The global indices for the subproblem
1280: - dm - The DM for the subproblem
1282: Level: intermediate
1284: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateFieldIS()
1285: @*/
1286: PetscErrorCode DMCreateSubDM(DM dm, PetscInt numFields, PetscInt fields[], IS *is, DM *subdm)
1287: {
1295: if (dm->ops->createsubdm) {
1296: (*dm->ops->createsubdm)(dm, numFields, fields, is, subdm);
1297: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "This type has no DMCreateSubDM implementation defined");
1298: return(0);
1299: }
1304: /*@C
1305: DMCreateDomainDecomposition - Returns lists of IS objects defining a decomposition of a problem into subproblems
1306: corresponding to restrictions to pairs nested subdomains: each IS contains the global
1307: indices of the dofs of the corresponding subdomains. The inner subdomains conceptually
1308: define a nonoverlapping covering, while outer subdomains can overlap.
1309: The optional list of DMs define the DM for each subproblem.
1311: Not collective
1313: Input Parameter:
1314: . dm - the DM object
1316: Output Parameters:
1317: + len - The number of subproblems in the domain decomposition (or NULL if not requested)
1318: . namelist - The name for each subdomain (or NULL if not requested)
1319: . innerislist - The global indices for each inner subdomain (or NULL, if not requested)
1320: . outerislist - The global indices for each outer subdomain (or NULL, if not requested)
1321: - dmlist - The DMs for each subdomain subproblem (or NULL, if not requested; if NULL is returned, no DMs are defined)
1323: Level: intermediate
1325: Notes:
1326: The user is responsible for freeing all requested arrays. In particular, every entry of names should be freed with
1327: PetscFree(), every entry of is should be destroyed with ISDestroy(), every entry of dm should be destroyed with DMDestroy(),
1328: and all of the arrays should be freed with PetscFree().
1330: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateDomainDecompositionDM(), DMCreateFieldDecomposition()
1331: @*/
1332: PetscErrorCode DMCreateDomainDecomposition(DM dm, PetscInt *len, char ***namelist, IS **innerislist, IS **outerislist, DM **dmlist)
1333: {
1334: PetscErrorCode ierr;
1335: DMSubDomainHookLink link;
1336: PetscInt i,l;
1345: /*
1346: Is it a good idea to apply the following check across all impls?
1347: Perhaps some impls can have a well-defined decomposition before DMSetUp?
1348: This, however, follows the general principle that accessors are not well-behaved until the object is set up.
1349: */
1350: if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE, "Decomposition defined only after DMSetUp");
1351: if (dm->ops->createdomaindecomposition) {
1352: (*dm->ops->createdomaindecomposition)(dm,&l,namelist,innerislist,outerislist,dmlist);
1353: /* copy subdomain hooks and context over to the subdomain DMs */
1354: if (dmlist) {
1355: for (i = 0; i < l; i++) {
1356: for (link=dm->subdomainhook; link; link=link->next) {
1357: if (link->ddhook) {(*link->ddhook)(dm,(*dmlist)[i],link->ctx);}
1358: }
1359: (*dmlist)[i]->ctx = dm->ctx;
1360: }
1361: }
1362: if (len) *len = l;
1363: }
1364: return(0);
1365: }
1370: /*@C
1371: DMCreateDomainDecompositionScatters - Returns scatters to the subdomain vectors from the global vector
1373: Not collective
1375: Input Parameters:
1376: + dm - the DM object
1377: . n - the number of subdomain scatters
1378: - subdms - the local subdomains
1380: Output Parameters:
1381: + n - the number of scatters returned
1382: . iscat - scatter from global vector to nonoverlapping global vector entries on subdomain
1383: . oscat - scatter from global vector to overlapping global vector entries on subdomain
1384: - gscat - scatter from global vector to local vector on subdomain (fills in ghosts)
1386: Notes: This is an alternative to the iis and ois arguments in DMCreateDomainDecomposition that allow for the solution
1387: of general nonlinear problems with overlapping subdomain methods. While merely having index sets that enable subsets
1388: of the residual equations to be created is fine for linear problems, nonlinear problems require local assembly of
1389: solution and residual data.
1391: Level: developer
1393: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateFieldIS()
1394: @*/
1395: PetscErrorCode DMCreateDomainDecompositionScatters(DM dm,PetscInt n,DM *subdms,VecScatter **iscat,VecScatter **oscat,VecScatter **gscat)
1396: {
1402: if (dm->ops->createddscatters) {
1403: (*dm->ops->createddscatters)(dm,n,subdms,iscat,oscat,gscat);
1404: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "This type has no DMCreateDomainDecompositionLocalScatter implementation defined");
1405: return(0);
1406: }
1410: /*@
1411: DMRefine - Refines a DM object
1413: Collective on DM
1415: Input Parameter:
1416: + dm - the DM object
1417: - comm - the communicator to contain the new DM object (or MPI_COMM_NULL)
1419: Output Parameter:
1420: . dmf - the refined DM, or NULL
1422: Note: If no refinement was done, the return value is NULL
1424: Level: developer
1426: .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation()
1427: @*/
1428: PetscErrorCode DMRefine(DM dm,MPI_Comm comm,DM *dmf)
1429: {
1430: PetscErrorCode ierr;
1431: DMRefineHookLink link;
1435: (*dm->ops->refine)(dm,comm,dmf);
1436: if (*dmf) {
1437: (*dmf)->ops->creatematrix = dm->ops->creatematrix;
1439: PetscObjectCopyFortranFunctionPointers((PetscObject)dm,(PetscObject)*dmf);
1441: (*dmf)->ctx = dm->ctx;
1442: (*dmf)->leveldown = dm->leveldown;
1443: (*dmf)->levelup = dm->levelup + 1;
1445: DMSetMatType(*dmf,dm->mattype);
1446: for (link=dm->refinehook; link; link=link->next) {
1447: if (link->refinehook) {
1448: (*link->refinehook)(dm,*dmf,link->ctx);
1449: }
1450: }
1451: }
1452: return(0);
1453: }
1457: /*@C
1458: DMRefineHookAdd - adds a callback to be run when interpolating a nonlinear problem to a finer grid
1460: Logically Collective
1462: Input Arguments:
1463: + coarse - nonlinear solver context on which to run a hook when restricting to a coarser level
1464: . refinehook - function to run when setting up a coarser level
1465: . interphook - function to run to update data on finer levels (once per SNESSolve())
1466: - ctx - [optional] user-defined context for provide data for the hooks (may be NULL)
1468: Calling sequence of refinehook:
1469: $ refinehook(DM coarse,DM fine,void *ctx);
1471: + coarse - coarse level DM
1472: . fine - fine level DM to interpolate problem to
1473: - ctx - optional user-defined function context
1475: Calling sequence for interphook:
1476: $ interphook(DM coarse,Mat interp,DM fine,void *ctx)
1478: + coarse - coarse level DM
1479: . interp - matrix interpolating a coarse-level solution to the finer grid
1480: . fine - fine level DM to update
1481: - ctx - optional user-defined function context
1483: Level: advanced
1485: Notes:
1486: This function is only needed if auxiliary data needs to be passed to fine grids while grid sequencing
1488: If this function is called multiple times, the hooks will be run in the order they are added.
1490: This function is currently not available from Fortran.
1492: .seealso: DMCoarsenHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
1493: @*/
1494: PetscErrorCode DMRefineHookAdd(DM coarse,PetscErrorCode (*refinehook)(DM,DM,void*),PetscErrorCode (*interphook)(DM,Mat,DM,void*),void *ctx)
1495: {
1496: PetscErrorCode ierr;
1497: DMRefineHookLink link,*p;
1501: for (p=&coarse->refinehook; *p; p=&(*p)->next) {} /* Scan to the end of the current list of hooks */
1502: PetscMalloc(sizeof(struct _DMRefineHookLink),&link);
1503: link->refinehook = refinehook;
1504: link->interphook = interphook;
1505: link->ctx = ctx;
1506: link->next = NULL;
1507: *p = link;
1508: return(0);
1509: }
1513: /*@
1514: DMInterpolate - interpolates user-defined problem data to a finer DM by running hooks registered by DMRefineHookAdd()
1516: Collective if any hooks are
1518: Input Arguments:
1519: + coarse - coarser DM to use as a base
1520: . restrct - interpolation matrix, apply using MatInterpolate()
1521: - fine - finer DM to update
1523: Level: developer
1525: .seealso: DMRefineHookAdd(), MatInterpolate()
1526: @*/
1527: PetscErrorCode DMInterpolate(DM coarse,Mat interp,DM fine)
1528: {
1529: PetscErrorCode ierr;
1530: DMRefineHookLink link;
1533: for (link=fine->refinehook; link; link=link->next) {
1534: if (link->interphook) {
1535: (*link->interphook)(coarse,interp,fine,link->ctx);
1536: }
1537: }
1538: return(0);
1539: }
1543: /*@
1544: DMGetRefineLevel - Get's the number of refinements that have generated this DM.
1546: Not Collective
1548: Input Parameter:
1549: . dm - the DM object
1551: Output Parameter:
1552: . level - number of refinements
1554: Level: developer
1556: .seealso DMCoarsen(), DMGetCoarsenLevel(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation()
1558: @*/
1559: PetscErrorCode DMGetRefineLevel(DM dm,PetscInt *level)
1560: {
1563: *level = dm->levelup;
1564: return(0);
1565: }
1569: /*@C
1570: DMGlobalToLocalHookAdd - adds a callback to be run when global to local is called
1572: Logically Collective
1574: Input Arguments:
1575: + dm - the DM
1576: . beginhook - function to run at the beginning of DMGlobalToLocalBegin()
1577: . endhook - function to run after DMGlobalToLocalEnd() has completed
1578: - ctx - [optional] user-defined context for provide data for the hooks (may be NULL)
1580: Calling sequence for beginhook:
1581: $ beginhook(DM fine,VecScatter out,VecScatter in,DM coarse,void *ctx)
1583: + dm - global DM
1584: . g - global vector
1585: . mode - mode
1586: . l - local vector
1587: - ctx - optional user-defined function context
1590: Calling sequence for endhook:
1591: $ endhook(DM fine,VecScatter out,VecScatter in,DM coarse,void *ctx)
1593: + global - global DM
1594: - ctx - optional user-defined function context
1596: Level: advanced
1598: .seealso: DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
1599: @*/
1600: PetscErrorCode DMGlobalToLocalHookAdd(DM dm,PetscErrorCode (*beginhook)(DM,Vec,InsertMode,Vec,void*),PetscErrorCode (*endhook)(DM,Vec,InsertMode,Vec,void*),void *ctx)
1601: {
1602: PetscErrorCode ierr;
1603: DMGlobalToLocalHookLink link,*p;
1607: for (p=&dm->gtolhook; *p; p=&(*p)->next) {} /* Scan to the end of the current list of hooks */
1608: PetscMalloc(sizeof(struct _DMGlobalToLocalHookLink),&link);
1609: link->beginhook = beginhook;
1610: link->endhook = endhook;
1611: link->ctx = ctx;
1612: link->next = NULL;
1613: *p = link;
1614: return(0);
1615: }
1619: /*@
1620: DMGlobalToLocalBegin - Begins updating local vectors from global vector
1622: Neighbor-wise Collective on DM
1624: Input Parameters:
1625: + dm - the DM object
1626: . g - the global vector
1627: . mode - INSERT_VALUES or ADD_VALUES
1628: - l - the local vector
1631: Level: beginner
1633: .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMGlobalToLocalEnd(), DMLocalToGlobalBegin()
1635: @*/
1636: PetscErrorCode DMGlobalToLocalBegin(DM dm,Vec g,InsertMode mode,Vec l)
1637: {
1638: PetscSF sf;
1639: PetscErrorCode ierr;
1640: DMGlobalToLocalHookLink link;
1644: for (link=dm->gtolhook; link; link=link->next) {
1645: if (link->beginhook) {
1646: (*link->beginhook)(dm,g,mode,l,link->ctx);
1647: }
1648: }
1649: DMGetDefaultSF(dm, &sf);
1650: if (sf) {
1651: PetscScalar *lArray, *gArray;
1653: if (mode == ADD_VALUES) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode);
1654: VecGetArray(l, &lArray);
1655: VecGetArray(g, &gArray);
1656: PetscSFBcastBegin(sf, MPIU_SCALAR, gArray, lArray);
1657: VecRestoreArray(l, &lArray);
1658: VecRestoreArray(g, &gArray);
1659: } else {
1660: (*dm->ops->globaltolocalbegin)(dm,g,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),l);
1661: }
1662: return(0);
1663: }
1667: /*@
1668: DMGlobalToLocalEnd - Ends updating local vectors from global vector
1670: Neighbor-wise Collective on DM
1672: Input Parameters:
1673: + dm - the DM object
1674: . g - the global vector
1675: . mode - INSERT_VALUES or ADD_VALUES
1676: - l - the local vector
1679: Level: beginner
1681: .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMGlobalToLocalEnd(), DMLocalToGlobalBegin()
1683: @*/
1684: PetscErrorCode DMGlobalToLocalEnd(DM dm,Vec g,InsertMode mode,Vec l)
1685: {
1686: PetscSF sf;
1687: PetscErrorCode ierr;
1688: PetscScalar *lArray, *gArray;
1689: DMGlobalToLocalHookLink link;
1693: DMGetDefaultSF(dm, &sf);
1694: if (sf) {
1695: if (mode == ADD_VALUES) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode);
1697: VecGetArray(l, &lArray);
1698: VecGetArray(g, &gArray);
1699: PetscSFBcastEnd(sf, MPIU_SCALAR, gArray, lArray);
1700: VecRestoreArray(l, &lArray);
1701: VecRestoreArray(g, &gArray);
1702: } else {
1703: (*dm->ops->globaltolocalend)(dm,g,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),l);
1704: }
1705: for (link=dm->gtolhook; link; link=link->next) {
1706: if (link->endhook) {(*link->endhook)(dm,g,mode,l,link->ctx);}
1707: }
1708: return(0);
1709: }
1713: /*@
1714: DMLocalToGlobalBegin - updates global vectors from local vectors
1716: Neighbor-wise Collective on DM
1718: Input Parameters:
1719: + dm - the DM object
1720: . l - the local vector
1721: . 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.
1722: - g - the global vector
1724: Notes: In the ADD_VALUES case you normally would zero the receiving vector before beginning this operation. If you would like to simply add the non-ghosted values in the local
1725: array into the global array you need to either (1) zero the ghosted locations and use ADD_VALUES or (2) use INSERT_VALUES into a work global array and then add the work
1726: global array to the final global array with VecAXPY().
1728: Level: beginner
1730: .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMGlobalToLocalEnd(), DMGlobalToLocalBegin()
1732: @*/
1733: PetscErrorCode DMLocalToGlobalBegin(DM dm,Vec l,InsertMode mode,Vec g)
1734: {
1735: PetscSF sf;
1740: DMGetDefaultSF(dm, &sf);
1741: if (sf) {
1742: MPI_Op op;
1743: PetscScalar *lArray, *gArray;
1745: switch (mode) {
1746: case INSERT_VALUES:
1747: case INSERT_ALL_VALUES:
1748: op = MPIU_REPLACE; break;
1749: case ADD_VALUES:
1750: case ADD_ALL_VALUES:
1751: op = MPI_SUM; break;
1752: default:
1753: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode);
1754: }
1755: VecGetArray(l, &lArray);
1756: VecGetArray(g, &gArray);
1757: PetscSFReduceBegin(sf, MPIU_SCALAR, lArray, gArray, op);
1758: VecRestoreArray(l, &lArray);
1759: VecRestoreArray(g, &gArray);
1760: } else {
1761: (*dm->ops->localtoglobalbegin)(dm,l,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),g);
1762: }
1763: return(0);
1764: }
1768: /*@
1769: DMLocalToGlobalEnd - updates global vectors from local vectors
1771: Neighbor-wise Collective on DM
1773: Input Parameters:
1774: + dm - the DM object
1775: . l - the local vector
1776: . mode - INSERT_VALUES or ADD_VALUES
1777: - g - the global vector
1780: Level: beginner
1782: .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMGlobalToLocalEnd(), DMGlobalToLocalEnd()
1784: @*/
1785: PetscErrorCode DMLocalToGlobalEnd(DM dm,Vec l,InsertMode mode,Vec g)
1786: {
1787: PetscSF sf;
1792: DMGetDefaultSF(dm, &sf);
1793: if (sf) {
1794: MPI_Op op;
1795: PetscScalar *lArray, *gArray;
1797: switch (mode) {
1798: case INSERT_VALUES:
1799: case INSERT_ALL_VALUES:
1800: op = MPIU_REPLACE; break;
1801: case ADD_VALUES:
1802: case ADD_ALL_VALUES:
1803: op = MPI_SUM; break;
1804: default:
1805: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode);
1806: }
1807: VecGetArray(l, &lArray);
1808: VecGetArray(g, &gArray);
1809: PetscSFReduceEnd(sf, MPIU_SCALAR, lArray, gArray, op);
1810: VecRestoreArray(l, &lArray);
1811: VecRestoreArray(g, &gArray);
1812: } else {
1813: (*dm->ops->localtoglobalend)(dm,l,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),g);
1814: }
1815: return(0);
1816: }
1820: /*@
1821: DMLocalToLocalBegin - Maps from a local vector (including ghost points
1822: that contain irrelevant values) to another local vector where the ghost
1823: points in the second are set correctly. Must be followed by DMLocalToLocalEnd().
1825: Neighbor-wise Collective on DM and Vec
1827: Input Parameters:
1828: + dm - the DM object
1829: . g - the original local vector
1830: - mode - one of INSERT_VALUES or ADD_VALUES
1832: Output Parameter:
1833: . l - the local vector with correct ghost values
1835: Level: intermediate
1837: Notes:
1838: The local vectors used here need not be the same as those
1839: obtained from DMCreateLocalVector(), BUT they
1840: must have the same parallel data layout; they could, for example, be
1841: obtained with VecDuplicate() from the DM originating vectors.
1843: .keywords: DM, local-to-local, begin
1844: .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateLocalVector(), DMCreateGlobalVector(), DMCreateInterpolation(), DMLocalToLocalEnd(), DMGlobalToLocalEnd(), DMLocalToGlobalBegin()
1846: @*/
1847: PetscErrorCode DMLocalToLocalBegin(DM dm,Vec g,InsertMode mode,Vec l)
1848: {
1849: PetscErrorCode ierr;
1853: (*dm->ops->localtolocalbegin)(dm,g,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),l);
1854: return(0);
1855: }
1859: /*@
1860: DMLocalToLocalEnd - Maps from a local vector (including ghost points
1861: that contain irrelevant values) to another local vector where the ghost
1862: points in the second are set correctly. Must be preceded by DMLocalToLocalBegin().
1864: Neighbor-wise Collective on DM and Vec
1866: Input Parameters:
1867: + da - the DM object
1868: . g - the original local vector
1869: - mode - one of INSERT_VALUES or ADD_VALUES
1871: Output Parameter:
1872: . l - the local vector with correct ghost values
1874: Level: intermediate
1876: Notes:
1877: The local vectors used here need not be the same as those
1878: obtained from DMCreateLocalVector(), BUT they
1879: must have the same parallel data layout; they could, for example, be
1880: obtained with VecDuplicate() from the DM originating vectors.
1882: .keywords: DM, local-to-local, end
1883: .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateLocalVector(), DMCreateGlobalVector(), DMCreateInterpolation(), DMLocalToLocalBegin(), DMGlobalToLocalEnd(), DMLocalToGlobalBegin()
1885: @*/
1886: PetscErrorCode DMLocalToLocalEnd(DM dm,Vec g,InsertMode mode,Vec l)
1887: {
1888: PetscErrorCode ierr;
1892: (*dm->ops->localtolocalend)(dm,g,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),l);
1893: return(0);
1894: }
1899: /*@
1900: DMCoarsen - Coarsens a DM object
1902: Collective on DM
1904: Input Parameter:
1905: + dm - the DM object
1906: - comm - the communicator to contain the new DM object (or MPI_COMM_NULL)
1908: Output Parameter:
1909: . dmc - the coarsened DM
1911: Level: developer
1913: .seealso DMRefine(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation()
1915: @*/
1916: PetscErrorCode DMCoarsen(DM dm, MPI_Comm comm, DM *dmc)
1917: {
1918: PetscErrorCode ierr;
1919: DMCoarsenHookLink link;
1923: (*dm->ops->coarsen)(dm, comm, dmc);
1924: (*dmc)->ops->creatematrix = dm->ops->creatematrix;
1925: PetscObjectCopyFortranFunctionPointers((PetscObject)dm,(PetscObject)*dmc);
1926: (*dmc)->ctx = dm->ctx;
1927: (*dmc)->levelup = dm->levelup;
1928: (*dmc)->leveldown = dm->leveldown + 1;
1929: DMSetMatType(*dmc,dm->mattype);
1930: for (link=dm->coarsenhook; link; link=link->next) {
1931: if (link->coarsenhook) {(*link->coarsenhook)(dm,*dmc,link->ctx);}
1932: }
1933: return(0);
1934: }
1938: /*@C
1939: DMCoarsenHookAdd - adds a callback to be run when restricting a nonlinear problem to the coarse grid
1941: Logically Collective
1943: Input Arguments:
1944: + fine - nonlinear solver context on which to run a hook when restricting to a coarser level
1945: . coarsenhook - function to run when setting up a coarser level
1946: . restricthook - function to run to update data on coarser levels (once per SNESSolve())
1947: - ctx - [optional] user-defined context for provide data for the hooks (may be NULL)
1949: Calling sequence of coarsenhook:
1950: $ coarsenhook(DM fine,DM coarse,void *ctx);
1952: + fine - fine level DM
1953: . coarse - coarse level DM to restrict problem to
1954: - ctx - optional user-defined function context
1956: Calling sequence for restricthook:
1957: $ restricthook(DM fine,Mat mrestrict,Vec rscale,Mat inject,DM coarse,void *ctx)
1959: + fine - fine level DM
1960: . mrestrict - matrix restricting a fine-level solution to the coarse grid
1961: . rscale - scaling vector for restriction
1962: . inject - matrix restricting by injection
1963: . coarse - coarse level DM to update
1964: - ctx - optional user-defined function context
1966: Level: advanced
1968: Notes:
1969: This function is only needed if auxiliary data needs to be set up on coarse grids.
1971: If this function is called multiple times, the hooks will be run in the order they are added.
1973: In order to compose with nonlinear preconditioning without duplicating storage, the hook should be implemented to
1974: extract the finest level information from its context (instead of from the SNES).
1976: This function is currently not available from Fortran.
1978: .seealso: DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
1979: @*/
1980: PetscErrorCode DMCoarsenHookAdd(DM fine,PetscErrorCode (*coarsenhook)(DM,DM,void*),PetscErrorCode (*restricthook)(DM,Mat,Vec,Mat,DM,void*),void *ctx)
1981: {
1982: PetscErrorCode ierr;
1983: DMCoarsenHookLink link,*p;
1987: for (p=&fine->coarsenhook; *p; p=&(*p)->next) {} /* Scan to the end of the current list of hooks */
1988: PetscMalloc(sizeof(struct _DMCoarsenHookLink),&link);
1989: link->coarsenhook = coarsenhook;
1990: link->restricthook = restricthook;
1991: link->ctx = ctx;
1992: link->next = NULL;
1993: *p = link;
1994: return(0);
1995: }
1999: /*@
2000: DMRestrict - restricts user-defined problem data to a coarser DM by running hooks registered by DMCoarsenHookAdd()
2002: Collective if any hooks are
2004: Input Arguments:
2005: + fine - finer DM to use as a base
2006: . restrct - restriction matrix, apply using MatRestrict()
2007: . inject - injection matrix, also use MatRestrict()
2008: - coarse - coarer DM to update
2010: Level: developer
2012: .seealso: DMCoarsenHookAdd(), MatRestrict()
2013: @*/
2014: PetscErrorCode DMRestrict(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse)
2015: {
2016: PetscErrorCode ierr;
2017: DMCoarsenHookLink link;
2020: for (link=fine->coarsenhook; link; link=link->next) {
2021: if (link->restricthook) {
2022: (*link->restricthook)(fine,restrct,rscale,inject,coarse,link->ctx);
2023: }
2024: }
2025: return(0);
2026: }
2030: /*@C
2031: DMSubDomainHookAdd - adds a callback to be run when restricting a problem to the coarse grid
2033: Logically Collective
2035: Input Arguments:
2036: + global - global DM
2037: . ddhook - function to run to pass data to the decomposition DM upon its creation
2038: . restricthook - function to run to update data on block solve (at the beginning of the block solve)
2039: - ctx - [optional] user-defined context for provide data for the hooks (may be NULL)
2042: Calling sequence for ddhook:
2043: $ ddhook(DM global,DM block,void *ctx)
2045: + global - global DM
2046: . block - block DM
2047: - ctx - optional user-defined function context
2049: Calling sequence for restricthook:
2050: $ restricthook(DM global,VecScatter out,VecScatter in,DM block,void *ctx)
2052: + global - global DM
2053: . out - scatter to the outer (with ghost and overlap points) block vector
2054: . in - scatter to block vector values only owned locally
2055: . block - block DM
2056: - ctx - optional user-defined function context
2058: Level: advanced
2060: Notes:
2061: This function is only needed if auxiliary data needs to be set up on subdomain DMs.
2063: If this function is called multiple times, the hooks will be run in the order they are added.
2065: In order to compose with nonlinear preconditioning without duplicating storage, the hook should be implemented to
2066: extract the global information from its context (instead of from the SNES).
2068: This function is currently not available from Fortran.
2070: .seealso: DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2071: @*/
2072: PetscErrorCode DMSubDomainHookAdd(DM global,PetscErrorCode (*ddhook)(DM,DM,void*),PetscErrorCode (*restricthook)(DM,VecScatter,VecScatter,DM,void*),void *ctx)
2073: {
2074: PetscErrorCode ierr;
2075: DMSubDomainHookLink link,*p;
2079: for (p=&global->subdomainhook; *p; p=&(*p)->next) {} /* Scan to the end of the current list of hooks */
2080: PetscMalloc(sizeof(struct _DMSubDomainHookLink),&link);
2081: link->restricthook = restricthook;
2082: link->ddhook = ddhook;
2083: link->ctx = ctx;
2084: link->next = NULL;
2085: *p = link;
2086: return(0);
2087: }
2091: /*@
2092: DMSubDomainRestrict - restricts user-defined problem data to a block DM by running hooks registered by DMSubDomainHookAdd()
2094: Collective if any hooks are
2096: Input Arguments:
2097: + fine - finer DM to use as a base
2098: . oscatter - scatter from domain global vector filling subdomain global vector with overlap
2099: . gscatter - scatter from domain global vector filling subdomain local vector with ghosts
2100: - coarse - coarer DM to update
2102: Level: developer
2104: .seealso: DMCoarsenHookAdd(), MatRestrict()
2105: @*/
2106: PetscErrorCode DMSubDomainRestrict(DM global,VecScatter oscatter,VecScatter gscatter,DM subdm)
2107: {
2108: PetscErrorCode ierr;
2109: DMSubDomainHookLink link;
2112: for (link=global->subdomainhook; link; link=link->next) {
2113: if (link->restricthook) {
2114: (*link->restricthook)(global,oscatter,gscatter,subdm,link->ctx);
2115: }
2116: }
2117: return(0);
2118: }
2122: /*@
2123: DMGetCoarsenLevel - Get's the number of coarsenings that have generated this DM.
2125: Not Collective
2127: Input Parameter:
2128: . dm - the DM object
2130: Output Parameter:
2131: . level - number of coarsenings
2133: Level: developer
2135: .seealso DMCoarsen(), DMGetRefineLevel(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation()
2137: @*/
2138: PetscErrorCode DMGetCoarsenLevel(DM dm,PetscInt *level)
2139: {
2142: *level = dm->leveldown;
2143: return(0);
2144: }
2150: /*@C
2151: DMRefineHierarchy - Refines a DM object, all levels at once
2153: Collective on DM
2155: Input Parameter:
2156: + dm - the DM object
2157: - nlevels - the number of levels of refinement
2159: Output Parameter:
2160: . dmf - the refined DM hierarchy
2162: Level: developer
2164: .seealso DMCoarsenHierarchy(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation()
2166: @*/
2167: PetscErrorCode DMRefineHierarchy(DM dm,PetscInt nlevels,DM dmf[])
2168: {
2173: if (nlevels < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"nlevels cannot be negative");
2174: if (nlevels == 0) return(0);
2175: if (dm->ops->refinehierarchy) {
2176: (*dm->ops->refinehierarchy)(dm,nlevels,dmf);
2177: } else if (dm->ops->refine) {
2178: PetscInt i;
2180: DMRefine(dm,PetscObjectComm((PetscObject)dm),&dmf[0]);
2181: for (i=1; i<nlevels; i++) {
2182: DMRefine(dmf[i-1],PetscObjectComm((PetscObject)dm),&dmf[i]);
2183: }
2184: } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"No RefineHierarchy for this DM yet");
2185: return(0);
2186: }
2190: /*@C
2191: DMCoarsenHierarchy - Coarsens a DM object, all levels at once
2193: Collective on DM
2195: Input Parameter:
2196: + dm - the DM object
2197: - nlevels - the number of levels of coarsening
2199: Output Parameter:
2200: . dmc - the coarsened DM hierarchy
2202: Level: developer
2204: .seealso DMRefineHierarchy(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation()
2206: @*/
2207: PetscErrorCode DMCoarsenHierarchy(DM dm, PetscInt nlevels, DM dmc[])
2208: {
2213: if (nlevels < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"nlevels cannot be negative");
2214: if (nlevels == 0) return(0);
2216: if (dm->ops->coarsenhierarchy) {
2217: (*dm->ops->coarsenhierarchy)(dm, nlevels, dmc);
2218: } else if (dm->ops->coarsen) {
2219: PetscInt i;
2221: DMCoarsen(dm,PetscObjectComm((PetscObject)dm),&dmc[0]);
2222: for (i=1; i<nlevels; i++) {
2223: DMCoarsen(dmc[i-1],PetscObjectComm((PetscObject)dm),&dmc[i]);
2224: }
2225: } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"No CoarsenHierarchy for this DM yet");
2226: return(0);
2227: }
2231: /*@
2232: DMCreateAggregates - Gets the aggregates that map between
2233: grids associated with two DMs.
2235: Collective on DM
2237: Input Parameters:
2238: + dmc - the coarse grid DM
2239: - dmf - the fine grid DM
2241: Output Parameters:
2242: . rest - the restriction matrix (transpose of the projection matrix)
2244: Level: intermediate
2246: .keywords: interpolation, restriction, multigrid
2248: .seealso: DMRefine(), DMCreateInjection(), DMCreateInterpolation()
2249: @*/
2250: PetscErrorCode DMCreateAggregates(DM dmc, DM dmf, Mat *rest)
2251: {
2257: (*dmc->ops->getaggregates)(dmc, dmf, rest);
2258: return(0);
2259: }
2263: /*@C
2264: DMSetApplicationContextDestroy - Sets a user function that will be called to destroy the application context when the DM is destroyed
2266: Not Collective
2268: Input Parameters:
2269: + dm - the DM object
2270: - destroy - the destroy function
2272: Level: intermediate
2274: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext()
2276: @*/
2277: PetscErrorCode DMSetApplicationContextDestroy(DM dm,PetscErrorCode (*destroy)(void**))
2278: {
2281: dm->ctxdestroy = destroy;
2282: return(0);
2283: }
2287: /*@
2288: DMSetApplicationContext - Set a user context into a DM object
2290: Not Collective
2292: Input Parameters:
2293: + dm - the DM object
2294: - ctx - the user context
2296: Level: intermediate
2298: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext()
2300: @*/
2301: PetscErrorCode DMSetApplicationContext(DM dm,void *ctx)
2302: {
2305: dm->ctx = ctx;
2306: return(0);
2307: }
2311: /*@
2312: DMGetApplicationContext - Gets a user context from a DM object
2314: Not Collective
2316: Input Parameter:
2317: . dm - the DM object
2319: Output Parameter:
2320: . ctx - the user context
2322: Level: intermediate
2324: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext()
2326: @*/
2327: PetscErrorCode DMGetApplicationContext(DM dm,void *ctx)
2328: {
2331: *(void**)ctx = dm->ctx;
2332: return(0);
2333: }
2337: /*@C
2338: DMSetVariableBounds - sets a function to compute the lower and upper bound vectors for SNESVI.
2340: Logically Collective on DM
2342: Input Parameter:
2343: + dm - the DM object
2344: - f - the function that computes variable bounds used by SNESVI (use NULL to cancel a previous function that was set)
2346: Level: intermediate
2348: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext(),
2349: DMSetJacobian()
2351: @*/
2352: PetscErrorCode DMSetVariableBounds(DM dm,PetscErrorCode (*f)(DM,Vec,Vec))
2353: {
2355: dm->ops->computevariablebounds = f;
2356: return(0);
2357: }
2361: /*@
2362: DMHasVariableBounds - does the DM object have a variable bounds function?
2364: Not Collective
2366: Input Parameter:
2367: . dm - the DM object to destroy
2369: Output Parameter:
2370: . flg - PETSC_TRUE if the variable bounds function exists
2372: Level: developer
2374: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext()
2376: @*/
2377: PetscErrorCode DMHasVariableBounds(DM dm,PetscBool *flg)
2378: {
2380: *flg = (dm->ops->computevariablebounds) ? PETSC_TRUE : PETSC_FALSE;
2381: return(0);
2382: }
2386: /*@C
2387: DMComputeVariableBounds - compute variable bounds used by SNESVI.
2389: Logically Collective on DM
2391: Input Parameters:
2392: + dm - the DM object to destroy
2393: - x - current solution at which the bounds are computed
2395: Output parameters:
2396: + xl - lower bound
2397: - xu - upper bound
2399: Level: intermediate
2401: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext()
2403: @*/
2404: PetscErrorCode DMComputeVariableBounds(DM dm, Vec xl, Vec xu)
2405: {
2411: if (dm->ops->computevariablebounds) {
2412: (*dm->ops->computevariablebounds)(dm, xl,xu);
2413: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "This DM is incapable of computing variable bounds.");
2414: return(0);
2415: }
2419: /*@
2420: DMHasColoring - does the DM object have a method of providing a coloring?
2422: Not Collective
2424: Input Parameter:
2425: . dm - the DM object
2427: Output Parameter:
2428: . flg - PETSC_TRUE if the DM has facilities for DMCreateColoring().
2430: Level: developer
2432: .seealso DMHasFunction(), DMCreateColoring()
2434: @*/
2435: PetscErrorCode DMHasColoring(DM dm,PetscBool *flg)
2436: {
2438: *flg = (dm->ops->getcoloring) ? PETSC_TRUE : PETSC_FALSE;
2439: return(0);
2440: }
2442: #undef __FUNCT__
2444: /*@C
2445: DMSetVec - set the vector at which to compute residual, Jacobian and VI bounds, if the problem is nonlinear.
2447: Collective on DM
2449: Input Parameter:
2450: + dm - the DM object
2451: - x - location to compute residual and Jacobian, if NULL is passed to those routines; will be NULL for linear problems.
2453: Level: developer
2455: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext()
2457: @*/
2458: PetscErrorCode DMSetVec(DM dm,Vec x)
2459: {
2463: if (x) {
2464: if (!dm->x) {
2465: DMCreateGlobalVector(dm,&dm->x);
2466: }
2467: VecCopy(x,dm->x);
2468: } else if (dm->x) {
2469: VecDestroy(&dm->x);
2470: }
2471: return(0);
2472: }
2474: PetscFunctionList DMList = NULL;
2475: PetscBool DMRegisterAllCalled = PETSC_FALSE;
2479: /*@C
2480: DMSetType - Builds a DM, for a particular DM implementation.
2482: Collective on DM
2484: Input Parameters:
2485: + dm - The DM object
2486: - method - The name of the DM type
2488: Options Database Key:
2489: . -dm_type <type> - Sets the DM type; use -help for a list of available types
2491: Notes:
2492: See "petsc/include/petscdm.h" for available DM types (for instance, DM1D, DM2D, or DM3D).
2494: Level: intermediate
2496: .keywords: DM, set, type
2497: .seealso: DMGetType(), DMCreate()
2498: @*/
2499: PetscErrorCode DMSetType(DM dm, DMType method)
2500: {
2501: PetscErrorCode (*r)(DM);
2502: PetscBool match;
2507: PetscObjectTypeCompare((PetscObject) dm, method, &match);
2508: if (match) return(0);
2510: if (!DMRegisterAllCalled) {DMRegisterAll();}
2511: PetscFunctionListFind(DMList,method,&r);
2512: if (!r) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown DM type: %s", method);
2514: if (dm->ops->destroy) {
2515: (*dm->ops->destroy)(dm);
2516: dm->ops->destroy = NULL;
2517: }
2518: (*r)(dm);
2519: PetscObjectChangeTypeName((PetscObject)dm,method);
2520: return(0);
2521: }
2525: /*@C
2526: DMGetType - Gets the DM type name (as a string) from the DM.
2528: Not Collective
2530: Input Parameter:
2531: . dm - The DM
2533: Output Parameter:
2534: . type - The DM type name
2536: Level: intermediate
2538: .keywords: DM, get, type, name
2539: .seealso: DMSetType(), DMCreate()
2540: @*/
2541: PetscErrorCode DMGetType(DM dm, DMType *type)
2542: {
2548: if (!DMRegisterAllCalled) {
2549: DMRegisterAll();
2550: }
2551: *type = ((PetscObject)dm)->type_name;
2552: return(0);
2553: }
2557: /*@C
2558: DMConvert - Converts a DM to another DM, either of the same or different type.
2560: Collective on DM
2562: Input Parameters:
2563: + dm - the DM
2564: - newtype - new DM type (use "same" for the same type)
2566: Output Parameter:
2567: . M - pointer to new DM
2569: Notes:
2570: Cannot be used to convert a sequential DM to parallel or parallel to sequential,
2571: the MPI communicator of the generated DM is always the same as the communicator
2572: of the input DM.
2574: Level: intermediate
2576: .seealso: DMCreate()
2577: @*/
2578: PetscErrorCode DMConvert(DM dm, DMType newtype, DM *M)
2579: {
2580: DM B;
2581: char convname[256];
2582: PetscBool sametype, issame;
2589: PetscObjectTypeCompare((PetscObject) dm, newtype, &sametype);
2590: PetscStrcmp(newtype, "same", &issame);
2591: {
2592: PetscErrorCode (*conv)(DM, DMType, DM*) = NULL;
2594: /*
2595: Order of precedence:
2596: 1) See if a specialized converter is known to the current DM.
2597: 2) See if a specialized converter is known to the desired DM class.
2598: 3) See if a good general converter is registered for the desired class
2599: 4) See if a good general converter is known for the current matrix.
2600: 5) Use a really basic converter.
2601: */
2603: /* 1) See if a specialized converter is known to the current DM and the desired class */
2604: PetscStrcpy(convname,"DMConvert_");
2605: PetscStrcat(convname,((PetscObject) dm)->type_name);
2606: PetscStrcat(convname,"_");
2607: PetscStrcat(convname,newtype);
2608: PetscStrcat(convname,"_C");
2609: PetscObjectQueryFunction((PetscObject)dm,convname,&conv);
2610: if (conv) goto foundconv;
2612: /* 2) See if a specialized converter is known to the desired DM class. */
2613: DMCreate(PetscObjectComm((PetscObject)dm), &B);
2614: DMSetType(B, newtype);
2615: PetscStrcpy(convname,"DMConvert_");
2616: PetscStrcat(convname,((PetscObject) dm)->type_name);
2617: PetscStrcat(convname,"_");
2618: PetscStrcat(convname,newtype);
2619: PetscStrcat(convname,"_C");
2620: PetscObjectQueryFunction((PetscObject)B,convname,&conv);
2621: if (conv) {
2622: DMDestroy(&B);
2623: goto foundconv;
2624: }
2626: #if 0
2627: /* 3) See if a good general converter is registered for the desired class */
2628: conv = B->ops->convertfrom;
2629: DMDestroy(&B);
2630: if (conv) goto foundconv;
2632: /* 4) See if a good general converter is known for the current matrix */
2633: if (dm->ops->convert) {
2634: conv = dm->ops->convert;
2635: }
2636: if (conv) goto foundconv;
2637: #endif
2639: /* 5) Use a really basic converter. */
2640: SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "No conversion possible between DM types %s and %s", ((PetscObject) dm)->type_name, newtype);
2642: foundconv:
2643: PetscLogEventBegin(DM_Convert,dm,0,0,0);
2644: (*conv)(dm,newtype,M);
2645: PetscLogEventEnd(DM_Convert,dm,0,0,0);
2646: }
2647: PetscObjectStateIncrease((PetscObject) *M);
2648: return(0);
2649: }
2651: /*--------------------------------------------------------------------------------------------------------------------*/
2655: /*@C
2656: DMRegister - Adds a new DM component implementation
2658: Not Collective
2660: Input Parameters:
2661: + name - The name of a new user-defined creation routine
2662: - create_func - The creation routine itself
2664: Notes:
2665: DMRegister() may be called multiple times to add several user-defined DMs
2668: Sample usage:
2669: .vb
2670: DMRegister("my_da", MyDMCreate);
2671: .ve
2673: Then, your DM type can be chosen with the procedural interface via
2674: .vb
2675: DMCreate(MPI_Comm, DM *);
2676: DMSetType(DM,"my_da");
2677: .ve
2678: or at runtime via the option
2679: .vb
2680: -da_type my_da
2681: .ve
2683: Level: advanced
2685: .keywords: DM, register
2686: .seealso: DMRegisterAll(), DMRegisterDestroy()
2688: @*/
2689: PetscErrorCode DMRegister(const char sname[],PetscErrorCode (*function)(DM))
2690: {
2694: PetscFunctionListAdd(&DMList,sname,function);
2695: return(0);
2696: }
2700: /*@C
2701: DMLoad - Loads a DM that has been stored in binary with DMView().
2703: Collective on PetscViewer
2705: Input Parameters:
2706: + newdm - the newly loaded DM, this needs to have been created with DMCreate() or
2707: some related function before a call to DMLoad().
2708: - viewer - binary file viewer, obtained from PetscViewerBinaryOpen() or
2709: HDF5 file viewer, obtained from PetscViewerHDF5Open()
2711: Level: intermediate
2713: Notes:
2714: The type is determined by the data in the file, any type set into the DM before this call is ignored.
2716: Notes for advanced users:
2717: Most users should not need to know the details of the binary storage
2718: format, since DMLoad() and DMView() completely hide these details.
2719: But for anyone who's interested, the standard binary matrix storage
2720: format is
2721: .vb
2722: has not yet been determined
2723: .ve
2725: .seealso: PetscViewerBinaryOpen(), DMView(), MatLoad(), VecLoad()
2726: @*/
2727: PetscErrorCode DMLoad(DM newdm, PetscViewer viewer)
2728: {
2729: PetscBool isbinary, ishdf5;
2735: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
2736: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);
2737: if (isbinary) {
2738: PetscInt classid;
2739: char type[256];
2741: PetscViewerBinaryRead(viewer,&classid,1,PETSC_INT);
2742: if (classid != DM_FILE_CLASSID) SETERRQ1(PetscObjectComm((PetscObject)newdm),PETSC_ERR_ARG_WRONG,"Not DM next in file, classid found %d",(int)classid);
2743: PetscViewerBinaryRead(viewer,type,256,PETSC_CHAR);
2744: DMSetType(newdm, type);
2745: if (newdm->ops->load) {(*newdm->ops->load)(newdm,viewer);}
2746: } else if (ishdf5) {
2747: if (newdm->ops->load) {(*newdm->ops->load)(newdm,viewer);}
2748: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid viewer; open viewer with PetscViewerBinaryOpen() or PetscViewerHDF5Open()");
2749: return(0);
2750: }
2752: /******************************** FEM Support **********************************/
2756: PetscErrorCode DMPrintCellVector(PetscInt c, const char name[], PetscInt len, const PetscScalar x[])
2757: {
2758: PetscInt f;
2762: PetscPrintf(PETSC_COMM_SELF, "Cell %D Element %s\n", c, name);
2763: for (f = 0; f < len; ++f) {
2764: PetscPrintf(PETSC_COMM_SELF, " | %g |\n", (double)PetscRealPart(x[f]));
2765: }
2766: return(0);
2767: }
2771: PetscErrorCode DMPrintCellMatrix(PetscInt c, const char name[], PetscInt rows, PetscInt cols, const PetscScalar A[])
2772: {
2773: PetscInt f, g;
2777: PetscPrintf(PETSC_COMM_SELF, "Cell %D Element %s\n", c, name);
2778: for (f = 0; f < rows; ++f) {
2779: PetscPrintf(PETSC_COMM_SELF, " |");
2780: for (g = 0; g < cols; ++g) {
2781: PetscPrintf(PETSC_COMM_SELF, " % 9.5g", PetscRealPart(A[f*cols+g]));
2782: }
2783: PetscPrintf(PETSC_COMM_SELF, " |\n");
2784: }
2785: return(0);
2786: }
2790: PetscErrorCode DMPrintLocalVec(DM dm, const char name[], PetscReal tol, Vec X)
2791: {
2792: PetscMPIInt rank, numProcs;
2793: PetscInt p;
2797: MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
2798: MPI_Comm_size(PetscObjectComm((PetscObject) dm), &numProcs);
2799: PetscPrintf(PetscObjectComm((PetscObject) dm), "%s:\n", name);
2800: for (p = 0; p < numProcs; ++p) {
2801: if (p == rank) {
2802: Vec x;
2804: VecDuplicate(X, &x);
2805: VecCopy(X, x);
2806: VecChop(x, tol);
2807: VecView(x, PETSC_VIEWER_STDOUT_SELF);
2808: VecDestroy(&x);
2809: PetscViewerFlush(PETSC_VIEWER_STDOUT_SELF);
2810: }
2811: PetscBarrier((PetscObject) dm);
2812: }
2813: return(0);
2814: }
2818: /*@
2819: DMGetDefaultSection - Get the PetscSection encoding the local data layout for the DM.
2821: Input Parameter:
2822: . dm - The DM
2824: Output Parameter:
2825: . section - The PetscSection
2827: Level: intermediate
2829: Note: This gets a borrowed reference, so the user should not destroy this PetscSection.
2831: .seealso: DMSetDefaultSection(), DMGetDefaultGlobalSection()
2832: @*/
2833: PetscErrorCode DMGetDefaultSection(DM dm, PetscSection *section)
2834: {
2840: if (!dm->defaultSection && dm->ops->createdefaultsection) {(*dm->ops->createdefaultsection)(dm);}
2841: *section = dm->defaultSection;
2842: return(0);
2843: }
2847: /*@
2848: DMSetDefaultSection - Set the PetscSection encoding the local data layout for the DM.
2850: Input Parameters:
2851: + dm - The DM
2852: - section - The PetscSection
2854: Level: intermediate
2856: Note: Any existing Section will be destroyed
2858: .seealso: DMSetDefaultSection(), DMGetDefaultGlobalSection()
2859: @*/
2860: PetscErrorCode DMSetDefaultSection(DM dm, PetscSection section)
2861: {
2862: PetscInt numFields;
2863: PetscInt f;
2869: PetscObjectReference((PetscObject)section);
2870: PetscSectionDestroy(&dm->defaultSection);
2871: dm->defaultSection = section;
2872: PetscSectionGetNumFields(dm->defaultSection, &numFields);
2873: if (numFields) {
2874: DMSetNumFields(dm, numFields);
2875: for (f = 0; f < numFields; ++f) {
2876: PetscObject disc;
2877: const char *name;
2879: PetscSectionGetFieldName(dm->defaultSection, f, &name);
2880: DMGetField(dm, f, &disc);
2881: PetscObjectSetName(disc, name);
2882: }
2883: }
2884: /* The global section will be rebuilt in the next call to DMGetDefaultGlobalSection(). */
2885: PetscSectionDestroy(&dm->defaultGlobalSection);
2886: return(0);
2887: }
2891: /*@
2892: DMGetDefaultGlobalSection - Get the PetscSection encoding the global data layout for the DM.
2894: Collective on DM
2896: Input Parameter:
2897: . dm - The DM
2899: Output Parameter:
2900: . section - The PetscSection
2902: Level: intermediate
2904: Note: This gets a borrowed reference, so the user should not destroy this PetscSection.
2906: .seealso: DMSetDefaultSection(), DMGetDefaultSection()
2907: @*/
2908: PetscErrorCode DMGetDefaultGlobalSection(DM dm, PetscSection *section)
2909: {
2915: if (!dm->defaultGlobalSection) {
2916: PetscSection s;
2918: DMGetDefaultSection(dm, &s);
2919: if (!s) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONGSTATE, "DM must have a default PetscSection in order to create a global PetscSection");
2920: if (!dm->sf) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DM must have a default PetscSF in order to create a global PetscSection");
2921: PetscSectionCreateGlobalSection(s, dm->sf, PETSC_FALSE, &dm->defaultGlobalSection);
2922: PetscLayoutDestroy(&dm->map);
2923: PetscSectionGetValueLayout(PetscObjectComm((PetscObject)dm), dm->defaultGlobalSection, &dm->map);
2924: }
2925: *section = dm->defaultGlobalSection;
2926: return(0);
2927: }
2931: /*@
2932: DMSetDefaultGlobalSection - Set the PetscSection encoding the global data layout for the DM.
2934: Input Parameters:
2935: + dm - The DM
2936: - section - The PetscSection, or NULL
2938: Level: intermediate
2940: Note: Any existing Section will be destroyed
2942: .seealso: DMGetDefaultGlobalSection(), DMSetDefaultSection()
2943: @*/
2944: PetscErrorCode DMSetDefaultGlobalSection(DM dm, PetscSection section)
2945: {
2951: PetscObjectReference((PetscObject)section);
2952: PetscSectionDestroy(&dm->defaultGlobalSection);
2953: dm->defaultGlobalSection = section;
2954: return(0);
2955: }
2959: /*@
2960: DMGetDefaultSF - Get the PetscSF encoding the parallel dof overlap for the DM. If it has not been set,
2961: it is created from the default PetscSection layouts in the DM.
2963: Input Parameter:
2964: . dm - The DM
2966: Output Parameter:
2967: . sf - The PetscSF
2969: Level: intermediate
2971: Note: This gets a borrowed reference, so the user should not destroy this PetscSF.
2973: .seealso: DMSetDefaultSF(), DMCreateDefaultSF()
2974: @*/
2975: PetscErrorCode DMGetDefaultSF(DM dm, PetscSF *sf)
2976: {
2977: PetscInt nroots;
2983: PetscSFGetGraph(dm->defaultSF, &nroots, NULL, NULL, NULL);
2984: if (nroots < 0) {
2985: PetscSection section, gSection;
2987: DMGetDefaultSection(dm, §ion);
2988: if (section) {
2989: DMGetDefaultGlobalSection(dm, &gSection);
2990: DMCreateDefaultSF(dm, section, gSection);
2991: } else {
2992: *sf = NULL;
2993: return(0);
2994: }
2995: }
2996: *sf = dm->defaultSF;
2997: return(0);
2998: }
3002: /*@
3003: DMSetDefaultSF - Set the PetscSF encoding the parallel dof overlap for the DM
3005: Input Parameters:
3006: + dm - The DM
3007: - sf - The PetscSF
3009: Level: intermediate
3011: Note: Any previous SF is destroyed
3013: .seealso: DMGetDefaultSF(), DMCreateDefaultSF()
3014: @*/
3015: PetscErrorCode DMSetDefaultSF(DM dm, PetscSF sf)
3016: {
3022: PetscSFDestroy(&dm->defaultSF);
3023: dm->defaultSF = sf;
3024: return(0);
3025: }
3029: /*@C
3030: DMCreateDefaultSF - Create the PetscSF encoding the parallel dof overlap for the DM based upon the PetscSections
3031: describing the data layout.
3033: Input Parameters:
3034: + dm - The DM
3035: . localSection - PetscSection describing the local data layout
3036: - globalSection - PetscSection describing the global data layout
3038: Level: intermediate
3040: .seealso: DMGetDefaultSF(), DMSetDefaultSF()
3041: @*/
3042: PetscErrorCode DMCreateDefaultSF(DM dm, PetscSection localSection, PetscSection globalSection)
3043: {
3044: MPI_Comm comm;
3045: PetscLayout layout;
3046: const PetscInt *ranges;
3047: PetscInt *local;
3048: PetscSFNode *remote;
3049: PetscInt pStart, pEnd, p, nroots, nleaves = 0, l;
3050: PetscMPIInt size, rank;
3054: PetscObjectGetComm((PetscObject)dm,&comm);
3056: MPI_Comm_size(comm, &size);
3057: MPI_Comm_rank(comm, &rank);
3058: PetscSectionGetChart(globalSection, &pStart, &pEnd);
3059: PetscSectionGetConstrainedStorageSize(globalSection, &nroots);
3060: PetscLayoutCreate(comm, &layout);
3061: PetscLayoutSetBlockSize(layout, 1);
3062: PetscLayoutSetLocalSize(layout, nroots);
3063: PetscLayoutSetUp(layout);
3064: PetscLayoutGetRanges(layout, &ranges);
3065: for (p = pStart; p < pEnd; ++p) {
3066: PetscInt gdof, gcdof;
3068: PetscSectionGetDof(globalSection, p, &gdof);
3069: PetscSectionGetConstraintDof(globalSection, p, &gcdof);
3070: 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));
3071: nleaves += gdof < 0 ? -(gdof+1)-gcdof : gdof-gcdof;
3072: }
3073: PetscMalloc1(nleaves, &local);
3074: PetscMalloc1(nleaves, &remote);
3075: for (p = pStart, l = 0; p < pEnd; ++p) {
3076: const PetscInt *cind;
3077: PetscInt dof, cdof, off, gdof, gcdof, goff, gsize, d, c;
3079: PetscSectionGetDof(localSection, p, &dof);
3080: PetscSectionGetOffset(localSection, p, &off);
3081: PetscSectionGetConstraintDof(localSection, p, &cdof);
3082: PetscSectionGetConstraintIndices(localSection, p, &cind);
3083: PetscSectionGetDof(globalSection, p, &gdof);
3084: PetscSectionGetConstraintDof(globalSection, p, &gcdof);
3085: PetscSectionGetOffset(globalSection, p, &goff);
3086: if (!gdof) continue; /* Censored point */
3087: gsize = gdof < 0 ? -(gdof+1)-gcdof : gdof-gcdof;
3088: if (gsize != dof-cdof) {
3089: 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);
3090: cdof = 0; /* Ignore constraints */
3091: }
3092: for (d = 0, c = 0; d < dof; ++d) {
3093: if ((c < cdof) && (cind[c] == d)) {++c; continue;}
3094: local[l+d-c] = off+d;
3095: }
3096: if (gdof < 0) {
3097: for (d = 0; d < gsize; ++d, ++l) {
3098: PetscInt offset = -(goff+1) + d, r;
3100: PetscFindInt(offset,size+1,ranges,&r);
3101: if (r < 0) r = -(r+2);
3102: 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);
3103: remote[l].rank = r;
3104: remote[l].index = offset - ranges[r];
3105: }
3106: } else {
3107: for (d = 0; d < gsize; ++d, ++l) {
3108: remote[l].rank = rank;
3109: remote[l].index = goff+d - ranges[rank];
3110: }
3111: }
3112: }
3113: if (l != nleaves) SETERRQ2(comm, PETSC_ERR_PLIB, "Iteration error, l %d != nleaves %d", l, nleaves);
3114: PetscLayoutDestroy(&layout);
3115: PetscSFSetGraph(dm->defaultSF, nroots, nleaves, local, PETSC_OWN_POINTER, remote, PETSC_OWN_POINTER);
3116: return(0);
3117: }
3121: /*@
3122: DMGetPointSF - Get the PetscSF encoding the parallel section point overlap for the DM.
3124: Input Parameter:
3125: . dm - The DM
3127: Output Parameter:
3128: . sf - The PetscSF
3130: Level: intermediate
3132: Note: This gets a borrowed reference, so the user should not destroy this PetscSF.
3134: .seealso: DMSetPointSF(), DMGetDefaultSF(), DMSetDefaultSF(), DMCreateDefaultSF()
3135: @*/
3136: PetscErrorCode DMGetPointSF(DM dm, PetscSF *sf)
3137: {
3141: *sf = dm->sf;
3142: return(0);
3143: }
3147: /*@
3148: DMSetPointSF - Set the PetscSF encoding the parallel section point overlap for the DM.
3150: Input Parameters:
3151: + dm - The DM
3152: - sf - The PetscSF
3154: Level: intermediate
3156: .seealso: DMGetPointSF(), DMGetDefaultSF(), DMSetDefaultSF(), DMCreateDefaultSF()
3157: @*/
3158: PetscErrorCode DMSetPointSF(DM dm, PetscSF sf)
3159: {
3165: PetscSFDestroy(&dm->sf);
3166: PetscObjectReference((PetscObject) sf);
3167: dm->sf = sf;
3168: return(0);
3169: }
3173: /*@
3174: DMGetDS - Get the PetscDS
3176: Input Parameter:
3177: . dm - The DM
3179: Output Parameter:
3180: . prob - The PetscDS
3182: Level: developer
3184: .seealso: DMSetDS()
3185: @*/
3186: PetscErrorCode DMGetDS(DM dm, PetscDS *prob)
3187: {
3191: *prob = dm->prob;
3192: return(0);
3193: }
3197: /*@
3198: DMSetDS - Set the PetscDS
3200: Input Parameters:
3201: + dm - The DM
3202: - prob - The PetscDS
3204: Level: developer
3206: .seealso: DMGetDS()
3207: @*/
3208: PetscErrorCode DMSetDS(DM dm, PetscDS prob)
3209: {
3215: PetscDSDestroy(&dm->prob);
3216: dm->prob = prob;
3217: PetscObjectReference((PetscObject) dm->prob);
3218: return(0);
3219: }
3223: PetscErrorCode DMGetNumFields(DM dm, PetscInt *numFields)
3224: {
3229: PetscDSGetNumFields(dm->prob, numFields);
3230: return(0);
3231: }
3235: PetscErrorCode DMSetNumFields(DM dm, PetscInt numFields)
3236: {
3237: PetscInt Nf, f;
3242: PetscDSGetNumFields(dm->prob, &Nf);
3243: for (f = Nf; f < numFields; ++f) {
3244: PetscContainer obj;
3246: PetscContainerCreate(PetscObjectComm((PetscObject) dm), &obj);
3247: PetscDSSetDiscretization(dm->prob, f, (PetscObject) obj);
3248: PetscContainerDestroy(&obj);
3249: }
3250: return(0);
3251: }
3255: /*@
3256: DMGetField - Return the discretization object for a given DM field
3258: Not collective
3260: Input Parameters:
3261: + dm - The DM
3262: - f - The field number
3264: Output Parameter:
3265: . field - The discretization object
3267: Level: developer
3269: .seealso: DMSetField()
3270: @*/
3271: PetscErrorCode DMGetField(DM dm, PetscInt f, PetscObject *field)
3272: {
3277: PetscDSGetDiscretization(dm->prob, f, field);
3278: return(0);
3279: }
3283: /*@
3284: DMSetField - Set the discretization object for a given DM field
3286: Logically collective on DM
3288: Input Parameters:
3289: + dm - The DM
3290: . f - The field number
3291: - field - The discretization object
3293: Level: developer
3295: .seealso: DMGetField()
3296: @*/
3297: PetscErrorCode DMSetField(DM dm, PetscInt f, PetscObject field)
3298: {
3303: PetscDSSetDiscretization(dm->prob, f, field);
3304: return(0);
3305: }
3309: PetscErrorCode DMRestrictHook_Coordinates(DM dm,DM dmc,void *ctx)
3310: {
3311: DM dm_coord,dmc_coord;
3313: Vec coords,ccoords;
3314: VecScatter scat;
3316: DMGetCoordinateDM(dm,&dm_coord);
3317: DMGetCoordinateDM(dmc,&dmc_coord);
3318: DMGetCoordinates(dm,&coords);
3319: DMGetCoordinates(dmc,&ccoords);
3320: if (coords && !ccoords) {
3321: DMCreateGlobalVector(dmc_coord,&ccoords);
3322: DMCreateInjection(dmc_coord,dm_coord,&scat);
3323: VecScatterBegin(scat,coords,ccoords,INSERT_VALUES,SCATTER_FORWARD);
3324: VecScatterEnd(scat,coords,ccoords,INSERT_VALUES,SCATTER_FORWARD);
3325: DMSetCoordinates(dmc,ccoords);
3326: VecScatterDestroy(&scat);
3327: VecDestroy(&ccoords);
3328: }
3329: return(0);
3330: }
3334: static PetscErrorCode DMSubDomainHook_Coordinates(DM dm,DM subdm,void *ctx)
3335: {
3336: DM dm_coord,subdm_coord;
3338: Vec coords,ccoords,clcoords;
3339: VecScatter *scat_i,*scat_g;
3341: DMGetCoordinateDM(dm,&dm_coord);
3342: DMGetCoordinateDM(subdm,&subdm_coord);
3343: DMGetCoordinates(dm,&coords);
3344: DMGetCoordinates(subdm,&ccoords);
3345: if (coords && !ccoords) {
3346: DMCreateGlobalVector(subdm_coord,&ccoords);
3347: DMCreateLocalVector(subdm_coord,&clcoords);
3348: DMCreateDomainDecompositionScatters(dm_coord,1,&subdm_coord,NULL,&scat_i,&scat_g);
3349: VecScatterBegin(scat_i[0],coords,ccoords,INSERT_VALUES,SCATTER_FORWARD);
3350: VecScatterBegin(scat_g[0],coords,clcoords,INSERT_VALUES,SCATTER_FORWARD);
3351: VecScatterEnd(scat_i[0],coords,ccoords,INSERT_VALUES,SCATTER_FORWARD);
3352: VecScatterEnd(scat_g[0],coords,clcoords,INSERT_VALUES,SCATTER_FORWARD);
3353: DMSetCoordinates(subdm,ccoords);
3354: DMSetCoordinatesLocal(subdm,clcoords);
3355: VecScatterDestroy(&scat_i[0]);
3356: VecScatterDestroy(&scat_g[0]);
3357: VecDestroy(&ccoords);
3358: VecDestroy(&clcoords);
3359: PetscFree(scat_i);
3360: PetscFree(scat_g);
3361: }
3362: return(0);
3363: }
3367: /*@
3368: DMSetCoordinates - Sets into the DM a global vector that holds the coordinates
3370: Collective on DM
3372: Input Parameters:
3373: + dm - the DM
3374: - c - coordinate vector
3376: Note:
3377: The coordinates do include those for ghost points, which are in the local vector
3379: Level: intermediate
3381: .keywords: distributed array, get, corners, nodes, local indices, coordinates
3382: .seealso: DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLoca(), DMGetCoordinateDM()
3383: @*/
3384: PetscErrorCode DMSetCoordinates(DM dm, Vec c)
3385: {
3391: PetscObjectReference((PetscObject) c);
3392: VecDestroy(&dm->coordinates);
3393: dm->coordinates = c;
3394: VecDestroy(&dm->coordinatesLocal);
3395: DMCoarsenHookAdd(dm,DMRestrictHook_Coordinates,NULL,NULL);
3396: DMSubDomainHookAdd(dm,DMSubDomainHook_Coordinates,NULL,NULL);
3397: return(0);
3398: }
3402: /*@
3403: DMSetCoordinatesLocal - Sets into the DM a local vector that holds the coordinates
3405: Collective on DM
3407: Input Parameters:
3408: + dm - the DM
3409: - c - coordinate vector
3411: Note:
3412: The coordinates of ghost points can be set using DMSetCoordinates()
3413: followed by DMGetCoordinatesLocal(). This is intended to enable the
3414: setting of ghost coordinates outside of the domain.
3416: Level: intermediate
3418: .keywords: distributed array, get, corners, nodes, local indices, coordinates
3419: .seealso: DMGetCoordinatesLocal(), DMSetCoordinates(), DMGetCoordinates(), DMGetCoordinateDM()
3420: @*/
3421: PetscErrorCode DMSetCoordinatesLocal(DM dm, Vec c)
3422: {
3428: PetscObjectReference((PetscObject) c);
3429: VecDestroy(&dm->coordinatesLocal);
3431: dm->coordinatesLocal = c;
3433: VecDestroy(&dm->coordinates);
3434: return(0);
3435: }
3439: /*@
3440: DMGetCoordinates - Gets a global vector with the coordinates associated with the DM.
3442: Not Collective
3444: Input Parameter:
3445: . dm - the DM
3447: Output Parameter:
3448: . c - global coordinate vector
3450: Note:
3451: This is a borrowed reference, so the user should NOT destroy this vector
3453: Each process has only the local coordinates (does NOT have the ghost coordinates).
3455: For DMDA, in two and three dimensions coordinates are interlaced (x_0,y_0,x_1,y_1,...)
3456: and (x_0,y_0,z_0,x_1,y_1,z_1...)
3458: Level: intermediate
3460: .keywords: distributed array, get, corners, nodes, local indices, coordinates
3461: .seealso: DMSetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM()
3462: @*/
3463: PetscErrorCode DMGetCoordinates(DM dm, Vec *c)
3464: {
3470: if (!dm->coordinates && dm->coordinatesLocal) {
3471: DM cdm = NULL;
3473: DMGetCoordinateDM(dm, &cdm);
3474: DMCreateGlobalVector(cdm, &dm->coordinates);
3475: PetscObjectSetName((PetscObject) dm->coordinates, "coordinates");
3476: DMLocalToGlobalBegin(cdm, dm->coordinatesLocal, INSERT_VALUES, dm->coordinates);
3477: DMLocalToGlobalEnd(cdm, dm->coordinatesLocal, INSERT_VALUES, dm->coordinates);
3478: }
3479: *c = dm->coordinates;
3480: return(0);
3481: }
3485: /*@
3486: DMGetCoordinatesLocal - Gets a local vector with the coordinates associated with the DM.
3488: Collective on DM
3490: Input Parameter:
3491: . dm - the DM
3493: Output Parameter:
3494: . c - coordinate vector
3496: Note:
3497: This is a borrowed reference, so the user should NOT destroy this vector
3499: Each process has the local and ghost coordinates
3501: For DMDA, in two and three dimensions coordinates are interlaced (x_0,y_0,x_1,y_1,...)
3502: and (x_0,y_0,z_0,x_1,y_1,z_1...)
3504: Level: intermediate
3506: .keywords: distributed array, get, corners, nodes, local indices, coordinates
3507: .seealso: DMSetCoordinatesLocal(), DMGetCoordinates(), DMSetCoordinates(), DMGetCoordinateDM()
3508: @*/
3509: PetscErrorCode DMGetCoordinatesLocal(DM dm, Vec *c)
3510: {
3516: if (!dm->coordinatesLocal && dm->coordinates) {
3517: DM cdm = NULL;
3519: DMGetCoordinateDM(dm, &cdm);
3520: DMCreateLocalVector(cdm, &dm->coordinatesLocal);
3521: PetscObjectSetName((PetscObject) dm->coordinatesLocal, "coordinates");
3522: DMGlobalToLocalBegin(cdm, dm->coordinates, INSERT_VALUES, dm->coordinatesLocal);
3523: DMGlobalToLocalEnd(cdm, dm->coordinates, INSERT_VALUES, dm->coordinatesLocal);
3524: }
3525: *c = dm->coordinatesLocal;
3526: return(0);
3527: }
3531: /*@
3532: DMGetCoordinateDM - Gets the DM that prescribes coordinate layout and scatters between global and local coordinates
3534: Collective on DM
3536: Input Parameter:
3537: . dm - the DM
3539: Output Parameter:
3540: . cdm - coordinate DM
3542: Level: intermediate
3544: .keywords: distributed array, get, corners, nodes, local indices, coordinates
3545: .seealso: DMSetCoordinateDM(), DMSetCoordinates(), DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLocal()
3546: @*/
3547: PetscErrorCode DMGetCoordinateDM(DM dm, DM *cdm)
3548: {
3554: if (!dm->coordinateDM) {
3555: if (!dm->ops->createcoordinatedm) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unable to create coordinates for this DM");
3556: (*dm->ops->createcoordinatedm)(dm, &dm->coordinateDM);
3557: }
3558: *cdm = dm->coordinateDM;
3559: return(0);
3560: }
3564: /*@
3565: DMSetCoordinateDM - Sets the DM that prescribes coordinate layout and scatters between global and local coordinates
3567: Logically Collective on DM
3569: Input Parameters:
3570: + dm - the DM
3571: - cdm - coordinate DM
3573: Level: intermediate
3575: .keywords: distributed array, get, corners, nodes, local indices, coordinates
3576: .seealso: DMGetCoordinateDM(), DMSetCoordinates(), DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLocal()
3577: @*/
3578: PetscErrorCode DMSetCoordinateDM(DM dm, DM cdm)
3579: {
3585: DMDestroy(&dm->coordinateDM);
3586: dm->coordinateDM = cdm;
3587: PetscObjectReference((PetscObject) dm->coordinateDM);
3588: return(0);
3589: }
3593: /*@
3594: DMGetCoordinateSection - Retrieve the layout of coordinate values over the mesh.
3596: Not Collective
3598: Input Parameter:
3599: . dm - The DM object
3601: Output Parameter:
3602: . section - The PetscSection object
3604: Level: intermediate
3606: .keywords: mesh, coordinates
3607: .seealso: DMGetCoordinateDM(), DMGetDefaultSection(), DMSetDefaultSection()
3608: @*/
3609: PetscErrorCode DMGetCoordinateSection(DM dm, PetscSection *section)
3610: {
3611: DM cdm;
3617: DMGetCoordinateDM(dm, &cdm);
3618: DMGetDefaultSection(cdm, section);
3619: return(0);
3620: }
3624: /*@
3625: DMSetCoordinateSection - Set the layout of coordinate values over the mesh.
3627: Not Collective
3629: Input Parameters:
3630: + dm - The DM object
3631: - section - The PetscSection object
3633: Level: intermediate
3635: .keywords: mesh, coordinates
3636: .seealso: DMGetCoordinateSection(), DMGetDefaultSection(), DMSetDefaultSection()
3637: @*/
3638: PetscErrorCode DMSetCoordinateSection(DM dm, PetscSection section)
3639: {
3640: DM cdm;
3646: DMGetCoordinateDM(dm, &cdm);
3647: DMSetDefaultSection(cdm, section);
3648: return(0);
3649: }
3653: PetscErrorCode DMGetPeriodicity(DM dm, const PetscReal **maxCell, const PetscReal **L)
3654: {
3657: if (L) *L = dm->L;
3658: if (maxCell) *maxCell = dm->maxCell;
3659: return(0);
3660: }
3664: PetscErrorCode DMSetPeriodicity(DM dm, const PetscReal maxCell[], const PetscReal L[])
3665: {
3666: Vec coordinates;
3667: PetscInt dim, d;
3672: DMGetCoordinatesLocal(dm, &coordinates);
3673: VecGetBlockSize(coordinates, &dim);
3674: PetscMalloc1(dim,&dm->L);
3675: PetscMalloc1(dim,&dm->maxCell);
3676: for (d = 0; d < dim; ++d) {dm->L[d] = L[d]; dm->maxCell[d] = maxCell[d];}
3677: return(0);
3678: }
3682: /*@
3683: DMLocatePoints - Locate the points in v in the mesh and return an IS of the containing cells
3685: Not collective
3687: Input Parameters:
3688: + dm - The DM
3689: - v - The Vec of points
3691: Output Parameter:
3692: . cells - The local cell numbers for cells which contain the points
3694: Level: developer
3696: .keywords: point location, mesh
3697: .seealso: DMSetCoordinates(), DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLocal()
3698: @*/
3699: PetscErrorCode DMLocatePoints(DM dm, Vec v, IS *cells)
3700: {
3707: if (dm->ops->locatepoints) {
3708: (*dm->ops->locatepoints)(dm,v,cells);
3709: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Point location not available for this DM");
3710: return(0);
3711: }
3715: /*@
3716: DMGetOutputDM - Retrieve the DM associated with the layout for output
3718: Input Parameter:
3719: . dm - The original DM
3721: Output Parameter:
3722: . odm - The DM which provides the layout for output
3724: Level: intermediate
3726: .seealso: VecView(), DMGetDefaultGlobalSection()
3727: @*/
3728: PetscErrorCode DMGetOutputDM(DM dm, DM *odm)
3729: {
3730: PetscSection section;
3731: PetscBool hasConstraints;
3737: DMGetDefaultSection(dm, §ion);
3738: PetscSectionHasConstraints(section, &hasConstraints);
3739: if (!hasConstraints) {
3740: *odm = dm;
3741: return(0);
3742: }
3743: if (!dm->dmBC) {
3744: PetscSection newSection, gsection;
3745: PetscSF sf;
3747: DMClone(dm, &dm->dmBC);
3748: PetscSectionClone(section, &newSection);
3749: DMSetDefaultSection(dm->dmBC, newSection);
3750: PetscSectionDestroy(&newSection);
3751: DMGetPointSF(dm->dmBC, &sf);
3752: PetscSectionCreateGlobalSection(section, sf, PETSC_TRUE, &gsection);
3753: DMSetDefaultGlobalSection(dm->dmBC, gsection);
3754: PetscSectionDestroy(&gsection);
3755: }
3756: *odm = dm->dmBC;
3757: return(0);
3758: }
3762: /*@
3763: DMGetOutputSequenceNumber - Retrieve the sequence number/value for output
3765: Input Parameter:
3766: . dm - The original DM
3768: Output Parameters:
3769: + num - The output sequence number
3770: - val - The output sequence value
3772: Level: intermediate
3774: Note: This is intended for output that should appear in sequence, for instance
3775: a set of timesteps in an HDF5 file, or a set of realizations of a stochastic system.
3777: .seealso: VecView()
3778: @*/
3779: PetscErrorCode DMGetOutputSequenceNumber(DM dm, PetscInt *num, PetscReal *val)
3780: {
3785: return(0);
3786: }
3790: /*@
3791: DMSetOutputSequenceNumber - Set the sequence number/value for output
3793: Input Parameters:
3794: + dm - The original DM
3795: . num - The output sequence number
3796: - val - The output sequence value
3798: Level: intermediate
3800: Note: This is intended for output that should appear in sequence, for instance
3801: a set of timesteps in an HDF5 file, or a set of realizations of a stochastic system.
3803: .seealso: VecView()
3804: @*/
3805: PetscErrorCode DMSetOutputSequenceNumber(DM dm, PetscInt num, PetscReal val)
3806: {
3809: dm->outputSequenceNum = num;
3810: dm->outputSequenceVal = val;
3811: return(0);
3812: }
3816: /*@C
3817: DMOutputSequenceLoad - Retrieve the sequence value from a Viewer
3819: Input Parameters:
3820: + dm - The original DM
3821: . name - The sequence name
3822: - num - The output sequence number
3824: Output Parameter:
3825: . val - The output sequence value
3827: Level: intermediate
3829: Note: This is intended for output that should appear in sequence, for instance
3830: a set of timesteps in an HDF5 file, or a set of realizations of a stochastic system.
3832: .seealso: DMGetOutputSequenceNumber(), DMSetOutputSequenceNumber(), VecView()
3833: @*/
3834: PetscErrorCode DMOutputSequenceLoad(DM dm, PetscViewer viewer, const char *name, PetscInt num, PetscReal *val)
3835: {
3836: PetscBool ishdf5;
3843: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERHDF5, &ishdf5);
3844: if (ishdf5) {
3845: #if defined(PETSC_HAVE_HDF5)
3846: PetscScalar value;
3848: DMSequenceLoad_HDF5(dm, name, num, &value, viewer);
3849: *val = PetscRealPart(value);
3850: #endif
3851: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid viewer; open viewer with PetscViewerHDF5Open()");
3852: return(0);
3853: }