61 #if !defined(VINLINE_VGRID) 63 return Vmem_bytes(thee->
mem);
66 #define IJK(i,j,k) (((k)*(nx)*(ny))+((j)*(nx))+(i)) 70 VPRIVATE
double Vcompare;
71 VPRIVATE
char Vprecision[26];
219 if ((ihi<nx) && (jhi<ny) && (khi<nz) &&
220 (ilo>=0) && (jlo>=0) && (klo>=0)) {
222 dx = ifloat - (double)(ilo);
223 dy = jfloat - (double)(jlo);
224 dz = kfloat - (double)(klo);
225 u = dx *dy *dz *(thee->data[IJK(ihi,jhi,khi)])
226 + dx *(1.0-dy)*dz *(thee->data[IJK(ihi,jlo,khi)])
227 + dx *dy *(1.0-dz)*(thee->data[IJK(ihi,jhi,klo)])
228 + dx *(1.0-dy)*(1.0-dz)*(thee->data[IJK(ihi,jlo,klo)])
229 + (1.0-dx)*dy *dz *(thee->data[IJK(ilo,jhi,khi)])
230 + (1.0-dx)*(1.0-dy)*dz *(thee->data[IJK(ilo,jlo,khi)])
231 + (1.0-dx)*dy *(1.0-dz)*(thee->data[IJK(ilo,jhi,klo)])
232 + (1.0-dx)*(1.0-dy)*(1.0-dz)*(thee->data[IJK(ilo,jlo,klo)]);
237 Vnm_print(2,
"Vgrid_value: Got NaN!\n");
238 Vnm_print(2,
"Vgrid_value: (x, y, z) = (%4.3f, %4.3f, %4.3f)\n",
239 pt[0], pt[1], pt[2]);
240 Vnm_print(2,
"Vgrid_value: (ihi, jhi, khi) = (%d, %d, %d)\n",
242 Vnm_print(2,
"Vgrid_value: (ilo, jlo, klo) = (%d, %d, %d)\n",
244 Vnm_print(2,
"Vgrid_value: (nx, ny, nz) = (%d, %d, %d)\n",
246 Vnm_print(2,
"Vgrid_value: (dx, dy, dz) = (%4.3f, %4.3f, %4.3f)\n",
248 Vnm_print(2,
"Vgrid_value: data[IJK(ihi,jhi,khi)] = %g\n",
249 thee->data[IJK(ihi,jhi,khi)]);
250 Vnm_print(2,
"Vgrid_value: data[IJK(ihi,jlo,khi)] = %g\n",
251 thee->data[IJK(ihi,jlo,khi)]);
252 Vnm_print(2,
"Vgrid_value: data[IJK(ihi,jhi,klo)] = %g\n",
253 thee->data[IJK(ihi,jhi,klo)]);
254 Vnm_print(2,
"Vgrid_value: data[IJK(ihi,jlo,klo)] = %g\n",
255 thee->data[IJK(ihi,jlo,klo)]);
256 Vnm_print(2,
"Vgrid_value: data[IJK(ilo,jhi,khi)] = %g\n",
257 thee->data[IJK(ilo,jhi,khi)]);
258 Vnm_print(2,
"Vgrid_value: data[IJK(ilo,jlo,khi)] = %g\n",
259 thee->data[IJK(ilo,jlo,khi)]);
260 Vnm_print(2,
"Vgrid_value: data[IJK(ilo,jhi,klo)] = %g\n",
261 thee->data[IJK(ilo,jhi,klo)]);
262 Vnm_print(2,
"Vgrid_value: data[IJK(ilo,jlo,klo)] = %g\n",
263 thee->data[IJK(ilo,jlo,klo)]);
316 testpt[0] = pt[0] - hx;
318 testpt[0] = pt[0] + hx;
322 dxx = (uright - 2*umid + uleft)/(hx*hx);
326 testpt[1] = pt[1] - hy;
328 testpt[1] = pt[1] + hy;
332 dyy = (uright - 2*umid + uleft)/(hy*hy);
336 testpt[2] = pt[2] - hzed;
338 testpt[2] = pt[2] + hzed;
341 dzz = (uright - 2*umid + uleft)/(hzed*hzed);
346 curv = ( curv > fabs(dyy) ) ? curv : fabs(dyy);
347 curv = ( curv > fabs(dzz) ) ? curv : fabs(dzz);
348 }
else if ( cflag == 1 ) {
349 curv = (dxx + dyy + dzz)/3.0;
351 Vnm_print(2,
"Vgrid_curvature: support for cflag = %d not available!\n", cflag);
392 testpt[0] = pt[0] - hx;
393 if (
Vgrid_value( thee, testpt, &uleft)) haveleft = 1;
395 testpt[0] = pt[0] + hx;
396 if (
Vgrid_value( thee, testpt, &uright)) haveright = 1;
398 if (haveright && haveleft) grad[0] = (uright - uleft)/(2*hx);
399 else if (haveright) grad[0] = (uright - umid)/hx;
400 else if (haveleft) grad[0] = (umid - uleft)/hx;
408 testpt[1] = pt[1] - hy;
409 if (
Vgrid_value( thee, testpt, &uleft)) haveleft = 1;
411 testpt[1] = pt[1] + hy;
412 if (
Vgrid_value( thee, testpt, &uright)) haveright = 1;
414 if (haveright && haveleft) grad[1] = (uright - uleft)/(2*hy);
415 else if (haveright) grad[1] = (uright - umid)/hy;
416 else if (haveleft) grad[1] = (umid - uleft)/hy;
424 testpt[2] = pt[2] - hzed;
425 if (
Vgrid_value( thee, testpt, &uleft)) haveleft = 1;
427 testpt[2] = pt[2] + hzed;
428 if (
Vgrid_value( thee, testpt, &uright)) haveright = 1;
430 if (haveright && haveleft) grad[2] = (uright - uleft)/(2*hzed);
431 else if (haveright) grad[2] = (uright - umid)/hzed;
432 else if (haveleft) grad[2] = (umid - uleft)/hzed;
465 if (thee->data != VNULL) {
466 Vnm_print(1,
"%s: destroying existing data!\n", __func__);
467 Vmem_free(thee->mem, thee->nx * thee->ny * thee->nz,
sizeof(
double),
468 (
void **)&(thee->data));
474 infile = gzopen(fname,
"rb");
475 if (infile == Z_NULL) {
476 Vnm_print(2,
"%s: Problem opening compressed file %s\n", __func__, fname);
486 if(gzgets(infile, line, VMAX_ARGLEN) == Z_NULL){
491 if(strncmp(line,
"#", 1) == 0)
continue;
492 if(line[0] ==
'\n')
continue;
496 sscanf(line,
"object 1 class gridpositions counts %d %d %d",
497 &(thee->nx),&(thee->ny),&(thee->nz));
500 sscanf(line,
"origin %lf %lf %lf",
501 &(thee->xmin),&(thee->ymin),&(thee->zmin));
506 sscanf(line,
"delta %lf %lf %lf",&dtmp1,&dtmp2,&dtmp3);
519 Vnm_print(0,
"%s: allocating %d x %d x %d doubles for storage\n",
520 __func__, thee->nx, thee->ny, thee->nz);
521 len = thee->nx * thee->ny * thee->nz;
524 thee->data = Vmem_malloc(thee->mem, len,
sizeof(
double));
525 if (thee->data == VNULL) {
526 Vnm_print(2,
"%s: Unable to allocate space for data!\n", __func__);
534 temp = (
double *)malloc(len * (2 *
sizeof(
double)));
536 for (i = 0; i < len; i += 3){
537 memset(&line, 0,
sizeof(line));
538 gzgets(infile, line, VMAX_ARGLEN);
539 sscanf(line,
"%lf %lf %lf", &temp[i], &temp[i+1], &temp[i+2]);
544 for (i=0; i<thee->nx; i++) {
545 for (j=0; j<thee->ny; j++) {
546 for (k=0; k<thee->nz; k++) {
547 u = k*(thee->nx)*(thee->ny)+j*(thee->nx)+i;
548 (thee->data)[u] = temp[incr++];
554 thee->xmax = thee->xmin + (thee->nx-1)*thee->hx;
555 thee->ymax = thee->ymin + (thee->ny-1)*thee->hy;
556 thee->zmax = thee->zmin + (thee->nz-1)*thee->hzed;
563 Vnm_print(0,
"WARNING\n");
564 Vnm_print(0,
"Vgrid_readGZ: gzip read/write support is disabled in this build\n");
565 Vnm_print(0,
"Vgrid_readGZ: configure and compile without the --disable-zlib flag.\n");
566 Vnm_print(0,
"WARNING\n");
588 char tok[VMAX_BUFSIZE];
592 if (thee->
data != VNULL) {
593 Vnm_print(1,
"Vgrid_readDX: destroying existing data!\n");
594 Vmem_free(thee->
mem, (thee->
nx*thee->
ny*thee->
nz),
sizeof(
double),
595 (
void **)&(thee->
data)); }
600 sock = Vio_ctor(iodev,iofmt,thost,fname,
"r");
602 Vnm_print(2,
"Vgrid_readDX: Problem opening virtual socket %s\n",
606 if (Vio_accept(sock, 0) < 0) {
607 Vnm_print(2,
"Vgrid_readDX: Problem accepting virtual socket %s\n",
612 Vio_setWhiteChars(sock, MCwhiteChars);
613 Vio_setCommChars(sock, MCcommChars);
617 VJMPERR2(1 == Vio_scanf(sock,
"%s", tok));
618 VJMPERR1(!strcmp(tok,
"object"));
620 VJMPERR2(1 == Vio_scanf(sock,
"%d", &itmp));
622 VJMPERR2(1 == Vio_scanf(sock,
"%s", tok));
623 VJMPERR1(!strcmp(tok,
"class"));
625 VJMPERR2(1 == Vio_scanf(sock,
"%s", tok));
626 VJMPERR1(!strcmp(tok,
"gridpositions"));
628 VJMPERR2(1 == Vio_scanf(sock,
"%s", tok));
629 VJMPERR1(!strcmp(tok,
"counts"));
631 VJMPERR2(1 == Vio_scanf(sock,
"%s", tok));
632 VJMPERR1(1 == sscanf(tok,
"%d", &(thee->
nx)));
634 VJMPERR2(1 == Vio_scanf(sock,
"%s", tok));
635 VJMPERR1(1 == sscanf(tok,
"%d", &(thee->
ny)));
637 VJMPERR2(1 == Vio_scanf(sock,
"%s", tok));
638 VJMPERR1(1 == sscanf(tok,
"%d", &(thee->
nz)));
639 Vnm_print(0,
"Vgrid_readDX: Grid dimensions %d x %d x %d grid\n",
640 thee->
nx, thee->
ny, thee->
nz);
642 VJMPERR2(1 == Vio_scanf(sock,
"%s", tok));
643 VJMPERR1(!strcmp(tok,
"origin"));
645 VJMPERR2(1 == Vio_scanf(sock,
"%s", tok));
646 VJMPERR1(1 == sscanf(tok,
"%lf", &(thee->
xmin)));
648 VJMPERR2(1 == Vio_scanf(sock,
"%s", tok));
649 VJMPERR1(1 == sscanf(tok,
"%lf", &(thee->
ymin)));
651 VJMPERR2(1 == Vio_scanf(sock,
"%s", tok));
652 VJMPERR1(1 == sscanf(tok,
"%lf", &(thee->
zmin)));
653 Vnm_print(0,
"Vgrid_readDX: Grid origin = (%g, %g, %g)\n",
656 VJMPERR2(1 == Vio_scanf(sock,
"%s", tok));
657 VJMPERR1(!strcmp(tok,
"delta"));
659 VJMPERR2(1 == Vio_scanf(sock,
"%s", tok));
660 VJMPERR1(1 == sscanf(tok,
"%lf", &(thee->
hx)));
662 VJMPERR2(1 == Vio_scanf(sock,
"%s", tok));
663 VJMPERR1(1 == sscanf(tok,
"%lf", &dtmp));
664 VJMPERR1(dtmp == 0.0);
666 VJMPERR2(1 == Vio_scanf(sock,
"%s", tok));
667 VJMPERR1(1 == sscanf(tok,
"%lf", &dtmp));
668 VJMPERR1(dtmp == 0.0);
670 VJMPERR2(1 == Vio_scanf(sock,
"%s", tok));
671 VJMPERR1(!strcmp(tok,
"delta"));
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(1 == sscanf(tok,
"%lf", &(thee->
hy)));
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(!strcmp(tok,
"delta"));
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(1 == sscanf(tok,
"%lf", &dtmp));
693 VJMPERR1(dtmp == 0.0);
695 VJMPERR2(1 == Vio_scanf(sock,
"%s", tok));
696 VJMPERR1(1 == sscanf(tok,
"%lf", &(thee->
hzed)));
697 Vnm_print(0,
"Vgrid_readDX: Grid spacings = (%g, %g, %g)\n",
700 VJMPERR2(1 == Vio_scanf(sock,
"%s", tok));
701 VJMPERR1(!strcmp(tok,
"object"));
703 VJMPERR2(1 == Vio_scanf(sock,
"%s", tok));
705 VJMPERR2(1 == Vio_scanf(sock,
"%s", tok));
706 VJMPERR1(!strcmp(tok,
"class"));
708 VJMPERR2(1 == Vio_scanf(sock,
"%s", tok));
709 VJMPERR1(!strcmp(tok,
"gridconnections"));
711 VJMPERR2(1 == Vio_scanf(sock,
"%s", tok));
712 VJMPERR1(!strcmp(tok,
"counts"));
714 VJMPERR2(1 == Vio_scanf(sock,
"%s", tok));
715 VJMPERR2(1 == Vio_scanf(sock,
"%s", tok));
716 VJMPERR2(1 == Vio_scanf(sock,
"%s", tok));
718 VJMPERR2(1 == Vio_scanf(sock,
"%s", tok));
719 VJMPERR1(!strcmp(tok,
"object"));
721 VJMPERR2(1 == Vio_scanf(sock,
"%s", tok));
723 VJMPERR2(1 == Vio_scanf(sock,
"%s", tok));
724 VJMPERR1(!strcmp(tok,
"class"));
726 VJMPERR2(1 == Vio_scanf(sock,
"%s", tok));
727 VJMPERR1(!strcmp(tok,
"array"));
729 VJMPERR2(1 == Vio_scanf(sock,
"%s", tok));
730 VJMPERR1(!strcmp(tok,
"type"));
732 VJMPERR2(1 == Vio_scanf(sock,
"%s", tok));
733 VJMPERR1(!strcmp(tok,
"double"));
735 VJMPERR2(1 == Vio_scanf(sock,
"%s", tok));
736 VJMPERR1(!strcmp(tok,
"rank"));
738 VJMPERR2(1 == Vio_scanf(sock,
"%s", tok));
740 VJMPERR2(1 == Vio_scanf(sock,
"%s", tok));
741 VJMPERR1(!strcmp(tok,
"items"));
743 VJMPERR2(1 == Vio_scanf(sock,
"%s", tok));
744 VJMPERR1(1 == sscanf(tok,
"%d", &itmp));
745 VJMPERR1(((thee->
nx)*(thee->
ny)*(thee->
nz)) == itmp);
747 VJMPERR2(1 == Vio_scanf(sock,
"%s", tok));
748 VJMPERR1(!strcmp(tok,
"data"));
750 VJMPERR2(1 == Vio_scanf(sock,
"%s", tok));
751 VJMPERR1(!strcmp(tok,
"follows"));
754 Vnm_print(0,
"Vgrid_readDX: allocating %d x %d x %d doubles for storage\n",
755 thee->
nx, thee->
ny, thee->
nz);
757 thee->
data = (
double*)Vmem_malloc(thee->
mem, (thee->
nx)*(thee->
ny)*(thee->
nz),
759 if (thee->
data == VNULL) {
760 Vnm_print(2,
"Vgrid_readDX: Unable to allocate space for data!\n");
764 for (i=0; i<thee->
nx; i++) {
765 for (j=0; j<thee->
ny; j++) {
766 for (k=0; k<thee->
nz; k++) {
767 u = k*(thee->
nx)*(thee->
ny)+j*(thee->
nx)+i;
768 VJMPERR2(1 == Vio_scanf(sock,
"%s", tok));
769 VJMPERR1(1 == sscanf(tok,
"%lf", &dtmp));
770 (thee->
data)[u] = dtmp;
781 Vio_acceptFree(sock);
788 Vnm_print(2,
"Vgrid_readDX: Format problem with input file <%s>\n",
794 Vnm_print(2,
"Vgrid_readDX: I/O problem with input file <%s>\n",
850 Vnm_print(0,
"Vgrid_writeGZ: Opening file...\n");
851 outfile = gzopen(fname,
"wb");
863 for (k=0; k<nz; k++) {
865 for (j=0; j<ny; j++) {
867 for (i=0; i<nx; i++) {
869 if (pvec[IJK(i,j,k)] > 0.0) {
870 if (x < xminPART) xminPART = x;
871 if (y < yminPART) yminPART = y;
872 if (z < zminPART) zminPART = z;
878 for (k=0; k<nz; k++) {
880 for (j=0; j<ny; j++) {
881 for (i=0; i<nx; i++) {
882 if (pvec[IJK(i,j,k)] > 0.0) {
892 for (j=0; j<ny; j++) {
894 for (k=0; k<nz; k++) {
895 for (i=0; i<nx; i++) {
896 if (pvec[IJK(i,j,k)] > 0.0) {
906 for (i=0; i<nx; i++) {
908 for (k=0; k<nz; k++) {
909 for (j=0; j<ny; j++) {
910 if (pvec[IJK(i,j,k)] > 0.0) {
920 if ((nxPART != nx) || (nyPART != ny) || (nzPART != nz)) {
921 Vnm_print(0,
"Vgrid_writeGZ: printing only subset of domain\n");
924 txyz = (nxPART*nyPART*nzPART);
944 "object 1 class gridpositions counts %i %i %i\n" \
945 "origin %12.6e %12.6e %12.6e\n" \
946 "delta %12.6e 0.000000e+00 0.000000e+00\n" \
947 "delta 0.000000e+00 %12.6e 0.000000e+00\n" \
948 "delta 0.000000e+00 0.000000e+00 %12.6e\n" \
949 "object 2 class gridconnections counts %i %i %i\n"\
950 "object 3 class array type double rank 0 items %i data follows\n",
951 PACKAGE_STRING,title,nx,ny,nz,txmin,tymin,tzmin,
952 hx,hy,hzed,nx,ny,nz,txyz);
953 gzwrite(outfile, header, strlen(header)*
sizeof(
char));
957 for (i=0; i<nx; i++) {
958 for (j=0; j<ny; j++) {
959 for (k=0; k<nz; k++) {
960 u = k*(nx)*(ny)+j*(nx)+i;
962 sprintf(line,
"%12.6e ", thee->data[u]);
963 gzwrite(outfile, line, strlen(line)*
sizeof(
char));
967 gzwrite(outfile, newline, strlen(newline)*
sizeof(
char));
974 char newline[] =
"\n";
975 gzwrite(outfile, newline, strlen(newline)*
sizeof(
char));
979 sprintf(footer,
"attribute \"dep\" string \"positions\"\n" \
980 "object \"regular positions regular connections\" class field\n" \
981 "component \"positions\" value 1\n" \
982 "component \"connections\" value 2\n" \
983 "component \"data\" value 3\n");
984 gzwrite(outfile, footer, strlen(footer)*
sizeof(
char));
989 Vnm_print(0,
"WARNING\n");
990 Vnm_print(0,
"Vgrid_readGZ: gzip read/write support is disabled in this build\n");
991 Vnm_print(0,
"Vgrid_readGZ: configure and compile without the --disable-zlib flag.\n");
992 Vnm_print(0,
"WARNING\n");
1034 Vnm_print(0,
"Vgrid_writeDX: Opening virtual socket...\n");
1035 sock = Vio_ctor(iodev,iofmt,thost,fname,
"w");
1036 if (sock == VNULL) {
1037 Vnm_print(2,
"Vgrid_writeDX: Problem opening virtual socket %s\n",
1041 if (Vio_connect(sock, 0) < 0) {
1042 Vnm_print(2,
"Vgrid_writeDX: Problem connecting virtual socket %s\n",
1047 Vio_setWhiteChars(sock, MCwhiteChars);
1048 Vio_setCommChars(sock, MCcommChars);
1050 Vnm_print(0,
"Vgrid_writeDX: Writing to virtual socket...\n");
1062 for (k=0; k<nz; k++) {
1064 for (j=0; j<ny; j++) {
1066 for (i=0; i<nx; i++) {
1068 if (pvec[IJK(i,j,k)] > 0.0) {
1069 if (x < xminPART) xminPART = x;
1070 if (y < yminPART) yminPART = y;
1071 if (z < zminPART) zminPART = z;
1077 for (k=0; k<nz; k++) {
1079 for (j=0; j<ny; j++) {
1080 for (i=0; i<nx; i++) {
1081 if (pvec[IJK(i,j,k)] > 0.0) {
1088 if (gotit) nzPART++;
1091 for (j=0; j<ny; j++) {
1093 for (k=0; k<nz; k++) {
1094 for (i=0; i<nx; i++) {
1095 if (pvec[IJK(i,j,k)] > 0.0) {
1102 if (gotit) nyPART++;
1105 for (i=0; i<nx; i++) {
1107 for (k=0; k<nz; k++) {
1108 for (j=0; j<ny; j++) {
1109 if (pvec[IJK(i,j,k)] > 0.0) {
1116 if (gotit) nxPART++;
1119 if ((nxPART != nx) || (nyPART != ny) || (nzPART != nz)) {
1120 Vnm_print(0,
"Vgrid_writeDX: printing only subset of domain\n");
1126 Vnm_print(0,
"Vgrid_writeDX: Skipping comments for XDR format.\n");
1128 Vnm_print(0,
"Vgrid_writeDX: Writing comments for %s format.\n",
1130 Vio_printf(sock,
"# Data from %s\n", PACKAGE_STRING);
1131 Vio_printf(sock,
"# \n");
1132 Vio_printf(sock,
"# %s\n", title);
1133 Vio_printf(sock,
"# \n");
1137 Vio_printf(sock,
"object 1 class gridpositions counts %d %d %d\n",
1138 nxPART, nyPART, nzPART);
1140 sprintf(precFormat, Vprecision, xminPART, yminPART, zminPART);
1141 Vio_printf(sock,
"origin %s\n", precFormat);
1142 sprintf(precFormat, Vprecision, hx, 0.0, 0.0);
1143 Vio_printf(sock,
"delta %s\n", precFormat);
1144 sprintf(precFormat, Vprecision, 0.0, hy, 0.0);
1145 Vio_printf(sock,
"delta %s\n", precFormat);
1146 sprintf(precFormat, Vprecision, 0.0, 0.0, hzed);
1147 Vio_printf(sock,
"delta %s\n", precFormat);
1150 Vio_printf(sock,
"object 2 class gridconnections counts %d %d %d\n",
1151 nxPART, nyPART, nzPART);
1154 Vio_printf(sock,
"object 3 class array type double rank 0 items %d \ 1155 data follows\n", (nxPART*nyPART*nzPART));
1157 for (i=0; i<nx; i++) {
1158 for (j=0; j<ny; j++) {
1159 for (k=0; k<nz; k++) {
1160 u = k*(nx)*(ny)+j*(nx)+i;
1161 if (pvec[u] > 0.0) {
1162 Vio_printf(sock,
"%12.6e ", thee->data[u]);
1166 Vio_printf(sock,
"\n");
1173 if (icol != 0) Vio_printf(sock,
"\n");
1176 Vio_printf(sock,
"attribute \"dep\" string \"positions\"\n");
1177 Vio_printf(sock,
"object \"regular positions regular connections\" \ 1179 Vio_printf(sock,
"component \"positions\" value 1\n");
1180 Vio_printf(sock,
"component \"connections\" value 2\n");
1181 Vio_printf(sock,
"component \"data\" value 3\n");
1186 Vnm_print(0,
"Vgrid_writeDX: Skipping comments for XDR format.\n");
1188 Vnm_print(0,
"Vgrid_writeDX: Writing comments for %s format.\n",
1190 Vio_printf(sock,
"# Data from %s\n", PACKAGE_STRING);
1191 Vio_printf(sock,
"# \n");
1192 Vio_printf(sock,
"# %s\n", title);
1193 Vio_printf(sock,
"# \n");
1198 Vio_printf(sock,
"object 1 class gridpositions counts %d %d %d\n",
1201 sprintf(precFormat, Vprecision, xmin, ymin, zmin);
1202 Vio_printf(sock,
"origin %s\n", precFormat);
1203 sprintf(precFormat, Vprecision, hx, 0.0, 0.0);
1204 Vio_printf(sock,
"delta %s\n", precFormat);
1205 sprintf(precFormat, Vprecision, 0.0, hy, 0.0);
1206 Vio_printf(sock,
"delta %s\n", precFormat);
1207 sprintf(precFormat, Vprecision, 0.0, 0.0, hzed);
1208 Vio_printf(sock,
"delta %s\n", precFormat);
1211 Vio_printf(sock,
"object 2 class gridconnections counts %d %d %d\n",
1215 Vio_printf(sock,
"object 3 class array type double rank 0 items %d \ 1216 data follows\n", (nx*ny*nz));
1218 for (i=0; i<nx; i++) {
1219 for (j=0; j<ny; j++) {
1220 for (k=0; k<nz; k++) {
1221 u = k*(nx)*(ny)+j*(nx)+i;
1222 Vio_printf(sock,
"%12.6e ", thee->data[u]);
1226 Vio_printf(sock,
"\n");
1231 if (icol != 0) Vio_printf(sock,
"\n");
1234 Vio_printf(sock,
"attribute \"dep\" string \"positions\"\n");
1235 Vio_printf(sock,
"object \"regular positions regular connections\" \ 1237 Vio_printf(sock,
"component \"positions\" value 1\n");
1238 Vio_printf(sock,
"component \"connections\" value 2\n");
1239 Vio_printf(sock,
"component \"data\" value 3\n");
1243 Vio_connectFree(sock);
1275 sock = Vio_ctor(iodev,iofmt,thost,fname,
"w");
1276 if (sock == VNULL) {
1277 Vnm_print(2,
"Vgrid_writeUHBD: Problem opening virtual socket %s\n",
1281 if (Vio_connect(sock, 0) < 0) {
1282 Vnm_print(2,
"Vgrid_writeUHBD: Problem connecting virtual socket %s\n",
1300 if (pvec != VNULL) {
1302 for (i=0; i<(nx*ny*nz); i++) {
1309 Vnm_print(2,
"Vgrid_writeUHBD: IGNORING PARTITION INFORMATION!\n");
1310 Vnm_print(2,
"Vgrid_writeUHBD: This means I/O from parallel runs \ 1311 will have significant overlap.\n");
1316 Vio_printf(sock,
"%72s\n", title);
1317 Vio_printf(sock,
"%12.5e%12.5e%7d%7d%7d%7d%7d\n", 1.0, 0.0, -1, 0,
1319 Vio_printf(sock,
"%7d%7d%7d%12.5e%12.5e%12.5e%12.5e\n", nx, ny, nz,
1320 hx, (xmin-hx), (ymin-hx), (zmin-hx));
1321 Vio_printf(sock,
"%12.5e%12.5e%12.5e%12.5e\n", 0.0, 0.0, 0.0, 0.0);
1322 Vio_printf(sock,
"%12.5e%12.5e%7d%7d", 0.0, 0.0, 0, 0);
1326 for (k=0; k<nz; k++) {
1327 Vio_printf(sock,
"\n%7d%7d%7d\n", k+1, thee->nx, thee->ny);
1329 for (j=0; j<ny; j++) {
1330 for (i=0; i<nx; i++) {
1331 u = k*(nx)*(ny)+j*(nx)+i;
1333 Vio_printf(sock,
" %12.5e", thee->data[u]);
1336 Vio_printf(sock,
"\n");
1341 if (icol != 0) Vio_printf(sock,
"\n");
1344 Vio_connectFree(sock);
1350 int i, j, k, nx, ny, nz;
1353 if (thee == VNULL) {
1354 Vnm_print(2,
"Vgrid_integrate: Got VNULL thee!\n");
1364 for (k=0; k<nz; k++) {
1366 if ((k==0) || (k==(nz-1))) w = w * 0.5;
1367 for (j=0; j<ny; j++) {
1369 if ((j==0) || (j==(ny-1))) w = w * 0.5;
1370 for (i=0; i<nx; i++) {
1372 if ((i==0) || (i==(nx-1))) w = w * 0.5;
1373 sum = sum + w*(thee->
data[IJK(i,j,k)]);
1378 sum = sum*(thee->
hx)*(thee->
hy)*(thee->
hzed);
1387 int i, j, k, nx, ny, nz;
1390 if (thee == VNULL) {
1391 Vnm_print(2,
"Vgrid_normL1: Got VNULL thee!\n");
1400 for (k=0; k<nz; k++) {
1401 for (j=0; j<ny; j++) {
1402 for (i=0; i<nx; i++) {
1403 sum = sum + VABS(thee->
data[IJK(i,j,k)]);
1408 sum = sum*(thee->
hx)*(thee->
hy)*(thee->
hzed);
1416 int i, j, k, nx, ny, nz;
1419 if (thee == VNULL) {
1420 Vnm_print(2,
"Vgrid_normL2: Got VNULL thee!\n");
1429 for (k=0; k<nz; k++) {
1430 for (j=0; j<ny; j++) {
1431 for (i=0; i<nx; i++) {
1432 sum = sum + VSQR(thee->
data[IJK(i,j,k)]);
1437 sum = sum*(thee->
hx)*(thee->
hy)*(thee->
hzed);
1445 int i, j, k, d, nx, ny, nz;
1446 double pt[3], grad[3], sum, hx, hy, hzed, xmin, ymin, zmin;
1448 if (thee == VNULL) {
1449 Vnm_print(2,
"Vgrid_seminormH1: Got VNULL thee!\n");
1464 for (k=0; k<nz; k++) {
1465 pt[2] = k*hzed + zmin;
1466 for (j=0; j<ny; j++) {
1467 pt[1] = j*hy + ymin;
1468 for (i=0; i<nx; i++) {
1469 pt[0] = i*hx + xmin;
1471 for (d=0; d<3; d++) sum = sum + VSQR(grad[d]);
1476 sum = sum*(hx)*(hy)*(hzed);
1478 if (VABS(sum) < VSMALL) sum = 0.0;
1479 else sum = VSQRT(sum);
1489 if (thee == VNULL) {
1490 Vnm_print(2,
"Vgrid_normH1: Got VNULL thee!\n");
1502 int i, j, k, nx, ny, nz, gotval;
1505 if (thee == VNULL) {
1506 Vnm_print(2,
"Vgrid_normLinf: Got VNULL thee!\n");
1516 for (k=0; k<nz; k++) {
1517 for (j=0; j<ny; j++) {
1518 for (i=0; i<nx; i++) {
1519 val = VABS(thee->
data[IJK(i,j,k)]);
1524 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.
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.