62 #if !defined(VINLINE_VGRID) 64 return Vmem_bytes(thee->
mem);
67 #define IJK(i,j,k) (((k)*(nx)*(ny))+((j)*(nx))+(i)) 69 #if defined(_WIN32) && (_MSC_VER < 1800) 79 VPRIVATE
double Vcompare;
80 VPRIVATE
char Vprecision[26];
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)]);
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);
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;
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",
811 const char *thost,
const char *fname) {
813 size_t i, j, k, itmp, u;
815 char tok[VMAX_BUFSIZE];
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)); }
828 FILE *fd = fopen(fname,
"rb");
830 printf(
"Vgrid_readDXBIN: Problem opening file %s\n",fname);
853 fgets(tok, VMAX_BUFSIZE, fd);
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");
863 printf(
"Vgrid_readDXBIN: Grid dimensions %d x %d x %d grid\n",thee->
nx, thee->
ny, thee->
nz);
893 if(fgets(tok,VMAX_BUFSIZE, fd) == NULL){
894 printf(
"Vgrid_readDXBIN: unexpected end of file.\n");
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");
903 printf(
"Vgrid_readDXBIN: Grid origin = (%g %g %g)\n",thee->
xmin, thee->
ymin, thee->
zmin);
906 if(fgets(tok,VMAX_BUFSIZE, fd) == NULL){
907 printf(
"Vgrid_readDXBIN: unexpected end of file.\n");
911 if(sscanf(tok,
"delta %lf %lf %lf",&(thee->
hx),&dtmp,&dtmp2)!=3){
912 printf(
"Vgrid_readDXBIN: Failed to read delta x data.\n");
917 if(fgets(tok,VMAX_BUFSIZE, fd) == NULL){
918 printf(
"Vgrid_readDXBIN: Unexpected end of file\n");
922 if(sscanf(tok,
"delta %lf %lf %lf",&dtmp,&(thee->
hy),&dtmp2)!=3){
923 printf(
"Vgrid_readDXBIN: Failed to read delta y data.\n");
928 if(fgets(tok,VMAX_BUFSIZE, fd) == NULL){
929 printf(
"Vgrid_readDXBIN: Unexpected end of file.\n");
933 if(sscanf(tok,
"delta %lf %lf %lf",&dtmp,&dtmp2,&(thee->
hzed))!=3){
934 printf(
"Vgrid_readDXBIN: Failed to read delta z data.\n");
938 printf(
"Vgrid_readDXBIN: Grid spacings = (%g, %g, %g)\n",thee->
hx, thee->
hy, thee->
hzed);
941 if(fgets(tok,VMAX_BUFSIZE, fd) == NULL){
942 printf(
"Vgrid_readDXBIN: Unexpected end of file.\n");
948 if(fgets(tok,VMAX_BUFSIZE, fd) == NULL){
949 printf(
"Vgrid_readDXBIN: Unexpected end of file.\n");
953 if(strstr(tok,
"binary")){
957 printf(
"Vgrid_readDXBIN: Binary tag not found. Will continue to try to read binary data.");
960 u = (size_t)thee->
nx * thee->
ny * thee->
nz;
961 int tot = thee->
nx * thee->
ny * thee->
nz;
964 printf(
"Vgrid_readDXBIN: allocating %d x %d x %d doubled for storage\n", thee->
nx, thee->
ny, thee->
nz);
966 thee->
data = (
double *)malloc(tot*
sizeof(
double));
968 if(thee->
data == NULL){
969 printf(
"Vgrid_readDXBIN: Unable to allocate space for data!\n");
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;
983 printf(
"Vgrid_readDXBIN: Failed to read doubles.\n");
992 printf(
"Vgrid_readDXBIN: Read double = %d not equal to items = %d\n",counter, tot);
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");
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",
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 \ 1361 data 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 \ 1422 data 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);
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");
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 \ 1753 will 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);
1796 if (thee == VNULL) {
1797 Vnm_print(2,
"Vgrid_integrate: Got VNULL thee!\n");
1807 for (k=0; k<nz; k++) {
1809 if ((k==0) || (k==(nz-1))) w = w * 0.5;
1810 for (j=0; j<ny; j++) {
1812 if ((j==0) || (j==(ny-1))) w = w * 0.5;
1813 for (i=0; i<nx; i++) {
1815 if ((i==0) || (i==(nx-1))) w = w * 0.5;
1816 sum = sum + w*(thee->
data[IJK(i,j,k)]);
1821 sum = sum*(thee->
hx)*(thee->
hy)*(thee->
hzed);
1834 if (thee == VNULL) {
1835 Vnm_print(2,
"Vgrid_normL1: Got VNULL thee!\n");
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)]);
1852 sum = sum*(thee->
hx)*(thee->
hy)*(thee->
hzed);
1864 if (thee == VNULL) {
1865 Vnm_print(2,
"Vgrid_normL2: Got VNULL thee!\n");
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)]);
1882 sum = sum*(thee->
hx)*(thee->
hy)*(thee->
hzed);
1892 double pt[3], grad[3], sum, hx, hy, hzed, xmin, ymin, zmin;
1894 if (thee == VNULL) {
1895 Vnm_print(2,
"Vgrid_seminormH1: Got VNULL thee!\n");
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;
1917 for (d=0; d<3; d++) sum = sum + VSQR(grad[d]);
1922 sum = sum*(hx)*(hy)*(hzed);
1924 if (VABS(sum) < VSMALL) sum = 0.0;
1925 else sum = VSQRT(sum);
1935 if (thee == VNULL) {
1936 Vnm_print(2,
"Vgrid_normH1: Got VNULL thee!\n");
1949 int nx, ny, nz, gotval;
1952 if (thee == VNULL) {
1953 Vnm_print(2,
"Vgrid_normLinf: Got VNULL thee!\n");
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)]);
1971 if (val > sum) sum = val;
VPUBLIC double Vgrid_normH1(Vgrid *thee)
Get the norm (or energy norm) of the data. This returns the integral: .
VPUBLIC double Vgrid_normL1(Vgrid *thee)
Get the norm of the data. This returns the integral: .
VPUBLIC double Vgrid_normLinf(Vgrid *thee)
Get the norm of the data. This returns the integral: .
Electrostatic potential oracle for Cartesian mesh data.
VPRIVATE char * MCcommChars
Comment characters for socket reads.
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.
Potential oracle for Cartesian mesh data.
VPUBLIC double Vgrid_integrate(Vgrid *thee)
Get the integral of the data.
#define VEMBED(rctag)
Allows embedding of RCS ID tags in object files.
VPUBLIC int Vstring_strcasecmp(const char *s1, const char *s2)
Case-insensitive string comparison (BSD standard)
VPUBLIC int Vgrid_gradient(Vgrid *thee, double pt[3], double grad[3])
Get first derivative values at a point.
VPUBLIC unsigned long int Vgrid_memChk(Vgrid *thee)
Return the memory used by this structure (and its contents) in bytes.
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.
VPUBLIC double Vgrid_normL2(Vgrid *thee)
Get the norm of the data. This returns the integral: .
VPUBLIC double Vgrid_seminormH1(Vgrid *thee)
Get the semi-norm of the data. This returns the integral: .
VPUBLIC int Vgrid_value(Vgrid *thee, double pt[3], double *value)
Get potential value (from mesh or approximation) at a point.
VPRIVATE char * MCwhiteChars
Whitespace characters for socket reads.