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