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