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