Actual source code: ex21.c

  1: static const char help[] = "Test DMCreateInjection() for mapping coordinates in 3D";

  3: #include <petscvec.h>
  4: #include <petscmat.h>
  5: #include <petscdm.h>
  6: #include <petscdmda.h>

  8: PetscErrorCode test1_DAInjection3d(PetscInt mx, PetscInt my, PetscInt mz)
  9: {
 10:   DM             dac, daf;
 11:   PetscViewer    vv;
 12:   Vec            ac, af;
 13:   PetscInt       periodicity;
 14:   DMBoundaryType bx, by, bz;

 16:   PetscFunctionBeginUser;
 17:   bx = DM_BOUNDARY_NONE;
 18:   by = DM_BOUNDARY_NONE;
 19:   bz = DM_BOUNDARY_NONE;

 21:   periodicity = 0;

 23:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-periodic", &periodicity, NULL));
 24:   if (periodicity == 1) {
 25:     bx = DM_BOUNDARY_PERIODIC;
 26:   } else if (periodicity == 2) {
 27:     by = DM_BOUNDARY_PERIODIC;
 28:   } else if (periodicity == 3) {
 29:     bz = DM_BOUNDARY_PERIODIC;
 30:   }

 32:   PetscCall(DMDACreate3d(PETSC_COMM_WORLD, bx, by, bz, DMDA_STENCIL_BOX, mx + 1, my + 1, mz + 1, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, 1, /* 1 dof */
 33:                          1, /* stencil = 1 */ NULL, NULL, NULL, &daf));
 34:   PetscCall(DMSetFromOptions(daf));
 35:   PetscCall(DMSetUp(daf));

 37:   PetscCall(DMCoarsen(daf, MPI_COMM_NULL, &dac));

 39:   PetscCall(DMDASetUniformCoordinates(dac, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0));
 40:   PetscCall(DMDASetUniformCoordinates(daf, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0));

 42:   {
 43:     DM         cdaf, cdac;
 44:     Vec        coordsc, coordsf, coordsf2;
 45:     Mat        inject;
 46:     VecScatter vscat;
 47:     Mat        interp;
 48:     PetscReal  norm;

 50:     PetscCall(DMGetCoordinateDM(dac, &cdac));
 51:     PetscCall(DMGetCoordinateDM(daf, &cdaf));

 53:     PetscCall(DMGetCoordinates(dac, &coordsc));
 54:     PetscCall(DMGetCoordinates(daf, &coordsf));

 56:     PetscCall(DMCreateInjection(cdac, cdaf, &inject));
 57:     PetscCall(MatScatterGetVecScatter(inject, &vscat));
 58:     PetscCall(VecScatterBegin(vscat, coordsf, coordsc, INSERT_VALUES, SCATTER_FORWARD));
 59:     PetscCall(VecScatterEnd(vscat, coordsf, coordsc, INSERT_VALUES, SCATTER_FORWARD));
 60:     PetscCall(MatDestroy(&inject));

 62:     PetscCall(DMCreateInterpolation(cdac, cdaf, &interp, NULL));
 63:     PetscCall(VecDuplicate(coordsf, &coordsf2));
 64:     PetscCall(MatInterpolate(interp, coordsc, coordsf2));
 65:     PetscCall(VecAXPY(coordsf2, -1.0, coordsf));
 66:     PetscCall(VecNorm(coordsf2, NORM_MAX, &norm));
 67:     /* The fine coordinates are only reproduced in certain cases */
 68:     if (!bx && !by && !bz && norm > PETSC_SQRT_MACHINE_EPSILON) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Norm %g\n", (double)norm));
 69:     PetscCall(VecDestroy(&coordsf2));
 70:     PetscCall(MatDestroy(&interp));
 71:   }

 73:   if (0) {
 74:     PetscCall(DMCreateGlobalVector(dac, &ac));
 75:     PetscCall(VecZeroEntries(ac));

 77:     PetscCall(DMCreateGlobalVector(daf, &af));
 78:     PetscCall(VecZeroEntries(af));

 80:     PetscCall(PetscViewerASCIIOpen(PETSC_COMM_WORLD, "dac_7.vtu", &vv));
 81:     PetscCall(VecView(ac, vv));
 82:     PetscCall(PetscViewerDestroy(&vv));

 84:     PetscCall(PetscViewerASCIIOpen(PETSC_COMM_WORLD, "daf_7.vtu", &vv));
 85:     PetscCall(VecView(af, vv));
 86:     PetscCall(PetscViewerDestroy(&vv));
 87:     PetscCall(VecDestroy(&ac));
 88:     PetscCall(VecDestroy(&af));
 89:   }
 90:   PetscCall(DMDestroy(&dac));
 91:   PetscCall(DMDestroy(&daf));
 92:   PetscFunctionReturn(PETSC_SUCCESS);
 93: }

 95: int main(int argc, char **argv)
 96: {
 97:   PetscInt mx, my, mz;

 99:   PetscFunctionBeginUser;
100:   PetscCall(PetscInitialize(&argc, &argv, 0, help));
101:   mx = 2;
102:   my = 2;
103:   mz = 2;
104:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-mx", &mx, 0));
105:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-my", &my, 0));
106:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-mz", &mz, 0));
107:   PetscCall(test1_DAInjection3d(mx, my, mz));
108:   PetscCall(PetscFinalize());
109:   return 0;
110: }

112: /*TEST

114:       test:
115:          nsize: 5
116:          args: -mx 30 -my 30 -mz 30 -periodic 0 -da_processors_x 5

118:       test:
119:          suffix: 2
120:          nsize: 5
121:          args: -mx 29 -my 30 -mz 30 -periodic 1 -da_processors_x 5

123:       test:
124:          suffix: 3
125:          nsize: 5
126:          args: -mx 30 -my 29 -mz 30 -periodic 2 -da_processors_x 5

128:       test:
129:          suffix: 4
130:          nsize: 5
131:          args: -mx 30 -my 30 -mz 29 -periodic 3 -da_processors_x 5

133:       test:
134:          suffix: 5
135:          nsize: 5
136:          args: -mx 30 -my 30 -mz 30 -periodic 0 -da_processors_y 5

138:       test:
139:          suffix: 6
140:          nsize: 5
141:          args: -mx 29 -my 30 -mz 30 -periodic 1 -da_processors_y 5

143:       test:
144:          suffix: 7
145:          nsize: 5
146:          args: -mx 30 -my 29 -mz 30 -periodic 2 -da_processors_y 5

148:       test:
149:          suffix: 8
150:          nsize: 5
151:          args: -mx 30 -my 30 -mz 29 -periodic 3 -da_processors_y 5

153:       test:
154:          suffix: 9
155:          nsize: 5
156:          args: -mx 30 -my 30 -mz 30 -periodic 0 -da_processors_z 5

158:       test:
159:          suffix: 10
160:          nsize: 5
161:          args: -mx 29 -my 30 -mz 30 -periodic 1 -da_processors_z 5

163:       test:
164:          suffix: 11
165:          nsize: 5
166:          args: -mx 30 -my 29 -mz 30 -periodic 2 -da_processors_z 5

168:       test:
169:          suffix: 12
170:          nsize: 5
171:          args: -mx 30 -my 30 -mz 29 -periodic 3 -da_processors_z 5

173:       test:
174:          suffix: 13
175:          nsize: 5
176:          args: -mx 30 -my 30 -mz 30 -periodic 0

178:       test:
179:          suffix: 14
180:          nsize: 5
181:          args: -mx 29 -my 30 -mz 30 -periodic 1

183:       test:
184:          suffix: 15
185:          nsize: 5
186:          args: -mx 30 -my 29 -mz 30 -periodic 2

188:       test:
189:          suffix: 16
190:          nsize: 5
191:          args: -mx 30 -my 30 -mz 29 -periodic 3

193: TEST*/