1 | /***************************************
2 | $Header: /cvsroot/petscgraphics/3dgf.c,v 1.13 2004/03/07 01:56:18 hazelsct Exp $
3 |
4 | This is a neat 3-D graphing application of Illuminator. Just put whatever
5 | function you like down in line 110 or so (or graph the examples provided), or
6 | run with -twodee and use PETSc's native 2-D graphics (though that would be
7 | BORING!). You might want to run it as:
8 | +latex+\begin{quote}{\tt
9 | +html+ <blockquote><tt>
10 | ./3dgf -da_grid_x 50 -da_grid_y 50 -da_grid_z 50
11 | +latex+}\end{quote}
12 | +html+ </tt></blockquote>
13 | and hit return to end.
14 | ***************************************/
15 |
16 |
17 | static char help[] = "A neat 3-D graphing application of Illuminator.\n\n\
18 | Use -da_grid_x etc. to set the discretization size.\n\
19 | Other options:\n\
20 | -twodee : (obvious)\n\
21 | So you would typically run this using:\n\
22 | 3dgf -da_grid_x 20 -da_grid_y 20 -da_grid_z 20\n";
23 |
24 |
25 | #include "illuminator.h"
26 |
27 |
28 | /* Set #if to 1 to debug, 0 otherwise */
29 | #undef DPRINTF
30 | #if 0
31 | #define DPRINTF(fmt, args...) \
32 | ierr = PetscSynchronizedPrintf (PETSC_COMM_WORLD, "[%d] %s: " fmt, rank, \
33 | __FUNCT__, args); CHKERRQ (ierr); \
34 | ierr = PetscSynchronizedFlush (PETSC_COMM_WORLD); CHKERRQ(ierr)
35 | #else
36 | #define DPRINTF(fmt, args...)
37 | #endif
38 |
39 |
40 | #undef __FUNCT__
41 | #define __FUNCT__ "function_3d"
42 |
43 | /*++++++++++++++++++++++++++++++++++++++
44 | This is where you put the 3-D function you'd like to graph, or the 2-D
45 | function you'd like to graph in 3-D using the zero contour of
46 | +latex+$f(x,y)-z$.
47 | +html+ <i>f</i>(<i>x</i>,<i>y</i>)-<i>z</i>.
48 |
49 | PetscScalar function_3d It returns the function value.
50 |
51 | PetscScalar x The
52 | +latex+$x$-coordinate
53 | +html+ <i>x</i>-coordinate
54 | at which to calculate the function value.
55 |
56 | PetscScalar y The
57 | +latex+$y$-coordinate
58 | +html+ <i>y</i>-coordinate
59 | at which to calculate the function value.
60 |
61 | PetscScalar z The
62 | +latex+$z$-coordinate
63 | +html+ <i>z</i>-coordinate
64 | at which to calculate the function value.
65 | ++++++++++++++++++++++++++++++++++++++*/
66 |
67 | static inline PetscScalar function_3d
68 | (PetscScalar x, PetscScalar y, PetscScalar z)
69 | {
70 | /* A real simple function for testing */
71 | /* return x+y+z; */
72 |
73 | /* The potential Green's function in 3-D: -1/(4 pi r), with one
74 | octant sliced out */
75 | if (x<0 || y<0 || z<0)
76 | return -.25/sqrt(x*x+y*y+z*z)/M_PI;
77 | else
78 | return 1000000000;
79 |
80 | /* The x-derivative of that Green's function: x/(4 pi r^3) */
81 | /* return .25*x/sqrt(x*x+y*y+z*z)/(x*x+y*y+z*z)/M_PI; */
82 |
83 | /* Demo of graphing z=f(x,y): the x-derivative of the 2-D Green's
84 | function; you need to modify the DATriangulate call below
85 | to plot one contour at f=0 */
86 | /* return x/(x*x+y*y)/2./M_PI - z; */
87 | }
88 |
89 |
90 | #undef __FUNCT__
91 | #define __FUNCT__ "function_2d"
92 |
93 | /*++++++++++++++++++++++++++++++++++++++
94 | This is where you put the 2-D function you'd like to graph using PETSc's
95 | native 2-D "contour" color graphics.
96 |
97 | PetscScalar function_2d It returns the function value.
98 |
99 | PetscScalar x The
100 | +latex+$x$-coordinate
101 | +html+ <i>x</i>-coordinate
102 | at which to calculate the function value.
103 |
104 | PetscScalar y The
105 | +latex+$y$-coordinate
106 | +html+ <i>y</i>-coordinate
107 | at which to calculate the function value.
108 | ++++++++++++++++++++++++++++++++++++++*/
109 |
110 | static inline PetscScalar function_2d (PetscScalar x, PetscScalar y)
111 | {
112 | /* The potential Green's function in 2-D: ln(r)/(2 pi) */
113 | return -log(1./sqrt(x*x+y*y))/2./M_PI;
114 |
115 | /* And its x-derivative: x/(2 pi r^2) */
116 | /* return x/(x*x+y*y)/2./M_PI; */
117 | }
118 |
119 |
120 | #undef __FUNCT__
121 | #define __FUNCT__ "main"
122 |
123 | /*++++++++++++++++++++++++++++++++++++++
124 | The usual main function.
125 |
126 | int main Returns 0 or error.
127 |
128 | int argc Number of args.
129 |
130 | char **argv The args.
131 | ++++++++++++++++++++++++++++++++++++++*/
132 |
133 | int main (int argc, char **argv)
134 | {
135 | DA da;
136 | Vec vector; /* solution vector */
137 | int xstart,ystart,zstart, xwidth,ywidth,zwidth, xglobal,yglobal,
138 | zglobal, i,j,k, rank, ierr;
139 | PetscScalar xmin=-.8,xmax=.8, ymin=-.8,ymax=.8, zmin=-.8,zmax=.8;
140 | PetscTruth threedee;
141 | PetscViewer theviewer;
142 |
143 | /*+The program first calls
144 | +latex+{\tt PetscInitialize()}
145 | +html+ <tt>PetscInitialize()</tt>
146 | and creates the distributed arrays. Note that even though this program
147 | doesn't do any communication between the CPUs, illuminator must do so in
148 | order to make the isoquants at CPU boundaries, so the stencil width must be
149 | at least one. +*/
150 | ierr = PetscInitialize (&argc, &argv, (char *)0, help); CHKERRQ (ierr);
151 | ierr = MPI_Comm_rank (PETSC_COMM_WORLD, &rank); CHKERRQ (ierr);
152 | DPRINTF ("Petsc initialized, moving forward\n", 0);
153 |
154 | /* Create DA */
155 | ierr = PetscOptionsHasName (PETSC_NULL, "-twodee", &threedee);
156 | threedee = !threedee;
157 | CHKERRQ (ierr);
158 |
159 | if (threedee)
160 | {
161 | ierr = DACreate3d
162 | (PETSC_COMM_WORLD, DA_NONPERIODIC, DA_STENCIL_BOX, PETSC_DECIDE,
163 | PETSC_DECIDE,PETSC_DECIDE, PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,
164 | 1, 1, PETSC_NULL,PETSC_NULL,PETSC_NULL, &da); CHKERRQ (ierr);
165 | }
166 | else
167 | {
168 | ierr = DACreate2d
169 | (PETSC_COMM_WORLD, DA_NONPERIODIC, DA_STENCIL_BOX, PETSC_DECIDE,
170 | PETSC_DECIDE, PETSC_DECIDE,PETSC_DECIDE, 1, 1, PETSC_NULL,PETSC_NULL,
171 | &da); CHKERRQ (ierr);
172 | }
173 |
174 | /*+Next it gets the distributed array's local corner and global size
175 | information. It gets the global vector, and loops over the part stored on
176 | this CPU to set all of the function values, using
177 | +latex+{\tt function\_3d()} or {\tt function\_2d()}
178 | +html+ <tt>function_3d()</tt> or <tt>function_2d()</tt>
179 | depending on whether the
180 | +latex+{\tt -twodee}
181 | +html+ <tt>-twodee</tt>
182 | command line switch was used at runtime. +*/
183 | ierr = DAGetCorners (da, &xstart, &ystart, &zstart, &xwidth, &ywidth,
184 | &zwidth); CHKERRQ (ierr);
185 | ierr = DAGetInfo (da, PETSC_NULL, &xglobal, &yglobal, &zglobal, PETSC_NULL,
186 | PETSC_NULL, PETSC_NULL, PETSC_NULL, PETSC_NULL, PETSC_NULL,
187 | PETSC_NULL); CHKERRQ (ierr);
188 | if (xglobal == 1 && yglobal == 1 && zglobal == 1)
189 | {
190 | ierr = PetscPrintf (PETSC_COMM_WORLD, "Grid dimensions 1x1x1 indicate you probably want to use -da_grid_x etc.\n");
191 | CHKERRQ (ierr);
192 | }
193 | ierr = DAGetGlobalVector (da, &vector); CHKERRQ (ierr);
194 |
195 | DPRINTF ("About to calculate function values\n", 0);
196 |
197 | if (threedee)
198 | {
199 | PetscScalar ***array3d;
200 |
201 | ierr = VecGetArray3d (vector, zwidth, ywidth, xwidth, zstart, ystart,
202 | xstart, &array3d); CHKERRQ (ierr);
203 | for (k=zstart; k<zstart+zwidth; k++)
204 | for (j=ystart; j<ystart+ywidth; j++)
205 | for (i=xstart; i<xstart+xwidth; i++)
206 | {
207 | PetscScalar x,y,z;
208 |
209 | x = xmin + (xmax-xmin) * (double) i/(xglobal-1);
210 | y = ymin + (ymax-ymin) * (double) j/(yglobal-1);
211 | z = zmin + (zmax-zmin) * (double) k/(zglobal-1);
212 |
213 | array3d [k][j][i] = function_3d (x, y, z);
214 | }
215 | ierr = VecRestoreArray3d (vector, zwidth, ywidth, xwidth, zstart, ystart,
216 | xstart, &array3d); CHKERRQ (ierr);
217 | }
218 | else
219 | {
220 | PetscScalar **array2d;
221 |
222 | ierr = VecGetArray2d (vector, ywidth, xwidth, ystart, xstart, &array2d);
223 | CHKERRQ (ierr);
224 | DASetFieldName (da, 0, "2-D Potential Green's Function");
225 | for (j=ystart; j<ystart+ywidth; j++)
226 | for (i=xstart; i<xstart+xwidth; i++)
227 | {
228 | double x,y;
229 |
230 | x = xmin + (xmax-xmin) * (double) i/(xglobal-1);
231 | y = ymin + (ymax-ymin) * (double) j/(yglobal-1);
232 |
233 | array2d [j][i] = function_2d (x,y);
234 | }
235 | ierr = VecRestoreArray2d (vector, ywidth, xwidth, ystart, xstart,
236 | &array2d); CHKERRQ (ierr);
237 | }
238 |
239 | /*+It then uses
240 | +latex+{\tt GeomviewBegin()} or {\tt PetscViewerDrawOpen()}
241 | +html+ <tt>GeomviewBegin()</tt> or <tt>PetscViewerDrawOpen()</tt>
242 | to start the viewer, and either
243 | +latex+{\tt DATriangulate()} and {\tt GeomviewDisplayTriangulation()} or
244 | +latex+{\tt VecView()}
245 | +html+ <tt>DATriangulate()</tt> and <tt>GeomviewDisplayTriangulation()</tt>
246 | +html+ or <tt>VecView()</tt>
247 | to display the solution. +*/
248 |
249 | DPRINTF ("About to open display\n", 0);
250 |
251 | if (threedee)
252 | {
253 | ierr = GeomviewBegin (PETSC_COMM_WORLD); CHKERRQ (ierr);
254 | }
255 | else
256 | {
257 | PetscDraw draw;
258 | ierr = PetscViewerDrawOpen (PETSC_COMM_WORLD, 0, "", PETSC_DECIDE,
259 | PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE,
260 | &theviewer); CHKERRQ (ierr);
261 | ierr = PetscViewerDrawGetDraw (theviewer, 0, &draw);
262 | CHKERRQ (ierr);
263 | ierr = PetscDrawSetDoubleBuffer (draw); CHKERRQ (ierr);
264 | }
265 |
266 | DPRINTF ("About to do the graphing\n", 0);
267 |
268 | if (threedee)
269 | {
270 | PetscScalar minmax[6] = { xmin,xmax, ymin,ymax, zmin,zmax };
271 | PetscScalar isovals[4] = { -.1, -.2, -.3, -.4 };
272 | PetscScalar colors[16] = { 1,0,0,.4, 1,1,0,.4, 0,1,0,.4, 0,0,1,.4 };
273 |
274 | DPRINTF ("Starting triangulation\n", 0);
275 | ierr = DATriangulate (da, vector, 0, minmax, 4, isovals, colors,
276 | PETSC_FALSE, PETSC_FALSE, PETSC_FALSE);
277 | CHKERRQ (ierr);
278 | DPRINTF ("Sending to Geomview\n", 0);
279 | ierr = GeomviewDisplayTriangulation
280 | (PETSC_COMM_WORLD, minmax, "Function-Contours", PETSC_TRUE);
281 | CHKERRQ (ierr);
282 | }
283 | else
284 | {
285 | ierr = VecView (vector, theviewer); CHKERRQ (ierr);
286 | }
287 |
288 | /*+ Finally, it prompts the user to hit <return> before wrapping up. +*/
289 |
290 | {
291 | char instring[100];
292 |
293 | ierr = PetscPrintf (PETSC_COMM_WORLD, "Press <return> to close up... ");
294 | CHKERRQ (ierr);
295 | ierr = PetscSynchronizedFGets (PETSC_COMM_WORLD, stdin, 99, instring);
296 | CHKERRQ (ierr);
297 | }
298 |
299 | DPRINTF ("Cleaning up\n", 0);
300 | if (threedee)
301 | {
302 | ierr = GeomviewEnd (PETSC_COMM_WORLD); CHKERRQ (ierr);
303 | }
304 | else
305 | {
306 | ierr = PetscViewerDestroy (theviewer); CHKERRQ (ierr);
307 | }
308 | ierr = DARestoreGlobalVector (da, &vector); CHKERRQ (ierr);
309 | ierr = DADestroy (da); CHKERRQ (ierr);
310 |
311 | PetscFinalize ();
312 | return 0;
313 | }