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  | }