Actual source code: pack.c
2: #include <../src/dm/impls/composite/packimpl.h>
3: #include <petsc/private/isimpl.h>
4: #include <petsc/private/glvisviewerimpl.h>
5: #include <petscds.h>
7: /*@C
8: DMCompositeSetCoupling - Sets user provided routines that compute the coupling between the
9: separate components (DMs) in a DMto build the correct matrix nonzero structure.
11: Logically Collective
13: Input Parameters:
14: + dm - the composite object
15: - formcouplelocations - routine to set the nonzero locations in the matrix
17: Level: advanced
19: Not available from Fortran
21: Notes:
22: See DMSetApplicationContext() and DMGetApplicationContext() for how to get user information into
23: this routine
25: @*/
26: PetscErrorCode DMCompositeSetCoupling(DM dm, PetscErrorCode (*FormCoupleLocations)(DM, Mat, PetscInt *, PetscInt *, PetscInt, PetscInt, PetscInt, PetscInt))
27: {
28: DM_Composite *com = (DM_Composite *)dm->data;
29: PetscBool flg;
31: PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg);
33: com->FormCoupleLocations = FormCoupleLocations;
34: return 0;
35: }
37: PetscErrorCode DMDestroy_Composite(DM dm)
38: {
39: struct DMCompositeLink *next, *prev;
40: DM_Composite *com = (DM_Composite *)dm->data;
42: next = com->next;
43: while (next) {
44: prev = next;
45: next = next->next;
46: DMDestroy(&prev->dm);
47: PetscFree(prev->grstarts);
48: PetscFree(prev);
49: }
50: PetscObjectComposeFunction((PetscObject)dm, "DMSetUpGLVisViewer_C", NULL);
51: /* This was originally freed in DMDestroy(), but that prevents reference counting of backend objects */
52: PetscFree(com);
53: return 0;
54: }
56: PetscErrorCode DMView_Composite(DM dm, PetscViewer v)
57: {
58: PetscBool iascii;
59: DM_Composite *com = (DM_Composite *)dm->data;
61: PetscObjectTypeCompare((PetscObject)v, PETSCVIEWERASCII, &iascii);
62: if (iascii) {
63: struct DMCompositeLink *lnk = com->next;
64: PetscInt i;
66: PetscViewerASCIIPrintf(v, "DM (%s)\n", ((PetscObject)dm)->prefix ? ((PetscObject)dm)->prefix : "no prefix");
67: PetscViewerASCIIPrintf(v, " contains %" PetscInt_FMT " DMs\n", com->nDM);
68: PetscViewerASCIIPushTab(v);
69: for (i = 0; lnk; lnk = lnk->next, i++) {
70: PetscViewerASCIIPrintf(v, "Link %" PetscInt_FMT ": DM of type %s\n", i, ((PetscObject)lnk->dm)->type_name);
71: PetscViewerASCIIPushTab(v);
72: DMView(lnk->dm, v);
73: PetscViewerASCIIPopTab(v);
74: }
75: PetscViewerASCIIPopTab(v);
76: }
77: return 0;
78: }
80: /* --------------------------------------------------------------------------------------*/
81: PetscErrorCode DMSetUp_Composite(DM dm)
82: {
83: PetscInt nprev = 0;
84: PetscMPIInt rank, size;
85: DM_Composite *com = (DM_Composite *)dm->data;
86: struct DMCompositeLink *next = com->next;
87: PetscLayout map;
90: PetscLayoutCreate(PetscObjectComm((PetscObject)dm), &map);
91: PetscLayoutSetLocalSize(map, com->n);
92: PetscLayoutSetSize(map, PETSC_DETERMINE);
93: PetscLayoutSetBlockSize(map, 1);
94: PetscLayoutSetUp(map);
95: PetscLayoutGetSize(map, &com->N);
96: PetscLayoutGetRange(map, &com->rstart, NULL);
97: PetscLayoutDestroy(&map);
99: /* now set the rstart for each linked vector */
100: MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
101: MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size);
102: while (next) {
103: next->rstart = nprev;
104: nprev += next->n;
105: next->grstart = com->rstart + next->rstart;
106: PetscMalloc1(size, &next->grstarts);
107: MPI_Allgather(&next->grstart, 1, MPIU_INT, next->grstarts, 1, MPIU_INT, PetscObjectComm((PetscObject)dm));
108: next = next->next;
109: }
110: com->setup = PETSC_TRUE;
111: return 0;
112: }
114: /* ----------------------------------------------------------------------------------*/
116: /*@
117: DMCompositeGetNumberDM - Get's the number of DM objects in the DMComposite
118: representation.
120: Not Collective
122: Input Parameter:
123: . dm - the packer object
125: Output Parameter:
126: . nDM - the number of DMs
128: Level: beginner
130: @*/
131: PetscErrorCode DMCompositeGetNumberDM(DM dm, PetscInt *nDM)
132: {
133: DM_Composite *com = (DM_Composite *)dm->data;
134: PetscBool flg;
137: PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg);
139: *nDM = com->nDM;
140: return 0;
141: }
143: /*@C
144: DMCompositeGetAccess - Allows one to access the individual packed vectors in their global
145: representation.
147: Collective on dm
149: Input Parameters:
150: + dm - the packer object
151: - gvec - the global vector
153: Output Parameters:
154: . Vec* ... - the packed parallel vectors, NULL for those that are not needed
156: Notes:
157: Use DMCompositeRestoreAccess() to return the vectors when you no longer need them
159: Fortran Notes:
161: Fortran callers must use numbered versions of this routine, e.g., DMCompositeGetAccess4(dm,gvec,vec1,vec2,vec3,vec4)
162: or use the alternative interface DMCompositeGetAccessArray().
164: Level: advanced
166: .seealso: `DMCompositeGetEntries()`, `DMCompositeScatter()`
167: @*/
168: PetscErrorCode DMCompositeGetAccess(DM dm, Vec gvec, ...)
169: {
170: va_list Argp;
171: struct DMCompositeLink *next;
172: DM_Composite *com = (DM_Composite *)dm->data;
173: PetscInt readonly;
174: PetscBool flg;
178: PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg);
180: next = com->next;
181: if (!com->setup) DMSetUp(dm);
183: VecLockGet(gvec, &readonly);
184: /* loop over packed objects, handling one at at time */
185: va_start(Argp, gvec);
186: while (next) {
187: Vec *vec;
188: vec = va_arg(Argp, Vec *);
189: if (vec) {
190: DMGetGlobalVector(next->dm, vec);
191: if (readonly) {
192: const PetscScalar *array;
193: VecGetArrayRead(gvec, &array);
194: VecPlaceArray(*vec, array + next->rstart);
195: VecLockReadPush(*vec);
196: VecRestoreArrayRead(gvec, &array);
197: } else {
198: PetscScalar *array;
199: VecGetArray(gvec, &array);
200: VecPlaceArray(*vec, array + next->rstart);
201: VecRestoreArray(gvec, &array);
202: }
203: }
204: next = next->next;
205: }
206: va_end(Argp);
207: return 0;
208: }
210: /*@C
211: DMCompositeGetAccessArray - Allows one to access the individual packed vectors in their global
212: representation.
214: Collective on dm
216: Input Parameters:
217: + dm - the packer object
218: . pvec - packed vector
219: . nwanted - number of vectors wanted
220: - wanted - sorted array of vectors wanted, or NULL to get all vectors
222: Output Parameters:
223: . vecs - array of requested global vectors (must be allocated)
225: Notes:
226: Use DMCompositeRestoreAccessArray() to return the vectors when you no longer need them
228: Level: advanced
230: .seealso: `DMCompositeGetAccess()`, `DMCompositeGetEntries()`, `DMCompositeScatter()`, `DMCompositeGather()`
231: @*/
232: PetscErrorCode DMCompositeGetAccessArray(DM dm, Vec pvec, PetscInt nwanted, const PetscInt *wanted, Vec *vecs)
233: {
234: struct DMCompositeLink *link;
235: PetscInt i, wnum;
236: DM_Composite *com = (DM_Composite *)dm->data;
237: PetscInt readonly;
238: PetscBool flg;
242: PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg);
244: if (!com->setup) DMSetUp(dm);
246: VecLockGet(pvec, &readonly);
247: for (i = 0, wnum = 0, link = com->next; link && wnum < nwanted; i++, link = link->next) {
248: if (!wanted || i == wanted[wnum]) {
249: Vec v;
250: DMGetGlobalVector(link->dm, &v);
251: if (readonly) {
252: const PetscScalar *array;
253: VecGetArrayRead(pvec, &array);
254: VecPlaceArray(v, array + link->rstart);
255: VecLockReadPush(v);
256: VecRestoreArrayRead(pvec, &array);
257: } else {
258: PetscScalar *array;
259: VecGetArray(pvec, &array);
260: VecPlaceArray(v, array + link->rstart);
261: VecRestoreArray(pvec, &array);
262: }
263: vecs[wnum++] = v;
264: }
265: }
266: return 0;
267: }
269: /*@C
270: DMCompositeGetLocalAccessArray - Allows one to access the individual
271: packed vectors in their local representation.
273: Collective on dm.
275: Input Parameters:
276: + dm - the packer object
277: . pvec - packed vector
278: . nwanted - number of vectors wanted
279: - wanted - sorted array of vectors wanted, or NULL to get all vectors
281: Output Parameters:
282: . vecs - array of requested local vectors (must be allocated)
284: Notes:
285: Use DMCompositeRestoreLocalAccessArray() to return the vectors
286: when you no longer need them.
288: Level: advanced
290: .seealso: `DMCompositeRestoreLocalAccessArray()`, `DMCompositeGetAccess()`,
291: `DMCompositeGetEntries()`, `DMCompositeScatter()`, `DMCompositeGather()`
292: @*/
293: PetscErrorCode DMCompositeGetLocalAccessArray(DM dm, Vec pvec, PetscInt nwanted, const PetscInt *wanted, Vec *vecs)
294: {
295: struct DMCompositeLink *link;
296: PetscInt i, wnum;
297: DM_Composite *com = (DM_Composite *)dm->data;
298: PetscInt readonly;
299: PetscInt nlocal = 0;
300: PetscBool flg;
304: PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg);
306: if (!com->setup) DMSetUp(dm);
308: VecLockGet(pvec, &readonly);
309: for (i = 0, wnum = 0, link = com->next; link && wnum < nwanted; i++, link = link->next) {
310: if (!wanted || i == wanted[wnum]) {
311: Vec v;
312: DMGetLocalVector(link->dm, &v);
313: if (readonly) {
314: const PetscScalar *array;
315: VecGetArrayRead(pvec, &array);
316: VecPlaceArray(v, array + nlocal);
317: // this method does not make sense. The local vectors are not updated with a global-to-local and the user can not do it because it is locked
318: VecLockReadPush(v);
319: VecRestoreArrayRead(pvec, &array);
320: } else {
321: PetscScalar *array;
322: VecGetArray(pvec, &array);
323: VecPlaceArray(v, array + nlocal);
324: VecRestoreArray(pvec, &array);
325: }
326: vecs[wnum++] = v;
327: }
329: nlocal += link->nlocal;
330: }
332: return 0;
333: }
335: /*@C
336: DMCompositeRestoreAccess - Returns the vectors obtained with DMCompositeGetAccess()
337: representation.
339: Collective on dm
341: Input Parameters:
342: + dm - the packer object
343: . gvec - the global vector
344: - Vec* ... - the individual parallel vectors, NULL for those that are not needed
346: Level: advanced
348: .seealso `DMCompositeAddDM()`, `DMCreateGlobalVector()`,
349: `DMCompositeGather()`, `DMCompositeCreate()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeScatter()`,
350: `DMCompositeRestoreAccess()`, `DMCompositeGetAccess()`
352: @*/
353: PetscErrorCode DMCompositeRestoreAccess(DM dm, Vec gvec, ...)
354: {
355: va_list Argp;
356: struct DMCompositeLink *next;
357: DM_Composite *com = (DM_Composite *)dm->data;
358: PetscInt readonly;
359: PetscBool flg;
363: PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg);
365: next = com->next;
366: if (!com->setup) DMSetUp(dm);
368: VecLockGet(gvec, &readonly);
369: /* loop over packed objects, handling one at at time */
370: va_start(Argp, gvec);
371: while (next) {
372: Vec *vec;
373: vec = va_arg(Argp, Vec *);
374: if (vec) {
375: VecResetArray(*vec);
376: if (readonly) VecLockReadPop(*vec);
377: DMRestoreGlobalVector(next->dm, vec);
378: }
379: next = next->next;
380: }
381: va_end(Argp);
382: return 0;
383: }
385: /*@C
386: DMCompositeRestoreAccessArray - Returns the vectors obtained with DMCompositeGetAccessArray()
388: Collective on dm
390: Input Parameters:
391: + dm - the packer object
392: . pvec - packed vector
393: . nwanted - number of vectors wanted
394: . wanted - sorted array of vectors wanted, or NULL to get all vectors
395: - vecs - array of global vectors to return
397: Level: advanced
399: .seealso: `DMCompositeRestoreAccess()`, `DMCompositeRestoreEntries()`, `DMCompositeScatter()`, `DMCompositeGather()`
400: @*/
401: PetscErrorCode DMCompositeRestoreAccessArray(DM dm, Vec pvec, PetscInt nwanted, const PetscInt *wanted, Vec *vecs)
402: {
403: struct DMCompositeLink *link;
404: PetscInt i, wnum;
405: DM_Composite *com = (DM_Composite *)dm->data;
406: PetscInt readonly;
407: PetscBool flg;
411: PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg);
413: if (!com->setup) DMSetUp(dm);
415: VecLockGet(pvec, &readonly);
416: for (i = 0, wnum = 0, link = com->next; link && wnum < nwanted; i++, link = link->next) {
417: if (!wanted || i == wanted[wnum]) {
418: VecResetArray(vecs[wnum]);
419: if (readonly) VecLockReadPop(vecs[wnum]);
420: DMRestoreGlobalVector(link->dm, &vecs[wnum]);
421: wnum++;
422: }
423: }
424: return 0;
425: }
427: /*@C
428: DMCompositeRestoreLocalAccessArray - Returns the vectors obtained with DMCompositeGetLocalAccessArray().
430: Collective on dm.
432: Input Parameters:
433: + dm - the packer object
434: . pvec - packed vector
435: . nwanted - number of vectors wanted
436: . wanted - sorted array of vectors wanted, or NULL to restore all vectors
437: - vecs - array of local vectors to return
439: Level: advanced
441: Notes:
442: nwanted and wanted must match the values given to DMCompositeGetLocalAccessArray()
443: otherwise the call will fail.
445: .seealso: `DMCompositeGetLocalAccessArray()`, `DMCompositeRestoreAccessArray()`,
446: `DMCompositeRestoreAccess()`, `DMCompositeRestoreEntries()`,
447: `DMCompositeScatter()`, `DMCompositeGather()`
448: @*/
449: PetscErrorCode DMCompositeRestoreLocalAccessArray(DM dm, Vec pvec, PetscInt nwanted, const PetscInt *wanted, Vec *vecs)
450: {
451: struct DMCompositeLink *link;
452: PetscInt i, wnum;
453: DM_Composite *com = (DM_Composite *)dm->data;
454: PetscInt readonly;
455: PetscBool flg;
459: PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg);
461: if (!com->setup) DMSetUp(dm);
463: VecLockGet(pvec, &readonly);
464: for (i = 0, wnum = 0, link = com->next; link && wnum < nwanted; i++, link = link->next) {
465: if (!wanted || i == wanted[wnum]) {
466: VecResetArray(vecs[wnum]);
467: if (readonly) VecLockReadPop(vecs[wnum]);
468: DMRestoreLocalVector(link->dm, &vecs[wnum]);
469: wnum++;
470: }
471: }
472: return 0;
473: }
475: /*@C
476: DMCompositeScatter - Scatters from a global packed vector into its individual local vectors
478: Collective on dm
480: Input Parameters:
481: + dm - the packer object
482: . gvec - the global vector
483: - Vec ... - the individual sequential vectors, NULL for those that are not needed
485: Level: advanced
487: Notes:
488: DMCompositeScatterArray() is a non-variadic alternative that is often more convenient for library callers and is
489: accessible from Fortran.
491: .seealso `DMDestroy()`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`,
492: `DMCompositeGather()`, `DMCompositeCreate()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeGetAccess()`,
493: `DMCompositeGetLocalVectors()`, `DMCompositeRestoreLocalVectors()`, `DMCompositeGetEntries()`
494: `DMCompositeScatterArray()`
496: @*/
497: PetscErrorCode DMCompositeScatter(DM dm, Vec gvec, ...)
498: {
499: va_list Argp;
500: struct DMCompositeLink *next;
501: PETSC_UNUSED PetscInt cnt;
502: DM_Composite *com = (DM_Composite *)dm->data;
503: PetscBool flg;
507: PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg);
509: if (!com->setup) DMSetUp(dm);
511: /* loop over packed objects, handling one at at time */
512: va_start(Argp, gvec);
513: for (cnt = 3, next = com->next; next; cnt++, next = next->next) {
514: Vec local;
515: local = va_arg(Argp, Vec);
516: if (local) {
517: Vec global;
518: const PetscScalar *array;
520: DMGetGlobalVector(next->dm, &global);
521: VecGetArrayRead(gvec, &array);
522: VecPlaceArray(global, array + next->rstart);
523: DMGlobalToLocalBegin(next->dm, global, INSERT_VALUES, local);
524: DMGlobalToLocalEnd(next->dm, global, INSERT_VALUES, local);
525: VecRestoreArrayRead(gvec, &array);
526: VecResetArray(global);
527: DMRestoreGlobalVector(next->dm, &global);
528: }
529: }
530: va_end(Argp);
531: return 0;
532: }
534: /*@
535: DMCompositeScatterArray - Scatters from a global packed vector into its individual local vectors
537: Collective on dm
539: Input Parameters:
540: + dm - the packer object
541: . gvec - the global vector
542: - lvecs - array of local vectors, NULL for any that are not needed
544: Level: advanced
546: Note:
547: This is a non-variadic alternative to DMCompositeScatter()
549: .seealso `DMDestroy()`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`
550: `DMCompositeGather()`, `DMCompositeCreate()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeGetAccess()`,
551: `DMCompositeGetLocalVectors()`, `DMCompositeRestoreLocalVectors()`, `DMCompositeGetEntries()`
553: @*/
554: PetscErrorCode DMCompositeScatterArray(DM dm, Vec gvec, Vec *lvecs)
555: {
556: struct DMCompositeLink *next;
557: PetscInt i;
558: DM_Composite *com = (DM_Composite *)dm->data;
559: PetscBool flg;
563: PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg);
565: if (!com->setup) DMSetUp(dm);
567: /* loop over packed objects, handling one at at time */
568: for (i = 0, next = com->next; next; next = next->next, i++) {
569: if (lvecs[i]) {
570: Vec global;
571: const PetscScalar *array;
573: DMGetGlobalVector(next->dm, &global);
574: VecGetArrayRead(gvec, &array);
575: VecPlaceArray(global, (PetscScalar *)array + next->rstart);
576: DMGlobalToLocalBegin(next->dm, global, INSERT_VALUES, lvecs[i]);
577: DMGlobalToLocalEnd(next->dm, global, INSERT_VALUES, lvecs[i]);
578: VecRestoreArrayRead(gvec, &array);
579: VecResetArray(global);
580: DMRestoreGlobalVector(next->dm, &global);
581: }
582: }
583: return 0;
584: }
586: /*@C
587: DMCompositeGather - Gathers into a global packed vector from its individual local vectors
589: Collective on dm
591: Input Parameters:
592: + dm - the packer object
593: . gvec - the global vector
594: . imode - INSERT_VALUES or ADD_VALUES
595: - Vec ... - the individual sequential vectors, NULL for any that are not needed
597: Level: advanced
599: Not available from Fortran, Fortran users can use DMCompositeGatherArray()
601: .seealso `DMDestroy()`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`,
602: `DMCompositeScatter()`, `DMCompositeCreate()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeGetAccess()`,
603: `DMCompositeGetLocalVectors()`, `DMCompositeRestoreLocalVectors()`, `DMCompositeGetEntries()`
605: @*/
606: PetscErrorCode DMCompositeGather(DM dm, InsertMode imode, Vec gvec, ...)
607: {
608: va_list Argp;
609: struct DMCompositeLink *next;
610: DM_Composite *com = (DM_Composite *)dm->data;
611: PETSC_UNUSED PetscInt cnt;
612: PetscBool flg;
616: PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg);
618: if (!com->setup) DMSetUp(dm);
620: /* loop over packed objects, handling one at at time */
621: va_start(Argp, gvec);
622: for (cnt = 3, next = com->next; next; cnt++, next = next->next) {
623: Vec local;
624: local = va_arg(Argp, Vec);
625: if (local) {
626: PetscScalar *array;
627: Vec global;
629: DMGetGlobalVector(next->dm, &global);
630: VecGetArray(gvec, &array);
631: VecPlaceArray(global, array + next->rstart);
632: DMLocalToGlobalBegin(next->dm, local, imode, global);
633: DMLocalToGlobalEnd(next->dm, local, imode, global);
634: VecRestoreArray(gvec, &array);
635: VecResetArray(global);
636: DMRestoreGlobalVector(next->dm, &global);
637: }
638: }
639: va_end(Argp);
640: return 0;
641: }
643: /*@
644: DMCompositeGatherArray - Gathers into a global packed vector from its individual local vectors
646: Collective on dm
648: Input Parameters:
649: + dm - the packer object
650: . gvec - the global vector
651: . imode - INSERT_VALUES or ADD_VALUES
652: - lvecs - the individual sequential vectors, NULL for any that are not needed
654: Level: advanced
656: Notes:
657: This is a non-variadic alternative to DMCompositeGather().
659: .seealso `DMDestroy()`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`,
660: `DMCompositeScatter()`, `DMCompositeCreate()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeGetAccess()`,
661: `DMCompositeGetLocalVectors()`, `DMCompositeRestoreLocalVectors()`, `DMCompositeGetEntries()`,
662: @*/
663: PetscErrorCode DMCompositeGatherArray(DM dm, InsertMode imode, Vec gvec, Vec *lvecs)
664: {
665: struct DMCompositeLink *next;
666: DM_Composite *com = (DM_Composite *)dm->data;
667: PetscInt i;
668: PetscBool flg;
672: PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg);
674: if (!com->setup) DMSetUp(dm);
676: /* loop over packed objects, handling one at at time */
677: for (next = com->next, i = 0; next; next = next->next, i++) {
678: if (lvecs[i]) {
679: PetscScalar *array;
680: Vec global;
682: DMGetGlobalVector(next->dm, &global);
683: VecGetArray(gvec, &array);
684: VecPlaceArray(global, array + next->rstart);
685: DMLocalToGlobalBegin(next->dm, lvecs[i], imode, global);
686: DMLocalToGlobalEnd(next->dm, lvecs[i], imode, global);
687: VecRestoreArray(gvec, &array);
688: VecResetArray(global);
689: DMRestoreGlobalVector(next->dm, &global);
690: }
691: }
692: return 0;
693: }
695: /*@
696: DMCompositeAddDM - adds a DM vector to a DMComposite
698: Collective on dm
700: Input Parameters:
701: + dmc - the DMComposite (packer) object
702: - dm - the DM object
704: Level: advanced
706: .seealso `DMDestroy()`, `DMCompositeGather()`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`,
707: `DMCompositeScatter()`, `DMCompositeCreate()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeGetAccess()`,
708: `DMCompositeGetLocalVectors()`, `DMCompositeRestoreLocalVectors()`, `DMCompositeGetEntries()`
710: @*/
711: PetscErrorCode DMCompositeAddDM(DM dmc, DM dm)
712: {
713: PetscInt n, nlocal;
714: struct DMCompositeLink *mine, *next;
715: Vec global, local;
716: DM_Composite *com = (DM_Composite *)dmc->data;
717: PetscBool flg;
721: PetscObjectTypeCompare((PetscObject)dmc, DMCOMPOSITE, &flg);
723: next = com->next;
726: /* create new link */
727: PetscNew(&mine);
728: PetscObjectReference((PetscObject)dm);
729: DMGetGlobalVector(dm, &global);
730: VecGetLocalSize(global, &n);
731: DMRestoreGlobalVector(dm, &global);
732: DMGetLocalVector(dm, &local);
733: VecGetSize(local, &nlocal);
734: DMRestoreLocalVector(dm, &local);
736: mine->n = n;
737: mine->nlocal = nlocal;
738: mine->dm = dm;
739: mine->next = NULL;
740: com->n += n;
741: com->nghost += nlocal;
743: /* add to end of list */
744: if (!next) com->next = mine;
745: else {
746: while (next->next) next = next->next;
747: next->next = mine;
748: }
749: com->nDM++;
750: com->nmine++;
751: return 0;
752: }
754: #include <petscdraw.h>
755: PETSC_EXTERN PetscErrorCode VecView_MPI(Vec, PetscViewer);
756: PetscErrorCode VecView_DMComposite(Vec gvec, PetscViewer viewer)
757: {
758: DM dm;
759: struct DMCompositeLink *next;
760: PetscBool isdraw;
761: DM_Composite *com;
763: VecGetDM(gvec, &dm);
765: com = (DM_Composite *)dm->data;
766: next = com->next;
768: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw);
769: if (!isdraw) {
770: /* do I really want to call this? */
771: VecView_MPI(gvec, viewer);
772: } else {
773: PetscInt cnt = 0;
775: /* loop over packed objects, handling one at at time */
776: while (next) {
777: Vec vec;
778: const PetscScalar *array;
779: PetscInt bs;
781: /* Should use VecGetSubVector() eventually, but would need to forward the DM for that to work */
782: DMGetGlobalVector(next->dm, &vec);
783: VecGetArrayRead(gvec, &array);
784: VecPlaceArray(vec, (PetscScalar *)array + next->rstart);
785: VecRestoreArrayRead(gvec, &array);
786: VecView(vec, viewer);
787: VecResetArray(vec);
788: VecGetBlockSize(vec, &bs);
789: DMRestoreGlobalVector(next->dm, &vec);
790: PetscViewerDrawBaseAdd(viewer, bs);
791: cnt += bs;
792: next = next->next;
793: }
794: PetscViewerDrawBaseAdd(viewer, -cnt);
795: }
796: return 0;
797: }
799: PetscErrorCode DMCreateGlobalVector_Composite(DM dm, Vec *gvec)
800: {
801: DM_Composite *com = (DM_Composite *)dm->data;
804: DMSetFromOptions(dm);
805: DMSetUp(dm);
806: VecCreate(PetscObjectComm((PetscObject)dm), gvec);
807: VecSetType(*gvec, dm->vectype);
808: VecSetSizes(*gvec, com->n, com->N);
809: VecSetDM(*gvec, dm);
810: VecSetOperation(*gvec, VECOP_VIEW, (void (*)(void))VecView_DMComposite);
811: return 0;
812: }
814: PetscErrorCode DMCreateLocalVector_Composite(DM dm, Vec *lvec)
815: {
816: DM_Composite *com = (DM_Composite *)dm->data;
819: if (!com->setup) {
820: DMSetFromOptions(dm);
821: DMSetUp(dm);
822: }
823: VecCreate(PETSC_COMM_SELF, lvec);
824: VecSetType(*lvec, dm->vectype);
825: VecSetSizes(*lvec, com->nghost, PETSC_DECIDE);
826: VecSetDM(*lvec, dm);
827: return 0;
828: }
830: /*@C
831: DMCompositeGetISLocalToGlobalMappings - gets an ISLocalToGlobalMapping for each DM in the DMComposite, maps to the composite global space
833: Collective on DM
835: Input Parameter:
836: . dm - the packer object
838: Output Parameters:
839: . ltogs - the individual mappings for each packed vector. Note that this includes
840: all the ghost points that individual ghosted DMDA's may have.
842: Level: advanced
844: Notes:
845: Each entry of ltogs should be destroyed with ISLocalToGlobalMappingDestroy(), the ltogs array should be freed with PetscFree().
847: Not available from Fortran
849: .seealso `DMDestroy()`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`,
850: `DMCompositeGather()`, `DMCompositeCreate()`, `DMCompositeGetAccess()`, `DMCompositeScatter()`,
851: `DMCompositeGetLocalVectors()`, `DMCompositeRestoreLocalVectors()`, `DMCompositeGetEntries()`
853: @*/
854: PetscErrorCode DMCompositeGetISLocalToGlobalMappings(DM dm, ISLocalToGlobalMapping **ltogs)
855: {
856: PetscInt i, *idx, n, cnt;
857: struct DMCompositeLink *next;
858: PetscMPIInt rank;
859: DM_Composite *com = (DM_Composite *)dm->data;
860: PetscBool flg;
863: PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg);
865: DMSetUp(dm);
866: PetscMalloc1(com->nDM, ltogs);
867: next = com->next;
868: MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
870: /* loop over packed objects, handling one at at time */
871: cnt = 0;
872: while (next) {
873: ISLocalToGlobalMapping ltog;
874: PetscMPIInt size;
875: const PetscInt *suboff, *indices;
876: Vec global;
878: /* Get sub-DM global indices for each local dof */
879: DMGetLocalToGlobalMapping(next->dm, <og);
880: ISLocalToGlobalMappingGetSize(ltog, &n);
881: ISLocalToGlobalMappingGetIndices(ltog, &indices);
882: PetscMalloc1(n, &idx);
884: /* Get the offsets for the sub-DM global vector */
885: DMGetGlobalVector(next->dm, &global);
886: VecGetOwnershipRanges(global, &suboff);
887: MPI_Comm_size(PetscObjectComm((PetscObject)global), &size);
889: /* Shift the sub-DM definition of the global space to the composite global space */
890: for (i = 0; i < n; i++) {
891: PetscInt subi = indices[i], lo = 0, hi = size, t;
892: /* There's no consensus on what a negative index means,
893: except for skipping when setting the values in vectors and matrices */
894: if (subi < 0) {
895: idx[i] = subi - next->grstarts[rank];
896: continue;
897: }
898: /* Binary search to find which rank owns subi */
899: while (hi - lo > 1) {
900: t = lo + (hi - lo) / 2;
901: if (suboff[t] > subi) hi = t;
902: else lo = t;
903: }
904: idx[i] = subi - suboff[lo] + next->grstarts[lo];
905: }
906: ISLocalToGlobalMappingRestoreIndices(ltog, &indices);
907: ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)dm), 1, n, idx, PETSC_OWN_POINTER, &(*ltogs)[cnt]);
908: DMRestoreGlobalVector(next->dm, &global);
909: next = next->next;
910: cnt++;
911: }
912: return 0;
913: }
915: /*@C
916: DMCompositeGetLocalISs - Gets index sets for each component of a composite local vector
918: Not Collective
920: Input Parameter:
921: . dm - composite DM
923: Output Parameter:
924: . is - array of serial index sets for each each component of the DMComposite
926: Level: intermediate
928: Notes:
929: At present, a composite local vector does not normally exist. This function is used to provide index sets for
930: MatGetLocalSubMatrix(). In the future, the scatters for each entry in the DMComposite may be be merged into a single
931: scatter to a composite local vector. The user should not typically need to know which is being done.
933: To get the composite global indices at all local points (including ghosts), use DMCompositeGetISLocalToGlobalMappings().
935: To get index sets for pieces of the composite global vector, use DMCompositeGetGlobalISs().
937: Each returned IS should be destroyed with ISDestroy(), the array should be freed with PetscFree().
939: Not available from Fortran
941: .seealso: `DMCompositeGetGlobalISs()`, `DMCompositeGetISLocalToGlobalMappings()`, `MatGetLocalSubMatrix()`, `MatCreateLocalRef()`
942: @*/
943: PetscErrorCode DMCompositeGetLocalISs(DM dm, IS **is)
944: {
945: DM_Composite *com = (DM_Composite *)dm->data;
946: struct DMCompositeLink *link;
947: PetscInt cnt, start;
948: PetscBool flg;
952: PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg);
954: PetscMalloc1(com->nmine, is);
955: for (cnt = 0, start = 0, link = com->next; link; start += link->nlocal, cnt++, link = link->next) {
956: PetscInt bs;
957: ISCreateStride(PETSC_COMM_SELF, link->nlocal, start, 1, &(*is)[cnt]);
958: DMGetBlockSize(link->dm, &bs);
959: ISSetBlockSize((*is)[cnt], bs);
960: }
961: return 0;
962: }
964: /*@C
965: DMCompositeGetGlobalISs - Gets the index sets for each composed object
967: Collective on dm
969: Input Parameter:
970: . dm - the packer object
972: Output Parameters:
973: . is - the array of index sets
975: Level: advanced
977: Notes:
978: The is entries should be destroyed with ISDestroy(), the is array should be freed with PetscFree()
980: These could be used to extract a subset of vector entries for a "multi-physics" preconditioner
982: Use DMCompositeGetLocalISs() for index sets in the packed local numbering, and
983: DMCompositeGetISLocalToGlobalMappings() for to map local sub-DM (including ghost) indices to packed global
984: indices.
986: Fortran Notes:
988: The output argument 'is' must be an allocated array of sufficient length, which can be learned using DMCompositeGetNumberDM().
990: .seealso `DMDestroy()`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`,
991: `DMCompositeGather()`, `DMCompositeCreate()`, `DMCompositeGetAccess()`, `DMCompositeScatter()`,
992: `DMCompositeGetLocalVectors()`, `DMCompositeRestoreLocalVectors()`, `DMCompositeGetEntries()`
994: @*/
995: PetscErrorCode DMCompositeGetGlobalISs(DM dm, IS *is[])
996: {
997: PetscInt cnt = 0;
998: struct DMCompositeLink *next;
999: PetscMPIInt rank;
1000: DM_Composite *com = (DM_Composite *)dm->data;
1001: PetscBool flg;
1004: PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg);
1007: PetscMalloc1(com->nDM, is);
1008: next = com->next;
1009: MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
1011: /* loop over packed objects, handling one at at time */
1012: while (next) {
1013: PetscDS prob;
1015: ISCreateStride(PetscObjectComm((PetscObject)dm), next->n, next->grstart, 1, &(*is)[cnt]);
1016: DMGetDS(dm, &prob);
1017: if (prob) {
1018: MatNullSpace space;
1019: Mat pmat;
1020: PetscObject disc;
1021: PetscInt Nf;
1023: PetscDSGetNumFields(prob, &Nf);
1024: if (cnt < Nf) {
1025: PetscDSGetDiscretization(prob, cnt, &disc);
1026: PetscObjectQuery(disc, "nullspace", (PetscObject *)&space);
1027: if (space) PetscObjectCompose((PetscObject)(*is)[cnt], "nullspace", (PetscObject)space);
1028: PetscObjectQuery(disc, "nearnullspace", (PetscObject *)&space);
1029: if (space) PetscObjectCompose((PetscObject)(*is)[cnt], "nearnullspace", (PetscObject)space);
1030: PetscObjectQuery(disc, "pmat", (PetscObject *)&pmat);
1031: if (pmat) PetscObjectCompose((PetscObject)(*is)[cnt], "pmat", (PetscObject)pmat);
1032: }
1033: }
1034: cnt++;
1035: next = next->next;
1036: }
1037: return 0;
1038: }
1040: PetscErrorCode DMCreateFieldIS_Composite(DM dm, PetscInt *numFields, char ***fieldNames, IS **fields)
1041: {
1042: PetscInt nDM;
1043: DM *dms;
1044: PetscInt i;
1046: DMCompositeGetNumberDM(dm, &nDM);
1047: if (numFields) *numFields = nDM;
1048: DMCompositeGetGlobalISs(dm, fields);
1049: if (fieldNames) {
1050: PetscMalloc1(nDM, &dms);
1051: PetscMalloc1(nDM, fieldNames);
1052: DMCompositeGetEntriesArray(dm, dms);
1053: for (i = 0; i < nDM; i++) {
1054: char buf[256];
1055: const char *splitname;
1057: /* Split naming precedence: object name, prefix, number */
1058: splitname = ((PetscObject)dm)->name;
1059: if (!splitname) {
1060: PetscObjectGetOptionsPrefix((PetscObject)dms[i], &splitname);
1061: if (splitname) {
1062: size_t len;
1063: PetscStrncpy(buf, splitname, sizeof(buf));
1064: buf[sizeof(buf) - 1] = 0;
1065: PetscStrlen(buf, &len);
1066: if (buf[len - 1] == '_') buf[len - 1] = 0; /* Remove trailing underscore if it was used */
1067: splitname = buf;
1068: }
1069: }
1070: if (!splitname) {
1071: PetscSNPrintf(buf, sizeof(buf), "%" PetscInt_FMT, i);
1072: splitname = buf;
1073: }
1074: PetscStrallocpy(splitname, &(*fieldNames)[i]);
1075: }
1076: PetscFree(dms);
1077: }
1078: return 0;
1079: }
1081: /*
1082: This could take over from DMCreateFieldIS(), as it is more general,
1083: making DMCreateFieldIS() a special case -- calling with dmlist == NULL;
1084: At this point it's probably best to be less intrusive, however.
1085: */
1086: PetscErrorCode DMCreateFieldDecomposition_Composite(DM dm, PetscInt *len, char ***namelist, IS **islist, DM **dmlist)
1087: {
1088: PetscInt nDM;
1089: PetscInt i;
1091: DMCreateFieldIS_Composite(dm, len, namelist, islist);
1092: if (dmlist) {
1093: DMCompositeGetNumberDM(dm, &nDM);
1094: PetscMalloc1(nDM, dmlist);
1095: DMCompositeGetEntriesArray(dm, *dmlist);
1096: for (i = 0; i < nDM; i++) PetscObjectReference((PetscObject)((*dmlist)[i]));
1097: }
1098: return 0;
1099: }
1101: /* -------------------------------------------------------------------------------------*/
1102: /*@C
1103: DMCompositeGetLocalVectors - Gets local vectors for each part of a DMComposite.
1104: Use DMCompositeRestoreLocalVectors() to return them.
1106: Not Collective
1108: Input Parameter:
1109: . dm - the packer object
1111: Output Parameter:
1112: . Vec ... - the individual sequential Vecs
1114: Level: advanced
1116: Not available from Fortran
1118: .seealso `DMDestroy()`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`,
1119: `DMCompositeGather()`, `DMCompositeCreate()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeGetAccess()`,
1120: `DMCompositeRestoreLocalVectors()`, `DMCompositeScatter()`, `DMCompositeGetEntries()`
1122: @*/
1123: PetscErrorCode DMCompositeGetLocalVectors(DM dm, ...)
1124: {
1125: va_list Argp;
1126: struct DMCompositeLink *next;
1127: DM_Composite *com = (DM_Composite *)dm->data;
1128: PetscBool flg;
1131: PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg);
1133: next = com->next;
1134: /* loop over packed objects, handling one at at time */
1135: va_start(Argp, dm);
1136: while (next) {
1137: Vec *vec;
1138: vec = va_arg(Argp, Vec *);
1139: if (vec) DMGetLocalVector(next->dm, vec);
1140: next = next->next;
1141: }
1142: va_end(Argp);
1143: return 0;
1144: }
1146: /*@C
1147: DMCompositeRestoreLocalVectors - Restores local vectors for each part of a DMComposite.
1149: Not Collective
1151: Input Parameter:
1152: . dm - the packer object
1154: Output Parameter:
1155: . Vec ... - the individual sequential Vecs
1157: Level: advanced
1159: Not available from Fortran
1161: .seealso `DMDestroy()`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`,
1162: `DMCompositeGather()`, `DMCompositeCreate()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeGetAccess()`,
1163: `DMCompositeGetLocalVectors()`, `DMCompositeScatter()`, `DMCompositeGetEntries()`
1165: @*/
1166: PetscErrorCode DMCompositeRestoreLocalVectors(DM dm, ...)
1167: {
1168: va_list Argp;
1169: struct DMCompositeLink *next;
1170: DM_Composite *com = (DM_Composite *)dm->data;
1171: PetscBool flg;
1174: PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg);
1176: next = com->next;
1177: /* loop over packed objects, handling one at at time */
1178: va_start(Argp, dm);
1179: while (next) {
1180: Vec *vec;
1181: vec = va_arg(Argp, Vec *);
1182: if (vec) DMRestoreLocalVector(next->dm, vec);
1183: next = next->next;
1184: }
1185: va_end(Argp);
1186: return 0;
1187: }
1189: /* -------------------------------------------------------------------------------------*/
1190: /*@C
1191: DMCompositeGetEntries - Gets the DM for each entry in a DMComposite.
1193: Not Collective
1195: Input Parameter:
1196: . dm - the packer object
1198: Output Parameter:
1199: . DM ... - the individual entries (DMs)
1201: Level: advanced
1203: Fortran Notes:
1204: Available as DMCompositeGetEntries() for one output DM, DMCompositeGetEntries2() for 2, etc
1206: .seealso `DMDestroy()`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`, `DMCompositeGetEntriesArray()`
1207: `DMCompositeGather()`, `DMCompositeCreate()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeGetAccess()`,
1208: `DMCompositeRestoreLocalVectors()`, `DMCompositeGetLocalVectors()`, `DMCompositeScatter()`,
1209: `DMCompositeGetLocalVectors()`, `DMCompositeRestoreLocalVectors()`
1211: @*/
1212: PetscErrorCode DMCompositeGetEntries(DM dm, ...)
1213: {
1214: va_list Argp;
1215: struct DMCompositeLink *next;
1216: DM_Composite *com = (DM_Composite *)dm->data;
1217: PetscBool flg;
1220: PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg);
1222: next = com->next;
1223: /* loop over packed objects, handling one at at time */
1224: va_start(Argp, dm);
1225: while (next) {
1226: DM *dmn;
1227: dmn = va_arg(Argp, DM *);
1228: if (dmn) *dmn = next->dm;
1229: next = next->next;
1230: }
1231: va_end(Argp);
1232: return 0;
1233: }
1235: /*@C
1236: DMCompositeGetEntriesArray - Gets the DM for each entry in a DMComposite.
1238: Not Collective
1240: Input Parameter:
1241: . dm - the packer object
1243: Output Parameter:
1244: . dms - array of sufficient length (see DMCompositeGetNumberDM()) to hold the individual DMs
1246: Level: advanced
1248: .seealso `DMDestroy()`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`, `DMCompositeGetEntries()`
1249: `DMCompositeGather()`, `DMCompositeCreate()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeGetAccess()`,
1250: `DMCompositeRestoreLocalVectors()`, `DMCompositeGetLocalVectors()`, `DMCompositeScatter()`,
1251: `DMCompositeGetLocalVectors()`, `DMCompositeRestoreLocalVectors()`
1253: @*/
1254: PetscErrorCode DMCompositeGetEntriesArray(DM dm, DM dms[])
1255: {
1256: struct DMCompositeLink *next;
1257: DM_Composite *com = (DM_Composite *)dm->data;
1258: PetscInt i;
1259: PetscBool flg;
1262: PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg);
1264: /* loop over packed objects, handling one at at time */
1265: for (next = com->next, i = 0; next; next = next->next, i++) dms[i] = next->dm;
1266: return 0;
1267: }
1269: typedef struct {
1270: DM dm;
1271: PetscViewer *subv;
1272: Vec *vecs;
1273: } GLVisViewerCtx;
1275: static PetscErrorCode DestroyGLVisViewerCtx_Private(void *vctx)
1276: {
1277: GLVisViewerCtx *ctx = (GLVisViewerCtx *)vctx;
1278: PetscInt i, n;
1280: DMCompositeGetNumberDM(ctx->dm, &n);
1281: for (i = 0; i < n; i++) PetscViewerDestroy(&ctx->subv[i]);
1282: PetscFree2(ctx->subv, ctx->vecs);
1283: DMDestroy(&ctx->dm);
1284: PetscFree(ctx);
1285: return 0;
1286: }
1288: static PetscErrorCode DMCompositeSampleGLVisFields_Private(PetscObject oX, PetscInt nf, PetscObject oXfield[], void *vctx)
1289: {
1290: Vec X = (Vec)oX;
1291: GLVisViewerCtx *ctx = (GLVisViewerCtx *)vctx;
1292: PetscInt i, n, cumf;
1294: DMCompositeGetNumberDM(ctx->dm, &n);
1295: DMCompositeGetAccessArray(ctx->dm, X, n, NULL, ctx->vecs);
1296: for (i = 0, cumf = 0; i < n; i++) {
1297: PetscErrorCode (*g2l)(PetscObject, PetscInt, PetscObject[], void *);
1298: void *fctx;
1299: PetscInt nfi;
1301: PetscViewerGLVisGetFields_Private(ctx->subv[i], &nfi, NULL, NULL, &g2l, NULL, &fctx);
1302: if (!nfi) continue;
1303: if (g2l) (*g2l)((PetscObject)ctx->vecs[i], nfi, oXfield + cumf, fctx);
1304: else VecCopy(ctx->vecs[i], (Vec)(oXfield[cumf]));
1305: cumf += nfi;
1306: }
1307: DMCompositeRestoreAccessArray(ctx->dm, X, n, NULL, ctx->vecs);
1308: return 0;
1309: }
1311: static PetscErrorCode DMSetUpGLVisViewer_Composite(PetscObject odm, PetscViewer viewer)
1312: {
1313: DM dm = (DM)odm, *dms;
1314: Vec *Ufds;
1315: GLVisViewerCtx *ctx;
1316: PetscInt i, n, tnf, *sdim;
1317: char **fecs;
1319: PetscNew(&ctx);
1320: PetscObjectReference((PetscObject)dm);
1321: ctx->dm = dm;
1322: DMCompositeGetNumberDM(dm, &n);
1323: PetscMalloc1(n, &dms);
1324: DMCompositeGetEntriesArray(dm, dms);
1325: PetscMalloc2(n, &ctx->subv, n, &ctx->vecs);
1326: for (i = 0, tnf = 0; i < n; i++) {
1327: PetscInt nf;
1329: PetscViewerCreate(PetscObjectComm(odm), &ctx->subv[i]);
1330: PetscViewerSetType(ctx->subv[i], PETSCVIEWERGLVIS);
1331: PetscViewerGLVisSetDM_Private(ctx->subv[i], (PetscObject)dms[i]);
1332: PetscViewerGLVisGetFields_Private(ctx->subv[i], &nf, NULL, NULL, NULL, NULL, NULL);
1333: tnf += nf;
1334: }
1335: PetscFree(dms);
1336: PetscMalloc3(tnf, &fecs, tnf, &sdim, tnf, &Ufds);
1337: for (i = 0, tnf = 0; i < n; i++) {
1338: PetscInt *sd, nf, f;
1339: const char **fec;
1340: Vec *Uf;
1342: PetscViewerGLVisGetFields_Private(ctx->subv[i], &nf, &fec, &sd, NULL, (PetscObject **)&Uf, NULL);
1343: for (f = 0; f < nf; f++) {
1344: PetscStrallocpy(fec[f], &fecs[tnf + f]);
1345: Ufds[tnf + f] = Uf[f];
1346: sdim[tnf + f] = sd[f];
1347: }
1348: tnf += nf;
1349: }
1350: PetscViewerGLVisSetFields(viewer, tnf, (const char **)fecs, sdim, DMCompositeSampleGLVisFields_Private, (PetscObject *)Ufds, ctx, DestroyGLVisViewerCtx_Private);
1351: for (i = 0; i < tnf; i++) PetscFree(fecs[i]);
1352: PetscFree3(fecs, sdim, Ufds);
1353: return 0;
1354: }
1356: PetscErrorCode DMRefine_Composite(DM dmi, MPI_Comm comm, DM *fine)
1357: {
1358: struct DMCompositeLink *next;
1359: DM_Composite *com = (DM_Composite *)dmi->data;
1360: DM dm;
1363: if (comm == MPI_COMM_NULL) PetscObjectGetComm((PetscObject)dmi, &comm);
1364: DMSetUp(dmi);
1365: next = com->next;
1366: DMCompositeCreate(comm, fine);
1368: /* loop over packed objects, handling one at at time */
1369: while (next) {
1370: DMRefine(next->dm, comm, &dm);
1371: DMCompositeAddDM(*fine, dm);
1372: PetscObjectDereference((PetscObject)dm);
1373: next = next->next;
1374: }
1375: return 0;
1376: }
1378: PetscErrorCode DMCoarsen_Composite(DM dmi, MPI_Comm comm, DM *fine)
1379: {
1380: struct DMCompositeLink *next;
1381: DM_Composite *com = (DM_Composite *)dmi->data;
1382: DM dm;
1385: DMSetUp(dmi);
1386: if (comm == MPI_COMM_NULL) PetscObjectGetComm((PetscObject)dmi, &comm);
1387: next = com->next;
1388: DMCompositeCreate(comm, fine);
1390: /* loop over packed objects, handling one at at time */
1391: while (next) {
1392: DMCoarsen(next->dm, comm, &dm);
1393: DMCompositeAddDM(*fine, dm);
1394: PetscObjectDereference((PetscObject)dm);
1395: next = next->next;
1396: }
1397: return 0;
1398: }
1400: PetscErrorCode DMCreateInterpolation_Composite(DM coarse, DM fine, Mat *A, Vec *v)
1401: {
1402: PetscInt m, n, M, N, nDM, i;
1403: struct DMCompositeLink *nextc;
1404: struct DMCompositeLink *nextf;
1405: Vec gcoarse, gfine, *vecs;
1406: DM_Composite *comcoarse = (DM_Composite *)coarse->data;
1407: DM_Composite *comfine = (DM_Composite *)fine->data;
1408: Mat *mats;
1412: DMSetUp(coarse);
1413: DMSetUp(fine);
1414: /* use global vectors only for determining matrix layout */
1415: DMGetGlobalVector(coarse, &gcoarse);
1416: DMGetGlobalVector(fine, &gfine);
1417: VecGetLocalSize(gcoarse, &n);
1418: VecGetLocalSize(gfine, &m);
1419: VecGetSize(gcoarse, &N);
1420: VecGetSize(gfine, &M);
1421: DMRestoreGlobalVector(coarse, &gcoarse);
1422: DMRestoreGlobalVector(fine, &gfine);
1424: nDM = comfine->nDM;
1426: PetscCalloc1(nDM * nDM, &mats);
1427: if (v) PetscCalloc1(nDM, &vecs);
1429: /* loop over packed objects, handling one at at time */
1430: for (nextc = comcoarse->next, nextf = comfine->next, i = 0; nextc; nextc = nextc->next, nextf = nextf->next, i++) {
1431: if (!v) DMCreateInterpolation(nextc->dm, nextf->dm, &mats[i * nDM + i], NULL);
1432: else DMCreateInterpolation(nextc->dm, nextf->dm, &mats[i * nDM + i], &vecs[i]);
1433: }
1434: MatCreateNest(PetscObjectComm((PetscObject)fine), nDM, NULL, nDM, NULL, mats, A);
1435: if (v) VecCreateNest(PetscObjectComm((PetscObject)fine), nDM, NULL, vecs, v);
1436: for (i = 0; i < nDM * nDM; i++) MatDestroy(&mats[i]);
1437: PetscFree(mats);
1438: if (v) {
1439: for (i = 0; i < nDM; i++) VecDestroy(&vecs[i]);
1440: PetscFree(vecs);
1441: }
1442: return 0;
1443: }
1445: static PetscErrorCode DMGetLocalToGlobalMapping_Composite(DM dm)
1446: {
1447: DM_Composite *com = (DM_Composite *)dm->data;
1448: ISLocalToGlobalMapping *ltogs;
1449: PetscInt i;
1451: /* Set the ISLocalToGlobalMapping on the new matrix */
1452: DMCompositeGetISLocalToGlobalMappings(dm, <ogs);
1453: ISLocalToGlobalMappingConcatenate(PetscObjectComm((PetscObject)dm), com->nDM, ltogs, &dm->ltogmap);
1454: for (i = 0; i < com->nDM; i++) ISLocalToGlobalMappingDestroy(<ogs[i]);
1455: PetscFree(ltogs);
1456: return 0;
1457: }
1459: PetscErrorCode DMCreateColoring_Composite(DM dm, ISColoringType ctype, ISColoring *coloring)
1460: {
1461: PetscInt n, i, cnt;
1462: ISColoringValue *colors;
1463: PetscBool dense = PETSC_FALSE;
1464: ISColoringValue maxcol = 0;
1465: DM_Composite *com = (DM_Composite *)dm->data;
1469: if (ctype == IS_COLORING_GLOBAL) {
1470: n = com->n;
1471: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Unknown ISColoringType");
1472: PetscMalloc1(n, &colors); /* freed in ISColoringDestroy() */
1474: PetscOptionsGetBool(((PetscObject)dm)->options, ((PetscObject)dm)->prefix, "-dmcomposite_dense_jacobian", &dense, NULL);
1475: if (dense) {
1476: for (i = 0; i < n; i++) colors[i] = (ISColoringValue)(com->rstart + i);
1477: maxcol = com->N;
1478: } else {
1479: struct DMCompositeLink *next = com->next;
1480: PetscMPIInt rank;
1482: MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
1483: cnt = 0;
1484: while (next) {
1485: ISColoring lcoloring;
1487: DMCreateColoring(next->dm, IS_COLORING_GLOBAL, &lcoloring);
1488: for (i = 0; i < lcoloring->N; i++) colors[cnt++] = maxcol + lcoloring->colors[i];
1489: maxcol += lcoloring->n;
1490: ISColoringDestroy(&lcoloring);
1491: next = next->next;
1492: }
1493: }
1494: ISColoringCreate(PetscObjectComm((PetscObject)dm), maxcol, n, colors, PETSC_OWN_POINTER, coloring);
1495: return 0;
1496: }
1498: PetscErrorCode DMGlobalToLocalBegin_Composite(DM dm, Vec gvec, InsertMode mode, Vec lvec)
1499: {
1500: struct DMCompositeLink *next;
1501: PetscScalar *garray, *larray;
1502: DM_Composite *com = (DM_Composite *)dm->data;
1507: if (!com->setup) DMSetUp(dm);
1509: VecGetArray(gvec, &garray);
1510: VecGetArray(lvec, &larray);
1512: /* loop over packed objects, handling one at at time */
1513: next = com->next;
1514: while (next) {
1515: Vec local, global;
1516: PetscInt N;
1518: DMGetGlobalVector(next->dm, &global);
1519: VecGetLocalSize(global, &N);
1520: VecPlaceArray(global, garray);
1521: DMGetLocalVector(next->dm, &local);
1522: VecPlaceArray(local, larray);
1523: DMGlobalToLocalBegin(next->dm, global, mode, local);
1524: DMGlobalToLocalEnd(next->dm, global, mode, local);
1525: VecResetArray(global);
1526: VecResetArray(local);
1527: DMRestoreGlobalVector(next->dm, &global);
1528: DMRestoreLocalVector(next->dm, &local);
1530: larray += next->nlocal;
1531: garray += next->n;
1532: next = next->next;
1533: }
1535: VecRestoreArray(gvec, NULL);
1536: VecRestoreArray(lvec, NULL);
1537: return 0;
1538: }
1540: PetscErrorCode DMGlobalToLocalEnd_Composite(DM dm, Vec gvec, InsertMode mode, Vec lvec)
1541: {
1545: return 0;
1546: }
1548: PetscErrorCode DMLocalToGlobalBegin_Composite(DM dm, Vec lvec, InsertMode mode, Vec gvec)
1549: {
1550: struct DMCompositeLink *next;
1551: PetscScalar *larray, *garray;
1552: DM_Composite *com = (DM_Composite *)dm->data;
1558: if (!com->setup) DMSetUp(dm);
1560: VecGetArray(lvec, &larray);
1561: VecGetArray(gvec, &garray);
1563: /* loop over packed objects, handling one at at time */
1564: next = com->next;
1565: while (next) {
1566: Vec global, local;
1568: DMGetLocalVector(next->dm, &local);
1569: VecPlaceArray(local, larray);
1570: DMGetGlobalVector(next->dm, &global);
1571: VecPlaceArray(global, garray);
1572: DMLocalToGlobalBegin(next->dm, local, mode, global);
1573: DMLocalToGlobalEnd(next->dm, local, mode, global);
1574: VecResetArray(local);
1575: VecResetArray(global);
1576: DMRestoreGlobalVector(next->dm, &global);
1577: DMRestoreLocalVector(next->dm, &local);
1579: garray += next->n;
1580: larray += next->nlocal;
1581: next = next->next;
1582: }
1584: VecRestoreArray(gvec, NULL);
1585: VecRestoreArray(lvec, NULL);
1587: return 0;
1588: }
1590: PetscErrorCode DMLocalToGlobalEnd_Composite(DM dm, Vec lvec, InsertMode mode, Vec gvec)
1591: {
1595: return 0;
1596: }
1598: PetscErrorCode DMLocalToLocalBegin_Composite(DM dm, Vec vec1, InsertMode mode, Vec vec2)
1599: {
1600: struct DMCompositeLink *next;
1601: PetscScalar *array1, *array2;
1602: DM_Composite *com = (DM_Composite *)dm->data;
1608: if (!com->setup) DMSetUp(dm);
1610: VecGetArray(vec1, &array1);
1611: VecGetArray(vec2, &array2);
1613: /* loop over packed objects, handling one at at time */
1614: next = com->next;
1615: while (next) {
1616: Vec local1, local2;
1618: DMGetLocalVector(next->dm, &local1);
1619: VecPlaceArray(local1, array1);
1620: DMGetLocalVector(next->dm, &local2);
1621: VecPlaceArray(local2, array2);
1622: DMLocalToLocalBegin(next->dm, local1, mode, local2);
1623: DMLocalToLocalEnd(next->dm, local1, mode, local2);
1624: VecResetArray(local2);
1625: DMRestoreLocalVector(next->dm, &local2);
1626: VecResetArray(local1);
1627: DMRestoreLocalVector(next->dm, &local1);
1629: array1 += next->nlocal;
1630: array2 += next->nlocal;
1631: next = next->next;
1632: }
1634: VecRestoreArray(vec1, NULL);
1635: VecRestoreArray(vec2, NULL);
1637: return 0;
1638: }
1640: PetscErrorCode DMLocalToLocalEnd_Composite(DM dm, Vec lvec, InsertMode mode, Vec gvec)
1641: {
1645: return 0;
1646: }
1648: /*MC
1649: DMCOMPOSITE = "composite" - A DM object that is used to manage data for a collection of DMs
1651: Level: intermediate
1653: .seealso: `DMType`, `DM`, `DMDACreate()`, `DMCreate()`, `DMSetType()`, `DMCompositeCreate()`
1654: M*/
1656: PETSC_EXTERN PetscErrorCode DMCreate_Composite(DM p)
1657: {
1658: DM_Composite *com;
1660: PetscNew(&com);
1661: p->data = com;
1662: com->n = 0;
1663: com->nghost = 0;
1664: com->next = NULL;
1665: com->nDM = 0;
1667: p->ops->createglobalvector = DMCreateGlobalVector_Composite;
1668: p->ops->createlocalvector = DMCreateLocalVector_Composite;
1669: p->ops->getlocaltoglobalmapping = DMGetLocalToGlobalMapping_Composite;
1670: p->ops->createfieldis = DMCreateFieldIS_Composite;
1671: p->ops->createfielddecomposition = DMCreateFieldDecomposition_Composite;
1672: p->ops->refine = DMRefine_Composite;
1673: p->ops->coarsen = DMCoarsen_Composite;
1674: p->ops->createinterpolation = DMCreateInterpolation_Composite;
1675: p->ops->creatematrix = DMCreateMatrix_Composite;
1676: p->ops->getcoloring = DMCreateColoring_Composite;
1677: p->ops->globaltolocalbegin = DMGlobalToLocalBegin_Composite;
1678: p->ops->globaltolocalend = DMGlobalToLocalEnd_Composite;
1679: p->ops->localtoglobalbegin = DMLocalToGlobalBegin_Composite;
1680: p->ops->localtoglobalend = DMLocalToGlobalEnd_Composite;
1681: p->ops->localtolocalbegin = DMLocalToLocalBegin_Composite;
1682: p->ops->localtolocalend = DMLocalToLocalEnd_Composite;
1683: p->ops->destroy = DMDestroy_Composite;
1684: p->ops->view = DMView_Composite;
1685: p->ops->setup = DMSetUp_Composite;
1687: PetscObjectComposeFunction((PetscObject)p, "DMSetUpGLVisViewer_C", DMSetUpGLVisViewer_Composite);
1688: return 0;
1689: }
1691: /*@
1692: DMCompositeCreate - Creates a vector packer, used to generate "composite"
1693: vectors made up of several subvectors.
1695: Collective
1697: Input Parameter:
1698: . comm - the processors that will share the global vector
1700: Output Parameters:
1701: . packer - the packer object
1703: Level: advanced
1705: .seealso `DMDestroy()`, `DMCompositeAddDM()`, `DMCompositeScatter()`, `DMCOMPOSITE`, `DMCreate()`
1706: `DMCompositeGather()`, `DMCreateGlobalVector()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeGetAccess()`
1707: `DMCompositeGetLocalVectors()`, `DMCompositeRestoreLocalVectors()`, `DMCompositeGetEntries()`
1709: @*/
1710: PetscErrorCode DMCompositeCreate(MPI_Comm comm, DM *packer)
1711: {
1713: DMCreate(comm, packer);
1714: DMSetType(*packer, DMCOMPOSITE);
1715: return 0;
1716: }