Actual source code: ex48.c
petsc-3.10.2 2018-10-09
1: static char help[] = "Test VTK structured (.vts) and rectilinear (.vtr) viewer support with multi-dof DMDAs\n\n";
3: #include <petscdm.h>
4: #include <petscdmda.h>
6: /*
7: Write 3D DMDA vector with coordinates in VTK format
8: */
9: PetscErrorCode test_3d(const char filename[],PetscInt dof)
10: {
11: MPI_Comm comm = MPI_COMM_WORLD;
12: const PetscInt M=10,N=15,P=30,sw=1;
13: const PetscScalar Lx=1.0,Ly=1.0,Lz=1.0;
14: DM da;
15: Vec v;
16: PetscViewer view;
17: DMDALocalInfo info;
18: PetscScalar ****va;
19: PetscInt i,j,k,c;
20: PetscErrorCode ierr;
22: DMDACreate3d(comm,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR, M,N,P,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,dof,sw,NULL,NULL,NULL,&da);
23: DMSetFromOptions(da);
24: DMSetUp(da);
26: DMDASetUniformCoordinates(da,0.0,Lx,0.0,Ly,0.0,Lz);
27: DMDAGetLocalInfo(da,&info);
28: DMCreateGlobalVector(da,&v);
29: DMDAVecGetArrayDOF(da,v,&va);
30: for (k=info.zs; k<info.zs+info.zm; k++) {
31: for (j=info.ys; j<info.ys+info.ym; j++) {
32: for (i=info.xs; i<info.xs+info.xm; i++) {
33: const PetscScalar x = (Lx*i)/M;
34: const PetscScalar y = (Ly*j)/N;
35: const PetscScalar z = (Lz*k)/P;
36: for (c=0; c<dof; ++ c) {
37: va[k][j][i][c] = PetscPowScalarInt(x-0.5*Lx,2)+PetscPowScalarInt(y-0.5*Ly,2)+PetscPowScalarInt(z-0.5*Lz,2) + 10.0*c;
38: }
39: }
40: }
41: }
42: DMDAVecRestoreArrayDOF(da,v,&va);
43: PetscViewerVTKOpen(comm,filename,FILE_MODE_WRITE,&view);
44: VecView(v,view);
45: PetscViewerDestroy(&view);
46: VecDestroy(&v);
47: DMDestroy(&da);
48: return 0;
49: }
52: /*
53: Write 2D DMDA vector with coordinates in VTK format
54: */
55: PetscErrorCode test_2d(const char filename[],PetscInt dof)
56: {
57: MPI_Comm comm = MPI_COMM_WORLD;
58: const PetscInt M=10,N=20,sw=1;
59: const PetscScalar Lx=1.0,Ly=1.0,Lz=1.0;
60: DM da;
61: Vec v;
62: PetscViewer view;
63: DMDALocalInfo info;
64: PetscScalar ***va;
65: PetscInt i,j,c;
66: PetscErrorCode ierr;
68: DMDACreate2d(comm,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR, M,N,PETSC_DECIDE,PETSC_DECIDE,dof,sw,NULL,NULL,&da);
69: DMSetFromOptions(da);
70: DMSetUp(da);
71: DMDASetUniformCoordinates(da,0.0,Lx,0.0,Ly,0.0,Lz);
72: DMDAGetLocalInfo(da,&info);
73: DMCreateGlobalVector(da,&v);
74: DMDAVecGetArrayDOF(da,v,&va);
75: for (j=info.ys; j<info.ys+info.ym; j++) {
76: for (i=info.xs; i<info.xs+info.xm; i++) {
77: const PetscScalar x = (Lx*i)/M;
78: const PetscScalar y = (Ly*j)/N;
79: for (c=0; c<dof; ++c) {
80: va[j][i][c] = PetscPowScalarInt(x-0.5*Lx,2)+PetscPowScalarInt(y-0.5*Ly,2) + 10.0*c;
81: }
82: }
83: }
84: DMDAVecRestoreArrayDOF(da,v,&va);
85: PetscViewerVTKOpen(comm,filename,FILE_MODE_WRITE,&view);
86: VecView(v,view);
87: PetscViewerDestroy(&view);
88: VecDestroy(&v);
89: DMDestroy(&da);
90: return 0;
91: }
93: /*
94: Write a scalar and a vector field from two compatible 3d DMDAs
95: */
96: PetscErrorCode test_3d_compat(const char filename[],PetscInt dof)
97: {
98: MPI_Comm comm = MPI_COMM_WORLD;
99: const PetscInt M=10,N=15,P=30,sw=1;
100: const PetscScalar Lx=1.0,Ly=1.0,Lz=1.0;
101: DM da,daVector;
102: Vec v,vVector;
103: PetscViewer view;
104: DMDALocalInfo info;
105: PetscScalar ***va,****vVectora;
106: PetscInt i,j,k,c;
107: PetscErrorCode ierr;
109: DMDACreate3d(comm,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR, M,N,P,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,/* dof:*/1,sw,NULL,NULL,NULL,&da);
110: DMSetFromOptions(da);
111: DMSetUp(da);
113: DMDASetUniformCoordinates(da,0.0,Lx,0.0,Ly,0.0,Lz);
114: DMDAGetLocalInfo(da,&info);
115: DMDACreateCompatibleDMDA(da,dof,&daVector);
116: DMCreateGlobalVector(da,&v);
117: DMCreateGlobalVector(daVector,&vVector);
118: DMDAVecGetArray(da,v,&va);
119: DMDAVecGetArrayDOF(daVector,vVector,&vVectora);
120: for (k=info.zs; k<info.zs+info.zm; k++) {
121: for (j=info.ys; j<info.ys+info.ym; j++) {
122: for (i=info.xs; i<info.xs+info.xm; i++) {
123: const PetscScalar x = (Lx*i)/M;
124: const PetscScalar y = (Ly*j)/N;
125: const PetscScalar z = (Lz*k)/P;
126: va[k][j][i] = PetscPowScalarInt(x-0.5*Lx,2)+PetscPowScalarInt(y-0.5*Ly,2)+PetscPowScalarInt(z-0.5*Lz,2);
127: for (c=0; c<dof; ++c) {
128: vVectora[k][j][i][c] = PetscPowScalarInt(x-0.5*Lx,2)+PetscPowScalarInt(y-0.5*Ly,2)+PetscPowScalarInt(z-0.5*Lz,2) + 10.0*c;
129: }
130: }
131: }
132: }
133: DMDAVecRestoreArray(da,v,&va);
134: DMDAVecRestoreArrayDOF(da,v,&vVectora);
135: PetscViewerVTKOpen(comm,filename,FILE_MODE_WRITE,&view);
136: VecView(v,view);
137: VecView(vVector,view);
138: PetscViewerDestroy(&view);
139: VecDestroy(&v);
140: VecDestroy(&vVector);
141: DMDestroy(&da);
142: DMDestroy(&daVector);
143: return 0;
144: }
146: /*
147: Write a scalar and a vector field from two compatible 2d DMDAs
148: */
149: PetscErrorCode test_2d_compat(const char filename[],PetscInt dof)
150: {
151: MPI_Comm comm = MPI_COMM_WORLD;
152: const PetscInt M=10,N=20,sw=1;
153: const PetscScalar Lx=1.0,Ly=1.0,Lz=1.0;
154: DM da, daVector;
155: Vec v,vVector;
156: PetscViewer view;
157: DMDALocalInfo info;
158: PetscScalar **va,***vVectora;
159: PetscInt i,j,c;
160: PetscErrorCode ierr;
162: DMDACreate2d(comm,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR, M,N,PETSC_DECIDE,PETSC_DECIDE,/* dof:*/ 1,sw,NULL,NULL,&da);
163: DMSetFromOptions(da);
164: DMSetUp(da);
165: DMDASetUniformCoordinates(da,0.0,Lx,0.0,Ly,0.0,Lz);
166: DMDACreateCompatibleDMDA(da,dof,&daVector);
167: DMDAGetLocalInfo(da,&info);
168: DMCreateGlobalVector(da,&v);
169: DMCreateGlobalVector(daVector,&vVector);
170: DMDAVecGetArray(da,v,&va);
171: DMDAVecGetArrayDOF(daVector,vVector,&vVectora);
172: for (j=info.ys; j<info.ys+info.ym; j++) {
173: for (i=info.xs; i<info.xs+info.xm; i++) {
174: const PetscScalar x = (Lx*i)/M;
175: const PetscScalar y = (Ly*j)/N;
176: va[j][i] = PetscPowScalarInt(x-0.5*Lx,2)+PetscPowScalarInt(y-0.5*Ly,2);
177: for (c=0; c<dof; ++c) {
178: vVectora[j][i][c] = PetscPowScalarInt(x-0.5*Lx,2)+PetscPowScalarInt(y-0.5*Ly,2) + 10.0*c;
179: }
180: }
181: }
182: DMDAVecRestoreArray(da,v,&va);
183: DMDAVecRestoreArrayDOF(daVector,vVector,&vVectora);
184: PetscViewerVTKOpen(comm,filename,FILE_MODE_WRITE,&view);
185: VecView(v,view);
186: VecView(vVector,view);
187: PetscViewerDestroy(&view);
188: VecDestroy(&v);
189: VecDestroy(&vVector);
190: DMDestroy(&da);
191: DMDestroy(&daVector);
192: return 0;
193: }
195: int main(int argc, char *argv[])
196: {
198: PetscInt dof;
200: PetscInitialize(&argc,&argv,0,help);if (ierr) return ierr;
201: dof = 2;
202: PetscOptionsGetInt(NULL,NULL,"-dof",&dof,NULL);
203: test_3d("3d.vtr",dof);
204: test_2d("2d.vtr",dof);
205: test_3d_compat("3d_compat.vtr",dof);
206: test_2d_compat("2d_compat.vtr",dof);
207: test_3d("3d.vts",dof);
208: test_2d("2d.vts",dof);
209: test_3d_compat("3d_compat.vts",dof);
210: test_2d_compat("2d_compat.vts",dof);
211: PetscFinalize();
212: return ierr;
213: }
215: /*TEST
217: build:
218: requires: !complex
220: test:
221: suffix: 1
222: nsize: 2
223: args: -dof 2
225: test:
226: suffix: 2
227: nsize: 2
228: args: -dof 2
230: test:
231: suffix: 3
232: nsize: 2
233: args: -dof 3
235: TEST*/