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_SINGLE_LIBRARY_INTERN PetscErrorCode VecView_MPI(Vec, PetscViewer);
760: static PetscErrorCode VecView_DMComposite(Vec gvec, PetscViewer viewer)
761: {
762: DM dm;
763: struct DMCompositeLink *next;
764: PetscBool isdraw;
765: DM_Composite *com;
767: PetscFunctionBegin;
768: PetscCall(VecGetDM(gvec, &dm));
769: PetscCheck(dm, PetscObjectComm((PetscObject)gvec), PETSC_ERR_ARG_WRONG, "Vector not generated from a DMComposite");
770: com = (DM_Composite *)dm->data;
771: next = com->next;
773: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw));
774: if (!isdraw) {
775: /* do I really want to call this? */
776: PetscCall(VecView_MPI(gvec, viewer));
777: } else {
778: PetscInt cnt = 0;
780: /* loop over packed objects, handling one at a time */
781: while (next) {
782: Vec vec;
783: const PetscScalar *array;
784: PetscInt bs;
786: /* Should use VecGetSubVector() eventually, but would need to forward the DM for that to work */
787: PetscCall(DMGetGlobalVector(next->dm, &vec));
788: PetscCall(VecGetArrayRead(gvec, &array));
789: PetscCall(VecPlaceArray(vec, (PetscScalar *)array + next->rstart));
790: PetscCall(VecRestoreArrayRead(gvec, &array));
791: PetscCall(VecView(vec, viewer));
792: PetscCall(VecResetArray(vec));
793: PetscCall(VecGetBlockSize(vec, &bs));
794: PetscCall(DMRestoreGlobalVector(next->dm, &vec));
795: PetscCall(PetscViewerDrawBaseAdd(viewer, bs));
796: cnt += bs;
797: next = next->next;
798: }
799: PetscCall(PetscViewerDrawBaseAdd(viewer, -cnt));
800: }
801: PetscFunctionReturn(PETSC_SUCCESS);
802: }
804: static PetscErrorCode DMCreateGlobalVector_Composite(DM dm, Vec *gvec)
805: {
806: DM_Composite *com = (DM_Composite *)dm->data;
808: PetscFunctionBegin;
810: PetscCall(DMSetFromOptions(dm));
811: PetscCall(DMSetUp(dm));
812: PetscCall(VecCreate(PetscObjectComm((PetscObject)dm), gvec));
813: PetscCall(VecSetType(*gvec, dm->vectype));
814: PetscCall(VecSetSizes(*gvec, com->n, com->N));
815: PetscCall(VecSetDM(*gvec, dm));
816: PetscCall(VecSetOperation(*gvec, VECOP_VIEW, (void (*)(void))VecView_DMComposite));
817: PetscFunctionReturn(PETSC_SUCCESS);
818: }
820: static PetscErrorCode DMCreateLocalVector_Composite(DM dm, Vec *lvec)
821: {
822: DM_Composite *com = (DM_Composite *)dm->data;
824: PetscFunctionBegin;
826: if (!com->setup) {
827: PetscCall(DMSetFromOptions(dm));
828: PetscCall(DMSetUp(dm));
829: }
830: PetscCall(VecCreate(PETSC_COMM_SELF, lvec));
831: PetscCall(VecSetType(*lvec, dm->vectype));
832: PetscCall(VecSetSizes(*lvec, com->nghost, PETSC_DECIDE));
833: PetscCall(VecSetDM(*lvec, dm));
834: PetscFunctionReturn(PETSC_SUCCESS);
835: }
837: /*@C
838: DMCompositeGetISLocalToGlobalMappings - gets an `ISLocalToGlobalMapping` for each `DM` in the `DMCOMPOSITE`, maps to the composite global space
840: Collective; No Fortran Support
842: Input Parameter:
843: . dm - the `DMCOMPOSITE` object
845: Output Parameter:
846: . ltogs - the individual mappings for each packed vector. Note that this includes
847: all the ghost points that individual ghosted `DMDA` may have.
849: Level: advanced
851: Note:
852: Each entry of `ltogs` should be destroyed with `ISLocalToGlobalMappingDestroy()`, `ltogs` should be freed with `PetscFree()`.
854: .seealso: `DMCOMPOSITE`, `DM`, `DMDestroy()`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`,
855: `DMCompositeGather()`, `DMCompositeCreate()`, `DMCompositeGetAccess()`, `DMCompositeScatter()`,
856: `DMCompositeGetLocalVectors()`, `DMCompositeRestoreLocalVectors()`, `DMCompositeGetEntries()`
857: @*/
858: PetscErrorCode DMCompositeGetISLocalToGlobalMappings(DM dm, ISLocalToGlobalMapping *ltogs[])
859: {
860: PetscInt i, *idx, n, cnt;
861: struct DMCompositeLink *next;
862: PetscMPIInt rank;
863: DM_Composite *com = (DM_Composite *)dm->data;
864: PetscBool flg;
866: PetscFunctionBegin;
868: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
869: PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
870: PetscCall(DMSetUp(dm));
871: PetscCall(PetscMalloc1(com->nDM, ltogs));
872: next = com->next;
873: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
875: /* loop over packed objects, handling one at a time */
876: cnt = 0;
877: while (next) {
878: ISLocalToGlobalMapping ltog;
879: PetscMPIInt size;
880: const PetscInt *suboff, *indices;
881: Vec global;
883: /* Get sub-DM global indices for each local dof */
884: PetscCall(DMGetLocalToGlobalMapping(next->dm, <og));
885: PetscCall(ISLocalToGlobalMappingGetSize(ltog, &n));
886: PetscCall(ISLocalToGlobalMappingGetIndices(ltog, &indices));
887: PetscCall(PetscMalloc1(n, &idx));
889: /* Get the offsets for the sub-DM global vector */
890: PetscCall(DMGetGlobalVector(next->dm, &global));
891: PetscCall(VecGetOwnershipRanges(global, &suboff));
892: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)global), &size));
894: /* Shift the sub-DM definition of the global space to the composite global space */
895: for (i = 0; i < n; i++) {
896: PetscInt subi = indices[i], lo = 0, hi = size, t;
897: /* There's no consensus on what a negative index means,
898: except for skipping when setting the values in vectors and matrices */
899: if (subi < 0) {
900: idx[i] = subi - next->grstarts[rank];
901: continue;
902: }
903: /* Binary search to find which rank owns subi */
904: while (hi - lo > 1) {
905: t = lo + (hi - lo) / 2;
906: if (suboff[t] > subi) hi = t;
907: else lo = t;
908: }
909: idx[i] = subi - suboff[lo] + next->grstarts[lo];
910: }
911: PetscCall(ISLocalToGlobalMappingRestoreIndices(ltog, &indices));
912: PetscCall(ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)dm), 1, n, idx, PETSC_OWN_POINTER, &(*ltogs)[cnt]));
913: PetscCall(DMRestoreGlobalVector(next->dm, &global));
914: next = next->next;
915: cnt++;
916: }
917: PetscFunctionReturn(PETSC_SUCCESS);
918: }
920: /*@C
921: DMCompositeGetLocalISs - Gets index sets for each component of a composite local vector
923: Not Collective; No Fortran Support
925: Input Parameter:
926: . dm - the `DMCOMPOSITE`
928: Output Parameter:
929: . is - array of serial index sets for each component of the `DMCOMPOSITE`
931: Level: intermediate
933: Notes:
934: At present, a composite local vector does not normally exist. This function is used to provide index sets for
935: `MatGetLocalSubMatrix()`. In the future, the scatters for each entry in the `DMCOMPOSITE` may be merged into a single
936: scatter to a composite local vector. The user should not typically need to know which is being done.
938: To get the composite global indices at all local points (including ghosts), use `DMCompositeGetISLocalToGlobalMappings()`.
940: To get index sets for pieces of the composite global vector, use `DMCompositeGetGlobalISs()`.
942: Each returned `IS` should be destroyed with `ISDestroy()`, the array should be freed with `PetscFree()`.
944: Fortran Note:
945: Use `DMCompositeRestoreLocalISs()` to release the `is`.
947: .seealso: `DMCOMPOSITE`, `DM`, `DMCompositeGetGlobalISs()`, `DMCompositeGetISLocalToGlobalMappings()`, `MatGetLocalSubMatrix()`,
948: `MatCreateLocalRef()`, `DMCompositeGetNumberDM()`
949: @*/
950: PetscErrorCode DMCompositeGetLocalISs(DM dm, IS *is[])
951: {
952: DM_Composite *com = (DM_Composite *)dm->data;
953: struct DMCompositeLink *link;
954: PetscInt cnt, start;
955: PetscBool flg;
957: PetscFunctionBegin;
959: PetscAssertPointer(is, 2);
960: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
961: PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
962: PetscCall(PetscMalloc1(com->nmine, is));
963: for (cnt = 0, start = 0, link = com->next; link; start += link->nlocal, cnt++, link = link->next) {
964: PetscInt bs;
965: PetscCall(ISCreateStride(PETSC_COMM_SELF, link->nlocal, start, 1, &(*is)[cnt]));
966: PetscCall(DMGetBlockSize(link->dm, &bs));
967: PetscCall(ISSetBlockSize((*is)[cnt], bs));
968: }
969: PetscFunctionReturn(PETSC_SUCCESS);
970: }
972: /*@C
973: DMCompositeGetGlobalISs - Gets the index sets for each composed object in a `DMCOMPOSITE`
975: Collective
977: Input Parameter:
978: . dm - the `DMCOMPOSITE` object
980: Output Parameter:
981: . is - the array of index sets
983: Level: advanced
985: Notes:
986: The `is` entries should be destroyed with `ISDestroy()`, `is` should be freed with `PetscFree()`
988: These could be used to extract a subset of vector entries for a "multi-physics" preconditioner
990: Use `DMCompositeGetLocalISs()` for index sets in the packed local numbering, and
991: `DMCompositeGetISLocalToGlobalMappings()` for to map local sub-`DM` (including ghost) indices to packed global
992: indices.
994: Fortran Note:
995: Use `DMCompositeRestoreGlobalISs()` to release the `is`.
997: .seealso: `DMCOMPOSITE`, `DM`, `DMDestroy()`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`,
998: `DMCompositeGather()`, `DMCompositeCreate()`, `DMCompositeGetAccess()`, `DMCompositeScatter()`,
999: `DMCompositeGetLocalVectors()`, `DMCompositeRestoreLocalVectors()`, `DMCompositeGetEntries()`
1000: @*/
1001: PetscErrorCode DMCompositeGetGlobalISs(DM dm, IS *is[])
1002: {
1003: PetscInt cnt = 0;
1004: struct DMCompositeLink *next;
1005: PetscMPIInt rank;
1006: DM_Composite *com = (DM_Composite *)dm->data;
1007: PetscBool flg;
1009: PetscFunctionBegin;
1011: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
1012: PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
1013: PetscCheck(dm->setupcalled, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "Must call DMSetUp() before");
1014: PetscCall(PetscMalloc1(com->nDM, is));
1015: next = com->next;
1016: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
1018: /* loop over packed objects, handling one at a time */
1019: while (next) {
1020: PetscDS prob;
1022: PetscCall(ISCreateStride(PetscObjectComm((PetscObject)dm), next->n, next->grstart, 1, &(*is)[cnt]));
1023: PetscCall(DMGetDS(dm, &prob));
1024: if (prob) {
1025: MatNullSpace space;
1026: Mat pmat;
1027: PetscObject disc;
1028: PetscInt Nf;
1030: PetscCall(PetscDSGetNumFields(prob, &Nf));
1031: if (cnt < Nf) {
1032: PetscCall(PetscDSGetDiscretization(prob, cnt, &disc));
1033: PetscCall(PetscObjectQuery(disc, "nullspace", (PetscObject *)&space));
1034: if (space) PetscCall(PetscObjectCompose((PetscObject)(*is)[cnt], "nullspace", (PetscObject)space));
1035: PetscCall(PetscObjectQuery(disc, "nearnullspace", (PetscObject *)&space));
1036: if (space) PetscCall(PetscObjectCompose((PetscObject)(*is)[cnt], "nearnullspace", (PetscObject)space));
1037: PetscCall(PetscObjectQuery(disc, "pmat", (PetscObject *)&pmat));
1038: if (pmat) PetscCall(PetscObjectCompose((PetscObject)(*is)[cnt], "pmat", (PetscObject)pmat));
1039: }
1040: }
1041: cnt++;
1042: next = next->next;
1043: }
1044: PetscFunctionReturn(PETSC_SUCCESS);
1045: }
1047: static PetscErrorCode DMCreateFieldIS_Composite(DM dm, PetscInt *numFields, char ***fieldNames, IS **fields)
1048: {
1049: PetscInt nDM;
1050: DM *dms;
1051: PetscInt i;
1053: PetscFunctionBegin;
1054: PetscCall(DMCompositeGetNumberDM(dm, &nDM));
1055: if (numFields) *numFields = nDM;
1056: PetscCall(DMCompositeGetGlobalISs(dm, fields));
1057: if (fieldNames) {
1058: PetscCall(PetscMalloc1(nDM, &dms));
1059: PetscCall(PetscMalloc1(nDM, fieldNames));
1060: PetscCall(DMCompositeGetEntriesArray(dm, dms));
1061: for (i = 0; i < nDM; i++) {
1062: char buf[256];
1063: const char *splitname;
1065: /* Split naming precedence: object name, prefix, number */
1066: splitname = ((PetscObject)dm)->name;
1067: if (!splitname) {
1068: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)dms[i], &splitname));
1069: if (splitname) {
1070: size_t len;
1071: PetscCall(PetscStrncpy(buf, splitname, sizeof(buf)));
1072: buf[sizeof(buf) - 1] = 0;
1073: PetscCall(PetscStrlen(buf, &len));
1074: if (buf[len - 1] == '_') buf[len - 1] = 0; /* Remove trailing underscore if it was used */
1075: splitname = buf;
1076: }
1077: }
1078: if (!splitname) {
1079: PetscCall(PetscSNPrintf(buf, sizeof(buf), "%" PetscInt_FMT, i));
1080: splitname = buf;
1081: }
1082: PetscCall(PetscStrallocpy(splitname, &(*fieldNames)[i]));
1083: }
1084: PetscCall(PetscFree(dms));
1085: }
1086: PetscFunctionReturn(PETSC_SUCCESS);
1087: }
1089: /*
1090: This could take over from DMCreateFieldIS(), as it is more general,
1091: making DMCreateFieldIS() a special case -- calling with dmlist == NULL;
1092: At this point it's probably best to be less intrusive, however.
1093: */
1094: static PetscErrorCode DMCreateFieldDecomposition_Composite(DM dm, PetscInt *len, char ***namelist, IS **islist, DM **dmlist)
1095: {
1096: PetscInt nDM;
1097: PetscInt i;
1099: PetscFunctionBegin;
1100: PetscCall(DMCreateFieldIS_Composite(dm, len, namelist, islist));
1101: if (dmlist) {
1102: PetscCall(DMCompositeGetNumberDM(dm, &nDM));
1103: PetscCall(PetscMalloc1(nDM, dmlist));
1104: PetscCall(DMCompositeGetEntriesArray(dm, *dmlist));
1105: for (i = 0; i < nDM; i++) PetscCall(PetscObjectReference((PetscObject)((*dmlist)[i])));
1106: }
1107: PetscFunctionReturn(PETSC_SUCCESS);
1108: }
1110: /*@C
1111: DMCompositeGetLocalVectors - Gets local vectors for each part of a `DMCOMPOSITE`
1112: Use `DMCompositeRestoreLocalVectors()` to return them.
1114: Not Collective; No Fortran Support
1116: Input Parameter:
1117: . dm - the `DMCOMPOSITE` object
1119: Output Parameter:
1120: . ... - the individual sequential `Vec`s
1122: Level: advanced
1124: .seealso: `DMCOMPOSITE`, `DM`, `DMDestroy()`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`,
1125: `DMCompositeGather()`, `DMCompositeCreate()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeGetAccess()`,
1126: `DMCompositeRestoreLocalVectors()`, `DMCompositeScatter()`, `DMCompositeGetEntries()`
1127: @*/
1128: PetscErrorCode DMCompositeGetLocalVectors(DM dm, ...)
1129: {
1130: va_list Argp;
1131: struct DMCompositeLink *next;
1132: DM_Composite *com = (DM_Composite *)dm->data;
1133: PetscBool flg;
1135: PetscFunctionBegin;
1137: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
1138: PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
1139: next = com->next;
1140: /* loop over packed objects, handling one at a time */
1141: va_start(Argp, dm);
1142: while (next) {
1143: Vec *vec;
1144: vec = va_arg(Argp, Vec *);
1145: if (vec) PetscCall(DMGetLocalVector(next->dm, vec));
1146: next = next->next;
1147: }
1148: va_end(Argp);
1149: PetscFunctionReturn(PETSC_SUCCESS);
1150: }
1152: /*@C
1153: DMCompositeRestoreLocalVectors - Restores local vectors for each part of a `DMCOMPOSITE`
1155: Not Collective; No Fortran Support
1157: Input Parameter:
1158: . dm - the `DMCOMPOSITE` object
1160: Output Parameter:
1161: . ... - the individual sequential `Vec`
1163: Level: advanced
1165: .seealso: `DMCOMPOSITE`, `DM`, `DMDestroy()`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`,
1166: `DMCompositeGather()`, `DMCompositeCreate()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeGetAccess()`,
1167: `DMCompositeGetLocalVectors()`, `DMCompositeScatter()`, `DMCompositeGetEntries()`
1168: @*/
1169: PetscErrorCode DMCompositeRestoreLocalVectors(DM dm, ...)
1170: {
1171: va_list Argp;
1172: struct DMCompositeLink *next;
1173: DM_Composite *com = (DM_Composite *)dm->data;
1174: PetscBool flg;
1176: PetscFunctionBegin;
1178: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
1179: PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
1180: next = com->next;
1181: /* loop over packed objects, handling one at a time */
1182: va_start(Argp, dm);
1183: while (next) {
1184: Vec *vec;
1185: vec = va_arg(Argp, Vec *);
1186: if (vec) PetscCall(DMRestoreLocalVector(next->dm, vec));
1187: next = next->next;
1188: }
1189: va_end(Argp);
1190: PetscFunctionReturn(PETSC_SUCCESS);
1191: }
1193: /*@C
1194: DMCompositeGetEntries - Gets the `DM` for each entry in a `DMCOMPOSITE`.
1196: Not Collective
1198: Input Parameter:
1199: . dm - the `DMCOMPOSITE` object
1201: Output Parameter:
1202: . ... - the individual entries `DM`
1204: Level: advanced
1206: Fortran Notes:
1207: Use `DMCompositeGetEntriesArray()`
1209: .seealso: `DMCOMPOSITE`, `DM`, `DMDestroy()`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`, `DMCompositeGetEntriesArray()`
1210: `DMCompositeGather()`, `DMCompositeCreate()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeGetAccess()`,
1211: `DMCompositeRestoreLocalVectors()`, `DMCompositeGetLocalVectors()`, `DMCompositeScatter()`
1212: @*/
1213: PetscErrorCode DMCompositeGetEntries(DM dm, ...)
1214: {
1215: va_list Argp;
1216: struct DMCompositeLink *next;
1217: DM_Composite *com = (DM_Composite *)dm->data;
1218: PetscBool flg;
1220: PetscFunctionBegin;
1222: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
1223: PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
1224: next = com->next;
1225: /* loop over packed objects, handling one at a time */
1226: va_start(Argp, dm);
1227: while (next) {
1228: DM *dmn;
1229: dmn = va_arg(Argp, DM *);
1230: if (dmn) *dmn = next->dm;
1231: next = next->next;
1232: }
1233: va_end(Argp);
1234: PetscFunctionReturn(PETSC_SUCCESS);
1235: }
1237: /*@
1238: DMCompositeGetEntriesArray - Gets the DM for each entry in a `DMCOMPOSITE`
1240: Not Collective
1242: Input Parameter:
1243: . dm - the `DMCOMPOSITE` object
1245: Output Parameter:
1246: . dms - array of sufficient length (see `DMCompositeGetNumberDM()`) to hold the individual `DM`
1248: Level: advanced
1250: .seealso: `DMCOMPOSITE`, `DM`, `DMDestroy()`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`, `DMCompositeGetEntries()`
1251: `DMCompositeGather()`, `DMCompositeCreate()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeGetAccess()`,
1252: `DMCompositeRestoreLocalVectors()`, `DMCompositeGetLocalVectors()`, `DMCompositeScatter()`
1253: @*/
1254: PetscErrorCode DMCompositeGetEntriesArray(DM dm, DM dms[])
1255: {
1256: struct DMCompositeLink *next;
1257: DM_Composite *com = (DM_Composite *)dm->data;
1258: PetscInt i;
1259: PetscBool flg;
1261: PetscFunctionBegin;
1263: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
1264: PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
1265: /* loop over packed objects, handling one at a time */
1266: for (next = com->next, i = 0; next; next = next->next, i++) dms[i] = next->dm;
1267: PetscFunctionReturn(PETSC_SUCCESS);
1268: }
1270: typedef struct {
1271: DM dm;
1272: PetscViewer *subv;
1273: Vec *vecs;
1274: } GLVisViewerCtx;
1276: static PetscErrorCode DestroyGLVisViewerCtx_Private(void *vctx)
1277: {
1278: GLVisViewerCtx *ctx = (GLVisViewerCtx *)vctx;
1279: PetscInt i, n;
1281: PetscFunctionBegin;
1282: PetscCall(DMCompositeGetNumberDM(ctx->dm, &n));
1283: for (i = 0; i < n; i++) PetscCall(PetscViewerDestroy(&ctx->subv[i]));
1284: PetscCall(PetscFree2(ctx->subv, ctx->vecs));
1285: PetscCall(DMDestroy(&ctx->dm));
1286: PetscCall(PetscFree(ctx));
1287: PetscFunctionReturn(PETSC_SUCCESS);
1288: }
1290: static PetscErrorCode DMCompositeSampleGLVisFields_Private(PetscObject oX, PetscInt nf, PetscObject oXfield[], void *vctx)
1291: {
1292: Vec X = (Vec)oX;
1293: GLVisViewerCtx *ctx = (GLVisViewerCtx *)vctx;
1294: PetscInt i, n, cumf;
1296: PetscFunctionBegin;
1297: PetscCall(DMCompositeGetNumberDM(ctx->dm, &n));
1298: PetscCall(DMCompositeGetAccessArray(ctx->dm, X, n, NULL, ctx->vecs));
1299: for (i = 0, cumf = 0; i < n; i++) {
1300: PetscErrorCode (*g2l)(PetscObject, PetscInt, PetscObject[], void *);
1301: void *fctx;
1302: PetscInt nfi;
1304: PetscCall(PetscViewerGLVisGetFields_Internal(ctx->subv[i], &nfi, NULL, NULL, &g2l, NULL, &fctx));
1305: if (!nfi) continue;
1306: if (g2l) PetscCall((*g2l)((PetscObject)ctx->vecs[i], nfi, oXfield + cumf, fctx));
1307: else PetscCall(VecCopy(ctx->vecs[i], (Vec)oXfield[cumf]));
1308: cumf += nfi;
1309: }
1310: PetscCall(DMCompositeRestoreAccessArray(ctx->dm, X, n, NULL, ctx->vecs));
1311: PetscFunctionReturn(PETSC_SUCCESS);
1312: }
1314: static PetscErrorCode DMSetUpGLVisViewer_Composite(PetscObject odm, PetscViewer viewer)
1315: {
1316: DM dm = (DM)odm, *dms;
1317: Vec *Ufds;
1318: GLVisViewerCtx *ctx;
1319: PetscInt i, n, tnf, *sdim;
1320: char **fecs;
1322: PetscFunctionBegin;
1323: PetscCall(PetscNew(&ctx));
1324: PetscCall(PetscObjectReference((PetscObject)dm));
1325: ctx->dm = dm;
1326: PetscCall(DMCompositeGetNumberDM(dm, &n));
1327: PetscCall(PetscMalloc1(n, &dms));
1328: PetscCall(DMCompositeGetEntriesArray(dm, dms));
1329: PetscCall(PetscMalloc2(n, &ctx->subv, n, &ctx->vecs));
1330: for (i = 0, tnf = 0; i < n; i++) {
1331: PetscInt nf;
1333: PetscCall(PetscViewerCreate(PetscObjectComm(odm), &ctx->subv[i]));
1334: PetscCall(PetscViewerSetType(ctx->subv[i], PETSCVIEWERGLVIS));
1335: PetscCall(PetscViewerGLVisSetDM_Internal(ctx->subv[i], (PetscObject)dms[i]));
1336: PetscCall(PetscViewerGLVisGetFields_Internal(ctx->subv[i], &nf, NULL, NULL, NULL, NULL, NULL));
1337: tnf += nf;
1338: }
1339: PetscCall(PetscFree(dms));
1340: PetscCall(PetscMalloc3(tnf, &fecs, tnf, &sdim, tnf, &Ufds));
1341: for (i = 0, tnf = 0; i < n; i++) {
1342: PetscInt *sd, nf, f;
1343: const char **fec;
1344: Vec *Uf;
1346: PetscCall(PetscViewerGLVisGetFields_Internal(ctx->subv[i], &nf, &fec, &sd, NULL, (PetscObject **)&Uf, NULL));
1347: for (f = 0; f < nf; f++) {
1348: PetscCall(PetscStrallocpy(fec[f], &fecs[tnf + f]));
1349: Ufds[tnf + f] = Uf[f];
1350: sdim[tnf + f] = sd[f];
1351: }
1352: tnf += nf;
1353: }
1354: PetscCall(PetscViewerGLVisSetFields(viewer, tnf, (const char **)fecs, sdim, DMCompositeSampleGLVisFields_Private, (PetscObject *)Ufds, ctx, DestroyGLVisViewerCtx_Private));
1355: for (i = 0; i < tnf; i++) PetscCall(PetscFree(fecs[i]));
1356: PetscCall(PetscFree3(fecs, sdim, Ufds));
1357: PetscFunctionReturn(PETSC_SUCCESS);
1358: }
1360: static PetscErrorCode DMRefine_Composite(DM dmi, MPI_Comm comm, DM *fine)
1361: {
1362: struct DMCompositeLink *next;
1363: DM_Composite *com = (DM_Composite *)dmi->data;
1364: DM dm;
1366: PetscFunctionBegin;
1368: if (comm == MPI_COMM_NULL) PetscCall(PetscObjectGetComm((PetscObject)dmi, &comm));
1369: PetscCall(DMSetUp(dmi));
1370: next = com->next;
1371: PetscCall(DMCompositeCreate(comm, fine));
1373: /* loop over packed objects, handling one at a time */
1374: while (next) {
1375: PetscCall(DMRefine(next->dm, comm, &dm));
1376: PetscCall(DMCompositeAddDM(*fine, dm));
1377: PetscCall(PetscObjectDereference((PetscObject)dm));
1378: next = next->next;
1379: }
1380: PetscFunctionReturn(PETSC_SUCCESS);
1381: }
1383: static PetscErrorCode DMCoarsen_Composite(DM dmi, MPI_Comm comm, DM *fine)
1384: {
1385: struct DMCompositeLink *next;
1386: DM_Composite *com = (DM_Composite *)dmi->data;
1387: DM dm;
1389: PetscFunctionBegin;
1391: PetscCall(DMSetUp(dmi));
1392: if (comm == MPI_COMM_NULL) PetscCall(PetscObjectGetComm((PetscObject)dmi, &comm));
1393: next = com->next;
1394: PetscCall(DMCompositeCreate(comm, fine));
1396: /* loop over packed objects, handling one at a time */
1397: while (next) {
1398: PetscCall(DMCoarsen(next->dm, comm, &dm));
1399: PetscCall(DMCompositeAddDM(*fine, dm));
1400: PetscCall(PetscObjectDereference((PetscObject)dm));
1401: next = next->next;
1402: }
1403: PetscFunctionReturn(PETSC_SUCCESS);
1404: }
1406: static PetscErrorCode DMCreateInterpolation_Composite(DM coarse, DM fine, Mat *A, Vec *v)
1407: {
1408: PetscInt m, n, M, N, nDM, i;
1409: struct DMCompositeLink *nextc;
1410: struct DMCompositeLink *nextf;
1411: Vec gcoarse, gfine, *vecs;
1412: DM_Composite *comcoarse = (DM_Composite *)coarse->data;
1413: DM_Composite *comfine = (DM_Composite *)fine->data;
1414: Mat *mats;
1416: PetscFunctionBegin;
1419: PetscCall(DMSetUp(coarse));
1420: PetscCall(DMSetUp(fine));
1421: /* use global vectors only for determining matrix layout */
1422: PetscCall(DMGetGlobalVector(coarse, &gcoarse));
1423: PetscCall(DMGetGlobalVector(fine, &gfine));
1424: PetscCall(VecGetLocalSize(gcoarse, &n));
1425: PetscCall(VecGetLocalSize(gfine, &m));
1426: PetscCall(VecGetSize(gcoarse, &N));
1427: PetscCall(VecGetSize(gfine, &M));
1428: PetscCall(DMRestoreGlobalVector(coarse, &gcoarse));
1429: PetscCall(DMRestoreGlobalVector(fine, &gfine));
1431: nDM = comfine->nDM;
1432: 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);
1433: PetscCall(PetscCalloc1(nDM * nDM, &mats));
1434: if (v) PetscCall(PetscCalloc1(nDM, &vecs));
1436: /* loop over packed objects, handling one at a time */
1437: for (nextc = comcoarse->next, nextf = comfine->next, i = 0; nextc; nextc = nextc->next, nextf = nextf->next, i++) {
1438: if (!v) PetscCall(DMCreateInterpolation(nextc->dm, nextf->dm, &mats[i * nDM + i], NULL));
1439: else PetscCall(DMCreateInterpolation(nextc->dm, nextf->dm, &mats[i * nDM + i], &vecs[i]));
1440: }
1441: PetscCall(MatCreateNest(PetscObjectComm((PetscObject)fine), nDM, NULL, nDM, NULL, mats, A));
1442: if (v) PetscCall(VecCreateNest(PetscObjectComm((PetscObject)fine), nDM, NULL, vecs, v));
1443: for (i = 0; i < nDM * nDM; i++) PetscCall(MatDestroy(&mats[i]));
1444: PetscCall(PetscFree(mats));
1445: if (v) {
1446: for (i = 0; i < nDM; i++) PetscCall(VecDestroy(&vecs[i]));
1447: PetscCall(PetscFree(vecs));
1448: }
1449: PetscFunctionReturn(PETSC_SUCCESS);
1450: }
1452: static PetscErrorCode DMGetLocalToGlobalMapping_Composite(DM dm)
1453: {
1454: DM_Composite *com = (DM_Composite *)dm->data;
1455: ISLocalToGlobalMapping *ltogs;
1456: PetscInt i;
1458: PetscFunctionBegin;
1459: /* Set the ISLocalToGlobalMapping on the new matrix */
1460: PetscCall(DMCompositeGetISLocalToGlobalMappings(dm, <ogs));
1461: PetscCall(ISLocalToGlobalMappingConcatenate(PetscObjectComm((PetscObject)dm), com->nDM, ltogs, &dm->ltogmap));
1462: for (i = 0; i < com->nDM; i++) PetscCall(ISLocalToGlobalMappingDestroy(<ogs[i]));
1463: PetscCall(PetscFree(ltogs));
1464: PetscFunctionReturn(PETSC_SUCCESS);
1465: }
1467: static PetscErrorCode DMCreateColoring_Composite(DM dm, ISColoringType ctype, ISColoring *coloring)
1468: {
1469: PetscInt n, i, cnt;
1470: ISColoringValue *colors;
1471: PetscBool dense = PETSC_FALSE;
1472: ISColoringValue maxcol = 0;
1473: DM_Composite *com = (DM_Composite *)dm->data;
1475: PetscFunctionBegin;
1477: PetscCheck(ctype != IS_COLORING_LOCAL, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Only global coloring supported");
1478: if (ctype == IS_COLORING_GLOBAL) {
1479: n = com->n;
1480: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Unknown ISColoringType");
1481: PetscCall(PetscMalloc1(n, &colors)); /* freed in ISColoringDestroy() */
1483: PetscCall(PetscOptionsGetBool(((PetscObject)dm)->options, ((PetscObject)dm)->prefix, "-dmcomposite_dense_jacobian", &dense, NULL));
1484: if (dense) {
1485: PetscCall(ISColoringValueCast(com->N, &maxcol));
1486: for (i = 0; i < n; i++) PetscCall(ISColoringValueCast(com->rstart + i, colors + i));
1487: } else {
1488: struct DMCompositeLink *next = com->next;
1489: PetscMPIInt rank;
1491: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
1492: cnt = 0;
1493: while (next) {
1494: ISColoring lcoloring;
1496: PetscCall(DMCreateColoring(next->dm, IS_COLORING_GLOBAL, &lcoloring));
1497: for (i = 0; i < lcoloring->N; i++) PetscCall(ISColoringValueCast(maxcol + lcoloring->colors[i], colors + cnt++));
1498: maxcol += lcoloring->n;
1499: PetscCall(ISColoringDestroy(&lcoloring));
1500: next = next->next;
1501: }
1502: }
1503: PetscCall(ISColoringCreate(PetscObjectComm((PetscObject)dm), maxcol, n, colors, PETSC_OWN_POINTER, coloring));
1504: PetscFunctionReturn(PETSC_SUCCESS);
1505: }
1507: static PetscErrorCode DMGlobalToLocalBegin_Composite(DM dm, Vec gvec, InsertMode mode, Vec lvec)
1508: {
1509: struct DMCompositeLink *next;
1510: PetscScalar *garray, *larray;
1511: DM_Composite *com = (DM_Composite *)dm->data;
1513: PetscFunctionBegin;
1517: if (!com->setup) PetscCall(DMSetUp(dm));
1519: PetscCall(VecGetArray(gvec, &garray));
1520: PetscCall(VecGetArray(lvec, &larray));
1522: /* loop over packed objects, handling one at a time */
1523: next = com->next;
1524: while (next) {
1525: Vec local, global;
1526: PetscInt N;
1528: PetscCall(DMGetGlobalVector(next->dm, &global));
1529: PetscCall(VecGetLocalSize(global, &N));
1530: PetscCall(VecPlaceArray(global, garray));
1531: PetscCall(DMGetLocalVector(next->dm, &local));
1532: PetscCall(VecPlaceArray(local, larray));
1533: PetscCall(DMGlobalToLocalBegin(next->dm, global, mode, local));
1534: PetscCall(DMGlobalToLocalEnd(next->dm, global, mode, local));
1535: PetscCall(VecResetArray(global));
1536: PetscCall(VecResetArray(local));
1537: PetscCall(DMRestoreGlobalVector(next->dm, &global));
1538: PetscCall(DMRestoreLocalVector(next->dm, &local));
1540: larray += next->nlocal;
1541: garray += next->n;
1542: next = next->next;
1543: }
1545: PetscCall(VecRestoreArray(gvec, NULL));
1546: PetscCall(VecRestoreArray(lvec, NULL));
1547: PetscFunctionReturn(PETSC_SUCCESS);
1548: }
1550: static PetscErrorCode DMGlobalToLocalEnd_Composite(DM dm, Vec gvec, InsertMode mode, Vec lvec)
1551: {
1552: PetscFunctionBegin;
1556: PetscFunctionReturn(PETSC_SUCCESS);
1557: }
1559: static PetscErrorCode DMLocalToGlobalBegin_Composite(DM dm, Vec lvec, InsertMode mode, Vec gvec)
1560: {
1561: struct DMCompositeLink *next;
1562: PetscScalar *larray, *garray;
1563: DM_Composite *com = (DM_Composite *)dm->data;
1565: PetscFunctionBegin;
1570: if (!com->setup) PetscCall(DMSetUp(dm));
1572: PetscCall(VecGetArray(lvec, &larray));
1573: PetscCall(VecGetArray(gvec, &garray));
1575: /* loop over packed objects, handling one at a time */
1576: next = com->next;
1577: while (next) {
1578: Vec global, local;
1580: PetscCall(DMGetLocalVector(next->dm, &local));
1581: PetscCall(VecPlaceArray(local, larray));
1582: PetscCall(DMGetGlobalVector(next->dm, &global));
1583: PetscCall(VecPlaceArray(global, garray));
1584: PetscCall(DMLocalToGlobalBegin(next->dm, local, mode, global));
1585: PetscCall(DMLocalToGlobalEnd(next->dm, local, mode, global));
1586: PetscCall(VecResetArray(local));
1587: PetscCall(VecResetArray(global));
1588: PetscCall(DMRestoreGlobalVector(next->dm, &global));
1589: PetscCall(DMRestoreLocalVector(next->dm, &local));
1591: garray += next->n;
1592: larray += next->nlocal;
1593: next = next->next;
1594: }
1596: PetscCall(VecRestoreArray(gvec, NULL));
1597: PetscCall(VecRestoreArray(lvec, NULL));
1598: PetscFunctionReturn(PETSC_SUCCESS);
1599: }
1601: static PetscErrorCode DMLocalToGlobalEnd_Composite(DM dm, Vec lvec, InsertMode mode, Vec gvec)
1602: {
1603: PetscFunctionBegin;
1607: PetscFunctionReturn(PETSC_SUCCESS);
1608: }
1610: static PetscErrorCode DMLocalToLocalBegin_Composite(DM dm, Vec vec1, InsertMode mode, Vec vec2)
1611: {
1612: struct DMCompositeLink *next;
1613: PetscScalar *array1, *array2;
1614: DM_Composite *com = (DM_Composite *)dm->data;
1616: PetscFunctionBegin;
1621: if (!com->setup) PetscCall(DMSetUp(dm));
1623: PetscCall(VecGetArray(vec1, &array1));
1624: PetscCall(VecGetArray(vec2, &array2));
1626: /* loop over packed objects, handling one at a time */
1627: next = com->next;
1628: while (next) {
1629: Vec local1, local2;
1631: PetscCall(DMGetLocalVector(next->dm, &local1));
1632: PetscCall(VecPlaceArray(local1, array1));
1633: PetscCall(DMGetLocalVector(next->dm, &local2));
1634: PetscCall(VecPlaceArray(local2, array2));
1635: PetscCall(DMLocalToLocalBegin(next->dm, local1, mode, local2));
1636: PetscCall(DMLocalToLocalEnd(next->dm, local1, mode, local2));
1637: PetscCall(VecResetArray(local2));
1638: PetscCall(DMRestoreLocalVector(next->dm, &local2));
1639: PetscCall(VecResetArray(local1));
1640: PetscCall(DMRestoreLocalVector(next->dm, &local1));
1642: array1 += next->nlocal;
1643: array2 += next->nlocal;
1644: next = next->next;
1645: }
1647: PetscCall(VecRestoreArray(vec1, NULL));
1648: PetscCall(VecRestoreArray(vec2, NULL));
1649: PetscFunctionReturn(PETSC_SUCCESS);
1650: }
1652: static PetscErrorCode DMLocalToLocalEnd_Composite(DM dm, Vec lvec, InsertMode mode, Vec gvec)
1653: {
1654: PetscFunctionBegin;
1658: PetscFunctionReturn(PETSC_SUCCESS);
1659: }
1661: /*MC
1662: DMCOMPOSITE = "composite" - A `DM` object that is used to manage data for a collection of `DM`
1664: Level: intermediate
1666: .seealso: `DMType`, `DM`, `DMDACreate()`, `DMCreate()`, `DMSetType()`, `DMCompositeCreate()`
1667: M*/
1669: PETSC_EXTERN PetscErrorCode DMCreate_Composite(DM p)
1670: {
1671: DM_Composite *com;
1673: PetscFunctionBegin;
1674: PetscCall(PetscNew(&com));
1675: p->data = com;
1676: com->n = 0;
1677: com->nghost = 0;
1678: com->next = NULL;
1679: com->nDM = 0;
1681: p->ops->createglobalvector = DMCreateGlobalVector_Composite;
1682: p->ops->createlocalvector = DMCreateLocalVector_Composite;
1683: p->ops->getlocaltoglobalmapping = DMGetLocalToGlobalMapping_Composite;
1684: p->ops->createfieldis = DMCreateFieldIS_Composite;
1685: p->ops->createfielddecomposition = DMCreateFieldDecomposition_Composite;
1686: p->ops->refine = DMRefine_Composite;
1687: p->ops->coarsen = DMCoarsen_Composite;
1688: p->ops->createinterpolation = DMCreateInterpolation_Composite;
1689: p->ops->creatematrix = DMCreateMatrix_Composite;
1690: p->ops->getcoloring = DMCreateColoring_Composite;
1691: p->ops->globaltolocalbegin = DMGlobalToLocalBegin_Composite;
1692: p->ops->globaltolocalend = DMGlobalToLocalEnd_Composite;
1693: p->ops->localtoglobalbegin = DMLocalToGlobalBegin_Composite;
1694: p->ops->localtoglobalend = DMLocalToGlobalEnd_Composite;
1695: p->ops->localtolocalbegin = DMLocalToLocalBegin_Composite;
1696: p->ops->localtolocalend = DMLocalToLocalEnd_Composite;
1697: p->ops->destroy = DMDestroy_Composite;
1698: p->ops->view = DMView_Composite;
1699: p->ops->setup = DMSetUp_Composite;
1701: PetscCall(PetscObjectComposeFunction((PetscObject)p, "DMSetUpGLVisViewer_C", DMSetUpGLVisViewer_Composite));
1702: PetscFunctionReturn(PETSC_SUCCESS);
1703: }
1705: /*@
1706: DMCompositeCreate - Creates a `DMCOMPOSITE`, used to generate "composite"
1707: vectors made up of several subvectors.
1709: Collective
1711: Input Parameter:
1712: . comm - the processors that will share the global vector
1714: Output Parameter:
1715: . packer - the `DMCOMPOSITE` object
1717: Level: advanced
1719: .seealso: `DMCOMPOSITE`, `DM`, `DMDestroy()`, `DMCompositeAddDM()`, `DMCompositeScatter()`, `DMCreate()`
1720: `DMCompositeGather()`, `DMCreateGlobalVector()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeGetAccess()`
1721: `DMCompositeGetLocalVectors()`, `DMCompositeRestoreLocalVectors()`, `DMCompositeGetEntries()`
1722: @*/
1723: PetscErrorCode DMCompositeCreate(MPI_Comm comm, DM *packer)
1724: {
1725: PetscFunctionBegin;
1726: PetscAssertPointer(packer, 2);
1727: PetscCall(DMCreate(comm, packer));
1728: PetscCall(DMSetType(*packer, DMCOMPOSITE));
1729: PetscFunctionReturn(PETSC_SUCCESS);
1730: }