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: /* --------------------------------------------------------------------------------------*/
100: static PetscErrorCode DMSetUp_Composite(DM dm)
101: {
102: PetscInt nprev = 0;
103: PetscMPIInt rank, size;
104: DM_Composite *com = (DM_Composite *)dm->data;
105: struct DMCompositeLink *next = com->next;
106: PetscLayout map;
108: PetscFunctionBegin;
109: PetscCheck(!com->setup, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "Packer has already been setup");
110: PetscCall(PetscLayoutCreate(PetscObjectComm((PetscObject)dm), &map));
111: PetscCall(PetscLayoutSetLocalSize(map, com->n));
112: PetscCall(PetscLayoutSetSize(map, PETSC_DETERMINE));
113: PetscCall(PetscLayoutSetBlockSize(map, 1));
114: PetscCall(PetscLayoutSetUp(map));
115: PetscCall(PetscLayoutGetSize(map, &com->N));
116: PetscCall(PetscLayoutGetRange(map, &com->rstart, NULL));
117: PetscCall(PetscLayoutDestroy(&map));
119: /* now set the rstart for each linked vector */
120: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
121: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size));
122: while (next) {
123: next->rstart = nprev;
124: nprev += next->n;
125: next->grstart = com->rstart + next->rstart;
126: PetscCall(PetscMalloc1(size, &next->grstarts));
127: PetscCallMPI(MPI_Allgather(&next->grstart, 1, MPIU_INT, next->grstarts, 1, MPIU_INT, PetscObjectComm((PetscObject)dm)));
128: next = next->next;
129: }
130: com->setup = PETSC_TRUE;
131: PetscFunctionReturn(PETSC_SUCCESS);
132: }
134: /*@
135: DMCompositeGetNumberDM - Gets the number of `DM` objects in the `DMCOMPOSITE`
136: representation.
138: Not Collective
140: Input Parameter:
141: . dm - the `DMCOMPOSITE` object
143: Output Parameter:
144: . nDM - the number of `DM`
146: Level: beginner
148: .seealso: `DMCOMPOSITE`, `DM`
149: @*/
150: PetscErrorCode DMCompositeGetNumberDM(DM dm, PetscInt *nDM)
151: {
152: DM_Composite *com = (DM_Composite *)dm->data;
153: PetscBool flg;
155: PetscFunctionBegin;
157: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
158: PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
159: *nDM = com->nDM;
160: PetscFunctionReturn(PETSC_SUCCESS);
161: }
163: /*@C
164: DMCompositeGetAccess - Allows one to access the individual packed vectors in their global
165: representation.
167: Collective
169: Input Parameters:
170: + dm - the `DMCOMPOSITE` object
171: - gvec - the global vector
173: Output Parameter:
174: . ... - the packed parallel vectors, `NULL` for those that are not needed
176: Level: advanced
178: Note:
179: Use `DMCompositeRestoreAccess()` to return the vectors when you no longer need them
181: .seealso: `DMCOMPOSITE`, `DM`, `DMCompositeGetEntries()`, `DMCompositeScatter()`
182: @*/
183: PetscErrorCode DMCompositeGetAccess(DM dm, Vec gvec, ...)
184: {
185: va_list Argp;
186: struct DMCompositeLink *next;
187: DM_Composite *com = (DM_Composite *)dm->data;
188: PetscInt readonly;
189: PetscBool flg;
191: PetscFunctionBegin;
194: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
195: PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
196: next = com->next;
197: if (!com->setup) PetscCall(DMSetUp(dm));
199: PetscCall(VecLockGet(gvec, &readonly));
200: /* loop over packed objects, handling one at a time */
201: va_start(Argp, gvec);
202: while (next) {
203: Vec *vec;
204: vec = va_arg(Argp, Vec *);
205: if (vec) {
206: PetscCall(DMGetGlobalVector(next->dm, vec));
207: if (readonly) {
208: const PetscScalar *array;
209: PetscCall(VecGetArrayRead(gvec, &array));
210: PetscCall(VecPlaceArray(*vec, array + next->rstart));
211: PetscCall(VecLockReadPush(*vec));
212: PetscCall(VecRestoreArrayRead(gvec, &array));
213: } else {
214: PetscScalar *array;
215: PetscCall(VecGetArray(gvec, &array));
216: PetscCall(VecPlaceArray(*vec, array + next->rstart));
217: PetscCall(VecRestoreArray(gvec, &array));
218: }
219: }
220: next = next->next;
221: }
222: va_end(Argp);
223: PetscFunctionReturn(PETSC_SUCCESS);
224: }
226: /*@
227: DMCompositeGetAccessArray - Allows one to access the individual packed vectors in their global
228: representation.
230: Collective
232: Input Parameters:
233: + dm - the `DMCOMPOSITE`
234: . pvec - packed vector
235: . nwanted - number of vectors wanted
236: - wanted - sorted array of integers indicating thde vectors wanted, or `NULL` to get all vectors, length `nwanted`
238: Output Parameter:
239: . vecs - array of requested global vectors (must be previously allocated and of length `nwanted`)
241: Level: advanced
243: Note:
244: Use `DMCompositeRestoreAccessArray()` to return the vectors when you no longer need them
246: .seealso: `DMCOMPOSITE`, `DM`, `DMCompositeGetAccess()`, `DMCompositeGetEntries()`, `DMCompositeScatter()`, `DMCompositeGather()`
247: @*/
248: PetscErrorCode DMCompositeGetAccessArray(DM dm, Vec pvec, PetscInt nwanted, const PetscInt wanted[], Vec vecs[])
249: {
250: struct DMCompositeLink *link;
251: PetscInt i, wnum;
252: DM_Composite *com = (DM_Composite *)dm->data;
253: PetscInt readonly;
254: PetscBool flg;
256: PetscFunctionBegin;
259: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
260: PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
261: if (!com->setup) PetscCall(DMSetUp(dm));
263: PetscCall(VecLockGet(pvec, &readonly));
264: for (i = 0, wnum = 0, link = com->next; link && wnum < nwanted; i++, link = link->next) {
265: if (!wanted || i == wanted[wnum]) {
266: Vec v;
267: PetscCall(DMGetGlobalVector(link->dm, &v));
268: if (readonly) {
269: const PetscScalar *array;
270: PetscCall(VecGetArrayRead(pvec, &array));
271: PetscCall(VecPlaceArray(v, array + link->rstart));
272: PetscCall(VecLockReadPush(v));
273: PetscCall(VecRestoreArrayRead(pvec, &array));
274: } else {
275: PetscScalar *array;
276: PetscCall(VecGetArray(pvec, &array));
277: PetscCall(VecPlaceArray(v, array + link->rstart));
278: PetscCall(VecRestoreArray(pvec, &array));
279: }
280: vecs[wnum++] = v;
281: }
282: }
283: PetscFunctionReturn(PETSC_SUCCESS);
284: }
286: /*@
287: DMCompositeGetLocalAccessArray - Allows one to access the individual
288: packed vectors in their local representation.
290: Collective
292: Input Parameters:
293: + dm - the `DMCOMPOSITE`
294: . pvec - packed vector
295: . nwanted - number of vectors wanted
296: - wanted - sorted array of vectors wanted, or `NULL` to get all vectors, length `nwanted`
298: Output Parameter:
299: . vecs - array of requested local vectors (must be allocated and of length `nwanted`)
301: Level: advanced
303: Note:
304: Use `DMCompositeRestoreLocalAccessArray()` to return the vectors
305: when you no longer need them.
307: .seealso: `DMCOMPOSITE`, `DM`, `DMCompositeRestoreLocalAccessArray()`, `DMCompositeGetAccess()`,
308: `DMCompositeGetEntries()`, `DMCompositeScatter()`, `DMCompositeGather()`
309: @*/
310: PetscErrorCode DMCompositeGetLocalAccessArray(DM dm, Vec pvec, PetscInt nwanted, const PetscInt wanted[], Vec vecs[])
311: {
312: struct DMCompositeLink *link;
313: PetscInt i, wnum;
314: DM_Composite *com = (DM_Composite *)dm->data;
315: PetscInt readonly;
316: PetscInt nlocal = 0;
317: PetscBool flg;
319: PetscFunctionBegin;
322: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
323: PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
324: if (!com->setup) PetscCall(DMSetUp(dm));
326: PetscCall(VecLockGet(pvec, &readonly));
327: for (i = 0, wnum = 0, link = com->next; link && wnum < nwanted; i++, link = link->next) {
328: if (!wanted || i == wanted[wnum]) {
329: Vec v;
330: PetscCall(DMGetLocalVector(link->dm, &v));
331: if (readonly) {
332: const PetscScalar *array;
333: PetscCall(VecGetArrayRead(pvec, &array));
334: PetscCall(VecPlaceArray(v, array + nlocal));
335: // 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
336: PetscCall(VecLockReadPush(v));
337: PetscCall(VecRestoreArrayRead(pvec, &array));
338: } else {
339: PetscScalar *array;
340: PetscCall(VecGetArray(pvec, &array));
341: PetscCall(VecPlaceArray(v, array + nlocal));
342: PetscCall(VecRestoreArray(pvec, &array));
343: }
344: vecs[wnum++] = v;
345: }
347: nlocal += link->nlocal;
348: }
349: PetscFunctionReturn(PETSC_SUCCESS);
350: }
352: /*@C
353: DMCompositeRestoreAccess - Returns the vectors obtained with `DMCompositeGetAccess()`
354: representation.
356: Collective
358: Input Parameters:
359: + dm - the `DMCOMPOSITE` object
360: . gvec - the global vector
361: - ... - the individual parallel vectors, `NULL` for those that are not needed
363: Level: advanced
365: .seealso: `DMCOMPOSITE`, `DM`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`,
366: `DMCompositeGather()`, `DMCompositeCreate()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeScatter()`,
367: `DMCompositeGetAccess()`
368: @*/
369: PetscErrorCode DMCompositeRestoreAccess(DM dm, Vec gvec, ...)
370: {
371: va_list Argp;
372: struct DMCompositeLink *next;
373: DM_Composite *com = (DM_Composite *)dm->data;
374: PetscInt readonly;
375: PetscBool flg;
377: PetscFunctionBegin;
380: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
381: PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
382: next = com->next;
383: if (!com->setup) PetscCall(DMSetUp(dm));
385: PetscCall(VecLockGet(gvec, &readonly));
386: /* loop over packed objects, handling one at a time */
387: va_start(Argp, gvec);
388: while (next) {
389: Vec *vec;
390: vec = va_arg(Argp, Vec *);
391: if (vec) {
392: PetscCall(VecResetArray(*vec));
393: if (readonly) PetscCall(VecLockReadPop(*vec));
394: PetscCall(DMRestoreGlobalVector(next->dm, vec));
395: }
396: next = next->next;
397: }
398: va_end(Argp);
399: PetscFunctionReturn(PETSC_SUCCESS);
400: }
402: /*@
403: DMCompositeRestoreAccessArray - Returns the vectors obtained with `DMCompositeGetAccessArray()`
405: Collective
407: Input Parameters:
408: + dm - the `DMCOMPOSITE` object
409: . pvec - packed vector
410: . nwanted - number of vectors wanted
411: . wanted - sorted array of vectors wanted, or `NULL` to restore all vectors
412: - vecs - array of global vectors
414: Level: advanced
416: .seealso: `DMCOMPOSITE`, `DM`, `DMCompositeRestoreAccess()`, `DMCompositeRestoreEntries()`, `DMCompositeScatter()`, `DMCompositeGather()`
417: @*/
418: PetscErrorCode DMCompositeRestoreAccessArray(DM dm, Vec pvec, PetscInt nwanted, const PetscInt wanted[], Vec vecs[])
419: {
420: struct DMCompositeLink *link;
421: PetscInt i, wnum;
422: DM_Composite *com = (DM_Composite *)dm->data;
423: PetscInt readonly;
424: PetscBool flg;
426: PetscFunctionBegin;
429: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
430: PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
431: if (!com->setup) PetscCall(DMSetUp(dm));
433: PetscCall(VecLockGet(pvec, &readonly));
434: for (i = 0, wnum = 0, link = com->next; link && wnum < nwanted; i++, link = link->next) {
435: if (!wanted || i == wanted[wnum]) {
436: PetscCall(VecResetArray(vecs[wnum]));
437: if (readonly) PetscCall(VecLockReadPop(vecs[wnum]));
438: PetscCall(DMRestoreGlobalVector(link->dm, &vecs[wnum]));
439: wnum++;
440: }
441: }
442: PetscFunctionReturn(PETSC_SUCCESS);
443: }
445: /*@
446: DMCompositeRestoreLocalAccessArray - Returns the vectors obtained with `DMCompositeGetLocalAccessArray()`.
448: Collective
450: Input Parameters:
451: + dm - the `DMCOMPOSITE` object
452: . pvec - packed vector
453: . nwanted - number of vectors wanted
454: . wanted - sorted array of vectors wanted, or `NULL` to restore all vectors
455: - vecs - array of local vectors
457: Level: advanced
459: Note:
460: `nwanted` and `wanted` must match the values given to `DMCompositeGetLocalAccessArray()`
461: otherwise the call will fail.
463: .seealso: `DMCOMPOSITE`, `DM`, `DMCompositeGetLocalAccessArray()`, `DMCompositeRestoreAccessArray()`,
464: `DMCompositeRestoreAccess()`, `DMCompositeRestoreEntries()`,
465: `DMCompositeScatter()`, `DMCompositeGather()`
466: @*/
467: PetscErrorCode DMCompositeRestoreLocalAccessArray(DM dm, Vec pvec, PetscInt nwanted, const PetscInt wanted[], Vec *vecs)
468: {
469: struct DMCompositeLink *link;
470: PetscInt i, wnum;
471: DM_Composite *com = (DM_Composite *)dm->data;
472: PetscInt readonly;
473: PetscBool flg;
475: PetscFunctionBegin;
478: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
479: PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
480: if (!com->setup) PetscCall(DMSetUp(dm));
482: PetscCall(VecLockGet(pvec, &readonly));
483: for (i = 0, wnum = 0, link = com->next; link && wnum < nwanted; i++, link = link->next) {
484: if (!wanted || i == wanted[wnum]) {
485: PetscCall(VecResetArray(vecs[wnum]));
486: if (readonly) PetscCall(VecLockReadPop(vecs[wnum]));
487: PetscCall(DMRestoreLocalVector(link->dm, &vecs[wnum]));
488: wnum++;
489: }
490: }
491: PetscFunctionReturn(PETSC_SUCCESS);
492: }
494: /*@C
495: DMCompositeScatter - Scatters from a global packed vector into its individual local vectors
497: Collective
499: Input Parameters:
500: + dm - the `DMCOMPOSITE` object
501: . gvec - the global vector
502: - ... - the individual sequential vectors, `NULL` for those that are not needed
504: Level: advanced
506: Note:
507: `DMCompositeScatterArray()` is a non-variadic alternative that is often more convenient for library callers and is
508: accessible from Fortran.
510: .seealso: `DMCOMPOSITE`, `DM`, `DMDestroy()`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`,
511: `DMCompositeGather()`, `DMCompositeCreate()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeGetAccess()`,
512: `DMCompositeGetLocalVectors()`, `DMCompositeRestoreLocalVectors()`, `DMCompositeGetEntries()`
513: `DMCompositeScatterArray()`
514: @*/
515: PetscErrorCode DMCompositeScatter(DM dm, Vec gvec, ...)
516: {
517: va_list Argp;
518: struct DMCompositeLink *next;
519: PETSC_UNUSED PetscInt cnt;
520: DM_Composite *com = (DM_Composite *)dm->data;
521: PetscBool flg;
523: PetscFunctionBegin;
526: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
527: PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
528: if (!com->setup) PetscCall(DMSetUp(dm));
530: /* loop over packed objects, handling one at a time */
531: va_start(Argp, gvec);
532: for (cnt = 3, next = com->next; next; cnt++, next = next->next) {
533: Vec local;
534: local = va_arg(Argp, Vec);
535: if (local) {
536: Vec global;
537: const PetscScalar *array;
539: PetscCall(DMGetGlobalVector(next->dm, &global));
540: PetscCall(VecGetArrayRead(gvec, &array));
541: PetscCall(VecPlaceArray(global, array + next->rstart));
542: PetscCall(DMGlobalToLocalBegin(next->dm, global, INSERT_VALUES, local));
543: PetscCall(DMGlobalToLocalEnd(next->dm, global, INSERT_VALUES, local));
544: PetscCall(VecRestoreArrayRead(gvec, &array));
545: PetscCall(VecResetArray(global));
546: PetscCall(DMRestoreGlobalVector(next->dm, &global));
547: }
548: }
549: va_end(Argp);
550: PetscFunctionReturn(PETSC_SUCCESS);
551: }
553: /*@
554: DMCompositeScatterArray - Scatters from a global packed vector into its individual local vectors
556: Collective
558: Input Parameters:
559: + dm - the `DMCOMPOSITE` object
560: . gvec - the global vector
561: - lvecs - array of local vectors, NULL for any that are not needed
563: Level: advanced
565: Note:
566: This is a non-variadic alternative to `DMCompositeScatter()`
568: .seealso: `DMCOMPOSITE`, `DM`, `DMDestroy()`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`
569: `DMCompositeGather()`, `DMCompositeCreate()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeGetAccess()`,
570: `DMCompositeGetLocalVectors()`, `DMCompositeRestoreLocalVectors()`, `DMCompositeGetEntries()`
571: @*/
572: PetscErrorCode DMCompositeScatterArray(DM dm, Vec gvec, Vec *lvecs)
573: {
574: struct DMCompositeLink *next;
575: PetscInt i;
576: DM_Composite *com = (DM_Composite *)dm->data;
577: PetscBool flg;
579: PetscFunctionBegin;
582: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
583: PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
584: if (!com->setup) PetscCall(DMSetUp(dm));
586: /* loop over packed objects, handling one at a time */
587: for (i = 0, next = com->next; next; next = next->next, i++) {
588: if (lvecs[i]) {
589: Vec global;
590: const PetscScalar *array;
592: PetscCall(DMGetGlobalVector(next->dm, &global));
593: PetscCall(VecGetArrayRead(gvec, &array));
594: PetscCall(VecPlaceArray(global, (PetscScalar *)array + next->rstart));
595: PetscCall(DMGlobalToLocalBegin(next->dm, global, INSERT_VALUES, lvecs[i]));
596: PetscCall(DMGlobalToLocalEnd(next->dm, global, INSERT_VALUES, lvecs[i]));
597: PetscCall(VecRestoreArrayRead(gvec, &array));
598: PetscCall(VecResetArray(global));
599: PetscCall(DMRestoreGlobalVector(next->dm, &global));
600: }
601: }
602: PetscFunctionReturn(PETSC_SUCCESS);
603: }
605: /*@C
606: DMCompositeGather - Gathers into a global packed vector from its individual local vectors
608: Collective
610: Input Parameters:
611: + dm - the `DMCOMPOSITE` object
612: . imode - `INSERT_VALUES` or `ADD_VALUES`
613: . gvec - the global vector
614: - ... - the individual sequential vectors, `NULL` for any that are not needed
616: Level: advanced
618: Fortran Notes:
619: Fortran users should use `DMCompositeGatherArray()`
621: .seealso: `DMCOMPOSITE`, `DM`, `DMDestroy()`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`,
622: `DMCompositeScatter()`, `DMCompositeCreate()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeGetAccess()`,
623: `DMCompositeGetLocalVectors()`, `DMCompositeRestoreLocalVectors()`, `DMCompositeGetEntries()`
624: @*/
625: PetscErrorCode DMCompositeGather(DM dm, InsertMode imode, Vec gvec, ...)
626: {
627: va_list Argp;
628: struct DMCompositeLink *next;
629: DM_Composite *com = (DM_Composite *)dm->data;
630: PETSC_UNUSED PetscInt cnt;
631: PetscBool flg;
633: PetscFunctionBegin;
636: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
637: PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
638: if (!com->setup) PetscCall(DMSetUp(dm));
640: /* loop over packed objects, handling one at a time */
641: va_start(Argp, gvec);
642: for (cnt = 3, next = com->next; next; cnt++, next = next->next) {
643: Vec local;
644: local = va_arg(Argp, Vec);
645: if (local) {
646: PetscScalar *array;
647: Vec global;
649: PetscCall(DMGetGlobalVector(next->dm, &global));
650: PetscCall(VecGetArray(gvec, &array));
651: PetscCall(VecPlaceArray(global, array + next->rstart));
652: PetscCall(DMLocalToGlobalBegin(next->dm, local, imode, global));
653: PetscCall(DMLocalToGlobalEnd(next->dm, local, imode, global));
654: PetscCall(VecRestoreArray(gvec, &array));
655: PetscCall(VecResetArray(global));
656: PetscCall(DMRestoreGlobalVector(next->dm, &global));
657: }
658: }
659: va_end(Argp);
660: PetscFunctionReturn(PETSC_SUCCESS);
661: }
663: /*@
664: DMCompositeGatherArray - Gathers into a global packed vector from its individual local vectors
666: Collective
668: Input Parameters:
669: + dm - the `DMCOMPOSITE` object
670: . gvec - the global vector
671: . imode - `INSERT_VALUES` or `ADD_VALUES`
672: - lvecs - the individual sequential vectors, NULL for any that are not needed
674: Level: advanced
676: Note:
677: This is a non-variadic alternative to `DMCompositeGather()`.
679: .seealso: `DMCOMPOSITE`, `DM`, `DMDestroy()`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`,
680: `DMCompositeScatter()`, `DMCompositeCreate()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeGetAccess()`,
681: `DMCompositeGetLocalVectors()`, `DMCompositeRestoreLocalVectors()`, `DMCompositeGetEntries()`,
682: @*/
683: PetscErrorCode DMCompositeGatherArray(DM dm, InsertMode imode, Vec gvec, Vec *lvecs)
684: {
685: struct DMCompositeLink *next;
686: DM_Composite *com = (DM_Composite *)dm->data;
687: PetscInt i;
688: PetscBool flg;
690: PetscFunctionBegin;
693: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
694: PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
695: if (!com->setup) PetscCall(DMSetUp(dm));
697: /* loop over packed objects, handling one at a time */
698: for (next = com->next, i = 0; next; next = next->next, i++) {
699: if (lvecs[i]) {
700: PetscScalar *array;
701: Vec global;
703: PetscCall(DMGetGlobalVector(next->dm, &global));
704: PetscCall(VecGetArray(gvec, &array));
705: PetscCall(VecPlaceArray(global, array + next->rstart));
706: PetscCall(DMLocalToGlobalBegin(next->dm, lvecs[i], imode, global));
707: PetscCall(DMLocalToGlobalEnd(next->dm, lvecs[i], imode, global));
708: PetscCall(VecRestoreArray(gvec, &array));
709: PetscCall(VecResetArray(global));
710: PetscCall(DMRestoreGlobalVector(next->dm, &global));
711: }
712: }
713: PetscFunctionReturn(PETSC_SUCCESS);
714: }
716: /*@
717: DMCompositeAddDM - adds a `DM` vector to a `DMCOMPOSITE`
719: Collective
721: Input Parameters:
722: + dmc - the `DMCOMPOSITE` object
723: - dm - the `DM` object
725: Level: advanced
727: .seealso: `DMCOMPOSITE`, `DM`, `DMDestroy()`, `DMCompositeGather()`, `DMCreateGlobalVector()`,
728: `DMCompositeScatter()`, `DMCompositeCreate()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeGetAccess()`,
729: `DMCompositeGetLocalVectors()`, `DMCompositeRestoreLocalVectors()`, `DMCompositeGetEntries()`
730: @*/
731: PetscErrorCode DMCompositeAddDM(DM dmc, DM dm)
732: {
733: PetscInt n, nlocal;
734: struct DMCompositeLink *mine, *next;
735: Vec global, local;
736: DM_Composite *com = (DM_Composite *)dmc->data;
737: PetscBool flg;
739: PetscFunctionBegin;
742: PetscCall(PetscObjectTypeCompare((PetscObject)dmc, DMCOMPOSITE, &flg));
743: PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
744: next = com->next;
745: PetscCheck(!com->setup, PetscObjectComm((PetscObject)dmc), PETSC_ERR_ARG_WRONGSTATE, "Cannot add a DM once you have used the DMComposite");
747: /* create new link */
748: PetscCall(PetscNew(&mine));
749: PetscCall(PetscObjectReference((PetscObject)dm));
750: PetscCall(DMCreateGlobalVector(dm, &global)); // Not using DMGetGlobalVector() since it will stash the vector with a type decided by dm,
751: PetscCall(VecGetLocalSize(global, &n)); // while we may want to set dmc's vectype later, say via DMSetFromOptions(dmc), and we
752: PetscCall(VecDestroy(&global)); // want to propagate the type to dm.
753: PetscCall(DMCreateLocalVector(dm, &local)); // Not using DMGetLocalVector(), same reason as above.
754: PetscCall(VecGetSize(local, &nlocal));
755: PetscCall(VecDestroy(&local));
757: mine->n = n;
758: mine->nlocal = nlocal;
759: mine->dm = dm;
760: mine->next = NULL;
761: com->n += n;
762: com->nghost += nlocal;
764: /* add to end of list */
765: if (!next) com->next = mine;
766: else {
767: while (next->next) next = next->next;
768: next->next = mine;
769: }
770: com->nDM++;
771: com->nmine++;
772: PetscFunctionReturn(PETSC_SUCCESS);
773: }
775: #include <petscdraw.h>
776: PETSC_SINGLE_LIBRARY_INTERN PetscErrorCode VecView_MPI(Vec, PetscViewer);
778: static PetscErrorCode VecView_DMComposite(Vec gvec, PetscViewer viewer)
779: {
780: DM dm;
781: struct DMCompositeLink *next;
782: PetscBool isdraw;
783: DM_Composite *com;
785: PetscFunctionBegin;
786: PetscCall(VecGetDM(gvec, &dm));
787: PetscCheck(dm, PetscObjectComm((PetscObject)gvec), PETSC_ERR_ARG_WRONG, "Vector not generated from a DMComposite");
788: com = (DM_Composite *)dm->data;
789: next = com->next;
791: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw));
792: if (!isdraw) {
793: /* do I really want to call this? */
794: PetscCall(VecView_MPI(gvec, viewer));
795: } else {
796: PetscInt cnt = 0;
798: /* loop over packed objects, handling one at a time */
799: while (next) {
800: Vec vec;
801: const PetscScalar *array;
802: PetscInt bs;
804: /* Should use VecGetSubVector() eventually, but would need to forward the DM for that to work */
805: PetscCall(DMGetGlobalVector(next->dm, &vec));
806: PetscCall(VecGetArrayRead(gvec, &array));
807: PetscCall(VecPlaceArray(vec, (PetscScalar *)array + next->rstart));
808: PetscCall(VecRestoreArrayRead(gvec, &array));
809: PetscCall(VecView(vec, viewer));
810: PetscCall(VecResetArray(vec));
811: PetscCall(VecGetBlockSize(vec, &bs));
812: PetscCall(DMRestoreGlobalVector(next->dm, &vec));
813: PetscCall(PetscViewerDrawBaseAdd(viewer, bs));
814: cnt += bs;
815: next = next->next;
816: }
817: PetscCall(PetscViewerDrawBaseAdd(viewer, -cnt));
818: }
819: PetscFunctionReturn(PETSC_SUCCESS);
820: }
822: static PetscErrorCode DMCreateGlobalVector_Composite(DM dm, Vec *gvec)
823: {
824: DM_Composite *com = (DM_Composite *)dm->data;
826: PetscFunctionBegin;
828: PetscCall(DMSetFromOptions(dm));
829: PetscCall(DMSetUp(dm));
830: PetscCall(VecCreate(PetscObjectComm((PetscObject)dm), gvec));
831: PetscCall(VecSetType(*gvec, dm->vectype));
832: PetscCall(VecSetSizes(*gvec, com->n, com->N));
833: PetscCall(VecSetDM(*gvec, dm));
834: PetscCall(VecSetOperation(*gvec, VECOP_VIEW, (PetscErrorCodeFn *)VecView_DMComposite));
835: PetscFunctionReturn(PETSC_SUCCESS);
836: }
838: static PetscErrorCode DMCreateLocalVector_Composite(DM dm, Vec *lvec)
839: {
840: DM_Composite *com = (DM_Composite *)dm->data;
842: PetscFunctionBegin;
844: if (!com->setup) {
845: PetscCall(DMSetFromOptions(dm));
846: PetscCall(DMSetUp(dm));
847: }
848: PetscCall(VecCreate(PETSC_COMM_SELF, lvec));
849: PetscCall(VecSetType(*lvec, dm->vectype));
850: PetscCall(VecSetSizes(*lvec, com->nghost, PETSC_DECIDE));
851: PetscCall(VecSetDM(*lvec, dm));
852: PetscFunctionReturn(PETSC_SUCCESS);
853: }
855: /*@C
856: DMCompositeGetISLocalToGlobalMappings - gets an `ISLocalToGlobalMapping` for each `DM` in the `DMCOMPOSITE`, maps to the composite global space
858: Collective; No Fortran Support
860: Input Parameter:
861: . dm - the `DMCOMPOSITE` object
863: Output Parameter:
864: . ltogs - the individual mappings for each packed vector. Note that this includes
865: all the ghost points that individual ghosted `DMDA` may have.
867: Level: advanced
869: Note:
870: Each entry of `ltogs` should be destroyed with `ISLocalToGlobalMappingDestroy()`, `ltogs` should be freed with `PetscFree()`.
872: .seealso: `DMCOMPOSITE`, `DM`, `DMDestroy()`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`,
873: `DMCompositeGather()`, `DMCompositeCreate()`, `DMCompositeGetAccess()`, `DMCompositeScatter()`,
874: `DMCompositeGetLocalVectors()`, `DMCompositeRestoreLocalVectors()`, `DMCompositeGetEntries()`
875: @*/
876: PetscErrorCode DMCompositeGetISLocalToGlobalMappings(DM dm, ISLocalToGlobalMapping *ltogs[])
877: {
878: PetscInt i, *idx, n, cnt;
879: struct DMCompositeLink *next;
880: PetscMPIInt rank;
881: DM_Composite *com = (DM_Composite *)dm->data;
882: PetscBool flg;
884: PetscFunctionBegin;
886: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
887: PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
888: PetscCall(DMSetUp(dm));
889: PetscCall(PetscMalloc1(com->nDM, ltogs));
890: next = com->next;
891: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
893: /* loop over packed objects, handling one at a time */
894: cnt = 0;
895: while (next) {
896: ISLocalToGlobalMapping ltog;
897: PetscMPIInt size;
898: const PetscInt *suboff, *indices;
899: Vec global;
901: /* Get sub-DM global indices for each local dof */
902: PetscCall(DMGetLocalToGlobalMapping(next->dm, <og));
903: PetscCall(ISLocalToGlobalMappingGetSize(ltog, &n));
904: PetscCall(ISLocalToGlobalMappingGetIndices(ltog, &indices));
905: PetscCall(PetscMalloc1(n, &idx));
907: /* Get the offsets for the sub-DM global vector */
908: PetscCall(DMGetGlobalVector(next->dm, &global));
909: PetscCall(VecGetOwnershipRanges(global, &suboff));
910: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)global), &size));
912: /* Shift the sub-DM definition of the global space to the composite global space */
913: for (i = 0; i < n; i++) {
914: PetscInt subi = indices[i], lo = 0, hi = size, t;
915: /* There's no consensus on what a negative index means,
916: except for skipping when setting the values in vectors and matrices */
917: if (subi < 0) {
918: idx[i] = subi - next->grstarts[rank];
919: continue;
920: }
921: /* Binary search to find which rank owns subi */
922: while (hi - lo > 1) {
923: t = lo + (hi - lo) / 2;
924: if (suboff[t] > subi) hi = t;
925: else lo = t;
926: }
927: idx[i] = subi - suboff[lo] + next->grstarts[lo];
928: }
929: PetscCall(ISLocalToGlobalMappingRestoreIndices(ltog, &indices));
930: PetscCall(ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)dm), 1, n, idx, PETSC_OWN_POINTER, &(*ltogs)[cnt]));
931: PetscCall(DMRestoreGlobalVector(next->dm, &global));
932: next = next->next;
933: cnt++;
934: }
935: PetscFunctionReturn(PETSC_SUCCESS);
936: }
938: /*@C
939: DMCompositeGetLocalISs - Gets index sets for each component of a composite local vector
941: Not Collective; No Fortran Support
943: Input Parameter:
944: . dm - the `DMCOMPOSITE`
946: Output Parameter:
947: . is - array of serial index sets for each component of the `DMCOMPOSITE`
949: Level: intermediate
951: Notes:
952: At present, a composite local vector does not normally exist. This function is used to provide index sets for
953: `MatGetLocalSubMatrix()`. In the future, the scatters for each entry in the `DMCOMPOSITE` may be merged into a single
954: scatter to a composite local vector. The user should not typically need to know which is being done.
956: To get the composite global indices at all local points (including ghosts), use `DMCompositeGetISLocalToGlobalMappings()`.
958: To get index sets for pieces of the composite global vector, use `DMCompositeGetGlobalISs()`.
960: Each returned `IS` should be destroyed with `ISDestroy()`, the array should be freed with `PetscFree()`.
962: Fortran Note:
963: Use `DMCompositeRestoreLocalISs()` to release the `is`.
965: .seealso: `DMCOMPOSITE`, `DM`, `DMCompositeGetGlobalISs()`, `DMCompositeGetISLocalToGlobalMappings()`, `MatGetLocalSubMatrix()`,
966: `MatCreateLocalRef()`, `DMCompositeGetNumberDM()`
967: @*/
968: PetscErrorCode DMCompositeGetLocalISs(DM dm, IS *is[])
969: {
970: DM_Composite *com = (DM_Composite *)dm->data;
971: struct DMCompositeLink *link;
972: PetscInt cnt, start;
973: PetscBool flg;
975: PetscFunctionBegin;
977: PetscAssertPointer(is, 2);
978: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
979: PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
980: PetscCall(PetscMalloc1(com->nmine, is));
981: for (cnt = 0, start = 0, link = com->next; link; start += link->nlocal, cnt++, link = link->next) {
982: PetscInt bs;
983: PetscCall(ISCreateStride(PETSC_COMM_SELF, link->nlocal, start, 1, &(*is)[cnt]));
984: PetscCall(DMGetBlockSize(link->dm, &bs));
985: PetscCall(ISSetBlockSize((*is)[cnt], bs));
986: }
987: PetscFunctionReturn(PETSC_SUCCESS);
988: }
990: /*@C
991: DMCompositeGetGlobalISs - Gets the index sets for each composed object in a `DMCOMPOSITE`
993: Collective
995: Input Parameter:
996: . dm - the `DMCOMPOSITE` object
998: Output Parameter:
999: . is - the array of index sets
1001: Level: advanced
1003: Notes:
1004: The `is` entries should be destroyed with `ISDestroy()`, `is` should be freed with `PetscFree()`
1006: These could be used to extract a subset of vector entries for a "multi-physics" preconditioner
1008: Use `DMCompositeGetLocalISs()` for index sets in the packed local numbering, and
1009: `DMCompositeGetISLocalToGlobalMappings()` for to map local sub-`DM` (including ghost) indices to packed global
1010: indices.
1012: Fortran Note:
1013: Use `DMCompositeRestoreGlobalISs()` to release the `is`.
1015: .seealso: `DMCOMPOSITE`, `DM`, `DMDestroy()`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`,
1016: `DMCompositeGather()`, `DMCompositeCreate()`, `DMCompositeGetAccess()`, `DMCompositeScatter()`,
1017: `DMCompositeGetLocalVectors()`, `DMCompositeRestoreLocalVectors()`, `DMCompositeGetEntries()`
1018: @*/
1019: PetscErrorCode DMCompositeGetGlobalISs(DM dm, IS *is[])
1020: {
1021: PetscInt cnt = 0;
1022: struct DMCompositeLink *next;
1023: PetscMPIInt rank;
1024: DM_Composite *com = (DM_Composite *)dm->data;
1025: PetscBool flg;
1027: PetscFunctionBegin;
1029: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
1030: PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
1031: PetscCheck(dm->setupcalled, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "Must call DMSetUp() before");
1032: PetscCall(PetscMalloc1(com->nDM, is));
1033: next = com->next;
1034: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
1036: /* loop over packed objects, handling one at a time */
1037: while (next) {
1038: PetscDS prob;
1040: PetscCall(ISCreateStride(PetscObjectComm((PetscObject)dm), next->n, next->grstart, 1, &(*is)[cnt]));
1041: PetscCall(DMGetDS(dm, &prob));
1042: if (prob) {
1043: MatNullSpace space;
1044: Mat pmat;
1045: PetscObject disc;
1046: PetscInt Nf;
1048: PetscCall(PetscDSGetNumFields(prob, &Nf));
1049: if (cnt < Nf) {
1050: PetscCall(PetscDSGetDiscretization(prob, cnt, &disc));
1051: PetscCall(PetscObjectQuery(disc, "nullspace", (PetscObject *)&space));
1052: if (space) PetscCall(PetscObjectCompose((PetscObject)(*is)[cnt], "nullspace", (PetscObject)space));
1053: PetscCall(PetscObjectQuery(disc, "nearnullspace", (PetscObject *)&space));
1054: if (space) PetscCall(PetscObjectCompose((PetscObject)(*is)[cnt], "nearnullspace", (PetscObject)space));
1055: PetscCall(PetscObjectQuery(disc, "pmat", (PetscObject *)&pmat));
1056: if (pmat) PetscCall(PetscObjectCompose((PetscObject)(*is)[cnt], "pmat", (PetscObject)pmat));
1057: }
1058: }
1059: cnt++;
1060: next = next->next;
1061: }
1062: PetscFunctionReturn(PETSC_SUCCESS);
1063: }
1065: static PetscErrorCode DMCreateFieldIS_Composite(DM dm, PetscInt *numFields, char ***fieldNames, IS **fields)
1066: {
1067: PetscInt nDM;
1068: DM *dms;
1069: PetscInt i;
1071: PetscFunctionBegin;
1072: PetscCall(DMCompositeGetNumberDM(dm, &nDM));
1073: if (numFields) *numFields = nDM;
1074: PetscCall(DMCompositeGetGlobalISs(dm, fields));
1075: if (fieldNames) {
1076: PetscCall(PetscMalloc1(nDM, &dms));
1077: PetscCall(PetscMalloc1(nDM, fieldNames));
1078: PetscCall(DMCompositeGetEntriesArray(dm, dms));
1079: for (i = 0; i < nDM; i++) {
1080: char buf[256];
1081: const char *splitname;
1083: /* Split naming precedence: object name, prefix, number */
1084: splitname = ((PetscObject)dm)->name;
1085: if (!splitname) {
1086: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)dms[i], &splitname));
1087: if (splitname) {
1088: size_t len;
1089: PetscCall(PetscStrncpy(buf, splitname, sizeof(buf)));
1090: buf[sizeof(buf) - 1] = 0;
1091: PetscCall(PetscStrlen(buf, &len));
1092: if (buf[len - 1] == '_') buf[len - 1] = 0; /* Remove trailing underscore if it was used */
1093: splitname = buf;
1094: }
1095: }
1096: if (!splitname) {
1097: PetscCall(PetscSNPrintf(buf, sizeof(buf), "%" PetscInt_FMT, i));
1098: splitname = buf;
1099: }
1100: PetscCall(PetscStrallocpy(splitname, &(*fieldNames)[i]));
1101: }
1102: PetscCall(PetscFree(dms));
1103: }
1104: PetscFunctionReturn(PETSC_SUCCESS);
1105: }
1107: /*
1108: This could take over from DMCreateFieldIS(), as it is more general,
1109: making DMCreateFieldIS() a special case -- calling with dmlist == NULL;
1110: At this point it's probably best to be less intrusive, however.
1111: */
1112: static PetscErrorCode DMCreateFieldDecomposition_Composite(DM dm, PetscInt *len, char ***namelist, IS **islist, DM **dmlist)
1113: {
1114: PetscInt nDM;
1115: PetscInt i;
1117: PetscFunctionBegin;
1118: PetscCall(DMCreateFieldIS_Composite(dm, len, namelist, islist));
1119: if (dmlist) {
1120: PetscCall(DMCompositeGetNumberDM(dm, &nDM));
1121: PetscCall(PetscMalloc1(nDM, dmlist));
1122: PetscCall(DMCompositeGetEntriesArray(dm, *dmlist));
1123: for (i = 0; i < nDM; i++) PetscCall(PetscObjectReference((PetscObject)((*dmlist)[i])));
1124: }
1125: PetscFunctionReturn(PETSC_SUCCESS);
1126: }
1128: /*@C
1129: DMCompositeGetLocalVectors - Gets local vectors for each part of a `DMCOMPOSITE`
1130: Use `DMCompositeRestoreLocalVectors()` to return them.
1132: Not Collective; No Fortran Support
1134: Input Parameter:
1135: . dm - the `DMCOMPOSITE` object
1137: Output Parameter:
1138: . ... - the individual sequential `Vec`s
1140: Level: advanced
1142: .seealso: `DMCOMPOSITE`, `DM`, `DMDestroy()`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`,
1143: `DMCompositeGather()`, `DMCompositeCreate()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeGetAccess()`,
1144: `DMCompositeRestoreLocalVectors()`, `DMCompositeScatter()`, `DMCompositeGetEntries()`
1145: @*/
1146: PetscErrorCode DMCompositeGetLocalVectors(DM dm, ...)
1147: {
1148: va_list Argp;
1149: struct DMCompositeLink *next;
1150: DM_Composite *com = (DM_Composite *)dm->data;
1151: PetscBool flg;
1153: PetscFunctionBegin;
1155: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
1156: PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
1157: next = com->next;
1158: /* loop over packed objects, handling one at a time */
1159: va_start(Argp, dm);
1160: while (next) {
1161: Vec *vec;
1162: vec = va_arg(Argp, Vec *);
1163: if (vec) PetscCall(DMGetLocalVector(next->dm, vec));
1164: next = next->next;
1165: }
1166: va_end(Argp);
1167: PetscFunctionReturn(PETSC_SUCCESS);
1168: }
1170: /*@C
1171: DMCompositeRestoreLocalVectors - Restores local vectors for each part of a `DMCOMPOSITE`
1173: Not Collective; No Fortran Support
1175: Input Parameter:
1176: . dm - the `DMCOMPOSITE` object
1178: Output Parameter:
1179: . ... - the individual sequential `Vec`
1181: Level: advanced
1183: .seealso: `DMCOMPOSITE`, `DM`, `DMDestroy()`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`,
1184: `DMCompositeGather()`, `DMCompositeCreate()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeGetAccess()`,
1185: `DMCompositeGetLocalVectors()`, `DMCompositeScatter()`, `DMCompositeGetEntries()`
1186: @*/
1187: PetscErrorCode DMCompositeRestoreLocalVectors(DM dm, ...)
1188: {
1189: va_list Argp;
1190: struct DMCompositeLink *next;
1191: DM_Composite *com = (DM_Composite *)dm->data;
1192: PetscBool flg;
1194: PetscFunctionBegin;
1196: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
1197: PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
1198: next = com->next;
1199: /* loop over packed objects, handling one at a time */
1200: va_start(Argp, dm);
1201: while (next) {
1202: Vec *vec;
1203: vec = va_arg(Argp, Vec *);
1204: if (vec) PetscCall(DMRestoreLocalVector(next->dm, vec));
1205: next = next->next;
1206: }
1207: va_end(Argp);
1208: PetscFunctionReturn(PETSC_SUCCESS);
1209: }
1211: /*@C
1212: DMCompositeGetEntries - Gets the `DM` for each entry in a `DMCOMPOSITE`.
1214: Not Collective
1216: Input Parameter:
1217: . dm - the `DMCOMPOSITE` object
1219: Output Parameter:
1220: . ... - the individual entries `DM`
1222: Level: advanced
1224: Fortran Notes:
1225: Use `DMCompositeGetEntriesArray()`
1227: .seealso: `DMCOMPOSITE`, `DM`, `DMDestroy()`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`, `DMCompositeGetEntriesArray()`
1228: `DMCompositeGather()`, `DMCompositeCreate()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeGetAccess()`,
1229: `DMCompositeRestoreLocalVectors()`, `DMCompositeGetLocalVectors()`, `DMCompositeScatter()`
1230: @*/
1231: PetscErrorCode DMCompositeGetEntries(DM dm, ...)
1232: {
1233: va_list Argp;
1234: struct DMCompositeLink *next;
1235: DM_Composite *com = (DM_Composite *)dm->data;
1236: PetscBool flg;
1238: PetscFunctionBegin;
1240: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
1241: PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
1242: next = com->next;
1243: /* loop over packed objects, handling one at a time */
1244: va_start(Argp, dm);
1245: while (next) {
1246: DM *dmn;
1247: dmn = va_arg(Argp, DM *);
1248: if (dmn) *dmn = next->dm;
1249: next = next->next;
1250: }
1251: va_end(Argp);
1252: PetscFunctionReturn(PETSC_SUCCESS);
1253: }
1255: /*@
1256: DMCompositeGetEntriesArray - Gets the DM for each entry in a `DMCOMPOSITE`
1258: Not Collective
1260: Input Parameter:
1261: . dm - the `DMCOMPOSITE` object
1263: Output Parameter:
1264: . dms - array of sufficient length (see `DMCompositeGetNumberDM()`) to hold the individual `DM`
1266: Level: advanced
1268: .seealso: `DMCOMPOSITE`, `DM`, `DMDestroy()`, `DMCompositeAddDM()`, `DMCreateGlobalVector()`, `DMCompositeGetEntries()`
1269: `DMCompositeGather()`, `DMCompositeCreate()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeGetAccess()`,
1270: `DMCompositeRestoreLocalVectors()`, `DMCompositeGetLocalVectors()`, `DMCompositeScatter()`
1271: @*/
1272: PetscErrorCode DMCompositeGetEntriesArray(DM dm, DM dms[])
1273: {
1274: struct DMCompositeLink *next;
1275: DM_Composite *com = (DM_Composite *)dm->data;
1276: PetscInt i;
1277: PetscBool flg;
1279: PetscFunctionBegin;
1281: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMCOMPOSITE, &flg));
1282: PetscCheck(flg, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for type %s", ((PetscObject)dm)->type_name);
1283: /* loop over packed objects, handling one at a time */
1284: for (next = com->next, i = 0; next; next = next->next, i++) dms[i] = next->dm;
1285: PetscFunctionReturn(PETSC_SUCCESS);
1286: }
1288: typedef struct {
1289: DM dm;
1290: PetscViewer *subv;
1291: Vec *vecs;
1292: } GLVisViewerCtx;
1294: static PetscErrorCode DestroyGLVisViewerCtx_Private(void **vctx)
1295: {
1296: GLVisViewerCtx *ctx = *(GLVisViewerCtx **)vctx;
1297: PetscInt i, n;
1299: PetscFunctionBegin;
1300: PetscCall(DMCompositeGetNumberDM(ctx->dm, &n));
1301: for (i = 0; i < n; i++) PetscCall(PetscViewerDestroy(&ctx->subv[i]));
1302: PetscCall(PetscFree2(ctx->subv, ctx->vecs));
1303: PetscCall(DMDestroy(&ctx->dm));
1304: PetscCall(PetscFree(ctx));
1305: PetscFunctionReturn(PETSC_SUCCESS);
1306: }
1308: static PetscErrorCode DMCompositeSampleGLVisFields_Private(PetscObject oX, PetscInt nf, PetscObject oXfield[], void *vctx)
1309: {
1310: Vec X = (Vec)oX;
1311: GLVisViewerCtx *ctx = (GLVisViewerCtx *)vctx;
1312: PetscInt i, n, cumf;
1314: PetscFunctionBegin;
1315: PetscCall(DMCompositeGetNumberDM(ctx->dm, &n));
1316: PetscCall(DMCompositeGetAccessArray(ctx->dm, X, n, NULL, ctx->vecs));
1317: for (i = 0, cumf = 0; i < n; i++) {
1318: PetscErrorCode (*g2l)(PetscObject, PetscInt, PetscObject[], void *);
1319: void *fctx;
1320: PetscInt nfi;
1322: PetscCall(PetscViewerGLVisGetFields_Internal(ctx->subv[i], &nfi, NULL, NULL, &g2l, NULL, &fctx));
1323: if (!nfi) continue;
1324: if (g2l) PetscCall((*g2l)((PetscObject)ctx->vecs[i], nfi, oXfield + cumf, fctx));
1325: else PetscCall(VecCopy(ctx->vecs[i], (Vec)oXfield[cumf]));
1326: cumf += nfi;
1327: }
1328: PetscCall(DMCompositeRestoreAccessArray(ctx->dm, X, n, NULL, ctx->vecs));
1329: PetscFunctionReturn(PETSC_SUCCESS);
1330: }
1332: static PetscErrorCode DMSetUpGLVisViewer_Composite(PetscObject odm, PetscViewer viewer)
1333: {
1334: DM dm = (DM)odm, *dms;
1335: Vec *Ufds;
1336: GLVisViewerCtx *ctx;
1337: PetscInt i, n, tnf, *sdim;
1338: char **fecs;
1340: PetscFunctionBegin;
1341: PetscCall(PetscNew(&ctx));
1342: PetscCall(PetscObjectReference((PetscObject)dm));
1343: ctx->dm = dm;
1344: PetscCall(DMCompositeGetNumberDM(dm, &n));
1345: PetscCall(PetscMalloc1(n, &dms));
1346: PetscCall(DMCompositeGetEntriesArray(dm, dms));
1347: PetscCall(PetscMalloc2(n, &ctx->subv, n, &ctx->vecs));
1348: for (i = 0, tnf = 0; i < n; i++) {
1349: PetscInt nf;
1351: PetscCall(PetscViewerCreate(PetscObjectComm(odm), &ctx->subv[i]));
1352: PetscCall(PetscViewerSetType(ctx->subv[i], PETSCVIEWERGLVIS));
1353: PetscCall(PetscViewerGLVisSetDM_Internal(ctx->subv[i], (PetscObject)dms[i]));
1354: PetscCall(PetscViewerGLVisGetFields_Internal(ctx->subv[i], &nf, NULL, NULL, NULL, NULL, NULL));
1355: tnf += nf;
1356: }
1357: PetscCall(PetscFree(dms));
1358: PetscCall(PetscMalloc3(tnf, &fecs, tnf, &sdim, tnf, &Ufds));
1359: for (i = 0, tnf = 0; i < n; i++) {
1360: PetscInt *sd, nf, f;
1361: const char **fec;
1362: Vec *Uf;
1364: PetscCall(PetscViewerGLVisGetFields_Internal(ctx->subv[i], &nf, &fec, &sd, NULL, (PetscObject **)&Uf, NULL));
1365: for (f = 0; f < nf; f++) {
1366: PetscCall(PetscStrallocpy(fec[f], &fecs[tnf + f]));
1367: Ufds[tnf + f] = Uf[f];
1368: sdim[tnf + f] = sd[f];
1369: }
1370: tnf += nf;
1371: }
1372: PetscCall(PetscViewerGLVisSetFields(viewer, tnf, (const char **)fecs, sdim, DMCompositeSampleGLVisFields_Private, (PetscObject *)Ufds, ctx, DestroyGLVisViewerCtx_Private));
1373: for (i = 0; i < tnf; i++) PetscCall(PetscFree(fecs[i]));
1374: PetscCall(PetscFree3(fecs, sdim, Ufds));
1375: PetscFunctionReturn(PETSC_SUCCESS);
1376: }
1378: static PetscErrorCode DMRefine_Composite(DM dmi, MPI_Comm comm, DM *fine)
1379: {
1380: struct DMCompositeLink *next;
1381: DM_Composite *com = (DM_Composite *)dmi->data;
1382: DM dm;
1384: PetscFunctionBegin;
1386: if (comm == MPI_COMM_NULL) PetscCall(PetscObjectGetComm((PetscObject)dmi, &comm));
1387: PetscCall(DMSetUp(dmi));
1388: next = com->next;
1389: PetscCall(DMCompositeCreate(comm, fine));
1391: /* loop over packed objects, handling one at a time */
1392: while (next) {
1393: PetscCall(DMRefine(next->dm, comm, &dm));
1394: PetscCall(DMCompositeAddDM(*fine, dm));
1395: PetscCall(PetscObjectDereference((PetscObject)dm));
1396: next = next->next;
1397: }
1398: PetscFunctionReturn(PETSC_SUCCESS);
1399: }
1401: static PetscErrorCode DMCoarsen_Composite(DM dmi, MPI_Comm comm, DM *fine)
1402: {
1403: struct DMCompositeLink *next;
1404: DM_Composite *com = (DM_Composite *)dmi->data;
1405: DM dm;
1407: PetscFunctionBegin;
1409: PetscCall(DMSetUp(dmi));
1410: if (comm == MPI_COMM_NULL) PetscCall(PetscObjectGetComm((PetscObject)dmi, &comm));
1411: next = com->next;
1412: PetscCall(DMCompositeCreate(comm, fine));
1414: /* loop over packed objects, handling one at a time */
1415: while (next) {
1416: PetscCall(DMCoarsen(next->dm, comm, &dm));
1417: PetscCall(DMCompositeAddDM(*fine, dm));
1418: PetscCall(PetscObjectDereference((PetscObject)dm));
1419: next = next->next;
1420: }
1421: PetscFunctionReturn(PETSC_SUCCESS);
1422: }
1424: static PetscErrorCode DMCreateInterpolation_Composite(DM coarse, DM fine, Mat *A, Vec *v)
1425: {
1426: PetscInt m, n, M, N, nDM, i;
1427: struct DMCompositeLink *nextc;
1428: struct DMCompositeLink *nextf;
1429: Vec gcoarse, gfine, *vecs;
1430: DM_Composite *comcoarse = (DM_Composite *)coarse->data;
1431: DM_Composite *comfine = (DM_Composite *)fine->data;
1432: Mat *mats;
1434: PetscFunctionBegin;
1437: PetscCall(DMSetUp(coarse));
1438: PetscCall(DMSetUp(fine));
1439: /* use global vectors only for determining matrix layout */
1440: PetscCall(DMGetGlobalVector(coarse, &gcoarse));
1441: PetscCall(DMGetGlobalVector(fine, &gfine));
1442: PetscCall(VecGetLocalSize(gcoarse, &n));
1443: PetscCall(VecGetLocalSize(gfine, &m));
1444: PetscCall(VecGetSize(gcoarse, &N));
1445: PetscCall(VecGetSize(gfine, &M));
1446: PetscCall(DMRestoreGlobalVector(coarse, &gcoarse));
1447: PetscCall(DMRestoreGlobalVector(fine, &gfine));
1449: nDM = comfine->nDM;
1450: 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);
1451: PetscCall(PetscCalloc1(nDM * nDM, &mats));
1452: if (v) PetscCall(PetscCalloc1(nDM, &vecs));
1454: /* loop over packed objects, handling one at a time */
1455: for (nextc = comcoarse->next, nextf = comfine->next, i = 0; nextc; nextc = nextc->next, nextf = nextf->next, i++) {
1456: if (!v) PetscCall(DMCreateInterpolation(nextc->dm, nextf->dm, &mats[i * nDM + i], NULL));
1457: else PetscCall(DMCreateInterpolation(nextc->dm, nextf->dm, &mats[i * nDM + i], &vecs[i]));
1458: }
1459: PetscCall(MatCreateNest(PetscObjectComm((PetscObject)fine), nDM, NULL, nDM, NULL, mats, A));
1460: if (v) PetscCall(VecCreateNest(PetscObjectComm((PetscObject)fine), nDM, NULL, vecs, v));
1461: for (i = 0; i < nDM * nDM; i++) PetscCall(MatDestroy(&mats[i]));
1462: PetscCall(PetscFree(mats));
1463: if (v) {
1464: for (i = 0; i < nDM; i++) PetscCall(VecDestroy(&vecs[i]));
1465: PetscCall(PetscFree(vecs));
1466: }
1467: PetscFunctionReturn(PETSC_SUCCESS);
1468: }
1470: static PetscErrorCode DMGetLocalToGlobalMapping_Composite(DM dm)
1471: {
1472: DM_Composite *com = (DM_Composite *)dm->data;
1473: ISLocalToGlobalMapping *ltogs;
1474: PetscInt i;
1476: PetscFunctionBegin;
1477: /* Set the ISLocalToGlobalMapping on the new matrix */
1478: PetscCall(DMCompositeGetISLocalToGlobalMappings(dm, <ogs));
1479: PetscCall(ISLocalToGlobalMappingConcatenate(PetscObjectComm((PetscObject)dm), com->nDM, ltogs, &dm->ltogmap));
1480: for (i = 0; i < com->nDM; i++) PetscCall(ISLocalToGlobalMappingDestroy(<ogs[i]));
1481: PetscCall(PetscFree(ltogs));
1482: PetscFunctionReturn(PETSC_SUCCESS);
1483: }
1485: static PetscErrorCode DMCreateColoring_Composite(DM dm, ISColoringType ctype, ISColoring *coloring)
1486: {
1487: PetscInt n, i, cnt;
1488: ISColoringValue *colors;
1489: PetscBool dense = PETSC_FALSE;
1490: ISColoringValue maxcol = 0;
1491: DM_Composite *com = (DM_Composite *)dm->data;
1493: PetscFunctionBegin;
1495: PetscCheck(ctype != IS_COLORING_LOCAL, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Only global coloring supported");
1496: PetscCheck(ctype == IS_COLORING_GLOBAL, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Unknown ISColoringType");
1497: n = com->n;
1498: PetscCall(PetscMalloc1(n, &colors)); /* freed in ISColoringDestroy() */
1500: PetscCall(PetscOptionsGetBool(((PetscObject)dm)->options, ((PetscObject)dm)->prefix, "-dmcomposite_dense_jacobian", &dense, NULL));
1501: if (dense) {
1502: PetscCall(ISColoringValueCast(com->N, &maxcol));
1503: for (i = 0; i < n; i++) PetscCall(ISColoringValueCast(com->rstart + i, colors + i));
1504: } else {
1505: struct DMCompositeLink *next = com->next;
1506: PetscMPIInt rank;
1508: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
1509: cnt = 0;
1510: while (next) {
1511: ISColoring lcoloring;
1513: PetscCall(DMCreateColoring(next->dm, IS_COLORING_GLOBAL, &lcoloring));
1514: for (i = 0; i < lcoloring->N; i++) PetscCall(ISColoringValueCast(maxcol + lcoloring->colors[i], colors + cnt++));
1515: maxcol += lcoloring->n;
1516: PetscCall(ISColoringDestroy(&lcoloring));
1517: next = next->next;
1518: }
1519: }
1520: PetscCall(ISColoringCreate(PetscObjectComm((PetscObject)dm), maxcol, n, colors, PETSC_OWN_POINTER, coloring));
1521: PetscFunctionReturn(PETSC_SUCCESS);
1522: }
1524: static PetscErrorCode DMGlobalToLocalBegin_Composite(DM dm, Vec gvec, InsertMode mode, Vec lvec)
1525: {
1526: struct DMCompositeLink *next;
1527: const PetscScalar *garray, *garrayhead;
1528: PetscScalar *larray, *larrayhead;
1529: DM_Composite *com = (DM_Composite *)dm->data;
1531: PetscFunctionBegin;
1535: if (!com->setup) PetscCall(DMSetUp(dm));
1537: PetscCall(VecGetArrayRead(gvec, &garray));
1538: garrayhead = garray;
1539: PetscCall(VecGetArray(lvec, &larray));
1540: larrayhead = larray;
1542: /* loop over packed objects, handling one at a time */
1543: next = com->next;
1544: while (next) {
1545: Vec local, global;
1546: PetscInt N;
1548: PetscCall(DMGetGlobalVector(next->dm, &global));
1549: PetscCall(VecGetLocalSize(global, &N));
1550: PetscCall(VecPlaceArray(global, garrayhead));
1551: PetscCall(DMGetLocalVector(next->dm, &local));
1552: PetscCall(VecPlaceArray(local, larrayhead));
1553: PetscCall(DMGlobalToLocalBegin(next->dm, global, mode, local));
1554: PetscCall(DMGlobalToLocalEnd(next->dm, global, mode, local));
1555: PetscCall(VecResetArray(global));
1556: PetscCall(VecResetArray(local));
1557: PetscCall(DMRestoreGlobalVector(next->dm, &global));
1558: PetscCall(DMRestoreLocalVector(next->dm, &local));
1560: larrayhead += next->nlocal;
1561: garrayhead += next->n;
1562: next = next->next;
1563: }
1564: PetscCall(VecRestoreArrayRead(gvec, &garray));
1565: PetscCall(VecRestoreArray(lvec, &larray));
1566: PetscFunctionReturn(PETSC_SUCCESS);
1567: }
1569: static PetscErrorCode DMGlobalToLocalEnd_Composite(DM dm, Vec gvec, InsertMode mode, Vec lvec)
1570: {
1571: PetscFunctionBegin;
1575: PetscFunctionReturn(PETSC_SUCCESS);
1576: }
1578: static PetscErrorCode DMLocalToGlobalBegin_Composite(DM dm, Vec lvec, InsertMode mode, Vec gvec)
1579: {
1580: struct DMCompositeLink *next;
1581: const PetscScalar *larray, *larrayhead;
1582: PetscScalar *garray, *garrayhead;
1583: DM_Composite *com = (DM_Composite *)dm->data;
1585: PetscFunctionBegin;
1590: if (!com->setup) PetscCall(DMSetUp(dm));
1592: PetscCall(VecGetArrayRead(lvec, &larray));
1593: larrayhead = larray;
1594: PetscCall(VecGetArray(gvec, &garray));
1595: garrayhead = garray;
1597: /* loop over packed objects, handling one at a time */
1598: next = com->next;
1599: while (next) {
1600: Vec global, local;
1602: PetscCall(DMGetLocalVector(next->dm, &local));
1603: PetscCall(VecPlaceArray(local, larrayhead));
1604: PetscCall(DMGetGlobalVector(next->dm, &global));
1605: PetscCall(VecPlaceArray(global, garrayhead));
1606: PetscCall(DMLocalToGlobalBegin(next->dm, local, mode, global));
1607: PetscCall(DMLocalToGlobalEnd(next->dm, local, mode, global));
1608: PetscCall(VecResetArray(local));
1609: PetscCall(VecResetArray(global));
1610: PetscCall(DMRestoreGlobalVector(next->dm, &global));
1611: PetscCall(DMRestoreLocalVector(next->dm, &local));
1613: garrayhead += next->n;
1614: larrayhead += next->nlocal;
1615: next = next->next;
1616: }
1618: PetscCall(VecRestoreArray(gvec, &garray));
1619: PetscCall(VecRestoreArrayRead(lvec, &larray));
1620: PetscFunctionReturn(PETSC_SUCCESS);
1621: }
1623: static PetscErrorCode DMLocalToGlobalEnd_Composite(DM dm, Vec lvec, InsertMode mode, Vec gvec)
1624: {
1625: PetscFunctionBegin;
1629: PetscFunctionReturn(PETSC_SUCCESS);
1630: }
1632: static PetscErrorCode DMLocalToLocalBegin_Composite(DM dm, Vec vec1, InsertMode mode, Vec vec2)
1633: {
1634: struct DMCompositeLink *next;
1635: const PetscScalar *array1, *array1head;
1636: PetscScalar *array2, *array2head;
1637: DM_Composite *com = (DM_Composite *)dm->data;
1639: PetscFunctionBegin;
1644: if (!com->setup) PetscCall(DMSetUp(dm));
1646: PetscCall(VecGetArrayRead(vec1, &array1));
1647: array1head = array1;
1648: PetscCall(VecGetArray(vec2, &array2));
1649: array2head = array2;
1651: /* loop over packed objects, handling one at a time */
1652: next = com->next;
1653: while (next) {
1654: Vec local1, local2;
1656: PetscCall(DMGetLocalVector(next->dm, &local1));
1657: PetscCall(VecPlaceArray(local1, array1head));
1658: PetscCall(DMGetLocalVector(next->dm, &local2));
1659: PetscCall(VecPlaceArray(local2, array2head));
1660: PetscCall(DMLocalToLocalBegin(next->dm, local1, mode, local2));
1661: PetscCall(DMLocalToLocalEnd(next->dm, local1, mode, local2));
1662: PetscCall(VecResetArray(local2));
1663: PetscCall(DMRestoreLocalVector(next->dm, &local2));
1664: PetscCall(VecResetArray(local1));
1665: PetscCall(DMRestoreLocalVector(next->dm, &local1));
1667: array1head += next->nlocal;
1668: array2head += next->nlocal;
1669: next = next->next;
1670: }
1672: PetscCall(VecRestoreArrayRead(vec1, &array1));
1673: PetscCall(VecRestoreArray(vec2, &array2));
1674: PetscFunctionReturn(PETSC_SUCCESS);
1675: }
1677: static PetscErrorCode DMLocalToLocalEnd_Composite(DM dm, Vec lvec, InsertMode mode, Vec gvec)
1678: {
1679: PetscFunctionBegin;
1683: PetscFunctionReturn(PETSC_SUCCESS);
1684: }
1686: /*MC
1687: DMCOMPOSITE = "composite" - A `DM` object that is used to manage data for a collection of `DM`
1689: Level: intermediate
1691: .seealso: `DMType`, `DM`, `DMDACreate()`, `DMCreate()`, `DMSetType()`, `DMCompositeCreate()`
1692: M*/
1694: PETSC_EXTERN PetscErrorCode DMCreate_Composite(DM p)
1695: {
1696: DM_Composite *com;
1698: PetscFunctionBegin;
1699: PetscCall(PetscNew(&com));
1700: p->data = com;
1701: com->n = 0;
1702: com->nghost = 0;
1703: com->next = NULL;
1704: com->nDM = 0;
1706: p->ops->createglobalvector = DMCreateGlobalVector_Composite;
1707: p->ops->createlocalvector = DMCreateLocalVector_Composite;
1708: p->ops->getlocaltoglobalmapping = DMGetLocalToGlobalMapping_Composite;
1709: p->ops->createfieldis = DMCreateFieldIS_Composite;
1710: p->ops->createfielddecomposition = DMCreateFieldDecomposition_Composite;
1711: p->ops->refine = DMRefine_Composite;
1712: p->ops->coarsen = DMCoarsen_Composite;
1713: p->ops->createinterpolation = DMCreateInterpolation_Composite;
1714: p->ops->creatematrix = DMCreateMatrix_Composite;
1715: p->ops->getcoloring = DMCreateColoring_Composite;
1716: p->ops->globaltolocalbegin = DMGlobalToLocalBegin_Composite;
1717: p->ops->globaltolocalend = DMGlobalToLocalEnd_Composite;
1718: p->ops->localtoglobalbegin = DMLocalToGlobalBegin_Composite;
1719: p->ops->localtoglobalend = DMLocalToGlobalEnd_Composite;
1720: p->ops->localtolocalbegin = DMLocalToLocalBegin_Composite;
1721: p->ops->localtolocalend = DMLocalToLocalEnd_Composite;
1722: p->ops->destroy = DMDestroy_Composite;
1723: p->ops->view = DMView_Composite;
1724: p->ops->setup = DMSetUp_Composite;
1726: PetscCall(PetscObjectComposeFunction((PetscObject)p, "DMSetUpGLVisViewer_C", DMSetUpGLVisViewer_Composite));
1727: PetscFunctionReturn(PETSC_SUCCESS);
1728: }
1730: /*@
1731: DMCompositeCreate - Creates a `DMCOMPOSITE`, used to generate "composite"
1732: vectors made up of several subvectors.
1734: Collective
1736: Input Parameter:
1737: . comm - the processors that will share the global vector
1739: Output Parameter:
1740: . packer - the `DMCOMPOSITE` object
1742: Level: advanced
1744: .seealso: `DMCOMPOSITE`, `DM`, `DMDestroy()`, `DMCompositeAddDM()`, `DMCompositeScatter()`, `DMCreate()`
1745: `DMCompositeGather()`, `DMCreateGlobalVector()`, `DMCompositeGetISLocalToGlobalMappings()`, `DMCompositeGetAccess()`
1746: `DMCompositeGetLocalVectors()`, `DMCompositeRestoreLocalVectors()`, `DMCompositeGetEntries()`
1747: @*/
1748: PetscErrorCode DMCompositeCreate(MPI_Comm comm, DM *packer)
1749: {
1750: PetscFunctionBegin;
1751: PetscAssertPointer(packer, 2);
1752: PetscCall(DMCreate(comm, packer));
1753: PetscCall(DMSetType(*packer, DMCOMPOSITE));
1754: PetscFunctionReturn(PETSC_SUCCESS);
1755: }