1    | /***************************************
2    |   $Header$
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  |   ISurface      Surf;
143  |   IDisplay      Disp;
144  | 
145  |   /*+The program first calls
146  |     +latex+{\tt PetscInitialize()}
147  |     +html+ <tt>PetscInitialize()</tt>
148  |     and creates the distributed arrays.  Note that even though this program
149  |     doesn't do any communication between the CPUs, illuminator must do so in
150  |     order to make the isoquants at CPU boundaries, so the stencil width must be
151  |     at least one. +*/
152  |   ierr = PetscInitialize (&argc, &argv, (char *)0, help); CHKERRQ (ierr);
153  |   ierr = MPI_Comm_rank (PETSC_COMM_WORLD, &rank); CHKERRQ (ierr);
154  |   DPRINTF ("Petsc initialized, moving forward\n", 0);
155  | 
156  |   /* Create DA */
157  |   ierr = PetscOptionsHasName (PETSC_NULL, "-twodee", &threedee);
158  |   threedee = !threedee;
159  |   CHKERRQ (ierr);
160  | 
161  |   if (threedee)
162  |     {
163  |       ierr = DACreate3d
164  | 	(PETSC_COMM_WORLD, DA_NONPERIODIC, DA_STENCIL_BOX, PETSC_DECIDE,
165  | 	 PETSC_DECIDE,PETSC_DECIDE, PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,
166  | 	 1, 1, PETSC_NULL,PETSC_NULL,PETSC_NULL, &da); CHKERRQ (ierr);
167  |     }
168  |   else
169  |     {
170  |       ierr = DACreate2d
171  | 	(PETSC_COMM_WORLD, DA_NONPERIODIC, DA_STENCIL_BOX, PETSC_DECIDE,
172  | 	 PETSC_DECIDE, PETSC_DECIDE,PETSC_DECIDE, 1, 1, PETSC_NULL,PETSC_NULL,
173  | 	 &da); CHKERRQ (ierr);
174  |     }
175  | 
176  |   /*+Next it gets the distributed array's local corner and global size
177  |     information.  It gets the global vector, and loops over the part stored on
178  |     this CPU to set all of the function values, using
179  |     +latex+{\tt function\_3d()} or {\tt function\_2d()}
180  |     +html+ <tt>function_3d()</tt> or <tt>function_2d()</tt>
181  |     depending on whether the
182  |     +latex+{\tt -twodee}
183  |     +html+ <tt>-twodee</tt>
184  |     command line switch was used at runtime. +*/
185  |   ierr = DAGetCorners (da, &xstart, &ystart, &zstart, &xwidth, &ywidth,
186  | 		       &zwidth); CHKERRQ (ierr);
187  |   ierr = DAGetInfo (da, PETSC_NULL, &xglobal, &yglobal, &zglobal, PETSC_NULL,
188  | 		    PETSC_NULL, PETSC_NULL, PETSC_NULL, PETSC_NULL, PETSC_NULL,
189  | 		    PETSC_NULL); CHKERRQ (ierr);
190  |   if (xglobal == 1 && yglobal == 1 && zglobal == 1)
191  |     {
192  |       ierr = PetscPrintf (PETSC_COMM_WORLD, "Grid dimensions 1x1x1 indicate you probably want to use -da_grid_x etc.\n");
193  |       CHKERRQ (ierr);
194  |     }
195  |   ierr = DAGetGlobalVector (da, &vector); CHKERRQ (ierr);
196  | 
197  |   DPRINTF ("About to calculate function values\n", 0);
198  | 
199  |   if (threedee)
200  |     {
201  |       PetscScalar ***array3d;
202  | 
203  |       ierr = VecGetArray3d (vector, zwidth, ywidth, xwidth, zstart, ystart,
204  | 			    xstart, &array3d); CHKERRQ (ierr);
205  |       for (k=zstart; k<zstart+zwidth; k++)
206  | 	for (j=ystart; j<ystart+ywidth; j++)
207  | 	  for (i=xstart; i<xstart+xwidth; i++)
208  | 	    {
209  | 	      PetscScalar x,y,z;
210  | 
211  | 	      x = xmin + (xmax-xmin) * (double) i/(xglobal-1);
212  | 	      y = ymin + (ymax-ymin) * (double) j/(yglobal-1);
213  | 	      z = zmin + (zmax-zmin) * (double) k/(zglobal-1);
214  | 
215  | 	      array3d [k][j][i] = function_3d (x, y, z);
216  | 	    }
217  |       ierr = VecRestoreArray3d (vector, zwidth, ywidth, xwidth, zstart, ystart,
218  | 				xstart, &array3d); CHKERRQ (ierr);
219  |     }
220  |   else
221  |     {
222  |       PetscScalar **array2d;
223  | 
224  |       ierr = VecGetArray2d (vector, ywidth, xwidth, ystart, xstart, &array2d);
225  |       CHKERRQ (ierr);
226  |       DASetFieldName (da, 0, "2-D Potential Green's Function");
227  |       for (j=ystart; j<ystart+ywidth; j++)
228  | 	for (i=xstart; i<xstart+xwidth; i++)
229  | 	  {
230  | 	    double x,y;
231  | 
232  | 	    x = xmin + (xmax-xmin) * (double) i/(xglobal-1);
233  | 	    y = ymin + (ymax-ymin) * (double) j/(yglobal-1);
234  | 
235  | 	    array2d [j][i] = function_2d (x,y);
236  | 	  }
237  |       ierr = VecRestoreArray2d (vector, ywidth, xwidth, ystart, xstart,
238  | 				&array2d); CHKERRQ (ierr);
239  |     }
240  | 
241  |   /*+It then uses
242  |     +latex+{\tt GeomviewBegin()} or {\tt PetscViewerDrawOpen()}
243  |     +html+ <tt>GeomviewBegin()</tt> or <tt>PetscViewerDrawOpen()</tt>
244  |     to start the viewer, and either
245  |     +latex+{\tt DATriangulate()} and {\tt GeomviewDisplayTriangulation()} or
246  |     +latex+{\tt VecView()}
247  |     +html+ <tt>DATriangulate()</tt> and <tt>GeomviewDisplayTriangulation()</tt>
248  |     +html+ or <tt>VecView()</tt>
249  |     to display the solution. +*/
250  | 
251  |   DPRINTF ("About to open display\n", 0);
252  | 
253  |   if (threedee)
254  |     {
255  |       ierr = SurfCreate (&Surf); CHKERRQ (ierr);
256  |       ierr = GeomviewBegin (PETSC_COMM_WORLD, &Disp); CHKERRQ (ierr);
257  |     }
258  |   else
259  |     {
260  |       PetscDraw draw;
261  |       ierr = PetscViewerDrawOpen (PETSC_COMM_WORLD, 0, "", PETSC_DECIDE,
262  | 				  PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE,
263  | 				  &theviewer); CHKERRQ (ierr);
264  |       ierr = PetscViewerDrawGetDraw (theviewer, 0, &draw);
265  |       CHKERRQ (ierr);
266  |       ierr = PetscDrawSetDoubleBuffer (draw); CHKERRQ (ierr);
267  |     }
268  | 
269  |   DPRINTF ("About to do the graphing\n", 0);
270  | 
271  |   if (threedee)
272  |     {
273  |       PetscScalar minmax[6] = { xmin,xmax, ymin,ymax, zmin,zmax };
274  |       PetscScalar isovals[4] = { -.1, -.2, -.3, -.4 };
275  |       PetscScalar colors[16] = { 1,0,0,.4, 1,1,0,.4, 0,1,0,.4, 0,0,1,.4 };
276  | 
277  |       DPRINTF ("Starting triangulation\n", 0);
278  |       ierr = DATriangulate (Surf, da, vector, 0, minmax, 4, isovals, colors,
279  | 			    PETSC_FALSE, PETSC_FALSE, PETSC_FALSE);
280  |       CHKERRQ (ierr);
281  | 
282  |       DPRINTF ("Sending to Geomview\n", 0);
283  |       ierr = GeomviewDisplayTriangulation
284  | 	(PETSC_COMM_WORLD, Surf, Disp, minmax, "Function-Contours",PETSC_TRUE);
285  |       CHKERRQ (ierr);
286  | 
287  |       ierr = SurfClear (Surf); CHKERRQ (ierr);
288  |     }
289  |   else
290  |     {
291  |       ierr = VecView (vector, theviewer); CHKERRQ (ierr);
292  |     }
293  | 
294  |   /*+ Finally, it prompts the user to hit <return> before wrapping up. +*/
295  | 
296  |   {
297  |     char instring[100];
298  | 
299  |     ierr = PetscPrintf (PETSC_COMM_WORLD, "Press <return> to close up... ");
300  |     CHKERRQ (ierr);
301  |     ierr = PetscSynchronizedFGets (PETSC_COMM_WORLD, stdin, 99, instring);
302  |     CHKERRQ (ierr);
303  |   }
304  | 
305  |   DPRINTF ("Cleaning up\n", 0);
306  |   if (threedee)
307  |     {
308  |       ierr = GeomviewEnd (PETSC_COMM_WORLD, Disp); CHKERRQ (ierr);
309  |       ierr = ISurfDestroy (Surf); CHKERRQ (ierr);
310  |     }
311  |   else
312  |     {
313  |       ierr = PetscViewerDestroy (theviewer); CHKERRQ (ierr);
314  |     }
315  |   ierr = DARestoreGlobalVector (da, &vector); CHKERRQ (ierr);
316  |   ierr = DADestroy (da); CHKERRQ (ierr);
317  | 
318  |   PetscFinalize ();
319  |   return 0;
320  | }