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