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;
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: }
153: /*@C
154: DMCompositeGetAccess - Allows one to access the individual packed vectors in their global
155: representation.
157: Collective on dm
159: Input Parameters:
160: + dm - the packer object
161: - gvec - the global vector
163: Output Parameters:
164: . Vec* ... - the packed parallel vectors, NULL for those that are not needed
166: Notes:
167: 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: VecLockReadPush(*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 dm
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:
240: Use DMCompositeRestoreAccessArray() to return the vectors when you no longer need them
242: Level: advanced
244: .seealso: DMCompositeGetAccess(), DMCompositeGetEntries(), DMCompositeScatter(), DMCompositeGather()
245: @*/
246: PetscErrorCode DMCompositeGetAccessArray(DM dm,Vec pvec,PetscInt nwanted,const PetscInt *wanted,Vec *vecs)
247: {
248: PetscErrorCode ierr;
249: struct DMCompositeLink *link;
250: PetscInt i,wnum;
251: DM_Composite *com = (DM_Composite*)dm->data;
252: PetscInt readonly;
253: PetscBool flg;
258: PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
259: if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
260: if (!com->setup) {
261: DMSetUp(dm);
262: }
264: VecLockGet(pvec,&readonly);
265: for (i=0,wnum=0,link=com->next; link && wnum<nwanted; i++,link=link->next) {
266: if (!wanted || i == wanted[wnum]) {
267: Vec v;
268: DMGetGlobalVector(link->dm,&v);
269: if (readonly) {
270: const PetscScalar *array;
271: VecGetArrayRead(pvec,&array);
272: VecPlaceArray(v,array+link->rstart);
273: VecLockReadPush(v);
274: VecRestoreArrayRead(pvec,&array);
275: } else {
276: PetscScalar *array;
277: VecGetArray(pvec,&array);
278: VecPlaceArray(v,array+link->rstart);
279: VecRestoreArray(pvec,&array);
280: }
281: vecs[wnum++] = v;
282: }
283: }
284: return(0);
285: }
287: /*@C
288: DMCompositeGetLocalAccessArray - Allows one to access the individual
289: packed vectors in their local representation.
291: Collective on dm.
293: Input Parameters:
294: + dm - the packer object
295: . pvec - packed vector
296: . nwanted - number of vectors wanted
297: - wanted - sorted array of vectors wanted, or NULL to get all vectors
299: Output Parameters:
300: . vecs - array of requested local vectors (must be allocated)
302: Notes:
303: Use DMCompositeRestoreLocalAccessArray() to return the vectors
304: when you no longer need them.
306: Level: advanced
308: .seealso: DMCompositeRestoreLocalAccessArray(), DMCompositeGetAccess(),
309: DMCompositeGetEntries(), DMCompositeScatter(), DMCompositeGather()
310: @*/
311: PetscErrorCode DMCompositeGetLocalAccessArray(DM dm,Vec pvec,PetscInt nwanted,const PetscInt *wanted,Vec *vecs)
312: {
313: PetscErrorCode ierr;
314: struct DMCompositeLink *link;
315: PetscInt i,wnum;
316: DM_Composite *com = (DM_Composite*)dm->data;
317: PetscInt readonly;
318: PetscInt nlocal = 0;
319: PetscBool flg;
324: PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
325: if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
326: if (!com->setup) {
327: DMSetUp(dm);
328: }
330: VecLockGet(pvec,&readonly);
331: for (i=0,wnum=0,link=com->next; link && wnum<nwanted; i++,link=link->next) {
332: if (!wanted || i == wanted[wnum]) {
333: Vec v;
334: DMGetLocalVector(link->dm,&v);
335: if (readonly) {
336: const PetscScalar *array;
337: VecGetArrayRead(pvec,&array);
338: VecPlaceArray(v,array+nlocal);
339: // 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
340: VecLockReadPush(v);
341: VecRestoreArrayRead(pvec,&array);
342: } else {
343: PetscScalar *array;
344: VecGetArray(pvec,&array);
345: VecPlaceArray(v,array+nlocal);
346: VecRestoreArray(pvec,&array);
347: }
348: vecs[wnum++] = v;
349: }
351: nlocal += link->nlocal;
352: }
354: return(0);
355: }
357: /*@C
358: DMCompositeRestoreAccess - Returns the vectors obtained with DMCompositeGetAccess()
359: representation.
361: Collective on dm
363: Input Parameters:
364: + dm - the packer object
365: . gvec - the global vector
366: - Vec* ... - the individual parallel vectors, NULL for those that are not needed
368: Level: advanced
370: .seealso DMCompositeAddDM(), DMCreateGlobalVector(),
371: DMCompositeGather(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeScatter(),
372: DMCompositeRestoreAccess(), DMCompositeGetAccess()
374: @*/
375: PetscErrorCode DMCompositeRestoreAccess(DM dm,Vec gvec,...)
376: {
377: va_list Argp;
378: PetscErrorCode ierr;
379: struct DMCompositeLink *next;
380: DM_Composite *com = (DM_Composite*)dm->data;
381: PetscInt readonly;
382: PetscBool flg;
387: PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
388: if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
389: next = com->next;
390: if (!com->setup) {
391: DMSetUp(dm);
392: }
394: VecLockGet(gvec,&readonly);
395: /* loop over packed objects, handling one at at time */
396: va_start(Argp,gvec);
397: while (next) {
398: Vec *vec;
399: vec = va_arg(Argp, Vec*);
400: if (vec) {
401: VecResetArray(*vec);
402: if (readonly) {
403: VecLockReadPop(*vec);
404: }
405: DMRestoreGlobalVector(next->dm,vec);
406: }
407: next = next->next;
408: }
409: va_end(Argp);
410: return(0);
411: }
413: /*@C
414: DMCompositeRestoreAccessArray - Returns the vectors obtained with DMCompositeGetAccessArray()
416: Collective on dm
418: Input Parameters:
419: + dm - the packer object
420: . pvec - packed vector
421: . nwanted - number of vectors wanted
422: . wanted - sorted array of vectors wanted, or NULL to get all vectors
423: - vecs - array of global vectors to return
425: Level: advanced
427: .seealso: DMCompositeRestoreAccess(), DMCompositeRestoreEntries(), DMCompositeScatter(), DMCompositeGather()
428: @*/
429: PetscErrorCode DMCompositeRestoreAccessArray(DM dm,Vec pvec,PetscInt nwanted,const PetscInt *wanted,Vec *vecs)
430: {
431: PetscErrorCode ierr;
432: struct DMCompositeLink *link;
433: PetscInt i,wnum;
434: DM_Composite *com = (DM_Composite*)dm->data;
435: PetscInt readonly;
436: PetscBool flg;
441: PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
442: if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
443: if (!com->setup) {
444: DMSetUp(dm);
445: }
447: VecLockGet(pvec,&readonly);
448: for (i=0,wnum=0,link=com->next; link && wnum<nwanted; i++,link=link->next) {
449: if (!wanted || i == wanted[wnum]) {
450: VecResetArray(vecs[wnum]);
451: if (readonly) {
452: VecLockReadPop(vecs[wnum]);
453: }
454: DMRestoreGlobalVector(link->dm,&vecs[wnum]);
455: wnum++;
456: }
457: }
458: return(0);
459: }
461: /*@C
462: DMCompositeRestoreLocalAccessArray - Returns the vectors obtained with DMCompositeGetLocalAccessArray().
464: Collective on dm.
466: Input Parameters:
467: + dm - the packer object
468: . pvec - packed vector
469: . nwanted - number of vectors wanted
470: . wanted - sorted array of vectors wanted, or NULL to restore all vectors
471: - vecs - array of local vectors to return
473: Level: advanced
475: Notes:
476: nwanted and wanted must match the values given to DMCompositeGetLocalAccessArray()
477: otherwise the call will fail.
479: .seealso: DMCompositeGetLocalAccessArray(), DMCompositeRestoreAccessArray(),
480: DMCompositeRestoreAccess(), DMCompositeRestoreEntries(),
481: DMCompositeScatter(), DMCompositeGather()
482: @*/
483: PetscErrorCode DMCompositeRestoreLocalAccessArray(DM dm,Vec pvec,PetscInt nwanted,const PetscInt *wanted,Vec *vecs)
484: {
485: PetscErrorCode ierr;
486: struct DMCompositeLink *link;
487: PetscInt i,wnum;
488: DM_Composite *com = (DM_Composite*)dm->data;
489: PetscInt readonly;
490: PetscBool flg;
495: PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
496: if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
497: if (!com->setup) {
498: DMSetUp(dm);
499: }
501: VecLockGet(pvec,&readonly);
502: for (i=0,wnum=0,link=com->next; link && wnum<nwanted; i++,link=link->next) {
503: if (!wanted || i == wanted[wnum]) {
504: VecResetArray(vecs[wnum]);
505: if (readonly) {
506: VecLockReadPop(vecs[wnum]);
507: }
508: DMRestoreLocalVector(link->dm,&vecs[wnum]);
509: wnum++;
510: }
511: }
512: return(0);
513: }
515: /*@C
516: DMCompositeScatter - Scatters from a global packed vector into its individual local vectors
518: Collective on dm
520: Input Parameters:
521: + dm - the packer object
522: . gvec - the global vector
523: - Vec ... - the individual sequential vectors, NULL for those that are not needed
525: Level: advanced
527: Notes:
528: DMCompositeScatterArray() is a non-variadic alternative that is often more convenient for library callers and is
529: accessible from Fortran.
531: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(),
532: DMCompositeGather(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
533: DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(), DMCompositeGetEntries()
534: DMCompositeScatterArray()
536: @*/
537: PetscErrorCode DMCompositeScatter(DM dm,Vec gvec,...)
538: {
539: va_list Argp;
540: PetscErrorCode ierr;
541: struct DMCompositeLink *next;
542: PetscInt cnt;
543: DM_Composite *com = (DM_Composite*)dm->data;
544: PetscBool flg;
549: PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
550: if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
551: if (!com->setup) {
552: DMSetUp(dm);
553: }
555: /* loop over packed objects, handling one at at time */
556: va_start(Argp,gvec);
557: for (cnt=3,next=com->next; next; cnt++,next=next->next) {
558: Vec local;
559: local = va_arg(Argp, Vec);
560: if (local) {
561: Vec global;
562: const PetscScalar *array;
564: DMGetGlobalVector(next->dm,&global);
565: VecGetArrayRead(gvec,&array);
566: VecPlaceArray(global,array+next->rstart);
567: DMGlobalToLocalBegin(next->dm,global,INSERT_VALUES,local);
568: DMGlobalToLocalEnd(next->dm,global,INSERT_VALUES,local);
569: VecRestoreArrayRead(gvec,&array);
570: VecResetArray(global);
571: DMRestoreGlobalVector(next->dm,&global);
572: }
573: }
574: va_end(Argp);
575: return(0);
576: }
578: /*@
579: DMCompositeScatterArray - Scatters from a global packed vector into its individual local vectors
581: Collective on dm
583: Input Parameters:
584: + dm - the packer object
585: . gvec - the global vector
586: - lvecs - array of local vectors, NULL for any that are not needed
588: Level: advanced
590: Note:
591: This is a non-variadic alternative to DMCompositeScatter()
593: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector()
594: DMCompositeGather(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
595: DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(), DMCompositeGetEntries()
597: @*/
598: PetscErrorCode DMCompositeScatterArray(DM dm,Vec gvec,Vec *lvecs)
599: {
600: PetscErrorCode ierr;
601: struct DMCompositeLink *next;
602: PetscInt i;
603: DM_Composite *com = (DM_Composite*)dm->data;
604: PetscBool flg;
609: PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
610: if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
611: if (!com->setup) {
612: DMSetUp(dm);
613: }
615: /* loop over packed objects, handling one at at time */
616: for (i=0,next=com->next; next; next=next->next,i++) {
617: if (lvecs[i]) {
618: Vec global;
619: const PetscScalar *array;
621: DMGetGlobalVector(next->dm,&global);
622: VecGetArrayRead(gvec,&array);
623: VecPlaceArray(global,(PetscScalar*)array+next->rstart);
624: DMGlobalToLocalBegin(next->dm,global,INSERT_VALUES,lvecs[i]);
625: DMGlobalToLocalEnd(next->dm,global,INSERT_VALUES,lvecs[i]);
626: VecRestoreArrayRead(gvec,&array);
627: VecResetArray(global);
628: DMRestoreGlobalVector(next->dm,&global);
629: }
630: }
631: return(0);
632: }
634: /*@C
635: DMCompositeGather - Gathers into a global packed vector from its individual local vectors
637: Collective on dm
639: Input Parameters:
640: + dm - the packer object
641: . gvec - the global vector
642: . imode - INSERT_VALUES or ADD_VALUES
643: - Vec ... - the individual sequential vectors, NULL for any that are not needed
645: Level: advanced
647: Not available from Fortran, Fortran users can use DMCompositeGatherArray()
649: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(),
650: DMCompositeScatter(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
651: DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(), DMCompositeGetEntries()
653: @*/
654: PetscErrorCode DMCompositeGather(DM dm,InsertMode imode,Vec gvec,...)
655: {
656: va_list Argp;
657: PetscErrorCode ierr;
658: struct DMCompositeLink *next;
659: DM_Composite *com = (DM_Composite*)dm->data;
660: PetscInt cnt;
661: PetscBool flg;
666: PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
667: if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
668: if (!com->setup) {
669: DMSetUp(dm);
670: }
672: /* loop over packed objects, handling one at at time */
673: va_start(Argp,gvec);
674: for (cnt=3,next=com->next; next; cnt++,next=next->next) {
675: Vec local;
676: local = va_arg(Argp, Vec);
677: if (local) {
678: PetscScalar *array;
679: Vec global;
681: DMGetGlobalVector(next->dm,&global);
682: VecGetArray(gvec,&array);
683: VecPlaceArray(global,array+next->rstart);
684: DMLocalToGlobalBegin(next->dm,local,imode,global);
685: DMLocalToGlobalEnd(next->dm,local,imode,global);
686: VecRestoreArray(gvec,&array);
687: VecResetArray(global);
688: DMRestoreGlobalVector(next->dm,&global);
689: }
690: }
691: va_end(Argp);
692: return(0);
693: }
695: /*@
696: DMCompositeGatherArray - Gathers into a global packed vector from its individual local vectors
698: Collective on dm
700: Input Parameters:
701: + dm - the packer object
702: . gvec - the global vector
703: . imode - INSERT_VALUES or ADD_VALUES
704: - lvecs - the individual sequential vectors, NULL for any that are not needed
706: Level: advanced
708: Notes:
709: This is a non-variadic alternative to DMCompositeGather().
711: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(),
712: DMCompositeScatter(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
713: DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(), DMCompositeGetEntries(),
714: @*/
715: PetscErrorCode DMCompositeGatherArray(DM dm,InsertMode imode,Vec gvec,Vec *lvecs)
716: {
717: PetscErrorCode ierr;
718: struct DMCompositeLink *next;
719: DM_Composite *com = (DM_Composite*)dm->data;
720: PetscInt i;
721: PetscBool flg;
726: PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
727: if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
728: if (!com->setup) {
729: DMSetUp(dm);
730: }
732: /* loop over packed objects, handling one at at time */
733: for (next=com->next,i=0; next; next=next->next,i++) {
734: if (lvecs[i]) {
735: PetscScalar *array;
736: Vec global;
738: DMGetGlobalVector(next->dm,&global);
739: VecGetArray(gvec,&array);
740: VecPlaceArray(global,array+next->rstart);
741: DMLocalToGlobalBegin(next->dm,lvecs[i],imode,global);
742: DMLocalToGlobalEnd(next->dm,lvecs[i],imode,global);
743: VecRestoreArray(gvec,&array);
744: VecResetArray(global);
745: DMRestoreGlobalVector(next->dm,&global);
746: }
747: }
748: return(0);
749: }
751: /*@
752: DMCompositeAddDM - adds a DM vector to a DMComposite
754: Collective on dm
756: Input Parameters:
757: + dmc - the DMComposite (packer) object
758: - dm - the DM object
760: Level: advanced
762: .seealso DMDestroy(), DMCompositeGather(), DMCompositeAddDM(), DMCreateGlobalVector(),
763: DMCompositeScatter(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
764: DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(), DMCompositeGetEntries()
766: @*/
767: PetscErrorCode DMCompositeAddDM(DM dmc,DM dm)
768: {
769: PetscErrorCode ierr;
770: PetscInt n,nlocal;
771: struct DMCompositeLink *mine,*next;
772: Vec global,local;
773: DM_Composite *com = (DM_Composite*)dmc->data;
774: PetscBool flg;
779: PetscObjectTypeCompare((PetscObject)dmc,DMCOMPOSITE,&flg);
780: if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
781: next = com->next;
782: if (com->setup) SETERRQ(PetscObjectComm((PetscObject)dmc),PETSC_ERR_ARG_WRONGSTATE,"Cannot add a DM once you have used the DMComposite");
784: /* create new link */
785: PetscNew(&mine);
786: PetscObjectReference((PetscObject)dm);
787: DMGetGlobalVector(dm,&global);
788: VecGetLocalSize(global,&n);
789: DMRestoreGlobalVector(dm,&global);
790: DMGetLocalVector(dm,&local);
791: VecGetSize(local,&nlocal);
792: DMRestoreLocalVector(dm,&local);
794: mine->n = n;
795: mine->nlocal = nlocal;
796: mine->dm = dm;
797: mine->next = NULL;
798: com->n += n;
799: com->nghost += nlocal;
801: /* add to end of list */
802: if (!next) com->next = mine;
803: else {
804: while (next->next) next = next->next;
805: next->next = mine;
806: }
807: com->nDM++;
808: com->nmine++;
809: return(0);
810: }
812: #include <petscdraw.h>
813: PETSC_EXTERN PetscErrorCode VecView_MPI(Vec,PetscViewer);
814: PetscErrorCode VecView_DMComposite(Vec gvec,PetscViewer viewer)
815: {
816: DM dm;
817: PetscErrorCode ierr;
818: struct DMCompositeLink *next;
819: PetscBool isdraw;
820: DM_Composite *com;
823: VecGetDM(gvec, &dm);
824: if (!dm) SETERRQ(PetscObjectComm((PetscObject)gvec),PETSC_ERR_ARG_WRONG,"Vector not generated from a DMComposite");
825: com = (DM_Composite*)dm->data;
826: next = com->next;
828: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
829: if (!isdraw) {
830: /* do I really want to call this? */
831: VecView_MPI(gvec,viewer);
832: } else {
833: PetscInt cnt = 0;
835: /* loop over packed objects, handling one at at time */
836: while (next) {
837: Vec vec;
838: const PetscScalar *array;
839: PetscInt bs;
841: /* Should use VecGetSubVector() eventually, but would need to forward the DM for that to work */
842: DMGetGlobalVector(next->dm,&vec);
843: VecGetArrayRead(gvec,&array);
844: VecPlaceArray(vec,(PetscScalar*)array+next->rstart);
845: VecRestoreArrayRead(gvec,&array);
846: VecView(vec,viewer);
847: VecResetArray(vec);
848: VecGetBlockSize(vec,&bs);
849: DMRestoreGlobalVector(next->dm,&vec);
850: PetscViewerDrawBaseAdd(viewer,bs);
851: cnt += bs;
852: next = next->next;
853: }
854: PetscViewerDrawBaseAdd(viewer,-cnt);
855: }
856: return(0);
857: }
859: PetscErrorCode DMCreateGlobalVector_Composite(DM dm,Vec *gvec)
860: {
862: DM_Composite *com = (DM_Composite*)dm->data;
866: DMSetFromOptions(dm);
867: DMSetUp(dm);
868: VecCreate(PetscObjectComm((PetscObject)dm),gvec);
869: VecSetType(*gvec,dm->vectype);
870: VecSetSizes(*gvec,com->n,com->N);
871: VecSetDM(*gvec, dm);
872: VecSetOperation(*gvec,VECOP_VIEW,(void (*)(void))VecView_DMComposite);
873: return(0);
874: }
876: PetscErrorCode DMCreateLocalVector_Composite(DM dm,Vec *lvec)
877: {
879: DM_Composite *com = (DM_Composite*)dm->data;
883: if (!com->setup) {
884: DMSetFromOptions(dm);
885: DMSetUp(dm);
886: }
887: VecCreate(PETSC_COMM_SELF,lvec);
888: VecSetType(*lvec,dm->vectype);
889: VecSetSizes(*lvec,com->nghost,PETSC_DECIDE);
890: VecSetDM(*lvec, dm);
891: return(0);
892: }
894: /*@C
895: DMCompositeGetISLocalToGlobalMappings - gets an ISLocalToGlobalMapping for each DM in the DMComposite, maps to the composite global space
897: Collective on DM
899: Input Parameter:
900: . dm - the packer object
902: Output Parameters:
903: . ltogs - the individual mappings for each packed vector. Note that this includes
904: all the ghost points that individual ghosted DMDA's may have.
906: Level: advanced
908: Notes:
909: Each entry of ltogs should be destroyed with ISLocalToGlobalMappingDestroy(), the ltogs array should be freed with PetscFree().
911: Not available from Fortran
913: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(),
914: DMCompositeGather(), DMCompositeCreate(), DMCompositeGetAccess(), DMCompositeScatter(),
915: DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(),DMCompositeGetEntries()
917: @*/
918: PetscErrorCode DMCompositeGetISLocalToGlobalMappings(DM dm,ISLocalToGlobalMapping **ltogs)
919: {
920: PetscErrorCode ierr;
921: PetscInt i,*idx,n,cnt;
922: struct DMCompositeLink *next;
923: PetscMPIInt rank;
924: DM_Composite *com = (DM_Composite*)dm->data;
925: PetscBool flg;
929: PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
930: if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
931: DMSetUp(dm);
932: PetscMalloc1(com->nDM,ltogs);
933: next = com->next;
934: MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);
936: /* loop over packed objects, handling one at at time */
937: cnt = 0;
938: while (next) {
939: ISLocalToGlobalMapping ltog;
940: PetscMPIInt size;
941: const PetscInt *suboff,*indices;
942: Vec global;
944: /* Get sub-DM global indices for each local dof */
945: DMGetLocalToGlobalMapping(next->dm,<og);
946: ISLocalToGlobalMappingGetSize(ltog,&n);
947: ISLocalToGlobalMappingGetIndices(ltog,&indices);
948: PetscMalloc1(n,&idx);
950: /* Get the offsets for the sub-DM global vector */
951: DMGetGlobalVector(next->dm,&global);
952: VecGetOwnershipRanges(global,&suboff);
953: MPI_Comm_size(PetscObjectComm((PetscObject)global),&size);
955: /* Shift the sub-DM definition of the global space to the composite global space */
956: for (i=0; i<n; i++) {
957: PetscInt subi = indices[i],lo = 0,hi = size,t;
958: /* There's no consensus on what a negative index means,
959: except for skipping when setting the values in vectors and matrices */
960: if (subi < 0) { idx[i] = subi - next->grstarts[rank]; continue; }
961: /* Binary search to find which rank owns subi */
962: while (hi-lo > 1) {
963: t = lo + (hi-lo)/2;
964: if (suboff[t] > subi) hi = t;
965: else lo = t;
966: }
967: idx[i] = subi - suboff[lo] + next->grstarts[lo];
968: }
969: ISLocalToGlobalMappingRestoreIndices(ltog,&indices);
970: ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)dm),1,n,idx,PETSC_OWN_POINTER,&(*ltogs)[cnt]);
971: DMRestoreGlobalVector(next->dm,&global);
972: next = next->next;
973: cnt++;
974: }
975: return(0);
976: }
978: /*@C
979: DMCompositeGetLocalISs - Gets index sets for each component of a composite local vector
981: Not Collective
983: Input Parameter:
984: . dm - composite DM
986: Output Parameter:
987: . is - array of serial index sets for each each component of the DMComposite
989: Level: intermediate
991: Notes:
992: At present, a composite local vector does not normally exist. This function is used to provide index sets for
993: MatGetLocalSubMatrix(). In the future, the scatters for each entry in the DMComposite may be be merged into a single
994: scatter to a composite local vector. The user should not typically need to know which is being done.
996: To get the composite global indices at all local points (including ghosts), use DMCompositeGetISLocalToGlobalMappings().
998: To get index sets for pieces of the composite global vector, use DMCompositeGetGlobalISs().
1000: Each returned IS should be destroyed with ISDestroy(), the array should be freed with PetscFree().
1002: Not available from Fortran
1004: .seealso: DMCompositeGetGlobalISs(), DMCompositeGetISLocalToGlobalMappings(), MatGetLocalSubMatrix(), MatCreateLocalRef()
1005: @*/
1006: PetscErrorCode DMCompositeGetLocalISs(DM dm,IS **is)
1007: {
1008: PetscErrorCode ierr;
1009: DM_Composite *com = (DM_Composite*)dm->data;
1010: struct DMCompositeLink *link;
1011: PetscInt cnt,start;
1012: PetscBool flg;
1017: PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
1018: if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
1019: PetscMalloc1(com->nmine,is);
1020: for (cnt=0,start=0,link=com->next; link; start+=link->nlocal,cnt++,link=link->next) {
1021: PetscInt bs;
1022: ISCreateStride(PETSC_COMM_SELF,link->nlocal,start,1,&(*is)[cnt]);
1023: DMGetBlockSize(link->dm,&bs);
1024: ISSetBlockSize((*is)[cnt],bs);
1025: }
1026: return(0);
1027: }
1029: /*@C
1030: DMCompositeGetGlobalISs - Gets the index sets for each composed object
1032: Collective on dm
1034: Input Parameter:
1035: . dm - the packer object
1037: Output Parameters:
1038: . is - the array of index sets
1040: Level: advanced
1042: Notes:
1043: The is entries should be destroyed with ISDestroy(), the is array should be freed with PetscFree()
1045: These could be used to extract a subset of vector entries for a "multi-physics" preconditioner
1047: Use DMCompositeGetLocalISs() for index sets in the packed local numbering, and
1048: DMCompositeGetISLocalToGlobalMappings() for to map local sub-DM (including ghost) indices to packed global
1049: indices.
1051: Fortran Notes:
1053: The output argument 'is' must be an allocated array of sufficient length, which can be learned using DMCompositeGetNumberDM().
1055: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(),
1056: DMCompositeGather(), DMCompositeCreate(), DMCompositeGetAccess(), DMCompositeScatter(),
1057: DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(),DMCompositeGetEntries()
1059: @*/
1060: PetscErrorCode DMCompositeGetGlobalISs(DM dm,IS *is[])
1061: {
1062: PetscErrorCode ierr;
1063: PetscInt cnt = 0;
1064: struct DMCompositeLink *next;
1065: PetscMPIInt rank;
1066: DM_Composite *com = (DM_Composite*)dm->data;
1067: PetscBool flg;
1071: PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
1072: if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
1073: if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Must call DMSetUp() before");
1074: PetscMalloc1(com->nDM,is);
1075: next = com->next;
1076: MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);
1078: /* loop over packed objects, handling one at at time */
1079: while (next) {
1080: PetscDS prob;
1082: ISCreateStride(PetscObjectComm((PetscObject)dm),next->n,next->grstart,1,&(*is)[cnt]);
1083: DMGetDS(dm, &prob);
1084: if (prob) {
1085: MatNullSpace space;
1086: Mat pmat;
1087: PetscObject disc;
1088: PetscInt Nf;
1090: PetscDSGetNumFields(prob, &Nf);
1091: if (cnt < Nf) {
1092: PetscDSGetDiscretization(prob, cnt, &disc);
1093: PetscObjectQuery(disc, "nullspace", (PetscObject*) &space);
1094: if (space) {PetscObjectCompose((PetscObject) (*is)[cnt], "nullspace", (PetscObject) space);}
1095: PetscObjectQuery(disc, "nearnullspace", (PetscObject*) &space);
1096: if (space) {PetscObjectCompose((PetscObject) (*is)[cnt], "nearnullspace", (PetscObject) space);}
1097: PetscObjectQuery(disc, "pmat", (PetscObject*) &pmat);
1098: if (pmat) {PetscObjectCompose((PetscObject) (*is)[cnt], "pmat", (PetscObject) pmat);}
1099: }
1100: }
1101: cnt++;
1102: next = next->next;
1103: }
1104: return(0);
1105: }
1107: PetscErrorCode DMCreateFieldIS_Composite(DM dm, PetscInt *numFields,char ***fieldNames, IS **fields)
1108: {
1109: PetscInt nDM;
1110: DM *dms;
1111: PetscInt i;
1115: DMCompositeGetNumberDM(dm, &nDM);
1116: if (numFields) *numFields = nDM;
1117: DMCompositeGetGlobalISs(dm, fields);
1118: if (fieldNames) {
1119: PetscMalloc1(nDM, &dms);
1120: PetscMalloc1(nDM, fieldNames);
1121: DMCompositeGetEntriesArray(dm, dms);
1122: for (i=0; i<nDM; i++) {
1123: char buf[256];
1124: const char *splitname;
1126: /* Split naming precedence: object name, prefix, number */
1127: splitname = ((PetscObject) dm)->name;
1128: if (!splitname) {
1129: PetscObjectGetOptionsPrefix((PetscObject)dms[i],&splitname);
1130: if (splitname) {
1131: size_t len;
1132: PetscStrncpy(buf,splitname,sizeof(buf));
1133: buf[sizeof(buf) - 1] = 0;
1134: PetscStrlen(buf,&len);
1135: if (buf[len-1] == '_') buf[len-1] = 0; /* Remove trailing underscore if it was used */
1136: splitname = buf;
1137: }
1138: }
1139: if (!splitname) {
1140: PetscSNPrintf(buf,sizeof(buf),"%D",i);
1141: splitname = buf;
1142: }
1143: PetscStrallocpy(splitname,&(*fieldNames)[i]);
1144: }
1145: PetscFree(dms);
1146: }
1147: return(0);
1148: }
1150: /*
1151: This could take over from DMCreateFieldIS(), as it is more general,
1152: making DMCreateFieldIS() a special case -- calling with dmlist == NULL;
1153: At this point it's probably best to be less intrusive, however.
1154: */
1155: PetscErrorCode DMCreateFieldDecomposition_Composite(DM dm, PetscInt *len,char ***namelist, IS **islist, DM **dmlist)
1156: {
1157: PetscInt nDM;
1158: PetscInt i;
1162: DMCreateFieldIS_Composite(dm, len, namelist, islist);
1163: if (dmlist) {
1164: DMCompositeGetNumberDM(dm, &nDM);
1165: PetscMalloc1(nDM, dmlist);
1166: DMCompositeGetEntriesArray(dm, *dmlist);
1167: for (i=0; i<nDM; i++) {
1168: PetscObjectReference((PetscObject)((*dmlist)[i]));
1169: }
1170: }
1171: return(0);
1172: }
1174: /* -------------------------------------------------------------------------------------*/
1175: /*@C
1176: DMCompositeGetLocalVectors - Gets local vectors for each part of a DMComposite.
1177: Use DMCompositeRestoreLocalVectors() to return them.
1179: Not Collective
1181: Input Parameter:
1182: . dm - the packer object
1184: Output Parameter:
1185: . Vec ... - the individual sequential Vecs
1187: Level: advanced
1189: Not available from Fortran
1191: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(),
1192: DMCompositeGather(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
1193: DMCompositeRestoreLocalVectors(), DMCompositeScatter(), DMCompositeGetEntries()
1195: @*/
1196: PetscErrorCode DMCompositeGetLocalVectors(DM dm,...)
1197: {
1198: va_list Argp;
1199: PetscErrorCode ierr;
1200: struct DMCompositeLink *next;
1201: DM_Composite *com = (DM_Composite*)dm->data;
1202: PetscBool flg;
1206: PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
1207: if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
1208: next = com->next;
1209: /* loop over packed objects, handling one at at time */
1210: va_start(Argp,dm);
1211: while (next) {
1212: Vec *vec;
1213: vec = va_arg(Argp, Vec*);
1214: if (vec) {DMGetLocalVector(next->dm,vec);}
1215: next = next->next;
1216: }
1217: va_end(Argp);
1218: return(0);
1219: }
1221: /*@C
1222: DMCompositeRestoreLocalVectors - Restores local vectors for each part of a DMComposite.
1224: Not Collective
1226: Input Parameter:
1227: . dm - the packer object
1229: Output Parameter:
1230: . Vec ... - the individual sequential Vecs
1232: Level: advanced
1234: Not available from Fortran
1236: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(),
1237: DMCompositeGather(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
1238: DMCompositeGetLocalVectors(), DMCompositeScatter(), DMCompositeGetEntries()
1240: @*/
1241: PetscErrorCode DMCompositeRestoreLocalVectors(DM dm,...)
1242: {
1243: va_list Argp;
1244: PetscErrorCode ierr;
1245: struct DMCompositeLink *next;
1246: DM_Composite *com = (DM_Composite*)dm->data;
1247: PetscBool flg;
1251: PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
1252: if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
1253: next = com->next;
1254: /* loop over packed objects, handling one at at time */
1255: va_start(Argp,dm);
1256: while (next) {
1257: Vec *vec;
1258: vec = va_arg(Argp, Vec*);
1259: if (vec) {DMRestoreLocalVector(next->dm,vec);}
1260: next = next->next;
1261: }
1262: va_end(Argp);
1263: return(0);
1264: }
1266: /* -------------------------------------------------------------------------------------*/
1267: /*@C
1268: DMCompositeGetEntries - Gets the DM for each entry in a DMComposite.
1270: Not Collective
1272: Input Parameter:
1273: . dm - the packer object
1275: Output Parameter:
1276: . DM ... - the individual entries (DMs)
1278: Level: advanced
1280: Fortran Notes:
1281: Available as DMCompositeGetEntries() for one output DM, DMCompositeGetEntries2() for 2, etc
1283: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(), DMCompositeGetEntriesArray()
1284: DMCompositeGather(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
1285: DMCompositeRestoreLocalVectors(), DMCompositeGetLocalVectors(), DMCompositeScatter(),
1286: DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors()
1288: @*/
1289: PetscErrorCode DMCompositeGetEntries(DM dm,...)
1290: {
1291: va_list Argp;
1292: struct DMCompositeLink *next;
1293: DM_Composite *com = (DM_Composite*)dm->data;
1294: PetscBool flg;
1295: PetscErrorCode ierr;
1299: PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
1300: if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
1301: next = com->next;
1302: /* loop over packed objects, handling one at at time */
1303: va_start(Argp,dm);
1304: while (next) {
1305: DM *dmn;
1306: dmn = va_arg(Argp, DM*);
1307: if (dmn) *dmn = next->dm;
1308: next = next->next;
1309: }
1310: va_end(Argp);
1311: return(0);
1312: }
1314: /*@C
1315: DMCompositeGetEntriesArray - Gets the DM for each entry in a DMComposite.
1317: Not Collective
1319: Input Parameter:
1320: . dm - the packer object
1322: Output Parameter:
1323: . dms - array of sufficient length (see DMCompositeGetNumberDM()) to hold the individual DMs
1325: Level: advanced
1327: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(), DMCompositeGetEntries()
1328: DMCompositeGather(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
1329: DMCompositeRestoreLocalVectors(), DMCompositeGetLocalVectors(), DMCompositeScatter(),
1330: DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors()
1332: @*/
1333: PetscErrorCode DMCompositeGetEntriesArray(DM dm,DM dms[])
1334: {
1335: struct DMCompositeLink *next;
1336: DM_Composite *com = (DM_Composite*)dm->data;
1337: PetscInt i;
1338: PetscBool flg;
1339: PetscErrorCode ierr;
1343: PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
1344: if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
1345: /* loop over packed objects, handling one at at time */
1346: for (next=com->next,i=0; next; next=next->next,i++) dms[i] = next->dm;
1347: return(0);
1348: }
1350: typedef struct {
1351: DM dm;
1352: PetscViewer *subv;
1353: Vec *vecs;
1354: } GLVisViewerCtx;
1356: static PetscErrorCode DestroyGLVisViewerCtx_Private(void *vctx)
1357: {
1358: GLVisViewerCtx *ctx = (GLVisViewerCtx*)vctx;
1359: PetscInt i,n;
1363: DMCompositeGetNumberDM(ctx->dm,&n);
1364: for (i = 0; i < n; i++) {
1365: PetscViewerDestroy(&ctx->subv[i]);
1366: }
1367: PetscFree2(ctx->subv,ctx->vecs);
1368: DMDestroy(&ctx->dm);
1369: PetscFree(ctx);
1370: return(0);
1371: }
1373: static PetscErrorCode DMCompositeSampleGLVisFields_Private(PetscObject oX, PetscInt nf, PetscObject oXfield[], void *vctx)
1374: {
1375: Vec X = (Vec)oX;
1376: GLVisViewerCtx *ctx = (GLVisViewerCtx*)vctx;
1377: PetscInt i,n,cumf;
1381: DMCompositeGetNumberDM(ctx->dm,&n);
1382: DMCompositeGetAccessArray(ctx->dm,X,n,NULL,ctx->vecs);
1383: for (i = 0, cumf = 0; i < n; i++) {
1384: PetscErrorCode (*g2l)(PetscObject,PetscInt,PetscObject[],void*);
1385: void *fctx;
1386: PetscInt nfi;
1388: PetscViewerGLVisGetFields_Private(ctx->subv[i],&nfi,NULL,NULL,&g2l,NULL,&fctx);
1389: if (!nfi) continue;
1390: if (g2l) {
1391: (*g2l)((PetscObject)ctx->vecs[i],nfi,oXfield+cumf,fctx);
1392: } else {
1393: VecCopy(ctx->vecs[i],(Vec)(oXfield[cumf]));
1394: }
1395: cumf += nfi;
1396: }
1397: DMCompositeRestoreAccessArray(ctx->dm,X,n,NULL,ctx->vecs);
1398: return(0);
1399: }
1401: static PetscErrorCode DMSetUpGLVisViewer_Composite(PetscObject odm, PetscViewer viewer)
1402: {
1403: DM dm = (DM)odm, *dms;
1404: Vec *Ufds;
1405: GLVisViewerCtx *ctx;
1406: PetscInt i,n,tnf,*sdim;
1407: char **fecs;
1411: PetscNew(&ctx);
1412: PetscObjectReference((PetscObject)dm);
1413: ctx->dm = dm;
1414: DMCompositeGetNumberDM(dm,&n);
1415: PetscMalloc1(n,&dms);
1416: DMCompositeGetEntriesArray(dm,dms);
1417: PetscMalloc2(n,&ctx->subv,n,&ctx->vecs);
1418: for (i = 0, tnf = 0; i < n; i++) {
1419: PetscInt nf;
1421: PetscViewerCreate(PetscObjectComm(odm),&ctx->subv[i]);
1422: PetscViewerSetType(ctx->subv[i],PETSCVIEWERGLVIS);
1423: PetscViewerGLVisSetDM_Private(ctx->subv[i],(PetscObject)dms[i]);
1424: PetscViewerGLVisGetFields_Private(ctx->subv[i],&nf,NULL,NULL,NULL,NULL,NULL);
1425: tnf += nf;
1426: }
1427: PetscFree(dms);
1428: PetscMalloc3(tnf,&fecs,tnf,&sdim,tnf,&Ufds);
1429: for (i = 0, tnf = 0; i < n; i++) {
1430: PetscInt *sd,nf,f;
1431: const char **fec;
1432: Vec *Uf;
1434: PetscViewerGLVisGetFields_Private(ctx->subv[i],&nf,&fec,&sd,NULL,(PetscObject**)&Uf,NULL);
1435: for (f = 0; f < nf; f++) {
1436: PetscStrallocpy(fec[f],&fecs[tnf+f]);
1437: Ufds[tnf+f] = Uf[f];
1438: sdim[tnf+f] = sd[f];
1439: }
1440: tnf += nf;
1441: }
1442: PetscViewerGLVisSetFields(viewer,tnf,(const char**)fecs,sdim,DMCompositeSampleGLVisFields_Private,(PetscObject*)Ufds,ctx,DestroyGLVisViewerCtx_Private);
1443: for (i = 0; i < tnf; i++) {
1444: PetscFree(fecs[i]);
1445: }
1446: PetscFree3(fecs,sdim,Ufds);
1447: return(0);
1448: }
1450: PetscErrorCode DMRefine_Composite(DM dmi,MPI_Comm comm,DM *fine)
1451: {
1452: PetscErrorCode ierr;
1453: struct DMCompositeLink *next;
1454: DM_Composite *com = (DM_Composite*)dmi->data;
1455: DM dm;
1459: if (comm == MPI_COMM_NULL) {
1460: PetscObjectGetComm((PetscObject)dmi,&comm);
1461: }
1462: DMSetUp(dmi);
1463: next = com->next;
1464: DMCompositeCreate(comm,fine);
1466: /* loop over packed objects, handling one at at time */
1467: while (next) {
1468: DMRefine(next->dm,comm,&dm);
1469: DMCompositeAddDM(*fine,dm);
1470: PetscObjectDereference((PetscObject)dm);
1471: next = next->next;
1472: }
1473: return(0);
1474: }
1476: PetscErrorCode DMCoarsen_Composite(DM dmi,MPI_Comm comm,DM *fine)
1477: {
1478: PetscErrorCode ierr;
1479: struct DMCompositeLink *next;
1480: DM_Composite *com = (DM_Composite*)dmi->data;
1481: DM dm;
1485: DMSetUp(dmi);
1486: if (comm == MPI_COMM_NULL) {
1487: PetscObjectGetComm((PetscObject)dmi,&comm);
1488: }
1489: next = com->next;
1490: DMCompositeCreate(comm,fine);
1492: /* loop over packed objects, handling one at at time */
1493: while (next) {
1494: DMCoarsen(next->dm,comm,&dm);
1495: DMCompositeAddDM(*fine,dm);
1496: PetscObjectDereference((PetscObject)dm);
1497: next = next->next;
1498: }
1499: return(0);
1500: }
1502: PetscErrorCode DMCreateInterpolation_Composite(DM coarse,DM fine,Mat *A,Vec *v)
1503: {
1504: PetscErrorCode ierr;
1505: PetscInt m,n,M,N,nDM,i;
1506: struct DMCompositeLink *nextc;
1507: struct DMCompositeLink *nextf;
1508: Vec gcoarse,gfine,*vecs;
1509: DM_Composite *comcoarse = (DM_Composite*)coarse->data;
1510: DM_Composite *comfine = (DM_Composite*)fine->data;
1511: Mat *mats;
1516: DMSetUp(coarse);
1517: DMSetUp(fine);
1518: /* use global vectors only for determining matrix layout */
1519: DMGetGlobalVector(coarse,&gcoarse);
1520: DMGetGlobalVector(fine,&gfine);
1521: VecGetLocalSize(gcoarse,&n);
1522: VecGetLocalSize(gfine,&m);
1523: VecGetSize(gcoarse,&N);
1524: VecGetSize(gfine,&M);
1525: DMRestoreGlobalVector(coarse,&gcoarse);
1526: DMRestoreGlobalVector(fine,&gfine);
1528: nDM = comfine->nDM;
1529: if (nDM != comcoarse->nDM) SETERRQ2(PetscObjectComm((PetscObject)fine),PETSC_ERR_ARG_INCOMP,"Fine DMComposite has %D entries, but coarse has %D",nDM,comcoarse->nDM);
1530: PetscCalloc1(nDM*nDM,&mats);
1531: if (v) {
1532: PetscCalloc1(nDM,&vecs);
1533: }
1535: /* loop over packed objects, handling one at at time */
1536: for (nextc=comcoarse->next,nextf=comfine->next,i=0; nextc; nextc=nextc->next,nextf=nextf->next,i++) {
1537: if (!v) {
1538: DMCreateInterpolation(nextc->dm,nextf->dm,&mats[i*nDM+i],NULL);
1539: } else {
1540: DMCreateInterpolation(nextc->dm,nextf->dm,&mats[i*nDM+i],&vecs[i]);
1541: }
1542: }
1543: MatCreateNest(PetscObjectComm((PetscObject)fine),nDM,NULL,nDM,NULL,mats,A);
1544: if (v) {
1545: VecCreateNest(PetscObjectComm((PetscObject)fine),nDM,NULL,vecs,v);
1546: }
1547: for (i=0; i<nDM*nDM; i++) {MatDestroy(&mats[i]);}
1548: PetscFree(mats);
1549: if (v) {
1550: for (i=0; i<nDM; i++) {VecDestroy(&vecs[i]);}
1551: PetscFree(vecs);
1552: }
1553: return(0);
1554: }
1556: static PetscErrorCode DMGetLocalToGlobalMapping_Composite(DM dm)
1557: {
1558: DM_Composite *com = (DM_Composite*)dm->data;
1559: ISLocalToGlobalMapping *ltogs;
1560: PetscInt i;
1561: PetscErrorCode ierr;
1564: /* Set the ISLocalToGlobalMapping on the new matrix */
1565: DMCompositeGetISLocalToGlobalMappings(dm,<ogs);
1566: ISLocalToGlobalMappingConcatenate(PetscObjectComm((PetscObject)dm),com->nDM,ltogs,&dm->ltogmap);
1567: for (i=0; i<com->nDM; i++) {ISLocalToGlobalMappingDestroy(<ogs[i]);}
1568: PetscFree(ltogs);
1569: return(0);
1570: }
1572: PetscErrorCode DMCreateColoring_Composite(DM dm,ISColoringType ctype,ISColoring *coloring)
1573: {
1574: PetscErrorCode ierr;
1575: PetscInt n,i,cnt;
1576: ISColoringValue *colors;
1577: PetscBool dense = PETSC_FALSE;
1578: ISColoringValue maxcol = 0;
1579: DM_Composite *com = (DM_Composite*)dm->data;
1583: if (ctype == IS_COLORING_LOCAL) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Only global coloring supported");
1584: else if (ctype == IS_COLORING_GLOBAL) {
1585: n = com->n;
1586: } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unknown ISColoringType");
1587: PetscMalloc1(n,&colors); /* freed in ISColoringDestroy() */
1589: PetscOptionsGetBool(((PetscObject)dm)->options,((PetscObject)dm)->prefix,"-dmcomposite_dense_jacobian",&dense,NULL);
1590: if (dense) {
1591: for (i=0; i<n; i++) {
1592: colors[i] = (ISColoringValue)(com->rstart + i);
1593: }
1594: maxcol = com->N;
1595: } else {
1596: struct DMCompositeLink *next = com->next;
1597: PetscMPIInt rank;
1599: MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);
1600: cnt = 0;
1601: while (next) {
1602: ISColoring lcoloring;
1604: DMCreateColoring(next->dm,IS_COLORING_GLOBAL,&lcoloring);
1605: for (i=0; i<lcoloring->N; i++) {
1606: colors[cnt++] = maxcol + lcoloring->colors[i];
1607: }
1608: maxcol += lcoloring->n;
1609: ISColoringDestroy(&lcoloring);
1610: next = next->next;
1611: }
1612: }
1613: ISColoringCreate(PetscObjectComm((PetscObject)dm),maxcol,n,colors,PETSC_OWN_POINTER,coloring);
1614: return(0);
1615: }
1617: PetscErrorCode DMGlobalToLocalBegin_Composite(DM dm,Vec gvec,InsertMode mode,Vec lvec)
1618: {
1619: PetscErrorCode ierr;
1620: struct DMCompositeLink *next;
1621: PetscScalar *garray,*larray;
1622: DM_Composite *com = (DM_Composite*)dm->data;
1628: if (!com->setup) {
1629: DMSetUp(dm);
1630: }
1632: VecGetArray(gvec,&garray);
1633: VecGetArray(lvec,&larray);
1635: /* loop over packed objects, handling one at at time */
1636: next = com->next;
1637: while (next) {
1638: Vec local,global;
1639: PetscInt N;
1641: DMGetGlobalVector(next->dm,&global);
1642: VecGetLocalSize(global,&N);
1643: VecPlaceArray(global,garray);
1644: DMGetLocalVector(next->dm,&local);
1645: VecPlaceArray(local,larray);
1646: DMGlobalToLocalBegin(next->dm,global,mode,local);
1647: DMGlobalToLocalEnd(next->dm,global,mode,local);
1648: VecResetArray(global);
1649: VecResetArray(local);
1650: DMRestoreGlobalVector(next->dm,&global);
1651: DMRestoreLocalVector(next->dm,&local);
1653: larray += next->nlocal;
1654: garray += next->n;
1655: next = next->next;
1656: }
1658: VecRestoreArray(gvec,NULL);
1659: VecRestoreArray(lvec,NULL);
1660: return(0);
1661: }
1663: PetscErrorCode DMGlobalToLocalEnd_Composite(DM dm,Vec gvec,InsertMode mode,Vec lvec)
1664: {
1669: return(0);
1670: }
1672: PetscErrorCode DMLocalToGlobalBegin_Composite(DM dm,Vec lvec,InsertMode mode,Vec gvec)
1673: {
1674: PetscErrorCode ierr;
1675: struct DMCompositeLink *next;
1676: PetscScalar *larray,*garray;
1677: DM_Composite *com = (DM_Composite*)dm->data;
1684: if (!com->setup) {
1685: DMSetUp(dm);
1686: }
1688: VecGetArray(lvec,&larray);
1689: VecGetArray(gvec,&garray);
1691: /* loop over packed objects, handling one at at time */
1692: next = com->next;
1693: while (next) {
1694: Vec global,local;
1696: DMGetLocalVector(next->dm,&local);
1697: VecPlaceArray(local,larray);
1698: DMGetGlobalVector(next->dm,&global);
1699: VecPlaceArray(global,garray);
1700: DMLocalToGlobalBegin(next->dm,local,mode,global);
1701: DMLocalToGlobalEnd(next->dm,local,mode,global);
1702: VecResetArray(local);
1703: VecResetArray(global);
1704: DMRestoreGlobalVector(next->dm,&global);
1705: DMRestoreLocalVector(next->dm,&local);
1707: garray += next->n;
1708: larray += next->nlocal;
1709: next = next->next;
1710: }
1712: VecRestoreArray(gvec,NULL);
1713: VecRestoreArray(lvec,NULL);
1715: return(0);
1716: }
1718: PetscErrorCode DMLocalToGlobalEnd_Composite(DM dm,Vec lvec,InsertMode mode,Vec gvec)
1719: {
1724: return(0);
1725: }
1727: PetscErrorCode DMLocalToLocalBegin_Composite(DM dm,Vec vec1,InsertMode mode,Vec vec2)
1728: {
1729: PetscErrorCode ierr;
1730: struct DMCompositeLink *next;
1731: PetscScalar *array1,*array2;
1732: DM_Composite *com = (DM_Composite*)dm->data;
1739: if (!com->setup) {
1740: DMSetUp(dm);
1741: }
1743: VecGetArray(vec1,&array1);
1744: VecGetArray(vec2,&array2);
1746: /* loop over packed objects, handling one at at time */
1747: next = com->next;
1748: while (next) {
1749: Vec local1,local2;
1751: DMGetLocalVector(next->dm,&local1);
1752: VecPlaceArray(local1,array1);
1753: DMGetLocalVector(next->dm,&local2);
1754: VecPlaceArray(local2,array2);
1755: DMLocalToLocalBegin(next->dm,local1,mode,local2);
1756: DMLocalToLocalEnd(next->dm,local1,mode,local2);
1757: VecResetArray(local2);
1758: DMRestoreLocalVector(next->dm,&local2);
1759: VecResetArray(local1);
1760: DMRestoreLocalVector(next->dm,&local1);
1762: array1 += next->nlocal;
1763: array2 += next->nlocal;
1764: next = next->next;
1765: }
1767: VecRestoreArray(vec1,NULL);
1768: VecRestoreArray(vec2,NULL);
1770: return(0);
1771: }
1773: PetscErrorCode DMLocalToLocalEnd_Composite(DM dm,Vec lvec,InsertMode mode,Vec gvec)
1774: {
1779: return(0);
1780: }
1782: /*MC
1783: DMCOMPOSITE = "composite" - A DM object that is used to manage data for a collection of DMs
1785: Level: intermediate
1787: .seealso: DMType, DM, DMDACreate(), DMCreate(), DMSetType(), DMCompositeCreate()
1788: M*/
1790: PETSC_EXTERN PetscErrorCode DMCreate_Composite(DM p)
1791: {
1793: DM_Composite *com;
1796: PetscNewLog(p,&com);
1797: p->data = com;
1798: com->n = 0;
1799: com->nghost = 0;
1800: com->next = NULL;
1801: com->nDM = 0;
1803: p->ops->createglobalvector = DMCreateGlobalVector_Composite;
1804: p->ops->createlocalvector = DMCreateLocalVector_Composite;
1805: p->ops->getlocaltoglobalmapping = DMGetLocalToGlobalMapping_Composite;
1806: p->ops->createfieldis = DMCreateFieldIS_Composite;
1807: p->ops->createfielddecomposition = DMCreateFieldDecomposition_Composite;
1808: p->ops->refine = DMRefine_Composite;
1809: p->ops->coarsen = DMCoarsen_Composite;
1810: p->ops->createinterpolation = DMCreateInterpolation_Composite;
1811: p->ops->creatematrix = DMCreateMatrix_Composite;
1812: p->ops->getcoloring = DMCreateColoring_Composite;
1813: p->ops->globaltolocalbegin = DMGlobalToLocalBegin_Composite;
1814: p->ops->globaltolocalend = DMGlobalToLocalEnd_Composite;
1815: p->ops->localtoglobalbegin = DMLocalToGlobalBegin_Composite;
1816: p->ops->localtoglobalend = DMLocalToGlobalEnd_Composite;
1817: p->ops->localtolocalbegin = DMLocalToLocalBegin_Composite;
1818: p->ops->localtolocalend = DMLocalToLocalEnd_Composite;
1819: p->ops->destroy = DMDestroy_Composite;
1820: p->ops->view = DMView_Composite;
1821: p->ops->setup = DMSetUp_Composite;
1823: PetscObjectComposeFunction((PetscObject)p,"DMSetUpGLVisViewer_C",DMSetUpGLVisViewer_Composite);
1824: return(0);
1825: }
1827: /*@
1828: DMCompositeCreate - Creates a vector packer, used to generate "composite"
1829: vectors made up of several subvectors.
1831: Collective
1833: Input Parameter:
1834: . comm - the processors that will share the global vector
1836: Output Parameters:
1837: . packer - the packer object
1839: Level: advanced
1841: .seealso DMDestroy(), DMCompositeAddDM(), DMCompositeScatter(), DMCOMPOSITE,DMCreate()
1842: DMCompositeGather(), DMCreateGlobalVector(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess()
1843: DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(), DMCompositeGetEntries()
1845: @*/
1846: PetscErrorCode DMCompositeCreate(MPI_Comm comm,DM *packer)
1847: {
1852: DMCreate(comm,packer);
1853: DMSetType(*packer,DMCOMPOSITE);
1854: return(0);
1855: }