Actual source code: pack.c
2: #include <../src/dm/impls/composite/packimpl.h>
3: #include <petsc/private/isimpl.h>
4: #include <petsc/private/glvisviewerimpl.h>
5: #include <petscds.h>
7: /*@C
8: DMCompositeSetCoupling - Sets user provided routines that compute the coupling between the
9: separate components (DMs) in a DMto build the correct matrix nonzero structure.
11: Logically Collective
13: Input Parameters:
14: + dm - the composite object
15: - formcouplelocations - routine to set the nonzero locations in the matrix
17: Level: advanced
19: Not available from Fortran
21: Notes:
22: See DMSetApplicationContext() and DMGetApplicationContext() for how to get user information into
23: this routine
25: @*/
26: PetscErrorCode DMCompositeSetCoupling(DM dm,PetscErrorCode (*FormCoupleLocations)(DM,Mat,PetscInt*,PetscInt*,PetscInt,PetscInt,PetscInt,PetscInt))
27: {
28: DM_Composite *com = (DM_Composite*)dm->data;
29: PetscBool flg;
31: PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
33: com->FormCoupleLocations = FormCoupleLocations;
34: return 0;
35: }
37: PetscErrorCode DMDestroy_Composite(DM dm)
38: {
39: struct DMCompositeLink *next, *prev;
40: DM_Composite *com = (DM_Composite*)dm->data;
42: next = com->next;
43: while (next) {
44: prev = next;
45: next = next->next;
46: DMDestroy(&prev->dm);
47: PetscFree(prev->grstarts);
48: PetscFree(prev);
49: }
50: PetscObjectComposeFunction((PetscObject)dm,"DMSetUpGLVisViewer_C",NULL);
51: /* This was originally freed in DMDestroy(), but that prevents reference counting of backend objects */
52: PetscFree(com);
53: return 0;
54: }
56: PetscErrorCode DMView_Composite(DM dm,PetscViewer v)
57: {
58: PetscBool iascii;
59: DM_Composite *com = (DM_Composite*)dm->data;
61: PetscObjectTypeCompare((PetscObject)v,PETSCVIEWERASCII,&iascii);
62: if (iascii) {
63: struct DMCompositeLink *lnk = com->next;
64: PetscInt i;
66: PetscViewerASCIIPrintf(v,"DM (%s)\n",((PetscObject)dm)->prefix ? ((PetscObject)dm)->prefix : "no prefix");
67: PetscViewerASCIIPrintf(v," contains %D DMs\n",com->nDM);
68: PetscViewerASCIIPushTab(v);
69: for (i=0; lnk; lnk=lnk->next,i++) {
70: PetscViewerASCIIPrintf(v,"Link %D: DM of type %s\n",i,((PetscObject)lnk->dm)->type_name);
71: PetscViewerASCIIPushTab(v);
72: DMView(lnk->dm,v);
73: PetscViewerASCIIPopTab(v);
74: }
75: PetscViewerASCIIPopTab(v);
76: }
77: return 0;
78: }
80: /* --------------------------------------------------------------------------------------*/
81: PetscErrorCode DMSetUp_Composite(DM dm)
82: {
83: PetscInt nprev = 0;
84: PetscMPIInt rank,size;
85: DM_Composite *com = (DM_Composite*)dm->data;
86: struct DMCompositeLink *next = com->next;
87: PetscLayout map;
90: PetscLayoutCreate(PetscObjectComm((PetscObject)dm),&map);
91: PetscLayoutSetLocalSize(map,com->n);
92: PetscLayoutSetSize(map,PETSC_DETERMINE);
93: PetscLayoutSetBlockSize(map,1);
94: PetscLayoutSetUp(map);
95: PetscLayoutGetSize(map,&com->N);
96: PetscLayoutGetRange(map,&com->rstart,NULL);
97: PetscLayoutDestroy(&map);
99: /* now set the rstart for each linked vector */
100: MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);
101: MPI_Comm_size(PetscObjectComm((PetscObject)dm),&size);
102: while (next) {
103: next->rstart = nprev;
104: nprev += next->n;
105: next->grstart = com->rstart + next->rstart;
106: PetscMalloc1(size,&next->grstarts);
107: MPI_Allgather(&next->grstart,1,MPIU_INT,next->grstarts,1,MPIU_INT,PetscObjectComm((PetscObject)dm));
108: next = next->next;
109: }
110: com->setup = PETSC_TRUE;
111: return 0;
112: }
114: /* ----------------------------------------------------------------------------------*/
116: /*@
117: DMCompositeGetNumberDM - Get's the number of DM objects in the DMComposite
118: representation.
120: Not Collective
122: Input Parameter:
123: . dm - the packer object
125: Output Parameter:
126: . nDM - the number of DMs
128: Level: beginner
130: @*/
131: PetscErrorCode DMCompositeGetNumberDM(DM dm,PetscInt *nDM)
132: {
133: DM_Composite *com = (DM_Composite*)dm->data;
134: PetscBool flg;
137: PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
139: *nDM = com->nDM;
140: return 0;
141: }
143: /*@C
144: DMCompositeGetAccess - Allows one to access the individual packed vectors in their global
145: representation.
147: Collective on dm
149: Input Parameters:
150: + dm - the packer object
151: - gvec - the global vector
153: Output Parameters:
154: . Vec* ... - the packed parallel vectors, NULL for those that are not needed
156: Notes:
157: Use DMCompositeRestoreAccess() to return the vectors when you no longer need them
159: Fortran Notes:
161: Fortran callers must use numbered versions of this routine, e.g., DMCompositeGetAccess4(dm,gvec,vec1,vec2,vec3,vec4)
162: or use the alternative interface DMCompositeGetAccessArray().
164: Level: advanced
166: .seealso: DMCompositeGetEntries(), DMCompositeScatter()
167: @*/
168: PetscErrorCode DMCompositeGetAccess(DM dm,Vec gvec,...)
169: {
170: va_list Argp;
171: struct DMCompositeLink *next;
172: DM_Composite *com = (DM_Composite*)dm->data;
173: PetscInt readonly;
174: PetscBool flg;
178: PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
180: next = com->next;
181: if (!com->setup) {
182: DMSetUp(dm);
183: }
185: VecLockGet(gvec,&readonly);
186: /* loop over packed objects, handling one at at time */
187: va_start(Argp,gvec);
188: while (next) {
189: Vec *vec;
190: vec = va_arg(Argp, Vec*);
191: if (vec) {
192: DMGetGlobalVector(next->dm,vec);
193: if (readonly) {
194: const PetscScalar *array;
195: VecGetArrayRead(gvec,&array);
196: VecPlaceArray(*vec,array+next->rstart);
197: VecLockReadPush(*vec);
198: VecRestoreArrayRead(gvec,&array);
199: } else {
200: PetscScalar *array;
201: VecGetArray(gvec,&array);
202: VecPlaceArray(*vec,array+next->rstart);
203: VecRestoreArray(gvec,&array);
204: }
205: }
206: next = next->next;
207: }
208: va_end(Argp);
209: return 0;
210: }
212: /*@C
213: DMCompositeGetAccessArray - Allows one to access the individual packed vectors in their global
214: representation.
216: Collective on dm
218: Input Parameters:
219: + dm - the packer object
220: . pvec - packed vector
221: . nwanted - number of vectors wanted
222: - wanted - sorted array of vectors wanted, or NULL to get all vectors
224: Output Parameters:
225: . vecs - array of requested global vectors (must be allocated)
227: Notes:
228: Use DMCompositeRestoreAccessArray() to return the vectors when you no longer need them
230: Level: advanced
232: .seealso: DMCompositeGetAccess(), DMCompositeGetEntries(), DMCompositeScatter(), DMCompositeGather()
233: @*/
234: PetscErrorCode DMCompositeGetAccessArray(DM dm,Vec pvec,PetscInt nwanted,const PetscInt *wanted,Vec *vecs)
235: {
236: struct DMCompositeLink *link;
237: PetscInt i,wnum;
238: DM_Composite *com = (DM_Composite*)dm->data;
239: PetscInt readonly;
240: PetscBool flg;
244: PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
246: if (!com->setup) {
247: DMSetUp(dm);
248: }
250: VecLockGet(pvec,&readonly);
251: for (i=0,wnum=0,link=com->next; link && wnum<nwanted; i++,link=link->next) {
252: if (!wanted || i == wanted[wnum]) {
253: Vec v;
254: DMGetGlobalVector(link->dm,&v);
255: if (readonly) {
256: const PetscScalar *array;
257: VecGetArrayRead(pvec,&array);
258: VecPlaceArray(v,array+link->rstart);
259: VecLockReadPush(v);
260: VecRestoreArrayRead(pvec,&array);
261: } else {
262: PetscScalar *array;
263: VecGetArray(pvec,&array);
264: VecPlaceArray(v,array+link->rstart);
265: VecRestoreArray(pvec,&array);
266: }
267: vecs[wnum++] = v;
268: }
269: }
270: return 0;
271: }
273: /*@C
274: DMCompositeGetLocalAccessArray - Allows one to access the individual
275: packed vectors in their local representation.
277: Collective on dm.
279: Input Parameters:
280: + dm - the packer object
281: . pvec - packed vector
282: . nwanted - number of vectors wanted
283: - wanted - sorted array of vectors wanted, or NULL to get all vectors
285: Output Parameters:
286: . vecs - array of requested local vectors (must be allocated)
288: Notes:
289: Use DMCompositeRestoreLocalAccessArray() to return the vectors
290: when you no longer need them.
292: Level: advanced
294: .seealso: DMCompositeRestoreLocalAccessArray(), DMCompositeGetAccess(),
295: DMCompositeGetEntries(), DMCompositeScatter(), DMCompositeGather()
296: @*/
297: PetscErrorCode DMCompositeGetLocalAccessArray(DM dm,Vec pvec,PetscInt nwanted,const PetscInt *wanted,Vec *vecs)
298: {
299: struct DMCompositeLink *link;
300: PetscInt i,wnum;
301: DM_Composite *com = (DM_Composite*)dm->data;
302: PetscInt readonly;
303: PetscInt nlocal = 0;
304: PetscBool flg;
308: PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
310: if (!com->setup) {
311: DMSetUp(dm);
312: }
314: VecLockGet(pvec,&readonly);
315: for (i=0,wnum=0,link=com->next; link && wnum<nwanted; i++,link=link->next) {
316: if (!wanted || i == wanted[wnum]) {
317: Vec v;
318: DMGetLocalVector(link->dm,&v);
319: if (readonly) {
320: const PetscScalar *array;
321: VecGetArrayRead(pvec,&array);
322: VecPlaceArray(v,array+nlocal);
323: // this method does not make sense. The local vectors are not updated with a global-to-local and the user can not do it because it is locked
324: VecLockReadPush(v);
325: VecRestoreArrayRead(pvec,&array);
326: } else {
327: PetscScalar *array;
328: VecGetArray(pvec,&array);
329: VecPlaceArray(v,array+nlocal);
330: VecRestoreArray(pvec,&array);
331: }
332: vecs[wnum++] = v;
333: }
335: nlocal += link->nlocal;
336: }
338: return 0;
339: }
341: /*@C
342: DMCompositeRestoreAccess - Returns the vectors obtained with DMCompositeGetAccess()
343: representation.
345: Collective on dm
347: Input Parameters:
348: + dm - the packer object
349: . gvec - the global vector
350: - Vec* ... - the individual parallel vectors, NULL for those that are not needed
352: Level: advanced
354: .seealso DMCompositeAddDM(), DMCreateGlobalVector(),
355: DMCompositeGather(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeScatter(),
356: DMCompositeRestoreAccess(), DMCompositeGetAccess()
358: @*/
359: PetscErrorCode DMCompositeRestoreAccess(DM dm,Vec gvec,...)
360: {
361: va_list Argp;
362: struct DMCompositeLink *next;
363: DM_Composite *com = (DM_Composite*)dm->data;
364: PetscInt readonly;
365: PetscBool flg;
369: PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
371: next = com->next;
372: if (!com->setup) {
373: DMSetUp(dm);
374: }
376: VecLockGet(gvec,&readonly);
377: /* loop over packed objects, handling one at at time */
378: va_start(Argp,gvec);
379: while (next) {
380: Vec *vec;
381: vec = va_arg(Argp, Vec*);
382: if (vec) {
383: VecResetArray(*vec);
384: if (readonly) {
385: VecLockReadPop(*vec);
386: }
387: DMRestoreGlobalVector(next->dm,vec);
388: }
389: next = next->next;
390: }
391: va_end(Argp);
392: return 0;
393: }
395: /*@C
396: DMCompositeRestoreAccessArray - Returns the vectors obtained with DMCompositeGetAccessArray()
398: Collective on dm
400: Input Parameters:
401: + dm - the packer object
402: . pvec - packed vector
403: . nwanted - number of vectors wanted
404: . wanted - sorted array of vectors wanted, or NULL to get all vectors
405: - vecs - array of global vectors to return
407: Level: advanced
409: .seealso: DMCompositeRestoreAccess(), DMCompositeRestoreEntries(), DMCompositeScatter(), DMCompositeGather()
410: @*/
411: PetscErrorCode DMCompositeRestoreAccessArray(DM dm,Vec pvec,PetscInt nwanted,const PetscInt *wanted,Vec *vecs)
412: {
413: struct DMCompositeLink *link;
414: PetscInt i,wnum;
415: DM_Composite *com = (DM_Composite*)dm->data;
416: PetscInt readonly;
417: PetscBool flg;
421: PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
423: if (!com->setup) {
424: DMSetUp(dm);
425: }
427: VecLockGet(pvec,&readonly);
428: for (i=0,wnum=0,link=com->next; link && wnum<nwanted; i++,link=link->next) {
429: if (!wanted || i == wanted[wnum]) {
430: VecResetArray(vecs[wnum]);
431: if (readonly) {
432: VecLockReadPop(vecs[wnum]);
433: }
434: DMRestoreGlobalVector(link->dm,&vecs[wnum]);
435: wnum++;
436: }
437: }
438: return 0;
439: }
441: /*@C
442: DMCompositeRestoreLocalAccessArray - Returns the vectors obtained with DMCompositeGetLocalAccessArray().
444: Collective on dm.
446: Input Parameters:
447: + dm - the packer object
448: . pvec - packed vector
449: . nwanted - number of vectors wanted
450: . wanted - sorted array of vectors wanted, or NULL to restore all vectors
451: - vecs - array of local vectors to return
453: Level: advanced
455: Notes:
456: nwanted and wanted must match the values given to DMCompositeGetLocalAccessArray()
457: otherwise the call will fail.
459: .seealso: DMCompositeGetLocalAccessArray(), DMCompositeRestoreAccessArray(),
460: DMCompositeRestoreAccess(), DMCompositeRestoreEntries(),
461: DMCompositeScatter(), DMCompositeGather()
462: @*/
463: PetscErrorCode DMCompositeRestoreLocalAccessArray(DM dm,Vec pvec,PetscInt nwanted,const PetscInt *wanted,Vec *vecs)
464: {
465: struct DMCompositeLink *link;
466: PetscInt i,wnum;
467: DM_Composite *com = (DM_Composite*)dm->data;
468: PetscInt readonly;
469: PetscBool flg;
473: PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
475: if (!com->setup) {
476: DMSetUp(dm);
477: }
479: VecLockGet(pvec,&readonly);
480: for (i=0,wnum=0,link=com->next; link && wnum<nwanted; i++,link=link->next) {
481: if (!wanted || i == wanted[wnum]) {
482: VecResetArray(vecs[wnum]);
483: if (readonly) {
484: VecLockReadPop(vecs[wnum]);
485: }
486: DMRestoreLocalVector(link->dm,&vecs[wnum]);
487: wnum++;
488: }
489: }
490: return 0;
491: }
493: /*@C
494: DMCompositeScatter - Scatters from a global packed vector into its individual local vectors
496: Collective on dm
498: Input Parameters:
499: + dm - the packer object
500: . gvec - the global vector
501: - Vec ... - the individual sequential vectors, NULL for those that are not needed
503: Level: advanced
505: Notes:
506: DMCompositeScatterArray() is a non-variadic alternative that is often more convenient for library callers and is
507: accessible from Fortran.
509: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(),
510: DMCompositeGather(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
511: DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(), DMCompositeGetEntries()
512: DMCompositeScatterArray()
514: @*/
515: PetscErrorCode DMCompositeScatter(DM dm,Vec gvec,...)
516: {
517: va_list Argp;
518: struct DMCompositeLink *next;
519: PetscInt cnt;
520: DM_Composite *com = (DM_Composite*)dm->data;
521: PetscBool flg;
525: PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
527: if (!com->setup) {
528: DMSetUp(dm);
529: }
531: /* loop over packed objects, handling one at at time */
532: va_start(Argp,gvec);
533: for (cnt=3,next=com->next; next; cnt++,next=next->next) {
534: Vec local;
535: local = va_arg(Argp, Vec);
536: if (local) {
537: Vec global;
538: const PetscScalar *array;
540: DMGetGlobalVector(next->dm,&global);
541: VecGetArrayRead(gvec,&array);
542: VecPlaceArray(global,array+next->rstart);
543: DMGlobalToLocalBegin(next->dm,global,INSERT_VALUES,local);
544: DMGlobalToLocalEnd(next->dm,global,INSERT_VALUES,local);
545: VecRestoreArrayRead(gvec,&array);
546: VecResetArray(global);
547: DMRestoreGlobalVector(next->dm,&global);
548: }
549: }
550: va_end(Argp);
551: return 0;
552: }
554: /*@
555: DMCompositeScatterArray - Scatters from a global packed vector into its individual local vectors
557: Collective on dm
559: Input Parameters:
560: + dm - the packer object
561: . gvec - the global vector
562: - lvecs - array of local vectors, NULL for any that are not needed
564: Level: advanced
566: Note:
567: This is a non-variadic alternative to DMCompositeScatter()
569: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector()
570: DMCompositeGather(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
571: DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(), DMCompositeGetEntries()
573: @*/
574: PetscErrorCode DMCompositeScatterArray(DM dm,Vec gvec,Vec *lvecs)
575: {
576: struct DMCompositeLink *next;
577: PetscInt i;
578: DM_Composite *com = (DM_Composite*)dm->data;
579: PetscBool flg;
583: PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
585: if (!com->setup) {
586: DMSetUp(dm);
587: }
589: /* loop over packed objects, handling one at at time */
590: for (i=0,next=com->next; next; next=next->next,i++) {
591: if (lvecs[i]) {
592: Vec global;
593: const PetscScalar *array;
595: DMGetGlobalVector(next->dm,&global);
596: VecGetArrayRead(gvec,&array);
597: VecPlaceArray(global,(PetscScalar*)array+next->rstart);
598: DMGlobalToLocalBegin(next->dm,global,INSERT_VALUES,lvecs[i]);
599: DMGlobalToLocalEnd(next->dm,global,INSERT_VALUES,lvecs[i]);
600: VecRestoreArrayRead(gvec,&array);
601: VecResetArray(global);
602: DMRestoreGlobalVector(next->dm,&global);
603: }
604: }
605: return 0;
606: }
608: /*@C
609: DMCompositeGather - Gathers into a global packed vector from its individual local vectors
611: Collective on dm
613: Input Parameters:
614: + dm - the packer object
615: . gvec - the global vector
616: . imode - INSERT_VALUES or ADD_VALUES
617: - Vec ... - the individual sequential vectors, NULL for any that are not needed
619: Level: advanced
621: Not available from Fortran, Fortran users can use DMCompositeGatherArray()
623: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(),
624: DMCompositeScatter(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
625: DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(), DMCompositeGetEntries()
627: @*/
628: PetscErrorCode DMCompositeGather(DM dm,InsertMode imode,Vec gvec,...)
629: {
630: va_list Argp;
631: struct DMCompositeLink *next;
632: DM_Composite *com = (DM_Composite*)dm->data;
633: PetscInt cnt;
634: PetscBool flg;
638: PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
640: if (!com->setup) {
641: DMSetUp(dm);
642: }
644: /* loop over packed objects, handling one at at time */
645: va_start(Argp,gvec);
646: for (cnt=3,next=com->next; next; cnt++,next=next->next) {
647: Vec local;
648: local = va_arg(Argp, Vec);
649: if (local) {
650: PetscScalar *array;
651: Vec global;
653: DMGetGlobalVector(next->dm,&global);
654: VecGetArray(gvec,&array);
655: VecPlaceArray(global,array+next->rstart);
656: DMLocalToGlobalBegin(next->dm,local,imode,global);
657: DMLocalToGlobalEnd(next->dm,local,imode,global);
658: VecRestoreArray(gvec,&array);
659: VecResetArray(global);
660: DMRestoreGlobalVector(next->dm,&global);
661: }
662: }
663: va_end(Argp);
664: return 0;
665: }
667: /*@
668: DMCompositeGatherArray - Gathers into a global packed vector from its individual local vectors
670: Collective on dm
672: Input Parameters:
673: + dm - the packer object
674: . gvec - the global vector
675: . imode - INSERT_VALUES or ADD_VALUES
676: - lvecs - the individual sequential vectors, NULL for any that are not needed
678: Level: advanced
680: Notes:
681: This is a non-variadic alternative to DMCompositeGather().
683: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(),
684: DMCompositeScatter(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
685: DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(), DMCompositeGetEntries(),
686: @*/
687: PetscErrorCode DMCompositeGatherArray(DM dm,InsertMode imode,Vec gvec,Vec *lvecs)
688: {
689: struct DMCompositeLink *next;
690: DM_Composite *com = (DM_Composite*)dm->data;
691: PetscInt i;
692: PetscBool flg;
696: PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
698: if (!com->setup) {
699: DMSetUp(dm);
700: }
702: /* loop over packed objects, handling one at at time */
703: for (next=com->next,i=0; next; next=next->next,i++) {
704: if (lvecs[i]) {
705: PetscScalar *array;
706: Vec global;
708: DMGetGlobalVector(next->dm,&global);
709: VecGetArray(gvec,&array);
710: VecPlaceArray(global,array+next->rstart);
711: DMLocalToGlobalBegin(next->dm,lvecs[i],imode,global);
712: DMLocalToGlobalEnd(next->dm,lvecs[i],imode,global);
713: VecRestoreArray(gvec,&array);
714: VecResetArray(global);
715: DMRestoreGlobalVector(next->dm,&global);
716: }
717: }
718: return 0;
719: }
721: /*@
722: DMCompositeAddDM - adds a DM vector to a DMComposite
724: Collective on dm
726: Input Parameters:
727: + dmc - the DMComposite (packer) object
728: - dm - the DM object
730: Level: advanced
732: .seealso DMDestroy(), DMCompositeGather(), DMCompositeAddDM(), DMCreateGlobalVector(),
733: DMCompositeScatter(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
734: DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(), DMCompositeGetEntries()
736: @*/
737: PetscErrorCode DMCompositeAddDM(DM dmc,DM dm)
738: {
739: PetscInt n,nlocal;
740: struct DMCompositeLink *mine,*next;
741: Vec global,local;
742: DM_Composite *com = (DM_Composite*)dmc->data;
743: PetscBool flg;
747: PetscObjectTypeCompare((PetscObject)dmc,DMCOMPOSITE,&flg);
749: next = com->next;
752: /* create new link */
753: PetscNew(&mine);
754: PetscObjectReference((PetscObject)dm);
755: DMGetGlobalVector(dm,&global);
756: VecGetLocalSize(global,&n);
757: DMRestoreGlobalVector(dm,&global);
758: DMGetLocalVector(dm,&local);
759: VecGetSize(local,&nlocal);
760: DMRestoreLocalVector(dm,&local);
762: mine->n = n;
763: mine->nlocal = nlocal;
764: mine->dm = dm;
765: mine->next = NULL;
766: com->n += n;
767: com->nghost += nlocal;
769: /* add to end of list */
770: if (!next) com->next = mine;
771: else {
772: while (next->next) next = next->next;
773: next->next = mine;
774: }
775: com->nDM++;
776: com->nmine++;
777: return 0;
778: }
780: #include <petscdraw.h>
781: PETSC_EXTERN PetscErrorCode VecView_MPI(Vec,PetscViewer);
782: PetscErrorCode VecView_DMComposite(Vec gvec,PetscViewer viewer)
783: {
784: DM dm;
785: struct DMCompositeLink *next;
786: PetscBool isdraw;
787: DM_Composite *com;
789: VecGetDM(gvec, &dm);
791: com = (DM_Composite*)dm->data;
792: next = com->next;
794: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
795: if (!isdraw) {
796: /* do I really want to call this? */
797: VecView_MPI(gvec,viewer);
798: } else {
799: PetscInt cnt = 0;
801: /* loop over packed objects, handling one at at time */
802: while (next) {
803: Vec vec;
804: const PetscScalar *array;
805: PetscInt bs;
807: /* Should use VecGetSubVector() eventually, but would need to forward the DM for that to work */
808: DMGetGlobalVector(next->dm,&vec);
809: VecGetArrayRead(gvec,&array);
810: VecPlaceArray(vec,(PetscScalar*)array+next->rstart);
811: VecRestoreArrayRead(gvec,&array);
812: VecView(vec,viewer);
813: VecResetArray(vec);
814: VecGetBlockSize(vec,&bs);
815: DMRestoreGlobalVector(next->dm,&vec);
816: PetscViewerDrawBaseAdd(viewer,bs);
817: cnt += bs;
818: next = next->next;
819: }
820: PetscViewerDrawBaseAdd(viewer,-cnt);
821: }
822: return 0;
823: }
825: PetscErrorCode DMCreateGlobalVector_Composite(DM dm,Vec *gvec)
826: {
827: DM_Composite *com = (DM_Composite*)dm->data;
830: DMSetFromOptions(dm);
831: DMSetUp(dm);
832: VecCreate(PetscObjectComm((PetscObject)dm),gvec);
833: VecSetType(*gvec,dm->vectype);
834: VecSetSizes(*gvec,com->n,com->N);
835: VecSetDM(*gvec, dm);
836: VecSetOperation(*gvec,VECOP_VIEW,(void (*)(void))VecView_DMComposite);
837: return 0;
838: }
840: PetscErrorCode DMCreateLocalVector_Composite(DM dm,Vec *lvec)
841: {
842: DM_Composite *com = (DM_Composite*)dm->data;
845: if (!com->setup) {
846: DMSetFromOptions(dm);
847: DMSetUp(dm);
848: }
849: VecCreate(PETSC_COMM_SELF,lvec);
850: VecSetType(*lvec,dm->vectype);
851: VecSetSizes(*lvec,com->nghost,PETSC_DECIDE);
852: VecSetDM(*lvec, dm);
853: return 0;
854: }
856: /*@C
857: DMCompositeGetISLocalToGlobalMappings - gets an ISLocalToGlobalMapping for each DM in the DMComposite, maps to the composite global space
859: Collective on DM
861: Input Parameter:
862: . dm - the packer object
864: Output Parameters:
865: . ltogs - the individual mappings for each packed vector. Note that this includes
866: all the ghost points that individual ghosted DMDA's may have.
868: Level: advanced
870: Notes:
871: Each entry of ltogs should be destroyed with ISLocalToGlobalMappingDestroy(), the ltogs array should be freed with PetscFree().
873: Not available from Fortran
875: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(),
876: DMCompositeGather(), DMCompositeCreate(), DMCompositeGetAccess(), DMCompositeScatter(),
877: DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(),DMCompositeGetEntries()
879: @*/
880: PetscErrorCode DMCompositeGetISLocalToGlobalMappings(DM dm,ISLocalToGlobalMapping **ltogs)
881: {
882: PetscInt i,*idx,n,cnt;
883: struct DMCompositeLink *next;
884: PetscMPIInt rank;
885: DM_Composite *com = (DM_Composite*)dm->data;
886: PetscBool flg;
889: PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
891: DMSetUp(dm);
892: PetscMalloc1(com->nDM,ltogs);
893: next = com->next;
894: MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);
896: /* loop over packed objects, handling one at at time */
897: cnt = 0;
898: while (next) {
899: ISLocalToGlobalMapping ltog;
900: PetscMPIInt size;
901: const PetscInt *suboff,*indices;
902: Vec global;
904: /* Get sub-DM global indices for each local dof */
905: DMGetLocalToGlobalMapping(next->dm,<og);
906: ISLocalToGlobalMappingGetSize(ltog,&n);
907: ISLocalToGlobalMappingGetIndices(ltog,&indices);
908: PetscMalloc1(n,&idx);
910: /* Get the offsets for the sub-DM global vector */
911: DMGetGlobalVector(next->dm,&global);
912: VecGetOwnershipRanges(global,&suboff);
913: MPI_Comm_size(PetscObjectComm((PetscObject)global),&size);
915: /* Shift the sub-DM definition of the global space to the composite global space */
916: for (i=0; i<n; i++) {
917: PetscInt subi = indices[i],lo = 0,hi = size,t;
918: /* There's no consensus on what a negative index means,
919: except for skipping when setting the values in vectors and matrices */
920: if (subi < 0) { idx[i] = subi - next->grstarts[rank]; continue; }
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: ISLocalToGlobalMappingRestoreIndices(ltog,&indices);
930: ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)dm),1,n,idx,PETSC_OWN_POINTER,&(*ltogs)[cnt]);
931: DMRestoreGlobalVector(next->dm,&global);
932: next = next->next;
933: cnt++;
934: }
935: return 0;
936: }
938: /*@C
939: DMCompositeGetLocalISs - Gets index sets for each component of a composite local vector
941: Not Collective
943: Input Parameter:
944: . dm - composite DM
946: Output Parameter:
947: . is - array of serial index sets for each 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 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: Not available from Fortran
964: .seealso: DMCompositeGetGlobalISs(), DMCompositeGetISLocalToGlobalMappings(), MatGetLocalSubMatrix(), MatCreateLocalRef()
965: @*/
966: PetscErrorCode DMCompositeGetLocalISs(DM dm,IS **is)
967: {
968: DM_Composite *com = (DM_Composite*)dm->data;
969: struct DMCompositeLink *link;
970: PetscInt cnt,start;
971: PetscBool flg;
975: PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
977: PetscMalloc1(com->nmine,is);
978: for (cnt=0,start=0,link=com->next; link; start+=link->nlocal,cnt++,link=link->next) {
979: PetscInt bs;
980: ISCreateStride(PETSC_COMM_SELF,link->nlocal,start,1,&(*is)[cnt]);
981: DMGetBlockSize(link->dm,&bs);
982: ISSetBlockSize((*is)[cnt],bs);
983: }
984: return 0;
985: }
987: /*@C
988: DMCompositeGetGlobalISs - Gets the index sets for each composed object
990: Collective on dm
992: Input Parameter:
993: . dm - the packer object
995: Output Parameters:
996: . is - the array of index sets
998: Level: advanced
1000: Notes:
1001: The is entries should be destroyed with ISDestroy(), the is array should be freed with PetscFree()
1003: These could be used to extract a subset of vector entries for a "multi-physics" preconditioner
1005: Use DMCompositeGetLocalISs() for index sets in the packed local numbering, and
1006: DMCompositeGetISLocalToGlobalMappings() for to map local sub-DM (including ghost) indices to packed global
1007: indices.
1009: Fortran Notes:
1011: The output argument 'is' must be an allocated array of sufficient length, which can be learned using DMCompositeGetNumberDM().
1013: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(),
1014: DMCompositeGather(), DMCompositeCreate(), DMCompositeGetAccess(), DMCompositeScatter(),
1015: 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;
1027: PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
1030: PetscMalloc1(com->nDM,is);
1031: next = com->next;
1032: MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);
1034: /* loop over packed objects, handling one at at time */
1035: while (next) {
1036: PetscDS prob;
1038: ISCreateStride(PetscObjectComm((PetscObject)dm),next->n,next->grstart,1,&(*is)[cnt]);
1039: DMGetDS(dm, &prob);
1040: if (prob) {
1041: MatNullSpace space;
1042: Mat pmat;
1043: PetscObject disc;
1044: PetscInt Nf;
1046: PetscDSGetNumFields(prob, &Nf);
1047: if (cnt < Nf) {
1048: PetscDSGetDiscretization(prob, cnt, &disc);
1049: PetscObjectQuery(disc, "nullspace", (PetscObject*) &space);
1050: if (space) PetscObjectCompose((PetscObject) (*is)[cnt], "nullspace", (PetscObject) space);
1051: PetscObjectQuery(disc, "nearnullspace", (PetscObject*) &space);
1052: if (space) PetscObjectCompose((PetscObject) (*is)[cnt], "nearnullspace", (PetscObject) space);
1053: PetscObjectQuery(disc, "pmat", (PetscObject*) &pmat);
1054: if (pmat) PetscObjectCompose((PetscObject) (*is)[cnt], "pmat", (PetscObject) pmat);
1055: }
1056: }
1057: cnt++;
1058: next = next->next;
1059: }
1060: return 0;
1061: }
1063: PetscErrorCode DMCreateFieldIS_Composite(DM dm, PetscInt *numFields,char ***fieldNames, IS **fields)
1064: {
1065: PetscInt nDM;
1066: DM *dms;
1067: PetscInt i;
1069: DMCompositeGetNumberDM(dm, &nDM);
1070: if (numFields) *numFields = nDM;
1071: DMCompositeGetGlobalISs(dm, fields);
1072: if (fieldNames) {
1073: PetscMalloc1(nDM, &dms);
1074: PetscMalloc1(nDM, fieldNames);
1075: DMCompositeGetEntriesArray(dm, dms);
1076: for (i=0; i<nDM; i++) {
1077: char buf[256];
1078: const char *splitname;
1080: /* Split naming precedence: object name, prefix, number */
1081: splitname = ((PetscObject) dm)->name;
1082: if (!splitname) {
1083: PetscObjectGetOptionsPrefix((PetscObject)dms[i],&splitname);
1084: if (splitname) {
1085: size_t len;
1086: PetscStrncpy(buf,splitname,sizeof(buf));
1087: buf[sizeof(buf) - 1] = 0;
1088: PetscStrlen(buf,&len);
1089: if (buf[len-1] == '_') buf[len-1] = 0; /* Remove trailing underscore if it was used */
1090: splitname = buf;
1091: }
1092: }
1093: if (!splitname) {
1094: PetscSNPrintf(buf,sizeof(buf),"%D",i);
1095: splitname = buf;
1096: }
1097: PetscStrallocpy(splitname,&(*fieldNames)[i]);
1098: }
1099: PetscFree(dms);
1100: }
1101: return 0;
1102: }
1104: /*
1105: This could take over from DMCreateFieldIS(), as it is more general,
1106: making DMCreateFieldIS() a special case -- calling with dmlist == NULL;
1107: At this point it's probably best to be less intrusive, however.
1108: */
1109: PetscErrorCode DMCreateFieldDecomposition_Composite(DM dm, PetscInt *len,char ***namelist, IS **islist, DM **dmlist)
1110: {
1111: PetscInt nDM;
1112: PetscInt i;
1114: DMCreateFieldIS_Composite(dm, len, namelist, islist);
1115: if (dmlist) {
1116: DMCompositeGetNumberDM(dm, &nDM);
1117: PetscMalloc1(nDM, dmlist);
1118: DMCompositeGetEntriesArray(dm, *dmlist);
1119: for (i=0; i<nDM; i++) {
1120: PetscObjectReference((PetscObject)((*dmlist)[i]));
1121: }
1122: }
1123: return 0;
1124: }
1126: /* -------------------------------------------------------------------------------------*/
1127: /*@C
1128: DMCompositeGetLocalVectors - Gets local vectors for each part of a DMComposite.
1129: Use DMCompositeRestoreLocalVectors() to return them.
1131: Not Collective
1133: Input Parameter:
1134: . dm - the packer object
1136: Output Parameter:
1137: . Vec ... - the individual sequential Vecs
1139: Level: advanced
1141: Not available from Fortran
1143: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(),
1144: DMCompositeGather(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
1145: DMCompositeRestoreLocalVectors(), DMCompositeScatter(), DMCompositeGetEntries()
1147: @*/
1148: PetscErrorCode DMCompositeGetLocalVectors(DM dm,...)
1149: {
1150: va_list Argp;
1151: struct DMCompositeLink *next;
1152: DM_Composite *com = (DM_Composite*)dm->data;
1153: PetscBool flg;
1156: PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
1158: next = com->next;
1159: /* loop over packed objects, handling one at at time */
1160: va_start(Argp,dm);
1161: while (next) {
1162: Vec *vec;
1163: vec = va_arg(Argp, Vec*);
1164: if (vec) DMGetLocalVector(next->dm,vec);
1165: next = next->next;
1166: }
1167: va_end(Argp);
1168: return 0;
1169: }
1171: /*@C
1172: DMCompositeRestoreLocalVectors - Restores local vectors for each part of a DMComposite.
1174: Not Collective
1176: Input Parameter:
1177: . dm - the packer object
1179: Output Parameter:
1180: . Vec ... - the individual sequential Vecs
1182: Level: advanced
1184: Not available from Fortran
1186: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(),
1187: DMCompositeGather(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
1188: DMCompositeGetLocalVectors(), DMCompositeScatter(), DMCompositeGetEntries()
1190: @*/
1191: PetscErrorCode DMCompositeRestoreLocalVectors(DM dm,...)
1192: {
1193: va_list Argp;
1194: struct DMCompositeLink *next;
1195: DM_Composite *com = (DM_Composite*)dm->data;
1196: PetscBool flg;
1199: PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
1201: next = com->next;
1202: /* loop over packed objects, handling one at at time */
1203: va_start(Argp,dm);
1204: while (next) {
1205: Vec *vec;
1206: vec = va_arg(Argp, Vec*);
1207: if (vec) DMRestoreLocalVector(next->dm,vec);
1208: next = next->next;
1209: }
1210: va_end(Argp);
1211: return 0;
1212: }
1214: /* -------------------------------------------------------------------------------------*/
1215: /*@C
1216: DMCompositeGetEntries - Gets the DM for each entry in a DMComposite.
1218: Not Collective
1220: Input Parameter:
1221: . dm - the packer object
1223: Output Parameter:
1224: . DM ... - the individual entries (DMs)
1226: Level: advanced
1228: Fortran Notes:
1229: Available as DMCompositeGetEntries() for one output DM, DMCompositeGetEntries2() for 2, etc
1231: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(), DMCompositeGetEntriesArray()
1232: DMCompositeGather(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
1233: DMCompositeRestoreLocalVectors(), DMCompositeGetLocalVectors(), DMCompositeScatter(),
1234: DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors()
1236: @*/
1237: PetscErrorCode DMCompositeGetEntries(DM dm,...)
1238: {
1239: va_list Argp;
1240: struct DMCompositeLink *next;
1241: DM_Composite *com = (DM_Composite*)dm->data;
1242: PetscBool flg;
1245: PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
1247: next = com->next;
1248: /* loop over packed objects, handling one at at time */
1249: va_start(Argp,dm);
1250: while (next) {
1251: DM *dmn;
1252: dmn = va_arg(Argp, DM*);
1253: if (dmn) *dmn = next->dm;
1254: next = next->next;
1255: }
1256: va_end(Argp);
1257: return 0;
1258: }
1260: /*@C
1261: DMCompositeGetEntriesArray - Gets the DM for each entry in a DMComposite.
1263: Not Collective
1265: Input Parameter:
1266: . dm - the packer object
1268: Output Parameter:
1269: . dms - array of sufficient length (see DMCompositeGetNumberDM()) to hold the individual DMs
1271: Level: advanced
1273: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(), DMCompositeGetEntries()
1274: DMCompositeGather(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
1275: DMCompositeRestoreLocalVectors(), DMCompositeGetLocalVectors(), DMCompositeScatter(),
1276: DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors()
1278: @*/
1279: PetscErrorCode DMCompositeGetEntriesArray(DM dm,DM dms[])
1280: {
1281: struct DMCompositeLink *next;
1282: DM_Composite *com = (DM_Composite*)dm->data;
1283: PetscInt i;
1284: PetscBool flg;
1287: PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
1289: /* loop over packed objects, handling one at at time */
1290: for (next=com->next,i=0; next; next=next->next,i++) dms[i] = next->dm;
1291: return 0;
1292: }
1294: typedef struct {
1295: DM dm;
1296: PetscViewer *subv;
1297: Vec *vecs;
1298: } GLVisViewerCtx;
1300: static PetscErrorCode DestroyGLVisViewerCtx_Private(void *vctx)
1301: {
1302: GLVisViewerCtx *ctx = (GLVisViewerCtx*)vctx;
1303: PetscInt i,n;
1305: DMCompositeGetNumberDM(ctx->dm,&n);
1306: for (i = 0; i < n; i++) {
1307: PetscViewerDestroy(&ctx->subv[i]);
1308: }
1309: PetscFree2(ctx->subv,ctx->vecs);
1310: DMDestroy(&ctx->dm);
1311: PetscFree(ctx);
1312: return 0;
1313: }
1315: static PetscErrorCode DMCompositeSampleGLVisFields_Private(PetscObject oX, PetscInt nf, PetscObject oXfield[], void *vctx)
1316: {
1317: Vec X = (Vec)oX;
1318: GLVisViewerCtx *ctx = (GLVisViewerCtx*)vctx;
1319: PetscInt i,n,cumf;
1321: DMCompositeGetNumberDM(ctx->dm,&n);
1322: DMCompositeGetAccessArray(ctx->dm,X,n,NULL,ctx->vecs);
1323: for (i = 0, cumf = 0; i < n; i++) {
1324: PetscErrorCode (*g2l)(PetscObject,PetscInt,PetscObject[],void*);
1325: void *fctx;
1326: PetscInt nfi;
1328: PetscViewerGLVisGetFields_Private(ctx->subv[i],&nfi,NULL,NULL,&g2l,NULL,&fctx);
1329: if (!nfi) continue;
1330: if (g2l) {
1331: (*g2l)((PetscObject)ctx->vecs[i],nfi,oXfield+cumf,fctx);
1332: } else {
1333: VecCopy(ctx->vecs[i],(Vec)(oXfield[cumf]));
1334: }
1335: cumf += nfi;
1336: }
1337: DMCompositeRestoreAccessArray(ctx->dm,X,n,NULL,ctx->vecs);
1338: return 0;
1339: }
1341: static PetscErrorCode DMSetUpGLVisViewer_Composite(PetscObject odm, PetscViewer viewer)
1342: {
1343: DM dm = (DM)odm, *dms;
1344: Vec *Ufds;
1345: GLVisViewerCtx *ctx;
1346: PetscInt i,n,tnf,*sdim;
1347: char **fecs;
1349: PetscNew(&ctx);
1350: PetscObjectReference((PetscObject)dm);
1351: ctx->dm = dm;
1352: DMCompositeGetNumberDM(dm,&n);
1353: PetscMalloc1(n,&dms);
1354: DMCompositeGetEntriesArray(dm,dms);
1355: PetscMalloc2(n,&ctx->subv,n,&ctx->vecs);
1356: for (i = 0, tnf = 0; i < n; i++) {
1357: PetscInt nf;
1359: PetscViewerCreate(PetscObjectComm(odm),&ctx->subv[i]);
1360: PetscViewerSetType(ctx->subv[i],PETSCVIEWERGLVIS);
1361: PetscViewerGLVisSetDM_Private(ctx->subv[i],(PetscObject)dms[i]);
1362: PetscViewerGLVisGetFields_Private(ctx->subv[i],&nf,NULL,NULL,NULL,NULL,NULL);
1363: tnf += nf;
1364: }
1365: PetscFree(dms);
1366: PetscMalloc3(tnf,&fecs,tnf,&sdim,tnf,&Ufds);
1367: for (i = 0, tnf = 0; i < n; i++) {
1368: PetscInt *sd,nf,f;
1369: const char **fec;
1370: Vec *Uf;
1372: PetscViewerGLVisGetFields_Private(ctx->subv[i],&nf,&fec,&sd,NULL,(PetscObject**)&Uf,NULL);
1373: for (f = 0; f < nf; f++) {
1374: PetscStrallocpy(fec[f],&fecs[tnf+f]);
1375: Ufds[tnf+f] = Uf[f];
1376: sdim[tnf+f] = sd[f];
1377: }
1378: tnf += nf;
1379: }
1380: PetscViewerGLVisSetFields(viewer,tnf,(const char**)fecs,sdim,DMCompositeSampleGLVisFields_Private,(PetscObject*)Ufds,ctx,DestroyGLVisViewerCtx_Private);
1381: for (i = 0; i < tnf; i++) {
1382: PetscFree(fecs[i]);
1383: }
1384: PetscFree3(fecs,sdim,Ufds);
1385: return 0;
1386: }
1388: PetscErrorCode DMRefine_Composite(DM dmi,MPI_Comm comm,DM *fine)
1389: {
1390: struct DMCompositeLink *next;
1391: DM_Composite *com = (DM_Composite*)dmi->data;
1392: DM dm;
1395: if (comm == MPI_COMM_NULL) {
1396: PetscObjectGetComm((PetscObject)dmi,&comm);
1397: }
1398: DMSetUp(dmi);
1399: next = com->next;
1400: DMCompositeCreate(comm,fine);
1402: /* loop over packed objects, handling one at at time */
1403: while (next) {
1404: DMRefine(next->dm,comm,&dm);
1405: DMCompositeAddDM(*fine,dm);
1406: PetscObjectDereference((PetscObject)dm);
1407: next = next->next;
1408: }
1409: return 0;
1410: }
1412: PetscErrorCode DMCoarsen_Composite(DM dmi,MPI_Comm comm,DM *fine)
1413: {
1414: struct DMCompositeLink *next;
1415: DM_Composite *com = (DM_Composite*)dmi->data;
1416: DM dm;
1419: DMSetUp(dmi);
1420: if (comm == MPI_COMM_NULL) {
1421: PetscObjectGetComm((PetscObject)dmi,&comm);
1422: }
1423: next = com->next;
1424: DMCompositeCreate(comm,fine);
1426: /* loop over packed objects, handling one at at time */
1427: while (next) {
1428: DMCoarsen(next->dm,comm,&dm);
1429: DMCompositeAddDM(*fine,dm);
1430: PetscObjectDereference((PetscObject)dm);
1431: next = next->next;
1432: }
1433: return 0;
1434: }
1436: PetscErrorCode DMCreateInterpolation_Composite(DM coarse,DM fine,Mat *A,Vec *v)
1437: {
1438: PetscInt m,n,M,N,nDM,i;
1439: struct DMCompositeLink *nextc;
1440: struct DMCompositeLink *nextf;
1441: Vec gcoarse,gfine,*vecs;
1442: DM_Composite *comcoarse = (DM_Composite*)coarse->data;
1443: DM_Composite *comfine = (DM_Composite*)fine->data;
1444: Mat *mats;
1448: DMSetUp(coarse);
1449: DMSetUp(fine);
1450: /* use global vectors only for determining matrix layout */
1451: DMGetGlobalVector(coarse,&gcoarse);
1452: DMGetGlobalVector(fine,&gfine);
1453: VecGetLocalSize(gcoarse,&n);
1454: VecGetLocalSize(gfine,&m);
1455: VecGetSize(gcoarse,&N);
1456: VecGetSize(gfine,&M);
1457: DMRestoreGlobalVector(coarse,&gcoarse);
1458: DMRestoreGlobalVector(fine,&gfine);
1460: nDM = comfine->nDM;
1462: PetscCalloc1(nDM*nDM,&mats);
1463: if (v) {
1464: PetscCalloc1(nDM,&vecs);
1465: }
1467: /* loop over packed objects, handling one at at time */
1468: for (nextc=comcoarse->next,nextf=comfine->next,i=0; nextc; nextc=nextc->next,nextf=nextf->next,i++) {
1469: if (!v) {
1470: DMCreateInterpolation(nextc->dm,nextf->dm,&mats[i*nDM+i],NULL);
1471: } else {
1472: DMCreateInterpolation(nextc->dm,nextf->dm,&mats[i*nDM+i],&vecs[i]);
1473: }
1474: }
1475: MatCreateNest(PetscObjectComm((PetscObject)fine),nDM,NULL,nDM,NULL,mats,A);
1476: if (v) {
1477: VecCreateNest(PetscObjectComm((PetscObject)fine),nDM,NULL,vecs,v);
1478: }
1479: for (i=0; i<nDM*nDM; i++) MatDestroy(&mats[i]);
1480: PetscFree(mats);
1481: if (v) {
1482: for (i=0; i<nDM; i++) VecDestroy(&vecs[i]);
1483: PetscFree(vecs);
1484: }
1485: return 0;
1486: }
1488: static PetscErrorCode DMGetLocalToGlobalMapping_Composite(DM dm)
1489: {
1490: DM_Composite *com = (DM_Composite*)dm->data;
1491: ISLocalToGlobalMapping *ltogs;
1492: PetscInt i;
1494: /* Set the ISLocalToGlobalMapping on the new matrix */
1495: DMCompositeGetISLocalToGlobalMappings(dm,<ogs);
1496: ISLocalToGlobalMappingConcatenate(PetscObjectComm((PetscObject)dm),com->nDM,ltogs,&dm->ltogmap);
1497: for (i=0; i<com->nDM; i++) ISLocalToGlobalMappingDestroy(<ogs[i]);
1498: PetscFree(ltogs);
1499: return 0;
1500: }
1502: PetscErrorCode DMCreateColoring_Composite(DM dm,ISColoringType ctype,ISColoring *coloring)
1503: {
1504: PetscInt n,i,cnt;
1505: ISColoringValue *colors;
1506: PetscBool dense = PETSC_FALSE;
1507: ISColoringValue maxcol = 0;
1508: DM_Composite *com = (DM_Composite*)dm->data;
1512: else if (ctype == IS_COLORING_GLOBAL) {
1513: n = com->n;
1514: } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unknown ISColoringType");
1515: PetscMalloc1(n,&colors); /* freed in ISColoringDestroy() */
1517: PetscOptionsGetBool(((PetscObject)dm)->options,((PetscObject)dm)->prefix,"-dmcomposite_dense_jacobian",&dense,NULL);
1518: if (dense) {
1519: for (i=0; i<n; i++) {
1520: colors[i] = (ISColoringValue)(com->rstart + i);
1521: }
1522: maxcol = com->N;
1523: } else {
1524: struct DMCompositeLink *next = com->next;
1525: PetscMPIInt rank;
1527: MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);
1528: cnt = 0;
1529: while (next) {
1530: ISColoring lcoloring;
1532: DMCreateColoring(next->dm,IS_COLORING_GLOBAL,&lcoloring);
1533: for (i=0; i<lcoloring->N; i++) {
1534: colors[cnt++] = maxcol + lcoloring->colors[i];
1535: }
1536: maxcol += lcoloring->n;
1537: ISColoringDestroy(&lcoloring);
1538: next = next->next;
1539: }
1540: }
1541: ISColoringCreate(PetscObjectComm((PetscObject)dm),maxcol,n,colors,PETSC_OWN_POINTER,coloring);
1542: return 0;
1543: }
1545: PetscErrorCode DMGlobalToLocalBegin_Composite(DM dm,Vec gvec,InsertMode mode,Vec lvec)
1546: {
1547: struct DMCompositeLink *next;
1548: PetscScalar *garray,*larray;
1549: DM_Composite *com = (DM_Composite*)dm->data;
1554: if (!com->setup) {
1555: DMSetUp(dm);
1556: }
1558: VecGetArray(gvec,&garray);
1559: VecGetArray(lvec,&larray);
1561: /* loop over packed objects, handling one at at time */
1562: next = com->next;
1563: while (next) {
1564: Vec local,global;
1565: PetscInt N;
1567: DMGetGlobalVector(next->dm,&global);
1568: VecGetLocalSize(global,&N);
1569: VecPlaceArray(global,garray);
1570: DMGetLocalVector(next->dm,&local);
1571: VecPlaceArray(local,larray);
1572: DMGlobalToLocalBegin(next->dm,global,mode,local);
1573: DMGlobalToLocalEnd(next->dm,global,mode,local);
1574: VecResetArray(global);
1575: VecResetArray(local);
1576: DMRestoreGlobalVector(next->dm,&global);
1577: DMRestoreLocalVector(next->dm,&local);
1579: larray += next->nlocal;
1580: garray += next->n;
1581: next = next->next;
1582: }
1584: VecRestoreArray(gvec,NULL);
1585: VecRestoreArray(lvec,NULL);
1586: return 0;
1587: }
1589: PetscErrorCode DMGlobalToLocalEnd_Composite(DM dm,Vec gvec,InsertMode mode,Vec lvec)
1590: {
1594: return 0;
1595: }
1597: PetscErrorCode DMLocalToGlobalBegin_Composite(DM dm,Vec lvec,InsertMode mode,Vec gvec)
1598: {
1599: struct DMCompositeLink *next;
1600: PetscScalar *larray,*garray;
1601: DM_Composite *com = (DM_Composite*)dm->data;
1607: if (!com->setup) {
1608: DMSetUp(dm);
1609: }
1611: VecGetArray(lvec,&larray);
1612: VecGetArray(gvec,&garray);
1614: /* loop over packed objects, handling one at at time */
1615: next = com->next;
1616: while (next) {
1617: Vec global,local;
1619: DMGetLocalVector(next->dm,&local);
1620: VecPlaceArray(local,larray);
1621: DMGetGlobalVector(next->dm,&global);
1622: VecPlaceArray(global,garray);
1623: DMLocalToGlobalBegin(next->dm,local,mode,global);
1624: DMLocalToGlobalEnd(next->dm,local,mode,global);
1625: VecResetArray(local);
1626: VecResetArray(global);
1627: DMRestoreGlobalVector(next->dm,&global);
1628: DMRestoreLocalVector(next->dm,&local);
1630: garray += next->n;
1631: larray += next->nlocal;
1632: next = next->next;
1633: }
1635: VecRestoreArray(gvec,NULL);
1636: VecRestoreArray(lvec,NULL);
1638: return 0;
1639: }
1641: PetscErrorCode DMLocalToGlobalEnd_Composite(DM dm,Vec lvec,InsertMode mode,Vec gvec)
1642: {
1646: return 0;
1647: }
1649: PetscErrorCode DMLocalToLocalBegin_Composite(DM dm,Vec vec1,InsertMode mode,Vec vec2)
1650: {
1651: struct DMCompositeLink *next;
1652: PetscScalar *array1,*array2;
1653: DM_Composite *com = (DM_Composite*)dm->data;
1659: if (!com->setup) {
1660: DMSetUp(dm);
1661: }
1663: VecGetArray(vec1,&array1);
1664: VecGetArray(vec2,&array2);
1666: /* loop over packed objects, handling one at at time */
1667: next = com->next;
1668: while (next) {
1669: Vec local1,local2;
1671: DMGetLocalVector(next->dm,&local1);
1672: VecPlaceArray(local1,array1);
1673: DMGetLocalVector(next->dm,&local2);
1674: VecPlaceArray(local2,array2);
1675: DMLocalToLocalBegin(next->dm,local1,mode,local2);
1676: DMLocalToLocalEnd(next->dm,local1,mode,local2);
1677: VecResetArray(local2);
1678: DMRestoreLocalVector(next->dm,&local2);
1679: VecResetArray(local1);
1680: DMRestoreLocalVector(next->dm,&local1);
1682: array1 += next->nlocal;
1683: array2 += next->nlocal;
1684: next = next->next;
1685: }
1687: VecRestoreArray(vec1,NULL);
1688: VecRestoreArray(vec2,NULL);
1690: return 0;
1691: }
1693: PetscErrorCode DMLocalToLocalEnd_Composite(DM dm,Vec lvec,InsertMode mode,Vec gvec)
1694: {
1698: return 0;
1699: }
1701: /*MC
1702: DMCOMPOSITE = "composite" - A DM object that is used to manage data for a collection of DMs
1704: Level: intermediate
1706: .seealso: DMType, DM, DMDACreate(), DMCreate(), DMSetType(), DMCompositeCreate()
1707: M*/
1709: PETSC_EXTERN PetscErrorCode DMCreate_Composite(DM p)
1710: {
1711: DM_Composite *com;
1713: PetscNewLog(p,&com);
1714: p->data = com;
1715: com->n = 0;
1716: com->nghost = 0;
1717: com->next = NULL;
1718: com->nDM = 0;
1720: p->ops->createglobalvector = DMCreateGlobalVector_Composite;
1721: p->ops->createlocalvector = DMCreateLocalVector_Composite;
1722: p->ops->getlocaltoglobalmapping = DMGetLocalToGlobalMapping_Composite;
1723: p->ops->createfieldis = DMCreateFieldIS_Composite;
1724: p->ops->createfielddecomposition = DMCreateFieldDecomposition_Composite;
1725: p->ops->refine = DMRefine_Composite;
1726: p->ops->coarsen = DMCoarsen_Composite;
1727: p->ops->createinterpolation = DMCreateInterpolation_Composite;
1728: p->ops->creatematrix = DMCreateMatrix_Composite;
1729: p->ops->getcoloring = DMCreateColoring_Composite;
1730: p->ops->globaltolocalbegin = DMGlobalToLocalBegin_Composite;
1731: p->ops->globaltolocalend = DMGlobalToLocalEnd_Composite;
1732: p->ops->localtoglobalbegin = DMLocalToGlobalBegin_Composite;
1733: p->ops->localtoglobalend = DMLocalToGlobalEnd_Composite;
1734: p->ops->localtolocalbegin = DMLocalToLocalBegin_Composite;
1735: p->ops->localtolocalend = DMLocalToLocalEnd_Composite;
1736: p->ops->destroy = DMDestroy_Composite;
1737: p->ops->view = DMView_Composite;
1738: p->ops->setup = DMSetUp_Composite;
1740: PetscObjectComposeFunction((PetscObject)p,"DMSetUpGLVisViewer_C",DMSetUpGLVisViewer_Composite);
1741: return 0;
1742: }
1744: /*@
1745: DMCompositeCreate - Creates a vector packer, used to generate "composite"
1746: vectors made up of several subvectors.
1748: Collective
1750: Input Parameter:
1751: . comm - the processors that will share the global vector
1753: Output Parameters:
1754: . packer - the packer object
1756: Level: advanced
1758: .seealso DMDestroy(), DMCompositeAddDM(), DMCompositeScatter(), DMCOMPOSITE,DMCreate()
1759: DMCompositeGather(), DMCreateGlobalVector(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess()
1760: DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(), DMCompositeGetEntries()
1762: @*/
1763: PetscErrorCode DMCompositeCreate(MPI_Comm comm,DM *packer)
1764: {
1766: DMCreate(comm,packer);
1767: DMSetType(*packer,DMCOMPOSITE);
1768: return 0;
1769: }