Actual source code: pack.c
1: #include <../src/dm/impls/composite/packimpl.h>
2: #include <petsc/private/isimpl.h>
3: #include <petsc/private/glvisviewerimpl.h>
4: #include <petscds.h>
6: /*@C
7: DMCompositeSetCoupling - Sets user provided routines that compute the coupling between the
8: separate components `DM` in a `DMCOMPOSITE` to build the correct matrix nonzero structure.
10: Logically Collective; No Fortran Support
12: Input Parameters:
13: + dm - the composite object
14: - FormCoupleLocations - routine to set the nonzero locations in the matrix
16: Level: advanced
18: Note:
19: See `DMSetApplicationContext()` and `DMGetApplicationContext()` for how to get user information into
20: this routine
22: .seealso: `DMCOMPOSITE`, `DM`
23: @*/
24: PetscErrorCode DMCompositeSetCoupling(DM dm, PetscErrorCode (*FormCoupleLocations)(DM, Mat, PetscInt *, PetscInt *, PetscInt, PetscInt, PetscInt, PetscInt))
25: {
26: DM_Composite *com = (DM_Composite *)dm->data;
27: PetscBool flg;
29: PetscFunctionBegin;
30: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
31: PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
32: com->FormCoupleLocations = FormCoupleLocations;
33: PetscFunctionReturn(PETSC_SUCCESS);
34: }
36: static PetscErrorCode DMDestroy_Composite(DM dm)
37: {
38: struct DMCompositeLink *next, *prev;
39: DM_Composite *com = (DM_Composite *)dm->data;
41: PetscFunctionBegin;
42: next = com->next;
43: while (next) {
44: prev = next;
45: next = next->next;
46: PetscCall(DMDestroy(&prev->dm));
47: PetscCall(PetscFree(prev->grstarts));
48: PetscCall(PetscFree(prev));
49: }
50: PetscCall(PetscObjectComposeFunction((PetscObject)dm, "DMSetUpGLVisViewer_C", NULL));
51: /* This was originally freed in DMDestroy(), but that prevents reference counting of backend objects */
52: PetscCall(PetscFree(com));
53: PetscFunctionReturn(PETSC_SUCCESS);
54: }
56: static PetscErrorCode DMView_Composite(DM dm, PetscViewer v)
57: {
58: PetscBool iascii;
59: DM_Composite *com = (DM_Composite *)dm->data;
61: PetscFunctionBegin;
62: PetscCall(PetscObjectTypeCompare((PetscObject)v, PETSCVIEWERASCII, &iascii));
63: if (iascii) {
64: struct DMCompositeLink *lnk = com->next;
65: PetscInt i;
67: PetscCall(PetscViewerASCIIPrintf(v, "DM (%s)\n", ((PetscObject)dm)->prefix ? ((PetscObject)dm)->prefix : "no prefix"));
68: PetscCall(PetscViewerASCIIPrintf(v, " contains %" PetscInt_FMT " DMs\n", com->nDM));
69: PetscCall(PetscViewerASCIIPushTab(v));
70: for (i = 0; lnk; lnk = lnk->next, i++) {
71: PetscCall(PetscViewerASCIIPrintf(v, "Link %" PetscInt_FMT ": DM of type %s\n", i, ((PetscObject)lnk->dm)->type_name));
72: PetscCall(PetscViewerASCIIPushTab(v));
73: PetscCall(DMView(lnk->dm, v));
74: PetscCall(PetscViewerASCIIPopTab(v));
75: }
76: PetscCall(PetscViewerASCIIPopTab(v));
77: }
78: PetscFunctionReturn(PETSC_SUCCESS);
79: }
81: /* --------------------------------------------------------------------------------------*/
82: static PetscErrorCode DMSetUp_Composite(DM dm)
83: {
84: PetscInt nprev = 0;
85: PetscMPIInt rank, size;
86: DM_Composite *com = (DM_Composite *)dm->data;
87: struct DMCompositeLink *next = com->next;
88: PetscLayout map;
90: PetscFunctionBegin;
91: PetscCheck(!com->setup, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "Packer has already been setup");
92: PetscCall(PetscLayoutCreate(PetscObjectComm((PetscObject)dm), &map));
93: PetscCall(PetscLayoutSetLocalSize(map, com->n));
94: PetscCall(PetscLayoutSetSize(map, PETSC_DETERMINE));
95: PetscCall(PetscLayoutSetBlockSize(map, 1));
96: PetscCall(PetscLayoutSetUp(map));
97: PetscCall(PetscLayoutGetSize(map, &com->N));
98: PetscCall(PetscLayoutGetRange(map, &com->rstart, NULL));
99: PetscCall(PetscLayoutDestroy(&map));
101: /* now set the rstart for each linked vector */
102: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
103: PetscCallMPI(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: PetscCall(PetscMalloc1(size, &next->grstarts));
109: PetscCallMPI(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: PetscFunctionReturn(PETSC_SUCCESS);
114: }
116: /*@
117: DMCompositeGetNumberDM - Gets the number of `DM` objects in the `DMCOMPOSITE`
118: representation.
120: Not Collective
122: Input Parameter:
123: . dm - the `DMCOMPOSITE` object
125: Output Parameter:
126: . nDM - the number of `DM`
128: Level: beginner
130: .seealso: `DMCOMPOSITE`, `DM`
131: @*/
132: PetscErrorCode DMCompositeGetNumberDM(DM dm, PetscInt *nDM)
133: {
134: DM_Composite *com = (DM_Composite *)dm->data;
135: PetscBool flg;
137: PetscFunctionBegin;
139: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
140: PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
141: *nDM = com->nDM;
142: PetscFunctionReturn(PETSC_SUCCESS);
143: }
145: /*@C
146: DMCompositeGetAccess - Allows one to access the individual packed vectors in their global
147: representation.
149: Collective
151: Input Parameters:
152: + dm - the `DMCOMPOSITE` object
153: - gvec - the global vector
155: Output Parameter:
156: . ... - the packed parallel vectors, `NULL` for those that are not needed
158: Level: advanced
160: Note:
161: Use `DMCompositeRestoreAccess()` to return the vectors when you no longer need them
163: .seealso: `DMCOMPOSITE`, `DM`, `DMCompositeGetEntries()`, `DMCompositeScatter()`
164: @*/
165: PetscErrorCode DMCompositeGetAccess(DM dm, Vec gvec, ...)
166: {
167: va_list Argp;
168: struct DMCompositeLink *next;
169: DM_Composite *com = (DM_Composite *)dm->data;
170: PetscInt readonly;
171: PetscBool flg;
173: PetscFunctionBegin;
176: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
177: PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
178: next = com->next;
179: if (!com->setup) PetscCall(DMSetUp(dm));
181: PetscCall(VecLockGet(gvec, &readonly));
182: /* loop over packed objects, handling one at a time */
183: va_start(Argp, gvec);
184: while (next) {
185: Vec *vec;
186: vec = va_arg(Argp, Vec *);
187: if (vec) {
188: PetscCall(DMGetGlobalVector(next->dm, vec));
189: if (readonly) {
190: const PetscScalar *array;
191: PetscCall(VecGetArrayRead(gvec, &array));
192: PetscCall(VecPlaceArray(*vec, array + next->rstart));
193: PetscCall(VecLockReadPush(*vec));
194: PetscCall(VecRestoreArrayRead(gvec, &array));
195: } else {
196: PetscScalar *array;
197: PetscCall(VecGetArray(gvec, &array));
198: PetscCall(VecPlaceArray(*vec, array + next->rstart));
199: PetscCall(VecRestoreArray(gvec, &array));
200: }
201: }
202: next = next->next;
203: }
204: va_end(Argp);
205: PetscFunctionReturn(PETSC_SUCCESS);
206: }
208: /*@
209: DMCompositeGetAccessArray - Allows one to access the individual packed vectors in their global
210: representation.
212: Collective
214: Input Parameters:
215: + dm - the `DMCOMPOSITE`
216: . pvec - packed vector
217: . nwanted - number of vectors wanted
218: - wanted - sorted array of integers indicating thde vectors wanted, or `NULL` to get all vectors, length `nwanted`
220: Output Parameter:
221: . vecs - array of requested global vectors (must be previously allocated and of length `nwanted`)
223: Level: advanced
225: Note:
226: Use `DMCompositeRestoreAccessArray()` to return the vectors when you no longer need them
228: .seealso: `DMCOMPOSITE`, `DM`, `DMCompositeGetAccess()`, `DMCompositeGetEntries()`, `DMCompositeScatter()`, `DMCompositeGather()`
229: @*/
230: PetscErrorCode DMCompositeGetAccessArray(DM dm, Vec pvec, PetscInt nwanted, const PetscInt wanted[], Vec vecs[])
231: {
232: struct DMCompositeLink *link;
233: PetscInt i, wnum;
234: DM_Composite *com = (DM_Composite *)dm->data;
235: PetscInt readonly;
236: PetscBool flg;
238: PetscFunctionBegin;
241: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
242: PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
243: if (!com->setup) PetscCall(DMSetUp(dm));
245: PetscCall(VecLockGet(pvec, &readonly));
246: for (i = 0, wnum = 0, link = com->next; link && wnum < nwanted; i++, link = link->next) {
247: if (!wanted || i == wanted[wnum]) {
248: Vec v;
249: PetscCall(DMGetGlobalVector(link->dm, &v));
250: if (readonly) {
251: const PetscScalar *array;
252: PetscCall(VecGetArrayRead(pvec, &array));
253: PetscCall(VecPlaceArray(v, array + link->rstart));
254: PetscCall(VecLockReadPush(v));
255: PetscCall(VecRestoreArrayRead(pvec, &array));
256: } else {
257: PetscScalar *array;
258: PetscCall(VecGetArray(pvec, &array));
259: PetscCall(VecPlaceArray(v, array + link->rstart));
260: PetscCall(VecRestoreArray(pvec, &array));
261: }
262: vecs[wnum++] = v;
263: }
264: }
265: PetscFunctionReturn(PETSC_SUCCESS);
266: }
268: /*@
269: DMCompositeGetLocalAccessArray - Allows one to access the individual
270: packed vectors in their local representation.
272: Collective
274: Input Parameters:
275: + dm - the `DMCOMPOSITE`
276: . pvec - packed vector
277: . nwanted - number of vectors wanted
278: - wanted - sorted array of vectors wanted, or `NULL` to get all vectors, length `nwanted`
280: Output Parameter:
281: . vecs - array of requested local vectors (must be allocated and of length `nwanted`)
283: Level: advanced
285: Note:
286: Use `DMCompositeRestoreLocalAccessArray()` to return the vectors
287: when you no longer need them.
289: .seealso: `DMCOMPOSITE`, `DM`, `DMCompositeRestoreLocalAccessArray()`, `DMCompositeGetAccess()`,
290: `DMCompositeGetEntries()`, `DMCompositeScatter()`, `DMCompositeGather()`
291: @*/
292: PetscErrorCode DMCompositeGetLocalAccessArray(DM dm, Vec pvec, PetscInt nwanted, const PetscInt wanted[], Vec vecs[])
293: {
294: struct DMCompositeLink *link;
295: PetscInt i, wnum;
296: DM_Composite *com = (DM_Composite *)dm->data;
297: PetscInt readonly;
298: PetscInt nlocal = 0;
299: PetscBool flg;
301: PetscFunctionBegin;
304: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
305: PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
306: if (!com->setup) PetscCall(DMSetUp(dm));
308: PetscCall(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: PetscCall(DMGetLocalVector(link->dm, &v));
313: if (readonly) {
314: const PetscScalar *array;
315: PetscCall(VecGetArrayRead(pvec, &array));
316: PetscCall(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: PetscCall(VecLockReadPush(v));
319: PetscCall(VecRestoreArrayRead(pvec, &array));
320: } else {
321: PetscScalar *array;
322: PetscCall(VecGetArray(pvec, &array));
323: PetscCall(VecPlaceArray(v, array + nlocal));
324: PetscCall(VecRestoreArray(pvec, &array));
325: }
326: vecs[wnum++] = v;
327: }
329: nlocal += link->nlocal;
330: }
331: PetscFunctionReturn(PETSC_SUCCESS);
332: }
334: /*@C
335: DMCompositeRestoreAccess - Returns the vectors obtained with `DMCompositeGetAccess()`
336: representation.
338: Collective
340: Input Parameters:
341: + dm - the `DMCOMPOSITE` object
342: . gvec - the global vector
343: - ... - the individual parallel vectors, `NULL` for those that are not needed
345: Level: advanced
347: .seealso: `DMCOMPOSITE`, `DM`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`,
348: `DMCompositeGather()`, `DMCompositeCreate()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeScatter()`,
349: `DMCompositeGetAccess()`
350: @*/
351: PetscErrorCode DMCompositeRestoreAccess(DM dm, Vec gvec, ...)
352: {
353: va_list Argp;
354: struct DMCompositeLink *next;
355: DM_Composite *com = (DM_Composite *)dm->data;
356: PetscInt readonly;
357: PetscBool flg;
359: PetscFunctionBegin;
362: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
363: PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
364: next = com->next;
365: if (!com->setup) PetscCall(DMSetUp(dm));
367: PetscCall(VecLockGet(gvec, &readonly));
368: /* loop over packed objects, handling one at a time */
369: va_start(Argp, gvec);
370: while (next) {
371: Vec *vec;
372: vec = va_arg(Argp, Vec *);
373: if (vec) {
374: PetscCall(VecResetArray(*vec));
375: if (readonly) PetscCall(VecLockReadPop(*vec));
376: PetscCall(DMRestoreGlobalVector(next->dm, vec));
377: }
378: next = next->next;
379: }
380: va_end(Argp);
381: PetscFunctionReturn(PETSC_SUCCESS);
382: }
384: /*@
385: DMCompositeRestoreAccessArray - Returns the vectors obtained with `DMCompositeGetAccessArray()`
387: Collective
389: Input Parameters:
390: + dm - the `DMCOMPOSITE` object
391: . pvec - packed vector
392: . nwanted - number of vectors wanted
393: . wanted - sorted array of vectors wanted, or `NULL` to restore all vectors
394: - vecs - array of global vectors
396: Level: advanced
398: .seealso: `DMCOMPOSITE`, `DM`, `DMCompositeRestoreAccess()`, `DMCompositeRestoreEntries()`, `DMCompositeScatter()`, `DMCompositeGather()`
399: @*/
400: PetscErrorCode DMCompositeRestoreAccessArray(DM dm, Vec pvec, PetscInt nwanted, const PetscInt wanted[], Vec vecs[])
401: {
402: struct DMCompositeLink *link;
403: PetscInt i, wnum;
404: DM_Composite *com = (DM_Composite *)dm->data;
405: PetscInt readonly;
406: PetscBool flg;
408: PetscFunctionBegin;
411: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
412: PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
413: if (!com->setup) PetscCall(DMSetUp(dm));
415: PetscCall(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: PetscCall(VecResetArray(vecs[wnum]));
419: if (readonly) PetscCall(VecLockReadPop(vecs[wnum]));
420: PetscCall(DMRestoreGlobalVector(link->dm, &vecs[wnum]));
421: wnum++;
422: }
423: }
424: PetscFunctionReturn(PETSC_SUCCESS);
425: }
427: /*@
428: DMCompositeRestoreLocalAccessArray - Returns the vectors obtained with `DMCompositeGetLocalAccessArray()`.
430: Collective
432: Input Parameters:
433: + dm - the `DMCOMPOSITE` 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
439: Level: advanced
441: Note:
442: `nwanted` and `wanted` must match the values given to `DMCompositeGetLocalAccessArray()`
443: otherwise the call will fail.
445: .seealso: `DMCOMPOSITE`, `DM`, `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;
457: PetscFunctionBegin;
460: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
461: PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
462: if (!com->setup) PetscCall(DMSetUp(dm));
464: PetscCall(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: PetscCall(VecResetArray(vecs[wnum]));
468: if (readonly) PetscCall(VecLockReadPop(vecs[wnum]));
469: PetscCall(DMRestoreLocalVector(link->dm, &vecs[wnum]));
470: wnum++;
471: }
472: }
473: PetscFunctionReturn(PETSC_SUCCESS);
474: }
476: /*@C
477: DMCompositeScatter - Scatters from a global packed vector into its individual local vectors
479: Collective
481: Input Parameters:
482: + dm - the `DMCOMPOSITE` object
483: . gvec - the global vector
484: - ... - 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;
505: PetscFunctionBegin;
508: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
509: PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
510: if (!com->setup) PetscCall(DMSetUp(dm));
512: /* loop over packed objects, handling one at a time */
513: va_start(Argp, gvec);
514: for (cnt = 3, next = com->next; next; cnt++, next = next->next) {
515: Vec local;
516: local = va_arg(Argp, Vec);
517: if (local) {
518: Vec global;
519: const PetscScalar *array;
521: PetscCall(DMGetGlobalVector(next->dm, &global));
522: PetscCall(VecGetArrayRead(gvec, &array));
523: PetscCall(VecPlaceArray(global, array + next->rstart));
524: PetscCall(DMGlobalToLocalBegin(next->dm, global, INSERT_VALUES, local));
525: PetscCall(DMGlobalToLocalEnd(next->dm, global, INSERT_VALUES, local));
526: PetscCall(VecRestoreArrayRead(gvec, &array));
527: PetscCall(VecResetArray(global));
528: PetscCall(DMRestoreGlobalVector(next->dm, &global));
529: }
530: }
531: va_end(Argp);
532: PetscFunctionReturn(PETSC_SUCCESS);
533: }
535: /*@
536: DMCompositeScatterArray - Scatters from a global packed vector into its individual local vectors
538: Collective
540: Input Parameters:
541: + dm - the `DMCOMPOSITE` object
542: . gvec - the global vector
543: - lvecs - array of local vectors, NULL for any that are not needed
545: Level: advanced
547: Note:
548: This is a non-variadic alternative to `DMCompositeScatter()`
550: .seealso: `DMCOMPOSITE`, `DM`, `DMDestroy()`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`
551: `DMCompositeGather()`, `DMCompositeCreate()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeGetAccess()`,
552: `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;
561: PetscFunctionBegin;
564: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
565: PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
566: if (!com->setup) PetscCall(DMSetUp(dm));
568: /* loop over packed objects, handling one at a time */
569: for (i = 0, next = com->next; next; next = next->next, i++) {
570: if (lvecs[i]) {
571: Vec global;
572: const PetscScalar *array;
574: PetscCall(DMGetGlobalVector(next->dm, &global));
575: PetscCall(VecGetArrayRead(gvec, &array));
576: PetscCall(VecPlaceArray(global, (PetscScalar *)array + next->rstart));
577: PetscCall(DMGlobalToLocalBegin(next->dm, global, INSERT_VALUES, lvecs[i]));
578: PetscCall(DMGlobalToLocalEnd(next->dm, global, INSERT_VALUES, lvecs[i]));
579: PetscCall(VecRestoreArrayRead(gvec, &array));
580: PetscCall(VecResetArray(global));
581: PetscCall(DMRestoreGlobalVector(next->dm, &global));
582: }
583: }
584: PetscFunctionReturn(PETSC_SUCCESS);
585: }
587: /*@C
588: DMCompositeGather - Gathers into a global packed vector from its individual local vectors
590: Collective
592: Input Parameters:
593: + dm - the `DMCOMPOSITE` object
594: . imode - `INSERT_VALUES` or `ADD_VALUES`
595: . gvec - the global vector
596: - ... - the individual sequential vectors, `NULL` for any that are not needed
598: Level: advanced
600: Fortran Notes:
601: Fortran users should use `DMCompositeGatherArray()`
603: .seealso: `DMCOMPOSITE`, `DM`, `DMDestroy()`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`,
604: `DMCompositeScatter()`, `DMCompositeCreate()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeGetAccess()`,
605: `DMCompositeGetLocalVectors()`, `DMCompositeRestoreLocalVectors()`, `DMCompositeGetEntries()`
606: @*/
607: PetscErrorCode DMCompositeGather(DM dm, InsertMode imode, Vec gvec, ...)
608: {
609: va_list Argp;
610: struct DMCompositeLink *next;
611: DM_Composite *com = (DM_Composite *)dm->data;
612: PETSC_UNUSED PetscInt cnt;
613: PetscBool flg;
615: PetscFunctionBegin;
618: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
619: PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
620: if (!com->setup) PetscCall(DMSetUp(dm));
622: /* loop over packed objects, handling one at a time */
623: va_start(Argp, gvec);
624: for (cnt = 3, next = com->next; next; cnt++, next = next->next) {
625: Vec local;
626: local = va_arg(Argp, Vec);
627: if (local) {
628: PetscScalar *array;
629: Vec global;
631: PetscCall(DMGetGlobalVector(next->dm, &global));
632: PetscCall(VecGetArray(gvec, &array));
633: PetscCall(VecPlaceArray(global, array + next->rstart));
634: PetscCall(DMLocalToGlobalBegin(next->dm, local, imode, global));
635: PetscCall(DMLocalToGlobalEnd(next->dm, local, imode, global));
636: PetscCall(VecRestoreArray(gvec, &array));
637: PetscCall(VecResetArray(global));
638: PetscCall(DMRestoreGlobalVector(next->dm, &global));
639: }
640: }
641: va_end(Argp);
642: PetscFunctionReturn(PETSC_SUCCESS);
643: }
645: /*@
646: DMCompositeGatherArray - Gathers into a global packed vector from its individual local vectors
648: Collective
650: Input Parameters:
651: + dm - the `DMCOMPOSITE` object
652: . gvec - the global vector
653: . imode - `INSERT_VALUES` or `ADD_VALUES`
654: - lvecs - the individual sequential vectors, NULL for any that are not needed
656: Level: advanced
658: Note:
659: This is a non-variadic alternative to `DMCompositeGather()`.
661: .seealso: `DMCOMPOSITE`, `DM`, `DMDestroy()`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`,
662: `DMCompositeScatter()`, `DMCompositeCreate()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeGetAccess()`,
663: `DMCompositeGetLocalVectors()`, `DMCompositeRestoreLocalVectors()`, `DMCompositeGetEntries()`,
664: @*/
665: PetscErrorCode DMCompositeGatherArray(DM dm, InsertMode imode, Vec gvec, Vec *lvecs)
666: {
667: struct DMCompositeLink *next;
668: DM_Composite *com = (DM_Composite *)dm->data;
669: PetscInt i;
670: PetscBool flg;
672: PetscFunctionBegin;
675: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
676: PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
677: if (!com->setup) PetscCall(DMSetUp(dm));
679: /* loop over packed objects, handling one at a time */
680: for (next = com->next, i = 0; next; next = next->next, i++) {
681: if (lvecs[i]) {
682: PetscScalar *array;
683: Vec global;
685: PetscCall(DMGetGlobalVector(next->dm, &global));
686: PetscCall(VecGetArray(gvec, &array));
687: PetscCall(VecPlaceArray(global, array + next->rstart));
688: PetscCall(DMLocalToGlobalBegin(next->dm, lvecs[i], imode, global));
689: PetscCall(DMLocalToGlobalEnd(next->dm, lvecs[i], imode, global));
690: PetscCall(VecRestoreArray(gvec, &array));
691: PetscCall(VecResetArray(global));
692: PetscCall(DMRestoreGlobalVector(next->dm, &global));
693: }
694: }
695: PetscFunctionReturn(PETSC_SUCCESS);
696: }
698: /*@
699: DMCompositeAddDM - adds a `DM` vector to a `DMCOMPOSITE`
701: Collective
703: Input Parameters:
704: + dmc - the `DMCOMPOSITE` object
705: - dm - the `DM` object
707: Level: advanced
709: .seealso: `DMCOMPOSITE`, `DM`, `DMDestroy()`, `DMCompositeGather()`, `DMCreateGlobalVector()`,
710: `DMCompositeScatter()`, `DMCompositeCreate()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeGetAccess()`,
711: `DMCompositeGetLocalVectors()`, `DMCompositeRestoreLocalVectors()`, `DMCompositeGetEntries()`
712: @*/
713: PetscErrorCode DMCompositeAddDM(DM dmc, DM dm)
714: {
715: PetscInt n, nlocal;
716: struct DMCompositeLink *mine, *next;
717: Vec global, local;
718: DM_Composite *com = (DM_Composite *)dmc->data;
719: PetscBool flg;
721: PetscFunctionBegin;
724: PetscCall(PetscObjectTypeCompare((PetscObject)dmc, DMCOMPOSITE, &flg));
725: PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
726: next = com->next;
727: PetscCheck(!com->setup, PetscObjectComm((PetscObject)dmc), PETSC_ERR_ARG_WRONGSTATE, "Cannot add a DM once you have used the DMComposite");
729: /* create new link */
730: PetscCall(PetscNew(&mine));
731: PetscCall(PetscObjectReference((PetscObject)dm));
732: PetscCall(DMCreateGlobalVector(dm, &global)); // Not using DMGetGlobalVector() since it will stash the vector with a type decided by dm,
733: PetscCall(VecGetLocalSize(global, &n)); // while we may want to set dmc's vectype later, say via DMSetFromOptions(dmc), and we
734: PetscCall(VecDestroy(&global)); // want to propagate the type to dm.
735: PetscCall(DMCreateLocalVector(dm, &local)); // Not using DMGetLocalVector(), same reason as above.
736: PetscCall(VecGetSize(local, &nlocal));
737: PetscCall(VecDestroy(&local));
739: mine->n = n;
740: mine->nlocal = nlocal;
741: mine->dm = dm;
742: mine->next = NULL;
743: com->n += n;
744: com->nghost += nlocal;
746: /* add to end of list */
747: if (!next) com->next = mine;
748: else {
749: while (next->next) next = next->next;
750: next->next = mine;
751: }
752: com->nDM++;
753: com->nmine++;
754: PetscFunctionReturn(PETSC_SUCCESS);
755: }
757: #include <petscdraw.h>
758: PETSC_EXTERN PetscErrorCode VecView_MPI(Vec, PetscViewer);
759: static PetscErrorCode VecView_DMComposite(Vec gvec, PetscViewer viewer)
760: {
761: DM dm;
762: struct DMCompositeLink *next;
763: PetscBool isdraw;
764: DM_Composite *com;
766: PetscFunctionBegin;
767: PetscCall(VecGetDM(gvec, &dm));
768: PetscCheck(dm, PetscObjectComm((PetscObject)gvec), PETSC_ERR_ARG_WRONG, "Vector not generated from a DMComposite");
769: com = (DM_Composite *)dm->data;
770: next = com->next;
772: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw));
773: if (!isdraw) {
774: /* do I really want to call this? */
775: PetscCall(VecView_MPI(gvec, viewer));
776: } else {
777: PetscInt cnt = 0;
779: /* loop over packed objects, handling one at a time */
780: while (next) {
781: Vec vec;
782: const PetscScalar *array;
783: PetscInt bs;
785: /* Should use VecGetSubVector() eventually, but would need to forward the DM for that to work */
786: PetscCall(DMGetGlobalVector(next->dm, &vec));
787: PetscCall(VecGetArrayRead(gvec, &array));
788: PetscCall(VecPlaceArray(vec, (PetscScalar *)array + next->rstart));
789: PetscCall(VecRestoreArrayRead(gvec, &array));
790: PetscCall(VecView(vec, viewer));
791: PetscCall(VecResetArray(vec));
792: PetscCall(VecGetBlockSize(vec, &bs));
793: PetscCall(DMRestoreGlobalVector(next->dm, &vec));
794: PetscCall(PetscViewerDrawBaseAdd(viewer, bs));
795: cnt += bs;
796: next = next->next;
797: }
798: PetscCall(PetscViewerDrawBaseAdd(viewer, -cnt));
799: }
800: PetscFunctionReturn(PETSC_SUCCESS);
801: }
803: static PetscErrorCode DMCreateGlobalVector_Composite(DM dm, Vec *gvec)
804: {
805: DM_Composite *com = (DM_Composite *)dm->data;
807: PetscFunctionBegin;
809: PetscCall(DMSetFromOptions(dm));
810: PetscCall(DMSetUp(dm));
811: PetscCall(VecCreate(PetscObjectComm((PetscObject)dm), gvec));
812: PetscCall(VecSetType(*gvec, dm->vectype));
813: PetscCall(VecSetSizes(*gvec, com->n, com->N));
814: PetscCall(VecSetDM(*gvec, dm));
815: PetscCall(VecSetOperation(*gvec, VECOP_VIEW, (void (*)(void))VecView_DMComposite));
816: PetscFunctionReturn(PETSC_SUCCESS);
817: }
819: static PetscErrorCode DMCreateLocalVector_Composite(DM dm, Vec *lvec)
820: {
821: DM_Composite *com = (DM_Composite *)dm->data;
823: PetscFunctionBegin;
825: if (!com->setup) {
826: PetscCall(DMSetFromOptions(dm));
827: PetscCall(DMSetUp(dm));
828: }
829: PetscCall(VecCreate(PETSC_COMM_SELF, lvec));
830: PetscCall(VecSetType(*lvec, dm->vectype));
831: PetscCall(VecSetSizes(*lvec, com->nghost, PETSC_DECIDE));
832: PetscCall(VecSetDM(*lvec, dm));
833: PetscFunctionReturn(PETSC_SUCCESS);
834: }
836: /*@C
837: DMCompositeGetISLocalToGlobalMappings - gets an `ISLocalToGlobalMapping` for each `DM` in the `DMCOMPOSITE`, maps to the composite global space
839: Collective; No Fortran Support
841: Input Parameter:
842: . dm - the `DMCOMPOSITE` object
844: Output Parameter:
845: . ltogs - the individual mappings for each packed vector. Note that this includes
846: all the ghost points that individual ghosted `DMDA` may have.
848: Level: advanced
850: Note:
851: Each entry of `ltogs` should be destroyed with `ISLocalToGlobalMappingDestroy()`, `ltogs` should be freed with `PetscFree()`.
853: .seealso: `DMCOMPOSITE`, `DM`, `DMDestroy()`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`,
854: `DMCompositeGather()`, `DMCompositeCreate()`, `DMCompositeGetAccess()`, `DMCompositeScatter()`,
855: `DMCompositeGetLocalVectors()`, `DMCompositeRestoreLocalVectors()`, `DMCompositeGetEntries()`
856: @*/
857: PetscErrorCode DMCompositeGetISLocalToGlobalMappings(DM dm, ISLocalToGlobalMapping *ltogs[])
858: {
859: PetscInt i, *idx, n, cnt;
860: struct DMCompositeLink *next;
861: PetscMPIInt rank;
862: DM_Composite *com = (DM_Composite *)dm->data;
863: PetscBool flg;
865: PetscFunctionBegin;
867: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
868: PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
869: PetscCall(DMSetUp(dm));
870: PetscCall(PetscMalloc1(com->nDM, ltogs));
871: next = com->next;
872: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
874: /* loop over packed objects, handling one at a time */
875: cnt = 0;
876: while (next) {
877: ISLocalToGlobalMapping ltog;
878: PetscMPIInt size;
879: const PetscInt *suboff, *indices;
880: Vec global;
882: /* Get sub-DM global indices for each local dof */
883: PetscCall(DMGetLocalToGlobalMapping(next->dm, <og));
884: PetscCall(ISLocalToGlobalMappingGetSize(ltog, &n));
885: PetscCall(ISLocalToGlobalMappingGetIndices(ltog, &indices));
886: PetscCall(PetscMalloc1(n, &idx));
888: /* Get the offsets for the sub-DM global vector */
889: PetscCall(DMGetGlobalVector(next->dm, &global));
890: PetscCall(VecGetOwnershipRanges(global, &suboff));
891: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)global), &size));
893: /* Shift the sub-DM definition of the global space to the composite global space */
894: for (i = 0; i < n; i++) {
895: PetscInt subi = indices[i], lo = 0, hi = size, t;
896: /* There's no consensus on what a negative index means,
897: except for skipping when setting the values in vectors and matrices */
898: if (subi < 0) {
899: idx[i] = subi - next->grstarts[rank];
900: continue;
901: }
902: /* Binary search to find which rank owns subi */
903: while (hi - lo > 1) {
904: t = lo + (hi - lo) / 2;
905: if (suboff[t] > subi) hi = t;
906: else lo = t;
907: }
908: idx[i] = subi - suboff[lo] + next->grstarts[lo];
909: }
910: PetscCall(ISLocalToGlobalMappingRestoreIndices(ltog, &indices));
911: PetscCall(ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)dm), 1, n, idx, PETSC_OWN_POINTER, &(*ltogs)[cnt]));
912: PetscCall(DMRestoreGlobalVector(next->dm, &global));
913: next = next->next;
914: cnt++;
915: }
916: PetscFunctionReturn(PETSC_SUCCESS);
917: }
919: /*@C
920: DMCompositeGetLocalISs - Gets index sets for each component of a composite local vector
922: Not Collective; No Fortran Support
924: Input Parameter:
925: . dm - the `DMCOMPOSITE`
927: Output Parameter:
928: . is - array of serial index sets for each component of the `DMCOMPOSITE`
930: Level: intermediate
932: Notes:
933: At present, a composite local vector does not normally exist. This function is used to provide index sets for
934: `MatGetLocalSubMatrix()`. In the future, the scatters for each entry in the `DMCOMPOSITE` may be merged into a single
935: scatter to a composite local vector. The user should not typically need to know which is being done.
937: To get the composite global indices at all local points (including ghosts), use `DMCompositeGetISLocalToGlobalMappings()`.
939: To get index sets for pieces of the composite global vector, use `DMCompositeGetGlobalISs()`.
941: Each returned `IS` should be destroyed with `ISDestroy()`, the array should be freed with `PetscFree()`.
943: Fortran Note:
944: Use `DMCompositeRestoreLocalISs()` to release the `is`.
946: .seealso: `DMCOMPOSITE`, `DM`, `DMCompositeGetGlobalISs()`, `DMCompositeGetISLocalToGlobalMappings()`, `MatGetLocalSubMatrix()`,
947: `MatCreateLocalRef()`, `DMCompositeGetNumberDM()`
948: @*/
949: PetscErrorCode DMCompositeGetLocalISs(DM dm, IS *is[])
950: {
951: DM_Composite *com = (DM_Composite *)dm->data;
952: struct DMCompositeLink *link;
953: PetscInt cnt, start;
954: PetscBool flg;
956: PetscFunctionBegin;
958: PetscAssertPointer(is, 2);
959: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
960: PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
961: PetscCall(PetscMalloc1(com->nmine, is));
962: for (cnt = 0, start = 0, link = com->next; link; start += link->nlocal, cnt++, link = link->next) {
963: PetscInt bs;
964: PetscCall(ISCreateStride(PETSC_COMM_SELF, link->nlocal, start, 1, &(*is)[cnt]));
965: PetscCall(DMGetBlockSize(link->dm, &bs));
966: PetscCall(ISSetBlockSize((*is)[cnt], bs));
967: }
968: PetscFunctionReturn(PETSC_SUCCESS);
969: }
971: /*@C
972: DMCompositeGetGlobalISs - Gets the index sets for each composed object in a `DMCOMPOSITE`
974: Collective
976: Input Parameter:
977: . dm - the `DMCOMPOSITE` object
979: Output Parameter:
980: . is - the array of index sets
982: Level: advanced
984: Notes:
985: The `is` entries should be destroyed with `ISDestroy()`, `is` should be freed with `PetscFree()`
987: These could be used to extract a subset of vector entries for a "multi-physics" preconditioner
989: Use `DMCompositeGetLocalISs()` for index sets in the packed local numbering, and
990: `DMCompositeGetISLocalToGlobalMappings()` for to map local sub-`DM` (including ghost) indices to packed global
991: indices.
993: Fortran Note:
994: Use `DMCompositeRestoreGlobalISs()` to release the `is`.
996: .seealso: `DMCOMPOSITE`, `DM`, `DMDestroy()`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`,
997: `DMCompositeGather()`, `DMCompositeCreate()`, `DMCompositeGetAccess()`, `DMCompositeScatter()`,
998: `DMCompositeGetLocalVectors()`, `DMCompositeRestoreLocalVectors()`, `DMCompositeGetEntries()`
999: @*/
1000: PetscErrorCode DMCompositeGetGlobalISs(DM dm, IS *is[])
1001: {
1002: PetscInt cnt = 0;
1003: struct DMCompositeLink *next;
1004: PetscMPIInt rank;
1005: DM_Composite *com = (DM_Composite *)dm->data;
1006: PetscBool flg;
1008: PetscFunctionBegin;
1010: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
1011: PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
1012: PetscCheck(dm->setupcalled, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "Must call DMSetUp() before");
1013: PetscCall(PetscMalloc1(com->nDM, is));
1014: next = com->next;
1015: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
1017: /* loop over packed objects, handling one at a time */
1018: while (next) {
1019: PetscDS prob;
1021: PetscCall(ISCreateStride(PetscObjectComm((PetscObject)dm), next->n, next->grstart, 1, &(*is)[cnt]));
1022: PetscCall(DMGetDS(dm, &prob));
1023: if (prob) {
1024: MatNullSpace space;
1025: Mat pmat;
1026: PetscObject disc;
1027: PetscInt Nf;
1029: PetscCall(PetscDSGetNumFields(prob, &Nf));
1030: if (cnt < Nf) {
1031: PetscCall(PetscDSGetDiscretization(prob, cnt, &disc));
1032: PetscCall(PetscObjectQuery(disc, "nullspace", (PetscObject *)&space));
1033: if (space) PetscCall(PetscObjectCompose((PetscObject)(*is)[cnt], "nullspace", (PetscObject)space));
1034: PetscCall(PetscObjectQuery(disc, "nearnullspace", (PetscObject *)&space));
1035: if (space) PetscCall(PetscObjectCompose((PetscObject)(*is)[cnt], "nearnullspace", (PetscObject)space));
1036: PetscCall(PetscObjectQuery(disc, "pmat", (PetscObject *)&pmat));
1037: if (pmat) PetscCall(PetscObjectCompose((PetscObject)(*is)[cnt], "pmat", (PetscObject)pmat));
1038: }
1039: }
1040: cnt++;
1041: next = next->next;
1042: }
1043: PetscFunctionReturn(PETSC_SUCCESS);
1044: }
1046: static PetscErrorCode DMCreateFieldIS_Composite(DM dm, PetscInt *numFields, char ***fieldNames, IS **fields)
1047: {
1048: PetscInt nDM;
1049: DM *dms;
1050: PetscInt i;
1052: PetscFunctionBegin;
1053: PetscCall(DMCompositeGetNumberDM(dm, &nDM));
1054: if (numFields) *numFields = nDM;
1055: PetscCall(DMCompositeGetGlobalISs(dm, fields));
1056: if (fieldNames) {
1057: PetscCall(PetscMalloc1(nDM, &dms));
1058: PetscCall(PetscMalloc1(nDM, fieldNames));
1059: PetscCall(DMCompositeGetEntriesArray(dm, dms));
1060: for (i = 0; i < nDM; i++) {
1061: char buf[256];
1062: const char *splitname;
1064: /* Split naming precedence: object name, prefix, number */
1065: splitname = ((PetscObject)dm)->name;
1066: if (!splitname) {
1067: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)dms[i], &splitname));
1068: if (splitname) {
1069: size_t len;
1070: PetscCall(PetscStrncpy(buf, splitname, sizeof(buf)));
1071: buf[sizeof(buf) - 1] = 0;
1072: PetscCall(PetscStrlen(buf, &len));
1073: if (buf[len - 1] == '_') buf[len - 1] = 0; /* Remove trailing underscore if it was used */
1074: splitname = buf;
1075: }
1076: }
1077: if (!splitname) {
1078: PetscCall(PetscSNPrintf(buf, sizeof(buf), "%" PetscInt_FMT, i));
1079: splitname = buf;
1080: }
1081: PetscCall(PetscStrallocpy(splitname, &(*fieldNames)[i]));
1082: }
1083: PetscCall(PetscFree(dms));
1084: }
1085: PetscFunctionReturn(PETSC_SUCCESS);
1086: }
1088: /*
1089: This could take over from DMCreateFieldIS(), as it is more general,
1090: making DMCreateFieldIS() a special case -- calling with dmlist == NULL;
1091: At this point it's probably best to be less intrusive, however.
1092: */
1093: static PetscErrorCode DMCreateFieldDecomposition_Composite(DM dm, PetscInt *len, char ***namelist, IS **islist, DM **dmlist)
1094: {
1095: PetscInt nDM;
1096: PetscInt i;
1098: PetscFunctionBegin;
1099: PetscCall(DMCreateFieldIS_Composite(dm, len, namelist, islist));
1100: if (dmlist) {
1101: PetscCall(DMCompositeGetNumberDM(dm, &nDM));
1102: PetscCall(PetscMalloc1(nDM, dmlist));
1103: PetscCall(DMCompositeGetEntriesArray(dm, *dmlist));
1104: for (i = 0; i < nDM; i++) PetscCall(PetscObjectReference((PetscObject)((*dmlist)[i])));
1105: }
1106: PetscFunctionReturn(PETSC_SUCCESS);
1107: }
1109: /*@C
1110: DMCompositeGetLocalVectors - Gets local vectors for each part of a `DMCOMPOSITE`
1111: Use `DMCompositeRestoreLocalVectors()` to return them.
1113: Not Collective; No Fortran Support
1115: Input Parameter:
1116: . dm - the `DMCOMPOSITE` object
1118: Output Parameter:
1119: . ... - the individual sequential `Vec`s
1121: Level: advanced
1123: .seealso: `DMCOMPOSITE`, `DM`, `DMDestroy()`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`,
1124: `DMCompositeGather()`, `DMCompositeCreate()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeGetAccess()`,
1125: `DMCompositeRestoreLocalVectors()`, `DMCompositeScatter()`, `DMCompositeGetEntries()`
1126: @*/
1127: PetscErrorCode DMCompositeGetLocalVectors(DM dm, ...)
1128: {
1129: va_list Argp;
1130: struct DMCompositeLink *next;
1131: DM_Composite *com = (DM_Composite *)dm->data;
1132: PetscBool flg;
1134: PetscFunctionBegin;
1136: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
1137: PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
1138: next = com->next;
1139: /* loop over packed objects, handling one at a time */
1140: va_start(Argp, dm);
1141: while (next) {
1142: Vec *vec;
1143: vec = va_arg(Argp, Vec *);
1144: if (vec) PetscCall(DMGetLocalVector(next->dm, vec));
1145: next = next->next;
1146: }
1147: va_end(Argp);
1148: PetscFunctionReturn(PETSC_SUCCESS);
1149: }
1151: /*@C
1152: DMCompositeRestoreLocalVectors - Restores local vectors for each part of a `DMCOMPOSITE`
1154: Not Collective; No Fortran Support
1156: Input Parameter:
1157: . dm - the `DMCOMPOSITE` object
1159: Output Parameter:
1160: . ... - the individual sequential `Vec`
1162: Level: advanced
1164: .seealso: `DMCOMPOSITE`, `DM`, `DMDestroy()`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`,
1165: `DMCompositeGather()`, `DMCompositeCreate()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeGetAccess()`,
1166: `DMCompositeGetLocalVectors()`, `DMCompositeScatter()`, `DMCompositeGetEntries()`
1167: @*/
1168: PetscErrorCode DMCompositeRestoreLocalVectors(DM dm, ...)
1169: {
1170: va_list Argp;
1171: struct DMCompositeLink *next;
1172: DM_Composite *com = (DM_Composite *)dm->data;
1173: PetscBool flg;
1175: PetscFunctionBegin;
1177: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
1178: PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
1179: next = com->next;
1180: /* loop over packed objects, handling one at a time */
1181: va_start(Argp, dm);
1182: while (next) {
1183: Vec *vec;
1184: vec = va_arg(Argp, Vec *);
1185: if (vec) PetscCall(DMRestoreLocalVector(next->dm, vec));
1186: next = next->next;
1187: }
1188: va_end(Argp);
1189: PetscFunctionReturn(PETSC_SUCCESS);
1190: }
1192: /*@C
1193: DMCompositeGetEntries - Gets the `DM` for each entry in a `DMCOMPOSITE`.
1195: Not Collective
1197: Input Parameter:
1198: . dm - the `DMCOMPOSITE` object
1200: Output Parameter:
1201: . ... - the individual entries `DM`
1203: Level: advanced
1205: Fortran Notes:
1206: Use `DMCompositeGetEntriesArray()`
1208: .seealso: `DMCOMPOSITE`, `DM`, `DMDestroy()`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`, `DMCompositeGetEntriesArray()`
1209: `DMCompositeGather()`, `DMCompositeCreate()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeGetAccess()`,
1210: `DMCompositeRestoreLocalVectors()`, `DMCompositeGetLocalVectors()`, `DMCompositeScatter()`
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;
1219: PetscFunctionBegin;
1221: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
1222: PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
1223: next = com->next;
1224: /* loop over packed objects, handling one at a time */
1225: va_start(Argp, dm);
1226: while (next) {
1227: DM *dmn;
1228: dmn = va_arg(Argp, DM *);
1229: if (dmn) *dmn = next->dm;
1230: next = next->next;
1231: }
1232: va_end(Argp);
1233: PetscFunctionReturn(PETSC_SUCCESS);
1234: }
1236: /*@
1237: DMCompositeGetEntriesArray - Gets the DM for each entry in a `DMCOMPOSITE`
1239: Not Collective
1241: Input Parameter:
1242: . dm - the `DMCOMPOSITE` object
1244: Output Parameter:
1245: . dms - array of sufficient length (see `DMCompositeGetNumberDM()`) to hold the individual `DM`
1247: Level: advanced
1249: .seealso: `DMCOMPOSITE`, `DM`, `DMDestroy()`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`, `DMCompositeGetEntries()`
1250: `DMCompositeGather()`, `DMCompositeCreate()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeGetAccess()`,
1251: `DMCompositeRestoreLocalVectors()`, `DMCompositeGetLocalVectors()`, `DMCompositeScatter()`
1252: @*/
1253: PetscErrorCode DMCompositeGetEntriesArray(DM dm, DM dms[])
1254: {
1255: struct DMCompositeLink *next;
1256: DM_Composite *com = (DM_Composite *)dm->data;
1257: PetscInt i;
1258: PetscBool flg;
1260: PetscFunctionBegin;
1262: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
1263: PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
1264: /* loop over packed objects, handling one at a time */
1265: for (next = com->next, i = 0; next; next = next->next, i++) dms[i] = next->dm;
1266: PetscFunctionReturn(PETSC_SUCCESS);
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: PetscFunctionBegin;
1281: PetscCall(DMCompositeGetNumberDM(ctx->dm, &n));
1282: for (i = 0; i < n; i++) PetscCall(PetscViewerDestroy(&ctx->subv[i]));
1283: PetscCall(PetscFree2(ctx->subv, ctx->vecs));
1284: PetscCall(DMDestroy(&ctx->dm));
1285: PetscCall(PetscFree(ctx));
1286: PetscFunctionReturn(PETSC_SUCCESS);
1287: }
1289: static PetscErrorCode DMCompositeSampleGLVisFields_Private(PetscObject oX, PetscInt nf, PetscObject oXfield[], void *vctx)
1290: {
1291: Vec X = (Vec)oX;
1292: GLVisViewerCtx *ctx = (GLVisViewerCtx *)vctx;
1293: PetscInt i, n, cumf;
1295: PetscFunctionBegin;
1296: PetscCall(DMCompositeGetNumberDM(ctx->dm, &n));
1297: PetscCall(DMCompositeGetAccessArray(ctx->dm, X, n, NULL, ctx->vecs));
1298: for (i = 0, cumf = 0; i < n; i++) {
1299: PetscErrorCode (*g2l)(PetscObject, PetscInt, PetscObject[], void *);
1300: void *fctx;
1301: PetscInt nfi;
1303: PetscCall(PetscViewerGLVisGetFields_Internal(ctx->subv[i], &nfi, NULL, NULL, &g2l, NULL, &fctx));
1304: if (!nfi) continue;
1305: if (g2l) PetscCall((*g2l)((PetscObject)ctx->vecs[i], nfi, oXfield + cumf, fctx));
1306: else PetscCall(VecCopy(ctx->vecs[i], (Vec)oXfield[cumf]));
1307: cumf += nfi;
1308: }
1309: PetscCall(DMCompositeRestoreAccessArray(ctx->dm, X, n, NULL, ctx->vecs));
1310: PetscFunctionReturn(PETSC_SUCCESS);
1311: }
1313: static PetscErrorCode DMSetUpGLVisViewer_Composite(PetscObject odm, PetscViewer viewer)
1314: {
1315: DM dm = (DM)odm, *dms;
1316: Vec *Ufds;
1317: GLVisViewerCtx *ctx;
1318: PetscInt i, n, tnf, *sdim;
1319: char **fecs;
1321: PetscFunctionBegin;
1322: PetscCall(PetscNew(&ctx));
1323: PetscCall(PetscObjectReference((PetscObject)dm));
1324: ctx->dm = dm;
1325: PetscCall(DMCompositeGetNumberDM(dm, &n));
1326: PetscCall(PetscMalloc1(n, &dms));
1327: PetscCall(DMCompositeGetEntriesArray(dm, dms));
1328: PetscCall(PetscMalloc2(n, &ctx->subv, n, &ctx->vecs));
1329: for (i = 0, tnf = 0; i < n; i++) {
1330: PetscInt nf;
1332: PetscCall(PetscViewerCreate(PetscObjectComm(odm), &ctx->subv[i]));
1333: PetscCall(PetscViewerSetType(ctx->subv[i], PETSCVIEWERGLVIS));
1334: PetscCall(PetscViewerGLVisSetDM_Internal(ctx->subv[i], (PetscObject)dms[i]));
1335: PetscCall(PetscViewerGLVisGetFields_Internal(ctx->subv[i], &nf, NULL, NULL, NULL, NULL, NULL));
1336: tnf += nf;
1337: }
1338: PetscCall(PetscFree(dms));
1339: PetscCall(PetscMalloc3(tnf, &fecs, tnf, &sdim, tnf, &Ufds));
1340: for (i = 0, tnf = 0; i < n; i++) {
1341: PetscInt *sd, nf, f;
1342: const char **fec;
1343: Vec *Uf;
1345: PetscCall(PetscViewerGLVisGetFields_Internal(ctx->subv[i], &nf, &fec, &sd, NULL, (PetscObject **)&Uf, NULL));
1346: for (f = 0; f < nf; f++) {
1347: PetscCall(PetscStrallocpy(fec[f], &fecs[tnf + f]));
1348: Ufds[tnf + f] = Uf[f];
1349: sdim[tnf + f] = sd[f];
1350: }
1351: tnf += nf;
1352: }
1353: PetscCall(PetscViewerGLVisSetFields(viewer, tnf, (const char **)fecs, sdim, DMCompositeSampleGLVisFields_Private, (PetscObject *)Ufds, ctx, DestroyGLVisViewerCtx_Private));
1354: for (i = 0; i < tnf; i++) PetscCall(PetscFree(fecs[i]));
1355: PetscCall(PetscFree3(fecs, sdim, Ufds));
1356: PetscFunctionReturn(PETSC_SUCCESS);
1357: }
1359: static PetscErrorCode DMRefine_Composite(DM dmi, MPI_Comm comm, DM *fine)
1360: {
1361: struct DMCompositeLink *next;
1362: DM_Composite *com = (DM_Composite *)dmi->data;
1363: DM dm;
1365: PetscFunctionBegin;
1367: if (comm == MPI_COMM_NULL) PetscCall(PetscObjectGetComm((PetscObject)dmi, &comm));
1368: PetscCall(DMSetUp(dmi));
1369: next = com->next;
1370: PetscCall(DMCompositeCreate(comm, fine));
1372: /* loop over packed objects, handling one at a time */
1373: while (next) {
1374: PetscCall(DMRefine(next->dm, comm, &dm));
1375: PetscCall(DMCompositeAddDM(*fine, dm));
1376: PetscCall(PetscObjectDereference((PetscObject)dm));
1377: next = next->next;
1378: }
1379: PetscFunctionReturn(PETSC_SUCCESS);
1380: }
1382: static PetscErrorCode DMCoarsen_Composite(DM dmi, MPI_Comm comm, DM *fine)
1383: {
1384: struct DMCompositeLink *next;
1385: DM_Composite *com = (DM_Composite *)dmi->data;
1386: DM dm;
1388: PetscFunctionBegin;
1390: PetscCall(DMSetUp(dmi));
1391: if (comm == MPI_COMM_NULL) PetscCall(PetscObjectGetComm((PetscObject)dmi, &comm));
1392: next = com->next;
1393: PetscCall(DMCompositeCreate(comm, fine));
1395: /* loop over packed objects, handling one at a time */
1396: while (next) {
1397: PetscCall(DMCoarsen(next->dm, comm, &dm));
1398: PetscCall(DMCompositeAddDM(*fine, dm));
1399: PetscCall(PetscObjectDereference((PetscObject)dm));
1400: next = next->next;
1401: }
1402: PetscFunctionReturn(PETSC_SUCCESS);
1403: }
1405: static PetscErrorCode DMCreateInterpolation_Composite(DM coarse, DM fine, Mat *A, Vec *v)
1406: {
1407: PetscInt m, n, M, N, nDM, i;
1408: struct DMCompositeLink *nextc;
1409: struct DMCompositeLink *nextf;
1410: Vec gcoarse, gfine, *vecs;
1411: DM_Composite *comcoarse = (DM_Composite *)coarse->data;
1412: DM_Composite *comfine = (DM_Composite *)fine->data;
1413: Mat *mats;
1415: PetscFunctionBegin;
1418: PetscCall(DMSetUp(coarse));
1419: PetscCall(DMSetUp(fine));
1420: /* use global vectors only for determining matrix layout */
1421: PetscCall(DMGetGlobalVector(coarse, &gcoarse));
1422: PetscCall(DMGetGlobalVector(fine, &gfine));
1423: PetscCall(VecGetLocalSize(gcoarse, &n));
1424: PetscCall(VecGetLocalSize(gfine, &m));
1425: PetscCall(VecGetSize(gcoarse, &N));
1426: PetscCall(VecGetSize(gfine, &M));
1427: PetscCall(DMRestoreGlobalVector(coarse, &gcoarse));
1428: PetscCall(DMRestoreGlobalVector(fine, &gfine));
1430: nDM = comfine->nDM;
1431: PetscCheck(nDM == comcoarse->nDM, PetscObjectComm((PetscObject)fine), PETSC_ERR_ARG_INCOMP, "Fine DMComposite has %" PetscInt_FMT " entries, but coarse has %" PetscInt_FMT, nDM, comcoarse->nDM);
1432: PetscCall(PetscCalloc1(nDM * nDM, &mats));
1433: if (v) PetscCall(PetscCalloc1(nDM, &vecs));
1435: /* loop over packed objects, handling one at a time */
1436: for (nextc = comcoarse->next, nextf = comfine->next, i = 0; nextc; nextc = nextc->next, nextf = nextf->next, i++) {
1437: if (!v) PetscCall(DMCreateInterpolation(nextc->dm, nextf->dm, &mats[i * nDM + i], NULL));
1438: else PetscCall(DMCreateInterpolation(nextc->dm, nextf->dm, &mats[i * nDM + i], &vecs[i]));
1439: }
1440: PetscCall(MatCreateNest(PetscObjectComm((PetscObject)fine), nDM, NULL, nDM, NULL, mats, A));
1441: if (v) PetscCall(VecCreateNest(PetscObjectComm((PetscObject)fine), nDM, NULL, vecs, v));
1442: for (i = 0; i < nDM * nDM; i++) PetscCall(MatDestroy(&mats[i]));
1443: PetscCall(PetscFree(mats));
1444: if (v) {
1445: for (i = 0; i < nDM; i++) PetscCall(VecDestroy(&vecs[i]));
1446: PetscCall(PetscFree(vecs));
1447: }
1448: PetscFunctionReturn(PETSC_SUCCESS);
1449: }
1451: static PetscErrorCode DMGetLocalToGlobalMapping_Composite(DM dm)
1452: {
1453: DM_Composite *com = (DM_Composite *)dm->data;
1454: ISLocalToGlobalMapping *ltogs;
1455: PetscInt i;
1457: PetscFunctionBegin;
1458: /* Set the ISLocalToGlobalMapping on the new matrix */
1459: PetscCall(DMCompositeGetISLocalToGlobalMappings(dm, <ogs));
1460: PetscCall(ISLocalToGlobalMappingConcatenate(PetscObjectComm((PetscObject)dm), com->nDM, ltogs, &dm->ltogmap));
1461: for (i = 0; i < com->nDM; i++) PetscCall(ISLocalToGlobalMappingDestroy(<ogs[i]));
1462: PetscCall(PetscFree(ltogs));
1463: PetscFunctionReturn(PETSC_SUCCESS);
1464: }
1466: static PetscErrorCode DMCreateColoring_Composite(DM dm, ISColoringType ctype, ISColoring *coloring)
1467: {
1468: PetscInt n, i, cnt;
1469: ISColoringValue *colors;
1470: PetscBool dense = PETSC_FALSE;
1471: ISColoringValue maxcol = 0;
1472: DM_Composite *com = (DM_Composite *)dm->data;
1474: PetscFunctionBegin;
1476: PetscCheck(ctype != IS_COLORING_LOCAL, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Only global coloring supported");
1477: if (ctype == IS_COLORING_GLOBAL) {
1478: n = com->n;
1479: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Unknown ISColoringType");
1480: PetscCall(PetscMalloc1(n, &colors)); /* freed in ISColoringDestroy() */
1482: PetscCall(PetscOptionsGetBool(((PetscObject)dm)->options, ((PetscObject)dm)->prefix, "-dmcomposite_dense_jacobian", &dense, NULL));
1483: if (dense) {
1484: PetscCall(ISColoringValueCast(com->N, &maxcol));
1485: for (i = 0; i < n; i++) PetscCall(ISColoringValueCast(com->rstart + i, colors + i));
1486: } else {
1487: struct DMCompositeLink *next = com->next;
1488: PetscMPIInt rank;
1490: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
1491: cnt = 0;
1492: while (next) {
1493: ISColoring lcoloring;
1495: PetscCall(DMCreateColoring(next->dm, IS_COLORING_GLOBAL, &lcoloring));
1496: for (i = 0; i < lcoloring->N; i++) PetscCall(ISColoringValueCast(maxcol + lcoloring->colors[i], colors + cnt++));
1497: maxcol += lcoloring->n;
1498: PetscCall(ISColoringDestroy(&lcoloring));
1499: next = next->next;
1500: }
1501: }
1502: PetscCall(ISColoringCreate(PetscObjectComm((PetscObject)dm), maxcol, n, colors, PETSC_OWN_POINTER, coloring));
1503: PetscFunctionReturn(PETSC_SUCCESS);
1504: }
1506: static PetscErrorCode DMGlobalToLocalBegin_Composite(DM dm, Vec gvec, InsertMode mode, Vec lvec)
1507: {
1508: struct DMCompositeLink *next;
1509: PetscScalar *garray, *larray;
1510: DM_Composite *com = (DM_Composite *)dm->data;
1512: PetscFunctionBegin;
1516: if (!com->setup) PetscCall(DMSetUp(dm));
1518: PetscCall(VecGetArray(gvec, &garray));
1519: PetscCall(VecGetArray(lvec, &larray));
1521: /* loop over packed objects, handling one at a time */
1522: next = com->next;
1523: while (next) {
1524: Vec local, global;
1525: PetscInt N;
1527: PetscCall(DMGetGlobalVector(next->dm, &global));
1528: PetscCall(VecGetLocalSize(global, &N));
1529: PetscCall(VecPlaceArray(global, garray));
1530: PetscCall(DMGetLocalVector(next->dm, &local));
1531: PetscCall(VecPlaceArray(local, larray));
1532: PetscCall(DMGlobalToLocalBegin(next->dm, global, mode, local));
1533: PetscCall(DMGlobalToLocalEnd(next->dm, global, mode, local));
1534: PetscCall(VecResetArray(global));
1535: PetscCall(VecResetArray(local));
1536: PetscCall(DMRestoreGlobalVector(next->dm, &global));
1537: PetscCall(DMRestoreLocalVector(next->dm, &local));
1539: larray += next->nlocal;
1540: garray += next->n;
1541: next = next->next;
1542: }
1544: PetscCall(VecRestoreArray(gvec, NULL));
1545: PetscCall(VecRestoreArray(lvec, NULL));
1546: PetscFunctionReturn(PETSC_SUCCESS);
1547: }
1549: static PetscErrorCode DMGlobalToLocalEnd_Composite(DM dm, Vec gvec, InsertMode mode, Vec lvec)
1550: {
1551: PetscFunctionBegin;
1555: PetscFunctionReturn(PETSC_SUCCESS);
1556: }
1558: static PetscErrorCode DMLocalToGlobalBegin_Composite(DM dm, Vec lvec, InsertMode mode, Vec gvec)
1559: {
1560: struct DMCompositeLink *next;
1561: PetscScalar *larray, *garray;
1562: DM_Composite *com = (DM_Composite *)dm->data;
1564: PetscFunctionBegin;
1569: if (!com->setup) PetscCall(DMSetUp(dm));
1571: PetscCall(VecGetArray(lvec, &larray));
1572: PetscCall(VecGetArray(gvec, &garray));
1574: /* loop over packed objects, handling one at a time */
1575: next = com->next;
1576: while (next) {
1577: Vec global, local;
1579: PetscCall(DMGetLocalVector(next->dm, &local));
1580: PetscCall(VecPlaceArray(local, larray));
1581: PetscCall(DMGetGlobalVector(next->dm, &global));
1582: PetscCall(VecPlaceArray(global, garray));
1583: PetscCall(DMLocalToGlobalBegin(next->dm, local, mode, global));
1584: PetscCall(DMLocalToGlobalEnd(next->dm, local, mode, global));
1585: PetscCall(VecResetArray(local));
1586: PetscCall(VecResetArray(global));
1587: PetscCall(DMRestoreGlobalVector(next->dm, &global));
1588: PetscCall(DMRestoreLocalVector(next->dm, &local));
1590: garray += next->n;
1591: larray += next->nlocal;
1592: next = next->next;
1593: }
1595: PetscCall(VecRestoreArray(gvec, NULL));
1596: PetscCall(VecRestoreArray(lvec, NULL));
1597: PetscFunctionReturn(PETSC_SUCCESS);
1598: }
1600: static PetscErrorCode DMLocalToGlobalEnd_Composite(DM dm, Vec lvec, InsertMode mode, Vec gvec)
1601: {
1602: PetscFunctionBegin;
1606: PetscFunctionReturn(PETSC_SUCCESS);
1607: }
1609: static PetscErrorCode DMLocalToLocalBegin_Composite(DM dm, Vec vec1, InsertMode mode, Vec vec2)
1610: {
1611: struct DMCompositeLink *next;
1612: PetscScalar *array1, *array2;
1613: DM_Composite *com = (DM_Composite *)dm->data;
1615: PetscFunctionBegin;
1620: if (!com->setup) PetscCall(DMSetUp(dm));
1622: PetscCall(VecGetArray(vec1, &array1));
1623: PetscCall(VecGetArray(vec2, &array2));
1625: /* loop over packed objects, handling one at a time */
1626: next = com->next;
1627: while (next) {
1628: Vec local1, local2;
1630: PetscCall(DMGetLocalVector(next->dm, &local1));
1631: PetscCall(VecPlaceArray(local1, array1));
1632: PetscCall(DMGetLocalVector(next->dm, &local2));
1633: PetscCall(VecPlaceArray(local2, array2));
1634: PetscCall(DMLocalToLocalBegin(next->dm, local1, mode, local2));
1635: PetscCall(DMLocalToLocalEnd(next->dm, local1, mode, local2));
1636: PetscCall(VecResetArray(local2));
1637: PetscCall(DMRestoreLocalVector(next->dm, &local2));
1638: PetscCall(VecResetArray(local1));
1639: PetscCall(DMRestoreLocalVector(next->dm, &local1));
1641: array1 += next->nlocal;
1642: array2 += next->nlocal;
1643: next = next->next;
1644: }
1646: PetscCall(VecRestoreArray(vec1, NULL));
1647: PetscCall(VecRestoreArray(vec2, NULL));
1648: PetscFunctionReturn(PETSC_SUCCESS);
1649: }
1651: static PetscErrorCode DMLocalToLocalEnd_Composite(DM dm, Vec lvec, InsertMode mode, Vec gvec)
1652: {
1653: PetscFunctionBegin;
1657: PetscFunctionReturn(PETSC_SUCCESS);
1658: }
1660: /*MC
1661: DMCOMPOSITE = "composite" - A `DM` object that is used to manage data for a collection of `DM`
1663: Level: intermediate
1665: .seealso: `DMType`, `DM`, `DMDACreate()`, `DMCreate()`, `DMSetType()`, `DMCompositeCreate()`
1666: M*/
1668: PETSC_EXTERN PetscErrorCode DMCreate_Composite(DM p)
1669: {
1670: DM_Composite *com;
1672: PetscFunctionBegin;
1673: PetscCall(PetscNew(&com));
1674: p->data = com;
1675: com->n = 0;
1676: com->nghost = 0;
1677: com->next = NULL;
1678: com->nDM = 0;
1680: p->ops->createglobalvector = DMCreateGlobalVector_Composite;
1681: p->ops->createlocalvector = DMCreateLocalVector_Composite;
1682: p->ops->getlocaltoglobalmapping = DMGetLocalToGlobalMapping_Composite;
1683: p->ops->createfieldis = DMCreateFieldIS_Composite;
1684: p->ops->createfielddecomposition = DMCreateFieldDecomposition_Composite;
1685: p->ops->refine = DMRefine_Composite;
1686: p->ops->coarsen = DMCoarsen_Composite;
1687: p->ops->createinterpolation = DMCreateInterpolation_Composite;
1688: p->ops->creatematrix = DMCreateMatrix_Composite;
1689: p->ops->getcoloring = DMCreateColoring_Composite;
1690: p->ops->globaltolocalbegin = DMGlobalToLocalBegin_Composite;
1691: p->ops->globaltolocalend = DMGlobalToLocalEnd_Composite;
1692: p->ops->localtoglobalbegin = DMLocalToGlobalBegin_Composite;
1693: p->ops->localtoglobalend = DMLocalToGlobalEnd_Composite;
1694: p->ops->localtolocalbegin = DMLocalToLocalBegin_Composite;
1695: p->ops->localtolocalend = DMLocalToLocalEnd_Composite;
1696: p->ops->destroy = DMDestroy_Composite;
1697: p->ops->view = DMView_Composite;
1698: p->ops->setup = DMSetUp_Composite;
1700: PetscCall(PetscObjectComposeFunction((PetscObject)p, "DMSetUpGLVisViewer_C", DMSetUpGLVisViewer_Composite));
1701: PetscFunctionReturn(PETSC_SUCCESS);
1702: }
1704: /*@
1705: DMCompositeCreate - Creates a `DMCOMPOSITE`, used to generate "composite"
1706: vectors made up of several subvectors.
1708: Collective
1710: Input Parameter:
1711: . comm - the processors that will share the global vector
1713: Output Parameter:
1714: . packer - the `DMCOMPOSITE` object
1716: Level: advanced
1718: .seealso: `DMCOMPOSITE`, `DM`, `DMDestroy()`, `DMCompositeAddDM()`, `DMCompositeScatter()`, `DMCreate()`
1719: `DMCompositeGather()`, `DMCreateGlobalVector()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeGetAccess()`
1720: `DMCompositeGetLocalVectors()`, `DMCompositeRestoreLocalVectors()`, `DMCompositeGetEntries()`
1721: @*/
1722: PetscErrorCode DMCompositeCreate(MPI_Comm comm, DM *packer)
1723: {
1724: PetscFunctionBegin;
1725: PetscAssertPointer(packer, 2);
1726: PetscCall(DMCreate(comm, packer));
1727: PetscCall(DMSetType(*packer, DMCOMPOSITE));
1728: PetscFunctionReturn(PETSC_SUCCESS);
1729: }