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;

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

 21:   periodicity = 0;

 23:   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:   DMSetFromOptions(daf);
 35:   DMSetUp(daf);

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

 39:   DMDASetUniformCoordinates(dac, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0);
 40:   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:     DMGetCoordinateDM(dac, &cdac);
 51:     DMGetCoordinateDM(daf, &cdaf);

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

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

 62:     DMCreateInterpolation(cdac, cdaf, &interp, NULL);
 63:     VecDuplicate(coordsf, &coordsf2);
 64:     MatInterpolate(interp, coordsc, coordsf2);
 65:     VecAXPY(coordsf2, -1.0, coordsf);
 66:     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) PetscPrintf(PETSC_COMM_WORLD, "Norm %g\n", (double)norm);
 69:     VecDestroy(&coordsf2);
 70:     MatDestroy(&interp);
 71:   }

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

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

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

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

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

100:   PetscInitialize(&argc, &argv, 0, help);
101:   mx = 2;
102:   my = 2;
103:   mz = 2;
104:   PetscOptionsGetInt(NULL, NULL, "-mx", &mx, 0);
105:   PetscOptionsGetInt(NULL, NULL, "-my", &my, 0);
106:   PetscOptionsGetInt(NULL, NULL, "-mz", &mz, 0);
107:   test1_DAInjection3d(mx, my, mz);
108:   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*/