#include <../src/dm/impls/composite/packimpl.h> /*I "petscdmcomposite.h" I*/ #include #include #include /*@C DMCompositeSetCoupling - Sets user provided routines that compute the coupling between the separate components (DMs) in a DMto build the correct matrix nonzero structure. Logically Collective Input Parameter: + dm - the composite object - formcouplelocations - routine to set the nonzero locations in the matrix Level: advanced Not available from Fortran Notes: See DMSetApplicationContext() and DMGetApplicationContext() for how to get user information into this routine @*/ PetscErrorCode DMCompositeSetCoupling(DM dm,PetscErrorCode (*FormCoupleLocations)(DM,Mat,PetscInt*,PetscInt*,PetscInt,PetscInt,PetscInt,PetscInt)) { DM_Composite *com = (DM_Composite*)dm->data; PetscBool flg; PetscErrorCode ierr; PetscFunctionBegin; ierr = PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);CHKERRQ(ierr); if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name); com->FormCoupleLocations = FormCoupleLocations; PetscFunctionReturn(0); } PetscErrorCode DMDestroy_Composite(DM dm) { PetscErrorCode ierr; struct DMCompositeLink *next, *prev; DM_Composite *com = (DM_Composite*)dm->data; PetscFunctionBegin; next = com->next; while (next) { prev = next; next = next->next; ierr = DMDestroy(&prev->dm);CHKERRQ(ierr); ierr = PetscFree(prev->grstarts);CHKERRQ(ierr); ierr = PetscFree(prev);CHKERRQ(ierr); } ierr = PetscObjectComposeFunction((PetscObject)dm,"DMSetUpGLVisViewer_C",NULL);CHKERRQ(ierr); /* This was originally freed in DMDestroy(), but that prevents reference counting of backend objects */ ierr = PetscFree(com);CHKERRQ(ierr); PetscFunctionReturn(0); } PetscErrorCode DMView_Composite(DM dm,PetscViewer v) { PetscErrorCode ierr; PetscBool iascii; DM_Composite *com = (DM_Composite*)dm->data; PetscFunctionBegin; ierr = PetscObjectTypeCompare((PetscObject)v,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); if (iascii) { struct DMCompositeLink *lnk = com->next; PetscInt i; ierr = PetscViewerASCIIPrintf(v,"DM (%s)\n",((PetscObject)dm)->prefix ? ((PetscObject)dm)->prefix : "no prefix");CHKERRQ(ierr); ierr = PetscViewerASCIIPrintf(v," contains %D DMs\n",com->nDM);CHKERRQ(ierr); ierr = PetscViewerASCIIPushTab(v);CHKERRQ(ierr); for (i=0; lnk; lnk=lnk->next,i++) { ierr = PetscViewerASCIIPrintf(v,"Link %D: DM of type %s\n",i,((PetscObject)lnk->dm)->type_name);CHKERRQ(ierr); ierr = PetscViewerASCIIPushTab(v);CHKERRQ(ierr); ierr = DMView(lnk->dm,v);CHKERRQ(ierr); ierr = PetscViewerASCIIPopTab(v);CHKERRQ(ierr); } ierr = PetscViewerASCIIPopTab(v);CHKERRQ(ierr); } PetscFunctionReturn(0); } /* --------------------------------------------------------------------------------------*/ PetscErrorCode DMSetUp_Composite(DM dm) { PetscErrorCode ierr; PetscInt nprev = 0; PetscMPIInt rank,size; DM_Composite *com = (DM_Composite*)dm->data; struct DMCompositeLink *next = com->next; PetscLayout map; PetscFunctionBegin; if (com->setup) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Packer has already been setup"); ierr = PetscLayoutCreate(PetscObjectComm((PetscObject)dm),&map);CHKERRQ(ierr); ierr = PetscLayoutSetLocalSize(map,com->n);CHKERRQ(ierr); ierr = PetscLayoutSetSize(map,PETSC_DETERMINE);CHKERRQ(ierr); ierr = PetscLayoutSetBlockSize(map,1);CHKERRQ(ierr); ierr = PetscLayoutSetUp(map);CHKERRQ(ierr); ierr = PetscLayoutGetSize(map,&com->N);CHKERRQ(ierr); ierr = PetscLayoutGetRange(map,&com->rstart,NULL);CHKERRQ(ierr); ierr = PetscLayoutDestroy(&map);CHKERRQ(ierr); /* now set the rstart for each linked vector */ ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);CHKERRQ(ierr); ierr = MPI_Comm_size(PetscObjectComm((PetscObject)dm),&size);CHKERRQ(ierr); while (next) { next->rstart = nprev; nprev += next->n; next->grstart = com->rstart + next->rstart; ierr = PetscMalloc1(size,&next->grstarts);CHKERRQ(ierr); ierr = MPI_Allgather(&next->grstart,1,MPIU_INT,next->grstarts,1,MPIU_INT,PetscObjectComm((PetscObject)dm));CHKERRQ(ierr); next = next->next; } com->setup = PETSC_TRUE; PetscFunctionReturn(0); } /* ----------------------------------------------------------------------------------*/ /*@ DMCompositeGetNumberDM - Get's the number of DM objects in the DMComposite representation. Not Collective Input Parameter: . dm - the packer object Output Parameter: . nDM - the number of DMs Level: beginner @*/ PetscErrorCode DMCompositeGetNumberDM(DM dm,PetscInt *nDM) { DM_Composite *com = (DM_Composite*)dm->data; PetscBool flg; PetscErrorCode ierr; PetscFunctionBegin; PetscValidHeaderSpecific(dm,DM_CLASSID,1); ierr = PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);CHKERRQ(ierr); if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name); *nDM = com->nDM; PetscFunctionReturn(0); } /*@C DMCompositeGetAccess - Allows one to access the individual packed vectors in their global representation. Collective on dm Input Parameters: + dm - the packer object - gvec - the global vector Output Parameters: . Vec* ... - the packed parallel vectors, NULL for those that are not needed Notes: Use DMCompositeRestoreAccess() to return the vectors when you no longer need them Fortran Notes: Fortran callers must use numbered versions of this routine, e.g., DMCompositeGetAccess4(dm,gvec,vec1,vec2,vec3,vec4) or use the alternative interface DMCompositeGetAccessArray(). Level: advanced .seealso: DMCompositeGetEntries(), DMCompositeScatter() @*/ PetscErrorCode DMCompositeGetAccess(DM dm,Vec gvec,...) { va_list Argp; PetscErrorCode ierr; struct DMCompositeLink *next; DM_Composite *com = (DM_Composite*)dm->data; PetscInt readonly; PetscBool flg; PetscFunctionBegin; PetscValidHeaderSpecific(dm,DM_CLASSID,1); PetscValidHeaderSpecific(gvec,VEC_CLASSID,2); ierr = PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);CHKERRQ(ierr); if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name); next = com->next; if (!com->setup) { ierr = DMSetUp(dm);CHKERRQ(ierr); } ierr = VecLockGet(gvec,&readonly);CHKERRQ(ierr); /* loop over packed objects, handling one at at time */ va_start(Argp,gvec); while (next) { Vec *vec; vec = va_arg(Argp, Vec*); if (vec) { ierr = DMGetGlobalVector(next->dm,vec);CHKERRQ(ierr); if (readonly) { const PetscScalar *array; ierr = VecGetArrayRead(gvec,&array);CHKERRQ(ierr); ierr = VecPlaceArray(*vec,array+next->rstart);CHKERRQ(ierr); ierr = VecLockReadPush(*vec);CHKERRQ(ierr); ierr = VecRestoreArrayRead(gvec,&array);CHKERRQ(ierr); } else { PetscScalar *array; ierr = VecGetArray(gvec,&array);CHKERRQ(ierr); ierr = VecPlaceArray(*vec,array+next->rstart);CHKERRQ(ierr); ierr = VecRestoreArray(gvec,&array);CHKERRQ(ierr); } } next = next->next; } va_end(Argp); PetscFunctionReturn(0); } /*@C DMCompositeGetAccessArray - Allows one to access the individual packed vectors in their global representation. Collective on dm Input Parameters: + dm - the packer object . pvec - packed vector . nwanted - number of vectors wanted - wanted - sorted array of vectors wanted, or NULL to get all vectors Output Parameters: . vecs - array of requested global vectors (must be allocated) Notes: Use DMCompositeRestoreAccessArray() to return the vectors when you no longer need them Level: advanced .seealso: DMCompositeGetAccess(), DMCompositeGetEntries(), DMCompositeScatter(), DMCompositeGather() @*/ PetscErrorCode DMCompositeGetAccessArray(DM dm,Vec pvec,PetscInt nwanted,const PetscInt *wanted,Vec *vecs) { PetscErrorCode ierr; struct DMCompositeLink *link; PetscInt i,wnum; DM_Composite *com = (DM_Composite*)dm->data; PetscInt readonly; PetscBool flg; PetscFunctionBegin; PetscValidHeaderSpecific(dm,DM_CLASSID,1); PetscValidHeaderSpecific(pvec,VEC_CLASSID,2); ierr = PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);CHKERRQ(ierr); if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name); if (!com->setup) { ierr = DMSetUp(dm);CHKERRQ(ierr); } ierr = VecLockGet(pvec,&readonly);CHKERRQ(ierr); for (i=0,wnum=0,link=com->next; link && wnumnext) { if (!wanted || i == wanted[wnum]) { Vec v; ierr = DMGetGlobalVector(link->dm,&v);CHKERRQ(ierr); if (readonly) { const PetscScalar *array; ierr = VecGetArrayRead(pvec,&array);CHKERRQ(ierr); ierr = VecPlaceArray(v,array+link->rstart);CHKERRQ(ierr); ierr = VecLockReadPush(v);CHKERRQ(ierr); ierr = VecRestoreArrayRead(pvec,&array);CHKERRQ(ierr); } else { PetscScalar *array; ierr = VecGetArray(pvec,&array);CHKERRQ(ierr); ierr = VecPlaceArray(v,array+link->rstart);CHKERRQ(ierr); ierr = VecRestoreArray(pvec,&array);CHKERRQ(ierr); } vecs[wnum++] = v; } } PetscFunctionReturn(0); } /*@C DMCompositeGetLocalAccessArray - Allows one to access the individual packed vectors in their local representation. Collective on dm. Input Parameters: + dm - the packer object . pvec - packed vector . nwanted - number of vectors wanted - wanted - sorted array of vectors wanted, or NULL to get all vectors Output Parameters: . vecs - array of requested local vectors (must be allocated) Notes: Use DMCompositeRestoreLocalAccessArray() to return the vectors when you no longer need them. Level: advanced .seealso: DMCompositeRestoreLocalAccessArray(), DMCompositeGetAccess(), DMCompositeGetEntries(), DMCompositeScatter(), DMCompositeGather() @*/ PetscErrorCode DMCompositeGetLocalAccessArray(DM dm,Vec pvec,PetscInt nwanted,const PetscInt *wanted,Vec *vecs) { PetscErrorCode ierr; struct DMCompositeLink *link; PetscInt i,wnum; DM_Composite *com = (DM_Composite*)dm->data; PetscInt readonly; PetscInt nlocal = 0; PetscBool flg; PetscFunctionBegin; PetscValidHeaderSpecific(dm,DM_CLASSID,1); PetscValidHeaderSpecific(pvec,VEC_CLASSID,2); ierr = PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);CHKERRQ(ierr); if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name); if (!com->setup) { ierr = DMSetUp(dm);CHKERRQ(ierr); } ierr = VecLockGet(pvec,&readonly);CHKERRQ(ierr); for (i=0,wnum=0,link=com->next; link && wnumnext) { if (!wanted || i == wanted[wnum]) { Vec v; ierr = DMGetLocalVector(link->dm,&v);CHKERRQ(ierr); if (readonly) { const PetscScalar *array; ierr = VecGetArrayRead(pvec,&array);CHKERRQ(ierr); ierr = VecPlaceArray(v,array+nlocal);CHKERRQ(ierr); ierr = VecLockReadPush(v);CHKERRQ(ierr); ierr = VecRestoreArrayRead(pvec,&array);CHKERRQ(ierr); } else { PetscScalar *array; ierr = VecGetArray(pvec,&array);CHKERRQ(ierr); ierr = VecPlaceArray(v,array+nlocal);CHKERRQ(ierr); ierr = VecRestoreArray(pvec,&array);CHKERRQ(ierr); } vecs[wnum++] = v; } nlocal += link->nlocal; } PetscFunctionReturn(0); } /*@C DMCompositeRestoreAccess - Returns the vectors obtained with DMCompositeGetAccess() representation. Collective on dm Input Parameters: + dm - the packer object . gvec - the global vector - Vec* ... - the individual parallel vectors, NULL for those that are not needed Level: advanced .seealso DMCompositeAddDM(), DMCreateGlobalVector(), DMCompositeGather(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeScatter(), DMCompositeRestoreAccess(), DMCompositeGetAccess() @*/ PetscErrorCode DMCompositeRestoreAccess(DM dm,Vec gvec,...) { va_list Argp; PetscErrorCode ierr; struct DMCompositeLink *next; DM_Composite *com = (DM_Composite*)dm->data; PetscInt readonly; PetscBool flg; PetscFunctionBegin; PetscValidHeaderSpecific(dm,DM_CLASSID,1); PetscValidHeaderSpecific(gvec,VEC_CLASSID,2); ierr = PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);CHKERRQ(ierr); if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name); next = com->next; if (!com->setup) { ierr = DMSetUp(dm);CHKERRQ(ierr); } ierr = VecLockGet(gvec,&readonly);CHKERRQ(ierr); /* loop over packed objects, handling one at at time */ va_start(Argp,gvec); while (next) { Vec *vec; vec = va_arg(Argp, Vec*); if (vec) { ierr = VecResetArray(*vec);CHKERRQ(ierr); if (readonly) { ierr = VecLockReadPop(*vec);CHKERRQ(ierr); } ierr = DMRestoreGlobalVector(next->dm,vec);CHKERRQ(ierr); } next = next->next; } va_end(Argp); PetscFunctionReturn(0); } /*@C DMCompositeRestoreAccessArray - Returns the vectors obtained with DMCompositeGetAccessArray() Collective on dm Input Parameters: + dm - the packer object . pvec - packed vector . nwanted - number of vectors wanted . wanted - sorted array of vectors wanted, or NULL to get all vectors - vecs - array of global vectors to return Level: advanced .seealso: DMCompositeRestoreAccess(), DMCompositeRestoreEntries(), DMCompositeScatter(), DMCompositeGather() @*/ PetscErrorCode DMCompositeRestoreAccessArray(DM dm,Vec pvec,PetscInt nwanted,const PetscInt *wanted,Vec *vecs) { PetscErrorCode ierr; struct DMCompositeLink *link; PetscInt i,wnum; DM_Composite *com = (DM_Composite*)dm->data; PetscInt readonly; PetscBool flg; PetscFunctionBegin; PetscValidHeaderSpecific(dm,DM_CLASSID,1); PetscValidHeaderSpecific(pvec,VEC_CLASSID,2); ierr = PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);CHKERRQ(ierr); if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name); if (!com->setup) { ierr = DMSetUp(dm);CHKERRQ(ierr); } ierr = VecLockGet(pvec,&readonly);CHKERRQ(ierr); for (i=0,wnum=0,link=com->next; link && wnumnext) { if (!wanted || i == wanted[wnum]) { ierr = VecResetArray(vecs[wnum]);CHKERRQ(ierr); if (readonly) { ierr = VecLockReadPop(vecs[wnum]);CHKERRQ(ierr); } ierr = DMRestoreGlobalVector(link->dm,&vecs[wnum]);CHKERRQ(ierr); wnum++; } } PetscFunctionReturn(0); } /*@C DMCompositeRestoreLocalAccessArray - Returns the vectors obtained with DMCompositeGetLocalAccessArray(). Collective on dm. Input Parameters: + dm - the packer object . pvec - packed vector . nwanted - number of vectors wanted . wanted - sorted array of vectors wanted, or NULL to restore all vectors - vecs - array of local vectors to return Level: advanced Notes: nwanted and wanted must match the values given to DMCompositeGetLocalAccessArray() otherwise the call will fail. .seealso: DMCompositeGetLocalAccessArray(), DMCompositeRestoreAccessArray(), DMCompositeRestoreAccess(), DMCompositeRestoreEntries(), DMCompositeScatter(), DMCompositeGather() @*/ PetscErrorCode DMCompositeRestoreLocalAccessArray(DM dm,Vec pvec,PetscInt nwanted,const PetscInt *wanted,Vec *vecs) { PetscErrorCode ierr; struct DMCompositeLink *link; PetscInt i,wnum; DM_Composite *com = (DM_Composite*)dm->data; PetscInt readonly; PetscBool flg; PetscFunctionBegin; PetscValidHeaderSpecific(dm,DM_CLASSID,1); PetscValidHeaderSpecific(pvec,VEC_CLASSID,2); ierr = PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);CHKERRQ(ierr); if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name); if (!com->setup) { ierr = DMSetUp(dm);CHKERRQ(ierr); } ierr = VecLockGet(pvec,&readonly);CHKERRQ(ierr); for (i=0,wnum=0,link=com->next; link && wnumnext) { if (!wanted || i == wanted[wnum]) { ierr = VecResetArray(vecs[wnum]);CHKERRQ(ierr); if (readonly) { ierr = VecLockReadPop(vecs[wnum]);CHKERRQ(ierr); } ierr = DMRestoreLocalVector(link->dm,&vecs[wnum]);CHKERRQ(ierr); wnum++; } } PetscFunctionReturn(0); } /*@C DMCompositeScatter - Scatters from a global packed vector into its individual local vectors Collective on dm Input Parameters: + dm - the packer object . gvec - the global vector - Vec ... - the individual sequential vectors, NULL for those that are not needed Level: advanced Notes: DMCompositeScatterArray() is a non-variadic alternative that is often more convenient for library callers and is accessible from Fortran. .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(), DMCompositeGather(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(), DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(), DMCompositeGetEntries() DMCompositeScatterArray() @*/ PetscErrorCode DMCompositeScatter(DM dm,Vec gvec,...) { va_list Argp; PetscErrorCode ierr; struct DMCompositeLink *next; PetscInt cnt; DM_Composite *com = (DM_Composite*)dm->data; PetscBool flg; PetscFunctionBegin; PetscValidHeaderSpecific(dm,DM_CLASSID,1); PetscValidHeaderSpecific(gvec,VEC_CLASSID,2); ierr = PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);CHKERRQ(ierr); if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name); if (!com->setup) { ierr = DMSetUp(dm);CHKERRQ(ierr); } /* loop over packed objects, handling one at at time */ va_start(Argp,gvec); for (cnt=3,next=com->next; next; cnt++,next=next->next) { Vec local; local = va_arg(Argp, Vec); if (local) { Vec global; const PetscScalar *array; PetscValidHeaderSpecific(local,VEC_CLASSID,cnt); ierr = DMGetGlobalVector(next->dm,&global);CHKERRQ(ierr); ierr = VecGetArrayRead(gvec,&array);CHKERRQ(ierr); ierr = VecPlaceArray(global,array+next->rstart);CHKERRQ(ierr); ierr = DMGlobalToLocalBegin(next->dm,global,INSERT_VALUES,local);CHKERRQ(ierr); ierr = DMGlobalToLocalEnd(next->dm,global,INSERT_VALUES,local);CHKERRQ(ierr); ierr = VecRestoreArrayRead(gvec,&array);CHKERRQ(ierr); ierr = VecResetArray(global);CHKERRQ(ierr); ierr = DMRestoreGlobalVector(next->dm,&global);CHKERRQ(ierr); } } va_end(Argp); PetscFunctionReturn(0); } /*@ DMCompositeScatterArray - Scatters from a global packed vector into its individual local vectors Collective on dm Input Parameters: + dm - the packer object . gvec - the global vector - lvecs - array of local vectors, NULL for any that are not needed Level: advanced Note: This is a non-variadic alternative to DMCompositeScatter() .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector() DMCompositeGather(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(), DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(), DMCompositeGetEntries() @*/ PetscErrorCode DMCompositeScatterArray(DM dm,Vec gvec,Vec *lvecs) { PetscErrorCode ierr; struct DMCompositeLink *next; PetscInt i; DM_Composite *com = (DM_Composite*)dm->data; PetscBool flg; PetscFunctionBegin; PetscValidHeaderSpecific(dm,DM_CLASSID,1); PetscValidHeaderSpecific(gvec,VEC_CLASSID,2); ierr = PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);CHKERRQ(ierr); if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name); if (!com->setup) { ierr = DMSetUp(dm);CHKERRQ(ierr); } /* loop over packed objects, handling one at at time */ for (i=0,next=com->next; next; next=next->next,i++) { if (lvecs[i]) { Vec global; const PetscScalar *array; PetscValidHeaderSpecific(lvecs[i],VEC_CLASSID,3); ierr = DMGetGlobalVector(next->dm,&global);CHKERRQ(ierr); ierr = VecGetArrayRead(gvec,&array);CHKERRQ(ierr); ierr = VecPlaceArray(global,(PetscScalar*)array+next->rstart);CHKERRQ(ierr); ierr = DMGlobalToLocalBegin(next->dm,global,INSERT_VALUES,lvecs[i]);CHKERRQ(ierr); ierr = DMGlobalToLocalEnd(next->dm,global,INSERT_VALUES,lvecs[i]);CHKERRQ(ierr); ierr = VecRestoreArrayRead(gvec,&array);CHKERRQ(ierr); ierr = VecResetArray(global);CHKERRQ(ierr); ierr = DMRestoreGlobalVector(next->dm,&global);CHKERRQ(ierr); } } PetscFunctionReturn(0); } /*@C DMCompositeGather - Gathers into a global packed vector from its individual local vectors Collective on dm Input Parameter: + dm - the packer object . gvec - the global vector . imode - INSERT_VALUES or ADD_VALUES - Vec ... - the individual sequential vectors, NULL for any that are not needed Level: advanced Not available from Fortran, Fortran users can use DMCompositeGatherArray() .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(), DMCompositeScatter(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(), DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(), DMCompositeGetEntries() @*/ PetscErrorCode DMCompositeGather(DM dm,InsertMode imode,Vec gvec,...) { va_list Argp; PetscErrorCode ierr; struct DMCompositeLink *next; DM_Composite *com = (DM_Composite*)dm->data; PetscInt cnt; PetscBool flg; PetscFunctionBegin; PetscValidHeaderSpecific(dm,DM_CLASSID,1); PetscValidHeaderSpecific(gvec,VEC_CLASSID,2); ierr = PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);CHKERRQ(ierr); if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name); if (!com->setup) { ierr = DMSetUp(dm);CHKERRQ(ierr); } /* loop over packed objects, handling one at at time */ va_start(Argp,gvec); for (cnt=3,next=com->next; next; cnt++,next=next->next) { Vec local; local = va_arg(Argp, Vec); if (local) { PetscScalar *array; Vec global; PetscValidHeaderSpecific(local,VEC_CLASSID,cnt); ierr = DMGetGlobalVector(next->dm,&global);CHKERRQ(ierr); ierr = VecGetArray(gvec,&array);CHKERRQ(ierr); ierr = VecPlaceArray(global,array+next->rstart);CHKERRQ(ierr); ierr = DMLocalToGlobalBegin(next->dm,local,imode,global);CHKERRQ(ierr); ierr = DMLocalToGlobalEnd(next->dm,local,imode,global);CHKERRQ(ierr); ierr = VecRestoreArray(gvec,&array);CHKERRQ(ierr); ierr = VecResetArray(global);CHKERRQ(ierr); ierr = DMRestoreGlobalVector(next->dm,&global);CHKERRQ(ierr); } } va_end(Argp); PetscFunctionReturn(0); } /*@ DMCompositeGatherArray - Gathers into a global packed vector from its individual local vectors Collective on dm Input Parameter: + dm - the packer object . gvec - the global vector . imode - INSERT_VALUES or ADD_VALUES - lvecs - the individual sequential vectors, NULL for any that are not needed Level: advanced Notes: This is a non-variadic alternative to DMCompositeGather(). .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(), DMCompositeScatter(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(), DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(), DMCompositeGetEntries(), @*/ PetscErrorCode DMCompositeGatherArray(DM dm,InsertMode imode,Vec gvec,Vec *lvecs) { PetscErrorCode ierr; struct DMCompositeLink *next; DM_Composite *com = (DM_Composite*)dm->data; PetscInt i; PetscBool flg; PetscFunctionBegin; PetscValidHeaderSpecific(dm,DM_CLASSID,1); PetscValidHeaderSpecific(gvec,VEC_CLASSID,2); ierr = PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);CHKERRQ(ierr); if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name); if (!com->setup) { ierr = DMSetUp(dm);CHKERRQ(ierr); } /* loop over packed objects, handling one at at time */ for (next=com->next,i=0; next; next=next->next,i++) { if (lvecs[i]) { PetscScalar *array; Vec global; PetscValidHeaderSpecific(lvecs[i],VEC_CLASSID,3); ierr = DMGetGlobalVector(next->dm,&global);CHKERRQ(ierr); ierr = VecGetArray(gvec,&array);CHKERRQ(ierr); ierr = VecPlaceArray(global,array+next->rstart);CHKERRQ(ierr); ierr = DMLocalToGlobalBegin(next->dm,lvecs[i],imode,global);CHKERRQ(ierr); ierr = DMLocalToGlobalEnd(next->dm,lvecs[i],imode,global);CHKERRQ(ierr); ierr = VecRestoreArray(gvec,&array);CHKERRQ(ierr); ierr = VecResetArray(global);CHKERRQ(ierr); ierr = DMRestoreGlobalVector(next->dm,&global);CHKERRQ(ierr); } } PetscFunctionReturn(0); } /*@ DMCompositeAddDM - adds a DM vector to a DMComposite Collective on dm Input Parameter: + dmc - the DMComposite (packer) object - dm - the DM object Level: advanced .seealso DMDestroy(), DMCompositeGather(), DMCompositeAddDM(), DMCreateGlobalVector(), DMCompositeScatter(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(), DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(), DMCompositeGetEntries() @*/ PetscErrorCode DMCompositeAddDM(DM dmc,DM dm) { PetscErrorCode ierr; PetscInt n,nlocal; struct DMCompositeLink *mine,*next; Vec global,local; DM_Composite *com = (DM_Composite*)dmc->data; PetscBool flg; PetscFunctionBegin; PetscValidHeaderSpecific(dmc,DM_CLASSID,1); PetscValidHeaderSpecific(dm,DM_CLASSID,2); ierr = PetscObjectTypeCompare((PetscObject)dmc,DMCOMPOSITE,&flg);CHKERRQ(ierr); if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name); next = com->next; if (com->setup) SETERRQ(PetscObjectComm((PetscObject)dmc),PETSC_ERR_ARG_WRONGSTATE,"Cannot add a DM once you have used the DMComposite"); /* create new link */ ierr = PetscNew(&mine);CHKERRQ(ierr); ierr = PetscObjectReference((PetscObject)dm);CHKERRQ(ierr); ierr = DMGetGlobalVector(dm,&global);CHKERRQ(ierr); ierr = VecGetLocalSize(global,&n);CHKERRQ(ierr); ierr = DMRestoreGlobalVector(dm,&global);CHKERRQ(ierr); ierr = DMGetLocalVector(dm,&local);CHKERRQ(ierr); ierr = VecGetSize(local,&nlocal);CHKERRQ(ierr); ierr = DMRestoreLocalVector(dm,&local);CHKERRQ(ierr); mine->n = n; mine->nlocal = nlocal; mine->dm = dm; mine->next = NULL; com->n += n; com->nghost += nlocal; /* add to end of list */ if (!next) com->next = mine; else { while (next->next) next = next->next; next->next = mine; } com->nDM++; com->nmine++; PetscFunctionReturn(0); } #include PETSC_EXTERN PetscErrorCode VecView_MPI(Vec,PetscViewer); PetscErrorCode VecView_DMComposite(Vec gvec,PetscViewer viewer) { DM dm; PetscErrorCode ierr; struct DMCompositeLink *next; PetscBool isdraw; DM_Composite *com; PetscFunctionBegin; ierr = VecGetDM(gvec, &dm);CHKERRQ(ierr); if (!dm) SETERRQ(PetscObjectComm((PetscObject)gvec),PETSC_ERR_ARG_WRONG,"Vector not generated from a DMComposite"); com = (DM_Composite*)dm->data; next = com->next; ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);CHKERRQ(ierr); if (!isdraw) { /* do I really want to call this? */ ierr = VecView_MPI(gvec,viewer);CHKERRQ(ierr); } else { PetscInt cnt = 0; /* loop over packed objects, handling one at at time */ while (next) { Vec vec; const PetscScalar *array; PetscInt bs; /* Should use VecGetSubVector() eventually, but would need to forward the DM for that to work */ ierr = DMGetGlobalVector(next->dm,&vec);CHKERRQ(ierr); ierr = VecGetArrayRead(gvec,&array);CHKERRQ(ierr); ierr = VecPlaceArray(vec,(PetscScalar*)array+next->rstart);CHKERRQ(ierr); ierr = VecRestoreArrayRead(gvec,&array);CHKERRQ(ierr); ierr = VecView(vec,viewer);CHKERRQ(ierr); ierr = VecResetArray(vec);CHKERRQ(ierr); ierr = VecGetBlockSize(vec,&bs);CHKERRQ(ierr); ierr = DMRestoreGlobalVector(next->dm,&vec);CHKERRQ(ierr); ierr = PetscViewerDrawBaseAdd(viewer,bs);CHKERRQ(ierr); cnt += bs; next = next->next; } ierr = PetscViewerDrawBaseAdd(viewer,-cnt);CHKERRQ(ierr); } PetscFunctionReturn(0); } PetscErrorCode DMCreateGlobalVector_Composite(DM dm,Vec *gvec) { PetscErrorCode ierr; DM_Composite *com = (DM_Composite*)dm->data; PetscFunctionBegin; PetscValidHeaderSpecific(dm,DM_CLASSID,1); ierr = DMSetFromOptions(dm);CHKERRQ(ierr); ierr = DMSetUp(dm);CHKERRQ(ierr); ierr = VecCreate(PetscObjectComm((PetscObject)dm),gvec);CHKERRQ(ierr); ierr = VecSetType(*gvec,dm->vectype);CHKERRQ(ierr); ierr = VecSetSizes(*gvec,com->n,com->N);CHKERRQ(ierr); ierr = VecSetDM(*gvec, dm);CHKERRQ(ierr); ierr = VecSetOperation(*gvec,VECOP_VIEW,(void (*)(void))VecView_DMComposite);CHKERRQ(ierr); PetscFunctionReturn(0); } PetscErrorCode DMCreateLocalVector_Composite(DM dm,Vec *lvec) { PetscErrorCode ierr; DM_Composite *com = (DM_Composite*)dm->data; PetscFunctionBegin; PetscValidHeaderSpecific(dm,DM_CLASSID,1); if (!com->setup) { ierr = DMSetFromOptions(dm);CHKERRQ(ierr); ierr = DMSetUp(dm);CHKERRQ(ierr); } ierr = VecCreate(PETSC_COMM_SELF,lvec);CHKERRQ(ierr); ierr = VecSetType(*lvec,dm->vectype);CHKERRQ(ierr); ierr = VecSetSizes(*lvec,com->nghost,PETSC_DECIDE);CHKERRQ(ierr); ierr = VecSetDM(*lvec, dm);CHKERRQ(ierr); PetscFunctionReturn(0); } /*@C DMCompositeGetISLocalToGlobalMappings - gets an ISLocalToGlobalMapping for each DM in the DMComposite, maps to the composite global space Collective on DM Input Parameter: . dm - the packer object Output Parameters: . ltogs - the individual mappings for each packed vector. Note that this includes all the ghost points that individual ghosted DMDA's may have. Level: advanced Notes: Each entry of ltogs should be destroyed with ISLocalToGlobalMappingDestroy(), the ltogs array should be freed with PetscFree(). Not available from Fortran .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(), DMCompositeGather(), DMCompositeCreate(), DMCompositeGetAccess(), DMCompositeScatter(), DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(),DMCompositeGetEntries() @*/ PetscErrorCode DMCompositeGetISLocalToGlobalMappings(DM dm,ISLocalToGlobalMapping **ltogs) { PetscErrorCode ierr; PetscInt i,*idx,n,cnt; struct DMCompositeLink *next; PetscMPIInt rank; DM_Composite *com = (DM_Composite*)dm->data; PetscBool flg; PetscFunctionBegin; PetscValidHeaderSpecific(dm,DM_CLASSID,1); ierr = PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);CHKERRQ(ierr); if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name); ierr = DMSetUp(dm);CHKERRQ(ierr); ierr = PetscMalloc1(com->nDM,ltogs);CHKERRQ(ierr); next = com->next; ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);CHKERRQ(ierr); /* loop over packed objects, handling one at at time */ cnt = 0; while (next) { ISLocalToGlobalMapping ltog; PetscMPIInt size; const PetscInt *suboff,*indices; Vec global; /* Get sub-DM global indices for each local dof */ ierr = DMGetLocalToGlobalMapping(next->dm,<og);CHKERRQ(ierr); ierr = ISLocalToGlobalMappingGetSize(ltog,&n);CHKERRQ(ierr); ierr = ISLocalToGlobalMappingGetIndices(ltog,&indices);CHKERRQ(ierr); ierr = PetscMalloc1(n,&idx);CHKERRQ(ierr); /* Get the offsets for the sub-DM global vector */ ierr = DMGetGlobalVector(next->dm,&global);CHKERRQ(ierr); ierr = VecGetOwnershipRanges(global,&suboff);CHKERRQ(ierr); ierr = MPI_Comm_size(PetscObjectComm((PetscObject)global),&size);CHKERRQ(ierr); /* Shift the sub-DM definition of the global space to the composite global space */ for (i=0; igrstarts[rank]; continue; } /* Binary search to find which rank owns subi */ while (hi-lo > 1) { t = lo + (hi-lo)/2; if (suboff[t] > subi) hi = t; else lo = t; } idx[i] = subi - suboff[lo] + next->grstarts[lo]; } ierr = ISLocalToGlobalMappingRestoreIndices(ltog,&indices);CHKERRQ(ierr); ierr = ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)dm),1,n,idx,PETSC_OWN_POINTER,&(*ltogs)[cnt]);CHKERRQ(ierr); ierr = DMRestoreGlobalVector(next->dm,&global);CHKERRQ(ierr); next = next->next; cnt++; } PetscFunctionReturn(0); } /*@C DMCompositeGetLocalISs - Gets index sets for each component of a composite local vector Not Collective Input Arguments: . dm - composite DM Output Arguments: . is - array of serial index sets for each each component of the DMComposite Level: intermediate Notes: At present, a composite local vector does not normally exist. This function is used to provide index sets for MatGetLocalSubMatrix(). In the future, the scatters for each entry in the DMComposite may be be merged into a single scatter to a composite local vector. The user should not typically need to know which is being done. To get the composite global indices at all local points (including ghosts), use DMCompositeGetISLocalToGlobalMappings(). To get index sets for pieces of the composite global vector, use DMCompositeGetGlobalISs(). Each returned IS should be destroyed with ISDestroy(), the array should be freed with PetscFree(). Not available from Fortran .seealso: DMCompositeGetGlobalISs(), DMCompositeGetISLocalToGlobalMappings(), MatGetLocalSubMatrix(), MatCreateLocalRef() @*/ PetscErrorCode DMCompositeGetLocalISs(DM dm,IS **is) { PetscErrorCode ierr; DM_Composite *com = (DM_Composite*)dm->data; struct DMCompositeLink *link; PetscInt cnt,start; PetscBool flg; PetscFunctionBegin; PetscValidHeaderSpecific(dm,DM_CLASSID,1); PetscValidPointer(is,2); ierr = PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);CHKERRQ(ierr); if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name); ierr = PetscMalloc1(com->nmine,is);CHKERRQ(ierr); for (cnt=0,start=0,link=com->next; link; start+=link->nlocal,cnt++,link=link->next) { PetscInt bs; ierr = ISCreateStride(PETSC_COMM_SELF,link->nlocal,start,1,&(*is)[cnt]);CHKERRQ(ierr); ierr = DMGetBlockSize(link->dm,&bs);CHKERRQ(ierr); ierr = ISSetBlockSize((*is)[cnt],bs);CHKERRQ(ierr); } PetscFunctionReturn(0); } /*@C DMCompositeGetGlobalISs - Gets the index sets for each composed object Collective on dm Input Parameter: . dm - the packer object Output Parameters: . is - the array of index sets Level: advanced Notes: The is entries should be destroyed with ISDestroy(), the is array should be freed with PetscFree() These could be used to extract a subset of vector entries for a "multi-physics" preconditioner Use DMCompositeGetLocalISs() for index sets in the packed local numbering, and DMCompositeGetISLocalToGlobalMappings() for to map local sub-DM (including ghost) indices to packed global indices. Fortran Notes: The output argument 'is' must be an allocated array of sufficient length, which can be learned using DMCompositeGetNumberDM(). .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(), DMCompositeGather(), DMCompositeCreate(), DMCompositeGetAccess(), DMCompositeScatter(), DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(),DMCompositeGetEntries() @*/ PetscErrorCode DMCompositeGetGlobalISs(DM dm,IS *is[]) { PetscErrorCode ierr; PetscInt cnt = 0; struct DMCompositeLink *next; PetscMPIInt rank; DM_Composite *com = (DM_Composite*)dm->data; PetscBool flg; PetscFunctionBegin; PetscValidHeaderSpecific(dm,DM_CLASSID,1); ierr = PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);CHKERRQ(ierr); if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name); if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Must call DMSetUp() before"); ierr = PetscMalloc1(com->nDM,is);CHKERRQ(ierr); next = com->next; ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);CHKERRQ(ierr); /* loop over packed objects, handling one at at time */ while (next) { PetscDS prob; ierr = ISCreateStride(PetscObjectComm((PetscObject)dm),next->n,next->grstart,1,&(*is)[cnt]);CHKERRQ(ierr); ierr = DMGetDS(dm, &prob);CHKERRQ(ierr); if (prob) { MatNullSpace space; Mat pmat; PetscObject disc; PetscInt Nf; ierr = PetscDSGetNumFields(prob, &Nf);CHKERRQ(ierr); if (cnt < Nf) { ierr = PetscDSGetDiscretization(prob, cnt, &disc);CHKERRQ(ierr); ierr = PetscObjectQuery(disc, "nullspace", (PetscObject*) &space);CHKERRQ(ierr); if (space) {ierr = PetscObjectCompose((PetscObject) (*is)[cnt], "nullspace", (PetscObject) space);CHKERRQ(ierr);} ierr = PetscObjectQuery(disc, "nearnullspace", (PetscObject*) &space);CHKERRQ(ierr); if (space) {ierr = PetscObjectCompose((PetscObject) (*is)[cnt], "nearnullspace", (PetscObject) space);CHKERRQ(ierr);} ierr = PetscObjectQuery(disc, "pmat", (PetscObject*) &pmat);CHKERRQ(ierr); if (pmat) {ierr = PetscObjectCompose((PetscObject) (*is)[cnt], "pmat", (PetscObject) pmat);CHKERRQ(ierr);} } } cnt++; next = next->next; } PetscFunctionReturn(0); } PetscErrorCode DMCreateFieldIS_Composite(DM dm, PetscInt *numFields,char ***fieldNames, IS **fields) { PetscInt nDM; DM *dms; PetscInt i; PetscErrorCode ierr; PetscFunctionBegin; ierr = DMCompositeGetNumberDM(dm, &nDM);CHKERRQ(ierr); if (numFields) *numFields = nDM; ierr = DMCompositeGetGlobalISs(dm, fields);CHKERRQ(ierr); if (fieldNames) { ierr = PetscMalloc1(nDM, &dms);CHKERRQ(ierr); ierr = PetscMalloc1(nDM, fieldNames);CHKERRQ(ierr); ierr = DMCompositeGetEntriesArray(dm, dms);CHKERRQ(ierr); for (i=0; iname; if (!splitname) { ierr = PetscObjectGetOptionsPrefix((PetscObject)dms[i],&splitname);CHKERRQ(ierr); if (splitname) { size_t len; ierr = PetscStrncpy(buf,splitname,sizeof(buf));CHKERRQ(ierr); buf[sizeof(buf) - 1] = 0; ierr = PetscStrlen(buf,&len);CHKERRQ(ierr); if (buf[len-1] == '_') buf[len-1] = 0; /* Remove trailing underscore if it was used */ splitname = buf; } } if (!splitname) { ierr = PetscSNPrintf(buf,sizeof(buf),"%D",i);CHKERRQ(ierr); splitname = buf; } ierr = PetscStrallocpy(splitname,&(*fieldNames)[i]);CHKERRQ(ierr); } ierr = PetscFree(dms);CHKERRQ(ierr); } PetscFunctionReturn(0); } /* This could take over from DMCreateFieldIS(), as it is more general, making DMCreateFieldIS() a special case -- calling with dmlist == NULL; At this point it's probably best to be less intrusive, however. */ PetscErrorCode DMCreateFieldDecomposition_Composite(DM dm, PetscInt *len,char ***namelist, IS **islist, DM **dmlist) { PetscInt nDM; PetscInt i; PetscErrorCode ierr; PetscFunctionBegin; ierr = DMCreateFieldIS_Composite(dm, len, namelist, islist);CHKERRQ(ierr); if (dmlist) { ierr = DMCompositeGetNumberDM(dm, &nDM);CHKERRQ(ierr); ierr = PetscMalloc1(nDM, dmlist);CHKERRQ(ierr); ierr = DMCompositeGetEntriesArray(dm, *dmlist);CHKERRQ(ierr); for (i=0; idata; PetscBool flg; PetscFunctionBegin; PetscValidHeaderSpecific(dm,DM_CLASSID,1); ierr = PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);CHKERRQ(ierr); if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name); next = com->next; /* loop over packed objects, handling one at at time */ va_start(Argp,dm); while (next) { Vec *vec; vec = va_arg(Argp, Vec*); if (vec) {ierr = DMGetLocalVector(next->dm,vec);CHKERRQ(ierr);} next = next->next; } va_end(Argp); PetscFunctionReturn(0); } /*@C DMCompositeRestoreLocalVectors - Restores local vectors for each part of a DMComposite. Not Collective Input Parameter: . dm - the packer object Output Parameter: . Vec ... - the individual sequential Vecs Level: advanced Not available from Fortran .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(), DMCompositeGather(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(), DMCompositeGetLocalVectors(), DMCompositeScatter(), DMCompositeGetEntries() @*/ PetscErrorCode DMCompositeRestoreLocalVectors(DM dm,...) { va_list Argp; PetscErrorCode ierr; struct DMCompositeLink *next; DM_Composite *com = (DM_Composite*)dm->data; PetscBool flg; PetscFunctionBegin; PetscValidHeaderSpecific(dm,DM_CLASSID,1); ierr = PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);CHKERRQ(ierr); if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name); next = com->next; /* loop over packed objects, handling one at at time */ va_start(Argp,dm); while (next) { Vec *vec; vec = va_arg(Argp, Vec*); if (vec) {ierr = DMRestoreLocalVector(next->dm,vec);CHKERRQ(ierr);} next = next->next; } va_end(Argp); PetscFunctionReturn(0); } /* -------------------------------------------------------------------------------------*/ /*@C DMCompositeGetEntries - Gets the DM for each entry in a DMComposite. Not Collective Input Parameter: . dm - the packer object Output Parameter: . DM ... - the individual entries (DMs) Level: advanced Fortran Notes: Available as DMCompositeGetEntries() for one output DM, DMCompositeGetEntries2() for 2, etc .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(), DMCompositeGetEntriesArray() DMCompositeGather(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(), DMCompositeRestoreLocalVectors(), DMCompositeGetLocalVectors(), DMCompositeScatter(), DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors() @*/ PetscErrorCode DMCompositeGetEntries(DM dm,...) { va_list Argp; struct DMCompositeLink *next; DM_Composite *com = (DM_Composite*)dm->data; PetscBool flg; PetscErrorCode ierr; PetscFunctionBegin; PetscValidHeaderSpecific(dm,DM_CLASSID,1); ierr = PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);CHKERRQ(ierr); if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name); next = com->next; /* loop over packed objects, handling one at at time */ va_start(Argp,dm); while (next) { DM *dmn; dmn = va_arg(Argp, DM*); if (dmn) *dmn = next->dm; next = next->next; } va_end(Argp); PetscFunctionReturn(0); } /*@C DMCompositeGetEntriesArray - Gets the DM for each entry in a DMComposite. Not Collective Input Parameter: . dm - the packer object Output Parameter: . dms - array of sufficient length (see DMCompositeGetNumberDM()) to hold the individual DMs Level: advanced .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(), DMCompositeGetEntries() DMCompositeGather(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(), DMCompositeRestoreLocalVectors(), DMCompositeGetLocalVectors(), DMCompositeScatter(), DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors() @*/ PetscErrorCode DMCompositeGetEntriesArray(DM dm,DM dms[]) { struct DMCompositeLink *next; DM_Composite *com = (DM_Composite*)dm->data; PetscInt i; PetscBool flg; PetscErrorCode ierr; PetscFunctionBegin; PetscValidHeaderSpecific(dm,DM_CLASSID,1); ierr = PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);CHKERRQ(ierr); if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name); /* loop over packed objects, handling one at at time */ for (next=com->next,i=0; next; next=next->next,i++) dms[i] = next->dm; PetscFunctionReturn(0); } typedef struct { DM dm; PetscViewer *subv; Vec *vecs; } GLVisViewerCtx; static PetscErrorCode DestroyGLVisViewerCtx_Private(void *vctx) { GLVisViewerCtx *ctx = (GLVisViewerCtx*)vctx; PetscInt i,n; PetscErrorCode ierr; PetscFunctionBegin; ierr = DMCompositeGetNumberDM(ctx->dm,&n);CHKERRQ(ierr); for (i = 0; i < n; i++) { ierr = PetscViewerDestroy(&ctx->subv[i]);CHKERRQ(ierr); } ierr = PetscFree2(ctx->subv,ctx->vecs);CHKERRQ(ierr); ierr = DMDestroy(&ctx->dm);CHKERRQ(ierr); ierr = PetscFree(ctx);CHKERRQ(ierr); PetscFunctionReturn(0); } static PetscErrorCode DMCompositeSampleGLVisFields_Private(PetscObject oX, PetscInt nf, PetscObject oXfield[], void *vctx) { Vec X = (Vec)oX; GLVisViewerCtx *ctx = (GLVisViewerCtx*)vctx; PetscInt i,n,cumf; PetscErrorCode ierr; PetscFunctionBegin; ierr = DMCompositeGetNumberDM(ctx->dm,&n);CHKERRQ(ierr); ierr = DMCompositeGetAccessArray(ctx->dm,X,n,NULL,ctx->vecs);CHKERRQ(ierr); for (i = 0, cumf = 0; i < n; i++) { PetscErrorCode (*g2l)(PetscObject,PetscInt,PetscObject[],void*); void *fctx; PetscInt nfi; ierr = PetscViewerGLVisGetFields_Private(ctx->subv[i],&nfi,NULL,NULL,&g2l,NULL,&fctx);CHKERRQ(ierr); if (!nfi) continue; if (g2l) { ierr = (*g2l)((PetscObject)ctx->vecs[i],nfi,oXfield+cumf,fctx);CHKERRQ(ierr); } else { ierr = VecCopy(ctx->vecs[i],(Vec)(oXfield[cumf]));CHKERRQ(ierr); } cumf += nfi; } ierr = DMCompositeRestoreAccessArray(ctx->dm,X,n,NULL,ctx->vecs);CHKERRQ(ierr); PetscFunctionReturn(0); } static PetscErrorCode DMSetUpGLVisViewer_Composite(PetscObject odm, PetscViewer viewer) { DM dm = (DM)odm, *dms; Vec *Ufds; GLVisViewerCtx *ctx; PetscInt i,n,tnf,*sdim; char **fecs; PetscErrorCode ierr; PetscFunctionBegin; ierr = PetscNew(&ctx);CHKERRQ(ierr); ierr = PetscObjectReference((PetscObject)dm);CHKERRQ(ierr); ctx->dm = dm; ierr = DMCompositeGetNumberDM(dm,&n);CHKERRQ(ierr); ierr = PetscMalloc1(n,&dms);CHKERRQ(ierr); ierr = DMCompositeGetEntriesArray(dm,dms);CHKERRQ(ierr); ierr = PetscMalloc2(n,&ctx->subv,n,&ctx->vecs);CHKERRQ(ierr); for (i = 0, tnf = 0; i < n; i++) { PetscInt nf; ierr = PetscViewerCreate(PetscObjectComm(odm),&ctx->subv[i]);CHKERRQ(ierr); ierr = PetscViewerSetType(ctx->subv[i],PETSCVIEWERGLVIS);CHKERRQ(ierr); ierr = PetscViewerGLVisSetDM_Private(ctx->subv[i],(PetscObject)dms[i]);CHKERRQ(ierr); ierr = PetscViewerGLVisGetFields_Private(ctx->subv[i],&nf,NULL,NULL,NULL,NULL,NULL);CHKERRQ(ierr); tnf += nf; } ierr = PetscFree(dms);CHKERRQ(ierr); ierr = PetscMalloc3(tnf,&fecs,tnf,&sdim,tnf,&Ufds);CHKERRQ(ierr); for (i = 0, tnf = 0; i < n; i++) { PetscInt *sd,nf,f; const char **fec; Vec *Uf; ierr = PetscViewerGLVisGetFields_Private(ctx->subv[i],&nf,&fec,&sd,NULL,(PetscObject**)&Uf,NULL);CHKERRQ(ierr); for (f = 0; f < nf; f++) { ierr = PetscStrallocpy(fec[f],&fecs[tnf+f]);CHKERRQ(ierr); Ufds[tnf+f] = Uf[f]; sdim[tnf+f] = sd[f]; } tnf += nf; } ierr = PetscViewerGLVisSetFields(viewer,tnf,(const char**)fecs,sdim,DMCompositeSampleGLVisFields_Private,(PetscObject*)Ufds,ctx,DestroyGLVisViewerCtx_Private);CHKERRQ(ierr); for (i = 0; i < tnf; i++) { ierr = PetscFree(fecs[i]);CHKERRQ(ierr); } ierr = PetscFree3(fecs,sdim,Ufds);CHKERRQ(ierr); PetscFunctionReturn(0); } PetscErrorCode DMRefine_Composite(DM dmi,MPI_Comm comm,DM *fine) { PetscErrorCode ierr; struct DMCompositeLink *next; DM_Composite *com = (DM_Composite*)dmi->data; DM dm; PetscFunctionBegin; PetscValidHeaderSpecific(dmi,DM_CLASSID,1); if (comm == MPI_COMM_NULL) { ierr = PetscObjectGetComm((PetscObject)dmi,&comm);CHKERRQ(ierr); } ierr = DMSetUp(dmi);CHKERRQ(ierr); next = com->next; ierr = DMCompositeCreate(comm,fine);CHKERRQ(ierr); /* loop over packed objects, handling one at at time */ while (next) { ierr = DMRefine(next->dm,comm,&dm);CHKERRQ(ierr); ierr = DMCompositeAddDM(*fine,dm);CHKERRQ(ierr); ierr = PetscObjectDereference((PetscObject)dm);CHKERRQ(ierr); next = next->next; } PetscFunctionReturn(0); } PetscErrorCode DMCoarsen_Composite(DM dmi,MPI_Comm comm,DM *fine) { PetscErrorCode ierr; struct DMCompositeLink *next; DM_Composite *com = (DM_Composite*)dmi->data; DM dm; PetscFunctionBegin; PetscValidHeaderSpecific(dmi,DM_CLASSID,1); ierr = DMSetUp(dmi);CHKERRQ(ierr); if (comm == MPI_COMM_NULL) { ierr = PetscObjectGetComm((PetscObject)dmi,&comm);CHKERRQ(ierr); } next = com->next; ierr = DMCompositeCreate(comm,fine);CHKERRQ(ierr); /* loop over packed objects, handling one at at time */ while (next) { ierr = DMCoarsen(next->dm,comm,&dm);CHKERRQ(ierr); ierr = DMCompositeAddDM(*fine,dm);CHKERRQ(ierr); ierr = PetscObjectDereference((PetscObject)dm);CHKERRQ(ierr); next = next->next; } PetscFunctionReturn(0); } PetscErrorCode DMCreateInterpolation_Composite(DM coarse,DM fine,Mat *A,Vec *v) { PetscErrorCode ierr; PetscInt m,n,M,N,nDM,i; struct DMCompositeLink *nextc; struct DMCompositeLink *nextf; Vec gcoarse,gfine,*vecs; DM_Composite *comcoarse = (DM_Composite*)coarse->data; DM_Composite *comfine = (DM_Composite*)fine->data; Mat *mats; PetscFunctionBegin; PetscValidHeaderSpecific(coarse,DM_CLASSID,1); PetscValidHeaderSpecific(fine,DM_CLASSID,2); ierr = DMSetUp(coarse);CHKERRQ(ierr); ierr = DMSetUp(fine);CHKERRQ(ierr); /* use global vectors only for determining matrix layout */ ierr = DMGetGlobalVector(coarse,&gcoarse);CHKERRQ(ierr); ierr = DMGetGlobalVector(fine,&gfine);CHKERRQ(ierr); ierr = VecGetLocalSize(gcoarse,&n);CHKERRQ(ierr); ierr = VecGetLocalSize(gfine,&m);CHKERRQ(ierr); ierr = VecGetSize(gcoarse,&N);CHKERRQ(ierr); ierr = VecGetSize(gfine,&M);CHKERRQ(ierr); ierr = DMRestoreGlobalVector(coarse,&gcoarse);CHKERRQ(ierr); ierr = DMRestoreGlobalVector(fine,&gfine);CHKERRQ(ierr); nDM = comfine->nDM; if (nDM != comcoarse->nDM) SETERRQ2(PetscObjectComm((PetscObject)fine),PETSC_ERR_ARG_INCOMP,"Fine DMComposite has %D entries, but coarse has %D",nDM,comcoarse->nDM); ierr = PetscCalloc1(nDM*nDM,&mats);CHKERRQ(ierr); if (v) { ierr = PetscCalloc1(nDM,&vecs);CHKERRQ(ierr); } /* loop over packed objects, handling one at at time */ for (nextc=comcoarse->next,nextf=comfine->next,i=0; nextc; nextc=nextc->next,nextf=nextf->next,i++) { if (!v) { ierr = DMCreateInterpolation(nextc->dm,nextf->dm,&mats[i*nDM+i],NULL);CHKERRQ(ierr); } else { ierr = DMCreateInterpolation(nextc->dm,nextf->dm,&mats[i*nDM+i],&vecs[i]);CHKERRQ(ierr); } } ierr = MatCreateNest(PetscObjectComm((PetscObject)fine),nDM,NULL,nDM,NULL,mats,A);CHKERRQ(ierr); if (v) { ierr = VecCreateNest(PetscObjectComm((PetscObject)fine),nDM,NULL,vecs,v);CHKERRQ(ierr); } for (i=0; idata; ISLocalToGlobalMapping *ltogs; PetscInt i; PetscErrorCode ierr; PetscFunctionBegin; /* Set the ISLocalToGlobalMapping on the new matrix */ ierr = DMCompositeGetISLocalToGlobalMappings(dm,<ogs);CHKERRQ(ierr); ierr = ISLocalToGlobalMappingConcatenate(PetscObjectComm((PetscObject)dm),com->nDM,ltogs,&dm->ltogmap);CHKERRQ(ierr); for (i=0; inDM; i++) {ierr = ISLocalToGlobalMappingDestroy(<ogs[i]);CHKERRQ(ierr);} ierr = PetscFree(ltogs);CHKERRQ(ierr); PetscFunctionReturn(0); } PetscErrorCode DMCreateColoring_Composite(DM dm,ISColoringType ctype,ISColoring *coloring) { PetscErrorCode ierr; PetscInt n,i,cnt; ISColoringValue *colors; PetscBool dense = PETSC_FALSE; ISColoringValue maxcol = 0; DM_Composite *com = (DM_Composite*)dm->data; PetscFunctionBegin; PetscValidHeaderSpecific(dm,DM_CLASSID,1); if (ctype == IS_COLORING_LOCAL) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Only global coloring supported"); else if (ctype == IS_COLORING_GLOBAL) { n = com->n; } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unknown ISColoringType"); ierr = PetscMalloc1(n,&colors);CHKERRQ(ierr); /* freed in ISColoringDestroy() */ ierr = PetscOptionsGetBool(((PetscObject)dm)->options,((PetscObject)dm)->prefix,"-dmcomposite_dense_jacobian",&dense,NULL);CHKERRQ(ierr); if (dense) { for (i=0; irstart + i); } maxcol = com->N; } else { struct DMCompositeLink *next = com->next; PetscMPIInt rank; ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);CHKERRQ(ierr); cnt = 0; while (next) { ISColoring lcoloring; ierr = DMCreateColoring(next->dm,IS_COLORING_GLOBAL,&lcoloring);CHKERRQ(ierr); for (i=0; iN; i++) { colors[cnt++] = maxcol + lcoloring->colors[i]; } maxcol += lcoloring->n; ierr = ISColoringDestroy(&lcoloring);CHKERRQ(ierr); next = next->next; } } ierr = ISColoringCreate(PetscObjectComm((PetscObject)dm),maxcol,n,colors,PETSC_OWN_POINTER,coloring);CHKERRQ(ierr); PetscFunctionReturn(0); } PetscErrorCode DMGlobalToLocalBegin_Composite(DM dm,Vec gvec,InsertMode mode,Vec lvec) { PetscErrorCode ierr; struct DMCompositeLink *next; PetscScalar *garray,*larray; DM_Composite *com = (DM_Composite*)dm->data; PetscFunctionBegin; PetscValidHeaderSpecific(dm,DM_CLASSID,1); PetscValidHeaderSpecific(gvec,VEC_CLASSID,2); if (!com->setup) { ierr = DMSetUp(dm);CHKERRQ(ierr); } ierr = VecGetArray(gvec,&garray);CHKERRQ(ierr); ierr = VecGetArray(lvec,&larray);CHKERRQ(ierr); /* loop over packed objects, handling one at at time */ next = com->next; while (next) { Vec local,global; PetscInt N; ierr = DMGetGlobalVector(next->dm,&global);CHKERRQ(ierr); ierr = VecGetLocalSize(global,&N);CHKERRQ(ierr); ierr = VecPlaceArray(global,garray);CHKERRQ(ierr); ierr = DMGetLocalVector(next->dm,&local);CHKERRQ(ierr); ierr = VecPlaceArray(local,larray);CHKERRQ(ierr); ierr = DMGlobalToLocalBegin(next->dm,global,mode,local);CHKERRQ(ierr); ierr = DMGlobalToLocalEnd(next->dm,global,mode,local);CHKERRQ(ierr); ierr = VecResetArray(global);CHKERRQ(ierr); ierr = VecResetArray(local);CHKERRQ(ierr); ierr = DMRestoreGlobalVector(next->dm,&global);CHKERRQ(ierr); ierr = DMRestoreLocalVector(next->dm,&local);CHKERRQ(ierr); larray += next->nlocal; garray += next->n; next = next->next; } ierr = VecRestoreArray(gvec,NULL);CHKERRQ(ierr); ierr = VecRestoreArray(lvec,NULL);CHKERRQ(ierr); PetscFunctionReturn(0); } PetscErrorCode DMGlobalToLocalEnd_Composite(DM dm,Vec gvec,InsertMode mode,Vec lvec) { PetscFunctionBegin; PetscValidHeaderSpecific(dm,DM_CLASSID,1); PetscValidHeaderSpecific(gvec,VEC_CLASSID,2); PetscValidHeaderSpecific(lvec,VEC_CLASSID,4); PetscFunctionReturn(0); } PetscErrorCode DMLocalToGlobalBegin_Composite(DM dm,Vec lvec,InsertMode mode,Vec gvec) { PetscErrorCode ierr; struct DMCompositeLink *next; PetscScalar *larray,*garray; DM_Composite *com = (DM_Composite*)dm->data; PetscFunctionBegin; PetscValidHeaderSpecific(dm,DM_CLASSID,1); PetscValidHeaderSpecific(lvec,VEC_CLASSID,2); PetscValidHeaderSpecific(gvec,VEC_CLASSID,4); if (!com->setup) { ierr = DMSetUp(dm);CHKERRQ(ierr); } ierr = VecGetArray(lvec,&larray);CHKERRQ(ierr); ierr = VecGetArray(gvec,&garray);CHKERRQ(ierr); /* loop over packed objects, handling one at at time */ next = com->next; while (next) { Vec global,local; ierr = DMGetLocalVector(next->dm,&local);CHKERRQ(ierr); ierr = VecPlaceArray(local,larray);CHKERRQ(ierr); ierr = DMGetGlobalVector(next->dm,&global);CHKERRQ(ierr); ierr = VecPlaceArray(global,garray);CHKERRQ(ierr); ierr = DMLocalToGlobalBegin(next->dm,local,mode,global);CHKERRQ(ierr); ierr = DMLocalToGlobalEnd(next->dm,local,mode,global);CHKERRQ(ierr); ierr = VecResetArray(local);CHKERRQ(ierr); ierr = VecResetArray(global);CHKERRQ(ierr); ierr = DMRestoreGlobalVector(next->dm,&global);CHKERRQ(ierr); ierr = DMRestoreLocalVector(next->dm,&local);CHKERRQ(ierr); garray += next->n; larray += next->nlocal; next = next->next; } ierr = VecRestoreArray(gvec,NULL);CHKERRQ(ierr); ierr = VecRestoreArray(lvec,NULL);CHKERRQ(ierr); PetscFunctionReturn(0); } PetscErrorCode DMLocalToGlobalEnd_Composite(DM dm,Vec lvec,InsertMode mode,Vec gvec) { PetscFunctionBegin; PetscValidHeaderSpecific(dm,DM_CLASSID,1); PetscValidHeaderSpecific(lvec,VEC_CLASSID,2); PetscValidHeaderSpecific(gvec,VEC_CLASSID,4); PetscFunctionReturn(0); } PetscErrorCode DMLocalToLocalBegin_Composite(DM dm,Vec vec1,InsertMode mode,Vec vec2) { PetscErrorCode ierr; struct DMCompositeLink *next; PetscScalar *array1,*array2; DM_Composite *com = (DM_Composite*)dm->data; PetscFunctionBegin; PetscValidHeaderSpecific(dm,DM_CLASSID,1); PetscValidHeaderSpecific(vec1,VEC_CLASSID,2); PetscValidHeaderSpecific(vec2,VEC_CLASSID,4); if (!com->setup) { ierr = DMSetUp(dm);CHKERRQ(ierr); } ierr = VecGetArray(vec1,&array1);CHKERRQ(ierr); ierr = VecGetArray(vec2,&array2);CHKERRQ(ierr); /* loop over packed objects, handling one at at time */ next = com->next; while (next) { Vec local1,local2; ierr = DMGetLocalVector(next->dm,&local1);CHKERRQ(ierr); ierr = VecPlaceArray(local1,array1);CHKERRQ(ierr); ierr = DMGetLocalVector(next->dm,&local2);CHKERRQ(ierr); ierr = VecPlaceArray(local2,array2);CHKERRQ(ierr); ierr = DMLocalToLocalBegin(next->dm,local1,mode,local2);CHKERRQ(ierr); ierr = DMLocalToLocalEnd(next->dm,local1,mode,local2);CHKERRQ(ierr); ierr = VecResetArray(local2);CHKERRQ(ierr); ierr = DMRestoreLocalVector(next->dm,&local2);CHKERRQ(ierr); ierr = VecResetArray(local1);CHKERRQ(ierr); ierr = DMRestoreLocalVector(next->dm,&local1);CHKERRQ(ierr); array1 += next->nlocal; array2 += next->nlocal; next = next->next; } ierr = VecRestoreArray(vec1,NULL);CHKERRQ(ierr); ierr = VecRestoreArray(vec2,NULL);CHKERRQ(ierr); PetscFunctionReturn(0); } PetscErrorCode DMLocalToLocalEnd_Composite(DM dm,Vec lvec,InsertMode mode,Vec gvec) { PetscFunctionBegin; PetscValidHeaderSpecific(dm,DM_CLASSID,1); PetscValidHeaderSpecific(lvec,VEC_CLASSID,2); PetscValidHeaderSpecific(gvec,VEC_CLASSID,4); PetscFunctionReturn(0); } /*MC DMCOMPOSITE = "composite" - A DM object that is used to manage data for a collection of DMs Level: intermediate .seealso: DMType, DM, DMDACreate(), DMCreate(), DMSetType(), DMCompositeCreate() M*/ PETSC_EXTERN PetscErrorCode DMCreate_Composite(DM p) { PetscErrorCode ierr; DM_Composite *com; PetscFunctionBegin; ierr = PetscNewLog(p,&com);CHKERRQ(ierr); p->data = com; com->n = 0; com->nghost = 0; com->next = NULL; com->nDM = 0; p->ops->createglobalvector = DMCreateGlobalVector_Composite; p->ops->createlocalvector = DMCreateLocalVector_Composite; p->ops->getlocaltoglobalmapping = DMGetLocalToGlobalMapping_Composite; p->ops->createfieldis = DMCreateFieldIS_Composite; p->ops->createfielddecomposition = DMCreateFieldDecomposition_Composite; p->ops->refine = DMRefine_Composite; p->ops->coarsen = DMCoarsen_Composite; p->ops->createinterpolation = DMCreateInterpolation_Composite; p->ops->creatematrix = DMCreateMatrix_Composite; p->ops->getcoloring = DMCreateColoring_Composite; p->ops->globaltolocalbegin = DMGlobalToLocalBegin_Composite; p->ops->globaltolocalend = DMGlobalToLocalEnd_Composite; p->ops->localtoglobalbegin = DMLocalToGlobalBegin_Composite; p->ops->localtoglobalend = DMLocalToGlobalEnd_Composite; p->ops->localtolocalbegin = DMLocalToLocalBegin_Composite; p->ops->localtolocalend = DMLocalToLocalEnd_Composite; p->ops->destroy = DMDestroy_Composite; p->ops->view = DMView_Composite; p->ops->setup = DMSetUp_Composite; ierr = PetscObjectComposeFunction((PetscObject)p,"DMSetUpGLVisViewer_C",DMSetUpGLVisViewer_Composite);CHKERRQ(ierr); PetscFunctionReturn(0); } /*@ DMCompositeCreate - Creates a vector packer, used to generate "composite" vectors made up of several subvectors. Collective Input Parameter: . comm - the processors that will share the global vector Output Parameters: . packer - the packer object Level: advanced .seealso DMDestroy(), DMCompositeAddDM(), DMCompositeScatter(), DMCOMPOSITE,DMCreate() DMCompositeGather(), DMCreateGlobalVector(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess() DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(), DMCompositeGetEntries() @*/ PetscErrorCode DMCompositeCreate(MPI_Comm comm,DM *packer) { PetscErrorCode ierr; PetscFunctionBegin; PetscValidPointer(packer,2); ierr = DMCreate(comm,packer);CHKERRQ(ierr); ierr = DMSetType(*packer,DMCOMPOSITE);CHKERRQ(ierr); PetscFunctionReturn(0); }