APBS 3.0.0
Loading...
Searching...
No Matches
vgrid.c
Go to the documentation of this file.
1
57#include "vgrid.h"
58#include <stdio.h>
59
60VEMBED(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>
71int isnan(double d)
72{
73 return _isnan(d);
74}
75#endif
76
77VPRIVATE char *MCwhiteChars = " =,;\t\n";
78VPRIVATE char *MCcommChars = "#%";
79VPRIVATE double Vcompare;
80VPRIVATE char Vprecision[26];
81
82/* ///////////////////////////////////////////////////////////////////////////
83// Routine: Vgrid_ctor
84// Author: Nathan Baker
86VPUBLIC 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
112VPUBLIC 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
152VPUBLIC 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
165VPUBLIC 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
179VPUBLIC 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
299VPUBLIC 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
379VPUBLIC 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
462VPUBLIC 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
586VPUBLIC 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
810VPUBLIC 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
1011VPUBLIC 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
1206VPUBLIC 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 \
1361data 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\" \
1384class 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 \
1422data 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\" \
1442class 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
1458VPUBLIC 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
1692VPUBLIC 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 \
1712spacing\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 \
1753will 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
1790VPUBLIC 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
1828VPUBLIC 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
1858VPUBLIC 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
1888VPUBLIC 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
1931VPUBLIC 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
1946VPUBLIC 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}
VPUBLIC unsigned long int Vgrid_memChk(Vgrid *thee)
Return the memory used by this structure (and its contents) in bytes.
Definition vgrid.c:63
VPUBLIC double Vgrid_normLinf(Vgrid *thee)
Get the norm of the data. This returns the integral:
Definition vgrid.c:1946
VPUBLIC double Vgrid_integrate(Vgrid *thee)
Get the integral of the data.
Definition vgrid.c:1790
VPUBLIC int Vgrid_gradient(Vgrid *thee, double pt[3], double grad[3])
Get first derivative values at a point.
Definition vgrid.c:379
VPUBLIC double Vgrid_normL1(Vgrid *thee)
Get the norm of the data. This returns the integral:
Definition vgrid.c:1828
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_normL2(Vgrid *thee)
Get the norm of the data. This returns the integral:
Definition vgrid.c:1858
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 int Vgrid_value(Vgrid *thee, double pt[3], double *value)
Get potential value (from mesh or approximation) at a point.
Definition vgrid.c:179
VPUBLIC double Vgrid_seminormH1(Vgrid *thee)
Get the semi-norm of the data. This returns the integral:
Definition vgrid.c:1888
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
#define VEMBED(rctag)
Allows embedding of RCS ID tags in object files.
Definition vhal.h:556
@ VRC_FAILURE
Definition vhal.h:69
@ VRC_SUCCESS
Definition vhal.h:70
VPRIVATE char * MCcommChars
Comment characters for socket reads.
Definition vparam.c:71
VPRIVATE char * MCwhiteChars
Whitespace characters for socket reads.
Definition vparam.c:65
VPUBLIC int Vstring_strcasecmp(const char *s1, const char *s2)
Case-insensitive string comparison (BSD standard)
Definition vstring.c:66
Electrostatic potential oracle for Cartesian mesh data.
Definition vgrid.h:81
int nx
Definition vgrid.h:83
double hx
Definition vgrid.h:86
double hy
Definition vgrid.h:87
Vmem * mem
Definition vgrid.h:99
double * data
Definition vgrid.h:95
double xmax
Definition vgrid.h:92
double ymin
Definition vgrid.h:90
int ny
Definition vgrid.h:84
double zmin
Definition vgrid.h:91
double xmin
Definition vgrid.h:89
int readdata
Definition vgrid.h:96
int nz
Definition vgrid.h:85
int ctordata
Definition vgrid.h:97
double zmax
Definition vgrid.h:94
double hzed
Definition vgrid.h:88
double ymax
Definition vgrid.h:93
Potential oracle for Cartesian mesh data.