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.
12: Logically Collective
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 dm
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: VecLockReadPush(*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 dm
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: VecLockReadPush(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 dm.
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: VecLockReadPush(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 dm
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: VecLockReadPop(*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 dm
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: VecLockReadPop(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 dm.
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: VecLockReadPop(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 dm
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 dm
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 dm
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 dm
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 dm
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: DMSetFromOptions(dm);
868: DMSetUp(dm);
869: VecCreate(PetscObjectComm((PetscObject)dm),gvec);
870: VecSetType(*gvec,dm->vectype);
871: VecSetSizes(*gvec,com->n,com->N);
872: VecSetDM(*gvec, dm);
873: VecSetOperation(*gvec,VECOP_VIEW,(void (*)(void))VecView_DMComposite);
874: return(0);
875: }
877: PetscErrorCode DMCreateLocalVector_Composite(DM dm,Vec *lvec)
878: {
880: DM_Composite *com = (DM_Composite*)dm->data;
884: if (!com->setup) {
885: DMSetFromOptions(dm);
886: DMSetUp(dm);
887: }
888: VecCreate(PETSC_COMM_SELF,lvec);
889: VecSetType(*lvec,dm->vectype);
890: VecSetSizes(*lvec,com->nghost,PETSC_DECIDE);
891: VecSetDM(*lvec, dm);
892: return(0);
893: }
895: /*@C
896: DMCompositeGetISLocalToGlobalMappings - gets an ISLocalToGlobalMapping for each DM in the DMComposite, maps to the composite global space
898: Collective on DM
900: Input Parameter:
901: . dm - the packer object
903: Output Parameters:
904: . ltogs - the individual mappings for each packed vector. Note that this includes
905: all the ghost points that individual ghosted DMDA's may have.
907: Level: advanced
909: Notes:
910: Each entry of ltogs should be destroyed with ISLocalToGlobalMappingDestroy(), the ltogs array should be freed with PetscFree().
912: Not available from Fortran
914: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(),
915: DMCompositeGather(), DMCompositeCreate(), DMCompositeGetAccess(), DMCompositeScatter(),
916: DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(),DMCompositeGetEntries()
918: @*/
919: PetscErrorCode DMCompositeGetISLocalToGlobalMappings(DM dm,ISLocalToGlobalMapping **ltogs)
920: {
921: PetscErrorCode ierr;
922: PetscInt i,*idx,n,cnt;
923: struct DMCompositeLink *next;
924: PetscMPIInt rank;
925: DM_Composite *com = (DM_Composite*)dm->data;
926: PetscBool flg;
930: PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
931: if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
932: DMSetUp(dm);
933: PetscMalloc1(com->nDM,ltogs);
934: next = com->next;
935: MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);
937: /* loop over packed objects, handling one at at time */
938: cnt = 0;
939: while (next) {
940: ISLocalToGlobalMapping ltog;
941: PetscMPIInt size;
942: const PetscInt *suboff,*indices;
943: Vec global;
945: /* Get sub-DM global indices for each local dof */
946: DMGetLocalToGlobalMapping(next->dm,<og);
947: ISLocalToGlobalMappingGetSize(ltog,&n);
948: ISLocalToGlobalMappingGetIndices(ltog,&indices);
949: PetscMalloc1(n,&idx);
951: /* Get the offsets for the sub-DM global vector */
952: DMGetGlobalVector(next->dm,&global);
953: VecGetOwnershipRanges(global,&suboff);
954: MPI_Comm_size(PetscObjectComm((PetscObject)global),&size);
956: /* Shift the sub-DM definition of the global space to the composite global space */
957: for (i=0; i<n; i++) {
958: PetscInt subi = indices[i],lo = 0,hi = size,t;
959: /* There's no consensus on what a negative index means,
960: except for skipping when setting the values in vectors and matrices */
961: if (subi < 0) { idx[i] = subi - next->grstarts[rank]; continue; }
962: /* Binary search to find which rank owns subi */
963: while (hi-lo > 1) {
964: t = lo + (hi-lo)/2;
965: if (suboff[t] > subi) hi = t;
966: else lo = t;
967: }
968: idx[i] = subi - suboff[lo] + next->grstarts[lo];
969: }
970: ISLocalToGlobalMappingRestoreIndices(ltog,&indices);
971: ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)dm),1,n,idx,PETSC_OWN_POINTER,&(*ltogs)[cnt]);
972: DMRestoreGlobalVector(next->dm,&global);
973: next = next->next;
974: cnt++;
975: }
976: return(0);
977: }
979: /*@C
980: DMCompositeGetLocalISs - Gets index sets for each component of a composite local vector
982: Not Collective
984: Input Arguments:
985: . dm - composite DM
987: Output Arguments:
988: . is - array of serial index sets for each each component of the DMComposite
990: Level: intermediate
992: Notes:
993: At present, a composite local vector does not normally exist. This function is used to provide index sets for
994: MatGetLocalSubMatrix(). In the future, the scatters for each entry in the DMComposite may be be merged into a single
995: scatter to a composite local vector. The user should not typically need to know which is being done.
997: To get the composite global indices at all local points (including ghosts), use DMCompositeGetISLocalToGlobalMappings().
999: To get index sets for pieces of the composite global vector, use DMCompositeGetGlobalISs().
1001: Each returned IS should be destroyed with ISDestroy(), the array should be freed with PetscFree().
1003: Not available from Fortran
1005: .seealso: DMCompositeGetGlobalISs(), DMCompositeGetISLocalToGlobalMappings(), MatGetLocalSubMatrix(), MatCreateLocalRef()
1006: @*/
1007: PetscErrorCode DMCompositeGetLocalISs(DM dm,IS **is)
1008: {
1009: PetscErrorCode ierr;
1010: DM_Composite *com = (DM_Composite*)dm->data;
1011: struct DMCompositeLink *link;
1012: PetscInt cnt,start;
1013: PetscBool flg;
1018: PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
1019: if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
1020: PetscMalloc1(com->nmine,is);
1021: for (cnt=0,start=0,link=com->next; link; start+=link->nlocal,cnt++,link=link->next) {
1022: PetscInt bs;
1023: ISCreateStride(PETSC_COMM_SELF,link->nlocal,start,1,&(*is)[cnt]);
1024: DMGetBlockSize(link->dm,&bs);
1025: ISSetBlockSize((*is)[cnt],bs);
1026: }
1027: return(0);
1028: }
1030: /*@C
1031: DMCompositeGetGlobalISs - Gets the index sets for each composed object
1033: Collective on dm
1035: Input Parameter:
1036: . dm - the packer object
1038: Output Parameters:
1039: . is - the array of index sets
1041: Level: advanced
1043: Notes:
1044: The is entries should be destroyed with ISDestroy(), the is array should be freed with PetscFree()
1046: These could be used to extract a subset of vector entries for a "multi-physics" preconditioner
1048: Use DMCompositeGetLocalISs() for index sets in the packed local numbering, and
1049: DMCompositeGetISLocalToGlobalMappings() for to map local sub-DM (including ghost) indices to packed global
1050: indices.
1052: Fortran Notes:
1054: The output argument 'is' must be an allocated array of sufficient length, which can be learned using DMCompositeGetNumberDM().
1056: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(),
1057: DMCompositeGather(), DMCompositeCreate(), DMCompositeGetAccess(), DMCompositeScatter(),
1058: DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(),DMCompositeGetEntries()
1060: @*/
1061: PetscErrorCode DMCompositeGetGlobalISs(DM dm,IS *is[])
1062: {
1063: PetscErrorCode ierr;
1064: PetscInt cnt = 0;
1065: struct DMCompositeLink *next;
1066: PetscMPIInt rank;
1067: DM_Composite *com = (DM_Composite*)dm->data;
1068: PetscBool flg;
1072: PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
1073: if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
1074: if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Must call DMSetUp() before");
1075: PetscMalloc1(com->nDM,is);
1076: next = com->next;
1077: MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);
1079: /* loop over packed objects, handling one at at time */
1080: while (next) {
1081: PetscDS prob;
1083: ISCreateStride(PetscObjectComm((PetscObject)dm),next->n,next->grstart,1,&(*is)[cnt]);
1084: DMGetDS(dm, &prob);
1085: if (prob) {
1086: MatNullSpace space;
1087: Mat pmat;
1088: PetscObject disc;
1089: PetscInt Nf;
1091: PetscDSGetNumFields(prob, &Nf);
1092: if (cnt < Nf) {
1093: PetscDSGetDiscretization(prob, cnt, &disc);
1094: PetscObjectQuery(disc, "nullspace", (PetscObject*) &space);
1095: if (space) {PetscObjectCompose((PetscObject) (*is)[cnt], "nullspace", (PetscObject) space);}
1096: PetscObjectQuery(disc, "nearnullspace", (PetscObject*) &space);
1097: if (space) {PetscObjectCompose((PetscObject) (*is)[cnt], "nearnullspace", (PetscObject) space);}
1098: PetscObjectQuery(disc, "pmat", (PetscObject*) &pmat);
1099: if (pmat) {PetscObjectCompose((PetscObject) (*is)[cnt], "pmat", (PetscObject) pmat);}
1100: }
1101: }
1102: cnt++;
1103: next = next->next;
1104: }
1105: return(0);
1106: }
1108: PetscErrorCode DMCreateFieldIS_Composite(DM dm, PetscInt *numFields,char ***fieldNames, IS **fields)
1109: {
1110: PetscInt nDM;
1111: DM *dms;
1112: PetscInt i;
1116: DMCompositeGetNumberDM(dm, &nDM);
1117: if (numFields) *numFields = nDM;
1118: DMCompositeGetGlobalISs(dm, fields);
1119: if (fieldNames) {
1120: PetscMalloc1(nDM, &dms);
1121: PetscMalloc1(nDM, fieldNames);
1122: DMCompositeGetEntriesArray(dm, dms);
1123: for (i=0; i<nDM; i++) {
1124: char buf[256];
1125: const char *splitname;
1127: /* Split naming precedence: object name, prefix, number */
1128: splitname = ((PetscObject) dm)->name;
1129: if (!splitname) {
1130: PetscObjectGetOptionsPrefix((PetscObject)dms[i],&splitname);
1131: if (splitname) {
1132: size_t len;
1133: PetscStrncpy(buf,splitname,sizeof(buf));
1134: buf[sizeof(buf) - 1] = 0;
1135: PetscStrlen(buf,&len);
1136: if (buf[len-1] == '_') buf[len-1] = 0; /* Remove trailing underscore if it was used */
1137: splitname = buf;
1138: }
1139: }
1140: if (!splitname) {
1141: PetscSNPrintf(buf,sizeof(buf),"%D",i);
1142: splitname = buf;
1143: }
1144: PetscStrallocpy(splitname,&(*fieldNames)[i]);
1145: }
1146: PetscFree(dms);
1147: }
1148: return(0);
1149: }
1151: /*
1152: This could take over from DMCreateFieldIS(), as it is more general,
1153: making DMCreateFieldIS() a special case -- calling with dmlist == NULL;
1154: At this point it's probably best to be less intrusive, however.
1155: */
1156: PetscErrorCode DMCreateFieldDecomposition_Composite(DM dm, PetscInt *len,char ***namelist, IS **islist, DM **dmlist)
1157: {
1158: PetscInt nDM;
1159: PetscInt i;
1163: DMCreateFieldIS_Composite(dm, len, namelist, islist);
1164: if (dmlist) {
1165: DMCompositeGetNumberDM(dm, &nDM);
1166: PetscMalloc1(nDM, dmlist);
1167: DMCompositeGetEntriesArray(dm, *dmlist);
1168: for (i=0; i<nDM; i++) {
1169: PetscObjectReference((PetscObject)((*dmlist)[i]));
1170: }
1171: }
1172: return(0);
1173: }
1177: /* -------------------------------------------------------------------------------------*/
1178: /*@C
1179: DMCompositeGetLocalVectors - Gets local vectors for each part of a DMComposite.
1180: Use DMCompositeRestoreLocalVectors() to return them.
1182: Not Collective
1184: Input Parameter:
1185: . dm - the packer object
1187: Output Parameter:
1188: . Vec ... - the individual sequential Vecs
1190: Level: advanced
1192: Not available from Fortran
1194: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(),
1195: DMCompositeGather(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
1196: DMCompositeRestoreLocalVectors(), DMCompositeScatter(), DMCompositeGetEntries()
1198: @*/
1199: PetscErrorCode DMCompositeGetLocalVectors(DM dm,...)
1200: {
1201: va_list Argp;
1202: PetscErrorCode ierr;
1203: struct DMCompositeLink *next;
1204: DM_Composite *com = (DM_Composite*)dm->data;
1205: PetscBool flg;
1209: PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
1210: if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
1211: next = com->next;
1212: /* loop over packed objects, handling one at at time */
1213: va_start(Argp,dm);
1214: while (next) {
1215: Vec *vec;
1216: vec = va_arg(Argp, Vec*);
1217: if (vec) {DMGetLocalVector(next->dm,vec);}
1218: next = next->next;
1219: }
1220: va_end(Argp);
1221: return(0);
1222: }
1224: /*@C
1225: DMCompositeRestoreLocalVectors - Restores local vectors for each part of a DMComposite.
1227: Not Collective
1229: Input Parameter:
1230: . dm - the packer object
1232: Output Parameter:
1233: . Vec ... - the individual sequential Vecs
1235: Level: advanced
1237: Not available from Fortran
1239: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(),
1240: DMCompositeGather(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
1241: DMCompositeGetLocalVectors(), DMCompositeScatter(), DMCompositeGetEntries()
1243: @*/
1244: PetscErrorCode DMCompositeRestoreLocalVectors(DM dm,...)
1245: {
1246: va_list Argp;
1247: PetscErrorCode ierr;
1248: struct DMCompositeLink *next;
1249: DM_Composite *com = (DM_Composite*)dm->data;
1250: PetscBool flg;
1254: PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
1255: if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
1256: next = com->next;
1257: /* loop over packed objects, handling one at at time */
1258: va_start(Argp,dm);
1259: while (next) {
1260: Vec *vec;
1261: vec = va_arg(Argp, Vec*);
1262: if (vec) {DMRestoreLocalVector(next->dm,vec);}
1263: next = next->next;
1264: }
1265: va_end(Argp);
1266: return(0);
1267: }
1269: /* -------------------------------------------------------------------------------------*/
1270: /*@C
1271: DMCompositeGetEntries - Gets the DM for each entry in a DMComposite.
1273: Not Collective
1275: Input Parameter:
1276: . dm - the packer object
1278: Output Parameter:
1279: . DM ... - the individual entries (DMs)
1281: Level: advanced
1283: Fortran Notes:
1284: Available as DMCompositeGetEntries() for one output DM, DMCompositeGetEntries2() for 2, etc
1286: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(), DMCompositeGetEntriesArray()
1287: DMCompositeGather(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
1288: DMCompositeRestoreLocalVectors(), DMCompositeGetLocalVectors(), DMCompositeScatter(),
1289: DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors()
1291: @*/
1292: PetscErrorCode DMCompositeGetEntries(DM dm,...)
1293: {
1294: va_list Argp;
1295: struct DMCompositeLink *next;
1296: DM_Composite *com = (DM_Composite*)dm->data;
1297: PetscBool flg;
1298: PetscErrorCode ierr;
1302: PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
1303: if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
1304: next = com->next;
1305: /* loop over packed objects, handling one at at time */
1306: va_start(Argp,dm);
1307: while (next) {
1308: DM *dmn;
1309: dmn = va_arg(Argp, DM*);
1310: if (dmn) *dmn = next->dm;
1311: next = next->next;
1312: }
1313: va_end(Argp);
1314: return(0);
1315: }
1317: /*@C
1318: DMCompositeGetEntriesArray - Gets the DM for each entry in a DMComposite.
1320: Not Collective
1322: Input Parameter:
1323: . dm - the packer object
1325: Output Parameter:
1326: . dms - array of sufficient length (see DMCompositeGetNumberDM()) to hold the individual DMs
1328: Level: advanced
1330: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(), DMCompositeGetEntries()
1331: DMCompositeGather(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
1332: DMCompositeRestoreLocalVectors(), DMCompositeGetLocalVectors(), DMCompositeScatter(),
1333: DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors()
1335: @*/
1336: PetscErrorCode DMCompositeGetEntriesArray(DM dm,DM dms[])
1337: {
1338: struct DMCompositeLink *next;
1339: DM_Composite *com = (DM_Composite*)dm->data;
1340: PetscInt i;
1341: PetscBool flg;
1342: PetscErrorCode ierr;
1346: PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
1347: if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
1348: /* loop over packed objects, handling one at at time */
1349: for (next=com->next,i=0; next; next=next->next,i++) dms[i] = next->dm;
1350: return(0);
1351: }
1353: typedef struct {
1354: DM dm;
1355: PetscViewer *subv;
1356: Vec *vecs;
1357: } GLVisViewerCtx;
1359: static PetscErrorCode DestroyGLVisViewerCtx_Private(void *vctx)
1360: {
1361: GLVisViewerCtx *ctx = (GLVisViewerCtx*)vctx;
1362: PetscInt i,n;
1366: DMCompositeGetNumberDM(ctx->dm,&n);
1367: for (i = 0; i < n; i++) {
1368: PetscViewerDestroy(&ctx->subv[i]);
1369: }
1370: PetscFree2(ctx->subv,ctx->vecs);
1371: DMDestroy(&ctx->dm);
1372: PetscFree(ctx);
1373: return(0);
1374: }
1376: static PetscErrorCode DMCompositeSampleGLVisFields_Private(PetscObject oX, PetscInt nf, PetscObject oXfield[], void *vctx)
1377: {
1378: Vec X = (Vec)oX;
1379: GLVisViewerCtx *ctx = (GLVisViewerCtx*)vctx;
1380: PetscInt i,n,cumf;
1384: DMCompositeGetNumberDM(ctx->dm,&n);
1385: DMCompositeGetAccessArray(ctx->dm,X,n,NULL,ctx->vecs);
1386: for (i = 0, cumf = 0; i < n; i++) {
1387: PetscErrorCode (*g2l)(PetscObject,PetscInt,PetscObject[],void*);
1388: void *fctx;
1389: PetscInt nfi;
1391: PetscViewerGLVisGetFields_Private(ctx->subv[i],&nfi,NULL,NULL,&g2l,NULL,&fctx);
1392: if (!nfi) continue;
1393: if (g2l) {
1394: (*g2l)((PetscObject)ctx->vecs[i],nfi,oXfield+cumf,fctx);
1395: } else {
1396: VecCopy(ctx->vecs[i],(Vec)(oXfield[cumf]));
1397: }
1398: cumf += nfi;
1399: }
1400: DMCompositeRestoreAccessArray(ctx->dm,X,n,NULL,ctx->vecs);
1401: return(0);
1402: }
1404: static PetscErrorCode DMSetUpGLVisViewer_Composite(PetscObject odm, PetscViewer viewer)
1405: {
1406: DM dm = (DM)odm, *dms;
1407: Vec *Ufds;
1408: GLVisViewerCtx *ctx;
1409: PetscInt i,n,tnf,*sdim;
1410: char **fecs;
1414: PetscNew(&ctx);
1415: PetscObjectReference((PetscObject)dm);
1416: ctx->dm = dm;
1417: DMCompositeGetNumberDM(dm,&n);
1418: PetscMalloc1(n,&dms);
1419: DMCompositeGetEntriesArray(dm,dms);
1420: PetscMalloc2(n,&ctx->subv,n,&ctx->vecs);
1421: for (i = 0, tnf = 0; i < n; i++) {
1422: PetscInt nf;
1424: PetscViewerCreate(PetscObjectComm(odm),&ctx->subv[i]);
1425: PetscViewerSetType(ctx->subv[i],PETSCVIEWERGLVIS);
1426: PetscViewerGLVisSetDM_Private(ctx->subv[i],(PetscObject)dms[i]);
1427: PetscViewerGLVisGetFields_Private(ctx->subv[i],&nf,NULL,NULL,NULL,NULL,NULL);
1428: tnf += nf;
1429: }
1430: PetscFree(dms);
1431: PetscMalloc3(tnf,&fecs,tnf,&sdim,tnf,&Ufds);
1432: for (i = 0, tnf = 0; i < n; i++) {
1433: PetscInt *sd,nf,f;
1434: const char **fec;
1435: Vec *Uf;
1437: PetscViewerGLVisGetFields_Private(ctx->subv[i],&nf,&fec,&sd,NULL,(PetscObject**)&Uf,NULL);
1438: for (f = 0; f < nf; f++) {
1439: PetscStrallocpy(fec[f],&fecs[tnf+f]);
1440: Ufds[tnf+f] = Uf[f];
1441: sdim[tnf+f] = sd[f];
1442: }
1443: tnf += nf;
1444: }
1445: PetscViewerGLVisSetFields(viewer,tnf,(const char**)fecs,sdim,DMCompositeSampleGLVisFields_Private,(PetscObject*)Ufds,ctx,DestroyGLVisViewerCtx_Private);
1446: for (i = 0; i < tnf; i++) {
1447: PetscFree(fecs[i]);
1448: }
1449: PetscFree3(fecs,sdim,Ufds);
1450: return(0);
1451: }
1453: PetscErrorCode DMRefine_Composite(DM dmi,MPI_Comm comm,DM *fine)
1454: {
1455: PetscErrorCode ierr;
1456: struct DMCompositeLink *next;
1457: DM_Composite *com = (DM_Composite*)dmi->data;
1458: DM dm;
1462: if (comm == MPI_COMM_NULL) {
1463: PetscObjectGetComm((PetscObject)dmi,&comm);
1464: }
1465: DMSetUp(dmi);
1466: next = com->next;
1467: DMCompositeCreate(comm,fine);
1469: /* loop over packed objects, handling one at at time */
1470: while (next) {
1471: DMRefine(next->dm,comm,&dm);
1472: DMCompositeAddDM(*fine,dm);
1473: PetscObjectDereference((PetscObject)dm);
1474: next = next->next;
1475: }
1476: return(0);
1477: }
1479: PetscErrorCode DMCoarsen_Composite(DM dmi,MPI_Comm comm,DM *fine)
1480: {
1481: PetscErrorCode ierr;
1482: struct DMCompositeLink *next;
1483: DM_Composite *com = (DM_Composite*)dmi->data;
1484: DM dm;
1488: DMSetUp(dmi);
1489: if (comm == MPI_COMM_NULL) {
1490: PetscObjectGetComm((PetscObject)dmi,&comm);
1491: }
1492: next = com->next;
1493: DMCompositeCreate(comm,fine);
1495: /* loop over packed objects, handling one at at time */
1496: while (next) {
1497: DMCoarsen(next->dm,comm,&dm);
1498: DMCompositeAddDM(*fine,dm);
1499: PetscObjectDereference((PetscObject)dm);
1500: next = next->next;
1501: }
1502: return(0);
1503: }
1505: PetscErrorCode DMCreateInterpolation_Composite(DM coarse,DM fine,Mat *A,Vec *v)
1506: {
1507: PetscErrorCode ierr;
1508: PetscInt m,n,M,N,nDM,i;
1509: struct DMCompositeLink *nextc;
1510: struct DMCompositeLink *nextf;
1511: Vec gcoarse,gfine,*vecs;
1512: DM_Composite *comcoarse = (DM_Composite*)coarse->data;
1513: DM_Composite *comfine = (DM_Composite*)fine->data;
1514: Mat *mats;
1519: DMSetUp(coarse);
1520: DMSetUp(fine);
1521: /* use global vectors only for determining matrix layout */
1522: DMGetGlobalVector(coarse,&gcoarse);
1523: DMGetGlobalVector(fine,&gfine);
1524: VecGetLocalSize(gcoarse,&n);
1525: VecGetLocalSize(gfine,&m);
1526: VecGetSize(gcoarse,&N);
1527: VecGetSize(gfine,&M);
1528: DMRestoreGlobalVector(coarse,&gcoarse);
1529: DMRestoreGlobalVector(fine,&gfine);
1531: nDM = comfine->nDM;
1532: if (nDM != comcoarse->nDM) SETERRQ2(PetscObjectComm((PetscObject)fine),PETSC_ERR_ARG_INCOMP,"Fine DMComposite has %D entries, but coarse has %D",nDM,comcoarse->nDM);
1533: PetscCalloc1(nDM*nDM,&mats);
1534: if (v) {
1535: PetscCalloc1(nDM,&vecs);
1536: }
1538: /* loop over packed objects, handling one at at time */
1539: for (nextc=comcoarse->next,nextf=comfine->next,i=0; nextc; nextc=nextc->next,nextf=nextf->next,i++) {
1540: if (!v) {
1541: DMCreateInterpolation(nextc->dm,nextf->dm,&mats[i*nDM+i],NULL);
1542: } else {
1543: DMCreateInterpolation(nextc->dm,nextf->dm,&mats[i*nDM+i],&vecs[i]);
1544: }
1545: }
1546: MatCreateNest(PetscObjectComm((PetscObject)fine),nDM,NULL,nDM,NULL,mats,A);
1547: if (v) {
1548: VecCreateNest(PetscObjectComm((PetscObject)fine),nDM,NULL,vecs,v);
1549: }
1550: for (i=0; i<nDM*nDM; i++) {MatDestroy(&mats[i]);}
1551: PetscFree(mats);
1552: if (v) {
1553: for (i=0; i<nDM; i++) {VecDestroy(&vecs[i]);}
1554: PetscFree(vecs);
1555: }
1556: return(0);
1557: }
1559: static PetscErrorCode DMGetLocalToGlobalMapping_Composite(DM dm)
1560: {
1561: DM_Composite *com = (DM_Composite*)dm->data;
1562: ISLocalToGlobalMapping *ltogs;
1563: PetscInt i;
1564: PetscErrorCode ierr;
1567: /* Set the ISLocalToGlobalMapping on the new matrix */
1568: DMCompositeGetISLocalToGlobalMappings(dm,<ogs);
1569: ISLocalToGlobalMappingConcatenate(PetscObjectComm((PetscObject)dm),com->nDM,ltogs,&dm->ltogmap);
1570: for (i=0; i<com->nDM; i++) {ISLocalToGlobalMappingDestroy(<ogs[i]);}
1571: PetscFree(ltogs);
1572: return(0);
1573: }
1576: PetscErrorCode DMCreateColoring_Composite(DM dm,ISColoringType ctype,ISColoring *coloring)
1577: {
1578: PetscErrorCode ierr;
1579: PetscInt n,i,cnt;
1580: ISColoringValue *colors;
1581: PetscBool dense = PETSC_FALSE;
1582: ISColoringValue maxcol = 0;
1583: DM_Composite *com = (DM_Composite*)dm->data;
1587: if (ctype == IS_COLORING_LOCAL) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Only global coloring supported");
1588: else if (ctype == IS_COLORING_GLOBAL) {
1589: n = com->n;
1590: } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unknown ISColoringType");
1591: PetscMalloc1(n,&colors); /* freed in ISColoringDestroy() */
1593: PetscOptionsGetBool(((PetscObject)dm)->options,((PetscObject)dm)->prefix,"-dmcomposite_dense_jacobian",&dense,NULL);
1594: if (dense) {
1595: for (i=0; i<n; i++) {
1596: colors[i] = (ISColoringValue)(com->rstart + i);
1597: }
1598: maxcol = com->N;
1599: } else {
1600: struct DMCompositeLink *next = com->next;
1601: PetscMPIInt rank;
1603: MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);
1604: cnt = 0;
1605: while (next) {
1606: ISColoring lcoloring;
1608: DMCreateColoring(next->dm,IS_COLORING_GLOBAL,&lcoloring);
1609: for (i=0; i<lcoloring->N; i++) {
1610: colors[cnt++] = maxcol + lcoloring->colors[i];
1611: }
1612: maxcol += lcoloring->n;
1613: ISColoringDestroy(&lcoloring);
1614: next = next->next;
1615: }
1616: }
1617: ISColoringCreate(PetscObjectComm((PetscObject)dm),maxcol,n,colors,PETSC_OWN_POINTER,coloring);
1618: return(0);
1619: }
1621: PetscErrorCode DMGlobalToLocalBegin_Composite(DM dm,Vec gvec,InsertMode mode,Vec lvec)
1622: {
1623: PetscErrorCode ierr;
1624: struct DMCompositeLink *next;
1625: PetscScalar *garray,*larray;
1626: DM_Composite *com = (DM_Composite*)dm->data;
1632: if (!com->setup) {
1633: DMSetUp(dm);
1634: }
1636: VecGetArray(gvec,&garray);
1637: VecGetArray(lvec,&larray);
1639: /* loop over packed objects, handling one at at time */
1640: next = com->next;
1641: while (next) {
1642: Vec local,global;
1643: PetscInt N;
1645: DMGetGlobalVector(next->dm,&global);
1646: VecGetLocalSize(global,&N);
1647: VecPlaceArray(global,garray);
1648: DMGetLocalVector(next->dm,&local);
1649: VecPlaceArray(local,larray);
1650: DMGlobalToLocalBegin(next->dm,global,mode,local);
1651: DMGlobalToLocalEnd(next->dm,global,mode,local);
1652: VecResetArray(global);
1653: VecResetArray(local);
1654: DMRestoreGlobalVector(next->dm,&global);
1655: DMRestoreLocalVector(next->dm,&local);
1657: larray += next->nlocal;
1658: garray += next->n;
1659: next = next->next;
1660: }
1662: VecRestoreArray(gvec,NULL);
1663: VecRestoreArray(lvec,NULL);
1664: return(0);
1665: }
1667: PetscErrorCode DMGlobalToLocalEnd_Composite(DM dm,Vec gvec,InsertMode mode,Vec lvec)
1668: {
1673: return(0);
1674: }
1676: PetscErrorCode DMLocalToGlobalBegin_Composite(DM dm,Vec lvec,InsertMode mode,Vec gvec)
1677: {
1678: PetscErrorCode ierr;
1679: struct DMCompositeLink *next;
1680: PetscScalar *larray,*garray;
1681: DM_Composite *com = (DM_Composite*)dm->data;
1688: if (!com->setup) {
1689: DMSetUp(dm);
1690: }
1692: VecGetArray(lvec,&larray);
1693: VecGetArray(gvec,&garray);
1695: /* loop over packed objects, handling one at at time */
1696: next = com->next;
1697: while (next) {
1698: Vec global,local;
1700: DMGetLocalVector(next->dm,&local);
1701: VecPlaceArray(local,larray);
1702: DMGetGlobalVector(next->dm,&global);
1703: VecPlaceArray(global,garray);
1704: DMLocalToGlobalBegin(next->dm,local,mode,global);
1705: DMLocalToGlobalEnd(next->dm,local,mode,global);
1706: VecResetArray(local);
1707: VecResetArray(global);
1708: DMRestoreGlobalVector(next->dm,&global);
1709: DMRestoreLocalVector(next->dm,&local);
1711: garray += next->n;
1712: larray += next->nlocal;
1713: next = next->next;
1714: }
1716: VecRestoreArray(gvec,NULL);
1717: VecRestoreArray(lvec,NULL);
1719: return(0);
1720: }
1722: PetscErrorCode DMLocalToGlobalEnd_Composite(DM dm,Vec lvec,InsertMode mode,Vec gvec)
1723: {
1728: return(0);
1729: }
1731: PetscErrorCode DMLocalToLocalBegin_Composite(DM dm,Vec vec1,InsertMode mode,Vec vec2)
1732: {
1733: PetscErrorCode ierr;
1734: struct DMCompositeLink *next;
1735: PetscScalar *array1,*array2;
1736: DM_Composite *com = (DM_Composite*)dm->data;
1743: if (!com->setup) {
1744: DMSetUp(dm);
1745: }
1747: VecGetArray(vec1,&array1);
1748: VecGetArray(vec2,&array2);
1750: /* loop over packed objects, handling one at at time */
1751: next = com->next;
1752: while (next) {
1753: Vec local1,local2;
1755: DMGetLocalVector(next->dm,&local1);
1756: VecPlaceArray(local1,array1);
1757: DMGetLocalVector(next->dm,&local2);
1758: VecPlaceArray(local2,array2);
1759: DMLocalToLocalBegin(next->dm,local1,mode,local2);
1760: DMLocalToLocalEnd(next->dm,local1,mode,local2);
1761: VecResetArray(local2);
1762: DMRestoreLocalVector(next->dm,&local2);
1763: VecResetArray(local1);
1764: DMRestoreLocalVector(next->dm,&local1);
1766: array1 += next->nlocal;
1767: array2 += next->nlocal;
1768: next = next->next;
1769: }
1771: VecRestoreArray(vec1,NULL);
1772: VecRestoreArray(vec2,NULL);
1774: return(0);
1775: }
1777: PetscErrorCode DMLocalToLocalEnd_Composite(DM dm,Vec lvec,InsertMode mode,Vec gvec)
1778: {
1783: return(0);
1784: }
1786: /*MC
1787: DMCOMPOSITE = "composite" - A DM object that is used to manage data for a collection of DMs
1789: Level: intermediate
1791: .seealso: DMType, DM, DMDACreate(), DMCreate(), DMSetType(), DMCompositeCreate()
1792: M*/
1795: PETSC_EXTERN PetscErrorCode DMCreate_Composite(DM p)
1796: {
1798: DM_Composite *com;
1801: PetscNewLog(p,&com);
1802: p->data = com;
1803: com->n = 0;
1804: com->nghost = 0;
1805: com->next = NULL;
1806: com->nDM = 0;
1808: p->ops->createglobalvector = DMCreateGlobalVector_Composite;
1809: p->ops->createlocalvector = DMCreateLocalVector_Composite;
1810: p->ops->getlocaltoglobalmapping = DMGetLocalToGlobalMapping_Composite;
1811: p->ops->createfieldis = DMCreateFieldIS_Composite;
1812: p->ops->createfielddecomposition = DMCreateFieldDecomposition_Composite;
1813: p->ops->refine = DMRefine_Composite;
1814: p->ops->coarsen = DMCoarsen_Composite;
1815: p->ops->createinterpolation = DMCreateInterpolation_Composite;
1816: p->ops->creatematrix = DMCreateMatrix_Composite;
1817: p->ops->getcoloring = DMCreateColoring_Composite;
1818: p->ops->globaltolocalbegin = DMGlobalToLocalBegin_Composite;
1819: p->ops->globaltolocalend = DMGlobalToLocalEnd_Composite;
1820: p->ops->localtoglobalbegin = DMLocalToGlobalBegin_Composite;
1821: p->ops->localtoglobalend = DMLocalToGlobalEnd_Composite;
1822: p->ops->localtolocalbegin = DMLocalToLocalBegin_Composite;
1823: p->ops->localtolocalend = DMLocalToLocalEnd_Composite;
1824: p->ops->destroy = DMDestroy_Composite;
1825: p->ops->view = DMView_Composite;
1826: p->ops->setup = DMSetUp_Composite;
1828: PetscObjectComposeFunction((PetscObject)p,"DMSetUpGLVisViewer_C",DMSetUpGLVisViewer_Composite);
1829: return(0);
1830: }
1832: /*@
1833: DMCompositeCreate - Creates a vector packer, used to generate "composite"
1834: vectors made up of several subvectors.
1836: Collective
1838: Input Parameter:
1839: . comm - the processors that will share the global vector
1841: Output Parameters:
1842: . packer - the packer object
1844: Level: advanced
1846: .seealso DMDestroy(), DMCompositeAddDM(), DMCompositeScatter(), DMCOMPOSITE,DMCreate()
1847: DMCompositeGather(), DMCreateGlobalVector(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess()
1848: DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(), DMCompositeGetEntries()
1850: @*/
1851: PetscErrorCode DMCompositeCreate(MPI_Comm comm,DM *packer)
1852: {
1857: DMCreate(comm,packer);
1858: DMSetType(*packer,DMCOMPOSITE);
1859: return(0);
1860: }