1    | /***************************************
2    |   $Header$
3    | 
4    |   This file has the Geomview interface, including the PETSc vector gather
5    |   operations to bring everything to CPU 0.
6    | ***************************************/
7    | 
8    | 
9    | #include <stdio.h>
10   | #include <petscvec.h>
11   | #include "config.h"      /* esp. for inline */
12   | #include "illuminator.h" /* Just to make sure the interface is "right" */
13   | #include "structures.h"     /* For isurface structure definition */
14   | #include <stdlib.h>      /* For malloc() */
15   | #include <unistd.h>      /* For execlp() */
16   | #include <glib.h>        /* For guint*, GUINT*_SWAP_LE_BE */
17   | 
18   | /* Build with -DDEBUG for debugging output */
19   | #undef DPRINTF
20   | #ifdef DEBUG
21   | #define DPRINTF(fmt, args...) PetscPrintf (PETSC_COMM_WORLD, "%s: " fmt, __FUNCT__, args)
22   | #else
23   | #define DPRINTF(fmt, args...)
24   | #endif
25   | 
26   | 
27   | #undef __FUNCT__
28   | #define __FUNCT__ "GeomviewBegin"
29   | 
30   | /*++++++++++++++++++++++++++++++++++++++
31   |   Spawn a new geomview process.  Most of this was shamelessly ripped from Ken
32   |   Brakke's Surface Evolver.
33   | 
34   |   int GeomviewBegin Returns 0 or an error code.
35   | 
36   |   MPI_Comm comm MPI communicator for rank information, if NULL it uses
37   |   PETSC_COMM_WORLD.
38   | 
39   |   IDisplay *newdisp Pointer in which to put a new IDisplay object.
40   |   ++++++++++++++++++++++++++++++++++++++*/
41   | 
42   | int GeomviewBegin(MPI_Comm comm, IDisplay *newdisp)
43   | {
44   |   int rank, ierr, gv_pid,gv_pipe[2], to_gv_pipe[2];
45   |   char gv_version[100];
46   |   struct idisplay *thedisp;
47   | 
48   |   if (!(thedisp = (struct idisplay *) malloc (sizeof (struct idisplay))))
49   |     SETERRQ (PETSC_ERR_MEM, "Unable to allocate memory for isurface object");
50   | 
51   |   /* Initialize display image to zero */
52   |   thedisp->rgb = NULL;
53   |   thedisp->rgb_width = thedisp->rgb_height = thedisp->rgb_rowskip =
54   |     thedisp->rgb_bpp = 0;
55   |   thedisp->zbuf = NULL;
56   |   thedisp->zbuf_width = thedisp->zbuf_height = thedisp->zbuf_depth =
57   |     thedisp->zbuf_rowskip = 0;
58   | 
59   |   if (!comm)
60   |     comm = PETSC_COMM_WORLD;
61   |   ierr = MPI_Comm_rank (comm,&rank); CHKERRQ(ierr);
62   |   if (!rank)
63   |     {
64   |       pipe (gv_pipe); /* from geomview stdout */
65   |       pipe (to_gv_pipe); /* to geomview stdin */
66   |       gv_pid = fork ();
67   | 
68   |       if(gv_pid==0) /* child */
69   | 	{
70   | 	  close (0);
71   | 	  dup (to_gv_pipe[0]);
72   | 	  close (to_gv_pipe[0]);
73   | 	  close (to_gv_pipe[1]);
74   | 	  close (1);
75   | 	  dup (gv_pipe[1]);
76   | 	  close (gv_pipe[0]);
77   | 	  close (gv_pipe[1]);
78   | 	  /* signal(SIGINT,SIG_IGN); */
79   | 	  execlp (GEOMVIEW,GEOMVIEW,"-c","(interest (pick world))","-",NULL);
80   | 	  perror (GEOMVIEW); /* only error gets here */
81   | 	  SETERRQ (PETSC_ERR_ARG_WRONGSTATE,"Geomview invocation failed.\n");
82   | 	}
83   | 
84   |       /* PETSc program execution resumes here */
85   |       close (gv_pipe[1]);
86   |       close (to_gv_pipe[0]);
87   |       thedisp->to_geomview = fdopen (to_gv_pipe[1], "w"); /* geomview stdin */
88   |       fprintf (thedisp->to_geomview, "(echo (geomview-version) \"\n\")\n");
89   |       fflush (thedisp->to_geomview);
90   |       read (gv_pipe[0], gv_version, sizeof (gv_version));
91   |     }
92   | 
93   |   *newdisp = (IDisplay) thedisp;
94   |   return 0;
95   | }
96   | 
97   | 
98   | #undef __FUNCT__
99   | #define __FUNCT__ "GeomviewEnd"
100  | 
101  | /*++++++++++++++++++++++++++++++++++++++
102  |   Exit the current running Geomview process and close its pipe.  Based in part
103  |   on Ken Brakke's Surface Evolver.
104  | 
105  |   int GeomviewEnd Returns 0 or an error code.
106  | 
107  |   MPI_Comm comm MPI communicator for rank information, if NULL it uses
108  |   PETSC_COMM_WORLD.
109  | 
110  |   IDisplay Disp IDisplay object whose geomview FILE we're closing.
111  |   ++++++++++++++++++++++++++++++++++++++*/
112  | 
113  | int GeomviewEnd (MPI_Comm comm, IDisplay Disp)
114  | {
115  |   int rank, ierr;
116  |   struct idisplay *thedisp = (struct idisplay *) Disp;
117  | 
118  |   if (!comm)
119  |     comm = PETSC_COMM_WORLD;
120  |   ierr = MPI_Comm_rank (comm, &rank); CHKERRQ (ierr);
121  |   if (!rank)
122  |     {
123  |       if (thedisp->to_geomview)
124  | 	{
125  | 	  fprintf (thedisp->to_geomview, "(exit)");
126  | 	  fclose (thedisp->to_geomview);
127  | 	  thedisp->to_geomview = NULL;
128  | 	}
129  |     }
130  |   return 0;
131  | }
132  | 
133  | 
134  | #ifdef WORDS_BIGENDIAN
135  | # define GEOMVIEW_SET_INT(j,x) gvint[j]=(guint32)(x)
136  | # define GEOMVIEW_SET_FLOAT(j,x) gvfloat[j]=(float)(x)
137  | #else
138  | # define GEOMVIEW_SET_INT(j,x) \
139  |   gvint[j]=GUINT32_SWAP_LE_BE_CONSTANT((guint32) (x))
140  | # define GEOMVIEW_SET_FLOAT(j,x) \
141  |   {gvfloat[j]=(float)(x); gvint[j]=GUINT32_SWAP_LE_BE_CONSTANT(gvint[j]);}
142  | #endif
143  | 
144  | /* #define GEOMVIEW_BINARY */
145  | 
146  | 
147  | #undef __FUNCT__
148  | #define __FUNCT__ "GeomviewDisplayTriangulation"
149  | 
150  | /*++++++++++++++++++++++++++++++++++++++
151  |   Pipe the current triangulation to Geomview for display.  Much of this is
152  |   based on Ken Brakke's Surface Evolver.
153  | 
154  |   int GeomviewDisplayTriangulation Returns 0 or an error code.
155  | 
156  |   MPI_Comm comm MPI communicator for rank information, if NULL it uses
157  |   PETSC_COMM_WORLD.
158  | 
159  |   ISurface Surf ISurface object containing triangles to render using Geomview.
160  | 
161  |   IDisplay Disp IDisplay object whose geomview FILE we're closing.
162  | 
163  |   PetscScalar *minmax Position of block corners: xmin, xmax, ymin, ymax, zmin, zmax.
164  | 
165  |   char *name Name to give the Geomview OOGL object which we create.
166  | 
167  |   PetscTruth transparent Geomview transparency flag.
168  |   ++++++++++++++++++++++++++++++++++++++*/
169  | 
170  | int GeomviewDisplayTriangulation
171  | (MPI_Comm comm, ISurface Surf, IDisplay Disp, PetscScalar *minmax, char *name,
172  |  PetscTruth transparent)
173  | {
174  |   int ierr, i, total_entries, start, end, rank;
175  |   PetscScalar *localvals;
176  |   Vec globalverts, localverts;
177  |   VecScatter vecscat;
178  |   IS islocal;
179  |   struct isurface *thesurf = (struct isurface *) Surf;
180  |   struct idisplay *thedisp = (struct idisplay *) Disp;
181  | 
182  |   /* Simple "assertion" check */
183  |   if (!comm)
184  |     comm = PETSC_COMM_WORLD;
185  |   ierr = MPI_Comm_rank (comm, &rank); CHKERRQ (ierr);
186  |   if (!rank && !thedisp->to_geomview)
187  |     SETERRQ (PETSC_ERR_ARG_WRONGSTATE,"Geomview display is not open!");
188  | 
189  |   /*+ First, this creates global and local vectors for all of the triangle
190  |     vertices. +*/
191  |   ierr = VecCreateMPIWithArray (comm, 13*thesurf->num_triangles,
192  | 				PETSC_DETERMINE,
193  | 				thesurf->vertices, &globalverts); CHKERRQ (ierr);
194  |   ierr = VecGetSize (globalverts, &total_entries); CHKERRQ (ierr);
195  |   ierr = VecCreateMPI (comm, (rank == 0) ? total_entries:0, total_entries,
196  | 		       &localverts); CHKERRQ (ierr);
197  |   DPRINTF ("Total triangles: %d\n", total_entries/13);
198  | 
199  |   /* Diagnostics which may be useful to somebody sometime */
200  |   /* ierr = VecGetOwnershipRange (globalverts, &start, &end); CHKERRQ (ierr);
201  |   ierr = PetscSynchronizedPrintf
202  |     (comm, "[%d] Global vector local size: %d, ownership from %d to %d\n",
203  |      rank, end-start, start, end); CHKERRQ (ierr);
204  |   ierr = PetscSynchronizedFlush (comm); CHKERRQ (ierr);
205  | 
206  |   ierr = VecGetOwnershipRange (localverts, &start, &end); CHKERRQ (ierr);
207  |   ierr = PetscSynchronizedPrintf
208  |     (comm, "[%d] Local vector local size: %d, ownership from %d to %d\n", rank,
209  |      end-start, start, end); CHKERRQ (ierr);
210  |   ierr = PetscSynchronizedFlush (comm); CHKERRQ (ierr);
211  |   ierr = PetscPrintf (comm, "Global vector size: %d\n", total_entries);
212  |   CHKERRQ (ierr); */
213  | 
214  |   /*+ It then gathers (``scatters'') all vertex data to processor zero, +*/
215  |   ierr = ISCreateStride (comm, total_entries, 0, 1, &islocal); CHKERRQ (ierr);
216  |   ierr = VecScatterCreate (globalverts, PETSC_NULL, localverts, islocal,
217  | 			   &vecscat); CHKERRQ (ierr);
218  |   DPRINTF ("Starting triangle scatter to head node\n",0);
219  |   ierr = VecScatterBegin (vecscat, globalverts, localverts, INSERT_VALUES,
220  | 			  SCATTER_FORWARD); CHKERRQ (ierr);
221  |   DPRINTF ("Triangle scatter to head node in progress...\n",0);
222  |   ierr = VecScatterEnd (vecscat, globalverts, localverts, INSERT_VALUES,
223  | 			SCATTER_FORWARD); CHKERRQ (ierr);
224  |   DPRINTF ("Triangle scatter to head node complete\n",0);
225  |   ierr = VecScatterDestroy (vecscat); CHKERRQ (ierr);
226  |   ierr = ISDestroy (islocal); CHKERRQ (ierr);
227  | 
228  |   /*+ and puts them in an array. +*/
229  |   ierr = VecGetArray (localverts, &localvals); CHKERRQ (ierr);
230  | 
231  |   /*+ Finally, it sends everything to Geomview, +*/
232  |   if (!rank) {
233  | #ifdef GEOMVIEW_BINARY
234  |     guint32 gvint[10];
235  |     float *gvfloat = (float *)gvint;
236  | #endif
237  | 
238  |     fprintf (thedisp->to_geomview, "(geometry \"%s\" { : bor })", name);
239  |     fprintf (thedisp->to_geomview, "(read geometry { define bor \n");
240  |     fprintf (thedisp->to_geomview, "appearance {%s}\nOFF%s\n",
241  | 	     transparent ? "+transparent" : "",
242  | #ifdef GEOMVIEW_BINARY
243  | 	     " BINARY"
244  | #else
245  | 	     ""
246  | #endif
247  | 	     );
248  | #ifndef GEOMVIEW_BINARY
249  |     fprintf (thedisp->to_geomview, "%d %d 0\n", 3*(total_entries/13) + 2,
250  | 	     total_entries/13);
251  | #else
252  |     GEOMVIEW_SET_INT (0, 3*(total_entries/13) + 2);
253  |     GEOMVIEW_SET_INT (1, total_entries/13);
254  |     GEOMVIEW_SET_INT (2, 0);
255  |     fwrite (gvint, sizeof (guint32), 3, thedisp->to_geomview);
256  | #endif
257  |     for (i=0; i<total_entries; i+=13)
258  |       {
259  | #ifndef GEOMVIEW_BINARY
260  | 	fprintf (thedisp->to_geomview, "%g %g %g\n%g %g %g\n%g %g %g\n",
261  | 		 localvals[i], localvals[i+1], localvals[i+2], localvals[i+3],
262  | 		 localvals[i+4], localvals[i+5], localvals[i+6],localvals[i+7],
263  | 		 localvals[i+8]);
264  | #else
265  | 	GEOMVIEW_SET_FLOAT(0, localvals[i]);
266  | 	GEOMVIEW_SET_FLOAT(1, localvals[i+1]);
267  | 	GEOMVIEW_SET_FLOAT(2, localvals[i+2]);
268  | 	GEOMVIEW_SET_FLOAT(3, localvals[i+3]);
269  | 	GEOMVIEW_SET_FLOAT(4, localvals[i+4]);
270  | 	GEOMVIEW_SET_FLOAT(5, localvals[i+5]);
271  | 	GEOMVIEW_SET_FLOAT(6, localvals[i+6]);
272  | 	GEOMVIEW_SET_FLOAT(7, localvals[i+7]);
273  | 	GEOMVIEW_SET_FLOAT(8, localvals[i+8]);
274  | 	fwrite (gvfloat, sizeof (float), 9, thedisp->to_geomview);
275  | #endif
276  |       }
277  | #ifndef GEOMVIEW_BINARY
278  |     fprintf (thedisp->to_geomview, "%g %g %g\n%g %g %g\n", minmax[0],minmax[2],
279  | 	     minmax[4], minmax[1], minmax[3], minmax[5]);
280  | #else
281  |     GEOMVIEW_SET_FLOAT(0, minmax[0]);
282  |     GEOMVIEW_SET_FLOAT(1, minmax[2]);
283  |     GEOMVIEW_SET_FLOAT(2, minmax[4]);
284  |     GEOMVIEW_SET_FLOAT(3, minmax[1]);
285  |     GEOMVIEW_SET_FLOAT(4, minmax[3]);
286  |     GEOMVIEW_SET_FLOAT(5, minmax[5]);
287  |     fwrite (gvfloat, sizeof (float), 6, thedisp->to_geomview);
288  | #endif
289  |     for (i=0; i<total_entries/13; i++)
290  |       {
291  | #ifndef GEOMVIEW_BINARY
292  |         fprintf (thedisp->to_geomview,"3 %d %d %d %g %g %g %g\n",
293  | 		 3*i, 3*i+1, 3*i+2,
294  | 		 localvals[13*i+9], localvals[13*i+10], localvals[13*i+11],
295  | 		 localvals[13*i+12]);
296  | #else
297  | 	GEOMVIEW_SET_INT(0, 3);
298  | 	GEOMVIEW_SET_INT(1, 3*i);
299  | 	GEOMVIEW_SET_INT(2, 3*i+1);
300  | 	GEOMVIEW_SET_INT(3, 3*i+2);
301  | 	GEOMVIEW_SET_INT(4, 4);
302  | 	fwrite (gvint, sizeof (guint32), 5, thedisp->to_geomview);
303  | 	GEOMVIEW_SET_FLOAT(0, localvals[13*i+9]);
304  | 	GEOMVIEW_SET_FLOAT(1, localvals[13*i+10]);
305  | 	GEOMVIEW_SET_FLOAT(2, localvals[13*i+11]);
306  | 	GEOMVIEW_SET_FLOAT(3, localvals[13*i+12]);
307  | 	fwrite (gvfloat, sizeof (float), 4, thedisp->to_geomview);
308  | #endif
309  |       }
310  |     fprintf (thedisp->to_geomview, "})\n");
311  |     fflush (thedisp->to_geomview);
312  |   }
313  | 
314  |   /*+ and cleans up the mess. +*/
315  |   ierr = VecRestoreArray (localverts, &localvals); CHKERRQ (ierr);
316  |   ierr = VecDestroy (localverts); CHKERRQ (ierr);
317  |   ierr = VecDestroy (globalverts); CHKERRQ (ierr);
318  | 
319  |   return 0;
320  | }