Actual source code: pack.c
petsc-3.11.4 2019-09-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: 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 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: 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 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: 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 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: 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 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: 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 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: 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 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: PetscDS prob;
1074: ISCreateStride(PetscObjectComm((PetscObject)dm),next->n,next->grstart,1,&(*is)[cnt]);
1075: DMGetDS(dm, &prob);
1076: if (prob) {
1077: MatNullSpace space;
1078: Mat pmat;
1079: PetscObject disc;
1080: PetscInt Nf;
1082: PetscDSGetNumFields(prob, &Nf);
1083: if (cnt < Nf) {
1084: PetscDSGetDiscretization(prob, cnt, &disc);
1085: PetscObjectQuery(disc, "nullspace", (PetscObject*) &space);
1086: if (space) {PetscObjectCompose((PetscObject) (*is)[cnt], "nullspace", (PetscObject) space);}
1087: PetscObjectQuery(disc, "nearnullspace", (PetscObject*) &space);
1088: if (space) {PetscObjectCompose((PetscObject) (*is)[cnt], "nearnullspace", (PetscObject) space);}
1089: PetscObjectQuery(disc, "pmat", (PetscObject*) &pmat);
1090: if (pmat) {PetscObjectCompose((PetscObject) (*is)[cnt], "pmat", (PetscObject) pmat);}
1091: }
1092: }
1093: cnt++;
1094: next = next->next;
1095: }
1096: return(0);
1097: }
1099: PetscErrorCode DMCreateFieldIS_Composite(DM dm, PetscInt *numFields,char ***fieldNames, IS **fields)
1100: {
1101: PetscInt nDM;
1102: DM *dms;
1103: PetscInt i;
1107: DMCompositeGetNumberDM(dm, &nDM);
1108: if (numFields) *numFields = nDM;
1109: DMCompositeGetGlobalISs(dm, fields);
1110: if (fieldNames) {
1111: PetscMalloc1(nDM, &dms);
1112: PetscMalloc1(nDM, fieldNames);
1113: DMCompositeGetEntriesArray(dm, dms);
1114: for (i=0; i<nDM; i++) {
1115: char buf[256];
1116: const char *splitname;
1118: /* Split naming precedence: object name, prefix, number */
1119: splitname = ((PetscObject) dm)->name;
1120: if (!splitname) {
1121: PetscObjectGetOptionsPrefix((PetscObject)dms[i],&splitname);
1122: if (splitname) {
1123: size_t len;
1124: PetscStrncpy(buf,splitname,sizeof(buf));
1125: buf[sizeof(buf) - 1] = 0;
1126: PetscStrlen(buf,&len);
1127: if (buf[len-1] == '_') buf[len-1] = 0; /* Remove trailing underscore if it was used */
1128: splitname = buf;
1129: }
1130: }
1131: if (!splitname) {
1132: PetscSNPrintf(buf,sizeof(buf),"%D",i);
1133: splitname = buf;
1134: }
1135: PetscStrallocpy(splitname,&(*fieldNames)[i]);
1136: }
1137: PetscFree(dms);
1138: }
1139: return(0);
1140: }
1142: /*
1143: This could take over from DMCreateFieldIS(), as it is more general,
1144: making DMCreateFieldIS() a special case -- calling with dmlist == NULL;
1145: At this point it's probably best to be less intrusive, however.
1146: */
1147: PetscErrorCode DMCreateFieldDecomposition_Composite(DM dm, PetscInt *len,char ***namelist, IS **islist, DM **dmlist)
1148: {
1149: PetscInt nDM;
1150: PetscInt i;
1154: DMCreateFieldIS_Composite(dm, len, namelist, islist);
1155: if (dmlist) {
1156: DMCompositeGetNumberDM(dm, &nDM);
1157: PetscMalloc1(nDM, dmlist);
1158: DMCompositeGetEntriesArray(dm, *dmlist);
1159: for (i=0; i<nDM; i++) {
1160: PetscObjectReference((PetscObject)((*dmlist)[i]));
1161: }
1162: }
1163: return(0);
1164: }
1168: /* -------------------------------------------------------------------------------------*/
1169: /*@C
1170: DMCompositeGetLocalVectors - Gets local vectors for each part of a DMComposite.
1171: Use DMCompositeRestoreLocalVectors() to return them.
1173: Not Collective
1175: Input Parameter:
1176: . dm - the packer object
1178: Output Parameter:
1179: . Vec ... - the individual sequential Vecs
1181: Level: advanced
1183: Not available from Fortran
1185: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(),
1186: DMCompositeGather(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
1187: DMCompositeRestoreLocalVectors(), DMCompositeScatter(), DMCompositeGetEntries()
1189: @*/
1190: PetscErrorCode DMCompositeGetLocalVectors(DM dm,...)
1191: {
1192: va_list Argp;
1193: PetscErrorCode ierr;
1194: struct DMCompositeLink *next;
1195: DM_Composite *com = (DM_Composite*)dm->data;
1196: PetscBool flg;
1200: PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
1201: if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
1202: next = com->next;
1203: /* loop over packed objects, handling one at at time */
1204: va_start(Argp,dm);
1205: while (next) {
1206: Vec *vec;
1207: vec = va_arg(Argp, Vec*);
1208: if (vec) {DMGetLocalVector(next->dm,vec);}
1209: next = next->next;
1210: }
1211: va_end(Argp);
1212: return(0);
1213: }
1215: /*@C
1216: DMCompositeRestoreLocalVectors - Restores local vectors for each part of a DMComposite.
1218: Not Collective
1220: Input Parameter:
1221: . dm - the packer object
1223: Output Parameter:
1224: . Vec ... - the individual sequential Vecs
1226: Level: advanced
1228: Not available from Fortran
1230: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(),
1231: DMCompositeGather(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
1232: DMCompositeGetLocalVectors(), DMCompositeScatter(), DMCompositeGetEntries()
1234: @*/
1235: PetscErrorCode DMCompositeRestoreLocalVectors(DM dm,...)
1236: {
1237: va_list Argp;
1238: PetscErrorCode ierr;
1239: struct DMCompositeLink *next;
1240: DM_Composite *com = (DM_Composite*)dm->data;
1241: PetscBool flg;
1245: PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
1246: if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
1247: next = com->next;
1248: /* loop over packed objects, handling one at at time */
1249: va_start(Argp,dm);
1250: while (next) {
1251: Vec *vec;
1252: vec = va_arg(Argp, Vec*);
1253: if (vec) {DMRestoreLocalVector(next->dm,vec);}
1254: next = next->next;
1255: }
1256: va_end(Argp);
1257: return(0);
1258: }
1260: /* -------------------------------------------------------------------------------------*/
1261: /*@C
1262: DMCompositeGetEntries - Gets the DM for each entry in a DMComposite.
1264: Not Collective
1266: Input Parameter:
1267: . dm - the packer object
1269: Output Parameter:
1270: . DM ... - the individual entries (DMs)
1272: Level: advanced
1274: Fortran Notes:
1275: Available as DMCompositeGetEntries() for one output DM, DMCompositeGetEntries2() for 2, etc
1277: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(), DMCompositeGetEntriesArray()
1278: DMCompositeGather(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
1279: DMCompositeRestoreLocalVectors(), DMCompositeGetLocalVectors(), DMCompositeScatter(),
1280: DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors()
1282: @*/
1283: PetscErrorCode DMCompositeGetEntries(DM dm,...)
1284: {
1285: va_list Argp;
1286: struct DMCompositeLink *next;
1287: DM_Composite *com = (DM_Composite*)dm->data;
1288: PetscBool flg;
1289: PetscErrorCode ierr;
1293: PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
1294: if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
1295: next = com->next;
1296: /* loop over packed objects, handling one at at time */
1297: va_start(Argp,dm);
1298: while (next) {
1299: DM *dmn;
1300: dmn = va_arg(Argp, DM*);
1301: if (dmn) *dmn = next->dm;
1302: next = next->next;
1303: }
1304: va_end(Argp);
1305: return(0);
1306: }
1308: /*@C
1309: DMCompositeGetEntriesArray - Gets the DM for each entry in a DMComposite.
1311: Not Collective
1313: Input Parameter:
1314: . dm - the packer object
1316: Output Parameter:
1317: . dms - array of sufficient length (see DMCompositeGetNumberDM()) to hold the individual DMs
1319: Level: advanced
1321: .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(), DMCompositeGetEntries()
1322: DMCompositeGather(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
1323: DMCompositeRestoreLocalVectors(), DMCompositeGetLocalVectors(), DMCompositeScatter(),
1324: DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors()
1326: @*/
1327: PetscErrorCode DMCompositeGetEntriesArray(DM dm,DM dms[])
1328: {
1329: struct DMCompositeLink *next;
1330: DM_Composite *com = (DM_Composite*)dm->data;
1331: PetscInt i;
1332: PetscBool flg;
1333: PetscErrorCode ierr;
1337: PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);
1338: if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
1339: /* loop over packed objects, handling one at at time */
1340: for (next=com->next,i=0; next; next=next->next,i++) dms[i] = next->dm;
1341: return(0);
1342: }
1344: typedef struct {
1345: DM dm;
1346: PetscViewer *subv;
1347: Vec *vecs;
1348: } GLVisViewerCtx;
1350: static PetscErrorCode DestroyGLVisViewerCtx_Private(void *vctx)
1351: {
1352: GLVisViewerCtx *ctx = (GLVisViewerCtx*)vctx;
1353: PetscInt i,n;
1357: DMCompositeGetNumberDM(ctx->dm,&n);
1358: for (i = 0; i < n; i++) {
1359: PetscViewerDestroy(&ctx->subv[i]);
1360: }
1361: PetscFree2(ctx->subv,ctx->vecs);
1362: DMDestroy(&ctx->dm);
1363: PetscFree(ctx);
1364: return(0);
1365: }
1367: static PetscErrorCode DMCompositeSampleGLVisFields_Private(PetscObject oX, PetscInt nf, PetscObject oXfield[], void *vctx)
1368: {
1369: Vec X = (Vec)oX;
1370: GLVisViewerCtx *ctx = (GLVisViewerCtx*)vctx;
1371: PetscInt i,n,cumf;
1375: DMCompositeGetNumberDM(ctx->dm,&n);
1376: DMCompositeGetAccessArray(ctx->dm,X,n,NULL,ctx->vecs);
1377: for (i = 0, cumf = 0; i < n; i++) {
1378: PetscErrorCode (*g2l)(PetscObject,PetscInt,PetscObject[],void*);
1379: void *fctx;
1380: PetscInt nfi;
1382: PetscViewerGLVisGetFields_Private(ctx->subv[i],&nfi,NULL,NULL,&g2l,NULL,&fctx);
1383: if (!nfi) continue;
1384: if (g2l) {
1385: (*g2l)((PetscObject)ctx->vecs[i],nfi,oXfield+cumf,fctx);
1386: } else {
1387: VecCopy(ctx->vecs[i],(Vec)(oXfield[cumf]));
1388: }
1389: cumf += nfi;
1390: }
1391: DMCompositeRestoreAccessArray(ctx->dm,X,n,NULL,ctx->vecs);
1392: return(0);
1393: }
1395: static PetscErrorCode DMSetUpGLVisViewer_Composite(PetscObject odm, PetscViewer viewer)
1396: {
1397: DM dm = (DM)odm, *dms;
1398: Vec *Ufds;
1399: GLVisViewerCtx *ctx;
1400: PetscInt i,n,tnf,*sdim;
1401: char **fecs;
1405: PetscNew(&ctx);
1406: PetscObjectReference((PetscObject)dm);
1407: ctx->dm = dm;
1408: DMCompositeGetNumberDM(dm,&n);
1409: PetscMalloc1(n,&dms);
1410: DMCompositeGetEntriesArray(dm,dms);
1411: PetscMalloc2(n,&ctx->subv,n,&ctx->vecs);
1412: for (i = 0, tnf = 0; i < n; i++) {
1413: PetscInt nf;
1415: PetscViewerCreate(PetscObjectComm(odm),&ctx->subv[i]);
1416: PetscViewerSetType(ctx->subv[i],PETSCVIEWERGLVIS);
1417: PetscViewerGLVisSetDM_Private(ctx->subv[i],(PetscObject)dms[i]);
1418: PetscViewerGLVisGetFields_Private(ctx->subv[i],&nf,NULL,NULL,NULL,NULL,NULL);
1419: tnf += nf;
1420: }
1421: PetscFree(dms);
1422: PetscMalloc3(tnf,&fecs,tnf,&sdim,tnf,&Ufds);
1423: for (i = 0, tnf = 0; i < n; i++) {
1424: PetscInt *sd,nf,f;
1425: const char **fec;
1426: Vec *Uf;
1428: PetscViewerGLVisGetFields_Private(ctx->subv[i],&nf,&fec,&sd,NULL,(PetscObject**)&Uf,NULL);
1429: for (f = 0; f < nf; f++) {
1430: PetscStrallocpy(fec[f],&fecs[tnf+f]);
1431: Ufds[tnf+f] = Uf[f];
1432: sdim[tnf+f] = sd[f];
1433: }
1434: tnf += nf;
1435: }
1436: PetscViewerGLVisSetFields(viewer,tnf,(const char**)fecs,sdim,DMCompositeSampleGLVisFields_Private,(PetscObject*)Ufds,ctx,DestroyGLVisViewerCtx_Private);
1437: for (i = 0; i < tnf; i++) {
1438: PetscFree(fecs[i]);
1439: }
1440: PetscFree3(fecs,sdim,Ufds);
1441: return(0);
1442: }
1444: PetscErrorCode DMRefine_Composite(DM dmi,MPI_Comm comm,DM *fine)
1445: {
1446: PetscErrorCode ierr;
1447: struct DMCompositeLink *next;
1448: DM_Composite *com = (DM_Composite*)dmi->data;
1449: DM dm;
1453: if (comm == MPI_COMM_NULL) {
1454: PetscObjectGetComm((PetscObject)dmi,&comm);
1455: }
1456: DMSetUp(dmi);
1457: next = com->next;
1458: DMCompositeCreate(comm,fine);
1460: /* loop over packed objects, handling one at at time */
1461: while (next) {
1462: DMRefine(next->dm,comm,&dm);
1463: DMCompositeAddDM(*fine,dm);
1464: PetscObjectDereference((PetscObject)dm);
1465: next = next->next;
1466: }
1467: return(0);
1468: }
1470: PetscErrorCode DMCoarsen_Composite(DM dmi,MPI_Comm comm,DM *fine)
1471: {
1472: PetscErrorCode ierr;
1473: struct DMCompositeLink *next;
1474: DM_Composite *com = (DM_Composite*)dmi->data;
1475: DM dm;
1479: DMSetUp(dmi);
1480: if (comm == MPI_COMM_NULL) {
1481: PetscObjectGetComm((PetscObject)dmi,&comm);
1482: }
1483: next = com->next;
1484: DMCompositeCreate(comm,fine);
1486: /* loop over packed objects, handling one at at time */
1487: while (next) {
1488: DMCoarsen(next->dm,comm,&dm);
1489: DMCompositeAddDM(*fine,dm);
1490: PetscObjectDereference((PetscObject)dm);
1491: next = next->next;
1492: }
1493: return(0);
1494: }
1496: PetscErrorCode DMCreateInterpolation_Composite(DM coarse,DM fine,Mat *A,Vec *v)
1497: {
1498: PetscErrorCode ierr;
1499: PetscInt m,n,M,N,nDM,i;
1500: struct DMCompositeLink *nextc;
1501: struct DMCompositeLink *nextf;
1502: Vec gcoarse,gfine,*vecs;
1503: DM_Composite *comcoarse = (DM_Composite*)coarse->data;
1504: DM_Composite *comfine = (DM_Composite*)fine->data;
1505: Mat *mats;
1510: DMSetUp(coarse);
1511: DMSetUp(fine);
1512: /* use global vectors only for determining matrix layout */
1513: DMGetGlobalVector(coarse,&gcoarse);
1514: DMGetGlobalVector(fine,&gfine);
1515: VecGetLocalSize(gcoarse,&n);
1516: VecGetLocalSize(gfine,&m);
1517: VecGetSize(gcoarse,&N);
1518: VecGetSize(gfine,&M);
1519: DMRestoreGlobalVector(coarse,&gcoarse);
1520: DMRestoreGlobalVector(fine,&gfine);
1522: nDM = comfine->nDM;
1523: if (nDM != comcoarse->nDM) SETERRQ2(PetscObjectComm((PetscObject)fine),PETSC_ERR_ARG_INCOMP,"Fine DMComposite has %D entries, but coarse has %D",nDM,comcoarse->nDM);
1524: PetscCalloc1(nDM*nDM,&mats);
1525: if (v) {
1526: PetscCalloc1(nDM,&vecs);
1527: }
1529: /* loop over packed objects, handling one at at time */
1530: for (nextc=comcoarse->next,nextf=comfine->next,i=0; nextc; nextc=nextc->next,nextf=nextf->next,i++) {
1531: if (!v) {
1532: DMCreateInterpolation(nextc->dm,nextf->dm,&mats[i*nDM+i],NULL);
1533: } else {
1534: DMCreateInterpolation(nextc->dm,nextf->dm,&mats[i*nDM+i],&vecs[i]);
1535: }
1536: }
1537: MatCreateNest(PetscObjectComm((PetscObject)fine),nDM,NULL,nDM,NULL,mats,A);
1538: if (v) {
1539: VecCreateNest(PetscObjectComm((PetscObject)fine),nDM,NULL,vecs,v);
1540: }
1541: for (i=0; i<nDM*nDM; i++) {MatDestroy(&mats[i]);}
1542: PetscFree(mats);
1543: if (v) {
1544: for (i=0; i<nDM; i++) {VecDestroy(&vecs[i]);}
1545: PetscFree(vecs);
1546: }
1547: return(0);
1548: }
1550: static PetscErrorCode DMGetLocalToGlobalMapping_Composite(DM dm)
1551: {
1552: DM_Composite *com = (DM_Composite*)dm->data;
1553: ISLocalToGlobalMapping *ltogs;
1554: PetscInt i;
1555: PetscErrorCode ierr;
1558: /* Set the ISLocalToGlobalMapping on the new matrix */
1559: DMCompositeGetISLocalToGlobalMappings(dm,<ogs);
1560: ISLocalToGlobalMappingConcatenate(PetscObjectComm((PetscObject)dm),com->nDM,ltogs,&dm->ltogmap);
1561: for (i=0; i<com->nDM; i++) {ISLocalToGlobalMappingDestroy(<ogs[i]);}
1562: PetscFree(ltogs);
1563: return(0);
1564: }
1567: PetscErrorCode DMCreateColoring_Composite(DM dm,ISColoringType ctype,ISColoring *coloring)
1568: {
1569: PetscErrorCode ierr;
1570: PetscInt n,i,cnt;
1571: ISColoringValue *colors;
1572: PetscBool dense = PETSC_FALSE;
1573: ISColoringValue maxcol = 0;
1574: DM_Composite *com = (DM_Composite*)dm->data;
1578: if (ctype == IS_COLORING_LOCAL) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Only global coloring supported");
1579: else if (ctype == IS_COLORING_GLOBAL) {
1580: n = com->n;
1581: } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unknown ISColoringType");
1582: PetscMalloc1(n,&colors); /* freed in ISColoringDestroy() */
1584: PetscOptionsGetBool(((PetscObject)dm)->options,((PetscObject)dm)->prefix,"-dmcomposite_dense_jacobian",&dense,NULL);
1585: if (dense) {
1586: for (i=0; i<n; i++) {
1587: colors[i] = (ISColoringValue)(com->rstart + i);
1588: }
1589: maxcol = com->N;
1590: } else {
1591: struct DMCompositeLink *next = com->next;
1592: PetscMPIInt rank;
1594: MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);
1595: cnt = 0;
1596: while (next) {
1597: ISColoring lcoloring;
1599: DMCreateColoring(next->dm,IS_COLORING_GLOBAL,&lcoloring);
1600: for (i=0; i<lcoloring->N; i++) {
1601: colors[cnt++] = maxcol + lcoloring->colors[i];
1602: }
1603: maxcol += lcoloring->n;
1604: ISColoringDestroy(&lcoloring);
1605: next = next->next;
1606: }
1607: }
1608: ISColoringCreate(PetscObjectComm((PetscObject)dm),maxcol,n,colors,PETSC_OWN_POINTER,coloring);
1609: return(0);
1610: }
1612: PetscErrorCode DMGlobalToLocalBegin_Composite(DM dm,Vec gvec,InsertMode mode,Vec lvec)
1613: {
1614: PetscErrorCode ierr;
1615: struct DMCompositeLink *next;
1616: PetscScalar *garray,*larray;
1617: DM_Composite *com = (DM_Composite*)dm->data;
1623: if (!com->setup) {
1624: DMSetUp(dm);
1625: }
1627: VecGetArray(gvec,&garray);
1628: VecGetArray(lvec,&larray);
1630: /* loop over packed objects, handling one at at time */
1631: next = com->next;
1632: while (next) {
1633: Vec local,global;
1634: PetscInt N;
1636: DMGetGlobalVector(next->dm,&global);
1637: VecGetLocalSize(global,&N);
1638: VecPlaceArray(global,garray);
1639: DMGetLocalVector(next->dm,&local);
1640: VecPlaceArray(local,larray);
1641: DMGlobalToLocalBegin(next->dm,global,mode,local);
1642: DMGlobalToLocalEnd(next->dm,global,mode,local);
1643: VecResetArray(global);
1644: VecResetArray(local);
1645: DMRestoreGlobalVector(next->dm,&global);
1646: DMRestoreLocalVector(next->dm,&local);
1648: larray += next->nlocal;
1649: garray += next->n;
1650: next = next->next;
1651: }
1653: VecRestoreArray(gvec,NULL);
1654: VecRestoreArray(lvec,NULL);
1655: return(0);
1656: }
1658: PetscErrorCode DMGlobalToLocalEnd_Composite(DM dm,Vec gvec,InsertMode mode,Vec lvec)
1659: {
1664: return(0);
1665: }
1667: PetscErrorCode DMLocalToGlobalBegin_Composite(DM dm,Vec lvec,InsertMode mode,Vec gvec)
1668: {
1669: PetscErrorCode ierr;
1670: struct DMCompositeLink *next;
1671: PetscScalar *larray,*garray;
1672: DM_Composite *com = (DM_Composite*)dm->data;
1679: if (!com->setup) {
1680: DMSetUp(dm);
1681: }
1683: VecGetArray(lvec,&larray);
1684: VecGetArray(gvec,&garray);
1686: /* loop over packed objects, handling one at at time */
1687: next = com->next;
1688: while (next) {
1689: Vec global,local;
1691: DMGetLocalVector(next->dm,&local);
1692: VecPlaceArray(local,larray);
1693: DMGetGlobalVector(next->dm,&global);
1694: VecPlaceArray(global,garray);
1695: DMLocalToGlobalBegin(next->dm,local,mode,global);
1696: DMLocalToGlobalEnd(next->dm,local,mode,global);
1697: VecResetArray(local);
1698: VecResetArray(global);
1699: DMRestoreGlobalVector(next->dm,&global);
1700: DMRestoreLocalVector(next->dm,&local);
1702: garray += next->n;
1703: larray += next->nlocal;
1704: next = next->next;
1705: }
1707: VecRestoreArray(gvec,NULL);
1708: VecRestoreArray(lvec,NULL);
1710: return(0);
1711: }
1713: PetscErrorCode DMLocalToGlobalEnd_Composite(DM dm,Vec lvec,InsertMode mode,Vec gvec)
1714: {
1719: return(0);
1720: }
1722: PetscErrorCode DMLocalToLocalBegin_Composite(DM dm,Vec vec1,InsertMode mode,Vec vec2)
1723: {
1724: PetscErrorCode ierr;
1725: struct DMCompositeLink *next;
1726: PetscScalar *array1,*array2;
1727: DM_Composite *com = (DM_Composite*)dm->data;
1734: if (!com->setup) {
1735: DMSetUp(dm);
1736: }
1738: VecGetArray(vec1,&array1);
1739: VecGetArray(vec2,&array2);
1741: /* loop over packed objects, handling one at at time */
1742: next = com->next;
1743: while (next) {
1744: Vec local1,local2;
1746: DMGetLocalVector(next->dm,&local1);
1747: VecPlaceArray(local1,array1);
1748: DMGetLocalVector(next->dm,&local2);
1749: VecPlaceArray(local2,array2);
1750: DMLocalToLocalBegin(next->dm,local1,mode,local2);
1751: DMLocalToLocalEnd(next->dm,local1,mode,local2);
1752: VecResetArray(local2);
1753: DMRestoreLocalVector(next->dm,&local2);
1754: VecResetArray(local1);
1755: DMRestoreLocalVector(next->dm,&local1);
1757: array1 += next->nlocal;
1758: array2 += next->nlocal;
1759: next = next->next;
1760: }
1762: VecRestoreArray(vec1,NULL);
1763: VecRestoreArray(vec2,NULL);
1765: return(0);
1766: }
1768: PetscErrorCode DMLocalToLocalEnd_Composite(DM dm,Vec lvec,InsertMode mode,Vec gvec)
1769: {
1774: return(0);
1775: }
1777: /*MC
1778: DMCOMPOSITE = "composite" - A DM object that is used to manage data for a collection of DMs
1780: Level: intermediate
1782: .seealso: DMType, DM, DMDACreate(), DMCreate(), DMSetType(), DMCompositeCreate()
1783: M*/
1786: PETSC_EXTERN PetscErrorCode DMCreate_Composite(DM p)
1787: {
1789: DM_Composite *com;
1792: PetscNewLog(p,&com);
1793: p->data = com;
1794: PetscObjectChangeTypeName((PetscObject)p,"DMComposite");
1795: com->n = 0;
1796: com->nghost = 0;
1797: com->next = NULL;
1798: com->nDM = 0;
1800: p->ops->createglobalvector = DMCreateGlobalVector_Composite;
1801: p->ops->createlocalvector = DMCreateLocalVector_Composite;
1802: p->ops->getlocaltoglobalmapping = DMGetLocalToGlobalMapping_Composite;
1803: p->ops->createfieldis = DMCreateFieldIS_Composite;
1804: p->ops->createfielddecomposition = DMCreateFieldDecomposition_Composite;
1805: p->ops->refine = DMRefine_Composite;
1806: p->ops->coarsen = DMCoarsen_Composite;
1807: p->ops->createinterpolation = DMCreateInterpolation_Composite;
1808: p->ops->creatematrix = DMCreateMatrix_Composite;
1809: p->ops->getcoloring = DMCreateColoring_Composite;
1810: p->ops->globaltolocalbegin = DMGlobalToLocalBegin_Composite;
1811: p->ops->globaltolocalend = DMGlobalToLocalEnd_Composite;
1812: p->ops->localtoglobalbegin = DMLocalToGlobalBegin_Composite;
1813: p->ops->localtoglobalend = DMLocalToGlobalEnd_Composite;
1814: p->ops->localtolocalbegin = DMLocalToLocalBegin_Composite;
1815: p->ops->localtolocalend = DMLocalToLocalEnd_Composite;
1816: p->ops->destroy = DMDestroy_Composite;
1817: p->ops->view = DMView_Composite;
1818: p->ops->setup = DMSetUp_Composite;
1820: PetscObjectComposeFunction((PetscObject)p,"DMSetUpGLVisViewer_C",DMSetUpGLVisViewer_Composite);
1821: return(0);
1822: }
1824: /*@
1825: DMCompositeCreate - Creates a vector packer, used to generate "composite"
1826: vectors made up of several subvectors.
1828: Collective on MPI_Comm
1830: Input Parameter:
1831: . comm - the processors that will share the global vector
1833: Output Parameters:
1834: . packer - the packer object
1836: Level: advanced
1838: .seealso DMDestroy(), DMCompositeAddDM(), DMCompositeScatter(), DMCOMPOSITE,DMCreate()
1839: DMCompositeGather(), DMCreateGlobalVector(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess()
1840: DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(), DMCompositeGetEntries()
1842: @*/
1843: PetscErrorCode DMCompositeCreate(MPI_Comm comm,DM *packer)
1844: {
1849: DMCreate(comm,packer);
1850: DMSetType(*packer,DMCOMPOSITE);
1851: return(0);
1852: }