APBS  1.5
vgrid.c
Go to the documentation of this file.
1 
57 #include "vgrid.h"
58 #include <stdio.h>
59 
60 VEMBED(rcsid="$Id$")
61 
62 #if !defined(VINLINE_VGRID)
63  VPUBLIC unsigned long int Vgrid_memChk(Vgrid *thee) {
64  return Vmem_bytes(thee->mem);
65  }
66 #endif
67 #define IJK(i,j,k) (((k)*(nx)*(ny))+((j)*(nx))+(i))
68 
69 #if defined(_WIN32) && (_MSC_VER < 1800)
70 #include <float.h>
71 int isnan(double d)
72 {
73  return _isnan(d);
74 }
75 #endif
76 
77 VPRIVATE char *MCwhiteChars = " =,;\t\n";
78 VPRIVATE char *MCcommChars = "#%";
79 VPRIVATE double Vcompare;
80 VPRIVATE char Vprecision[26];
81 
82 /* ///////////////////////////////////////////////////////////////////////////
83 // Routine: Vgrid_ctor
84 // Author: Nathan Baker
86 VPUBLIC Vgrid* Vgrid_ctor(int nx,
87  int ny,
88  int nz,
89  double hx,
90  double hy,
91  double hzed,
92  double xmin,
93  double ymin,
94  double zmin,
95  double *data
96  ) {
97 
98  Vgrid *thee = VNULL;
99 
100  thee = (Vgrid*)Vmem_malloc(VNULL, 1, sizeof(Vgrid));
101  VASSERT(thee != VNULL);
102  VASSERT(Vgrid_ctor2(thee, nx, ny, nz, hx, hy, hzed,
103  xmin, ymin, zmin, data));
104 
105  return thee;
106 }
107 
108 /* ///////////////////////////////////////////////////////////////////////////
109 // Routine: Vgrid_ctor2
110 // Author: Nathan Baker
112 VPUBLIC int Vgrid_ctor2(Vgrid *thee, int nx, int ny, int nz,
113  double hx, double hy, double hzed,
114  double xmin, double ymin, double zmin,
115  double *data) {
116 
117  if (thee == VNULL) return 0;
118  thee->nx = nx;
119  thee->ny = ny;
120  thee->nz = nz;
121  thee->hx = hx;
122  thee->hy = hy;
123  thee->hzed = hzed;
124  thee->xmin = xmin;
125  thee->xmax = xmin + (nx-1)*hx;
126  thee->ymin = ymin;
127  thee->ymax = ymin + (ny-1)*hy;
128  thee->zmin = zmin;
129  thee->zmax = zmin + (nz-1)*hzed;
130  if (data == VNULL) {
131  thee->ctordata = 0;
132  thee->readdata = 0;
133  } else {
134  thee->ctordata = 1;
135  thee->readdata = 0;
136  thee->data = data;
137  }
138 
139  thee->mem = Vmem_ctor("APBS:VGRID");
140 
141  Vcompare = pow(10,-1*(VGRID_DIGITS - 2));
142  sprintf(Vprecision,"%%12.%de %%12.%de %%12.%de", VGRID_DIGITS,
143  VGRID_DIGITS, VGRID_DIGITS);
144 
145  return 1;
146 }
147 
148 /* ///////////////////////////////////////////////////////////////////////////
149 // Routine: Vgrid_dtor
150 // Author: Nathan Baker
152 VPUBLIC void Vgrid_dtor(Vgrid **thee) {
153 
154  if ((*thee) != VNULL) {
155  Vgrid_dtor2(*thee);
156  Vmem_free(VNULL, 1, sizeof(Vgrid), (void **)thee);
157  (*thee) = VNULL;
158  }
159 }
160 
161 /* ///////////////////////////////////////////////////////////////////////////
162 // Routine: Vgrid_dtor2
163 // Author: Nathan Baker
165 VPUBLIC void Vgrid_dtor2(Vgrid *thee) {
166 
167  if (thee->readdata) {
168  Vmem_free(thee->mem, (thee->nx*thee->ny*thee->nz), sizeof(double),
169  (void **)&(thee->data));
170  }
171  Vmem_dtor(&(thee->mem));
172 
173 }
174 
175 /* ///////////////////////////////////////////////////////////////////////////
176 // Routine: Vgrid_value
177 // Author: Nathan Baker
179 VPUBLIC int Vgrid_value(Vgrid *thee, double pt[3], double *value) {
180 
181  int nx, ny, nz;
182  size_t ihi, jhi, khi, ilo, jlo, klo;
183  double hx, hy, hzed, xmin, ymin, zmin, ifloat, jfloat, kfloat;
184  double xmax, ymax, zmax;
185  double u, dx, dy, dz;
186 
187  if (thee == VNULL) {
188  Vnm_print(2, "Vgrid_value: Error -- got VNULL thee!\n");
189  VASSERT(0);
190  }
191  if (!(thee->ctordata || thee->readdata)) {
192  Vnm_print(2, "Vgrid_value: Error -- no data available!\n");
193  VASSERT(0);
194  }
195 
196  nx = thee->nx;
197  ny = thee->ny;
198  nz = thee->nz;
199  hx = thee->hx;
200  hy = thee->hy;
201  hzed = thee->hzed;
202  xmin = thee->xmin;
203  ymin = thee->ymin;
204  zmin = thee->zmin;
205  xmax = thee->xmax;
206  ymax = thee->ymax;
207  zmax = thee->zmax;
208 
209  u = 0;
210 
211  ifloat = (pt[0] - xmin)/hx;
212  jfloat = (pt[1] - ymin)/hy;
213  kfloat = (pt[2] - zmin)/hzed;
214 
215  ihi = (int)ceil(ifloat);
216  jhi = (int)ceil(jfloat);
217  khi = (int)ceil(kfloat);
218  ilo = (int)floor(ifloat);
219  jlo = (int)floor(jfloat);
220  klo = (int)floor(kfloat);
221  if (VABS(pt[0] - xmin) < Vcompare) ilo = 0;
222  if (VABS(pt[1] - ymin) < Vcompare) jlo = 0;
223  if (VABS(pt[2] - zmin) < Vcompare) klo = 0;
224  if (VABS(pt[0] - xmax) < Vcompare) ihi = nx-1;
225  if (VABS(pt[1] - ymax) < Vcompare) jhi = ny-1;
226  if (VABS(pt[2] - zmax) < Vcompare) khi = nz-1;
227 
228  /* See if we're on the mesh */
229  /*the condions starting with ilo>=0 seem unnecessary since they are of type size_t*/
230  if ((ihi<nx) && (jhi<ny) && (khi<nz) /*&&
231  (ilo>=0) && (jlo>=0) && (klo>=0)*/) {
232 
233  dx = ifloat - (double)(ilo);
234  dy = jfloat - (double)(jlo);
235  dz = kfloat - (double)(klo);
236  u = dx *dy *dz *(thee->data[IJK(ihi,jhi,khi)])
237  + dx *(1.0-dy)*dz *(thee->data[IJK(ihi,jlo,khi)])
238  + dx *dy *(1.0-dz)*(thee->data[IJK(ihi,jhi,klo)])
239  + dx *(1.0-dy)*(1.0-dz)*(thee->data[IJK(ihi,jlo,klo)])
240  + (1.0-dx)*dy *dz *(thee->data[IJK(ilo,jhi,khi)])
241  + (1.0-dx)*(1.0-dy)*dz *(thee->data[IJK(ilo,jlo,khi)])
242  + (1.0-dx)*dy *(1.0-dz)*(thee->data[IJK(ilo,jhi,klo)])
243  + (1.0-dx)*(1.0-dy)*(1.0-dz)*(thee->data[IJK(ilo,jlo,klo)]);
244 
245  *value = u;
246 
247  if (isnan(u)) {
248  Vnm_print(2, "Vgrid_value: Got NaN!\n");
249  Vnm_print(2, "Vgrid_value: (x, y, z) = (%4.3f, %4.3f, %4.3f)\n",
250  pt[0], pt[1], pt[2]);
251  Vnm_print(2, "Vgrid_value: (ihi, jhi, khi) = (%d, %d, %d)\n",
252  ihi, jhi, khi);
253  Vnm_print(2, "Vgrid_value: (ilo, jlo, klo) = (%d, %d, %d)\n",
254  ilo, jlo, klo);
255  Vnm_print(2, "Vgrid_value: (nx, ny, nz) = (%d, %d, %d)\n",
256  nx, ny, nz);
257  Vnm_print(2, "Vgrid_value: (dx, dy, dz) = (%4.3f, %4.3f, %4.3f)\n",
258  dx, dy, dz);
259  Vnm_print(2, "Vgrid_value: data[IJK(ihi,jhi,khi)] = %g\n",
260  thee->data[IJK(ihi,jhi,khi)]);
261  Vnm_print(2, "Vgrid_value: data[IJK(ihi,jlo,khi)] = %g\n",
262  thee->data[IJK(ihi,jlo,khi)]);
263  Vnm_print(2, "Vgrid_value: data[IJK(ihi,jhi,klo)] = %g\n",
264  thee->data[IJK(ihi,jhi,klo)]);
265  Vnm_print(2, "Vgrid_value: data[IJK(ihi,jlo,klo)] = %g\n",
266  thee->data[IJK(ihi,jlo,klo)]);
267  Vnm_print(2, "Vgrid_value: data[IJK(ilo,jhi,khi)] = %g\n",
268  thee->data[IJK(ilo,jhi,khi)]);
269  Vnm_print(2, "Vgrid_value: data[IJK(ilo,jlo,khi)] = %g\n",
270  thee->data[IJK(ilo,jlo,khi)]);
271  Vnm_print(2, "Vgrid_value: data[IJK(ilo,jhi,klo)] = %g\n",
272  thee->data[IJK(ilo,jhi,klo)]);
273  Vnm_print(2, "Vgrid_value: data[IJK(ilo,jlo,klo)] = %g\n",
274  thee->data[IJK(ilo,jlo,klo)]);
275  }
276  return 1;
277 
278  } else {
279 
280  *value = 0;
281  return 0;
282 
283  }
284 
285  return 0;
286 
287 }
288 
289 /* ///////////////////////////////////////////////////////////////////////////
290 // Routine: Vgrid_curvature
291 //
292 // Notes: cflag=0 ==> Reduced Maximal Curvature
293 // cflag=1 ==> Mean Curvature (Laplace)
294 // cflag=2 ==> Gauss Curvature
295 // cflag=3 ==> True Maximal Curvature
296 //
297 // Authors: Stephen Bond and Nathan Baker
299 VPUBLIC int Vgrid_curvature(Vgrid *thee, double pt[3], int cflag,
300  double *value) {
301 
302  double hx, hy, hzed, curv;
303  double dxx, dyy, dzz;
304  double uleft, umid, uright, testpt[3];
305 
306  if (thee == VNULL) {
307  Vnm_print(2, "Vgrid_curvature: Error -- got VNULL thee!\n");
308  VASSERT(0);
309  }
310  if (!(thee->ctordata || thee->readdata)) {
311  Vnm_print(2, "Vgrid_curvature: Error -- no data available!\n");
312  VASSERT(0);
313  }
314 
315  hx = thee->hx;
316  hy = thee->hy;
317  hzed = thee->hzed;
318 
319  curv = 0.0;
320 
321  testpt[0] = pt[0];
322  testpt[1] = pt[1];
323  testpt[2] = pt[2];
324 
325  /* Compute 2nd derivative in the x-direction */
326  VJMPERR1(Vgrid_value( thee, testpt, &umid));
327  testpt[0] = pt[0] - hx;
328  VJMPERR1(Vgrid_value( thee, testpt, &uleft));
329  testpt[0] = pt[0] + hx;
330  VJMPERR1(Vgrid_value( thee, testpt, &uright));
331  testpt[0] = pt[0];
332 
333  dxx = (uright - 2*umid + uleft)/(hx*hx);
334 
335  /* Compute 2nd derivative in the y-direction */
336  VJMPERR1(Vgrid_value( thee, testpt, &umid));
337  testpt[1] = pt[1] - hy;
338  VJMPERR1(Vgrid_value( thee, testpt, &uleft));
339  testpt[1] = pt[1] + hy;
340  VJMPERR1(Vgrid_value( thee, testpt, &uright));
341  testpt[1] = pt[1];
342 
343  dyy = (uright - 2*umid + uleft)/(hy*hy);
344 
345  /* Compute 2nd derivative in the z-direction */
346  VJMPERR1(Vgrid_value( thee, testpt, &umid));
347  testpt[2] = pt[2] - hzed;
348  VJMPERR1(Vgrid_value( thee, testpt, &uleft));
349  testpt[2] = pt[2] + hzed;
350  VJMPERR1(Vgrid_value( thee, testpt, &uright));
351 
352  dzz = (uright - 2*umid + uleft)/(hzed*hzed);
353 
354 
355  if ( cflag == 0 ) {
356  curv = fabs(dxx);
357  curv = ( curv > fabs(dyy) ) ? curv : fabs(dyy);
358  curv = ( curv > fabs(dzz) ) ? curv : fabs(dzz);
359  } else if ( cflag == 1 ) {
360  curv = (dxx + dyy + dzz)/3.0;
361  } else {
362  Vnm_print(2, "Vgrid_curvature: support for cflag = %d not available!\n", cflag);
363  VASSERT( 0 ); /* Feature Not Coded Yet! */
364  }
365 
366  *value = curv;
367  return 1;
368 
369  VERROR1:
370  return 0;
371 
372 }
373 
374 /* ///////////////////////////////////////////////////////////////////////////
375 // Routine: Vgrid_gradient
376 //
377 // Authors: Nathan Baker and Stephen Bond
379 VPUBLIC int Vgrid_gradient(Vgrid *thee, double pt[3], double grad[3]) {
380 
381  double hx, hy, hzed;
382  double uleft, umid, uright, testpt[3];
383  int haveleft, haveright;
384 
385  if (thee == VNULL) {
386  Vnm_print(2, "Vgrid_gradient: Error -- got VNULL thee!\n");
387  VASSERT(0);
388  }
389  if (!(thee->ctordata || thee->readdata)) {
390  Vnm_print(2, "Vgrid_gradient: Error -- no data available!\n");
391  VASSERT(0);
392  }
393 
394  hx = thee->hx;
395  hy = thee->hy;
396  hzed = thee->hzed;
397 
398  /* Compute derivative in the x-direction */
399  testpt[0] = pt[0];
400  testpt[1] = pt[1];
401  testpt[2] = pt[2];
402  VJMPERR1( Vgrid_value( thee, testpt, &umid));
403  testpt[0] = pt[0] - hx;
404  if (Vgrid_value( thee, testpt, &uleft)) haveleft = 1;
405  else haveleft = 0;
406  testpt[0] = pt[0] + hx;
407  if (Vgrid_value( thee, testpt, &uright)) haveright = 1;
408  else haveright = 0;
409  if (haveright && haveleft) grad[0] = (uright - uleft)/(2*hx);
410  else if (haveright) grad[0] = (uright - umid)/hx;
411  else if (haveleft) grad[0] = (umid - uleft)/hx;
412  else VJMPERR1(0);
413 
414  /* Compute derivative in the y-direction */
415  testpt[0] = pt[0];
416  testpt[1] = pt[1];
417  testpt[2] = pt[2];
418  VJMPERR1(Vgrid_value(thee, testpt, &umid));
419  testpt[1] = pt[1] - hy;
420  if (Vgrid_value( thee, testpt, &uleft)) haveleft = 1;
421  else haveleft = 0;
422  testpt[1] = pt[1] + hy;
423  if (Vgrid_value( thee, testpt, &uright)) haveright = 1;
424  else haveright = 0;
425  if (haveright && haveleft) grad[1] = (uright - uleft)/(2*hy);
426  else if (haveright) grad[1] = (uright - umid)/hy;
427  else if (haveleft) grad[1] = (umid - uleft)/hy;
428  else VJMPERR1(0);
429 
430  /* Compute derivative in the z-direction */
431  testpt[0] = pt[0];
432  testpt[1] = pt[1];
433  testpt[2] = pt[2];
434  VJMPERR1(Vgrid_value(thee, testpt, &umid));
435  testpt[2] = pt[2] - hzed;
436  if (Vgrid_value( thee, testpt, &uleft)) haveleft = 1;
437  else haveleft = 0;
438  testpt[2] = pt[2] + hzed;
439  if (Vgrid_value( thee, testpt, &uright)) haveright = 1;
440  else haveright = 0;
441  if (haveright && haveleft) grad[2] = (uright - uleft)/(2*hzed);
442  else if (haveright) grad[2] = (uright - umid)/hzed;
443  else if (haveleft) grad[2] = (umid - uleft)/hzed;
444  else VJMPERR1(0);
445 
446  return 1;
447 
448  VERROR1:
449  return 0;
450 
451 }
452 
453 /* ///////////////////////////////////////////////////////////////////////////
454  // Routine: Vgrid_readGZ
455  //
456  // Author: David Gohara
458 #ifdef HAVE_ZLIB
459 #define off_t long
460 #include "zlib.h"
461 #endif
462 VPUBLIC int Vgrid_readGZ(Vgrid *thee, const char *fname) {
463 
464 #ifdef HAVE_ZLIB
465  size_t i, j, k, u;
466  size_t len; // Temporary counter variable for loop conditionals
467  size_t header, incr;
468  double *temp;
469  double dtmp1, dtmp2, dtmp3;
470  gzFile infile;
471  char line[VMAX_ARGLEN];
472 
473  header = 0;
474 
475  /* Check to see if the existing data is null and, if not, clear it out */
476  if (thee->data != VNULL) {
477  Vnm_print(1, "%s: destroying existing data!\n", __func__);
478  Vmem_free(thee->mem, thee->nx * thee->ny * thee->nz, sizeof(double),
479  (void **)&(thee->data));
480  }
481 
482  thee->readdata = 1;
483  thee->ctordata = 0;
484 
485  infile = gzopen(fname, "rb");
486  if (infile == Z_NULL) {
487  Vnm_print(2, "%s: Problem opening compressed file %s\n", __func__, fname);
488  return VRC_FAILURE;
489  }
490 
491  thee->hx = 0.0;
492  thee->hy = 0.0;
493  thee->hzed = 0.0;
494 
495  //read data here
496  while (header < 7) {
497  if(gzgets(infile, line, VMAX_ARGLEN) == Z_NULL){
498  return VRC_FAILURE;
499  }
500 
501  // Skip comments and newlines
502  if(strncmp(line, "#", 1) == 0) continue;
503  if(line[0] == '\n') continue;
504 
505  switch (header) {
506  case 0:
507  sscanf(line, "object 1 class gridpositions counts %d %d %d",
508  &(thee->nx),&(thee->ny),&(thee->nz));
509  break;
510  case 1:
511  sscanf(line, "origin %lf %lf %lf",
512  &(thee->xmin),&(thee->ymin),&(thee->zmin));
513  break;
514  case 2:
515  case 3:
516  case 4:
517  sscanf(line, "delta %lf %lf %lf",&dtmp1,&dtmp2,&dtmp3);
518  thee->hx += dtmp1;
519  thee->hy += dtmp2;
520  thee->hzed += dtmp3;
521  break;
522  default:
523  break;
524  }
525 
526  header++;
527  }
528 
529  /* Allocate space for the data */
530  Vnm_print(0, "%s: allocating %d x %d x %d doubles for storage\n",
531  __func__, thee->nx, thee->ny, thee->nz);
532  len = thee->nx * thee->ny * thee->nz;
533 
534  thee->data = VNULL;
535  thee->data = Vmem_malloc(thee->mem, len, sizeof(double));
536  if (thee->data == VNULL) {
537  Vnm_print(2, "%s: Unable to allocate space for data!\n", __func__);
538  return 0;
539  }
540 
541  /* Allocate a temporary buffer to store the compressed
542  * data into (column major order). Add 2 to ensure the buffer is
543  * big enough to take extra data on the final read loop.
544  */
545  temp = (double *)malloc(len * (2 * sizeof(double)));
546 
547  for (i = 0; i < len; i += 3){
548  memset(&line, 0, sizeof(line));
549  gzgets(infile, line, VMAX_ARGLEN);
550  sscanf(line, "%lf %lf %lf", &temp[i], &temp[i+1], &temp[i+2]);
551  }
552 
553  /* Now move the data to row major order */
554  incr = 0;
555  for (i=0; i<thee->nx; i++) {
556  for (j=0; j<thee->ny; j++) {
557  for (k=0; k<thee->nz; k++) {
558  u = k*(thee->nx)*(thee->ny)+j*(thee->nx)+i;
559  (thee->data)[u] = temp[incr++];
560  }
561  }
562  }
563 
564  /* calculate grid maxima */
565  thee->xmax = thee->xmin + (thee->nx-1)*thee->hx;
566  thee->ymax = thee->ymin + (thee->ny-1)*thee->hy;
567  thee->zmax = thee->zmin + (thee->nz-1)*thee->hzed;
568 
569  /* Close off the socket */
570  gzclose(infile);
571  free(temp);
572 #else
573 
574  Vnm_print(0, "WARNING\n");
575  Vnm_print(0, "Vgrid_readGZ: gzip read/write support is disabled in this build\n");
576  Vnm_print(0, "Vgrid_readGZ: configure and compile without the --disable-zlib flag.\n");
577  Vnm_print(0, "WARNING\n");
578 #endif
579  return VRC_SUCCESS;
580 }
581 
586 VPUBLIC int Vgrid_readDX(Vgrid *thee,
587  const char *iodev,
588  const char *iofmt,
589  const char *thost,
590  const char *fname
591  ) {
592 
593  size_t i, j, k, itmp, u;
594  double dtmp;
595  char tok[VMAX_BUFSIZE];
596  Vio *sock;
597 
598  /* Check to see if the existing data is null and, if not, clear it out */
599  if (thee->data != VNULL) {
600  Vnm_print(1, "Vgrid_readDX: destroying existing data!\n");
601  Vmem_free(thee->mem, (thee->nx*thee->ny*thee->nz), sizeof(double),
602  (void **)&(thee->data)); }
603  thee->readdata = 1;
604  thee->ctordata = 0;
605 
606  /* Set up the virtual socket */
607  sock = Vio_ctor(iodev,iofmt,thost,fname,"r");
608  if (sock == VNULL) {
609  Vnm_print(2, "Vgrid_readDX: Problem opening virtual socket %s\n",
610  fname);
611  return 0;
612  }
613  if (Vio_accept(sock, 0) < 0) {
614  Vnm_print(2, "Vgrid_readDX: Problem accepting virtual socket %s\n",
615  fname);
616  return 0;
617  }
618 
619  Vio_setWhiteChars(sock, MCwhiteChars);
620  Vio_setCommChars(sock, MCcommChars);
621 
622  /* Read in the DX regular positions */
623  /* Get "object" */
624  VJMPERR2(1 == Vio_scanf(sock, "%s", tok));
625  VJMPERR1(!strcmp(tok, "object"));
626  /* Get "1" */
627  VJMPERR2(1 == Vio_scanf(sock, "%d", &itmp));
628  /* Get "class" */
629  VJMPERR2(1 == Vio_scanf(sock, "%s", tok));
630  VJMPERR1(!strcmp(tok, "class"));
631  /* Get "gridpositions" */
632  VJMPERR2(1 == Vio_scanf(sock, "%s", tok));
633  VJMPERR1(!strcmp(tok, "gridpositions"));
634  /* Get "counts" */
635  VJMPERR2(1 == Vio_scanf(sock, "%s", tok));
636  VJMPERR1(!strcmp(tok, "counts"));
637  /* Get nx */
638  VJMPERR2(1 == Vio_scanf(sock, "%s", tok));
639  VJMPERR1(1 == sscanf(tok, "%d", &(thee->nx)));
640  /* Get ny */
641  VJMPERR2(1 == Vio_scanf(sock, "%s", tok));
642  VJMPERR1(1 == sscanf(tok, "%d", &(thee->ny)));
643  /* Get nz */
644  VJMPERR2(1 == Vio_scanf(sock, "%s", tok));
645  VJMPERR1(1 == sscanf(tok, "%d", &(thee->nz)));
646  Vnm_print(0, "Vgrid_readDX: Grid dimensions %d x %d x %d grid\n",
647  thee->nx, thee->ny, thee->nz);
648  /* Get "origin" */
649  VJMPERR2(1 == Vio_scanf(sock, "%s", tok));
650  VJMPERR1(!strcmp(tok, "origin"));
651  /* Get xmin */
652  VJMPERR2(1 == Vio_scanf(sock, "%s", tok));
653  VJMPERR1(1 == sscanf(tok, "%lf", &(thee->xmin)));
654  /* Get ymin */
655  VJMPERR2(1 == Vio_scanf(sock, "%s", tok));
656  VJMPERR1(1 == sscanf(tok, "%lf", &(thee->ymin)));
657  /* Get zmin */
658  VJMPERR2(1 == Vio_scanf(sock, "%s", tok));
659  VJMPERR1(1 == sscanf(tok, "%lf", &(thee->zmin)));
660  Vnm_print(0, "Vgrid_readDX: Grid origin = (%g, %g, %g)\n",
661  thee->xmin, thee->ymin, thee->zmin);
662  /* Get "delta" */
663  VJMPERR2(1 == Vio_scanf(sock, "%s", tok));
664  VJMPERR1(!strcmp(tok, "delta"));
665  /* Get hx */
666  VJMPERR2(1 == Vio_scanf(sock, "%s", tok));
667  VJMPERR1(1 == sscanf(tok, "%lf", &(thee->hx)));
668  /* Get 0.0 */
669  VJMPERR2(1 == Vio_scanf(sock, "%s", tok));
670  VJMPERR1(1 == sscanf(tok, "%lf", &dtmp));
671  VJMPERR1(dtmp == 0.0);
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 "delta" */
677  VJMPERR2(1 == Vio_scanf(sock, "%s", tok));
678  VJMPERR1(!strcmp(tok, "delta"));
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 hy */
684  VJMPERR2(1 == Vio_scanf(sock, "%s", tok));
685  VJMPERR1(1 == sscanf(tok, "%lf", &(thee->hy)));
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 "delta" */
691  VJMPERR2(1 == Vio_scanf(sock, "%s", tok));
692  VJMPERR1(!strcmp(tok, "delta"));
693  /* Get 0.0 */
694  VJMPERR2(1 == Vio_scanf(sock, "%s", tok));
695  VJMPERR1(1 == sscanf(tok, "%lf", &dtmp));
696  VJMPERR1(dtmp == 0.0);
697  /* Get 0.0 */
698  VJMPERR2(1 == Vio_scanf(sock, "%s", tok));
699  VJMPERR1(1 == sscanf(tok, "%lf", &dtmp));
700  VJMPERR1(dtmp == 0.0);
701  /* Get hz */
702  VJMPERR2(1 == Vio_scanf(sock, "%s", tok));
703  VJMPERR1(1 == sscanf(tok, "%lf", &(thee->hzed)));
704  Vnm_print(0, "Vgrid_readDX: Grid spacings = (%g, %g, %g)\n",
705  thee->hx, thee->hy, thee->hzed);
706  /* Get "object" */
707  VJMPERR2(1 == Vio_scanf(sock, "%s", tok));
708  VJMPERR1(!strcmp(tok, "object"));
709  /* Get "2" */
710  VJMPERR2(1 == Vio_scanf(sock, "%s", tok));
711  /* Get "class" */
712  VJMPERR2(1 == Vio_scanf(sock, "%s", tok));
713  VJMPERR1(!strcmp(tok, "class"));
714  /* Get "gridconnections" */
715  VJMPERR2(1 == Vio_scanf(sock, "%s", tok));
716  VJMPERR1(!strcmp(tok, "gridconnections"));
717  /* Get "counts" */
718  VJMPERR2(1 == Vio_scanf(sock, "%s", tok));
719  VJMPERR1(!strcmp(tok, "counts"));
720  /* Get the dimensions again */
721  VJMPERR2(1 == Vio_scanf(sock, "%s", tok));
722  VJMPERR2(1 == Vio_scanf(sock, "%s", tok));
723  VJMPERR2(1 == Vio_scanf(sock, "%s", tok));
724  /* Get "object" */
725  VJMPERR2(1 == Vio_scanf(sock, "%s", tok));
726  VJMPERR1(!strcmp(tok, "object"));
727  /* Get # */
728  VJMPERR2(1 == Vio_scanf(sock, "%s", tok));
729  /* Get "class" */
730  VJMPERR2(1 == Vio_scanf(sock, "%s", tok));
731  VJMPERR1(!strcmp(tok, "class"));
732  /* Get "array" */
733  VJMPERR2(1 == Vio_scanf(sock, "%s", tok));
734  VJMPERR1(!strcmp(tok, "array"));
735  /* Get "type" */
736  VJMPERR2(1 == Vio_scanf(sock, "%s", tok));
737  VJMPERR1(!strcmp(tok, "type"));
738  /* Get "double" */
739  VJMPERR2(1 == Vio_scanf(sock, "%s", tok));
740  VJMPERR1(!strcmp(tok, "double"));
741  /* Get "rank" */
742  VJMPERR2(1 == Vio_scanf(sock, "%s", tok));
743  VJMPERR1(!strcmp(tok, "rank"));
744  /* Get # */
745  VJMPERR2(1 == Vio_scanf(sock, "%s", tok));
746  /* Get "items" */
747  VJMPERR2(1 == Vio_scanf(sock, "%s", tok));
748  VJMPERR1(!strcmp(tok, "items"));
749  /* Get # */
750  VJMPERR2(1 == Vio_scanf(sock, "%s", tok));
751  VJMPERR1(1 == sscanf(tok, "%lu", &itmp));
752  u = (size_t)thee->nx * thee->ny * thee->nz;
753  VJMPERR1(u == itmp);
754  /* Get "data" */
755  VJMPERR2(1 == Vio_scanf(sock, "%s", tok));
756  VJMPERR1(!strcmp(tok, "data"));
757  /* Get "follows" */
758  VJMPERR2(1 == Vio_scanf(sock, "%s", tok));
759  VJMPERR1(!strcmp(tok, "follows"));
760 
761  /* Allocate space for the data */
762  Vnm_print(0, "Vgrid_readDX: allocating %d x %d x %d doubles for storage\n",
763  thee->nx, thee->ny, thee->nz);
764  thee->data = VNULL;
765  thee->data = (double*)Vmem_malloc(thee->mem, u, sizeof(double));
766  if (thee->data == VNULL) {
767  Vnm_print(2, "Vgrid_readDX: Unable to allocate space for data!\n");
768  return 0;
769  }
770 
771  for (i=0; i<thee->nx; i++) {
772  for (j=0; j<thee->ny; j++) {
773  for (k=0; k<thee->nz; k++) {
774  u = k*(thee->nx)*(thee->ny)+j*(thee->nx)+i;
775  VJMPERR2(1 == Vio_scanf(sock, "%s", tok));
776  VJMPERR1(1 == sscanf(tok, "%lf", &dtmp));
777  (thee->data)[u] = dtmp;
778  }
779  }
780  }
781 
782  /* calculate grid maxima */
783  thee->xmax = thee->xmin + (thee->nx-1)*thee->hx;
784  thee->ymax = thee->ymin + (thee->ny-1)*thee->hy;
785  thee->zmax = thee->zmin + (thee->nz-1)*thee->hzed;
786 
787  /* Close off the socket */
788  Vio_acceptFree(sock);
789  Vio_dtor(&sock);
790 
791  return 1;
792 
793  VERROR1:
794  Vio_dtor(&sock);
795  Vnm_print(2, "Vgrid_readDX: Format problem with input file <%s>\n",
796  fname);
797  return 0;
798 
799  VERROR2:
800  Vio_dtor(&sock);
801  Vnm_print(2, "Vgrid_readDX: I/O problem with input file <%s>\n",
802  fname);
803  return 0;
804 }
805 
810 VPUBLIC int Vgrid_readDXBIN(Vgrid *thee, const char *iodev, const char *iofmt,
811  const char *thost, const char *fname) {
812 
813  size_t i, j, k, itmp, u;
814  double dtmp, dtmp2;
815  char tok[VMAX_BUFSIZE];
816  int isBinary = 0;
817  //Vio *sock;
818 
819  /* Check to see if the existing data is null and, if not, clear it out */
820  if (thee->data != VNULL) {
821  Vnm_print(1, "Vgrid_readDXBIN: destroying existing data!\n");
822  Vmem_free(thee->mem, (thee->nx*thee->ny*thee->nz), sizeof(double),
823  (void **)&(thee->data)); }
824  thee->readdata = 1;
825  thee->ctordata = 0;
826 
827  /*Opend file fd for binary reading*/
828  FILE *fd = fopen(fname,"rb");
829  if(fd == NULL){
830  printf("Vgrid_readDXBIN: Problem opening file %s\n",fname);
831  fclose(fd);
832  return 0;
833  }
834 
835  /* Set up the virtual socket */
836 // sock = Vio_ctor(iodev,iofmt,thost,fname,"r");
837 // if (sock == VNULL) {
838 // Vnm_print(2, "Vgrid_readDX: Problem opening virtual socket %s\n",
839 // fname);
840 // return 0;
841 // }
842 // if (Vio_accept(sock, 0) < 0) {
843 // Vnm_print(2, "Vgrid_readDX: Problem accepting virtual socket %s\n",
844 // fname);
845 // return 0;
846 // }
847 //
848 // Vio_setWhiteChars(sock, MCwhiteChars);
849 // Vio_setCommChars(sock, MCcommChars);
850 
851  //skip comments
852  do{
853  fgets(tok, VMAX_BUFSIZE, fd);
854  }
855  while(tok[0]=='#');
856 
857  //get counts
858  if(sscanf(tok,"object 1 class gridpositions counts %i %i %i\n",&(thee->nx),&(thee->ny),&(thee->nz))!= 3){
859  printf("Vgrid_readDXBIN: Failed to read dimensions.\n");
860  fclose(fd);
861  return 0;
862  }
863  printf("Vgrid_readDXBIN: Grid dimensions %d x %d x %d grid\n",thee->nx, thee->ny, thee->nz);
864 
865  /* Read in the DX regular positions */
866 
867  /* Get "object" */
868 // VJMPERR2(1 == Vio_scanf(sock, "%s", tok));
869 // VJMPERR1(!strcmp(tok, "object"));
870 // /* Get "1" */
871 // VJMPERR2(1 == Vio_scanf(sock, "%d", &itmp));
872 // /* Get "class" */
873 // VJMPERR2(1 == Vio_scanf(sock, "%s", tok));
874 // VJMPERR1(!strcmp(tok, "class"));
875 // /* Get "gridpositions" */
876 // VJMPERR2(1 == Vio_scanf(sock, "%s", tok));
877 // VJMPERR1(!strcmp(tok, "gridpositions"));
878 // /* Get "counts" */
879 // VJMPERR2(1 == Vio_scanf(sock, "%s", tok));
880 // VJMPERR1(!strcmp(tok, "counts"));
881 // /* Get nx */
882 // VJMPERR2(1 == Vio_scanf(sock, "%s", tok));
883 // VJMPERR1(1 == sscanf(tok, "%d", &(thee->nx)));
884 // /* Get ny */
885 // VJMPERR2(1 == Vio_scanf(sock, "%s", tok));
886 // VJMPERR1(1 == sscanf(tok, "%d", &(thee->ny)));
887 // /* Get nz */
888 // VJMPERR2(1 == Vio_scanf(sock, "%s", tok));
889 // VJMPERR1(1 == sscanf(tok, "%d", &(thee->nz)));
890 // Vnm_print(0, "Vgrid_readDX: Grid dimensions %d x %d x %d grid\n",
891 // thee->nx, thee->ny, thee->nz);
892 
893  if(fgets(tok,VMAX_BUFSIZE, fd) == NULL){
894  printf("Vgrid_readDXBIN: unexpected end of file.\n");
895  fclose(fd);
896  return 0;
897  }
898  if(sscanf(tok,"origin %lf %lf %lf",&(thee->xmin),&(thee->ymin),&(thee->zmin))!=3){
899  printf("Vgrid_readDXBIN: Failed to read origin cell data.\n");
900  fclose(fd);
901  return 0;
902  }
903  printf("Vgrid_readDXBIN: Grid origin = (%g %g %g)\n",thee->xmin, thee->ymin, thee->zmin);
904 
905  //get Delta x
906  if(fgets(tok,VMAX_BUFSIZE, fd) == NULL){
907  printf("Vgrid_readDXBIN: unexpected end of file.\n");
908  fclose(fd);
909  return 0;
910  }
911  if(sscanf(tok,"delta %lf %lf %lf",&(thee->hx),&dtmp,&dtmp2)!=3){
912  printf("Vgrid_readDXBIN: Failed to read delta x data.\n");
913  fclose(fd);
914  return 0;
915  }
916  //get Delta y
917  if(fgets(tok,VMAX_BUFSIZE, fd) == NULL){
918  printf("Vgrid_readDXBIN: Unexpected end of file\n");
919  fclose(fd);
920  return 0;
921  }
922  if(sscanf(tok,"delta %lf %lf %lf",&dtmp,&(thee->hy),&dtmp2)!=3){
923  printf("Vgrid_readDXBIN: Failed to read delta y data.\n");
924  fclose(fd);
925  return 0;
926  }
927  //get Delta z
928  if(fgets(tok,VMAX_BUFSIZE, fd) == NULL){
929  printf("Vgrid_readDXBIN: Unexpected end of file.\n");
930  fclose(fd);
931  return 0;
932  }
933  if(sscanf(tok,"delta %lf %lf %lf",&dtmp,&dtmp2,&(thee->hzed))!=3){
934  printf("Vgrid_readDXBIN: Failed to read delta z data.\n");
935  fclose(fd);
936  return 0;
937  }
938  printf("Vgrid_readDXBIN: Grid spacings = (%g, %g, %g)\n",thee->hx, thee->hy, thee->hzed);
939 
940  //skip a line
941  if(fgets(tok,VMAX_BUFSIZE, fd) == NULL){
942  printf("Vgrid_readDXBIN: Unexpected end of file.\n");
943  fclose(fd);
944  return 0;
945  }
946 
947  //scan the buffer for the word binary
948  if(fgets(tok,VMAX_BUFSIZE, fd) == NULL){
949  printf("Vgrid_readDXBIN: Unexpected end of file.\n");
950  fclose(fd);
951  return 0;
952  }
953  if(strstr(tok,"binary")){
954  isBinary = 1;
955  }
956  else{
957  printf("Vgrid_readDXBIN: Binary tag not found. Will continue to try to read binary data.");
958  }
959 
960  u = (size_t)thee->nx * thee->ny * thee->nz;
961  int tot = thee->nx * thee->ny * thee->nz;
962 
963  /*Allocate space for the data*/
964  printf("Vgrid_readDXBIN: allocating %d x %d x %d doubled for storage\n", thee->nx, thee->ny, thee->nz);
965  thee->data = NULL;
966  thee->data = (double *)malloc(tot*sizeof(double));
967 
968  if(thee->data == NULL){
969  printf("Vgrid_readDXBIN: Unable to allocate space for data!\n");
970  fclose(fd);
971  return 0;
972  }
973 
974  int counter = 0, r;
975 
976  for (i=0; i<thee->nx; i++) {
977  for (j=0; j<thee->ny; j++) {
978  for (k=0; k<thee->nz; k++) {
979  u = k*(thee->nx)*(thee->ny)+j*(thee->nx)+i;
980  r = fread(&dtmp,sizeof(double),1,fd);
981  (thee->data)[u] = dtmp;
982  if(r!= 1){
983  printf("Vgrid_readDXBIN: Failed to read doubles.\n");
984  return 0;
985  }
986  counter++;
987  }
988  }
989  }
990 
991  if(counter!=tot){
992  printf("Vgrid_readDXBIN: Read double = %d not equal to items = %d\n",counter, tot);
993  }
994 
995  /* calculate grid maxima */
996  thee->xmax = thee->xmin + (thee->nx-1)*thee->hx;
997  thee->ymax = thee->ymin + (thee->ny-1)*thee->hy;
998  thee->zmax = thee->zmin + (thee->nz-1)*thee->hzed;
999 
1000  fclose(fd);
1001 
1002  return 1;
1003 }
1004 
1005 
1006 /* ///////////////////////////////////////////////////////////////////////////
1007  // Routine: Vgrid_writeGZ
1008  //
1009  // Author: Nathan Baker
1011 VPUBLIC void Vgrid_writeGZ(Vgrid *thee, const char *iodev, const char *iofmt,
1012  const char *thost, const char *fname, char *title, double *pvec) {
1013 
1014 #ifdef HAVE_ZLIB
1015  double xmin, ymin, zmin, hx, hy, hzed;
1016 
1017  int nx, ny, nz, nxPART, nyPART, nzPART;
1018  int usepart, gotit;
1019  size_t icol, i, j, k, u;
1020  double x, y, z, xminPART, yminPART, zminPART;
1021 
1022  size_t txyz;
1023  double txmin, tymin, tzmin;
1024 
1025  char header[8196];
1026  char footer[8196];
1027  char line[80];
1028  char newline[] = "\n";
1029  gzFile outfile;
1030  char precFormat[VMAX_BUFSIZE];
1031 
1032  if (thee == VNULL) {
1033  Vnm_print(2, "Vgrid_writeGZ: Error -- got VNULL thee!\n");
1034  VASSERT(0);
1035  }
1036  if (!(thee->ctordata || thee->readdata)) {
1037  Vnm_print(2, "Vgrid_writeGZ: Error -- no data available!\n");
1038  VASSERT(0);
1039  }
1040 
1041  hx = thee->hx;
1042  hy = thee->hy;
1043  hzed = thee->hzed;
1044  nx = thee->nx;
1045  ny = thee->ny;
1046  nz = thee->nz;
1047  xmin = thee->xmin;
1048  ymin = thee->ymin;
1049  zmin = thee->zmin;
1050 
1051  if (pvec == VNULL) usepart = 0;
1052  else usepart = 1;
1053 
1054  /* Set up the virtual socket */
1055  Vnm_print(0, "Vgrid_writeGZ: Opening file...\n");
1056  outfile = gzopen(fname, "wb");
1057 
1058  if (usepart) {
1059  /* Get the lower corner and number of grid points for the local
1060  * partition */
1061  xminPART = VLARGE;
1062  yminPART = VLARGE;
1063  zminPART = VLARGE;
1064  nxPART = 0;
1065  nyPART = 0;
1066  nzPART = 0;
1067  /* First, search for the lower corner */
1068  for (k=0; k<nz; k++) {
1069  z = k*hzed + zmin;
1070  for (j=0; j<ny; j++) {
1071  y = j*hy + ymin;
1072  for (i=0; i<nx; i++) {
1073  x = i*hx + xmin;
1074  if (pvec[IJK(i,j,k)] > 0.0) {
1075  if (x < xminPART) xminPART = x;
1076  if (y < yminPART) yminPART = y;
1077  if (z < zminPART) zminPART = z;
1078  }
1079  }
1080  }
1081  }
1082  /* Now search for the number of grid points in the z direction */
1083  for (k=0; k<nz; k++) {
1084  gotit = 0;
1085  for (j=0; j<ny; j++) {
1086  for (i=0; i<nx; i++) {
1087  if (pvec[IJK(i,j,k)] > 0.0) {
1088  gotit = 1;
1089  break;
1090  }
1091  }
1092  if (gotit) break;
1093  }
1094  if (gotit) nzPART++;
1095  }
1096  /* Now search for the number of grid points in the y direction */
1097  for (j=0; j<ny; j++) {
1098  gotit = 0;
1099  for (k=0; k<nz; k++) {
1100  for (i=0; i<nx; i++) {
1101  if (pvec[IJK(i,j,k)] > 0.0) {
1102  gotit = 1;
1103  break;
1104  }
1105  }
1106  if (gotit) break;
1107  }
1108  if (gotit) nyPART++;
1109  }
1110  /* Now search for the number of grid points in the x direction */
1111  for (i=0; i<nx; i++) {
1112  gotit = 0;
1113  for (k=0; k<nz; k++) {
1114  for (j=0; j<ny; j++) {
1115  if (pvec[IJK(i,j,k)] > 0.0) {
1116  gotit = 1;
1117  break;
1118  }
1119  }
1120  if (gotit) break;
1121  }
1122  if (gotit) nxPART++;
1123  }
1124 
1125  if ((nxPART != nx) || (nyPART != ny) || (nzPART != nz)) {
1126  Vnm_print(0, "Vgrid_writeGZ: printing only subset of domain\n");
1127  }
1128 
1129  txyz = (nxPART*nyPART*nzPART);
1130  txmin = xminPART;
1131  tymin = yminPART;
1132  tzmin = zminPART;
1133 
1134  }else {
1135 
1136  txyz = (nx*ny*nz);
1137  txmin = xmin;
1138  tymin = ymin;
1139  tzmin = zmin;
1140 
1141  }
1142 
1143  /* Write off the title (if we're not XDR) */
1144  sprintf(header,
1145  "# Data from %s\n" \
1146  "# \n" \
1147  "# %s\n" \
1148  "# \n" \
1149  "object 1 class gridpositions counts %i %i %i\n" \
1150  "origin %12.6e %12.6e %12.6e\n" \
1151  "delta %12.6e 0.000000e+00 0.000000e+00\n" \
1152  "delta 0.000000e+00 %12.6e 0.000000e+00\n" \
1153  "delta 0.000000e+00 0.000000e+00 %12.6e\n" \
1154  "object 2 class gridconnections counts %i %i %i\n"\
1155  "object 3 class array type double rank 0 items %lu data follows\n",
1156  PACKAGE_STRING,title,nx,ny,nz,txmin,tymin,tzmin,
1157  hx,hy,hzed,nx,ny,nz,txyz);
1158  gzwrite(outfile, header, strlen(header)*sizeof(char));
1159 
1160  /* Now write the data */
1161  icol = 0;
1162  for (i=0; i<nx; i++) {
1163  for (j=0; j<ny; j++) {
1164  for (k=0; k<nz; k++) {
1165  u = k*(nx)*(ny)+j*(nx)+i;
1166  if (pvec[u] > 0.0) {
1167  sprintf(line, "%12.6e ", thee->data[u]);
1168  gzwrite(outfile, line, strlen(line)*sizeof(char));
1169  icol++;
1170  if (icol == 3) {
1171  icol = 0;
1172  gzwrite(outfile, newline, strlen(newline)*sizeof(char));
1173  }
1174  }
1175  }
1176  }
1177  }
1178  if(icol < 3){
1179  char newline[] = "\n";
1180  gzwrite(outfile, newline, strlen(newline)*sizeof(char));
1181  }
1182 
1183  /* Create the field */
1184  sprintf(footer, "attribute \"dep\" string \"positions\"\n" \
1185  "object \"regular positions regular connections\" class field\n" \
1186  "component \"positions\" value 1\n" \
1187  "component \"connections\" value 2\n" \
1188  "component \"data\" value 3\n");
1189  gzwrite(outfile, footer, strlen(footer)*sizeof(char));
1190 
1191  gzclose(outfile);
1192 #else
1193 
1194  Vnm_print(0, "WARNING\n");
1195  Vnm_print(0, "Vgrid_readGZ: gzip read/write support is disabled in this build\n");
1196  Vnm_print(0, "Vgrid_readGZ: configure and compile without the --disable-zlib flag.\n");
1197  Vnm_print(0, "WARNING\n");
1198 #endif
1199 }
1200 
1201 /* ///////////////////////////////////////////////////////////////////////////
1202 // Routine: Vgrid_writeDX
1203 //
1204 // Author: Nathan Baker
1206 VPUBLIC void Vgrid_writeDX(Vgrid *thee, const char *iodev, const char *iofmt,
1207  const char *thost, const char *fname, char *title, double *pvec) {
1208 
1209  double xmin, ymin, zmin, hx, hy, hzed;
1210  int nx, ny, nz, nxPART, nyPART, nzPART;
1211  int usepart, gotit;
1212  size_t icol, i, j, k, u;
1213  double x, y, z, xminPART, yminPART, zminPART;
1214  Vio *sock;
1215  char precFormat[VMAX_BUFSIZE];
1216 
1217  if (thee == VNULL) {
1218  Vnm_print(2, "Vgrid_writeDX: Error -- got VNULL thee!\n");
1219  VASSERT(0);
1220  }
1221  if (!(thee->ctordata || thee->readdata)) {
1222  Vnm_print(2, "Vgrid_writeDX: Error -- no data available!\n");
1223  VASSERT(0);
1224  }
1225 
1226  hx = thee->hx;
1227  hy = thee->hy;
1228  hzed = thee->hzed;
1229  nx = thee->nx;
1230  ny = thee->ny;
1231  nz = thee->nz;
1232  xmin = thee->xmin;
1233  ymin = thee->ymin;
1234  zmin = thee->zmin;
1235 
1236  if (pvec == VNULL) usepart = 0;
1237  else usepart = 1;
1238 
1239  /* Set up the virtual socket */
1240  Vnm_print(0, "Vgrid_writeDX: Opening virtual socket...\n");
1241  sock = Vio_ctor(iodev,iofmt,thost,fname,"w");
1242  if (sock == VNULL) {
1243  Vnm_print(2, "Vgrid_writeDX: Problem opening virtual socket %s\n",
1244  fname);
1245  return;
1246  }
1247  if (Vio_connect(sock, 0) < 0) {
1248  Vnm_print(2, "Vgrid_writeDX: Problem connecting virtual socket %s\n",
1249  fname);
1250  return;
1251  }
1252 
1253  Vio_setWhiteChars(sock, MCwhiteChars);
1254  Vio_setCommChars(sock, MCcommChars);
1255 
1256  Vnm_print(0, "Vgrid_writeDX: Writing to virtual socket...\n");
1257 
1258  if (usepart) {
1259  /* Get the lower corner and number of grid points for the local
1260  * partition */
1261  xminPART = VLARGE;
1262  yminPART = VLARGE;
1263  zminPART = VLARGE;
1264  nxPART = 0;
1265  nyPART = 0;
1266  nzPART = 0;
1267  /* First, search for the lower corner */
1268  for (k=0; k<nz; k++) {
1269  z = k*hzed + zmin;
1270  for (j=0; j<ny; j++) {
1271  y = j*hy + ymin;
1272  for (i=0; i<nx; i++) {
1273  x = i*hx + xmin;
1274  if (pvec[IJK(i,j,k)] > 0.0) {
1275  if (x < xminPART) xminPART = x;
1276  if (y < yminPART) yminPART = y;
1277  if (z < zminPART) zminPART = z;
1278  }
1279  }
1280  }
1281  }
1282  /* Now search for the number of grid points in the z direction */
1283  for (k=0; k<nz; k++) {
1284  gotit = 0;
1285  for (j=0; j<ny; j++) {
1286  for (i=0; i<nx; i++) {
1287  if (pvec[IJK(i,j,k)] > 0.0) {
1288  gotit = 1;
1289  break;
1290  }
1291  }
1292  if (gotit) break;
1293  }
1294  if (gotit) nzPART++;
1295  }
1296  /* Now search for the number of grid points in the y direction */
1297  for (j=0; j<ny; j++) {
1298  gotit = 0;
1299  for (k=0; k<nz; k++) {
1300  for (i=0; i<nx; i++) {
1301  if (pvec[IJK(i,j,k)] > 0.0) {
1302  gotit = 1;
1303  break;
1304  }
1305  }
1306  if (gotit) break;
1307  }
1308  if (gotit) nyPART++;
1309  }
1310  /* Now search for the number of grid points in the x direction */
1311  for (i=0; i<nx; i++) {
1312  gotit = 0;
1313  for (k=0; k<nz; k++) {
1314  for (j=0; j<ny; j++) {
1315  if (pvec[IJK(i,j,k)] > 0.0) {
1316  gotit = 1;
1317  break;
1318  }
1319  }
1320  if (gotit) break;
1321  }
1322  if (gotit) nxPART++;
1323  }
1324 
1325  if ((nxPART != nx) || (nyPART != ny) || (nzPART != nz)) {
1326  Vnm_print(0, "Vgrid_writeDX: printing only subset of domain\n");
1327  }
1328 
1329 
1330  /* Write off the title (if we're not XDR) */
1331  if (Vstring_strcasecmp(iofmt, "XDR") == 0) {
1332  Vnm_print(0, "Vgrid_writeDX: Skipping comments for XDR format.\n");
1333  } else {
1334  Vnm_print(0, "Vgrid_writeDX: Writing comments for %s format.\n",
1335  iofmt);
1336  Vio_printf(sock, "# Data from %s\n", PACKAGE_STRING);
1337  Vio_printf(sock, "# \n");
1338  Vio_printf(sock, "# %s\n", title);
1339  Vio_printf(sock, "# \n");
1340  }
1341 
1342  /* Write off the DX regular positions */
1343  Vio_printf(sock, "object 1 class gridpositions counts %d %d %d\n",
1344  nxPART, nyPART, nzPART);
1345 
1346  sprintf(precFormat, Vprecision, xminPART, yminPART, zminPART);
1347  Vio_printf(sock, "origin %s\n", precFormat);
1348  sprintf(precFormat, Vprecision, hx, 0.0, 0.0);
1349  Vio_printf(sock, "delta %s\n", precFormat);
1350  sprintf(precFormat, Vprecision, 0.0, hy, 0.0);
1351  Vio_printf(sock, "delta %s\n", precFormat);
1352  sprintf(precFormat, Vprecision, 0.0, 0.0, hzed);
1353  Vio_printf(sock, "delta %s\n", precFormat);
1354 
1355  /* Write off the DX regular connections */
1356  Vio_printf(sock, "object 2 class gridconnections counts %d %d %d\n",
1357  nxPART, nyPART, nzPART);
1358 
1359  /* Write off the DX data */
1360  Vio_printf(sock, "object 3 class array type double rank 0 items %lu \
1361 data follows\n", (nxPART*nyPART*nzPART));
1362  icol = 0;
1363  for (i=0; i<nx; i++) {
1364  for (j=0; j<ny; j++) {
1365  for (k=0; k<nz; k++) {
1366  u = k*(nx)*(ny)+j*(nx)+i;
1367  if (pvec[u] > 0.0) {
1368  Vio_printf(sock, "%12.6e ", thee->data[u]);
1369  icol++;
1370  if (icol == 3) {
1371  icol = 0;
1372  Vio_printf(sock, "\n");
1373  }
1374  }
1375  }
1376  }
1377  }
1378 
1379  if (icol != 0) Vio_printf(sock, "\n");
1380 
1381  /* Create the field */
1382  Vio_printf(sock, "attribute \"dep\" string \"positions\"\n");
1383  Vio_printf(sock, "object \"regular positions regular connections\" \
1384 class field\n");
1385  Vio_printf(sock, "component \"positions\" value 1\n");
1386  Vio_printf(sock, "component \"connections\" value 2\n");
1387  Vio_printf(sock, "component \"data\" value 3\n");
1388 
1389  } else {
1390  /* Write off the title (if we're not XDR) */
1391  if (Vstring_strcasecmp(iofmt, "XDR") == 0) {
1392  Vnm_print(0, "Vgrid_writeDX: Skipping comments for XDR format.\n");
1393  } else {
1394  Vnm_print(0, "Vgrid_writeDX: Writing comments for %s format.\n",
1395  iofmt);
1396  Vio_printf(sock, "# Data from %s\n", PACKAGE_STRING);
1397  Vio_printf(sock, "# \n");
1398  Vio_printf(sock, "# %s\n", title);
1399  Vio_printf(sock, "# \n");
1400  }
1401 
1402 
1403  /* Write off the DX regular positions */
1404  Vio_printf(sock, "object 1 class gridpositions counts %d %d %d\n",
1405  nx, ny, nz);
1406 
1407  sprintf(precFormat, Vprecision, xmin, ymin, zmin);
1408  Vio_printf(sock, "origin %s\n", precFormat);
1409  sprintf(precFormat, Vprecision, hx, 0.0, 0.0);
1410  Vio_printf(sock, "delta %s\n", precFormat);
1411  sprintf(precFormat, Vprecision, 0.0, hy, 0.0);
1412  Vio_printf(sock, "delta %s\n", precFormat);
1413  sprintf(precFormat, Vprecision, 0.0, 0.0, hzed);
1414  Vio_printf(sock, "delta %s\n", precFormat);
1415 
1416  /* Write off the DX regular connections */
1417  Vio_printf(sock, "object 2 class gridconnections counts %d %d %d\n",
1418  nx, ny, nz);
1419 
1420  /* Write off the DX data */
1421  Vio_printf(sock, "object 3 class array type double rank 0 items %lu \
1422 data follows\n", (nx*ny*nz));
1423  icol = 0;
1424  for (i=0; i<nx; i++) {
1425  for (j=0; j<ny; j++) {
1426  for (k=0; k<nz; k++) {
1427  u = k*(nx)*(ny)+j*(nx)+i;
1428  Vio_printf(sock, "%12.6e ", thee->data[u]);
1429  icol++;
1430  if (icol == 3) {
1431  icol = 0;
1432  Vio_printf(sock, "\n");
1433  }
1434  }
1435  }
1436  }
1437  if (icol != 0) Vio_printf(sock, "\n");
1438 
1439  /* Create the field */
1440  Vio_printf(sock, "attribute \"dep\" string \"positions\"\n");
1441  Vio_printf(sock, "object \"regular positions regular connections\" \
1442 class field\n");
1443  Vio_printf(sock, "component \"positions\" value 1\n");
1444  Vio_printf(sock, "component \"connections\" value 2\n");
1445  Vio_printf(sock, "component \"data\" value 3\n");
1446  }
1447 
1448  /* Close off the socket */
1449  Vio_connectFree(sock);
1450  Vio_dtor(&sock);
1451 }
1452 
1453 /* ///////////////////////////////////////////////////////////////////////////
1454 // Routine: Vgrid_writeDXBIN
1455 //
1456 // Author: Juan Brandi
1458 VPUBLIC void Vgrid_writeDXBIN(Vgrid *thee, const char *iodev, const char *iofmt,
1459  const char *thost, const char *fname, char *title, double *pvec){
1460 
1461  double xmin, ymin, zmin, hx, hy, hzed;
1462  int nx, ny, nz, nxPART, nyPART, nzPART;
1463  int usepart, gotit;
1464  size_t icol, i, j, k, u;
1465  double x, y, z, xminPART, yminPART, zminPART;
1466  //Vio *sock;
1467  char precFormat[VMAX_BUFSIZE];
1468 
1469  if (thee == VNULL) {
1470  Vnm_print(2, "Vgrid_writeDXBIN: Error -- got VNULL thee!\n");
1471  VASSERT(0);
1472  }
1473  if (!(thee->ctordata || thee->readdata)) {
1474  Vnm_print(2, "Vgrid_writeDXBIN: Error -- no data available!\n");
1475  VASSERT(0);
1476  }
1477 
1478 
1479  hx = thee->hx;
1480  hy = thee->hy;
1481  hzed = thee->hzed;
1482  nx = thee->nx;
1483  ny = thee->ny;
1484  nz = thee->nz;
1485  xmin = thee->xmin;
1486  ymin = thee->ymin;
1487  zmin = thee->zmin;
1488 
1489  if (pvec == VNULL) usepart = 0;
1490  else usepart = 1;
1491 
1492  /*will not use vio methods to try to avoid using malloc.*/
1493  FILE *fd = fopen(fname,"wb");
1494 
1495  //check to se if the file was created/open successfully.
1496  if(fd == NULL){
1497  printf("Vgrid_writeDXBIN: Problem opening file %s for writing.\n", fname);
1498  return;
1499  }
1500 
1501  printf("Vgrid_writeDXBIN: Writing to file...\n");
1502 
1503  if (usepart) {
1504  /* Get the lower corner and number of grid points for the local
1505  * partition */
1506  xminPART = VLARGE;
1507  yminPART = VLARGE;
1508  zminPART = VLARGE;
1509  nxPART = 0;
1510  nyPART = 0;
1511  nzPART = 0;
1512  /* First, search for the lower corner */
1513  for (k=0; k<nz; k++) {
1514  z = k*hzed + zmin;
1515  for (j=0; j<ny; j++) {
1516  y = j*hy + ymin;
1517  for (i=0; i<nx; i++) {
1518  x = i*hx + xmin;
1519  if (pvec[IJK(i,j,k)] > 0.0) {
1520  if (x < xminPART) xminPART = x;
1521  if (y < yminPART) yminPART = y;
1522  if (z < zminPART) zminPART = z;
1523  }
1524  }
1525  }
1526  }
1527  /* Now search for the number of grid points in the z direction */
1528  for (k=0; k<nz; k++) {
1529  gotit = 0;
1530  for (j=0; j<ny; j++) {
1531  for (i=0; i<nx; i++) {
1532  if (pvec[IJK(i,j,k)] > 0.0) {
1533  gotit = 1;
1534  break;
1535  }
1536  }
1537  if (gotit) break;
1538  }
1539  if (gotit) nzPART++;
1540  }
1541  /* Now search for the number of grid points in the y direction */
1542  for (j=0; j<ny; j++) {
1543  gotit = 0;
1544  for (k=0; k<nz; k++) {
1545  for (i=0; i<nx; i++) {
1546  if (pvec[IJK(i,j,k)] > 0.0) {
1547  gotit = 1;
1548  break;
1549  }
1550  }
1551  if (gotit) break;
1552  }
1553  if (gotit) nyPART++;
1554  }
1555  /* Now search for the number of grid points in the x direction */
1556  for (i=0; i<nx; i++) {
1557  gotit = 0;
1558  for (k=0; k<nz; k++) {
1559  for (j=0; j<ny; j++) {
1560  if (pvec[IJK(i,j,k)] > 0.0) {
1561  gotit = 1;
1562  break;
1563  }
1564  }
1565  if (gotit) break;
1566  }
1567  if (gotit) nxPART++;
1568  }
1569 
1570  if ((nxPART != nx) || (nyPART != ny) || (nzPART != nz)) {
1571  Vnm_print(0, "Vgrid_writeDXBIN: printing only subset of domain\n");
1572  }
1573 
1574  /* Write title (we're in XDR and "wb") */
1575  //Vnm_print(0, "Vgrid_writeDXBIN: Writing comments for dxbin format.\n");
1576  printf("Vgrid_writeDXBIN: Writing comments for dxbin format\n");
1577  fprintf(fd, "# Data from %s\n",PACKAGE_STRING);
1578  fprintf(fd, "# \n");
1579  fprintf(fd, "# %s\n",title);
1580  fprintf(fd, "# \n");
1581 
1582  /* Write off the DX regular positions */
1583  fprintf(fd, "object 1 class gridpositions counts %d %d %d\n",nxPART, nyPART, nzPART);
1584 
1585  sprintf(precFormat, Vprecision, xminPART, yminPART, zminPART);
1586  fprintf(fd, "origin %s\n",precFormat);
1587 
1588  sprintf(precFormat, Vprecision, hx, 0.0, 0.0);
1589  fprintf(fd, "delta %s\n",precFormat);
1590 
1591  sprintf(precFormat, Vprecision, 0.0, hy, 0.0);
1592  fprintf(fd, "delta %s\n", precFormat);
1593 
1594  sprintf(precFormat, Vprecision, 0.0, 0.0, hzed);
1595  fprintf(fd, "delta %s\n", precFormat);
1596 
1597  /* Write off the DX regular connections */
1598  fprintf(fd, "object 2 class gridconnections counts %d %d %d\n", nxPART, nyPART, nzPART);
1599 
1600  /* Write off the DX data */
1601  fprintf(fd, "object 3 class array type double rank 0 items %d binary data follows\n",(nxPART*nyPART*nzPART));
1602 
1603  icol = 0;
1604  for (i=0; i<nx; i++) {
1605  for (j=0; j<ny; j++) {
1606  for (k=0; k<nz; k++) {
1607  u = k*(nx)*(ny)+j*(nx)+i;
1608  if (pvec[u] > 0.0) {
1609  fwrite(&(thee->data)[u],sizeof(double),1,fd);
1610  icol++;
1611  /*don't need the column formating to write binary doubles.*/
1612  if (icol == 3) {
1613  icol = 0;
1614  }
1615  }
1616  }
1617  }
1618  }
1619 
1620  fprintf(fd,"\n");
1621 
1622  /* Create the field */
1623  fprintf(fd, "attribute \"dep\" string \"positions\"\n");
1624  fprintf(fd, "object \"regular positions regular connections\" class field\n");
1625  fprintf(fd, "component \"positions\" value 1\n");
1626  fprintf(fd, "component \"connections\" value 2\n");
1627  fprintf(fd, "component \"data\" value 3\n");
1628 
1629  fclose(fd);
1630 
1631  } else {
1632  /*write dx format title*/
1633  printf("Vgrid_writeDXBIN: Writing comments for %s format.\n",iofmt);
1634  fprintf(fd,"# Data from %s\n", PACKAGE_STRING);
1635  fprintf(fd,"# \n");
1636  fprintf(fd, "# %s\n", title);
1637  fprintf(fd, "# \n");
1638 
1639  /* Write off the DX regular positions */
1640  fprintf(fd, "object 1 class gridpositions counts %d %d %d\n", nx, ny, nz);
1641 
1642  sprintf(precFormat, Vprecision, xmin, ymin, zmin);
1643  fprintf(fd, "origin %s\n", precFormat);
1644 
1645  sprintf(precFormat, Vprecision, hx, 0.0, 0.0);
1646  fprintf(fd, "delta %s\n", precFormat);
1647 
1648  sprintf(precFormat, Vprecision, 0.0, hy, 0.0);
1649  fprintf(fd, "delta %s\n", precFormat);
1650 
1651  sprintf(precFormat, Vprecision, 0.0, 0.0, hzed);
1652  fprintf(fd, "delta %s\n", precFormat);
1653 
1654  /* Write off the DX regular connections */
1655  fprintf(fd, "object 2 class gridconnections counts %d %d %d\n", nx, ny, nz);
1656 
1657  /* Write off the DX data */
1658  fprintf(fd, "object 3 class array type double rank 0 items %d binary data follows\n", (nx*ny*nz));
1659 
1660  icol = 0;
1661  for (i=0; i<nx; i++) {
1662  for (j=0; j<ny; j++) {
1663  for (k=0; k<nz; k++) {
1664  u = k*(nx)*(ny)+j*(nx)+i;
1665  fwrite(&(thee->data)[u],sizeof(double),1,fd);
1666  icol++;
1667  if (icol == 3) {
1668  icol = 0;
1669  }
1670  }
1671  }
1672  }
1673 
1674  fprintf(fd, "\n");
1675 
1676  /* Create the field */
1677  fprintf(fd, "attribute \"dep\" string \"positions\"\n");
1678  fprintf(fd, "object \"regular positions regular connections\" class field\n");
1679  fprintf(fd, "component \"positions\" value 1\n");
1680  fprintf(fd, "component \"connections\" value 2\n");
1681  fprintf(fd, "component \"data\" value 3\n");
1682 
1683  fclose(fd);
1684  }
1685 }
1686 
1687 
1688 /* ///////////////////////////////////////////////////////////////////////////
1689 // Routine: Vgrid_writeUHBD
1690 // Author: Nathan Baker
1692 VPUBLIC void Vgrid_writeUHBD(Vgrid *thee, const char *iodev, const char *iofmt,
1693  const char *thost, const char *fname, char *title, double *pvec) {
1694 
1695  size_t u, icol, i, j, k;
1696  size_t gotit, nx, ny, nz;
1697  double xmin, ymin, zmin, hzed, hy, hx;
1698  Vio *sock;
1699 
1700  if (thee == VNULL) {
1701  Vnm_print(2, "Vgrid_writeUHBD: Error -- got VNULL thee!\n");
1702  VASSERT(0);
1703  }
1704  if (!(thee->ctordata || thee->readdata)) {
1705  Vnm_print(2, "Vgrid_writeUHBD: Error -- no data available!\n");
1706  VASSERT(0);
1707  }
1708 
1709  if ((thee->hx!=thee->hy) || (thee->hy!=thee->hzed)
1710  || (thee->hx!=thee->hzed)) {
1711  Vnm_print(2, "Vgrid_writeUHBD: can't write UHBD mesh with non-uniform \
1712 spacing\n");
1713  return;
1714  }
1715 
1716  /* Set up the virtual socket */
1717  sock = Vio_ctor(iodev,iofmt,thost,fname,"w");
1718  if (sock == VNULL) {
1719  Vnm_print(2, "Vgrid_writeUHBD: Problem opening virtual socket %s\n",
1720  fname);
1721  return;
1722  }
1723  if (Vio_connect(sock, 0) < 0) {
1724  Vnm_print(2, "Vgrid_writeUHBD: Problem connecting virtual socket %s\n",
1725  fname);
1726  return;
1727  }
1728 
1729  /* Get the lower corner and number of grid points for the local
1730  * partition */
1731  hx = thee->hx;
1732  hy = thee->hy;
1733  hzed = thee->hzed;
1734  nx = thee->nx;
1735  ny = thee->ny;
1736  nz = thee->nz;
1737  xmin = thee->xmin;
1738  ymin = thee->ymin;
1739  zmin = thee->zmin;
1740 
1741  /* Let interested folks know that partition information is ignored */
1742  if (pvec != VNULL) {
1743  gotit = 0;
1744  for (i=0; i<(nx*ny*nz); i++) {
1745  if (pvec[i] == 0) {
1746  gotit = 1;
1747  break;
1748  }
1749  }
1750  if (gotit) {
1751  Vnm_print(2, "Vgrid_writeUHBD: IGNORING PARTITION INFORMATION!\n");
1752  Vnm_print(2, "Vgrid_writeUHBD: This means I/O from parallel runs \
1753 will have significant overlap.\n");
1754  }
1755  }
1756 
1757  /* Write out the header */
1758  Vio_printf(sock, "%72s\n", title);
1759  Vio_printf(sock, "%12.5e%12.5e%7d%7d%7d%7d%7d\n", 1.0, 0.0, -1, 0,
1760  nz, 1, nz);
1761  Vio_printf(sock, "%7d%7d%7d%12.5e%12.5e%12.5e%12.5e\n", nx, ny, nz,
1762  hx, (xmin-hx), (ymin-hx), (zmin-hx));
1763  Vio_printf(sock, "%12.5e%12.5e%12.5e%12.5e\n", 0.0, 0.0, 0.0, 0.0);
1764  Vio_printf(sock, "%12.5e%12.5e%7d%7d", 0.0, 0.0, 0, 0);
1765 
1766  /* Write out the entries */
1767  icol = 0;
1768  for (k=0; k<nz; k++) {
1769  Vio_printf(sock, "\n%7d%7d%7d\n", k+1, thee->nx, thee->ny);
1770  icol = 0;
1771  for (j=0; j<ny; j++) {
1772  for (i=0; i<nx; i++) {
1773  u = k*(nx)*(ny)+j*(nx)+i;
1774  icol++;
1775  Vio_printf(sock, " %12.5e", thee->data[u]);
1776  if (icol == 6) {
1777  icol = 0;
1778  Vio_printf(sock, "\n");
1779  }
1780  }
1781  }
1782  }
1783  if (icol != 0) Vio_printf(sock, "\n");
1784 
1785  /* Close off the socket */
1786  Vio_connectFree(sock);
1787  Vio_dtor(&sock);
1788 }
1789 
1790 VPUBLIC double Vgrid_integrate(Vgrid *thee) {
1791 
1792  size_t i, j, k;
1793  int nx, ny, nz;
1794  double sum, w;
1795 
1796  if (thee == VNULL) {
1797  Vnm_print(2, "Vgrid_integrate: Got VNULL thee!\n");
1798  VASSERT(0);
1799  }
1800 
1801  nx = thee->nx;
1802  ny = thee->ny;
1803  nz = thee->nz;
1804 
1805  sum = 0.0;
1806 
1807  for (k=0; k<nz; k++) {
1808  w = 1.0;
1809  if ((k==0) || (k==(nz-1))) w = w * 0.5;
1810  for (j=0; j<ny; j++) {
1811  w = 1.0;
1812  if ((j==0) || (j==(ny-1))) w = w * 0.5;
1813  for (i=0; i<nx; i++) {
1814  w = 1.0;
1815  if ((i==0) || (i==(nx-1))) w = w * 0.5;
1816  sum = sum + w*(thee->data[IJK(i,j,k)]);
1817  }
1818  }
1819  }
1820 
1821  sum = sum*(thee->hx)*(thee->hy)*(thee->hzed);
1822 
1823  return sum;
1824 
1825 }
1826 
1827 
1828 VPUBLIC double Vgrid_normL1(Vgrid *thee) {
1829 
1830  size_t i, j, k;
1831  int nx, ny, nz;
1832  double sum;
1833 
1834  if (thee == VNULL) {
1835  Vnm_print(2, "Vgrid_normL1: Got VNULL thee!\n");
1836  VASSERT(0);
1837  }
1838 
1839  nx = thee->nx;
1840  ny = thee->ny;
1841  nz = thee->nz;
1842 
1843  sum = 0.0;
1844  for (k=0; k<nz; k++) {
1845  for (j=0; j<ny; j++) {
1846  for (i=0; i<nx; i++) {
1847  sum = sum + VABS(thee->data[IJK(i,j,k)]);
1848  }
1849  }
1850  }
1851 
1852  sum = sum*(thee->hx)*(thee->hy)*(thee->hzed);
1853 
1854  return sum;
1855 
1856 }
1857 
1858 VPUBLIC double Vgrid_normL2(Vgrid *thee) {
1859 
1860  size_t i, j, k;
1861  int nx, ny, nz;
1862  double sum;
1863 
1864  if (thee == VNULL) {
1865  Vnm_print(2, "Vgrid_normL2: Got VNULL thee!\n");
1866  VASSERT(0);
1867  }
1868 
1869  nx = thee->nx;
1870  ny = thee->ny;
1871  nz = thee->nz;
1872 
1873  sum = 0.0;
1874  for (k=0; k<nz; k++) {
1875  for (j=0; j<ny; j++) {
1876  for (i=0; i<nx; i++) {
1877  sum = sum + VSQR(thee->data[IJK(i,j,k)]);
1878  }
1879  }
1880  }
1881 
1882  sum = sum*(thee->hx)*(thee->hy)*(thee->hzed);
1883 
1884  return VSQRT(sum);
1885 
1886 }
1887 
1888 VPUBLIC double Vgrid_seminormH1(Vgrid *thee) {
1889 
1890  size_t i, j, k;
1891  int nx, ny, nz, d;
1892  double pt[3], grad[3], sum, hx, hy, hzed, xmin, ymin, zmin;
1893 
1894  if (thee == VNULL) {
1895  Vnm_print(2, "Vgrid_seminormH1: Got VNULL thee!\n");
1896  VASSERT(0);
1897  }
1898 
1899  nx = thee->nx;
1900  ny = thee->ny;
1901  nz = thee->nz;
1902  hx = thee->hx;
1903  hy = thee->hy;
1904  hzed = thee->hzed;
1905  xmin = thee->xmin;
1906  ymin = thee->ymin;
1907  zmin = thee->zmin;
1908 
1909  sum = 0.0;
1910  for (k=0; k<nz; k++) {
1911  pt[2] = k*hzed + zmin;
1912  for (j=0; j<ny; j++) {
1913  pt[1] = j*hy + ymin;
1914  for (i=0; i<nx; i++) {
1915  pt[0] = i*hx + xmin;
1916  VASSERT(Vgrid_gradient(thee, pt, grad));
1917  for (d=0; d<3; d++) sum = sum + VSQR(grad[d]);
1918  }
1919  }
1920  }
1921 
1922  sum = sum*(hx)*(hy)*(hzed);
1923 
1924  if (VABS(sum) < VSMALL) sum = 0.0;
1925  else sum = VSQRT(sum);
1926 
1927  return sum;
1928 
1929 }
1930 
1931 VPUBLIC double Vgrid_normH1(Vgrid *thee) {
1932 
1933  double sum = 0.0;
1934 
1935  if (thee == VNULL) {
1936  Vnm_print(2, "Vgrid_normH1: Got VNULL thee!\n");
1937  VASSERT(0);
1938  }
1939 
1940  sum = VSQR(Vgrid_seminormH1(thee)) + VSQR(Vgrid_normL2(thee));
1941 
1942  return VSQRT(sum);
1943 
1944 }
1945 
1946 VPUBLIC double Vgrid_normLinf(Vgrid *thee) {
1947 
1948  size_t i, j, k;
1949  int nx, ny, nz, gotval;
1950  double sum, val;
1951 
1952  if (thee == VNULL) {
1953  Vnm_print(2, "Vgrid_normLinf: Got VNULL thee!\n");
1954  VASSERT(0);
1955  }
1956 
1957  nx = thee->nx;
1958  ny = thee->ny;
1959  nz = thee->nz;
1960 
1961  sum = 0.0;
1962  gotval = 0;
1963  for (k=0; k<nz; k++) {
1964  for (j=0; j<ny; j++) {
1965  for (i=0; i<nx; i++) {
1966  val = VABS(thee->data[IJK(i,j,k)]);
1967  if (!gotval) {
1968  gotval = 1;
1969  sum = val;
1970  }
1971  if (val > sum) sum = val;
1972  }
1973  }
1974  }
1975 
1976  return sum;
1977 
1978 }
1979 
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:1931
VPUBLIC double Vgrid_normL1(Vgrid *thee)
Get the norm of the data. This returns the integral: .
Definition: vgrid.c:1828
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:1946
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
VPUBLIC int Vgrid_readDXBIN(Vgrid *thee, const char *iodev, const char *iofmt, const char *thost, const char *fname)
Read in binary data in OpenDX grid format.
Definition: vgrid.c:810
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:1790
#define VEMBED(rctag)
Allows embedding of RCS ID tags in object files.
Definition: vhal.h:556
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:379
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:63
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:586
VPUBLIC double Vgrid_normL2(Vgrid *thee)
Get the norm of the data. This returns the integral: .
Definition: vgrid.c:1858
VPUBLIC double Vgrid_seminormH1(Vgrid *thee)
Get the semi-norm of the data. This returns the integral: .
Definition: vgrid.c:1888
VPUBLIC int Vgrid_value(Vgrid *thee, double pt[3], double *value)
Get potential value (from mesh or approximation) at a point.
Definition: vgrid.c:179
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