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 */
DMDAGetNatural_Private(DM da,PetscInt * outNlocal,IS * isnatural)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 @*/
DMDASetAOType(DM da,AOType aotype)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 @*/
DMDAGetAO(DM da,AO * ao)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