Actual source code: dm.c
petsc-3.10.5 2019-03-28
1: #include <petsc/private/dmimpl.h>
2: #include <petsc/private/dmlabelimpl.h>
3: #include <petsc/private/petscdsimpl.h>
4: #include <petscdmplex.h>
5: #include <petscdmfield.h>
6: #include <petscsf.h>
7: #include <petscds.h>
9: PetscClassId DM_CLASSID;
10: PetscLogEvent DM_Convert, DM_GlobalToLocal, DM_LocalToGlobal, DM_LocalToLocal, DM_LocatePoints, DM_Coarsen, DM_Refine, DM_CreateInterpolation, DM_CreateRestriction;
12: const char *const DMBoundaryTypes[] = {"NONE","GHOSTED","MIRROR","PERIODIC","TWIST","DMBoundaryType","DM_BOUNDARY_",0};
14: static PetscErrorCode DMHasCreateInjection_Default(DM dm, PetscBool *flg)
15: {
19: *flg = PETSC_FALSE;
20: return(0);
21: }
23: /*@
24: DMCreate - Creates an empty DM object. The type can then be set with DMSetType().
26: If you never call DMSetType() it will generate an
27: error when you try to use the vector.
29: Collective on MPI_Comm
31: Input Parameter:
32: . comm - The communicator for the DM object
34: Output Parameter:
35: . dm - The DM object
37: Level: beginner
39: .seealso: DMSetType(), DMDA, DMSLICED, DMCOMPOSITE, DMPLEX, DMMOAB, DMNETWORK
40: @*/
41: PetscErrorCode DMCreate(MPI_Comm comm,DM *dm)
42: {
43: DM v;
48: *dm = NULL;
49: DMInitializePackage();
51: PetscHeaderCreate(v, DM_CLASSID, "DM", "Distribution Manager", "DM", comm, DMDestroy, DMView);
53: v->ltogmap = NULL;
54: v->bs = 1;
55: v->coloringtype = IS_COLORING_GLOBAL;
56: PetscSFCreate(comm, &v->sf);
57: PetscSFCreate(comm, &v->defaultSF);
58: v->labels = NULL;
59: v->depthLabel = NULL;
60: v->defaultSection = NULL;
61: v->defaultGlobalSection = NULL;
62: v->defaultConstraintSection = NULL;
63: v->defaultConstraintMat = NULL;
64: v->L = NULL;
65: v->maxCell = NULL;
66: v->bdtype = NULL;
67: v->dimEmbed = PETSC_DEFAULT;
68: v->dim = PETSC_DETERMINE;
69: {
70: PetscInt i;
71: for (i = 0; i < 10; ++i) {
72: v->nullspaceConstructors[i] = NULL;
73: }
74: }
75: PetscDSCreate(comm, &v->prob);
76: v->dmBC = NULL;
77: v->coarseMesh = NULL;
78: v->outputSequenceNum = -1;
79: v->outputSequenceVal = 0.0;
80: DMSetVecType(v,VECSTANDARD);
81: DMSetMatType(v,MATAIJ);
82: PetscNew(&(v->labels));
83: v->labels->refct = 1;
85: v->ops->hascreateinjection = DMHasCreateInjection_Default;
87: *dm = v;
88: return(0);
89: }
91: /*@
92: DMClone - Creates a DM object with the same topology as the original.
94: Collective on MPI_Comm
96: Input Parameter:
97: . dm - The original DM object
99: Output Parameter:
100: . newdm - The new DM object
102: Level: beginner
104: .keywords: DM, topology, create
105: @*/
106: PetscErrorCode DMClone(DM dm, DM *newdm)
107: {
108: PetscSF sf;
109: Vec coords;
110: void *ctx;
111: PetscInt dim, cdim;
117: DMCreate(PetscObjectComm((PetscObject) dm), newdm);
118: PetscFree((*newdm)->labels);
119: dm->labels->refct++;
120: (*newdm)->labels = dm->labels;
121: (*newdm)->depthLabel = dm->depthLabel;
122: DMGetDimension(dm, &dim);
123: DMSetDimension(*newdm, dim);
124: if (dm->ops->clone) {
125: (*dm->ops->clone)(dm, newdm);
126: }
127: (*newdm)->setupcalled = dm->setupcalled;
128: DMGetPointSF(dm, &sf);
129: DMSetPointSF(*newdm, sf);
130: DMGetApplicationContext(dm, &ctx);
131: DMSetApplicationContext(*newdm, ctx);
132: if (dm->coordinateDM) {
133: DM ncdm;
134: PetscSection cs;
135: PetscInt pEnd = -1, pEndMax = -1;
137: DMGetSection(dm->coordinateDM, &cs);
138: if (cs) {PetscSectionGetChart(cs, NULL, &pEnd);}
139: MPI_Allreduce(&pEnd,&pEndMax,1,MPIU_INT,MPI_MAX,PetscObjectComm((PetscObject)dm));
140: if (pEndMax >= 0) {
141: DMClone(dm->coordinateDM, &ncdm);
142: DMSetSection(ncdm, cs);
143: DMSetCoordinateDM(*newdm, ncdm);
144: DMDestroy(&ncdm);
145: }
146: }
147: DMGetCoordinateDim(dm, &cdim);
148: DMSetCoordinateDim(*newdm, cdim);
149: DMGetCoordinatesLocal(dm, &coords);
150: if (coords) {
151: DMSetCoordinatesLocal(*newdm, coords);
152: } else {
153: DMGetCoordinates(dm, &coords);
154: if (coords) {DMSetCoordinates(*newdm, coords);}
155: }
156: {
157: PetscBool isper;
158: const PetscReal *maxCell, *L;
159: const DMBoundaryType *bd;
160: DMGetPeriodicity(dm, &isper, &maxCell, &L, &bd);
161: DMSetPeriodicity(*newdm, isper, maxCell, L, bd);
162: }
163: return(0);
164: }
166: /*@C
167: DMSetVecType - Sets the type of vector created with DMCreateLocalVector() and DMCreateGlobalVector()
169: Logically Collective on DM
171: Input Parameter:
172: + da - initial distributed array
173: . ctype - the vector type, currently either VECSTANDARD, VECCUDA, or VECVIENNACL
175: Options Database:
176: . -dm_vec_type ctype
178: Level: intermediate
180: .seealso: DMCreate(), DMDestroy(), DM, DMDAInterpolationType, VecType, DMGetVecType()
181: @*/
182: PetscErrorCode DMSetVecType(DM da,VecType ctype)
183: {
188: PetscFree(da->vectype);
189: PetscStrallocpy(ctype,(char**)&da->vectype);
190: return(0);
191: }
193: /*@C
194: DMGetVecType - Gets the type of vector created with DMCreateLocalVector() and DMCreateGlobalVector()
196: Logically Collective on DM
198: Input Parameter:
199: . da - initial distributed array
201: Output Parameter:
202: . ctype - the vector type
204: Level: intermediate
206: .seealso: DMCreate(), DMDestroy(), DM, DMDAInterpolationType, VecType
207: @*/
208: PetscErrorCode DMGetVecType(DM da,VecType *ctype)
209: {
212: *ctype = da->vectype;
213: return(0);
214: }
216: /*@
217: VecGetDM - Gets the DM defining the data layout of the vector
219: Not collective
221: Input Parameter:
222: . v - The Vec
224: Output Parameter:
225: . dm - The DM
227: Level: intermediate
229: .seealso: VecSetDM(), DMGetLocalVector(), DMGetGlobalVector(), DMSetVecType()
230: @*/
231: PetscErrorCode VecGetDM(Vec v, DM *dm)
232: {
238: PetscObjectQuery((PetscObject) v, "__PETSc_dm", (PetscObject*) dm);
239: return(0);
240: }
242: /*@
243: VecSetDM - Sets the DM defining the data layout of the vector.
245: Not collective
247: Input Parameters:
248: + v - The Vec
249: - dm - The DM
251: Note: This is NOT the same as DMCreateGlobalVector() since it does not change the view methods or perform other customization, but merely sets the DM member.
253: Level: intermediate
255: .seealso: VecGetDM(), DMGetLocalVector(), DMGetGlobalVector(), DMSetVecType()
256: @*/
257: PetscErrorCode VecSetDM(Vec v, DM dm)
258: {
264: PetscObjectCompose((PetscObject) v, "__PETSc_dm", (PetscObject) dm);
265: return(0);
266: }
268: /*@C
269: DMSetISColoringType - Sets the type of coloring, global or local, that is created by the DM
271: Logically Collective on DM
273: Input Parameters:
274: + dm - the DM context
275: - ctype - the matrix type
277: Options Database:
278: . -dm_is_coloring_type - global or local
280: Level: intermediate
282: .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMCreateMatrix(), DMSetMatrixPreallocateOnly(), MatType, DMGetMatType(),
283: DMGetISColoringType()
284: @*/
285: PetscErrorCode DMSetISColoringType(DM dm,ISColoringType ctype)
286: {
289: dm->coloringtype = ctype;
290: return(0);
291: }
293: /*@C
294: DMGetISColoringType - Gets the type of coloring, global or local, that is created by the DM
296: Logically Collective on DM
298: Input Parameter:
299: . dm - the DM context
301: Output Parameter:
302: . ctype - the matrix type
304: Options Database:
305: . -dm_is_coloring_type - global or local
307: Level: intermediate
309: .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMCreateMatrix(), DMSetMatrixPreallocateOnly(), MatType, DMGetMatType(),
310: DMGetISColoringType()
311: @*/
312: PetscErrorCode DMGetISColoringType(DM dm,ISColoringType *ctype)
313: {
316: *ctype = dm->coloringtype;
317: return(0);
318: }
320: /*@C
321: DMSetMatType - Sets the type of matrix created with DMCreateMatrix()
323: Logically Collective on DM
325: Input Parameters:
326: + dm - the DM context
327: - ctype - the matrix type
329: Options Database:
330: . -dm_mat_type ctype
332: Level: intermediate
334: .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMCreateMatrix(), DMSetMatrixPreallocateOnly(), MatType, DMGetMatType()
335: @*/
336: PetscErrorCode DMSetMatType(DM dm,MatType ctype)
337: {
342: PetscFree(dm->mattype);
343: PetscStrallocpy(ctype,(char**)&dm->mattype);
344: return(0);
345: }
347: /*@C
348: DMGetMatType - Gets the type of matrix created with DMCreateMatrix()
350: Logically Collective on DM
352: Input Parameter:
353: . dm - the DM context
355: Output Parameter:
356: . ctype - the matrix type
358: Options Database:
359: . -dm_mat_type ctype
361: Level: intermediate
363: .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMCreateMatrix(), DMSetMatrixPreallocateOnly(), MatType, DMSetMatType()
364: @*/
365: PetscErrorCode DMGetMatType(DM dm,MatType *ctype)
366: {
369: *ctype = dm->mattype;
370: return(0);
371: }
373: /*@
374: MatGetDM - Gets the DM defining the data layout of the matrix
376: Not collective
378: Input Parameter:
379: . A - The Mat
381: Output Parameter:
382: . dm - The DM
384: Level: intermediate
386: Developer Note: Since the Mat class doesn't know about the DM class the DM object is associated with
387: the Mat through a PetscObjectCompose() operation
389: .seealso: MatSetDM(), DMCreateMatrix(), DMSetMatType()
390: @*/
391: PetscErrorCode MatGetDM(Mat A, DM *dm)
392: {
398: PetscObjectQuery((PetscObject) A, "__PETSc_dm", (PetscObject*) dm);
399: return(0);
400: }
402: /*@
403: MatSetDM - Sets the DM defining the data layout of the matrix
405: Not collective
407: Input Parameters:
408: + A - The Mat
409: - dm - The DM
411: Level: intermediate
413: Developer Note: Since the Mat class doesn't know about the DM class the DM object is associated with
414: the Mat through a PetscObjectCompose() operation
417: .seealso: MatGetDM(), DMCreateMatrix(), DMSetMatType()
418: @*/
419: PetscErrorCode MatSetDM(Mat A, DM dm)
420: {
426: PetscObjectCompose((PetscObject) A, "__PETSc_dm", (PetscObject) dm);
427: return(0);
428: }
430: /*@C
431: DMSetOptionsPrefix - Sets the prefix used for searching for all
432: DM options in the database.
434: Logically Collective on DM
436: Input Parameter:
437: + da - the DM context
438: - prefix - the prefix to prepend to all option names
440: Notes:
441: A hyphen (-) must NOT be given at the beginning of the prefix name.
442: The first character of all runtime options is AUTOMATICALLY the hyphen.
444: Level: advanced
446: .keywords: DM, set, options, prefix, database
448: .seealso: DMSetFromOptions()
449: @*/
450: PetscErrorCode DMSetOptionsPrefix(DM dm,const char prefix[])
451: {
456: PetscObjectSetOptionsPrefix((PetscObject)dm,prefix);
457: if (dm->sf) {
458: PetscObjectSetOptionsPrefix((PetscObject)dm->sf,prefix);
459: }
460: if (dm->defaultSF) {
461: PetscObjectSetOptionsPrefix((PetscObject)dm->defaultSF,prefix);
462: }
463: return(0);
464: }
466: /*@C
467: DMAppendOptionsPrefix - Appends to the prefix used for searching for all
468: DM options in the database.
470: Logically Collective on DM
472: Input Parameters:
473: + dm - the DM context
474: - prefix - the prefix string to prepend to all DM option requests
476: Notes:
477: A hyphen (-) must NOT be given at the beginning of the prefix name.
478: The first character of all runtime options is AUTOMATICALLY the hyphen.
480: Level: advanced
482: .keywords: DM, append, options, prefix, database
484: .seealso: DMSetOptionsPrefix(), DMGetOptionsPrefix()
485: @*/
486: PetscErrorCode DMAppendOptionsPrefix(DM dm,const char prefix[])
487: {
492: PetscObjectAppendOptionsPrefix((PetscObject)dm,prefix);
493: return(0);
494: }
496: /*@C
497: DMGetOptionsPrefix - Gets the prefix used for searching for all
498: DM options in the database.
500: Not Collective
502: Input Parameters:
503: . dm - the DM context
505: Output Parameters:
506: . prefix - pointer to the prefix string used is returned
508: Notes:
509: On the fortran side, the user should pass in a string 'prefix' of
510: sufficient length to hold the prefix.
512: Level: advanced
514: .keywords: DM, set, options, prefix, database
516: .seealso: DMSetOptionsPrefix(), DMAppendOptionsPrefix()
517: @*/
518: PetscErrorCode DMGetOptionsPrefix(DM dm,const char *prefix[])
519: {
524: PetscObjectGetOptionsPrefix((PetscObject)dm,prefix);
525: return(0);
526: }
528: static PetscErrorCode DMCountNonCyclicReferences(DM dm, PetscBool recurseCoarse, PetscBool recurseFine, PetscInt *ncrefct)
529: {
530: PetscInt i, refct = ((PetscObject) dm)->refct;
531: DMNamedVecLink nlink;
535: *ncrefct = 0;
536: /* count all the circular references of DM and its contained Vecs */
537: for (i=0; i<DM_MAX_WORK_VECTORS; i++) {
538: if (dm->localin[i]) refct--;
539: if (dm->globalin[i]) refct--;
540: }
541: for (nlink=dm->namedglobal; nlink; nlink=nlink->next) refct--;
542: for (nlink=dm->namedlocal; nlink; nlink=nlink->next) refct--;
543: if (dm->x) {
544: DM obj;
545: VecGetDM(dm->x, &obj);
546: if (obj == dm) refct--;
547: }
548: if (dm->coarseMesh && dm->coarseMesh->fineMesh == dm) {
549: refct--;
550: if (recurseCoarse) {
551: PetscInt coarseCount;
553: DMCountNonCyclicReferences(dm->coarseMesh, PETSC_TRUE, PETSC_FALSE,&coarseCount);
554: refct += coarseCount;
555: }
556: }
557: if (dm->fineMesh && dm->fineMesh->coarseMesh == dm) {
558: refct--;
559: if (recurseFine) {
560: PetscInt fineCount;
562: DMCountNonCyclicReferences(dm->fineMesh, PETSC_FALSE, PETSC_TRUE,&fineCount);
563: refct += fineCount;
564: }
565: }
566: *ncrefct = refct;
567: return(0);
568: }
570: PetscErrorCode DMDestroyLabelLinkList(DM dm)
571: {
575: if (!--(dm->labels->refct)) {
576: DMLabelLink next = dm->labels->next;
578: /* destroy the labels */
579: while (next) {
580: DMLabelLink tmp = next->next;
582: DMLabelDestroy(&next->label);
583: PetscFree(next);
584: next = tmp;
585: }
586: PetscFree(dm->labels);
587: }
588: return(0);
589: }
591: /*@
592: DMDestroy - Destroys a vector packer or DM.
594: Collective on DM
596: Input Parameter:
597: . dm - the DM object to destroy
599: Level: developer
601: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix()
603: @*/
604: PetscErrorCode DMDestroy(DM *dm)
605: {
606: PetscInt i, cnt;
607: DMNamedVecLink nlink,nnext;
611: if (!*dm) return(0);
614: /* count all non-cyclic references in the doubly-linked list of coarse<->fine meshes */
615: DMCountNonCyclicReferences(*dm,PETSC_TRUE,PETSC_TRUE,&cnt);
616: --((PetscObject)(*dm))->refct;
617: if (--cnt > 0) {*dm = 0; return(0);}
618: /*
619: Need this test because the dm references the vectors that
620: reference the dm, so destroying the dm calls destroy on the
621: vectors that cause another destroy on the dm
622: */
623: if (((PetscObject)(*dm))->refct < 0) return(0);
624: ((PetscObject) (*dm))->refct = 0;
625: for (i=0; i<DM_MAX_WORK_VECTORS; i++) {
626: if ((*dm)->localout[i]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Destroying a DM that has a local vector obtained with DMGetLocalVector()");
627: VecDestroy(&(*dm)->localin[i]);
628: }
629: nnext=(*dm)->namedglobal;
630: (*dm)->namedglobal = NULL;
631: for (nlink=nnext; nlink; nlink=nnext) { /* Destroy the named vectors */
632: nnext = nlink->next;
633: if (nlink->status != DMVEC_STATUS_IN) SETERRQ1(((PetscObject)*dm)->comm,PETSC_ERR_ARG_WRONGSTATE,"DM still has Vec named '%s' checked out",nlink->name);
634: PetscFree(nlink->name);
635: VecDestroy(&nlink->X);
636: PetscFree(nlink);
637: }
638: nnext=(*dm)->namedlocal;
639: (*dm)->namedlocal = NULL;
640: for (nlink=nnext; nlink; nlink=nnext) { /* Destroy the named local vectors */
641: nnext = nlink->next;
642: if (nlink->status != DMVEC_STATUS_IN) SETERRQ1(((PetscObject)*dm)->comm,PETSC_ERR_ARG_WRONGSTATE,"DM still has Vec named '%s' checked out",nlink->name);
643: PetscFree(nlink->name);
644: VecDestroy(&nlink->X);
645: PetscFree(nlink);
646: }
648: /* Destroy the list of hooks */
649: {
650: DMCoarsenHookLink link,next;
651: for (link=(*dm)->coarsenhook; link; link=next) {
652: next = link->next;
653: PetscFree(link);
654: }
655: (*dm)->coarsenhook = NULL;
656: }
657: {
658: DMRefineHookLink link,next;
659: for (link=(*dm)->refinehook; link; link=next) {
660: next = link->next;
661: PetscFree(link);
662: }
663: (*dm)->refinehook = NULL;
664: }
665: {
666: DMSubDomainHookLink link,next;
667: for (link=(*dm)->subdomainhook; link; link=next) {
668: next = link->next;
669: PetscFree(link);
670: }
671: (*dm)->subdomainhook = NULL;
672: }
673: {
674: DMGlobalToLocalHookLink link,next;
675: for (link=(*dm)->gtolhook; link; link=next) {
676: next = link->next;
677: PetscFree(link);
678: }
679: (*dm)->gtolhook = NULL;
680: }
681: {
682: DMLocalToGlobalHookLink link,next;
683: for (link=(*dm)->ltoghook; link; link=next) {
684: next = link->next;
685: PetscFree(link);
686: }
687: (*dm)->ltoghook = NULL;
688: }
689: /* Destroy the work arrays */
690: {
691: DMWorkLink link,next;
692: if ((*dm)->workout) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Work array still checked out");
693: for (link=(*dm)->workin; link; link=next) {
694: next = link->next;
695: PetscFree(link->mem);
696: PetscFree(link);
697: }
698: (*dm)->workin = NULL;
699: }
700: if (!--((*dm)->labels->refct)) {
701: DMLabelLink next = (*dm)->labels->next;
703: /* destroy the labels */
704: while (next) {
705: DMLabelLink tmp = next->next;
707: DMLabelDestroy(&next->label);
708: PetscFree(next);
709: next = tmp;
710: }
711: PetscFree((*dm)->labels);
712: }
713: {
714: DMBoundary next = (*dm)->boundary;
715: while (next) {
716: DMBoundary b = next;
718: next = b->next;
719: PetscFree(b);
720: }
721: }
723: PetscObjectDestroy(&(*dm)->dmksp);
724: PetscObjectDestroy(&(*dm)->dmsnes);
725: PetscObjectDestroy(&(*dm)->dmts);
727: if ((*dm)->ctx && (*dm)->ctxdestroy) {
728: (*(*dm)->ctxdestroy)(&(*dm)->ctx);
729: }
730: VecDestroy(&(*dm)->x);
731: MatFDColoringDestroy(&(*dm)->fd);
732: DMClearGlobalVectors(*dm);
733: ISLocalToGlobalMappingDestroy(&(*dm)->ltogmap);
734: PetscFree((*dm)->vectype);
735: PetscFree((*dm)->mattype);
737: PetscSectionDestroy(&(*dm)->defaultSection);
738: PetscSectionDestroy(&(*dm)->defaultGlobalSection);
739: PetscLayoutDestroy(&(*dm)->map);
740: PetscSectionDestroy(&(*dm)->defaultConstraintSection);
741: MatDestroy(&(*dm)->defaultConstraintMat);
742: PetscSFDestroy(&(*dm)->sf);
743: PetscSFDestroy(&(*dm)->defaultSF);
744: if ((*dm)->useNatural) {
745: if ((*dm)->sfNatural) {
746: PetscSFDestroy(&(*dm)->sfNatural);
747: }
748: PetscObjectDereference((PetscObject) (*dm)->sfMigration);
749: }
750: if ((*dm)->coarseMesh && (*dm)->coarseMesh->fineMesh == *dm) {
751: DMSetFineDM((*dm)->coarseMesh,NULL);
752: }
753: DMDestroy(&(*dm)->coarseMesh);
754: if ((*dm)->fineMesh && (*dm)->fineMesh->coarseMesh == *dm) {
755: DMSetCoarseDM((*dm)->fineMesh,NULL);
756: }
757: DMDestroy(&(*dm)->fineMesh);
758: DMFieldDestroy(&(*dm)->coordinateField);
759: DMDestroy(&(*dm)->coordinateDM);
760: VecDestroy(&(*dm)->coordinates);
761: VecDestroy(&(*dm)->coordinatesLocal);
762: PetscFree3((*dm)->L,(*dm)->maxCell,(*dm)->bdtype);
764: PetscDSDestroy(&(*dm)->prob);
765: DMDestroy(&(*dm)->dmBC);
766: /* if memory was published with SAWs then destroy it */
767: PetscObjectSAWsViewOff((PetscObject)*dm);
769: if ((*dm)->ops->destroy) {
770: (*(*dm)->ops->destroy)(*dm);
771: }
772: /* We do not destroy (*dm)->data here so that we can reference count backend objects */
773: PetscHeaderDestroy(dm);
774: return(0);
775: }
777: /*@
778: DMSetUp - sets up the data structures inside a DM object
780: Collective on DM
782: Input Parameter:
783: . dm - the DM object to setup
785: Level: developer
787: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix()
789: @*/
790: PetscErrorCode DMSetUp(DM dm)
791: {
796: if (dm->setupcalled) return(0);
797: if (dm->ops->setup) {
798: (*dm->ops->setup)(dm);
799: }
800: dm->setupcalled = PETSC_TRUE;
801: return(0);
802: }
804: /*@
805: DMSetFromOptions - sets parameters in a DM from the options database
807: Collective on DM
809: Input Parameter:
810: . dm - the DM object to set options for
812: Options Database:
813: + -dm_preallocate_only - Only preallocate the matrix for DMCreateMatrix(), but do not fill it with zeros
814: . -dm_vec_type <type> - type of vector to create inside DM
815: . -dm_mat_type <type> - type of matrix to create inside DM
816: - -dm_is_coloring_type - <global or local>
818: Level: developer
820: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix()
822: @*/
823: PetscErrorCode DMSetFromOptions(DM dm)
824: {
825: char typeName[256];
826: PetscBool flg;
831: if (dm->prob) {
832: PetscDSSetFromOptions(dm->prob);
833: }
834: if (dm->sf) {
835: PetscSFSetFromOptions(dm->sf);
836: }
837: if (dm->defaultSF) {
838: PetscSFSetFromOptions(dm->defaultSF);
839: }
840: PetscObjectOptionsBegin((PetscObject)dm);
841: PetscOptionsBool("-dm_preallocate_only","only preallocate matrix, but do not set column indices","DMSetMatrixPreallocateOnly",dm->prealloc_only,&dm->prealloc_only,NULL);
842: PetscOptionsFList("-dm_vec_type","Vector type used for created vectors","DMSetVecType",VecList,dm->vectype,typeName,256,&flg);
843: if (flg) {
844: DMSetVecType(dm,typeName);
845: }
846: PetscOptionsFList("-dm_mat_type","Matrix type used for created matrices","DMSetMatType",MatList,dm->mattype ? dm->mattype : typeName,typeName,sizeof(typeName),&flg);
847: if (flg) {
848: DMSetMatType(dm,typeName);
849: }
850: PetscOptionsEnum("-dm_is_coloring_type","Global or local coloring of Jacobian","DMSetISColoringType",ISColoringTypes,(PetscEnum)dm->coloringtype,(PetscEnum*)&dm->coloringtype,NULL);
851: if (dm->ops->setfromoptions) {
852: (*dm->ops->setfromoptions)(PetscOptionsObject,dm);
853: }
854: /* process any options handlers added with PetscObjectAddOptionsHandler() */
855: PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject) dm);
856: PetscOptionsEnd();
857: return(0);
858: }
860: /*@C
861: DMView - Views a DM
863: Collective on DM
865: Input Parameter:
866: + dm - the DM object to view
867: - v - the viewer
869: Level: beginner
871: .seealso DMDestroy(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix()
873: @*/
874: PetscErrorCode DMView(DM dm,PetscViewer v)
875: {
876: PetscErrorCode ierr;
877: PetscBool isbinary;
878: PetscMPIInt size;
879: PetscViewerFormat format;
883: if (!v) {
884: PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)dm),&v);
885: }
886: PetscViewerGetFormat(v,&format);
887: MPI_Comm_size(PetscObjectComm((PetscObject)dm),&size);
888: if (size == 1 && format == PETSC_VIEWER_LOAD_BALANCE) return(0);
889: PetscObjectPrintClassNamePrefixType((PetscObject)dm,v);
890: PetscObjectTypeCompare((PetscObject)v,PETSCVIEWERBINARY,&isbinary);
891: if (isbinary) {
892: PetscInt classid = DM_FILE_CLASSID;
893: char type[256];
895: PetscViewerBinaryWrite(v,&classid,1,PETSC_INT,PETSC_FALSE);
896: PetscStrncpy(type,((PetscObject)dm)->type_name,256);
897: PetscViewerBinaryWrite(v,type,256,PETSC_CHAR,PETSC_FALSE);
898: }
899: if (dm->ops->view) {
900: (*dm->ops->view)(dm,v);
901: }
902: return(0);
903: }
905: /*@
906: DMCreateGlobalVector - Creates a global vector from a DM object
908: Collective on DM
910: Input Parameter:
911: . dm - the DM object
913: Output Parameter:
914: . vec - the global vector
916: Level: beginner
918: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix()
920: @*/
921: PetscErrorCode DMCreateGlobalVector(DM dm,Vec *vec)
922: {
927: (*dm->ops->createglobalvector)(dm,vec);
928: return(0);
929: }
931: /*@
932: DMCreateLocalVector - Creates a local vector from a DM object
934: Not Collective
936: Input Parameter:
937: . dm - the DM object
939: Output Parameter:
940: . vec - the local vector
942: Level: beginner
944: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix()
946: @*/
947: PetscErrorCode DMCreateLocalVector(DM dm,Vec *vec)
948: {
953: (*dm->ops->createlocalvector)(dm,vec);
954: return(0);
955: }
957: /*@
958: DMGetLocalToGlobalMapping - Accesses the local-to-global mapping in a DM.
960: Collective on DM
962: Input Parameter:
963: . dm - the DM that provides the mapping
965: Output Parameter:
966: . ltog - the mapping
968: Level: intermediate
970: Notes:
971: This mapping can then be used by VecSetLocalToGlobalMapping() or
972: MatSetLocalToGlobalMapping().
974: .seealso: DMCreateLocalVector()
975: @*/
976: PetscErrorCode DMGetLocalToGlobalMapping(DM dm,ISLocalToGlobalMapping *ltog)
977: {
978: PetscInt bs = -1, bsLocal[2], bsMinMax[2];
984: if (!dm->ltogmap) {
985: PetscSection section, sectionGlobal;
987: DMGetSection(dm, §ion);
988: if (section) {
989: const PetscInt *cdofs;
990: PetscInt *ltog;
991: PetscInt pStart, pEnd, n, p, k, l;
993: DMGetGlobalSection(dm, §ionGlobal);
994: PetscSectionGetChart(section, &pStart, &pEnd);
995: PetscSectionGetStorageSize(section, &n);
996: PetscMalloc1(n, <og); /* We want the local+overlap size */
997: for (p = pStart, l = 0; p < pEnd; ++p) {
998: PetscInt bdof, cdof, dof, off, c, cind = 0;
1000: /* Should probably use constrained dofs */
1001: PetscSectionGetDof(section, p, &dof);
1002: PetscSectionGetConstraintDof(section, p, &cdof);
1003: PetscSectionGetConstraintIndices(section, p, &cdofs);
1004: PetscSectionGetOffset(sectionGlobal, p, &off);
1005: /* If you have dofs, and constraints, and they are unequal, we set the blocksize to 1 */
1006: bdof = cdof && (dof-cdof) ? 1 : dof;
1007: if (dof) {
1008: if (bs < 0) {bs = bdof;}
1009: else if (bs != bdof) {bs = 1;}
1010: }
1011: for (c = 0; c < dof; ++c, ++l) {
1012: if ((cind < cdof) && (c == cdofs[cind])) ltog[l] = off < 0 ? off-c : off+c;
1013: else ltog[l] = (off < 0 ? -(off+1) : off) + c;
1014: }
1015: }
1016: /* Must have same blocksize on all procs (some might have no points) */
1017: bsLocal[0] = bs < 0 ? PETSC_MAX_INT : bs; bsLocal[1] = bs;
1018: PetscGlobalMinMaxInt(PetscObjectComm((PetscObject) dm), bsLocal, bsMinMax);
1019: if (bsMinMax[0] != bsMinMax[1]) {bs = 1;}
1020: else {bs = bsMinMax[0];}
1021: bs = bs < 0 ? 1 : bs;
1022: /* Must reduce indices by blocksize */
1023: if (bs > 1) {
1024: for (l = 0, k = 0; l < n; l += bs, ++k) ltog[k] = ltog[l]/bs;
1025: n /= bs;
1026: }
1027: ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)dm), bs, n, ltog, PETSC_OWN_POINTER, &dm->ltogmap);
1028: PetscLogObjectParent((PetscObject)dm, (PetscObject)dm->ltogmap);
1029: } else {
1030: if (!dm->ops->getlocaltoglobalmapping) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM can not create LocalToGlobalMapping");
1031: (*dm->ops->getlocaltoglobalmapping)(dm);
1032: }
1033: }
1034: *ltog = dm->ltogmap;
1035: return(0);
1036: }
1038: /*@
1039: DMGetBlockSize - Gets the inherent block size associated with a DM
1041: Not Collective
1043: Input Parameter:
1044: . dm - the DM with block structure
1046: Output Parameter:
1047: . bs - the block size, 1 implies no exploitable block structure
1049: Level: intermediate
1051: .seealso: ISCreateBlock(), VecSetBlockSize(), MatSetBlockSize(), DMGetLocalToGlobalMapping()
1052: @*/
1053: PetscErrorCode DMGetBlockSize(DM dm,PetscInt *bs)
1054: {
1058: if (dm->bs < 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"DM does not have enough information to provide a block size yet");
1059: *bs = dm->bs;
1060: return(0);
1061: }
1063: /*@
1064: DMCreateInterpolation - Gets interpolation matrix between two DM objects
1066: Collective on DM
1068: Input Parameter:
1069: + dm1 - the DM object
1070: - dm2 - the second, finer DM object
1072: Output Parameter:
1073: + mat - the interpolation
1074: - vec - the scaling (optional)
1076: Level: developer
1078: Notes:
1079: For DMDA objects this only works for "uniform refinement", that is the refined mesh was obtained DMRefine() or the coarse mesh was obtained by
1080: DMCoarsen(). The coordinates set into the DMDA are completely ignored in computing the interpolation.
1082: For DMDA objects you can use this interpolation (more precisely the interpolation from the DMGetCoordinateDM()) to interpolate the mesh coordinate vectors
1083: EXCEPT in the periodic case where it does not make sense since the coordinate vectors are not periodic.
1086: .seealso DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateColoring(), DMCreateMatrix(), DMRefine(), DMCoarsen(), DMCreateRestriction()
1088: @*/
1089: PetscErrorCode DMCreateInterpolation(DM dm1,DM dm2,Mat *mat,Vec *vec)
1090: {
1096: PetscLogEventBegin(DM_CreateInterpolation,dm1,dm2,0,0);
1097: (*dm1->ops->createinterpolation)(dm1,dm2,mat,vec);
1098: PetscLogEventEnd(DM_CreateInterpolation,dm1,dm2,0,0);
1099: return(0);
1100: }
1102: /*@
1103: DMCreateRestriction - Gets restriction matrix between two DM objects
1105: Collective on DM
1107: Input Parameter:
1108: + dm1 - the DM object
1109: - dm2 - the second, finer DM object
1111: Output Parameter:
1112: . mat - the restriction
1115: Level: developer
1117: Notes:
1118: For DMDA objects this only works for "uniform refinement", that is the refined mesh was obtained DMRefine() or the coarse mesh was obtained by
1119: DMCoarsen(). The coordinates set into the DMDA are completely ignored in computing the interpolation.
1122: .seealso DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateColoring(), DMCreateMatrix(), DMRefine(), DMCoarsen(), DMCreateInterpolation()
1124: @*/
1125: PetscErrorCode DMCreateRestriction(DM dm1,DM dm2,Mat *mat)
1126: {
1132: PetscLogEventBegin(DM_CreateRestriction,dm1,dm2,0,0);
1133: if (!dm1->ops->createrestriction) SETERRQ(PetscObjectComm((PetscObject)dm1),PETSC_ERR_SUP,"DMCreateRestriction not implemented for this type");
1134: (*dm1->ops->createrestriction)(dm1,dm2,mat);
1135: PetscLogEventEnd(DM_CreateRestriction,dm1,dm2,0,0);
1136: return(0);
1137: }
1139: /*@
1140: DMCreateInjection - Gets injection matrix between two DM objects
1142: Collective on DM
1144: Input Parameter:
1145: + dm1 - the DM object
1146: - dm2 - the second, finer DM object
1148: Output Parameter:
1149: . mat - the injection
1151: Level: developer
1153: Notes:
1154: For DMDA objects this only works for "uniform refinement", that is the refined mesh was obtained DMRefine() or the coarse mesh was obtained by
1155: DMCoarsen(). The coordinates set into the DMDA are completely ignored in computing the injection.
1157: .seealso DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateColoring(), DMCreateMatrix(), DMCreateInterpolation()
1159: @*/
1160: PetscErrorCode DMCreateInjection(DM dm1,DM dm2,Mat *mat)
1161: {
1167: if (!dm1->ops->getinjection) SETERRQ(PetscObjectComm((PetscObject)dm1),PETSC_ERR_SUP,"DMCreateInjection not implemented for this type");
1168: (*dm1->ops->getinjection)(dm1,dm2,mat);
1169: return(0);
1170: }
1172: /*@
1173: DMCreateMassMatrix - Gets mass matrix between two DM objects, M_ij = \int \phi_i \psi_j
1175: Collective on DM
1177: Input Parameter:
1178: + dm1 - the DM object
1179: - dm2 - the second, finer DM object
1181: Output Parameter:
1182: . mat - the interpolation
1184: Level: developer
1186: .seealso DMCreateMatrix(), DMRefine(), DMCoarsen(), DMCreateRestriction(), DMCreateInterpolation(), DMCreateInjection()
1187: @*/
1188: PetscErrorCode DMCreateMassMatrix(DM dm1, DM dm2, Mat *mat)
1189: {
1195: (*dm1->ops->createmassmatrix)(dm1, dm2, mat);
1196: return(0);
1197: }
1199: /*@
1200: DMCreateColoring - Gets coloring for a DM
1202: Collective on DM
1204: Input Parameter:
1205: + dm - the DM object
1206: - ctype - IS_COLORING_LOCAL or IS_COLORING_GLOBAL
1208: Output Parameter:
1209: . coloring - the coloring
1211: Level: developer
1213: .seealso DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateMatrix(), DMSetMatType()
1215: @*/
1216: PetscErrorCode DMCreateColoring(DM dm,ISColoringType ctype,ISColoring *coloring)
1217: {
1222: if (!dm->ops->getcoloring) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"No coloring for this type of DM yet");
1223: (*dm->ops->getcoloring)(dm,ctype,coloring);
1224: return(0);
1225: }
1227: /*@
1228: DMCreateMatrix - Gets empty Jacobian for a DM
1230: Collective on DM
1232: Input Parameter:
1233: . dm - the DM object
1235: Output Parameter:
1236: . mat - the empty Jacobian
1238: Level: beginner
1240: Notes:
1241: This properly preallocates the number of nonzeros in the sparse matrix so you
1242: do not need to do it yourself.
1244: By default it also sets the nonzero structure and puts in the zero entries. To prevent setting
1245: the nonzero pattern call DMSetMatrixPreallocateOnly()
1247: For structured grid problems, when you call MatView() on this matrix it is displayed using the global natural ordering, NOT in the ordering used
1248: internally by PETSc.
1250: For structured grid problems, in general it is easiest to use MatSetValuesStencil() or MatSetValuesLocal() to put values into the matrix because MatSetValues() requires
1251: the indices for the global numbering for DMDAs which is complicated.
1253: .seealso DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMSetMatType()
1255: @*/
1256: PetscErrorCode DMCreateMatrix(DM dm,Mat *mat)
1257: {
1262: MatInitializePackage();
1265: (*dm->ops->creatematrix)(dm,mat);
1266: return(0);
1267: }
1269: /*@
1270: DMSetMatrixPreallocateOnly - When DMCreateMatrix() is called the matrix will be properly
1271: preallocated but the nonzero structure and zero values will not be set.
1273: Logically Collective on DM
1275: Input Parameter:
1276: + dm - the DM
1277: - only - PETSC_TRUE if only want preallocation
1279: Level: developer
1280: .seealso DMCreateMatrix(), DMSetMatrixStructureOnly()
1281: @*/
1282: PetscErrorCode DMSetMatrixPreallocateOnly(DM dm, PetscBool only)
1283: {
1286: dm->prealloc_only = only;
1287: return(0);
1288: }
1290: /*@
1291: DMSetMatrixStructureOnly - When DMCreateMatrix() is called, the matrix structure will be created
1292: but the array for values will not be allocated.
1294: Logically Collective on DM
1296: Input Parameter:
1297: + dm - the DM
1298: - only - PETSC_TRUE if only want matrix stucture
1300: Level: developer
1301: .seealso DMCreateMatrix(), DMSetMatrixPreallocateOnly()
1302: @*/
1303: PetscErrorCode DMSetMatrixStructureOnly(DM dm, PetscBool only)
1304: {
1307: dm->structure_only = only;
1308: return(0);
1309: }
1311: /*@C
1312: DMGetWorkArray - Gets a work array guaranteed to be at least the input size, restore with DMRestoreWorkArray()
1314: Not Collective
1316: Input Parameters:
1317: + dm - the DM object
1318: . count - The minium size
1319: - dtype - MPI data type, often MPIU_REAL, MPIU_SCALAR, MPIU_INT)
1321: Output Parameter:
1322: . array - the work array
1324: Level: developer
1326: .seealso DMDestroy(), DMCreate()
1327: @*/
1328: PetscErrorCode DMGetWorkArray(DM dm,PetscInt count,MPI_Datatype dtype,void *mem)
1329: {
1331: DMWorkLink link;
1332: PetscMPIInt dsize;
1337: if (dm->workin) {
1338: link = dm->workin;
1339: dm->workin = dm->workin->next;
1340: } else {
1341: PetscNewLog(dm,&link);
1342: }
1343: MPI_Type_size(dtype,&dsize);
1344: if (((size_t)dsize*count) > link->bytes) {
1345: PetscFree(link->mem);
1346: PetscMalloc(dsize*count,&link->mem);
1347: link->bytes = dsize*count;
1348: }
1349: link->next = dm->workout;
1350: dm->workout = link;
1351: *(void**)mem = link->mem;
1352: return(0);
1353: }
1355: /*@C
1356: DMRestoreWorkArray - Restores a work array guaranteed to be at least the input size, restore with DMRestoreWorkArray()
1358: Not Collective
1360: Input Parameters:
1361: + dm - the DM object
1362: . count - The minium size
1363: - dtype - MPI data type, often MPIU_REAL, MPIU_SCALAR, MPIU_INT
1365: Output Parameter:
1366: . array - the work array
1368: Level: developer
1370: Developer Notes:
1371: count and dtype are ignored, they are only needed for DMGetWorkArray()
1372: .seealso DMDestroy(), DMCreate()
1373: @*/
1374: PetscErrorCode DMRestoreWorkArray(DM dm,PetscInt count,MPI_Datatype dtype,void *mem)
1375: {
1376: DMWorkLink *p,link;
1381: for (p=&dm->workout; (link=*p); p=&link->next) {
1382: if (link->mem == *(void**)mem) {
1383: *p = link->next;
1384: link->next = dm->workin;
1385: dm->workin = link;
1386: *(void**)mem = NULL;
1387: return(0);
1388: }
1389: }
1390: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Array was not checked out");
1391: }
1393: PetscErrorCode DMSetNullSpaceConstructor(DM dm, PetscInt field, PetscErrorCode (*nullsp)(DM dm, PetscInt field, MatNullSpace *nullSpace))
1394: {
1397: if (field >= 10) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle %d >= 10 fields", field);
1398: dm->nullspaceConstructors[field] = nullsp;
1399: return(0);
1400: }
1402: PetscErrorCode DMGetNullSpaceConstructor(DM dm, PetscInt field, PetscErrorCode (**nullsp)(DM dm, PetscInt field, MatNullSpace *nullSpace))
1403: {
1406: if (field >= 10) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle %d >= 10 fields", field);
1407: *nullsp = dm->nullspaceConstructors[field];
1408: return(0);
1409: }
1411: /*@C
1412: DMCreateFieldIS - Creates a set of IS objects with the global indices of dofs for each field
1414: Not collective
1416: Input Parameter:
1417: . dm - the DM object
1419: Output Parameters:
1420: + numFields - The number of fields (or NULL if not requested)
1421: . fieldNames - The name for each field (or NULL if not requested)
1422: - fields - The global indices for each field (or NULL if not requested)
1424: Level: intermediate
1426: Notes:
1427: The user is responsible for freeing all requested arrays. In particular, every entry of names should be freed with
1428: PetscFree(), every entry of fields should be destroyed with ISDestroy(), and both arrays should be freed with
1429: PetscFree().
1431: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix()
1432: @*/
1433: PetscErrorCode DMCreateFieldIS(DM dm, PetscInt *numFields, char ***fieldNames, IS **fields)
1434: {
1435: PetscSection section, sectionGlobal;
1440: if (numFields) {
1442: *numFields = 0;
1443: }
1444: if (fieldNames) {
1446: *fieldNames = NULL;
1447: }
1448: if (fields) {
1450: *fields = NULL;
1451: }
1452: DMGetSection(dm, §ion);
1453: if (section) {
1454: PetscInt *fieldSizes, **fieldIndices;
1455: PetscInt nF, f, pStart, pEnd, p;
1457: DMGetGlobalSection(dm, §ionGlobal);
1458: PetscSectionGetNumFields(section, &nF);
1459: PetscMalloc2(nF,&fieldSizes,nF,&fieldIndices);
1460: PetscSectionGetChart(sectionGlobal, &pStart, &pEnd);
1461: for (f = 0; f < nF; ++f) {
1462: fieldSizes[f] = 0;
1463: }
1464: for (p = pStart; p < pEnd; ++p) {
1465: PetscInt gdof;
1467: PetscSectionGetDof(sectionGlobal, p, &gdof);
1468: if (gdof > 0) {
1469: for (f = 0; f < nF; ++f) {
1470: PetscInt fdof, fcdof;
1472: PetscSectionGetFieldDof(section, p, f, &fdof);
1473: PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);
1474: fieldSizes[f] += fdof-fcdof;
1475: }
1476: }
1477: }
1478: for (f = 0; f < nF; ++f) {
1479: PetscMalloc1(fieldSizes[f], &fieldIndices[f]);
1480: fieldSizes[f] = 0;
1481: }
1482: for (p = pStart; p < pEnd; ++p) {
1483: PetscInt gdof, goff;
1485: PetscSectionGetDof(sectionGlobal, p, &gdof);
1486: if (gdof > 0) {
1487: PetscSectionGetOffset(sectionGlobal, p, &goff);
1488: for (f = 0; f < nF; ++f) {
1489: PetscInt fdof, fcdof, fc;
1491: PetscSectionGetFieldDof(section, p, f, &fdof);
1492: PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);
1493: for (fc = 0; fc < fdof-fcdof; ++fc, ++fieldSizes[f]) {
1494: fieldIndices[f][fieldSizes[f]] = goff++;
1495: }
1496: }
1497: }
1498: }
1499: if (numFields) *numFields = nF;
1500: if (fieldNames) {
1501: PetscMalloc1(nF, fieldNames);
1502: for (f = 0; f < nF; ++f) {
1503: const char *fieldName;
1505: PetscSectionGetFieldName(section, f, &fieldName);
1506: PetscStrallocpy(fieldName, (char**) &(*fieldNames)[f]);
1507: }
1508: }
1509: if (fields) {
1510: PetscMalloc1(nF, fields);
1511: for (f = 0; f < nF; ++f) {
1512: ISCreateGeneral(PetscObjectComm((PetscObject)dm), fieldSizes[f], fieldIndices[f], PETSC_OWN_POINTER, &(*fields)[f]);
1513: }
1514: }
1515: PetscFree2(fieldSizes,fieldIndices);
1516: } else if (dm->ops->createfieldis) {
1517: (*dm->ops->createfieldis)(dm, numFields, fieldNames, fields);
1518: }
1519: return(0);
1520: }
1523: /*@C
1524: DMCreateFieldDecomposition - Returns a list of IS objects defining a decomposition of a problem into subproblems
1525: corresponding to different fields: each IS contains the global indices of the dofs of the
1526: corresponding field. The optional list of DMs define the DM for each subproblem.
1527: Generalizes DMCreateFieldIS().
1529: Not collective
1531: Input Parameter:
1532: . dm - the DM object
1534: Output Parameters:
1535: + len - The number of subproblems in the field decomposition (or NULL if not requested)
1536: . namelist - The name for each field (or NULL if not requested)
1537: . islist - The global indices for each field (or NULL if not requested)
1538: - dmlist - The DMs for each field subproblem (or NULL, if not requested; if NULL is returned, no DMs are defined)
1540: Level: intermediate
1542: Notes:
1543: The user is responsible for freeing all requested arrays. In particular, every entry of names should be freed with
1544: PetscFree(), every entry of is should be destroyed with ISDestroy(), every entry of dm should be destroyed with DMDestroy(),
1545: and all of the arrays should be freed with PetscFree().
1547: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateFieldIS()
1548: @*/
1549: PetscErrorCode DMCreateFieldDecomposition(DM dm, PetscInt *len, char ***namelist, IS **islist, DM **dmlist)
1550: {
1555: if (len) {
1557: *len = 0;
1558: }
1559: if (namelist) {
1561: *namelist = 0;
1562: }
1563: if (islist) {
1565: *islist = 0;
1566: }
1567: if (dmlist) {
1569: *dmlist = 0;
1570: }
1571: /*
1572: Is it a good idea to apply the following check across all impls?
1573: Perhaps some impls can have a well-defined decomposition before DMSetUp?
1574: This, however, follows the general principle that accessors are not well-behaved until the object is set up.
1575: */
1576: if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE, "Decomposition defined only after DMSetUp");
1577: if (!dm->ops->createfielddecomposition) {
1578: PetscSection section;
1579: PetscInt numFields, f;
1581: DMGetSection(dm, §ion);
1582: if (section) {PetscSectionGetNumFields(section, &numFields);}
1583: if (section && numFields && dm->ops->createsubdm) {
1584: if (len) *len = numFields;
1585: if (namelist) {PetscMalloc1(numFields,namelist);}
1586: if (islist) {PetscMalloc1(numFields,islist);}
1587: if (dmlist) {PetscMalloc1(numFields,dmlist);}
1588: for (f = 0; f < numFields; ++f) {
1589: const char *fieldName;
1591: DMCreateSubDM(dm, 1, &f, islist ? &(*islist)[f] : NULL, dmlist ? &(*dmlist)[f] : NULL);
1592: if (namelist) {
1593: PetscSectionGetFieldName(section, f, &fieldName);
1594: PetscStrallocpy(fieldName, (char**) &(*namelist)[f]);
1595: }
1596: }
1597: } else {
1598: DMCreateFieldIS(dm, len, namelist, islist);
1599: /* By default there are no DMs associated with subproblems. */
1600: if (dmlist) *dmlist = NULL;
1601: }
1602: } else {
1603: (*dm->ops->createfielddecomposition)(dm,len,namelist,islist,dmlist);
1604: }
1605: return(0);
1606: }
1608: /*@
1609: DMCreateSubDM - Returns an IS and DM encapsulating a subproblem defined by the fields passed in.
1610: The fields are defined by DMCreateFieldIS().
1612: Not collective
1614: Input Parameters:
1615: + dm - The DM object
1616: . numFields - The number of fields in this subproblem
1617: - fields - The field numbers of the selected fields
1619: Output Parameters:
1620: + is - The global indices for the subproblem
1621: - subdm - The DM for the subproblem
1623: Level: intermediate
1625: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateFieldIS()
1626: @*/
1627: PetscErrorCode DMCreateSubDM(DM dm, PetscInt numFields, const PetscInt fields[], IS *is, DM *subdm)
1628: {
1636: if (dm->ops->createsubdm) {
1637: (*dm->ops->createsubdm)(dm, numFields, fields, is, subdm);
1638: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "This type has no DMCreateSubDM implementation defined");
1639: return(0);
1640: }
1642: /*@C
1643: DMCreateSuperDM - Returns an arrays of ISes and DM encapsulating a superproblem defined by the DMs passed in.
1645: Not collective
1647: Input Parameter:
1648: + dms - The DM objects
1649: - len - The number of DMs
1651: Output Parameters:
1652: + is - The global indices for the subproblem, or NULL
1653: - superdm - The DM for the superproblem
1655: Level: intermediate
1657: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateFieldIS()
1658: @*/
1659: PetscErrorCode DMCreateSuperDM(DM dms[], PetscInt len, IS **is, DM *superdm)
1660: {
1661: PetscInt i;
1669: if (len < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of DMs must be nonnegative: %D", len);
1670: if (len) {
1671: if (dms[0]->ops->createsuperdm) {(*dms[0]->ops->createsuperdm)(dms, len, is, superdm);}
1672: else SETERRQ(PetscObjectComm((PetscObject) dms[0]), PETSC_ERR_SUP, "This type has no DMCreateSuperDM implementation defined");
1673: }
1674: return(0);
1675: }
1678: /*@C
1679: DMCreateDomainDecomposition - Returns lists of IS objects defining a decomposition of a problem into subproblems
1680: corresponding to restrictions to pairs nested subdomains: each IS contains the global
1681: indices of the dofs of the corresponding subdomains. The inner subdomains conceptually
1682: define a nonoverlapping covering, while outer subdomains can overlap.
1683: The optional list of DMs define the DM for each subproblem.
1685: Not collective
1687: Input Parameter:
1688: . dm - the DM object
1690: Output Parameters:
1691: + len - The number of subproblems in the domain decomposition (or NULL if not requested)
1692: . namelist - The name for each subdomain (or NULL if not requested)
1693: . innerislist - The global indices for each inner subdomain (or NULL, if not requested)
1694: . outerislist - The global indices for each outer subdomain (or NULL, if not requested)
1695: - dmlist - The DMs for each subdomain subproblem (or NULL, if not requested; if NULL is returned, no DMs are defined)
1697: Level: intermediate
1699: Notes:
1700: The user is responsible for freeing all requested arrays. In particular, every entry of names should be freed with
1701: PetscFree(), every entry of is should be destroyed with ISDestroy(), every entry of dm should be destroyed with DMDestroy(),
1702: and all of the arrays should be freed with PetscFree().
1704: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateDomainDecompositionDM(), DMCreateFieldDecomposition()
1705: @*/
1706: PetscErrorCode DMCreateDomainDecomposition(DM dm, PetscInt *len, char ***namelist, IS **innerislist, IS **outerislist, DM **dmlist)
1707: {
1708: PetscErrorCode ierr;
1709: DMSubDomainHookLink link;
1710: PetscInt i,l;
1719: /*
1720: Is it a good idea to apply the following check across all impls?
1721: Perhaps some impls can have a well-defined decomposition before DMSetUp?
1722: This, however, follows the general principle that accessors are not well-behaved until the object is set up.
1723: */
1724: if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE, "Decomposition defined only after DMSetUp");
1725: if (dm->ops->createdomaindecomposition) {
1726: (*dm->ops->createdomaindecomposition)(dm,&l,namelist,innerislist,outerislist,dmlist);
1727: /* copy subdomain hooks and context over to the subdomain DMs */
1728: if (dmlist && *dmlist) {
1729: for (i = 0; i < l; i++) {
1730: for (link=dm->subdomainhook; link; link=link->next) {
1731: if (link->ddhook) {(*link->ddhook)(dm,(*dmlist)[i],link->ctx);}
1732: }
1733: if (dm->ctx) (*dmlist)[i]->ctx = dm->ctx;
1734: }
1735: }
1736: if (len) *len = l;
1737: }
1738: return(0);
1739: }
1742: /*@C
1743: DMCreateDomainDecompositionScatters - Returns scatters to the subdomain vectors from the global vector
1745: Not collective
1747: Input Parameters:
1748: + dm - the DM object
1749: . n - the number of subdomain scatters
1750: - subdms - the local subdomains
1752: Output Parameters:
1753: + n - the number of scatters returned
1754: . iscat - scatter from global vector to nonoverlapping global vector entries on subdomain
1755: . oscat - scatter from global vector to overlapping global vector entries on subdomain
1756: - gscat - scatter from global vector to local vector on subdomain (fills in ghosts)
1758: Notes:
1759: This is an alternative to the iis and ois arguments in DMCreateDomainDecomposition that allow for the solution
1760: of general nonlinear problems with overlapping subdomain methods. While merely having index sets that enable subsets
1761: of the residual equations to be created is fine for linear problems, nonlinear problems require local assembly of
1762: solution and residual data.
1764: Level: developer
1766: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateFieldIS()
1767: @*/
1768: PetscErrorCode DMCreateDomainDecompositionScatters(DM dm,PetscInt n,DM *subdms,VecScatter **iscat,VecScatter **oscat,VecScatter **gscat)
1769: {
1775: if (dm->ops->createddscatters) {
1776: (*dm->ops->createddscatters)(dm,n,subdms,iscat,oscat,gscat);
1777: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "This type has no DMCreateDomainDecompositionScatter implementation defined");
1778: return(0);
1779: }
1781: /*@
1782: DMRefine - Refines a DM object
1784: Collective on DM
1786: Input Parameter:
1787: + dm - the DM object
1788: - comm - the communicator to contain the new DM object (or MPI_COMM_NULL)
1790: Output Parameter:
1791: . dmf - the refined DM, or NULL
1793: Note: If no refinement was done, the return value is NULL
1795: Level: developer
1797: .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation()
1798: @*/
1799: PetscErrorCode DMRefine(DM dm,MPI_Comm comm,DM *dmf)
1800: {
1801: PetscErrorCode ierr;
1802: DMRefineHookLink link;
1806: PetscLogEventBegin(DM_Refine,dm,0,0,0);
1807: if (!dm->ops->refine) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"This DM cannot refine");
1808: (*dm->ops->refine)(dm,comm,dmf);
1809: if (*dmf) {
1810: (*dmf)->ops->creatematrix = dm->ops->creatematrix;
1812: PetscObjectCopyFortranFunctionPointers((PetscObject)dm,(PetscObject)*dmf);
1814: (*dmf)->ctx = dm->ctx;
1815: (*dmf)->leveldown = dm->leveldown;
1816: (*dmf)->levelup = dm->levelup + 1;
1818: DMSetMatType(*dmf,dm->mattype);
1819: for (link=dm->refinehook; link; link=link->next) {
1820: if (link->refinehook) {
1821: (*link->refinehook)(dm,*dmf,link->ctx);
1822: }
1823: }
1824: }
1825: PetscLogEventEnd(DM_Refine,dm,0,0,0);
1826: return(0);
1827: }
1829: /*@C
1830: DMRefineHookAdd - adds a callback to be run when interpolating a nonlinear problem to a finer grid
1832: Logically Collective
1834: Input Arguments:
1835: + coarse - nonlinear solver context on which to run a hook when restricting to a coarser level
1836: . refinehook - function to run when setting up a coarser level
1837: . interphook - function to run to update data on finer levels (once per SNESSolve())
1838: - ctx - [optional] user-defined context for provide data for the hooks (may be NULL)
1840: Calling sequence of refinehook:
1841: $ refinehook(DM coarse,DM fine,void *ctx);
1843: + coarse - coarse level DM
1844: . fine - fine level DM to interpolate problem to
1845: - ctx - optional user-defined function context
1847: Calling sequence for interphook:
1848: $ interphook(DM coarse,Mat interp,DM fine,void *ctx)
1850: + coarse - coarse level DM
1851: . interp - matrix interpolating a coarse-level solution to the finer grid
1852: . fine - fine level DM to update
1853: - ctx - optional user-defined function context
1855: Level: advanced
1857: Notes:
1858: This function is only needed if auxiliary data needs to be passed to fine grids while grid sequencing
1860: If this function is called multiple times, the hooks will be run in the order they are added.
1862: This function is currently not available from Fortran.
1864: .seealso: DMCoarsenHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
1865: @*/
1866: PetscErrorCode DMRefineHookAdd(DM coarse,PetscErrorCode (*refinehook)(DM,DM,void*),PetscErrorCode (*interphook)(DM,Mat,DM,void*),void *ctx)
1867: {
1868: PetscErrorCode ierr;
1869: DMRefineHookLink link,*p;
1873: for (p=&coarse->refinehook; *p; p=&(*p)->next) { /* Scan to the end of the current list of hooks */
1874: if ((*p)->refinehook == refinehook && (*p)->interphook == interphook && (*p)->ctx == ctx) return(0);
1875: }
1876: PetscNew(&link);
1877: link->refinehook = refinehook;
1878: link->interphook = interphook;
1879: link->ctx = ctx;
1880: link->next = NULL;
1881: *p = link;
1882: return(0);
1883: }
1885: /*@C
1886: DMRefineHookRemove - remove a callback from the list of hooks to be run when interpolating a nonlinear problem to a finer grid
1888: Logically Collective
1890: Input Arguments:
1891: + coarse - nonlinear solver context on which to run a hook when restricting to a coarser level
1892: . refinehook - function to run when setting up a coarser level
1893: . interphook - function to run to update data on finer levels (once per SNESSolve())
1894: - ctx - [optional] user-defined context for provide data for the hooks (may be NULL)
1896: Level: advanced
1898: Notes:
1899: This function does nothing if the hook is not in the list.
1901: This function is currently not available from Fortran.
1903: .seealso: DMCoarsenHookRemove(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
1904: @*/
1905: PetscErrorCode DMRefineHookRemove(DM coarse,PetscErrorCode (*refinehook)(DM,DM,void*),PetscErrorCode (*interphook)(DM,Mat,DM,void*),void *ctx)
1906: {
1907: PetscErrorCode ierr;
1908: DMRefineHookLink link,*p;
1912: for (p=&coarse->refinehook; *p; p=&(*p)->next) { /* Search the list of current hooks */
1913: if ((*p)->refinehook == refinehook && (*p)->interphook == interphook && (*p)->ctx == ctx) {
1914: link = *p;
1915: *p = link->next;
1916: PetscFree(link);
1917: break;
1918: }
1919: }
1920: return(0);
1921: }
1923: /*@
1924: DMInterpolate - interpolates user-defined problem data to a finer DM by running hooks registered by DMRefineHookAdd()
1926: Collective if any hooks are
1928: Input Arguments:
1929: + coarse - coarser DM to use as a base
1930: . interp - interpolation matrix, apply using MatInterpolate()
1931: - fine - finer DM to update
1933: Level: developer
1935: .seealso: DMRefineHookAdd(), MatInterpolate()
1936: @*/
1937: PetscErrorCode DMInterpolate(DM coarse,Mat interp,DM fine)
1938: {
1939: PetscErrorCode ierr;
1940: DMRefineHookLink link;
1943: for (link=fine->refinehook; link; link=link->next) {
1944: if (link->interphook) {
1945: (*link->interphook)(coarse,interp,fine,link->ctx);
1946: }
1947: }
1948: return(0);
1949: }
1951: /*@
1952: DMGetRefineLevel - Get's the number of refinements that have generated this DM.
1954: Not Collective
1956: Input Parameter:
1957: . dm - the DM object
1959: Output Parameter:
1960: . level - number of refinements
1962: Level: developer
1964: .seealso DMCoarsen(), DMGetCoarsenLevel(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation()
1966: @*/
1967: PetscErrorCode DMGetRefineLevel(DM dm,PetscInt *level)
1968: {
1971: *level = dm->levelup;
1972: return(0);
1973: }
1975: /*@
1976: DMSetRefineLevel - Set's the number of refinements that have generated this DM.
1978: Not Collective
1980: Input Parameter:
1981: + dm - the DM object
1982: - level - number of refinements
1984: Level: advanced
1986: Notes:
1987: This value is used by PCMG to determine how many multigrid levels to use
1989: .seealso DMCoarsen(), DMGetCoarsenLevel(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation()
1991: @*/
1992: PetscErrorCode DMSetRefineLevel(DM dm,PetscInt level)
1993: {
1996: dm->levelup = level;
1997: return(0);
1998: }
2000: /*@C
2001: DMGlobalToLocalHookAdd - adds a callback to be run when global to local is called
2003: Logically Collective
2005: Input Arguments:
2006: + dm - the DM
2007: . beginhook - function to run at the beginning of DMGlobalToLocalBegin()
2008: . endhook - function to run after DMGlobalToLocalEnd() has completed
2009: - ctx - [optional] user-defined context for provide data for the hooks (may be NULL)
2011: Calling sequence for beginhook:
2012: $ beginhook(DM fine,VecScatter out,VecScatter in,DM coarse,void *ctx)
2014: + dm - global DM
2015: . g - global vector
2016: . mode - mode
2017: . l - local vector
2018: - ctx - optional user-defined function context
2021: Calling sequence for endhook:
2022: $ endhook(DM fine,VecScatter out,VecScatter in,DM coarse,void *ctx)
2024: + global - global DM
2025: - ctx - optional user-defined function context
2027: Level: advanced
2029: .seealso: DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2030: @*/
2031: PetscErrorCode DMGlobalToLocalHookAdd(DM dm,PetscErrorCode (*beginhook)(DM,Vec,InsertMode,Vec,void*),PetscErrorCode (*endhook)(DM,Vec,InsertMode,Vec,void*),void *ctx)
2032: {
2033: PetscErrorCode ierr;
2034: DMGlobalToLocalHookLink link,*p;
2038: for (p=&dm->gtolhook; *p; p=&(*p)->next) {} /* Scan to the end of the current list of hooks */
2039: PetscNew(&link);
2040: link->beginhook = beginhook;
2041: link->endhook = endhook;
2042: link->ctx = ctx;
2043: link->next = NULL;
2044: *p = link;
2045: return(0);
2046: }
2048: static PetscErrorCode DMGlobalToLocalHook_Constraints(DM dm, Vec g, InsertMode mode, Vec l, void *ctx)
2049: {
2050: Mat cMat;
2051: Vec cVec;
2052: PetscSection section, cSec;
2053: PetscInt pStart, pEnd, p, dof;
2058: DMGetDefaultConstraints(dm,&cSec,&cMat);
2059: if (cMat && (mode == INSERT_VALUES || mode == INSERT_ALL_VALUES || mode == INSERT_BC_VALUES)) {
2060: PetscInt nRows;
2062: MatGetSize(cMat,&nRows,NULL);
2063: if (nRows <= 0) return(0);
2064: DMGetSection(dm,§ion);
2065: MatCreateVecs(cMat,NULL,&cVec);
2066: MatMult(cMat,l,cVec);
2067: PetscSectionGetChart(cSec,&pStart,&pEnd);
2068: for (p = pStart; p < pEnd; p++) {
2069: PetscSectionGetDof(cSec,p,&dof);
2070: if (dof) {
2071: PetscScalar *vals;
2072: VecGetValuesSection(cVec,cSec,p,&vals);
2073: VecSetValuesSection(l,section,p,vals,INSERT_ALL_VALUES);
2074: }
2075: }
2076: VecDestroy(&cVec);
2077: }
2078: return(0);
2079: }
2081: /*@
2082: DMGlobalToLocalBegin - Begins updating local vectors from global vector
2084: Neighbor-wise Collective on DM
2086: Input Parameters:
2087: + dm - the DM object
2088: . g - the global vector
2089: . mode - INSERT_VALUES or ADD_VALUES
2090: - l - the local vector
2093: Level: beginner
2095: .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMGlobalToLocalEnd(), DMLocalToGlobalBegin()
2097: @*/
2098: PetscErrorCode DMGlobalToLocalBegin(DM dm,Vec g,InsertMode mode,Vec l)
2099: {
2100: PetscSF sf;
2101: PetscErrorCode ierr;
2102: DMGlobalToLocalHookLink link;
2106: for (link=dm->gtolhook; link; link=link->next) {
2107: if (link->beginhook) {
2108: (*link->beginhook)(dm,g,mode,l,link->ctx);
2109: }
2110: }
2111: DMGetDefaultSF(dm, &sf);
2112: if (sf) {
2113: const PetscScalar *gArray;
2114: PetscScalar *lArray;
2116: if (mode == ADD_VALUES) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode);
2117: VecGetArray(l, &lArray);
2118: VecGetArrayRead(g, &gArray);
2119: PetscSFBcastBegin(sf, MPIU_SCALAR, gArray, lArray);
2120: VecRestoreArray(l, &lArray);
2121: VecRestoreArrayRead(g, &gArray);
2122: } else {
2123: (*dm->ops->globaltolocalbegin)(dm,g,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),l);
2124: }
2125: return(0);
2126: }
2128: /*@
2129: DMGlobalToLocalEnd - Ends updating local vectors from global vector
2131: Neighbor-wise Collective on DM
2133: Input Parameters:
2134: + dm - the DM object
2135: . g - the global vector
2136: . mode - INSERT_VALUES or ADD_VALUES
2137: - l - the local vector
2140: Level: beginner
2142: .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMGlobalToLocalEnd(), DMLocalToGlobalBegin()
2144: @*/
2145: PetscErrorCode DMGlobalToLocalEnd(DM dm,Vec g,InsertMode mode,Vec l)
2146: {
2147: PetscSF sf;
2148: PetscErrorCode ierr;
2149: const PetscScalar *gArray;
2150: PetscScalar *lArray;
2151: DMGlobalToLocalHookLink link;
2155: DMGetDefaultSF(dm, &sf);
2156: if (sf) {
2157: if (mode == ADD_VALUES) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode);
2159: VecGetArray(l, &lArray);
2160: VecGetArrayRead(g, &gArray);
2161: PetscSFBcastEnd(sf, MPIU_SCALAR, gArray, lArray);
2162: VecRestoreArray(l, &lArray);
2163: VecRestoreArrayRead(g, &gArray);
2164: } else {
2165: (*dm->ops->globaltolocalend)(dm,g,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),l);
2166: }
2167: DMGlobalToLocalHook_Constraints(dm,g,mode,l,NULL);
2168: for (link=dm->gtolhook; link; link=link->next) {
2169: if (link->endhook) {(*link->endhook)(dm,g,mode,l,link->ctx);}
2170: }
2171: return(0);
2172: }
2174: /*@C
2175: DMLocalToGlobalHookAdd - adds a callback to be run when a local to global is called
2177: Logically Collective
2179: Input Arguments:
2180: + dm - the DM
2181: . beginhook - function to run at the beginning of DMLocalToGlobalBegin()
2182: . endhook - function to run after DMLocalToGlobalEnd() has completed
2183: - ctx - [optional] user-defined context for provide data for the hooks (may be NULL)
2185: Calling sequence for beginhook:
2186: $ beginhook(DM fine,Vec l,InsertMode mode,Vec g,void *ctx)
2188: + dm - global DM
2189: . l - local vector
2190: . mode - mode
2191: . g - global vector
2192: - ctx - optional user-defined function context
2195: Calling sequence for endhook:
2196: $ endhook(DM fine,Vec l,InsertMode mode,Vec g,void *ctx)
2198: + global - global DM
2199: . l - local vector
2200: . mode - mode
2201: . g - global vector
2202: - ctx - optional user-defined function context
2204: Level: advanced
2206: .seealso: DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2207: @*/
2208: PetscErrorCode DMLocalToGlobalHookAdd(DM dm,PetscErrorCode (*beginhook)(DM,Vec,InsertMode,Vec,void*),PetscErrorCode (*endhook)(DM,Vec,InsertMode,Vec,void*),void *ctx)
2209: {
2210: PetscErrorCode ierr;
2211: DMLocalToGlobalHookLink link,*p;
2215: for (p=&dm->ltoghook; *p; p=&(*p)->next) {} /* Scan to the end of the current list of hooks */
2216: PetscNew(&link);
2217: link->beginhook = beginhook;
2218: link->endhook = endhook;
2219: link->ctx = ctx;
2220: link->next = NULL;
2221: *p = link;
2222: return(0);
2223: }
2225: static PetscErrorCode DMLocalToGlobalHook_Constraints(DM dm, Vec l, InsertMode mode, Vec g, void *ctx)
2226: {
2227: Mat cMat;
2228: Vec cVec;
2229: PetscSection section, cSec;
2230: PetscInt pStart, pEnd, p, dof;
2235: DMGetDefaultConstraints(dm,&cSec,&cMat);
2236: if (cMat && (mode == ADD_VALUES || mode == ADD_ALL_VALUES || mode == ADD_BC_VALUES)) {
2237: PetscInt nRows;
2239: MatGetSize(cMat,&nRows,NULL);
2240: if (nRows <= 0) return(0);
2241: DMGetSection(dm,§ion);
2242: MatCreateVecs(cMat,NULL,&cVec);
2243: PetscSectionGetChart(cSec,&pStart,&pEnd);
2244: for (p = pStart; p < pEnd; p++) {
2245: PetscSectionGetDof(cSec,p,&dof);
2246: if (dof) {
2247: PetscInt d;
2248: PetscScalar *vals;
2249: VecGetValuesSection(l,section,p,&vals);
2250: VecSetValuesSection(cVec,cSec,p,vals,mode);
2251: /* for this to be the true transpose, we have to zero the values that
2252: * we just extracted */
2253: for (d = 0; d < dof; d++) {
2254: vals[d] = 0.;
2255: }
2256: }
2257: }
2258: MatMultTransposeAdd(cMat,cVec,l,l);
2259: VecDestroy(&cVec);
2260: }
2261: return(0);
2262: }
2264: /*@
2265: DMLocalToGlobalBegin - updates global vectors from local vectors
2267: Neighbor-wise Collective on DM
2269: Input Parameters:
2270: + dm - the DM object
2271: . l - the local vector
2272: . mode - if INSERT_VALUES then no parallel communication is used, if ADD_VALUES then all ghost points from the same base point accumulate into that base point.
2273: - g - the global vector
2275: Notes:
2276: In the ADD_VALUES case you normally would zero the receiving vector before beginning this operation.
2277: INSERT_VALUES is not supported for DMDA, in that case simply compute the values directly into a global vector instead of a local one.
2279: Level: beginner
2281: .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMGlobalToLocalEnd(), DMGlobalToLocalBegin()
2283: @*/
2284: PetscErrorCode DMLocalToGlobalBegin(DM dm,Vec l,InsertMode mode,Vec g)
2285: {
2286: PetscSF sf;
2287: PetscSection s, gs;
2288: DMLocalToGlobalHookLink link;
2289: const PetscScalar *lArray;
2290: PetscScalar *gArray;
2291: PetscBool isInsert;
2292: PetscErrorCode ierr;
2296: for (link=dm->ltoghook; link; link=link->next) {
2297: if (link->beginhook) {
2298: (*link->beginhook)(dm,l,mode,g,link->ctx);
2299: }
2300: }
2301: DMLocalToGlobalHook_Constraints(dm,l,mode,g,NULL);
2302: DMGetDefaultSF(dm, &sf);
2303: DMGetSection(dm, &s);
2304: switch (mode) {
2305: case INSERT_VALUES:
2306: case INSERT_ALL_VALUES:
2307: case INSERT_BC_VALUES:
2308: isInsert = PETSC_TRUE; break;
2309: case ADD_VALUES:
2310: case ADD_ALL_VALUES:
2311: case ADD_BC_VALUES:
2312: isInsert = PETSC_FALSE; break;
2313: default:
2314: SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode);
2315: }
2316: if (sf && !isInsert) {
2317: VecGetArrayRead(l, &lArray);
2318: VecGetArray(g, &gArray);
2319: PetscSFReduceBegin(sf, MPIU_SCALAR, lArray, gArray, MPIU_SUM);
2320: VecRestoreArrayRead(l, &lArray);
2321: VecRestoreArray(g, &gArray);
2322: } else if (s && isInsert) {
2323: PetscInt gStart, pStart, pEnd, p;
2325: DMGetGlobalSection(dm, &gs);
2326: PetscSectionGetChart(s, &pStart, &pEnd);
2327: VecGetOwnershipRange(g, &gStart, NULL);
2328: VecGetArrayRead(l, &lArray);
2329: VecGetArray(g, &gArray);
2330: for (p = pStart; p < pEnd; ++p) {
2331: PetscInt dof, gdof, cdof, gcdof, off, goff, d, e;
2333: PetscSectionGetDof(s, p, &dof);
2334: PetscSectionGetDof(gs, p, &gdof);
2335: PetscSectionGetConstraintDof(s, p, &cdof);
2336: PetscSectionGetConstraintDof(gs, p, &gcdof);
2337: PetscSectionGetOffset(s, p, &off);
2338: PetscSectionGetOffset(gs, p, &goff);
2339: /* Ignore off-process data and points with no global data */
2340: if (!gdof || goff < 0) continue;
2341: if (dof != gdof) SETERRQ5(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Inconsistent sizes, p: %d dof: %d gdof: %d cdof: %d gcdof: %d", p, dof, gdof, cdof, gcdof);
2342: /* If no constraints are enforced in the global vector */
2343: if (!gcdof) {
2344: for (d = 0; d < dof; ++d) gArray[goff-gStart+d] = lArray[off+d];
2345: /* If constraints are enforced in the global vector */
2346: } else if (cdof == gcdof) {
2347: const PetscInt *cdofs;
2348: PetscInt cind = 0;
2350: PetscSectionGetConstraintIndices(s, p, &cdofs);
2351: for (d = 0, e = 0; d < dof; ++d) {
2352: if ((cind < cdof) && (d == cdofs[cind])) {++cind; continue;}
2353: gArray[goff-gStart+e++] = lArray[off+d];
2354: }
2355: } else SETERRQ5(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Inconsistent sizes, p: %d dof: %d gdof: %d cdof: %d gcdof: %d", p, dof, gdof, cdof, gcdof);
2356: }
2357: VecRestoreArrayRead(l, &lArray);
2358: VecRestoreArray(g, &gArray);
2359: } else {
2360: (*dm->ops->localtoglobalbegin)(dm,l,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),g);
2361: }
2362: return(0);
2363: }
2365: /*@
2366: DMLocalToGlobalEnd - updates global vectors from local vectors
2368: Neighbor-wise Collective on DM
2370: Input Parameters:
2371: + dm - the DM object
2372: . l - the local vector
2373: . mode - INSERT_VALUES or ADD_VALUES
2374: - g - the global vector
2377: Level: beginner
2379: .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMGlobalToLocalEnd(), DMGlobalToLocalEnd()
2381: @*/
2382: PetscErrorCode DMLocalToGlobalEnd(DM dm,Vec l,InsertMode mode,Vec g)
2383: {
2384: PetscSF sf;
2385: PetscSection s;
2386: DMLocalToGlobalHookLink link;
2387: PetscBool isInsert;
2388: PetscErrorCode ierr;
2392: DMGetDefaultSF(dm, &sf);
2393: DMGetSection(dm, &s);
2394: switch (mode) {
2395: case INSERT_VALUES:
2396: case INSERT_ALL_VALUES:
2397: isInsert = PETSC_TRUE; break;
2398: case ADD_VALUES:
2399: case ADD_ALL_VALUES:
2400: isInsert = PETSC_FALSE; break;
2401: default:
2402: SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode);
2403: }
2404: if (sf && !isInsert) {
2405: const PetscScalar *lArray;
2406: PetscScalar *gArray;
2408: VecGetArrayRead(l, &lArray);
2409: VecGetArray(g, &gArray);
2410: PetscSFReduceEnd(sf, MPIU_SCALAR, lArray, gArray, MPIU_SUM);
2411: VecRestoreArrayRead(l, &lArray);
2412: VecRestoreArray(g, &gArray);
2413: } else if (s && isInsert) {
2414: } else {
2415: (*dm->ops->localtoglobalend)(dm,l,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),g);
2416: }
2417: for (link=dm->ltoghook; link; link=link->next) {
2418: if (link->endhook) {(*link->endhook)(dm,g,mode,l,link->ctx);}
2419: }
2420: return(0);
2421: }
2423: /*@
2424: DMLocalToLocalBegin - Maps from a local vector (including ghost points
2425: that contain irrelevant values) to another local vector where the ghost
2426: points in the second are set correctly. Must be followed by DMLocalToLocalEnd().
2428: Neighbor-wise Collective on DM and Vec
2430: Input Parameters:
2431: + dm - the DM object
2432: . g - the original local vector
2433: - mode - one of INSERT_VALUES or ADD_VALUES
2435: Output Parameter:
2436: . l - the local vector with correct ghost values
2438: Level: intermediate
2440: Notes:
2441: The local vectors used here need not be the same as those
2442: obtained from DMCreateLocalVector(), BUT they
2443: must have the same parallel data layout; they could, for example, be
2444: obtained with VecDuplicate() from the DM originating vectors.
2446: .keywords: DM, local-to-local, begin
2447: .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateLocalVector(), DMCreateGlobalVector(), DMCreateInterpolation(), DMLocalToLocalEnd(), DMGlobalToLocalEnd(), DMLocalToGlobalBegin()
2449: @*/
2450: PetscErrorCode DMLocalToLocalBegin(DM dm,Vec g,InsertMode mode,Vec l)
2451: {
2452: PetscErrorCode ierr;
2456: if (!dm->ops->localtolocalbegin) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"This DM does not support local to local maps");
2457: (*dm->ops->localtolocalbegin)(dm,g,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),l);
2458: return(0);
2459: }
2461: /*@
2462: DMLocalToLocalEnd - Maps from a local vector (including ghost points
2463: that contain irrelevant values) to another local vector where the ghost
2464: points in the second are set correctly. Must be preceded by DMLocalToLocalBegin().
2466: Neighbor-wise Collective on DM and Vec
2468: Input Parameters:
2469: + da - the DM object
2470: . g - the original local vector
2471: - mode - one of INSERT_VALUES or ADD_VALUES
2473: Output Parameter:
2474: . l - the local vector with correct ghost values
2476: Level: intermediate
2478: Notes:
2479: The local vectors used here need not be the same as those
2480: obtained from DMCreateLocalVector(), BUT they
2481: must have the same parallel data layout; they could, for example, be
2482: obtained with VecDuplicate() from the DM originating vectors.
2484: .keywords: DM, local-to-local, end
2485: .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateLocalVector(), DMCreateGlobalVector(), DMCreateInterpolation(), DMLocalToLocalBegin(), DMGlobalToLocalEnd(), DMLocalToGlobalBegin()
2487: @*/
2488: PetscErrorCode DMLocalToLocalEnd(DM dm,Vec g,InsertMode mode,Vec l)
2489: {
2490: PetscErrorCode ierr;
2494: if (!dm->ops->localtolocalend) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"This DM does not support local to local maps");
2495: (*dm->ops->localtolocalend)(dm,g,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),l);
2496: return(0);
2497: }
2500: /*@
2501: DMCoarsen - Coarsens a DM object
2503: Collective on DM
2505: Input Parameter:
2506: + dm - the DM object
2507: - comm - the communicator to contain the new DM object (or MPI_COMM_NULL)
2509: Output Parameter:
2510: . dmc - the coarsened DM
2512: Level: developer
2514: .seealso DMRefine(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation()
2516: @*/
2517: PetscErrorCode DMCoarsen(DM dm, MPI_Comm comm, DM *dmc)
2518: {
2519: PetscErrorCode ierr;
2520: DMCoarsenHookLink link;
2524: if (!dm->ops->coarsen) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"This DM cannot coarsen");
2525: PetscLogEventBegin(DM_Coarsen,dm,0,0,0);
2526: (*dm->ops->coarsen)(dm, comm, dmc);
2527: if (!(*dmc)) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "NULL coarse mesh produced");
2528: DMSetCoarseDM(dm,*dmc);
2529: (*dmc)->ops->creatematrix = dm->ops->creatematrix;
2530: PetscObjectCopyFortranFunctionPointers((PetscObject)dm,(PetscObject)*dmc);
2531: (*dmc)->ctx = dm->ctx;
2532: (*dmc)->levelup = dm->levelup;
2533: (*dmc)->leveldown = dm->leveldown + 1;
2534: DMSetMatType(*dmc,dm->mattype);
2535: for (link=dm->coarsenhook; link; link=link->next) {
2536: if (link->coarsenhook) {(*link->coarsenhook)(dm,*dmc,link->ctx);}
2537: }
2538: PetscLogEventEnd(DM_Coarsen,dm,0,0,0);
2539: return(0);
2540: }
2542: /*@C
2543: DMCoarsenHookAdd - adds a callback to be run when restricting a nonlinear problem to the coarse grid
2545: Logically Collective
2547: Input Arguments:
2548: + fine - nonlinear solver context on which to run a hook when restricting to a coarser level
2549: . coarsenhook - function to run when setting up a coarser level
2550: . restricthook - function to run to update data on coarser levels (once per SNESSolve())
2551: - ctx - [optional] user-defined context for provide data for the hooks (may be NULL)
2553: Calling sequence of coarsenhook:
2554: $ coarsenhook(DM fine,DM coarse,void *ctx);
2556: + fine - fine level DM
2557: . coarse - coarse level DM to restrict problem to
2558: - ctx - optional user-defined function context
2560: Calling sequence for restricthook:
2561: $ restricthook(DM fine,Mat mrestrict,Vec rscale,Mat inject,DM coarse,void *ctx)
2563: + fine - fine level DM
2564: . mrestrict - matrix restricting a fine-level solution to the coarse grid
2565: . rscale - scaling vector for restriction
2566: . inject - matrix restricting by injection
2567: . coarse - coarse level DM to update
2568: - ctx - optional user-defined function context
2570: Level: advanced
2572: Notes:
2573: This function is only needed if auxiliary data needs to be set up on coarse grids.
2575: If this function is called multiple times, the hooks will be run in the order they are added.
2577: In order to compose with nonlinear preconditioning without duplicating storage, the hook should be implemented to
2578: extract the finest level information from its context (instead of from the SNES).
2580: This function is currently not available from Fortran.
2582: .seealso: DMCoarsenHookRemove(), DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2583: @*/
2584: PetscErrorCode DMCoarsenHookAdd(DM fine,PetscErrorCode (*coarsenhook)(DM,DM,void*),PetscErrorCode (*restricthook)(DM,Mat,Vec,Mat,DM,void*),void *ctx)
2585: {
2586: PetscErrorCode ierr;
2587: DMCoarsenHookLink link,*p;
2591: for (p=&fine->coarsenhook; *p; p=&(*p)->next) { /* Scan to the end of the current list of hooks */
2592: if ((*p)->coarsenhook == coarsenhook && (*p)->restricthook == restricthook && (*p)->ctx == ctx) return(0);
2593: }
2594: PetscNew(&link);
2595: link->coarsenhook = coarsenhook;
2596: link->restricthook = restricthook;
2597: link->ctx = ctx;
2598: link->next = NULL;
2599: *p = link;
2600: return(0);
2601: }
2603: /*@C
2604: DMCoarsenHookRemove - remove a callback from the list of hooks to be run when restricting a nonlinear problem to the coarse grid
2606: Logically Collective
2608: Input Arguments:
2609: + fine - nonlinear solver context on which to run a hook when restricting to a coarser level
2610: . coarsenhook - function to run when setting up a coarser level
2611: . restricthook - function to run to update data on coarser levels (once per SNESSolve())
2612: - ctx - [optional] user-defined context for provide data for the hooks (may be NULL)
2614: Level: advanced
2616: Notes:
2617: This function does nothing if the hook is not in the list.
2619: This function is currently not available from Fortran.
2621: .seealso: DMCoarsenHookAdd(), DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2622: @*/
2623: PetscErrorCode DMCoarsenHookRemove(DM fine,PetscErrorCode (*coarsenhook)(DM,DM,void*),PetscErrorCode (*restricthook)(DM,Mat,Vec,Mat,DM,void*),void *ctx)
2624: {
2625: PetscErrorCode ierr;
2626: DMCoarsenHookLink link,*p;
2630: for (p=&fine->coarsenhook; *p; p=&(*p)->next) { /* Search the list of current hooks */
2631: if ((*p)->coarsenhook == coarsenhook && (*p)->restricthook == restricthook && (*p)->ctx == ctx) {
2632: link = *p;
2633: *p = link->next;
2634: PetscFree(link);
2635: break;
2636: }
2637: }
2638: return(0);
2639: }
2642: /*@
2643: DMRestrict - restricts user-defined problem data to a coarser DM by running hooks registered by DMCoarsenHookAdd()
2645: Collective if any hooks are
2647: Input Arguments:
2648: + fine - finer DM to use as a base
2649: . restrct - restriction matrix, apply using MatRestrict()
2650: . rscale - scaling vector for restriction
2651: . inject - injection matrix, also use MatRestrict()
2652: - coarse - coarser DM to update
2654: Level: developer
2656: .seealso: DMCoarsenHookAdd(), MatRestrict()
2657: @*/
2658: PetscErrorCode DMRestrict(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse)
2659: {
2660: PetscErrorCode ierr;
2661: DMCoarsenHookLink link;
2664: for (link=fine->coarsenhook; link; link=link->next) {
2665: if (link->restricthook) {
2666: (*link->restricthook)(fine,restrct,rscale,inject,coarse,link->ctx);
2667: }
2668: }
2669: return(0);
2670: }
2672: /*@C
2673: DMSubDomainHookAdd - adds a callback to be run when restricting a problem to the coarse grid
2675: Logically Collective
2677: Input Arguments:
2678: + global - global DM
2679: . ddhook - function to run to pass data to the decomposition DM upon its creation
2680: . restricthook - function to run to update data on block solve (at the beginning of the block solve)
2681: - ctx - [optional] user-defined context for provide data for the hooks (may be NULL)
2684: Calling sequence for ddhook:
2685: $ ddhook(DM global,DM block,void *ctx)
2687: + global - global DM
2688: . block - block DM
2689: - ctx - optional user-defined function context
2691: Calling sequence for restricthook:
2692: $ restricthook(DM global,VecScatter out,VecScatter in,DM block,void *ctx)
2694: + global - global DM
2695: . out - scatter to the outer (with ghost and overlap points) block vector
2696: . in - scatter to block vector values only owned locally
2697: . block - block DM
2698: - ctx - optional user-defined function context
2700: Level: advanced
2702: Notes:
2703: This function is only needed if auxiliary data needs to be set up on subdomain DMs.
2705: If this function is called multiple times, the hooks will be run in the order they are added.
2707: In order to compose with nonlinear preconditioning without duplicating storage, the hook should be implemented to
2708: extract the global information from its context (instead of from the SNES).
2710: This function is currently not available from Fortran.
2712: .seealso: DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2713: @*/
2714: PetscErrorCode DMSubDomainHookAdd(DM global,PetscErrorCode (*ddhook)(DM,DM,void*),PetscErrorCode (*restricthook)(DM,VecScatter,VecScatter,DM,void*),void *ctx)
2715: {
2716: PetscErrorCode ierr;
2717: DMSubDomainHookLink link,*p;
2721: for (p=&global->subdomainhook; *p; p=&(*p)->next) { /* Scan to the end of the current list of hooks */
2722: if ((*p)->ddhook == ddhook && (*p)->restricthook == restricthook && (*p)->ctx == ctx) return(0);
2723: }
2724: PetscNew(&link);
2725: link->restricthook = restricthook;
2726: link->ddhook = ddhook;
2727: link->ctx = ctx;
2728: link->next = NULL;
2729: *p = link;
2730: return(0);
2731: }
2733: /*@C
2734: DMSubDomainHookRemove - remove a callback from the list to be run when restricting a problem to the coarse grid
2736: Logically Collective
2738: Input Arguments:
2739: + global - global DM
2740: . ddhook - function to run to pass data to the decomposition DM upon its creation
2741: . restricthook - function to run to update data on block solve (at the beginning of the block solve)
2742: - ctx - [optional] user-defined context for provide data for the hooks (may be NULL)
2744: Level: advanced
2746: Notes:
2748: This function is currently not available from Fortran.
2750: .seealso: DMSubDomainHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2751: @*/
2752: PetscErrorCode DMSubDomainHookRemove(DM global,PetscErrorCode (*ddhook)(DM,DM,void*),PetscErrorCode (*restricthook)(DM,VecScatter,VecScatter,DM,void*),void *ctx)
2753: {
2754: PetscErrorCode ierr;
2755: DMSubDomainHookLink link,*p;
2759: for (p=&global->subdomainhook; *p; p=&(*p)->next) { /* Search the list of current hooks */
2760: if ((*p)->ddhook == ddhook && (*p)->restricthook == restricthook && (*p)->ctx == ctx) {
2761: link = *p;
2762: *p = link->next;
2763: PetscFree(link);
2764: break;
2765: }
2766: }
2767: return(0);
2768: }
2770: /*@
2771: DMSubDomainRestrict - restricts user-defined problem data to a block DM by running hooks registered by DMSubDomainHookAdd()
2773: Collective if any hooks are
2775: Input Arguments:
2776: + fine - finer DM to use as a base
2777: . oscatter - scatter from domain global vector filling subdomain global vector with overlap
2778: . gscatter - scatter from domain global vector filling subdomain local vector with ghosts
2779: - coarse - coarer DM to update
2781: Level: developer
2783: .seealso: DMCoarsenHookAdd(), MatRestrict()
2784: @*/
2785: PetscErrorCode DMSubDomainRestrict(DM global,VecScatter oscatter,VecScatter gscatter,DM subdm)
2786: {
2787: PetscErrorCode ierr;
2788: DMSubDomainHookLink link;
2791: for (link=global->subdomainhook; link; link=link->next) {
2792: if (link->restricthook) {
2793: (*link->restricthook)(global,oscatter,gscatter,subdm,link->ctx);
2794: }
2795: }
2796: return(0);
2797: }
2799: /*@
2800: DMGetCoarsenLevel - Get's the number of coarsenings that have generated this DM.
2802: Not Collective
2804: Input Parameter:
2805: . dm - the DM object
2807: Output Parameter:
2808: . level - number of coarsenings
2810: Level: developer
2812: .seealso DMCoarsen(), DMGetRefineLevel(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation()
2814: @*/
2815: PetscErrorCode DMGetCoarsenLevel(DM dm,PetscInt *level)
2816: {
2819: *level = dm->leveldown;
2820: return(0);
2821: }
2825: /*@C
2826: DMRefineHierarchy - Refines a DM object, all levels at once
2828: Collective on DM
2830: Input Parameter:
2831: + dm - the DM object
2832: - nlevels - the number of levels of refinement
2834: Output Parameter:
2835: . dmf - the refined DM hierarchy
2837: Level: developer
2839: .seealso DMCoarsenHierarchy(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation()
2841: @*/
2842: PetscErrorCode DMRefineHierarchy(DM dm,PetscInt nlevels,DM dmf[])
2843: {
2848: if (nlevels < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"nlevels cannot be negative");
2849: if (nlevels == 0) return(0);
2850: if (dm->ops->refinehierarchy) {
2851: (*dm->ops->refinehierarchy)(dm,nlevels,dmf);
2852: } else if (dm->ops->refine) {
2853: PetscInt i;
2855: DMRefine(dm,PetscObjectComm((PetscObject)dm),&dmf[0]);
2856: for (i=1; i<nlevels; i++) {
2857: DMRefine(dmf[i-1],PetscObjectComm((PetscObject)dm),&dmf[i]);
2858: }
2859: } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"No RefineHierarchy for this DM yet");
2860: return(0);
2861: }
2863: /*@C
2864: DMCoarsenHierarchy - Coarsens a DM object, all levels at once
2866: Collective on DM
2868: Input Parameter:
2869: + dm - the DM object
2870: - nlevels - the number of levels of coarsening
2872: Output Parameter:
2873: . dmc - the coarsened DM hierarchy
2875: Level: developer
2877: .seealso DMRefineHierarchy(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation()
2879: @*/
2880: PetscErrorCode DMCoarsenHierarchy(DM dm, PetscInt nlevels, DM dmc[])
2881: {
2886: if (nlevels < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"nlevels cannot be negative");
2887: if (nlevels == 0) return(0);
2889: if (dm->ops->coarsenhierarchy) {
2890: (*dm->ops->coarsenhierarchy)(dm, nlevels, dmc);
2891: } else if (dm->ops->coarsen) {
2892: PetscInt i;
2894: DMCoarsen(dm,PetscObjectComm((PetscObject)dm),&dmc[0]);
2895: for (i=1; i<nlevels; i++) {
2896: DMCoarsen(dmc[i-1],PetscObjectComm((PetscObject)dm),&dmc[i]);
2897: }
2898: } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"No CoarsenHierarchy for this DM yet");
2899: return(0);
2900: }
2902: /*@
2903: DMCreateAggregates - Gets the aggregates that map between
2904: grids associated with two DMs.
2906: Collective on DM
2908: Input Parameters:
2909: + dmc - the coarse grid DM
2910: - dmf - the fine grid DM
2912: Output Parameters:
2913: . rest - the restriction matrix (transpose of the projection matrix)
2915: Level: intermediate
2917: .keywords: interpolation, restriction, multigrid
2919: .seealso: DMRefine(), DMCreateInjection(), DMCreateInterpolation()
2920: @*/
2921: PetscErrorCode DMCreateAggregates(DM dmc, DM dmf, Mat *rest)
2922: {
2928: (*dmc->ops->getaggregates)(dmc, dmf, rest);
2929: return(0);
2930: }
2932: /*@C
2933: DMSetApplicationContextDestroy - Sets a user function that will be called to destroy the application context when the DM is destroyed
2935: Not Collective
2937: Input Parameters:
2938: + dm - the DM object
2939: - destroy - the destroy function
2941: Level: intermediate
2943: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext()
2945: @*/
2946: PetscErrorCode DMSetApplicationContextDestroy(DM dm,PetscErrorCode (*destroy)(void**))
2947: {
2950: dm->ctxdestroy = destroy;
2951: return(0);
2952: }
2954: /*@
2955: DMSetApplicationContext - Set a user context into a DM object
2957: Not Collective
2959: Input Parameters:
2960: + dm - the DM object
2961: - ctx - the user context
2963: Level: intermediate
2965: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext()
2967: @*/
2968: PetscErrorCode DMSetApplicationContext(DM dm,void *ctx)
2969: {
2972: dm->ctx = ctx;
2973: return(0);
2974: }
2976: /*@
2977: DMGetApplicationContext - Gets a user context from a DM object
2979: Not Collective
2981: Input Parameter:
2982: . dm - the DM object
2984: Output Parameter:
2985: . ctx - the user context
2987: Level: intermediate
2989: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext()
2991: @*/
2992: PetscErrorCode DMGetApplicationContext(DM dm,void *ctx)
2993: {
2996: *(void**)ctx = dm->ctx;
2997: return(0);
2998: }
3000: /*@C
3001: DMSetVariableBounds - sets a function to compute the lower and upper bound vectors for SNESVI.
3003: Logically Collective on DM
3005: Input Parameter:
3006: + dm - the DM object
3007: - f - the function that computes variable bounds used by SNESVI (use NULL to cancel a previous function that was set)
3009: Level: intermediate
3011: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext(),
3012: DMSetJacobian()
3014: @*/
3015: PetscErrorCode DMSetVariableBounds(DM dm,PetscErrorCode (*f)(DM,Vec,Vec))
3016: {
3018: dm->ops->computevariablebounds = f;
3019: return(0);
3020: }
3022: /*@
3023: DMHasVariableBounds - does the DM object have a variable bounds function?
3025: Not Collective
3027: Input Parameter:
3028: . dm - the DM object to destroy
3030: Output Parameter:
3031: . flg - PETSC_TRUE if the variable bounds function exists
3033: Level: developer
3035: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext()
3037: @*/
3038: PetscErrorCode DMHasVariableBounds(DM dm,PetscBool *flg)
3039: {
3041: *flg = (dm->ops->computevariablebounds) ? PETSC_TRUE : PETSC_FALSE;
3042: return(0);
3043: }
3045: /*@C
3046: DMComputeVariableBounds - compute variable bounds used by SNESVI.
3048: Logically Collective on DM
3050: Input Parameters:
3051: . dm - the DM object
3053: Output parameters:
3054: + xl - lower bound
3055: - xu - upper bound
3057: Level: advanced
3059: Notes:
3060: This is generally not called by users. It calls the function provided by the user with DMSetVariableBounds()
3062: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext()
3064: @*/
3065: PetscErrorCode DMComputeVariableBounds(DM dm, Vec xl, Vec xu)
3066: {
3072: if (dm->ops->computevariablebounds) {
3073: (*dm->ops->computevariablebounds)(dm, xl,xu);
3074: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "This DM is incapable of computing variable bounds.");
3075: return(0);
3076: }
3078: /*@
3079: DMHasColoring - does the DM object have a method of providing a coloring?
3081: Not Collective
3083: Input Parameter:
3084: . dm - the DM object
3086: Output Parameter:
3087: . flg - PETSC_TRUE if the DM has facilities for DMCreateColoring().
3089: Level: developer
3091: .seealso DMHasFunction(), DMCreateColoring()
3093: @*/
3094: PetscErrorCode DMHasColoring(DM dm,PetscBool *flg)
3095: {
3097: *flg = (dm->ops->getcoloring) ? PETSC_TRUE : PETSC_FALSE;
3098: return(0);
3099: }
3101: /*@
3102: DMHasCreateRestriction - does the DM object have a method of providing a restriction?
3104: Not Collective
3106: Input Parameter:
3107: . dm - the DM object
3109: Output Parameter:
3110: . flg - PETSC_TRUE if the DM has facilities for DMCreateRestriction().
3112: Level: developer
3114: .seealso DMHasFunction(), DMCreateRestriction()
3116: @*/
3117: PetscErrorCode DMHasCreateRestriction(DM dm,PetscBool *flg)
3118: {
3120: *flg = (dm->ops->createrestriction) ? PETSC_TRUE : PETSC_FALSE;
3121: return(0);
3122: }
3125: /*@
3126: DMHasCreateInjection - does the DM object have a method of providing an injection?
3128: Not Collective
3130: Input Parameter:
3131: . dm - the DM object
3133: Output Parameter:
3134: . flg - PETSC_TRUE if the DM has facilities for DMCreateInjection().
3136: Level: developer
3138: .seealso DMHasFunction(), DMCreateInjection()
3140: @*/
3141: PetscErrorCode DMHasCreateInjection(DM dm,PetscBool *flg)
3142: {
3145: if (!dm->ops->hascreateinjection) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DMHasCreateInjection not implemented for this type");
3146: (*dm->ops->hascreateinjection)(dm,flg);
3147: return(0);
3148: }
3150: /*@C
3151: DMSetVec - set the vector at which to compute residual, Jacobian and VI bounds, if the problem is nonlinear.
3153: Collective on DM
3155: Input Parameter:
3156: + dm - the DM object
3157: - x - location to compute residual and Jacobian, if NULL is passed to those routines; will be NULL for linear problems.
3159: Level: developer
3161: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext()
3163: @*/
3164: PetscErrorCode DMSetVec(DM dm,Vec x)
3165: {
3169: if (x) {
3170: if (!dm->x) {
3171: DMCreateGlobalVector(dm,&dm->x);
3172: }
3173: VecCopy(x,dm->x);
3174: } else if (dm->x) {
3175: VecDestroy(&dm->x);
3176: }
3177: return(0);
3178: }
3180: PetscFunctionList DMList = NULL;
3181: PetscBool DMRegisterAllCalled = PETSC_FALSE;
3183: /*@C
3184: DMSetType - Builds a DM, for a particular DM implementation.
3186: Collective on DM
3188: Input Parameters:
3189: + dm - The DM object
3190: - method - The name of the DM type
3192: Options Database Key:
3193: . -dm_type <type> - Sets the DM type; use -help for a list of available types
3195: Notes:
3196: See "petsc/include/petscdm.h" for available DM types (for instance, DM1D, DM2D, or DM3D).
3198: Level: intermediate
3200: .keywords: DM, set, type
3201: .seealso: DMGetType(), DMCreate()
3202: @*/
3203: PetscErrorCode DMSetType(DM dm, DMType method)
3204: {
3205: PetscErrorCode (*r)(DM);
3206: PetscBool match;
3211: PetscObjectTypeCompare((PetscObject) dm, method, &match);
3212: if (match) return(0);
3214: DMRegisterAll();
3215: PetscFunctionListFind(DMList,method,&r);
3216: if (!r) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown DM type: %s", method);
3218: if (dm->ops->destroy) {
3219: (*dm->ops->destroy)(dm);
3220: dm->ops->destroy = NULL;
3221: }
3222: (*r)(dm);
3223: PetscObjectChangeTypeName((PetscObject)dm,method);
3224: return(0);
3225: }
3227: /*@C
3228: DMGetType - Gets the DM type name (as a string) from the DM.
3230: Not Collective
3232: Input Parameter:
3233: . dm - The DM
3235: Output Parameter:
3236: . type - The DM type name
3238: Level: intermediate
3240: .keywords: DM, get, type, name
3241: .seealso: DMSetType(), DMCreate()
3242: @*/
3243: PetscErrorCode DMGetType(DM dm, DMType *type)
3244: {
3250: DMRegisterAll();
3251: *type = ((PetscObject)dm)->type_name;
3252: return(0);
3253: }
3255: /*@C
3256: DMConvert - Converts a DM to another DM, either of the same or different type.
3258: Collective on DM
3260: Input Parameters:
3261: + dm - the DM
3262: - newtype - new DM type (use "same" for the same type)
3264: Output Parameter:
3265: . M - pointer to new DM
3267: Notes:
3268: Cannot be used to convert a sequential DM to parallel or parallel to sequential,
3269: the MPI communicator of the generated DM is always the same as the communicator
3270: of the input DM.
3272: Level: intermediate
3274: .seealso: DMCreate()
3275: @*/
3276: PetscErrorCode DMConvert(DM dm, DMType newtype, DM *M)
3277: {
3278: DM B;
3279: char convname[256];
3280: PetscBool sametype/*, issame */;
3287: PetscObjectTypeCompare((PetscObject) dm, newtype, &sametype);
3288: /* PetscStrcmp(newtype, "same", &issame); */
3289: if (sametype) {
3290: *M = dm;
3291: PetscObjectReference((PetscObject) dm);
3292: return(0);
3293: } else {
3294: PetscErrorCode (*conv)(DM, DMType, DM*) = NULL;
3296: /*
3297: Order of precedence:
3298: 1) See if a specialized converter is known to the current DM.
3299: 2) See if a specialized converter is known to the desired DM class.
3300: 3) See if a good general converter is registered for the desired class
3301: 4) See if a good general converter is known for the current matrix.
3302: 5) Use a really basic converter.
3303: */
3305: /* 1) See if a specialized converter is known to the current DM and the desired class */
3306: PetscStrncpy(convname,"DMConvert_",sizeof(convname));
3307: PetscStrlcat(convname,((PetscObject) dm)->type_name,sizeof(convname));
3308: PetscStrlcat(convname,"_",sizeof(convname));
3309: PetscStrlcat(convname,newtype,sizeof(convname));
3310: PetscStrlcat(convname,"_C",sizeof(convname));
3311: PetscObjectQueryFunction((PetscObject)dm,convname,&conv);
3312: if (conv) goto foundconv;
3314: /* 2) See if a specialized converter is known to the desired DM class. */
3315: DMCreate(PetscObjectComm((PetscObject)dm), &B);
3316: DMSetType(B, newtype);
3317: PetscStrncpy(convname,"DMConvert_",sizeof(convname));
3318: PetscStrlcat(convname,((PetscObject) dm)->type_name,sizeof(convname));
3319: PetscStrlcat(convname,"_",sizeof(convname));
3320: PetscStrlcat(convname,newtype,sizeof(convname));
3321: PetscStrlcat(convname,"_C",sizeof(convname));
3322: PetscObjectQueryFunction((PetscObject)B,convname,&conv);
3323: if (conv) {
3324: DMDestroy(&B);
3325: goto foundconv;
3326: }
3328: #if 0
3329: /* 3) See if a good general converter is registered for the desired class */
3330: conv = B->ops->convertfrom;
3331: DMDestroy(&B);
3332: if (conv) goto foundconv;
3334: /* 4) See if a good general converter is known for the current matrix */
3335: if (dm->ops->convert) {
3336: conv = dm->ops->convert;
3337: }
3338: if (conv) goto foundconv;
3339: #endif
3341: /* 5) Use a really basic converter. */
3342: SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "No conversion possible between DM types %s and %s", ((PetscObject) dm)->type_name, newtype);
3344: foundconv:
3345: PetscLogEventBegin(DM_Convert,dm,0,0,0);
3346: (*conv)(dm,newtype,M);
3347: /* Things that are independent of DM type: We should consult DMClone() here */
3348: {
3349: PetscBool isper;
3350: const PetscReal *maxCell, *L;
3351: const DMBoundaryType *bd;
3352: DMGetPeriodicity(dm, &isper, &maxCell, &L, &bd);
3353: DMSetPeriodicity(*M, isper, maxCell, L, bd);
3354: }
3355: PetscLogEventEnd(DM_Convert,dm,0,0,0);
3356: }
3357: PetscObjectStateIncrease((PetscObject) *M);
3358: return(0);
3359: }
3361: /*--------------------------------------------------------------------------------------------------------------------*/
3363: /*@C
3364: DMRegister - Adds a new DM component implementation
3366: Not Collective
3368: Input Parameters:
3369: + name - The name of a new user-defined creation routine
3370: - create_func - The creation routine itself
3372: Notes:
3373: DMRegister() may be called multiple times to add several user-defined DMs
3376: Sample usage:
3377: .vb
3378: DMRegister("my_da", MyDMCreate);
3379: .ve
3381: Then, your DM type can be chosen with the procedural interface via
3382: .vb
3383: DMCreate(MPI_Comm, DM *);
3384: DMSetType(DM,"my_da");
3385: .ve
3386: or at runtime via the option
3387: .vb
3388: -da_type my_da
3389: .ve
3391: Level: advanced
3393: .keywords: DM, register
3394: .seealso: DMRegisterAll(), DMRegisterDestroy()
3396: @*/
3397: PetscErrorCode DMRegister(const char sname[],PetscErrorCode (*function)(DM))
3398: {
3402: DMInitializePackage();
3403: PetscFunctionListAdd(&DMList,sname,function);
3404: return(0);
3405: }
3407: /*@C
3408: DMLoad - Loads a DM that has been stored in binary with DMView().
3410: Collective on PetscViewer
3412: Input Parameters:
3413: + newdm - the newly loaded DM, this needs to have been created with DMCreate() or
3414: some related function before a call to DMLoad().
3415: - viewer - binary file viewer, obtained from PetscViewerBinaryOpen() or
3416: HDF5 file viewer, obtained from PetscViewerHDF5Open()
3418: Level: intermediate
3420: Notes:
3421: The type is determined by the data in the file, any type set into the DM before this call is ignored.
3423: Notes for advanced users:
3424: Most users should not need to know the details of the binary storage
3425: format, since DMLoad() and DMView() completely hide these details.
3426: But for anyone who's interested, the standard binary matrix storage
3427: format is
3428: .vb
3429: has not yet been determined
3430: .ve
3432: .seealso: PetscViewerBinaryOpen(), DMView(), MatLoad(), VecLoad()
3433: @*/
3434: PetscErrorCode DMLoad(DM newdm, PetscViewer viewer)
3435: {
3436: PetscBool isbinary, ishdf5;
3442: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
3443: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);
3444: if (isbinary) {
3445: PetscInt classid;
3446: char type[256];
3448: PetscViewerBinaryRead(viewer,&classid,1,NULL,PETSC_INT);
3449: if (classid != DM_FILE_CLASSID) SETERRQ1(PetscObjectComm((PetscObject)newdm),PETSC_ERR_ARG_WRONG,"Not DM next in file, classid found %d",(int)classid);
3450: PetscViewerBinaryRead(viewer,type,256,NULL,PETSC_CHAR);
3451: DMSetType(newdm, type);
3452: if (newdm->ops->load) {(*newdm->ops->load)(newdm,viewer);}
3453: } else if (ishdf5) {
3454: if (newdm->ops->load) {(*newdm->ops->load)(newdm,viewer);}
3455: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid viewer; open viewer with PetscViewerBinaryOpen() or PetscViewerHDF5Open()");
3456: return(0);
3457: }
3459: /******************************** FEM Support **********************************/
3461: PetscErrorCode DMPrintCellVector(PetscInt c, const char name[], PetscInt len, const PetscScalar x[])
3462: {
3463: PetscInt f;
3467: PetscPrintf(PETSC_COMM_SELF, "Cell %D Element %s\n", c, name);
3468: for (f = 0; f < len; ++f) {
3469: PetscPrintf(PETSC_COMM_SELF, " | %g |\n", (double)PetscRealPart(x[f]));
3470: }
3471: return(0);
3472: }
3474: PetscErrorCode DMPrintCellMatrix(PetscInt c, const char name[], PetscInt rows, PetscInt cols, const PetscScalar A[])
3475: {
3476: PetscInt f, g;
3480: PetscPrintf(PETSC_COMM_SELF, "Cell %D Element %s\n", c, name);
3481: for (f = 0; f < rows; ++f) {
3482: PetscPrintf(PETSC_COMM_SELF, " |");
3483: for (g = 0; g < cols; ++g) {
3484: PetscPrintf(PETSC_COMM_SELF, " % 9.5g", PetscRealPart(A[f*cols+g]));
3485: }
3486: PetscPrintf(PETSC_COMM_SELF, " |\n");
3487: }
3488: return(0);
3489: }
3491: PetscErrorCode DMPrintLocalVec(DM dm, const char name[], PetscReal tol, Vec X)
3492: {
3493: PetscInt localSize, bs;
3494: PetscMPIInt size;
3495: Vec x, xglob;
3496: const PetscScalar *xarray;
3497: PetscErrorCode ierr;
3500: MPI_Comm_size(PetscObjectComm((PetscObject) dm),&size);
3501: VecDuplicate(X, &x);
3502: VecCopy(X, x);
3503: VecChop(x, tol);
3504: PetscPrintf(PetscObjectComm((PetscObject) dm),"%s:\n",name);
3505: if (size > 1) {
3506: VecGetLocalSize(x,&localSize);
3507: VecGetArrayRead(x,&xarray);
3508: VecGetBlockSize(x,&bs);
3509: VecCreateMPIWithArray(PetscObjectComm((PetscObject) dm),bs,localSize,PETSC_DETERMINE,xarray,&xglob);
3510: } else {
3511: xglob = x;
3512: }
3513: VecView(xglob,PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject) dm)));
3514: if (size > 1) {
3515: VecDestroy(&xglob);
3516: VecRestoreArrayRead(x,&xarray);
3517: }
3518: VecDestroy(&x);
3519: return(0);
3520: }
3522: /*@
3523: DMGetSection - Get the PetscSection encoding the local data layout for the DM.
3525: Input Parameter:
3526: . dm - The DM
3528: Output Parameter:
3529: . section - The PetscSection
3531: Options Database Keys:
3532: . -dm_petscsection_view - View the Section created by the DM
3534: Level: intermediate
3536: Note: This gets a borrowed reference, so the user should not destroy this PetscSection.
3538: .seealso: DMSetSection(), DMGetGlobalSection()
3539: @*/
3540: PetscErrorCode DMGetSection(DM dm, PetscSection *section)
3541: {
3547: if (!dm->defaultSection && dm->ops->createdefaultsection) {
3548: (*dm->ops->createdefaultsection)(dm);
3549: if (dm->defaultSection) {PetscObjectViewFromOptions((PetscObject) dm->defaultSection, NULL, "-dm_petscsection_view");}
3550: }
3551: *section = dm->defaultSection;
3552: return(0);
3553: }
3555: /*@
3556: DMSetSection - Set the PetscSection encoding the local data layout for the DM.
3558: Input Parameters:
3559: + dm - The DM
3560: - section - The PetscSection
3562: Level: intermediate
3564: Note: Any existing Section will be destroyed
3566: .seealso: DMSetSection(), DMGetGlobalSection()
3567: @*/
3568: PetscErrorCode DMSetSection(DM dm, PetscSection section)
3569: {
3570: PetscInt numFields = 0;
3571: PetscInt f;
3576: if (section) {
3578: PetscObjectReference((PetscObject)section);
3579: }
3580: PetscSectionDestroy(&dm->defaultSection);
3581: dm->defaultSection = section;
3582: if (section) {PetscSectionGetNumFields(dm->defaultSection, &numFields);}
3583: if (numFields) {
3584: DMSetNumFields(dm, numFields);
3585: for (f = 0; f < numFields; ++f) {
3586: PetscObject disc;
3587: const char *name;
3589: PetscSectionGetFieldName(dm->defaultSection, f, &name);
3590: DMGetField(dm, f, &disc);
3591: PetscObjectSetName(disc, name);
3592: }
3593: }
3594: /* The global section will be rebuilt in the next call to DMGetGlobalSection(). */
3595: PetscSectionDestroy(&dm->defaultGlobalSection);
3596: return(0);
3597: }
3599: /*@
3600: DMGetDefaultConstraints - Get the PetscSection and Mat the specify the local constraint interpolation. See DMSetDefaultConstraints() for a description of the purpose of constraint interpolation.
3602: not collective
3604: Input Parameter:
3605: . dm - The DM
3607: Output Parameter:
3608: + section - The PetscSection describing the range of the constraint matrix: relates rows of the constraint matrix to dofs of the default section. Returns NULL if there are no local constraints.
3609: - mat - The Mat that interpolates local constraints: its width should be the layout size of the default section. Returns NULL if there are no local constraints.
3611: Level: advanced
3613: Note: This gets borrowed references, so the user should not destroy the PetscSection or the Mat.
3615: .seealso: DMSetDefaultConstraints()
3616: @*/
3617: PetscErrorCode DMGetDefaultConstraints(DM dm, PetscSection *section, Mat *mat)
3618: {
3623: if (!dm->defaultConstraintSection && !dm->defaultConstraintMat && dm->ops->createdefaultconstraints) {(*dm->ops->createdefaultconstraints)(dm);}
3624: if (section) {*section = dm->defaultConstraintSection;}
3625: if (mat) {*mat = dm->defaultConstraintMat;}
3626: return(0);
3627: }
3629: /*@
3630: DMSetDefaultConstraints - Set the PetscSection and Mat the specify the local constraint interpolation.
3632: If a constraint matrix is specified, then it is applied during DMGlobalToLocalEnd() when mode is INSERT_VALUES, INSERT_BC_VALUES, or INSERT_ALL_VALUES. Without a constraint matrix, the local vector l returned by DMGlobalToLocalEnd() contains values that have been scattered from a global vector without modification; with a constraint matrix A, l is modified by computing c = A * l, l[s[i]] = c[i], where the scatter s is defined by the PetscSection returned by DMGetDefaultConstraintMatrix().
3634: If a constraint matrix is specified, then its adjoint is applied during DMLocalToGlobalBegin() when mode is ADD_VALUES, ADD_BC_VALUES, or ADD_ALL_VALUES. Without a constraint matrix, the local vector l is accumulated into a global vector without modification; with a constraint matrix A, l is first modified by computing c[i] = l[s[i]], l[s[i]] = 0, l = l + A'*c, which is the adjoint of the operation described above.
3636: collective on dm
3638: Input Parameters:
3639: + dm - The DM
3640: + section - The PetscSection describing the range of the constraint matrix: relates rows of the constraint matrix to dofs of the default section. Must have a local communicator (PETSC_COMM_SELF or derivative).
3641: - mat - The Mat that interpolates local constraints: its width should be the layout size of the default section: NULL indicates no constraints. Must have a local communicator (PETSC_COMM_SELF or derivative).
3643: Level: advanced
3645: Note: This increments the references of the PetscSection and the Mat, so they user can destroy them
3647: .seealso: DMGetDefaultConstraints()
3648: @*/
3649: PetscErrorCode DMSetDefaultConstraints(DM dm, PetscSection section, Mat mat)
3650: {
3651: PetscMPIInt result;
3656: if (section) {
3658: MPI_Comm_compare(PETSC_COMM_SELF,PetscObjectComm((PetscObject)section),&result);
3659: if (result != MPI_CONGRUENT && result != MPI_IDENT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMECOMM,"constraint section must have local communicator");
3660: }
3661: if (mat) {
3663: MPI_Comm_compare(PETSC_COMM_SELF,PetscObjectComm((PetscObject)mat),&result);
3664: if (result != MPI_CONGRUENT && result != MPI_IDENT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMECOMM,"constraint matrix must have local communicator");
3665: }
3666: PetscObjectReference((PetscObject)section);
3667: PetscSectionDestroy(&dm->defaultConstraintSection);
3668: dm->defaultConstraintSection = section;
3669: PetscObjectReference((PetscObject)mat);
3670: MatDestroy(&dm->defaultConstraintMat);
3671: dm->defaultConstraintMat = mat;
3672: return(0);
3673: }
3675: #if defined(PETSC_USE_DEBUG)
3676: /*
3677: DMDefaultSectionCheckConsistency - Check the consistentcy of the global and local sections.
3679: Input Parameters:
3680: + dm - The DM
3681: . localSection - PetscSection describing the local data layout
3682: - globalSection - PetscSection describing the global data layout
3684: Level: intermediate
3686: .seealso: DMGetDefaultSF(), DMSetDefaultSF()
3687: */
3688: static PetscErrorCode DMDefaultSectionCheckConsistency_Internal(DM dm, PetscSection localSection, PetscSection globalSection)
3689: {
3690: MPI_Comm comm;
3691: PetscLayout layout;
3692: const PetscInt *ranges;
3693: PetscInt pStart, pEnd, p, nroots;
3694: PetscMPIInt size, rank;
3695: PetscBool valid = PETSC_TRUE, gvalid;
3696: PetscErrorCode ierr;
3699: PetscObjectGetComm((PetscObject)dm,&comm);
3701: MPI_Comm_size(comm, &size);
3702: MPI_Comm_rank(comm, &rank);
3703: PetscSectionGetChart(globalSection, &pStart, &pEnd);
3704: PetscSectionGetConstrainedStorageSize(globalSection, &nroots);
3705: PetscLayoutCreate(comm, &layout);
3706: PetscLayoutSetBlockSize(layout, 1);
3707: PetscLayoutSetLocalSize(layout, nroots);
3708: PetscLayoutSetUp(layout);
3709: PetscLayoutGetRanges(layout, &ranges);
3710: for (p = pStart; p < pEnd; ++p) {
3711: PetscInt dof, cdof, off, gdof, gcdof, goff, gsize, d;
3713: PetscSectionGetDof(localSection, p, &dof);
3714: PetscSectionGetOffset(localSection, p, &off);
3715: PetscSectionGetConstraintDof(localSection, p, &cdof);
3716: PetscSectionGetDof(globalSection, p, &gdof);
3717: PetscSectionGetConstraintDof(globalSection, p, &gcdof);
3718: PetscSectionGetOffset(globalSection, p, &goff);
3719: if (!gdof) continue; /* Censored point */
3720: if ((gdof < 0 ? -(gdof+1) : gdof) != dof) {PetscSynchronizedPrintf(comm, "[%d]Global dof %d for point %d not equal to local dof %d\n", rank, gdof, p, dof); valid = PETSC_FALSE;}
3721: if (gcdof && (gcdof != cdof)) {PetscSynchronizedPrintf(comm, "[%d]Global constraints %d for point %d not equal to local constraints %d\n", rank, gcdof, p, cdof); valid = PETSC_FALSE;}
3722: if (gdof < 0) {
3723: gsize = gdof < 0 ? -(gdof+1)-gcdof : gdof-gcdof;
3724: for (d = 0; d < gsize; ++d) {
3725: PetscInt offset = -(goff+1) + d, r;
3727: PetscFindInt(offset,size+1,ranges,&r);
3728: if (r < 0) r = -(r+2);
3729: if ((r < 0) || (r >= size)) {PetscSynchronizedPrintf(comm, "[%d]Point %d mapped to invalid process %d (%d, %d)\n", rank, p, r, gdof, goff); valid = PETSC_FALSE;break;}
3730: }
3731: }
3732: }
3733: PetscLayoutDestroy(&layout);
3734: PetscSynchronizedFlush(comm, NULL);
3735: MPIU_Allreduce(&valid, &gvalid, 1, MPIU_BOOL, MPI_LAND, comm);
3736: if (!gvalid) {
3737: DMView(dm, NULL);
3738: SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Inconsistent local and global sections");
3739: }
3740: return(0);
3741: }
3742: #endif
3744: /*@
3745: DMGetGlobalSection - Get the PetscSection encoding the global data layout for the DM.
3747: Collective on DM
3749: Input Parameter:
3750: . dm - The DM
3752: Output Parameter:
3753: . section - The PetscSection
3755: Level: intermediate
3757: Note: This gets a borrowed reference, so the user should not destroy this PetscSection.
3759: .seealso: DMSetSection(), DMGetSection()
3760: @*/
3761: PetscErrorCode DMGetGlobalSection(DM dm, PetscSection *section)
3762: {
3768: if (!dm->defaultGlobalSection) {
3769: PetscSection s;
3771: DMGetSection(dm, &s);
3772: if (!s) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONGSTATE, "DM must have a default PetscSection in order to create a global PetscSection");
3773: if (!dm->sf) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DM must have a default PetscSF in order to create a global PetscSection");
3774: PetscSectionCreateGlobalSection(s, dm->sf, PETSC_FALSE, PETSC_FALSE, &dm->defaultGlobalSection);
3775: PetscLayoutDestroy(&dm->map);
3776: PetscSectionGetValueLayout(PetscObjectComm((PetscObject)dm), dm->defaultGlobalSection, &dm->map);
3777: PetscSectionViewFromOptions(dm->defaultGlobalSection, NULL, "-global_section_view");
3778: }
3779: *section = dm->defaultGlobalSection;
3780: return(0);
3781: }
3783: /*@
3784: DMSetGlobalSection - Set the PetscSection encoding the global data layout for the DM.
3786: Input Parameters:
3787: + dm - The DM
3788: - section - The PetscSection, or NULL
3790: Level: intermediate
3792: Note: Any existing Section will be destroyed
3794: .seealso: DMGetGlobalSection(), DMSetSection()
3795: @*/
3796: PetscErrorCode DMSetGlobalSection(DM dm, PetscSection section)
3797: {
3803: PetscObjectReference((PetscObject)section);
3804: PetscSectionDestroy(&dm->defaultGlobalSection);
3805: dm->defaultGlobalSection = section;
3806: #if defined(PETSC_USE_DEBUG)
3807: if (section) {DMDefaultSectionCheckConsistency_Internal(dm, dm->defaultSection, section);}
3808: #endif
3809: return(0);
3810: }
3812: /*@
3813: DMGetDefaultSF - Get the PetscSF encoding the parallel dof overlap for the DM. If it has not been set,
3814: it is created from the default PetscSection layouts in the DM.
3816: Input Parameter:
3817: . dm - The DM
3819: Output Parameter:
3820: . sf - The PetscSF
3822: Level: intermediate
3824: Note: This gets a borrowed reference, so the user should not destroy this PetscSF.
3826: .seealso: DMSetDefaultSF(), DMCreateDefaultSF()
3827: @*/
3828: PetscErrorCode DMGetDefaultSF(DM dm, PetscSF *sf)
3829: {
3830: PetscInt nroots;
3836: PetscSFGetGraph(dm->defaultSF, &nroots, NULL, NULL, NULL);
3837: if (nroots < 0) {
3838: PetscSection section, gSection;
3840: DMGetSection(dm, §ion);
3841: if (section) {
3842: DMGetGlobalSection(dm, &gSection);
3843: DMCreateDefaultSF(dm, section, gSection);
3844: } else {
3845: *sf = NULL;
3846: return(0);
3847: }
3848: }
3849: *sf = dm->defaultSF;
3850: return(0);
3851: }
3853: /*@
3854: DMSetDefaultSF - Set the PetscSF encoding the parallel dof overlap for the DM
3856: Input Parameters:
3857: + dm - The DM
3858: - sf - The PetscSF
3860: Level: intermediate
3862: Note: Any previous SF is destroyed
3864: .seealso: DMGetDefaultSF(), DMCreateDefaultSF()
3865: @*/
3866: PetscErrorCode DMSetDefaultSF(DM dm, PetscSF sf)
3867: {
3873: PetscSFDestroy(&dm->defaultSF);
3874: dm->defaultSF = sf;
3875: return(0);
3876: }
3878: /*@C
3879: DMCreateDefaultSF - Create the PetscSF encoding the parallel dof overlap for the DM based upon the PetscSections
3880: describing the data layout.
3882: Input Parameters:
3883: + dm - The DM
3884: . localSection - PetscSection describing the local data layout
3885: - globalSection - PetscSection describing the global data layout
3887: Level: intermediate
3889: .seealso: DMGetDefaultSF(), DMSetDefaultSF()
3890: @*/
3891: PetscErrorCode DMCreateDefaultSF(DM dm, PetscSection localSection, PetscSection globalSection)
3892: {
3893: MPI_Comm comm;
3894: PetscLayout layout;
3895: const PetscInt *ranges;
3896: PetscInt *local;
3897: PetscSFNode *remote;
3898: PetscInt pStart, pEnd, p, nroots, nleaves = 0, l;
3899: PetscMPIInt size, rank;
3903: PetscObjectGetComm((PetscObject)dm,&comm);
3905: MPI_Comm_size(comm, &size);
3906: MPI_Comm_rank(comm, &rank);
3907: PetscSectionGetChart(globalSection, &pStart, &pEnd);
3908: PetscSectionGetConstrainedStorageSize(globalSection, &nroots);
3909: PetscLayoutCreate(comm, &layout);
3910: PetscLayoutSetBlockSize(layout, 1);
3911: PetscLayoutSetLocalSize(layout, nroots);
3912: PetscLayoutSetUp(layout);
3913: PetscLayoutGetRanges(layout, &ranges);
3914: for (p = pStart; p < pEnd; ++p) {
3915: PetscInt gdof, gcdof;
3917: PetscSectionGetDof(globalSection, p, &gdof);
3918: PetscSectionGetConstraintDof(globalSection, p, &gcdof);
3919: if (gcdof > (gdof < 0 ? -(gdof+1) : gdof)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %d has %d constraints > %d dof", p, gcdof, (gdof < 0 ? -(gdof+1) : gdof));
3920: nleaves += gdof < 0 ? -(gdof+1)-gcdof : gdof-gcdof;
3921: }
3922: PetscMalloc1(nleaves, &local);
3923: PetscMalloc1(nleaves, &remote);
3924: for (p = pStart, l = 0; p < pEnd; ++p) {
3925: const PetscInt *cind;
3926: PetscInt dof, cdof, off, gdof, gcdof, goff, gsize, d, c;
3928: PetscSectionGetDof(localSection, p, &dof);
3929: PetscSectionGetOffset(localSection, p, &off);
3930: PetscSectionGetConstraintDof(localSection, p, &cdof);
3931: PetscSectionGetConstraintIndices(localSection, p, &cind);
3932: PetscSectionGetDof(globalSection, p, &gdof);
3933: PetscSectionGetConstraintDof(globalSection, p, &gcdof);
3934: PetscSectionGetOffset(globalSection, p, &goff);
3935: if (!gdof) continue; /* Censored point */
3936: gsize = gdof < 0 ? -(gdof+1)-gcdof : gdof-gcdof;
3937: if (gsize != dof-cdof) {
3938: 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);
3939: cdof = 0; /* Ignore constraints */
3940: }
3941: for (d = 0, c = 0; d < dof; ++d) {
3942: if ((c < cdof) && (cind[c] == d)) {++c; continue;}
3943: local[l+d-c] = off+d;
3944: }
3945: if (gdof < 0) {
3946: for (d = 0; d < gsize; ++d, ++l) {
3947: PetscInt offset = -(goff+1) + d, r;
3949: PetscFindInt(offset,size+1,ranges,&r);
3950: if (r < 0) r = -(r+2);
3951: if ((r < 0) || (r >= size)) SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %d mapped to invalid process %d (%d, %d)", p, r, gdof, goff);
3952: remote[l].rank = r;
3953: remote[l].index = offset - ranges[r];
3954: }
3955: } else {
3956: for (d = 0; d < gsize; ++d, ++l) {
3957: remote[l].rank = rank;
3958: remote[l].index = goff+d - ranges[rank];
3959: }
3960: }
3961: }
3962: if (l != nleaves) SETERRQ2(comm, PETSC_ERR_PLIB, "Iteration error, l %d != nleaves %d", l, nleaves);
3963: PetscLayoutDestroy(&layout);
3964: PetscSFSetGraph(dm->defaultSF, nroots, nleaves, local, PETSC_OWN_POINTER, remote, PETSC_OWN_POINTER);
3965: return(0);
3966: }
3968: /*@
3969: DMGetPointSF - Get the PetscSF encoding the parallel section point overlap for the DM.
3971: Input Parameter:
3972: . dm - The DM
3974: Output Parameter:
3975: . sf - The PetscSF
3977: Level: intermediate
3979: Note: This gets a borrowed reference, so the user should not destroy this PetscSF.
3981: .seealso: DMSetPointSF(), DMGetDefaultSF(), DMSetDefaultSF(), DMCreateDefaultSF()
3982: @*/
3983: PetscErrorCode DMGetPointSF(DM dm, PetscSF *sf)
3984: {
3988: *sf = dm->sf;
3989: return(0);
3990: }
3992: /*@
3993: DMSetPointSF - Set the PetscSF encoding the parallel section point overlap for the DM.
3995: Input Parameters:
3996: + dm - The DM
3997: - sf - The PetscSF
3999: Level: intermediate
4001: .seealso: DMGetPointSF(), DMGetDefaultSF(), DMSetDefaultSF(), DMCreateDefaultSF()
4002: @*/
4003: PetscErrorCode DMSetPointSF(DM dm, PetscSF sf)
4004: {
4010: PetscSFDestroy(&dm->sf);
4011: PetscObjectReference((PetscObject) sf);
4012: dm->sf = sf;
4013: return(0);
4014: }
4016: /*@
4017: DMGetDS - Get the PetscDS
4019: Input Parameter:
4020: . dm - The DM
4022: Output Parameter:
4023: . prob - The PetscDS
4025: Level: developer
4027: .seealso: DMSetDS()
4028: @*/
4029: PetscErrorCode DMGetDS(DM dm, PetscDS *prob)
4030: {
4034: *prob = dm->prob;
4035: return(0);
4036: }
4038: /*@
4039: DMSetDS - Set the PetscDS
4041: Input Parameters:
4042: + dm - The DM
4043: - prob - The PetscDS
4045: Level: developer
4047: .seealso: DMGetDS()
4048: @*/
4049: PetscErrorCode DMSetDS(DM dm, PetscDS prob)
4050: {
4051: PetscInt dimEmbed;
4057: PetscObjectReference((PetscObject) prob);
4058: PetscDSDestroy(&dm->prob);
4059: dm->prob = prob;
4060: DMGetCoordinateDim(dm, &dimEmbed);
4061: PetscDSSetCoordinateDimension(prob, dimEmbed);
4062: return(0);
4063: }
4065: PetscErrorCode DMGetNumFields(DM dm, PetscInt *numFields)
4066: {
4071: PetscDSGetNumFields(dm->prob, numFields);
4072: return(0);
4073: }
4075: PetscErrorCode DMSetNumFields(DM dm, PetscInt numFields)
4076: {
4077: PetscInt Nf, f;
4082: PetscDSGetNumFields(dm->prob, &Nf);
4083: for (f = Nf; f < numFields; ++f) {
4084: PetscContainer obj;
4086: PetscContainerCreate(PetscObjectComm((PetscObject) dm), &obj);
4087: PetscDSSetDiscretization(dm->prob, f, (PetscObject) obj);
4088: PetscContainerDestroy(&obj);
4089: }
4090: return(0);
4091: }
4093: /*@
4094: DMGetField - Return the discretization object for a given DM field
4096: Not collective
4098: Input Parameters:
4099: + dm - The DM
4100: - f - The field number
4102: Output Parameter:
4103: . field - The discretization object
4105: Level: developer
4107: .seealso: DMSetField()
4108: @*/
4109: PetscErrorCode DMGetField(DM dm, PetscInt f, PetscObject *field)
4110: {
4115: PetscDSGetDiscretization(dm->prob, f, field);
4116: return(0);
4117: }
4119: /*@
4120: DMSetField - Set the discretization object for a given DM field
4122: Logically collective on DM
4124: Input Parameters:
4125: + dm - The DM
4126: . f - The field number
4127: - field - The discretization object
4129: Level: developer
4131: .seealso: DMGetField()
4132: @*/
4133: PetscErrorCode DMSetField(DM dm, PetscInt f, PetscObject field)
4134: {
4139: PetscDSSetDiscretization(dm->prob, f, field);
4140: return(0);
4141: }
4143: PetscErrorCode DMRestrictHook_Coordinates(DM dm,DM dmc,void *ctx)
4144: {
4145: DM dm_coord,dmc_coord;
4147: Vec coords,ccoords;
4148: Mat inject;
4150: DMGetCoordinateDM(dm,&dm_coord);
4151: DMGetCoordinateDM(dmc,&dmc_coord);
4152: DMGetCoordinates(dm,&coords);
4153: DMGetCoordinates(dmc,&ccoords);
4154: if (coords && !ccoords) {
4155: DMCreateGlobalVector(dmc_coord,&ccoords);
4156: PetscObjectSetName((PetscObject)ccoords,"coordinates");
4157: DMCreateInjection(dmc_coord,dm_coord,&inject);
4158: MatRestrict(inject,coords,ccoords);
4159: MatDestroy(&inject);
4160: DMSetCoordinates(dmc,ccoords);
4161: VecDestroy(&ccoords);
4162: }
4163: return(0);
4164: }
4166: static PetscErrorCode DMSubDomainHook_Coordinates(DM dm,DM subdm,void *ctx)
4167: {
4168: DM dm_coord,subdm_coord;
4170: Vec coords,ccoords,clcoords;
4171: VecScatter *scat_i,*scat_g;
4173: DMGetCoordinateDM(dm,&dm_coord);
4174: DMGetCoordinateDM(subdm,&subdm_coord);
4175: DMGetCoordinates(dm,&coords);
4176: DMGetCoordinates(subdm,&ccoords);
4177: if (coords && !ccoords) {
4178: DMCreateGlobalVector(subdm_coord,&ccoords);
4179: PetscObjectSetName((PetscObject)ccoords,"coordinates");
4180: DMCreateLocalVector(subdm_coord,&clcoords);
4181: PetscObjectSetName((PetscObject)clcoords,"coordinates");
4182: DMCreateDomainDecompositionScatters(dm_coord,1,&subdm_coord,NULL,&scat_i,&scat_g);
4183: VecScatterBegin(scat_i[0],coords,ccoords,INSERT_VALUES,SCATTER_FORWARD);
4184: VecScatterBegin(scat_g[0],coords,clcoords,INSERT_VALUES,SCATTER_FORWARD);
4185: VecScatterEnd(scat_i[0],coords,ccoords,INSERT_VALUES,SCATTER_FORWARD);
4186: VecScatterEnd(scat_g[0],coords,clcoords,INSERT_VALUES,SCATTER_FORWARD);
4187: DMSetCoordinates(subdm,ccoords);
4188: DMSetCoordinatesLocal(subdm,clcoords);
4189: VecScatterDestroy(&scat_i[0]);
4190: VecScatterDestroy(&scat_g[0]);
4191: VecDestroy(&ccoords);
4192: VecDestroy(&clcoords);
4193: PetscFree(scat_i);
4194: PetscFree(scat_g);
4195: }
4196: return(0);
4197: }
4199: /*@
4200: DMGetDimension - Return the topological dimension of the DM
4202: Not collective
4204: Input Parameter:
4205: . dm - The DM
4207: Output Parameter:
4208: . dim - The topological dimension
4210: Level: beginner
4212: .seealso: DMSetDimension(), DMCreate()
4213: @*/
4214: PetscErrorCode DMGetDimension(DM dm, PetscInt *dim)
4215: {
4219: *dim = dm->dim;
4220: return(0);
4221: }
4223: /*@
4224: DMSetDimension - Set the topological dimension of the DM
4226: Collective on dm
4228: Input Parameters:
4229: + dm - The DM
4230: - dim - The topological dimension
4232: Level: beginner
4234: .seealso: DMGetDimension(), DMCreate()
4235: @*/
4236: PetscErrorCode DMSetDimension(DM dm, PetscInt dim)
4237: {
4243: dm->dim = dim;
4244: if (dm->prob->dimEmbed < 0) {PetscDSSetCoordinateDimension(dm->prob, dm->dim);}
4245: return(0);
4246: }
4248: /*@
4249: DMGetDimPoints - Get the half-open interval for all points of a given dimension
4251: Collective on DM
4253: Input Parameters:
4254: + dm - the DM
4255: - dim - the dimension
4257: Output Parameters:
4258: + pStart - The first point of the given dimension
4259: . pEnd - The first point following points of the given dimension
4261: Note:
4262: The points are vertices in the Hasse diagram encoding the topology. This is explained in
4263: http://arxiv.org/abs/0908.4427. If not points exist of this dimension in the storage scheme,
4264: then the interval is empty.
4266: Level: intermediate
4268: .keywords: point, Hasse Diagram, dimension
4269: .seealso: DMPLEX, DMPlexGetDepthStratum(), DMPlexGetHeightStratum()
4270: @*/
4271: PetscErrorCode DMGetDimPoints(DM dm, PetscInt dim, PetscInt *pStart, PetscInt *pEnd)
4272: {
4273: PetscInt d;
4278: DMGetDimension(dm, &d);
4279: if ((dim < 0) || (dim > d)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension %d 1", dim, d);
4280: (*dm->ops->getdimpoints)(dm, dim, pStart, pEnd);
4281: return(0);
4282: }
4284: /*@
4285: DMSetCoordinates - Sets into the DM a global vector that holds the coordinates
4287: Collective on DM
4289: Input Parameters:
4290: + dm - the DM
4291: - c - coordinate vector
4293: Notes:
4294: The coordinates do include those for ghost points, which are in the local vector.
4296: The vector c should be destroyed by the caller.
4298: Level: intermediate
4300: .keywords: distributed array, get, corners, nodes, local indices, coordinates
4301: .seealso: DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM()
4302: @*/
4303: PetscErrorCode DMSetCoordinates(DM dm, Vec c)
4304: {
4310: PetscObjectReference((PetscObject) c);
4311: VecDestroy(&dm->coordinates);
4312: dm->coordinates = c;
4313: VecDestroy(&dm->coordinatesLocal);
4314: DMCoarsenHookAdd(dm,DMRestrictHook_Coordinates,NULL,NULL);
4315: DMSubDomainHookAdd(dm,DMSubDomainHook_Coordinates,NULL,NULL);
4316: return(0);
4317: }
4319: /*@
4320: DMSetCoordinatesLocal - Sets into the DM a local vector that holds the coordinates
4322: Collective on DM
4324: Input Parameters:
4325: + dm - the DM
4326: - c - coordinate vector
4328: Notes:
4329: The coordinates of ghost points can be set using DMSetCoordinates()
4330: followed by DMGetCoordinatesLocal(). This is intended to enable the
4331: setting of ghost coordinates outside of the domain.
4333: The vector c should be destroyed by the caller.
4335: Level: intermediate
4337: .keywords: distributed array, get, corners, nodes, local indices, coordinates
4338: .seealso: DMGetCoordinatesLocal(), DMSetCoordinates(), DMGetCoordinates(), DMGetCoordinateDM()
4339: @*/
4340: PetscErrorCode DMSetCoordinatesLocal(DM dm, Vec c)
4341: {
4347: PetscObjectReference((PetscObject) c);
4348: VecDestroy(&dm->coordinatesLocal);
4350: dm->coordinatesLocal = c;
4352: VecDestroy(&dm->coordinates);
4353: return(0);
4354: }
4356: /*@
4357: DMGetCoordinates - Gets a global vector with the coordinates associated with the DM.
4359: Not Collective
4361: Input Parameter:
4362: . dm - the DM
4364: Output Parameter:
4365: . c - global coordinate vector
4367: Note:
4368: This is a borrowed reference, so the user should NOT destroy this vector
4370: Each process has only the local coordinates (does NOT have the ghost coordinates).
4372: For DMDA, in two and three dimensions coordinates are interlaced (x_0,y_0,x_1,y_1,...)
4373: and (x_0,y_0,z_0,x_1,y_1,z_1...)
4375: Level: intermediate
4377: .keywords: distributed array, get, corners, nodes, local indices, coordinates
4378: .seealso: DMSetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM()
4379: @*/
4380: PetscErrorCode DMGetCoordinates(DM dm, Vec *c)
4381: {
4387: if (!dm->coordinates && dm->coordinatesLocal) {
4388: DM cdm = NULL;
4390: DMGetCoordinateDM(dm, &cdm);
4391: DMCreateGlobalVector(cdm, &dm->coordinates);
4392: PetscObjectSetName((PetscObject) dm->coordinates, "coordinates");
4393: DMLocalToGlobalBegin(cdm, dm->coordinatesLocal, INSERT_VALUES, dm->coordinates);
4394: DMLocalToGlobalEnd(cdm, dm->coordinatesLocal, INSERT_VALUES, dm->coordinates);
4395: }
4396: *c = dm->coordinates;
4397: return(0);
4398: }
4400: /*@
4401: DMGetCoordinatesLocal - Gets a local vector with the coordinates associated with the DM.
4403: Collective on DM
4405: Input Parameter:
4406: . dm - the DM
4408: Output Parameter:
4409: . c - coordinate vector
4411: Note:
4412: This is a borrowed reference, so the user should NOT destroy this vector
4414: Each process has the local and ghost coordinates
4416: For DMDA, in two and three dimensions coordinates are interlaced (x_0,y_0,x_1,y_1,...)
4417: and (x_0,y_0,z_0,x_1,y_1,z_1...)
4419: Level: intermediate
4421: .keywords: distributed array, get, corners, nodes, local indices, coordinates
4422: .seealso: DMSetCoordinatesLocal(), DMGetCoordinates(), DMSetCoordinates(), DMGetCoordinateDM()
4423: @*/
4424: PetscErrorCode DMGetCoordinatesLocal(DM dm, Vec *c)
4425: {
4431: if (!dm->coordinatesLocal && dm->coordinates) {
4432: DM cdm = NULL;
4434: DMGetCoordinateDM(dm, &cdm);
4435: DMCreateLocalVector(cdm, &dm->coordinatesLocal);
4436: PetscObjectSetName((PetscObject) dm->coordinatesLocal, "coordinates");
4437: DMGlobalToLocalBegin(cdm, dm->coordinates, INSERT_VALUES, dm->coordinatesLocal);
4438: DMGlobalToLocalEnd(cdm, dm->coordinates, INSERT_VALUES, dm->coordinatesLocal);
4439: }
4440: *c = dm->coordinatesLocal;
4441: return(0);
4442: }
4444: PetscErrorCode DMGetCoordinateField(DM dm, DMField *field)
4445: {
4451: if (!dm->coordinateField) {
4452: if (dm->ops->createcoordinatefield) {
4453: (*dm->ops->createcoordinatefield)(dm,&dm->coordinateField);
4454: }
4455: }
4456: *field = dm->coordinateField;
4457: return(0);
4458: }
4460: PetscErrorCode DMSetCoordinateField(DM dm, DMField field)
4461: {
4467: PetscObjectReference((PetscObject)field);
4468: DMFieldDestroy(&dm->coordinateField);
4469: dm->coordinateField = field;
4470: return(0);
4471: }
4473: /*@
4474: DMGetCoordinateDM - Gets the DM that prescribes coordinate layout and scatters between global and local coordinates
4476: Collective on DM
4478: Input Parameter:
4479: . dm - the DM
4481: Output Parameter:
4482: . cdm - coordinate DM
4484: Level: intermediate
4486: .keywords: distributed array, get, corners, nodes, local indices, coordinates
4487: .seealso: DMSetCoordinateDM(), DMSetCoordinates(), DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLocal()
4488: @*/
4489: PetscErrorCode DMGetCoordinateDM(DM dm, DM *cdm)
4490: {
4496: if (!dm->coordinateDM) {
4497: if (!dm->ops->createcoordinatedm) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unable to create coordinates for this DM");
4498: (*dm->ops->createcoordinatedm)(dm, &dm->coordinateDM);
4499: }
4500: *cdm = dm->coordinateDM;
4501: return(0);
4502: }
4504: /*@
4505: DMSetCoordinateDM - Sets the DM that prescribes coordinate layout and scatters between global and local coordinates
4507: Logically Collective on DM
4509: Input Parameters:
4510: + dm - the DM
4511: - cdm - coordinate DM
4513: Level: intermediate
4515: .keywords: distributed array, get, corners, nodes, local indices, coordinates
4516: .seealso: DMGetCoordinateDM(), DMSetCoordinates(), DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLocal()
4517: @*/
4518: PetscErrorCode DMSetCoordinateDM(DM dm, DM cdm)
4519: {
4525: PetscObjectReference((PetscObject)cdm);
4526: DMDestroy(&dm->coordinateDM);
4527: dm->coordinateDM = cdm;
4528: return(0);
4529: }
4531: /*@
4532: DMGetCoordinateDim - Retrieve the dimension of embedding space for coordinate values.
4534: Not Collective
4536: Input Parameter:
4537: . dm - The DM object
4539: Output Parameter:
4540: . dim - The embedding dimension
4542: Level: intermediate
4544: .keywords: mesh, coordinates
4545: .seealso: DMSetCoordinateDim(), DMGetCoordinateSection(), DMGetCoordinateDM(), DMGetSection(), DMSetSection()
4546: @*/
4547: PetscErrorCode DMGetCoordinateDim(DM dm, PetscInt *dim)
4548: {
4552: if (dm->dimEmbed == PETSC_DEFAULT) {
4553: dm->dimEmbed = dm->dim;
4554: }
4555: *dim = dm->dimEmbed;
4556: return(0);
4557: }
4559: /*@
4560: DMSetCoordinateDim - Set the dimension of the embedding space for coordinate values.
4562: Not Collective
4564: Input Parameters:
4565: + dm - The DM object
4566: - dim - The embedding dimension
4568: Level: intermediate
4570: .keywords: mesh, coordinates
4571: .seealso: DMGetCoordinateDim(), DMSetCoordinateSection(), DMGetCoordinateSection(), DMGetSection(), DMSetSection()
4572: @*/
4573: PetscErrorCode DMSetCoordinateDim(DM dm, PetscInt dim)
4574: {
4579: dm->dimEmbed = dim;
4580: PetscDSSetCoordinateDimension(dm->prob, dm->dimEmbed);
4581: return(0);
4582: }
4584: /*@
4585: DMGetCoordinateSection - Retrieve the layout of coordinate values over the mesh.
4587: Not Collective
4589: Input Parameter:
4590: . dm - The DM object
4592: Output Parameter:
4593: . section - The PetscSection object
4595: Level: intermediate
4597: .keywords: mesh, coordinates
4598: .seealso: DMGetCoordinateDM(), DMGetSection(), DMSetSection()
4599: @*/
4600: PetscErrorCode DMGetCoordinateSection(DM dm, PetscSection *section)
4601: {
4602: DM cdm;
4608: DMGetCoordinateDM(dm, &cdm);
4609: DMGetSection(cdm, section);
4610: return(0);
4611: }
4613: /*@
4614: DMSetCoordinateSection - Set the layout of coordinate values over the mesh.
4616: Not Collective
4618: Input Parameters:
4619: + dm - The DM object
4620: . dim - The embedding dimension, or PETSC_DETERMINE
4621: - section - The PetscSection object
4623: Level: intermediate
4625: .keywords: mesh, coordinates
4626: .seealso: DMGetCoordinateSection(), DMGetSection(), DMSetSection()
4627: @*/
4628: PetscErrorCode DMSetCoordinateSection(DM dm, PetscInt dim, PetscSection section)
4629: {
4630: DM cdm;
4636: DMGetCoordinateDM(dm, &cdm);
4637: DMSetSection(cdm, section);
4638: if (dim == PETSC_DETERMINE) {
4639: PetscInt d = PETSC_DEFAULT;
4640: PetscInt pStart, pEnd, vStart, vEnd, v, dd;
4642: PetscSectionGetChart(section, &pStart, &pEnd);
4643: DMGetDimPoints(dm, 0, &vStart, &vEnd);
4644: pStart = PetscMax(vStart, pStart);
4645: pEnd = PetscMin(vEnd, pEnd);
4646: for (v = pStart; v < pEnd; ++v) {
4647: PetscSectionGetDof(section, v, &dd);
4648: if (dd) {d = dd; break;}
4649: }
4650: if (d < 0) d = PETSC_DEFAULT;
4651: DMSetCoordinateDim(dm, d);
4652: }
4653: return(0);
4654: }
4656: /*@C
4657: DMGetPeriodicity - Get the description of mesh periodicity
4659: Input Parameters:
4660: . dm - The DM object
4662: Output Parameters:
4663: + per - Whether the DM is periodic or not
4664: . maxCell - Over distances greater than this, we can assume a point has crossed over to another sheet, when trying to localize cell coordinates
4665: . L - If we assume the mesh is a torus, this is the length of each coordinate
4666: - bd - This describes the type of periodicity in each topological dimension
4668: Level: developer
4670: .seealso: DMGetPeriodicity()
4671: @*/
4672: PetscErrorCode DMGetPeriodicity(DM dm, PetscBool *per, const PetscReal **maxCell, const PetscReal **L, const DMBoundaryType **bd)
4673: {
4676: if (per) *per = dm->periodic;
4677: if (L) *L = dm->L;
4678: if (maxCell) *maxCell = dm->maxCell;
4679: if (bd) *bd = dm->bdtype;
4680: return(0);
4681: }
4683: /*@C
4684: DMSetPeriodicity - Set the description of mesh periodicity
4686: Input Parameters:
4687: + dm - The DM object
4688: . per - Whether the DM is periodic or not. If maxCell is not provided, coordinates need to be localized
4689: . maxCell - Over distances greater than this, we can assume a point has crossed over to another sheet, when trying to localize cell coordinates
4690: . L - If we assume the mesh is a torus, this is the length of each coordinate
4691: - bd - This describes the type of periodicity in each topological dimension
4693: Level: developer
4695: .seealso: DMGetPeriodicity()
4696: @*/
4697: PetscErrorCode DMSetPeriodicity(DM dm, PetscBool per, const PetscReal maxCell[], const PetscReal L[], const DMBoundaryType bd[])
4698: {
4699: PetscInt dim, d;
4705: if (maxCell) {
4709: }
4710: PetscFree3(dm->L,dm->maxCell,dm->bdtype);
4711: DMGetDimension(dm, &dim);
4712: if (maxCell) {
4713: PetscMalloc3(dim,&dm->L,dim,&dm->maxCell,dim,&dm->bdtype);
4714: for (d = 0; d < dim; ++d) {dm->L[d] = L[d]; dm->maxCell[d] = maxCell[d]; dm->bdtype[d] = bd[d];}
4715: dm->periodic = PETSC_TRUE;
4716: } else {
4717: dm->periodic = per;
4718: }
4719: return(0);
4720: }
4722: /*@
4723: DMLocalizeCoordinate - If a mesh is periodic (a torus with lengths L_i, some of which can be infinite), project the coordinate onto [0, L_i) in each dimension.
4725: Input Parameters:
4726: + dm - The DM
4727: . in - The input coordinate point (dim numbers)
4728: - endpoint - Include the endpoint L_i
4730: Output Parameter:
4731: . out - The localized coordinate point
4733: Level: developer
4735: .seealso: DMLocalizeCoordinates(), DMLocalizeAddCoordinate()
4736: @*/
4737: PetscErrorCode DMLocalizeCoordinate(DM dm, const PetscScalar in[], PetscBool endpoint, PetscScalar out[])
4738: {
4739: PetscInt dim, d;
4743: DMGetCoordinateDim(dm, &dim);
4744: if (!dm->maxCell) {
4745: for (d = 0; d < dim; ++d) out[d] = in[d];
4746: } else {
4747: if (endpoint) {
4748: for (d = 0; d < dim; ++d) {
4749: if ((PetscAbsReal(PetscRealPart(in[d])/dm->L[d] - PetscFloorReal(PetscRealPart(in[d])/dm->L[d])) < PETSC_SMALL) && (PetscRealPart(in[d])/dm->L[d] > PETSC_SMALL)) {
4750: out[d] = in[d] - dm->L[d]*(PetscFloorReal(PetscRealPart(in[d])/dm->L[d]) - 1);
4751: } else {
4752: out[d] = in[d] - dm->L[d]*PetscFloorReal(PetscRealPart(in[d])/dm->L[d]);
4753: }
4754: }
4755: } else {
4756: for (d = 0; d < dim; ++d) {
4757: out[d] = in[d] - dm->L[d]*PetscFloorReal(PetscRealPart(in[d])/dm->L[d]);
4758: }
4759: }
4760: }
4761: return(0);
4762: }
4764: /*
4765: DMLocalizeCoordinate_Internal - If a mesh is periodic, and the input point is far from the anchor, pick the coordinate sheet of the torus which moves it closer.
4767: Input Parameters:
4768: + dm - The DM
4769: . dim - The spatial dimension
4770: . anchor - The anchor point, the input point can be no more than maxCell away from it
4771: - in - The input coordinate point (dim numbers)
4773: Output Parameter:
4774: . out - The localized coordinate point
4776: Level: developer
4778: Note: This is meant to get a set of coordinates close to each other, as in a cell. The anchor is usually the one of the vertices on a containing cell
4780: .seealso: DMLocalizeCoordinates(), DMLocalizeAddCoordinate()
4781: */
4782: PetscErrorCode DMLocalizeCoordinate_Internal(DM dm, PetscInt dim, const PetscScalar anchor[], const PetscScalar in[], PetscScalar out[])
4783: {
4784: PetscInt d;
4787: if (!dm->maxCell) {
4788: for (d = 0; d < dim; ++d) out[d] = in[d];
4789: } else {
4790: for (d = 0; d < dim; ++d) {
4791: if ((dm->bdtype[d] != DM_BOUNDARY_NONE) && (PetscAbsScalar(anchor[d] - in[d]) > dm->maxCell[d])) {
4792: out[d] = PetscRealPart(anchor[d]) > PetscRealPart(in[d]) ? dm->L[d] + in[d] : in[d] - dm->L[d];
4793: } else {
4794: out[d] = in[d];
4795: }
4796: }
4797: }
4798: return(0);
4799: }
4800: PetscErrorCode DMLocalizeCoordinateReal_Internal(DM dm, PetscInt dim, const PetscReal anchor[], const PetscReal in[], PetscReal out[])
4801: {
4802: PetscInt d;
4805: if (!dm->maxCell) {
4806: for (d = 0; d < dim; ++d) out[d] = in[d];
4807: } else {
4808: for (d = 0; d < dim; ++d) {
4809: if ((dm->bdtype[d] != DM_BOUNDARY_NONE) && (PetscAbsReal(anchor[d] - in[d]) > dm->maxCell[d])) {
4810: out[d] = anchor[d] > in[d] ? dm->L[d] + in[d] : in[d] - dm->L[d];
4811: } else {
4812: out[d] = in[d];
4813: }
4814: }
4815: }
4816: return(0);
4817: }
4819: /*
4820: DMLocalizeAddCoordinate_Internal - If a mesh is periodic, and the input point is far from the anchor, pick the coordinate sheet of the torus which moves it closer.
4822: Input Parameters:
4823: + dm - The DM
4824: . dim - The spatial dimension
4825: . anchor - The anchor point, the input point can be no more than maxCell away from it
4826: . in - The input coordinate delta (dim numbers)
4827: - out - The input coordinate point (dim numbers)
4829: Output Parameter:
4830: . out - The localized coordinate in + out
4832: Level: developer
4834: Note: This is meant to get a set of coordinates close to each other, as in a cell. The anchor is usually the one of the vertices on a containing cell
4836: .seealso: DMLocalizeCoordinates(), DMLocalizeCoordinate()
4837: */
4838: PetscErrorCode DMLocalizeAddCoordinate_Internal(DM dm, PetscInt dim, const PetscScalar anchor[], const PetscScalar in[], PetscScalar out[])
4839: {
4840: PetscInt d;
4843: if (!dm->maxCell) {
4844: for (d = 0; d < dim; ++d) out[d] += in[d];
4845: } else {
4846: for (d = 0; d < dim; ++d) {
4847: if ((dm->bdtype[d] != DM_BOUNDARY_NONE) && (PetscAbsScalar(anchor[d] - in[d]) > dm->maxCell[d])) {
4848: out[d] += PetscRealPart(anchor[d]) > PetscRealPart(in[d]) ? dm->L[d] + in[d] : in[d] - dm->L[d];
4849: } else {
4850: out[d] += in[d];
4851: }
4852: }
4853: }
4854: return(0);
4855: }
4857: /*@
4858: DMGetCoordinatesLocalized - Check if the DM coordinates have been localized for cells
4860: Input Parameter:
4861: . dm - The DM
4863: Output Parameter:
4864: areLocalized - True if localized
4866: Level: developer
4868: .seealso: DMLocalizeCoordinates()
4869: @*/
4870: PetscErrorCode DMGetCoordinatesLocalized(DM dm,PetscBool *areLocalized)
4871: {
4872: DM cdm;
4873: PetscSection coordSection;
4874: PetscInt cStart, cEnd, sStart, sEnd, c, dof;
4875: PetscBool isPlex, alreadyLocalized;
4882: *areLocalized = PETSC_FALSE;
4883: if (!dm->periodic) return(0); /* This is a hideous optimization hack! */
4885: /* We need some generic way of refering to cells/vertices */
4886: DMGetCoordinateDM(dm, &cdm);
4887: DMGetCoordinateSection(dm, &coordSection);
4888: PetscObjectTypeCompare((PetscObject) cdm, DMPLEX, &isPlex);
4889: if (!isPlex) SETERRQ(PetscObjectComm((PetscObject) cdm), PETSC_ERR_ARG_WRONG, "Coordinate localization requires a DMPLEX coordinate DM");
4891: DMPlexGetHeightStratum(cdm, 0, &cStart, &cEnd);
4892: PetscSectionGetChart(coordSection, &sStart, &sEnd);
4893: alreadyLocalized = PETSC_FALSE;
4894: for (c = cStart; c < cEnd; ++c) {
4895: if (c < sStart || c >= sEnd) continue;
4896: PetscSectionGetDof(coordSection, c, &dof);
4897: if (dof) { alreadyLocalized = PETSC_TRUE; break; }
4898: }
4899: MPIU_Allreduce(&alreadyLocalized,areLocalized,1,MPIU_BOOL,MPI_LOR,PetscObjectComm((PetscObject)dm));
4900: return(0);
4901: }
4904: /*@
4905: DMLocalizeCoordinates - If a mesh is periodic, create local coordinates for cells having periodic faces
4907: Input Parameter:
4908: . dm - The DM
4910: Level: developer
4912: .seealso: DMLocalizeCoordinate(), DMLocalizeAddCoordinate()
4913: @*/
4914: PetscErrorCode DMLocalizeCoordinates(DM dm)
4915: {
4916: DM cdm;
4917: PetscSection coordSection, cSection;
4918: Vec coordinates, cVec;
4919: PetscScalar *coords, *coords2, *anchor, *localized;
4920: PetscInt Nc, vStart, vEnd, v, sStart, sEnd, newStart = PETSC_MAX_INT, newEnd = PETSC_MIN_INT, dof, d, off, off2, bs, coordSize;
4921: PetscBool alreadyLocalized, alreadyLocalizedGlobal;
4922: PetscInt maxHeight = 0, h;
4923: PetscInt *pStart = NULL, *pEnd = NULL;
4928: if (!dm->periodic) return(0);
4929: DMGetCoordinatesLocalized(dm, &alreadyLocalized);
4930: if (alreadyLocalized) return(0);
4932: /* We need some generic way of refering to cells/vertices */
4933: DMGetCoordinateDM(dm, &cdm);
4934: {
4935: PetscBool isplex;
4937: PetscObjectTypeCompare((PetscObject) cdm, DMPLEX, &isplex);
4938: if (isplex) {
4939: DMPlexGetDepthStratum(cdm, 0, &vStart, &vEnd);
4940: DMPlexGetMaxProjectionHeight(cdm,&maxHeight);
4941: DMGetWorkArray(dm,2*(maxHeight + 1),MPIU_INT,&pStart);
4942: pEnd = &pStart[maxHeight + 1];
4943: newStart = vStart;
4944: newEnd = vEnd;
4945: for (h = 0; h <= maxHeight; h++) {
4946: DMPlexGetHeightStratum(cdm, h, &pStart[h], &pEnd[h]);
4947: newStart = PetscMin(newStart,pStart[h]);
4948: newEnd = PetscMax(newEnd,pEnd[h]);
4949: }
4950: } else SETERRQ(PetscObjectComm((PetscObject) cdm), PETSC_ERR_ARG_WRONG, "Coordinate localization requires a DMPLEX coordinate DM");
4951: }
4952: DMGetCoordinatesLocal(dm, &coordinates);
4953: if (!coordinates) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Missing local coordinates vector");
4954: DMGetCoordinateSection(dm, &coordSection);
4955: VecGetBlockSize(coordinates, &bs);
4956: PetscSectionGetChart(coordSection,&sStart,&sEnd);
4958: PetscSectionCreate(PetscObjectComm((PetscObject) dm), &cSection);
4959: PetscSectionSetNumFields(cSection, 1);
4960: PetscSectionGetFieldComponents(coordSection, 0, &Nc);
4961: PetscSectionSetFieldComponents(cSection, 0, Nc);
4962: PetscSectionSetChart(cSection, newStart, newEnd);
4964: DMGetWorkArray(dm, 2 * bs, MPIU_SCALAR, &anchor);
4965: localized = &anchor[bs];
4966: alreadyLocalized = alreadyLocalizedGlobal = PETSC_TRUE;
4967: for (h = 0; h <= maxHeight; h++) {
4968: PetscInt cStart = pStart[h], cEnd = pEnd[h], c;
4970: for (c = cStart; c < cEnd; ++c) {
4971: PetscScalar *cellCoords = NULL;
4972: PetscInt b;
4974: if (c < sStart || c >= sEnd) alreadyLocalized = PETSC_FALSE;
4975: DMPlexVecGetClosure(cdm, coordSection, coordinates, c, &dof, &cellCoords);
4976: for (b = 0; b < bs; ++b) anchor[b] = cellCoords[b];
4977: for (d = 0; d < dof/bs; ++d) {
4978: DMLocalizeCoordinate_Internal(dm, bs, anchor, &cellCoords[d*bs], localized);
4979: for (b = 0; b < bs; b++) {
4980: if (cellCoords[d*bs + b] != localized[b]) break;
4981: }
4982: if (b < bs) break;
4983: }
4984: if (d < dof/bs) {
4985: if (c >= sStart && c < sEnd) {
4986: PetscInt cdof;
4988: PetscSectionGetDof(coordSection, c, &cdof);
4989: if (cdof != dof) alreadyLocalized = PETSC_FALSE;
4990: }
4991: PetscSectionSetDof(cSection, c, dof);
4992: PetscSectionSetFieldDof(cSection, c, 0, dof);
4993: }
4994: DMPlexVecRestoreClosure(cdm, coordSection, coordinates, c, &dof, &cellCoords);
4995: }
4996: }
4997: MPI_Allreduce(&alreadyLocalized,&alreadyLocalizedGlobal,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)dm));
4998: if (alreadyLocalizedGlobal) {
4999: DMRestoreWorkArray(dm, 2 * bs, MPIU_SCALAR, &anchor);
5000: PetscSectionDestroy(&cSection);
5001: DMRestoreWorkArray(dm,2*(maxHeight + 1),MPIU_INT,&pStart);
5002: return(0);
5003: }
5004: for (v = vStart; v < vEnd; ++v) {
5005: PetscSectionGetDof(coordSection, v, &dof);
5006: PetscSectionSetDof(cSection, v, dof);
5007: PetscSectionSetFieldDof(cSection, v, 0, dof);
5008: }
5009: PetscSectionSetUp(cSection);
5010: PetscSectionGetStorageSize(cSection, &coordSize);
5011: VecCreate(PETSC_COMM_SELF, &cVec);
5012: PetscObjectSetName((PetscObject)cVec,"coordinates");
5013: VecSetBlockSize(cVec, bs);
5014: VecSetSizes(cVec, coordSize, PETSC_DETERMINE);
5015: VecSetType(cVec, VECSTANDARD);
5016: VecGetArrayRead(coordinates, (const PetscScalar**)&coords);
5017: VecGetArray(cVec, &coords2);
5018: for (v = vStart; v < vEnd; ++v) {
5019: PetscSectionGetDof(coordSection, v, &dof);
5020: PetscSectionGetOffset(coordSection, v, &off);
5021: PetscSectionGetOffset(cSection, v, &off2);
5022: for (d = 0; d < dof; ++d) coords2[off2+d] = coords[off+d];
5023: }
5024: for (h = 0; h <= maxHeight; h++) {
5025: PetscInt cStart = pStart[h], cEnd = pEnd[h], c;
5027: for (c = cStart; c < cEnd; ++c) {
5028: PetscScalar *cellCoords = NULL;
5029: PetscInt b, cdof;
5031: PetscSectionGetDof(cSection,c,&cdof);
5032: if (!cdof) continue;
5033: DMPlexVecGetClosure(cdm, coordSection, coordinates, c, &dof, &cellCoords);
5034: PetscSectionGetOffset(cSection, c, &off2);
5035: for (b = 0; b < bs; ++b) anchor[b] = cellCoords[b];
5036: for (d = 0; d < dof/bs; ++d) {DMLocalizeCoordinate_Internal(dm, bs, anchor, &cellCoords[d*bs], &coords2[off2+d*bs]);}
5037: DMPlexVecRestoreClosure(cdm, coordSection, coordinates, c, &dof, &cellCoords);
5038: }
5039: }
5040: DMRestoreWorkArray(dm, 2 * bs, MPIU_SCALAR, &anchor);
5041: DMRestoreWorkArray(dm,2*(maxHeight + 1),MPIU_INT,&pStart);
5042: VecRestoreArrayRead(coordinates, (const PetscScalar**)&coords);
5043: VecRestoreArray(cVec, &coords2);
5044: DMSetCoordinateSection(dm, PETSC_DETERMINE, cSection);
5045: DMSetCoordinatesLocal(dm, cVec);
5046: VecDestroy(&cVec);
5047: PetscSectionDestroy(&cSection);
5048: return(0);
5049: }
5051: /*@
5052: DMLocatePoints - Locate the points in v in the mesh and return a PetscSF of the containing cells
5054: Collective on Vec v (see explanation below)
5056: Input Parameters:
5057: + dm - The DM
5058: . v - The Vec of points
5059: . ltype - The type of point location, e.g. DM_POINTLOCATION_NONE or DM_POINTLOCATION_NEAREST
5060: - cells - Points to either NULL, or a PetscSF with guesses for which cells contain each point.
5062: Output Parameter:
5063: + v - The Vec of points, which now contains the nearest mesh points to the given points if DM_POINTLOCATION_NEAREST is used
5064: - cells - The PetscSF containing the ranks and local indices of the containing points.
5067: Level: developer
5069: Notes:
5070: To do a search of the local cells of the mesh, v should have PETSC_COMM_SELF as its communicator.
5071: To do a search of all the cells in the distributed mesh, v should have the same communicator as dm.
5073: If *cellSF is NULL on input, a PetscSF will be created.
5074: If *cellSF is not NULL on input, it should point to an existing PetscSF, whose graph will be used as initial guesses.
5076: An array that maps each point to its containing cell can be obtained with
5078: $ const PetscSFNode *cells;
5079: $ PetscInt nFound;
5080: $ const PetscInt *found;
5081: $
5082: $ PetscSFGetGraph(cellSF,NULL,&nFound,&found,&cells);
5084: Where cells[i].rank is the rank of the cell containing point found[i] (or i if found == NULL), and cells[i].index is
5085: the index of the cell in its rank's local numbering.
5087: .keywords: point location, mesh
5088: .seealso: DMSetCoordinates(), DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLocal(), DMPointLocationType
5089: @*/
5090: PetscErrorCode DMLocatePoints(DM dm, Vec v, DMPointLocationType ltype, PetscSF *cellSF)
5091: {
5098: if (*cellSF) {
5099: PetscMPIInt result;
5102: MPI_Comm_compare(PetscObjectComm((PetscObject)v),PetscObjectComm((PetscObject)*cellSF),&result);
5103: if (result != MPI_IDENT && result != MPI_CONGRUENT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"cellSF must have a communicator congruent to v's");
5104: } else {
5105: PetscSFCreate(PetscObjectComm((PetscObject)v),cellSF);
5106: }
5107: PetscLogEventBegin(DM_LocatePoints,dm,0,0,0);
5108: if (dm->ops->locatepoints) {
5109: (*dm->ops->locatepoints)(dm,v,ltype,*cellSF);
5110: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Point location not available for this DM");
5111: PetscLogEventEnd(DM_LocatePoints,dm,0,0,0);
5112: return(0);
5113: }
5115: /*@
5116: DMGetOutputDM - Retrieve the DM associated with the layout for output
5118: Input Parameter:
5119: . dm - The original DM
5121: Output Parameter:
5122: . odm - The DM which provides the layout for output
5124: Level: intermediate
5126: .seealso: VecView(), DMGetGlobalSection()
5127: @*/
5128: PetscErrorCode DMGetOutputDM(DM dm, DM *odm)
5129: {
5130: PetscSection section;
5131: PetscBool hasConstraints, ghasConstraints;
5137: DMGetSection(dm, §ion);
5138: PetscSectionHasConstraints(section, &hasConstraints);
5139: MPI_Allreduce(&hasConstraints, &ghasConstraints, 1, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject) dm));
5140: if (!ghasConstraints) {
5141: *odm = dm;
5142: return(0);
5143: }
5144: if (!dm->dmBC) {
5145: PetscDS ds;
5146: PetscSection newSection, gsection;
5147: PetscSF sf;
5149: DMClone(dm, &dm->dmBC);
5150: DMGetDS(dm, &ds);
5151: DMSetDS(dm->dmBC, ds);
5152: PetscSectionClone(section, &newSection);
5153: DMSetSection(dm->dmBC, newSection);
5154: PetscSectionDestroy(&newSection);
5155: DMGetPointSF(dm->dmBC, &sf);
5156: PetscSectionCreateGlobalSection(section, sf, PETSC_TRUE, PETSC_FALSE, &gsection);
5157: DMSetGlobalSection(dm->dmBC, gsection);
5158: PetscSectionDestroy(&gsection);
5159: }
5160: *odm = dm->dmBC;
5161: return(0);
5162: }
5164: /*@
5165: DMGetOutputSequenceNumber - Retrieve the sequence number/value for output
5167: Input Parameter:
5168: . dm - The original DM
5170: Output Parameters:
5171: + num - The output sequence number
5172: - val - The output sequence value
5174: Level: intermediate
5176: Note: This is intended for output that should appear in sequence, for instance
5177: a set of timesteps in an HDF5 file, or a set of realizations of a stochastic system.
5179: .seealso: VecView()
5180: @*/
5181: PetscErrorCode DMGetOutputSequenceNumber(DM dm, PetscInt *num, PetscReal *val)
5182: {
5187: return(0);
5188: }
5190: /*@
5191: DMSetOutputSequenceNumber - Set the sequence number/value for output
5193: Input Parameters:
5194: + dm - The original DM
5195: . num - The output sequence number
5196: - val - The output sequence value
5198: Level: intermediate
5200: Note: This is intended for output that should appear in sequence, for instance
5201: a set of timesteps in an HDF5 file, or a set of realizations of a stochastic system.
5203: .seealso: VecView()
5204: @*/
5205: PetscErrorCode DMSetOutputSequenceNumber(DM dm, PetscInt num, PetscReal val)
5206: {
5209: dm->outputSequenceNum = num;
5210: dm->outputSequenceVal = val;
5211: return(0);
5212: }
5214: /*@C
5215: DMOutputSequenceLoad - Retrieve the sequence value from a Viewer
5217: Input Parameters:
5218: + dm - The original DM
5219: . name - The sequence name
5220: - num - The output sequence number
5222: Output Parameter:
5223: . val - The output sequence value
5225: Level: intermediate
5227: Note: This is intended for output that should appear in sequence, for instance
5228: a set of timesteps in an HDF5 file, or a set of realizations of a stochastic system.
5230: .seealso: DMGetOutputSequenceNumber(), DMSetOutputSequenceNumber(), VecView()
5231: @*/
5232: PetscErrorCode DMOutputSequenceLoad(DM dm, PetscViewer viewer, const char *name, PetscInt num, PetscReal *val)
5233: {
5234: PetscBool ishdf5;
5241: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERHDF5, &ishdf5);
5242: if (ishdf5) {
5243: #if defined(PETSC_HAVE_HDF5)
5244: PetscScalar value;
5246: DMSequenceLoad_HDF5_Internal(dm, name, num, &value, viewer);
5247: *val = PetscRealPart(value);
5248: #endif
5249: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid viewer; open viewer with PetscViewerHDF5Open()");
5250: return(0);
5251: }
5253: /*@
5254: DMGetUseNatural - Get the flag for creating a mapping to the natural order on distribution
5256: Not collective
5258: Input Parameter:
5259: . dm - The DM
5261: Output Parameter:
5262: . useNatural - The flag to build the mapping to a natural order during distribution
5264: Level: beginner
5266: .seealso: DMSetUseNatural(), DMCreate()
5267: @*/
5268: PetscErrorCode DMGetUseNatural(DM dm, PetscBool *useNatural)
5269: {
5273: *useNatural = dm->useNatural;
5274: return(0);
5275: }
5277: /*@
5278: DMSetUseNatural - Set the flag for creating a mapping to the natural order on distribution
5280: Collective on dm
5282: Input Parameters:
5283: + dm - The DM
5284: - useNatural - The flag to build the mapping to a natural order during distribution
5286: Level: beginner
5288: .seealso: DMGetUseNatural(), DMCreate()
5289: @*/
5290: PetscErrorCode DMSetUseNatural(DM dm, PetscBool useNatural)
5291: {
5295: dm->useNatural = useNatural;
5296: return(0);
5297: }
5300: /*@C
5301: DMCreateLabel - Create a label of the given name if it does not already exist
5303: Not Collective
5305: Input Parameters:
5306: + dm - The DM object
5307: - name - The label name
5309: Level: intermediate
5311: .keywords: mesh
5312: .seealso: DMLabelCreate(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5313: @*/
5314: PetscErrorCode DMCreateLabel(DM dm, const char name[])
5315: {
5316: DMLabelLink next = dm->labels->next;
5317: PetscBool flg = PETSC_FALSE;
5323: while (next) {
5324: PetscStrcmp(name, next->label->name, &flg);
5325: if (flg) break;
5326: next = next->next;
5327: }
5328: if (!flg) {
5329: DMLabelLink tmpLabel;
5331: PetscCalloc1(1, &tmpLabel);
5332: DMLabelCreate(name, &tmpLabel->label);
5333: tmpLabel->output = PETSC_TRUE;
5334: tmpLabel->next = dm->labels->next;
5335: dm->labels->next = tmpLabel;
5336: }
5337: return(0);
5338: }
5340: /*@C
5341: DMGetLabelValue - Get the value in a Sieve Label for the given point, with 0 as the default
5343: Not Collective
5345: Input Parameters:
5346: + dm - The DM object
5347: . name - The label name
5348: - point - The mesh point
5350: Output Parameter:
5351: . value - The label value for this point, or -1 if the point is not in the label
5353: Level: beginner
5355: .keywords: mesh
5356: .seealso: DMLabelGetValue(), DMSetLabelValue(), DMGetStratumIS()
5357: @*/
5358: PetscErrorCode DMGetLabelValue(DM dm, const char name[], PetscInt point, PetscInt *value)
5359: {
5360: DMLabel label;
5366: DMGetLabel(dm, name, &label);
5367: if (!label) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No label named %s was found", name);
5368: DMLabelGetValue(label, point, value);
5369: return(0);
5370: }
5372: /*@C
5373: DMSetLabelValue - Add a point to a Sieve Label with given value
5375: Not Collective
5377: Input Parameters:
5378: + dm - The DM object
5379: . name - The label name
5380: . point - The mesh point
5381: - value - The label value for this point
5383: Output Parameter:
5385: Level: beginner
5387: .keywords: mesh
5388: .seealso: DMLabelSetValue(), DMGetStratumIS(), DMClearLabelValue()
5389: @*/
5390: PetscErrorCode DMSetLabelValue(DM dm, const char name[], PetscInt point, PetscInt value)
5391: {
5392: DMLabel label;
5398: DMGetLabel(dm, name, &label);
5399: if (!label) {
5400: DMCreateLabel(dm, name);
5401: DMGetLabel(dm, name, &label);
5402: }
5403: DMLabelSetValue(label, point, value);
5404: return(0);
5405: }
5407: /*@C
5408: DMClearLabelValue - Remove a point from a Sieve Label with given value
5410: Not Collective
5412: Input Parameters:
5413: + dm - The DM object
5414: . name - The label name
5415: . point - The mesh point
5416: - value - The label value for this point
5418: Output Parameter:
5420: Level: beginner
5422: .keywords: mesh
5423: .seealso: DMLabelClearValue(), DMSetLabelValue(), DMGetStratumIS()
5424: @*/
5425: PetscErrorCode DMClearLabelValue(DM dm, const char name[], PetscInt point, PetscInt value)
5426: {
5427: DMLabel label;
5433: DMGetLabel(dm, name, &label);
5434: if (!label) return(0);
5435: DMLabelClearValue(label, point, value);
5436: return(0);
5437: }
5439: /*@C
5440: DMGetLabelSize - Get the number of different integer ids in a Label
5442: Not Collective
5444: Input Parameters:
5445: + dm - The DM object
5446: - name - The label name
5448: Output Parameter:
5449: . size - The number of different integer ids, or 0 if the label does not exist
5451: Level: beginner
5453: .keywords: mesh
5454: .seealso: DMLabeGetNumValues(), DMSetLabelValue()
5455: @*/
5456: PetscErrorCode DMGetLabelSize(DM dm, const char name[], PetscInt *size)
5457: {
5458: DMLabel label;
5465: DMGetLabel(dm, name, &label);
5466: *size = 0;
5467: if (!label) return(0);
5468: DMLabelGetNumValues(label, size);
5469: return(0);
5470: }
5472: /*@C
5473: DMGetLabelIdIS - Get the integer ids in a label
5475: Not Collective
5477: Input Parameters:
5478: + mesh - The DM object
5479: - name - The label name
5481: Output Parameter:
5482: . ids - The integer ids, or NULL if the label does not exist
5484: Level: beginner
5486: .keywords: mesh
5487: .seealso: DMLabelGetValueIS(), DMGetLabelSize()
5488: @*/
5489: PetscErrorCode DMGetLabelIdIS(DM dm, const char name[], IS *ids)
5490: {
5491: DMLabel label;
5498: DMGetLabel(dm, name, &label);
5499: *ids = NULL;
5500: if (label) {
5501: DMLabelGetValueIS(label, ids);
5502: } else {
5503: /* returning an empty IS */
5504: ISCreateGeneral(PETSC_COMM_SELF,0,NULL,PETSC_USE_POINTER,ids);
5505: }
5506: return(0);
5507: }
5509: /*@C
5510: DMGetStratumSize - Get the number of points in a label stratum
5512: Not Collective
5514: Input Parameters:
5515: + dm - The DM object
5516: . name - The label name
5517: - value - The stratum value
5519: Output Parameter:
5520: . size - The stratum size
5522: Level: beginner
5524: .keywords: mesh
5525: .seealso: DMLabelGetStratumSize(), DMGetLabelSize(), DMGetLabelIds()
5526: @*/
5527: PetscErrorCode DMGetStratumSize(DM dm, const char name[], PetscInt value, PetscInt *size)
5528: {
5529: DMLabel label;
5536: DMGetLabel(dm, name, &label);
5537: *size = 0;
5538: if (!label) return(0);
5539: DMLabelGetStratumSize(label, value, size);
5540: return(0);
5541: }
5543: /*@C
5544: DMGetStratumIS - Get the points in a label stratum
5546: Not Collective
5548: Input Parameters:
5549: + dm - The DM object
5550: . name - The label name
5551: - value - The stratum value
5553: Output Parameter:
5554: . points - The stratum points, or NULL if the label does not exist or does not have that value
5556: Level: beginner
5558: .keywords: mesh
5559: .seealso: DMLabelGetStratumIS(), DMGetStratumSize()
5560: @*/
5561: PetscErrorCode DMGetStratumIS(DM dm, const char name[], PetscInt value, IS *points)
5562: {
5563: DMLabel label;
5570: DMGetLabel(dm, name, &label);
5571: *points = NULL;
5572: if (!label) return(0);
5573: DMLabelGetStratumIS(label, value, points);
5574: return(0);
5575: }
5577: /*@C
5578: DMSetStratumIS - Set the points in a label stratum
5580: Not Collective
5582: Input Parameters:
5583: + dm - The DM object
5584: . name - The label name
5585: . value - The stratum value
5586: - points - The stratum points
5588: Level: beginner
5590: .keywords: mesh
5591: .seealso: DMLabelSetStratumIS(), DMGetStratumSize()
5592: @*/
5593: PetscErrorCode DMSetStratumIS(DM dm, const char name[], PetscInt value, IS points)
5594: {
5595: DMLabel label;
5602: DMGetLabel(dm, name, &label);
5603: if (!label) return(0);
5604: DMLabelSetStratumIS(label, value, points);
5605: return(0);
5606: }
5608: /*@C
5609: DMClearLabelStratum - Remove all points from a stratum from a Sieve Label
5611: Not Collective
5613: Input Parameters:
5614: + dm - The DM object
5615: . name - The label name
5616: - value - The label value for this point
5618: Output Parameter:
5620: Level: beginner
5622: .keywords: mesh
5623: .seealso: DMLabelClearStratum(), DMSetLabelValue(), DMGetStratumIS(), DMClearLabelValue()
5624: @*/
5625: PetscErrorCode DMClearLabelStratum(DM dm, const char name[], PetscInt value)
5626: {
5627: DMLabel label;
5633: DMGetLabel(dm, name, &label);
5634: if (!label) return(0);
5635: DMLabelClearStratum(label, value);
5636: return(0);
5637: }
5639: /*@
5640: DMGetNumLabels - Return the number of labels defined by the mesh
5642: Not Collective
5644: Input Parameter:
5645: . dm - The DM object
5647: Output Parameter:
5648: . numLabels - the number of Labels
5650: Level: intermediate
5652: .keywords: mesh
5653: .seealso: DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5654: @*/
5655: PetscErrorCode DMGetNumLabels(DM dm, PetscInt *numLabels)
5656: {
5657: DMLabelLink next = dm->labels->next;
5658: PetscInt n = 0;
5663: while (next) {++n; next = next->next;}
5664: *numLabels = n;
5665: return(0);
5666: }
5668: /*@C
5669: DMGetLabelName - Return the name of nth label
5671: Not Collective
5673: Input Parameters:
5674: + dm - The DM object
5675: - n - the label number
5677: Output Parameter:
5678: . name - the label name
5680: Level: intermediate
5682: .keywords: mesh
5683: .seealso: DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5684: @*/
5685: PetscErrorCode DMGetLabelName(DM dm, PetscInt n, const char **name)
5686: {
5687: DMLabelLink next = dm->labels->next;
5688: PetscInt l = 0;
5693: while (next) {
5694: if (l == n) {
5695: *name = next->label->name;
5696: return(0);
5697: }
5698: ++l;
5699: next = next->next;
5700: }
5701: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label %D does not exist in this DM", n);
5702: }
5704: /*@C
5705: DMHasLabel - Determine whether the mesh has a label of a given name
5707: Not Collective
5709: Input Parameters:
5710: + dm - The DM object
5711: - name - The label name
5713: Output Parameter:
5714: . hasLabel - PETSC_TRUE if the label is present
5716: Level: intermediate
5718: .keywords: mesh
5719: .seealso: DMCreateLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5720: @*/
5721: PetscErrorCode DMHasLabel(DM dm, const char name[], PetscBool *hasLabel)
5722: {
5723: DMLabelLink next = dm->labels->next;
5730: *hasLabel = PETSC_FALSE;
5731: while (next) {
5732: PetscStrcmp(name, next->label->name, hasLabel);
5733: if (*hasLabel) break;
5734: next = next->next;
5735: }
5736: return(0);
5737: }
5739: /*@C
5740: DMGetLabel - Return the label of a given name, or NULL
5742: Not Collective
5744: Input Parameters:
5745: + dm - The DM object
5746: - name - The label name
5748: Output Parameter:
5749: . label - The DMLabel, or NULL if the label is absent
5751: Level: intermediate
5753: .keywords: mesh
5754: .seealso: DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5755: @*/
5756: PetscErrorCode DMGetLabel(DM dm, const char name[], DMLabel *label)
5757: {
5758: DMLabelLink next = dm->labels->next;
5759: PetscBool hasLabel;
5766: *label = NULL;
5767: while (next) {
5768: PetscStrcmp(name, next->label->name, &hasLabel);
5769: if (hasLabel) {
5770: *label = next->label;
5771: break;
5772: }
5773: next = next->next;
5774: }
5775: return(0);
5776: }
5778: /*@C
5779: DMGetLabelByNum - Return the nth label
5781: Not Collective
5783: Input Parameters:
5784: + dm - The DM object
5785: - n - the label number
5787: Output Parameter:
5788: . label - the label
5790: Level: intermediate
5792: .keywords: mesh
5793: .seealso: DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5794: @*/
5795: PetscErrorCode DMGetLabelByNum(DM dm, PetscInt n, DMLabel *label)
5796: {
5797: DMLabelLink next = dm->labels->next;
5798: PetscInt l = 0;
5803: while (next) {
5804: if (l == n) {
5805: *label = next->label;
5806: return(0);
5807: }
5808: ++l;
5809: next = next->next;
5810: }
5811: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label %D does not exist in this DM", n);
5812: }
5814: /*@C
5815: DMAddLabel - Add the label to this mesh
5817: Not Collective
5819: Input Parameters:
5820: + dm - The DM object
5821: - label - The DMLabel
5823: Level: developer
5825: .keywords: mesh
5826: .seealso: DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5827: @*/
5828: PetscErrorCode DMAddLabel(DM dm, DMLabel label)
5829: {
5830: DMLabelLink tmpLabel;
5831: PetscBool hasLabel;
5836: DMHasLabel(dm, label->name, &hasLabel);
5837: if (hasLabel) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label %s already exists in this DM", label->name);
5838: PetscCalloc1(1, &tmpLabel);
5839: tmpLabel->label = label;
5840: tmpLabel->output = PETSC_TRUE;
5841: tmpLabel->next = dm->labels->next;
5842: dm->labels->next = tmpLabel;
5843: return(0);
5844: }
5846: /*@C
5847: DMRemoveLabel - Remove the label from this mesh
5849: Not Collective
5851: Input Parameters:
5852: + dm - The DM object
5853: - name - The label name
5855: Output Parameter:
5856: . label - The DMLabel, or NULL if the label is absent
5858: Level: developer
5860: .keywords: mesh
5861: .seealso: DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5862: @*/
5863: PetscErrorCode DMRemoveLabel(DM dm, const char name[], DMLabel *label)
5864: {
5865: DMLabelLink next = dm->labels->next;
5866: DMLabelLink last = NULL;
5867: PetscBool hasLabel;
5872: DMHasLabel(dm, name, &hasLabel);
5873: *label = NULL;
5874: if (!hasLabel) return(0);
5875: while (next) {
5876: PetscStrcmp(name, next->label->name, &hasLabel);
5877: if (hasLabel) {
5878: if (last) last->next = next->next;
5879: else dm->labels->next = next->next;
5880: next->next = NULL;
5881: *label = next->label;
5882: PetscStrcmp(name, "depth", &hasLabel);
5883: if (hasLabel) {
5884: dm->depthLabel = NULL;
5885: }
5886: PetscFree(next);
5887: break;
5888: }
5889: last = next;
5890: next = next->next;
5891: }
5892: return(0);
5893: }
5895: /*@C
5896: DMGetLabelOutput - Get the output flag for a given label
5898: Not Collective
5900: Input Parameters:
5901: + dm - The DM object
5902: - name - The label name
5904: Output Parameter:
5905: . output - The flag for output
5907: Level: developer
5909: .keywords: mesh
5910: .seealso: DMSetLabelOutput(), DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5911: @*/
5912: PetscErrorCode DMGetLabelOutput(DM dm, const char name[], PetscBool *output)
5913: {
5914: DMLabelLink next = dm->labels->next;
5921: while (next) {
5922: PetscBool flg;
5924: PetscStrcmp(name, next->label->name, &flg);
5925: if (flg) {*output = next->output; return(0);}
5926: next = next->next;
5927: }
5928: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No label named %s was present in this dm", name);
5929: }
5931: /*@C
5932: DMSetLabelOutput - Set the output flag for a given label
5934: Not Collective
5936: Input Parameters:
5937: + dm - The DM object
5938: . name - The label name
5939: - output - The flag for output
5941: Level: developer
5943: .keywords: mesh
5944: .seealso: DMGetLabelOutput(), DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5945: @*/
5946: PetscErrorCode DMSetLabelOutput(DM dm, const char name[], PetscBool output)
5947: {
5948: DMLabelLink next = dm->labels->next;
5954: while (next) {
5955: PetscBool flg;
5957: PetscStrcmp(name, next->label->name, &flg);
5958: if (flg) {next->output = output; return(0);}
5959: next = next->next;
5960: }
5961: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No label named %s was present in this dm", name);
5962: }
5965: /*@
5966: DMCopyLabels - Copy labels from one mesh to another with a superset of the points
5968: Collective on DM
5970: Input Parameter:
5971: . dmA - The DM object with initial labels
5973: Output Parameter:
5974: . dmB - The DM object with copied labels
5976: Level: intermediate
5978: Note: This is typically used when interpolating or otherwise adding to a mesh
5980: .keywords: mesh
5981: .seealso: DMCopyCoordinates(), DMGetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM(), DMGetCoordinateSection()
5982: @*/
5983: PetscErrorCode DMCopyLabels(DM dmA, DM dmB)
5984: {
5985: PetscInt numLabels, l;
5989: if (dmA == dmB) return(0);
5990: DMGetNumLabels(dmA, &numLabels);
5991: for (l = 0; l < numLabels; ++l) {
5992: DMLabel label, labelNew;
5993: const char *name;
5994: PetscBool flg;
5996: DMGetLabelName(dmA, l, &name);
5997: PetscStrcmp(name, "depth", &flg);
5998: if (flg) continue;
5999: PetscStrcmp(name, "dim", &flg);
6000: if (flg) continue;
6001: DMGetLabel(dmA, name, &label);
6002: DMLabelDuplicate(label, &labelNew);
6003: DMAddLabel(dmB, labelNew);
6004: }
6005: return(0);
6006: }
6008: /*@
6009: DMGetCoarseDM - Get the coarse mesh from which this was obtained by refinement
6011: Input Parameter:
6012: . dm - The DM object
6014: Output Parameter:
6015: . cdm - The coarse DM
6017: Level: intermediate
6019: .seealso: DMSetCoarseDM()
6020: @*/
6021: PetscErrorCode DMGetCoarseDM(DM dm, DM *cdm)
6022: {
6026: *cdm = dm->coarseMesh;
6027: return(0);
6028: }
6030: /*@
6031: DMSetCoarseDM - Set the coarse mesh from which this was obtained by refinement
6033: Input Parameters:
6034: + dm - The DM object
6035: - cdm - The coarse DM
6037: Level: intermediate
6039: .seealso: DMGetCoarseDM()
6040: @*/
6041: PetscErrorCode DMSetCoarseDM(DM dm, DM cdm)
6042: {
6048: PetscObjectReference((PetscObject)cdm);
6049: DMDestroy(&dm->coarseMesh);
6050: dm->coarseMesh = cdm;
6051: return(0);
6052: }
6054: /*@
6055: DMGetFineDM - Get the fine mesh from which this was obtained by refinement
6057: Input Parameter:
6058: . dm - The DM object
6060: Output Parameter:
6061: . fdm - The fine DM
6063: Level: intermediate
6065: .seealso: DMSetFineDM()
6066: @*/
6067: PetscErrorCode DMGetFineDM(DM dm, DM *fdm)
6068: {
6072: *fdm = dm->fineMesh;
6073: return(0);
6074: }
6076: /*@
6077: DMSetFineDM - Set the fine mesh from which this was obtained by refinement
6079: Input Parameters:
6080: + dm - The DM object
6081: - fdm - The fine DM
6083: Level: intermediate
6085: .seealso: DMGetFineDM()
6086: @*/
6087: PetscErrorCode DMSetFineDM(DM dm, DM fdm)
6088: {
6094: PetscObjectReference((PetscObject)fdm);
6095: DMDestroy(&dm->fineMesh);
6096: dm->fineMesh = fdm;
6097: return(0);
6098: }
6100: /*=== DMBoundary code ===*/
6102: PetscErrorCode DMCopyBoundary(DM dm, DM dmNew)
6103: {
6107: PetscDSCopyBoundary(dm->prob,dmNew->prob);
6108: return(0);
6109: }
6111: /*@C
6112: DMAddBoundary - Add a boundary condition to the model
6114: Input Parameters:
6115: + dm - The DM, with a PetscDS that matches the problem being constrained
6116: . type - The type of condition, e.g. DM_BC_ESSENTIAL_ANALYTIC/DM_BC_ESSENTIAL_FIELD (Dirichlet), or DM_BC_NATURAL (Neumann)
6117: . name - The BC name
6118: . labelname - The label defining constrained points
6119: . field - The field to constrain
6120: . numcomps - The number of constrained field components
6121: . comps - An array of constrained component numbers
6122: . bcFunc - A pointwise function giving boundary values
6123: . numids - The number of DMLabel ids for constrained points
6124: . ids - An array of ids for constrained points
6125: - ctx - An optional user context for bcFunc
6127: Options Database Keys:
6128: + -bc_<boundary name> <num> - Overrides the boundary ids
6129: - -bc_<boundary name>_comp <num> - Overrides the boundary components
6131: Level: developer
6133: .seealso: DMGetBoundary()
6134: @*/
6135: PetscErrorCode DMAddBoundary(DM dm, DMBoundaryConditionType type, const char name[], const char labelname[], PetscInt field, PetscInt numcomps, const PetscInt *comps, void (*bcFunc)(void), PetscInt numids, const PetscInt *ids, void *ctx)
6136: {
6141: PetscDSAddBoundary(dm->prob,type,name,labelname,field,numcomps,comps,bcFunc,numids,ids,ctx);
6142: return(0);
6143: }
6145: /*@
6146: DMGetNumBoundary - Get the number of registered BC
6148: Input Parameters:
6149: . dm - The mesh object
6151: Output Parameters:
6152: . numBd - The number of BC
6154: Level: intermediate
6156: .seealso: DMAddBoundary(), DMGetBoundary()
6157: @*/
6158: PetscErrorCode DMGetNumBoundary(DM dm, PetscInt *numBd)
6159: {
6164: PetscDSGetNumBoundary(dm->prob,numBd);
6165: return(0);
6166: }
6168: /*@C
6169: DMGetBoundary - Get a model boundary condition
6171: Input Parameters:
6172: + dm - The mesh object
6173: - bd - The BC number
6175: Output Parameters:
6176: + type - The type of condition, e.g. DM_BC_ESSENTIAL_ANALYTIC/DM_BC_ESSENTIAL_FIELD (Dirichlet), or DM_BC_NATURAL (Neumann)
6177: . name - The BC name
6178: . labelname - The label defining constrained points
6179: . field - The field to constrain
6180: . numcomps - The number of constrained field components
6181: . comps - An array of constrained component numbers
6182: . bcFunc - A pointwise function giving boundary values
6183: . numids - The number of DMLabel ids for constrained points
6184: . ids - An array of ids for constrained points
6185: - ctx - An optional user context for bcFunc
6187: Options Database Keys:
6188: + -bc_<boundary name> <num> - Overrides the boundary ids
6189: - -bc_<boundary name>_comp <num> - Overrides the boundary components
6191: Level: developer
6193: .seealso: DMAddBoundary()
6194: @*/
6195: PetscErrorCode DMGetBoundary(DM dm, PetscInt bd, DMBoundaryConditionType *type, const char **name, const char **labelname, PetscInt *field, PetscInt *numcomps, const PetscInt **comps, void (**func)(void), PetscInt *numids, const PetscInt **ids, void **ctx)
6196: {
6201: PetscDSGetBoundary(dm->prob,bd,type,name,labelname,field,numcomps,comps,func,numids,ids,ctx);
6202: return(0);
6203: }
6205: static PetscErrorCode DMPopulateBoundary(DM dm)
6206: {
6207: DMBoundary *lastnext;
6208: DSBoundary dsbound;
6212: dsbound = dm->prob->boundary;
6213: if (dm->boundary) {
6214: DMBoundary next = dm->boundary;
6216: /* quick check to see if the PetscDS has changed */
6217: if (next->dsboundary == dsbound) return(0);
6218: /* the PetscDS has changed: tear down and rebuild */
6219: while (next) {
6220: DMBoundary b = next;
6222: next = b->next;
6223: PetscFree(b);
6224: }
6225: dm->boundary = NULL;
6226: }
6228: lastnext = &(dm->boundary);
6229: while (dsbound) {
6230: DMBoundary dmbound;
6232: PetscNew(&dmbound);
6233: dmbound->dsboundary = dsbound;
6234: DMGetLabel(dm, dsbound->labelname, &(dmbound->label));
6235: if (!dmbound->label) {PetscInfo2(dm, "DSBoundary %s wants label %s, which is not in this dm.\n",dsbound->name,dsbound->labelname);}
6236: /* push on the back instead of the front so that it is in the same order as in the PetscDS */
6237: *lastnext = dmbound;
6238: lastnext = &(dmbound->next);
6239: dsbound = dsbound->next;
6240: }
6241: return(0);
6242: }
6244: PetscErrorCode DMIsBoundaryPoint(DM dm, PetscInt point, PetscBool *isBd)
6245: {
6246: DMBoundary b;
6252: *isBd = PETSC_FALSE;
6253: DMPopulateBoundary(dm);
6254: b = dm->boundary;
6255: while (b && !(*isBd)) {
6256: DMLabel label = b->label;
6257: DSBoundary dsb = b->dsboundary;
6259: if (label) {
6260: PetscInt i;
6262: for (i = 0; i < dsb->numids && !(*isBd); ++i) {
6263: DMLabelStratumHasPoint(label, dsb->ids[i], point, isBd);
6264: }
6265: }
6266: b = b->next;
6267: }
6268: return(0);
6269: }
6271: /*@C
6272: DMProjectFunction - This projects the given function into the function space provided.
6274: Input Parameters:
6275: + dm - The DM
6276: . time - The time
6277: . funcs - The coordinate functions to evaluate, one per field
6278: . ctxs - Optional array of contexts to pass to each coordinate function. ctxs itself may be null.
6279: - mode - The insertion mode for values
6281: Output Parameter:
6282: . X - vector
6284: Calling sequence of func:
6285: $ func(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar u[], void *ctx);
6287: + dim - The spatial dimension
6288: . x - The coordinates
6289: . Nf - The number of fields
6290: . u - The output field values
6291: - ctx - optional user-defined function context
6293: Level: developer
6295: .seealso: DMComputeL2Diff()
6296: @*/
6297: PetscErrorCode DMProjectFunction(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec X)
6298: {
6299: Vec localX;
6304: DMGetLocalVector(dm, &localX);
6305: DMProjectFunctionLocal(dm, time, funcs, ctxs, mode, localX);
6306: DMLocalToGlobalBegin(dm, localX, mode, X);
6307: DMLocalToGlobalEnd(dm, localX, mode, X);
6308: DMRestoreLocalVector(dm, &localX);
6309: return(0);
6310: }
6312: PetscErrorCode DMProjectFunctionLocal(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX)
6313: {
6319: if (!dm->ops->projectfunctionlocal) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMProjectFunctionLocal",((PetscObject)dm)->type_name);
6320: (dm->ops->projectfunctionlocal) (dm, time, funcs, ctxs, mode, localX);
6321: return(0);
6322: }
6324: PetscErrorCode DMProjectFunctionLabel(DM dm, PetscReal time, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscInt Nc, const PetscInt comps[], PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec X)
6325: {
6326: Vec localX;
6331: DMGetLocalVector(dm, &localX);
6332: DMProjectFunctionLabelLocal(dm, time, label, numIds, ids, Nc, comps, funcs, ctxs, mode, localX);
6333: DMLocalToGlobalBegin(dm, localX, mode, X);
6334: DMLocalToGlobalEnd(dm, localX, mode, X);
6335: DMRestoreLocalVector(dm, &localX);
6336: return(0);
6337: }
6339: PetscErrorCode DMProjectFunctionLabelLocal(DM dm, PetscReal time, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscInt Nc, const PetscInt comps[], PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX)
6340: {
6346: if (!dm->ops->projectfunctionlabellocal) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMProjectFunctionLabelLocal",((PetscObject)dm)->type_name);
6347: (dm->ops->projectfunctionlabellocal) (dm, time, label, numIds, ids, Nc, comps, funcs, ctxs, mode, localX);
6348: return(0);
6349: }
6351: PetscErrorCode DMProjectFieldLocal(DM dm, PetscReal time, Vec localU,
6352: void (**funcs)(PetscInt, PetscInt, PetscInt,
6353: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
6354: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
6355: PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]),
6356: InsertMode mode, Vec localX)
6357: {
6364: if (!dm->ops->projectfieldlocal) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMProjectFieldLocal",((PetscObject)dm)->type_name);
6365: (dm->ops->projectfieldlocal) (dm, time, localU, funcs, mode, localX);
6366: return(0);
6367: }
6369: PetscErrorCode DMProjectFieldLabelLocal(DM dm, PetscReal time, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscInt Nc, const PetscInt comps[], Vec localU,
6370: void (**funcs)(PetscInt, PetscInt, PetscInt,
6371: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
6372: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
6373: PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]),
6374: InsertMode mode, Vec localX)
6375: {
6382: if (!dm->ops->projectfieldlabellocal) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMProjectFieldLocal",((PetscObject)dm)->type_name);
6383: (dm->ops->projectfieldlabellocal)(dm, time, label, numIds, ids, Nc, comps, localU, funcs, mode, localX);
6384: return(0);
6385: }
6387: /*@C
6388: DMComputeL2Diff - This function computes the L_2 difference between a function u and an FEM interpolant solution u_h.
6390: Input Parameters:
6391: + dm - The DM
6392: . time - The time
6393: . funcs - The functions to evaluate for each field component
6394: . ctxs - Optional array of contexts to pass to each function, or NULL.
6395: - X - The coefficient vector u_h, a global vector
6397: Output Parameter:
6398: . diff - The diff ||u - u_h||_2
6400: Level: developer
6402: .seealso: DMProjectFunction(), DMComputeL2FieldDiff(), DMComputeL2GradientDiff()
6403: @*/
6404: PetscErrorCode DMComputeL2Diff(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff)
6405: {
6411: if (!dm->ops->computel2diff) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMComputeL2Diff",((PetscObject)dm)->type_name);
6412: (dm->ops->computel2diff)(dm,time,funcs,ctxs,X,diff);
6413: return(0);
6414: }
6416: /*@C
6417: DMComputeL2GradientDiff - This function computes the L_2 difference between the gradient of a function u and an FEM interpolant solution grad u_h.
6419: Input Parameters:
6420: + dm - The DM
6421: , time - The time
6422: . funcs - The gradient functions to evaluate for each field component
6423: . ctxs - Optional array of contexts to pass to each function, or NULL.
6424: . X - The coefficient vector u_h, a global vector
6425: - n - The vector to project along
6427: Output Parameter:
6428: . diff - The diff ||(grad u - grad u_h) . n||_2
6430: Level: developer
6432: .seealso: DMProjectFunction(), DMComputeL2Diff()
6433: @*/
6434: PetscErrorCode DMComputeL2GradientDiff(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], const PetscReal[], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, const PetscReal n[], PetscReal *diff)
6435: {
6441: if (!dm->ops->computel2gradientdiff) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMComputeL2GradientDiff",((PetscObject)dm)->type_name);
6442: (dm->ops->computel2gradientdiff)(dm,time,funcs,ctxs,X,n,diff);
6443: return(0);
6444: }
6446: /*@C
6447: DMComputeL2FieldDiff - This function computes the L_2 difference between a function u and an FEM interpolant solution u_h, separated into field components.
6449: Input Parameters:
6450: + dm - The DM
6451: . time - The time
6452: . funcs - The functions to evaluate for each field component
6453: . ctxs - Optional array of contexts to pass to each function, or NULL.
6454: - X - The coefficient vector u_h, a global vector
6456: Output Parameter:
6457: . diff - The array of differences, ||u^f - u^f_h||_2
6459: Level: developer
6461: .seealso: DMProjectFunction(), DMComputeL2FieldDiff(), DMComputeL2GradientDiff()
6462: @*/
6463: PetscErrorCode DMComputeL2FieldDiff(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal diff[])
6464: {
6470: if (!dm->ops->computel2fielddiff) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMComputeL2FieldDiff",((PetscObject)dm)->type_name);
6471: (dm->ops->computel2fielddiff)(dm,time,funcs,ctxs,X,diff);
6472: return(0);
6473: }
6475: /*@C
6476: DMAdaptLabel - Adapt a dm based on a label with values interpreted as coarsening and refining flags. Specific implementations of DM maybe have
6477: specialized flags, but all implementations should accept flag values DM_ADAPT_DETERMINE, DM_ADAPT_KEEP, DM_ADAPT_REFINE, and DM_ADAPT_COARSEN.
6479: Collective on dm
6481: Input parameters:
6482: + dm - the pre-adaptation DM object
6483: - label - label with the flags
6485: Output parameters:
6486: . dmAdapt - the adapted DM object: may be NULL if an adapted DM could not be produced.
6488: Level: intermediate
6490: .seealso: DMAdaptMetric(), DMCoarsen(), DMRefine()
6491: @*/
6492: PetscErrorCode DMAdaptLabel(DM dm, DMLabel label, DM *dmAdapt)
6493: {
6500: *dmAdapt = NULL;
6501: if (!dm->ops->adaptlabel) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMAdaptLabel",((PetscObject)dm)->type_name);
6502: (dm->ops->adaptlabel)(dm, label, dmAdapt);
6503: return(0);
6504: }
6506: /*@C
6507: DMAdaptMetric - Generates a mesh adapted to the specified metric field using the pragmatic library.
6509: Input Parameters:
6510: + dm - The DM object
6511: . metric - The metric to which the mesh is adapted, defined vertex-wise.
6512: - bdLabel - Label for boundary tags, which will be preserved in the output mesh. bdLabel should be NULL if there is no such label, and should be different from "_boundary_".
6514: Output Parameter:
6515: . dmAdapt - Pointer to the DM object containing the adapted mesh
6517: Note: The label in the adapted mesh will be registered under the name of the input DMLabel object
6519: Level: advanced
6521: .seealso: DMAdaptLabel(), DMCoarsen(), DMRefine()
6522: @*/
6523: PetscErrorCode DMAdaptMetric(DM dm, Vec metric, DMLabel bdLabel, DM *dmAdapt)
6524: {
6532: *dmAdapt = NULL;
6533: if (!dm->ops->adaptmetric) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMAdaptMetric",((PetscObject)dm)->type_name);
6534: (dm->ops->adaptmetric)(dm, metric, bdLabel, dmAdapt);
6535: return(0);
6536: }
6538: /*@C
6539: DMGetNeighbors - Gets an array containing the MPI rank of all the processes neighbors
6541: Not Collective
6543: Input Parameter:
6544: . dm - The DM
6546: Output Parameter:
6547: . nranks - the number of neighbours
6548: . ranks - the neighbors ranks
6550: Notes:
6551: Do not free the array, it is freed when the DM is destroyed.
6553: Level: beginner
6555: .seealso: DMDAGetNeighbors(), PetscSFGetRanks()
6556: @*/
6557: PetscErrorCode DMGetNeighbors(DM dm,PetscInt *nranks,const PetscMPIInt *ranks[])
6558: {
6563: if (!dm->ops->getneighbors) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMGetNeighbors",((PetscObject)dm)->type_name);
6564: (dm->ops->getneighbors)(dm,nranks,ranks);
6565: return(0);
6566: }
6568: #include <petsc/private/matimpl.h> /* Needed because of coloring->ctype below */
6570: /*
6571: Converts the input vector to a ghosted vector and then calls the standard coloring code.
6572: This has be a different function because it requires DM which is not defined in the Mat library
6573: */
6574: PetscErrorCode MatFDColoringApply_AIJDM(Mat J,MatFDColoring coloring,Vec x1,void *sctx)
6575: {
6579: if (coloring->ctype == IS_COLORING_LOCAL) {
6580: Vec x1local;
6581: DM dm;
6582: MatGetDM(J,&dm);
6583: if (!dm) SETERRQ(PetscObjectComm((PetscObject)J),PETSC_ERR_ARG_INCOMP,"IS_COLORING_LOCAL requires a DM");
6584: DMGetLocalVector(dm,&x1local);
6585: DMGlobalToLocalBegin(dm,x1,INSERT_VALUES,x1local);
6586: DMGlobalToLocalEnd(dm,x1,INSERT_VALUES,x1local);
6587: x1 = x1local;
6588: }
6589: MatFDColoringApply_AIJ(J,coloring,x1,sctx);
6590: if (coloring->ctype == IS_COLORING_LOCAL) {
6591: DM dm;
6592: MatGetDM(J,&dm);
6593: DMRestoreLocalVector(dm,&x1);
6594: }
6595: return(0);
6596: }
6598: /*@
6599: MatFDColoringUseDM - allows a MatFDColoring object to use the DM associated with the matrix to use a IS_COLORING_LOCAL coloring
6601: Input Parameter:
6602: . coloring - the MatFDColoring object
6604: Developer Notes:
6605: this routine exists because the PETSc Mat library does not know about the DM objects
6607: Level: advanced
6609: .seealso: MatFDColoring, MatFDColoringCreate(), ISColoringType
6610: @*/
6611: PetscErrorCode MatFDColoringUseDM(Mat coloring,MatFDColoring fdcoloring)
6612: {
6614: coloring->ops->fdcoloringapply = MatFDColoringApply_AIJDM;
6615: return(0);
6616: }
6618: /*@
6619: DMGetCompatibility - determine if two DMs are compatible
6621: Collective
6623: Input Parameters:
6624: + dm - the first DM
6625: - dm2 - the second DM
6627: Output Parameters:
6628: + compatible - whether or not the two DMs are compatible
6629: - set - whether or not the compatible value was set
6631: Notes:
6632: Two DMs are deemed compatible if they represent the same parallel decomposition
6633: of the same topology. This implies that the the section (field data) on one
6634: "makes sense" with respect to the topology and parallel decomposition of the other.
6635: Loosely speaking, compatibile DMs represent the same domain, with the same parallel
6636: decomposition, with different data.
6638: Typically, one would confirm compatibility if intending to simultaneously iterate
6639: over a pair of vectors obtained from different DMs.
6641: For example, two DMDA objects are compatible if they have the same local
6642: and global sizes and the same stencil width. They can have different numbers
6643: of degrees of freedom per node. Thus, one could use the node numbering from
6644: either DM in bounds for a loop over vectors derived from either DM.
6646: Consider the operation of summing data living on a 2-dof DMDA to data living
6647: on a 1-dof DMDA, which should be compatible, as in the following snippet.
6648: .vb
6649: ...
6650: DMGetCompatibility(da1,da2,&compatible,&set);
6651: if (set && compatible) {
6652: DMDAVecGetArrayDOF(da1,vec1,&arr1);
6653: DMDAVecGetArrayDOF(da2,vec2,&arr2);
6654: DMDAGetCorners(da1,&x,&y,NULL,&m,&n);
6655: for (j=y; j<y+n; ++j) {
6656: for (i=x; i<x+m, ++i) {
6657: arr1[j][i][0] = arr2[j][i][0] + arr2[j][i][1];
6658: }
6659: }
6660: DMDAVecRestoreArrayDOF(da1,vec1,&arr1);
6661: DMDAVecRestoreArrayDOF(da2,vec2,&arr2);
6662: } else {
6663: SETERRQ(PetscObjectComm((PetscObject)da1,PETSC_ERR_ARG_INCOMP,"DMDA objects incompatible");
6664: }
6665: ...
6666: .ve
6668: Checking compatibility might be expensive for a given implementation of DM,
6669: or might be impossible to unambiguously confirm or deny. For this reason,
6670: this function may decline to determine compatibility, and hence users should
6671: always check the "set" output parameter.
6673: A DM is always compatible with itself.
6675: In the current implementation, DMs which live on "unequal" communicators
6676: (MPI_UNEQUAL in the terminology of MPI_Comm_compare()) are always deemed
6677: incompatible.
6679: This function is labeled "Collective," as information about all subdomains
6680: is required on each rank. However, in DM implementations which store all this
6681: information locally, this function may be merely "Logically Collective".
6683: Developer Notes:
6684: Compatibility is assumed to be a symmetric concept; if DM A is compatible with DM B,
6685: the DM B is compatible with DM A. Thus, this function checks the implementations
6686: of both dm and dm2 (if they are of different types), attempting to determine
6687: compatibility. It is left to DM implementers to ensure that symmetry is
6688: preserved. The simplest way to do this is, when implementing type-specific
6689: logic for this function, to check for existing logic in the implementation
6690: of other DM types and let *set = PETSC_FALSE if found; the logic of this
6691: function will then call that logic.
6693: Level: advanced
6695: .seealso: DM, DMDACreateCompatibleDMDA()
6696: @*/
6698: PetscErrorCode DMGetCompatibility(DM dm,DM dm2,PetscBool *compatible,PetscBool *set)
6699: {
6701: PetscMPIInt compareResult;
6702: DMType type,type2;
6703: PetscBool sameType;
6709: /* Declare a DM compatible with itself */
6710: if (dm == dm2) {
6711: *set = PETSC_TRUE;
6712: *compatible = PETSC_TRUE;
6713: return(0);
6714: }
6716: /* Declare a DM incompatible with a DM that lives on an "unequal"
6717: communicator. Note that this does not preclude compatibility with
6718: DMs living on "congruent" or "similar" communicators, but this must be
6719: determined by the implementation-specific logic */
6720: MPI_Comm_compare(PetscObjectComm((PetscObject)dm),PetscObjectComm((PetscObject)dm2),&compareResult);
6721: if (compareResult == MPI_UNEQUAL) {
6722: *set = PETSC_TRUE;
6723: *compatible = PETSC_FALSE;
6724: return(0);
6725: }
6727: /* Pass to the implementation-specific routine, if one exists. */
6728: if (dm->ops->getcompatibility) {
6729: (*dm->ops->getcompatibility)(dm,dm2,compatible,set);
6730: if (*set) {
6731: return(0);
6732: }
6733: }
6735: /* If dm and dm2 are of different types, then attempt to check compatibility
6736: with an implementation of this function from dm2 */
6737: DMGetType(dm,&type);
6738: DMGetType(dm2,&type2);
6739: PetscStrcmp(type,type2,&sameType);
6740: if (!sameType && dm2->ops->getcompatibility) {
6741: (*dm2->ops->getcompatibility)(dm2,dm,compatible,set); /* Note argument order */
6742: } else {
6743: *set = PETSC_FALSE;
6744: }
6745: return(0);
6746: }