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,&ltog);
    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,&ltogs);
    1569:   ISLocalToGlobalMappingConcatenate(PetscObjectComm((PetscObject)dm),com->nDM,ltogs,&dm->ltogmap);
    1570:   for (i=0; i<com->nDM; i++) {ISLocalToGlobalMappingDestroy(&ltogs[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: }