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 `DM` in a `DMCOMPOSITE` to 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: Note:
20: See `DMSetApplicationContext()` and `DMGetApplicationContext()` for how to get user information into
21: this routine
23: Fortran Note:
24: Not available from Fortran
26: .seealso: `DMCOMPOSITE`, `DM`
27: @*/
28: PetscErrorCode DMCompositeSetCoupling(DM dm, PetscErrorCode (*FormCoupleLocations)(DM, Mat, PetscInt *, PetscInt *, PetscInt, PetscInt, PetscInt, PetscInt))
29: {
30: DM_Composite *com = (DM_Composite *)dm->data;
31: PetscBool flg;
33: PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg);
35: com->FormCoupleLocations = FormCoupleLocations;
36: return 0;
37: }
39: PetscErrorCode DMDestroy_Composite(DM dm)
40: {
41: struct DMCompositeLink *next, *prev;
42: DM_Composite *com = (DM_Composite *)dm->data;
44: next = com->next;
45: while (next) {
46: prev = next;
47: next = next->next;
48: DMDestroy(&prev->dm);
49: PetscFree(prev->grstarts);
50: PetscFree(prev);
51: }
52: PetscObjectComposeFunction((PetscObject)dm, "DMSetUpGLVisViewer_C", NULL);
53: /* This was originally freed in DMDestroy(), but that prevents reference counting of backend objects */
54: PetscFree(com);
55: return 0;
56: }
58: PetscErrorCode DMView_Composite(DM dm, PetscViewer v)
59: {
60: PetscBool iascii;
61: DM_Composite *com = (DM_Composite *)dm->data;
63: PetscObjectTypeCompare((PetscObject)v, PETSCVIEWERASCII, &iascii);
64: if (iascii) {
65: struct DMCompositeLink *lnk = com->next;
66: PetscInt i;
68: PetscViewerASCIIPrintf(v, "DM (%s)\n", ((PetscObject)dm)->prefix ? ((PetscObject)dm)->prefix : "no prefix");
69: PetscViewerASCIIPrintf(v, " contains %" PetscInt_FMT " DMs\n", com->nDM);
70: PetscViewerASCIIPushTab(v);
71: for (i = 0; lnk; lnk = lnk->next, i++) {
72: PetscViewerASCIIPrintf(v, "Link %" PetscInt_FMT ": DM of type %s\n", i, ((PetscObject)lnk->dm)->type_name);
73: PetscViewerASCIIPushTab(v);
74: DMView(lnk->dm, v);
75: PetscViewerASCIIPopTab(v);
76: }
77: PetscViewerASCIIPopTab(v);
78: }
79: return 0;
80: }
82: /* --------------------------------------------------------------------------------------*/
83: PetscErrorCode DMSetUp_Composite(DM dm)
84: {
85: PetscInt nprev = 0;
86: PetscMPIInt rank, size;
87: DM_Composite *com = (DM_Composite *)dm->data;
88: struct DMCompositeLink *next = com->next;
89: PetscLayout map;
92: PetscLayoutCreate(PetscObjectComm((PetscObject)dm), &map);
93: PetscLayoutSetLocalSize(map, com->n);
94: PetscLayoutSetSize(map, PETSC_DETERMINE);
95: PetscLayoutSetBlockSize(map, 1);
96: PetscLayoutSetUp(map);
97: PetscLayoutGetSize(map, &com->N);
98: PetscLayoutGetRange(map, &com->rstart, NULL);
99: PetscLayoutDestroy(&map);
101: /* now set the rstart for each linked vector */
102: MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
103: MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size);
104: while (next) {
105: next->rstart = nprev;
106: nprev += next->n;
107: next->grstart = com->rstart + next->rstart;
108: PetscMalloc1(size, &next->grstarts);
109: MPI_Allgather(&next->grstart, 1, MPIU_INT, next->grstarts, 1, MPIU_INT, PetscObjectComm((PetscObject)dm));
110: next = next->next;
111: }
112: com->setup = PETSC_TRUE;
113: return 0;
114: }
116: /* ----------------------------------------------------------------------------------*/
118: /*@
119: DMCompositeGetNumberDM - Get's the number of `DM` objects in the `DMCOMPOSITE`
120: representation.
122: Not Collective
124: Input Parameter:
125: . dm - the `DMCOMPOSITE` object
127: Output Parameter:
128: . nDM - the number of `DM`
130: Level: beginner
132: .seealso: `DMCOMPOSITE`, `DM`
133: @*/
134: PetscErrorCode DMCompositeGetNumberDM(DM dm, PetscInt *nDM)
135: {
136: DM_Composite *com = (DM_Composite *)dm->data;
137: PetscBool flg;
140: PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg);
142: *nDM = com->nDM;
143: return 0;
144: }
146: /*@C
147: DMCompositeGetAccess - Allows one to access the individual packed vectors in their global
148: representation.
150: Collective on dm
152: Input Parameters:
153: + dm - the `DMCOMPOSITE` object
154: - gvec - the global vector
156: Output Parameters:
157: . Vec* ... - the packed parallel vectors, NULL for those that are not needed
159: Level: advanced
161: Note:
162: Use `DMCompositeRestoreAccess()` to return the vectors when you no longer need them
164: Fortran Note:
165: Fortran callers must use numbered versions of this routine, e.g., DMCompositeGetAccess4(dm,gvec,vec1,vec2,vec3,vec4)
166: or use the alternative interface `DMCompositeGetAccessArray()`.
168: .seealso: `DMCOMPOSITE`, `DM`, `DMCompositeGetEntries()`, `DMCompositeScatter()`
169: @*/
170: PetscErrorCode DMCompositeGetAccess(DM dm, Vec gvec, ...)
171: {
172: va_list Argp;
173: struct DMCompositeLink *next;
174: DM_Composite *com = (DM_Composite *)dm->data;
175: PetscInt readonly;
176: PetscBool flg;
180: PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg);
182: next = com->next;
183: if (!com->setup) DMSetUp(dm);
185: VecLockGet(gvec, &readonly);
186: /* loop over packed objects, handling one at at time */
187: va_start(Argp, gvec);
188: while (next) {
189: Vec *vec;
190: vec = va_arg(Argp, Vec *);
191: if (vec) {
192: DMGetGlobalVector(next->dm, vec);
193: if (readonly) {
194: const PetscScalar *array;
195: VecGetArrayRead(gvec, &array);
196: VecPlaceArray(*vec, array + next->rstart);
197: VecLockReadPush(*vec);
198: VecRestoreArrayRead(gvec, &array);
199: } else {
200: PetscScalar *array;
201: VecGetArray(gvec, &array);
202: VecPlaceArray(*vec, array + next->rstart);
203: VecRestoreArray(gvec, &array);
204: }
205: }
206: next = next->next;
207: }
208: va_end(Argp);
209: return 0;
210: }
212: /*@C
213: DMCompositeGetAccessArray - Allows one to access the individual packed vectors in their global
214: representation.
216: Collective on dm
218: Input Parameters:
219: + dm - the `DMCOMPOSITE`
220: . pvec - packed vector
221: . nwanted - number of vectors wanted
222: - wanted - sorted array of vectors wanted, or NULL to get all vectors
224: Output Parameters:
225: . vecs - array of requested global vectors (must be allocated)
227: Level: advanced
229: Note:
230: Use `DMCompositeRestoreAccessArray()` to return the vectors when you no longer need them
232: .seealso: `DMCOMPOSITE`, `DM`, `DMCompositeGetAccess()`, `DMCompositeGetEntries()`, `DMCompositeScatter()`, `DMCompositeGather()`
233: @*/
234: PetscErrorCode DMCompositeGetAccessArray(DM dm, Vec pvec, PetscInt nwanted, const PetscInt *wanted, Vec *vecs)
235: {
236: struct DMCompositeLink *link;
237: PetscInt i, wnum;
238: DM_Composite *com = (DM_Composite *)dm->data;
239: PetscInt readonly;
240: PetscBool flg;
244: PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg);
246: if (!com->setup) DMSetUp(dm);
248: VecLockGet(pvec, &readonly);
249: for (i = 0, wnum = 0, link = com->next; link && wnum < nwanted; i++, link = link->next) {
250: if (!wanted || i == wanted[wnum]) {
251: Vec v;
252: DMGetGlobalVector(link->dm, &v);
253: if (readonly) {
254: const PetscScalar *array;
255: VecGetArrayRead(pvec, &array);
256: VecPlaceArray(v, array + link->rstart);
257: VecLockReadPush(v);
258: VecRestoreArrayRead(pvec, &array);
259: } else {
260: PetscScalar *array;
261: VecGetArray(pvec, &array);
262: VecPlaceArray(v, array + link->rstart);
263: VecRestoreArray(pvec, &array);
264: }
265: vecs[wnum++] = v;
266: }
267: }
268: return 0;
269: }
271: /*@C
272: DMCompositeGetLocalAccessArray - Allows one to access the individual
273: packed vectors in their local representation.
275: Collective on dm.
277: Input Parameters:
278: + dm - the `DMCOMPOSITE`
279: . pvec - packed vector
280: . nwanted - number of vectors wanted
281: - wanted - sorted array of vectors wanted, or NULL to get all vectors
283: Output Parameters:
284: . vecs - array of requested local vectors (must be allocated)
286: Level: advanced
288: Note:
289: Use `DMCompositeRestoreLocalAccessArray()` to return the vectors
290: when you no longer need them.
292: .seealso: `DMCOMPOSITE`, `DM`, `DMCompositeRestoreLocalAccessArray()`, `DMCompositeGetAccess()`,
293: `DMCompositeGetEntries()`, `DMCompositeScatter()`, `DMCompositeGather()`
294: @*/
295: PetscErrorCode DMCompositeGetLocalAccessArray(DM dm, Vec pvec, PetscInt nwanted, const PetscInt *wanted, Vec *vecs)
296: {
297: struct DMCompositeLink *link;
298: PetscInt i, wnum;
299: DM_Composite *com = (DM_Composite *)dm->data;
300: PetscInt readonly;
301: PetscInt nlocal = 0;
302: PetscBool flg;
306: PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg);
308: if (!com->setup) DMSetUp(dm);
310: VecLockGet(pvec, &readonly);
311: for (i = 0, wnum = 0, link = com->next; link && wnum < nwanted; i++, link = link->next) {
312: if (!wanted || i == wanted[wnum]) {
313: Vec v;
314: DMGetLocalVector(link->dm, &v);
315: if (readonly) {
316: const PetscScalar *array;
317: VecGetArrayRead(pvec, &array);
318: VecPlaceArray(v, array + nlocal);
319: // 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
320: VecLockReadPush(v);
321: VecRestoreArrayRead(pvec, &array);
322: } else {
323: PetscScalar *array;
324: VecGetArray(pvec, &array);
325: VecPlaceArray(v, array + nlocal);
326: VecRestoreArray(pvec, &array);
327: }
328: vecs[wnum++] = v;
329: }
331: nlocal += link->nlocal;
332: }
334: return 0;
335: }
337: /*@C
338: DMCompositeRestoreAccess - Returns the vectors obtained with `DMCompositeGetAccess()`
339: representation.
341: Collective on dm
343: Input Parameters:
344: + dm - the `DMCOMPOSITE` object
345: . gvec - the global vector
346: - Vec* ... - the individual parallel vectors, NULL for those that are not needed
348: Level: advanced
350: .seealso: `DMCOMPOSITE`, `DM`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`,
351: `DMCompositeGather()`, `DMCompositeCreate()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeScatter()`,
352: `DMCompositeRestoreAccess()`, `DMCompositeGetAccess()`
353: @*/
354: PetscErrorCode DMCompositeRestoreAccess(DM dm, Vec gvec, ...)
355: {
356: va_list Argp;
357: struct DMCompositeLink *next;
358: DM_Composite *com = (DM_Composite *)dm->data;
359: PetscInt readonly;
360: PetscBool flg;
364: PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg);
366: next = com->next;
367: if (!com->setup) DMSetUp(dm);
369: VecLockGet(gvec, &readonly);
370: /* loop over packed objects, handling one at at time */
371: va_start(Argp, gvec);
372: while (next) {
373: Vec *vec;
374: vec = va_arg(Argp, Vec *);
375: if (vec) {
376: VecResetArray(*vec);
377: if (readonly) VecLockReadPop(*vec);
378: DMRestoreGlobalVector(next->dm, vec);
379: }
380: next = next->next;
381: }
382: va_end(Argp);
383: return 0;
384: }
386: /*@C
387: DMCompositeRestoreAccessArray - Returns the vectors obtained with `DMCompositeGetAccessArray()`
389: Collective on dm
391: Input Parameters:
392: + dm - the `DMCOMPOSITE` object
393: . pvec - packed vector
394: . nwanted - number of vectors wanted
395: . wanted - sorted array of vectors wanted, or NULL to get all vectors
396: - vecs - array of global vectors to return
398: Level: advanced
400: .seealso: `DMCOMPOSITE`, `DM`, `DMCompositeRestoreAccess()`, `DMCompositeRestoreEntries()`, `DMCompositeScatter()`, `DMCompositeGather()`
401: @*/
402: PetscErrorCode DMCompositeRestoreAccessArray(DM dm, Vec pvec, PetscInt nwanted, const PetscInt *wanted, Vec *vecs)
403: {
404: struct DMCompositeLink *link;
405: PetscInt i, wnum;
406: DM_Composite *com = (DM_Composite *)dm->data;
407: PetscInt readonly;
408: PetscBool flg;
412: PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg);
414: if (!com->setup) DMSetUp(dm);
416: VecLockGet(pvec, &readonly);
417: for (i = 0, wnum = 0, link = com->next; link && wnum < nwanted; i++, link = link->next) {
418: if (!wanted || i == wanted[wnum]) {
419: VecResetArray(vecs[wnum]);
420: if (readonly) VecLockReadPop(vecs[wnum]);
421: DMRestoreGlobalVector(link->dm, &vecs[wnum]);
422: wnum++;
423: }
424: }
425: return 0;
426: }
428: /*@C
429: DMCompositeRestoreLocalAccessArray - Returns the vectors obtained with `DMCompositeGetLocalAccessArray()`.
431: Collective on dm.
433: Input Parameters:
434: + dm - the `DMCOMPOSITE` object
435: . pvec - packed vector
436: . nwanted - number of vectors wanted
437: . wanted - sorted array of vectors wanted, or NULL to restore all vectors
438: - vecs - array of local vectors to return
440: Level: advanced
442: Note:
443: nwanted and wanted must match the values given to `DMCompositeGetLocalAccessArray()`
444: otherwise the call will fail.
446: .seealso: `DMCOMPOSITE`, `DM`, `DMCompositeGetLocalAccessArray()`, `DMCompositeRestoreAccessArray()`,
447: `DMCompositeRestoreAccess()`, `DMCompositeRestoreEntries()`,
448: `DMCompositeScatter()`, `DMCompositeGather()`
449: @*/
450: PetscErrorCode DMCompositeRestoreLocalAccessArray(DM dm, Vec pvec, PetscInt nwanted, const PetscInt *wanted, Vec *vecs)
451: {
452: struct DMCompositeLink *link;
453: PetscInt i, wnum;
454: DM_Composite *com = (DM_Composite *)dm->data;
455: PetscInt readonly;
456: PetscBool flg;
460: PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg);
462: if (!com->setup) DMSetUp(dm);
464: VecLockGet(pvec, &readonly);
465: for (i = 0, wnum = 0, link = com->next; link && wnum < nwanted; i++, link = link->next) {
466: if (!wanted || i == wanted[wnum]) {
467: VecResetArray(vecs[wnum]);
468: if (readonly) VecLockReadPop(vecs[wnum]);
469: DMRestoreLocalVector(link->dm, &vecs[wnum]);
470: wnum++;
471: }
472: }
473: return 0;
474: }
476: /*@C
477: DMCompositeScatter - Scatters from a global packed vector into its individual local vectors
479: Collective on dm
481: Input Parameters:
482: + dm - the `DMCOMPOSITE` object
483: . gvec - the global vector
484: - Vec ... - the individual sequential vectors, NULL for those that are not needed
486: Level: advanced
488: Note:
489: `DMCompositeScatterArray()` is a non-variadic alternative that is often more convenient for library callers and is
490: accessible from Fortran.
492: .seealso: `DMCOMPOSITE`, `DM`, `DMDestroy()`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`,
493: `DMCompositeGather()`, `DMCompositeCreate()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeGetAccess()`,
494: `DMCompositeGetLocalVectors()`, `DMCompositeRestoreLocalVectors()`, `DMCompositeGetEntries()`
495: `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 `DMCOMPOSITE` 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: `DMCOMPOSITE`, `DM`, `DMDestroy()`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`
550: `DMCompositeGather()`, `DMCompositeCreate()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeGetAccess()`,
551: `DMCompositeGetLocalVectors()`, `DMCompositeRestoreLocalVectors()`, `DMCompositeGetEntries()`
552: @*/
553: PetscErrorCode DMCompositeScatterArray(DM dm, Vec gvec, Vec *lvecs)
554: {
555: struct DMCompositeLink *next;
556: PetscInt i;
557: DM_Composite *com = (DM_Composite *)dm->data;
558: PetscBool flg;
562: PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg);
564: if (!com->setup) DMSetUp(dm);
566: /* loop over packed objects, handling one at at time */
567: for (i = 0, next = com->next; next; next = next->next, i++) {
568: if (lvecs[i]) {
569: Vec global;
570: const PetscScalar *array;
572: DMGetGlobalVector(next->dm, &global);
573: VecGetArrayRead(gvec, &array);
574: VecPlaceArray(global, (PetscScalar *)array + next->rstart);
575: DMGlobalToLocalBegin(next->dm, global, INSERT_VALUES, lvecs[i]);
576: DMGlobalToLocalEnd(next->dm, global, INSERT_VALUES, lvecs[i]);
577: VecRestoreArrayRead(gvec, &array);
578: VecResetArray(global);
579: DMRestoreGlobalVector(next->dm, &global);
580: }
581: }
582: return 0;
583: }
585: /*@C
586: DMCompositeGather - Gathers into a global packed vector from its individual local vectors
588: Collective on dm
590: Input Parameters:
591: + dm - the `DMCOMPOSITE` object
592: . gvec - the global vector
593: . imode - `INSERT_VALUES` or `ADD_VALUES`
594: - Vec ... - the individual sequential vectors, NULL for any that are not needed
596: Level: advanced
598: Fortran Note:
599: Fortran users should use `DMCompositeGatherArray()`
601: .seealso: `DMCOMPOSITE`, `DM`, `DMDestroy()`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`,
602: `DMCompositeScatter()`, `DMCompositeCreate()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeGetAccess()`,
603: `DMCompositeGetLocalVectors()`, `DMCompositeRestoreLocalVectors()`, `DMCompositeGetEntries()`
604: @*/
605: PetscErrorCode DMCompositeGather(DM dm, InsertMode imode, Vec gvec, ...)
606: {
607: va_list Argp;
608: struct DMCompositeLink *next;
609: DM_Composite *com = (DM_Composite *)dm->data;
610: PETSC_UNUSED PetscInt cnt;
611: PetscBool flg;
615: PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg);
617: if (!com->setup) DMSetUp(dm);
619: /* loop over packed objects, handling one at at time */
620: va_start(Argp, gvec);
621: for (cnt = 3, next = com->next; next; cnt++, next = next->next) {
622: Vec local;
623: local = va_arg(Argp, Vec);
624: if (local) {
625: PetscScalar *array;
626: Vec global;
628: DMGetGlobalVector(next->dm, &global);
629: VecGetArray(gvec, &array);
630: VecPlaceArray(global, array + next->rstart);
631: DMLocalToGlobalBegin(next->dm, local, imode, global);
632: DMLocalToGlobalEnd(next->dm, local, imode, global);
633: VecRestoreArray(gvec, &array);
634: VecResetArray(global);
635: DMRestoreGlobalVector(next->dm, &global);
636: }
637: }
638: va_end(Argp);
639: return 0;
640: }
642: /*@
643: DMCompositeGatherArray - Gathers into a global packed vector from its individual local vectors
645: Collective on dm
647: Input Parameters:
648: + dm - the `DMCOMPOSITE` object
649: . gvec - the global vector
650: . imode - `INSERT_VALUES` or `ADD_VALUES`
651: - lvecs - the individual sequential vectors, NULL for any that are not needed
653: Level: advanced
655: Note:
656: This is a non-variadic alternative to `DMCompositeGather()`.
658: .seealso: `DMCOMPOSITE`, `DM`, `DMDestroy()`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`,
659: `DMCompositeScatter()`, `DMCompositeCreate()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeGetAccess()`,
660: `DMCompositeGetLocalVectors()`, `DMCompositeRestoreLocalVectors()`, `DMCompositeGetEntries()`,
661: @*/
662: PetscErrorCode DMCompositeGatherArray(DM dm, InsertMode imode, Vec gvec, Vec *lvecs)
663: {
664: struct DMCompositeLink *next;
665: DM_Composite *com = (DM_Composite *)dm->data;
666: PetscInt i;
667: PetscBool flg;
671: PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg);
673: if (!com->setup) DMSetUp(dm);
675: /* loop over packed objects, handling one at at time */
676: for (next = com->next, i = 0; next; next = next->next, i++) {
677: if (lvecs[i]) {
678: PetscScalar *array;
679: Vec global;
681: DMGetGlobalVector(next->dm, &global);
682: VecGetArray(gvec, &array);
683: VecPlaceArray(global, array + next->rstart);
684: DMLocalToGlobalBegin(next->dm, lvecs[i], imode, global);
685: DMLocalToGlobalEnd(next->dm, lvecs[i], imode, global);
686: VecRestoreArray(gvec, &array);
687: VecResetArray(global);
688: DMRestoreGlobalVector(next->dm, &global);
689: }
690: }
691: return 0;
692: }
694: /*@
695: DMCompositeAddDM - adds a `DM` vector to a `DMCOMPOSITE`
697: Collective on dm
699: Input Parameters:
700: + dmc - the `DMCOMPOSITE` object
701: - dm - the `DM` object
703: Level: advanced
705: .seealso: `DMCOMPOSITE`, `DM`, `DMDestroy()`, `DMCompositeGather()`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`,
706: `DMCompositeScatter()`, `DMCompositeCreate()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeGetAccess()`,
707: `DMCompositeGetLocalVectors()`, `DMCompositeRestoreLocalVectors()`, `DMCompositeGetEntries()`
708: @*/
709: PetscErrorCode DMCompositeAddDM(DM dmc, DM dm)
710: {
711: PetscInt n, nlocal;
712: struct DMCompositeLink *mine, *next;
713: Vec global, local;
714: DM_Composite *com = (DM_Composite *)dmc->data;
715: PetscBool flg;
719: PetscObjectTypeCompare((PetscObject)dmc, DMCOMPOSITE, &flg);
721: next = com->next;
724: /* create new link */
725: PetscNew(&mine);
726: PetscObjectReference((PetscObject)dm);
727: DMGetGlobalVector(dm, &global);
728: VecGetLocalSize(global, &n);
729: DMRestoreGlobalVector(dm, &global);
730: DMGetLocalVector(dm, &local);
731: VecGetSize(local, &nlocal);
732: DMRestoreLocalVector(dm, &local);
734: mine->n = n;
735: mine->nlocal = nlocal;
736: mine->dm = dm;
737: mine->next = NULL;
738: com->n += n;
739: com->nghost += nlocal;
741: /* add to end of list */
742: if (!next) com->next = mine;
743: else {
744: while (next->next) next = next->next;
745: next->next = mine;
746: }
747: com->nDM++;
748: com->nmine++;
749: return 0;
750: }
752: #include <petscdraw.h>
753: PETSC_EXTERN PetscErrorCode VecView_MPI(Vec, PetscViewer);
754: PetscErrorCode VecView_DMComposite(Vec gvec, PetscViewer viewer)
755: {
756: DM dm;
757: struct DMCompositeLink *next;
758: PetscBool isdraw;
759: DM_Composite *com;
761: VecGetDM(gvec, &dm);
763: com = (DM_Composite *)dm->data;
764: next = com->next;
766: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw);
767: if (!isdraw) {
768: /* do I really want to call this? */
769: VecView_MPI(gvec, viewer);
770: } else {
771: PetscInt cnt = 0;
773: /* loop over packed objects, handling one at at time */
774: while (next) {
775: Vec vec;
776: const PetscScalar *array;
777: PetscInt bs;
779: /* Should use VecGetSubVector() eventually, but would need to forward the DM for that to work */
780: DMGetGlobalVector(next->dm, &vec);
781: VecGetArrayRead(gvec, &array);
782: VecPlaceArray(vec, (PetscScalar *)array + next->rstart);
783: VecRestoreArrayRead(gvec, &array);
784: VecView(vec, viewer);
785: VecResetArray(vec);
786: VecGetBlockSize(vec, &bs);
787: DMRestoreGlobalVector(next->dm, &vec);
788: PetscViewerDrawBaseAdd(viewer, bs);
789: cnt += bs;
790: next = next->next;
791: }
792: PetscViewerDrawBaseAdd(viewer, -cnt);
793: }
794: return 0;
795: }
797: PetscErrorCode DMCreateGlobalVector_Composite(DM dm, Vec *gvec)
798: {
799: DM_Composite *com = (DM_Composite *)dm->data;
802: DMSetFromOptions(dm);
803: DMSetUp(dm);
804: VecCreate(PetscObjectComm((PetscObject)dm), gvec);
805: VecSetType(*gvec, dm->vectype);
806: VecSetSizes(*gvec, com->n, com->N);
807: VecSetDM(*gvec, dm);
808: VecSetOperation(*gvec, VECOP_VIEW, (void (*)(void))VecView_DMComposite);
809: return 0;
810: }
812: PetscErrorCode DMCreateLocalVector_Composite(DM dm, Vec *lvec)
813: {
814: DM_Composite *com = (DM_Composite *)dm->data;
817: if (!com->setup) {
818: DMSetFromOptions(dm);
819: DMSetUp(dm);
820: }
821: VecCreate(PETSC_COMM_SELF, lvec);
822: VecSetType(*lvec, dm->vectype);
823: VecSetSizes(*lvec, com->nghost, PETSC_DECIDE);
824: VecSetDM(*lvec, dm);
825: return 0;
826: }
828: /*@C
829: DMCompositeGetISLocalToGlobalMappings - gets an `ISLocalToGlobalMapping` for each `DM` in the `DMCOMPOSITE`, maps to the composite global space
831: Collective on dm
833: Input Parameter:
834: . dm - the `DMCOMPOSITE` object
836: Output Parameters:
837: . ltogs - the individual mappings for each packed vector. Note that this includes
838: all the ghost points that individual ghosted `DMDA` may have.
840: Level: advanced
842: Note:
843: Each entry of ltogs should be destroyed with `ISLocalToGlobalMappingDestroy()`, the ltogs array should be freed with `PetscFree()`.
845: Fortran Note:
846: Not available from Fortran
848: .seealso: `DMCOMPOSITE`, `DM`, `DMDestroy()`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`,
849: `DMCompositeGather()`, `DMCompositeCreate()`, `DMCompositeGetAccess()`, `DMCompositeScatter()`,
850: `DMCompositeGetLocalVectors()`, `DMCompositeRestoreLocalVectors()`, `DMCompositeGetEntries()`
851: @*/
852: PetscErrorCode DMCompositeGetISLocalToGlobalMappings(DM dm, ISLocalToGlobalMapping **ltogs)
853: {
854: PetscInt i, *idx, n, cnt;
855: struct DMCompositeLink *next;
856: PetscMPIInt rank;
857: DM_Composite *com = (DM_Composite *)dm->data;
858: PetscBool flg;
861: PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg);
863: DMSetUp(dm);
864: PetscMalloc1(com->nDM, ltogs);
865: next = com->next;
866: MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
868: /* loop over packed objects, handling one at at time */
869: cnt = 0;
870: while (next) {
871: ISLocalToGlobalMapping ltog;
872: PetscMPIInt size;
873: const PetscInt *suboff, *indices;
874: Vec global;
876: /* Get sub-DM global indices for each local dof */
877: DMGetLocalToGlobalMapping(next->dm, <og);
878: ISLocalToGlobalMappingGetSize(ltog, &n);
879: ISLocalToGlobalMappingGetIndices(ltog, &indices);
880: PetscMalloc1(n, &idx);
882: /* Get the offsets for the sub-DM global vector */
883: DMGetGlobalVector(next->dm, &global);
884: VecGetOwnershipRanges(global, &suboff);
885: MPI_Comm_size(PetscObjectComm((PetscObject)global), &size);
887: /* Shift the sub-DM definition of the global space to the composite global space */
888: for (i = 0; i < n; i++) {
889: PetscInt subi = indices[i], lo = 0, hi = size, t;
890: /* There's no consensus on what a negative index means,
891: except for skipping when setting the values in vectors and matrices */
892: if (subi < 0) {
893: idx[i] = subi - next->grstarts[rank];
894: continue;
895: }
896: /* Binary search to find which rank owns subi */
897: while (hi - lo > 1) {
898: t = lo + (hi - lo) / 2;
899: if (suboff[t] > subi) hi = t;
900: else lo = t;
901: }
902: idx[i] = subi - suboff[lo] + next->grstarts[lo];
903: }
904: ISLocalToGlobalMappingRestoreIndices(ltog, &indices);
905: ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)dm), 1, n, idx, PETSC_OWN_POINTER, &(*ltogs)[cnt]);
906: DMRestoreGlobalVector(next->dm, &global);
907: next = next->next;
908: cnt++;
909: }
910: return 0;
911: }
913: /*@C
914: DMCompositeGetLocalISs - Gets index sets for each component of a composite local vector
916: Not Collective
918: Input Parameter:
919: . dm - the `DMCOMPOSITE`
921: Output Parameter:
922: . is - array of serial index sets for each each component of the `DMCOMPOSITE`
924: Level: intermediate
926: Notes:
927: At present, a composite local vector does not normally exist. This function is used to provide index sets for
928: `MatGetLocalSubMatrix()`. In the future, the scatters for each entry in the `DMCOMPOSITE` may be be merged into a single
929: scatter to a composite local vector. The user should not typically need to know which is being done.
931: To get the composite global indices at all local points (including ghosts), use `DMCompositeGetISLocalToGlobalMappings()`.
933: To get index sets for pieces of the composite global vector, use `DMCompositeGetGlobalISs()`.
935: Each returned `IS` should be destroyed with `ISDestroy()`, the array should be freed with `PetscFree()`.
937: Fortran Note:
938: Not available from Fortran
940: .seealso: `DMCOMPOSITE`, `DM`, `DMCompositeGetGlobalISs()`, `DMCompositeGetISLocalToGlobalMappings()`, `MatGetLocalSubMatrix()`, `MatCreateLocalRef()`
941: @*/
942: PetscErrorCode DMCompositeGetLocalISs(DM dm, IS **is)
943: {
944: DM_Composite *com = (DM_Composite *)dm->data;
945: struct DMCompositeLink *link;
946: PetscInt cnt, start;
947: PetscBool flg;
951: PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg);
953: PetscMalloc1(com->nmine, is);
954: for (cnt = 0, start = 0, link = com->next; link; start += link->nlocal, cnt++, link = link->next) {
955: PetscInt bs;
956: ISCreateStride(PETSC_COMM_SELF, link->nlocal, start, 1, &(*is)[cnt]);
957: DMGetBlockSize(link->dm, &bs);
958: ISSetBlockSize((*is)[cnt], bs);
959: }
960: return 0;
961: }
963: /*@C
964: DMCompositeGetGlobalISs - Gets the index sets for each composed object in a `DMCOMPOSITE`
966: Collective on dm
968: Input Parameter:
969: . dm - the `DMCOMPOSITE` object
971: Output Parameters:
972: . is - the array of index sets
974: Level: advanced
976: Notes:
977: The is entries should be destroyed with `ISDestroy()`, the is array should be freed with `PetscFree()`
979: These could be used to extract a subset of vector entries for a "multi-physics" preconditioner
981: Use `DMCompositeGetLocalISs()` for index sets in the packed local numbering, and
982: `DMCompositeGetISLocalToGlobalMappings()` for to map local sub-`DM` (including ghost) indices to packed global
983: indices.
985: Fortran Note:
986: The output argument 'is' must be an allocated array of sufficient length, which can be learned using `DMCompositeGetNumberDM()`.
988: .seealso: `DMCOMPOSITE`, `DM`, `DMDestroy()`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`,
989: `DMCompositeGather()`, `DMCompositeCreate()`, `DMCompositeGetAccess()`, `DMCompositeScatter()`,
990: `DMCompositeGetLocalVectors()`, `DMCompositeRestoreLocalVectors()`, `DMCompositeGetEntries()`
991: @*/
992: PetscErrorCode DMCompositeGetGlobalISs(DM dm, IS *is[])
993: {
994: PetscInt cnt = 0;
995: struct DMCompositeLink *next;
996: PetscMPIInt rank;
997: DM_Composite *com = (DM_Composite *)dm->data;
998: PetscBool flg;
1001: PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg);
1004: PetscMalloc1(com->nDM, is);
1005: next = com->next;
1006: MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
1008: /* loop over packed objects, handling one at at time */
1009: while (next) {
1010: PetscDS prob;
1012: ISCreateStride(PetscObjectComm((PetscObject)dm), next->n, next->grstart, 1, &(*is)[cnt]);
1013: DMGetDS(dm, &prob);
1014: if (prob) {
1015: MatNullSpace space;
1016: Mat pmat;
1017: PetscObject disc;
1018: PetscInt Nf;
1020: PetscDSGetNumFields(prob, &Nf);
1021: if (cnt < Nf) {
1022: PetscDSGetDiscretization(prob, cnt, &disc);
1023: PetscObjectQuery(disc, "nullspace", (PetscObject *)&space);
1024: if (space) PetscObjectCompose((PetscObject)(*is)[cnt], "nullspace", (PetscObject)space);
1025: PetscObjectQuery(disc, "nearnullspace", (PetscObject *)&space);
1026: if (space) PetscObjectCompose((PetscObject)(*is)[cnt], "nearnullspace", (PetscObject)space);
1027: PetscObjectQuery(disc, "pmat", (PetscObject *)&pmat);
1028: if (pmat) PetscObjectCompose((PetscObject)(*is)[cnt], "pmat", (PetscObject)pmat);
1029: }
1030: }
1031: cnt++;
1032: next = next->next;
1033: }
1034: return 0;
1035: }
1037: PetscErrorCode DMCreateFieldIS_Composite(DM dm, PetscInt *numFields, char ***fieldNames, IS **fields)
1038: {
1039: PetscInt nDM;
1040: DM *dms;
1041: PetscInt i;
1043: DMCompositeGetNumberDM(dm, &nDM);
1044: if (numFields) *numFields = nDM;
1045: DMCompositeGetGlobalISs(dm, fields);
1046: if (fieldNames) {
1047: PetscMalloc1(nDM, &dms);
1048: PetscMalloc1(nDM, fieldNames);
1049: DMCompositeGetEntriesArray(dm, dms);
1050: for (i = 0; i < nDM; i++) {
1051: char buf[256];
1052: const char *splitname;
1054: /* Split naming precedence: object name, prefix, number */
1055: splitname = ((PetscObject)dm)->name;
1056: if (!splitname) {
1057: PetscObjectGetOptionsPrefix((PetscObject)dms[i], &splitname);
1058: if (splitname) {
1059: size_t len;
1060: PetscStrncpy(buf, splitname, sizeof(buf));
1061: buf[sizeof(buf) - 1] = 0;
1062: PetscStrlen(buf, &len);
1063: if (buf[len - 1] == '_') buf[len - 1] = 0; /* Remove trailing underscore if it was used */
1064: splitname = buf;
1065: }
1066: }
1067: if (!splitname) {
1068: PetscSNPrintf(buf, sizeof(buf), "%" PetscInt_FMT, i);
1069: splitname = buf;
1070: }
1071: PetscStrallocpy(splitname, &(*fieldNames)[i]);
1072: }
1073: PetscFree(dms);
1074: }
1075: return 0;
1076: }
1078: /*
1079: This could take over from DMCreateFieldIS(), as it is more general,
1080: making DMCreateFieldIS() a special case -- calling with dmlist == NULL;
1081: At this point it's probably best to be less intrusive, however.
1082: */
1083: PetscErrorCode DMCreateFieldDecomposition_Composite(DM dm, PetscInt *len, char ***namelist, IS **islist, DM **dmlist)
1084: {
1085: PetscInt nDM;
1086: PetscInt i;
1088: DMCreateFieldIS_Composite(dm, len, namelist, islist);
1089: if (dmlist) {
1090: DMCompositeGetNumberDM(dm, &nDM);
1091: PetscMalloc1(nDM, dmlist);
1092: DMCompositeGetEntriesArray(dm, *dmlist);
1093: for (i = 0; i < nDM; i++) PetscObjectReference((PetscObject)((*dmlist)[i]));
1094: }
1095: return 0;
1096: }
1098: /* -------------------------------------------------------------------------------------*/
1099: /*@C
1100: DMCompositeGetLocalVectors - Gets local vectors for each part of a `DMCOMPOSITE`
1101: Use `DMCompositeRestoreLocalVectors()` to return them.
1103: Not Collective
1105: Input Parameter:
1106: . dm - the `DMCOMPOSITE` object
1108: Output Parameter:
1109: . Vec ... - the individual sequential Vecs
1111: Level: advanced
1113: Fortran Note:
1114: Not available from Fortran
1116: .seealso: `DMCOMPOSITE`, `DM`, `DMDestroy()`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`,
1117: `DMCompositeGather()`, `DMCompositeCreate()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeGetAccess()`,
1118: `DMCompositeRestoreLocalVectors()`, `DMCompositeScatter()`, `DMCompositeGetEntries()`
1119: @*/
1120: PetscErrorCode DMCompositeGetLocalVectors(DM dm, ...)
1121: {
1122: va_list Argp;
1123: struct DMCompositeLink *next;
1124: DM_Composite *com = (DM_Composite *)dm->data;
1125: PetscBool flg;
1128: PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg);
1130: next = com->next;
1131: /* loop over packed objects, handling one at at time */
1132: va_start(Argp, dm);
1133: while (next) {
1134: Vec *vec;
1135: vec = va_arg(Argp, Vec *);
1136: if (vec) DMGetLocalVector(next->dm, vec);
1137: next = next->next;
1138: }
1139: va_end(Argp);
1140: return 0;
1141: }
1143: /*@C
1144: DMCompositeRestoreLocalVectors - Restores local vectors for each part of a `DMCOMPOSITE`
1146: Not Collective
1148: Input Parameter:
1149: . dm - the `DMCOMPOSITE` object
1151: Output Parameter:
1152: . Vec ... - the individual sequential `Vec`
1154: Level: advanced
1156: Fortran Note:
1157: Not available from Fortran
1159: .seealso: `DMCOMPOSITE`, `DM`, `DMDestroy()`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`,
1160: `DMCompositeGather()`, `DMCompositeCreate()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeGetAccess()`,
1161: `DMCompositeGetLocalVectors()`, `DMCompositeScatter()`, `DMCompositeGetEntries()`
1162: @*/
1163: PetscErrorCode DMCompositeRestoreLocalVectors(DM dm, ...)
1164: {
1165: va_list Argp;
1166: struct DMCompositeLink *next;
1167: DM_Composite *com = (DM_Composite *)dm->data;
1168: PetscBool flg;
1171: PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg);
1173: next = com->next;
1174: /* loop over packed objects, handling one at at time */
1175: va_start(Argp, dm);
1176: while (next) {
1177: Vec *vec;
1178: vec = va_arg(Argp, Vec *);
1179: if (vec) DMRestoreLocalVector(next->dm, vec);
1180: next = next->next;
1181: }
1182: va_end(Argp);
1183: return 0;
1184: }
1186: /* -------------------------------------------------------------------------------------*/
1187: /*@C
1188: DMCompositeGetEntries - Gets the `DM` for each entry in a `DMCOMPOSITE`.
1190: Not Collective
1192: Input Parameter:
1193: . dm - the `DMCOMPOSITE` object
1195: Output Parameter:
1196: . DM ... - the individual entries `DM`
1198: Level: advanced
1200: Fortran Note:
1201: Available as `DMCompositeGetEntries()` for one output `DM`, DMCompositeGetEntries2() for 2, etc
1203: .seealso: `DMCOMPOSITE`, `DM`, `DMDestroy()`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`, `DMCompositeGetEntriesArray()`
1204: `DMCompositeGather()`, `DMCompositeCreate()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeGetAccess()`,
1205: `DMCompositeRestoreLocalVectors()`, `DMCompositeGetLocalVectors()`, `DMCompositeScatter()`,
1206: `DMCompositeGetLocalVectors()`, `DMCompositeRestoreLocalVectors()`
1207: @*/
1208: PetscErrorCode DMCompositeGetEntries(DM dm, ...)
1209: {
1210: va_list Argp;
1211: struct DMCompositeLink *next;
1212: DM_Composite *com = (DM_Composite *)dm->data;
1213: PetscBool flg;
1216: PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg);
1218: next = com->next;
1219: /* loop over packed objects, handling one at at time */
1220: va_start(Argp, dm);
1221: while (next) {
1222: DM *dmn;
1223: dmn = va_arg(Argp, DM *);
1224: if (dmn) *dmn = next->dm;
1225: next = next->next;
1226: }
1227: va_end(Argp);
1228: return 0;
1229: }
1231: /*@C
1232: DMCompositeGetEntriesArray - Gets the DM for each entry in a `DMCOMPOSITE`
1234: Not Collective
1236: Input Parameter:
1237: . dm - the `DMCOMPOSITE` object
1239: Output Parameter:
1240: . dms - array of sufficient length (see `DMCompositeGetNumberDM()`) to hold the individual `DM`
1242: Level: advanced
1244: .seealso: `DMCOMPOSITE`, `DM`, `DMDestroy()`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`, `DMCompositeGetEntries()`
1245: `DMCompositeGather()`, `DMCompositeCreate()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeGetAccess()`,
1246: `DMCompositeRestoreLocalVectors()`, `DMCompositeGetLocalVectors()`, `DMCompositeScatter()`,
1247: `DMCompositeGetLocalVectors()`, `DMCompositeRestoreLocalVectors()`
1248: @*/
1249: PetscErrorCode DMCompositeGetEntriesArray(DM dm, DM dms[])
1250: {
1251: struct DMCompositeLink *next;
1252: DM_Composite *com = (DM_Composite *)dm->data;
1253: PetscInt i;
1254: PetscBool flg;
1257: PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg);
1259: /* loop over packed objects, handling one at at time */
1260: for (next = com->next, i = 0; next; next = next->next, i++) dms[i] = next->dm;
1261: return 0;
1262: }
1264: typedef struct {
1265: DM dm;
1266: PetscViewer *subv;
1267: Vec *vecs;
1268: } GLVisViewerCtx;
1270: static PetscErrorCode DestroyGLVisViewerCtx_Private(void *vctx)
1271: {
1272: GLVisViewerCtx *ctx = (GLVisViewerCtx *)vctx;
1273: PetscInt i, n;
1275: DMCompositeGetNumberDM(ctx->dm, &n);
1276: for (i = 0; i < n; i++) PetscViewerDestroy(&ctx->subv[i]);
1277: PetscFree2(ctx->subv, ctx->vecs);
1278: DMDestroy(&ctx->dm);
1279: PetscFree(ctx);
1280: return 0;
1281: }
1283: static PetscErrorCode DMCompositeSampleGLVisFields_Private(PetscObject oX, PetscInt nf, PetscObject oXfield[], void *vctx)
1284: {
1285: Vec X = (Vec)oX;
1286: GLVisViewerCtx *ctx = (GLVisViewerCtx *)vctx;
1287: PetscInt i, n, cumf;
1289: DMCompositeGetNumberDM(ctx->dm, &n);
1290: DMCompositeGetAccessArray(ctx->dm, X, n, NULL, ctx->vecs);
1291: for (i = 0, cumf = 0; i < n; i++) {
1292: PetscErrorCode (*g2l)(PetscObject, PetscInt, PetscObject[], void *);
1293: void *fctx;
1294: PetscInt nfi;
1296: PetscViewerGLVisGetFields_Private(ctx->subv[i], &nfi, NULL, NULL, &g2l, NULL, &fctx);
1297: if (!nfi) continue;
1298: if (g2l) (*g2l)((PetscObject)ctx->vecs[i], nfi, oXfield + cumf, fctx);
1299: else VecCopy(ctx->vecs[i], (Vec)(oXfield[cumf]));
1300: cumf += nfi;
1301: }
1302: DMCompositeRestoreAccessArray(ctx->dm, X, n, NULL, ctx->vecs);
1303: return 0;
1304: }
1306: static PetscErrorCode DMSetUpGLVisViewer_Composite(PetscObject odm, PetscViewer viewer)
1307: {
1308: DM dm = (DM)odm, *dms;
1309: Vec *Ufds;
1310: GLVisViewerCtx *ctx;
1311: PetscInt i, n, tnf, *sdim;
1312: char **fecs;
1314: PetscNew(&ctx);
1315: PetscObjectReference((PetscObject)dm);
1316: ctx->dm = dm;
1317: DMCompositeGetNumberDM(dm, &n);
1318: PetscMalloc1(n, &dms);
1319: DMCompositeGetEntriesArray(dm, dms);
1320: PetscMalloc2(n, &ctx->subv, n, &ctx->vecs);
1321: for (i = 0, tnf = 0; i < n; i++) {
1322: PetscInt nf;
1324: PetscViewerCreate(PetscObjectComm(odm), &ctx->subv[i]);
1325: PetscViewerSetType(ctx->subv[i], PETSCVIEWERGLVIS);
1326: PetscViewerGLVisSetDM_Private(ctx->subv[i], (PetscObject)dms[i]);
1327: PetscViewerGLVisGetFields_Private(ctx->subv[i], &nf, NULL, NULL, NULL, NULL, NULL);
1328: tnf += nf;
1329: }
1330: PetscFree(dms);
1331: PetscMalloc3(tnf, &fecs, tnf, &sdim, tnf, &Ufds);
1332: for (i = 0, tnf = 0; i < n; i++) {
1333: PetscInt *sd, nf, f;
1334: const char **fec;
1335: Vec *Uf;
1337: PetscViewerGLVisGetFields_Private(ctx->subv[i], &nf, &fec, &sd, NULL, (PetscObject **)&Uf, NULL);
1338: for (f = 0; f < nf; f++) {
1339: PetscStrallocpy(fec[f], &fecs[tnf + f]);
1340: Ufds[tnf + f] = Uf[f];
1341: sdim[tnf + f] = sd[f];
1342: }
1343: tnf += nf;
1344: }
1345: PetscViewerGLVisSetFields(viewer, tnf, (const char **)fecs, sdim, DMCompositeSampleGLVisFields_Private, (PetscObject *)Ufds, ctx, DestroyGLVisViewerCtx_Private);
1346: for (i = 0; i < tnf; i++) PetscFree(fecs[i]);
1347: PetscFree3(fecs, sdim, Ufds);
1348: return 0;
1349: }
1351: PetscErrorCode DMRefine_Composite(DM dmi, MPI_Comm comm, DM *fine)
1352: {
1353: struct DMCompositeLink *next;
1354: DM_Composite *com = (DM_Composite *)dmi->data;
1355: DM dm;
1358: if (comm == MPI_COMM_NULL) PetscObjectGetComm((PetscObject)dmi, &comm);
1359: DMSetUp(dmi);
1360: next = com->next;
1361: DMCompositeCreate(comm, fine);
1363: /* loop over packed objects, handling one at at time */
1364: while (next) {
1365: DMRefine(next->dm, comm, &dm);
1366: DMCompositeAddDM(*fine, dm);
1367: PetscObjectDereference((PetscObject)dm);
1368: next = next->next;
1369: }
1370: return 0;
1371: }
1373: PetscErrorCode DMCoarsen_Composite(DM dmi, MPI_Comm comm, DM *fine)
1374: {
1375: struct DMCompositeLink *next;
1376: DM_Composite *com = (DM_Composite *)dmi->data;
1377: DM dm;
1380: DMSetUp(dmi);
1381: if (comm == MPI_COMM_NULL) PetscObjectGetComm((PetscObject)dmi, &comm);
1382: next = com->next;
1383: DMCompositeCreate(comm, fine);
1385: /* loop over packed objects, handling one at at time */
1386: while (next) {
1387: DMCoarsen(next->dm, comm, &dm);
1388: DMCompositeAddDM(*fine, dm);
1389: PetscObjectDereference((PetscObject)dm);
1390: next = next->next;
1391: }
1392: return 0;
1393: }
1395: PetscErrorCode DMCreateInterpolation_Composite(DM coarse, DM fine, Mat *A, Vec *v)
1396: {
1397: PetscInt m, n, M, N, nDM, i;
1398: struct DMCompositeLink *nextc;
1399: struct DMCompositeLink *nextf;
1400: Vec gcoarse, gfine, *vecs;
1401: DM_Composite *comcoarse = (DM_Composite *)coarse->data;
1402: DM_Composite *comfine = (DM_Composite *)fine->data;
1403: Mat *mats;
1407: DMSetUp(coarse);
1408: DMSetUp(fine);
1409: /* use global vectors only for determining matrix layout */
1410: DMGetGlobalVector(coarse, &gcoarse);
1411: DMGetGlobalVector(fine, &gfine);
1412: VecGetLocalSize(gcoarse, &n);
1413: VecGetLocalSize(gfine, &m);
1414: VecGetSize(gcoarse, &N);
1415: VecGetSize(gfine, &M);
1416: DMRestoreGlobalVector(coarse, &gcoarse);
1417: DMRestoreGlobalVector(fine, &gfine);
1419: nDM = comfine->nDM;
1421: PetscCalloc1(nDM * nDM, &mats);
1422: if (v) PetscCalloc1(nDM, &vecs);
1424: /* loop over packed objects, handling one at at time */
1425: for (nextc = comcoarse->next, nextf = comfine->next, i = 0; nextc; nextc = nextc->next, nextf = nextf->next, i++) {
1426: if (!v) DMCreateInterpolation(nextc->dm, nextf->dm, &mats[i * nDM + i], NULL);
1427: else DMCreateInterpolation(nextc->dm, nextf->dm, &mats[i * nDM + i], &vecs[i]);
1428: }
1429: MatCreateNest(PetscObjectComm((PetscObject)fine), nDM, NULL, nDM, NULL, mats, A);
1430: if (v) VecCreateNest(PetscObjectComm((PetscObject)fine), nDM, NULL, vecs, v);
1431: for (i = 0; i < nDM * nDM; i++) MatDestroy(&mats[i]);
1432: PetscFree(mats);
1433: if (v) {
1434: for (i = 0; i < nDM; i++) VecDestroy(&vecs[i]);
1435: PetscFree(vecs);
1436: }
1437: return 0;
1438: }
1440: static PetscErrorCode DMGetLocalToGlobalMapping_Composite(DM dm)
1441: {
1442: DM_Composite *com = (DM_Composite *)dm->data;
1443: ISLocalToGlobalMapping *ltogs;
1444: PetscInt i;
1446: /* Set the ISLocalToGlobalMapping on the new matrix */
1447: DMCompositeGetISLocalToGlobalMappings(dm, <ogs);
1448: ISLocalToGlobalMappingConcatenate(PetscObjectComm((PetscObject)dm), com->nDM, ltogs, &dm->ltogmap);
1449: for (i = 0; i < com->nDM; i++) ISLocalToGlobalMappingDestroy(<ogs[i]);
1450: PetscFree(ltogs);
1451: return 0;
1452: }
1454: PetscErrorCode DMCreateColoring_Composite(DM dm, ISColoringType ctype, ISColoring *coloring)
1455: {
1456: PetscInt n, i, cnt;
1457: ISColoringValue *colors;
1458: PetscBool dense = PETSC_FALSE;
1459: ISColoringValue maxcol = 0;
1460: DM_Composite *com = (DM_Composite *)dm->data;
1464: if (ctype == IS_COLORING_GLOBAL) {
1465: n = com->n;
1466: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Unknown ISColoringType");
1467: PetscMalloc1(n, &colors); /* freed in ISColoringDestroy() */
1469: PetscOptionsGetBool(((PetscObject)dm)->options, ((PetscObject)dm)->prefix, "-dmcomposite_dense_jacobian", &dense, NULL);
1470: if (dense) {
1471: for (i = 0; i < n; i++) colors[i] = (ISColoringValue)(com->rstart + i);
1472: maxcol = com->N;
1473: } else {
1474: struct DMCompositeLink *next = com->next;
1475: PetscMPIInt rank;
1477: MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
1478: cnt = 0;
1479: while (next) {
1480: ISColoring lcoloring;
1482: DMCreateColoring(next->dm, IS_COLORING_GLOBAL, &lcoloring);
1483: for (i = 0; i < lcoloring->N; i++) colors[cnt++] = maxcol + lcoloring->colors[i];
1484: maxcol += lcoloring->n;
1485: ISColoringDestroy(&lcoloring);
1486: next = next->next;
1487: }
1488: }
1489: ISColoringCreate(PetscObjectComm((PetscObject)dm), maxcol, n, colors, PETSC_OWN_POINTER, coloring);
1490: return 0;
1491: }
1493: PetscErrorCode DMGlobalToLocalBegin_Composite(DM dm, Vec gvec, InsertMode mode, Vec lvec)
1494: {
1495: struct DMCompositeLink *next;
1496: PetscScalar *garray, *larray;
1497: DM_Composite *com = (DM_Composite *)dm->data;
1502: if (!com->setup) DMSetUp(dm);
1504: VecGetArray(gvec, &garray);
1505: VecGetArray(lvec, &larray);
1507: /* loop over packed objects, handling one at at time */
1508: next = com->next;
1509: while (next) {
1510: Vec local, global;
1511: PetscInt N;
1513: DMGetGlobalVector(next->dm, &global);
1514: VecGetLocalSize(global, &N);
1515: VecPlaceArray(global, garray);
1516: DMGetLocalVector(next->dm, &local);
1517: VecPlaceArray(local, larray);
1518: DMGlobalToLocalBegin(next->dm, global, mode, local);
1519: DMGlobalToLocalEnd(next->dm, global, mode, local);
1520: VecResetArray(global);
1521: VecResetArray(local);
1522: DMRestoreGlobalVector(next->dm, &global);
1523: DMRestoreLocalVector(next->dm, &local);
1525: larray += next->nlocal;
1526: garray += next->n;
1527: next = next->next;
1528: }
1530: VecRestoreArray(gvec, NULL);
1531: VecRestoreArray(lvec, NULL);
1532: return 0;
1533: }
1535: PetscErrorCode DMGlobalToLocalEnd_Composite(DM dm, Vec gvec, InsertMode mode, Vec lvec)
1536: {
1540: return 0;
1541: }
1543: PetscErrorCode DMLocalToGlobalBegin_Composite(DM dm, Vec lvec, InsertMode mode, Vec gvec)
1544: {
1545: struct DMCompositeLink *next;
1546: PetscScalar *larray, *garray;
1547: DM_Composite *com = (DM_Composite *)dm->data;
1553: if (!com->setup) DMSetUp(dm);
1555: VecGetArray(lvec, &larray);
1556: VecGetArray(gvec, &garray);
1558: /* loop over packed objects, handling one at at time */
1559: next = com->next;
1560: while (next) {
1561: Vec global, local;
1563: DMGetLocalVector(next->dm, &local);
1564: VecPlaceArray(local, larray);
1565: DMGetGlobalVector(next->dm, &global);
1566: VecPlaceArray(global, garray);
1567: DMLocalToGlobalBegin(next->dm, local, mode, global);
1568: DMLocalToGlobalEnd(next->dm, local, mode, global);
1569: VecResetArray(local);
1570: VecResetArray(global);
1571: DMRestoreGlobalVector(next->dm, &global);
1572: DMRestoreLocalVector(next->dm, &local);
1574: garray += next->n;
1575: larray += next->nlocal;
1576: next = next->next;
1577: }
1579: VecRestoreArray(gvec, NULL);
1580: VecRestoreArray(lvec, NULL);
1582: return 0;
1583: }
1585: PetscErrorCode DMLocalToGlobalEnd_Composite(DM dm, Vec lvec, InsertMode mode, Vec gvec)
1586: {
1590: return 0;
1591: }
1593: PetscErrorCode DMLocalToLocalBegin_Composite(DM dm, Vec vec1, InsertMode mode, Vec vec2)
1594: {
1595: struct DMCompositeLink *next;
1596: PetscScalar *array1, *array2;
1597: DM_Composite *com = (DM_Composite *)dm->data;
1603: if (!com->setup) DMSetUp(dm);
1605: VecGetArray(vec1, &array1);
1606: VecGetArray(vec2, &array2);
1608: /* loop over packed objects, handling one at at time */
1609: next = com->next;
1610: while (next) {
1611: Vec local1, local2;
1613: DMGetLocalVector(next->dm, &local1);
1614: VecPlaceArray(local1, array1);
1615: DMGetLocalVector(next->dm, &local2);
1616: VecPlaceArray(local2, array2);
1617: DMLocalToLocalBegin(next->dm, local1, mode, local2);
1618: DMLocalToLocalEnd(next->dm, local1, mode, local2);
1619: VecResetArray(local2);
1620: DMRestoreLocalVector(next->dm, &local2);
1621: VecResetArray(local1);
1622: DMRestoreLocalVector(next->dm, &local1);
1624: array1 += next->nlocal;
1625: array2 += next->nlocal;
1626: next = next->next;
1627: }
1629: VecRestoreArray(vec1, NULL);
1630: VecRestoreArray(vec2, NULL);
1632: return 0;
1633: }
1635: PetscErrorCode DMLocalToLocalEnd_Composite(DM dm, Vec lvec, InsertMode mode, Vec gvec)
1636: {
1640: return 0;
1641: }
1643: /*MC
1644: DMCOMPOSITE = "composite" - A `DM` object that is used to manage data for a collection of `DM`
1646: Level: intermediate
1648: .seealso: `DMType`, `DM`, `DMDACreate()`, `DMCreate()`, `DMSetType()`, `DMCompositeCreate()`
1649: M*/
1651: PETSC_EXTERN PetscErrorCode DMCreate_Composite(DM p)
1652: {
1653: DM_Composite *com;
1655: PetscNew(&com);
1656: p->data = com;
1657: com->n = 0;
1658: com->nghost = 0;
1659: com->next = NULL;
1660: com->nDM = 0;
1662: p->ops->createglobalvector = DMCreateGlobalVector_Composite;
1663: p->ops->createlocalvector = DMCreateLocalVector_Composite;
1664: p->ops->getlocaltoglobalmapping = DMGetLocalToGlobalMapping_Composite;
1665: p->ops->createfieldis = DMCreateFieldIS_Composite;
1666: p->ops->createfielddecomposition = DMCreateFieldDecomposition_Composite;
1667: p->ops->refine = DMRefine_Composite;
1668: p->ops->coarsen = DMCoarsen_Composite;
1669: p->ops->createinterpolation = DMCreateInterpolation_Composite;
1670: p->ops->creatematrix = DMCreateMatrix_Composite;
1671: p->ops->getcoloring = DMCreateColoring_Composite;
1672: p->ops->globaltolocalbegin = DMGlobalToLocalBegin_Composite;
1673: p->ops->globaltolocalend = DMGlobalToLocalEnd_Composite;
1674: p->ops->localtoglobalbegin = DMLocalToGlobalBegin_Composite;
1675: p->ops->localtoglobalend = DMLocalToGlobalEnd_Composite;
1676: p->ops->localtolocalbegin = DMLocalToLocalBegin_Composite;
1677: p->ops->localtolocalend = DMLocalToLocalEnd_Composite;
1678: p->ops->destroy = DMDestroy_Composite;
1679: p->ops->view = DMView_Composite;
1680: p->ops->setup = DMSetUp_Composite;
1682: PetscObjectComposeFunction((PetscObject)p, "DMSetUpGLVisViewer_C", DMSetUpGLVisViewer_Composite);
1683: return 0;
1684: }
1686: /*@
1687: DMCompositeCreate - Creates a `DMCOMPOSITE`, used to generate "composite"
1688: vectors made up of several subvectors.
1690: Collective
1692: Input Parameter:
1693: . comm - the processors that will share the global vector
1695: Output Parameters:
1696: . packer - the `DMCOMPOSITE` object
1698: Level: advanced
1700: .seealso: `DMCOMPOSITE`, `DM`, `DMDestroy()`, `DMCompositeAddDM()`, `DMCompositeScatter()`, `DMCOMPOSITE`, `DMCreate()`
1701: `DMCompositeGather()`, `DMCreateGlobalVector()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeGetAccess()`
1702: `DMCompositeGetLocalVectors()`, `DMCompositeRestoreLocalVectors()`, `DMCompositeGetEntries()`
1703: @*/
1704: PetscErrorCode DMCompositeCreate(MPI_Comm comm, DM *packer)
1705: {
1707: DMCreate(comm, packer);
1708: DMSetType(*packer, DMCOMPOSITE);
1709: return 0;
1710: }