Actual source code: pack.c
petsc-3.9.4 2018-09-11
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.
12: Logically Collective on MPI_Comm
14: Input Parameter:
15: + dm - the composite object
16: - formcouplelocations - routine to set the nonzero locations in the matrix
18: Level: advanced
20: Not available from Fortran
22: Notes: 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;
33: PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
34: if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
35: com->FormCoupleLocations = FormCoupleLocations;
36: return(0);
37: }
39: PetscErrorCode DMDestroy_Composite(DM dm)
40: {
41: PetscErrorCode ierr;
42: struct DMCompositeLink *next, *prev;
43: DM_Composite *com = (DM_Composite*)dm->data;
46: next = com->next;
47: while (next) {
48: prev = next;
49: next = next->next;
50: DMDestroy(&prev->dm);
51: PetscFree(prev->grstarts);
52: PetscFree(prev);
53: }
54: PetscObjectComposeFunction((PetscObject)dm,"DMSetUpGLVisViewer_C",NULL);
55: /* This was originally freed in DMDestroy(), but that prevents reference counting of backend objects */
56: PetscFree(com);
57: return(0);
58: }
60: PetscErrorCode DMView_Composite(DM dm,PetscViewer v)
61: {
63: PetscBool iascii;
64: DM_Composite *com = (DM_Composite*)dm->data;
67: PetscObjectTypeCompare((PetscObject)v,PETSCVIEWERASCII,&iascii);
68: if (iascii) {
69: struct DMCompositeLink *lnk = com->next;
70: PetscInt i;
72: PetscViewerASCIIPrintf(v,"DM (%s)\n",((PetscObject)dm)->prefix ? ((PetscObject)dm)->prefix : "no prefix");
73: PetscViewerASCIIPrintf(v," contains %D DMs\n",com->nDM);
74: PetscViewerASCIIPushTab(v);
75: for (i=0; lnk; lnk=lnk->next,i++) {
76: PetscViewerASCIIPrintf(v,"Link %D: DM of type %s\n",i,((PetscObject)lnk->dm)->type_name);
77: PetscViewerASCIIPushTab(v);
78: DMView(lnk->dm,v);
79: PetscViewerASCIIPopTab(v);
80: }
81: PetscViewerASCIIPopTab(v);
82: }
83: return(0);
84: }
86: /* --------------------------------------------------------------------------------------*/
87: PetscErrorCode DMSetUp_Composite(DM dm)
88: {
89: PetscErrorCode ierr;
90: PetscInt nprev = 0;
91: PetscMPIInt rank,size;
92: DM_Composite *com = (DM_Composite*)dm->data;
93: struct DMCompositeLink *next = com->next;
94: PetscLayout map;
97: if (com->setup) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Packer has already been setup");
98: PetscLayoutCreate(PetscObjectComm((PetscObject)dm),&map);
99: PetscLayoutSetLocalSize(map,com->n);
100: PetscLayoutSetSize(map,PETSC_DETERMINE);
101: PetscLayoutSetBlockSize(map,1);
102: PetscLayoutSetUp(map);
103: PetscLayoutGetSize(map,&com->N);
104: PetscLayoutGetRange(map,&com->rstart,NULL);
105: PetscLayoutDestroy(&map);
107: /* now set the rstart for each linked vector */
108: MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);
109: MPI_Comm_size(PetscObjectComm((PetscObject)dm),&size);
110: while (next) {
111: next->rstart = nprev;
112: nprev += next->n;
113: next->grstart = com->rstart + next->rstart;
114: PetscMalloc1(size,&next->grstarts);
115: MPI_Allgather(&next->grstart,1,MPIU_INT,next->grstarts,1,MPIU_INT,PetscObjectComm((PetscObject)dm));
116: next = next->next;
117: }
118: com->setup = PETSC_TRUE;
119: return(0);
120: }
122: /* ----------------------------------------------------------------------------------*/
124: /*@
125: DMCompositeGetNumberDM - Get's the number of DM objects in the DMComposite
126: representation.
128: Not Collective
130: Input Parameter:
131: . dm - the packer object
133: Output Parameter:
134: . nDM - the number of DMs
136: Level: beginner
138: @*/
139: PetscErrorCode DMCompositeGetNumberDM(DM dm,PetscInt *nDM)
140: {
141: DM_Composite *com = (DM_Composite*)dm->data;
142: PetscBool flg;
147: PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
148: if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
149: *nDM = com->nDM;
150: return(0);
151: }
154: /*@C
155: DMCompositeGetAccess - Allows one to access the individual packed vectors in their global
156: representation.
158: Collective on DMComposite
160: Input Parameters:
161: + dm - the packer object
162: - gvec - the global vector
164: Output Parameters:
165: . Vec* ... - the packed parallel vectors, NULL for those that are not needed
167: Notes: Use DMCompositeRestoreAccess() to return the vectors when you no longer need them
169: Fortran Notes:
171: Fortran callers must use numbered versions of this routine, e.g., DMCompositeGetAccess4(dm,gvec,vec1,vec2,vec3,vec4)
172: or use the alternative interface DMCompositeGetAccessArray().
174: Level: advanced
176: .seealso: DMCompositeGetEntries(), DMCompositeScatter()
177: @*/
178: PetscErrorCode DMCompositeGetAccess(DM dm,Vec gvec,...)
179: {
180: va_list Argp;
181: PetscErrorCode ierr;
182: struct DMCompositeLink *next;
183: DM_Composite *com = (DM_Composite*)dm->data;
184: PetscInt readonly;
185: PetscBool flg;
190: PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
191: if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
192: next = com->next;
193: if (!com->setup) {
194: DMSetUp(dm);
195: }
197: VecLockGet(gvec,&readonly);
198: /* loop over packed objects, handling one at at time */
199: va_start(Argp,gvec);
200: while (next) {
201: Vec *vec;
202: vec = va_arg(Argp, Vec*);
203: if (vec) {
204: DMGetGlobalVector(next->dm,vec);
205: if (readonly) {
206: const PetscScalar *array;
207: VecGetArrayRead(gvec,&array);
208: VecPlaceArray(*vec,array+next->rstart);
209: VecLockPush(*vec);
210: VecRestoreArrayRead(gvec,&array);
211: } else {
212: PetscScalar *array;
213: VecGetArray(gvec,&array);
214: VecPlaceArray(*vec,array+next->rstart);
215: VecRestoreArray(gvec,&array);
216: }
217: }
218: next = next->next;
219: }
220: va_end(Argp);
221: return(0);
222: }
224: /*@C
225: DMCompositeGetAccessArray - Allows one to access the individual packed vectors in their global
226: representation.
228: Collective on DMComposite
230: Input Parameters:
231: + dm - the packer object
232: . pvec - packed vector
233: . nwanted - number of vectors wanted
234: - wanted - sorted array of vectors wanted, or NULL to get all vectors
236: Output Parameters:
237: . vecs - array of requested global vectors (must be allocated)
239: Notes: Use DMCompositeRestoreAccessArray() to return the vectors when you no longer need them
241: Level: advanced
243: .seealso: DMCompositeGetAccess(), DMCompositeGetEntries(), DMCompositeScatter(), DMCompositeGather()
244: @*/
245: PetscErrorCode DMCompositeGetAccessArray(DM dm,Vec pvec,PetscInt nwanted,const PetscInt *wanted,Vec *vecs)
246: {
247: PetscErrorCode ierr;
248: struct DMCompositeLink *link;
249: PetscInt i,wnum;
250: DM_Composite *com = (DM_Composite*)dm->data;
251: PetscInt readonly;
252: PetscBool flg;
257: PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
258: if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
259: if (!com->setup) {
260: DMSetUp(dm);
261: }
263: 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: DMGetGlobalVector(link->dm,&v);
268: if (readonly) {
269: const PetscScalar *array;
270: VecGetArrayRead(pvec,&array);
271: VecPlaceArray(v,array+link->rstart);
272: VecLockPush(v);
273: VecRestoreArrayRead(pvec,&array);
274: } else {
275: PetscScalar *array;
276: VecGetArray(pvec,&array);
277: VecPlaceArray(v,array+link->rstart);
278: VecRestoreArray(pvec,&array);
279: }
280: vecs[wnum++] = v;
281: }
282: }
283: return(0);
284: }
286: /*@C
287: DMCompositeGetLocalAccessArray - Allows one to access the individual
288: packed vectors in their local representation.
290: Collective on DMComposite.
292: Input Parameters:
293: + dm - the packer object
294: . pvec - packed vector
295: . nwanted - number of vectors wanted
296: - wanted - sorted array of vectors wanted, or NULL to get all vectors
298: Output Parameters:
299: . vecs - array of requested local vectors (must be allocated)
301: Notes: Use DMCompositeRestoreLocalAccessArray() to return the vectors
302: when you no longer need them.
304: Level: advanced
306: .seealso: DMCompositeRestoreLocalAccessArray(), DMCompositeGetAccess(),
307: DMCompositeGetEntries(), DMCompositeScatter(), DMCompositeGather()
308: @*/
309: PetscErrorCode DMCompositeGetLocalAccessArray(DM dm,Vec pvec,PetscInt nwanted,const PetscInt *wanted,Vec *vecs)
310: {
311: PetscErrorCode ierr;
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;
322: PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
323: if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
324: if (!com->setup) {
325: DMSetUp(dm);
326: }
328: VecLockGet(pvec,&readonly);
329: for (i=0,wnum=0,link=com->next; link && wnum<nwanted; i++,link=link->next) {
330: if (!wanted || i == wanted[wnum]) {
331: Vec v;
332: DMGetLocalVector(link->dm,&v);
333: if (readonly) {
334: const PetscScalar *array;
335: VecGetArrayRead(pvec,&array);
336: VecPlaceArray(v,array+nlocal);
337: VecLockPush(v);
338: VecRestoreArrayRead(pvec,&array);
339: } else {
340: PetscScalar *array;
341: VecGetArray(pvec,&array);
342: VecPlaceArray(v,array+nlocal);
343: VecRestoreArray(pvec,&array);
344: }
345: vecs[wnum++] = v;
346: }
348: nlocal += link->nlocal;
349: }
351: return(0);
352: }
354: /*@C
355: DMCompositeRestoreAccess - Returns the vectors obtained with DMCompositeGetAccess()
356: representation.
358: Collective on DMComposite
360: Input Parameters:
361: + dm - the packer object
362: . gvec - the global vector
363: - Vec* ... - the individual parallel vectors, NULL for those that are not needed
365: Level: advanced
367: .seealso DMCompositeAddDM(), DMCreateGlobalVector(),
368: DMCompositeGather(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeScatter(),
369: DMCompositeRestoreAccess(), DMCompositeGetAccess()
371: @*/
372: PetscErrorCode DMCompositeRestoreAccess(DM dm,Vec gvec,...)
373: {
374: va_list Argp;
375: PetscErrorCode ierr;
376: struct DMCompositeLink *next;
377: DM_Composite *com = (DM_Composite*)dm->data;
378: PetscInt readonly;
379: PetscBool flg;
384: PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
385: if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
386: next = com->next;
387: if (!com->setup) {
388: DMSetUp(dm);
389: }
391: VecLockGet(gvec,&readonly);
392: /* loop over packed objects, handling one at at time */
393: va_start(Argp,gvec);
394: while (next) {
395: Vec *vec;
396: vec = va_arg(Argp, Vec*);
397: if (vec) {
398: VecResetArray(*vec);
399: if (readonly) {
400: VecLockPop(*vec);
401: }
402: DMRestoreGlobalVector(next->dm,vec);
403: }
404: next = next->next;
405: }
406: va_end(Argp);
407: return(0);
408: }
410: /*@C
411: DMCompositeRestoreAccessArray - Returns the vectors obtained with DMCompositeGetAccessArray()
413: Collective on DMComposite
415: Input Parameters:
416: + dm - the packer object
417: . pvec - packed vector
418: . nwanted - number of vectors wanted
419: . wanted - sorted array of vectors wanted, or NULL to get all vectors
420: - vecs - array of global vectors to return
422: Level: advanced
424: .seealso: DMCompositeRestoreAccess(), DMCompositeRestoreEntries(), DMCompositeScatter(), DMCompositeGather()
425: @*/
426: PetscErrorCode DMCompositeRestoreAccessArray(DM dm,Vec pvec,PetscInt nwanted,const PetscInt *wanted,Vec *vecs)
427: {
428: PetscErrorCode ierr;
429: struct DMCompositeLink *link;
430: PetscInt i,wnum;
431: DM_Composite *com = (DM_Composite*)dm->data;
432: PetscInt readonly;
433: PetscBool flg;
438: PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
439: if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
440: if (!com->setup) {
441: DMSetUp(dm);
442: }
444: VecLockGet(pvec,&readonly);
445: for (i=0,wnum=0,link=com->next; link && wnum<nwanted; i++,link=link->next) {
446: if (!wanted || i == wanted[wnum]) {
447: VecResetArray(vecs[wnum]);
448: if (readonly) {
449: VecLockPop(vecs[wnum]);
450: }
451: DMRestoreGlobalVector(link->dm,&vecs[wnum]);
452: wnum++;
453: }
454: }
455: return(0);
456: }
458: /*@C
459: DMCompositeRestoreLocalAccessArray - Returns the vectors obtained with DMCompositeGetLocalAccessArray().
461: Collective on DMComposite.
463: Input Parameters:
464: + dm - the packer object
465: . pvec - packed vector
466: . nwanted - number of vectors wanted
467: . wanted - sorted array of vectors wanted, or NULL to restore all vectors
468: - vecs - array of local vectors to return
470: Level: advanced
472: Notes:
473: nwanted and wanted must match the values given to DMCompositeGetLocalAccessArray()
474: otherwise the call will fail.
476: .seealso: DMCompositeGetLocalAccessArray(), DMCompositeRestoreAccessArray(),
477: DMCompositeRestoreAccess(), DMCompositeRestoreEntries(),
478: DMCompositeScatter(), DMCompositeGather()
479: @*/
480: PetscErrorCode DMCompositeRestoreLocalAccessArray(DM dm,Vec pvec,PetscInt nwanted,const PetscInt *wanted,Vec *vecs)
481: {
482: PetscErrorCode ierr;
483: struct DMCompositeLink *link;
484: PetscInt i,wnum;
485: DM_Composite *com = (DM_Composite*)dm->data;
486: PetscInt readonly;
487: PetscBool flg;
492: PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
493: if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
494: if (!com->setup) {
495: DMSetUp(dm);
496: }
498: VecLockGet(pvec,&readonly);
499: for (i=0,wnum=0,link=com->next; link && wnum<nwanted; i++,link=link->next) {
500: if (!wanted || i == wanted[wnum]) {
501: VecResetArray(vecs[wnum]);
502: if (readonly) {
503: VecLockPop(vecs[wnum]);
504: }
505: DMRestoreLocalVector(link->dm,&vecs[wnum]);
506: wnum++;
507: }
508: }
509: return(0);
510: }
512: /*@C
513: DMCompositeScatter - Scatters from a global packed vector into its individual local vectors
515: Collective on DMComposite
517: Input Parameters:
518: + dm - the packer object
519: . gvec - the global vector
520: - Vec ... - the individual sequential vectors, NULL for those that are not needed
522: Level: advanced
524: Notes:
525: DMCompositeScatterArray() is a non-variadic alternative that is often more convenient for library callers and is
526: accessible from Fortran.
528: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(),
529: DMCompositeGather(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
530: DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(), DMCompositeGetEntries()
531: DMCompositeScatterArray()
533: @*/
534: PetscErrorCode DMCompositeScatter(DM dm,Vec gvec,...)
535: {
536: va_list Argp;
537: PetscErrorCode ierr;
538: struct DMCompositeLink *next;
539: PetscInt cnt;
540: DM_Composite *com = (DM_Composite*)dm->data;
541: PetscBool flg;
546: PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
547: if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
548: if (!com->setup) {
549: DMSetUp(dm);
550: }
552: /* loop over packed objects, handling one at at time */
553: va_start(Argp,gvec);
554: for (cnt=3,next=com->next; next; cnt++,next=next->next) {
555: Vec local;
556: local = va_arg(Argp, Vec);
557: if (local) {
558: Vec global;
559: const PetscScalar *array;
561: DMGetGlobalVector(next->dm,&global);
562: VecGetArrayRead(gvec,&array);
563: VecPlaceArray(global,array+next->rstart);
564: DMGlobalToLocalBegin(next->dm,global,INSERT_VALUES,local);
565: DMGlobalToLocalEnd(next->dm,global,INSERT_VALUES,local);
566: VecRestoreArrayRead(gvec,&array);
567: VecResetArray(global);
568: DMRestoreGlobalVector(next->dm,&global);
569: }
570: }
571: va_end(Argp);
572: return(0);
573: }
575: /*@
576: DMCompositeScatterArray - Scatters from a global packed vector into its individual local vectors
578: Collective on DMComposite
580: Input Parameters:
581: + dm - the packer object
582: . gvec - the global vector
583: . lvecs - array of local vectors, NULL for any that are not needed
585: Level: advanced
587: Note:
588: This is a non-variadic alternative to DMCompositeScatter()
590: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector()
591: DMCompositeGather(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
592: DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(), DMCompositeGetEntries()
594: @*/
595: PetscErrorCode DMCompositeScatterArray(DM dm,Vec gvec,Vec *lvecs)
596: {
597: PetscErrorCode ierr;
598: struct DMCompositeLink *next;
599: PetscInt i;
600: DM_Composite *com = (DM_Composite*)dm->data;
601: PetscBool flg;
606: PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
607: if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
608: if (!com->setup) {
609: DMSetUp(dm);
610: }
612: /* loop over packed objects, handling one at at time */
613: for (i=0,next=com->next; next; next=next->next,i++) {
614: if (lvecs[i]) {
615: Vec global;
616: const PetscScalar *array;
618: DMGetGlobalVector(next->dm,&global);
619: VecGetArrayRead(gvec,&array);
620: VecPlaceArray(global,(PetscScalar*)array+next->rstart);
621: DMGlobalToLocalBegin(next->dm,global,INSERT_VALUES,lvecs[i]);
622: DMGlobalToLocalEnd(next->dm,global,INSERT_VALUES,lvecs[i]);
623: VecRestoreArrayRead(gvec,&array);
624: VecResetArray(global);
625: DMRestoreGlobalVector(next->dm,&global);
626: }
627: }
628: return(0);
629: }
631: /*@C
632: DMCompositeGather - Gathers into a global packed vector from its individual local vectors
634: Collective on DMComposite
636: Input Parameter:
637: + dm - the packer object
638: . gvec - the global vector
639: . imode - INSERT_VALUES or ADD_VALUES
640: - Vec ... - the individual sequential vectors, NULL for any that are not needed
642: Level: advanced
644: Not available from Fortran, Fortran users can use DMCompositeGatherArray()
646: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(),
647: DMCompositeScatter(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
648: DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(), DMCompositeGetEntries()
650: @*/
651: PetscErrorCode DMCompositeGather(DM dm,InsertMode imode,Vec gvec,...)
652: {
653: va_list Argp;
654: PetscErrorCode ierr;
655: struct DMCompositeLink *next;
656: DM_Composite *com = (DM_Composite*)dm->data;
657: PetscInt cnt;
658: PetscBool flg;
663: PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
664: if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
665: if (!com->setup) {
666: DMSetUp(dm);
667: }
669: /* loop over packed objects, handling one at at time */
670: va_start(Argp,gvec);
671: for (cnt=3,next=com->next; next; cnt++,next=next->next) {
672: Vec local;
673: local = va_arg(Argp, Vec);
674: if (local) {
675: PetscScalar *array;
676: Vec global;
678: DMGetGlobalVector(next->dm,&global);
679: VecGetArray(gvec,&array);
680: VecPlaceArray(global,array+next->rstart);
681: DMLocalToGlobalBegin(next->dm,local,imode,global);
682: DMLocalToGlobalEnd(next->dm,local,imode,global);
683: VecRestoreArray(gvec,&array);
684: VecResetArray(global);
685: DMRestoreGlobalVector(next->dm,&global);
686: }
687: }
688: va_end(Argp);
689: return(0);
690: }
692: /*@
693: DMCompositeGatherArray - Gathers into a global packed vector from its individual local vectors
695: Collective on DMComposite
697: Input Parameter:
698: + dm - the packer object
699: . gvec - the global vector
700: . imode - INSERT_VALUES or ADD_VALUES
701: - lvecs - the individual sequential vectors, NULL for any that are not needed
703: Level: advanced
705: Notes:
706: This is a non-variadic alternative to DMCompositeGather().
708: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(),
709: DMCompositeScatter(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
710: DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(), DMCompositeGetEntries(),
711: @*/
712: PetscErrorCode DMCompositeGatherArray(DM dm,InsertMode imode,Vec gvec,Vec *lvecs)
713: {
714: PetscErrorCode ierr;
715: struct DMCompositeLink *next;
716: DM_Composite *com = (DM_Composite*)dm->data;
717: PetscInt i;
718: PetscBool flg;
723: PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
724: if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
725: if (!com->setup) {
726: DMSetUp(dm);
727: }
729: /* loop over packed objects, handling one at at time */
730: for (next=com->next,i=0; next; next=next->next,i++) {
731: if (lvecs[i]) {
732: PetscScalar *array;
733: Vec global;
735: DMGetGlobalVector(next->dm,&global);
736: VecGetArray(gvec,&array);
737: VecPlaceArray(global,array+next->rstart);
738: DMLocalToGlobalBegin(next->dm,lvecs[i],imode,global);
739: DMLocalToGlobalEnd(next->dm,lvecs[i],imode,global);
740: VecRestoreArray(gvec,&array);
741: VecResetArray(global);
742: DMRestoreGlobalVector(next->dm,&global);
743: }
744: }
745: return(0);
746: }
748: /*@
749: DMCompositeAddDM - adds a DM vector to a DMComposite
751: Collective on DMComposite
753: Input Parameter:
754: + dmc - the DMComposite (packer) object
755: - dm - the DM object
757: Level: advanced
759: .seealso DMDestroy(), DMCompositeGather(), DMCompositeAddDM(), DMCreateGlobalVector(),
760: DMCompositeScatter(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
761: DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(), DMCompositeGetEntries()
763: @*/
764: PetscErrorCode DMCompositeAddDM(DM dmc,DM dm)
765: {
766: PetscErrorCode ierr;
767: PetscInt n,nlocal;
768: struct DMCompositeLink *mine,*next;
769: Vec global,local;
770: DM_Composite *com = (DM_Composite*)dmc->data;
771: PetscBool flg;
776: PetscObjectTypeCompare((PetscObject)dmc,DMCOMPOSITE,&flg);
777: if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
778: next = com->next;
779: if (com->setup) SETERRQ(PetscObjectComm((PetscObject)dmc),PETSC_ERR_ARG_WRONGSTATE,"Cannot add a DM once you have used the DMComposite");
781: /* create new link */
782: PetscNew(&mine);
783: PetscObjectReference((PetscObject)dm);
784: DMGetGlobalVector(dm,&global);
785: VecGetLocalSize(global,&n);
786: DMRestoreGlobalVector(dm,&global);
787: DMGetLocalVector(dm,&local);
788: VecGetSize(local,&nlocal);
789: DMRestoreLocalVector(dm,&local);
791: mine->n = n;
792: mine->nlocal = nlocal;
793: mine->dm = dm;
794: mine->next = NULL;
795: com->n += n;
796: com->nghost += nlocal;
798: /* add to end of list */
799: if (!next) com->next = mine;
800: else {
801: while (next->next) next = next->next;
802: next->next = mine;
803: }
804: com->nDM++;
805: com->nmine++;
806: return(0);
807: }
809: #include <petscdraw.h>
810: PETSC_EXTERN PetscErrorCode VecView_MPI(Vec,PetscViewer);
811: PetscErrorCode VecView_DMComposite(Vec gvec,PetscViewer viewer)
812: {
813: DM dm;
814: PetscErrorCode ierr;
815: struct DMCompositeLink *next;
816: PetscBool isdraw;
817: DM_Composite *com;
820: VecGetDM(gvec, &dm);
821: if (!dm) SETERRQ(PetscObjectComm((PetscObject)gvec),PETSC_ERR_ARG_WRONG,"Vector not generated from a DMComposite");
822: com = (DM_Composite*)dm->data;
823: next = com->next;
825: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
826: if (!isdraw) {
827: /* do I really want to call this? */
828: VecView_MPI(gvec,viewer);
829: } else {
830: PetscInt cnt = 0;
832: /* loop over packed objects, handling one at at time */
833: while (next) {
834: Vec vec;
835: const PetscScalar *array;
836: PetscInt bs;
838: /* Should use VecGetSubVector() eventually, but would need to forward the DM for that to work */
839: DMGetGlobalVector(next->dm,&vec);
840: VecGetArrayRead(gvec,&array);
841: VecPlaceArray(vec,(PetscScalar*)array+next->rstart);
842: VecRestoreArrayRead(gvec,&array);
843: VecView(vec,viewer);
844: VecResetArray(vec);
845: VecGetBlockSize(vec,&bs);
846: DMRestoreGlobalVector(next->dm,&vec);
847: PetscViewerDrawBaseAdd(viewer,bs);
848: cnt += bs;
849: next = next->next;
850: }
851: PetscViewerDrawBaseAdd(viewer,-cnt);
852: }
853: return(0);
854: }
856: PetscErrorCode DMCreateGlobalVector_Composite(DM dm,Vec *gvec)
857: {
859: DM_Composite *com = (DM_Composite*)dm->data;
863: DMSetUp(dm);
864: VecCreateMPI(PetscObjectComm((PetscObject)dm),com->n,com->N,gvec);
865: VecSetDM(*gvec, dm);
866: VecSetOperation(*gvec,VECOP_VIEW,(void (*)(void))VecView_DMComposite);
867: return(0);
868: }
870: PetscErrorCode DMCreateLocalVector_Composite(DM dm,Vec *lvec)
871: {
873: DM_Composite *com = (DM_Composite*)dm->data;
877: if (!com->setup) {
878: DMSetUp(dm);
879: }
880: VecCreateSeq(PETSC_COMM_SELF,com->nghost,lvec);
881: VecSetDM(*lvec, dm);
882: return(0);
883: }
885: /*@C
886: DMCompositeGetISLocalToGlobalMappings - gets an ISLocalToGlobalMapping for each DM in the DMComposite, maps to the composite global space
888: Collective on DM
890: Input Parameter:
891: . dm - the packer object
893: Output Parameters:
894: . ltogs - the individual mappings for each packed vector. Note that this includes
895: all the ghost points that individual ghosted DMDA's may have.
897: Level: advanced
899: Notes:
900: Each entry of ltogs should be destroyed with ISLocalToGlobalMappingDestroy(), the ltogs array should be freed with PetscFree().
902: Not available from Fortran
904: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(),
905: DMCompositeGather(), DMCompositeCreate(), DMCompositeGetAccess(), DMCompositeScatter(),
906: DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(),DMCompositeGetEntries()
908: @*/
909: PetscErrorCode DMCompositeGetISLocalToGlobalMappings(DM dm,ISLocalToGlobalMapping **ltogs)
910: {
911: PetscErrorCode ierr;
912: PetscInt i,*idx,n,cnt;
913: struct DMCompositeLink *next;
914: PetscMPIInt rank;
915: DM_Composite *com = (DM_Composite*)dm->data;
916: PetscBool flg;
920: PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
921: if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
922: DMSetUp(dm);
923: PetscMalloc1(com->nDM,ltogs);
924: next = com->next;
925: MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);
927: /* loop over packed objects, handling one at at time */
928: cnt = 0;
929: while (next) {
930: ISLocalToGlobalMapping ltog;
931: PetscMPIInt size;
932: const PetscInt *suboff,*indices;
933: Vec global;
935: /* Get sub-DM global indices for each local dof */
936: DMGetLocalToGlobalMapping(next->dm,<og);
937: ISLocalToGlobalMappingGetSize(ltog,&n);
938: ISLocalToGlobalMappingGetIndices(ltog,&indices);
939: PetscMalloc1(n,&idx);
941: /* Get the offsets for the sub-DM global vector */
942: DMGetGlobalVector(next->dm,&global);
943: VecGetOwnershipRanges(global,&suboff);
944: MPI_Comm_size(PetscObjectComm((PetscObject)global),&size);
946: /* Shift the sub-DM definition of the global space to the composite global space */
947: for (i=0; i<n; i++) {
948: PetscInt subi = indices[i],lo = 0,hi = size,t;
949: /* Binary search to find which rank owns subi */
950: while (hi-lo > 1) {
951: t = lo + (hi-lo)/2;
952: if (suboff[t] > subi) hi = t;
953: else lo = t;
954: }
955: idx[i] = subi - suboff[lo] + next->grstarts[lo];
956: }
957: ISLocalToGlobalMappingRestoreIndices(ltog,&indices);
958: ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)dm),1,n,idx,PETSC_OWN_POINTER,&(*ltogs)[cnt]);
959: DMRestoreGlobalVector(next->dm,&global);
960: next = next->next;
961: cnt++;
962: }
963: return(0);
964: }
966: /*@C
967: DMCompositeGetLocalISs - Gets index sets for each component of a composite local vector
969: Not Collective
971: Input Arguments:
972: . dm - composite DM
974: Output Arguments:
975: . is - array of serial index sets for each each component of the DMComposite
977: Level: intermediate
979: Notes:
980: At present, a composite local vector does not normally exist. This function is used to provide index sets for
981: MatGetLocalSubMatrix(). In the future, the scatters for each entry in the DMComposite may be be merged into a single
982: scatter to a composite local vector. The user should not typically need to know which is being done.
984: To get the composite global indices at all local points (including ghosts), use DMCompositeGetISLocalToGlobalMappings().
986: To get index sets for pieces of the composite global vector, use DMCompositeGetGlobalISs().
988: Each returned IS should be destroyed with ISDestroy(), the array should be freed with PetscFree().
990: Not available from Fortran
992: .seealso: DMCompositeGetGlobalISs(), DMCompositeGetISLocalToGlobalMappings(), MatGetLocalSubMatrix(), MatCreateLocalRef()
993: @*/
994: PetscErrorCode DMCompositeGetLocalISs(DM dm,IS **is)
995: {
996: PetscErrorCode ierr;
997: DM_Composite *com = (DM_Composite*)dm->data;
998: struct DMCompositeLink *link;
999: PetscInt cnt,start;
1000: PetscBool flg;
1005: PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
1006: if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
1007: PetscMalloc1(com->nmine,is);
1008: for (cnt=0,start=0,link=com->next; link; start+=link->nlocal,cnt++,link=link->next) {
1009: PetscInt bs;
1010: ISCreateStride(PETSC_COMM_SELF,link->nlocal,start,1,&(*is)[cnt]);
1011: DMGetBlockSize(link->dm,&bs);
1012: ISSetBlockSize((*is)[cnt],bs);
1013: }
1014: return(0);
1015: }
1017: /*@C
1018: DMCompositeGetGlobalISs - Gets the index sets for each composed object
1020: Collective on DMComposite
1022: Input Parameter:
1023: . dm - the packer object
1025: Output Parameters:
1026: . is - the array of index sets
1028: Level: advanced
1030: Notes:
1031: The is entries should be destroyed with ISDestroy(), the is array should be freed with PetscFree()
1033: These could be used to extract a subset of vector entries for a "multi-physics" preconditioner
1035: Use DMCompositeGetLocalISs() for index sets in the packed local numbering, and
1036: DMCompositeGetISLocalToGlobalMappings() for to map local sub-DM (including ghost) indices to packed global
1037: indices.
1039: Fortran Notes:
1041: The output argument 'is' must be an allocated array of sufficient length, which can be learned using DMCompositeGetNumberDM().
1043: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(),
1044: DMCompositeGather(), DMCompositeCreate(), DMCompositeGetAccess(), DMCompositeScatter(),
1045: DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(),DMCompositeGetEntries()
1047: @*/
1048: PetscErrorCode DMCompositeGetGlobalISs(DM dm,IS *is[])
1049: {
1050: PetscErrorCode ierr;
1051: PetscInt cnt = 0;
1052: struct DMCompositeLink *next;
1053: PetscMPIInt rank;
1054: DM_Composite *com = (DM_Composite*)dm->data;
1055: PetscBool flg;
1059: PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
1060: if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
1061: if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Must call DMSetUp() before");
1062: PetscMalloc1(com->nDM,is);
1063: next = com->next;
1064: MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);
1066: /* loop over packed objects, handling one at at time */
1067: while (next) {
1068: ISCreateStride(PetscObjectComm((PetscObject)dm),next->n,next->grstart,1,&(*is)[cnt]);
1069: if (dm->prob) {
1070: MatNullSpace space;
1071: Mat pmat;
1072: PetscObject disc;
1073: PetscInt Nf;
1075: PetscDSGetNumFields(dm->prob, &Nf);
1076: if (cnt < Nf) {
1077: PetscDSGetDiscretization(dm->prob, cnt, &disc);
1078: PetscObjectQuery(disc, "nullspace", (PetscObject*) &space);
1079: if (space) {PetscObjectCompose((PetscObject) (*is)[cnt], "nullspace", (PetscObject) space);}
1080: PetscObjectQuery(disc, "nearnullspace", (PetscObject*) &space);
1081: if (space) {PetscObjectCompose((PetscObject) (*is)[cnt], "nearnullspace", (PetscObject) space);}
1082: PetscObjectQuery(disc, "pmat", (PetscObject*) &pmat);
1083: if (pmat) {PetscObjectCompose((PetscObject) (*is)[cnt], "pmat", (PetscObject) pmat);}
1084: }
1085: }
1086: cnt++;
1087: next = next->next;
1088: }
1089: return(0);
1090: }
1092: PetscErrorCode DMCreateFieldIS_Composite(DM dm, PetscInt *numFields,char ***fieldNames, IS **fields)
1093: {
1094: PetscInt nDM;
1095: DM *dms;
1096: PetscInt i;
1100: DMCompositeGetNumberDM(dm, &nDM);
1101: if (numFields) *numFields = nDM;
1102: DMCompositeGetGlobalISs(dm, fields);
1103: if (fieldNames) {
1104: PetscMalloc1(nDM, &dms);
1105: PetscMalloc1(nDM, fieldNames);
1106: DMCompositeGetEntriesArray(dm, dms);
1107: for (i=0; i<nDM; i++) {
1108: char buf[256];
1109: const char *splitname;
1111: /* Split naming precedence: object name, prefix, number */
1112: splitname = ((PetscObject) dm)->name;
1113: if (!splitname) {
1114: PetscObjectGetOptionsPrefix((PetscObject)dms[i],&splitname);
1115: if (splitname) {
1116: size_t len;
1117: PetscStrncpy(buf,splitname,sizeof(buf));
1118: buf[sizeof(buf) - 1] = 0;
1119: PetscStrlen(buf,&len);
1120: if (buf[len-1] == '_') buf[len-1] = 0; /* Remove trailing underscore if it was used */
1121: splitname = buf;
1122: }
1123: }
1124: if (!splitname) {
1125: PetscSNPrintf(buf,sizeof(buf),"%D",i);
1126: splitname = buf;
1127: }
1128: PetscStrallocpy(splitname,&(*fieldNames)[i]);
1129: }
1130: PetscFree(dms);
1131: }
1132: return(0);
1133: }
1135: /*
1136: This could take over from DMCreateFieldIS(), as it is more general,
1137: making DMCreateFieldIS() a special case -- calling with dmlist == NULL;
1138: At this point it's probably best to be less intrusive, however.
1139: */
1140: PetscErrorCode DMCreateFieldDecomposition_Composite(DM dm, PetscInt *len,char ***namelist, IS **islist, DM **dmlist)
1141: {
1142: PetscInt nDM;
1143: PetscInt i;
1147: DMCreateFieldIS_Composite(dm, len, namelist, islist);
1148: if (dmlist) {
1149: DMCompositeGetNumberDM(dm, &nDM);
1150: PetscMalloc1(nDM, dmlist);
1151: DMCompositeGetEntriesArray(dm, *dmlist);
1152: for (i=0; i<nDM; i++) {
1153: PetscObjectReference((PetscObject)((*dmlist)[i]));
1154: }
1155: }
1156: return(0);
1157: }
1161: /* -------------------------------------------------------------------------------------*/
1162: /*@C
1163: DMCompositeGetLocalVectors - Gets local vectors for each part of a DMComposite.
1164: Use DMCompositeRestoreLocalVectors() to return them.
1166: Not Collective
1168: Input Parameter:
1169: . dm - the packer object
1171: Output Parameter:
1172: . Vec ... - the individual sequential Vecs
1174: Level: advanced
1176: Not available from Fortran
1178: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(),
1179: DMCompositeGather(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
1180: DMCompositeRestoreLocalVectors(), DMCompositeScatter(), DMCompositeGetEntries()
1182: @*/
1183: PetscErrorCode DMCompositeGetLocalVectors(DM dm,...)
1184: {
1185: va_list Argp;
1186: PetscErrorCode ierr;
1187: struct DMCompositeLink *next;
1188: DM_Composite *com = (DM_Composite*)dm->data;
1189: PetscBool flg;
1193: PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
1194: if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
1195: next = com->next;
1196: /* loop over packed objects, handling one at at time */
1197: va_start(Argp,dm);
1198: while (next) {
1199: Vec *vec;
1200: vec = va_arg(Argp, Vec*);
1201: if (vec) {DMGetLocalVector(next->dm,vec);}
1202: next = next->next;
1203: }
1204: va_end(Argp);
1205: return(0);
1206: }
1208: /*@C
1209: DMCompositeRestoreLocalVectors - Restores local vectors for each part of a DMComposite.
1211: Not Collective
1213: Input Parameter:
1214: . dm - the packer object
1216: Output Parameter:
1217: . Vec ... - the individual sequential Vecs
1219: Level: advanced
1221: Not available from Fortran
1223: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(),
1224: DMCompositeGather(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
1225: DMCompositeGetLocalVectors(), DMCompositeScatter(), DMCompositeGetEntries()
1227: @*/
1228: PetscErrorCode DMCompositeRestoreLocalVectors(DM dm,...)
1229: {
1230: va_list Argp;
1231: PetscErrorCode ierr;
1232: struct DMCompositeLink *next;
1233: DM_Composite *com = (DM_Composite*)dm->data;
1234: PetscBool flg;
1238: PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
1239: if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
1240: next = com->next;
1241: /* loop over packed objects, handling one at at time */
1242: va_start(Argp,dm);
1243: while (next) {
1244: Vec *vec;
1245: vec = va_arg(Argp, Vec*);
1246: if (vec) {DMRestoreLocalVector(next->dm,vec);}
1247: next = next->next;
1248: }
1249: va_end(Argp);
1250: return(0);
1251: }
1253: /* -------------------------------------------------------------------------------------*/
1254: /*@C
1255: DMCompositeGetEntries - Gets the DM for each entry in a DMComposite.
1257: Not Collective
1259: Input Parameter:
1260: . dm - the packer object
1262: Output Parameter:
1263: . DM ... - the individual entries (DMs)
1265: Level: advanced
1267: Fortran Notes: Available as DMCompositeGetEntries() for one output DM, DMCompositeGetEntries2() for 2, etc
1269: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(), DMCompositeGetEntriesArray()
1270: DMCompositeGather(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
1271: DMCompositeRestoreLocalVectors(), DMCompositeGetLocalVectors(), DMCompositeScatter(),
1272: DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors()
1274: @*/
1275: PetscErrorCode DMCompositeGetEntries(DM dm,...)
1276: {
1277: va_list Argp;
1278: struct DMCompositeLink *next;
1279: DM_Composite *com = (DM_Composite*)dm->data;
1280: PetscBool flg;
1281: PetscErrorCode ierr;
1285: PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
1286: if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
1287: next = com->next;
1288: /* loop over packed objects, handling one at at time */
1289: va_start(Argp,dm);
1290: while (next) {
1291: DM *dmn;
1292: dmn = va_arg(Argp, DM*);
1293: if (dmn) *dmn = next->dm;
1294: next = next->next;
1295: }
1296: va_end(Argp);
1297: return(0);
1298: }
1300: /*@C
1301: DMCompositeGetEntriesArray - Gets the DM for each entry in a DMComposite.
1303: Not Collective
1305: Input Parameter:
1306: . dm - the packer object
1308: Output Parameter:
1309: . dms - array of sufficient length (see DMCompositeGetNumberDM()) to hold the individual DMs
1311: Level: advanced
1313: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(), DMCompositeGetEntries()
1314: DMCompositeGather(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
1315: DMCompositeRestoreLocalVectors(), DMCompositeGetLocalVectors(), DMCompositeScatter(),
1316: DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors()
1318: @*/
1319: PetscErrorCode DMCompositeGetEntriesArray(DM dm,DM dms[])
1320: {
1321: struct DMCompositeLink *next;
1322: DM_Composite *com = (DM_Composite*)dm->data;
1323: PetscInt i;
1324: PetscBool flg;
1325: PetscErrorCode ierr;
1329: PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
1330: if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
1331: /* loop over packed objects, handling one at at time */
1332: for (next=com->next,i=0; next; next=next->next,i++) dms[i] = next->dm;
1333: return(0);
1334: }
1336: typedef struct {
1337: DM dm;
1338: PetscViewer *subv;
1339: Vec *vecs;
1340: } GLVisViewerCtx;
1342: static PetscErrorCode DestroyGLVisViewerCtx_Private(void *vctx)
1343: {
1344: GLVisViewerCtx *ctx = (GLVisViewerCtx*)vctx;
1345: PetscInt i,n;
1349: DMCompositeGetNumberDM(ctx->dm,&n);
1350: for (i = 0; i < n; i++) {
1351: PetscViewerDestroy(&ctx->subv[i]);
1352: }
1353: PetscFree2(ctx->subv,ctx->vecs);
1354: DMDestroy(&ctx->dm);
1355: PetscFree(ctx);
1356: return(0);
1357: }
1359: static PetscErrorCode DMCompositeSampleGLVisFields_Private(PetscObject oX, PetscInt nf, PetscObject oXfield[], void *vctx)
1360: {
1361: Vec X = (Vec)oX;
1362: GLVisViewerCtx *ctx = (GLVisViewerCtx*)vctx;
1363: PetscInt i,n,cumf;
1367: DMCompositeGetNumberDM(ctx->dm,&n);
1368: DMCompositeGetAccessArray(ctx->dm,X,n,NULL,ctx->vecs);
1369: for (i = 0, cumf = 0; i < n; i++) {
1370: PetscErrorCode (*g2l)(PetscObject,PetscInt,PetscObject[],void*);
1371: void *fctx;
1372: PetscInt nfi;
1374: PetscViewerGLVisGetFields_Private(ctx->subv[i],&nfi,NULL,NULL,&g2l,NULL,&fctx);
1375: if (!nfi) continue;
1376: if (g2l) {
1377: (*g2l)((PetscObject)ctx->vecs[i],nfi,oXfield+cumf,fctx);
1378: } else {
1379: VecCopy(ctx->vecs[i],(Vec)(oXfield[cumf]));
1380: }
1381: cumf += nfi;
1382: }
1383: DMCompositeRestoreAccessArray(ctx->dm,X,n,NULL,ctx->vecs);
1384: return(0);
1385: }
1387: static PetscErrorCode DMSetUpGLVisViewer_Composite(PetscObject odm, PetscViewer viewer)
1388: {
1389: DM dm = (DM)odm, *dms;
1390: Vec *Ufds;
1391: GLVisViewerCtx *ctx;
1392: PetscInt i,n,tnf,*sdim;
1393: char **fecs;
1397: PetscNew(&ctx);
1398: PetscObjectReference((PetscObject)dm);
1399: ctx->dm = dm;
1400: DMCompositeGetNumberDM(dm,&n);
1401: PetscMalloc1(n,&dms);
1402: DMCompositeGetEntriesArray(dm,dms);
1403: PetscMalloc2(n,&ctx->subv,n,&ctx->vecs);
1404: for (i = 0, tnf = 0; i < n; i++) {
1405: PetscInt nf;
1407: PetscViewerCreate(PetscObjectComm(odm),&ctx->subv[i]);
1408: PetscViewerSetType(ctx->subv[i],PETSCVIEWERGLVIS);
1409: PetscViewerGLVisSetDM_Private(ctx->subv[i],(PetscObject)dms[i]);
1410: PetscViewerGLVisGetFields_Private(ctx->subv[i],&nf,NULL,NULL,NULL,NULL,NULL);
1411: tnf += nf;
1412: }
1413: PetscFree(dms);
1414: PetscMalloc3(tnf,&fecs,tnf,&sdim,tnf,&Ufds);
1415: for (i = 0, tnf = 0; i < n; i++) {
1416: PetscInt *sd,nf,f;
1417: const char **fec;
1418: Vec *Uf;
1420: PetscViewerGLVisGetFields_Private(ctx->subv[i],&nf,&fec,&sd,NULL,(PetscObject**)&Uf,NULL);
1421: for (f = 0; f < nf; f++) {
1422: PetscStrallocpy(fec[f],&fecs[tnf+f]);
1423: Ufds[tnf+f] = Uf[f];
1424: sdim[tnf+f] = sd[f];
1425: }
1426: tnf += nf;
1427: }
1428: PetscViewerGLVisSetFields(viewer,tnf,(const char**)fecs,sdim,DMCompositeSampleGLVisFields_Private,(PetscObject*)Ufds,ctx,DestroyGLVisViewerCtx_Private);
1429: for (i = 0; i < tnf; i++) {
1430: PetscFree(fecs[i]);
1431: }
1432: PetscFree3(fecs,sdim,Ufds);
1433: return(0);
1434: }
1436: PetscErrorCode DMRefine_Composite(DM dmi,MPI_Comm comm,DM *fine)
1437: {
1438: PetscErrorCode ierr;
1439: struct DMCompositeLink *next;
1440: DM_Composite *com = (DM_Composite*)dmi->data;
1441: DM dm;
1445: if (comm == MPI_COMM_NULL) {
1446: PetscObjectGetComm((PetscObject)dmi,&comm);
1447: }
1448: DMSetUp(dmi);
1449: next = com->next;
1450: DMCompositeCreate(comm,fine);
1452: /* loop over packed objects, handling one at at time */
1453: while (next) {
1454: DMRefine(next->dm,comm,&dm);
1455: DMCompositeAddDM(*fine,dm);
1456: PetscObjectDereference((PetscObject)dm);
1457: next = next->next;
1458: }
1459: return(0);
1460: }
1462: PetscErrorCode DMCoarsen_Composite(DM dmi,MPI_Comm comm,DM *fine)
1463: {
1464: PetscErrorCode ierr;
1465: struct DMCompositeLink *next;
1466: DM_Composite *com = (DM_Composite*)dmi->data;
1467: DM dm;
1471: DMSetUp(dmi);
1472: if (comm == MPI_COMM_NULL) {
1473: PetscObjectGetComm((PetscObject)dmi,&comm);
1474: }
1475: next = com->next;
1476: DMCompositeCreate(comm,fine);
1478: /* loop over packed objects, handling one at at time */
1479: while (next) {
1480: DMCoarsen(next->dm,comm,&dm);
1481: DMCompositeAddDM(*fine,dm);
1482: PetscObjectDereference((PetscObject)dm);
1483: next = next->next;
1484: }
1485: return(0);
1486: }
1488: PetscErrorCode DMCreateInterpolation_Composite(DM coarse,DM fine,Mat *A,Vec *v)
1489: {
1490: PetscErrorCode ierr;
1491: PetscInt m,n,M,N,nDM,i;
1492: struct DMCompositeLink *nextc;
1493: struct DMCompositeLink *nextf;
1494: Vec gcoarse,gfine,*vecs;
1495: DM_Composite *comcoarse = (DM_Composite*)coarse->data;
1496: DM_Composite *comfine = (DM_Composite*)fine->data;
1497: Mat *mats;
1502: DMSetUp(coarse);
1503: DMSetUp(fine);
1504: /* use global vectors only for determining matrix layout */
1505: DMGetGlobalVector(coarse,&gcoarse);
1506: DMGetGlobalVector(fine,&gfine);
1507: VecGetLocalSize(gcoarse,&n);
1508: VecGetLocalSize(gfine,&m);
1509: VecGetSize(gcoarse,&N);
1510: VecGetSize(gfine,&M);
1511: DMRestoreGlobalVector(coarse,&gcoarse);
1512: DMRestoreGlobalVector(fine,&gfine);
1514: nDM = comfine->nDM;
1515: if (nDM != comcoarse->nDM) SETERRQ2(PetscObjectComm((PetscObject)fine),PETSC_ERR_ARG_INCOMP,"Fine DMComposite has %D entries, but coarse has %D",nDM,comcoarse->nDM);
1516: PetscCalloc1(nDM*nDM,&mats);
1517: if (v) {
1518: PetscCalloc1(nDM,&vecs);
1519: }
1521: /* loop over packed objects, handling one at at time */
1522: for (nextc=comcoarse->next,nextf=comfine->next,i=0; nextc; nextc=nextc->next,nextf=nextf->next,i++) {
1523: if (!v) {
1524: DMCreateInterpolation(nextc->dm,nextf->dm,&mats[i*nDM+i],NULL);
1525: } else {
1526: DMCreateInterpolation(nextc->dm,nextf->dm,&mats[i*nDM+i],&vecs[i]);
1527: }
1528: }
1529: MatCreateNest(PetscObjectComm((PetscObject)fine),nDM,NULL,nDM,NULL,mats,A);
1530: if (v) {
1531: VecCreateNest(PetscObjectComm((PetscObject)fine),nDM,NULL,vecs,v);
1532: }
1533: for (i=0; i<nDM*nDM; i++) {MatDestroy(&mats[i]);}
1534: PetscFree(mats);
1535: if (v) {
1536: for (i=0; i<nDM; i++) {VecDestroy(&vecs[i]);}
1537: PetscFree(vecs);
1538: }
1539: return(0);
1540: }
1542: static PetscErrorCode DMGetLocalToGlobalMapping_Composite(DM dm)
1543: {
1544: DM_Composite *com = (DM_Composite*)dm->data;
1545: ISLocalToGlobalMapping *ltogs;
1546: PetscInt i;
1547: PetscErrorCode ierr;
1550: /* Set the ISLocalToGlobalMapping on the new matrix */
1551: DMCompositeGetISLocalToGlobalMappings(dm,<ogs);
1552: ISLocalToGlobalMappingConcatenate(PetscObjectComm((PetscObject)dm),com->nDM,ltogs,&dm->ltogmap);
1553: for (i=0; i<com->nDM; i++) {ISLocalToGlobalMappingDestroy(<ogs[i]);}
1554: PetscFree(ltogs);
1555: return(0);
1556: }
1559: PetscErrorCode DMCreateColoring_Composite(DM dm,ISColoringType ctype,ISColoring *coloring)
1560: {
1561: PetscErrorCode ierr;
1562: PetscInt n,i,cnt;
1563: ISColoringValue *colors;
1564: PetscBool dense = PETSC_FALSE;
1565: ISColoringValue maxcol = 0;
1566: DM_Composite *com = (DM_Composite*)dm->data;
1570: if (ctype == IS_COLORING_LOCAL) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Only global coloring supported");
1571: else if (ctype == IS_COLORING_GLOBAL) {
1572: n = com->n;
1573: } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unknown ISColoringType");
1574: PetscMalloc1(n,&colors); /* freed in ISColoringDestroy() */
1576: PetscOptionsGetBool(((PetscObject)dm)->options,((PetscObject)dm)->prefix,"-dmcomposite_dense_jacobian",&dense,NULL);
1577: if (dense) {
1578: for (i=0; i<n; i++) {
1579: colors[i] = (ISColoringValue)(com->rstart + i);
1580: }
1581: maxcol = com->N;
1582: } else {
1583: struct DMCompositeLink *next = com->next;
1584: PetscMPIInt rank;
1586: MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);
1587: cnt = 0;
1588: while (next) {
1589: ISColoring lcoloring;
1591: DMCreateColoring(next->dm,IS_COLORING_GLOBAL,&lcoloring);
1592: for (i=0; i<lcoloring->N; i++) {
1593: colors[cnt++] = maxcol + lcoloring->colors[i];
1594: }
1595: maxcol += lcoloring->n;
1596: ISColoringDestroy(&lcoloring);
1597: next = next->next;
1598: }
1599: }
1600: ISColoringCreate(PetscObjectComm((PetscObject)dm),maxcol,n,colors,PETSC_OWN_POINTER,coloring);
1601: return(0);
1602: }
1604: PetscErrorCode DMGlobalToLocalBegin_Composite(DM dm,Vec gvec,InsertMode mode,Vec lvec)
1605: {
1606: PetscErrorCode ierr;
1607: struct DMCompositeLink *next;
1608: PetscScalar *garray,*larray;
1609: DM_Composite *com = (DM_Composite*)dm->data;
1615: if (!com->setup) {
1616: DMSetUp(dm);
1617: }
1619: VecGetArray(gvec,&garray);
1620: VecGetArray(lvec,&larray);
1622: /* loop over packed objects, handling one at at time */
1623: next = com->next;
1624: while (next) {
1625: Vec local,global;
1626: PetscInt N;
1628: DMGetGlobalVector(next->dm,&global);
1629: VecGetLocalSize(global,&N);
1630: VecPlaceArray(global,garray);
1631: DMGetLocalVector(next->dm,&local);
1632: VecPlaceArray(local,larray);
1633: DMGlobalToLocalBegin(next->dm,global,mode,local);
1634: DMGlobalToLocalEnd(next->dm,global,mode,local);
1635: VecResetArray(global);
1636: VecResetArray(local);
1637: DMRestoreGlobalVector(next->dm,&global);
1638: DMRestoreLocalVector(next->dm,&local);
1640: larray += next->nlocal;
1641: garray += next->n;
1642: next = next->next;
1643: }
1645: VecRestoreArray(gvec,NULL);
1646: VecRestoreArray(lvec,NULL);
1647: return(0);
1648: }
1650: PetscErrorCode DMGlobalToLocalEnd_Composite(DM dm,Vec gvec,InsertMode mode,Vec lvec)
1651: {
1656: return(0);
1657: }
1659: PetscErrorCode DMLocalToGlobalBegin_Composite(DM dm,Vec lvec,InsertMode mode,Vec gvec)
1660: {
1661: PetscErrorCode ierr;
1662: struct DMCompositeLink *next;
1663: PetscScalar *larray,*garray;
1664: DM_Composite *com = (DM_Composite*)dm->data;
1671: if (!com->setup) {
1672: DMSetUp(dm);
1673: }
1675: VecGetArray(lvec,&larray);
1676: VecGetArray(gvec,&garray);
1678: /* loop over packed objects, handling one at at time */
1679: next = com->next;
1680: while (next) {
1681: Vec global,local;
1683: DMGetLocalVector(next->dm,&local);
1684: VecPlaceArray(local,larray);
1685: DMGetGlobalVector(next->dm,&global);
1686: VecPlaceArray(global,garray);
1687: DMLocalToGlobalBegin(next->dm,local,mode,global);
1688: DMLocalToGlobalEnd(next->dm,local,mode,global);
1689: VecResetArray(local);
1690: VecResetArray(global);
1691: DMRestoreGlobalVector(next->dm,&global);
1692: DMRestoreLocalVector(next->dm,&local);
1694: garray += next->n;
1695: larray += next->nlocal;
1696: next = next->next;
1697: }
1699: VecRestoreArray(gvec,NULL);
1700: VecRestoreArray(lvec,NULL);
1702: return(0);
1703: }
1705: PetscErrorCode DMLocalToGlobalEnd_Composite(DM dm,Vec lvec,InsertMode mode,Vec gvec)
1706: {
1711: return(0);
1712: }
1714: PetscErrorCode DMLocalToLocalBegin_Composite(DM dm,Vec vec1,InsertMode mode,Vec vec2)
1715: {
1716: PetscErrorCode ierr;
1717: struct DMCompositeLink *next;
1718: PetscScalar *array1,*array2;
1719: DM_Composite *com = (DM_Composite*)dm->data;
1726: if (!com->setup) {
1727: DMSetUp(dm);
1728: }
1730: VecGetArray(vec1,&array1);
1731: VecGetArray(vec2,&array2);
1733: /* loop over packed objects, handling one at at time */
1734: next = com->next;
1735: while (next) {
1736: Vec local1,local2;
1738: DMGetLocalVector(next->dm,&local1);
1739: VecPlaceArray(local1,array1);
1740: DMGetLocalVector(next->dm,&local2);
1741: VecPlaceArray(local2,array2);
1742: DMLocalToLocalBegin(next->dm,local1,mode,local2);
1743: DMLocalToLocalEnd(next->dm,local1,mode,local2);
1744: VecResetArray(local2);
1745: DMRestoreLocalVector(next->dm,&local2);
1746: VecResetArray(local1);
1747: DMRestoreLocalVector(next->dm,&local1);
1749: array1 += next->nlocal;
1750: array2 += next->nlocal;
1751: next = next->next;
1752: }
1754: VecRestoreArray(vec1,NULL);
1755: VecRestoreArray(vec2,NULL);
1757: return(0);
1758: }
1760: PetscErrorCode DMLocalToLocalEnd_Composite(DM dm,Vec lvec,InsertMode mode,Vec gvec)
1761: {
1766: return(0);
1767: }
1769: /*MC
1770: DMCOMPOSITE = "composite" - A DM object that is used to manage data for a collection of DMs
1772: Level: intermediate
1774: .seealso: DMType, DM, DMDACreate(), DMCreate(), DMSetType(), DMCompositeCreate()
1775: M*/
1778: PETSC_EXTERN PetscErrorCode DMCreate_Composite(DM p)
1779: {
1781: DM_Composite *com;
1784: PetscNewLog(p,&com);
1785: p->data = com;
1786: PetscObjectChangeTypeName((PetscObject)p,"DMComposite");
1787: com->n = 0;
1788: com->nghost = 0;
1789: com->next = NULL;
1790: com->nDM = 0;
1792: p->ops->createglobalvector = DMCreateGlobalVector_Composite;
1793: p->ops->createlocalvector = DMCreateLocalVector_Composite;
1794: p->ops->getlocaltoglobalmapping = DMGetLocalToGlobalMapping_Composite;
1795: p->ops->createfieldis = DMCreateFieldIS_Composite;
1796: p->ops->createfielddecomposition = DMCreateFieldDecomposition_Composite;
1797: p->ops->refine = DMRefine_Composite;
1798: p->ops->coarsen = DMCoarsen_Composite;
1799: p->ops->createinterpolation = DMCreateInterpolation_Composite;
1800: p->ops->creatematrix = DMCreateMatrix_Composite;
1801: p->ops->getcoloring = DMCreateColoring_Composite;
1802: p->ops->globaltolocalbegin = DMGlobalToLocalBegin_Composite;
1803: p->ops->globaltolocalend = DMGlobalToLocalEnd_Composite;
1804: p->ops->localtoglobalbegin = DMLocalToGlobalBegin_Composite;
1805: p->ops->localtoglobalend = DMLocalToGlobalEnd_Composite;
1806: p->ops->localtolocalbegin = DMLocalToLocalBegin_Composite;
1807: p->ops->localtolocalend = DMLocalToLocalEnd_Composite;
1808: p->ops->destroy = DMDestroy_Composite;
1809: p->ops->view = DMView_Composite;
1810: p->ops->setup = DMSetUp_Composite;
1812: PetscObjectComposeFunction((PetscObject)p,"DMSetUpGLVisViewer_C",DMSetUpGLVisViewer_Composite);
1813: return(0);
1814: }
1816: /*@
1817: DMCompositeCreate - Creates a vector packer, used to generate "composite"
1818: vectors made up of several subvectors.
1820: Collective on MPI_Comm
1822: Input Parameter:
1823: . comm - the processors that will share the global vector
1825: Output Parameters:
1826: . packer - the packer object
1828: Level: advanced
1830: .seealso DMDestroy(), DMCompositeAddDM(), DMCompositeScatter(), DMCOMPOSITE,DMCreate()
1831: DMCompositeGather(), DMCreateGlobalVector(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess()
1832: DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(), DMCompositeGetEntries()
1834: @*/
1835: PetscErrorCode DMCompositeCreate(MPI_Comm comm,DM *packer)
1836: {
1841: DMCreate(comm,packer);
1842: DMSetType(*packer,DMCOMPOSITE);
1843: return(0);
1844: }