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: Notes:
19: See `DMSetApplicationContext()` and `DMGetApplicationContext()` for how to get user information into
20: this routine
22: The provided function should have a signature matching
23: .vb
24: PetscErrorCode your_form_couple_method(DM dm, Mat J, PetscInt *dnz, PetscInt *onz, PetscInt rstart, PetscInt nrows, PetscInt start, PetscInt end);
25: .ve
26: where
27: + dm - the composite object
28: . J - the constructed matrix, or `NULL`. If provided, the function should fill the existing nonzero pattern with zeros (only `dm` and `rstart` are valid in this case).
29: . dnz - array counting the number of on-diagonal non-zero entries per row, where on-diagonal means that this process owns both the row and column
30: . onz - array counting the number of off-diagonal non-zero entries per row, where off-diagonal means that this process owns the row
31: . rstart - offset into `*nz` arrays, for local row index `r`, update `onz[r - rstart]` or `dnz[r - rstart]`
32: . nrows - number of owned global rows
33: . start - the first owned global index
34: - end - the last owned global index + 1
36: If `J` is not `NULL`, then the only other valid parameter is `rstart`
38: The user coupling function has a weird and poorly documented interface and is not tested, it should be removed
40: .seealso: `DMCOMPOSITE`, `DM`
41: @*/
42: PetscErrorCode DMCompositeSetCoupling(DM dm, PetscErrorCode (*FormCoupleLocations)(DM, Mat, PetscInt *, PetscInt *, PetscInt, PetscInt, PetscInt, PetscInt))
43: {
44: DM_Composite *com = (DM_Composite *)dm->data;
45: PetscBool flg;
47: PetscFunctionBegin;
48: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
49: PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
50: com->FormCoupleLocations = FormCoupleLocations;
51: PetscFunctionReturn(PETSC_SUCCESS);
52: }
54: static PetscErrorCode DMDestroy_Composite(DM dm)
55: {
56: struct DMCompositeLink *next, *prev;
57: DM_Composite *com = (DM_Composite *)dm->data;
59: PetscFunctionBegin;
60: next = com->next;
61: while (next) {
62: prev = next;
63: next = next->next;
64: PetscCall(DMDestroy(&prev->dm));
65: PetscCall(PetscFree(prev->grstarts));
66: PetscCall(PetscFree(prev));
67: }
68: PetscCall(PetscObjectComposeFunction((PetscObject)dm, "DMSetUpGLVisViewer_C", NULL));
69: /* This was originally freed in DMDestroy(), but that prevents reference counting of backend objects */
70: PetscCall(PetscFree(com));
71: PetscFunctionReturn(PETSC_SUCCESS);
72: }
74: static PetscErrorCode DMView_Composite(DM dm, PetscViewer v)
75: {
76: PetscBool isascii;
77: DM_Composite *com = (DM_Composite *)dm->data;
79: PetscFunctionBegin;
80: PetscCall(PetscObjectTypeCompare((PetscObject)v, PETSCVIEWERASCII, &isascii));
81: if (isascii) {
82: struct DMCompositeLink *lnk = com->next;
83: PetscInt i;
85: PetscCall(PetscViewerASCIIPrintf(v, "DM (%s)\n", ((PetscObject)dm)->prefix ? ((PetscObject)dm)->prefix : "no prefix"));
86: PetscCall(PetscViewerASCIIPrintf(v, " contains %" PetscInt_FMT " DMs\n", com->nDM));
87: PetscCall(PetscViewerASCIIPushTab(v));
88: for (i = 0; lnk; lnk = lnk->next, i++) {
89: PetscCall(PetscViewerASCIIPrintf(v, "Link %" PetscInt_FMT ": DM of type %s\n", i, ((PetscObject)lnk->dm)->type_name));
90: PetscCall(PetscViewerASCIIPushTab(v));
91: PetscCall(DMView(lnk->dm, v));
92: PetscCall(PetscViewerASCIIPopTab(v));
93: }
94: PetscCall(PetscViewerASCIIPopTab(v));
95: }
96: PetscFunctionReturn(PETSC_SUCCESS);
97: }
99: static PetscErrorCode DMSetUp_Composite(DM dm)
100: {
101: PetscInt nprev = 0;
102: PetscMPIInt rank, size;
103: DM_Composite *com = (DM_Composite *)dm->data;
104: struct DMCompositeLink *next = com->next;
105: PetscLayout map;
107: PetscFunctionBegin;
108: PetscCheck(!com->setup, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "Packer has already been setup");
109: PetscCall(PetscLayoutCreate(PetscObjectComm((PetscObject)dm), &map));
110: PetscCall(PetscLayoutSetLocalSize(map, com->n));
111: PetscCall(PetscLayoutSetSize(map, PETSC_DETERMINE));
112: PetscCall(PetscLayoutSetBlockSize(map, 1));
113: PetscCall(PetscLayoutSetUp(map));
114: PetscCall(PetscLayoutGetSize(map, &com->N));
115: PetscCall(PetscLayoutGetRange(map, &com->rstart, NULL));
116: PetscCall(PetscLayoutDestroy(&map));
118: /* now set the rstart for each linked vector */
119: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
120: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size));
121: while (next) {
122: next->rstart = nprev;
123: nprev += next->n;
124: next->grstart = com->rstart + next->rstart;
125: PetscCall(PetscMalloc1(size, &next->grstarts));
126: PetscCallMPI(MPI_Allgather(&next->grstart, 1, MPIU_INT, next->grstarts, 1, MPIU_INT, PetscObjectComm((PetscObject)dm)));
127: next = next->next;
128: }
129: com->setup = PETSC_TRUE;
130: PetscFunctionReturn(PETSC_SUCCESS);
131: }
133: /*@
134: DMCompositeGetNumberDM - Gets the number of `DM` objects in the `DMCOMPOSITE`
135: representation.
137: Not Collective
139: Input Parameter:
140: . dm - the `DMCOMPOSITE` object
142: Output Parameter:
143: . nDM - the number of `DM`
145: Level: beginner
147: .seealso: `DMCOMPOSITE`, `DM`
148: @*/
149: PetscErrorCode DMCompositeGetNumberDM(DM dm, PetscInt *nDM)
150: {
151: DM_Composite *com = (DM_Composite *)dm->data;
152: PetscBool flg;
154: PetscFunctionBegin;
156: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
157: PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
158: *nDM = com->nDM;
159: PetscFunctionReturn(PETSC_SUCCESS);
160: }
162: /*@C
163: DMCompositeGetAccess - Allows one to access the individual packed vectors in their global
164: representation.
166: Collective
168: Input Parameters:
169: + dm - the `DMCOMPOSITE` object
170: - gvec - the global vector
172: Output Parameter:
173: . ... - the packed parallel vectors, `NULL` for those that are not needed
175: Level: advanced
177: Note:
178: Use `DMCompositeRestoreAccess()` to return the vectors when you no longer need them
180: .seealso: `DMCOMPOSITE`, `DM`, `DMCompositeGetEntries()`, `DMCompositeScatter()`
181: @*/
182: PetscErrorCode DMCompositeGetAccess(DM dm, Vec gvec, ...)
183: {
184: va_list Argp;
185: struct DMCompositeLink *next;
186: DM_Composite *com = (DM_Composite *)dm->data;
187: PetscInt readonly;
188: PetscBool flg;
190: PetscFunctionBegin;
193: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
194: PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
195: next = com->next;
196: if (!com->setup) PetscCall(DMSetUp(dm));
198: PetscCall(VecLockGet(gvec, &readonly));
199: /* loop over packed objects, handling one at a time */
200: va_start(Argp, gvec);
201: while (next) {
202: Vec *vec;
203: vec = va_arg(Argp, Vec *);
204: if (vec) {
205: PetscCall(DMGetGlobalVector(next->dm, vec));
206: if (readonly) {
207: const PetscScalar *array;
208: PetscCall(VecGetArrayRead(gvec, &array));
209: PetscCall(VecPlaceArray(*vec, array + next->rstart));
210: PetscCall(VecLockReadPush(*vec));
211: PetscCall(VecRestoreArrayRead(gvec, &array));
212: } else {
213: PetscScalar *array;
214: PetscCall(VecGetArray(gvec, &array));
215: PetscCall(VecPlaceArray(*vec, array + next->rstart));
216: PetscCall(VecRestoreArray(gvec, &array));
217: }
218: }
219: next = next->next;
220: }
221: va_end(Argp);
222: PetscFunctionReturn(PETSC_SUCCESS);
223: }
225: /*@
226: DMCompositeGetAccessArray - Allows one to access the individual packed vectors in their global
227: representation.
229: Collective
231: Input Parameters:
232: + dm - the `DMCOMPOSITE`
233: . pvec - packed vector
234: . nwanted - number of vectors wanted
235: - wanted - sorted array of integers indicating thde vectors wanted, or `NULL` to get all vectors, length `nwanted`
237: Output Parameter:
238: . vecs - array of requested global vectors (must be previously allocated and of length `nwanted`)
240: Level: advanced
242: Note:
243: Use `DMCompositeRestoreAccessArray()` to return the vectors when you no longer need them
245: .seealso: `DMCOMPOSITE`, `DM`, `DMCompositeGetAccess()`, `DMCompositeGetEntries()`, `DMCompositeScatter()`, `DMCompositeGather()`
246: @*/
247: PetscErrorCode DMCompositeGetAccessArray(DM dm, Vec pvec, PetscInt nwanted, const PetscInt wanted[], Vec vecs[])
248: {
249: struct DMCompositeLink *link;
250: PetscInt i, wnum;
251: DM_Composite *com = (DM_Composite *)dm->data;
252: PetscInt readonly;
253: PetscBool flg;
255: PetscFunctionBegin;
258: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
259: PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
260: if (!com->setup) PetscCall(DMSetUp(dm));
262: PetscCall(VecLockGet(pvec, &readonly));
263: for (i = 0, wnum = 0, link = com->next; link && wnum < nwanted; i++, link = link->next) {
264: if (!wanted || i == wanted[wnum]) {
265: Vec v;
266: PetscCall(DMGetGlobalVector(link->dm, &v));
267: if (readonly) {
268: const PetscScalar *array;
269: PetscCall(VecGetArrayRead(pvec, &array));
270: PetscCall(VecPlaceArray(v, array + link->rstart));
271: PetscCall(VecLockReadPush(v));
272: PetscCall(VecRestoreArrayRead(pvec, &array));
273: } else {
274: PetscScalar *array;
275: PetscCall(VecGetArray(pvec, &array));
276: PetscCall(VecPlaceArray(v, array + link->rstart));
277: PetscCall(VecRestoreArray(pvec, &array));
278: }
279: vecs[wnum++] = v;
280: }
281: }
282: PetscFunctionReturn(PETSC_SUCCESS);
283: }
285: /*@
286: DMCompositeGetLocalAccessArray - Allows one to access the individual
287: packed vectors in their local representation.
289: Collective
291: Input Parameters:
292: + dm - the `DMCOMPOSITE`
293: . pvec - packed vector
294: . nwanted - number of vectors wanted
295: - wanted - sorted array of vectors wanted, or `NULL` to get all vectors, length `nwanted`
297: Output Parameter:
298: . vecs - array of requested local vectors (must be allocated and of length `nwanted`)
300: Level: advanced
302: Note:
303: Use `DMCompositeRestoreLocalAccessArray()` to return the vectors
304: when you no longer need them.
306: .seealso: `DMCOMPOSITE`, `DM`, `DMCompositeRestoreLocalAccessArray()`, `DMCompositeGetAccess()`,
307: `DMCompositeGetEntries()`, `DMCompositeScatter()`, `DMCompositeGather()`
308: @*/
309: PetscErrorCode DMCompositeGetLocalAccessArray(DM dm, Vec pvec, PetscInt nwanted, const PetscInt wanted[], Vec vecs[])
310: {
311: struct DMCompositeLink *link;
312: PetscInt i, wnum;
313: DM_Composite *com = (DM_Composite *)dm->data;
314: PetscInt readonly;
315: PetscInt nlocal = 0;
316: PetscBool flg;
318: PetscFunctionBegin;
321: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
322: PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
323: if (!com->setup) PetscCall(DMSetUp(dm));
325: PetscCall(VecLockGet(pvec, &readonly));
326: for (i = 0, wnum = 0, link = com->next; link && wnum < nwanted; i++, link = link->next) {
327: if (!wanted || i == wanted[wnum]) {
328: Vec v;
329: PetscCall(DMGetLocalVector(link->dm, &v));
330: if (readonly) {
331: const PetscScalar *array;
332: PetscCall(VecGetArrayRead(pvec, &array));
333: PetscCall(VecPlaceArray(v, array + nlocal));
334: // 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
335: PetscCall(VecLockReadPush(v));
336: PetscCall(VecRestoreArrayRead(pvec, &array));
337: } else {
338: PetscScalar *array;
339: PetscCall(VecGetArray(pvec, &array));
340: PetscCall(VecPlaceArray(v, array + nlocal));
341: PetscCall(VecRestoreArray(pvec, &array));
342: }
343: vecs[wnum++] = v;
344: }
346: nlocal += link->nlocal;
347: }
348: PetscFunctionReturn(PETSC_SUCCESS);
349: }
351: /*@C
352: DMCompositeRestoreAccess - Returns the vectors obtained with `DMCompositeGetAccess()`
353: representation.
355: Collective
357: Input Parameters:
358: + dm - the `DMCOMPOSITE` object
359: . gvec - the global vector
360: - ... - the individual parallel vectors, `NULL` for those that are not needed
362: Level: advanced
364: .seealso: `DMCOMPOSITE`, `DM`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`,
365: `DMCompositeGather()`, `DMCompositeCreate()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeScatter()`,
366: `DMCompositeGetAccess()`
367: @*/
368: PetscErrorCode DMCompositeRestoreAccess(DM dm, Vec gvec, ...)
369: {
370: va_list Argp;
371: struct DMCompositeLink *next;
372: DM_Composite *com = (DM_Composite *)dm->data;
373: PetscInt readonly;
374: PetscBool flg;
376: PetscFunctionBegin;
379: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
380: PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
381: next = com->next;
382: if (!com->setup) PetscCall(DMSetUp(dm));
384: PetscCall(VecLockGet(gvec, &readonly));
385: /* loop over packed objects, handling one at a time */
386: va_start(Argp, gvec);
387: while (next) {
388: Vec *vec;
389: vec = va_arg(Argp, Vec *);
390: if (vec) {
391: PetscCall(VecResetArray(*vec));
392: if (readonly) PetscCall(VecLockReadPop(*vec));
393: PetscCall(DMRestoreGlobalVector(next->dm, vec));
394: }
395: next = next->next;
396: }
397: va_end(Argp);
398: PetscFunctionReturn(PETSC_SUCCESS);
399: }
401: /*@
402: DMCompositeRestoreAccessArray - Returns the vectors obtained with `DMCompositeGetAccessArray()`
404: Collective
406: Input Parameters:
407: + dm - the `DMCOMPOSITE` object
408: . pvec - packed vector
409: . nwanted - number of vectors wanted
410: . wanted - sorted array of vectors wanted, or `NULL` to restore all vectors
411: - vecs - array of global vectors
413: Level: advanced
415: .seealso: `DMCOMPOSITE`, `DM`, `DMCompositeRestoreAccess()`, `DMCompositeRestoreEntries()`, `DMCompositeScatter()`, `DMCompositeGather()`
416: @*/
417: PetscErrorCode DMCompositeRestoreAccessArray(DM dm, Vec pvec, PetscInt nwanted, const PetscInt wanted[], Vec vecs[])
418: {
419: struct DMCompositeLink *link;
420: PetscInt i, wnum;
421: DM_Composite *com = (DM_Composite *)dm->data;
422: PetscInt readonly;
423: PetscBool flg;
425: PetscFunctionBegin;
428: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
429: PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
430: if (!com->setup) PetscCall(DMSetUp(dm));
432: PetscCall(VecLockGet(pvec, &readonly));
433: for (i = 0, wnum = 0, link = com->next; link && wnum < nwanted; i++, link = link->next) {
434: if (!wanted || i == wanted[wnum]) {
435: PetscCall(VecResetArray(vecs[wnum]));
436: if (readonly) PetscCall(VecLockReadPop(vecs[wnum]));
437: PetscCall(DMRestoreGlobalVector(link->dm, &vecs[wnum]));
438: wnum++;
439: }
440: }
441: PetscFunctionReturn(PETSC_SUCCESS);
442: }
444: /*@
445: DMCompositeRestoreLocalAccessArray - Returns the vectors obtained with `DMCompositeGetLocalAccessArray()`.
447: Collective
449: Input Parameters:
450: + dm - the `DMCOMPOSITE` object
451: . pvec - packed vector
452: . nwanted - number of vectors wanted
453: . wanted - sorted array of vectors wanted, or `NULL` to restore all vectors
454: - vecs - array of local vectors
456: Level: advanced
458: Note:
459: `nwanted` and `wanted` must match the values given to `DMCompositeGetLocalAccessArray()`
460: otherwise the call will fail.
462: .seealso: `DMCOMPOSITE`, `DM`, `DMCompositeGetLocalAccessArray()`, `DMCompositeRestoreAccessArray()`,
463: `DMCompositeRestoreAccess()`, `DMCompositeRestoreEntries()`,
464: `DMCompositeScatter()`, `DMCompositeGather()`
465: @*/
466: PetscErrorCode DMCompositeRestoreLocalAccessArray(DM dm, Vec pvec, PetscInt nwanted, const PetscInt wanted[], Vec *vecs)
467: {
468: struct DMCompositeLink *link;
469: PetscInt i, wnum;
470: DM_Composite *com = (DM_Composite *)dm->data;
471: PetscInt readonly;
472: PetscBool flg;
474: PetscFunctionBegin;
477: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
478: PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
479: if (!com->setup) PetscCall(DMSetUp(dm));
481: PetscCall(VecLockGet(pvec, &readonly));
482: for (i = 0, wnum = 0, link = com->next; link && wnum < nwanted; i++, link = link->next) {
483: if (!wanted || i == wanted[wnum]) {
484: PetscCall(VecResetArray(vecs[wnum]));
485: if (readonly) PetscCall(VecLockReadPop(vecs[wnum]));
486: PetscCall(DMRestoreLocalVector(link->dm, &vecs[wnum]));
487: wnum++;
488: }
489: }
490: PetscFunctionReturn(PETSC_SUCCESS);
491: }
493: /*@C
494: DMCompositeScatter - Scatters from a global packed vector into its individual local vectors
496: Collective
498: Input Parameters:
499: + dm - the `DMCOMPOSITE` object
500: . gvec - the global vector
501: - ... - the individual sequential vectors, `NULL` for those that are not needed
503: Level: advanced
505: Note:
506: `DMCompositeScatterArray()` is a non-variadic alternative that is often more convenient for library callers and is
507: accessible from Fortran.
509: .seealso: `DMCOMPOSITE`, `DM`, `DMDestroy()`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`,
510: `DMCompositeGather()`, `DMCompositeCreate()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeGetAccess()`,
511: `DMCompositeGetLocalVectors()`, `DMCompositeRestoreLocalVectors()`, `DMCompositeGetEntries()`
512: `DMCompositeScatterArray()`
513: @*/
514: PetscErrorCode DMCompositeScatter(DM dm, Vec gvec, ...)
515: {
516: va_list Argp;
517: struct DMCompositeLink *next;
518: PETSC_UNUSED PetscInt cnt;
519: DM_Composite *com = (DM_Composite *)dm->data;
520: PetscBool flg;
522: PetscFunctionBegin;
525: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
526: PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
527: if (!com->setup) PetscCall(DMSetUp(dm));
529: /* loop over packed objects, handling one at a time */
530: va_start(Argp, gvec);
531: for (cnt = 3, next = com->next; next; cnt++, next = next->next) {
532: Vec local;
533: local = va_arg(Argp, Vec);
534: if (local) {
535: Vec global;
536: const PetscScalar *array;
538: PetscCall(DMGetGlobalVector(next->dm, &global));
539: PetscCall(VecGetArrayRead(gvec, &array));
540: PetscCall(VecPlaceArray(global, array + next->rstart));
541: PetscCall(DMGlobalToLocalBegin(next->dm, global, INSERT_VALUES, local));
542: PetscCall(DMGlobalToLocalEnd(next->dm, global, INSERT_VALUES, local));
543: PetscCall(VecRestoreArrayRead(gvec, &array));
544: PetscCall(VecResetArray(global));
545: PetscCall(DMRestoreGlobalVector(next->dm, &global));
546: }
547: }
548: va_end(Argp);
549: PetscFunctionReturn(PETSC_SUCCESS);
550: }
552: /*@
553: DMCompositeScatterArray - Scatters from a global packed vector into its individual local vectors
555: Collective
557: Input Parameters:
558: + dm - the `DMCOMPOSITE` object
559: . gvec - the global vector
560: - lvecs - array of local vectors, NULL for any that are not needed
562: Level: advanced
564: Note:
565: This is a non-variadic alternative to `DMCompositeScatter()`
567: .seealso: `DMCOMPOSITE`, `DM`, `DMDestroy()`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`
568: `DMCompositeGather()`, `DMCompositeCreate()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeGetAccess()`,
569: `DMCompositeGetLocalVectors()`, `DMCompositeRestoreLocalVectors()`, `DMCompositeGetEntries()`
570: @*/
571: PetscErrorCode DMCompositeScatterArray(DM dm, Vec gvec, Vec *lvecs)
572: {
573: struct DMCompositeLink *next;
574: PetscInt i;
575: DM_Composite *com = (DM_Composite *)dm->data;
576: PetscBool flg;
578: PetscFunctionBegin;
581: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
582: PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
583: if (!com->setup) PetscCall(DMSetUp(dm));
585: /* loop over packed objects, handling one at a time */
586: for (i = 0, next = com->next; next; next = next->next, i++) {
587: if (lvecs[i]) {
588: Vec global;
589: const PetscScalar *array;
591: PetscCall(DMGetGlobalVector(next->dm, &global));
592: PetscCall(VecGetArrayRead(gvec, &array));
593: PetscCall(VecPlaceArray(global, (PetscScalar *)array + next->rstart));
594: PetscCall(DMGlobalToLocalBegin(next->dm, global, INSERT_VALUES, lvecs[i]));
595: PetscCall(DMGlobalToLocalEnd(next->dm, global, INSERT_VALUES, lvecs[i]));
596: PetscCall(VecRestoreArrayRead(gvec, &array));
597: PetscCall(VecResetArray(global));
598: PetscCall(DMRestoreGlobalVector(next->dm, &global));
599: }
600: }
601: PetscFunctionReturn(PETSC_SUCCESS);
602: }
604: /*@C
605: DMCompositeGather - Gathers into a global packed vector from its individual local vectors
607: Collective
609: Input Parameters:
610: + dm - the `DMCOMPOSITE` object
611: . imode - `INSERT_VALUES` or `ADD_VALUES`
612: . gvec - the global vector
613: - ... - the individual sequential vectors, `NULL` for any that are not needed
615: Level: advanced
617: Fortran Notes:
618: Fortran users should use `DMCompositeGatherArray()`
620: .seealso: `DMCOMPOSITE`, `DM`, `DMDestroy()`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`,
621: `DMCompositeScatter()`, `DMCompositeCreate()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeGetAccess()`,
622: `DMCompositeGetLocalVectors()`, `DMCompositeRestoreLocalVectors()`, `DMCompositeGetEntries()`
623: @*/
624: PetscErrorCode DMCompositeGather(DM dm, InsertMode imode, Vec gvec, ...)
625: {
626: va_list Argp;
627: struct DMCompositeLink *next;
628: DM_Composite *com = (DM_Composite *)dm->data;
629: PETSC_UNUSED PetscInt cnt;
630: PetscBool flg;
632: PetscFunctionBegin;
635: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
636: PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
637: if (!com->setup) PetscCall(DMSetUp(dm));
639: /* loop over packed objects, handling one at a time */
640: va_start(Argp, gvec);
641: for (cnt = 3, next = com->next; next; cnt++, next = next->next) {
642: Vec local;
643: local = va_arg(Argp, Vec);
644: if (local) {
645: PetscScalar *array;
646: Vec global;
648: PetscCall(DMGetGlobalVector(next->dm, &global));
649: PetscCall(VecGetArray(gvec, &array));
650: PetscCall(VecPlaceArray(global, array + next->rstart));
651: PetscCall(DMLocalToGlobalBegin(next->dm, local, imode, global));
652: PetscCall(DMLocalToGlobalEnd(next->dm, local, imode, global));
653: PetscCall(VecRestoreArray(gvec, &array));
654: PetscCall(VecResetArray(global));
655: PetscCall(DMRestoreGlobalVector(next->dm, &global));
656: }
657: }
658: va_end(Argp);
659: PetscFunctionReturn(PETSC_SUCCESS);
660: }
662: /*@
663: DMCompositeGatherArray - Gathers into a global packed vector from its individual local vectors
665: Collective
667: Input Parameters:
668: + dm - the `DMCOMPOSITE` object
669: . gvec - the global vector
670: . imode - `INSERT_VALUES` or `ADD_VALUES`
671: - lvecs - the individual sequential vectors, NULL for any that are not needed
673: Level: advanced
675: Note:
676: This is a non-variadic alternative to `DMCompositeGather()`.
678: .seealso: `DMCOMPOSITE`, `DM`, `DMDestroy()`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`,
679: `DMCompositeScatter()`, `DMCompositeCreate()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeGetAccess()`,
680: `DMCompositeGetLocalVectors()`, `DMCompositeRestoreLocalVectors()`, `DMCompositeGetEntries()`,
681: @*/
682: PetscErrorCode DMCompositeGatherArray(DM dm, InsertMode imode, Vec gvec, Vec *lvecs)
683: {
684: struct DMCompositeLink *next;
685: DM_Composite *com = (DM_Composite *)dm->data;
686: PetscInt i;
687: PetscBool flg;
689: PetscFunctionBegin;
692: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
693: PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
694: if (!com->setup) PetscCall(DMSetUp(dm));
696: /* loop over packed objects, handling one at a time */
697: for (next = com->next, i = 0; next; next = next->next, i++) {
698: if (lvecs[i]) {
699: PetscScalar *array;
700: Vec global;
702: PetscCall(DMGetGlobalVector(next->dm, &global));
703: PetscCall(VecGetArray(gvec, &array));
704: PetscCall(VecPlaceArray(global, array + next->rstart));
705: PetscCall(DMLocalToGlobalBegin(next->dm, lvecs[i], imode, global));
706: PetscCall(DMLocalToGlobalEnd(next->dm, lvecs[i], imode, global));
707: PetscCall(VecRestoreArray(gvec, &array));
708: PetscCall(VecResetArray(global));
709: PetscCall(DMRestoreGlobalVector(next->dm, &global));
710: }
711: }
712: PetscFunctionReturn(PETSC_SUCCESS);
713: }
715: /*@
716: DMCompositeAddDM - adds a `DM` vector to a `DMCOMPOSITE`
718: Collective
720: Input Parameters:
721: + dmc - the `DMCOMPOSITE` object
722: - dm - the `DM` object
724: Level: advanced
726: .seealso: `DMCOMPOSITE`, `DM`, `DMDestroy()`, `DMCompositeGather()`, `DMCreateGlobalVector()`,
727: `DMCompositeScatter()`, `DMCompositeCreate()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeGetAccess()`,
728: `DMCompositeGetLocalVectors()`, `DMCompositeRestoreLocalVectors()`, `DMCompositeGetEntries()`
729: @*/
730: PetscErrorCode DMCompositeAddDM(DM dmc, DM dm)
731: {
732: PetscInt n, nlocal;
733: struct DMCompositeLink *mine, *next;
734: Vec global, local;
735: DM_Composite *com = (DM_Composite *)dmc->data;
736: PetscBool flg;
738: PetscFunctionBegin;
741: PetscCall(PetscObjectTypeCompare((PetscObject)dmc, DMCOMPOSITE, &flg));
742: PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
743: next = com->next;
744: PetscCheck(!com->setup, PetscObjectComm((PetscObject)dmc), PETSC_ERR_ARG_WRONGSTATE, "Cannot add a DM once you have used the DMComposite");
746: /* create new link */
747: PetscCall(PetscNew(&mine));
748: PetscCall(PetscObjectReference((PetscObject)dm));
749: PetscCall(DMCreateGlobalVector(dm, &global)); // Not using DMGetGlobalVector() since it will stash the vector with a type decided by dm,
750: PetscCall(VecGetLocalSize(global, &n)); // while we may want to set dmc's vectype later, say via DMSetFromOptions(dmc), and we
751: PetscCall(VecDestroy(&global)); // want to propagate the type to dm.
752: PetscCall(DMCreateLocalVector(dm, &local)); // Not using DMGetLocalVector(), same reason as above.
753: PetscCall(VecGetSize(local, &nlocal));
754: PetscCall(VecDestroy(&local));
756: mine->n = n;
757: mine->nlocal = nlocal;
758: mine->dm = dm;
759: mine->next = NULL;
760: com->n += n;
761: com->nghost += nlocal;
763: /* add to end of list */
764: if (!next) com->next = mine;
765: else {
766: while (next->next) next = next->next;
767: next->next = mine;
768: }
769: com->nDM++;
770: com->nmine++;
771: PetscFunctionReturn(PETSC_SUCCESS);
772: }
774: #include <petscdraw.h>
775: PETSC_SINGLE_LIBRARY_INTERN PetscErrorCode VecView_MPI(Vec, PetscViewer);
777: static PetscErrorCode VecView_DMComposite(Vec gvec, PetscViewer viewer)
778: {
779: DM dm;
780: struct DMCompositeLink *next;
781: PetscBool isdraw;
782: DM_Composite *com;
784: PetscFunctionBegin;
785: PetscCall(VecGetDM(gvec, &dm));
786: PetscCheck(dm, PetscObjectComm((PetscObject)gvec), PETSC_ERR_ARG_WRONG, "Vector not generated from a DMComposite");
787: com = (DM_Composite *)dm->data;
788: next = com->next;
790: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw));
791: if (!isdraw) {
792: /* do I really want to call this? */
793: PetscCall(VecView_MPI(gvec, viewer));
794: } else {
795: PetscInt cnt = 0;
797: /* loop over packed objects, handling one at a time */
798: while (next) {
799: Vec vec;
800: const PetscScalar *array;
801: PetscInt bs;
803: /* Should use VecGetSubVector() eventually, but would need to forward the DM for that to work */
804: PetscCall(DMGetGlobalVector(next->dm, &vec));
805: PetscCall(VecGetArrayRead(gvec, &array));
806: PetscCall(VecPlaceArray(vec, (PetscScalar *)array + next->rstart));
807: PetscCall(VecRestoreArrayRead(gvec, &array));
808: PetscCall(VecView(vec, viewer));
809: PetscCall(VecResetArray(vec));
810: PetscCall(VecGetBlockSize(vec, &bs));
811: PetscCall(DMRestoreGlobalVector(next->dm, &vec));
812: PetscCall(PetscViewerDrawBaseAdd(viewer, bs));
813: cnt += bs;
814: next = next->next;
815: }
816: PetscCall(PetscViewerDrawBaseAdd(viewer, -cnt));
817: }
818: PetscFunctionReturn(PETSC_SUCCESS);
819: }
821: static PetscErrorCode DMCreateGlobalVector_Composite(DM dm, Vec *gvec)
822: {
823: DM_Composite *com = (DM_Composite *)dm->data;
825: PetscFunctionBegin;
827: PetscCall(DMSetFromOptions(dm));
828: PetscCall(DMSetUp(dm));
829: PetscCall(VecCreate(PetscObjectComm((PetscObject)dm), gvec));
830: PetscCall(VecSetType(*gvec, dm->vectype));
831: PetscCall(VecSetSizes(*gvec, com->n, com->N));
832: PetscCall(VecSetDM(*gvec, dm));
833: PetscCall(VecSetOperation(*gvec, VECOP_VIEW, (PetscErrorCodeFn *)VecView_DMComposite));
834: PetscFunctionReturn(PETSC_SUCCESS);
835: }
837: static PetscErrorCode DMCreateLocalVector_Composite(DM dm, Vec *lvec)
838: {
839: DM_Composite *com = (DM_Composite *)dm->data;
841: PetscFunctionBegin;
843: if (!com->setup) {
844: PetscCall(DMSetFromOptions(dm));
845: PetscCall(DMSetUp(dm));
846: }
847: PetscCall(VecCreate(PETSC_COMM_SELF, lvec));
848: PetscCall(VecSetType(*lvec, dm->vectype));
849: PetscCall(VecSetSizes(*lvec, com->nghost, PETSC_DECIDE));
850: PetscCall(VecSetDM(*lvec, dm));
851: PetscFunctionReturn(PETSC_SUCCESS);
852: }
854: /*@C
855: DMCompositeGetISLocalToGlobalMappings - gets an `ISLocalToGlobalMapping` for each `DM` in the `DMCOMPOSITE`, maps to the composite global space
857: Collective; No Fortran Support
859: Input Parameter:
860: . dm - the `DMCOMPOSITE` object
862: Output Parameter:
863: . ltogs - the individual mappings for each packed vector. Note that this includes
864: all the ghost points that individual ghosted `DMDA` may have.
866: Level: advanced
868: Note:
869: Each entry of `ltogs` should be destroyed with `ISLocalToGlobalMappingDestroy()`, `ltogs` should be freed with `PetscFree()`.
871: .seealso: `DMCOMPOSITE`, `DM`, `DMDestroy()`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`,
872: `DMCompositeGather()`, `DMCompositeCreate()`, `DMCompositeGetAccess()`, `DMCompositeScatter()`,
873: `DMCompositeGetLocalVectors()`, `DMCompositeRestoreLocalVectors()`, `DMCompositeGetEntries()`
874: @*/
875: PetscErrorCode DMCompositeGetISLocalToGlobalMappings(DM dm, ISLocalToGlobalMapping *ltogs[])
876: {
877: PetscInt i, *idx, n, cnt;
878: struct DMCompositeLink *next;
879: PetscMPIInt rank;
880: DM_Composite *com = (DM_Composite *)dm->data;
881: PetscBool flg;
883: PetscFunctionBegin;
885: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
886: PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
887: PetscCall(DMSetUp(dm));
888: PetscCall(PetscMalloc1(com->nDM, ltogs));
889: next = com->next;
890: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
892: /* loop over packed objects, handling one at a time */
893: cnt = 0;
894: while (next) {
895: ISLocalToGlobalMapping ltog;
896: PetscMPIInt size;
897: const PetscInt *suboff, *indices;
898: Vec global;
900: /* Get sub-DM global indices for each local dof */
901: PetscCall(DMGetLocalToGlobalMapping(next->dm, <og));
902: PetscCall(ISLocalToGlobalMappingGetSize(ltog, &n));
903: PetscCall(ISLocalToGlobalMappingGetIndices(ltog, &indices));
904: PetscCall(PetscMalloc1(n, &idx));
906: /* Get the offsets for the sub-DM global vector */
907: PetscCall(DMGetGlobalVector(next->dm, &global));
908: PetscCall(VecGetOwnershipRanges(global, &suboff));
909: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)global), &size));
911: /* Shift the sub-DM definition of the global space to the composite global space */
912: for (i = 0; i < n; i++) {
913: PetscInt subi = indices[i], lo = 0, hi = size, t;
914: /* There's no consensus on what a negative index means,
915: except for skipping when setting the values in vectors and matrices */
916: if (subi < 0) {
917: idx[i] = subi - next->grstarts[rank];
918: continue;
919: }
920: /* Binary search to find which rank owns subi */
921: while (hi - lo > 1) {
922: t = lo + (hi - lo) / 2;
923: if (suboff[t] > subi) hi = t;
924: else lo = t;
925: }
926: idx[i] = subi - suboff[lo] + next->grstarts[lo];
927: }
928: PetscCall(ISLocalToGlobalMappingRestoreIndices(ltog, &indices));
929: PetscCall(ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)dm), 1, n, idx, PETSC_OWN_POINTER, &(*ltogs)[cnt]));
930: PetscCall(DMRestoreGlobalVector(next->dm, &global));
931: next = next->next;
932: cnt++;
933: }
934: PetscFunctionReturn(PETSC_SUCCESS);
935: }
937: /*@C
938: DMCompositeGetLocalISs - Gets index sets for each component of a composite local vector
940: Not Collective; No Fortran Support
942: Input Parameter:
943: . dm - the `DMCOMPOSITE`
945: Output Parameter:
946: . is - array of serial index sets for each component of the `DMCOMPOSITE`
948: Level: intermediate
950: Notes:
951: At present, a composite local vector does not normally exist. This function is used to provide index sets for
952: `MatGetLocalSubMatrix()`. In the future, the scatters for each entry in the `DMCOMPOSITE` may be merged into a single
953: scatter to a composite local vector. The user should not typically need to know which is being done.
955: To get the composite global indices at all local points (including ghosts), use `DMCompositeGetISLocalToGlobalMappings()`.
957: To get index sets for pieces of the composite global vector, use `DMCompositeGetGlobalISs()`.
959: Each returned `IS` should be destroyed with `ISDestroy()`, the array should be freed with `PetscFree()`.
961: Fortran Note:
962: Use `DMCompositeRestoreLocalISs()` to release the `is`.
964: .seealso: `DMCOMPOSITE`, `DM`, `DMCompositeGetGlobalISs()`, `DMCompositeGetISLocalToGlobalMappings()`, `MatGetLocalSubMatrix()`,
965: `MatCreateLocalRef()`, `DMCompositeGetNumberDM()`
966: @*/
967: PetscErrorCode DMCompositeGetLocalISs(DM dm, IS *is[])
968: {
969: DM_Composite *com = (DM_Composite *)dm->data;
970: struct DMCompositeLink *link;
971: PetscInt cnt, start;
972: PetscBool flg;
974: PetscFunctionBegin;
976: PetscAssertPointer(is, 2);
977: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
978: PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
979: PetscCall(PetscMalloc1(com->nmine, is));
980: for (cnt = 0, start = 0, link = com->next; link; start += link->nlocal, cnt++, link = link->next) {
981: PetscInt bs;
982: PetscCall(ISCreateStride(PETSC_COMM_SELF, link->nlocal, start, 1, &(*is)[cnt]));
983: PetscCall(DMGetBlockSize(link->dm, &bs));
984: PetscCall(ISSetBlockSize((*is)[cnt], bs));
985: }
986: PetscFunctionReturn(PETSC_SUCCESS);
987: }
989: /*@C
990: DMCompositeGetGlobalISs - Gets the index sets for each composed object in a `DMCOMPOSITE`
992: Collective
994: Input Parameter:
995: . dm - the `DMCOMPOSITE` object
997: Output Parameter:
998: . is - the array of index sets
1000: Level: advanced
1002: Notes:
1003: The `is` entries should be destroyed with `ISDestroy()`, `is` should be freed with `PetscFree()`
1005: These could be used to extract a subset of vector entries for a "multi-physics" preconditioner
1007: Use `DMCompositeGetLocalISs()` for index sets in the packed local numbering, and
1008: `DMCompositeGetISLocalToGlobalMappings()` for to map local sub-`DM` (including ghost) indices to packed global
1009: indices.
1011: Fortran Note:
1012: Use `DMCompositeRestoreGlobalISs()` to release the `is`.
1014: .seealso: `DMCOMPOSITE`, `DM`, `DMDestroy()`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`,
1015: `DMCompositeGather()`, `DMCompositeCreate()`, `DMCompositeGetAccess()`, `DMCompositeScatter()`,
1016: `DMCompositeGetLocalVectors()`, `DMCompositeRestoreLocalVectors()`, `DMCompositeGetEntries()`
1017: @*/
1018: PetscErrorCode DMCompositeGetGlobalISs(DM dm, IS *is[])
1019: {
1020: PetscInt cnt = 0;
1021: struct DMCompositeLink *next;
1022: PetscMPIInt rank;
1023: DM_Composite *com = (DM_Composite *)dm->data;
1024: PetscBool flg;
1026: PetscFunctionBegin;
1028: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
1029: PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
1030: PetscCheck(dm->setupcalled, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "Must call DMSetUp() before");
1031: PetscCall(PetscMalloc1(com->nDM, is));
1032: next = com->next;
1033: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
1035: /* loop over packed objects, handling one at a time */
1036: while (next) {
1037: PetscDS prob;
1039: PetscCall(ISCreateStride(PetscObjectComm((PetscObject)dm), next->n, next->grstart, 1, &(*is)[cnt]));
1040: PetscCall(DMGetDS(dm, &prob));
1041: if (prob) {
1042: MatNullSpace space;
1043: Mat pmat;
1044: PetscObject disc;
1045: PetscInt Nf;
1047: PetscCall(PetscDSGetNumFields(prob, &Nf));
1048: if (cnt < Nf) {
1049: PetscCall(PetscDSGetDiscretization(prob, cnt, &disc));
1050: PetscCall(PetscObjectQuery(disc, "nullspace", (PetscObject *)&space));
1051: if (space) PetscCall(PetscObjectCompose((PetscObject)(*is)[cnt], "nullspace", (PetscObject)space));
1052: PetscCall(PetscObjectQuery(disc, "nearnullspace", (PetscObject *)&space));
1053: if (space) PetscCall(PetscObjectCompose((PetscObject)(*is)[cnt], "nearnullspace", (PetscObject)space));
1054: PetscCall(PetscObjectQuery(disc, "pmat", (PetscObject *)&pmat));
1055: if (pmat) PetscCall(PetscObjectCompose((PetscObject)(*is)[cnt], "pmat", (PetscObject)pmat));
1056: }
1057: }
1058: cnt++;
1059: next = next->next;
1060: }
1061: PetscFunctionReturn(PETSC_SUCCESS);
1062: }
1064: static PetscErrorCode DMCreateFieldIS_Composite(DM dm, PetscInt *numFields, char ***fieldNames, IS **fields)
1065: {
1066: PetscInt nDM;
1067: DM *dms;
1068: PetscInt i;
1070: PetscFunctionBegin;
1071: PetscCall(DMCompositeGetNumberDM(dm, &nDM));
1072: if (numFields) *numFields = nDM;
1073: PetscCall(DMCompositeGetGlobalISs(dm, fields));
1074: if (fieldNames) {
1075: PetscCall(PetscMalloc1(nDM, &dms));
1076: PetscCall(PetscMalloc1(nDM, fieldNames));
1077: PetscCall(DMCompositeGetEntriesArray(dm, dms));
1078: for (i = 0; i < nDM; i++) {
1079: char buf[256];
1080: const char *splitname;
1082: /* Split naming precedence: object name, prefix, number */
1083: splitname = ((PetscObject)dm)->name;
1084: if (!splitname) {
1085: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)dms[i], &splitname));
1086: if (splitname) {
1087: size_t len;
1088: PetscCall(PetscStrncpy(buf, splitname, sizeof(buf)));
1089: buf[sizeof(buf) - 1] = 0;
1090: PetscCall(PetscStrlen(buf, &len));
1091: if (buf[len - 1] == '_') buf[len - 1] = 0; /* Remove trailing underscore if it was used */
1092: splitname = buf;
1093: }
1094: }
1095: if (!splitname) {
1096: PetscCall(PetscSNPrintf(buf, sizeof(buf), "%" PetscInt_FMT, i));
1097: splitname = buf;
1098: }
1099: PetscCall(PetscStrallocpy(splitname, &(*fieldNames)[i]));
1100: }
1101: PetscCall(PetscFree(dms));
1102: }
1103: PetscFunctionReturn(PETSC_SUCCESS);
1104: }
1106: /*
1107: This could take over from DMCreateFieldIS(), as it is more general,
1108: making DMCreateFieldIS() a special case -- calling with dmlist == NULL;
1109: At this point it's probably best to be less intrusive, however.
1110: */
1111: static PetscErrorCode DMCreateFieldDecomposition_Composite(DM dm, PetscInt *len, char ***namelist, IS **islist, DM **dmlist)
1112: {
1113: PetscInt nDM;
1114: PetscInt i;
1116: PetscFunctionBegin;
1117: PetscCall(DMCreateFieldIS_Composite(dm, len, namelist, islist));
1118: if (dmlist) {
1119: PetscCall(DMCompositeGetNumberDM(dm, &nDM));
1120: PetscCall(PetscMalloc1(nDM, dmlist));
1121: PetscCall(DMCompositeGetEntriesArray(dm, *dmlist));
1122: for (i = 0; i < nDM; i++) PetscCall(PetscObjectReference((PetscObject)((*dmlist)[i])));
1123: }
1124: PetscFunctionReturn(PETSC_SUCCESS);
1125: }
1127: /*@C
1128: DMCompositeGetLocalVectors - Gets local vectors for each part of a `DMCOMPOSITE`
1129: Use `DMCompositeRestoreLocalVectors()` to return them.
1131: Not Collective; No Fortran Support
1133: Input Parameter:
1134: . dm - the `DMCOMPOSITE` object
1136: Output Parameter:
1137: . ... - the individual sequential `Vec`s
1139: Level: advanced
1141: .seealso: `DMCOMPOSITE`, `DM`, `DMDestroy()`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`,
1142: `DMCompositeGather()`, `DMCompositeCreate()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeGetAccess()`,
1143: `DMCompositeRestoreLocalVectors()`, `DMCompositeScatter()`, `DMCompositeGetEntries()`
1144: @*/
1145: PetscErrorCode DMCompositeGetLocalVectors(DM dm, ...)
1146: {
1147: va_list Argp;
1148: struct DMCompositeLink *next;
1149: DM_Composite *com = (DM_Composite *)dm->data;
1150: PetscBool flg;
1152: PetscFunctionBegin;
1154: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
1155: PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
1156: next = com->next;
1157: /* loop over packed objects, handling one at a time */
1158: va_start(Argp, dm);
1159: while (next) {
1160: Vec *vec;
1161: vec = va_arg(Argp, Vec *);
1162: if (vec) PetscCall(DMGetLocalVector(next->dm, vec));
1163: next = next->next;
1164: }
1165: va_end(Argp);
1166: PetscFunctionReturn(PETSC_SUCCESS);
1167: }
1169: /*@C
1170: DMCompositeRestoreLocalVectors - Restores local vectors for each part of a `DMCOMPOSITE`
1172: Not Collective; No Fortran Support
1174: Input Parameter:
1175: . dm - the `DMCOMPOSITE` object
1177: Output Parameter:
1178: . ... - the individual sequential `Vec`
1180: Level: advanced
1182: .seealso: `DMCOMPOSITE`, `DM`, `DMDestroy()`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`,
1183: `DMCompositeGather()`, `DMCompositeCreate()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeGetAccess()`,
1184: `DMCompositeGetLocalVectors()`, `DMCompositeScatter()`, `DMCompositeGetEntries()`
1185: @*/
1186: PetscErrorCode DMCompositeRestoreLocalVectors(DM dm, ...)
1187: {
1188: va_list Argp;
1189: struct DMCompositeLink *next;
1190: DM_Composite *com = (DM_Composite *)dm->data;
1191: PetscBool flg;
1193: PetscFunctionBegin;
1195: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
1196: PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
1197: next = com->next;
1198: /* loop over packed objects, handling one at a time */
1199: va_start(Argp, dm);
1200: while (next) {
1201: Vec *vec;
1202: vec = va_arg(Argp, Vec *);
1203: if (vec) PetscCall(DMRestoreLocalVector(next->dm, vec));
1204: next = next->next;
1205: }
1206: va_end(Argp);
1207: PetscFunctionReturn(PETSC_SUCCESS);
1208: }
1210: /*@C
1211: DMCompositeGetEntries - Gets the `DM` for each entry in a `DMCOMPOSITE`.
1213: Not Collective
1215: Input Parameter:
1216: . dm - the `DMCOMPOSITE` object
1218: Output Parameter:
1219: . ... - the individual entries `DM`
1221: Level: advanced
1223: Fortran Notes:
1224: Use `DMCompositeGetEntriesArray()`
1226: .seealso: `DMCOMPOSITE`, `DM`, `DMDestroy()`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`, `DMCompositeGetEntriesArray()`
1227: `DMCompositeGather()`, `DMCompositeCreate()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeGetAccess()`,
1228: `DMCompositeRestoreLocalVectors()`, `DMCompositeGetLocalVectors()`, `DMCompositeScatter()`
1229: @*/
1230: PetscErrorCode DMCompositeGetEntries(DM dm, ...)
1231: {
1232: va_list Argp;
1233: struct DMCompositeLink *next;
1234: DM_Composite *com = (DM_Composite *)dm->data;
1235: PetscBool flg;
1237: PetscFunctionBegin;
1239: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
1240: PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
1241: next = com->next;
1242: /* loop over packed objects, handling one at a time */
1243: va_start(Argp, dm);
1244: while (next) {
1245: DM *dmn;
1246: dmn = va_arg(Argp, DM *);
1247: if (dmn) *dmn = next->dm;
1248: next = next->next;
1249: }
1250: va_end(Argp);
1251: PetscFunctionReturn(PETSC_SUCCESS);
1252: }
1254: /*@
1255: DMCompositeGetEntriesArray - Gets the DM for each entry in a `DMCOMPOSITE`
1257: Not Collective
1259: Input Parameter:
1260: . dm - the `DMCOMPOSITE` object
1262: Output Parameter:
1263: . dms - array of sufficient length (see `DMCompositeGetNumberDM()`) to hold the individual `DM`
1265: Level: advanced
1267: .seealso: `DMCOMPOSITE`, `DM`, `DMDestroy()`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`, `DMCompositeGetEntries()`
1268: `DMCompositeGather()`, `DMCompositeCreate()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeGetAccess()`,
1269: `DMCompositeRestoreLocalVectors()`, `DMCompositeGetLocalVectors()`, `DMCompositeScatter()`
1270: @*/
1271: PetscErrorCode DMCompositeGetEntriesArray(DM dm, DM dms[])
1272: {
1273: struct DMCompositeLink *next;
1274: DM_Composite *com = (DM_Composite *)dm->data;
1275: PetscInt i;
1276: PetscBool flg;
1278: PetscFunctionBegin;
1280: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
1281: PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
1282: /* loop over packed objects, handling one at a time */
1283: for (next = com->next, i = 0; next; next = next->next, i++) dms[i] = next->dm;
1284: PetscFunctionReturn(PETSC_SUCCESS);
1285: }
1287: typedef struct {
1288: DM dm;
1289: PetscViewer *subv;
1290: Vec *vecs;
1291: } GLVisViewerCtx;
1293: static PetscErrorCode DestroyGLVisViewerCtx_Private(PetscCtxRt vctx)
1294: {
1295: GLVisViewerCtx *ctx = *(GLVisViewerCtx **)vctx;
1296: PetscInt i, n;
1298: PetscFunctionBegin;
1299: PetscCall(DMCompositeGetNumberDM(ctx->dm, &n));
1300: for (i = 0; i < n; i++) PetscCall(PetscViewerDestroy(&ctx->subv[i]));
1301: PetscCall(PetscFree2(ctx->subv, ctx->vecs));
1302: PetscCall(DMDestroy(&ctx->dm));
1303: PetscCall(PetscFree(ctx));
1304: PetscFunctionReturn(PETSC_SUCCESS);
1305: }
1307: static PetscErrorCode DMCompositeSampleGLVisFields_Private(PetscObject oX, PetscInt nf, PetscObject oXfield[], void *vctx)
1308: {
1309: Vec X = (Vec)oX;
1310: GLVisViewerCtx *ctx = (GLVisViewerCtx *)vctx;
1311: PetscInt i, n, cumf;
1313: PetscFunctionBegin;
1314: PetscCall(DMCompositeGetNumberDM(ctx->dm, &n));
1315: PetscCall(DMCompositeGetAccessArray(ctx->dm, X, n, NULL, ctx->vecs));
1316: for (i = 0, cumf = 0; i < n; i++) {
1317: PetscErrorCode (*g2l)(PetscObject, PetscInt, PetscObject[], void *);
1318: void *fctx;
1319: PetscInt nfi;
1321: PetscCall(PetscViewerGLVisGetFields_Internal(ctx->subv[i], &nfi, NULL, NULL, &g2l, NULL, &fctx));
1322: if (!nfi) continue;
1323: if (g2l) PetscCall((*g2l)((PetscObject)ctx->vecs[i], nfi, oXfield + cumf, fctx));
1324: else PetscCall(VecCopy(ctx->vecs[i], (Vec)oXfield[cumf]));
1325: cumf += nfi;
1326: }
1327: PetscCall(DMCompositeRestoreAccessArray(ctx->dm, X, n, NULL, ctx->vecs));
1328: PetscFunctionReturn(PETSC_SUCCESS);
1329: }
1331: static PetscErrorCode DMSetUpGLVisViewer_Composite(PetscObject odm, PetscViewer viewer)
1332: {
1333: DM dm = (DM)odm, *dms;
1334: Vec *Ufds;
1335: GLVisViewerCtx *ctx;
1336: PetscInt i, n, tnf, *sdim;
1337: char **fecs;
1339: PetscFunctionBegin;
1340: PetscCall(PetscNew(&ctx));
1341: PetscCall(PetscObjectReference((PetscObject)dm));
1342: ctx->dm = dm;
1343: PetscCall(DMCompositeGetNumberDM(dm, &n));
1344: PetscCall(PetscMalloc1(n, &dms));
1345: PetscCall(DMCompositeGetEntriesArray(dm, dms));
1346: PetscCall(PetscMalloc2(n, &ctx->subv, n, &ctx->vecs));
1347: for (i = 0, tnf = 0; i < n; i++) {
1348: PetscInt nf;
1350: PetscCall(PetscViewerCreate(PetscObjectComm(odm), &ctx->subv[i]));
1351: PetscCall(PetscViewerSetType(ctx->subv[i], PETSCVIEWERGLVIS));
1352: PetscCall(PetscViewerGLVisSetDM_Internal(ctx->subv[i], (PetscObject)dms[i]));
1353: PetscCall(PetscViewerGLVisGetFields_Internal(ctx->subv[i], &nf, NULL, NULL, NULL, NULL, NULL));
1354: tnf += nf;
1355: }
1356: PetscCall(PetscFree(dms));
1357: PetscCall(PetscMalloc3(tnf, &fecs, tnf, &sdim, tnf, &Ufds));
1358: for (i = 0, tnf = 0; i < n; i++) {
1359: PetscInt *sd, nf, f;
1360: const char **fec;
1361: Vec *Uf;
1363: PetscCall(PetscViewerGLVisGetFields_Internal(ctx->subv[i], &nf, &fec, &sd, NULL, (PetscObject **)&Uf, NULL));
1364: for (f = 0; f < nf; f++) {
1365: PetscCall(PetscStrallocpy(fec[f], &fecs[tnf + f]));
1366: Ufds[tnf + f] = Uf[f];
1367: sdim[tnf + f] = sd[f];
1368: }
1369: tnf += nf;
1370: }
1371: PetscCall(PetscViewerGLVisSetFields(viewer, tnf, (const char **)fecs, sdim, DMCompositeSampleGLVisFields_Private, (PetscObject *)Ufds, ctx, DestroyGLVisViewerCtx_Private));
1372: for (i = 0; i < tnf; i++) PetscCall(PetscFree(fecs[i]));
1373: PetscCall(PetscFree3(fecs, sdim, Ufds));
1374: PetscFunctionReturn(PETSC_SUCCESS);
1375: }
1377: static PetscErrorCode DMRefine_Composite(DM dmi, MPI_Comm comm, DM *fine)
1378: {
1379: struct DMCompositeLink *next;
1380: DM_Composite *com = (DM_Composite *)dmi->data;
1381: DM dm;
1383: PetscFunctionBegin;
1385: if (comm == MPI_COMM_NULL) PetscCall(PetscObjectGetComm((PetscObject)dmi, &comm));
1386: PetscCall(DMSetUp(dmi));
1387: next = com->next;
1388: PetscCall(DMCompositeCreate(comm, fine));
1390: /* loop over packed objects, handling one at a time */
1391: while (next) {
1392: PetscCall(DMRefine(next->dm, comm, &dm));
1393: PetscCall(DMCompositeAddDM(*fine, dm));
1394: PetscCall(PetscObjectDereference((PetscObject)dm));
1395: next = next->next;
1396: }
1397: PetscFunctionReturn(PETSC_SUCCESS);
1398: }
1400: static PetscErrorCode DMCoarsen_Composite(DM dmi, MPI_Comm comm, DM *fine)
1401: {
1402: struct DMCompositeLink *next;
1403: DM_Composite *com = (DM_Composite *)dmi->data;
1404: DM dm;
1406: PetscFunctionBegin;
1408: PetscCall(DMSetUp(dmi));
1409: if (comm == MPI_COMM_NULL) PetscCall(PetscObjectGetComm((PetscObject)dmi, &comm));
1410: next = com->next;
1411: PetscCall(DMCompositeCreate(comm, fine));
1413: /* loop over packed objects, handling one at a time */
1414: while (next) {
1415: PetscCall(DMCoarsen(next->dm, comm, &dm));
1416: PetscCall(DMCompositeAddDM(*fine, dm));
1417: PetscCall(PetscObjectDereference((PetscObject)dm));
1418: next = next->next;
1419: }
1420: PetscFunctionReturn(PETSC_SUCCESS);
1421: }
1423: static PetscErrorCode DMCreateInterpolation_Composite(DM coarse, DM fine, Mat *A, Vec *v)
1424: {
1425: PetscInt m, n, M, N, nDM, i;
1426: struct DMCompositeLink *nextc;
1427: struct DMCompositeLink *nextf;
1428: Vec gcoarse, gfine, *vecs;
1429: DM_Composite *comcoarse = (DM_Composite *)coarse->data;
1430: DM_Composite *comfine = (DM_Composite *)fine->data;
1431: Mat *mats;
1433: PetscFunctionBegin;
1436: PetscCall(DMSetUp(coarse));
1437: PetscCall(DMSetUp(fine));
1438: /* use global vectors only for determining matrix layout */
1439: PetscCall(DMGetGlobalVector(coarse, &gcoarse));
1440: PetscCall(DMGetGlobalVector(fine, &gfine));
1441: PetscCall(VecGetLocalSize(gcoarse, &n));
1442: PetscCall(VecGetLocalSize(gfine, &m));
1443: PetscCall(VecGetSize(gcoarse, &N));
1444: PetscCall(VecGetSize(gfine, &M));
1445: PetscCall(DMRestoreGlobalVector(coarse, &gcoarse));
1446: PetscCall(DMRestoreGlobalVector(fine, &gfine));
1448: nDM = comfine->nDM;
1449: 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);
1450: PetscCall(PetscCalloc1(nDM * nDM, &mats));
1451: if (v) PetscCall(PetscCalloc1(nDM, &vecs));
1453: /* loop over packed objects, handling one at a time */
1454: for (nextc = comcoarse->next, nextf = comfine->next, i = 0; nextc; nextc = nextc->next, nextf = nextf->next, i++) {
1455: if (!v) PetscCall(DMCreateInterpolation(nextc->dm, nextf->dm, &mats[i * nDM + i], NULL));
1456: else PetscCall(DMCreateInterpolation(nextc->dm, nextf->dm, &mats[i * nDM + i], &vecs[i]));
1457: }
1458: PetscCall(MatCreateNest(PetscObjectComm((PetscObject)fine), nDM, NULL, nDM, NULL, mats, A));
1459: if (v) PetscCall(VecCreateNest(PetscObjectComm((PetscObject)fine), nDM, NULL, vecs, v));
1460: for (i = 0; i < nDM * nDM; i++) PetscCall(MatDestroy(&mats[i]));
1461: PetscCall(PetscFree(mats));
1462: if (v) {
1463: for (i = 0; i < nDM; i++) PetscCall(VecDestroy(&vecs[i]));
1464: PetscCall(PetscFree(vecs));
1465: }
1466: PetscFunctionReturn(PETSC_SUCCESS);
1467: }
1469: static PetscErrorCode DMGetLocalToGlobalMapping_Composite(DM dm)
1470: {
1471: DM_Composite *com = (DM_Composite *)dm->data;
1472: ISLocalToGlobalMapping *ltogs;
1473: PetscInt i;
1475: PetscFunctionBegin;
1476: /* Set the ISLocalToGlobalMapping on the new matrix */
1477: PetscCall(DMCompositeGetISLocalToGlobalMappings(dm, <ogs));
1478: PetscCall(ISLocalToGlobalMappingConcatenate(PetscObjectComm((PetscObject)dm), com->nDM, ltogs, &dm->ltogmap));
1479: for (i = 0; i < com->nDM; i++) PetscCall(ISLocalToGlobalMappingDestroy(<ogs[i]));
1480: PetscCall(PetscFree(ltogs));
1481: PetscFunctionReturn(PETSC_SUCCESS);
1482: }
1484: static PetscErrorCode DMCreateColoring_Composite(DM dm, ISColoringType ctype, ISColoring *coloring)
1485: {
1486: PetscInt n, i, cnt;
1487: ISColoringValue *colors;
1488: PetscBool dense = PETSC_FALSE;
1489: ISColoringValue maxcol = 0;
1490: DM_Composite *com = (DM_Composite *)dm->data;
1492: PetscFunctionBegin;
1494: PetscCheck(ctype != IS_COLORING_LOCAL, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Only global coloring supported");
1495: PetscCheck(ctype == IS_COLORING_GLOBAL, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Unknown ISColoringType");
1496: n = com->n;
1497: PetscCall(PetscMalloc1(n, &colors)); /* freed in ISColoringDestroy() */
1499: PetscCall(PetscOptionsGetBool(((PetscObject)dm)->options, ((PetscObject)dm)->prefix, "-dmcomposite_dense_jacobian", &dense, NULL));
1500: if (dense) {
1501: PetscCall(ISColoringValueCast(com->N, &maxcol));
1502: for (i = 0; i < n; i++) PetscCall(ISColoringValueCast(com->rstart + i, colors + i));
1503: } else {
1504: struct DMCompositeLink *next = com->next;
1505: PetscMPIInt rank;
1507: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
1508: cnt = 0;
1509: while (next) {
1510: ISColoring lcoloring;
1512: PetscCall(DMCreateColoring(next->dm, IS_COLORING_GLOBAL, &lcoloring));
1513: for (i = 0; i < lcoloring->N; i++) PetscCall(ISColoringValueCast(maxcol + lcoloring->colors[i], colors + cnt++));
1514: maxcol += lcoloring->n;
1515: PetscCall(ISColoringDestroy(&lcoloring));
1516: next = next->next;
1517: }
1518: }
1519: PetscCall(ISColoringCreate(PetscObjectComm((PetscObject)dm), maxcol, n, colors, PETSC_OWN_POINTER, coloring));
1520: PetscFunctionReturn(PETSC_SUCCESS);
1521: }
1523: static PetscErrorCode DMGlobalToLocalBegin_Composite(DM dm, Vec gvec, InsertMode mode, Vec lvec)
1524: {
1525: struct DMCompositeLink *next;
1526: const PetscScalar *garray, *garrayhead;
1527: PetscScalar *larray, *larrayhead;
1528: DM_Composite *com = (DM_Composite *)dm->data;
1530: PetscFunctionBegin;
1534: if (!com->setup) PetscCall(DMSetUp(dm));
1536: PetscCall(VecGetArrayRead(gvec, &garray));
1537: garrayhead = garray;
1538: PetscCall(VecGetArray(lvec, &larray));
1539: larrayhead = larray;
1541: /* loop over packed objects, handling one at a time */
1542: next = com->next;
1543: while (next) {
1544: Vec local, global;
1545: PetscInt N;
1547: PetscCall(DMGetGlobalVector(next->dm, &global));
1548: PetscCall(VecGetLocalSize(global, &N));
1549: PetscCall(VecPlaceArray(global, garrayhead));
1550: PetscCall(DMGetLocalVector(next->dm, &local));
1551: PetscCall(VecPlaceArray(local, larrayhead));
1552: PetscCall(DMGlobalToLocalBegin(next->dm, global, mode, local));
1553: PetscCall(DMGlobalToLocalEnd(next->dm, global, mode, local));
1554: PetscCall(VecResetArray(global));
1555: PetscCall(VecResetArray(local));
1556: PetscCall(DMRestoreGlobalVector(next->dm, &global));
1557: PetscCall(DMRestoreLocalVector(next->dm, &local));
1559: larrayhead += next->nlocal;
1560: garrayhead += next->n;
1561: next = next->next;
1562: }
1563: PetscCall(VecRestoreArrayRead(gvec, &garray));
1564: PetscCall(VecRestoreArray(lvec, &larray));
1565: PetscFunctionReturn(PETSC_SUCCESS);
1566: }
1568: static PetscErrorCode DMGlobalToLocalEnd_Composite(DM dm, Vec gvec, InsertMode mode, Vec lvec)
1569: {
1570: PetscFunctionBegin;
1574: PetscFunctionReturn(PETSC_SUCCESS);
1575: }
1577: static PetscErrorCode DMLocalToGlobalBegin_Composite(DM dm, Vec lvec, InsertMode mode, Vec gvec)
1578: {
1579: struct DMCompositeLink *next;
1580: const PetscScalar *larray, *larrayhead;
1581: PetscScalar *garray, *garrayhead;
1582: DM_Composite *com = (DM_Composite *)dm->data;
1584: PetscFunctionBegin;
1589: if (!com->setup) PetscCall(DMSetUp(dm));
1591: PetscCall(VecGetArrayRead(lvec, &larray));
1592: larrayhead = larray;
1593: PetscCall(VecGetArray(gvec, &garray));
1594: garrayhead = garray;
1596: /* loop over packed objects, handling one at a time */
1597: next = com->next;
1598: while (next) {
1599: Vec global, local;
1601: PetscCall(DMGetLocalVector(next->dm, &local));
1602: PetscCall(VecPlaceArray(local, larrayhead));
1603: PetscCall(DMGetGlobalVector(next->dm, &global));
1604: PetscCall(VecPlaceArray(global, garrayhead));
1605: PetscCall(DMLocalToGlobalBegin(next->dm, local, mode, global));
1606: PetscCall(DMLocalToGlobalEnd(next->dm, local, mode, global));
1607: PetscCall(VecResetArray(local));
1608: PetscCall(VecResetArray(global));
1609: PetscCall(DMRestoreGlobalVector(next->dm, &global));
1610: PetscCall(DMRestoreLocalVector(next->dm, &local));
1612: garrayhead += next->n;
1613: larrayhead += next->nlocal;
1614: next = next->next;
1615: }
1617: PetscCall(VecRestoreArray(gvec, &garray));
1618: PetscCall(VecRestoreArrayRead(lvec, &larray));
1619: PetscFunctionReturn(PETSC_SUCCESS);
1620: }
1622: static PetscErrorCode DMLocalToGlobalEnd_Composite(DM dm, Vec lvec, InsertMode mode, Vec gvec)
1623: {
1624: PetscFunctionBegin;
1628: PetscFunctionReturn(PETSC_SUCCESS);
1629: }
1631: static PetscErrorCode DMLocalToLocalBegin_Composite(DM dm, Vec vec1, InsertMode mode, Vec vec2)
1632: {
1633: struct DMCompositeLink *next;
1634: const PetscScalar *array1, *array1head;
1635: PetscScalar *array2, *array2head;
1636: DM_Composite *com = (DM_Composite *)dm->data;
1638: PetscFunctionBegin;
1643: if (!com->setup) PetscCall(DMSetUp(dm));
1645: PetscCall(VecGetArrayRead(vec1, &array1));
1646: array1head = array1;
1647: PetscCall(VecGetArray(vec2, &array2));
1648: array2head = array2;
1650: /* loop over packed objects, handling one at a time */
1651: next = com->next;
1652: while (next) {
1653: Vec local1, local2;
1655: PetscCall(DMGetLocalVector(next->dm, &local1));
1656: PetscCall(VecPlaceArray(local1, array1head));
1657: PetscCall(DMGetLocalVector(next->dm, &local2));
1658: PetscCall(VecPlaceArray(local2, array2head));
1659: PetscCall(DMLocalToLocalBegin(next->dm, local1, mode, local2));
1660: PetscCall(DMLocalToLocalEnd(next->dm, local1, mode, local2));
1661: PetscCall(VecResetArray(local2));
1662: PetscCall(DMRestoreLocalVector(next->dm, &local2));
1663: PetscCall(VecResetArray(local1));
1664: PetscCall(DMRestoreLocalVector(next->dm, &local1));
1666: array1head += next->nlocal;
1667: array2head += next->nlocal;
1668: next = next->next;
1669: }
1671: PetscCall(VecRestoreArrayRead(vec1, &array1));
1672: PetscCall(VecRestoreArray(vec2, &array2));
1673: PetscFunctionReturn(PETSC_SUCCESS);
1674: }
1676: static PetscErrorCode DMLocalToLocalEnd_Composite(DM dm, Vec lvec, InsertMode mode, Vec gvec)
1677: {
1678: PetscFunctionBegin;
1682: PetscFunctionReturn(PETSC_SUCCESS);
1683: }
1685: /*MC
1686: DMCOMPOSITE = "composite" - A `DM` object that is used to manage data for a collection of `DM`
1688: Level: intermediate
1690: .seealso: `DMType`, `DM`, `DMDACreate()`, `DMCreate()`, `DMSetType()`, `DMCompositeCreate()`
1691: M*/
1693: PETSC_EXTERN PetscErrorCode DMCreate_Composite(DM p)
1694: {
1695: DM_Composite *com;
1697: PetscFunctionBegin;
1698: PetscCall(PetscNew(&com));
1699: p->data = com;
1700: com->n = 0;
1701: com->nghost = 0;
1702: com->next = NULL;
1703: com->nDM = 0;
1705: p->ops->createglobalvector = DMCreateGlobalVector_Composite;
1706: p->ops->createlocalvector = DMCreateLocalVector_Composite;
1707: p->ops->getlocaltoglobalmapping = DMGetLocalToGlobalMapping_Composite;
1708: p->ops->createfieldis = DMCreateFieldIS_Composite;
1709: p->ops->createfielddecomposition = DMCreateFieldDecomposition_Composite;
1710: p->ops->refine = DMRefine_Composite;
1711: p->ops->coarsen = DMCoarsen_Composite;
1712: p->ops->createinterpolation = DMCreateInterpolation_Composite;
1713: p->ops->creatematrix = DMCreateMatrix_Composite;
1714: p->ops->getcoloring = DMCreateColoring_Composite;
1715: p->ops->globaltolocalbegin = DMGlobalToLocalBegin_Composite;
1716: p->ops->globaltolocalend = DMGlobalToLocalEnd_Composite;
1717: p->ops->localtoglobalbegin = DMLocalToGlobalBegin_Composite;
1718: p->ops->localtoglobalend = DMLocalToGlobalEnd_Composite;
1719: p->ops->localtolocalbegin = DMLocalToLocalBegin_Composite;
1720: p->ops->localtolocalend = DMLocalToLocalEnd_Composite;
1721: p->ops->destroy = DMDestroy_Composite;
1722: p->ops->view = DMView_Composite;
1723: p->ops->setup = DMSetUp_Composite;
1725: PetscCall(PetscObjectComposeFunction((PetscObject)p, "DMSetUpGLVisViewer_C", DMSetUpGLVisViewer_Composite));
1726: PetscFunctionReturn(PETSC_SUCCESS);
1727: }
1729: /*@
1730: DMCompositeCreate - Creates a `DMCOMPOSITE`, used to generate "composite"
1731: vectors made up of several subvectors.
1733: Collective
1735: Input Parameter:
1736: . comm - the processors that will share the global vector
1738: Output Parameter:
1739: . packer - the `DMCOMPOSITE` object
1741: Level: advanced
1743: .seealso: `DMCOMPOSITE`, `DM`, `DMDestroy()`, `DMCompositeAddDM()`, `DMCompositeScatter()`, `DMCreate()`
1744: `DMCompositeGather()`, `DMCreateGlobalVector()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeGetAccess()`
1745: `DMCompositeGetLocalVectors()`, `DMCompositeRestoreLocalVectors()`, `DMCompositeGetEntries()`
1746: @*/
1747: PetscErrorCode DMCompositeCreate(MPI_Comm comm, DM *packer)
1748: {
1749: PetscFunctionBegin;
1750: PetscAssertPointer(packer, 2);
1751: PetscCall(DMCreate(comm, packer));
1752: PetscCall(DMSetType(*packer, DMCOMPOSITE));
1753: PetscFunctionReturn(PETSC_SUCCESS);
1754: }