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: /* ----------------------------------------------------------------------------------*/
118: /*@
119: DMCompositeGetNumberDM - Gets the number of `DM` objects in the `DMCOMPOSITE`
120: representation.
122: Not Collective
124: Input Parameter:
125: . dm - the `DMCOMPOSITE` object
127: Output Parameter:
128: . nDM - the number of `DM`
130: Level: beginner
132: .seealso: `DMCOMPOSITE`, `DM`
133: @*/
134: PetscErrorCode DMCompositeGetNumberDM(DM dm, PetscInt *nDM)
135: {
136: DM_Composite *com = (DM_Composite *)dm->data;
137: PetscBool flg;
139: PetscFunctionBegin;
141: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
142: PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
143: *nDM = com->nDM;
144: PetscFunctionReturn(PETSC_SUCCESS);
145: }
147: /*@C
148: DMCompositeGetAccess - Allows one to access the individual packed vectors in their global
149: representation.
151: Collective
153: Input Parameters:
154: + dm - the `DMCOMPOSITE` object
155: - gvec - the global vector
157: Output Parameter:
158: . ... - the packed parallel vectors, `NULL` for those that are not needed
160: Level: advanced
162: Note:
163: Use `DMCompositeRestoreAccess()` to return the vectors when you no longer need them
165: Fortran Notes:
166: Fortran callers must use numbered versions of this routine, e.g., DMCompositeGetAccess4(dm,gvec,vec1,vec2,vec3,vec4)
167: or use the alternative interface `DMCompositeGetAccessArray()`.
169: .seealso: `DMCOMPOSITE`, `DM`, `DMCompositeGetEntries()`, `DMCompositeScatter()`
170: @*/
171: PetscErrorCode DMCompositeGetAccess(DM dm, Vec gvec, ...)
172: {
173: va_list Argp;
174: struct DMCompositeLink *next;
175: DM_Composite *com = (DM_Composite *)dm->data;
176: PetscInt readonly;
177: PetscBool flg;
179: PetscFunctionBegin;
182: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
183: PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
184: next = com->next;
185: if (!com->setup) PetscCall(DMSetUp(dm));
187: PetscCall(VecLockGet(gvec, &readonly));
188: /* loop over packed objects, handling one at a time */
189: va_start(Argp, gvec);
190: while (next) {
191: Vec *vec;
192: vec = va_arg(Argp, Vec *);
193: if (vec) {
194: PetscCall(DMGetGlobalVector(next->dm, vec));
195: if (readonly) {
196: const PetscScalar *array;
197: PetscCall(VecGetArrayRead(gvec, &array));
198: PetscCall(VecPlaceArray(*vec, array + next->rstart));
199: PetscCall(VecLockReadPush(*vec));
200: PetscCall(VecRestoreArrayRead(gvec, &array));
201: } else {
202: PetscScalar *array;
203: PetscCall(VecGetArray(gvec, &array));
204: PetscCall(VecPlaceArray(*vec, array + next->rstart));
205: PetscCall(VecRestoreArray(gvec, &array));
206: }
207: }
208: next = next->next;
209: }
210: va_end(Argp);
211: PetscFunctionReturn(PETSC_SUCCESS);
212: }
214: /*@C
215: DMCompositeGetAccessArray - Allows one to access the individual packed vectors in their global
216: representation.
218: Collective
220: Input Parameters:
221: + dm - the `DMCOMPOSITE`
222: . pvec - packed vector
223: . nwanted - number of vectors wanted
224: - wanted - sorted array of vectors wanted, or `NULL` to get all vectors, length `nwanted`
226: Output Parameter:
227: . vecs - array of requested global vectors (must be previously allocated and of length `nwanted`)
229: Level: advanced
231: Note:
232: Use `DMCompositeRestoreAccessArray()` to return the vectors when you no longer need them
234: .seealso: `DMCOMPOSITE`, `DM`, `DMCompositeGetAccess()`, `DMCompositeGetEntries()`, `DMCompositeScatter()`, `DMCompositeGather()`
235: @*/
236: PetscErrorCode DMCompositeGetAccessArray(DM dm, Vec pvec, PetscInt nwanted, const PetscInt *wanted, Vec *vecs)
237: {
238: struct DMCompositeLink *link;
239: PetscInt i, wnum;
240: DM_Composite *com = (DM_Composite *)dm->data;
241: PetscInt readonly;
242: PetscBool flg;
244: PetscFunctionBegin;
247: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
248: PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
249: if (!com->setup) PetscCall(DMSetUp(dm));
251: PetscCall(VecLockGet(pvec, &readonly));
252: for (i = 0, wnum = 0, link = com->next; link && wnum < nwanted; i++, link = link->next) {
253: if (!wanted || i == wanted[wnum]) {
254: Vec v;
255: PetscCall(DMGetGlobalVector(link->dm, &v));
256: if (readonly) {
257: const PetscScalar *array;
258: PetscCall(VecGetArrayRead(pvec, &array));
259: PetscCall(VecPlaceArray(v, array + link->rstart));
260: PetscCall(VecLockReadPush(v));
261: PetscCall(VecRestoreArrayRead(pvec, &array));
262: } else {
263: PetscScalar *array;
264: PetscCall(VecGetArray(pvec, &array));
265: PetscCall(VecPlaceArray(v, array + link->rstart));
266: PetscCall(VecRestoreArray(pvec, &array));
267: }
268: vecs[wnum++] = v;
269: }
270: }
271: PetscFunctionReturn(PETSC_SUCCESS);
272: }
274: /*@C
275: DMCompositeGetLocalAccessArray - Allows one to access the individual
276: packed vectors in their local representation.
278: Collective
280: Input Parameters:
281: + dm - the `DMCOMPOSITE`
282: . pvec - packed vector
283: . nwanted - number of vectors wanted
284: - wanted - sorted array of vectors wanted, or `NULL` to get all vectors, length `nwanted`
286: Output Parameter:
287: . vecs - array of requested local vectors (must be allocated and of length `nwanted`)
289: Level: advanced
291: Note:
292: Use `DMCompositeRestoreLocalAccessArray()` to return the vectors
293: when you no longer need them.
295: .seealso: `DMCOMPOSITE`, `DM`, `DMCompositeRestoreLocalAccessArray()`, `DMCompositeGetAccess()`,
296: `DMCompositeGetEntries()`, `DMCompositeScatter()`, `DMCompositeGather()`
297: @*/
298: PetscErrorCode DMCompositeGetLocalAccessArray(DM dm, Vec pvec, PetscInt nwanted, const PetscInt *wanted, Vec *vecs)
299: {
300: struct DMCompositeLink *link;
301: PetscInt i, wnum;
302: DM_Composite *com = (DM_Composite *)dm->data;
303: PetscInt readonly;
304: PetscInt nlocal = 0;
305: PetscBool flg;
307: PetscFunctionBegin;
310: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
311: PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
312: if (!com->setup) PetscCall(DMSetUp(dm));
314: PetscCall(VecLockGet(pvec, &readonly));
315: for (i = 0, wnum = 0, link = com->next; link && wnum < nwanted; i++, link = link->next) {
316: if (!wanted || i == wanted[wnum]) {
317: Vec v;
318: PetscCall(DMGetLocalVector(link->dm, &v));
319: if (readonly) {
320: const PetscScalar *array;
321: PetscCall(VecGetArrayRead(pvec, &array));
322: PetscCall(VecPlaceArray(v, array + nlocal));
323: // 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
324: PetscCall(VecLockReadPush(v));
325: PetscCall(VecRestoreArrayRead(pvec, &array));
326: } else {
327: PetscScalar *array;
328: PetscCall(VecGetArray(pvec, &array));
329: PetscCall(VecPlaceArray(v, array + nlocal));
330: PetscCall(VecRestoreArray(pvec, &array));
331: }
332: vecs[wnum++] = v;
333: }
335: nlocal += link->nlocal;
336: }
337: PetscFunctionReturn(PETSC_SUCCESS);
338: }
340: /*@C
341: DMCompositeRestoreAccess - Returns the vectors obtained with `DMCompositeGetAccess()`
342: representation.
344: Collective
346: Input Parameters:
347: + dm - the `DMCOMPOSITE` object
348: . gvec - the global vector
349: - ... - the individual parallel vectors, `NULL` for those that are not needed
351: Level: advanced
353: .seealso: `DMCOMPOSITE`, `DM`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`,
354: `DMCompositeGather()`, `DMCompositeCreate()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeScatter()`,
355: `DMCompositeGetAccess()`
356: @*/
357: PetscErrorCode DMCompositeRestoreAccess(DM dm, Vec gvec, ...)
358: {
359: va_list Argp;
360: struct DMCompositeLink *next;
361: DM_Composite *com = (DM_Composite *)dm->data;
362: PetscInt readonly;
363: PetscBool flg;
365: PetscFunctionBegin;
368: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
369: PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
370: next = com->next;
371: if (!com->setup) PetscCall(DMSetUp(dm));
373: PetscCall(VecLockGet(gvec, &readonly));
374: /* loop over packed objects, handling one at a time */
375: va_start(Argp, gvec);
376: while (next) {
377: Vec *vec;
378: vec = va_arg(Argp, Vec *);
379: if (vec) {
380: PetscCall(VecResetArray(*vec));
381: if (readonly) PetscCall(VecLockReadPop(*vec));
382: PetscCall(DMRestoreGlobalVector(next->dm, vec));
383: }
384: next = next->next;
385: }
386: va_end(Argp);
387: PetscFunctionReturn(PETSC_SUCCESS);
388: }
390: /*@C
391: DMCompositeRestoreAccessArray - Returns the vectors obtained with `DMCompositeGetAccessArray()`
393: Collective
395: Input Parameters:
396: + dm - the `DMCOMPOSITE` object
397: . pvec - packed vector
398: . nwanted - number of vectors wanted
399: . wanted - sorted array of vectors wanted, or `NULL` to restore all vectors
400: - vecs - array of global vectors
402: Level: advanced
404: .seealso: `DMCOMPOSITE`, `DM`, `DMCompositeRestoreAccess()`, `DMCompositeRestoreEntries()`, `DMCompositeScatter()`, `DMCompositeGather()`
405: @*/
406: PetscErrorCode DMCompositeRestoreAccessArray(DM dm, Vec pvec, PetscInt nwanted, const PetscInt *wanted, Vec *vecs)
407: {
408: struct DMCompositeLink *link;
409: PetscInt i, wnum;
410: DM_Composite *com = (DM_Composite *)dm->data;
411: PetscInt readonly;
412: PetscBool flg;
414: PetscFunctionBegin;
417: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
418: PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
419: if (!com->setup) PetscCall(DMSetUp(dm));
421: PetscCall(VecLockGet(pvec, &readonly));
422: for (i = 0, wnum = 0, link = com->next; link && wnum < nwanted; i++, link = link->next) {
423: if (!wanted || i == wanted[wnum]) {
424: PetscCall(VecResetArray(vecs[wnum]));
425: if (readonly) PetscCall(VecLockReadPop(vecs[wnum]));
426: PetscCall(DMRestoreGlobalVector(link->dm, &vecs[wnum]));
427: wnum++;
428: }
429: }
430: PetscFunctionReturn(PETSC_SUCCESS);
431: }
433: /*@C
434: DMCompositeRestoreLocalAccessArray - Returns the vectors obtained with `DMCompositeGetLocalAccessArray()`.
436: Collective
438: Input Parameters:
439: + dm - the `DMCOMPOSITE` object
440: . pvec - packed vector
441: . nwanted - number of vectors wanted
442: . wanted - sorted array of vectors wanted, or `NULL` to restore all vectors
443: - vecs - array of local vectors
445: Level: advanced
447: Note:
448: `nwanted` and `wanted` must match the values given to `DMCompositeGetLocalAccessArray()`
449: otherwise the call will fail.
451: .seealso: `DMCOMPOSITE`, `DM`, `DMCompositeGetLocalAccessArray()`, `DMCompositeRestoreAccessArray()`,
452: `DMCompositeRestoreAccess()`, `DMCompositeRestoreEntries()`,
453: `DMCompositeScatter()`, `DMCompositeGather()`
454: @*/
455: PetscErrorCode DMCompositeRestoreLocalAccessArray(DM dm, Vec pvec, PetscInt nwanted, const PetscInt *wanted, Vec *vecs)
456: {
457: struct DMCompositeLink *link;
458: PetscInt i, wnum;
459: DM_Composite *com = (DM_Composite *)dm->data;
460: PetscInt readonly;
461: PetscBool flg;
463: PetscFunctionBegin;
466: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
467: PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
468: if (!com->setup) PetscCall(DMSetUp(dm));
470: PetscCall(VecLockGet(pvec, &readonly));
471: for (i = 0, wnum = 0, link = com->next; link && wnum < nwanted; i++, link = link->next) {
472: if (!wanted || i == wanted[wnum]) {
473: PetscCall(VecResetArray(vecs[wnum]));
474: if (readonly) PetscCall(VecLockReadPop(vecs[wnum]));
475: PetscCall(DMRestoreLocalVector(link->dm, &vecs[wnum]));
476: wnum++;
477: }
478: }
479: PetscFunctionReturn(PETSC_SUCCESS);
480: }
482: /*@C
483: DMCompositeScatter - Scatters from a global packed vector into its individual local vectors
485: Collective
487: Input Parameters:
488: + dm - the `DMCOMPOSITE` object
489: . gvec - the global vector
490: - ... - the individual sequential vectors, `NULL` for those that are not needed
492: Level: advanced
494: Note:
495: `DMCompositeScatterArray()` is a non-variadic alternative that is often more convenient for library callers and is
496: accessible from Fortran.
498: .seealso: `DMCOMPOSITE`, `DM`, `DMDestroy()`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`,
499: `DMCompositeGather()`, `DMCompositeCreate()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeGetAccess()`,
500: `DMCompositeGetLocalVectors()`, `DMCompositeRestoreLocalVectors()`, `DMCompositeGetEntries()`
501: `DMCompositeScatterArray()`
502: @*/
503: PetscErrorCode DMCompositeScatter(DM dm, Vec gvec, ...)
504: {
505: va_list Argp;
506: struct DMCompositeLink *next;
507: PETSC_UNUSED PetscInt cnt;
508: DM_Composite *com = (DM_Composite *)dm->data;
509: PetscBool flg;
511: PetscFunctionBegin;
514: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
515: PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
516: if (!com->setup) PetscCall(DMSetUp(dm));
518: /* loop over packed objects, handling one at a time */
519: va_start(Argp, gvec);
520: for (cnt = 3, next = com->next; next; cnt++, next = next->next) {
521: Vec local;
522: local = va_arg(Argp, Vec);
523: if (local) {
524: Vec global;
525: const PetscScalar *array;
527: PetscCall(DMGetGlobalVector(next->dm, &global));
528: PetscCall(VecGetArrayRead(gvec, &array));
529: PetscCall(VecPlaceArray(global, array + next->rstart));
530: PetscCall(DMGlobalToLocalBegin(next->dm, global, INSERT_VALUES, local));
531: PetscCall(DMGlobalToLocalEnd(next->dm, global, INSERT_VALUES, local));
532: PetscCall(VecRestoreArrayRead(gvec, &array));
533: PetscCall(VecResetArray(global));
534: PetscCall(DMRestoreGlobalVector(next->dm, &global));
535: }
536: }
537: va_end(Argp);
538: PetscFunctionReturn(PETSC_SUCCESS);
539: }
541: /*@
542: DMCompositeScatterArray - Scatters from a global packed vector into its individual local vectors
544: Collective
546: Input Parameters:
547: + dm - the `DMCOMPOSITE` object
548: . gvec - the global vector
549: - lvecs - array of local vectors, NULL for any that are not needed
551: Level: advanced
553: Note:
554: This is a non-variadic alternative to `DMCompositeScatter()`
556: .seealso: `DMCOMPOSITE`, `DM`, `DMDestroy()`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`
557: `DMCompositeGather()`, `DMCompositeCreate()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeGetAccess()`,
558: `DMCompositeGetLocalVectors()`, `DMCompositeRestoreLocalVectors()`, `DMCompositeGetEntries()`
559: @*/
560: PetscErrorCode DMCompositeScatterArray(DM dm, Vec gvec, Vec *lvecs)
561: {
562: struct DMCompositeLink *next;
563: PetscInt i;
564: DM_Composite *com = (DM_Composite *)dm->data;
565: PetscBool flg;
567: PetscFunctionBegin;
570: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
571: PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
572: if (!com->setup) PetscCall(DMSetUp(dm));
574: /* loop over packed objects, handling one at a time */
575: for (i = 0, next = com->next; next; next = next->next, i++) {
576: if (lvecs[i]) {
577: Vec global;
578: const PetscScalar *array;
580: PetscCall(DMGetGlobalVector(next->dm, &global));
581: PetscCall(VecGetArrayRead(gvec, &array));
582: PetscCall(VecPlaceArray(global, (PetscScalar *)array + next->rstart));
583: PetscCall(DMGlobalToLocalBegin(next->dm, global, INSERT_VALUES, lvecs[i]));
584: PetscCall(DMGlobalToLocalEnd(next->dm, global, INSERT_VALUES, lvecs[i]));
585: PetscCall(VecRestoreArrayRead(gvec, &array));
586: PetscCall(VecResetArray(global));
587: PetscCall(DMRestoreGlobalVector(next->dm, &global));
588: }
589: }
590: PetscFunctionReturn(PETSC_SUCCESS);
591: }
593: /*@C
594: DMCompositeGather - Gathers into a global packed vector from its individual local vectors
596: Collective
598: Input Parameters:
599: + dm - the `DMCOMPOSITE` object
600: . imode - `INSERT_VALUES` or `ADD_VALUES`
601: . gvec - the global vector
602: - ... - the individual sequential vectors, `NULL` for any that are not needed
604: Level: advanced
606: Fortran Notes:
607: Fortran users should use `DMCompositeGatherArray()`
609: .seealso: `DMCOMPOSITE`, `DM`, `DMDestroy()`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`,
610: `DMCompositeScatter()`, `DMCompositeCreate()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeGetAccess()`,
611: `DMCompositeGetLocalVectors()`, `DMCompositeRestoreLocalVectors()`, `DMCompositeGetEntries()`
612: @*/
613: PetscErrorCode DMCompositeGather(DM dm, InsertMode imode, Vec gvec, ...)
614: {
615: va_list Argp;
616: struct DMCompositeLink *next;
617: DM_Composite *com = (DM_Composite *)dm->data;
618: PETSC_UNUSED PetscInt cnt;
619: PetscBool flg;
621: PetscFunctionBegin;
624: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
625: PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
626: if (!com->setup) PetscCall(DMSetUp(dm));
628: /* loop over packed objects, handling one at a time */
629: va_start(Argp, gvec);
630: for (cnt = 3, next = com->next; next; cnt++, next = next->next) {
631: Vec local;
632: local = va_arg(Argp, Vec);
633: if (local) {
634: PetscScalar *array;
635: Vec global;
637: PetscCall(DMGetGlobalVector(next->dm, &global));
638: PetscCall(VecGetArray(gvec, &array));
639: PetscCall(VecPlaceArray(global, array + next->rstart));
640: PetscCall(DMLocalToGlobalBegin(next->dm, local, imode, global));
641: PetscCall(DMLocalToGlobalEnd(next->dm, local, imode, global));
642: PetscCall(VecRestoreArray(gvec, &array));
643: PetscCall(VecResetArray(global));
644: PetscCall(DMRestoreGlobalVector(next->dm, &global));
645: }
646: }
647: va_end(Argp);
648: PetscFunctionReturn(PETSC_SUCCESS);
649: }
651: /*@
652: DMCompositeGatherArray - Gathers into a global packed vector from its individual local vectors
654: Collective
656: Input Parameters:
657: + dm - the `DMCOMPOSITE` object
658: . gvec - the global vector
659: . imode - `INSERT_VALUES` or `ADD_VALUES`
660: - lvecs - the individual sequential vectors, NULL for any that are not needed
662: Level: advanced
664: Note:
665: This is a non-variadic alternative to `DMCompositeGather()`.
667: .seealso: `DMCOMPOSITE`, `DM`, `DMDestroy()`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`,
668: `DMCompositeScatter()`, `DMCompositeCreate()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeGetAccess()`,
669: `DMCompositeGetLocalVectors()`, `DMCompositeRestoreLocalVectors()`, `DMCompositeGetEntries()`,
670: @*/
671: PetscErrorCode DMCompositeGatherArray(DM dm, InsertMode imode, Vec gvec, Vec *lvecs)
672: {
673: struct DMCompositeLink *next;
674: DM_Composite *com = (DM_Composite *)dm->data;
675: PetscInt i;
676: PetscBool flg;
678: PetscFunctionBegin;
681: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
682: PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
683: if (!com->setup) PetscCall(DMSetUp(dm));
685: /* loop over packed objects, handling one at a time */
686: for (next = com->next, i = 0; next; next = next->next, i++) {
687: if (lvecs[i]) {
688: PetscScalar *array;
689: Vec global;
691: PetscCall(DMGetGlobalVector(next->dm, &global));
692: PetscCall(VecGetArray(gvec, &array));
693: PetscCall(VecPlaceArray(global, array + next->rstart));
694: PetscCall(DMLocalToGlobalBegin(next->dm, lvecs[i], imode, global));
695: PetscCall(DMLocalToGlobalEnd(next->dm, lvecs[i], imode, global));
696: PetscCall(VecRestoreArray(gvec, &array));
697: PetscCall(VecResetArray(global));
698: PetscCall(DMRestoreGlobalVector(next->dm, &global));
699: }
700: }
701: PetscFunctionReturn(PETSC_SUCCESS);
702: }
704: /*@
705: DMCompositeAddDM - adds a `DM` vector to a `DMCOMPOSITE`
707: Collective
709: Input Parameters:
710: + dmc - the `DMCOMPOSITE` object
711: - dm - the `DM` object
713: Level: advanced
715: .seealso: `DMCOMPOSITE`, `DM`, `DMDestroy()`, `DMCompositeGather()`, `DMCreateGlobalVector()`,
716: `DMCompositeScatter()`, `DMCompositeCreate()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeGetAccess()`,
717: `DMCompositeGetLocalVectors()`, `DMCompositeRestoreLocalVectors()`, `DMCompositeGetEntries()`
718: @*/
719: PetscErrorCode DMCompositeAddDM(DM dmc, DM dm)
720: {
721: PetscInt n, nlocal;
722: struct DMCompositeLink *mine, *next;
723: Vec global, local;
724: DM_Composite *com = (DM_Composite *)dmc->data;
725: PetscBool flg;
727: PetscFunctionBegin;
730: PetscCall(PetscObjectTypeCompare((PetscObject)dmc, DMCOMPOSITE, &flg));
731: PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
732: next = com->next;
733: PetscCheck(!com->setup, PetscObjectComm((PetscObject)dmc), PETSC_ERR_ARG_WRONGSTATE, "Cannot add a DM once you have used the DMComposite");
735: /* create new link */
736: PetscCall(PetscNew(&mine));
737: PetscCall(PetscObjectReference((PetscObject)dm));
738: PetscCall(DMGetGlobalVector(dm, &global));
739: PetscCall(VecGetLocalSize(global, &n));
740: PetscCall(DMRestoreGlobalVector(dm, &global));
741: PetscCall(DMGetLocalVector(dm, &local));
742: PetscCall(VecGetSize(local, &nlocal));
743: PetscCall(DMRestoreLocalVector(dm, &local));
745: mine->n = n;
746: mine->nlocal = nlocal;
747: mine->dm = dm;
748: mine->next = NULL;
749: com->n += n;
750: com->nghost += nlocal;
752: /* add to end of list */
753: if (!next) com->next = mine;
754: else {
755: while (next->next) next = next->next;
756: next->next = mine;
757: }
758: com->nDM++;
759: com->nmine++;
760: PetscFunctionReturn(PETSC_SUCCESS);
761: }
763: #include <petscdraw.h>
764: PETSC_EXTERN PetscErrorCode VecView_MPI(Vec, PetscViewer);
765: static PetscErrorCode VecView_DMComposite(Vec gvec, PetscViewer viewer)
766: {
767: DM dm;
768: struct DMCompositeLink *next;
769: PetscBool isdraw;
770: DM_Composite *com;
772: PetscFunctionBegin;
773: PetscCall(VecGetDM(gvec, &dm));
774: PetscCheck(dm, PetscObjectComm((PetscObject)gvec), PETSC_ERR_ARG_WRONG, "Vector not generated from a DMComposite");
775: com = (DM_Composite *)dm->data;
776: next = com->next;
778: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw));
779: if (!isdraw) {
780: /* do I really want to call this? */
781: PetscCall(VecView_MPI(gvec, viewer));
782: } else {
783: PetscInt cnt = 0;
785: /* loop over packed objects, handling one at a time */
786: while (next) {
787: Vec vec;
788: const PetscScalar *array;
789: PetscInt bs;
791: /* Should use VecGetSubVector() eventually, but would need to forward the DM for that to work */
792: PetscCall(DMGetGlobalVector(next->dm, &vec));
793: PetscCall(VecGetArrayRead(gvec, &array));
794: PetscCall(VecPlaceArray(vec, (PetscScalar *)array + next->rstart));
795: PetscCall(VecRestoreArrayRead(gvec, &array));
796: PetscCall(VecView(vec, viewer));
797: PetscCall(VecResetArray(vec));
798: PetscCall(VecGetBlockSize(vec, &bs));
799: PetscCall(DMRestoreGlobalVector(next->dm, &vec));
800: PetscCall(PetscViewerDrawBaseAdd(viewer, bs));
801: cnt += bs;
802: next = next->next;
803: }
804: PetscCall(PetscViewerDrawBaseAdd(viewer, -cnt));
805: }
806: PetscFunctionReturn(PETSC_SUCCESS);
807: }
809: static PetscErrorCode DMCreateGlobalVector_Composite(DM dm, Vec *gvec)
810: {
811: DM_Composite *com = (DM_Composite *)dm->data;
813: PetscFunctionBegin;
815: PetscCall(DMSetFromOptions(dm));
816: PetscCall(DMSetUp(dm));
817: PetscCall(VecCreate(PetscObjectComm((PetscObject)dm), gvec));
818: PetscCall(VecSetType(*gvec, dm->vectype));
819: PetscCall(VecSetSizes(*gvec, com->n, com->N));
820: PetscCall(VecSetDM(*gvec, dm));
821: PetscCall(VecSetOperation(*gvec, VECOP_VIEW, (void (*)(void))VecView_DMComposite));
822: PetscFunctionReturn(PETSC_SUCCESS);
823: }
825: static PetscErrorCode DMCreateLocalVector_Composite(DM dm, Vec *lvec)
826: {
827: DM_Composite *com = (DM_Composite *)dm->data;
829: PetscFunctionBegin;
831: if (!com->setup) {
832: PetscCall(DMSetFromOptions(dm));
833: PetscCall(DMSetUp(dm));
834: }
835: PetscCall(VecCreate(PETSC_COMM_SELF, lvec));
836: PetscCall(VecSetType(*lvec, dm->vectype));
837: PetscCall(VecSetSizes(*lvec, com->nghost, PETSC_DECIDE));
838: PetscCall(VecSetDM(*lvec, dm));
839: PetscFunctionReturn(PETSC_SUCCESS);
840: }
842: /*@C
843: DMCompositeGetISLocalToGlobalMappings - gets an `ISLocalToGlobalMapping` for each `DM` in the `DMCOMPOSITE`, maps to the composite global space
845: Collective; No Fortran Support
847: Input Parameter:
848: . dm - the `DMCOMPOSITE` object
850: Output Parameter:
851: . ltogs - the individual mappings for each packed vector. Note that this includes
852: all the ghost points that individual ghosted `DMDA` may have.
854: Level: advanced
856: Note:
857: Each entry of `ltogs` should be destroyed with `ISLocalToGlobalMappingDestroy()`, `ltogs` should be freed with `PetscFree()`.
859: .seealso: `DMCOMPOSITE`, `DM`, `DMDestroy()`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`,
860: `DMCompositeGather()`, `DMCompositeCreate()`, `DMCompositeGetAccess()`, `DMCompositeScatter()`,
861: `DMCompositeGetLocalVectors()`, `DMCompositeRestoreLocalVectors()`, `DMCompositeGetEntries()`
862: @*/
863: PetscErrorCode DMCompositeGetISLocalToGlobalMappings(DM dm, ISLocalToGlobalMapping **ltogs)
864: {
865: PetscInt i, *idx, n, cnt;
866: struct DMCompositeLink *next;
867: PetscMPIInt rank;
868: DM_Composite *com = (DM_Composite *)dm->data;
869: PetscBool flg;
871: PetscFunctionBegin;
873: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
874: PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
875: PetscCall(DMSetUp(dm));
876: PetscCall(PetscMalloc1(com->nDM, ltogs));
877: next = com->next;
878: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
880: /* loop over packed objects, handling one at a time */
881: cnt = 0;
882: while (next) {
883: ISLocalToGlobalMapping ltog;
884: PetscMPIInt size;
885: const PetscInt *suboff, *indices;
886: Vec global;
888: /* Get sub-DM global indices for each local dof */
889: PetscCall(DMGetLocalToGlobalMapping(next->dm, <og));
890: PetscCall(ISLocalToGlobalMappingGetSize(ltog, &n));
891: PetscCall(ISLocalToGlobalMappingGetIndices(ltog, &indices));
892: PetscCall(PetscMalloc1(n, &idx));
894: /* Get the offsets for the sub-DM global vector */
895: PetscCall(DMGetGlobalVector(next->dm, &global));
896: PetscCall(VecGetOwnershipRanges(global, &suboff));
897: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)global), &size));
899: /* Shift the sub-DM definition of the global space to the composite global space */
900: for (i = 0; i < n; i++) {
901: PetscInt subi = indices[i], lo = 0, hi = size, t;
902: /* There's no consensus on what a negative index means,
903: except for skipping when setting the values in vectors and matrices */
904: if (subi < 0) {
905: idx[i] = subi - next->grstarts[rank];
906: continue;
907: }
908: /* Binary search to find which rank owns subi */
909: while (hi - lo > 1) {
910: t = lo + (hi - lo) / 2;
911: if (suboff[t] > subi) hi = t;
912: else lo = t;
913: }
914: idx[i] = subi - suboff[lo] + next->grstarts[lo];
915: }
916: PetscCall(ISLocalToGlobalMappingRestoreIndices(ltog, &indices));
917: PetscCall(ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)dm), 1, n, idx, PETSC_OWN_POINTER, &(*ltogs)[cnt]));
918: PetscCall(DMRestoreGlobalVector(next->dm, &global));
919: next = next->next;
920: cnt++;
921: }
922: PetscFunctionReturn(PETSC_SUCCESS);
923: }
925: /*@C
926: DMCompositeGetLocalISs - Gets index sets for each component of a composite local vector
928: Not Collective; No Fortran Support
930: Input Parameter:
931: . dm - the `DMCOMPOSITE`
933: Output Parameter:
934: . is - array of serial index sets for each component of the `DMCOMPOSITE`
936: Level: intermediate
938: Notes:
939: At present, a composite local vector does not normally exist. This function is used to provide index sets for
940: `MatGetLocalSubMatrix()`. In the future, the scatters for each entry in the `DMCOMPOSITE` may be merged into a single
941: scatter to a composite local vector. The user should not typically need to know which is being done.
943: To get the composite global indices at all local points (including ghosts), use `DMCompositeGetISLocalToGlobalMappings()`.
945: To get index sets for pieces of the composite global vector, use `DMCompositeGetGlobalISs()`.
947: Each returned `IS` should be destroyed with `ISDestroy()`, the array should be freed with `PetscFree()`.
949: Fortran Note:
950: Pass in an array long enough to hold all the `IS`, see `DMCompositeGetNumberDM()`
952: .seealso: `DMCOMPOSITE`, `DM`, `DMCompositeGetGlobalISs()`, `DMCompositeGetISLocalToGlobalMappings()`, `MatGetLocalSubMatrix()`,
953: `MatCreateLocalRef()`, `DMCompositeGetNumberDM()`
954: @*/
955: PetscErrorCode DMCompositeGetLocalISs(DM dm, IS **is)
956: {
957: DM_Composite *com = (DM_Composite *)dm->data;
958: struct DMCompositeLink *link;
959: PetscInt cnt, start;
960: PetscBool flg;
962: PetscFunctionBegin;
964: PetscAssertPointer(is, 2);
965: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
966: PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
967: PetscCall(PetscMalloc1(com->nmine, is));
968: for (cnt = 0, start = 0, link = com->next; link; start += link->nlocal, cnt++, link = link->next) {
969: PetscInt bs;
970: PetscCall(ISCreateStride(PETSC_COMM_SELF, link->nlocal, start, 1, &(*is)[cnt]));
971: PetscCall(DMGetBlockSize(link->dm, &bs));
972: PetscCall(ISSetBlockSize((*is)[cnt], bs));
973: }
974: PetscFunctionReturn(PETSC_SUCCESS);
975: }
977: /*@C
978: DMCompositeGetGlobalISs - Gets the index sets for each composed object in a `DMCOMPOSITE`
980: Collective
982: Input Parameter:
983: . dm - the `DMCOMPOSITE` object
985: Output Parameter:
986: . is - the array of index sets
988: Level: advanced
990: Notes:
991: The `is` entries should be destroyed with `ISDestroy()`, `is` should be freed with `PetscFree()`
993: These could be used to extract a subset of vector entries for a "multi-physics" preconditioner
995: Use `DMCompositeGetLocalISs()` for index sets in the packed local numbering, and
996: `DMCompositeGetISLocalToGlobalMappings()` for to map local sub-`DM` (including ghost) indices to packed global
997: indices.
999: Fortran Notes:
1000: The output argument 'is' must be an allocated array of sufficient length, which can be learned using `DMCompositeGetNumberDM()`.
1002: .seealso: `DMCOMPOSITE`, `DM`, `DMDestroy()`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`,
1003: `DMCompositeGather()`, `DMCompositeCreate()`, `DMCompositeGetAccess()`, `DMCompositeScatter()`,
1004: `DMCompositeGetLocalVectors()`, `DMCompositeRestoreLocalVectors()`, `DMCompositeGetEntries()`
1005: @*/
1006: PetscErrorCode DMCompositeGetGlobalISs(DM dm, IS *is[])
1007: {
1008: PetscInt cnt = 0;
1009: struct DMCompositeLink *next;
1010: PetscMPIInt rank;
1011: DM_Composite *com = (DM_Composite *)dm->data;
1012: PetscBool flg;
1014: PetscFunctionBegin;
1016: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
1017: PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
1018: PetscCheck(dm->setupcalled, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "Must call DMSetUp() before");
1019: PetscCall(PetscMalloc1(com->nDM, is));
1020: next = com->next;
1021: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
1023: /* loop over packed objects, handling one at a time */
1024: while (next) {
1025: PetscDS prob;
1027: PetscCall(ISCreateStride(PetscObjectComm((PetscObject)dm), next->n, next->grstart, 1, &(*is)[cnt]));
1028: PetscCall(DMGetDS(dm, &prob));
1029: if (prob) {
1030: MatNullSpace space;
1031: Mat pmat;
1032: PetscObject disc;
1033: PetscInt Nf;
1035: PetscCall(PetscDSGetNumFields(prob, &Nf));
1036: if (cnt < Nf) {
1037: PetscCall(PetscDSGetDiscretization(prob, cnt, &disc));
1038: PetscCall(PetscObjectQuery(disc, "nullspace", (PetscObject *)&space));
1039: if (space) PetscCall(PetscObjectCompose((PetscObject)(*is)[cnt], "nullspace", (PetscObject)space));
1040: PetscCall(PetscObjectQuery(disc, "nearnullspace", (PetscObject *)&space));
1041: if (space) PetscCall(PetscObjectCompose((PetscObject)(*is)[cnt], "nearnullspace", (PetscObject)space));
1042: PetscCall(PetscObjectQuery(disc, "pmat", (PetscObject *)&pmat));
1043: if (pmat) PetscCall(PetscObjectCompose((PetscObject)(*is)[cnt], "pmat", (PetscObject)pmat));
1044: }
1045: }
1046: cnt++;
1047: next = next->next;
1048: }
1049: PetscFunctionReturn(PETSC_SUCCESS);
1050: }
1052: static PetscErrorCode DMCreateFieldIS_Composite(DM dm, PetscInt *numFields, char ***fieldNames, IS **fields)
1053: {
1054: PetscInt nDM;
1055: DM *dms;
1056: PetscInt i;
1058: PetscFunctionBegin;
1059: PetscCall(DMCompositeGetNumberDM(dm, &nDM));
1060: if (numFields) *numFields = nDM;
1061: PetscCall(DMCompositeGetGlobalISs(dm, fields));
1062: if (fieldNames) {
1063: PetscCall(PetscMalloc1(nDM, &dms));
1064: PetscCall(PetscMalloc1(nDM, fieldNames));
1065: PetscCall(DMCompositeGetEntriesArray(dm, dms));
1066: for (i = 0; i < nDM; i++) {
1067: char buf[256];
1068: const char *splitname;
1070: /* Split naming precedence: object name, prefix, number */
1071: splitname = ((PetscObject)dm)->name;
1072: if (!splitname) {
1073: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)dms[i], &splitname));
1074: if (splitname) {
1075: size_t len;
1076: PetscCall(PetscStrncpy(buf, splitname, sizeof(buf)));
1077: buf[sizeof(buf) - 1] = 0;
1078: PetscCall(PetscStrlen(buf, &len));
1079: if (buf[len - 1] == '_') buf[len - 1] = 0; /* Remove trailing underscore if it was used */
1080: splitname = buf;
1081: }
1082: }
1083: if (!splitname) {
1084: PetscCall(PetscSNPrintf(buf, sizeof(buf), "%" PetscInt_FMT, i));
1085: splitname = buf;
1086: }
1087: PetscCall(PetscStrallocpy(splitname, &(*fieldNames)[i]));
1088: }
1089: PetscCall(PetscFree(dms));
1090: }
1091: PetscFunctionReturn(PETSC_SUCCESS);
1092: }
1094: /*
1095: This could take over from DMCreateFieldIS(), as it is more general,
1096: making DMCreateFieldIS() a special case -- calling with dmlist == NULL;
1097: At this point it's probably best to be less intrusive, however.
1098: */
1099: static PetscErrorCode DMCreateFieldDecomposition_Composite(DM dm, PetscInt *len, char ***namelist, IS **islist, DM **dmlist)
1100: {
1101: PetscInt nDM;
1102: PetscInt i;
1104: PetscFunctionBegin;
1105: PetscCall(DMCreateFieldIS_Composite(dm, len, namelist, islist));
1106: if (dmlist) {
1107: PetscCall(DMCompositeGetNumberDM(dm, &nDM));
1108: PetscCall(PetscMalloc1(nDM, dmlist));
1109: PetscCall(DMCompositeGetEntriesArray(dm, *dmlist));
1110: for (i = 0; i < nDM; i++) PetscCall(PetscObjectReference((PetscObject)((*dmlist)[i])));
1111: }
1112: PetscFunctionReturn(PETSC_SUCCESS);
1113: }
1115: /* -------------------------------------------------------------------------------------*/
1116: /*@C
1117: DMCompositeGetLocalVectors - Gets local vectors for each part of a `DMCOMPOSITE`
1118: Use `DMCompositeRestoreLocalVectors()` to return them.
1120: Not Collective; No Fortran Support
1122: Input Parameter:
1123: . dm - the `DMCOMPOSITE` object
1125: Output Parameter:
1126: . ... - the individual sequential `Vec`s
1128: Level: advanced
1130: .seealso: `DMCOMPOSITE`, `DM`, `DMDestroy()`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`,
1131: `DMCompositeGather()`, `DMCompositeCreate()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeGetAccess()`,
1132: `DMCompositeRestoreLocalVectors()`, `DMCompositeScatter()`, `DMCompositeGetEntries()`
1133: @*/
1134: PetscErrorCode DMCompositeGetLocalVectors(DM dm, ...)
1135: {
1136: va_list Argp;
1137: struct DMCompositeLink *next;
1138: DM_Composite *com = (DM_Composite *)dm->data;
1139: PetscBool flg;
1141: PetscFunctionBegin;
1143: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
1144: PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
1145: next = com->next;
1146: /* loop over packed objects, handling one at a time */
1147: va_start(Argp, dm);
1148: while (next) {
1149: Vec *vec;
1150: vec = va_arg(Argp, Vec *);
1151: if (vec) PetscCall(DMGetLocalVector(next->dm, vec));
1152: next = next->next;
1153: }
1154: va_end(Argp);
1155: PetscFunctionReturn(PETSC_SUCCESS);
1156: }
1158: /*@C
1159: DMCompositeRestoreLocalVectors - Restores local vectors for each part of a `DMCOMPOSITE`
1161: Not Collective; No Fortran Support
1163: Input Parameter:
1164: . dm - the `DMCOMPOSITE` object
1166: Output Parameter:
1167: . ... - the individual sequential `Vec`
1169: Level: advanced
1171: .seealso: `DMCOMPOSITE`, `DM`, `DMDestroy()`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`,
1172: `DMCompositeGather()`, `DMCompositeCreate()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeGetAccess()`,
1173: `DMCompositeGetLocalVectors()`, `DMCompositeScatter()`, `DMCompositeGetEntries()`
1174: @*/
1175: PetscErrorCode DMCompositeRestoreLocalVectors(DM dm, ...)
1176: {
1177: va_list Argp;
1178: struct DMCompositeLink *next;
1179: DM_Composite *com = (DM_Composite *)dm->data;
1180: PetscBool flg;
1182: PetscFunctionBegin;
1184: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
1185: PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
1186: next = com->next;
1187: /* loop over packed objects, handling one at a time */
1188: va_start(Argp, dm);
1189: while (next) {
1190: Vec *vec;
1191: vec = va_arg(Argp, Vec *);
1192: if (vec) PetscCall(DMRestoreLocalVector(next->dm, vec));
1193: next = next->next;
1194: }
1195: va_end(Argp);
1196: PetscFunctionReturn(PETSC_SUCCESS);
1197: }
1199: /* -------------------------------------------------------------------------------------*/
1200: /*@C
1201: DMCompositeGetEntries - Gets the `DM` for each entry in a `DMCOMPOSITE`.
1203: Not Collective
1205: Input Parameter:
1206: . dm - the `DMCOMPOSITE` object
1208: Output Parameter:
1209: . ... - the individual entries `DM`
1211: Level: advanced
1213: Fortran Notes:
1214: Available as `DMCompositeGetEntries()` for one output `DM`, DMCompositeGetEntries2() for 2, etc
1216: .seealso: `DMCOMPOSITE`, `DM`, `DMDestroy()`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`, `DMCompositeGetEntriesArray()`
1217: `DMCompositeGather()`, `DMCompositeCreate()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeGetAccess()`,
1218: `DMCompositeRestoreLocalVectors()`, `DMCompositeGetLocalVectors()`, `DMCompositeScatter()`
1219: @*/
1220: PetscErrorCode DMCompositeGetEntries(DM dm, ...)
1221: {
1222: va_list Argp;
1223: struct DMCompositeLink *next;
1224: DM_Composite *com = (DM_Composite *)dm->data;
1225: PetscBool flg;
1227: PetscFunctionBegin;
1229: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
1230: PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
1231: next = com->next;
1232: /* loop over packed objects, handling one at a time */
1233: va_start(Argp, dm);
1234: while (next) {
1235: DM *dmn;
1236: dmn = va_arg(Argp, DM *);
1237: if (dmn) *dmn = next->dm;
1238: next = next->next;
1239: }
1240: va_end(Argp);
1241: PetscFunctionReturn(PETSC_SUCCESS);
1242: }
1244: /*@C
1245: DMCompositeGetEntriesArray - Gets the DM for each entry in a `DMCOMPOSITE`
1247: Not Collective
1249: Input Parameter:
1250: . dm - the `DMCOMPOSITE` object
1252: Output Parameter:
1253: . dms - array of sufficient length (see `DMCompositeGetNumberDM()`) to hold the individual `DM`
1255: Level: advanced
1257: .seealso: `DMCOMPOSITE`, `DM`, `DMDestroy()`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`, `DMCompositeGetEntries()`
1258: `DMCompositeGather()`, `DMCompositeCreate()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeGetAccess()`,
1259: `DMCompositeRestoreLocalVectors()`, `DMCompositeGetLocalVectors()`, `DMCompositeScatter()`
1260: @*/
1261: PetscErrorCode DMCompositeGetEntriesArray(DM dm, DM dms[])
1262: {
1263: struct DMCompositeLink *next;
1264: DM_Composite *com = (DM_Composite *)dm->data;
1265: PetscInt i;
1266: PetscBool flg;
1268: PetscFunctionBegin;
1270: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
1271: PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
1272: /* loop over packed objects, handling one at a time */
1273: for (next = com->next, i = 0; next; next = next->next, i++) dms[i] = next->dm;
1274: PetscFunctionReturn(PETSC_SUCCESS);
1275: }
1277: typedef struct {
1278: DM dm;
1279: PetscViewer *subv;
1280: Vec *vecs;
1281: } GLVisViewerCtx;
1283: static PetscErrorCode DestroyGLVisViewerCtx_Private(void *vctx)
1284: {
1285: GLVisViewerCtx *ctx = (GLVisViewerCtx *)vctx;
1286: PetscInt i, n;
1288: PetscFunctionBegin;
1289: PetscCall(DMCompositeGetNumberDM(ctx->dm, &n));
1290: for (i = 0; i < n; i++) PetscCall(PetscViewerDestroy(&ctx->subv[i]));
1291: PetscCall(PetscFree2(ctx->subv, ctx->vecs));
1292: PetscCall(DMDestroy(&ctx->dm));
1293: PetscCall(PetscFree(ctx));
1294: PetscFunctionReturn(PETSC_SUCCESS);
1295: }
1297: static PetscErrorCode DMCompositeSampleGLVisFields_Private(PetscObject oX, PetscInt nf, PetscObject oXfield[], void *vctx)
1298: {
1299: Vec X = (Vec)oX;
1300: GLVisViewerCtx *ctx = (GLVisViewerCtx *)vctx;
1301: PetscInt i, n, cumf;
1303: PetscFunctionBegin;
1304: PetscCall(DMCompositeGetNumberDM(ctx->dm, &n));
1305: PetscCall(DMCompositeGetAccessArray(ctx->dm, X, n, NULL, ctx->vecs));
1306: for (i = 0, cumf = 0; i < n; i++) {
1307: PetscErrorCode (*g2l)(PetscObject, PetscInt, PetscObject[], void *);
1308: void *fctx;
1309: PetscInt nfi;
1311: PetscCall(PetscViewerGLVisGetFields_Internal(ctx->subv[i], &nfi, NULL, NULL, &g2l, NULL, &fctx));
1312: if (!nfi) continue;
1313: if (g2l) PetscCall((*g2l)((PetscObject)ctx->vecs[i], nfi, oXfield + cumf, fctx));
1314: else PetscCall(VecCopy(ctx->vecs[i], (Vec)oXfield[cumf]));
1315: cumf += nfi;
1316: }
1317: PetscCall(DMCompositeRestoreAccessArray(ctx->dm, X, n, NULL, ctx->vecs));
1318: PetscFunctionReturn(PETSC_SUCCESS);
1319: }
1321: static PetscErrorCode DMSetUpGLVisViewer_Composite(PetscObject odm, PetscViewer viewer)
1322: {
1323: DM dm = (DM)odm, *dms;
1324: Vec *Ufds;
1325: GLVisViewerCtx *ctx;
1326: PetscInt i, n, tnf, *sdim;
1327: char **fecs;
1329: PetscFunctionBegin;
1330: PetscCall(PetscNew(&ctx));
1331: PetscCall(PetscObjectReference((PetscObject)dm));
1332: ctx->dm = dm;
1333: PetscCall(DMCompositeGetNumberDM(dm, &n));
1334: PetscCall(PetscMalloc1(n, &dms));
1335: PetscCall(DMCompositeGetEntriesArray(dm, dms));
1336: PetscCall(PetscMalloc2(n, &ctx->subv, n, &ctx->vecs));
1337: for (i = 0, tnf = 0; i < n; i++) {
1338: PetscInt nf;
1340: PetscCall(PetscViewerCreate(PetscObjectComm(odm), &ctx->subv[i]));
1341: PetscCall(PetscViewerSetType(ctx->subv[i], PETSCVIEWERGLVIS));
1342: PetscCall(PetscViewerGLVisSetDM_Internal(ctx->subv[i], (PetscObject)dms[i]));
1343: PetscCall(PetscViewerGLVisGetFields_Internal(ctx->subv[i], &nf, NULL, NULL, NULL, NULL, NULL));
1344: tnf += nf;
1345: }
1346: PetscCall(PetscFree(dms));
1347: PetscCall(PetscMalloc3(tnf, &fecs, tnf, &sdim, tnf, &Ufds));
1348: for (i = 0, tnf = 0; i < n; i++) {
1349: PetscInt *sd, nf, f;
1350: const char **fec;
1351: Vec *Uf;
1353: PetscCall(PetscViewerGLVisGetFields_Internal(ctx->subv[i], &nf, &fec, &sd, NULL, (PetscObject **)&Uf, NULL));
1354: for (f = 0; f < nf; f++) {
1355: PetscCall(PetscStrallocpy(fec[f], &fecs[tnf + f]));
1356: Ufds[tnf + f] = Uf[f];
1357: sdim[tnf + f] = sd[f];
1358: }
1359: tnf += nf;
1360: }
1361: PetscCall(PetscViewerGLVisSetFields(viewer, tnf, (const char **)fecs, sdim, DMCompositeSampleGLVisFields_Private, (PetscObject *)Ufds, ctx, DestroyGLVisViewerCtx_Private));
1362: for (i = 0; i < tnf; i++) PetscCall(PetscFree(fecs[i]));
1363: PetscCall(PetscFree3(fecs, sdim, Ufds));
1364: PetscFunctionReturn(PETSC_SUCCESS);
1365: }
1367: static PetscErrorCode DMRefine_Composite(DM dmi, MPI_Comm comm, DM *fine)
1368: {
1369: struct DMCompositeLink *next;
1370: DM_Composite *com = (DM_Composite *)dmi->data;
1371: DM dm;
1373: PetscFunctionBegin;
1375: if (comm == MPI_COMM_NULL) PetscCall(PetscObjectGetComm((PetscObject)dmi, &comm));
1376: PetscCall(DMSetUp(dmi));
1377: next = com->next;
1378: PetscCall(DMCompositeCreate(comm, fine));
1380: /* loop over packed objects, handling one at a time */
1381: while (next) {
1382: PetscCall(DMRefine(next->dm, comm, &dm));
1383: PetscCall(DMCompositeAddDM(*fine, dm));
1384: PetscCall(PetscObjectDereference((PetscObject)dm));
1385: next = next->next;
1386: }
1387: PetscFunctionReturn(PETSC_SUCCESS);
1388: }
1390: static PetscErrorCode DMCoarsen_Composite(DM dmi, MPI_Comm comm, DM *fine)
1391: {
1392: struct DMCompositeLink *next;
1393: DM_Composite *com = (DM_Composite *)dmi->data;
1394: DM dm;
1396: PetscFunctionBegin;
1398: PetscCall(DMSetUp(dmi));
1399: if (comm == MPI_COMM_NULL) PetscCall(PetscObjectGetComm((PetscObject)dmi, &comm));
1400: next = com->next;
1401: PetscCall(DMCompositeCreate(comm, fine));
1403: /* loop over packed objects, handling one at a time */
1404: while (next) {
1405: PetscCall(DMCoarsen(next->dm, comm, &dm));
1406: PetscCall(DMCompositeAddDM(*fine, dm));
1407: PetscCall(PetscObjectDereference((PetscObject)dm));
1408: next = next->next;
1409: }
1410: PetscFunctionReturn(PETSC_SUCCESS);
1411: }
1413: static PetscErrorCode DMCreateInterpolation_Composite(DM coarse, DM fine, Mat *A, Vec *v)
1414: {
1415: PetscInt m, n, M, N, nDM, i;
1416: struct DMCompositeLink *nextc;
1417: struct DMCompositeLink *nextf;
1418: Vec gcoarse, gfine, *vecs;
1419: DM_Composite *comcoarse = (DM_Composite *)coarse->data;
1420: DM_Composite *comfine = (DM_Composite *)fine->data;
1421: Mat *mats;
1423: PetscFunctionBegin;
1426: PetscCall(DMSetUp(coarse));
1427: PetscCall(DMSetUp(fine));
1428: /* use global vectors only for determining matrix layout */
1429: PetscCall(DMGetGlobalVector(coarse, &gcoarse));
1430: PetscCall(DMGetGlobalVector(fine, &gfine));
1431: PetscCall(VecGetLocalSize(gcoarse, &n));
1432: PetscCall(VecGetLocalSize(gfine, &m));
1433: PetscCall(VecGetSize(gcoarse, &N));
1434: PetscCall(VecGetSize(gfine, &M));
1435: PetscCall(DMRestoreGlobalVector(coarse, &gcoarse));
1436: PetscCall(DMRestoreGlobalVector(fine, &gfine));
1438: nDM = comfine->nDM;
1439: 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);
1440: PetscCall(PetscCalloc1(nDM * nDM, &mats));
1441: if (v) PetscCall(PetscCalloc1(nDM, &vecs));
1443: /* loop over packed objects, handling one at a time */
1444: for (nextc = comcoarse->next, nextf = comfine->next, i = 0; nextc; nextc = nextc->next, nextf = nextf->next, i++) {
1445: if (!v) PetscCall(DMCreateInterpolation(nextc->dm, nextf->dm, &mats[i * nDM + i], NULL));
1446: else PetscCall(DMCreateInterpolation(nextc->dm, nextf->dm, &mats[i * nDM + i], &vecs[i]));
1447: }
1448: PetscCall(MatCreateNest(PetscObjectComm((PetscObject)fine), nDM, NULL, nDM, NULL, mats, A));
1449: if (v) PetscCall(VecCreateNest(PetscObjectComm((PetscObject)fine), nDM, NULL, vecs, v));
1450: for (i = 0; i < nDM * nDM; i++) PetscCall(MatDestroy(&mats[i]));
1451: PetscCall(PetscFree(mats));
1452: if (v) {
1453: for (i = 0; i < nDM; i++) PetscCall(VecDestroy(&vecs[i]));
1454: PetscCall(PetscFree(vecs));
1455: }
1456: PetscFunctionReturn(PETSC_SUCCESS);
1457: }
1459: static PetscErrorCode DMGetLocalToGlobalMapping_Composite(DM dm)
1460: {
1461: DM_Composite *com = (DM_Composite *)dm->data;
1462: ISLocalToGlobalMapping *ltogs;
1463: PetscInt i;
1465: PetscFunctionBegin;
1466: /* Set the ISLocalToGlobalMapping on the new matrix */
1467: PetscCall(DMCompositeGetISLocalToGlobalMappings(dm, <ogs));
1468: PetscCall(ISLocalToGlobalMappingConcatenate(PetscObjectComm((PetscObject)dm), com->nDM, ltogs, &dm->ltogmap));
1469: for (i = 0; i < com->nDM; i++) PetscCall(ISLocalToGlobalMappingDestroy(<ogs[i]));
1470: PetscCall(PetscFree(ltogs));
1471: PetscFunctionReturn(PETSC_SUCCESS);
1472: }
1474: static PetscErrorCode DMCreateColoring_Composite(DM dm, ISColoringType ctype, ISColoring *coloring)
1475: {
1476: PetscInt n, i, cnt;
1477: ISColoringValue *colors;
1478: PetscBool dense = PETSC_FALSE;
1479: ISColoringValue maxcol = 0;
1480: DM_Composite *com = (DM_Composite *)dm->data;
1482: PetscFunctionBegin;
1484: PetscCheck(ctype != IS_COLORING_LOCAL, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Only global coloring supported");
1485: if (ctype == IS_COLORING_GLOBAL) {
1486: n = com->n;
1487: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Unknown ISColoringType");
1488: PetscCall(PetscMalloc1(n, &colors)); /* freed in ISColoringDestroy() */
1490: PetscCall(PetscOptionsGetBool(((PetscObject)dm)->options, ((PetscObject)dm)->prefix, "-dmcomposite_dense_jacobian", &dense, NULL));
1491: if (dense) {
1492: for (i = 0; i < n; i++) colors[i] = (ISColoringValue)(com->rstart + i);
1493: maxcol = com->N;
1494: } else {
1495: struct DMCompositeLink *next = com->next;
1496: PetscMPIInt rank;
1498: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
1499: cnt = 0;
1500: while (next) {
1501: ISColoring lcoloring;
1503: PetscCall(DMCreateColoring(next->dm, IS_COLORING_GLOBAL, &lcoloring));
1504: for (i = 0; i < lcoloring->N; i++) colors[cnt++] = maxcol + lcoloring->colors[i];
1505: maxcol += lcoloring->n;
1506: PetscCall(ISColoringDestroy(&lcoloring));
1507: next = next->next;
1508: }
1509: }
1510: PetscCall(ISColoringCreate(PetscObjectComm((PetscObject)dm), maxcol, n, colors, PETSC_OWN_POINTER, coloring));
1511: PetscFunctionReturn(PETSC_SUCCESS);
1512: }
1514: static PetscErrorCode DMGlobalToLocalBegin_Composite(DM dm, Vec gvec, InsertMode mode, Vec lvec)
1515: {
1516: struct DMCompositeLink *next;
1517: PetscScalar *garray, *larray;
1518: DM_Composite *com = (DM_Composite *)dm->data;
1520: PetscFunctionBegin;
1524: if (!com->setup) PetscCall(DMSetUp(dm));
1526: PetscCall(VecGetArray(gvec, &garray));
1527: PetscCall(VecGetArray(lvec, &larray));
1529: /* loop over packed objects, handling one at a time */
1530: next = com->next;
1531: while (next) {
1532: Vec local, global;
1533: PetscInt N;
1535: PetscCall(DMGetGlobalVector(next->dm, &global));
1536: PetscCall(VecGetLocalSize(global, &N));
1537: PetscCall(VecPlaceArray(global, garray));
1538: PetscCall(DMGetLocalVector(next->dm, &local));
1539: PetscCall(VecPlaceArray(local, larray));
1540: PetscCall(DMGlobalToLocalBegin(next->dm, global, mode, local));
1541: PetscCall(DMGlobalToLocalEnd(next->dm, global, mode, local));
1542: PetscCall(VecResetArray(global));
1543: PetscCall(VecResetArray(local));
1544: PetscCall(DMRestoreGlobalVector(next->dm, &global));
1545: PetscCall(DMRestoreLocalVector(next->dm, &local));
1547: larray += next->nlocal;
1548: garray += next->n;
1549: next = next->next;
1550: }
1552: PetscCall(VecRestoreArray(gvec, NULL));
1553: PetscCall(VecRestoreArray(lvec, NULL));
1554: PetscFunctionReturn(PETSC_SUCCESS);
1555: }
1557: static PetscErrorCode DMGlobalToLocalEnd_Composite(DM dm, Vec gvec, InsertMode mode, Vec lvec)
1558: {
1559: PetscFunctionBegin;
1563: PetscFunctionReturn(PETSC_SUCCESS);
1564: }
1566: static PetscErrorCode DMLocalToGlobalBegin_Composite(DM dm, Vec lvec, InsertMode mode, Vec gvec)
1567: {
1568: struct DMCompositeLink *next;
1569: PetscScalar *larray, *garray;
1570: DM_Composite *com = (DM_Composite *)dm->data;
1572: PetscFunctionBegin;
1577: if (!com->setup) PetscCall(DMSetUp(dm));
1579: PetscCall(VecGetArray(lvec, &larray));
1580: PetscCall(VecGetArray(gvec, &garray));
1582: /* loop over packed objects, handling one at a time */
1583: next = com->next;
1584: while (next) {
1585: Vec global, local;
1587: PetscCall(DMGetLocalVector(next->dm, &local));
1588: PetscCall(VecPlaceArray(local, larray));
1589: PetscCall(DMGetGlobalVector(next->dm, &global));
1590: PetscCall(VecPlaceArray(global, garray));
1591: PetscCall(DMLocalToGlobalBegin(next->dm, local, mode, global));
1592: PetscCall(DMLocalToGlobalEnd(next->dm, local, mode, global));
1593: PetscCall(VecResetArray(local));
1594: PetscCall(VecResetArray(global));
1595: PetscCall(DMRestoreGlobalVector(next->dm, &global));
1596: PetscCall(DMRestoreLocalVector(next->dm, &local));
1598: garray += next->n;
1599: larray += next->nlocal;
1600: next = next->next;
1601: }
1603: PetscCall(VecRestoreArray(gvec, NULL));
1604: PetscCall(VecRestoreArray(lvec, NULL));
1605: PetscFunctionReturn(PETSC_SUCCESS);
1606: }
1608: static PetscErrorCode DMLocalToGlobalEnd_Composite(DM dm, Vec lvec, InsertMode mode, Vec gvec)
1609: {
1610: PetscFunctionBegin;
1614: PetscFunctionReturn(PETSC_SUCCESS);
1615: }
1617: static PetscErrorCode DMLocalToLocalBegin_Composite(DM dm, Vec vec1, InsertMode mode, Vec vec2)
1618: {
1619: struct DMCompositeLink *next;
1620: PetscScalar *array1, *array2;
1621: DM_Composite *com = (DM_Composite *)dm->data;
1623: PetscFunctionBegin;
1628: if (!com->setup) PetscCall(DMSetUp(dm));
1630: PetscCall(VecGetArray(vec1, &array1));
1631: PetscCall(VecGetArray(vec2, &array2));
1633: /* loop over packed objects, handling one at a time */
1634: next = com->next;
1635: while (next) {
1636: Vec local1, local2;
1638: PetscCall(DMGetLocalVector(next->dm, &local1));
1639: PetscCall(VecPlaceArray(local1, array1));
1640: PetscCall(DMGetLocalVector(next->dm, &local2));
1641: PetscCall(VecPlaceArray(local2, array2));
1642: PetscCall(DMLocalToLocalBegin(next->dm, local1, mode, local2));
1643: PetscCall(DMLocalToLocalEnd(next->dm, local1, mode, local2));
1644: PetscCall(VecResetArray(local2));
1645: PetscCall(DMRestoreLocalVector(next->dm, &local2));
1646: PetscCall(VecResetArray(local1));
1647: PetscCall(DMRestoreLocalVector(next->dm, &local1));
1649: array1 += next->nlocal;
1650: array2 += next->nlocal;
1651: next = next->next;
1652: }
1654: PetscCall(VecRestoreArray(vec1, NULL));
1655: PetscCall(VecRestoreArray(vec2, NULL));
1656: PetscFunctionReturn(PETSC_SUCCESS);
1657: }
1659: static PetscErrorCode DMLocalToLocalEnd_Composite(DM dm, Vec lvec, InsertMode mode, Vec gvec)
1660: {
1661: PetscFunctionBegin;
1665: PetscFunctionReturn(PETSC_SUCCESS);
1666: }
1668: /*MC
1669: DMCOMPOSITE = "composite" - A `DM` object that is used to manage data for a collection of `DM`
1671: Level: intermediate
1673: .seealso: `DMType`, `DM`, `DMDACreate()`, `DMCreate()`, `DMSetType()`, `DMCompositeCreate()`
1674: M*/
1676: PETSC_EXTERN PetscErrorCode DMCreate_Composite(DM p)
1677: {
1678: DM_Composite *com;
1680: PetscFunctionBegin;
1681: PetscCall(PetscNew(&com));
1682: p->data = com;
1683: com->n = 0;
1684: com->nghost = 0;
1685: com->next = NULL;
1686: com->nDM = 0;
1688: p->ops->createglobalvector = DMCreateGlobalVector_Composite;
1689: p->ops->createlocalvector = DMCreateLocalVector_Composite;
1690: p->ops->getlocaltoglobalmapping = DMGetLocalToGlobalMapping_Composite;
1691: p->ops->createfieldis = DMCreateFieldIS_Composite;
1692: p->ops->createfielddecomposition = DMCreateFieldDecomposition_Composite;
1693: p->ops->refine = DMRefine_Composite;
1694: p->ops->coarsen = DMCoarsen_Composite;
1695: p->ops->createinterpolation = DMCreateInterpolation_Composite;
1696: p->ops->creatematrix = DMCreateMatrix_Composite;
1697: p->ops->getcoloring = DMCreateColoring_Composite;
1698: p->ops->globaltolocalbegin = DMGlobalToLocalBegin_Composite;
1699: p->ops->globaltolocalend = DMGlobalToLocalEnd_Composite;
1700: p->ops->localtoglobalbegin = DMLocalToGlobalBegin_Composite;
1701: p->ops->localtoglobalend = DMLocalToGlobalEnd_Composite;
1702: p->ops->localtolocalbegin = DMLocalToLocalBegin_Composite;
1703: p->ops->localtolocalend = DMLocalToLocalEnd_Composite;
1704: p->ops->destroy = DMDestroy_Composite;
1705: p->ops->view = DMView_Composite;
1706: p->ops->setup = DMSetUp_Composite;
1708: PetscCall(PetscObjectComposeFunction((PetscObject)p, "DMSetUpGLVisViewer_C", DMSetUpGLVisViewer_Composite));
1709: PetscFunctionReturn(PETSC_SUCCESS);
1710: }
1712: /*@
1713: DMCompositeCreate - Creates a `DMCOMPOSITE`, used to generate "composite"
1714: vectors made up of several subvectors.
1716: Collective
1718: Input Parameter:
1719: . comm - the processors that will share the global vector
1721: Output Parameter:
1722: . packer - the `DMCOMPOSITE` object
1724: Level: advanced
1726: .seealso: `DMCOMPOSITE`, `DM`, `DMDestroy()`, `DMCompositeAddDM()`, `DMCompositeScatter()`, `DMCreate()`
1727: `DMCompositeGather()`, `DMCreateGlobalVector()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeGetAccess()`
1728: `DMCompositeGetLocalVectors()`, `DMCompositeRestoreLocalVectors()`, `DMCompositeGetEntries()`
1729: @*/
1730: PetscErrorCode DMCompositeCreate(MPI_Comm comm, DM *packer)
1731: {
1732: PetscFunctionBegin;
1733: PetscAssertPointer(packer, 2);
1734: PetscCall(DMCreate(comm, packer));
1735: PetscCall(DMSetType(*packer, DMCOMPOSITE));
1736: PetscFunctionReturn(PETSC_SUCCESS);
1737: }