Actual source code: pack.c
petsc-3.10.5 2019-03-28
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:
23: See DMSetApplicationContext() and DMGetApplicationContext() for how to get user information into
24: this routine
26: @*/
27: PetscErrorCode DMCompositeSetCoupling(DM dm,PetscErrorCode (*FormCoupleLocations)(DM,Mat,PetscInt*,PetscInt*,PetscInt,PetscInt,PetscInt,PetscInt))
28: {
29: DM_Composite *com = (DM_Composite*)dm->data;
30: PetscBool flg;
34: PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
35: if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
36: com->FormCoupleLocations = FormCoupleLocations;
37: return(0);
38: }
40: PetscErrorCode DMDestroy_Composite(DM dm)
41: {
42: PetscErrorCode ierr;
43: struct DMCompositeLink *next, *prev;
44: DM_Composite *com = (DM_Composite*)dm->data;
47: next = com->next;
48: while (next) {
49: prev = next;
50: next = next->next;
51: DMDestroy(&prev->dm);
52: PetscFree(prev->grstarts);
53: PetscFree(prev);
54: }
55: PetscObjectComposeFunction((PetscObject)dm,"DMSetUpGLVisViewer_C",NULL);
56: /* This was originally freed in DMDestroy(), but that prevents reference counting of backend objects */
57: PetscFree(com);
58: return(0);
59: }
61: PetscErrorCode DMView_Composite(DM dm,PetscViewer v)
62: {
64: PetscBool iascii;
65: DM_Composite *com = (DM_Composite*)dm->data;
68: PetscObjectTypeCompare((PetscObject)v,PETSCVIEWERASCII,&iascii);
69: if (iascii) {
70: struct DMCompositeLink *lnk = com->next;
71: PetscInt i;
73: PetscViewerASCIIPrintf(v,"DM (%s)\n",((PetscObject)dm)->prefix ? ((PetscObject)dm)->prefix : "no prefix");
74: PetscViewerASCIIPrintf(v," contains %D DMs\n",com->nDM);
75: PetscViewerASCIIPushTab(v);
76: for (i=0; lnk; lnk=lnk->next,i++) {
77: PetscViewerASCIIPrintf(v,"Link %D: DM of type %s\n",i,((PetscObject)lnk->dm)->type_name);
78: PetscViewerASCIIPushTab(v);
79: DMView(lnk->dm,v);
80: PetscViewerASCIIPopTab(v);
81: }
82: PetscViewerASCIIPopTab(v);
83: }
84: return(0);
85: }
87: /* --------------------------------------------------------------------------------------*/
88: PetscErrorCode DMSetUp_Composite(DM dm)
89: {
90: PetscErrorCode ierr;
91: PetscInt nprev = 0;
92: PetscMPIInt rank,size;
93: DM_Composite *com = (DM_Composite*)dm->data;
94: struct DMCompositeLink *next = com->next;
95: PetscLayout map;
98: if (com->setup) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Packer has already been setup");
99: PetscLayoutCreate(PetscObjectComm((PetscObject)dm),&map);
100: PetscLayoutSetLocalSize(map,com->n);
101: PetscLayoutSetSize(map,PETSC_DETERMINE);
102: PetscLayoutSetBlockSize(map,1);
103: PetscLayoutSetUp(map);
104: PetscLayoutGetSize(map,&com->N);
105: PetscLayoutGetRange(map,&com->rstart,NULL);
106: PetscLayoutDestroy(&map);
108: /* now set the rstart for each linked vector */
109: MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);
110: MPI_Comm_size(PetscObjectComm((PetscObject)dm),&size);
111: while (next) {
112: next->rstart = nprev;
113: nprev += next->n;
114: next->grstart = com->rstart + next->rstart;
115: PetscMalloc1(size,&next->grstarts);
116: MPI_Allgather(&next->grstart,1,MPIU_INT,next->grstarts,1,MPIU_INT,PetscObjectComm((PetscObject)dm));
117: next = next->next;
118: }
119: com->setup = PETSC_TRUE;
120: return(0);
121: }
123: /* ----------------------------------------------------------------------------------*/
125: /*@
126: DMCompositeGetNumberDM - Get's the number of DM objects in the DMComposite
127: representation.
129: Not Collective
131: Input Parameter:
132: . dm - the packer object
134: Output Parameter:
135: . nDM - the number of DMs
137: Level: beginner
139: @*/
140: PetscErrorCode DMCompositeGetNumberDM(DM dm,PetscInt *nDM)
141: {
142: DM_Composite *com = (DM_Composite*)dm->data;
143: PetscBool flg;
148: PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
149: if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
150: *nDM = com->nDM;
151: return(0);
152: }
155: /*@C
156: DMCompositeGetAccess - Allows one to access the individual packed vectors in their global
157: representation.
159: Collective on DMComposite
161: Input Parameters:
162: + dm - the packer object
163: - gvec - the global vector
165: Output Parameters:
166: . Vec* ... - the packed parallel vectors, NULL for those that are not needed
168: Notes:
169: Use DMCompositeRestoreAccess() to return the vectors when you no longer need them
171: Fortran Notes:
173: Fortran callers must use numbered versions of this routine, e.g., DMCompositeGetAccess4(dm,gvec,vec1,vec2,vec3,vec4)
174: or use the alternative interface DMCompositeGetAccessArray().
176: Level: advanced
178: .seealso: DMCompositeGetEntries(), DMCompositeScatter()
179: @*/
180: PetscErrorCode DMCompositeGetAccess(DM dm,Vec gvec,...)
181: {
182: va_list Argp;
183: PetscErrorCode ierr;
184: struct DMCompositeLink *next;
185: DM_Composite *com = (DM_Composite*)dm->data;
186: PetscInt readonly;
187: PetscBool flg;
192: PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
193: if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
194: next = com->next;
195: if (!com->setup) {
196: DMSetUp(dm);
197: }
199: VecLockGet(gvec,&readonly);
200: /* loop over packed objects, handling one at at time */
201: va_start(Argp,gvec);
202: while (next) {
203: Vec *vec;
204: vec = va_arg(Argp, Vec*);
205: if (vec) {
206: DMGetGlobalVector(next->dm,vec);
207: if (readonly) {
208: const PetscScalar *array;
209: VecGetArrayRead(gvec,&array);
210: VecPlaceArray(*vec,array+next->rstart);
211: VecLockPush(*vec);
212: VecRestoreArrayRead(gvec,&array);
213: } else {
214: PetscScalar *array;
215: VecGetArray(gvec,&array);
216: VecPlaceArray(*vec,array+next->rstart);
217: VecRestoreArray(gvec,&array);
218: }
219: }
220: next = next->next;
221: }
222: va_end(Argp);
223: return(0);
224: }
226: /*@C
227: DMCompositeGetAccessArray - Allows one to access the individual packed vectors in their global
228: representation.
230: Collective on DMComposite
232: Input Parameters:
233: + dm - the packer object
234: . pvec - packed vector
235: . nwanted - number of vectors wanted
236: - wanted - sorted array of vectors wanted, or NULL to get all vectors
238: Output Parameters:
239: . vecs - array of requested global vectors (must be allocated)
241: Notes:
242: Use DMCompositeRestoreAccessArray() to return the vectors when you no longer need them
244: Level: advanced
246: .seealso: DMCompositeGetAccess(), DMCompositeGetEntries(), DMCompositeScatter(), DMCompositeGather()
247: @*/
248: PetscErrorCode DMCompositeGetAccessArray(DM dm,Vec pvec,PetscInt nwanted,const PetscInt *wanted,Vec *vecs)
249: {
250: PetscErrorCode ierr;
251: struct DMCompositeLink *link;
252: PetscInt i,wnum;
253: DM_Composite *com = (DM_Composite*)dm->data;
254: PetscInt readonly;
255: PetscBool flg;
260: PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
261: if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
262: if (!com->setup) {
263: DMSetUp(dm);
264: }
266: VecLockGet(pvec,&readonly);
267: for (i=0,wnum=0,link=com->next; link && wnum<nwanted; i++,link=link->next) {
268: if (!wanted || i == wanted[wnum]) {
269: Vec v;
270: DMGetGlobalVector(link->dm,&v);
271: if (readonly) {
272: const PetscScalar *array;
273: VecGetArrayRead(pvec,&array);
274: VecPlaceArray(v,array+link->rstart);
275: VecLockPush(v);
276: VecRestoreArrayRead(pvec,&array);
277: } else {
278: PetscScalar *array;
279: VecGetArray(pvec,&array);
280: VecPlaceArray(v,array+link->rstart);
281: VecRestoreArray(pvec,&array);
282: }
283: vecs[wnum++] = v;
284: }
285: }
286: return(0);
287: }
289: /*@C
290: DMCompositeGetLocalAccessArray - Allows one to access the individual
291: packed vectors in their local representation.
293: Collective on DMComposite.
295: Input Parameters:
296: + dm - the packer object
297: . pvec - packed vector
298: . nwanted - number of vectors wanted
299: - wanted - sorted array of vectors wanted, or NULL to get all vectors
301: Output Parameters:
302: . vecs - array of requested local vectors (must be allocated)
304: Notes:
305: Use DMCompositeRestoreLocalAccessArray() to return the vectors
306: when you no longer need them.
308: Level: advanced
310: .seealso: DMCompositeRestoreLocalAccessArray(), DMCompositeGetAccess(),
311: DMCompositeGetEntries(), DMCompositeScatter(), DMCompositeGather()
312: @*/
313: PetscErrorCode DMCompositeGetLocalAccessArray(DM dm,Vec pvec,PetscInt nwanted,const PetscInt *wanted,Vec *vecs)
314: {
315: PetscErrorCode ierr;
316: struct DMCompositeLink *link;
317: PetscInt i,wnum;
318: DM_Composite *com = (DM_Composite*)dm->data;
319: PetscInt readonly;
320: PetscInt nlocal = 0;
321: PetscBool flg;
326: PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
327: if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
328: if (!com->setup) {
329: DMSetUp(dm);
330: }
332: VecLockGet(pvec,&readonly);
333: for (i=0,wnum=0,link=com->next; link && wnum<nwanted; i++,link=link->next) {
334: if (!wanted || i == wanted[wnum]) {
335: Vec v;
336: DMGetLocalVector(link->dm,&v);
337: if (readonly) {
338: const PetscScalar *array;
339: VecGetArrayRead(pvec,&array);
340: VecPlaceArray(v,array+nlocal);
341: VecLockPush(v);
342: VecRestoreArrayRead(pvec,&array);
343: } else {
344: PetscScalar *array;
345: VecGetArray(pvec,&array);
346: VecPlaceArray(v,array+nlocal);
347: VecRestoreArray(pvec,&array);
348: }
349: vecs[wnum++] = v;
350: }
352: nlocal += link->nlocal;
353: }
355: return(0);
356: }
358: /*@C
359: DMCompositeRestoreAccess - Returns the vectors obtained with DMCompositeGetAccess()
360: representation.
362: Collective on DMComposite
364: Input Parameters:
365: + dm - the packer object
366: . gvec - the global vector
367: - Vec* ... - the individual parallel vectors, NULL for those that are not needed
369: Level: advanced
371: .seealso DMCompositeAddDM(), DMCreateGlobalVector(),
372: DMCompositeGather(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeScatter(),
373: DMCompositeRestoreAccess(), DMCompositeGetAccess()
375: @*/
376: PetscErrorCode DMCompositeRestoreAccess(DM dm,Vec gvec,...)
377: {
378: va_list Argp;
379: PetscErrorCode ierr;
380: struct DMCompositeLink *next;
381: DM_Composite *com = (DM_Composite*)dm->data;
382: PetscInt readonly;
383: PetscBool flg;
388: PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
389: if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
390: next = com->next;
391: if (!com->setup) {
392: DMSetUp(dm);
393: }
395: VecLockGet(gvec,&readonly);
396: /* loop over packed objects, handling one at at time */
397: va_start(Argp,gvec);
398: while (next) {
399: Vec *vec;
400: vec = va_arg(Argp, Vec*);
401: if (vec) {
402: VecResetArray(*vec);
403: if (readonly) {
404: VecLockPop(*vec);
405: }
406: DMRestoreGlobalVector(next->dm,vec);
407: }
408: next = next->next;
409: }
410: va_end(Argp);
411: return(0);
412: }
414: /*@C
415: DMCompositeRestoreAccessArray - Returns the vectors obtained with DMCompositeGetAccessArray()
417: Collective on DMComposite
419: Input Parameters:
420: + dm - the packer object
421: . pvec - packed vector
422: . nwanted - number of vectors wanted
423: . wanted - sorted array of vectors wanted, or NULL to get all vectors
424: - vecs - array of global vectors to return
426: Level: advanced
428: .seealso: DMCompositeRestoreAccess(), DMCompositeRestoreEntries(), DMCompositeScatter(), DMCompositeGather()
429: @*/
430: PetscErrorCode DMCompositeRestoreAccessArray(DM dm,Vec pvec,PetscInt nwanted,const PetscInt *wanted,Vec *vecs)
431: {
432: PetscErrorCode ierr;
433: struct DMCompositeLink *link;
434: PetscInt i,wnum;
435: DM_Composite *com = (DM_Composite*)dm->data;
436: PetscInt readonly;
437: PetscBool flg;
442: PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
443: if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
444: if (!com->setup) {
445: DMSetUp(dm);
446: }
448: VecLockGet(pvec,&readonly);
449: for (i=0,wnum=0,link=com->next; link && wnum<nwanted; i++,link=link->next) {
450: if (!wanted || i == wanted[wnum]) {
451: VecResetArray(vecs[wnum]);
452: if (readonly) {
453: VecLockPop(vecs[wnum]);
454: }
455: DMRestoreGlobalVector(link->dm,&vecs[wnum]);
456: wnum++;
457: }
458: }
459: return(0);
460: }
462: /*@C
463: DMCompositeRestoreLocalAccessArray - Returns the vectors obtained with DMCompositeGetLocalAccessArray().
465: Collective on DMComposite.
467: Input Parameters:
468: + dm - the packer object
469: . pvec - packed vector
470: . nwanted - number of vectors wanted
471: . wanted - sorted array of vectors wanted, or NULL to restore all vectors
472: - vecs - array of local vectors to return
474: Level: advanced
476: Notes:
477: nwanted and wanted must match the values given to DMCompositeGetLocalAccessArray()
478: otherwise the call will fail.
480: .seealso: DMCompositeGetLocalAccessArray(), DMCompositeRestoreAccessArray(),
481: DMCompositeRestoreAccess(), DMCompositeRestoreEntries(),
482: DMCompositeScatter(), DMCompositeGather()
483: @*/
484: PetscErrorCode DMCompositeRestoreLocalAccessArray(DM dm,Vec pvec,PetscInt nwanted,const PetscInt *wanted,Vec *vecs)
485: {
486: PetscErrorCode ierr;
487: struct DMCompositeLink *link;
488: PetscInt i,wnum;
489: DM_Composite *com = (DM_Composite*)dm->data;
490: PetscInt readonly;
491: PetscBool flg;
496: PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
497: if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
498: if (!com->setup) {
499: DMSetUp(dm);
500: }
502: VecLockGet(pvec,&readonly);
503: for (i=0,wnum=0,link=com->next; link && wnum<nwanted; i++,link=link->next) {
504: if (!wanted || i == wanted[wnum]) {
505: VecResetArray(vecs[wnum]);
506: if (readonly) {
507: VecLockPop(vecs[wnum]);
508: }
509: DMRestoreLocalVector(link->dm,&vecs[wnum]);
510: wnum++;
511: }
512: }
513: return(0);
514: }
516: /*@C
517: DMCompositeScatter - Scatters from a global packed vector into its individual local vectors
519: Collective on DMComposite
521: Input Parameters:
522: + dm - the packer object
523: . gvec - the global vector
524: - Vec ... - the individual sequential vectors, NULL for those that are not needed
526: Level: advanced
528: Notes:
529: DMCompositeScatterArray() is a non-variadic alternative that is often more convenient for library callers and is
530: accessible from Fortran.
532: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(),
533: DMCompositeGather(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
534: DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(), DMCompositeGetEntries()
535: DMCompositeScatterArray()
537: @*/
538: PetscErrorCode DMCompositeScatter(DM dm,Vec gvec,...)
539: {
540: va_list Argp;
541: PetscErrorCode ierr;
542: struct DMCompositeLink *next;
543: PetscInt cnt;
544: DM_Composite *com = (DM_Composite*)dm->data;
545: PetscBool flg;
550: PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
551: if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
552: if (!com->setup) {
553: DMSetUp(dm);
554: }
556: /* loop over packed objects, handling one at at time */
557: va_start(Argp,gvec);
558: for (cnt=3,next=com->next; next; cnt++,next=next->next) {
559: Vec local;
560: local = va_arg(Argp, Vec);
561: if (local) {
562: Vec global;
563: const PetscScalar *array;
565: DMGetGlobalVector(next->dm,&global);
566: VecGetArrayRead(gvec,&array);
567: VecPlaceArray(global,array+next->rstart);
568: DMGlobalToLocalBegin(next->dm,global,INSERT_VALUES,local);
569: DMGlobalToLocalEnd(next->dm,global,INSERT_VALUES,local);
570: VecRestoreArrayRead(gvec,&array);
571: VecResetArray(global);
572: DMRestoreGlobalVector(next->dm,&global);
573: }
574: }
575: va_end(Argp);
576: return(0);
577: }
579: /*@
580: DMCompositeScatterArray - Scatters from a global packed vector into its individual local vectors
582: Collective on DMComposite
584: Input Parameters:
585: + dm - the packer object
586: . gvec - the global vector
587: . lvecs - array of local vectors, NULL for any that are not needed
589: Level: advanced
591: Note:
592: This is a non-variadic alternative to DMCompositeScatter()
594: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector()
595: DMCompositeGather(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
596: DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(), DMCompositeGetEntries()
598: @*/
599: PetscErrorCode DMCompositeScatterArray(DM dm,Vec gvec,Vec *lvecs)
600: {
601: PetscErrorCode ierr;
602: struct DMCompositeLink *next;
603: PetscInt i;
604: DM_Composite *com = (DM_Composite*)dm->data;
605: PetscBool flg;
610: PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
611: if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
612: if (!com->setup) {
613: DMSetUp(dm);
614: }
616: /* loop over packed objects, handling one at at time */
617: for (i=0,next=com->next; next; next=next->next,i++) {
618: if (lvecs[i]) {
619: Vec global;
620: const PetscScalar *array;
622: DMGetGlobalVector(next->dm,&global);
623: VecGetArrayRead(gvec,&array);
624: VecPlaceArray(global,(PetscScalar*)array+next->rstart);
625: DMGlobalToLocalBegin(next->dm,global,INSERT_VALUES,lvecs[i]);
626: DMGlobalToLocalEnd(next->dm,global,INSERT_VALUES,lvecs[i]);
627: VecRestoreArrayRead(gvec,&array);
628: VecResetArray(global);
629: DMRestoreGlobalVector(next->dm,&global);
630: }
631: }
632: return(0);
633: }
635: /*@C
636: DMCompositeGather - Gathers into a global packed vector from its individual local vectors
638: Collective on DMComposite
640: Input Parameter:
641: + dm - the packer object
642: . gvec - the global vector
643: . imode - INSERT_VALUES or ADD_VALUES
644: - Vec ... - the individual sequential vectors, NULL for any that are not needed
646: Level: advanced
648: Not available from Fortran, Fortran users can use DMCompositeGatherArray()
650: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(),
651: DMCompositeScatter(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
652: DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(), DMCompositeGetEntries()
654: @*/
655: PetscErrorCode DMCompositeGather(DM dm,InsertMode imode,Vec gvec,...)
656: {
657: va_list Argp;
658: PetscErrorCode ierr;
659: struct DMCompositeLink *next;
660: DM_Composite *com = (DM_Composite*)dm->data;
661: PetscInt cnt;
662: PetscBool flg;
667: PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
668: if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
669: if (!com->setup) {
670: DMSetUp(dm);
671: }
673: /* loop over packed objects, handling one at at time */
674: va_start(Argp,gvec);
675: for (cnt=3,next=com->next; next; cnt++,next=next->next) {
676: Vec local;
677: local = va_arg(Argp, Vec);
678: if (local) {
679: PetscScalar *array;
680: Vec global;
682: DMGetGlobalVector(next->dm,&global);
683: VecGetArray(gvec,&array);
684: VecPlaceArray(global,array+next->rstart);
685: DMLocalToGlobalBegin(next->dm,local,imode,global);
686: DMLocalToGlobalEnd(next->dm,local,imode,global);
687: VecRestoreArray(gvec,&array);
688: VecResetArray(global);
689: DMRestoreGlobalVector(next->dm,&global);
690: }
691: }
692: va_end(Argp);
693: return(0);
694: }
696: /*@
697: DMCompositeGatherArray - Gathers into a global packed vector from its individual local vectors
699: Collective on DMComposite
701: Input Parameter:
702: + dm - the packer object
703: . gvec - the global vector
704: . imode - INSERT_VALUES or ADD_VALUES
705: - lvecs - the individual sequential vectors, NULL for any that are not needed
707: Level: advanced
709: Notes:
710: This is a non-variadic alternative to DMCompositeGather().
712: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(),
713: DMCompositeScatter(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
714: DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(), DMCompositeGetEntries(),
715: @*/
716: PetscErrorCode DMCompositeGatherArray(DM dm,InsertMode imode,Vec gvec,Vec *lvecs)
717: {
718: PetscErrorCode ierr;
719: struct DMCompositeLink *next;
720: DM_Composite *com = (DM_Composite*)dm->data;
721: PetscInt i;
722: PetscBool flg;
727: PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
728: if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
729: if (!com->setup) {
730: DMSetUp(dm);
731: }
733: /* loop over packed objects, handling one at at time */
734: for (next=com->next,i=0; next; next=next->next,i++) {
735: if (lvecs[i]) {
736: PetscScalar *array;
737: Vec global;
739: DMGetGlobalVector(next->dm,&global);
740: VecGetArray(gvec,&array);
741: VecPlaceArray(global,array+next->rstart);
742: DMLocalToGlobalBegin(next->dm,lvecs[i],imode,global);
743: DMLocalToGlobalEnd(next->dm,lvecs[i],imode,global);
744: VecRestoreArray(gvec,&array);
745: VecResetArray(global);
746: DMRestoreGlobalVector(next->dm,&global);
747: }
748: }
749: return(0);
750: }
752: /*@
753: DMCompositeAddDM - adds a DM vector to a DMComposite
755: Collective on DMComposite
757: Input Parameter:
758: + dmc - the DMComposite (packer) object
759: - dm - the DM object
761: Level: advanced
763: .seealso DMDestroy(), DMCompositeGather(), DMCompositeAddDM(), DMCreateGlobalVector(),
764: DMCompositeScatter(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
765: DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(), DMCompositeGetEntries()
767: @*/
768: PetscErrorCode DMCompositeAddDM(DM dmc,DM dm)
769: {
770: PetscErrorCode ierr;
771: PetscInt n,nlocal;
772: struct DMCompositeLink *mine,*next;
773: Vec global,local;
774: DM_Composite *com = (DM_Composite*)dmc->data;
775: PetscBool flg;
780: PetscObjectTypeCompare((PetscObject)dmc,DMCOMPOSITE,&flg);
781: if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
782: next = com->next;
783: if (com->setup) SETERRQ(PetscObjectComm((PetscObject)dmc),PETSC_ERR_ARG_WRONGSTATE,"Cannot add a DM once you have used the DMComposite");
785: /* create new link */
786: PetscNew(&mine);
787: PetscObjectReference((PetscObject)dm);
788: DMGetGlobalVector(dm,&global);
789: VecGetLocalSize(global,&n);
790: DMRestoreGlobalVector(dm,&global);
791: DMGetLocalVector(dm,&local);
792: VecGetSize(local,&nlocal);
793: DMRestoreLocalVector(dm,&local);
795: mine->n = n;
796: mine->nlocal = nlocal;
797: mine->dm = dm;
798: mine->next = NULL;
799: com->n += n;
800: com->nghost += nlocal;
802: /* add to end of list */
803: if (!next) com->next = mine;
804: else {
805: while (next->next) next = next->next;
806: next->next = mine;
807: }
808: com->nDM++;
809: com->nmine++;
810: return(0);
811: }
813: #include <petscdraw.h>
814: PETSC_EXTERN PetscErrorCode VecView_MPI(Vec,PetscViewer);
815: PetscErrorCode VecView_DMComposite(Vec gvec,PetscViewer viewer)
816: {
817: DM dm;
818: PetscErrorCode ierr;
819: struct DMCompositeLink *next;
820: PetscBool isdraw;
821: DM_Composite *com;
824: VecGetDM(gvec, &dm);
825: if (!dm) SETERRQ(PetscObjectComm((PetscObject)gvec),PETSC_ERR_ARG_WRONG,"Vector not generated from a DMComposite");
826: com = (DM_Composite*)dm->data;
827: next = com->next;
829: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
830: if (!isdraw) {
831: /* do I really want to call this? */
832: VecView_MPI(gvec,viewer);
833: } else {
834: PetscInt cnt = 0;
836: /* loop over packed objects, handling one at at time */
837: while (next) {
838: Vec vec;
839: const PetscScalar *array;
840: PetscInt bs;
842: /* Should use VecGetSubVector() eventually, but would need to forward the DM for that to work */
843: DMGetGlobalVector(next->dm,&vec);
844: VecGetArrayRead(gvec,&array);
845: VecPlaceArray(vec,(PetscScalar*)array+next->rstart);
846: VecRestoreArrayRead(gvec,&array);
847: VecView(vec,viewer);
848: VecResetArray(vec);
849: VecGetBlockSize(vec,&bs);
850: DMRestoreGlobalVector(next->dm,&vec);
851: PetscViewerDrawBaseAdd(viewer,bs);
852: cnt += bs;
853: next = next->next;
854: }
855: PetscViewerDrawBaseAdd(viewer,-cnt);
856: }
857: return(0);
858: }
860: PetscErrorCode DMCreateGlobalVector_Composite(DM dm,Vec *gvec)
861: {
863: DM_Composite *com = (DM_Composite*)dm->data;
867: DMSetUp(dm);
868: VecCreateMPI(PetscObjectComm((PetscObject)dm),com->n,com->N,gvec);
869: VecSetDM(*gvec, dm);
870: VecSetOperation(*gvec,VECOP_VIEW,(void (*)(void))VecView_DMComposite);
871: return(0);
872: }
874: PetscErrorCode DMCreateLocalVector_Composite(DM dm,Vec *lvec)
875: {
877: DM_Composite *com = (DM_Composite*)dm->data;
881: if (!com->setup) {
882: DMSetUp(dm);
883: }
884: VecCreateSeq(PETSC_COMM_SELF,com->nghost,lvec);
885: VecSetDM(*lvec, dm);
886: return(0);
887: }
889: /*@C
890: DMCompositeGetISLocalToGlobalMappings - gets an ISLocalToGlobalMapping for each DM in the DMComposite, maps to the composite global space
892: Collective on DM
894: Input Parameter:
895: . dm - the packer object
897: Output Parameters:
898: . ltogs - the individual mappings for each packed vector. Note that this includes
899: all the ghost points that individual ghosted DMDA's may have.
901: Level: advanced
903: Notes:
904: Each entry of ltogs should be destroyed with ISLocalToGlobalMappingDestroy(), the ltogs array should be freed with PetscFree().
906: Not available from Fortran
908: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(),
909: DMCompositeGather(), DMCompositeCreate(), DMCompositeGetAccess(), DMCompositeScatter(),
910: DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(),DMCompositeGetEntries()
912: @*/
913: PetscErrorCode DMCompositeGetISLocalToGlobalMappings(DM dm,ISLocalToGlobalMapping **ltogs)
914: {
915: PetscErrorCode ierr;
916: PetscInt i,*idx,n,cnt;
917: struct DMCompositeLink *next;
918: PetscMPIInt rank;
919: DM_Composite *com = (DM_Composite*)dm->data;
920: PetscBool flg;
924: PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
925: if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
926: DMSetUp(dm);
927: PetscMalloc1(com->nDM,ltogs);
928: next = com->next;
929: MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);
931: /* loop over packed objects, handling one at at time */
932: cnt = 0;
933: while (next) {
934: ISLocalToGlobalMapping ltog;
935: PetscMPIInt size;
936: const PetscInt *suboff,*indices;
937: Vec global;
939: /* Get sub-DM global indices for each local dof */
940: DMGetLocalToGlobalMapping(next->dm,<og);
941: ISLocalToGlobalMappingGetSize(ltog,&n);
942: ISLocalToGlobalMappingGetIndices(ltog,&indices);
943: PetscMalloc1(n,&idx);
945: /* Get the offsets for the sub-DM global vector */
946: DMGetGlobalVector(next->dm,&global);
947: VecGetOwnershipRanges(global,&suboff);
948: MPI_Comm_size(PetscObjectComm((PetscObject)global),&size);
950: /* Shift the sub-DM definition of the global space to the composite global space */
951: for (i=0; i<n; i++) {
952: PetscInt subi = indices[i],lo = 0,hi = size,t;
953: /* Binary search to find which rank owns subi */
954: while (hi-lo > 1) {
955: t = lo + (hi-lo)/2;
956: if (suboff[t] > subi) hi = t;
957: else lo = t;
958: }
959: idx[i] = subi - suboff[lo] + next->grstarts[lo];
960: }
961: ISLocalToGlobalMappingRestoreIndices(ltog,&indices);
962: ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)dm),1,n,idx,PETSC_OWN_POINTER,&(*ltogs)[cnt]);
963: DMRestoreGlobalVector(next->dm,&global);
964: next = next->next;
965: cnt++;
966: }
967: return(0);
968: }
970: /*@C
971: DMCompositeGetLocalISs - Gets index sets for each component of a composite local vector
973: Not Collective
975: Input Arguments:
976: . dm - composite DM
978: Output Arguments:
979: . is - array of serial index sets for each each component of the DMComposite
981: Level: intermediate
983: Notes:
984: At present, a composite local vector does not normally exist. This function is used to provide index sets for
985: MatGetLocalSubMatrix(). In the future, the scatters for each entry in the DMComposite may be be merged into a single
986: scatter to a composite local vector. The user should not typically need to know which is being done.
988: To get the composite global indices at all local points (including ghosts), use DMCompositeGetISLocalToGlobalMappings().
990: To get index sets for pieces of the composite global vector, use DMCompositeGetGlobalISs().
992: Each returned IS should be destroyed with ISDestroy(), the array should be freed with PetscFree().
994: Not available from Fortran
996: .seealso: DMCompositeGetGlobalISs(), DMCompositeGetISLocalToGlobalMappings(), MatGetLocalSubMatrix(), MatCreateLocalRef()
997: @*/
998: PetscErrorCode DMCompositeGetLocalISs(DM dm,IS **is)
999: {
1000: PetscErrorCode ierr;
1001: DM_Composite *com = (DM_Composite*)dm->data;
1002: struct DMCompositeLink *link;
1003: PetscInt cnt,start;
1004: PetscBool flg;
1009: PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
1010: if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
1011: PetscMalloc1(com->nmine,is);
1012: for (cnt=0,start=0,link=com->next; link; start+=link->nlocal,cnt++,link=link->next) {
1013: PetscInt bs;
1014: ISCreateStride(PETSC_COMM_SELF,link->nlocal,start,1,&(*is)[cnt]);
1015: DMGetBlockSize(link->dm,&bs);
1016: ISSetBlockSize((*is)[cnt],bs);
1017: }
1018: return(0);
1019: }
1021: /*@C
1022: DMCompositeGetGlobalISs - Gets the index sets for each composed object
1024: Collective on DMComposite
1026: Input Parameter:
1027: . dm - the packer object
1029: Output Parameters:
1030: . is - the array of index sets
1032: Level: advanced
1034: Notes:
1035: The is entries should be destroyed with ISDestroy(), the is array should be freed with PetscFree()
1037: These could be used to extract a subset of vector entries for a "multi-physics" preconditioner
1039: Use DMCompositeGetLocalISs() for index sets in the packed local numbering, and
1040: DMCompositeGetISLocalToGlobalMappings() for to map local sub-DM (including ghost) indices to packed global
1041: indices.
1043: Fortran Notes:
1045: The output argument 'is' must be an allocated array of sufficient length, which can be learned using DMCompositeGetNumberDM().
1047: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(),
1048: DMCompositeGather(), DMCompositeCreate(), DMCompositeGetAccess(), DMCompositeScatter(),
1049: DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(),DMCompositeGetEntries()
1051: @*/
1052: PetscErrorCode DMCompositeGetGlobalISs(DM dm,IS *is[])
1053: {
1054: PetscErrorCode ierr;
1055: PetscInt cnt = 0;
1056: struct DMCompositeLink *next;
1057: PetscMPIInt rank;
1058: DM_Composite *com = (DM_Composite*)dm->data;
1059: PetscBool flg;
1063: PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
1064: if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
1065: if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Must call DMSetUp() before");
1066: PetscMalloc1(com->nDM,is);
1067: next = com->next;
1068: MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);
1070: /* loop over packed objects, handling one at at time */
1071: while (next) {
1072: ISCreateStride(PetscObjectComm((PetscObject)dm),next->n,next->grstart,1,&(*is)[cnt]);
1073: if (dm->prob) {
1074: MatNullSpace space;
1075: Mat pmat;
1076: PetscObject disc;
1077: PetscInt Nf;
1079: PetscDSGetNumFields(dm->prob, &Nf);
1080: if (cnt < Nf) {
1081: PetscDSGetDiscretization(dm->prob, cnt, &disc);
1082: PetscObjectQuery(disc, "nullspace", (PetscObject*) &space);
1083: if (space) {PetscObjectCompose((PetscObject) (*is)[cnt], "nullspace", (PetscObject) space);}
1084: PetscObjectQuery(disc, "nearnullspace", (PetscObject*) &space);
1085: if (space) {PetscObjectCompose((PetscObject) (*is)[cnt], "nearnullspace", (PetscObject) space);}
1086: PetscObjectQuery(disc, "pmat", (PetscObject*) &pmat);
1087: if (pmat) {PetscObjectCompose((PetscObject) (*is)[cnt], "pmat", (PetscObject) pmat);}
1088: }
1089: }
1090: cnt++;
1091: next = next->next;
1092: }
1093: return(0);
1094: }
1096: PetscErrorCode DMCreateFieldIS_Composite(DM dm, PetscInt *numFields,char ***fieldNames, IS **fields)
1097: {
1098: PetscInt nDM;
1099: DM *dms;
1100: PetscInt i;
1104: DMCompositeGetNumberDM(dm, &nDM);
1105: if (numFields) *numFields = nDM;
1106: DMCompositeGetGlobalISs(dm, fields);
1107: if (fieldNames) {
1108: PetscMalloc1(nDM, &dms);
1109: PetscMalloc1(nDM, fieldNames);
1110: DMCompositeGetEntriesArray(dm, dms);
1111: for (i=0; i<nDM; i++) {
1112: char buf[256];
1113: const char *splitname;
1115: /* Split naming precedence: object name, prefix, number */
1116: splitname = ((PetscObject) dm)->name;
1117: if (!splitname) {
1118: PetscObjectGetOptionsPrefix((PetscObject)dms[i],&splitname);
1119: if (splitname) {
1120: size_t len;
1121: PetscStrncpy(buf,splitname,sizeof(buf));
1122: buf[sizeof(buf) - 1] = 0;
1123: PetscStrlen(buf,&len);
1124: if (buf[len-1] == '_') buf[len-1] = 0; /* Remove trailing underscore if it was used */
1125: splitname = buf;
1126: }
1127: }
1128: if (!splitname) {
1129: PetscSNPrintf(buf,sizeof(buf),"%D",i);
1130: splitname = buf;
1131: }
1132: PetscStrallocpy(splitname,&(*fieldNames)[i]);
1133: }
1134: PetscFree(dms);
1135: }
1136: return(0);
1137: }
1139: /*
1140: This could take over from DMCreateFieldIS(), as it is more general,
1141: making DMCreateFieldIS() a special case -- calling with dmlist == NULL;
1142: At this point it's probably best to be less intrusive, however.
1143: */
1144: PetscErrorCode DMCreateFieldDecomposition_Composite(DM dm, PetscInt *len,char ***namelist, IS **islist, DM **dmlist)
1145: {
1146: PetscInt nDM;
1147: PetscInt i;
1151: DMCreateFieldIS_Composite(dm, len, namelist, islist);
1152: if (dmlist) {
1153: DMCompositeGetNumberDM(dm, &nDM);
1154: PetscMalloc1(nDM, dmlist);
1155: DMCompositeGetEntriesArray(dm, *dmlist);
1156: for (i=0; i<nDM; i++) {
1157: PetscObjectReference((PetscObject)((*dmlist)[i]));
1158: }
1159: }
1160: return(0);
1161: }
1165: /* -------------------------------------------------------------------------------------*/
1166: /*@C
1167: DMCompositeGetLocalVectors - Gets local vectors for each part of a DMComposite.
1168: Use DMCompositeRestoreLocalVectors() to return them.
1170: Not Collective
1172: Input Parameter:
1173: . dm - the packer object
1175: Output Parameter:
1176: . Vec ... - the individual sequential Vecs
1178: Level: advanced
1180: Not available from Fortran
1182: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(),
1183: DMCompositeGather(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
1184: DMCompositeRestoreLocalVectors(), DMCompositeScatter(), DMCompositeGetEntries()
1186: @*/
1187: PetscErrorCode DMCompositeGetLocalVectors(DM dm,...)
1188: {
1189: va_list Argp;
1190: PetscErrorCode ierr;
1191: struct DMCompositeLink *next;
1192: DM_Composite *com = (DM_Composite*)dm->data;
1193: PetscBool flg;
1197: PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
1198: if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
1199: next = com->next;
1200: /* loop over packed objects, handling one at at time */
1201: va_start(Argp,dm);
1202: while (next) {
1203: Vec *vec;
1204: vec = va_arg(Argp, Vec*);
1205: if (vec) {DMGetLocalVector(next->dm,vec);}
1206: next = next->next;
1207: }
1208: va_end(Argp);
1209: return(0);
1210: }
1212: /*@C
1213: DMCompositeRestoreLocalVectors - Restores local vectors for each part of a DMComposite.
1215: Not Collective
1217: Input Parameter:
1218: . dm - the packer object
1220: Output Parameter:
1221: . Vec ... - the individual sequential Vecs
1223: Level: advanced
1225: Not available from Fortran
1227: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(),
1228: DMCompositeGather(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
1229: DMCompositeGetLocalVectors(), DMCompositeScatter(), DMCompositeGetEntries()
1231: @*/
1232: PetscErrorCode DMCompositeRestoreLocalVectors(DM dm,...)
1233: {
1234: va_list Argp;
1235: PetscErrorCode ierr;
1236: struct DMCompositeLink *next;
1237: DM_Composite *com = (DM_Composite*)dm->data;
1238: PetscBool flg;
1242: PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
1243: if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
1244: next = com->next;
1245: /* loop over packed objects, handling one at at time */
1246: va_start(Argp,dm);
1247: while (next) {
1248: Vec *vec;
1249: vec = va_arg(Argp, Vec*);
1250: if (vec) {DMRestoreLocalVector(next->dm,vec);}
1251: next = next->next;
1252: }
1253: va_end(Argp);
1254: return(0);
1255: }
1257: /* -------------------------------------------------------------------------------------*/
1258: /*@C
1259: DMCompositeGetEntries - Gets the DM for each entry in a DMComposite.
1261: Not Collective
1263: Input Parameter:
1264: . dm - the packer object
1266: Output Parameter:
1267: . DM ... - the individual entries (DMs)
1269: Level: advanced
1271: Fortran Notes:
1272: Available as DMCompositeGetEntries() for one output DM, DMCompositeGetEntries2() for 2, etc
1274: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(), DMCompositeGetEntriesArray()
1275: DMCompositeGather(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
1276: DMCompositeRestoreLocalVectors(), DMCompositeGetLocalVectors(), DMCompositeScatter(),
1277: DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors()
1279: @*/
1280: PetscErrorCode DMCompositeGetEntries(DM dm,...)
1281: {
1282: va_list Argp;
1283: struct DMCompositeLink *next;
1284: DM_Composite *com = (DM_Composite*)dm->data;
1285: PetscBool flg;
1286: PetscErrorCode ierr;
1290: PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
1291: if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
1292: next = com->next;
1293: /* loop over packed objects, handling one at at time */
1294: va_start(Argp,dm);
1295: while (next) {
1296: DM *dmn;
1297: dmn = va_arg(Argp, DM*);
1298: if (dmn) *dmn = next->dm;
1299: next = next->next;
1300: }
1301: va_end(Argp);
1302: return(0);
1303: }
1305: /*@C
1306: DMCompositeGetEntriesArray - Gets the DM for each entry in a DMComposite.
1308: Not Collective
1310: Input Parameter:
1311: . dm - the packer object
1313: Output Parameter:
1314: . dms - array of sufficient length (see DMCompositeGetNumberDM()) to hold the individual DMs
1316: Level: advanced
1318: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(), DMCompositeGetEntries()
1319: DMCompositeGather(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
1320: DMCompositeRestoreLocalVectors(), DMCompositeGetLocalVectors(), DMCompositeScatter(),
1321: DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors()
1323: @*/
1324: PetscErrorCode DMCompositeGetEntriesArray(DM dm,DM dms[])
1325: {
1326: struct DMCompositeLink *next;
1327: DM_Composite *com = (DM_Composite*)dm->data;
1328: PetscInt i;
1329: PetscBool flg;
1330: PetscErrorCode ierr;
1334: PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
1335: if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
1336: /* loop over packed objects, handling one at at time */
1337: for (next=com->next,i=0; next; next=next->next,i++) dms[i] = next->dm;
1338: return(0);
1339: }
1341: typedef struct {
1342: DM dm;
1343: PetscViewer *subv;
1344: Vec *vecs;
1345: } GLVisViewerCtx;
1347: static PetscErrorCode DestroyGLVisViewerCtx_Private(void *vctx)
1348: {
1349: GLVisViewerCtx *ctx = (GLVisViewerCtx*)vctx;
1350: PetscInt i,n;
1354: DMCompositeGetNumberDM(ctx->dm,&n);
1355: for (i = 0; i < n; i++) {
1356: PetscViewerDestroy(&ctx->subv[i]);
1357: }
1358: PetscFree2(ctx->subv,ctx->vecs);
1359: DMDestroy(&ctx->dm);
1360: PetscFree(ctx);
1361: return(0);
1362: }
1364: static PetscErrorCode DMCompositeSampleGLVisFields_Private(PetscObject oX, PetscInt nf, PetscObject oXfield[], void *vctx)
1365: {
1366: Vec X = (Vec)oX;
1367: GLVisViewerCtx *ctx = (GLVisViewerCtx*)vctx;
1368: PetscInt i,n,cumf;
1372: DMCompositeGetNumberDM(ctx->dm,&n);
1373: DMCompositeGetAccessArray(ctx->dm,X,n,NULL,ctx->vecs);
1374: for (i = 0, cumf = 0; i < n; i++) {
1375: PetscErrorCode (*g2l)(PetscObject,PetscInt,PetscObject[],void*);
1376: void *fctx;
1377: PetscInt nfi;
1379: PetscViewerGLVisGetFields_Private(ctx->subv[i],&nfi,NULL,NULL,&g2l,NULL,&fctx);
1380: if (!nfi) continue;
1381: if (g2l) {
1382: (*g2l)((PetscObject)ctx->vecs[i],nfi,oXfield+cumf,fctx);
1383: } else {
1384: VecCopy(ctx->vecs[i],(Vec)(oXfield[cumf]));
1385: }
1386: cumf += nfi;
1387: }
1388: DMCompositeRestoreAccessArray(ctx->dm,X,n,NULL,ctx->vecs);
1389: return(0);
1390: }
1392: static PetscErrorCode DMSetUpGLVisViewer_Composite(PetscObject odm, PetscViewer viewer)
1393: {
1394: DM dm = (DM)odm, *dms;
1395: Vec *Ufds;
1396: GLVisViewerCtx *ctx;
1397: PetscInt i,n,tnf,*sdim;
1398: char **fecs;
1402: PetscNew(&ctx);
1403: PetscObjectReference((PetscObject)dm);
1404: ctx->dm = dm;
1405: DMCompositeGetNumberDM(dm,&n);
1406: PetscMalloc1(n,&dms);
1407: DMCompositeGetEntriesArray(dm,dms);
1408: PetscMalloc2(n,&ctx->subv,n,&ctx->vecs);
1409: for (i = 0, tnf = 0; i < n; i++) {
1410: PetscInt nf;
1412: PetscViewerCreate(PetscObjectComm(odm),&ctx->subv[i]);
1413: PetscViewerSetType(ctx->subv[i],PETSCVIEWERGLVIS);
1414: PetscViewerGLVisSetDM_Private(ctx->subv[i],(PetscObject)dms[i]);
1415: PetscViewerGLVisGetFields_Private(ctx->subv[i],&nf,NULL,NULL,NULL,NULL,NULL);
1416: tnf += nf;
1417: }
1418: PetscFree(dms);
1419: PetscMalloc3(tnf,&fecs,tnf,&sdim,tnf,&Ufds);
1420: for (i = 0, tnf = 0; i < n; i++) {
1421: PetscInt *sd,nf,f;
1422: const char **fec;
1423: Vec *Uf;
1425: PetscViewerGLVisGetFields_Private(ctx->subv[i],&nf,&fec,&sd,NULL,(PetscObject**)&Uf,NULL);
1426: for (f = 0; f < nf; f++) {
1427: PetscStrallocpy(fec[f],&fecs[tnf+f]);
1428: Ufds[tnf+f] = Uf[f];
1429: sdim[tnf+f] = sd[f];
1430: }
1431: tnf += nf;
1432: }
1433: PetscViewerGLVisSetFields(viewer,tnf,(const char**)fecs,sdim,DMCompositeSampleGLVisFields_Private,(PetscObject*)Ufds,ctx,DestroyGLVisViewerCtx_Private);
1434: for (i = 0; i < tnf; i++) {
1435: PetscFree(fecs[i]);
1436: }
1437: PetscFree3(fecs,sdim,Ufds);
1438: return(0);
1439: }
1441: PetscErrorCode DMRefine_Composite(DM dmi,MPI_Comm comm,DM *fine)
1442: {
1443: PetscErrorCode ierr;
1444: struct DMCompositeLink *next;
1445: DM_Composite *com = (DM_Composite*)dmi->data;
1446: DM dm;
1450: if (comm == MPI_COMM_NULL) {
1451: PetscObjectGetComm((PetscObject)dmi,&comm);
1452: }
1453: DMSetUp(dmi);
1454: next = com->next;
1455: DMCompositeCreate(comm,fine);
1457: /* loop over packed objects, handling one at at time */
1458: while (next) {
1459: DMRefine(next->dm,comm,&dm);
1460: DMCompositeAddDM(*fine,dm);
1461: PetscObjectDereference((PetscObject)dm);
1462: next = next->next;
1463: }
1464: return(0);
1465: }
1467: PetscErrorCode DMCoarsen_Composite(DM dmi,MPI_Comm comm,DM *fine)
1468: {
1469: PetscErrorCode ierr;
1470: struct DMCompositeLink *next;
1471: DM_Composite *com = (DM_Composite*)dmi->data;
1472: DM dm;
1476: DMSetUp(dmi);
1477: if (comm == MPI_COMM_NULL) {
1478: PetscObjectGetComm((PetscObject)dmi,&comm);
1479: }
1480: next = com->next;
1481: DMCompositeCreate(comm,fine);
1483: /* loop over packed objects, handling one at at time */
1484: while (next) {
1485: DMCoarsen(next->dm,comm,&dm);
1486: DMCompositeAddDM(*fine,dm);
1487: PetscObjectDereference((PetscObject)dm);
1488: next = next->next;
1489: }
1490: return(0);
1491: }
1493: PetscErrorCode DMCreateInterpolation_Composite(DM coarse,DM fine,Mat *A,Vec *v)
1494: {
1495: PetscErrorCode ierr;
1496: PetscInt m,n,M,N,nDM,i;
1497: struct DMCompositeLink *nextc;
1498: struct DMCompositeLink *nextf;
1499: Vec gcoarse,gfine,*vecs;
1500: DM_Composite *comcoarse = (DM_Composite*)coarse->data;
1501: DM_Composite *comfine = (DM_Composite*)fine->data;
1502: Mat *mats;
1507: DMSetUp(coarse);
1508: DMSetUp(fine);
1509: /* use global vectors only for determining matrix layout */
1510: DMGetGlobalVector(coarse,&gcoarse);
1511: DMGetGlobalVector(fine,&gfine);
1512: VecGetLocalSize(gcoarse,&n);
1513: VecGetLocalSize(gfine,&m);
1514: VecGetSize(gcoarse,&N);
1515: VecGetSize(gfine,&M);
1516: DMRestoreGlobalVector(coarse,&gcoarse);
1517: DMRestoreGlobalVector(fine,&gfine);
1519: nDM = comfine->nDM;
1520: if (nDM != comcoarse->nDM) SETERRQ2(PetscObjectComm((PetscObject)fine),PETSC_ERR_ARG_INCOMP,"Fine DMComposite has %D entries, but coarse has %D",nDM,comcoarse->nDM);
1521: PetscCalloc1(nDM*nDM,&mats);
1522: if (v) {
1523: PetscCalloc1(nDM,&vecs);
1524: }
1526: /* loop over packed objects, handling one at at time */
1527: for (nextc=comcoarse->next,nextf=comfine->next,i=0; nextc; nextc=nextc->next,nextf=nextf->next,i++) {
1528: if (!v) {
1529: DMCreateInterpolation(nextc->dm,nextf->dm,&mats[i*nDM+i],NULL);
1530: } else {
1531: DMCreateInterpolation(nextc->dm,nextf->dm,&mats[i*nDM+i],&vecs[i]);
1532: }
1533: }
1534: MatCreateNest(PetscObjectComm((PetscObject)fine),nDM,NULL,nDM,NULL,mats,A);
1535: if (v) {
1536: VecCreateNest(PetscObjectComm((PetscObject)fine),nDM,NULL,vecs,v);
1537: }
1538: for (i=0; i<nDM*nDM; i++) {MatDestroy(&mats[i]);}
1539: PetscFree(mats);
1540: if (v) {
1541: for (i=0; i<nDM; i++) {VecDestroy(&vecs[i]);}
1542: PetscFree(vecs);
1543: }
1544: return(0);
1545: }
1547: static PetscErrorCode DMGetLocalToGlobalMapping_Composite(DM dm)
1548: {
1549: DM_Composite *com = (DM_Composite*)dm->data;
1550: ISLocalToGlobalMapping *ltogs;
1551: PetscInt i;
1552: PetscErrorCode ierr;
1555: /* Set the ISLocalToGlobalMapping on the new matrix */
1556: DMCompositeGetISLocalToGlobalMappings(dm,<ogs);
1557: ISLocalToGlobalMappingConcatenate(PetscObjectComm((PetscObject)dm),com->nDM,ltogs,&dm->ltogmap);
1558: for (i=0; i<com->nDM; i++) {ISLocalToGlobalMappingDestroy(<ogs[i]);}
1559: PetscFree(ltogs);
1560: return(0);
1561: }
1564: PetscErrorCode DMCreateColoring_Composite(DM dm,ISColoringType ctype,ISColoring *coloring)
1565: {
1566: PetscErrorCode ierr;
1567: PetscInt n,i,cnt;
1568: ISColoringValue *colors;
1569: PetscBool dense = PETSC_FALSE;
1570: ISColoringValue maxcol = 0;
1571: DM_Composite *com = (DM_Composite*)dm->data;
1575: if (ctype == IS_COLORING_LOCAL) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Only global coloring supported");
1576: else if (ctype == IS_COLORING_GLOBAL) {
1577: n = com->n;
1578: } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unknown ISColoringType");
1579: PetscMalloc1(n,&colors); /* freed in ISColoringDestroy() */
1581: PetscOptionsGetBool(((PetscObject)dm)->options,((PetscObject)dm)->prefix,"-dmcomposite_dense_jacobian",&dense,NULL);
1582: if (dense) {
1583: for (i=0; i<n; i++) {
1584: colors[i] = (ISColoringValue)(com->rstart + i);
1585: }
1586: maxcol = com->N;
1587: } else {
1588: struct DMCompositeLink *next = com->next;
1589: PetscMPIInt rank;
1591: MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);
1592: cnt = 0;
1593: while (next) {
1594: ISColoring lcoloring;
1596: DMCreateColoring(next->dm,IS_COLORING_GLOBAL,&lcoloring);
1597: for (i=0; i<lcoloring->N; i++) {
1598: colors[cnt++] = maxcol + lcoloring->colors[i];
1599: }
1600: maxcol += lcoloring->n;
1601: ISColoringDestroy(&lcoloring);
1602: next = next->next;
1603: }
1604: }
1605: ISColoringCreate(PetscObjectComm((PetscObject)dm),maxcol,n,colors,PETSC_OWN_POINTER,coloring);
1606: return(0);
1607: }
1609: PetscErrorCode DMGlobalToLocalBegin_Composite(DM dm,Vec gvec,InsertMode mode,Vec lvec)
1610: {
1611: PetscErrorCode ierr;
1612: struct DMCompositeLink *next;
1613: PetscScalar *garray,*larray;
1614: DM_Composite *com = (DM_Composite*)dm->data;
1620: if (!com->setup) {
1621: DMSetUp(dm);
1622: }
1624: VecGetArray(gvec,&garray);
1625: VecGetArray(lvec,&larray);
1627: /* loop over packed objects, handling one at at time */
1628: next = com->next;
1629: while (next) {
1630: Vec local,global;
1631: PetscInt N;
1633: DMGetGlobalVector(next->dm,&global);
1634: VecGetLocalSize(global,&N);
1635: VecPlaceArray(global,garray);
1636: DMGetLocalVector(next->dm,&local);
1637: VecPlaceArray(local,larray);
1638: DMGlobalToLocalBegin(next->dm,global,mode,local);
1639: DMGlobalToLocalEnd(next->dm,global,mode,local);
1640: VecResetArray(global);
1641: VecResetArray(local);
1642: DMRestoreGlobalVector(next->dm,&global);
1643: DMRestoreLocalVector(next->dm,&local);
1645: larray += next->nlocal;
1646: garray += next->n;
1647: next = next->next;
1648: }
1650: VecRestoreArray(gvec,NULL);
1651: VecRestoreArray(lvec,NULL);
1652: return(0);
1653: }
1655: PetscErrorCode DMGlobalToLocalEnd_Composite(DM dm,Vec gvec,InsertMode mode,Vec lvec)
1656: {
1661: return(0);
1662: }
1664: PetscErrorCode DMLocalToGlobalBegin_Composite(DM dm,Vec lvec,InsertMode mode,Vec gvec)
1665: {
1666: PetscErrorCode ierr;
1667: struct DMCompositeLink *next;
1668: PetscScalar *larray,*garray;
1669: DM_Composite *com = (DM_Composite*)dm->data;
1676: if (!com->setup) {
1677: DMSetUp(dm);
1678: }
1680: VecGetArray(lvec,&larray);
1681: VecGetArray(gvec,&garray);
1683: /* loop over packed objects, handling one at at time */
1684: next = com->next;
1685: while (next) {
1686: Vec global,local;
1688: DMGetLocalVector(next->dm,&local);
1689: VecPlaceArray(local,larray);
1690: DMGetGlobalVector(next->dm,&global);
1691: VecPlaceArray(global,garray);
1692: DMLocalToGlobalBegin(next->dm,local,mode,global);
1693: DMLocalToGlobalEnd(next->dm,local,mode,global);
1694: VecResetArray(local);
1695: VecResetArray(global);
1696: DMRestoreGlobalVector(next->dm,&global);
1697: DMRestoreLocalVector(next->dm,&local);
1699: garray += next->n;
1700: larray += next->nlocal;
1701: next = next->next;
1702: }
1704: VecRestoreArray(gvec,NULL);
1705: VecRestoreArray(lvec,NULL);
1707: return(0);
1708: }
1710: PetscErrorCode DMLocalToGlobalEnd_Composite(DM dm,Vec lvec,InsertMode mode,Vec gvec)
1711: {
1716: return(0);
1717: }
1719: PetscErrorCode DMLocalToLocalBegin_Composite(DM dm,Vec vec1,InsertMode mode,Vec vec2)
1720: {
1721: PetscErrorCode ierr;
1722: struct DMCompositeLink *next;
1723: PetscScalar *array1,*array2;
1724: DM_Composite *com = (DM_Composite*)dm->data;
1731: if (!com->setup) {
1732: DMSetUp(dm);
1733: }
1735: VecGetArray(vec1,&array1);
1736: VecGetArray(vec2,&array2);
1738: /* loop over packed objects, handling one at at time */
1739: next = com->next;
1740: while (next) {
1741: Vec local1,local2;
1743: DMGetLocalVector(next->dm,&local1);
1744: VecPlaceArray(local1,array1);
1745: DMGetLocalVector(next->dm,&local2);
1746: VecPlaceArray(local2,array2);
1747: DMLocalToLocalBegin(next->dm,local1,mode,local2);
1748: DMLocalToLocalEnd(next->dm,local1,mode,local2);
1749: VecResetArray(local2);
1750: DMRestoreLocalVector(next->dm,&local2);
1751: VecResetArray(local1);
1752: DMRestoreLocalVector(next->dm,&local1);
1754: array1 += next->nlocal;
1755: array2 += next->nlocal;
1756: next = next->next;
1757: }
1759: VecRestoreArray(vec1,NULL);
1760: VecRestoreArray(vec2,NULL);
1762: return(0);
1763: }
1765: PetscErrorCode DMLocalToLocalEnd_Composite(DM dm,Vec lvec,InsertMode mode,Vec gvec)
1766: {
1771: return(0);
1772: }
1774: /*MC
1775: DMCOMPOSITE = "composite" - A DM object that is used to manage data for a collection of DMs
1777: Level: intermediate
1779: .seealso: DMType, DM, DMDACreate(), DMCreate(), DMSetType(), DMCompositeCreate()
1780: M*/
1783: PETSC_EXTERN PetscErrorCode DMCreate_Composite(DM p)
1784: {
1786: DM_Composite *com;
1789: PetscNewLog(p,&com);
1790: p->data = com;
1791: PetscObjectChangeTypeName((PetscObject)p,"DMComposite");
1792: com->n = 0;
1793: com->nghost = 0;
1794: com->next = NULL;
1795: com->nDM = 0;
1797: p->ops->createglobalvector = DMCreateGlobalVector_Composite;
1798: p->ops->createlocalvector = DMCreateLocalVector_Composite;
1799: p->ops->getlocaltoglobalmapping = DMGetLocalToGlobalMapping_Composite;
1800: p->ops->createfieldis = DMCreateFieldIS_Composite;
1801: p->ops->createfielddecomposition = DMCreateFieldDecomposition_Composite;
1802: p->ops->refine = DMRefine_Composite;
1803: p->ops->coarsen = DMCoarsen_Composite;
1804: p->ops->createinterpolation = DMCreateInterpolation_Composite;
1805: p->ops->creatematrix = DMCreateMatrix_Composite;
1806: p->ops->getcoloring = DMCreateColoring_Composite;
1807: p->ops->globaltolocalbegin = DMGlobalToLocalBegin_Composite;
1808: p->ops->globaltolocalend = DMGlobalToLocalEnd_Composite;
1809: p->ops->localtoglobalbegin = DMLocalToGlobalBegin_Composite;
1810: p->ops->localtoglobalend = DMLocalToGlobalEnd_Composite;
1811: p->ops->localtolocalbegin = DMLocalToLocalBegin_Composite;
1812: p->ops->localtolocalend = DMLocalToLocalEnd_Composite;
1813: p->ops->destroy = DMDestroy_Composite;
1814: p->ops->view = DMView_Composite;
1815: p->ops->setup = DMSetUp_Composite;
1817: PetscObjectComposeFunction((PetscObject)p,"DMSetUpGLVisViewer_C",DMSetUpGLVisViewer_Composite);
1818: return(0);
1819: }
1821: /*@
1822: DMCompositeCreate - Creates a vector packer, used to generate "composite"
1823: vectors made up of several subvectors.
1825: Collective on MPI_Comm
1827: Input Parameter:
1828: . comm - the processors that will share the global vector
1830: Output Parameters:
1831: . packer - the packer object
1833: Level: advanced
1835: .seealso DMDestroy(), DMCompositeAddDM(), DMCompositeScatter(), DMCOMPOSITE,DMCreate()
1836: DMCompositeGather(), DMCreateGlobalVector(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess()
1837: DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(), DMCompositeGetEntries()
1839: @*/
1840: PetscErrorCode DMCompositeCreate(MPI_Comm comm,DM *packer)
1841: {
1846: DMCreate(comm,packer);
1847: DMSetType(*packer,DMCOMPOSITE);
1848: return(0);
1849: }