APBS  1.4.1
vgrid.c
Go to the documentation of this file.
1 
57 #include "vgrid.h"
58 
59 VEMBED(rcsid="$Id$")
60 
61 #if !defined(VINLINE_VGRID)
62  VPUBLIC unsigned long int Vgrid_memChk(Vgrid *thee) {
63  return Vmem_bytes(thee->mem);
64  }
65 #endif
66 #define IJK(i,j,k) (((k)*(nx)*(ny))+((j)*(nx))+(i))
67 
68 VPRIVATE char *MCwhiteChars = " =,;\t\n";
69 VPRIVATE char *MCcommChars = "#%";
70 VPRIVATE double Vcompare;
71 VPRIVATE char Vprecision[26];
72 
73 /* ///////////////////////////////////////////////////////////////////////////
74 // Routine: Vgrid_ctor
75 // Author: Nathan Baker
77 VPUBLIC Vgrid* Vgrid_ctor(int nx,
78  int ny,
79  int nz,
80  double hx,
81  double hy,
82  double hzed,
83  double xmin,
84  double ymin,
85  double zmin,
86  double *data
87  ) {
88 
89  Vgrid *thee = VNULL;
90 
91  thee = (Vgrid*)Vmem_malloc(VNULL, 1, sizeof(Vgrid));
92  VASSERT(thee != VNULL);
93  VASSERT(Vgrid_ctor2(thee, nx, ny, nz, hx, hy, hzed,
94  xmin, ymin, zmin, data));
95 
96  return thee;
97 }
98 
99 /* ///////////////////////////////////////////////////////////////////////////
100 // Routine: Vgrid_ctor2
101 // Author: Nathan Baker
103 VPUBLIC int Vgrid_ctor2(Vgrid *thee, int nx, int ny, int nz,
104  double hx, double hy, double hzed,
105  double xmin, double ymin, double zmin,
106  double *data) {
107 
108  if (thee == VNULL) return 0;
109  thee->nx = nx;
110  thee->ny = ny;
111  thee->nz = nz;
112  thee->hx = hx;
113  thee->hy = hy;
114  thee->hzed = hzed;
115  thee->xmin = xmin;
116  thee->xmax = xmin + (nx-1)*hx;
117  thee->ymin = ymin;
118  thee->ymax = ymin + (ny-1)*hy;
119  thee->zmin = zmin;
120  thee->zmax = zmin + (nz-1)*hzed;
121  if (data == VNULL) {
122  thee->ctordata = 0;
123  thee->readdata = 0;
124  } else {
125  thee->ctordata = 1;
126  thee->readdata = 0;
127  thee->data = data;
128  }
129 
130  thee->mem = Vmem_ctor("APBS:VGRID");
131 
132  Vcompare = pow(10,-1*(VGRID_DIGITS - 2));
133  sprintf(Vprecision,"%%12.%de %%12.%de %%12.%de", VGRID_DIGITS,
134  VGRID_DIGITS, VGRID_DIGITS);
135 
136  return 1;
137 }
138 
139 /* ///////////////////////////////////////////////////////////////////////////
140 // Routine: Vgrid_dtor
141 // Author: Nathan Baker
143 VPUBLIC void Vgrid_dtor(Vgrid **thee) {
144 
145  if ((*thee) != VNULL) {
146  Vgrid_dtor2(*thee);
147  Vmem_free(VNULL, 1, sizeof(Vgrid), (void **)thee);
148  (*thee) = VNULL;
149  }
150 }
151 
152 /* ///////////////////////////////////////////////////////////////////////////
153 // Routine: Vgrid_dtor2
154 // Author: Nathan Baker
156 VPUBLIC void Vgrid_dtor2(Vgrid *thee) {
157 
158  if (thee->readdata) {
159  Vmem_free(thee->mem, (thee->nx*thee->ny*thee->nz), sizeof(double),
160  (void **)&(thee->data));
161  }
162  Vmem_dtor(&(thee->mem));
163 
164 }
165 
166 /* ///////////////////////////////////////////////////////////////////////////
167 // Routine: Vgrid_value
168 // Author: Nathan Baker
170 VPUBLIC int Vgrid_value(Vgrid *thee, double pt[3], double *value) {
171 
172  int nx, ny, nz, ihi, jhi, khi, ilo, jlo, klo;
173  double hx, hy, hzed, xmin, ymin, zmin, ifloat, jfloat, kfloat;
174  double xmax, ymax, zmax;
175  double u, dx, dy, dz;
176 
177  if (thee == VNULL) {
178  Vnm_print(2, "Vgrid_value: Error -- got VNULL thee!\n");
179  VASSERT(0);
180  }
181  if (!(thee->ctordata || thee->readdata)) {
182  Vnm_print(2, "Vgrid_value: Error -- no data available!\n");
183  VASSERT(0);
184  }
185 
186  nx = thee->nx;
187  ny = thee->ny;
188  nz = thee->nz;
189  hx = thee->hx;
190  hy = thee->hy;
191  hzed = thee->hzed;
192  xmin = thee->xmin;
193  ymin = thee->ymin;
194  zmin = thee->zmin;
195  xmax = thee->xmax;
196  ymax = thee->ymax;
197  zmax = thee->zmax;
198 
199  u = 0;
200 
201  ifloat = (pt[0] - xmin)/hx;
202  jfloat = (pt[1] - ymin)/hy;
203  kfloat = (pt[2] - zmin)/hzed;
204 
205  ihi = (int)ceil(ifloat);
206  jhi = (int)ceil(jfloat);
207  khi = (int)ceil(kfloat);
208  ilo = (int)floor(ifloat);
209  jlo = (int)floor(jfloat);
210  klo = (int)floor(kfloat);
211  if (VABS(pt[0] - xmin) < Vcompare) ilo = 0;
212  if (VABS(pt[1] - ymin) < Vcompare) jlo = 0;
213  if (VABS(pt[2] - zmin) < Vcompare) klo = 0;
214  if (VABS(pt[0] - xmax) < Vcompare) ihi = nx-1;
215  if (VABS(pt[1] - ymax) < Vcompare) jhi = ny-1;
216  if (VABS(pt[2] - zmax) < Vcompare) khi = nz-1;
217 
218  /* See if we're on the mesh */
219  if ((ihi<nx) && (jhi<ny) && (khi<nz) &&
220  (ilo>=0) && (jlo>=0) && (klo>=0)) {
221 
222  dx = ifloat - (double)(ilo);
223  dy = jfloat - (double)(jlo);
224  dz = kfloat - (double)(klo);
225  u = dx *dy *dz *(thee->data[IJK(ihi,jhi,khi)])
226  + dx *(1.0-dy)*dz *(thee->data[IJK(ihi,jlo,khi)])
227  + dx *dy *(1.0-dz)*(thee->data[IJK(ihi,jhi,klo)])
228  + dx *(1.0-dy)*(1.0-dz)*(thee->data[IJK(ihi,jlo,klo)])
229  + (1.0-dx)*dy *dz *(thee->data[IJK(ilo,jhi,khi)])
230  + (1.0-dx)*(1.0-dy)*dz *(thee->data[IJK(ilo,jlo,khi)])
231  + (1.0-dx)*dy *(1.0-dz)*(thee->data[IJK(ilo,jhi,klo)])
232  + (1.0-dx)*(1.0-dy)*(1.0-dz)*(thee->data[IJK(ilo,jlo,klo)]);
233 
234  *value = u;
235 
236  if (isnan(u)) {
237  Vnm_print(2, "Vgrid_value: Got NaN!\n");
238  Vnm_print(2, "Vgrid_value: (x, y, z) = (%4.3f, %4.3f, %4.3f)\n",
239  pt[0], pt[1], pt[2]);
240  Vnm_print(2, "Vgrid_value: (ihi, jhi, khi) = (%d, %d, %d)\n",
241  ihi, jhi, khi);
242  Vnm_print(2, "Vgrid_value: (ilo, jlo, klo) = (%d, %d, %d)\n",
243  ilo, jlo, klo);
244  Vnm_print(2, "Vgrid_value: (nx, ny, nz) = (%d, %d, %d)\n",
245  nx, ny, nz);
246  Vnm_print(2, "Vgrid_value: (dx, dy, dz) = (%4.3f, %4.3f, %4.3f)\n",
247  dx, dy, dz);
248  Vnm_print(2, "Vgrid_value: data[IJK(ihi,jhi,khi)] = %g\n",
249  thee->data[IJK(ihi,jhi,khi)]);
250  Vnm_print(2, "Vgrid_value: data[IJK(ihi,jlo,khi)] = %g\n",
251  thee->data[IJK(ihi,jlo,khi)]);
252  Vnm_print(2, "Vgrid_value: data[IJK(ihi,jhi,klo)] = %g\n",
253  thee->data[IJK(ihi,jhi,klo)]);
254  Vnm_print(2, "Vgrid_value: data[IJK(ihi,jlo,klo)] = %g\n",
255  thee->data[IJK(ihi,jlo,klo)]);
256  Vnm_print(2, "Vgrid_value: data[IJK(ilo,jhi,khi)] = %g\n",
257  thee->data[IJK(ilo,jhi,khi)]);
258  Vnm_print(2, "Vgrid_value: data[IJK(ilo,jlo,khi)] = %g\n",
259  thee->data[IJK(ilo,jlo,khi)]);
260  Vnm_print(2, "Vgrid_value: data[IJK(ilo,jhi,klo)] = %g\n",
261  thee->data[IJK(ilo,jhi,klo)]);
262  Vnm_print(2, "Vgrid_value: data[IJK(ilo,jlo,klo)] = %g\n",
263  thee->data[IJK(ilo,jlo,klo)]);
264  }
265  return 1;
266 
267  } else {
268 
269  *value = 0;
270  return 0;
271 
272  }
273 
274  return 0;
275 
276 }
277 
278 /* ///////////////////////////////////////////////////////////////////////////
279 // Routine: Vgrid_curvature
280 //
281 // Notes: cflag=0 ==> Reduced Maximal Curvature
282 // cflag=1 ==> Mean Curvature (Laplace)
283 // cflag=2 ==> Gauss Curvature
284 // cflag=3 ==> True Maximal Curvature
285 //
286 // Authors: Stephen Bond and Nathan Baker
288 VPUBLIC int Vgrid_curvature(Vgrid *thee, double pt[3], int cflag,
289  double *value) {
290 
291  double hx, hy, hzed, curv;
292  double dxx, dyy, dzz;
293  double uleft, umid, uright, testpt[3];
294 
295  if (thee == VNULL) {
296  Vnm_print(2, "Vgrid_curvature: Error -- got VNULL thee!\n");
297  VASSERT(0);
298  }
299  if (!(thee->ctordata || thee->readdata)) {
300  Vnm_print(2, "Vgrid_curvature: Error -- no data available!\n");
301  VASSERT(0);
302  }
303 
304  hx = thee->hx;
305  hy = thee->hy;
306  hzed = thee->hzed;
307 
308  curv = 0.0;
309 
310  testpt[0] = pt[0];
311  testpt[1] = pt[1];
312  testpt[2] = pt[2];
313 
314  /* Compute 2nd derivative in the x-direction */
315  VJMPERR1(Vgrid_value( thee, testpt, &umid));
316  testpt[0] = pt[0] - hx;
317  VJMPERR1(Vgrid_value( thee, testpt, &uleft));
318  testpt[0] = pt[0] + hx;
319  VJMPERR1(Vgrid_value( thee, testpt, &uright));
320  testpt[0] = pt[0];
321 
322  dxx = (uright - 2*umid + uleft)/(hx*hx);
323 
324  /* Compute 2nd derivative in the y-direction */
325  VJMPERR1(Vgrid_value( thee, testpt, &umid));
326  testpt[1] = pt[1] - hy;
327  VJMPERR1(Vgrid_value( thee, testpt, &uleft));
328  testpt[1] = pt[1] + hy;
329  VJMPERR1(Vgrid_value( thee, testpt, &uright));
330  testpt[1] = pt[1];
331 
332  dyy = (uright - 2*umid + uleft)/(hy*hy);
333 
334  /* Compute 2nd derivative in the z-direction */
335  VJMPERR1(Vgrid_value( thee, testpt, &umid));
336  testpt[2] = pt[2] - hzed;
337  VJMPERR1(Vgrid_value( thee, testpt, &uleft));
338  testpt[2] = pt[2] + hzed;
339  VJMPERR1(Vgrid_value( thee, testpt, &uright));
340 
341  dzz = (uright - 2*umid + uleft)/(hzed*hzed);
342 
343 
344  if ( cflag == 0 ) {
345  curv = fabs(dxx);
346  curv = ( curv > fabs(dyy) ) ? curv : fabs(dyy);
347  curv = ( curv > fabs(dzz) ) ? curv : fabs(dzz);
348  } else if ( cflag == 1 ) {
349  curv = (dxx + dyy + dzz)/3.0;
350  } else {
351  Vnm_print(2, "Vgrid_curvature: support for cflag = %d not available!\n", cflag);
352  VASSERT( 0 ); /* Feature Not Coded Yet! */
353  }
354 
355  *value = curv;
356  return 1;
357 
358  VERROR1:
359  return 0;
360 
361 }
362 
363 /* ///////////////////////////////////////////////////////////////////////////
364 // Routine: Vgrid_gradient
365 //
366 // Authors: Nathan Baker and Stephen Bond
368 VPUBLIC int Vgrid_gradient(Vgrid *thee, double pt[3], double grad[3]) {
369 
370  double hx, hy, hzed;
371  double uleft, umid, uright, testpt[3];
372  int haveleft, haveright;
373 
374  if (thee == VNULL) {
375  Vnm_print(2, "Vgrid_gradient: Error -- got VNULL thee!\n");
376  VASSERT(0);
377  }
378  if (!(thee->ctordata || thee->readdata)) {
379  Vnm_print(2, "Vgrid_gradient: Error -- no data available!\n");
380  VASSERT(0);
381  }
382 
383  hx = thee->hx;
384  hy = thee->hy;
385  hzed = thee->hzed;
386 
387  /* Compute derivative in the x-direction */
388  testpt[0] = pt[0];
389  testpt[1] = pt[1];
390  testpt[2] = pt[2];
391  VJMPERR1( Vgrid_value( thee, testpt, &umid));
392  testpt[0] = pt[0] - hx;
393  if (Vgrid_value( thee, testpt, &uleft)) haveleft = 1;
394  else haveleft = 0;
395  testpt[0] = pt[0] + hx;
396  if (Vgrid_value( thee, testpt, &uright)) haveright = 1;
397  else haveright = 0;
398  if (haveright && haveleft) grad[0] = (uright - uleft)/(2*hx);
399  else if (haveright) grad[0] = (uright - umid)/hx;
400  else if (haveleft) grad[0] = (umid - uleft)/hx;
401  else VJMPERR1(0);
402 
403  /* Compute derivative in the y-direction */
404  testpt[0] = pt[0];
405  testpt[1] = pt[1];
406  testpt[2] = pt[2];
407  VJMPERR1(Vgrid_value(thee, testpt, &umid));
408  testpt[1] = pt[1] - hy;
409  if (Vgrid_value( thee, testpt, &uleft)) haveleft = 1;
410  else haveleft = 0;
411  testpt[1] = pt[1] + hy;
412  if (Vgrid_value( thee, testpt, &uright)) haveright = 1;
413  else haveright = 0;
414  if (haveright && haveleft) grad[1] = (uright - uleft)/(2*hy);
415  else if (haveright) grad[1] = (uright - umid)/hy;
416  else if (haveleft) grad[1] = (umid - uleft)/hy;
417  else VJMPERR1(0);
418 
419  /* Compute derivative in the z-direction */
420  testpt[0] = pt[0];
421  testpt[1] = pt[1];
422  testpt[2] = pt[2];
423  VJMPERR1(Vgrid_value(thee, testpt, &umid));
424  testpt[2] = pt[2] - hzed;
425  if (Vgrid_value( thee, testpt, &uleft)) haveleft = 1;
426  else haveleft = 0;
427  testpt[2] = pt[2] + hzed;
428  if (Vgrid_value( thee, testpt, &uright)) haveright = 1;
429  else haveright = 0;
430  if (haveright && haveleft) grad[2] = (uright - uleft)/(2*hzed);
431  else if (haveright) grad[2] = (uright - umid)/hzed;
432  else if (haveleft) grad[2] = (umid - uleft)/hzed;
433  else VJMPERR1(0);
434 
435  return 1;
436 
437  VERROR1:
438  return 0;
439 
440 }
441 
442 /* ///////////////////////////////////////////////////////////////////////////
443  // Routine: Vgrid_readGZ
444  //
445  // Author: David Gohara
447 #ifdef HAVE_ZLIB
448 #define off_t long
449 #include "zlib.h"
450 #endif
451 VPUBLIC int Vgrid_readGZ(Vgrid *thee, const char *fname) {
452 
453 #ifdef HAVE_ZLIB
454  int i, j, k;
455  int len; // Temporary counter variable for loop conditionals
456  int q, itmp, u, header, incr;
457  double *temp;
458  double dtmp1, dtmp2, dtmp3;
459  gzFile infile;
460  char line[VMAX_ARGLEN];
461 
462  header = 0;
463 
464  /* Check to see if the existing data is null and, if not, clear it out */
465  if (thee->data != VNULL) {
466  Vnm_print(1, "%s: destroying existing data!\n", __func__);
467  Vmem_free(thee->mem, thee->nx * thee->ny * thee->nz, sizeof(double),
468  (void **)&(thee->data));
469  }
470 
471  thee->readdata = 1;
472  thee->ctordata = 0;
473 
474  infile = gzopen(fname, "rb");
475  if (infile == Z_NULL) {
476  Vnm_print(2, "%s: Problem opening compressed file %s\n", __func__, fname);
477  return VRC_FAILURE;
478  }
479 
480  thee->hx = 0.0;
481  thee->hy = 0.0;
482  thee->hzed = 0.0;
483 
484  //read data here
485  while (header < 7) {
486  if(gzgets(infile, line, VMAX_ARGLEN) == Z_NULL){
487  return VRC_FAILURE;
488  }
489 
490  // Skip comments and newlines
491  if(strncmp(line, "#", 1) == 0) continue;
492  if(line[0] == '\n') continue;
493 
494  switch (header) {
495  case 0:
496  sscanf(line, "object 1 class gridpositions counts %d %d %d",
497  &(thee->nx),&(thee->ny),&(thee->nz));
498  break;
499  case 1:
500  sscanf(line, "origin %lf %lf %lf",
501  &(thee->xmin),&(thee->ymin),&(thee->zmin));
502  break;
503  case 2:
504  case 3:
505  case 4:
506  sscanf(line, "delta %lf %lf %lf",&dtmp1,&dtmp2,&dtmp3);
507  thee->hx += dtmp1;
508  thee->hy += dtmp2;
509  thee->hzed += dtmp3;
510  break;
511  default:
512  break;
513  }
514 
515  header++;
516  }
517 
518  /* Allocate space for the data */
519  Vnm_print(0, "%s: allocating %d x %d x %d doubles for storage\n",
520  __func__, thee->nx, thee->ny, thee->nz);
521  len = thee->nx * thee->ny * thee->nz;
522 
523  thee->data = VNULL;
524  thee->data = Vmem_malloc(thee->mem, len, sizeof(double));
525  if (thee->data == VNULL) {
526  Vnm_print(2, "%s: Unable to allocate space for data!\n", __func__);
527  return 0;
528  }
529 
530  /* Allocate a temporary buffer to store the compressed
531  * data into (column major order). Add 2 to ensure the buffer is
532  * big enough to take extra data on the final read loop.
533  */
534  temp = (double *)malloc(len * (2 * sizeof(double)));
535 
536  for (i = 0; i < len; i += 3){
537  memset(&line, 0, sizeof(line));
538  gzgets(infile, line, VMAX_ARGLEN);
539  sscanf(line, "%lf %lf %lf", &temp[i], &temp[i+1], &temp[i+2]);
540  }
541 
542  /* Now move the data to row major order */
543  incr = 0;
544  for (i=0; i<thee->nx; i++) {
545  for (j=0; j<thee->ny; j++) {
546  for (k=0; k<thee->nz; k++) {
547  u = k*(thee->nx)*(thee->ny)+j*(thee->nx)+i;
548  (thee->data)[u] = temp[incr++];
549  }
550  }
551  }
552 
553  /* calculate grid maxima */
554  thee->xmax = thee->xmin + (thee->nx-1)*thee->hx;
555  thee->ymax = thee->ymin + (thee->ny-1)*thee->hy;
556  thee->zmax = thee->zmin + (thee->nz-1)*thee->hzed;
557 
558  /* Close off the socket */
559  gzclose(infile);
560  free(temp);
561 #else
562 
563  Vnm_print(0, "WARNING\n");
564  Vnm_print(0, "Vgrid_readGZ: gzip read/write support is disabled in this build\n");
565  Vnm_print(0, "Vgrid_readGZ: configure and compile without the --disable-zlib flag.\n");
566  Vnm_print(0, "WARNING\n");
567 #endif
568  return VRC_SUCCESS;
569 }
570 
575 VPUBLIC int Vgrid_readDX(Vgrid *thee,
576  const char *iodev,
577  const char *iofmt,
578  const char *thost,
579  const char *fname
580  ) {
581 
582  int i,
583  j,
584  k,
585  itmp,
586  u;
587  double dtmp;
588  char tok[VMAX_BUFSIZE];
589  Vio *sock;
590 
591  /* Check to see if the existing data is null and, if not, clear it out */
592  if (thee->data != VNULL) {
593  Vnm_print(1, "Vgrid_readDX: destroying existing data!\n");
594  Vmem_free(thee->mem, (thee->nx*thee->ny*thee->nz), sizeof(double),
595  (void **)&(thee->data)); }
596  thee->readdata = 1;
597  thee->ctordata = 0;
598 
599  /* Set up the virtual socket */
600  sock = Vio_ctor(iodev,iofmt,thost,fname,"r");
601  if (sock == VNULL) {
602  Vnm_print(2, "Vgrid_readDX: Problem opening virtual socket %s\n",
603  fname);
604  return 0;
605  }
606  if (Vio_accept(sock, 0) < 0) {
607  Vnm_print(2, "Vgrid_readDX: Problem accepting virtual socket %s\n",
608  fname);
609  return 0;
610  }
611 
612  Vio_setWhiteChars(sock, MCwhiteChars);
613  Vio_setCommChars(sock, MCcommChars);
614 
615  /* Read in the DX regular positions */
616  /* Get "object" */
617  VJMPERR2(1 == Vio_scanf(sock, "%s", tok));
618  VJMPERR1(!strcmp(tok, "object"));
619  /* Get "1" */
620  VJMPERR2(1 == Vio_scanf(sock, "%d", &itmp));
621  /* Get "class" */
622  VJMPERR2(1 == Vio_scanf(sock, "%s", tok));
623  VJMPERR1(!strcmp(tok, "class"));
624  /* Get "gridpositions" */
625  VJMPERR2(1 == Vio_scanf(sock, "%s", tok));
626  VJMPERR1(!strcmp(tok, "gridpositions"));
627  /* Get "counts" */
628  VJMPERR2(1 == Vio_scanf(sock, "%s", tok));
629  VJMPERR1(!strcmp(tok, "counts"));
630  /* Get nx */
631  VJMPERR2(1 == Vio_scanf(sock, "%s", tok));
632  VJMPERR1(1 == sscanf(tok, "%d", &(thee->nx)));
633  /* Get ny */
634  VJMPERR2(1 == Vio_scanf(sock, "%s", tok));
635  VJMPERR1(1 == sscanf(tok, "%d", &(thee->ny)));
636  /* Get nz */
637  VJMPERR2(1 == Vio_scanf(sock, "%s", tok));
638  VJMPERR1(1 == sscanf(tok, "%d", &(thee->nz)));
639  Vnm_print(0, "Vgrid_readDX: Grid dimensions %d x %d x %d grid\n",
640  thee->nx, thee->ny, thee->nz);
641  /* Get "origin" */
642  VJMPERR2(1 == Vio_scanf(sock, "%s", tok));
643  VJMPERR1(!strcmp(tok, "origin"));
644  /* Get xmin */
645  VJMPERR2(1 == Vio_scanf(sock, "%s", tok));
646  VJMPERR1(1 == sscanf(tok, "%lf", &(thee->xmin)));
647  /* Get ymin */
648  VJMPERR2(1 == Vio_scanf(sock, "%s", tok));
649  VJMPERR1(1 == sscanf(tok, "%lf", &(thee->ymin)));
650  /* Get zmin */
651  VJMPERR2(1 == Vio_scanf(sock, "%s", tok));
652  VJMPERR1(1 == sscanf(tok, "%lf", &(thee->zmin)));
653  Vnm_print(0, "Vgrid_readDX: Grid origin = (%g, %g, %g)\n",
654  thee->xmin, thee->ymin, thee->zmin);
655  /* Get "delta" */
656  VJMPERR2(1 == Vio_scanf(sock, "%s", tok));
657  VJMPERR1(!strcmp(tok, "delta"));
658  /* Get hx */
659  VJMPERR2(1 == Vio_scanf(sock, "%s", tok));
660  VJMPERR1(1 == sscanf(tok, "%lf", &(thee->hx)));
661  /* Get 0.0 */
662  VJMPERR2(1 == Vio_scanf(sock, "%s", tok));
663  VJMPERR1(1 == sscanf(tok, "%lf", &dtmp));
664  VJMPERR1(dtmp == 0.0);
665  /* Get 0.0 */
666  VJMPERR2(1 == Vio_scanf(sock, "%s", tok));
667  VJMPERR1(1 == sscanf(tok, "%lf", &dtmp));
668  VJMPERR1(dtmp == 0.0);
669  /* Get "delta" */
670  VJMPERR2(1 == Vio_scanf(sock, "%s", tok));
671  VJMPERR1(!strcmp(tok, "delta"));
672  /* Get 0.0 */
673  VJMPERR2(1 == Vio_scanf(sock, "%s", tok));
674  VJMPERR1(1 == sscanf(tok, "%lf", &dtmp));
675  VJMPERR1(dtmp == 0.0);
676  /* Get hy */
677  VJMPERR2(1 == Vio_scanf(sock, "%s", tok));
678  VJMPERR1(1 == sscanf(tok, "%lf", &(thee->hy)));
679  /* Get 0.0 */
680  VJMPERR2(1 == Vio_scanf(sock, "%s", tok));
681  VJMPERR1(1 == sscanf(tok, "%lf", &dtmp));
682  VJMPERR1(dtmp == 0.0);
683  /* Get "delta" */
684  VJMPERR2(1 == Vio_scanf(sock, "%s", tok));
685  VJMPERR1(!strcmp(tok, "delta"));
686  /* Get 0.0 */
687  VJMPERR2(1 == Vio_scanf(sock, "%s", tok));
688  VJMPERR1(1 == sscanf(tok, "%lf", &dtmp));
689  VJMPERR1(dtmp == 0.0);
690  /* Get 0.0 */
691  VJMPERR2(1 == Vio_scanf(sock, "%s", tok));
692  VJMPERR1(1 == sscanf(tok, "%lf", &dtmp));
693  VJMPERR1(dtmp == 0.0);
694  /* Get hz */
695  VJMPERR2(1 == Vio_scanf(sock, "%s", tok));
696  VJMPERR1(1 == sscanf(tok, "%lf", &(thee->hzed)));
697  Vnm_print(0, "Vgrid_readDX: Grid spacings = (%g, %g, %g)\n",
698  thee->hx, thee->hy, thee->hzed);
699  /* Get "object" */
700  VJMPERR2(1 == Vio_scanf(sock, "%s", tok));
701  VJMPERR1(!strcmp(tok, "object"));
702  /* Get "2" */
703  VJMPERR2(1 == Vio_scanf(sock, "%s", tok));
704  /* Get "class" */
705  VJMPERR2(1 == Vio_scanf(sock, "%s", tok));
706  VJMPERR1(!strcmp(tok, "class"));
707  /* Get "gridconnections" */
708  VJMPERR2(1 == Vio_scanf(sock, "%s", tok));
709  VJMPERR1(!strcmp(tok, "gridconnections"));
710  /* Get "counts" */
711  VJMPERR2(1 == Vio_scanf(sock, "%s", tok));
712  VJMPERR1(!strcmp(tok, "counts"));
713  /* Get the dimensions again */
714  VJMPERR2(1 == Vio_scanf(sock, "%s", tok));
715  VJMPERR2(1 == Vio_scanf(sock, "%s", tok));
716  VJMPERR2(1 == Vio_scanf(sock, "%s", tok));
717  /* Get "object" */
718  VJMPERR2(1 == Vio_scanf(sock, "%s", tok));
719  VJMPERR1(!strcmp(tok, "object"));
720  /* Get # */
721  VJMPERR2(1 == Vio_scanf(sock, "%s", tok));
722  /* Get "class" */
723  VJMPERR2(1 == Vio_scanf(sock, "%s", tok));
724  VJMPERR1(!strcmp(tok, "class"));
725  /* Get "array" */
726  VJMPERR2(1 == Vio_scanf(sock, "%s", tok));
727  VJMPERR1(!strcmp(tok, "array"));
728  /* Get "type" */
729  VJMPERR2(1 == Vio_scanf(sock, "%s", tok));
730  VJMPERR1(!strcmp(tok, "type"));
731  /* Get "double" */
732  VJMPERR2(1 == Vio_scanf(sock, "%s", tok));
733  VJMPERR1(!strcmp(tok, "double"));
734  /* Get "rank" */
735  VJMPERR2(1 == Vio_scanf(sock, "%s", tok));
736  VJMPERR1(!strcmp(tok, "rank"));
737  /* Get # */
738  VJMPERR2(1 == Vio_scanf(sock, "%s", tok));
739  /* Get "items" */
740  VJMPERR2(1 == Vio_scanf(sock, "%s", tok));
741  VJMPERR1(!strcmp(tok, "items"));
742  /* Get # */
743  VJMPERR2(1 == Vio_scanf(sock, "%s", tok));
744  VJMPERR1(1 == sscanf(tok, "%d", &itmp));
745  VJMPERR1(((thee->nx)*(thee->ny)*(thee->nz)) == itmp);
746  /* Get "data" */
747  VJMPERR2(1 == Vio_scanf(sock, "%s", tok));
748  VJMPERR1(!strcmp(tok, "data"));
749  /* Get "follows" */
750  VJMPERR2(1 == Vio_scanf(sock, "%s", tok));
751  VJMPERR1(!strcmp(tok, "follows"));
752 
753  /* Allocate space for the data */
754  Vnm_print(0, "Vgrid_readDX: allocating %d x %d x %d doubles for storage\n",
755  thee->nx, thee->ny, thee->nz);
756  thee->data = VNULL;
757  thee->data = (double*)Vmem_malloc(thee->mem, (thee->nx)*(thee->ny)*(thee->nz),
758  sizeof(double));
759  if (thee->data == VNULL) {
760  Vnm_print(2, "Vgrid_readDX: Unable to allocate space for data!\n");
761  return 0;
762  }
763 
764  for (i=0; i<thee->nx; i++) {
765  for (j=0; j<thee->ny; j++) {
766  for (k=0; k<thee->nz; k++) {
767  u = k*(thee->nx)*(thee->ny)+j*(thee->nx)+i;
768  VJMPERR2(1 == Vio_scanf(sock, "%s", tok));
769  VJMPERR1(1 == sscanf(tok, "%lf", &dtmp));
770  (thee->data)[u] = dtmp;
771  }
772  }
773  }
774 
775  /* calculate grid maxima */
776  thee->xmax = thee->xmin + (thee->nx-1)*thee->hx;
777  thee->ymax = thee->ymin + (thee->ny-1)*thee->hy;
778  thee->zmax = thee->zmin + (thee->nz-1)*thee->hzed;
779 
780  /* Close off the socket */
781  Vio_acceptFree(sock);
782  Vio_dtor(&sock);
783 
784  return 1;
785 
786  VERROR1:
787  Vio_dtor(&sock);
788  Vnm_print(2, "Vgrid_readDX: Format problem with input file <%s>\n",
789  fname);
790  return 0;
791 
792  VERROR2:
793  Vio_dtor(&sock);
794  Vnm_print(2, "Vgrid_readDX: I/O problem with input file <%s>\n",
795  fname);
796  return 0;
797 
798 
799 
800 }
801 
802 /* ///////////////////////////////////////////////////////////////////////////
803  // Routine: Vgrid_writeGZ
804  //
805  // Author: Nathan Baker
807 VPUBLIC void Vgrid_writeGZ(Vgrid *thee, const char *iodev, const char *iofmt,
808  const char *thost, const char *fname, char *title, double *pvec) {
809 
810 #ifdef HAVE_ZLIB
811  double xmin, ymin, zmin, hx, hy, hzed;
812 
813  int nx, ny, nz;
814  int icol, i, j, k, u, usepart, nxPART, nyPART, nzPART, gotit;
815  double x, y, z, xminPART, yminPART, zminPART;
816 
817  int txyz;
818  double txmin, tymin, tzmin;
819 
820  char header[8196];
821  char footer[8196];
822  char line[80];
823  char newline[] = "\n";
824  gzFile outfile;
825  char precFormat[VMAX_BUFSIZE];
826 
827  if (thee == VNULL) {
828  Vnm_print(2, "Vgrid_writeGZ: Error -- got VNULL thee!\n");
829  VASSERT(0);
830  }
831  if (!(thee->ctordata || thee->readdata)) {
832  Vnm_print(2, "Vgrid_writeGZ: Error -- no data available!\n");
833  VASSERT(0);
834  }
835 
836  hx = thee->hx;
837  hy = thee->hy;
838  hzed = thee->hzed;
839  nx = thee->nx;
840  ny = thee->ny;
841  nz = thee->nz;
842  xmin = thee->xmin;
843  ymin = thee->ymin;
844  zmin = thee->zmin;
845 
846  if (pvec == VNULL) usepart = 0;
847  else usepart = 1;
848 
849  /* Set up the virtual socket */
850  Vnm_print(0, "Vgrid_writeGZ: Opening file...\n");
851  outfile = gzopen(fname, "wb");
852 
853  if (usepart) {
854  /* Get the lower corner and number of grid points for the local
855  * partition */
856  xminPART = VLARGE;
857  yminPART = VLARGE;
858  zminPART = VLARGE;
859  nxPART = 0;
860  nyPART = 0;
861  nzPART = 0;
862  /* First, search for the lower corner */
863  for (k=0; k<nz; k++) {
864  z = k*hzed + zmin;
865  for (j=0; j<ny; j++) {
866  y = j*hy + ymin;
867  for (i=0; i<nx; i++) {
868  x = i*hx + xmin;
869  if (pvec[IJK(i,j,k)] > 0.0) {
870  if (x < xminPART) xminPART = x;
871  if (y < yminPART) yminPART = y;
872  if (z < zminPART) zminPART = z;
873  }
874  }
875  }
876  }
877  /* Now search for the number of grid points in the z direction */
878  for (k=0; k<nz; k++) {
879  gotit = 0;
880  for (j=0; j<ny; j++) {
881  for (i=0; i<nx; i++) {
882  if (pvec[IJK(i,j,k)] > 0.0) {
883  gotit = 1;
884  break;
885  }
886  }
887  if (gotit) break;
888  }
889  if (gotit) nzPART++;
890  }
891  /* Now search for the number of grid points in the y direction */
892  for (j=0; j<ny; j++) {
893  gotit = 0;
894  for (k=0; k<nz; k++) {
895  for (i=0; i<nx; i++) {
896  if (pvec[IJK(i,j,k)] > 0.0) {
897  gotit = 1;
898  break;
899  }
900  }
901  if (gotit) break;
902  }
903  if (gotit) nyPART++;
904  }
905  /* Now search for the number of grid points in the x direction */
906  for (i=0; i<nx; i++) {
907  gotit = 0;
908  for (k=0; k<nz; k++) {
909  for (j=0; j<ny; j++) {
910  if (pvec[IJK(i,j,k)] > 0.0) {
911  gotit = 1;
912  break;
913  }
914  }
915  if (gotit) break;
916  }
917  if (gotit) nxPART++;
918  }
919 
920  if ((nxPART != nx) || (nyPART != ny) || (nzPART != nz)) {
921  Vnm_print(0, "Vgrid_writeGZ: printing only subset of domain\n");
922  }
923 
924  txyz = (nxPART*nyPART*nzPART);
925  txmin = xminPART;
926  tymin = yminPART;
927  tzmin = zminPART;
928 
929  }else {
930 
931  txyz = (nx*ny*nz);
932  txmin = xmin;
933  tymin = ymin;
934  tzmin = zmin;
935 
936  }
937 
938  /* Write off the title (if we're not XDR) */
939  sprintf(header,
940  "# Data from %s\n" \
941  "# \n" \
942  "# %s\n" \
943  "# \n" \
944  "object 1 class gridpositions counts %i %i %i\n" \
945  "origin %12.6e %12.6e %12.6e\n" \
946  "delta %12.6e 0.000000e+00 0.000000e+00\n" \
947  "delta 0.000000e+00 %12.6e 0.000000e+00\n" \
948  "delta 0.000000e+00 0.000000e+00 %12.6e\n" \
949  "object 2 class gridconnections counts %i %i %i\n"\
950  "object 3 class array type double rank 0 items %i data follows\n",
951  PACKAGE_STRING,title,nx,ny,nz,txmin,tymin,tzmin,
952  hx,hy,hzed,nx,ny,nz,txyz);
953  gzwrite(outfile, header, strlen(header)*sizeof(char));
954 
955  /* Now write the data */
956  icol = 0;
957  for (i=0; i<nx; i++) {
958  for (j=0; j<ny; j++) {
959  for (k=0; k<nz; k++) {
960  u = k*(nx)*(ny)+j*(nx)+i;
961  if (pvec[u] > 0.0) {
962  sprintf(line, "%12.6e ", thee->data[u]);
963  gzwrite(outfile, line, strlen(line)*sizeof(char));
964  icol++;
965  if (icol == 3) {
966  icol = 0;
967  gzwrite(outfile, newline, strlen(newline)*sizeof(char));
968  }
969  }
970  }
971  }
972  }
973  if(icol < 3){
974  char newline[] = "\n";
975  gzwrite(outfile, newline, strlen(newline)*sizeof(char));
976  }
977 
978  /* Create the field */
979  sprintf(footer, "attribute \"dep\" string \"positions\"\n" \
980  "object \"regular positions regular connections\" class field\n" \
981  "component \"positions\" value 1\n" \
982  "component \"connections\" value 2\n" \
983  "component \"data\" value 3\n");
984  gzwrite(outfile, footer, strlen(footer)*sizeof(char));
985 
986  gzclose(outfile);
987 #else
988 
989  Vnm_print(0, "WARNING\n");
990  Vnm_print(0, "Vgrid_readGZ: gzip read/write support is disabled in this build\n");
991  Vnm_print(0, "Vgrid_readGZ: configure and compile without the --disable-zlib flag.\n");
992  Vnm_print(0, "WARNING\n");
993 #endif
994 }
995 
996 /* ///////////////////////////////////////////////////////////////////////////
997 // Routine: Vgrid_writeDX
998 //
999 // Author: Nathan Baker
1001 VPUBLIC void Vgrid_writeDX(Vgrid *thee, const char *iodev, const char *iofmt,
1002  const char *thost, const char *fname, char *title, double *pvec) {
1003 
1004  double xmin, ymin, zmin, hx, hy, hzed;
1005  int nx, ny, nz;
1006  int icol, i, j, k, u, usepart, nxPART, nyPART, nzPART, gotit;
1007  double x, y, z, xminPART, yminPART, zminPART;
1008  Vio *sock;
1009  char precFormat[VMAX_BUFSIZE];
1010 
1011  if (thee == VNULL) {
1012  Vnm_print(2, "Vgrid_writeDX: Error -- got VNULL thee!\n");
1013  VASSERT(0);
1014  }
1015  if (!(thee->ctordata || thee->readdata)) {
1016  Vnm_print(2, "Vgrid_writeDX: Error -- no data available!\n");
1017  VASSERT(0);
1018  }
1019 
1020  hx = thee->hx;
1021  hy = thee->hy;
1022  hzed = thee->hzed;
1023  nx = thee->nx;
1024  ny = thee->ny;
1025  nz = thee->nz;
1026  xmin = thee->xmin;
1027  ymin = thee->ymin;
1028  zmin = thee->zmin;
1029 
1030  if (pvec == VNULL) usepart = 0;
1031  else usepart = 1;
1032 
1033  /* Set up the virtual socket */
1034  Vnm_print(0, "Vgrid_writeDX: Opening virtual socket...\n");
1035  sock = Vio_ctor(iodev,iofmt,thost,fname,"w");
1036  if (sock == VNULL) {
1037  Vnm_print(2, "Vgrid_writeDX: Problem opening virtual socket %s\n",
1038  fname);
1039  return;
1040  }
1041  if (Vio_connect(sock, 0) < 0) {
1042  Vnm_print(2, "Vgrid_writeDX: Problem connecting virtual socket %s\n",
1043  fname);
1044  return;
1045  }
1046 
1047  Vio_setWhiteChars(sock, MCwhiteChars);
1048  Vio_setCommChars(sock, MCcommChars);
1049 
1050  Vnm_print(0, "Vgrid_writeDX: Writing to virtual socket...\n");
1051 
1052  if (usepart) {
1053  /* Get the lower corner and number of grid points for the local
1054  * partition */
1055  xminPART = VLARGE;
1056  yminPART = VLARGE;
1057  zminPART = VLARGE;
1058  nxPART = 0;
1059  nyPART = 0;
1060  nzPART = 0;
1061  /* First, search for the lower corner */
1062  for (k=0; k<nz; k++) {
1063  z = k*hzed + zmin;
1064  for (j=0; j<ny; j++) {
1065  y = j*hy + ymin;
1066  for (i=0; i<nx; i++) {
1067  x = i*hx + xmin;
1068  if (pvec[IJK(i,j,k)] > 0.0) {
1069  if (x < xminPART) xminPART = x;
1070  if (y < yminPART) yminPART = y;
1071  if (z < zminPART) zminPART = z;
1072  }
1073  }
1074  }
1075  }
1076  /* Now search for the number of grid points in the z direction */
1077  for (k=0; k<nz; k++) {
1078  gotit = 0;
1079  for (j=0; j<ny; j++) {
1080  for (i=0; i<nx; i++) {
1081  if (pvec[IJK(i,j,k)] > 0.0) {
1082  gotit = 1;
1083  break;
1084  }
1085  }
1086  if (gotit) break;
1087  }
1088  if (gotit) nzPART++;
1089  }
1090  /* Now search for the number of grid points in the y direction */
1091  for (j=0; j<ny; j++) {
1092  gotit = 0;
1093  for (k=0; k<nz; k++) {
1094  for (i=0; i<nx; i++) {
1095  if (pvec[IJK(i,j,k)] > 0.0) {
1096  gotit = 1;
1097  break;
1098  }
1099  }
1100  if (gotit) break;
1101  }
1102  if (gotit) nyPART++;
1103  }
1104  /* Now search for the number of grid points in the x direction */
1105  for (i=0; i<nx; i++) {
1106  gotit = 0;
1107  for (k=0; k<nz; k++) {
1108  for (j=0; j<ny; j++) {
1109  if (pvec[IJK(i,j,k)] > 0.0) {
1110  gotit = 1;
1111  break;
1112  }
1113  }
1114  if (gotit) break;
1115  }
1116  if (gotit) nxPART++;
1117  }
1118 
1119  if ((nxPART != nx) || (nyPART != ny) || (nzPART != nz)) {
1120  Vnm_print(0, "Vgrid_writeDX: printing only subset of domain\n");
1121  }
1122 
1123 
1124  /* Write off the title (if we're not XDR) */
1125  if (Vstring_strcasecmp(iofmt, "XDR") == 0) {
1126  Vnm_print(0, "Vgrid_writeDX: Skipping comments for XDR format.\n");
1127  } else {
1128  Vnm_print(0, "Vgrid_writeDX: Writing comments for %s format.\n",
1129  iofmt);
1130  Vio_printf(sock, "# Data from %s\n", PACKAGE_STRING);
1131  Vio_printf(sock, "# \n");
1132  Vio_printf(sock, "# %s\n", title);
1133  Vio_printf(sock, "# \n");
1134  }
1135 
1136  /* Write off the DX regular positions */
1137  Vio_printf(sock, "object 1 class gridpositions counts %d %d %d\n",
1138  nxPART, nyPART, nzPART);
1139 
1140  sprintf(precFormat, Vprecision, xminPART, yminPART, zminPART);
1141  Vio_printf(sock, "origin %s\n", precFormat);
1142  sprintf(precFormat, Vprecision, hx, 0.0, 0.0);
1143  Vio_printf(sock, "delta %s\n", precFormat);
1144  sprintf(precFormat, Vprecision, 0.0, hy, 0.0);
1145  Vio_printf(sock, "delta %s\n", precFormat);
1146  sprintf(precFormat, Vprecision, 0.0, 0.0, hzed);
1147  Vio_printf(sock, "delta %s\n", precFormat);
1148 
1149  /* Write off the DX regular connections */
1150  Vio_printf(sock, "object 2 class gridconnections counts %d %d %d\n",
1151  nxPART, nyPART, nzPART);
1152 
1153  /* Write off the DX data */
1154  Vio_printf(sock, "object 3 class array type double rank 0 items %d \
1155 data follows\n", (nxPART*nyPART*nzPART));
1156  icol = 0;
1157  for (i=0; i<nx; i++) {
1158  for (j=0; j<ny; j++) {
1159  for (k=0; k<nz; k++) {
1160  u = k*(nx)*(ny)+j*(nx)+i;
1161  if (pvec[u] > 0.0) {
1162  Vio_printf(sock, "%12.6e ", thee->data[u]);
1163  icol++;
1164  if (icol == 3) {
1165  icol = 0;
1166  Vio_printf(sock, "\n");
1167  }
1168  }
1169  }
1170  }
1171  }
1172 
1173  if (icol != 0) Vio_printf(sock, "\n");
1174 
1175  /* Create the field */
1176  Vio_printf(sock, "attribute \"dep\" string \"positions\"\n");
1177  Vio_printf(sock, "object \"regular positions regular connections\" \
1178 class field\n");
1179  Vio_printf(sock, "component \"positions\" value 1\n");
1180  Vio_printf(sock, "component \"connections\" value 2\n");
1181  Vio_printf(sock, "component \"data\" value 3\n");
1182 
1183  } else {
1184  /* Write off the title (if we're not XDR) */
1185  if (Vstring_strcasecmp(iofmt, "XDR") == 0) {
1186  Vnm_print(0, "Vgrid_writeDX: Skipping comments for XDR format.\n");
1187  } else {
1188  Vnm_print(0, "Vgrid_writeDX: Writing comments for %s format.\n",
1189  iofmt);
1190  Vio_printf(sock, "# Data from %s\n", PACKAGE_STRING);
1191  Vio_printf(sock, "# \n");
1192  Vio_printf(sock, "# %s\n", title);
1193  Vio_printf(sock, "# \n");
1194  }
1195 
1196 
1197  /* Write off the DX regular positions */
1198  Vio_printf(sock, "object 1 class gridpositions counts %d %d %d\n",
1199  nx, ny, nz);
1200 
1201  sprintf(precFormat, Vprecision, xmin, ymin, zmin);
1202  Vio_printf(sock, "origin %s\n", precFormat);
1203  sprintf(precFormat, Vprecision, hx, 0.0, 0.0);
1204  Vio_printf(sock, "delta %s\n", precFormat);
1205  sprintf(precFormat, Vprecision, 0.0, hy, 0.0);
1206  Vio_printf(sock, "delta %s\n", precFormat);
1207  sprintf(precFormat, Vprecision, 0.0, 0.0, hzed);
1208  Vio_printf(sock, "delta %s\n", precFormat);
1209 
1210  /* Write off the DX regular connections */
1211  Vio_printf(sock, "object 2 class gridconnections counts %d %d %d\n",
1212  nx, ny, nz);
1213 
1214  /* Write off the DX data */
1215  Vio_printf(sock, "object 3 class array type double rank 0 items %d \
1216 data follows\n", (nx*ny*nz));
1217  icol = 0;
1218  for (i=0; i<nx; i++) {
1219  for (j=0; j<ny; j++) {
1220  for (k=0; k<nz; k++) {
1221  u = k*(nx)*(ny)+j*(nx)+i;
1222  Vio_printf(sock, "%12.6e ", thee->data[u]);
1223  icol++;
1224  if (icol == 3) {
1225  icol = 0;
1226  Vio_printf(sock, "\n");
1227  }
1228  }
1229  }
1230  }
1231  if (icol != 0) Vio_printf(sock, "\n");
1232 
1233  /* Create the field */
1234  Vio_printf(sock, "attribute \"dep\" string \"positions\"\n");
1235  Vio_printf(sock, "object \"regular positions regular connections\" \
1236 class field\n");
1237  Vio_printf(sock, "component \"positions\" value 1\n");
1238  Vio_printf(sock, "component \"connections\" value 2\n");
1239  Vio_printf(sock, "component \"data\" value 3\n");
1240  }
1241 
1242  /* Close off the socket */
1243  Vio_connectFree(sock);
1244  Vio_dtor(&sock);
1245 }
1246 
1247 /* ///////////////////////////////////////////////////////////////////////////
1248 // Routine: Vgrid_writeUHBD
1249 // Author: Nathan Baker
1251 VPUBLIC void Vgrid_writeUHBD(Vgrid *thee, const char *iodev, const char *iofmt,
1252  const char *thost, const char *fname, char *title, double *pvec) {
1253 
1254  int icol, i, j, k, u, nx, ny, nz, gotit;
1255  double xmin, ymin, zmin, hzed, hy, hx;
1256  Vio *sock;
1257 
1258  if (thee == VNULL) {
1259  Vnm_print(2, "Vgrid_writeUHBD: Error -- got VNULL thee!\n");
1260  VASSERT(0);
1261  }
1262  if (!(thee->ctordata || thee->readdata)) {
1263  Vnm_print(2, "Vgrid_writeUHBD: Error -- no data available!\n");
1264  VASSERT(0);
1265  }
1266 
1267  if ((thee->hx!=thee->hy) || (thee->hy!=thee->hzed)
1268  || (thee->hx!=thee->hzed)) {
1269  Vnm_print(2, "Vgrid_writeUHBD: can't write UHBD mesh with non-uniform \
1270 spacing\n");
1271  return;
1272  }
1273 
1274  /* Set up the virtual socket */
1275  sock = Vio_ctor(iodev,iofmt,thost,fname,"w");
1276  if (sock == VNULL) {
1277  Vnm_print(2, "Vgrid_writeUHBD: Problem opening virtual socket %s\n",
1278  fname);
1279  return;
1280  }
1281  if (Vio_connect(sock, 0) < 0) {
1282  Vnm_print(2, "Vgrid_writeUHBD: Problem connecting virtual socket %s\n",
1283  fname);
1284  return;
1285  }
1286 
1287  /* Get the lower corner and number of grid points for the local
1288  * partition */
1289  hx = thee->hx;
1290  hy = thee->hy;
1291  hzed = thee->hzed;
1292  nx = thee->nx;
1293  ny = thee->ny;
1294  nz = thee->nz;
1295  xmin = thee->xmin;
1296  ymin = thee->ymin;
1297  zmin = thee->zmin;
1298 
1299  /* Let interested folks know that partition information is ignored */
1300  if (pvec != VNULL) {
1301  gotit = 0;
1302  for (i=0; i<(nx*ny*nz); i++) {
1303  if (pvec[i] == 0) {
1304  gotit = 1;
1305  break;
1306  }
1307  }
1308  if (gotit) {
1309  Vnm_print(2, "Vgrid_writeUHBD: IGNORING PARTITION INFORMATION!\n");
1310  Vnm_print(2, "Vgrid_writeUHBD: This means I/O from parallel runs \
1311 will have significant overlap.\n");
1312  }
1313  }
1314 
1315  /* Write out the header */
1316  Vio_printf(sock, "%72s\n", title);
1317  Vio_printf(sock, "%12.5e%12.5e%7d%7d%7d%7d%7d\n", 1.0, 0.0, -1, 0,
1318  nz, 1, nz);
1319  Vio_printf(sock, "%7d%7d%7d%12.5e%12.5e%12.5e%12.5e\n", nx, ny, nz,
1320  hx, (xmin-hx), (ymin-hx), (zmin-hx));
1321  Vio_printf(sock, "%12.5e%12.5e%12.5e%12.5e\n", 0.0, 0.0, 0.0, 0.0);
1322  Vio_printf(sock, "%12.5e%12.5e%7d%7d", 0.0, 0.0, 0, 0);
1323 
1324  /* Write out the entries */
1325  icol = 0;
1326  for (k=0; k<nz; k++) {
1327  Vio_printf(sock, "\n%7d%7d%7d\n", k+1, thee->nx, thee->ny);
1328  icol = 0;
1329  for (j=0; j<ny; j++) {
1330  for (i=0; i<nx; i++) {
1331  u = k*(nx)*(ny)+j*(nx)+i;
1332  icol++;
1333  Vio_printf(sock, " %12.5e", thee->data[u]);
1334  if (icol == 6) {
1335  icol = 0;
1336  Vio_printf(sock, "\n");
1337  }
1338  }
1339  }
1340  }
1341  if (icol != 0) Vio_printf(sock, "\n");
1342 
1343  /* Close off the socket */
1344  Vio_connectFree(sock);
1345  Vio_dtor(&sock);
1346 }
1347 
1348 VPUBLIC double Vgrid_integrate(Vgrid *thee) {
1349 
1350  int i, j, k, nx, ny, nz;
1351  double sum, w;
1352 
1353  if (thee == VNULL) {
1354  Vnm_print(2, "Vgrid_integrate: Got VNULL thee!\n");
1355  VASSERT(0);
1356  }
1357 
1358  nx = thee->nx;
1359  ny = thee->ny;
1360  nz = thee->nz;
1361 
1362  sum = 0.0;
1363 
1364  for (k=0; k<nz; k++) {
1365  w = 1.0;
1366  if ((k==0) || (k==(nz-1))) w = w * 0.5;
1367  for (j=0; j<ny; j++) {
1368  w = 1.0;
1369  if ((j==0) || (j==(ny-1))) w = w * 0.5;
1370  for (i=0; i<nx; i++) {
1371  w = 1.0;
1372  if ((i==0) || (i==(nx-1))) w = w * 0.5;
1373  sum = sum + w*(thee->data[IJK(i,j,k)]);
1374  }
1375  }
1376  }
1377 
1378  sum = sum*(thee->hx)*(thee->hy)*(thee->hzed);
1379 
1380  return sum;
1381 
1382 }
1383 
1384 
1385 VPUBLIC double Vgrid_normL1(Vgrid *thee) {
1386 
1387  int i, j, k, nx, ny, nz;
1388  double sum;
1389 
1390  if (thee == VNULL) {
1391  Vnm_print(2, "Vgrid_normL1: Got VNULL thee!\n");
1392  VASSERT(0);
1393  }
1394 
1395  nx = thee->nx;
1396  ny = thee->ny;
1397  nz = thee->nz;
1398 
1399  sum = 0.0;
1400  for (k=0; k<nz; k++) {
1401  for (j=0; j<ny; j++) {
1402  for (i=0; i<nx; i++) {
1403  sum = sum + VABS(thee->data[IJK(i,j,k)]);
1404  }
1405  }
1406  }
1407 
1408  sum = sum*(thee->hx)*(thee->hy)*(thee->hzed);
1409 
1410  return sum;
1411 
1412 }
1413 
1414 VPUBLIC double Vgrid_normL2(Vgrid *thee) {
1415 
1416  int i, j, k, nx, ny, nz;
1417  double sum;
1418 
1419  if (thee == VNULL) {
1420  Vnm_print(2, "Vgrid_normL2: Got VNULL thee!\n");
1421  VASSERT(0);
1422  }
1423 
1424  nx = thee->nx;
1425  ny = thee->ny;
1426  nz = thee->nz;
1427 
1428  sum = 0.0;
1429  for (k=0; k<nz; k++) {
1430  for (j=0; j<ny; j++) {
1431  for (i=0; i<nx; i++) {
1432  sum = sum + VSQR(thee->data[IJK(i,j,k)]);
1433  }
1434  }
1435  }
1436 
1437  sum = sum*(thee->hx)*(thee->hy)*(thee->hzed);
1438 
1439  return VSQRT(sum);
1440 
1441 }
1442 
1443 VPUBLIC double Vgrid_seminormH1(Vgrid *thee) {
1444 
1445  int i, j, k, d, nx, ny, nz;
1446  double pt[3], grad[3], sum, hx, hy, hzed, xmin, ymin, zmin;
1447 
1448  if (thee == VNULL) {
1449  Vnm_print(2, "Vgrid_seminormH1: Got VNULL thee!\n");
1450  VASSERT(0);
1451  }
1452 
1453  nx = thee->nx;
1454  ny = thee->ny;
1455  nz = thee->nz;
1456  hx = thee->hx;
1457  hy = thee->hy;
1458  hzed = thee->hzed;
1459  xmin = thee->xmin;
1460  ymin = thee->ymin;
1461  zmin = thee->zmin;
1462 
1463  sum = 0.0;
1464  for (k=0; k<nz; k++) {
1465  pt[2] = k*hzed + zmin;
1466  for (j=0; j<ny; j++) {
1467  pt[1] = j*hy + ymin;
1468  for (i=0; i<nx; i++) {
1469  pt[0] = i*hx + xmin;
1470  VASSERT(Vgrid_gradient(thee, pt, grad));
1471  for (d=0; d<3; d++) sum = sum + VSQR(grad[d]);
1472  }
1473  }
1474  }
1475 
1476  sum = sum*(hx)*(hy)*(hzed);
1477 
1478  if (VABS(sum) < VSMALL) sum = 0.0;
1479  else sum = VSQRT(sum);
1480 
1481  return sum;
1482 
1483 }
1484 
1485 VPUBLIC double Vgrid_normH1(Vgrid *thee) {
1486 
1487  double sum = 0.0;
1488 
1489  if (thee == VNULL) {
1490  Vnm_print(2, "Vgrid_normH1: Got VNULL thee!\n");
1491  VASSERT(0);
1492  }
1493 
1494  sum = VSQR(Vgrid_seminormH1(thee)) + VSQR(Vgrid_normL2(thee));
1495 
1496  return VSQRT(sum);
1497 
1498 }
1499 
1500 VPUBLIC double Vgrid_normLinf(Vgrid *thee) {
1501 
1502  int i, j, k, nx, ny, nz, gotval;
1503  double sum, val;
1504 
1505  if (thee == VNULL) {
1506  Vnm_print(2, "Vgrid_normLinf: Got VNULL thee!\n");
1507  VASSERT(0);
1508  }
1509 
1510  nx = thee->nx;
1511  ny = thee->ny;
1512  nz = thee->nz;
1513 
1514  sum = 0.0;
1515  gotval = 0;
1516  for (k=0; k<nz; k++) {
1517  for (j=0; j<ny; j++) {
1518  for (i=0; i<nx; i++) {
1519  val = VABS(thee->data[IJK(i,j,k)]);
1520  if (!gotval) {
1521  gotval = 1;
1522  sum = val;
1523  }
1524  if (val > sum) sum = val;
1525  }
1526  }
1527  }
1528 
1529  return sum;
1530 
1531 }
1532 
double xmin
Definition: vgrid.h:89
VPUBLIC double Vgrid_normH1(Vgrid *thee)
Get the norm (or energy norm) of the data. This returns the integral: .
Definition: vgrid.c:1485
VPUBLIC double Vgrid_normL1(Vgrid *thee)
Get the norm of the data. This returns the integral: .
Definition: vgrid.c:1385
double ymax
Definition: vgrid.h:93
VPUBLIC double Vgrid_normLinf(Vgrid *thee)
Get the norm of the data. This returns the integral: .
Definition: vgrid.c:1500
Electrostatic potential oracle for Cartesian mesh data.
Definition: vgrid.h:81
int ny
Definition: vgrid.h:84
VPRIVATE char * MCcommChars
Comment characters for socket reads.
Definition: vparam.c:71
int nz
Definition: vgrid.h:85
double * data
Definition: vgrid.h:95
double xmax
Definition: vgrid.h:92
Potential oracle for Cartesian mesh data.
int nx
Definition: vgrid.h:83
VPUBLIC double Vgrid_integrate(Vgrid *thee)
Get the integral of the data.
Definition: vgrid.c:1348
#define VEMBED(rctag)
Allows embedding of RCS ID tags in object files.
Definition: vhal.h:563
int ctordata
Definition: vgrid.h:97
VPUBLIC int Vstring_strcasecmp(const char *s1, const char *s2)
Case-insensitive string comparison (BSD standard)
Definition: vstring.c:66
VPUBLIC int Vgrid_gradient(Vgrid *thee, double pt[3], double grad[3])
Get first derivative values at a point.
Definition: vgrid.c:368
Vmem * mem
Definition: vgrid.h:99
VPUBLIC unsigned long int Vgrid_memChk(Vgrid *thee)
Return the memory used by this structure (and its contents) in bytes.
Definition: vgrid.c:62
double zmin
Definition: vgrid.h:91
int readdata
Definition: vgrid.h:96
double ymin
Definition: vgrid.h:90
VPUBLIC int Vgrid_readDX(Vgrid *thee, const char *iodev, const char *iofmt, const char *thost, const char *fname)
Read in data in OpenDX grid format.
Definition: vgrid.c:575
VPUBLIC double Vgrid_normL2(Vgrid *thee)
Get the norm of the data. This returns the integral: .
Definition: vgrid.c:1414
VPUBLIC double Vgrid_seminormH1(Vgrid *thee)
Get the semi-norm of the data. This returns the integral: .
Definition: vgrid.c:1443
VPUBLIC int Vgrid_value(Vgrid *thee, double pt[3], double *value)
Get potential value (from mesh or approximation) at a point.
Definition: vgrid.c:170
double hy
Definition: vgrid.h:87
double hzed
Definition: vgrid.h:88
double hx
Definition: vgrid.h:86
VPRIVATE char * MCwhiteChars
Whitespace characters for socket reads.
Definition: vparam.c:65
double zmax
Definition: vgrid.h:94