Actual source code: dm.c
petsc-3.4.2 2013-07-02
1: #include <petsc-private/dmimpl.h> /*I "petscdm.h" I*/
2: #include <petscsf.h>
4: PetscClassId DM_CLASSID;
5: PetscLogEvent DM_Convert, DM_GlobalToLocal, DM_LocalToGlobal;
9: /*
10: DMViewFromOptions - Processes command line options to determine if/how a DM is to be viewed.
12: Collective on Vec
14: Input Parameters:
15: + dm - the DM
16: . prefix - prefix to use for viewing, or NULL to use prefix of 'rnd'
17: - optionname - option to activate viewing
19: Level: intermediate
21: .keywords: DM, view, options, database
22: .seealso: VecViewFromOptions(), MatViewFromOptions()
23: */
24: PetscErrorCode DMViewFromOptions(DM dm,const char prefix[],const char optionname[])
25: {
26: PetscErrorCode ierr;
27: PetscBool flg;
28: PetscViewer viewer;
29: PetscViewerFormat format;
32: if (prefix) {
33: PetscOptionsGetViewer(PetscObjectComm((PetscObject)dm),prefix,optionname,&viewer,&format,&flg);
34: } else {
35: PetscOptionsGetViewer(PetscObjectComm((PetscObject)dm),((PetscObject)dm)->prefix,optionname,&viewer,&format,&flg);
36: }
37: if (flg) {
38: PetscViewerPushFormat(viewer,format);
39: DMView(dm,viewer);
40: PetscViewerPopFormat(viewer);
41: PetscViewerDestroy(&viewer);
42: }
43: return(0);
44: }
49: /*@
50: DMCreate - Creates an empty DM object. The type can then be set with DMSetType().
52: If you never call DMSetType() it will generate an
53: error when you try to use the vector.
55: Collective on MPI_Comm
57: Input Parameter:
58: . comm - The communicator for the DM object
60: Output Parameter:
61: . dm - The DM object
63: Level: beginner
65: .seealso: DMSetType(), DMDA, DMSLICED, DMCOMPOSITE
66: @*/
67: PetscErrorCode DMCreate(MPI_Comm comm,DM *dm)
68: {
69: DM v;
74: *dm = NULL;
75: #if !defined(PETSC_USE_DYNAMIC_LIBRARIES)
76: VecInitializePackage();
77: MatInitializePackage();
78: DMInitializePackage();
79: #endif
81: PetscHeaderCreate(v, _p_DM, struct _DMOps, DM_CLASSID, "DM", "Distribution Manager", "DM", comm, DMDestroy, DMView);
82: PetscMemzero(v->ops, sizeof(struct _DMOps));
85: v->ltogmap = NULL;
86: v->ltogmapb = NULL;
87: v->bs = 1;
88: v->coloringtype = IS_COLORING_GLOBAL;
89: PetscSFCreate(comm, &v->sf);
90: PetscSFCreate(comm, &v->defaultSF);
91: v->defaultSection = NULL;
92: v->defaultGlobalSection = NULL;
93: {
94: PetscInt i;
95: for (i = 0; i < 10; ++i) {
96: v->nullspaceConstructors[i] = NULL;
97: }
98: }
99: v->numFields = 0;
100: v->fields = NULL;
102: *dm = v;
103: return(0);
104: }
108: /*@C
109: DMSetVecType - Sets the type of vector created with DMCreateLocalVector() and DMCreateGlobalVector()
111: Logically Collective on DMDA
113: Input Parameter:
114: + da - initial distributed array
115: . ctype - the vector type, currently either VECSTANDARD or VECCUSP
117: Options Database:
118: . -dm_vec_type ctype
120: Level: intermediate
122: .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMDestroy(), DMDA, DMDAInterpolationType, VecType
123: @*/
124: PetscErrorCode DMSetVecType(DM da,VecType ctype)
125: {
130: PetscFree(da->vectype);
131: PetscStrallocpy(ctype,(char**)&da->vectype);
132: return(0);
133: }
137: /*@
138: VecGetDM - Gets the DM defining the data layout of the vector
140: Not collective
142: Input Parameter:
143: . v - The Vec
145: Output Parameter:
146: . dm - The DM
148: Level: intermediate
150: .seealso: VecSetDM(), DMGetLocalVector(), DMGetGlobalVector(), DMSetVecType()
151: @*/
152: PetscErrorCode VecGetDM(Vec v, DM *dm)
153: {
159: PetscObjectQuery((PetscObject) v, "__PETSc_dm", (PetscObject*) dm);
160: return(0);
161: }
165: /*@
166: VecSetDM - Sets the DM defining the data layout of the vector
168: Not collective
170: Input Parameters:
171: + v - The Vec
172: - dm - The DM
174: Level: intermediate
176: .seealso: VecGetDM(), DMGetLocalVector(), DMGetGlobalVector(), DMSetVecType()
177: @*/
178: PetscErrorCode VecSetDM(Vec v, DM dm)
179: {
185: PetscObjectCompose((PetscObject) v, "__PETSc_dm", (PetscObject) dm);
186: return(0);
187: }
191: /*@C
192: DMSetMatType - Sets the type of matrix created with DMCreateMatrix()
194: Logically Collective on DM
196: Input Parameter:
197: + dm - the DM context
198: . ctype - the matrix type
200: Options Database:
201: . -dm_mat_type ctype
203: Level: intermediate
205: .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMCreateMatrix(), DMSetMatrixPreallocateOnly(), MatType
206: @*/
207: PetscErrorCode DMSetMatType(DM dm,MatType ctype)
208: {
213: PetscFree(dm->mattype);
214: PetscStrallocpy(ctype,(char**)&dm->mattype);
215: return(0);
216: }
220: /*@
221: MatGetDM - Gets the DM defining the data layout of the matrix
223: Not collective
225: Input Parameter:
226: . A - The Mat
228: Output Parameter:
229: . dm - The DM
231: Level: intermediate
233: .seealso: MatSetDM(), DMCreateMatrix(), DMSetMatType()
234: @*/
235: PetscErrorCode MatGetDM(Mat A, DM *dm)
236: {
242: PetscObjectQuery((PetscObject) A, "__PETSc_dm", (PetscObject*) dm);
243: return(0);
244: }
248: /*@
249: MatSetDM - Sets the DM defining the data layout of the matrix
251: Not collective
253: Input Parameters:
254: + A - The Mat
255: - dm - The DM
257: Level: intermediate
259: .seealso: MatGetDM(), DMCreateMatrix(), DMSetMatType()
260: @*/
261: PetscErrorCode MatSetDM(Mat A, DM dm)
262: {
268: PetscObjectCompose((PetscObject) A, "__PETSc_dm", (PetscObject) dm);
269: return(0);
270: }
274: /*@C
275: DMSetOptionsPrefix - Sets the prefix used for searching for all
276: DMDA options in the database.
278: Logically Collective on DMDA
280: Input Parameter:
281: + da - the DMDA context
282: - prefix - the prefix to prepend to all option names
284: Notes:
285: A hyphen (-) must NOT be given at the beginning of the prefix name.
286: The first character of all runtime options is AUTOMATICALLY the hyphen.
288: Level: advanced
290: .keywords: DMDA, set, options, prefix, database
292: .seealso: DMSetFromOptions()
293: @*/
294: PetscErrorCode DMSetOptionsPrefix(DM dm,const char prefix[])
295: {
300: PetscObjectSetOptionsPrefix((PetscObject)dm,prefix);
301: return(0);
302: }
306: /*@
307: DMDestroy - Destroys a vector packer or DMDA.
309: Collective on DM
311: Input Parameter:
312: . dm - the DM object to destroy
314: Level: developer
316: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix()
318: @*/
319: PetscErrorCode DMDestroy(DM *dm)
320: {
321: PetscInt i, cnt = 0, f;
322: DMNamedVecLink nlink,nnext;
326: if (!*dm) return(0);
329: /* count all the circular references of DM and its contained Vecs */
330: for (i=0; i<DM_MAX_WORK_VECTORS; i++) {
331: if ((*dm)->localin[i]) cnt++;
332: if ((*dm)->globalin[i]) cnt++;
333: }
334: for (nlink=(*dm)->namedglobal; nlink; nlink=nlink->next) cnt++;
335: if ((*dm)->x) {
336: DM obj;
337: VecGetDM((*dm)->x, &obj);
338: if (obj == *dm) cnt++;
339: }
340: for (nlink=(*dm)->namedlocal; nlink; nlink=nlink->next) cnt++;
341: if ((*dm)->x) {
342: DM obj;
343: VecGetDM((*dm)->x, &obj);
344: if (obj == *dm) cnt++;
345: }
347: if (--((PetscObject)(*dm))->refct - cnt > 0) {*dm = 0; return(0);}
348: /*
349: Need this test because the dm references the vectors that
350: reference the dm, so destroying the dm calls destroy on the
351: vectors that cause another destroy on the dm
352: */
353: if (((PetscObject)(*dm))->refct < 0) return(0);
354: ((PetscObject) (*dm))->refct = 0;
355: for (i=0; i<DM_MAX_WORK_VECTORS; i++) {
356: if ((*dm)->localout[i]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Destroying a DM that has a local vector obtained with DMGetLocalVector()");
357: VecDestroy(&(*dm)->localin[i]);
358: }
359: for (nlink=(*dm)->namedglobal; nlink; nlink=nnext) { /* Destroy the named vectors */
360: nnext = nlink->next;
361: if (nlink->status != DMVEC_STATUS_IN) SETERRQ1(((PetscObject)*dm)->comm,PETSC_ERR_ARG_WRONGSTATE,"DM still has Vec named '%s' checked out",nlink->name);
362: PetscFree(nlink->name);
363: VecDestroy(&nlink->X);
364: PetscFree(nlink);
365: }
366: (*dm)->namedglobal = NULL;
368: for (nlink=(*dm)->namedlocal; nlink; nlink=nnext) { /* Destroy the named local vectors */
369: nnext = nlink->next;
370: if (nlink->status != DMVEC_STATUS_IN) SETERRQ1(((PetscObject)*dm)->comm,PETSC_ERR_ARG_WRONGSTATE,"DM still has Vec named '%s' checked out",nlink->name);
371: PetscFree(nlink->name);
372: VecDestroy(&nlink->X);
373: PetscFree(nlink);
374: }
375: (*dm)->namedlocal = NULL;
377: /* Destroy the list of hooks */
378: {
379: DMCoarsenHookLink link,next;
380: for (link=(*dm)->coarsenhook; link; link=next) {
381: next = link->next;
382: PetscFree(link);
383: }
384: (*dm)->coarsenhook = NULL;
385: }
386: {
387: DMRefineHookLink link,next;
388: for (link=(*dm)->refinehook; link; link=next) {
389: next = link->next;
390: PetscFree(link);
391: }
392: (*dm)->refinehook = NULL;
393: }
394: {
395: DMSubDomainHookLink link,next;
396: for (link=(*dm)->subdomainhook; link; link=next) {
397: next = link->next;
398: PetscFree(link);
399: }
400: (*dm)->subdomainhook = NULL;
401: }
402: {
403: DMGlobalToLocalHookLink link,next;
404: for (link=(*dm)->gtolhook; link; link=next) {
405: next = link->next;
406: PetscFree(link);
407: }
408: (*dm)->gtolhook = NULL;
409: }
410: /* Destroy the work arrays */
411: {
412: DMWorkLink link,next;
413: if ((*dm)->workout) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Work array still checked out");
414: for (link=(*dm)->workin; link; link=next) {
415: next = link->next;
416: PetscFree(link->mem);
417: PetscFree(link);
418: }
419: (*dm)->workin = NULL;
420: }
422: PetscObjectDestroy(&(*dm)->dmksp);
423: PetscObjectDestroy(&(*dm)->dmsnes);
424: PetscObjectDestroy(&(*dm)->dmts);
426: if ((*dm)->ctx && (*dm)->ctxdestroy) {
427: (*(*dm)->ctxdestroy)(&(*dm)->ctx);
428: }
429: VecDestroy(&(*dm)->x);
430: MatFDColoringDestroy(&(*dm)->fd);
431: DMClearGlobalVectors(*dm);
432: ISLocalToGlobalMappingDestroy(&(*dm)->ltogmap);
433: ISLocalToGlobalMappingDestroy(&(*dm)->ltogmapb);
434: PetscFree((*dm)->vectype);
435: PetscFree((*dm)->mattype);
437: PetscSectionDestroy(&(*dm)->defaultSection);
438: PetscSectionDestroy(&(*dm)->defaultGlobalSection);
439: PetscLayoutDestroy(&(*dm)->map);
440: PetscSFDestroy(&(*dm)->sf);
441: PetscSFDestroy(&(*dm)->defaultSF);
443: DMDestroy(&(*dm)->coordinateDM);
444: VecDestroy(&(*dm)->coordinates);
445: VecDestroy(&(*dm)->coordinatesLocal);
447: for (f = 0; f < (*dm)->numFields; ++f) {
448: PetscObjectDestroy(&(*dm)->fields[f]);
449: }
450: PetscFree((*dm)->fields);
451: /* if memory was published with AMS then destroy it */
452: PetscObjectAMSViewOff((PetscObject)*dm);
454: (*(*dm)->ops->destroy)(*dm);
455: /* We do not destroy (*dm)->data here so that we can reference count backend objects */
456: PetscHeaderDestroy(dm);
457: return(0);
458: }
462: /*@
463: DMSetUp - sets up the data structures inside a DM object
465: Collective on DM
467: Input Parameter:
468: . dm - the DM object to setup
470: Level: developer
472: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix()
474: @*/
475: PetscErrorCode DMSetUp(DM dm)
476: {
481: if (dm->setupcalled) return(0);
482: if (dm->ops->setup) {
483: (*dm->ops->setup)(dm);
484: }
485: dm->setupcalled = PETSC_TRUE;
486: return(0);
487: }
491: /*@
492: DMSetFromOptions - sets parameters in a DM from the options database
494: Collective on DM
496: Input Parameter:
497: . dm - the DM object to set options for
499: Options Database:
500: + -dm_preallocate_only: Only preallocate the matrix for DMCreateMatrix(), but do not fill it with zeros
501: . -dm_vec_type <type> type of vector to create inside DM
502: . -dm_mat_type <type> type of matrix to create inside DM
503: - -dm_coloring_type <global or ghosted>
505: Level: developer
507: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix()
509: @*/
510: PetscErrorCode DMSetFromOptions(DM dm)
511: {
512: char typeName[256] = MATAIJ;
513: PetscBool flg;
518: PetscObjectOptionsBegin((PetscObject)dm);
519: PetscOptionsBool("-dm_preallocate_only","only preallocate matrix, but do not set column indices","DMSetMatrixPreallocateOnly",dm->prealloc_only,&dm->prealloc_only,NULL);
520: PetscOptionsList("-dm_vec_type","Vector type used for created vectors","DMSetVecType",VecList,dm->vectype,typeName,256,&flg);
521: if (flg) {
522: DMSetVecType(dm,typeName);
523: }
524: PetscOptionsList("-dm_mat_type","Matrix type used for created matrices","DMSetMatType",MatList,dm->mattype ? dm->mattype : typeName,typeName,sizeof(typeName),&flg);
525: if (flg) {
526: DMSetMatType(dm,typeName);
527: }
528: PetscOptionsEnum("-dm_is_coloring_type","Global or local coloring of Jacobian","ISColoringType",ISColoringTypes,(PetscEnum)dm->coloringtype,(PetscEnum*)&dm->coloringtype,NULL);
529: if (dm->ops->setfromoptions) {
530: (*dm->ops->setfromoptions)(dm);
531: }
532: /* process any options handlers added with PetscObjectAddOptionsHandler() */
533: PetscObjectProcessOptionsHandlers((PetscObject) dm);
534: PetscOptionsEnd();
535: DMViewFromOptions(dm,NULL,"-dm_view");
536: return(0);
537: }
541: /*@C
542: DMView - Views a vector packer or DMDA.
544: Collective on DM
546: Input Parameter:
547: + dm - the DM object to view
548: - v - the viewer
550: Level: developer
552: .seealso DMDestroy(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix()
554: @*/
555: PetscErrorCode DMView(DM dm,PetscViewer v)
556: {
558: PetscBool isbinary;
562: if (!v) {
563: PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)dm),&v);
564: }
565: PetscObjectTypeCompare((PetscObject)v,PETSCVIEWERBINARY,&isbinary);
566: if (isbinary) {
567: PetscInt classid = DM_FILE_CLASSID;
568: char type[256];
570: PetscViewerBinaryWrite(v,&classid,1,PETSC_INT,PETSC_FALSE);
571: PetscStrncpy(type,((PetscObject)dm)->type_name,256);
572: PetscViewerBinaryWrite(v,type,256,PETSC_CHAR,PETSC_FALSE);
573: }
574: if (dm->ops->view) {
575: (*dm->ops->view)(dm,v);
576: }
577: return(0);
578: }
582: /*@
583: DMCreateGlobalVector - Creates a global vector from a DMDA or DMComposite object
585: Collective on DM
587: Input Parameter:
588: . dm - the DM object
590: Output Parameter:
591: . vec - the global vector
593: Level: beginner
595: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix()
597: @*/
598: PetscErrorCode DMCreateGlobalVector(DM dm,Vec *vec)
599: {
604: (*dm->ops->createglobalvector)(dm,vec);
605: return(0);
606: }
610: /*@
611: DMCreateLocalVector - Creates a local vector from a DMDA or DMComposite object
613: Not Collective
615: Input Parameter:
616: . dm - the DM object
618: Output Parameter:
619: . vec - the local vector
621: Level: beginner
623: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix()
625: @*/
626: PetscErrorCode DMCreateLocalVector(DM dm,Vec *vec)
627: {
632: (*dm->ops->createlocalvector)(dm,vec);
633: return(0);
634: }
638: /*@
639: DMGetLocalToGlobalMapping - Accesses the local-to-global mapping in a DM.
641: Collective on DM
643: Input Parameter:
644: . dm - the DM that provides the mapping
646: Output Parameter:
647: . ltog - the mapping
649: Level: intermediate
651: Notes:
652: This mapping can then be used by VecSetLocalToGlobalMapping() or
653: MatSetLocalToGlobalMapping().
655: .seealso: DMCreateLocalVector(), DMGetLocalToGlobalMappingBlock()
656: @*/
657: PetscErrorCode DMGetLocalToGlobalMapping(DM dm,ISLocalToGlobalMapping *ltog)
658: {
664: if (!dm->ltogmap) {
665: PetscSection section, sectionGlobal;
667: DMGetDefaultSection(dm, §ion);
668: if (section) {
669: PetscInt *ltog;
670: PetscInt pStart, pEnd, size, p, l;
672: DMGetDefaultGlobalSection(dm, §ionGlobal);
673: PetscSectionGetChart(section, &pStart, &pEnd);
674: PetscSectionGetStorageSize(section, &size);
675: PetscMalloc(size * sizeof(PetscInt), <og); /* We want the local+overlap size */
676: for (p = pStart, l = 0; p < pEnd; ++p) {
677: PetscInt dof, off, c;
679: /* Should probably use constrained dofs */
680: PetscSectionGetDof(section, p, &dof);
681: PetscSectionGetOffset(sectionGlobal, p, &off);
682: for (c = 0; c < dof; ++c, ++l) {
683: ltog[l] = off+c;
684: }
685: }
686: ISLocalToGlobalMappingCreate(PETSC_COMM_SELF, size, ltog, PETSC_OWN_POINTER, &dm->ltogmap);
687: PetscLogObjectParent(dm, dm->ltogmap);
688: } else {
689: if (!dm->ops->getlocaltoglobalmapping) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM can not create LocalToGlobalMapping");
690: (*dm->ops->getlocaltoglobalmapping)(dm);
691: }
692: }
693: *ltog = dm->ltogmap;
694: return(0);
695: }
699: /*@
700: DMGetLocalToGlobalMappingBlock - Accesses the blocked local-to-global mapping in a DM.
702: Collective on DM
704: Input Parameter:
705: . da - the distributed array that provides the mapping
707: Output Parameter:
708: . ltog - the block mapping
710: Level: intermediate
712: Notes:
713: This mapping can then be used by VecSetLocalToGlobalMappingBlock() or
714: MatSetLocalToGlobalMappingBlock().
716: .seealso: DMCreateLocalVector(), DMGetLocalToGlobalMapping(), DMGetBlockSize(), VecSetBlockSize(), MatSetBlockSize()
717: @*/
718: PetscErrorCode DMGetLocalToGlobalMappingBlock(DM dm,ISLocalToGlobalMapping *ltog)
719: {
725: if (!dm->ltogmapb) {
726: PetscInt bs;
727: DMGetBlockSize(dm,&bs);
728: if (bs > 1) {
729: if (!dm->ops->getlocaltoglobalmappingblock) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM can not create LocalToGlobalMappingBlock");
730: (*dm->ops->getlocaltoglobalmappingblock)(dm);
731: } else {
732: DMGetLocalToGlobalMapping(dm,&dm->ltogmapb);
733: PetscObjectReference((PetscObject)dm->ltogmapb);
734: }
735: }
736: *ltog = dm->ltogmapb;
737: return(0);
738: }
742: /*@
743: DMGetBlockSize - Gets the inherent block size associated with a DM
745: Not Collective
747: Input Parameter:
748: . dm - the DM with block structure
750: Output Parameter:
751: . bs - the block size, 1 implies no exploitable block structure
753: Level: intermediate
755: .seealso: ISCreateBlock(), VecSetBlockSize(), MatSetBlockSize(), DMGetLocalToGlobalMappingBlock()
756: @*/
757: PetscErrorCode DMGetBlockSize(DM dm,PetscInt *bs)
758: {
762: if (dm->bs < 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"DM does not have enough information to provide a block size yet");
763: *bs = dm->bs;
764: return(0);
765: }
769: /*@
770: DMCreateInterpolation - Gets interpolation matrix between two DMDA or DMComposite objects
772: Collective on DM
774: Input Parameter:
775: + dm1 - the DM object
776: - dm2 - the second, finer DM object
778: Output Parameter:
779: + mat - the interpolation
780: - vec - the scaling (optional)
782: Level: developer
784: 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
785: DMCoarsen(). The coordinates set into the DMDA are completely ignored in computing the interpolation.
787: For DMDA objects you can use this interpolation (more precisely the interpolation from the DMGetCoordinateDM()) to interpolate the mesh coordinate vectors
788: EXCEPT in the periodic case where it does not make sense since the coordinate vectors are not periodic.
791: .seealso DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateColoring(), DMCreateMatrix(), DMRefine(), DMCoarsen()
793: @*/
794: PetscErrorCode DMCreateInterpolation(DM dm1,DM dm2,Mat *mat,Vec *vec)
795: {
801: (*dm1->ops->createinterpolation)(dm1,dm2,mat,vec);
802: return(0);
803: }
807: /*@
808: DMCreateInjection - Gets injection matrix between two DMDA or DMComposite objects
810: Collective on DM
812: Input Parameter:
813: + dm1 - the DM object
814: - dm2 - the second, finer DM object
816: Output Parameter:
817: . ctx - the injection
819: Level: developer
821: 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
822: DMCoarsen(). The coordinates set into the DMDA are completely ignored in computing the injection.
824: .seealso DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateColoring(), DMCreateMatrix(), DMCreateInterpolation()
826: @*/
827: PetscErrorCode DMCreateInjection(DM dm1,DM dm2,VecScatter *ctx)
828: {
834: (*dm1->ops->getinjection)(dm1,dm2,ctx);
835: return(0);
836: }
840: /*@C
841: DMCreateColoring - Gets coloring for a DMDA or DMComposite
843: Collective on DM
845: Input Parameter:
846: + dm - the DM object
847: . ctype - IS_COLORING_GHOSTED or IS_COLORING_GLOBAL
848: - matype - either MATAIJ or MATBAIJ
850: Output Parameter:
851: . coloring - the coloring
853: Level: developer
855: .seealso DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateMatrix()
857: @*/
858: PetscErrorCode DMCreateColoring(DM dm,ISColoringType ctype,MatType mtype,ISColoring *coloring)
859: {
864: if (!dm->ops->getcoloring) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"No coloring for this type of DM yet");
865: (*dm->ops->getcoloring)(dm,ctype,mtype,coloring);
866: return(0);
867: }
871: /*@C
872: DMCreateMatrix - Gets empty Jacobian for a DMDA or DMComposite
874: Collective on DM
876: Input Parameter:
877: + dm - the DM object
878: - mtype - Supported types are MATSEQAIJ, MATMPIAIJ, MATSEQBAIJ, MATMPIBAIJ, or
879: any type which inherits from one of these (such as MATAIJ)
881: Output Parameter:
882: . mat - the empty Jacobian
884: Level: beginner
886: Notes: This properly preallocates the number of nonzeros in the sparse matrix so you
887: do not need to do it yourself.
889: By default it also sets the nonzero structure and puts in the zero entries. To prevent setting
890: the nonzero pattern call DMDASetMatPreallocateOnly()
892: For structured grid problems, when you call MatView() on this matrix it is displayed using the global natural ordering, NOT in the ordering used
893: internally by PETSc.
895: For structured grid problems, in general it is easiest to use MatSetValuesStencil() or MatSetValuesLocal() to put values into the matrix because MatSetValues() requires
896: the indices for the global numbering for DMDAs which is complicated.
898: .seealso DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation()
900: @*/
901: PetscErrorCode DMCreateMatrix(DM dm,MatType mtype,Mat *mat)
902: {
907: #if !defined(PETSC_USE_DYNAMIC_LIBRARIES)
908: MatInitializePackage();
909: #endif
912: if (dm->mattype) {
913: (*dm->ops->creatematrix)(dm,dm->mattype,mat);
914: } else {
915: (*dm->ops->creatematrix)(dm,mtype,mat);
916: }
917: return(0);
918: }
922: /*@
923: DMSetMatrixPreallocateOnly - When DMCreateMatrix() is called the matrix will be properly
924: preallocated but the nonzero structure and zero values will not be set.
926: Logically Collective on DMDA
928: Input Parameter:
929: + dm - the DM
930: - only - PETSC_TRUE if only want preallocation
932: Level: developer
933: .seealso DMCreateMatrix()
934: @*/
935: PetscErrorCode DMSetMatrixPreallocateOnly(DM dm, PetscBool only)
936: {
939: dm->prealloc_only = only;
940: return(0);
941: }
945: /*@C
946: DMGetWorkArray - Gets a work array guaranteed to be at least the input size, restore with DMRestoreWorkArray()
948: Not Collective
950: Input Parameters:
951: + dm - the DM object
952: . count - The minium size
953: - dtype - data type (PETSC_REAL, PETSC_SCALAR, PETSC_INT)
955: Output Parameter:
956: . array - the work array
958: Level: developer
960: .seealso DMDestroy(), DMCreate()
961: @*/
962: PetscErrorCode DMGetWorkArray(DM dm,PetscInt count,PetscDataType dtype,void *mem)
963: {
965: DMWorkLink link;
966: size_t size;
971: if (dm->workin) {
972: link = dm->workin;
973: dm->workin = dm->workin->next;
974: } else {
975: PetscNewLog(dm,struct _DMWorkLink,&link);
976: }
977: PetscDataTypeGetSize(dtype,&size);
978: if (size*count > link->bytes) {
979: PetscFree(link->mem);
980: PetscMalloc(size*count,&link->mem);
981: link->bytes = size*count;
982: }
983: link->next = dm->workout;
984: dm->workout = link;
985: *(void**)mem = link->mem;
986: return(0);
987: }
991: /*@C
992: DMRestoreWorkArray - Restores a work array guaranteed to be at least the input size, restore with DMRestoreWorkArray()
994: Not Collective
996: Input Parameters:
997: + dm - the DM object
998: . count - The minium size
999: - dtype - data type (PETSC_REAL, PETSC_SCALAR, PETSC_INT)
1001: Output Parameter:
1002: . array - the work array
1004: Level: developer
1006: .seealso DMDestroy(), DMCreate()
1007: @*/
1008: PetscErrorCode DMRestoreWorkArray(DM dm,PetscInt count,PetscDataType dtype,void *mem)
1009: {
1010: DMWorkLink *p,link;
1015: for (p=&dm->workout; (link=*p); p=&link->next) {
1016: if (link->mem == *(void**)mem) {
1017: *p = link->next;
1018: link->next = dm->workin;
1019: dm->workin = link;
1020: *(void**)mem = NULL;
1021: return(0);
1022: }
1023: }
1024: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Array was not checked out");
1025: return(0);
1026: }
1030: PetscErrorCode DMSetNullSpaceConstructor(DM dm, PetscInt field, PetscErrorCode (*nullsp)(DM dm, PetscInt field, MatNullSpace *nullSpace))
1031: {
1034: if (field >= 10) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle %d >= 10 fields", field);
1035: dm->nullspaceConstructors[field] = nullsp;
1036: return(0);
1037: }
1041: /*@C
1042: DMCreateFieldIS - Creates a set of IS objects with the global indices of dofs for each field
1044: Not collective
1046: Input Parameter:
1047: . dm - the DM object
1049: Output Parameters:
1050: + numFields - The number of fields (or NULL if not requested)
1051: . fieldNames - The name for each field (or NULL if not requested)
1052: - fields - The global indices for each field (or NULL if not requested)
1054: Level: intermediate
1056: Notes:
1057: The user is responsible for freeing all requested arrays. In particular, every entry of names should be freed with
1058: PetscFree(), every entry of fields should be destroyed with ISDestroy(), and both arrays should be freed with
1059: PetscFree().
1061: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix()
1062: @*/
1063: PetscErrorCode DMCreateFieldIS(DM dm, PetscInt *numFields, char ***fieldNames, IS **fields)
1064: {
1065: PetscSection section, sectionGlobal;
1070: if (numFields) {
1072: *numFields = 0;
1073: }
1074: if (fieldNames) {
1076: *fieldNames = NULL;
1077: }
1078: if (fields) {
1080: *fields = NULL;
1081: }
1082: DMGetDefaultSection(dm, §ion);
1083: if (section) {
1084: PetscInt *fieldSizes, **fieldIndices;
1085: PetscInt nF, f, pStart, pEnd, p;
1087: DMGetDefaultGlobalSection(dm, §ionGlobal);
1088: PetscSectionGetNumFields(section, &nF);
1089: PetscMalloc2(nF,PetscInt,&fieldSizes,nF,PetscInt*,&fieldIndices);
1090: PetscSectionGetChart(sectionGlobal, &pStart, &pEnd);
1091: for (f = 0; f < nF; ++f) {
1092: fieldSizes[f] = 0;
1093: }
1094: for (p = pStart; p < pEnd; ++p) {
1095: PetscInt gdof;
1097: PetscSectionGetDof(sectionGlobal, p, &gdof);
1098: if (gdof > 0) {
1099: for (f = 0; f < nF; ++f) {
1100: PetscInt fdof, fcdof;
1102: PetscSectionGetFieldDof(section, p, f, &fdof);
1103: PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);
1104: fieldSizes[f] += fdof-fcdof;
1105: }
1106: }
1107: }
1108: for (f = 0; f < nF; ++f) {
1109: PetscMalloc(fieldSizes[f] * sizeof(PetscInt), &fieldIndices[f]);
1110: fieldSizes[f] = 0;
1111: }
1112: for (p = pStart; p < pEnd; ++p) {
1113: PetscInt gdof, goff;
1115: PetscSectionGetDof(sectionGlobal, p, &gdof);
1116: if (gdof > 0) {
1117: PetscSectionGetOffset(sectionGlobal, p, &goff);
1118: for (f = 0; f < nF; ++f) {
1119: PetscInt fdof, fcdof, fc;
1121: PetscSectionGetFieldDof(section, p, f, &fdof);
1122: PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);
1123: for (fc = 0; fc < fdof-fcdof; ++fc, ++fieldSizes[f]) {
1124: fieldIndices[f][fieldSizes[f]] = goff++;
1125: }
1126: }
1127: }
1128: }
1129: if (numFields) *numFields = nF;
1130: if (fieldNames) {
1131: PetscMalloc(nF * sizeof(char*), fieldNames);
1132: for (f = 0; f < nF; ++f) {
1133: const char *fieldName;
1135: PetscSectionGetFieldName(section, f, &fieldName);
1136: PetscStrallocpy(fieldName, (char**) &(*fieldNames)[f]);
1137: }
1138: }
1139: if (fields) {
1140: PetscMalloc(nF * sizeof(IS), fields);
1141: for (f = 0; f < nF; ++f) {
1142: ISCreateGeneral(PetscObjectComm((PetscObject)dm), fieldSizes[f], fieldIndices[f], PETSC_OWN_POINTER, &(*fields)[f]);
1143: }
1144: }
1145: PetscFree2(fieldSizes,fieldIndices);
1146: } else if (dm->ops->createfieldis) {
1147: (*dm->ops->createfieldis)(dm, numFields, fieldNames, fields);
1148: }
1149: return(0);
1150: }
1155: /*@C
1156: DMCreateFieldDecomposition - Returns a list of IS objects defining a decomposition of a problem into subproblems
1157: corresponding to different fields: each IS contains the global indices of the dofs of the
1158: corresponding field. The optional list of DMs define the DM for each subproblem.
1159: Generalizes DMCreateFieldIS().
1161: Not collective
1163: Input Parameter:
1164: . dm - the DM object
1166: Output Parameters:
1167: + len - The number of subproblems in the field decomposition (or NULL if not requested)
1168: . namelist - The name for each field (or NULL if not requested)
1169: . islist - The global indices for each field (or NULL if not requested)
1170: - dmlist - The DMs for each field subproblem (or NULL, if not requested; if NULL is returned, no DMs are defined)
1172: Level: intermediate
1174: Notes:
1175: The user is responsible for freeing all requested arrays. In particular, every entry of names should be freed with
1176: PetscFree(), every entry of is should be destroyed with ISDestroy(), every entry of dm should be destroyed with DMDestroy(),
1177: and all of the arrays should be freed with PetscFree().
1179: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateFieldIS()
1180: @*/
1181: PetscErrorCode DMCreateFieldDecomposition(DM dm, PetscInt *len, char ***namelist, IS **islist, DM **dmlist)
1182: {
1187: if (len) {
1189: *len = 0;
1190: }
1191: if (namelist) {
1193: *namelist = 0;
1194: }
1195: if (islist) {
1197: *islist = 0;
1198: }
1199: if (dmlist) {
1201: *dmlist = 0;
1202: }
1203: /*
1204: Is it a good idea to apply the following check across all impls?
1205: Perhaps some impls can have a well-defined decomposition before DMSetUp?
1206: This, however, follows the general principle that accessors are not well-behaved until the object is set up.
1207: */
1208: if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE, "Decomposition defined only after DMSetUp");
1209: if (!dm->ops->createfielddecomposition) {
1210: PetscSection section;
1211: PetscInt numFields, f;
1213: DMGetDefaultSection(dm, §ion);
1214: if (section) {PetscSectionGetNumFields(section, &numFields);}
1215: if (section && numFields && dm->ops->createsubdm) {
1216: *len = numFields;
1217: PetscMalloc3(numFields,char*,namelist,numFields,IS,islist,numFields,DM,dmlist);
1218: for (f = 0; f < numFields; ++f) {
1219: const char *fieldName;
1221: DMCreateSubDM(dm, 1, &f, &(*islist)[f], &(*dmlist)[f]);
1222: PetscSectionGetFieldName(section, f, &fieldName);
1223: PetscStrallocpy(fieldName, (char**) &(*namelist)[f]);
1224: }
1225: } else {
1226: DMCreateFieldIS(dm, len, namelist, islist);
1227: /* By default there are no DMs associated with subproblems. */
1228: if (dmlist) *dmlist = NULL;
1229: }
1230: } else {
1231: (*dm->ops->createfielddecomposition)(dm,len,namelist,islist,dmlist);
1232: }
1233: return(0);
1234: }
1238: /*@C
1239: DMCreateSubDM - Returns an IS and DM encapsulating a subproblem defined by the fields passed in.
1240: The fields are defined by DMCreateFieldIS().
1242: Not collective
1244: Input Parameters:
1245: + dm - the DM object
1246: . numFields - number of fields in this subproblem
1247: - len - The number of subproblems in the decomposition (or NULL if not requested)
1249: Output Parameters:
1250: . is - The global indices for the subproblem
1251: - dm - The DM for the subproblem
1253: Level: intermediate
1255: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateFieldIS()
1256: @*/
1257: PetscErrorCode DMCreateSubDM(DM dm, PetscInt numFields, PetscInt fields[], IS *is, DM *subdm)
1258: {
1266: if (dm->ops->createsubdm) {
1267: (*dm->ops->createsubdm)(dm, numFields, fields, is, subdm);
1268: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "This type has no DMCreateSubDM implementation defined");
1269: return(0);
1270: }
1275: /*@C
1276: DMCreateDomainDecomposition - Returns lists of IS objects defining a decomposition of a problem into subproblems
1277: corresponding to restrictions to pairs nested subdomains: each IS contains the global
1278: indices of the dofs of the corresponding subdomains. The inner subdomains conceptually
1279: define a nonoverlapping covering, while outer subdomains can overlap.
1280: The optional list of DMs define the DM for each subproblem.
1282: Not collective
1284: Input Parameter:
1285: . dm - the DM object
1287: Output Parameters:
1288: + len - The number of subproblems in the domain decomposition (or NULL if not requested)
1289: . namelist - The name for each subdomain (or NULL if not requested)
1290: . innerislist - The global indices for each inner subdomain (or NULL, if not requested)
1291: . outerislist - The global indices for each outer subdomain (or NULL, if not requested)
1292: - dmlist - The DMs for each subdomain subproblem (or NULL, if not requested; if NULL is returned, no DMs are defined)
1294: Level: intermediate
1296: Notes:
1297: The user is responsible for freeing all requested arrays. In particular, every entry of names should be freed with
1298: PetscFree(), every entry of is should be destroyed with ISDestroy(), every entry of dm should be destroyed with DMDestroy(),
1299: and all of the arrays should be freed with PetscFree().
1301: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateDomainDecompositionDM(), DMCreateFieldDecomposition()
1302: @*/
1303: PetscErrorCode DMCreateDomainDecomposition(DM dm, PetscInt *len, char ***namelist, IS **innerislist, IS **outerislist, DM **dmlist)
1304: {
1305: PetscErrorCode ierr;
1306: DMSubDomainHookLink link;
1307: PetscInt i,l;
1316: /*
1317: Is it a good idea to apply the following check across all impls?
1318: Perhaps some impls can have a well-defined decomposition before DMSetUp?
1319: This, however, follows the general principle that accessors are not well-behaved until the object is set up.
1320: */
1321: if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE, "Decomposition defined only after DMSetUp");
1322: if (dm->ops->createdomaindecomposition) {
1323: (*dm->ops->createdomaindecomposition)(dm,&l,namelist,innerislist,outerislist,dmlist);
1324: /* copy subdomain hooks and context over to the subdomain DMs */
1325: if (dmlist) {
1326: for (i = 0; i < l; i++) {
1327: for (link=dm->subdomainhook; link; link=link->next) {
1328: if (link->ddhook) {(*link->ddhook)(dm,(*dmlist)[i],link->ctx);}
1329: }
1330: (*dmlist)[i]->ctx = dm->ctx;
1331: }
1332: }
1333: if (len) *len = l;
1334: }
1335: return(0);
1336: }
1341: /*@C
1342: DMCreateDomainDecompositionScatters - Returns scatters to the subdomain vectors from the global vector
1344: Not collective
1346: Input Parameters:
1347: + dm - the DM object
1348: . n - the number of subdomain scatters
1349: - subdms - the local subdomains
1351: Output Parameters:
1352: + n - the number of scatters returned
1353: . iscat - scatter from global vector to nonoverlapping global vector entries on subdomain
1354: . oscat - scatter from global vector to overlapping global vector entries on subdomain
1355: - gscat - scatter from global vector to local vector on subdomain (fills in ghosts)
1357: Notes: This is an alternative to the iis and ois arguments in DMCreateDomainDecomposition that allow for the solution
1358: of general nonlinear problems with overlapping subdomain methods. While merely having index sets that enable subsets
1359: of the residual equations to be created is fine for linear problems, nonlinear problems require local assembly of
1360: solution and residual data.
1362: Level: developer
1364: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateFieldIS()
1365: @*/
1366: PetscErrorCode DMCreateDomainDecompositionScatters(DM dm,PetscInt n,DM *subdms,VecScatter **iscat,VecScatter **oscat,VecScatter **gscat)
1367: {
1373: if (dm->ops->createddscatters) {
1374: (*dm->ops->createddscatters)(dm,n,subdms,iscat,oscat,gscat);
1375: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "This type has no DMCreateDomainDecompositionLocalScatter implementation defined");
1376: return(0);
1377: }
1381: /*@
1382: DMRefine - Refines a DM object
1384: Collective on DM
1386: Input Parameter:
1387: + dm - the DM object
1388: - comm - the communicator to contain the new DM object (or MPI_COMM_NULL)
1390: Output Parameter:
1391: . dmf - the refined DM, or NULL
1393: Note: If no refinement was done, the return value is NULL
1395: Level: developer
1397: .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation()
1398: @*/
1399: PetscErrorCode DMRefine(DM dm,MPI_Comm comm,DM *dmf)
1400: {
1401: PetscErrorCode ierr;
1402: DMRefineHookLink link;
1406: (*dm->ops->refine)(dm,comm,dmf);
1407: if (*dmf) {
1408: (*dmf)->ops->creatematrix = dm->ops->creatematrix;
1410: PetscObjectCopyFortranFunctionPointers((PetscObject)dm,(PetscObject)*dmf);
1412: (*dmf)->ctx = dm->ctx;
1413: (*dmf)->leveldown = dm->leveldown;
1414: (*dmf)->levelup = dm->levelup + 1;
1416: DMSetMatType(*dmf,dm->mattype);
1417: for (link=dm->refinehook; link; link=link->next) {
1418: if (link->refinehook) {
1419: (*link->refinehook)(dm,*dmf,link->ctx);
1420: }
1421: }
1422: }
1423: return(0);
1424: }
1428: /*@C
1429: DMRefineHookAdd - adds a callback to be run when interpolating a nonlinear problem to a finer grid
1431: Logically Collective
1433: Input Arguments:
1434: + coarse - nonlinear solver context on which to run a hook when restricting to a coarser level
1435: . refinehook - function to run when setting up a coarser level
1436: . interphook - function to run to update data on finer levels (once per SNESSolve())
1437: - ctx - [optional] user-defined context for provide data for the hooks (may be NULL)
1439: Calling sequence of refinehook:
1440: $ refinehook(DM coarse,DM fine,void *ctx);
1442: + coarse - coarse level DM
1443: . fine - fine level DM to interpolate problem to
1444: - ctx - optional user-defined function context
1446: Calling sequence for interphook:
1447: $ interphook(DM coarse,Mat interp,DM fine,void *ctx)
1449: + coarse - coarse level DM
1450: . interp - matrix interpolating a coarse-level solution to the finer grid
1451: . fine - fine level DM to update
1452: - ctx - optional user-defined function context
1454: Level: advanced
1456: Notes:
1457: This function is only needed if auxiliary data needs to be passed to fine grids while grid sequencing
1459: If this function is called multiple times, the hooks will be run in the order they are added.
1461: This function is currently not available from Fortran.
1463: .seealso: DMCoarsenHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
1464: @*/
1465: PetscErrorCode DMRefineHookAdd(DM coarse,PetscErrorCode (*refinehook)(DM,DM,void*),PetscErrorCode (*interphook)(DM,Mat,DM,void*),void *ctx)
1466: {
1467: PetscErrorCode ierr;
1468: DMRefineHookLink link,*p;
1472: for (p=&coarse->refinehook; *p; p=&(*p)->next) {} /* Scan to the end of the current list of hooks */
1473: PetscMalloc(sizeof(struct _DMRefineHookLink),&link);
1474: link->refinehook = refinehook;
1475: link->interphook = interphook;
1476: link->ctx = ctx;
1477: link->next = NULL;
1478: *p = link;
1479: return(0);
1480: }
1484: /*@
1485: DMInterpolate - interpolates user-defined problem data to a finer DM by running hooks registered by DMRefineHookAdd()
1487: Collective if any hooks are
1489: Input Arguments:
1490: + coarse - coarser DM to use as a base
1491: . restrct - interpolation matrix, apply using MatInterpolate()
1492: - fine - finer DM to update
1494: Level: developer
1496: .seealso: DMRefineHookAdd(), MatInterpolate()
1497: @*/
1498: PetscErrorCode DMInterpolate(DM coarse,Mat interp,DM fine)
1499: {
1500: PetscErrorCode ierr;
1501: DMRefineHookLink link;
1504: for (link=fine->refinehook; link; link=link->next) {
1505: if (link->interphook) {
1506: (*link->interphook)(coarse,interp,fine,link->ctx);
1507: }
1508: }
1509: return(0);
1510: }
1514: /*@
1515: DMGetRefineLevel - Get's the number of refinements that have generated this DM.
1517: Not Collective
1519: Input Parameter:
1520: . dm - the DM object
1522: Output Parameter:
1523: . level - number of refinements
1525: Level: developer
1527: .seealso DMCoarsen(), DMGetCoarsenLevel(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation()
1529: @*/
1530: PetscErrorCode DMGetRefineLevel(DM dm,PetscInt *level)
1531: {
1534: *level = dm->levelup;
1535: return(0);
1536: }
1540: /*@C
1541: DMGlobalToLocalHookAdd - adds a callback to be run when global to local is called
1543: Logically Collective
1545: Input Arguments:
1546: + dm - the DM
1547: . beginhook - function to run at the beginning of DMGlobalToLocalBegin()
1548: . endhook - function to run after DMGlobalToLocalEnd() has completed
1549: - ctx - [optional] user-defined context for provide data for the hooks (may be NULL)
1551: Calling sequence for beginhook:
1552: $ beginhook(DM fine,VecScatter out,VecScatter in,DM coarse,void *ctx)
1554: + dm - global DM
1555: . g - global vector
1556: . mode - mode
1557: . l - local vector
1558: - ctx - optional user-defined function context
1561: Calling sequence for endhook:
1562: $ endhook(DM fine,VecScatter out,VecScatter in,DM coarse,void *ctx)
1564: + global - global DM
1565: - ctx - optional user-defined function context
1567: Level: advanced
1569: .seealso: DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
1570: @*/
1571: PetscErrorCode DMGlobalToLocalHookAdd(DM dm,PetscErrorCode (*beginhook)(DM,Vec,InsertMode,Vec,void*),PetscErrorCode (*endhook)(DM,Vec,InsertMode,Vec,void*),void *ctx)
1572: {
1573: PetscErrorCode ierr;
1574: DMGlobalToLocalHookLink link,*p;
1578: for (p=&dm->gtolhook; *p; p=&(*p)->next) {} /* Scan to the end of the current list of hooks */
1579: PetscMalloc(sizeof(struct _DMGlobalToLocalHookLink),&link);
1580: link->beginhook = beginhook;
1581: link->endhook = endhook;
1582: link->ctx = ctx;
1583: link->next = NULL;
1584: *p = link;
1585: return(0);
1586: }
1590: /*@
1591: DMGlobalToLocalBegin - Begins updating local vectors from global vector
1593: Neighbor-wise Collective on DM
1595: Input Parameters:
1596: + dm - the DM object
1597: . g - the global vector
1598: . mode - INSERT_VALUES or ADD_VALUES
1599: - l - the local vector
1602: Level: beginner
1604: .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMGlobalToLocalEnd(), DMLocalToGlobalBegin()
1606: @*/
1607: PetscErrorCode DMGlobalToLocalBegin(DM dm,Vec g,InsertMode mode,Vec l)
1608: {
1609: PetscSF sf;
1610: PetscErrorCode ierr;
1611: DMGlobalToLocalHookLink link;
1615: for (link=dm->gtolhook; link; link=link->next) {
1616: if (link->beginhook) {
1617: (*link->beginhook)(dm,g,mode,l,link->ctx);
1618: }
1619: }
1620: DMGetDefaultSF(dm, &sf);
1621: if (sf) {
1622: PetscScalar *lArray, *gArray;
1624: if (mode == ADD_VALUES) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode);
1625: VecGetArray(l, &lArray);
1626: VecGetArray(g, &gArray);
1627: PetscSFBcastBegin(sf, MPIU_SCALAR, gArray, lArray);
1628: VecRestoreArray(l, &lArray);
1629: VecRestoreArray(g, &gArray);
1630: } else {
1631: (*dm->ops->globaltolocalbegin)(dm,g,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),l);
1632: }
1633: return(0);
1634: }
1638: /*@
1639: DMGlobalToLocalEnd - Ends updating local vectors from global vector
1641: Neighbor-wise Collective on DM
1643: Input Parameters:
1644: + dm - the DM object
1645: . g - the global vector
1646: . mode - INSERT_VALUES or ADD_VALUES
1647: - l - the local vector
1650: Level: beginner
1652: .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMGlobalToLocalEnd(), DMLocalToGlobalBegin()
1654: @*/
1655: PetscErrorCode DMGlobalToLocalEnd(DM dm,Vec g,InsertMode mode,Vec l)
1656: {
1657: PetscSF sf;
1658: PetscErrorCode ierr;
1659: PetscScalar *lArray, *gArray;
1660: DMGlobalToLocalHookLink link;
1664: DMGetDefaultSF(dm, &sf);
1665: if (sf) {
1666: if (mode == ADD_VALUES) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode);
1668: VecGetArray(l, &lArray);
1669: VecGetArray(g, &gArray);
1670: PetscSFBcastEnd(sf, MPIU_SCALAR, gArray, lArray);
1671: VecRestoreArray(l, &lArray);
1672: VecRestoreArray(g, &gArray);
1673: } else {
1674: (*dm->ops->globaltolocalend)(dm,g,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),l);
1675: }
1676: for (link=dm->gtolhook; link; link=link->next) {
1677: if (link->endhook) {(*link->endhook)(dm,g,mode,l,link->ctx);}
1678: }
1679: return(0);
1680: }
1684: /*@
1685: DMLocalToGlobalBegin - updates global vectors from local vectors
1687: Neighbor-wise Collective on DM
1689: Input Parameters:
1690: + dm - the DM object
1691: . l - the local vector
1692: . 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
1693: base point.
1694: - - the global vector
1696: 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
1697: 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
1698: global array to the final global array with VecAXPY().
1700: Level: beginner
1702: .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMGlobalToLocalEnd(), DMGlobalToLocalBegin()
1704: @*/
1705: PetscErrorCode DMLocalToGlobalBegin(DM dm,Vec l,InsertMode mode,Vec g)
1706: {
1707: PetscSF sf;
1712: DMGetDefaultSF(dm, &sf);
1713: if (sf) {
1714: MPI_Op op;
1715: PetscScalar *lArray, *gArray;
1717: switch (mode) {
1718: case INSERT_VALUES:
1719: case INSERT_ALL_VALUES:
1720: op = MPIU_REPLACE; break;
1721: case ADD_VALUES:
1722: case ADD_ALL_VALUES:
1723: op = MPI_SUM; break;
1724: default:
1725: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode);
1726: }
1727: VecGetArray(l, &lArray);
1728: VecGetArray(g, &gArray);
1729: PetscSFReduceBegin(sf, MPIU_SCALAR, lArray, gArray, op);
1730: VecRestoreArray(l, &lArray);
1731: VecRestoreArray(g, &gArray);
1732: } else {
1733: (*dm->ops->localtoglobalbegin)(dm,l,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),g);
1734: }
1735: return(0);
1736: }
1740: /*@
1741: DMLocalToGlobalEnd - updates global vectors from local vectors
1743: Neighbor-wise Collective on DM
1745: Input Parameters:
1746: + dm - the DM object
1747: . l - the local vector
1748: . mode - INSERT_VALUES or ADD_VALUES
1749: - g - the global vector
1752: Level: beginner
1754: .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMGlobalToLocalEnd(), DMGlobalToLocalEnd()
1756: @*/
1757: PetscErrorCode DMLocalToGlobalEnd(DM dm,Vec l,InsertMode mode,Vec g)
1758: {
1759: PetscSF sf;
1764: DMGetDefaultSF(dm, &sf);
1765: if (sf) {
1766: MPI_Op op;
1767: PetscScalar *lArray, *gArray;
1769: switch (mode) {
1770: case INSERT_VALUES:
1771: case INSERT_ALL_VALUES:
1772: op = MPIU_REPLACE; break;
1773: case ADD_VALUES:
1774: case ADD_ALL_VALUES:
1775: op = MPI_SUM; break;
1776: default:
1777: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode);
1778: }
1779: VecGetArray(l, &lArray);
1780: VecGetArray(g, &gArray);
1781: PetscSFReduceEnd(sf, MPIU_SCALAR, lArray, gArray, op);
1782: VecRestoreArray(l, &lArray);
1783: VecRestoreArray(g, &gArray);
1784: } else {
1785: (*dm->ops->localtoglobalend)(dm,l,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),g);
1786: }
1787: return(0);
1788: }
1792: /*@
1793: DMCoarsen - Coarsens a DM object
1795: Collective on DM
1797: Input Parameter:
1798: + dm - the DM object
1799: - comm - the communicator to contain the new DM object (or MPI_COMM_NULL)
1801: Output Parameter:
1802: . dmc - the coarsened DM
1804: Level: developer
1806: .seealso DMRefine(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation()
1808: @*/
1809: PetscErrorCode DMCoarsen(DM dm, MPI_Comm comm, DM *dmc)
1810: {
1811: PetscErrorCode ierr;
1812: DMCoarsenHookLink link;
1816: (*dm->ops->coarsen)(dm, comm, dmc);
1817: (*dmc)->ops->creatematrix = dm->ops->creatematrix;
1818: PetscObjectCopyFortranFunctionPointers((PetscObject)dm,(PetscObject)*dmc);
1819: (*dmc)->ctx = dm->ctx;
1820: (*dmc)->levelup = dm->levelup;
1821: (*dmc)->leveldown = dm->leveldown + 1;
1822: DMSetMatType(*dmc,dm->mattype);
1823: for (link=dm->coarsenhook; link; link=link->next) {
1824: if (link->coarsenhook) {(*link->coarsenhook)(dm,*dmc,link->ctx);}
1825: }
1826: return(0);
1827: }
1831: /*@C
1832: DMCoarsenHookAdd - adds a callback to be run when restricting a nonlinear problem to the coarse grid
1834: Logically Collective
1836: Input Arguments:
1837: + fine - nonlinear solver context on which to run a hook when restricting to a coarser level
1838: . coarsenhook - function to run when setting up a coarser level
1839: . restricthook - function to run to update data on coarser levels (once per SNESSolve())
1840: - ctx - [optional] user-defined context for provide data for the hooks (may be NULL)
1842: Calling sequence of coarsenhook:
1843: $ coarsenhook(DM fine,DM coarse,void *ctx);
1845: + fine - fine level DM
1846: . coarse - coarse level DM to restrict problem to
1847: - ctx - optional user-defined function context
1849: Calling sequence for restricthook:
1850: $ restricthook(DM fine,Mat mrestrict,Vec rscale,Mat inject,DM coarse,void *ctx)
1852: + fine - fine level DM
1853: . mrestrict - matrix restricting a fine-level solution to the coarse grid
1854: . rscale - scaling vector for restriction
1855: . inject - matrix restricting by injection
1856: . coarse - coarse level DM to update
1857: - ctx - optional user-defined function context
1859: Level: advanced
1861: Notes:
1862: This function is only needed if auxiliary data needs to be set up on coarse grids.
1864: If this function is called multiple times, the hooks will be run in the order they are added.
1866: In order to compose with nonlinear preconditioning without duplicating storage, the hook should be implemented to
1867: extract the finest level information from its context (instead of from the SNES).
1869: This function is currently not available from Fortran.
1871: .seealso: DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
1872: @*/
1873: PetscErrorCode DMCoarsenHookAdd(DM fine,PetscErrorCode (*coarsenhook)(DM,DM,void*),PetscErrorCode (*restricthook)(DM,Mat,Vec,Mat,DM,void*),void *ctx)
1874: {
1875: PetscErrorCode ierr;
1876: DMCoarsenHookLink link,*p;
1880: for (p=&fine->coarsenhook; *p; p=&(*p)->next) {} /* Scan to the end of the current list of hooks */
1881: PetscMalloc(sizeof(struct _DMCoarsenHookLink),&link);
1882: link->coarsenhook = coarsenhook;
1883: link->restricthook = restricthook;
1884: link->ctx = ctx;
1885: link->next = NULL;
1886: *p = link;
1887: return(0);
1888: }
1892: /*@
1893: DMRestrict - restricts user-defined problem data to a coarser DM by running hooks registered by DMCoarsenHookAdd()
1895: Collective if any hooks are
1897: Input Arguments:
1898: + fine - finer DM to use as a base
1899: . restrct - restriction matrix, apply using MatRestrict()
1900: . inject - injection matrix, also use MatRestrict()
1901: - coarse - coarer DM to update
1903: Level: developer
1905: .seealso: DMCoarsenHookAdd(), MatRestrict()
1906: @*/
1907: PetscErrorCode DMRestrict(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse)
1908: {
1909: PetscErrorCode ierr;
1910: DMCoarsenHookLink link;
1913: for (link=fine->coarsenhook; link; link=link->next) {
1914: if (link->restricthook) {
1915: (*link->restricthook)(fine,restrct,rscale,inject,coarse,link->ctx);
1916: }
1917: }
1918: return(0);
1919: }
1923: /*@C
1924: DMSubDomainHookAdd - adds a callback to be run when restricting a problem to the coarse grid
1926: Logically Collective
1928: Input Arguments:
1929: + global - global DM
1930: . ddhook - function to run to pass data to the decomposition DM upon its creation
1931: . restricthook - function to run to update data on block solve (at the beginning of the block solve)
1932: - ctx - [optional] user-defined context for provide data for the hooks (may be NULL)
1935: Calling sequence for ddhook:
1936: $ ddhook(DM global,DM block,void *ctx)
1938: + global - global DM
1939: . block - block DM
1940: - ctx - optional user-defined function context
1942: Calling sequence for restricthook:
1943: $ restricthook(DM global,VecScatter out,VecScatter in,DM block,void *ctx)
1945: + global - global DM
1946: . out - scatter to the outer (with ghost and overlap points) block vector
1947: . in - scatter to block vector values only owned locally
1948: . block - block DM
1949: - ctx - optional user-defined function context
1951: Level: advanced
1953: Notes:
1954: This function is only needed if auxiliary data needs to be set up on subdomain DMs.
1956: If this function is called multiple times, the hooks will be run in the order they are added.
1958: In order to compose with nonlinear preconditioning without duplicating storage, the hook should be implemented to
1959: extract the global information from its context (instead of from the SNES).
1961: This function is currently not available from Fortran.
1963: .seealso: DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
1964: @*/
1965: PetscErrorCode DMSubDomainHookAdd(DM global,PetscErrorCode (*ddhook)(DM,DM,void*),PetscErrorCode (*restricthook)(DM,VecScatter,VecScatter,DM,void*),void *ctx)
1966: {
1967: PetscErrorCode ierr;
1968: DMSubDomainHookLink link,*p;
1972: for (p=&global->subdomainhook; *p; p=&(*p)->next) {} /* Scan to the end of the current list of hooks */
1973: PetscMalloc(sizeof(struct _DMSubDomainHookLink),&link);
1974: link->restricthook = restricthook;
1975: link->ddhook = ddhook;
1976: link->ctx = ctx;
1977: link->next = NULL;
1978: *p = link;
1979: return(0);
1980: }
1984: /*@
1985: DMSubDomainRestrict - restricts user-defined problem data to a block DM by running hooks registered by DMSubDomainHookAdd()
1987: Collective if any hooks are
1989: Input Arguments:
1990: + fine - finer DM to use as a base
1991: . oscatter - scatter from domain global vector filling subdomain global vector with overlap
1992: . gscatter - scatter from domain global vector filling subdomain local vector with ghosts
1993: - coarse - coarer DM to update
1995: Level: developer
1997: .seealso: DMCoarsenHookAdd(), MatRestrict()
1998: @*/
1999: PetscErrorCode DMSubDomainRestrict(DM global,VecScatter oscatter,VecScatter gscatter,DM subdm)
2000: {
2001: PetscErrorCode ierr;
2002: DMSubDomainHookLink link;
2005: for (link=global->subdomainhook; link; link=link->next) {
2006: if (link->restricthook) {
2007: (*link->restricthook)(global,oscatter,gscatter,subdm,link->ctx);
2008: }
2009: }
2010: return(0);
2011: }
2015: /*@
2016: DMGetCoarsenLevel - Get's the number of coarsenings that have generated this DM.
2018: Not Collective
2020: Input Parameter:
2021: . dm - the DM object
2023: Output Parameter:
2024: . level - number of coarsenings
2026: Level: developer
2028: .seealso DMCoarsen(), DMGetRefineLevel(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation()
2030: @*/
2031: PetscErrorCode DMGetCoarsenLevel(DM dm,PetscInt *level)
2032: {
2035: *level = dm->leveldown;
2036: return(0);
2037: }
2043: /*@C
2044: DMRefineHierarchy - Refines a DM object, all levels at once
2046: Collective on DM
2048: Input Parameter:
2049: + dm - the DM object
2050: - nlevels - the number of levels of refinement
2052: Output Parameter:
2053: . dmf - the refined DM hierarchy
2055: Level: developer
2057: .seealso DMCoarsenHierarchy(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation()
2059: @*/
2060: PetscErrorCode DMRefineHierarchy(DM dm,PetscInt nlevels,DM dmf[])
2061: {
2066: if (nlevels < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"nlevels cannot be negative");
2067: if (nlevels == 0) return(0);
2068: if (dm->ops->refinehierarchy) {
2069: (*dm->ops->refinehierarchy)(dm,nlevels,dmf);
2070: } else if (dm->ops->refine) {
2071: PetscInt i;
2073: DMRefine(dm,PetscObjectComm((PetscObject)dm),&dmf[0]);
2074: for (i=1; i<nlevels; i++) {
2075: DMRefine(dmf[i-1],PetscObjectComm((PetscObject)dm),&dmf[i]);
2076: }
2077: } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"No RefineHierarchy for this DM yet");
2078: return(0);
2079: }
2083: /*@C
2084: DMCoarsenHierarchy - Coarsens a DM object, all levels at once
2086: Collective on DM
2088: Input Parameter:
2089: + dm - the DM object
2090: - nlevels - the number of levels of coarsening
2092: Output Parameter:
2093: . dmc - the coarsened DM hierarchy
2095: Level: developer
2097: .seealso DMRefineHierarchy(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation()
2099: @*/
2100: PetscErrorCode DMCoarsenHierarchy(DM dm, PetscInt nlevels, DM dmc[])
2101: {
2106: if (nlevels < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"nlevels cannot be negative");
2107: if (nlevels == 0) return(0);
2109: if (dm->ops->coarsenhierarchy) {
2110: (*dm->ops->coarsenhierarchy)(dm, nlevels, dmc);
2111: } else if (dm->ops->coarsen) {
2112: PetscInt i;
2114: DMCoarsen(dm,PetscObjectComm((PetscObject)dm),&dmc[0]);
2115: for (i=1; i<nlevels; i++) {
2116: DMCoarsen(dmc[i-1],PetscObjectComm((PetscObject)dm),&dmc[i]);
2117: }
2118: } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"No CoarsenHierarchy for this DM yet");
2119: return(0);
2120: }
2124: /*@
2125: DMCreateAggregates - Gets the aggregates that map between
2126: grids associated with two DMs.
2128: Collective on DM
2130: Input Parameters:
2131: + dmc - the coarse grid DM
2132: - dmf - the fine grid DM
2134: Output Parameters:
2135: . rest - the restriction matrix (transpose of the projection matrix)
2137: Level: intermediate
2139: .keywords: interpolation, restriction, multigrid
2141: .seealso: DMRefine(), DMCreateInjection(), DMCreateInterpolation()
2142: @*/
2143: PetscErrorCode DMCreateAggregates(DM dmc, DM dmf, Mat *rest)
2144: {
2150: (*dmc->ops->getaggregates)(dmc, dmf, rest);
2151: return(0);
2152: }
2156: /*@C
2157: DMSetApplicationContextDestroy - Sets a user function that will be called to destroy the application context when the DM is destroyed
2159: Not Collective
2161: Input Parameters:
2162: + dm - the DM object
2163: - destroy - the destroy function
2165: Level: intermediate
2167: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext()
2169: @*/
2170: PetscErrorCode DMSetApplicationContextDestroy(DM dm,PetscErrorCode (*destroy)(void**))
2171: {
2174: dm->ctxdestroy = destroy;
2175: return(0);
2176: }
2180: /*@
2181: DMSetApplicationContext - Set a user context into a DM object
2183: Not Collective
2185: Input Parameters:
2186: + dm - the DM object
2187: - ctx - the user context
2189: Level: intermediate
2191: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext()
2193: @*/
2194: PetscErrorCode DMSetApplicationContext(DM dm,void *ctx)
2195: {
2198: dm->ctx = ctx;
2199: return(0);
2200: }
2204: /*@
2205: DMGetApplicationContext - Gets a user context from a DM object
2207: Not Collective
2209: Input Parameter:
2210: . dm - the DM object
2212: Output Parameter:
2213: . ctx - the user context
2215: Level: intermediate
2217: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext()
2219: @*/
2220: PetscErrorCode DMGetApplicationContext(DM dm,void *ctx)
2221: {
2224: *(void**)ctx = dm->ctx;
2225: return(0);
2226: }
2230: /*@C
2231: DMSetVariableBounds - sets a function to compute the the lower and upper bound vectors for SNESVI.
2233: Logically Collective on DM
2235: Input Parameter:
2236: + dm - the DM object
2237: - f - the function that computes variable bounds used by SNESVI (use NULL to cancel a previous function that was set)
2239: Level: intermediate
2241: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext(),
2242: DMSetJacobian()
2244: @*/
2245: PetscErrorCode DMSetVariableBounds(DM dm,PetscErrorCode (*f)(DM,Vec,Vec))
2246: {
2248: dm->ops->computevariablebounds = f;
2249: return(0);
2250: }
2254: /*@
2255: DMHasVariableBounds - does the DM object have a variable bounds function?
2257: Not Collective
2259: Input Parameter:
2260: . dm - the DM object to destroy
2262: Output Parameter:
2263: . flg - PETSC_TRUE if the variable bounds function exists
2265: Level: developer
2267: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext()
2269: @*/
2270: PetscErrorCode DMHasVariableBounds(DM dm,PetscBool *flg)
2271: {
2273: *flg = (dm->ops->computevariablebounds) ? PETSC_TRUE : PETSC_FALSE;
2274: return(0);
2275: }
2279: /*@C
2280: DMComputeVariableBounds - compute variable bounds used by SNESVI.
2282: Logically Collective on DM
2284: Input Parameters:
2285: + dm - the DM object to destroy
2286: - x - current solution at which the bounds are computed
2288: Output parameters:
2289: + xl - lower bound
2290: - xu - upper bound
2292: Level: intermediate
2294: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext()
2296: @*/
2297: PetscErrorCode DMComputeVariableBounds(DM dm, Vec xl, Vec xu)
2298: {
2304: if (dm->ops->computevariablebounds) {
2305: (*dm->ops->computevariablebounds)(dm, xl,xu);
2306: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "This DM is incapable of computing variable bounds.");
2307: return(0);
2308: }
2312: /*@
2313: DMHasColoring - does the DM object have a method of providing a coloring?
2315: Not Collective
2317: Input Parameter:
2318: . dm - the DM object
2320: Output Parameter:
2321: . flg - PETSC_TRUE if the DM has facilities for DMCreateColoring().
2323: Level: developer
2325: .seealso DMHasFunction(), DMCreateColoring()
2327: @*/
2328: PetscErrorCode DMHasColoring(DM dm,PetscBool *flg)
2329: {
2331: *flg = (dm->ops->getcoloring) ? PETSC_TRUE : PETSC_FALSE;
2332: return(0);
2333: }
2335: #undef __FUNCT__
2337: /*@C
2338: DMSetVec - set the vector at which to compute residual, Jacobian and VI bounds, if the problem is nonlinear.
2340: Collective on DM
2342: Input Parameter:
2343: + dm - the DM object
2344: - x - location to compute residual and Jacobian, if NULL is passed to those routines; will be NULL for linear problems.
2346: Level: developer
2348: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext()
2350: @*/
2351: PetscErrorCode DMSetVec(DM dm,Vec x)
2352: {
2356: if (x) {
2357: if (!dm->x) {
2358: DMCreateGlobalVector(dm,&dm->x);
2359: }
2360: VecCopy(x,dm->x);
2361: } else if (dm->x) {
2362: VecDestroy(&dm->x);
2363: }
2364: return(0);
2365: }
2367: PetscFunctionList DMList = NULL;
2368: PetscBool DMRegisterAllCalled = PETSC_FALSE;
2372: /*@C
2373: DMSetType - Builds a DM, for a particular DM implementation.
2375: Collective on DM
2377: Input Parameters:
2378: + dm - The DM object
2379: - method - The name of the DM type
2381: Options Database Key:
2382: . -dm_type <type> - Sets the DM type; use -help for a list of available types
2384: Notes:
2385: See "petsc/include/petscdm.h" for available DM types (for instance, DM1D, DM2D, or DM3D).
2387: Level: intermediate
2389: .keywords: DM, set, type
2390: .seealso: DMGetType(), DMCreate()
2391: @*/
2392: PetscErrorCode DMSetType(DM dm, DMType method)
2393: {
2394: PetscErrorCode (*r)(DM);
2395: PetscBool match;
2400: PetscObjectTypeCompare((PetscObject) dm, method, &match);
2401: if (match) return(0);
2403: if (!DMRegisterAllCalled) {DMRegisterAll();}
2404: PetscFunctionListFind(DMList,method,&r);
2405: if (!r) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown DM type: %s", method);
2407: if (dm->ops->destroy) {
2408: (*dm->ops->destroy)(dm);
2409: dm->ops->destroy = NULL;
2410: }
2411: (*r)(dm);
2412: PetscObjectChangeTypeName((PetscObject)dm,method);
2413: return(0);
2414: }
2418: /*@C
2419: DMGetType - Gets the DM type name (as a string) from the DM.
2421: Not Collective
2423: Input Parameter:
2424: . dm - The DM
2426: Output Parameter:
2427: . type - The DM type name
2429: Level: intermediate
2431: .keywords: DM, get, type, name
2432: .seealso: DMSetType(), DMCreate()
2433: @*/
2434: PetscErrorCode DMGetType(DM dm, DMType *type)
2435: {
2441: if (!DMRegisterAllCalled) {
2442: DMRegisterAll();
2443: }
2444: *type = ((PetscObject)dm)->type_name;
2445: return(0);
2446: }
2450: /*@C
2451: DMConvert - Converts a DM to another DM, either of the same or different type.
2453: Collective on DM
2455: Input Parameters:
2456: + dm - the DM
2457: - newtype - new DM type (use "same" for the same type)
2459: Output Parameter:
2460: . M - pointer to new DM
2462: Notes:
2463: Cannot be used to convert a sequential DM to parallel or parallel to sequential,
2464: the MPI communicator of the generated DM is always the same as the communicator
2465: of the input DM.
2467: Level: intermediate
2469: .seealso: DMCreate()
2470: @*/
2471: PetscErrorCode DMConvert(DM dm, DMType newtype, DM *M)
2472: {
2473: DM B;
2474: char convname[256];
2475: PetscBool sametype, issame;
2482: PetscObjectTypeCompare((PetscObject) dm, newtype, &sametype);
2483: PetscStrcmp(newtype, "same", &issame);
2484: {
2485: PetscErrorCode (*conv)(DM, DMType, DM*) = NULL;
2487: /*
2488: Order of precedence:
2489: 1) See if a specialized converter is known to the current DM.
2490: 2) See if a specialized converter is known to the desired DM class.
2491: 3) See if a good general converter is registered for the desired class
2492: 4) See if a good general converter is known for the current matrix.
2493: 5) Use a really basic converter.
2494: */
2496: /* 1) See if a specialized converter is known to the current DM and the desired class */
2497: PetscStrcpy(convname,"DMConvert_");
2498: PetscStrcat(convname,((PetscObject) dm)->type_name);
2499: PetscStrcat(convname,"_");
2500: PetscStrcat(convname,newtype);
2501: PetscStrcat(convname,"_C");
2502: PetscObjectQueryFunction((PetscObject)dm,convname,&conv);
2503: if (conv) goto foundconv;
2505: /* 2) See if a specialized converter is known to the desired DM class. */
2506: DMCreate(PetscObjectComm((PetscObject)dm), &B);
2507: DMSetType(B, newtype);
2508: PetscStrcpy(convname,"DMConvert_");
2509: PetscStrcat(convname,((PetscObject) dm)->type_name);
2510: PetscStrcat(convname,"_");
2511: PetscStrcat(convname,newtype);
2512: PetscStrcat(convname,"_C");
2513: PetscObjectQueryFunction((PetscObject)B,convname,&conv);
2514: if (conv) {
2515: DMDestroy(&B);
2516: goto foundconv;
2517: }
2519: #if 0
2520: /* 3) See if a good general converter is registered for the desired class */
2521: conv = B->ops->convertfrom;
2522: DMDestroy(&B);
2523: if (conv) goto foundconv;
2525: /* 4) See if a good general converter is known for the current matrix */
2526: if (dm->ops->convert) {
2527: conv = dm->ops->convert;
2528: }
2529: if (conv) goto foundconv;
2530: #endif
2532: /* 5) Use a really basic converter. */
2533: SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "No conversion possible between DM types %s and %s", ((PetscObject) dm)->type_name, newtype);
2535: foundconv:
2536: PetscLogEventBegin(DM_Convert,dm,0,0,0);
2537: (*conv)(dm,newtype,M);
2538: PetscLogEventEnd(DM_Convert,dm,0,0,0);
2539: }
2540: PetscObjectStateIncrease((PetscObject) *M);
2541: return(0);
2542: }
2544: /*--------------------------------------------------------------------------------------------------------------------*/
2548: /*@C
2549: DMRegister - Adds a new DM component implementation
2551: Not Collective
2553: Input Parameters:
2554: + name - The name of a new user-defined creation routine
2555: - create_func - The creation routine itself
2557: Notes:
2558: DMRegister() may be called multiple times to add several user-defined DMs
2561: Sample usage:
2562: .vb
2563: DMRegister("my_da", MyDMCreate);
2564: .ve
2566: Then, your DM type can be chosen with the procedural interface via
2567: .vb
2568: DMCreate(MPI_Comm, DM *);
2569: DMSetType(DM,"my_da");
2570: .ve
2571: or at runtime via the option
2572: .vb
2573: -da_type my_da
2574: .ve
2576: Level: advanced
2578: .keywords: DM, register
2579: .seealso: DMRegisterAll(), DMRegisterDestroy()
2581: @*/
2582: PetscErrorCode DMRegister(const char sname[],PetscErrorCode (*function)(DM))
2583: {
2587: PetscFunctionListAdd(&DMList,sname,function);
2588: return(0);
2589: }
2593: /*@C
2594: DMLoad - Loads a DM that has been stored in binary with DMView().
2596: Collective on PetscViewer
2598: Input Parameters:
2599: + newdm - the newly loaded DM, this needs to have been created with DMCreate() or
2600: some related function before a call to DMLoad().
2601: - viewer - binary file viewer, obtained from PetscViewerBinaryOpen() or
2602: HDF5 file viewer, obtained from PetscViewerHDF5Open()
2604: Level: intermediate
2606: Notes:
2607: The type is determined by the data in the file, any type set into the DM before this call is ignored.
2609: Notes for advanced users:
2610: Most users should not need to know the details of the binary storage
2611: format, since DMLoad() and DMView() completely hide these details.
2612: But for anyone who's interested, the standard binary matrix storage
2613: format is
2614: .vb
2615: has not yet been determined
2616: .ve
2618: .seealso: PetscViewerBinaryOpen(), DMView(), MatLoad(), VecLoad()
2619: @*/
2620: PetscErrorCode DMLoad(DM newdm, PetscViewer viewer)
2621: {
2623: PetscBool isbinary;
2624: PetscInt classid;
2625: char type[256];
2630: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
2631: if (!isbinary) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid viewer; open viewer with PetscViewerBinaryOpen()");
2633: PetscViewerBinaryRead(viewer,&classid,1,PETSC_INT);
2634: if (classid != DM_FILE_CLASSID) SETERRQ(PetscObjectComm((PetscObject)newdm),PETSC_ERR_ARG_WRONG,"Not DM next in file");
2635: PetscViewerBinaryRead(viewer,type,256,PETSC_CHAR);
2636: DMSetType(newdm, type);
2637: if (newdm->ops->load) {
2638: (*newdm->ops->load)(newdm,viewer);
2639: }
2640: return(0);
2641: }
2643: /******************************** FEM Support **********************************/
2647: PetscErrorCode DMPrintCellVector(PetscInt c, const char name[], PetscInt len, const PetscScalar x[])
2648: {
2649: PetscInt f;
2653: PetscPrintf(PETSC_COMM_SELF, "Cell %D Element %s\n", c, name);
2654: for (f = 0; f < len; ++f) {
2655: PetscPrintf(PETSC_COMM_SELF, " | %G |\n", PetscRealPart(x[f]));
2656: }
2657: return(0);
2658: }
2662: PetscErrorCode DMPrintCellMatrix(PetscInt c, const char name[], PetscInt rows, PetscInt cols, const PetscScalar A[])
2663: {
2664: PetscInt f, g;
2668: PetscPrintf(PETSC_COMM_SELF, "Cell %D Element %s\n", c, name);
2669: for (f = 0; f < rows; ++f) {
2670: PetscPrintf(PETSC_COMM_SELF, " |");
2671: for (g = 0; g < cols; ++g) {
2672: PetscPrintf(PETSC_COMM_SELF, " % 9.5G", PetscRealPart(A[f*cols+g]));
2673: }
2674: PetscPrintf(PETSC_COMM_SELF, " |\n");
2675: }
2676: return(0);
2677: }
2681: /*@
2682: DMGetDefaultSection - Get the PetscSection encoding the local data layout for the DM.
2684: Input Parameter:
2685: . dm - The DM
2687: Output Parameter:
2688: . section - The PetscSection
2690: Level: intermediate
2692: Note: This gets a borrowed reference, so the user should not destroy this PetscSection.
2694: .seealso: DMSetDefaultSection(), DMGetDefaultGlobalSection()
2695: @*/
2696: PetscErrorCode DMGetDefaultSection(DM dm, PetscSection *section)
2697: {
2701: *section = dm->defaultSection;
2702: return(0);
2703: }
2707: /*@
2708: DMSetDefaultSection - Set the PetscSection encoding the local data layout for the DM.
2710: Input Parameters:
2711: + dm - The DM
2712: - section - The PetscSection
2714: Level: intermediate
2716: Note: Any existing Section will be destroyed
2718: .seealso: DMSetDefaultSection(), DMGetDefaultGlobalSection()
2719: @*/
2720: PetscErrorCode DMSetDefaultSection(DM dm, PetscSection section)
2721: {
2722: PetscInt numFields;
2723: PetscInt f;
2729: PetscObjectReference((PetscObject)section);
2730: PetscSectionDestroy(&dm->defaultSection);
2731: dm->defaultSection = section;
2732: PetscSectionGetNumFields(dm->defaultSection, &numFields);
2733: if (numFields) {
2734: DMSetNumFields(dm, numFields);
2735: for (f = 0; f < numFields; ++f) {
2736: const char *name;
2738: PetscSectionGetFieldName(dm->defaultSection, f, &name);
2739: PetscObjectSetName(dm->fields[f], name);
2740: }
2741: }
2742: /* The global section will be rebuilt in the next call to DMGetDefaultGlobalSection(). */
2743: PetscSectionDestroy(&dm->defaultGlobalSection);
2744: return(0);
2745: }
2749: /*@
2750: DMGetDefaultGlobalSection - Get the PetscSection encoding the global data layout for the DM.
2752: Collective on DM
2754: Input Parameter:
2755: . dm - The DM
2757: Output Parameter:
2758: . section - The PetscSection
2760: Level: intermediate
2762: Note: This gets a borrowed reference, so the user should not destroy this PetscSection.
2764: .seealso: DMSetDefaultSection(), DMGetDefaultSection()
2765: @*/
2766: PetscErrorCode DMGetDefaultGlobalSection(DM dm, PetscSection *section)
2767: {
2773: if (!dm->defaultGlobalSection) {
2774: if (!dm->defaultSection || !dm->sf) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DM must have a default PetscSection and PetscSF in order to create a global PetscSection");
2775: PetscSectionCreateGlobalSection(dm->defaultSection, dm->sf, PETSC_FALSE, &dm->defaultGlobalSection);
2776: PetscLayoutDestroy(&dm->map);
2777: PetscSectionGetValueLayout(PetscObjectComm((PetscObject)dm),dm->defaultGlobalSection,&dm->map);
2778: }
2779: *section = dm->defaultGlobalSection;
2780: return(0);
2781: }
2785: /*@
2786: DMSetDefaultGlobalSection - Set the PetscSection encoding the global data layout for the DM.
2788: Input Parameters:
2789: + dm - The DM
2790: - section - The PetscSection, or NULL
2792: Level: intermediate
2794: Note: Any existing Section will be destroyed
2796: .seealso: DMGetDefaultGlobalSection(), DMSetDefaultSection()
2797: @*/
2798: PetscErrorCode DMSetDefaultGlobalSection(DM dm, PetscSection section)
2799: {
2805: PetscObjectReference((PetscObject)section);
2806: PetscSectionDestroy(&dm->defaultGlobalSection);
2807: dm->defaultGlobalSection = section;
2808: return(0);
2809: }
2813: /*@
2814: DMGetDefaultSF - Get the PetscSF encoding the parallel dof overlap for the DM. If it has not been set,
2815: it is created from the default PetscSection layouts in the DM.
2817: Input Parameter:
2818: . dm - The DM
2820: Output Parameter:
2821: . sf - The PetscSF
2823: Level: intermediate
2825: Note: This gets a borrowed reference, so the user should not destroy this PetscSF.
2827: .seealso: DMSetDefaultSF(), DMCreateDefaultSF()
2828: @*/
2829: PetscErrorCode DMGetDefaultSF(DM dm, PetscSF *sf)
2830: {
2831: PetscInt nroots;
2837: PetscSFGetGraph(dm->defaultSF, &nroots, NULL, NULL, NULL);
2838: if (nroots < 0) {
2839: PetscSection section, gSection;
2841: DMGetDefaultSection(dm, §ion);
2842: if (section) {
2843: DMGetDefaultGlobalSection(dm, &gSection);
2844: DMCreateDefaultSF(dm, section, gSection);
2845: } else {
2846: *sf = NULL;
2847: return(0);
2848: }
2849: }
2850: *sf = dm->defaultSF;
2851: return(0);
2852: }
2856: /*@
2857: DMSetDefaultSF - Set the PetscSF encoding the parallel dof overlap for the DM
2859: Input Parameters:
2860: + dm - The DM
2861: - sf - The PetscSF
2863: Level: intermediate
2865: Note: Any previous SF is destroyed
2867: .seealso: DMGetDefaultSF(), DMCreateDefaultSF()
2868: @*/
2869: PetscErrorCode DMSetDefaultSF(DM dm, PetscSF sf)
2870: {
2876: PetscSFDestroy(&dm->defaultSF);
2877: dm->defaultSF = sf;
2878: return(0);
2879: }
2883: /*@C
2884: DMCreateDefaultSF - Create the PetscSF encoding the parallel dof overlap for the DM based upon the PetscSections
2885: describing the data layout.
2887: Input Parameters:
2888: + dm - The DM
2889: . localSection - PetscSection describing the local data layout
2890: - globalSection - PetscSection describing the global data layout
2892: Level: intermediate
2894: .seealso: DMGetDefaultSF(), DMSetDefaultSF()
2895: @*/
2896: PetscErrorCode DMCreateDefaultSF(DM dm, PetscSection localSection, PetscSection globalSection)
2897: {
2898: MPI_Comm comm;
2899: PetscLayout layout;
2900: const PetscInt *ranges;
2901: PetscInt *local;
2902: PetscSFNode *remote;
2903: PetscInt pStart, pEnd, p, nroots, nleaves, l;
2904: PetscMPIInt size, rank;
2908: PetscObjectGetComm((PetscObject)dm,&comm);
2910: MPI_Comm_size(comm, &size);
2911: MPI_Comm_rank(comm, &rank);
2912: PetscSectionGetChart(globalSection, &pStart, &pEnd);
2913: PetscSectionGetConstrainedStorageSize(globalSection, &nroots);
2914: PetscLayoutCreate(comm, &layout);
2915: PetscLayoutSetBlockSize(layout, 1);
2916: PetscLayoutSetLocalSize(layout, nroots);
2917: PetscLayoutSetUp(layout);
2918: PetscLayoutGetRanges(layout, &ranges);
2919: for (p = pStart, nleaves = 0; p < pEnd; ++p) {
2920: PetscInt gdof, gcdof;
2922: PetscSectionGetDof(globalSection, p, &gdof);
2923: PetscSectionGetConstraintDof(globalSection, p, &gcdof);
2924: nleaves += gdof < 0 ? -(gdof+1)-gcdof : gdof-gcdof;
2925: }
2926: PetscMalloc(nleaves * sizeof(PetscInt), &local);
2927: PetscMalloc(nleaves * sizeof(PetscSFNode), &remote);
2928: for (p = pStart, l = 0; p < pEnd; ++p) {
2929: const PetscInt *cind;
2930: PetscInt dof, cdof, off, gdof, gcdof, goff, gsize, d, c;
2932: PetscSectionGetDof(localSection, p, &dof);
2933: PetscSectionGetOffset(localSection, p, &off);
2934: PetscSectionGetConstraintDof(localSection, p, &cdof);
2935: PetscSectionGetConstraintIndices(localSection, p, &cind);
2936: PetscSectionGetDof(globalSection, p, &gdof);
2937: PetscSectionGetConstraintDof(globalSection, p, &gcdof);
2938: PetscSectionGetOffset(globalSection, p, &goff);
2939: if (!gdof) continue; /* Censored point */
2940: gsize = gdof < 0 ? -(gdof+1)-gcdof : gdof-gcdof;
2941: if (gsize != dof-cdof) {
2942: 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);
2943: cdof = 0; /* Ignore constraints */
2944: }
2945: for (d = 0, c = 0; d < dof; ++d) {
2946: if ((c < cdof) && (cind[c] == d)) {++c; continue;}
2947: local[l+d-c] = off+d;
2948: }
2949: if (gdof < 0) {
2950: for (d = 0; d < gsize; ++d, ++l) {
2951: PetscInt offset = -(goff+1) + d, r;
2953: PetscFindInt(offset,size,ranges,&r);
2954: if (r < 0) r = -(r+2);
2955: remote[l].rank = r;
2956: remote[l].index = offset - ranges[r];
2957: }
2958: } else {
2959: for (d = 0; d < gsize; ++d, ++l) {
2960: remote[l].rank = rank;
2961: remote[l].index = goff+d - ranges[rank];
2962: }
2963: }
2964: }
2965: if (l != nleaves) SETERRQ2(comm, PETSC_ERR_PLIB, "Iteration error, l %d != nleaves %d", l, nleaves);
2966: PetscLayoutDestroy(&layout);
2967: PetscSFSetGraph(dm->defaultSF, nroots, nleaves, local, PETSC_OWN_POINTER, remote, PETSC_OWN_POINTER);
2968: return(0);
2969: }
2973: /*@
2974: DMGetPointSF - Get the PetscSF encoding the parallel section point overlap for the DM.
2976: Input Parameter:
2977: . dm - The DM
2979: Output Parameter:
2980: . sf - The PetscSF
2982: Level: intermediate
2984: Note: This gets a borrowed reference, so the user should not destroy this PetscSF.
2986: .seealso: DMSetPointSF(), DMGetDefaultSF(), DMSetDefaultSF(), DMCreateDefaultSF()
2987: @*/
2988: PetscErrorCode DMGetPointSF(DM dm, PetscSF *sf)
2989: {
2993: *sf = dm->sf;
2994: return(0);
2995: }
2999: /*@
3000: DMSetPointSF - Set the PetscSF encoding the parallel section point overlap for the DM.
3002: Input Parameters:
3003: + dm - The DM
3004: - sf - The PetscSF
3006: Level: intermediate
3008: .seealso: DMGetPointSF(), DMGetDefaultSF(), DMSetDefaultSF(), DMCreateDefaultSF()
3009: @*/
3010: PetscErrorCode DMSetPointSF(DM dm, PetscSF sf)
3011: {
3017: PetscSFDestroy(&dm->sf);
3018: PetscObjectReference((PetscObject) sf);
3019: dm->sf = sf;
3020: return(0);
3021: }
3025: PetscErrorCode DMGetNumFields(DM dm, PetscInt *numFields)
3026: {
3030: *numFields = dm->numFields;
3031: return(0);
3032: }
3036: PetscErrorCode DMSetNumFields(DM dm, PetscInt numFields)
3037: {
3038: PetscInt f;
3043: for (f = 0; f < dm->numFields; ++f) {
3044: PetscObjectDestroy(&dm->fields[f]);
3045: }
3046: PetscFree(dm->fields);
3047: dm->numFields = numFields;
3048: PetscMalloc(dm->numFields * sizeof(PetscObject), &dm->fields);
3049: for (f = 0; f < dm->numFields; ++f) {
3050: PetscContainerCreate(PetscObjectComm((PetscObject)dm), (PetscContainer*) &dm->fields[f]);
3051: }
3052: return(0);
3053: }
3057: PetscErrorCode DMGetField(DM dm, PetscInt f, PetscObject *field)
3058: {
3062: if (!dm->fields) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "Fields have not been setup in this DM. Call DMSetNumFields()");
3063: if ((f < 0) || (f >= dm->numFields)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Field %d should be in [%d,%d)", f, 0, dm->numFields);
3064: *field = dm->fields[f];
3065: return(0);
3066: }
3070: /*@
3071: DMSetCoordinates - Sets into the DM a global vector that holds the coordinates
3073: Collective on DM
3075: Input Parameters:
3076: + dm - the DM
3077: - c - coordinate vector
3079: Note:
3080: The coordinates do include those for ghost points, which are in the local vector
3082: Level: intermediate
3084: .keywords: distributed array, get, corners, nodes, local indices, coordinates
3085: .seealso: DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLoca(), DMGetCoordinateDM()
3086: @*/
3087: PetscErrorCode DMSetCoordinates(DM dm, Vec c)
3088: {
3094: PetscObjectReference((PetscObject) c);
3095: VecDestroy(&dm->coordinates);
3096: dm->coordinates = c;
3097: VecDestroy(&dm->coordinatesLocal);
3098: return(0);
3099: }
3103: /*@
3104: DMSetCoordinatesLocal - Sets into the DM a local vector that holds the coordinates
3106: Collective on DM
3108: Input Parameters:
3109: + dm - the DM
3110: - c - coordinate vector
3112: Note:
3113: The coordinates of ghost points can be set using DMSetCoordinates()
3114: followed by DMGetCoordinatesLocal(). This is intended to enable the
3115: setting of ghost coordinates outside of the domain.
3117: Level: intermediate
3119: .keywords: distributed array, get, corners, nodes, local indices, coordinates
3120: .seealso: DMGetCoordinatesLocal(), DMSetCoordinates(), DMGetCoordinates(), DMGetCoordinateDM()
3121: @*/
3122: PetscErrorCode DMSetCoordinatesLocal(DM dm, Vec c)
3123: {
3129: PetscObjectReference((PetscObject) c);
3130: VecDestroy(&dm->coordinatesLocal);
3132: dm->coordinatesLocal = c;
3134: VecDestroy(&dm->coordinates);
3135: return(0);
3136: }
3140: /*@
3141: DMGetCoordinates - Gets a global vector with the coordinates associated with the DM.
3143: Not Collective
3145: Input Parameter:
3146: . dm - the DM
3148: Output Parameter:
3149: . c - global coordinate vector
3151: Note:
3152: This is a borrowed reference, so the user should NOT destroy this vector
3154: Each process has only the local coordinates (does NOT have the ghost coordinates).
3156: For DMDA, in two and three dimensions coordinates are interlaced (x_0,y_0,x_1,y_1,...)
3157: and (x_0,y_0,z_0,x_1,y_1,z_1...)
3159: Level: intermediate
3161: .keywords: distributed array, get, corners, nodes, local indices, coordinates
3162: .seealso: DMSetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM()
3163: @*/
3164: PetscErrorCode DMGetCoordinates(DM dm, Vec *c)
3165: {
3171: if (!dm->coordinates && dm->coordinatesLocal) {
3172: DM cdm = NULL;
3174: DMGetCoordinateDM(dm, &cdm);
3175: DMCreateGlobalVector(cdm, &dm->coordinates);
3176: PetscObjectSetName((PetscObject) dm->coordinates, "coordinates");
3177: DMLocalToGlobalBegin(cdm, dm->coordinatesLocal, INSERT_VALUES, dm->coordinates);
3178: DMLocalToGlobalEnd(cdm, dm->coordinatesLocal, INSERT_VALUES, dm->coordinates);
3179: }
3180: *c = dm->coordinates;
3181: return(0);
3182: }
3186: /*@
3187: DMGetCoordinatesLocal - Gets a local vector with the coordinates associated with the DM.
3189: Collective on DM
3191: Input Parameter:
3192: . dm - the DM
3194: Output Parameter:
3195: . c - coordinate vector
3197: Note:
3198: This is a borrowed reference, so the user should NOT destroy this vector
3200: Each process has the local and ghost coordinates
3202: For DMDA, in two and three dimensions coordinates are interlaced (x_0,y_0,x_1,y_1,...)
3203: and (x_0,y_0,z_0,x_1,y_1,z_1...)
3205: Level: intermediate
3207: .keywords: distributed array, get, corners, nodes, local indices, coordinates
3208: .seealso: DMSetCoordinatesLocal(), DMGetCoordinates(), DMSetCoordinates(), DMGetCoordinateDM()
3209: @*/
3210: PetscErrorCode DMGetCoordinatesLocal(DM dm, Vec *c)
3211: {
3217: if (!dm->coordinatesLocal && dm->coordinates) {
3218: DM cdm = NULL;
3220: DMGetCoordinateDM(dm, &cdm);
3221: DMCreateLocalVector(cdm, &dm->coordinatesLocal);
3222: PetscObjectSetName((PetscObject) dm->coordinatesLocal, "coordinates");
3223: DMGlobalToLocalBegin(cdm, dm->coordinates, INSERT_VALUES, dm->coordinatesLocal);
3224: DMGlobalToLocalEnd(cdm, dm->coordinates, INSERT_VALUES, dm->coordinatesLocal);
3225: }
3226: *c = dm->coordinatesLocal;
3227: return(0);
3228: }
3232: /*@
3233: DMGetCoordinateDM - Gets the DM that scatters between global and local coordinates
3235: Collective on DM
3237: Input Parameter:
3238: . dm - the DM
3240: Output Parameter:
3241: . cdm - coordinate DM
3243: Level: intermediate
3245: .keywords: distributed array, get, corners, nodes, local indices, coordinates
3246: .seealso: DMSetCoordinates(), DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLocal()
3247: @*/
3248: PetscErrorCode DMGetCoordinateDM(DM dm, DM *cdm)
3249: {
3255: if (!dm->coordinateDM) {
3256: if (!dm->ops->createcoordinatedm) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unable to create coordinates for this DM");
3257: (*dm->ops->createcoordinatedm)(dm, &dm->coordinateDM);
3258: }
3259: *cdm = dm->coordinateDM;
3260: return(0);
3261: }
3265: /*@
3266: DMLocatePoints - Locate the points in v in the mesh and return an IS of the containing cells
3268: Not collective
3270: Input Parameters:
3271: + dm - The DM
3272: - v - The Vec of points
3274: Output Parameter:
3275: . cells - The local cell numbers for cells which contain the points
3277: Level: developer
3279: .keywords: point location, mesh
3280: .seealso: DMSetCoordinates(), DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLocal()
3281: @*/
3282: PetscErrorCode DMLocatePoints(DM dm, Vec v, IS *cells)
3283: {
3290: if (dm->ops->locatepoints) {
3291: (*dm->ops->locatepoints)(dm,v,cells);
3292: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Point location not available for this DM");
3293: return(0);
3294: }