1 | /*************************************** 2 | $Header: /cvsroot/petscgraphics/3dgf.c,v 1.12 2003/04/29 13:38:37 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 | CHKERRQ (ierr); 277 | DPRINTF ("Sending to Geomview\n", 0); 278 | ierr = GeomviewDisplayTriangulation 279 | (PETSC_COMM_WORLD, minmax, "Function-Contours", PETSC_TRUE); 280 | CHKERRQ (ierr); 281 | } 282 | else 283 | { 284 | ierr = VecView (vector, theviewer); CHKERRQ (ierr); 285 | } 286 | 287 | /*+ Finally, it prompts the user to hit <return> before wrapping up. +*/ 288 | 289 | { 290 | char instring[100]; 291 | 292 | ierr = PetscPrintf (PETSC_COMM_WORLD, "Press <return> to close up... "); 293 | CHKERRQ (ierr); 294 | ierr = PetscSynchronizedFGets (PETSC_COMM_WORLD, stdin, 99, instring); 295 | CHKERRQ (ierr); 296 | } 297 | 298 | DPRINTF ("Cleaning up\n", 0); 299 | if (threedee) 300 | { 301 | ierr = GeomviewEnd (PETSC_COMM_WORLD); CHKERRQ (ierr); 302 | } 303 | else 304 | { 305 | ierr = PetscViewerDestroy (theviewer); CHKERRQ (ierr); 306 | } 307 | ierr = DARestoreGlobalVector (da, &vector); CHKERRQ (ierr); 308 | ierr = DADestroy (da); CHKERRQ (ierr); 309 | 310 | PetscFinalize (); 311 | return 0; 312 | }