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 DMDA 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 .keywords: distributed array, get, global, indices, local-to-global 69 70 .seealso: DMDACreate2d(), DMDAGetAO(), DMDAGetGhostCorners(), DMDAGetCorners(), DMDALocalToGlocal() 71 DMGlobalToLocalBegin(), DMGlobalToLocalEnd(), DMLocalToLocalBegin(), DMLocalToLocalEnd(), DMDAGetGlobalIndices(), DMDAGetOwnershipRanges(), 72 AO, AOPetscToApplication(), AOApplicationToPetsc() 73 @*/ 74 PetscErrorCode DMDASetAOType(DM da,AOType aotype) 75 { 76 DM_DA *dd; 77 PetscBool isdmda; 78 PetscErrorCode ierr; 79 80 PetscFunctionBegin; 81 PetscValidHeaderSpecificType(da,DM_CLASSID,1,DMDA); 82 ierr = PetscObjectTypeCompare((PetscObject)da,DMDA,&isdmda);CHKERRQ(ierr); 83 if (!isdmda) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Requires a DMDA as input"); 84 /* now we can safely dereference */ 85 dd = (DM_DA*)da->data; 86 if (dd->ao) { /* check if the already computed AO has the same type as requested */ 87 PetscBool match; 88 ierr = PetscObjectTypeCompare((PetscObject)dd->ao,aotype,&match);CHKERRQ(ierr); 89 if (!match) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Cannot change AO type"); 90 PetscFunctionReturn(0); 91 } 92 ierr = PetscFree(dd->aotype);CHKERRQ(ierr); 93 ierr = PetscStrallocpy(aotype,(char**)&dd->aotype);CHKERRQ(ierr); 94 PetscFunctionReturn(0); 95 } 96 97 /*@ 98 DMDAGetAO - Gets the application ordering context for a distributed array. 99 100 Collective on DMDA 101 102 Input Parameter: 103 . da - the distributed array 104 105 Output Parameters: 106 . ao - the application ordering context for DMDAs 107 108 Level: intermediate 109 110 Notes: 111 In this case, the AO maps to the natural grid ordering that would be used 112 for the DMDA if only 1 processor were employed (ordering most rapidly in the 113 x-direction, then y, then z). Multiple degrees of freedom are numbered 114 for each node (rather than 1 component for the whole grid, then the next 115 component, etc.) 116 117 Do NOT call AODestroy() on the ao returned by this function. 118 119 .keywords: distributed array, get, global, indices, local-to-global 120 121 .seealso: DMDACreate2d(), DMDASetAOType(), DMDAGetGhostCorners(), DMDAGetCorners(), DMDALocalToGlocal() 122 DMGlobalToLocalBegin(), DMGlobalToLocalEnd(), DMLocalToLocalBegin(), DMLocalToLocalEnd(), DMDAGetOwnershipRanges(), 123 AO, AOPetscToApplication(), AOApplicationToPetsc() 124 @*/ 125 PetscErrorCode DMDAGetAO(DM da,AO *ao) 126 { 127 DM_DA *dd; 128 PetscBool isdmda; 129 PetscErrorCode ierr; 130 131 PetscFunctionBegin; 132 PetscValidHeaderSpecificType(da,DM_CLASSID,1,DMDA); 133 PetscValidPointer(ao,2); 134 ierr = PetscObjectTypeCompare((PetscObject)da,DMDA,&isdmda);CHKERRQ(ierr); 135 if (!isdmda) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Requires a DMDA as input"); 136 /* now we can safely dereference */ 137 dd = (DM_DA*)da->data; 138 139 /* 140 Build the natural ordering to PETSc ordering mappings. 141 */ 142 if (!dd->ao) { 143 IS ispetsc,isnatural; 144 PetscErrorCode ierr; 145 PetscInt Nlocal; 146 147 ierr = DMDAGetNatural_Private(da,&Nlocal,&isnatural);CHKERRQ(ierr); 148 ierr = ISCreateStride(PetscObjectComm((PetscObject)da),Nlocal,dd->base,1,&ispetsc);CHKERRQ(ierr); 149 ierr = AOCreate(PetscObjectComm((PetscObject)da),&dd->ao);CHKERRQ(ierr); 150 ierr = AOSetIS(dd->ao,isnatural,ispetsc);CHKERRQ(ierr); 151 ierr = AOSetType(dd->ao,dd->aotype);CHKERRQ(ierr); 152 ierr = PetscLogObjectParent((PetscObject)da,(PetscObject)dd->ao);CHKERRQ(ierr); 153 ierr = ISDestroy(&ispetsc);CHKERRQ(ierr); 154 ierr = ISDestroy(&isnatural);CHKERRQ(ierr); 155 } 156 *ao = dd->ao; 157 PetscFunctionReturn(0); 158 } 159 160 161