xref: /petsc/src/dm/tests/ex21.c (revision 503c0ea9b45bcfbcebbb1ea5341243bbc69f0bea)
1 static const char help[] = "Test DMCreateInjection() for mapping coordinates in 3D";
2 
3 #include <petscvec.h>
4 #include <petscmat.h>
5 #include <petscdm.h>
6 #include <petscdmda.h>
7 
8 PetscErrorCode test1_DAInjection3d(PetscInt mx, PetscInt my, PetscInt mz)
9 {
10   PetscErrorCode   ierr;
11   DM               dac,daf;
12   PetscViewer      vv;
13   Vec              ac,af;
14   PetscInt         periodicity;
15   DMBoundaryType   bx,by,bz;
16 
17   PetscFunctionBeginUser;
18   bx = DM_BOUNDARY_NONE;
19   by = DM_BOUNDARY_NONE;
20   bz = DM_BOUNDARY_NONE;
21 
22   periodicity = 0;
23 
24   PetscCall(PetscOptionsGetInt(NULL,NULL,"-periodic", &periodicity, NULL));
25   if (periodicity==1) {
26     bx = DM_BOUNDARY_PERIODIC;
27   } else if (periodicity==2) {
28     by = DM_BOUNDARY_PERIODIC;
29   } else if (periodicity==3) {
30     bz = DM_BOUNDARY_PERIODIC;
31   }
32 
33   ierr = 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 */
34                       1, /* stencil = 1 */NULL,NULL,NULL,&daf);PetscCall(ierr);
35   PetscCall(DMSetFromOptions(daf));
36   PetscCall(DMSetUp(daf));
37 
38   PetscCall(DMCoarsen(daf,MPI_COMM_NULL,&dac));
39 
40   PetscCall(DMDASetUniformCoordinates(dac, -1.0,1.0, -1.0,1.0, -1.0,1.0));
41   PetscCall(DMDASetUniformCoordinates(daf, -1.0,1.0, -1.0,1.0, -1.0,1.0));
42 
43   {
44     DM         cdaf,cdac;
45     Vec        coordsc,coordsf,coordsf2;
46     Mat        inject;
47     VecScatter vscat;
48     Mat        interp;
49     PetscReal  norm;
50 
51     PetscCall(DMGetCoordinateDM(dac,&cdac));
52     PetscCall(DMGetCoordinateDM(daf,&cdaf));
53 
54     PetscCall(DMGetCoordinates(dac,&coordsc));
55     PetscCall(DMGetCoordinates(daf,&coordsf));
56 
57     PetscCall(DMCreateInjection(cdac,cdaf,&inject));
58     PetscCall(MatScatterGetVecScatter(inject,&vscat));
59     PetscCall(VecScatterBegin(vscat,coordsf,coordsc,INSERT_VALUES,SCATTER_FORWARD));
60     PetscCall(VecScatterEnd(vscat  ,coordsf,coordsc,INSERT_VALUES,SCATTER_FORWARD));
61     PetscCall(MatDestroy(&inject));
62 
63     PetscCall(DMCreateInterpolation(cdac,cdaf,&interp,NULL));
64     PetscCall(VecDuplicate(coordsf,&coordsf2));
65     PetscCall(MatInterpolate(interp,coordsc,coordsf2));
66     PetscCall(VecAXPY(coordsf2,-1.0,coordsf));
67     PetscCall(VecNorm(coordsf2,NORM_MAX,&norm));
68     /* The fine coordinates are only reproduced in certain cases */
69     if (!bx && !by && !bz && norm > PETSC_SQRT_MACHINE_EPSILON) PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Norm %g\n",(double)norm));
70     PetscCall(VecDestroy(&coordsf2));
71     PetscCall(MatDestroy(&interp));
72   }
73 
74   if (0) {
75     PetscCall(DMCreateGlobalVector(dac,&ac));
76     PetscCall(VecZeroEntries(ac));
77 
78     PetscCall(DMCreateGlobalVector(daf,&af));
79     PetscCall(VecZeroEntries(af));
80 
81     PetscCall(PetscViewerASCIIOpen(PETSC_COMM_WORLD, "dac_7.vtu", &vv));
82     PetscCall(VecView(ac, vv));
83     PetscCall(PetscViewerDestroy(&vv));
84 
85     PetscCall(PetscViewerASCIIOpen(PETSC_COMM_WORLD, "daf_7.vtu", &vv));
86     PetscCall(VecView(af, vv));
87     PetscCall(PetscViewerDestroy(&vv));
88     PetscCall(VecDestroy(&ac));
89     PetscCall(VecDestroy(&af));
90   }
91   PetscCall(DMDestroy(&dac));
92   PetscCall(DMDestroy(&daf));
93   PetscFunctionReturn(0);
94 }
95 
96 int main(int argc,char **argv)
97 {
98   PetscInt       mx,my,mz;
99 
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 }
111 
112 /*TEST
113 
114       test:
115          nsize: 5
116          args: -mx 30 -my 30 -mz 30 -periodic 0 -da_processors_x 5
117 
118       test:
119          suffix: 2
120          nsize: 5
121          args: -mx 29 -my 30 -mz 30 -periodic 1 -da_processors_x 5
122 
123       test:
124          suffix: 3
125          nsize: 5
126          args: -mx 30 -my 29 -mz 30 -periodic 2 -da_processors_x 5
127 
128       test:
129          suffix: 4
130          nsize: 5
131          args: -mx 30 -my 30 -mz 29 -periodic 3 -da_processors_x 5
132 
133       test:
134          suffix: 5
135          nsize: 5
136          args: -mx 30 -my 30 -mz 30 -periodic 0 -da_processors_y 5
137 
138       test:
139          suffix: 6
140          nsize: 5
141          args: -mx 29 -my 30 -mz 30 -periodic 1 -da_processors_y 5
142 
143       test:
144          suffix: 7
145          nsize: 5
146          args: -mx 30 -my 29 -mz 30 -periodic 2 -da_processors_y 5
147 
148       test:
149          suffix: 8
150          nsize: 5
151          args: -mx 30 -my 30 -mz 29 -periodic 3 -da_processors_y 5
152 
153       test:
154          suffix: 9
155          nsize: 5
156          args: -mx 30 -my 30 -mz 30 -periodic 0 -da_processors_z 5
157 
158       test:
159          suffix: 10
160          nsize: 5
161          args: -mx 29 -my 30 -mz 30 -periodic 1 -da_processors_z 5
162 
163       test:
164          suffix: 11
165          nsize: 5
166          args: -mx 30 -my 29 -mz 30 -periodic 2 -da_processors_z 5
167 
168       test:
169          suffix: 12
170          nsize: 5
171          args: -mx 30 -my 30 -mz 29 -periodic 3 -da_processors_z 5
172 
173       test:
174          suffix: 13
175          nsize: 5
176          args: -mx 30 -my 30 -mz 30 -periodic 0
177 
178       test:
179          suffix: 14
180          nsize: 5
181          args: -mx 29 -my 30 -mz 30 -periodic 1
182 
183       test:
184          suffix: 15
185          nsize: 5
186          args: -mx 30 -my 29 -mz 30 -periodic 2
187 
188       test:
189          suffix: 16
190          nsize: 5
191          args: -mx 30 -my 30 -mz 29 -periodic 3
192 
193 TEST*/
194