1    | /***************************************
2    |   $Header$
3    | 
4    |   This is the petsc.c main file.  It has all of the PETSc-dependent functions.
5    | ***************************************/
6    | 
7    | 
8    | #include <stdlib.h>      /* For exit() */
9    | #include <petscda.h>
10   | #include "config.h"      /* esp. for inline */
11   | #include "illuminator.h" /* Just to make sure the interface is "right" */
12   | #include "structures.h"     /* For the isurface structure definition */
13   | 
14   | 
15   | #undef __FUNCT__
16   | #define __FUNCT__ "DATriangulateRange"
17   | 
18   | /*++++++++++++++++++++++++++++++++++++++
19   |   Calculate vertices of isoquant triangles in a 3-D distributed array.  This
20   |   takes a PETSc DA object, does some sanity checks, calculates array sizes,
21   |   gets the local vector and array, and then calls
22   |   +latex+{\tt DATriangulateLocal()}
23   |   +html+ <tt>DATriangulateLocal()</tt>
24   |   to do the rest.  Note that global array access (i.e. this function) is
25   |   necessary for using default isoquant values, since we need to be able to
26   |   calculate the maximum and minimum on the global array.
27   | 
28   |   int DATriangulateRange Returns 0 or an error code.
29   | 
30   |   ISurface Surf ISurface object into which to draw this DA's isoquants.
31   | 
32   |   DA theda The PETSc distributed array object.
33   | 
34   |   Vec globalX PETSc global vector object associated with the DA with data we'd
35   |   like to graph.
36   | 
37   |   int this Index of the field we'd like to draw.
38   | 
39   |   PetscScalar *minmax Position of block corners: xmin, xmax, ymin, ymax, zmin,
40   |   zmax.
41   | 
42   |   int n_quants Number of isoquant surfaces to draw (isoquant values), or
43   |   +latex+{\tt PETSC\_DECIDE}
44   |   +html+ <tt>PETSC_DECIDE</tt>
45   |   to use red, yellow, green, blue at 0.2, 0.4, 0.6 and 0.8 between the vector's
46   |   global minimum and maximum values.
47   | 
48   |   PetscScalar *isoquants Array of function values at which to draw isoquants,
49   |   +latex+or {\tt PETSC\_NULL} if {\tt n\_quants=PETSC\_DECIDE}.
50   |   +html+ or <tt>PETSC_NULL</tt> if <tt>n_quants=PETSC\_DECIDE</tt>.
51   | 
52   |   PetscScalar *colors Array of color R,G,B,A quads for each isoquant, or
53   |   +latex+{\tt PETSC\_NULL} if {\tt n\_quants=PETSC\_DECIDE}.
54   |   +html+ <tt>PETSC_NULL</tt> if <tt>n_quants=PETSC\_DECIDE</tt>.
55   | 
56   |   int xmin Smallest grid x-coordinate to render.
57   | 
58   |   int xmax Largest grid x-coordinate to render, -1 goes to full x maximum, -2
59   |   in periodic systems goes to one short of x maximum.
60   | 
61   |   int ymin Smallest grid y-coordinate to render.
62   | 
63   |   int ymax Largest grid y-coordinate to render, -1 goes to full y maximum, -2
64   |   in periodic systems goes to one short of y maximum.
65   | 
66   |   int zmin Smallest grid z-coordinate to render.
67   | 
68   |   int zmax Largest grid z-coordinate to render, -1 goes to full z maximum, -2
69   |   in periodic systems goes to one short of z maximum.
70   |   ++++++++++++++++++++++++++++++++++++++*/
71   | 
72   | int DATriangulateRange
73   | (ISurface Surf, DA theda, Vec globalX, int this, PetscScalar *minmax,
74   |  int n_quants, PetscScalar *isoquants, PetscScalar *colors, int xmin,int xmax,
75   |  int ymin,int ymax, int zmin,int zmax)
76   | {
77   |   Vec localX;
78   |   PetscScalar isoquantdefaults[4],
79   |     colordefaults[16] = { 1,0,0,.5, 1,1,0,.5, 0,1,0,.5, 0,0,1,.5 };
80   |   PetscReal themax, themin;
81   |   int ierr;
82   | 
83   |   /* Default isoquants */
84   |   if (n_quants == PETSC_DECIDE) {
85   |     ierr = VecStrideMin (globalX, this, PETSC_NULL, &themin); CHKERRQ (ierr);
86   |     ierr = VecStrideMax (globalX, this, PETSC_NULL, &themax); CHKERRQ (ierr);
87   |     /* Red, yellow, green, blue at 0.2, 0.4, 0.6, 0.8, all with alpha=0.5 */
88   |     n_quants = 4;
89   |     isoquantdefaults[0] = themin + 0.2*(themax-themin);
90   |     isoquantdefaults[1] = themin + 0.4*(themax-themin);
91   |     isoquantdefaults[2] = themin + 0.6*(themax-themin);
92   |     isoquantdefaults[3] = themin + 0.8*(themax-themin);
93   |     isoquants = isoquantdefaults;
94   |     colors = colordefaults;
95   |   }
96   | 
97   |   /* Get the local array of points, with ghosts */
98   |   ierr = DACreateLocalVector (theda, &localX); CHKERRQ (ierr);
99   |   ierr = DAGlobalToLocalBegin (theda, globalX, INSERT_VALUES, localX); CHKERRQ(ierr);
100  |   ierr = DAGlobalToLocalEnd (theda, globalX, INSERT_VALUES, localX); CHKERRQ (ierr);
101  | 
102  |   /* Use DATriangulateLocalRange() to do the work */
103  |   ierr = DATriangulateLocalRange (Surf, theda, localX, this, minmax, n_quants,
104  | 				  isoquants, colors, xmin,xmax, ymin,ymax,
105  | 				  zmin,zmax); CHKERRQ (ierr);
106  | 
107  |   ierr = VecDestroy (localX); CHKERRQ (ierr);
108  | 
109  |   return 0;
110  | }
111  | 
112  | 
113  | #undef __FUNCT__
114  | #define __FUNCT__ "DATriangulateLocalRange"
115  | 
116  | /*++++++++++++++++++++++++++++++++++++++
117  |   Calculate vertices of isoquant triangles in a 3-D distributed array.  This
118  |   takes a PETSc DA object, does some sanity checks, calculates array sizes, and
119  |   then gets array and passes it to Draw3DBlock for triangulation.
120  | 
121  |   int DATriangulateLocalRange Returns 0 or an error code.
122  | 
123  |   ISurface Surf ISurface object into which to draw this DA's isoquants.
124  | 
125  |   DA theda The PETSc distributed array object.
126  | 
127  |   Vec localX PETSc local vector object associated with the DA with data we'd
128  |   like to graph.
129  | 
130  |   int this Index of the field we'd like to draw.
131  | 
132  |   PetscScalar *minmax Position of block corners: xmin, xmax, ymin, ymax, zmin,
133  |   zmax.
134  | 
135  |   int n_quants Number of isoquant surfaces to draw (isoquant values).  Note
136  |   +latex+{\tt PETSC\_DECIDE}
137  |   +html+ <tt>PETSC_DECIDE</tt>
138  |   is not a valid option here, because it's impossible to know the global
139  |   maximum and minimum and have consistent contours without user-supplied
140  |   information.
141  | 
142  |   PetscScalar *isoquants Array of function values at which to draw isoquants.
143  | 
144  |   PetscScalar *colors Array of color R,G,B,A quads for each isoquant.
145  | 
146  |   int xmin Smallest grid x-coordinate to render.
147  | 
148  |   int xmax Largest grid x-coordinate to render, -1 goes to full x maximum, -2
149  |   in periodic systems goes to one short of x maximum.
150  | 
151  |   int ymin Smallest grid y-coordinate to render.
152  | 
153  |   int ymax Largest grid y-coordinate to render, -1 goes to full y maximum, -2
154  |   in periodic systems goes to one short of y maximum.
155  | 
156  |   int zmin Smallest grid z-coordinate to render.
157  | 
158  |   int zmax Largest grid z-coordinate to render, -1 goes to full z maximum, -2
159  |   in periodic systems goes to one short of z maximum.
160  |   ++++++++++++++++++++++++++++++++++++++*/
161  | 
162  | int DATriangulateLocalRange
163  | (ISurface Surf, DA theda, Vec localX, int this, PetscScalar *minmax,
164  |  int n_quants, PetscScalar *isoquants, PetscScalar *colors, int xmin,int xmax,
165  |  int ymin,int ymax, int zmin,int zmax)
166  | {
167  |   PetscScalar *x, isoquantdefaults[4], localminmax[6],
168  |     colordefaults[16] = { 1,0,0,.5, 1,1,0,.5, 0,1,0,.5, 0,0,1,.5 };
169  |   DAStencilType stencil;
170  |   int i, ierr, fields, xw,yw,zw, xs,ys,zs, xm,ym,zm, gxs,gys,gzs, gxm,gym,gzm;
171  |   struct isurface *thesurf = (struct isurface *) Surf;
172  | 
173  |   /* Default isoquant error */
174  |   if (n_quants == PETSC_DECIDE)
175  |     SETERRQ (PETSC_ERR_ARG_WRONGSTATE, "Can't get default isoquants for local array");
176  | 
177  |   /* Get global and local grid boundaries */
178  |   ierr = DAGetInfo (theda, &i, &xw,&yw,&zw, PETSC_NULL,PETSC_NULL,PETSC_NULL,
179  | 		    &fields, PETSC_NULL, PETSC_NULL, &stencil);CHKERRQ(ierr);
180  |   if (i!=3)
181  |     SETERRQ (PETSC_ERR_ARG_WRONGSTATE, "DA must be 3-dimensional");
182  |   if (stencil!=DA_STENCIL_BOX)
183  |     SETERRQ (PETSC_ERR_ARG_WRONGSTATE, "DA must have a box stencil");
184  | 
185  |   ierr = DAGetCorners (theda, &xs,&ys,&zs, &xm,&ym,&zm); CHKERRQ (ierr);
186  |   ierr = DAGetGhostCorners (theda, &gxs,&gys,&gzs, &gxm,&gym,&gzm);
187  |   CHKERRQ (ierr);
188  | 
189  |   /* Get the local array of points, with ghosts */
190  |   ierr = VecGetArray (localX, &x); CHKERRQ (ierr);
191  | 
192  |   /* If the array is too big, cut it down to size */
193  |   if (gxm <= xs-gxs+xm)
194  |     xm = gxm-xs+gxs-1;
195  |   if (gym <= ys-gys+ym)
196  |     ym = gym-ys+gys-1;
197  |   if (gzm <= zs-gzs+zm)
198  |     zm = gzm-zs+gzs-1;
199  |   /* Eliminate final rows/planes if *cut and periodic. */
200  |   if (xmax == -2 && xs+xm>=xw)
201  |     xm--;
202  |   if (ymax == -2 && ys+ym>=yw)
203  |     ym--;
204  |   if (zmax == -2 && zs+zm>=zw)
205  |     zm--;
206  | 
207  |   /* Reframe to range */
208  |   xs = PetscMax (xs, xmin);
209  |   ys = PetscMax (ys, xmin);
210  |   zs = PetscMax (zs, xmin);
211  |   xm = (xmax > 0) ? PetscMin (xm, xmax-xs) : xm;
212  |   ym = (ymax > 0) ? PetscMin (ym, ymax-ys) : ym;
213  |   zm = (zmax > 0) ? PetscMin (zm, zmax-zs) : zm;
214  | 
215  |   /* Calculate local physical size */
216  |   localminmax[0] = minmax[0] + (minmax[1]-minmax[0])*xs/xw;
217  |   localminmax[1] = minmax[2] + (minmax[1]-minmax[0])*(xs+xm)/xw;
218  |   localminmax[2] = minmax[2] + (minmax[3]-minmax[2])*ys/yw;
219  |   localminmax[3] = minmax[2] + (minmax[3]-minmax[2])*(ys+ym)/yw;
220  |   localminmax[4] = minmax[4] + (minmax[5]-minmax[4])*zs/zw;
221  |   localminmax[5] = minmax[4] + (minmax[5]-minmax[4])*(zs+zm)/zw;
222  | 
223  |   /* Let 'er rip! */
224  |   ierr = Draw3DBlock (Surf, gxm,gym,gzm, xs-gxs,ys-gys,zs-gzs, xm,ym,zm,
225  | 		      localminmax, x+this, fields, n_quants,isoquants, colors);
226  |   CHKERRQ (ierr);
227  |   ierr = VecRestoreArray (localX, &x); CHKERRQ (ierr);
228  | 
229  |   /* Prints the number of triangles on all CPUs */
230  | #ifdef DEBUG
231  |   {
232  |     int rank;
233  |     ierr = MPI_Comm_rank (PETSC_COMM_WORLD, &rank); CHKERRQ (ierr);
234  |     ierr = PetscSynchronizedPrintf
235  |       (PETSC_COMM_WORLD, "[%d] zs=%d, zm=%d, zmin=%g, zmax=%g\n",
236  |        rank, zs, zm, localminmax[4], localminmax[5]); CHKERRQ (ierr);
237  |     ierr = PetscSynchronizedFlush (PETSC_COMM_WORLD); CHKERRQ (ierr);
238  |     ierr = PetscSynchronizedPrintf (PETSC_COMM_WORLD, "[%d] Triangles: %d\n",
239  | 				    rank, thesurf->num_triangles);
240  |     CHKERRQ (ierr);
241  |     ierr = PetscSynchronizedFlush (PETSC_COMM_WORLD); CHKERRQ (ierr);
242  |   }
243  | #endif
244  | 
245  |   return 0;
246  | }
247  | 
248  | 
249  | #undef __FUNCT__
250  | #define __FUNCT__ "IllErrorHandler"
251  | 
252  | /*++++++++++++++++++++++++++++++++++++++
253  |   Handle errors, in this case the PETSc way.
254  | 
255  |   int IllErrorHandler Returns the error code supplied.
256  | 
257  |   int id Index of the error, defined in petscerror.h.
258  | 
259  |   char *message Text of the error message.
260  |   ++++++++++++++++++++++++++++++++++++++*/
261  | 
262  | int IllErrorHandler (int id, char *message)
263  | {
264  |   int ierr;
265  |   ierr = PetscPrintf (PETSC_COMM_WORLD, "Warning: IllErrorHandler is broken and deprecated, use PETSc SETERRQ instead.\n");
266  |   SETERRQ (id, message);
267  |   exit (1);
268  | }