Actual source code: ex21.c
petsc-3.3-p6 2013-02-11
1: static const char help[] = "Test DMCreateInjection for mapping coordinates in 3D";
3: #include <petscvec.h>
4: #include <petscmat.h>
5: #include <petscdmda.h>
9: PetscErrorCode test1_DAInjection3d( PetscInt mx, PetscInt my, PetscInt mz )
10: {
12: DM dac,daf;
13: PetscViewer vv;
14: Vec ac,af;
15: PetscInt periodicity;
16: DMDABoundaryType bx,by,bz;
19: bx = DMDA_BOUNDARY_NONE;
20: by = DMDA_BOUNDARY_NONE;
21: bz = DMDA_BOUNDARY_NONE;
22: periodicity = 0;
23: PetscOptionsGetInt(PETSC_NULL,"-periodic", &periodicity, PETSC_NULL);
24: if (periodicity==1) {
25: bx = DMDA_BOUNDARY_PERIODIC;
26: } else if (periodicity==2) {
27: by = DMDA_BOUNDARY_PERIODIC;
28: } else if (periodicity==3) {
29: bz = DMDA_BOUNDARY_PERIODIC;
30: }
32: DMDACreate3d(PETSC_COMM_WORLD, bx,by,bz, DMDA_STENCIL_BOX,
33: mx+1, my+1,mz+1,
34: PETSC_DECIDE, PETSC_DECIDE,PETSC_DECIDE,
35: 1, /* 1 dof */
36: 1, /* stencil = 1 */
37: PETSC_NULL,PETSC_NULL,PETSC_NULL,
38: &daf);
40: DMSetFromOptions(daf);
42: DMCoarsen(daf,MPI_COMM_NULL,&dac);
44: DMDASetUniformCoordinates(dac, -1.0,1.0, -1.0,1.0, -1.0,1.0 );
45: DMDASetUniformCoordinates(daf, -1.0,1.0, -1.0,1.0, -1.0,1.0 );
47: {
48: DM cdaf,cdac;
49: Vec coordsc,coordsf,coordsf2;
50: VecScatter inject;
51: Mat interp;
52: PetscReal norm;
54: DMDAGetCoordinateDA(dac,&cdac);
55: DMDAGetCoordinateDA(daf,&cdaf);
57: DMDAGetCoordinates(dac,&coordsc);
58: DMDAGetCoordinates(daf,&coordsf);
60: DMCreateInjection(cdac,cdaf,&inject);
62: VecScatterBegin(inject,coordsf,coordsc,INSERT_VALUES,SCATTER_FORWARD);
63: VecScatterEnd(inject ,coordsf,coordsc,INSERT_VALUES,SCATTER_FORWARD);
64: VecScatterDestroy(&inject);
66: DMCreateInterpolation(cdac,cdaf,&interp,PETSC_NULL);
67: VecDuplicate(coordsf,&coordsf2);
68: MatInterpolate(interp,coordsc,coordsf2);
69: VecAXPY(coordsf2,-1.0,coordsf);
70: VecNorm(coordsf2,NORM_MAX,&norm);
71: /* The fine coordinates are only reproduced in certain cases */
72: if (!bx && !by && !bz && norm > 1.e-10) {PetscPrintf(PETSC_COMM_WORLD,"Norm %G\n",norm);}
73: VecDestroy(&coordsf2);
74: MatDestroy(&interp);
75: }
77: if (0) {
78: DMCreateGlobalVector(dac,&ac);
79: VecZeroEntries(ac);
81: DMCreateGlobalVector(daf,&af);
82: VecZeroEntries(af);
84: PetscViewerASCIIOpen(PETSC_COMM_WORLD, "dac_7.vtk", &vv);
85: PetscViewerSetFormat(vv, PETSC_VIEWER_ASCII_VTK);
86: DMView(dac, vv);
87: VecView(ac, vv);
88: PetscViewerDestroy(&vv);
90: PetscViewerASCIIOpen(PETSC_COMM_WORLD, "daf_7.vtk", &vv);
91: PetscViewerSetFormat(vv, PETSC_VIEWER_ASCII_VTK);
92: DMView(daf, vv);
93: VecView(af, vv);
94: PetscViewerDestroy(&vv);
95: VecDestroy(&ac);
96: VecDestroy(&af);
97: }
98: DMDestroy(&dac);
99: DMDestroy(&daf);
100: return(0);
101: }
105: int main(int argc,char **argv)
106: {
108: PetscInt mx,my,mz;
110: PetscInitialize(&argc,&argv,0,0);
111: mx = 2;
112: my = 2;
113: mz = 2;
114: PetscOptionsGetInt(PETSC_NULL,"-mx", &mx, 0 );
115: PetscOptionsGetInt(PETSC_NULL,"-my", &my, 0 );
116: PetscOptionsGetInt(PETSC_NULL,"-mz", &mz, 0 );
118: test1_DAInjection3d(mx,my,mz);
120: PetscFinalize();
121: return 0;
122: }