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,
117 if (thee == VNULL) return 0;
125 thee->xmax = xmin + (nx-1)*hx;
127 thee->ymax = ymin + (ny-1)*hy;
129 thee->zmax = zmin + (nz-1)*hzed;
139 thee->mem = Vmem_ctor("APBS:VGRID");
141 Vcompare = pow(10,-1*(VGRID_DIGITS - 2));
142 sprintf(Vprecision,"%%12.%de %%12.%de %%12.%de", VGRID_DIGITS,
143 VGRID_DIGITS, VGRID_DIGITS);
179VPUBLIC int Vgrid_value(Vgrid *thee, double pt[3], double *value) {
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;
188 Vnm_print(2, "Vgrid_value: Error -- got VNULL thee!\n");
191 if (!(thee->ctordata || thee->readdata)) {
192 Vnm_print(2, "Vgrid_value: Error -- no data available!\n");
211 ifloat = (pt[0] - xmin)/hx;
212 jfloat = (pt[1] - ymin)/hy;
213 kfloat = (pt[2] - zmin)/hzed;
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;
228 /* See if we're on the mesh */
230 if ((ihi<nx) && (jhi<ny) && (khi<nz)
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)]);
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",
253 Vnm_print(2,
"Vgrid_value: (ilo, jlo, klo) = (%d, %d, %d)\n",
255 Vnm_print(2,
"Vgrid_value: (nx, ny, nz) = (%d, %d, %d)\n",
257 Vnm_print(2,
"Vgrid_value: (dx, dy, dz) = (%4.3f, %4.3f, %4.3f)\n",
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)]);
299VPUBLIC int Vgrid_curvature(Vgrid *thee, double pt[3], int cflag,
302 double hx, hy, hzed, curv;
303 double dxx, dyy, dzz;
304 double uleft, umid, uright, testpt[3];
307 Vnm_print(2, "Vgrid_curvature: Error -- got VNULL thee!\n");
310 if (!(thee->ctordata || thee->readdata)) {
311 Vnm_print(2, "Vgrid_curvature: Error -- no data available!\n");
325 /* Compute 2nd derivative in the x-direction */
327 testpt[0] = pt[0] - hx;
329 testpt[0] = pt[0] + hx;
333 dxx = (uright - 2*umid + uleft)/(hx*hx);
337 testpt[1] = pt[1] - hy;
339 testpt[1] = pt[1] + hy;
343 dyy = (uright - 2*umid + uleft)/(hy*hy);
347 testpt[2] = pt[2] - hzed;
349 testpt[2] = pt[2] + hzed;
352 dzz = (uright - 2*umid + uleft)/(hzed*hzed);
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;
362 Vnm_print(2,
"Vgrid_curvature: support for cflag = %d not available!\n", cflag);
379VPUBLIC int Vgrid_gradient(Vgrid *thee, double pt[3], double grad[3]) {
382 double uleft, umid, uright, testpt[3];
383 int haveleft, haveright;
386 Vnm_print(2, "Vgrid_gradient: Error -- got VNULL thee!\n");
389 if (!(thee->ctordata || thee->readdata)) {
390 Vnm_print(2, "Vgrid_gradient: Error -- no data available!\n");
398 /* Compute derivative in the x-direction */
403 testpt[0] = pt[0] - hx;
404 if (
Vgrid_value( thee, testpt, &uleft)) haveleft = 1;
406 testpt[0] = pt[0] + hx;
407 if (
Vgrid_value( thee, testpt, &uright)) haveright = 1;
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;
419 testpt[1] = pt[1] - hy;
420 if (
Vgrid_value( thee, testpt, &uleft)) haveleft = 1;
422 testpt[1] = pt[1] + hy;
423 if (
Vgrid_value( thee, testpt, &uright)) haveright = 1;
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;
435 testpt[2] = pt[2] - hzed;
436 if (
Vgrid_value( thee, testpt, &uleft)) haveleft = 1;
438 testpt[2] = pt[2] + hzed;
439 if (
Vgrid_value( thee, testpt, &uright)) haveright = 1;
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;
462VPUBLIC int Vgrid_readGZ(Vgrid *thee, const char *fname) {
466 size_t len; // Temporary counter variable for loop conditionals
469 double dtmp1, dtmp2, dtmp3;
471 char line[VMAX_ARGLEN];
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));
485 infile = gzopen(fname,
"rb");
486 if (infile == Z_NULL) {
487 Vnm_print(2,
"%s: Problem opening compressed file %s\n", __func__, fname);
497 if(gzgets(infile, line, VMAX_ARGLEN) == Z_NULL){
502 if(strncmp(line,
"#", 1) == 0)
continue;
503 if(line[0] ==
'\n')
continue;
507 sscanf(line,
"object 1 class gridpositions counts %d %d %d",
508 &(thee->nx),&(thee->ny),&(thee->nz));
511 sscanf(line,
"origin %lf %lf %lf",
512 &(thee->xmin),&(thee->ymin),&(thee->zmin));
517 sscanf(line,
"delta %lf %lf %lf",&dtmp1,&dtmp2,&dtmp3);
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;
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__);
545 temp = (
double *)malloc(len * (2 *
sizeof(
double)));
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]);
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++];
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;
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");
593 size_t i, j, k, itmp, u;
595 char tok[VMAX_BUFSIZE];
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)); }
607 sock = Vio_ctor(iodev,iofmt,thost,fname,
"r");
609 Vnm_print(2,
"Vgrid_readDX: Problem opening virtual socket %s\n",
613 if (Vio_accept(sock, 0) < 0) {
614 Vnm_print(2,
"Vgrid_readDX: Problem accepting virtual socket %s\n",
624 VJMPERR2(1 == Vio_scanf(sock,
"%s", tok));
625 VJMPERR1(!strcmp(tok,
"object"));
627 VJMPERR2(1 == Vio_scanf(sock,
"%d", &itmp));
629 VJMPERR2(1 == Vio_scanf(sock,
"%s", tok));
630 VJMPERR1(!strcmp(tok,
"class"));
632 VJMPERR2(1 == Vio_scanf(sock,
"%s", tok));
633 VJMPERR1(!strcmp(tok,
"gridpositions"));
635 VJMPERR2(1 == Vio_scanf(sock,
"%s", tok));
636 VJMPERR1(!strcmp(tok,
"counts"));
638 VJMPERR2(1 == Vio_scanf(sock,
"%s", tok));
639 VJMPERR1(1 == sscanf(tok,
"%d", &(thee->
nx)));
641 VJMPERR2(1 == Vio_scanf(sock,
"%s", tok));
642 VJMPERR1(1 == sscanf(tok,
"%d", &(thee->
ny)));
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);
649 VJMPERR2(1 == Vio_scanf(sock,
"%s", tok));
650 VJMPERR1(!strcmp(tok,
"origin"));
652 VJMPERR2(1 == Vio_scanf(sock,
"%s", tok));
653 VJMPERR1(1 == sscanf(tok,
"%lf", &(thee->
xmin)));
655 VJMPERR2(1 == Vio_scanf(sock,
"%s", tok));
656 VJMPERR1(1 == sscanf(tok,
"%lf", &(thee->
ymin)));
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",
663 VJMPERR2(1 == Vio_scanf(sock,
"%s", tok));
664 VJMPERR1(!strcmp(tok,
"delta"));
666 VJMPERR2(1 == Vio_scanf(sock,
"%s", tok));
667 VJMPERR1(1 == sscanf(tok,
"%lf", &(thee->
hx)));
669 VJMPERR2(1 == Vio_scanf(sock,
"%s", tok));
670 VJMPERR1(1 == sscanf(tok,
"%lf", &dtmp));
671 VJMPERR1(dtmp == 0.0);
673 VJMPERR2(1 == Vio_scanf(sock,
"%s", tok));
674 VJMPERR1(1 == sscanf(tok,
"%lf", &dtmp));
675 VJMPERR1(dtmp == 0.0);
677 VJMPERR2(1 == Vio_scanf(sock,
"%s", tok));
678 VJMPERR1(!strcmp(tok,
"delta"));
680 VJMPERR2(1 == Vio_scanf(sock,
"%s", tok));
681 VJMPERR1(1 == sscanf(tok,
"%lf", &dtmp));
682 VJMPERR1(dtmp == 0.0);
684 VJMPERR2(1 == Vio_scanf(sock,
"%s", tok));
685 VJMPERR1(1 == sscanf(tok,
"%lf", &(thee->
hy)));
687 VJMPERR2(1 == Vio_scanf(sock,
"%s", tok));
688 VJMPERR1(1 == sscanf(tok,
"%lf", &dtmp));
689 VJMPERR1(dtmp == 0.0);
691 VJMPERR2(1 == Vio_scanf(sock,
"%s", tok));
692 VJMPERR1(!strcmp(tok,
"delta"));
694 VJMPERR2(1 == Vio_scanf(sock,
"%s", tok));
695 VJMPERR1(1 == sscanf(tok,
"%lf", &dtmp));
696 VJMPERR1(dtmp == 0.0);
698 VJMPERR2(1 == Vio_scanf(sock,
"%s", tok));
699 VJMPERR1(1 == sscanf(tok,
"%lf", &dtmp));
700 VJMPERR1(dtmp == 0.0);
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",
707 VJMPERR2(1 == Vio_scanf(sock,
"%s", tok));
708 VJMPERR1(!strcmp(tok,
"object"));
710 VJMPERR2(1 == Vio_scanf(sock,
"%s", tok));
712 VJMPERR2(1 == Vio_scanf(sock,
"%s", tok));
713 VJMPERR1(!strcmp(tok,
"class"));
715 VJMPERR2(1 == Vio_scanf(sock,
"%s", tok));
716 VJMPERR1(!strcmp(tok,
"gridconnections"));
718 VJMPERR2(1 == Vio_scanf(sock,
"%s", tok));
719 VJMPERR1(!strcmp(tok,
"counts"));
721 VJMPERR2(1 == Vio_scanf(sock,
"%s", tok));
722 VJMPERR2(1 == Vio_scanf(sock,
"%s", tok));
723 VJMPERR2(1 == Vio_scanf(sock,
"%s", tok));
725 VJMPERR2(1 == Vio_scanf(sock,
"%s", tok));
726 VJMPERR1(!strcmp(tok,
"object"));
728 VJMPERR2(1 == Vio_scanf(sock,
"%s", tok));
730 VJMPERR2(1 == Vio_scanf(sock,
"%s", tok));
731 VJMPERR1(!strcmp(tok,
"class"));
733 VJMPERR2(1 == Vio_scanf(sock,
"%s", tok));
734 VJMPERR1(!strcmp(tok,
"array"));
736 VJMPERR2(1 == Vio_scanf(sock,
"%s", tok));
737 VJMPERR1(!strcmp(tok,
"type"));
739 VJMPERR2(1 == Vio_scanf(sock,
"%s", tok));
740 VJMPERR1(!strcmp(tok,
"double"));
742 VJMPERR2(1 == Vio_scanf(sock,
"%s", tok));
743 VJMPERR1(!strcmp(tok,
"rank"));
745 VJMPERR2(1 == Vio_scanf(sock,
"%s", tok));
747 VJMPERR2(1 == Vio_scanf(sock,
"%s", tok));
748 VJMPERR1(!strcmp(tok,
"items"));
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;
755 VJMPERR2(1 == Vio_scanf(sock,
"%s", tok));
756 VJMPERR1(!strcmp(tok,
"data"));
758 VJMPERR2(1 == Vio_scanf(sock,
"%s", tok));
759 VJMPERR1(!strcmp(tok,
"follows"));
762 Vnm_print(0,
"Vgrid_readDX: allocating %d x %d x %d doubles for storage\n",
763 thee->
nx, thee->
ny, thee->
nz);
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");
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;
788 Vio_acceptFree(sock);
795 Vnm_print(2,
"Vgrid_readDX: Format problem with input file <%s>\n",
801 Vnm_print(2,
"Vgrid_readDX: I/O problem with input file <%s>\n",
1011VPUBLIC void Vgrid_writeGZ(Vgrid *thee, const char *iodev, const char *iofmt,
1012 const char *thost, const char *fname, char *title, double *pvec) {
1015 double xmin, ymin, zmin, hx, hy, hzed;
1017 int nx, ny, nz, nxPART, nyPART, nzPART;
1019 size_t icol, i, j, k, u;
1020 double x, y, z, xminPART, yminPART, zminPART;
1023 double txmin, tymin, tzmin;
1028 char newline[] = "\n";
1030 char precFormat[VMAX_BUFSIZE];
1032 if (thee == VNULL) {
1033 Vnm_print(2, "Vgrid_writeGZ: Error -- got VNULL thee!\n");
1036 if (!(thee->ctordata || thee->readdata)) {
1037 Vnm_print(2, "Vgrid_writeGZ: Error -- no data available!\n");
1051 if (pvec == VNULL) usepart = 0;
1054 /* Set up the virtual socket */
1055 Vnm_print(0,
"Vgrid_writeGZ: Opening file...\n");
1056 outfile = gzopen(fname,
"wb");
1068 for (k=0; k<nz; k++) {
1070 for (j=0; j<ny; j++) {
1072 for (i=0; i<nx; i++) {
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;
1083 for (k=0; k<nz; k++) {
1085 for (j=0; j<ny; j++) {
1086 for (i=0; i<nx; i++) {
1087 if (pvec[IJK(i,j,k)] > 0.0) {
1094 if (gotit) nzPART++;
1097 for (j=0; j<ny; j++) {
1099 for (k=0; k<nz; k++) {
1100 for (i=0; i<nx; i++) {
1101 if (pvec[IJK(i,j,k)] > 0.0) {
1108 if (gotit) nyPART++;
1111 for (i=0; i<nx; i++) {
1113 for (k=0; k<nz; k++) {
1114 for (j=0; j<ny; j++) {
1115 if (pvec[IJK(i,j,k)] > 0.0) {
1122 if (gotit) nxPART++;
1125 if ((nxPART != nx) || (nyPART != ny) || (nzPART != nz)) {
1126 Vnm_print(0,
"Vgrid_writeGZ: printing only subset of domain\n");
1129 txyz = (nxPART*nyPART*nzPART);
1145 "# Data from %s\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));
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));
1172 gzwrite(outfile, newline, strlen(newline)*
sizeof(
char));
1179 char newline[] =
"\n";
1180 gzwrite(outfile, newline, strlen(newline)*
sizeof(
char));
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));
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");
1206VPUBLIC void Vgrid_writeDX(Vgrid *thee, const char *iodev, const char *iofmt,
1207 const char *thost, const char *fname, char *title, double *pvec) {
1209 double xmin, ymin, zmin, hx, hy, hzed;
1210 int nx, ny, nz, nxPART, nyPART, nzPART;
1212 size_t icol, i, j, k, u;
1213 double x, y, z, xminPART, yminPART, zminPART;
1215 char precFormat[VMAX_BUFSIZE];
1217 if (thee == VNULL) {
1218 Vnm_print(2, "Vgrid_writeDX: Error -- got VNULL thee!\n");
1221 if (!(thee->ctordata || thee->readdata)) {
1222 Vnm_print(2, "Vgrid_writeDX: Error -- no data available!\n");
1236 if (pvec == VNULL) usepart = 0;
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",
1247 if (Vio_connect(sock, 0) < 0) {
1248 Vnm_print(2,
"Vgrid_writeDX: Problem connecting virtual socket %s\n",
1253 Vio_setWhiteChars(sock, MCwhiteChars);
1254 Vio_setCommChars(sock, MCcommChars);
1256 Vnm_print(0,
"Vgrid_writeDX: Writing to virtual socket...\n");
1268 for (k=0; k<nz; k++) {
1270 for (j=0; j<ny; j++) {
1272 for (i=0; i<nx; i++) {
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;
1283 for (k=0; k<nz; k++) {
1285 for (j=0; j<ny; j++) {
1286 for (i=0; i<nx; i++) {
1287 if (pvec[IJK(i,j,k)] > 0.0) {
1294 if (gotit) nzPART++;
1297 for (j=0; j<ny; j++) {
1299 for (k=0; k<nz; k++) {
1300 for (i=0; i<nx; i++) {
1301 if (pvec[IJK(i,j,k)] > 0.0) {
1308 if (gotit) nyPART++;
1311 for (i=0; i<nx; i++) {
1313 for (k=0; k<nz; k++) {
1314 for (j=0; j<ny; j++) {
1315 if (pvec[IJK(i,j,k)] > 0.0) {
1322 if (gotit) nxPART++;
1325 if ((nxPART != nx) || (nyPART != ny) || (nzPART != nz)) {
1326 Vnm_print(0,
"Vgrid_writeDX: printing only subset of domain\n");
1332 Vnm_print(0,
"Vgrid_writeDX: Skipping comments for XDR format.\n");
1334 Vnm_print(0,
"Vgrid_writeDX: Writing comments for %s format.\n",
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");
1343 Vio_printf(sock,
"object 1 class gridpositions counts %d %d %d\n",
1344 nxPART, nyPART, nzPART);
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);
1356 Vio_printf(sock,
"object 2 class gridconnections counts %d %d %d\n",
1357 nxPART, nyPART, nzPART);
1360 Vio_printf(sock,
"object 3 class array type double rank 0 items %lu \
1361data follows\n", (nxPART*nyPART*nzPART));
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]);
1372 Vio_printf(sock,
"\n");
1379 if (icol != 0) Vio_printf(sock,
"\n");
1382 Vio_printf(sock,
"attribute \"dep\" string \"positions\"\n");
1383 Vio_printf(sock,
"object \"regular positions regular connections\" \
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");
1392 Vnm_print(0,
"Vgrid_writeDX: Skipping comments for XDR format.\n");
1394 Vnm_print(0,
"Vgrid_writeDX: Writing comments for %s format.\n",
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");
1404 Vio_printf(sock,
"object 1 class gridpositions counts %d %d %d\n",
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);
1417 Vio_printf(sock,
"object 2 class gridconnections counts %d %d %d\n",
1421 Vio_printf(sock,
"object 3 class array type double rank 0 items %lu \
1422data follows\n", (nx*ny*nz));
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]);
1432 Vio_printf(sock,
"\n");
1437 if (icol != 0) Vio_printf(sock,
"\n");
1440 Vio_printf(sock,
"attribute \"dep\" string \"positions\"\n");
1441 Vio_printf(sock,
"object \"regular positions regular connections\" \
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");
1449 Vio_connectFree(sock);
1458VPUBLIC void Vgrid_writeDXBIN(Vgrid *thee, const char *iodev, const char *iofmt,
1459 const char *thost, const char *fname, char *title, double *pvec){
1461 double xmin, ymin, zmin, hx, hy, hzed;
1462 int nx, ny, nz, nxPART, nyPART, nzPART;
1464 size_t icol, i, j, k, u;
1465 double x, y, z, xminPART, yminPART, zminPART;
1467 char precFormat[VMAX_BUFSIZE];
1469 if (thee == VNULL) {
1470 Vnm_print(2, "Vgrid_writeDXBIN: Error -- got VNULL thee!\n");
1473 if (!(thee->ctordata || thee->readdata)) {
1474 Vnm_print(2, "Vgrid_writeDXBIN: Error -- no data available!\n");
1489 if (pvec == VNULL) usepart = 0;
1492 /*will not use vio methods to try to avoid using malloc.*/
1493 FILE *fd = fopen(fname,
"wb");
1497 printf(
"Vgrid_writeDXBIN: Problem opening file %s for writing.\n", fname);
1501 printf(
"Vgrid_writeDXBIN: Writing to file...\n");
1513 for (k=0; k<nz; k++) {
1515 for (j=0; j<ny; j++) {
1517 for (i=0; i<nx; i++) {
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;
1528 for (k=0; k<nz; k++) {
1530 for (j=0; j<ny; j++) {
1531 for (i=0; i<nx; i++) {
1532 if (pvec[IJK(i,j,k)] > 0.0) {
1539 if (gotit) nzPART++;
1542 for (j=0; j<ny; j++) {
1544 for (k=0; k<nz; k++) {
1545 for (i=0; i<nx; i++) {
1546 if (pvec[IJK(i,j,k)] > 0.0) {
1553 if (gotit) nyPART++;
1556 for (i=0; i<nx; i++) {
1558 for (k=0; k<nz; k++) {
1559 for (j=0; j<ny; j++) {
1560 if (pvec[IJK(i,j,k)] > 0.0) {
1567 if (gotit) nxPART++;
1570 if ((nxPART != nx) || (nyPART != ny) || (nzPART != nz)) {
1571 Vnm_print(0,
"Vgrid_writeDXBIN: printing only subset of domain\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");
1583 fprintf(fd,
"object 1 class gridpositions counts %d %d %d\n",nxPART, nyPART, nzPART);
1585 sprintf(precFormat, Vprecision, xminPART, yminPART, zminPART);
1586 fprintf(fd,
"origin %s\n",precFormat);
1588 sprintf(precFormat, Vprecision, hx, 0.0, 0.0);
1589 fprintf(fd,
"delta %s\n",precFormat);
1591 sprintf(precFormat, Vprecision, 0.0, hy, 0.0);
1592 fprintf(fd,
"delta %s\n", precFormat);
1594 sprintf(precFormat, Vprecision, 0.0, 0.0, hzed);
1595 fprintf(fd,
"delta %s\n", precFormat);
1598 fprintf(fd,
"object 2 class gridconnections counts %d %d %d\n", nxPART, nyPART, nzPART);
1601 fprintf(fd,
"object 3 class array type double rank 0 items %d binary data follows\n",(nxPART*nyPART*nzPART));
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);
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");
1633 printf(
"Vgrid_writeDXBIN: Writing comments for %s format.\n",iofmt);
1634 fprintf(fd,
"# Data from %s\n", PACKAGE_STRING);
1636 fprintf(fd,
"# %s\n", title);
1637 fprintf(fd,
"# \n");
1640 fprintf(fd,
"object 1 class gridpositions counts %d %d %d\n", nx, ny, nz);
1642 sprintf(precFormat, Vprecision, xmin, ymin, zmin);
1643 fprintf(fd,
"origin %s\n", precFormat);
1645 sprintf(precFormat, Vprecision, hx, 0.0, 0.0);
1646 fprintf(fd,
"delta %s\n", precFormat);
1648 sprintf(precFormat, Vprecision, 0.0, hy, 0.0);
1649 fprintf(fd,
"delta %s\n", precFormat);
1651 sprintf(precFormat, Vprecision, 0.0, 0.0, hzed);
1652 fprintf(fd,
"delta %s\n", precFormat);
1655 fprintf(fd,
"object 2 class gridconnections counts %d %d %d\n", nx, ny, nz);
1658 fprintf(fd,
"object 3 class array type double rank 0 items %d binary data follows\n", (nx*ny*nz));
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);
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");
1692VPUBLIC void Vgrid_writeUHBD(Vgrid *thee, const char *iodev, const char *iofmt,
1693 const char *thost, const char *fname, char *title, double *pvec) {
1695 size_t u, icol, i, j, k;
1696 size_t gotit, nx, ny, nz;
1697 double xmin, ymin, zmin, hzed, hy, hx;
1700 if (thee == VNULL) {
1701 Vnm_print(2, "Vgrid_writeUHBD: Error -- got VNULL thee!\n");
1704 if (!(thee->ctordata || thee->readdata)) {
1705 Vnm_print(2, "Vgrid_writeUHBD: Error -- no data available!\n");
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 \
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",
1723 if (Vio_connect(sock, 0) < 0) {
1724 Vnm_print(2,
"Vgrid_writeUHBD: Problem connecting virtual socket %s\n",
1742 if (pvec != VNULL) {
1744 for (i=0; i<(nx*ny*nz); i++) {
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");
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,
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);
1768 for (k=0; k<nz; k++) {
1769 Vio_printf(sock,
"\n%7d%7d%7d\n", k+1, thee->nx, thee->ny);
1771 for (j=0; j<ny; j++) {
1772 for (i=0; i<nx; i++) {
1773 u = k*(nx)*(ny)+j*(nx)+i;
1775 Vio_printf(sock,
" %12.5e", thee->data[u]);
1778 Vio_printf(sock,
"\n");
1783 if (icol != 0) Vio_printf(sock,
"\n");
1786 Vio_connectFree(sock);