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