Actual source code: dm.c
petsc-3.3-p7 2013-05-11
1: #include <petscsnes.h>
2: #include <petsc-private/dmimpl.h> /*I "petscdm.h" I*/
4: PetscClassId DM_CLASSID;
5: PetscLogEvent DM_Convert, DM_GlobalToLocal, DM_LocalToGlobal;
9: /*@
10: DMCreate - Creates an empty DM object. The type can then be set with DMSetType().
12: If you never call DMSetType() it will generate an
13: error when you try to use the vector.
15: Collective on MPI_Comm
17: Input Parameter:
18: . comm - The communicator for the DM object
20: Output Parameter:
21: . dm - The DM object
23: Level: beginner
25: .seealso: DMSetType(), DMDA, DMSLICED, DMCOMPOSITE
26: @*/
27: PetscErrorCode DMCreate(MPI_Comm comm,DM *dm)
28: {
29: DM v;
34: *dm = PETSC_NULL;
35: #ifndef PETSC_USE_DYNAMIC_LIBRARIES
36: VecInitializePackage(PETSC_NULL);
37: MatInitializePackage(PETSC_NULL);
38: DMInitializePackage(PETSC_NULL);
39: #endif
41: PetscHeaderCreate(v, _p_DM, struct _DMOps, DM_CLASSID, -1, "DM", "Distribution Manager", "DM", comm, DMDestroy, DMView);
42: PetscMemzero(v->ops, sizeof(struct _DMOps));
45: v->workSize = 0;
46: v->workArray = PETSC_NULL;
47: v->ltogmap = PETSC_NULL;
48: v->ltogmapb = PETSC_NULL;
49: v->bs = 1;
50: v->coloringtype = IS_COLORING_GLOBAL;
51: v->lf = PETSC_NULL;
52: v->lj = PETSC_NULL;
53: PetscSFCreate(comm, &v->sf);
54: PetscSFCreate(comm, &v->defaultSF);
55: v->defaultSection = PETSC_NULL;
56: v->defaultGlobalSection = PETSC_NULL;
58: *dm = v;
59: return(0);
60: }
65: /*@C
66: DMSetVecType - Sets the type of vector created with DMCreateLocalVector() and DMCreateGlobalVector()
68: Logically Collective on DMDA
70: Input Parameter:
71: + da - initial distributed array
72: . ctype - the vector type, currently either VECSTANDARD or VECCUSP
74: Options Database:
75: . -dm_vec_type ctype
77: Level: intermediate
79: .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMDestroy(), DMDA, DMDAInterpolationType, VecType
80: @*/
81: PetscErrorCode DMSetVecType(DM da,const VecType ctype)
82: {
87: PetscFree(da->vectype);
88: PetscStrallocpy(ctype,&da->vectype);
89: return(0);
90: }
94: /*@C
95: DMSetMatType - Sets the type of matrix created with DMCreateMatrix()
97: Logically Collective on DM
99: Input Parameter:
100: + dm - the DM context
101: . ctype - the matrix type
103: Options Database:
104: . -dm_mat_type ctype
106: Level: intermediate
108: .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMCreateMatrix(), DMSetMatrixPreallocateOnly(), MatType
109: @*/
110: PetscErrorCode DMSetMatType(DM dm,const MatType ctype)
111: {
115: PetscFree(dm->mattype);
116: PetscStrallocpy(ctype,&dm->mattype);
117: return(0);
118: }
122: /*@C
123: DMSetOptionsPrefix - Sets the prefix used for searching for all
124: DMDA options in the database.
126: Logically Collective on DMDA
128: Input Parameter:
129: + da - the DMDA context
130: - prefix - the prefix to prepend to all option names
132: Notes:
133: A hyphen (-) must NOT be given at the beginning of the prefix name.
134: The first character of all runtime options is AUTOMATICALLY the hyphen.
136: Level: advanced
138: .keywords: DMDA, set, options, prefix, database
140: .seealso: DMSetFromOptions()
141: @*/
142: PetscErrorCode DMSetOptionsPrefix(DM dm,const char prefix[])
143: {
148: PetscObjectSetOptionsPrefix((PetscObject)dm,prefix);
149: return(0);
150: }
154: /*@
155: DMDestroy - Destroys a vector packer or DMDA.
157: Collective on DM
159: Input Parameter:
160: . dm - the DM object to destroy
162: Level: developer
164: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix()
166: @*/
167: PetscErrorCode DMDestroy(DM *dm)
168: {
169: PetscInt i, cnt = 0;
170: DMNamedVecLink nlink,nnext;
174: if (!*dm) return(0);
177: /* count all the circular references of DM and its contained Vecs */
178: for (i=0; i<DM_MAX_WORK_VECTORS; i++) {
179: if ((*dm)->localin[i]) {cnt++;}
180: if ((*dm)->globalin[i]) {cnt++;}
181: }
182: for (nlink=(*dm)->namedglobal; nlink; nlink=nlink->next) cnt++;
183: if ((*dm)->x) {
184: PetscObject obj;
185: PetscObjectQuery((PetscObject)(*dm)->x,"DM",&obj);
186: if (obj == (PetscObject)*dm) cnt++;
187: }
189: if (--((PetscObject)(*dm))->refct - cnt > 0) {*dm = 0; return(0);}
190: /*
191: Need this test because the dm references the vectors that
192: reference the dm, so destroying the dm calls destroy on the
193: vectors that cause another destroy on the dm
194: */
195: if (((PetscObject)(*dm))->refct < 0) return(0);
196: ((PetscObject) (*dm))->refct = 0;
197: for (i=0; i<DM_MAX_WORK_VECTORS; i++) {
198: if ((*dm)->localout[i]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Destroying a DM that has a local vector obtained with DMGetLocalVector()");
199: VecDestroy(&(*dm)->localin[i]);
200: }
201: for (nlink=(*dm)->namedglobal; nlink; nlink=nnext) { /* Destroy the named vectors */
202: nnext = nlink->next;
203: if (nlink->status != DMVEC_STATUS_IN) SETERRQ1(((PetscObject)*dm)->comm,PETSC_ERR_ARG_WRONGSTATE,"DM still has Vec named '%s' checked out",nlink->name);
204: PetscFree(nlink->name);
205: VecDestroy(&nlink->X);
206: PetscFree(nlink);
207: }
208: (*dm)->namedglobal = PETSC_NULL;
210: /* Destroy the list of hooks */
211: {
212: DMCoarsenHookLink link,next;
213: for (link=(*dm)->coarsenhook; link; link=next) {
214: next = link->next;
215: PetscFree(link);
216: }
217: (*dm)->coarsenhook = PETSC_NULL;
218: }
219: {
220: DMRefineHookLink link,next;
221: for (link=(*dm)->refinehook; link; link=next) {
222: next = link->next;
223: PetscFree(link);
224: }
225: (*dm)->refinehook = PETSC_NULL;
226: }
228: if ((*dm)->ctx && (*dm)->ctxdestroy) {
229: (*(*dm)->ctxdestroy)(&(*dm)->ctx);
230: }
231: VecDestroy(&(*dm)->x);
232: MatFDColoringDestroy(&(*dm)->fd);
233: DMClearGlobalVectors(*dm);
234: ISLocalToGlobalMappingDestroy(&(*dm)->ltogmap);
235: ISLocalToGlobalMappingDestroy(&(*dm)->ltogmapb);
236: PetscFree((*dm)->vectype);
237: PetscFree((*dm)->mattype);
238: PetscFree((*dm)->workArray);
240: PetscSectionDestroy(&(*dm)->defaultSection);
241: PetscSectionDestroy(&(*dm)->defaultGlobalSection);
242: PetscSFDestroy(&(*dm)->sf);
243: PetscSFDestroy(&(*dm)->defaultSF);
244: /* if memory was published with AMS then destroy it */
245: PetscObjectDepublish(*dm);
247: (*(*dm)->ops->destroy)(*dm);
248: PetscFree((*dm)->data);
249: PetscHeaderDestroy(dm);
250: return(0);
251: }
255: /*@
256: DMSetUp - sets up the data structures inside a DM object
258: Collective on DM
260: Input Parameter:
261: . dm - the DM object to setup
263: Level: developer
265: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix()
267: @*/
268: PetscErrorCode DMSetUp(DM dm)
269: {
274: if (dm->setupcalled) return(0);
275: if (dm->ops->setup) {
276: (*dm->ops->setup)(dm);
277: }
278: dm->setupcalled = PETSC_TRUE;
279: return(0);
280: }
284: /*@
285: DMSetFromOptions - sets parameters in a DM from the options database
287: Collective on DM
289: Input Parameter:
290: . dm - the DM object to set options for
292: Options Database:
293: + -dm_preallocate_only: Only preallocate the matrix for DMCreateMatrix(), but do not fill it with zeros
294: . -dm_vec_type <type> type of vector to create inside DM
295: . -dm_mat_type <type> type of matrix to create inside DM
296: - -dm_coloring_type <global or ghosted>
298: Level: developer
300: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix()
302: @*/
303: PetscErrorCode DMSetFromOptions(DM dm)
304: {
305: PetscBool flg1 = PETSC_FALSE,flg2 = PETSC_FALSE,flg3 = PETSC_FALSE,flg4 = PETSC_FALSE,flg;
307: char typeName[256] = MATAIJ;
311: PetscObjectOptionsBegin((PetscObject)dm);
312: PetscOptionsBool("-dm_view", "Information on DM", "DMView", flg1, &flg1, PETSC_NULL);
313: PetscOptionsBool("-dm_view_detail", "Exhaustive mesh description", "DMView", flg2, &flg2, PETSC_NULL);
314: PetscOptionsBool("-dm_view_vtk", "Output mesh in VTK format", "DMView", flg3, &flg3, PETSC_NULL);
315: PetscOptionsBool("-dm_view_latex", "Output mesh in LaTeX TikZ format", "DMView", flg4, &flg4, PETSC_NULL);
316: PetscOptionsBool("-dm_preallocate_only","only preallocate matrix, but do not set column indices","DMSetMatrixPreallocateOnly",dm->prealloc_only,&dm->prealloc_only,PETSC_NULL);
317: PetscOptionsList("-dm_vec_type","Vector type used for created vectors","DMSetVecType",VecList,dm->vectype,typeName,256,&flg);
318: if (flg) {
319: DMSetVecType(dm,typeName);
320: }
321: PetscOptionsList("-dm_mat_type","Matrix type used for created matrices","DMSetMatType",MatList,dm->mattype?dm->mattype:typeName,typeName,sizeof typeName,&flg);
322: if (flg) {
323: DMSetMatType(dm,typeName);
324: }
325: PetscOptionsEnum("-dm_is_coloring_type","Global or local coloring of Jacobian","ISColoringType",ISColoringTypes,(PetscEnum)dm->coloringtype,(PetscEnum*)&dm->coloringtype,PETSC_NULL);
326: if (dm->ops->setfromoptions) {
327: (*dm->ops->setfromoptions)(dm);
328: }
329: /* process any options handlers added with PetscObjectAddOptionsHandler() */
330: PetscObjectProcessOptionsHandlers((PetscObject) dm);
331: PetscOptionsEnd();
332: if (flg1) {
333: DMView(dm, PETSC_VIEWER_STDOUT_WORLD);
334: }
335: if (flg2) {
336: PetscViewer viewer;
338: PetscViewerCreate(((PetscObject) dm)->comm, &viewer);
339: PetscViewerSetType(viewer, PETSCVIEWERASCII);
340: PetscViewerSetFormat(viewer, PETSC_VIEWER_ASCII_INFO_DETAIL);
341: DMView(dm, viewer);
342: PetscViewerDestroy(&viewer);
343: }
344: if (flg3) {
345: PetscViewer viewer;
347: PetscViewerCreate(((PetscObject) dm)->comm, &viewer);
348: PetscViewerSetType(viewer, PETSCVIEWERASCII);
349: PetscViewerSetFormat(viewer, PETSC_VIEWER_ASCII_VTK);
350: PetscViewerFileSetName(viewer, "mesh.vtk");
351: DMView(dm, viewer);
352: PetscViewerDestroy(&viewer);
353: }
354: if (flg4) {
355: PetscViewer viewer;
357: PetscViewerCreate(((PetscObject) dm)->comm, &viewer);
358: PetscViewerSetType(viewer, PETSCVIEWERASCII);
359: PetscViewerSetFormat(viewer, PETSC_VIEWER_ASCII_LATEX);
360: PetscViewerFileSetName(viewer, "mesh.tex");
361: DMView(dm, viewer);
362: PetscViewerDestroy(&viewer);
363: }
364: return(0);
365: }
369: /*@C
370: DMView - Views a vector packer or DMDA.
372: Collective on DM
374: Input Parameter:
375: + dm - the DM object to view
376: - v - the viewer
378: Level: developer
380: .seealso DMDestroy(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix()
382: @*/
383: PetscErrorCode DMView(DM dm,PetscViewer v)
384: {
389: if (!v) {
390: PetscViewerASCIIGetStdout(((PetscObject)dm)->comm,&v);
391: }
392: if (dm->ops->view) {
393: (*dm->ops->view)(dm,v);
394: }
395: return(0);
396: }
400: /*@
401: DMCreateGlobalVector - Creates a global vector from a DMDA or DMComposite object
403: Collective on DM
405: Input Parameter:
406: . dm - the DM object
408: Output Parameter:
409: . vec - the global vector
411: Level: beginner
413: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix()
415: @*/
416: PetscErrorCode DMCreateGlobalVector(DM dm,Vec *vec)
417: {
422: if (dm->defaultSection) {
423: PetscSection gSection;
424: PetscInt localSize;
426: DMGetDefaultGlobalSection(dm, &gSection);
427: PetscSectionGetConstrainedStorageSize(dm->defaultGlobalSection, &localSize);
428: VecCreate(((PetscObject) dm)->comm, vec);
429: VecSetSizes(*vec, localSize, PETSC_DETERMINE);
430: /* VecSetType(*vec, dm->vectype); */
431: VecSetFromOptions(*vec);
432: PetscObjectCompose((PetscObject) *vec, "DM", (PetscObject) dm);
433: /* VecSetLocalToGlobalMapping(*vec, dm->ltogmap); */
434: /* VecSetLocalToGlobalMappingBlock(*vec, dm->ltogmapb); */
435: /* VecSetOperation(*vec, VECOP_DUPLICATE, (void(*)(void)) VecDuplicate_MPI_DM); */
436: } else {
437: (*dm->ops->createglobalvector)(dm,vec);
438: }
439: return(0);
440: }
444: /*@
445: DMCreateLocalVector - Creates a local vector from a DMDA or DMComposite object
447: Not Collective
449: Input Parameter:
450: . dm - the DM object
452: Output Parameter:
453: . vec - the local vector
455: Level: beginner
457: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix()
459: @*/
460: PetscErrorCode DMCreateLocalVector(DM dm,Vec *vec)
461: {
466: if (dm->defaultSection) {
467: PetscInt localSize;
469: PetscSectionGetStorageSize(dm->defaultSection, &localSize);
470: VecCreate(PETSC_COMM_SELF, vec);
471: VecSetSizes(*vec, localSize, localSize);
472: VecSetFromOptions(*vec);
473: PetscObjectCompose((PetscObject) *vec, "DM", (PetscObject) dm);
474: } else {
475: (*dm->ops->createlocalvector)(dm,vec);
476: }
477: return(0);
478: }
482: /*@
483: DMGetLocalToGlobalMapping - Accesses the local-to-global mapping in a DM.
485: Collective on DM
487: Input Parameter:
488: . dm - the DM that provides the mapping
490: Output Parameter:
491: . ltog - the mapping
493: Level: intermediate
495: Notes:
496: This mapping can then be used by VecSetLocalToGlobalMapping() or
497: MatSetLocalToGlobalMapping().
499: .seealso: DMCreateLocalVector(), DMGetLocalToGlobalMappingBlock()
500: @*/
501: PetscErrorCode DMGetLocalToGlobalMapping(DM dm,ISLocalToGlobalMapping *ltog)
502: {
508: if (!dm->ltogmap) {
509: PetscSection section, sectionGlobal;
511: DMGetDefaultSection(dm, §ion);
512: if (section) {
513: PetscInt *ltog;
514: PetscInt pStart, pEnd, size, p, l;
516: DMGetDefaultGlobalSection(dm, §ionGlobal);
517: PetscSectionGetChart(section, &pStart, &pEnd);
518: PetscSectionGetStorageSize(section, &size);
519: PetscMalloc(size * sizeof(PetscInt), <og); /* We want the local+overlap size */
520: for(p = pStart, l = 0; p < pEnd; ++p) {
521: PetscInt dof, off, c;
523: /* Should probably use constrained dofs */
524: PetscSectionGetDof(section, p, &dof);
525: PetscSectionGetOffset(sectionGlobal, p, &off);
526: for(c = 0; c < dof; ++c, ++l) {
527: ltog[l] = off+c;
528: }
529: }
530: ISLocalToGlobalMappingCreate(PETSC_COMM_SELF, size, ltog, PETSC_OWN_POINTER, &dm->ltogmap);
531: PetscLogObjectParent(dm, dm->ltogmap);
532: } else {
533: if (!dm->ops->createlocaltoglobalmapping) SETERRQ(((PetscObject)dm)->comm,PETSC_ERR_SUP,"DM can not create LocalToGlobalMapping");
534: (*dm->ops->createlocaltoglobalmapping)(dm);
535: }
536: }
537: *ltog = dm->ltogmap;
538: return(0);
539: }
543: /*@
544: DMGetLocalToGlobalMappingBlock - Accesses the blocked local-to-global mapping in a DM.
546: Collective on DM
548: Input Parameter:
549: . da - the distributed array that provides the mapping
551: Output Parameter:
552: . ltog - the block mapping
554: Level: intermediate
556: Notes:
557: This mapping can then be used by VecSetLocalToGlobalMappingBlock() or
558: MatSetLocalToGlobalMappingBlock().
560: .seealso: DMCreateLocalVector(), DMGetLocalToGlobalMapping(), DMGetBlockSize(), VecSetBlockSize(), MatSetBlockSize()
561: @*/
562: PetscErrorCode DMGetLocalToGlobalMappingBlock(DM dm,ISLocalToGlobalMapping *ltog)
563: {
569: if (!dm->ltogmapb) {
570: PetscInt bs;
571: DMGetBlockSize(dm,&bs);
572: if (bs > 1) {
573: if (!dm->ops->createlocaltoglobalmappingblock) SETERRQ(((PetscObject)dm)->comm,PETSC_ERR_SUP,"DM can not create LocalToGlobalMappingBlock");
574: (*dm->ops->createlocaltoglobalmappingblock)(dm);
575: } else {
576: DMGetLocalToGlobalMapping(dm,&dm->ltogmapb);
577: PetscObjectReference((PetscObject)dm->ltogmapb);
578: }
579: }
580: *ltog = dm->ltogmapb;
581: return(0);
582: }
586: /*@
587: DMGetBlockSize - Gets the inherent block size associated with a DM
589: Not Collective
591: Input Parameter:
592: . dm - the DM with block structure
594: Output Parameter:
595: . bs - the block size, 1 implies no exploitable block structure
597: Level: intermediate
599: .seealso: ISCreateBlock(), VecSetBlockSize(), MatSetBlockSize(), DMGetLocalToGlobalMappingBlock()
600: @*/
601: PetscErrorCode DMGetBlockSize(DM dm,PetscInt *bs)
602: {
606: if (dm->bs < 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"DM does not have enough information to provide a block size yet");
607: *bs = dm->bs;
608: return(0);
609: }
613: /*@
614: DMCreateInterpolation - Gets interpolation matrix between two DMDA or DMComposite objects
616: Collective on DM
618: Input Parameter:
619: + dm1 - the DM object
620: - dm2 - the second, finer DM object
622: Output Parameter:
623: + mat - the interpolation
624: - vec - the scaling (optional)
626: Level: developer
628: 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
629: DMCoarsen(). The coordinates set into the DMDA are completely ignored in computing the interpolation.
631: For DMDA objects you can use this interpolation (more precisely the interpolation from the DMDAGetCoordinateDA()) to interpolate the mesh coordinate vectors
632: EXCEPT in the periodic case where it does not make sense since the coordinate vectors are not periodic.
633:
635: .seealso DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateColoring(), DMCreateMatrix(), DMRefine(), DMCoarsen()
637: @*/
638: PetscErrorCode DMCreateInterpolation(DM dm1,DM dm2,Mat *mat,Vec *vec)
639: {
645: (*dm1->ops->createinterpolation)(dm1,dm2,mat,vec);
646: return(0);
647: }
651: /*@
652: DMCreateInjection - Gets injection matrix between two DMDA or DMComposite objects
654: Collective on DM
656: Input Parameter:
657: + dm1 - the DM object
658: - dm2 - the second, finer DM object
660: Output Parameter:
661: . ctx - the injection
663: Level: developer
665: 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
666: DMCoarsen(). The coordinates set into the DMDA are completely ignored in computing the injection.
668: .seealso DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateColoring(), DMCreateMatrix(), DMCreateInterpolation()
670: @*/
671: PetscErrorCode DMCreateInjection(DM dm1,DM dm2,VecScatter *ctx)
672: {
678: (*dm1->ops->getinjection)(dm1,dm2,ctx);
679: return(0);
680: }
684: /*@C
685: DMCreateColoring - Gets coloring for a DMDA or DMComposite
687: Collective on DM
689: Input Parameter:
690: + dm - the DM object
691: . ctype - IS_COLORING_GHOSTED or IS_COLORING_GLOBAL
692: - matype - either MATAIJ or MATBAIJ
694: Output Parameter:
695: . coloring - the coloring
697: Level: developer
699: .seealso DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateMatrix()
701: @*/
702: PetscErrorCode DMCreateColoring(DM dm,ISColoringType ctype,const MatType mtype,ISColoring *coloring)
703: {
708: if (!dm->ops->getcoloring) SETERRQ(((PetscObject)dm)->comm,PETSC_ERR_SUP,"No coloring for this type of DM yet");
709: (*dm->ops->getcoloring)(dm,ctype,mtype,coloring);
710: return(0);
711: }
715: /*@C
716: DMCreateMatrix - Gets empty Jacobian for a DMDA or DMComposite
718: Collective on DM
720: Input Parameter:
721: + dm - the DM object
722: - mtype - Supported types are MATSEQAIJ, MATMPIAIJ, MATSEQBAIJ, MATMPIBAIJ, or
723: any type which inherits from one of these (such as MATAIJ)
725: Output Parameter:
726: . mat - the empty Jacobian
728: Level: beginner
730: Notes: This properly preallocates the number of nonzeros in the sparse matrix so you
731: do not need to do it yourself.
733: By default it also sets the nonzero structure and puts in the zero entries. To prevent setting
734: the nonzero pattern call DMDASetMatPreallocateOnly()
736: For structured grid problems, when you call MatView() on this matrix it is displayed using the global natural ordering, NOT in the ordering used
737: internally by PETSc.
739: For structured grid problems, in general it is easiest to use MatSetValuesStencil() or MatSetValuesLocal() to put values into the matrix because MatSetValues() requires
740: the indices for the global numbering for DMDAs which is complicated.
742: .seealso DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation()
744: @*/
745: PetscErrorCode DMCreateMatrix(DM dm,const MatType mtype,Mat *mat)
746: {
751: #ifndef PETSC_USE_DYNAMIC_LIBRARIES
752: MatInitializePackage(PETSC_NULL);
753: #endif
756: if (dm->mattype) {
757: (*dm->ops->creatematrix)(dm,dm->mattype,mat);
758: } else {
759: (*dm->ops->creatematrix)(dm,mtype,mat);
760: }
761: return(0);
762: }
766: /*@
767: DMSetMatrixPreallocateOnly - When DMCreateMatrix() is called the matrix will be properly
768: preallocated but the nonzero structure and zero values will not be set.
770: Logically Collective on DMDA
772: Input Parameter:
773: + dm - the DM
774: - only - PETSC_TRUE if only want preallocation
776: Level: developer
777: .seealso DMCreateMatrix()
778: @*/
779: PetscErrorCode DMSetMatrixPreallocateOnly(DM dm, PetscBool only)
780: {
783: dm->prealloc_only = only;
784: return(0);
785: }
789: /*@C
790: DMGetWorkArray - Gets a work array guaranteed to be at least the input size
792: Not Collective
794: Input Parameters:
795: + dm - the DM object
796: - size - The minium size
798: Output Parameter:
799: . array - the work array
801: Level: developer
803: .seealso DMDestroy(), DMCreate()
804: @*/
805: PetscErrorCode DMGetWorkArray(DM dm,PetscInt size,PetscScalar **array)
806: {
812: if (size > dm->workSize) {
813: dm->workSize = size;
814: PetscFree(dm->workArray);
815: PetscMalloc(dm->workSize * sizeof(PetscScalar), &dm->workArray);
816: }
817: *array = dm->workArray;
818: return(0);
819: }
824: /*@C
825: DMCreateFieldIS - Creates a set of IS objects with the global indices of dofs for each field
827: Not collective
829: Input Parameter:
830: . dm - the DM object
832: Output Parameters:
833: + numFields - The number of fields (or PETSC_NULL if not requested)
834: . fieldNames - The name for each field (or PETSC_NULL if not requested)
835: - fields - The global indices for each field (or PETSC_NULL if not requested)
837: Level: intermediate
839: Notes:
840: The user is responsible for freeing all requested arrays. In particular, every entry of names should be freed with
841: PetscFree(), every entry of fields should be destroyed with ISDestroy(), and both arrays should be freed with
842: PetscFree().
844: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix()
845: @*/
846: PetscErrorCode DMCreateFieldIS(DM dm, PetscInt *numFields, char ***fieldNames, IS **fields)
847: {
848: PetscSection section, sectionGlobal;
853: if (numFields) {
855: *numFields = 0;
856: }
857: if (fieldNames) {
859: *fieldNames = PETSC_NULL;
860: }
861: if (fields) {
863: *fields = PETSC_NULL;
864: }
865: DMGetDefaultSection(dm, §ion);
866: if (section) {
867: PetscInt *fieldSizes, **fieldIndices;
868: PetscInt nF, f, pStart, pEnd, p;
870: DMGetDefaultGlobalSection(dm, §ionGlobal);
871: PetscSectionGetNumFields(section, &nF);
872: PetscMalloc2(nF,PetscInt,&fieldSizes,nF,PetscInt *,&fieldIndices);
873: PetscSectionGetChart(sectionGlobal, &pStart, &pEnd);
874: for(f = 0; f < nF; ++f) {
875: fieldSizes[f] = 0;
876: }
877: for(p = pStart; p < pEnd; ++p) {
878: PetscInt gdof;
880: PetscSectionGetDof(sectionGlobal, p, &gdof);
881: if (gdof > 0) {
882: for(f = 0; f < nF; ++f) {
883: PetscInt fdof, fcdof;
885: PetscSectionGetFieldDof(section, p, f, &fdof);
886: PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);
887: fieldSizes[f] += fdof-fcdof;
888: }
889: }
890: }
891: for(f = 0; f < nF; ++f) {
892: PetscMalloc(fieldSizes[f] * sizeof(PetscInt), &fieldIndices[f]);
893: fieldSizes[f] = 0;
894: }
895: for(p = pStart; p < pEnd; ++p) {
896: PetscInt gdof, goff;
898: PetscSectionGetDof(sectionGlobal, p, &gdof);
899: if (gdof > 0) {
900: PetscSectionGetOffset(sectionGlobal, p, &goff);
901: for(f = 0; f < nF; ++f) {
902: PetscInt fdof, fcdof, fc;
904: PetscSectionGetFieldDof(section, p, f, &fdof);
905: PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);
906: for(fc = 0; fc < fdof-fcdof; ++fc, ++fieldSizes[f]) {
907: fieldIndices[f][fieldSizes[f]] = goff++;
908: }
909: }
910: }
911: }
912: if (numFields) {*numFields = nF;}
913: if (fieldNames) {
914: PetscMalloc(nF * sizeof(char *), fieldNames);
915: for(f = 0; f < nF; ++f) {
916: const char *fieldName;
918: PetscSectionGetFieldName(section, f, &fieldName);
919: PetscStrallocpy(fieldName, (char **) &(*fieldNames)[f]);
920: }
921: }
922: if (fields) {
923: PetscMalloc(nF * sizeof(IS), fields);
924: for(f = 0; f < nF; ++f) {
925: ISCreateGeneral(((PetscObject) dm)->comm, fieldSizes[f], fieldIndices[f], PETSC_OWN_POINTER, &(*fields)[f]);
926: }
927: }
928: PetscFree2(fieldSizes,fieldIndices);
929: } else {
930: if(dm->ops->createfieldis) {(*dm->ops->createfieldis)(dm, numFields, fieldNames, fields);}
931: }
932: return(0);
933: }
938: /*@C
939: DMCreateFieldDecompositionDM - creates a DM that encapsulates a decomposition of the original DM into fields.
941: Not Collective
943: Input Parameters:
944: + dm - the DM object
945: - name - the name of the field decomposition
947: Output Parameter:
948: . ddm - the field decomposition DM (PETSC_NULL, if no such decomposition is known)
950: Level: advanced
952: .seealso DMDestroy(), DMCreate(), DMCreateFieldDecomposition(), DMCreateDomainDecompositionDM()
953: @*/
954: PetscErrorCode DMCreateFieldDecompositionDM(DM dm, const char* name, DM *ddm)
955: {
962: *ddm = PETSC_NULL;
963: if(dm->ops->createfielddecompositiondm) {
964: (*dm->ops->createfielddecompositiondm)(dm,name,ddm);
965: }
966: return(0);
967: }
972: /*@C
973: DMCreateFieldDecomposition - Returns a list of IS objects defining a decomposition of a problem into subproblems
974: corresponding to different fields: each IS contains the global indices of the dofs of the
975: corresponding field. The optional list of DMs define the DM for each subproblem.
976: Generalizes DMCreateFieldIS().
978: Not collective
980: Input Parameter:
981: . dm - the DM object
983: Output Parameters:
984: + len - The number of subproblems in the field decomposition (or PETSC_NULL if not requested)
985: . namelist - The name for each field (or PETSC_NULL if not requested)
986: . islist - The global indices for each field (or PETSC_NULL if not requested)
987: - dmlist - The DMs for each field subproblem (or PETSC_NULL, if not requested; if PETSC_NULL is returned, no DMs are defined)
989: Level: intermediate
991: Notes:
992: The user is responsible for freeing all requested arrays. In particular, every entry of names should be freed with
993: PetscFree(), every entry of is should be destroyed with ISDestroy(), every entry of dm should be destroyed with DMDestroy(),
994: and all of the arrays should be freed with PetscFree().
996: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateFieldIS()
997: @*/
998: PetscErrorCode DMCreateFieldDecomposition(DM dm, PetscInt *len, char ***namelist, IS **islist, DM **dmlist)
999: {
1008: if(!dm->ops->createfielddecomposition) {
1009: DMCreateFieldIS(dm, len, namelist, islist);
1010: /* By default there are no DMs associated with subproblems. */
1011: if(dmlist) *dmlist = PETSC_NULL;
1012: }
1013: else {
1014: (*dm->ops->createfielddecomposition)(dm,len,namelist,islist,dmlist);
1015: }
1016: return(0);
1017: }
1021: /*@C
1022: DMCreateDomainDecompositionDM - creates a DM that encapsulates a decomposition of the original DM into subdomains.
1024: Not Collective
1026: Input Parameters:
1027: + dm - the DM object
1028: - name - the name of the subdomain decomposition
1030: Output Parameter:
1031: . ddm - the subdomain decomposition DM (PETSC_NULL, if no such decomposition is known)
1033: Level: advanced
1035: .seealso DMDestroy(), DMCreate(), DMCreateFieldDecomposition(), DMCreateDomainDecompositionDM()
1036: @*/
1037: PetscErrorCode DMCreateDomainDecompositionDM(DM dm, const char* name, DM *ddm)
1038: {
1045: *ddm = PETSC_NULL;
1046: if(dm->ops->createdomaindecompositiondm) {
1047: (*dm->ops->createdomaindecompositiondm)(dm,name,ddm);
1048: }
1049: return(0);
1050: }
1055: /*@C
1056: DMCreateDomainDecomposition - Returns lists of IS objects defining a decomposition of a problem into subproblems
1057: corresponding to restrictions to pairs nested subdomains: each IS contains the global
1058: indices of the dofs of the corresponding subdomains. The inner subdomains conceptually
1059: define a nonoverlapping covering, while outer subdomains can overlap.
1060: The optional list of DMs define the DM for each subproblem.
1062: Not collective
1064: Input Parameter:
1065: . dm - the DM object
1067: Output Parameters:
1068: + len - The number of subproblems in the domain decomposition (or PETSC_NULL if not requested)
1069: . namelist - The name for each subdomain (or PETSC_NULL if not requested)
1070: . innerislist - The global indices for each inner subdomain (or PETSC_NULL, if not requested)
1071: . outerislist - The global indices for each outer subdomain (or PETSC_NULL, if not requested)
1072: - dmlist - The DMs for each subdomain subproblem (or PETSC_NULL, if not requested; if PETSC_NULL is returned, no DMs are defined)
1074: Level: intermediate
1076: Notes:
1077: The user is responsible for freeing all requested arrays. In particular, every entry of names should be freed with
1078: PetscFree(), every entry of is should be destroyed with ISDestroy(), every entry of dm should be destroyed with DMDestroy(),
1079: and all of the arrays should be freed with PetscFree().
1081: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateDomainDecompositionDM(), DMCreateFieldDecomposition()
1082: @*/
1083: PetscErrorCode DMCreateDomainDecomposition(DM dm, PetscInt *len, char ***namelist, IS **innerislist, IS **outerislist, DM **dmlist)
1084: {
1094: if(dm->ops->createdomaindecomposition) {
1095: (*dm->ops->createdomaindecomposition)(dm,len,namelist,innerislist,outerislist,dmlist);
1096: }
1097: return(0);
1098: }
1102: /*@
1103: DMRefine - Refines a DM object
1105: Collective on DM
1107: Input Parameter:
1108: + dm - the DM object
1109: - comm - the communicator to contain the new DM object (or MPI_COMM_NULL)
1111: Output Parameter:
1112: . dmf - the refined DM, or PETSC_NULL
1114: Note: If no refinement was done, the return value is PETSC_NULL
1116: Level: developer
1118: .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation()
1119: @*/
1120: PetscErrorCode DMRefine(DM dm,MPI_Comm comm,DM *dmf)
1121: {
1123: DMRefineHookLink link;
1127: (*dm->ops->refine)(dm,comm,dmf);
1128: if (*dmf) {
1129: (*dmf)->ops->creatematrix = dm->ops->creatematrix;
1130: (*dmf)->ops->initialguess = dm->ops->initialguess;
1131: (*dmf)->ops->function = dm->ops->function;
1132: (*dmf)->ops->functionj = dm->ops->functionj;
1133: if (dm->ops->jacobian != DMComputeJacobianDefault) {
1134: (*dmf)->ops->jacobian = dm->ops->jacobian;
1135: }
1136: PetscObjectCopyFortranFunctionPointers((PetscObject)dm,(PetscObject)*dmf);
1137: (*dmf)->ctx = dm->ctx;
1138: (*dmf)->leveldown = dm->leveldown;
1139: (*dmf)->levelup = dm->levelup + 1;
1140: DMSetMatType(*dmf,dm->mattype);
1141: for (link=dm->refinehook; link; link=link->next) {
1142: if (link->refinehook) {(*link->refinehook)(dm,*dmf,link->ctx);}
1143: }
1144: }
1145: return(0);
1146: }
1150: /*@
1151: DMRefineHookAdd - adds a callback to be run when interpolating a nonlinear problem to a finer grid
1153: Logically Collective
1155: Input Arguments:
1156: + coarse - nonlinear solver context on which to run a hook when restricting to a coarser level
1157: . refinehook - function to run when setting up a coarser level
1158: . interphook - function to run to update data on finer levels (once per SNESSolve())
1159: - ctx - [optional] user-defined context for provide data for the hooks (may be PETSC_NULL)
1161: Calling sequence of refinehook:
1162: $ refinehook(DM coarse,DM fine,void *ctx);
1164: + coarse - coarse level DM
1165: . fine - fine level DM to interpolate problem to
1166: - ctx - optional user-defined function context
1168: Calling sequence for interphook:
1169: $ interphook(DM coarse,Mat interp,DM fine,void *ctx)
1171: + coarse - coarse level DM
1172: . interp - matrix interpolating a coarse-level solution to the finer grid
1173: . fine - fine level DM to update
1174: - ctx - optional user-defined function context
1176: Level: advanced
1178: Notes:
1179: This function is only needed if auxiliary data needs to be passed to fine grids while grid sequencing
1181: If this function is called multiple times, the hooks will be run in the order they are added.
1183: .seealso: DMCoarsenHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
1184: @*/
1185: PetscErrorCode DMRefineHookAdd(DM coarse,PetscErrorCode (*refinehook)(DM,DM,void*),PetscErrorCode (*interphook)(DM,Mat,DM,void*),void *ctx)
1186: {
1188: DMRefineHookLink link,*p;
1192: for (p=&coarse->refinehook; *p; p=&(*p)->next) {} /* Scan to the end of the current list of hooks */
1193: PetscMalloc(sizeof(struct _DMRefineHookLink),&link);
1194: link->refinehook = refinehook;
1195: link->interphook = interphook;
1196: link->ctx = ctx;
1197: link->next = PETSC_NULL;
1198: *p = link;
1199: return(0);
1200: }
1204: /*@
1205: DMInterpolate - interpolates user-defined problem data to a finer DM by running hooks registered by DMRefineHookAdd()
1207: Collective if any hooks are
1209: Input Arguments:
1210: + coarse - coarser DM to use as a base
1211: . restrct - interpolation matrix, apply using MatInterpolate()
1212: - fine - finer DM to update
1214: Level: developer
1216: .seealso: DMRefineHookAdd(), MatInterpolate()
1217: @*/
1218: PetscErrorCode DMInterpolate(DM coarse,Mat interp,DM fine)
1219: {
1221: DMRefineHookLink link;
1224: for (link=fine->refinehook; link; link=link->next) {
1225: if (link->interphook) {(*link->interphook)(coarse,interp,fine,link->ctx);}
1226: }
1227: return(0);
1228: }
1232: /*@
1233: DMGetRefineLevel - Get's the number of refinements that have generated this DM.
1235: Not Collective
1237: Input Parameter:
1238: . dm - the DM object
1240: Output Parameter:
1241: . level - number of refinements
1243: Level: developer
1245: .seealso DMCoarsen(), DMGetCoarsenLevel(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation()
1247: @*/
1248: PetscErrorCode DMGetRefineLevel(DM dm,PetscInt *level)
1249: {
1252: *level = dm->levelup;
1253: return(0);
1254: }
1258: /*@
1259: DMGlobalToLocalBegin - Begins updating local vectors from global vector
1261: Neighbor-wise Collective on DM
1263: Input Parameters:
1264: + dm - the DM object
1265: . g - the global vector
1266: . mode - INSERT_VALUES or ADD_VALUES
1267: - l - the local vector
1270: Level: beginner
1272: .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMGlobalToLocalEnd(), DMLocalToGlobalBegin()
1274: @*/
1275: PetscErrorCode DMGlobalToLocalBegin(DM dm,Vec g,InsertMode mode,Vec l)
1276: {
1277: PetscSF sf;
1282: DMGetDefaultSF(dm, &sf);
1283: if (sf) {
1284: PetscScalar *lArray, *gArray;
1286: if (mode == ADD_VALUES) SETERRQ1(((PetscObject) dm)->comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode);
1287: VecGetArray(l, &lArray);
1288: VecGetArray(g, &gArray);
1289: PetscSFBcastBegin(sf, MPIU_SCALAR, gArray, lArray);
1290: VecRestoreArray(l, &lArray);
1291: VecRestoreArray(g, &gArray);
1292: } else {
1293: (*dm->ops->globaltolocalbegin)(dm,g,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),l);
1294: }
1295: return(0);
1296: }
1300: /*@
1301: DMGlobalToLocalEnd - Ends updating local vectors from global vector
1303: Neighbor-wise Collective on DM
1305: Input Parameters:
1306: + dm - the DM object
1307: . g - the global vector
1308: . mode - INSERT_VALUES or ADD_VALUES
1309: - l - the local vector
1312: Level: beginner
1314: .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMGlobalToLocalEnd(), DMLocalToGlobalBegin()
1316: @*/
1317: PetscErrorCode DMGlobalToLocalEnd(DM dm,Vec g,InsertMode mode,Vec l)
1318: {
1319: PetscSF sf;
1321: PetscScalar *lArray, *gArray;
1325: DMGetDefaultSF(dm, &sf);
1326: if (sf) {
1327: if (mode == ADD_VALUES) SETERRQ1(((PetscObject) dm)->comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode);
1329: VecGetArray(l, &lArray);
1330: VecGetArray(g, &gArray);
1331: PetscSFBcastEnd(sf, MPIU_SCALAR, gArray, lArray);
1332: VecRestoreArray(l, &lArray);
1333: VecRestoreArray(g, &gArray);
1334: } else {
1335: (*dm->ops->globaltolocalend)(dm,g,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),l);
1336: }
1337: return(0);
1338: }
1342: /*@
1343: DMLocalToGlobalBegin - updates global vectors from local vectors
1345: Neighbor-wise Collective on DM
1347: Input Parameters:
1348: + dm - the DM object
1349: . l - the local vector
1350: . 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
1351: base point.
1352: - - the global vector
1354: 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
1355: 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
1356: global array to the final global array with VecAXPY().
1358: Level: beginner
1360: .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMGlobalToLocalEnd(), DMGlobalToLocalBegin()
1362: @*/
1363: PetscErrorCode DMLocalToGlobalBegin(DM dm,Vec l,InsertMode mode,Vec g)
1364: {
1365: PetscSF sf;
1370: DMGetDefaultSF(dm, &sf);
1371: if (sf) {
1372: MPI_Op op;
1373: PetscScalar *lArray, *gArray;
1375: switch(mode) {
1376: case INSERT_VALUES:
1377: case INSERT_ALL_VALUES:
1378: #if defined(PETSC_HAVE_MPI_REPLACE)
1379: op = MPI_REPLACE; break;
1380: #else
1381: SETERRQ(((PetscObject)dm)->comm,PETSC_ERR_SUP,"No support for INSERT_VALUES without an MPI-2 implementation");
1382: #endif
1383: case ADD_VALUES:
1384: case ADD_ALL_VALUES:
1385: op = MPI_SUM; break;
1386: default:
1387: SETERRQ1(((PetscObject) dm)->comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode);
1388: }
1389: VecGetArray(l, &lArray);
1390: VecGetArray(g, &gArray);
1391: PetscSFReduceBegin(sf, MPIU_SCALAR, lArray, gArray, op);
1392: VecRestoreArray(l, &lArray);
1393: VecRestoreArray(g, &gArray);
1394: } else {
1395: (*dm->ops->localtoglobalbegin)(dm,l,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),g);
1396: }
1397: return(0);
1398: }
1402: /*@
1403: DMLocalToGlobalEnd - updates global vectors from local vectors
1405: Neighbor-wise Collective on DM
1407: Input Parameters:
1408: + dm - the DM object
1409: . l - the local vector
1410: . mode - INSERT_VALUES or ADD_VALUES
1411: - g - the global vector
1414: Level: beginner
1416: .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMGlobalToLocalEnd(), DMGlobalToLocalEnd()
1418: @*/
1419: PetscErrorCode DMLocalToGlobalEnd(DM dm,Vec l,InsertMode mode,Vec g)
1420: {
1421: PetscSF sf;
1426: DMGetDefaultSF(dm, &sf);
1427: if (sf) {
1428: MPI_Op op;
1429: PetscScalar *lArray, *gArray;
1431: switch(mode) {
1432: case INSERT_VALUES:
1433: case INSERT_ALL_VALUES:
1434: #if defined(PETSC_HAVE_MPI_REPLACE)
1435: op = MPI_REPLACE; break;
1436: #else
1437: SETERRQ(((PetscObject)dm)->comm,PETSC_ERR_SUP,"No support for INSERT_VALUES without an MPI-2 implementation");
1438: #endif
1439: case ADD_VALUES:
1440: case ADD_ALL_VALUES:
1441: op = MPI_SUM; break;
1442: default:
1443: SETERRQ1(((PetscObject) dm)->comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode);
1444: }
1445: VecGetArray(l, &lArray);
1446: VecGetArray(g, &gArray);
1447: PetscSFReduceEnd(sf, MPIU_SCALAR, lArray, gArray, op);
1448: VecRestoreArray(l, &lArray);
1449: VecRestoreArray(g, &gArray);
1450: } else {
1451: (*dm->ops->localtoglobalend)(dm,l,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),g);
1452: }
1453: return(0);
1454: }
1458: /*@
1459: DMComputeJacobianDefault - computes the Jacobian using the DMComputeFunction() if Jacobian computer is not provided
1461: Collective on DM
1463: Input Parameter:
1464: + dm - the DM object
1465: . x - location to compute Jacobian at; may be ignored for linear problems
1466: . A - matrix that defines the operator for the linear solve
1467: - B - the matrix used to construct the preconditioner
1469: Level: developer
1471: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext(), DMSetInitialGuess(),
1472: DMSetFunction()
1474: @*/
1475: PetscErrorCode DMComputeJacobianDefault(DM dm,Vec x,Mat A,Mat B,MatStructure *stflag)
1476: {
1481: *stflag = SAME_NONZERO_PATTERN;
1482: MatFDColoringApply(B,dm->fd,x,stflag,dm);
1483: if (A != B) {
1484: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
1485: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
1486: }
1487: return(0);
1488: }
1492: /*@
1493: DMCoarsen - Coarsens a DM object
1495: Collective on DM
1497: Input Parameter:
1498: + dm - the DM object
1499: - comm - the communicator to contain the new DM object (or MPI_COMM_NULL)
1501: Output Parameter:
1502: . dmc - the coarsened DM
1504: Level: developer
1506: .seealso DMRefine(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation()
1508: @*/
1509: PetscErrorCode DMCoarsen(DM dm, MPI_Comm comm, DM *dmc)
1510: {
1512: DMCoarsenHookLink link;
1516: (*dm->ops->coarsen)(dm, comm, dmc);
1517: (*dmc)->ops->creatematrix = dm->ops->creatematrix;
1518: (*dmc)->ops->initialguess = dm->ops->initialguess;
1519: (*dmc)->ops->function = dm->ops->function;
1520: (*dmc)->ops->functionj = dm->ops->functionj;
1521: if (dm->ops->jacobian != DMComputeJacobianDefault) {
1522: (*dmc)->ops->jacobian = dm->ops->jacobian;
1523: }
1524: PetscObjectCopyFortranFunctionPointers((PetscObject)dm,(PetscObject)*dmc);
1525: (*dmc)->ctx = dm->ctx;
1526: (*dmc)->levelup = dm->levelup;
1527: (*dmc)->leveldown = dm->leveldown + 1;
1528: DMSetMatType(*dmc,dm->mattype);
1529: for (link=dm->coarsenhook; link; link=link->next) {
1530: if (link->coarsenhook) {(*link->coarsenhook)(dm,*dmc,link->ctx);}
1531: }
1532: return(0);
1533: }
1537: /*@
1538: DMCoarsenHookAdd - adds a callback to be run when restricting a nonlinear problem to the coarse grid
1540: Logically Collective
1542: Input Arguments:
1543: + fine - nonlinear solver context on which to run a hook when restricting to a coarser level
1544: . coarsenhook - function to run when setting up a coarser level
1545: . restricthook - function to run to update data on coarser levels (once per SNESSolve())
1546: - ctx - [optional] user-defined context for provide data for the hooks (may be PETSC_NULL)
1548: Calling sequence of coarsenhook:
1549: $ coarsenhook(DM fine,DM coarse,void *ctx);
1551: + fine - fine level DM
1552: . coarse - coarse level DM to restrict problem to
1553: - ctx - optional user-defined function context
1555: Calling sequence for restricthook:
1556: $ restricthook(DM fine,Mat mrestrict,Vec rscale,Mat inject,DM coarse,void *ctx)
1558: + fine - fine level DM
1559: . mrestrict - matrix restricting a fine-level solution to the coarse grid
1560: . rscale - scaling vector for restriction
1561: . inject - matrix restricting by injection
1562: . coarse - coarse level DM to update
1563: - ctx - optional user-defined function context
1565: Level: advanced
1567: Notes:
1568: This function is only needed if auxiliary data needs to be set up on coarse grids.
1570: If this function is called multiple times, the hooks will be run in the order they are added.
1572: In order to compose with nonlinear preconditioning without duplicating storage, the hook should be implemented to
1573: extract the finest level information from its context (instead of from the SNES).
1575: .seealso: DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
1576: @*/
1577: PetscErrorCode DMCoarsenHookAdd(DM fine,PetscErrorCode (*coarsenhook)(DM,DM,void*),PetscErrorCode (*restricthook)(DM,Mat,Vec,Mat,DM,void*),void *ctx)
1578: {
1580: DMCoarsenHookLink link,*p;
1584: for (p=&fine->coarsenhook; *p; p=&(*p)->next) {} /* Scan to the end of the current list of hooks */
1585: PetscMalloc(sizeof(struct _DMCoarsenHookLink),&link);
1586: link->coarsenhook = coarsenhook;
1587: link->restricthook = restricthook;
1588: link->ctx = ctx;
1589: link->next = PETSC_NULL;
1590: *p = link;
1591: return(0);
1592: }
1596: /*@
1597: DMRestrict - restricts user-defined problem data to a coarser DM by running hooks registered by DMCoarsenHookAdd()
1599: Collective if any hooks are
1601: Input Arguments:
1602: + fine - finer DM to use as a base
1603: . restrct - restriction matrix, apply using MatRestrict()
1604: . inject - injection matrix, also use MatRestrict()
1605: - coarse - coarer DM to update
1607: Level: developer
1609: .seealso: DMCoarsenHookAdd(), MatRestrict()
1610: @*/
1611: PetscErrorCode DMRestrict(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse)
1612: {
1614: DMCoarsenHookLink link;
1617: for (link=fine->coarsenhook; link; link=link->next) {
1618: if (link->restricthook) {(*link->restricthook)(fine,restrct,rscale,inject,coarse,link->ctx);}
1619: }
1620: return(0);
1621: }
1625: /*@
1626: DMGetCoarsenLevel - Get's the number of coarsenings that have generated this DM.
1628: Not Collective
1630: Input Parameter:
1631: . dm - the DM object
1633: Output Parameter:
1634: . level - number of coarsenings
1636: Level: developer
1638: .seealso DMCoarsen(), DMGetRefineLevel(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation()
1640: @*/
1641: PetscErrorCode DMGetCoarsenLevel(DM dm,PetscInt *level)
1642: {
1645: *level = dm->leveldown;
1646: return(0);
1647: }
1653: /*@C
1654: DMRefineHierarchy - Refines a DM object, all levels at once
1656: Collective on DM
1658: Input Parameter:
1659: + dm - the DM object
1660: - nlevels - the number of levels of refinement
1662: Output Parameter:
1663: . dmf - the refined DM hierarchy
1665: Level: developer
1667: .seealso DMCoarsenHierarchy(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation()
1669: @*/
1670: PetscErrorCode DMRefineHierarchy(DM dm,PetscInt nlevels,DM dmf[])
1671: {
1676: if (nlevels < 0) SETERRQ(((PetscObject)dm)->comm,PETSC_ERR_ARG_OUTOFRANGE,"nlevels cannot be negative");
1677: if (nlevels == 0) return(0);
1678: if (dm->ops->refinehierarchy) {
1679: (*dm->ops->refinehierarchy)(dm,nlevels,dmf);
1680: } else if (dm->ops->refine) {
1681: PetscInt i;
1683: DMRefine(dm,((PetscObject)dm)->comm,&dmf[0]);
1684: for (i=1; i<nlevels; i++) {
1685: DMRefine(dmf[i-1],((PetscObject)dm)->comm,&dmf[i]);
1686: }
1687: } else {
1688: SETERRQ(((PetscObject)dm)->comm,PETSC_ERR_SUP,"No RefineHierarchy for this DM yet");
1689: }
1690: return(0);
1691: }
1695: /*@C
1696: DMCoarsenHierarchy - Coarsens a DM object, all levels at once
1698: Collective on DM
1700: Input Parameter:
1701: + dm - the DM object
1702: - nlevels - the number of levels of coarsening
1704: Output Parameter:
1705: . dmc - the coarsened DM hierarchy
1707: Level: developer
1709: .seealso DMRefineHierarchy(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation()
1711: @*/
1712: PetscErrorCode DMCoarsenHierarchy(DM dm, PetscInt nlevels, DM dmc[])
1713: {
1718: if (nlevels < 0) SETERRQ(((PetscObject)dm)->comm,PETSC_ERR_ARG_OUTOFRANGE,"nlevels cannot be negative");
1719: if (nlevels == 0) return(0);
1721: if (dm->ops->coarsenhierarchy) {
1722: (*dm->ops->coarsenhierarchy)(dm, nlevels, dmc);
1723: } else if (dm->ops->coarsen) {
1724: PetscInt i;
1726: DMCoarsen(dm,((PetscObject)dm)->comm,&dmc[0]);
1727: for (i=1; i<nlevels; i++) {
1728: DMCoarsen(dmc[i-1],((PetscObject)dm)->comm,&dmc[i]);
1729: }
1730: } else {
1731: SETERRQ(((PetscObject)dm)->comm,PETSC_ERR_SUP,"No CoarsenHierarchy for this DM yet");
1732: }
1733: return(0);
1734: }
1738: /*@
1739: DMCreateAggregates - Gets the aggregates that map between
1740: grids associated with two DMs.
1742: Collective on DM
1744: Input Parameters:
1745: + dmc - the coarse grid DM
1746: - dmf - the fine grid DM
1748: Output Parameters:
1749: . rest - the restriction matrix (transpose of the projection matrix)
1751: Level: intermediate
1753: .keywords: interpolation, restriction, multigrid
1755: .seealso: DMRefine(), DMCreateInjection(), DMCreateInterpolation()
1756: @*/
1757: PetscErrorCode DMCreateAggregates(DM dmc, DM dmf, Mat *rest)
1758: {
1764: (*dmc->ops->getaggregates)(dmc, dmf, rest);
1765: return(0);
1766: }
1770: /*@C
1771: DMSetApplicationContextDestroy - Sets a user function that will be called to destroy the application context when the DM is destroyed
1773: Not Collective
1775: Input Parameters:
1776: + dm - the DM object
1777: - destroy - the destroy function
1779: Level: intermediate
1781: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext()
1783: @*/
1784: PetscErrorCode DMSetApplicationContextDestroy(DM dm,PetscErrorCode (*destroy)(void**))
1785: {
1788: dm->ctxdestroy = destroy;
1789: return(0);
1790: }
1794: /*@
1795: DMSetApplicationContext - Set a user context into a DM object
1797: Not Collective
1799: Input Parameters:
1800: + dm - the DM object
1801: - ctx - the user context
1803: Level: intermediate
1805: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext()
1807: @*/
1808: PetscErrorCode DMSetApplicationContext(DM dm,void *ctx)
1809: {
1812: dm->ctx = ctx;
1813: return(0);
1814: }
1818: /*@
1819: DMGetApplicationContext - Gets a user context from a DM object
1821: Not Collective
1823: Input Parameter:
1824: . dm - the DM object
1826: Output Parameter:
1827: . ctx - the user context
1829: Level: intermediate
1831: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext()
1833: @*/
1834: PetscErrorCode DMGetApplicationContext(DM dm,void *ctx)
1835: {
1838: *(void**)ctx = dm->ctx;
1839: return(0);
1840: }
1844: /*@C
1845: DMSetInitialGuess - sets a function to compute an initial guess vector entries for the solvers
1847: Logically Collective on DM
1849: Input Parameter:
1850: + dm - the DM object to destroy
1851: - f - the function to compute the initial guess
1853: Level: intermediate
1855: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext(), DMSetFunction(), DMSetJacobian()
1857: @*/
1858: PetscErrorCode DMSetInitialGuess(DM dm,PetscErrorCode (*f)(DM,Vec))
1859: {
1862: dm->ops->initialguess = f;
1863: return(0);
1864: }
1868: /*@C
1869: DMSetFunction - sets a function to compute the right hand side vector entries for the KSP solver or nonlinear function for SNES
1871: Logically Collective on DM
1873: Input Parameter:
1874: + dm - the DM object
1875: - f - the function to compute (use PETSC_NULL to cancel a previous function that was set)
1877: Level: intermediate
1879: Notes: This sets both the function for function evaluations and the function used to compute Jacobians via finite differences if no Jacobian
1880: computer is provided with DMSetJacobian(). Canceling cancels the function, but not the function used to compute the Jacobian.
1882: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext(), DMSetInitialGuess(),
1883: DMSetJacobian()
1885: @*/
1886: PetscErrorCode DMSetFunction(DM dm,PetscErrorCode (*f)(DM,Vec,Vec))
1887: {
1890: dm->ops->function = f;
1891: if (f) {
1892: dm->ops->functionj = f;
1893: }
1894: return(0);
1895: }
1899: /*@C
1900: DMSetJacobian - sets a function to compute the matrix entries for the KSP solver or Jacobian for SNES
1902: Logically Collective on DM
1904: Input Parameter:
1905: + dm - the DM object to destroy
1906: - f - the function to compute the matrix entries
1908: Level: intermediate
1910: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext(), DMSetInitialGuess(),
1911: DMSetFunction()
1913: @*/
1914: PetscErrorCode DMSetJacobian(DM dm,PetscErrorCode (*f)(DM,Vec,Mat,Mat,MatStructure*))
1915: {
1918: dm->ops->jacobian = f;
1919: return(0);
1920: }
1924: /*@C
1925: DMSetVariableBounds - sets a function to compute the the lower and upper bound vectors for SNESVI.
1927: Logically Collective on DM
1929: Input Parameter:
1930: + dm - the DM object
1931: - f - the function that computes variable bounds used by SNESVI (use PETSC_NULL to cancel a previous function that was set)
1933: Level: intermediate
1935: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext(), DMSetInitialGuess(),
1936: DMSetJacobian()
1938: @*/
1939: PetscErrorCode DMSetVariableBounds(DM dm,PetscErrorCode (*f)(DM,Vec,Vec))
1940: {
1942: dm->ops->computevariablebounds = f;
1943: return(0);
1944: }
1948: /*@
1949: DMHasVariableBounds - does the DM object have a variable bounds function?
1951: Not Collective
1953: Input Parameter:
1954: . dm - the DM object to destroy
1956: Output Parameter:
1957: . flg - PETSC_TRUE if the variable bounds function exists
1959: Level: developer
1961: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext(), DMSetFunction(), DMSetJacobian()
1963: @*/
1964: PetscErrorCode DMHasVariableBounds(DM dm,PetscBool *flg)
1965: {
1967: *flg = (dm->ops->computevariablebounds) ? PETSC_TRUE : PETSC_FALSE;
1968: return(0);
1969: }
1973: /*@C
1974: DMComputeVariableBounds - compute variable bounds used by SNESVI.
1976: Logically Collective on DM
1978: Input Parameters:
1979: + dm - the DM object to destroy
1980: - x - current solution at which the bounds are computed
1982: Output parameters:
1983: + xl - lower bound
1984: - xu - upper bound
1986: Level: intermediate
1988: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext(), DMSetInitialGuess(),
1989: DMSetFunction(), DMSetVariableBounds()
1991: @*/
1992: PetscErrorCode DMComputeVariableBounds(DM dm, Vec xl, Vec xu)
1993: {
1998: if(dm->ops->computevariablebounds) {
1999: (*dm->ops->computevariablebounds)(dm, xl,xu);
2000: }
2001: else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "This DM is incapable of computing variable bounds.");
2002: return(0);
2003: }
2007: /*@
2008: DMComputeInitialGuess - computes an initial guess vector entries for the KSP solvers
2010: Collective on DM
2012: Input Parameter:
2013: + dm - the DM object to destroy
2014: - x - the vector to hold the initial guess values
2016: Level: developer
2018: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext(), DMSetRhs(), DMSetMat()
2020: @*/
2021: PetscErrorCode DMComputeInitialGuess(DM dm,Vec x)
2022: {
2026: if (!dm->ops->initialguess) SETERRQ(((PetscObject)dm)->comm,PETSC_ERR_ARG_WRONGSTATE,"Need to provide function with DMSetInitialGuess()");
2027: (*dm->ops->initialguess)(dm,x);
2028: return(0);
2029: }
2033: /*@
2034: DMHasInitialGuess - does the DM object have an initial guess function
2036: Not Collective
2038: Input Parameter:
2039: . dm - the DM object to destroy
2041: Output Parameter:
2042: . flg - PETSC_TRUE if function exists
2044: Level: developer
2046: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext(), DMSetFunction(), DMSetJacobian()
2048: @*/
2049: PetscErrorCode DMHasInitialGuess(DM dm,PetscBool *flg)
2050: {
2053: *flg = (dm->ops->initialguess) ? PETSC_TRUE : PETSC_FALSE;
2054: return(0);
2055: }
2059: /*@
2060: DMHasFunction - does the DM object have a function
2062: Not Collective
2064: Input Parameter:
2065: . dm - the DM object to destroy
2067: Output Parameter:
2068: . flg - PETSC_TRUE if function exists
2070: Level: developer
2072: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext(), DMSetFunction(), DMSetJacobian()
2074: @*/
2075: PetscErrorCode DMHasFunction(DM dm,PetscBool *flg)
2076: {
2079: *flg = (dm->ops->function) ? PETSC_TRUE : PETSC_FALSE;
2080: return(0);
2081: }
2085: /*@
2086: DMHasJacobian - does the DM object have a matrix function
2088: Not Collective
2090: Input Parameter:
2091: . dm - the DM object to destroy
2093: Output Parameter:
2094: . flg - PETSC_TRUE if function exists
2096: Level: developer
2098: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext(), DMSetFunction(), DMSetJacobian()
2100: @*/
2101: PetscErrorCode DMHasJacobian(DM dm,PetscBool *flg)
2102: {
2105: *flg = (dm->ops->jacobian) ? PETSC_TRUE : PETSC_FALSE;
2106: return(0);
2107: }
2109: #undef __FUNCT__
2111: /*@C
2112: DMSetVec - set the vector at which to compute residual, Jacobian and VI bounds, if the problem is nonlinear.
2114: Collective on DM
2116: Input Parameter:
2117: + dm - the DM object
2118: - x - location to compute residual and Jacobian, if PETSC_NULL is passed to those routines; will be PETSC_NULL for linear problems.
2120: Level: developer
2122: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext(), DMSetInitialGuess(),
2123: DMSetFunction(), DMSetJacobian(), DMSetVariableBounds()
2125: @*/
2126: PetscErrorCode DMSetVec(DM dm,Vec x)
2127: {
2130: if (x) {
2131: if (!dm->x) {
2132: DMCreateGlobalVector(dm,&dm->x);
2133: }
2134: VecCopy(x,dm->x);
2135: }
2136: else if(dm->x) {
2137: VecDestroy(&dm->x);
2138: }
2139: return(0);
2140: }
2145: /*@
2146: DMComputeFunction - computes the right hand side vector entries for the KSP solver or nonlinear function for SNES
2148: Collective on DM
2150: Input Parameter:
2151: + dm - the DM object to destroy
2152: . x - the location where the function is evaluationed, may be ignored for linear problems
2153: - b - the vector to hold the right hand side entries
2155: Level: developer
2157: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext(), DMSetInitialGuess(),
2158: DMSetJacobian()
2160: @*/
2161: PetscErrorCode DMComputeFunction(DM dm,Vec x,Vec b)
2162: {
2166: if (!dm->ops->function) SETERRQ(((PetscObject)dm)->comm,PETSC_ERR_ARG_WRONGSTATE,"Need to provide function with DMSetFunction()");
2167: PetscStackPush("DM user function");
2168: (*dm->ops->function)(dm,x,b);
2169: PetscStackPop;
2170: return(0);
2171: }
2177: /*@
2178: DMComputeJacobian - compute the matrix entries for the solver
2180: Collective on DM
2182: Input Parameter:
2183: + dm - the DM object
2184: . x - location to compute Jacobian at; will be PETSC_NULL for linear problems, for nonlinear problems if not provided then pulled from DM
2185: . A - matrix that defines the operator for the linear solve
2186: - B - the matrix used to construct the preconditioner
2188: Level: developer
2190: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext(), DMSetInitialGuess(),
2191: DMSetFunction()
2193: @*/
2194: PetscErrorCode DMComputeJacobian(DM dm,Vec x,Mat A,Mat B,MatStructure *stflag)
2195: {
2200: if (!dm->ops->jacobian) {
2201: ISColoring coloring;
2202: MatFDColoring fd;
2203: const MatType mtype;
2204:
2205: PetscObjectGetType((PetscObject)B,&mtype);
2206: DMCreateColoring(dm,dm->coloringtype,mtype,&coloring);
2207: MatFDColoringCreate(B,coloring,&fd);
2208: ISColoringDestroy(&coloring);
2209: MatFDColoringSetFunction(fd,(PetscErrorCode (*)(void))dm->ops->functionj,dm);
2210: PetscObjectSetOptionsPrefix((PetscObject)fd,((PetscObject)dm)->prefix);
2211: MatFDColoringSetFromOptions(fd);
2213: dm->fd = fd;
2214: dm->ops->jacobian = DMComputeJacobianDefault;
2216: /* don't know why this is needed */
2217: PetscObjectDereference((PetscObject)dm);
2218: }
2219: if (!x) x = dm->x;
2220: (*dm->ops->jacobian)(dm,x,A,B,stflag);
2222: /* if matrix depends on x; i.e. nonlinear problem, keep copy of input vector since needed by multigrid methods to generate coarse grid matrices */
2223: if (x) {
2224: if (!dm->x) {
2225: DMCreateGlobalVector(dm,&dm->x);
2226: }
2227: VecCopy(x,dm->x);
2228: }
2229: if (A != B) {
2230: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
2231: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
2232: }
2233: return(0);
2234: }
2237: PetscFList DMList = PETSC_NULL;
2238: PetscBool DMRegisterAllCalled = PETSC_FALSE;
2242: /*@C
2243: DMSetType - Builds a DM, for a particular DM implementation.
2245: Collective on DM
2247: Input Parameters:
2248: + dm - The DM object
2249: - method - The name of the DM type
2251: Options Database Key:
2252: . -dm_type <type> - Sets the DM type; use -help for a list of available types
2254: Notes:
2255: See "petsc/include/petscdm.h" for available DM types (for instance, DM1D, DM2D, or DM3D).
2257: Level: intermediate
2259: .keywords: DM, set, type
2260: .seealso: DMGetType(), DMCreate()
2261: @*/
2262: PetscErrorCode DMSetType(DM dm, const DMType method)
2263: {
2264: PetscErrorCode (*r)(DM);
2265: PetscBool match;
2270: PetscObjectTypeCompare((PetscObject) dm, method, &match);
2271: if (match) return(0);
2273: if (!DMRegisterAllCalled) {DMRegisterAll(PETSC_NULL);}
2274: PetscFListFind(DMList, ((PetscObject)dm)->comm, method,PETSC_TRUE,(void (**)(void)) &r);
2275: if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown DM type: %s", method);
2277: if (dm->ops->destroy) {
2278: (*dm->ops->destroy)(dm);
2279: dm->ops->destroy = PETSC_NULL;
2280: }
2281: (*r)(dm);
2282: PetscObjectChangeTypeName((PetscObject)dm,method);
2283: return(0);
2284: }
2288: /*@C
2289: DMGetType - Gets the DM type name (as a string) from the DM.
2291: Not Collective
2293: Input Parameter:
2294: . dm - The DM
2296: Output Parameter:
2297: . type - The DM type name
2299: Level: intermediate
2301: .keywords: DM, get, type, name
2302: .seealso: DMSetType(), DMCreate()
2303: @*/
2304: PetscErrorCode DMGetType(DM dm, const DMType *type)
2305: {
2311: if (!DMRegisterAllCalled) {
2312: DMRegisterAll(PETSC_NULL);
2313: }
2314: *type = ((PetscObject)dm)->type_name;
2315: return(0);
2316: }
2320: /*@C
2321: DMConvert - Converts a DM to another DM, either of the same or different type.
2323: Collective on DM
2325: Input Parameters:
2326: + dm - the DM
2327: - newtype - new DM type (use "same" for the same type)
2329: Output Parameter:
2330: . M - pointer to new DM
2332: Notes:
2333: Cannot be used to convert a sequential DM to parallel or parallel to sequential,
2334: the MPI communicator of the generated DM is always the same as the communicator
2335: of the input DM.
2337: Level: intermediate
2339: .seealso: DMCreate()
2340: @*/
2341: PetscErrorCode DMConvert(DM dm, const DMType newtype, DM *M)
2342: {
2343: DM B;
2344: char convname[256];
2345: PetscBool sametype, issame;
2352: PetscObjectTypeCompare((PetscObject) dm, newtype, &sametype);
2353: PetscStrcmp(newtype, "same", &issame);
2354: {
2355: PetscErrorCode (*conv)(DM, const DMType, DM *) = PETSC_NULL;
2357: /*
2358: Order of precedence:
2359: 1) See if a specialized converter is known to the current DM.
2360: 2) See if a specialized converter is known to the desired DM class.
2361: 3) See if a good general converter is registered for the desired class
2362: 4) See if a good general converter is known for the current matrix.
2363: 5) Use a really basic converter.
2364: */
2366: /* 1) See if a specialized converter is known to the current DM and the desired class */
2367: PetscStrcpy(convname,"DMConvert_");
2368: PetscStrcat(convname,((PetscObject) dm)->type_name);
2369: PetscStrcat(convname,"_");
2370: PetscStrcat(convname,newtype);
2371: PetscStrcat(convname,"_C");
2372: PetscObjectQueryFunction((PetscObject)dm,convname,(void (**)(void))&conv);
2373: if (conv) goto foundconv;
2375: /* 2) See if a specialized converter is known to the desired DM class. */
2376: DMCreate(((PetscObject) dm)->comm, &B);
2377: DMSetType(B, newtype);
2378: PetscStrcpy(convname,"DMConvert_");
2379: PetscStrcat(convname,((PetscObject) dm)->type_name);
2380: PetscStrcat(convname,"_");
2381: PetscStrcat(convname,newtype);
2382: PetscStrcat(convname,"_C");
2383: PetscObjectQueryFunction((PetscObject)B,convname,(void (**)(void))&conv);
2384: if (conv) {
2385: DMDestroy(&B);
2386: goto foundconv;
2387: }
2389: #if 0
2390: /* 3) See if a good general converter is registered for the desired class */
2391: conv = B->ops->convertfrom;
2392: DMDestroy(&B);
2393: if (conv) goto foundconv;
2395: /* 4) See if a good general converter is known for the current matrix */
2396: if (dm->ops->convert) {
2397: conv = dm->ops->convert;
2398: }
2399: if (conv) goto foundconv;
2400: #endif
2402: /* 5) Use a really basic converter. */
2403: SETERRQ2(((PetscObject) dm)->comm, PETSC_ERR_SUP, "No conversion possible between DM types %s and %s", ((PetscObject) dm)->type_name, newtype);
2405: foundconv:
2406: PetscLogEventBegin(DM_Convert,dm,0,0,0);
2407: (*conv)(dm,newtype,M);
2408: PetscLogEventEnd(DM_Convert,dm,0,0,0);
2409: }
2410: PetscObjectStateIncrease((PetscObject) *M);
2411: return(0);
2412: }
2414: /*--------------------------------------------------------------------------------------------------------------------*/
2418: /*@C
2419: DMRegister - See DMRegisterDynamic()
2421: Level: advanced
2422: @*/
2423: PetscErrorCode DMRegister(const char sname[], const char path[], const char name[], PetscErrorCode (*function)(DM))
2424: {
2425: char fullname[PETSC_MAX_PATH_LEN];
2429: PetscStrcpy(fullname, path);
2430: PetscStrcat(fullname, ":");
2431: PetscStrcat(fullname, name);
2432: PetscFListAdd(&DMList, sname, fullname, (void (*)(void)) function);
2433: return(0);
2434: }
2437: /*--------------------------------------------------------------------------------------------------------------------*/
2440: /*@C
2441: DMRegisterDestroy - Frees the list of DM methods that were registered by DMRegister()/DMRegisterDynamic().
2443: Not Collective
2445: Level: advanced
2447: .keywords: DM, register, destroy
2448: .seealso: DMRegister(), DMRegisterAll(), DMRegisterDynamic()
2449: @*/
2450: PetscErrorCode DMRegisterDestroy(void)
2451: {
2455: PetscFListDestroy(&DMList);
2456: DMRegisterAllCalled = PETSC_FALSE;
2457: return(0);
2458: }
2460: #if defined(PETSC_HAVE_MATLAB_ENGINE)
2461: #include <mex.h>
2463: typedef struct {char *funcname; char *jacname; mxArray *ctx;} DMMatlabContext;
2467: /*
2468: DMComputeFunction_Matlab - Calls the function that has been set with
2469: DMSetFunctionMatlab().
2471: For linear problems x is null
2472:
2473: .seealso: DMSetFunction(), DMGetFunction()
2474: */
2475: PetscErrorCode DMComputeFunction_Matlab(DM dm,Vec x,Vec y)
2476: {
2477: PetscErrorCode ierr;
2478: DMMatlabContext *sctx;
2479: int nlhs = 1,nrhs = 4;
2480: mxArray *plhs[1],*prhs[4];
2481: long long int lx = 0,ly = 0,ls = 0;
2482:
2488: /* call Matlab function in ctx with arguments x and y */
2489: DMGetApplicationContext(dm,&sctx);
2490: PetscMemcpy(&ls,&dm,sizeof(dm));
2491: PetscMemcpy(&lx,&x,sizeof(x));
2492: PetscMemcpy(&ly,&y,sizeof(y));
2493: prhs[0] = mxCreateDoubleScalar((double)ls);
2494: prhs[1] = mxCreateDoubleScalar((double)lx);
2495: prhs[2] = mxCreateDoubleScalar((double)ly);
2496: prhs[3] = mxCreateString(sctx->funcname);
2497: mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscDMComputeFunctionInternal");
2498: mxGetScalar(plhs[0]);
2499: mxDestroyArray(prhs[0]);
2500: mxDestroyArray(prhs[1]);
2501: mxDestroyArray(prhs[2]);
2502: mxDestroyArray(prhs[3]);
2503: mxDestroyArray(plhs[0]);
2504: return(0);
2505: }
2510: /*
2511: DMSetFunctionMatlab - Sets the function evaluation routine
2513: */
2514: PetscErrorCode DMSetFunctionMatlab(DM dm,const char *func)
2515: {
2516: PetscErrorCode ierr;
2517: DMMatlabContext *sctx;
2521: /* currently sctx is memory bleed */
2522: DMGetApplicationContext(dm,&sctx);
2523: if (!sctx) {
2524: PetscMalloc(sizeof(DMMatlabContext),&sctx);
2525: }
2526: PetscStrallocpy(func,&sctx->funcname);
2527: DMSetApplicationContext(dm,sctx);
2528: DMSetFunction(dm,DMComputeFunction_Matlab);
2529: return(0);
2530: }
2534: /*
2535: DMComputeJacobian_Matlab - Calls the function that has been set with
2536: DMSetJacobianMatlab().
2538: For linear problems x is null
2539:
2540: .seealso: DMSetFunction(), DMGetFunction()
2541: */
2542: PetscErrorCode DMComputeJacobian_Matlab(DM dm,Vec x,Mat A,Mat B,MatStructure *str)
2543: {
2544: PetscErrorCode ierr;
2545: DMMatlabContext *sctx;
2546: int nlhs = 2,nrhs = 5;
2547: mxArray *plhs[2],*prhs[5];
2548: long long int lx = 0,lA = 0,lB = 0,ls = 0;
2549:
2554: /* call MATLAB function in ctx with arguments x, A, and B */
2555: DMGetApplicationContext(dm,&sctx);
2556: PetscMemcpy(&ls,&dm,sizeof(dm));
2557: PetscMemcpy(&lx,&x,sizeof(x));
2558: PetscMemcpy(&lA,&A,sizeof(A));
2559: PetscMemcpy(&lB,&B,sizeof(B));
2560: prhs[0] = mxCreateDoubleScalar((double)ls);
2561: prhs[1] = mxCreateDoubleScalar((double)lx);
2562: prhs[2] = mxCreateDoubleScalar((double)lA);
2563: prhs[3] = mxCreateDoubleScalar((double)lB);
2564: prhs[4] = mxCreateString(sctx->jacname);
2565: mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscDMComputeJacobianInternal");
2566: *str = (MatStructure) mxGetScalar(plhs[0]);
2567: (PetscInt) mxGetScalar(plhs[1]);
2568: mxDestroyArray(prhs[0]);
2569: mxDestroyArray(prhs[1]);
2570: mxDestroyArray(prhs[2]);
2571: mxDestroyArray(prhs[3]);
2572: mxDestroyArray(prhs[4]);
2573: mxDestroyArray(plhs[0]);
2574: mxDestroyArray(plhs[1]);
2575: return(0);
2576: }
2581: /*
2582: DMSetJacobianMatlab - Sets the Jacobian function evaluation routine
2584: */
2585: PetscErrorCode DMSetJacobianMatlab(DM dm,const char *func)
2586: {
2587: PetscErrorCode ierr;
2588: DMMatlabContext *sctx;
2592: /* currently sctx is memory bleed */
2593: DMGetApplicationContext(dm,&sctx);
2594: if (!sctx) {
2595: PetscMalloc(sizeof(DMMatlabContext),&sctx);
2596: }
2597: PetscStrallocpy(func,&sctx->jacname);
2598: DMSetApplicationContext(dm,sctx);
2599: DMSetJacobian(dm,DMComputeJacobian_Matlab);
2600: return(0);
2601: }
2602: #endif
2606: /*@C
2607: DMLoad - Loads a DM that has been stored in binary or HDF5 format
2608: with DMView().
2610: Collective on PetscViewer
2612: Input Parameters:
2613: + newdm - the newly loaded DM, this needs to have been created with DMCreate() or
2614: some related function before a call to DMLoad().
2615: - viewer - binary file viewer, obtained from PetscViewerBinaryOpen() or
2616: HDF5 file viewer, obtained from PetscViewerHDF5Open()
2618: Level: intermediate
2620: Notes:
2621: Defaults to the DM DA.
2623: Notes for advanced users:
2624: Most users should not need to know the details of the binary storage
2625: format, since DMLoad() and DMView() completely hide these details.
2626: But for anyone who's interested, the standard binary matrix storage
2627: format is
2628: .vb
2629: has not yet been determined
2630: .ve
2632: In addition, PETSc automatically does the byte swapping for
2633: machines that store the bytes reversed, e.g. DEC alpha, freebsd,
2634: linux, Windows and the paragon; thus if you write your own binary
2635: read/write routines you have to swap the bytes; see PetscBinaryRead()
2636: and PetscBinaryWrite() to see how this may be done.
2638: Concepts: vector^loading from file
2640: .seealso: PetscViewerBinaryOpen(), DMView(), MatLoad(), VecLoad()
2641: @*/
2642: PetscErrorCode DMLoad(DM newdm, PetscViewer viewer)
2643: {
2650: if (!((PetscObject)newdm)->type_name) {
2651: DMSetType(newdm, DMDA);
2652: }
2653: (*newdm->ops->load)(newdm,viewer);
2654: return(0);
2655: }
2657: /******************************** FEM Support **********************************/
2661: PetscErrorCode DMPrintCellVector(PetscInt c, const char name[], PetscInt len, const PetscScalar x[]) {
2662: PetscInt f;
2666: PetscPrintf(PETSC_COMM_SELF, "Cell %D Element %s\n", c, name);
2667: for(f = 0; f < len; ++f) {
2668: PetscPrintf(PETSC_COMM_SELF, " | %G |\n", PetscRealPart(x[f]));
2669: }
2670: return(0);
2671: }
2675: PetscErrorCode DMPrintCellMatrix(PetscInt c, const char name[], PetscInt rows, PetscInt cols, const PetscScalar A[]) {
2676: PetscInt f, g;
2680: PetscPrintf(PETSC_COMM_SELF, "Cell %D Element %s\n", c, name);
2681: for(f = 0; f < rows; ++f) {
2682: PetscPrintf(PETSC_COMM_SELF, " |");
2683: for(g = 0; g < cols; ++g) {
2684: PetscPrintf(PETSC_COMM_SELF, " % 9.5G", PetscRealPart(A[f*cols+g]));
2685: }
2686: PetscPrintf(PETSC_COMM_SELF, " |\n");
2687: }
2688: return(0);
2689: }
2693: /*@C
2694: DMGetLocalFunction - Get the local residual function from this DM
2696: Not collective
2698: Input Parameter:
2699: . dm - The DM
2701: Output Parameter:
2702: . lf - The local residual function
2704: Calling sequence of lf:
2705: $ lf (SNES snes, Vec x, Vec f, void *ctx);
2707: + snes - the SNES context
2708: . x - local vector with the state at which to evaluate residual
2709: . f - local vector to put residual in
2710: - ctx - optional user-defined function context
2712: Level: intermediate
2714: .seealso DMSetLocalFunction(), DMGetLocalJacobian(), DMSetLocalJacobian()
2715: @*/
2716: PetscErrorCode DMGetLocalFunction(DM dm, PetscErrorCode (**lf)(DM, Vec, Vec, void *))
2717: {
2720: if (lf) *lf = dm->lf;
2721: return(0);
2722: }
2726: /*@C
2727: DMSetLocalFunction - Set the local residual function from this DM
2729: Not collective
2731: Input Parameters:
2732: + dm - The DM
2733: - lf - The local residual function
2735: Calling sequence of lf:
2736: $ lf (SNES snes, Vec x, Vec f, void *ctx);
2738: + snes - the SNES context
2739: . x - local vector with the state at which to evaluate residual
2740: . f - local vector to put residual in
2741: - ctx - optional user-defined function context
2743: Level: intermediate
2745: .seealso DMGetLocalFunction(), DMGetLocalJacobian(), DMSetLocalJacobian()
2746: @*/
2747: PetscErrorCode DMSetLocalFunction(DM dm, PetscErrorCode (*lf)(DM, Vec, Vec, void *))
2748: {
2751: dm->lf = lf;
2752: return(0);
2753: }
2757: /*@C
2758: DMGetLocalJacobian - Get the local Jacobian function from this DM
2760: Not collective
2762: Input Parameter:
2763: . dm - The DM
2765: Output Parameter:
2766: . lj - The local Jacobian function
2768: Calling sequence of lj:
2769: $ lj (SNES snes, Vec x, Mat J, Mat M, void *ctx);
2771: + snes - the SNES context
2772: . x - local vector with the state at which to evaluate residual
2773: . J - matrix to put Jacobian in
2774: . M - matrix to use for defining Jacobian preconditioner
2775: - ctx - optional user-defined function context
2777: Level: intermediate
2779: .seealso DMSetLocalJacobian(), DMGetLocalFunction(), DMSetLocalFunction()
2780: @*/
2781: PetscErrorCode DMGetLocalJacobian(DM dm, PetscErrorCode (**lj)(DM, Vec, Mat, Mat, void *))
2782: {
2785: if (lj) *lj = dm->lj;
2786: return(0);
2787: }
2791: /*@C
2792: DMSetLocalJacobian - Set the local Jacobian function from this DM
2794: Not collective
2796: Input Parameters:
2797: + dm - The DM
2798: - lj - The local Jacobian function
2800: Calling sequence of lj:
2801: $ lj (SNES snes, Vec x, Mat J, Mat M, void *ctx);
2803: + snes - the SNES context
2804: . x - local vector with the state at which to evaluate residual
2805: . J - matrix to put Jacobian in
2806: . M - matrix to use for defining Jacobian preconditioner
2807: - ctx - optional user-defined function context
2809: Level: intermediate
2811: .seealso DMGetLocalJacobian(), DMGetLocalFunction(), DMSetLocalFunction()
2812: @*/
2813: PetscErrorCode DMSetLocalJacobian(DM dm, PetscErrorCode (*lj)(DM, Vec, Mat, Mat, void *))
2814: {
2817: dm->lj = lj;
2818: return(0);
2819: }
2823: /*@
2824: DMGetDefaultSection - Get the PetscSection encoding the local data layout for the DM.
2826: Input Parameter:
2827: . dm - The DM
2829: Output Parameter:
2830: . section - The PetscSection
2832: Level: intermediate
2834: Note: This gets a borrowed reference, so the user should not destroy this PetscSection.
2836: .seealso: DMSetDefaultSection(), DMGetDefaultGlobalSection()
2837: @*/
2838: PetscErrorCode DMGetDefaultSection(DM dm, PetscSection *section) {
2842: *section = dm->defaultSection;
2843: return(0);
2844: }
2848: /*@
2849: DMSetDefaultSection - Set the PetscSection encoding the local data layout for the DM.
2851: Input Parameters:
2852: + dm - The DM
2853: - section - The PetscSection
2855: Level: intermediate
2857: Note: Any existing Section will be destroyed
2859: .seealso: DMSetDefaultSection(), DMGetDefaultGlobalSection()
2860: @*/
2861: PetscErrorCode DMSetDefaultSection(DM dm, PetscSection section) {
2866: PetscSectionDestroy(&dm->defaultSection);
2867: PetscSectionDestroy(&dm->defaultGlobalSection);
2868: dm->defaultSection = section;
2869: return(0);
2870: }
2874: /*@
2875: DMGetDefaultGlobalSection - Get the PetscSection encoding the global data layout for the DM.
2877: Input Parameter:
2878: . dm - The DM
2880: Output Parameter:
2881: . section - The PetscSection
2883: Level: intermediate
2885: Note: This gets a borrowed reference, so the user should not destroy this PetscSection.
2887: .seealso: DMSetDefaultSection(), DMGetDefaultSection()
2888: @*/
2889: PetscErrorCode DMGetDefaultGlobalSection(DM dm, PetscSection *section) {
2895: if (!dm->defaultGlobalSection) {
2896: PetscSectionCreateGlobalSection(dm->defaultSection, dm->sf, &dm->defaultGlobalSection);
2897: }
2898: *section = dm->defaultGlobalSection;
2899: return(0);
2900: }
2904: /*@
2905: DMGetDefaultSF - Get the PetscSF encoding the parallel dof overlap for the DM. If it has not been set,
2906: it is created from the default PetscSection layouts in the DM.
2908: Input Parameter:
2909: . dm - The DM
2911: Output Parameter:
2912: . sf - The PetscSF
2914: Level: intermediate
2916: Note: This gets a borrowed reference, so the user should not destroy this PetscSF.
2918: .seealso: DMSetDefaultSF(), DMCreateDefaultSF()
2919: @*/
2920: PetscErrorCode DMGetDefaultSF(DM dm, PetscSF *sf) {
2921: PetscInt nroots;
2927: PetscSFGetGraph(dm->defaultSF, &nroots, PETSC_NULL, PETSC_NULL, PETSC_NULL);
2928: if (nroots < 0) {
2929: PetscSection section, gSection;
2931: DMGetDefaultSection(dm, §ion);
2932: if (section) {
2933: DMGetDefaultGlobalSection(dm, &gSection);
2934: DMCreateDefaultSF(dm, section, gSection);
2935: } else {
2936: *sf = PETSC_NULL;
2937: return(0);
2938: }
2939: }
2940: *sf = dm->defaultSF;
2941: return(0);
2942: }
2946: /*@
2947: DMSetDefaultSF - Set the PetscSF encoding the parallel dof overlap for the DM
2949: Input Parameters:
2950: + dm - The DM
2951: - sf - The PetscSF
2953: Level: intermediate
2955: Note: Any previous SF is destroyed
2957: .seealso: DMGetDefaultSF(), DMCreateDefaultSF()
2958: @*/
2959: PetscErrorCode DMSetDefaultSF(DM dm, PetscSF sf) {
2965: PetscSFDestroy(&dm->defaultSF);
2966: dm->defaultSF = sf;
2967: return(0);
2968: }
2972: /*@C
2973: DMCreateDefaultSF - Create the PetscSF encoding the parallel dof overlap for the DM based upon the PetscSections
2974: describing the data layout.
2976: Input Parameters:
2977: + dm - The DM
2978: . localSection - PetscSection describing the local data layout
2979: - globalSection - PetscSection describing the global data layout
2981: Level: intermediate
2983: .seealso: DMGetDefaultSF(), DMSetDefaultSF()
2984: @*/
2985: PetscErrorCode DMCreateDefaultSF(DM dm, PetscSection localSection, PetscSection globalSection)
2986: {
2987: MPI_Comm comm = ((PetscObject) dm)->comm;
2988: PetscLayout layout;
2989: const PetscInt *ranges;
2990: PetscInt *local;
2991: PetscSFNode *remote;
2992: PetscInt pStart, pEnd, p, nroots, nleaves, l;
2993: PetscMPIInt size, rank;
2994: PetscErrorCode ierr;
2998: MPI_Comm_size(comm, &size);
2999: MPI_Comm_rank(comm, &rank);
3000: PetscSectionGetChart(globalSection, &pStart, &pEnd);
3001: PetscSectionGetConstrainedStorageSize(globalSection, &nroots);
3002: PetscLayoutCreate(comm, &layout);
3003: PetscLayoutSetBlockSize(layout, 1);
3004: PetscLayoutSetLocalSize(layout, nroots);
3005: PetscLayoutSetUp(layout);
3006: PetscLayoutGetRanges(layout, &ranges);
3007: for(p = pStart, nleaves = 0; p < pEnd; ++p) {
3008: PetscInt dof, cdof;
3010: PetscSectionGetDof(globalSection, p, &dof);
3011: PetscSectionGetConstraintDof(globalSection, p, &cdof);
3012: nleaves += dof < 0 ? -(dof+1)-cdof : dof-cdof;
3013: }
3014: PetscMalloc(nleaves * sizeof(PetscInt), &local);
3015: PetscMalloc(nleaves * sizeof(PetscSFNode), &remote);
3016: for(p = pStart, l = 0; p < pEnd; ++p) {
3017: PetscInt *cind;
3018: PetscInt dof, gdof, cdof, dim, off, goff, d, c;
3020: PetscSectionGetDof(localSection, p, &dof);
3021: PetscSectionGetOffset(localSection, p, &off);
3022: PetscSectionGetConstraintDof(localSection, p, &cdof);
3023: PetscSectionGetConstraintIndices(localSection, p, &cind);
3024: PetscSectionGetDof(globalSection, p, &gdof);
3025: PetscSectionGetOffset(globalSection, p, &goff);
3026: dim = dof-cdof;
3027: for(d = 0, c = 0; d < dof; ++d) {
3028: if ((c < cdof) && (cind[c] == d)) {++c; continue;}
3029: local[l+d-c] = off+d;
3030: }
3031: if (gdof < 0) {
3032: for(d = 0; d < dim; ++d, ++l) {
3033: PetscInt offset = -(goff+1) + d, r;
3035: for(r = 0; r < size; ++r) {
3036: if ((offset >= ranges[r]) && (offset < ranges[r+1])) break;
3037: }
3038: remote[l].rank = r;
3039: remote[l].index = offset - ranges[r];
3040: }
3041: } else {
3042: for(d = 0; d < dim; ++d, ++l) {
3043: remote[l].rank = rank;
3044: remote[l].index = goff+d - ranges[rank];
3045: }
3046: }
3047: }
3048: PetscLayoutDestroy(&layout);
3049: PetscSFSetGraph(dm->defaultSF, nroots, nleaves, local, PETSC_OWN_POINTER, remote, PETSC_OWN_POINTER);
3050: return(0);
3051: }