1 2 /* 3 Code for manipulating distributed regular arrays in parallel. 4 */ 5 6 #include <petsc/private/dmdaimpl.h> /*I "petscdmda.h" I*/ 7 8 /* 9 Gets the natural number for each global number on the process. 10 11 Used by DMDAGetAO() and DMDAGlobalToNatural_Create() 12 */ 13 PetscErrorCode DMDAGetNatural_Private(DM da,PetscInt *outNlocal,IS *isnatural) 14 { 15 PetscErrorCode ierr; 16 PetscInt Nlocal,i,j,k,*lidx,lict = 0,dim = da->dim; 17 DM_DA *dd = (DM_DA*)da->data; 18 19 PetscFunctionBegin; 20 Nlocal = (dd->xe-dd->xs); 21 if (dim > 1) Nlocal *= (dd->ye-dd->ys); 22 if (dim > 2) Nlocal *= (dd->ze-dd->zs); 23 24 ierr = PetscMalloc1(Nlocal,&lidx);CHKERRQ(ierr); 25 26 if (dim == 1) { 27 for (i=dd->xs; i<dd->xe; i++) { 28 /* global number in natural ordering */ 29 lidx[lict++] = i; 30 } 31 } else if (dim == 2) { 32 for (j=dd->ys; j<dd->ye; j++) { 33 for (i=dd->xs; i<dd->xe; i++) { 34 /* global number in natural ordering */ 35 lidx[lict++] = i + j*dd->M*dd->w; 36 } 37 } 38 } else if (dim == 3) { 39 for (k=dd->zs; k<dd->ze; k++) { 40 for (j=dd->ys; j<dd->ye; j++) { 41 for (i=dd->xs; i<dd->xe; i++) { 42 lidx[lict++] = i + j*dd->M*dd->w + k*dd->M*dd->N*dd->w; 43 } 44 } 45 } 46 } 47 *outNlocal = Nlocal; 48 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)da),Nlocal,lidx,PETSC_OWN_POINTER,isnatural);CHKERRQ(ierr); 49 PetscFunctionReturn(0); 50 } 51 52 /*@C 53 DMDASetAOType - Sets the type of application ordering for a distributed array. 54 55 Collective on da 56 57 Input Parameter: 58 + da - the distributed array 59 - aotype - type of AO 60 61 Output Parameters: 62 63 Level: intermediate 64 65 Notes: 66 It will generate and error if an AO has already been obtained with a call to DMDAGetAO and the user sets a different AOType 67 68 .seealso: DMDACreate2d(), DMDAGetAO(), DMDAGetGhostCorners(), DMDAGetCorners(), DMLocalToGlobal() 69 DMGlobalToLocalBegin(), DMGlobalToLocalEnd(), DMLocalToLocalBegin(), DMLocalToLocalEnd(), DMDAGetGlobalIndices(), DMDAGetOwnershipRanges(), 70 AO, AOPetscToApplication(), AOApplicationToPetsc() 71 @*/ 72 PetscErrorCode DMDASetAOType(DM da,AOType aotype) 73 { 74 DM_DA *dd; 75 PetscBool isdmda; 76 PetscErrorCode ierr; 77 78 PetscFunctionBegin; 79 PetscValidHeaderSpecificType(da,DM_CLASSID,1,DMDA); 80 ierr = PetscObjectTypeCompare((PetscObject)da,DMDA,&isdmda);CHKERRQ(ierr); 81 if (!isdmda) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Requires a DMDA as input"); 82 /* now we can safely dereference */ 83 dd = (DM_DA*)da->data; 84 if (dd->ao) { /* check if the already computed AO has the same type as requested */ 85 PetscBool match; 86 ierr = PetscObjectTypeCompare((PetscObject)dd->ao,aotype,&match);CHKERRQ(ierr); 87 if (!match) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Cannot change AO type"); 88 PetscFunctionReturn(0); 89 } 90 ierr = PetscFree(dd->aotype);CHKERRQ(ierr); 91 ierr = PetscStrallocpy(aotype,(char**)&dd->aotype);CHKERRQ(ierr); 92 PetscFunctionReturn(0); 93 } 94 95 /*@ 96 DMDAGetAO - Gets the application ordering context for a distributed array. 97 98 Collective on da 99 100 Input Parameter: 101 . da - the distributed array 102 103 Output Parameters: 104 . ao - the application ordering context for DMDAs 105 106 Level: intermediate 107 108 Notes: 109 In this case, the AO maps to the natural grid ordering that would be used 110 for the DMDA if only 1 processor were employed (ordering most rapidly in the 111 x-direction, then y, then z). Multiple degrees of freedom are numbered 112 for each node (rather than 1 component for the whole grid, then the next 113 component, etc.) 114 115 Do NOT call AODestroy() on the ao returned by this function. 116 117 .seealso: DMDACreate2d(), DMDASetAOType(), DMDAGetGhostCorners(), DMDAGetCorners(), DMLocalToGlobal() 118 DMGlobalToLocalBegin(), DMGlobalToLocalEnd(), DMLocalToLocalBegin(), DMLocalToLocalEnd(), DMDAGetOwnershipRanges(), 119 AO, AOPetscToApplication(), AOApplicationToPetsc() 120 @*/ 121 PetscErrorCode DMDAGetAO(DM da,AO *ao) 122 { 123 DM_DA *dd; 124 PetscBool isdmda; 125 PetscErrorCode ierr; 126 127 PetscFunctionBegin; 128 PetscValidHeaderSpecificType(da,DM_CLASSID,1,DMDA); 129 PetscValidPointer(ao,2); 130 ierr = PetscObjectTypeCompare((PetscObject)da,DMDA,&isdmda);CHKERRQ(ierr); 131 if (!isdmda) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Requires a DMDA as input"); 132 /* now we can safely dereference */ 133 dd = (DM_DA*)da->data; 134 135 /* 136 Build the natural ordering to PETSc ordering mappings. 137 */ 138 if (!dd->ao) { 139 IS ispetsc,isnatural; 140 PetscErrorCode ierr; 141 PetscInt Nlocal; 142 143 ierr = DMDAGetNatural_Private(da,&Nlocal,&isnatural);CHKERRQ(ierr); 144 ierr = ISCreateStride(PetscObjectComm((PetscObject)da),Nlocal,dd->base,1,&ispetsc);CHKERRQ(ierr); 145 ierr = AOCreate(PetscObjectComm((PetscObject)da),&dd->ao);CHKERRQ(ierr); 146 ierr = AOSetIS(dd->ao,isnatural,ispetsc);CHKERRQ(ierr); 147 ierr = AOSetType(dd->ao,dd->aotype);CHKERRQ(ierr); 148 ierr = PetscLogObjectParent((PetscObject)da,(PetscObject)dd->ao);CHKERRQ(ierr); 149 ierr = ISDestroy(&ispetsc);CHKERRQ(ierr); 150 ierr = ISDestroy(&isnatural);CHKERRQ(ierr); 151 } 152 *ao = dd->ao; 153 PetscFunctionReturn(0); 154 } 155 156 157