/* // d2p.c // created by F.Mizutani 960313 // d2p002 Plus minus iso value // d2p003 Data convert, orbital data support // d2p004 Moleculer data // d2p005 Auto output filename // d2p006 unix,dos support // d2p007 double,triple polygon // d2p011 BUG fix: normal vector // d2p012 Volume count 2001/08/02 */ /* unix header file */ #include #include #include #include #include #include #include #define DOSCOD 0x1111 #define UNXCOD 0x2222 #ifdef DOS #define DCODE DOSCOD #else #define DCODE UNXCOD #endif #define PCODE 11002 #define GCODE 10004 /* program code */ #define PMAX 32 /* parameters */ #define DIV 50 /* find neighbor polygon < step/DIV */ #define OBTMAX 64 #define FALSE 0 #define TRUE !FALSE float viso; /* iso vale */ float vvvv; float xmax, ymax, zmax, datamax; float xmin, ymin, zmin, datamin; typedef struct { float x; float y; float z; float data; int flag; }IDST; IDST *ids; typedef struct { float x; float y; float z; } XYZ; typedef struct { float x[3]; float y[3]; float z[3]; float vx[3]; float vy[3]; float vz[3]; float vnx, vny, vnz; } PGD; typedef struct { long mn; float mw; float x; float y; float z; } MOL; PGD *pgd; int pgdct; int pgdmx; long idsz, /* Mesh data size. */ idmx, /* Mesh size for x direction. */ idmy, /* Mesh size for y direction. */ idmz; /* Mesh size for z direction. */ float mbgx, mbgy, mbgz; float mstx, msty, mstz; float mdvx, mdvy, mdvz; long dtype; long oct; long obt[OBTMAX]; MOL *mol; long molct; int sfg; int dble = 1; int cmode; long incount[12]; long opencount = 0; void swapper(p, n) char *p; int n; { char d[4]; int i, j; for (j = 0; j < n; j++){ for (i = 0; i < 4; i++) d[3 - i] = p[i + 4 * j]; for (i = 0; i < 4; i++) p[i + 4 * j] = d[i]; } } void m_alloc(size) int size; { if ((ids = (IDST *)malloc(size * sizeof(IDST)))==NULL){ perror("ids malloc error "); exit(1); } pgdmx = size / 3; if ((pgd = (PGD *)malloc(pgdmx * sizeof(PGD))) == NULL){ perror("pgd malloc error "); exit(1); } } /* * read mesh data from input file * idsz : Mesh data size * idmx : Mesh size for x direction * idmy : Mesh size for y direction * idmz : Mesh size for z direction * ids : Mesh struct data */ void InputData(fp) FILE *fp; { long i, j, k, n, dcode, pcode; float x, y, z, data, *dt; static float smax = (float)(1.0e30); static float smin = - (float)(1.0e30); printf("ReadFile ... "); fread(&dcode, sizeof(long), 1, fp); if (dcode != DOSCOD && dcode != UNXCOD){ fclose(fp); printf("This data is not made by isoconv.\n"); fclose(fp); } printf("Done.\n"); if (dcode == DCODE) sfg = FALSE; else sfg = TRUE; fread(&pcode, sizeof(long), 1, fp); if (sfg) swapper((char *)&pcode, 1); if (pcode != PCODE){ fclose(fp); printf("This data is not made by isoconv.\n"); exit(1); } fread(&dtype, sizeof(long), 1, fp); fread(&datamax, sizeof(float), 1, fp); fread(&datamin, sizeof(float), 1, fp); if (sfg) swapper((char *)&dtype, 1); if (sfg) swapper((char *)&datamax, 1); if (sfg) swapper((char *)&datamin, 1); fread(&idmx, sizeof(long), 1, fp); fread(&idmy, sizeof(long), 1, fp); fread(&idmz, sizeof(long), 1, fp); if (sfg) swapper((char *)&idmx, 1); if (sfg) swapper((char *)&idmy, 1); if (sfg) swapper((char *)&idmz, 1); idsz = idmx * idmy * idmz; printf("Resolution X=%d Y=%d Z=%d TOTAL=%d\n",idmx,idmy,idmz,idsz); m_alloc(idsz); fread(&mbgx, sizeof(float), 1, fp); fread(&mbgy, sizeof(float), 1, fp); fread(&mbgz, sizeof(float), 1, fp); fread(&mstx, sizeof(float), 1, fp); fread(&msty, sizeof(float), 1, fp); fread(&mstz, sizeof(float), 1, fp); if (sfg) swapper((char *)&mbgx, 1); if (sfg) swapper((char *)&mbgy, 1); if (sfg) swapper((char *)&mbgz, 1); if (sfg) swapper((char *)&mstx, 1); if (sfg) swapper((char *)&msty, 1); if (sfg) swapper((char *)&mstz, 1); printf("Begin AXIS X=%+f Y=%+f Z=%+f\n",mbgx,mbgy,mbgz); printf("Step AXIS X=%+f Y=%+f Z=%+f\n",mstx,msty,mstz); mdvx = mstx / DIV; mdvy = msty / DIV; mdvz = mstz / DIV; fread(&molct, sizeof(long), 1, fp); if (sfg) swapper((char *)&molct, 1); mol = (MOL *)malloc(molct * sizeof(MOL)); fread(mol, sizeof(MOL), molct, fp); if (sfg) swapper((char *)mol, sizeof(MOL)/sizeof(long)); for (i = 0; i < molct; i++){ printf("%3d: No.%3d W=%+f X=%+f Y=%+f Z=%+f\n", i+1, mol[i].mn, mol[i].mw, mol[i].x, mol[i].y, mol[i].z); } if (dtype == 1){ printf("Orbital data\n"); fread(&oct, sizeof(long), 1, fp); if (sfg) swapper((char *)&oct, 1); fread(obt, sizeof(long), oct, fp); if (sfg) swapper((char *)obt, oct); } else printf("Density data\n"); dt = (float *)malloc(idmz * sizeof(float)); xmax = smin; ymax = smin; zmax = smin; /*datamax = smin;*/ xmin = smax; ymin = smax; zmin = smax; /*datamin = smax;*/ printf("X loop "); for(i = 0; i < idmx; i++){ printf("[%d]", i + 1); fflush(stdout); for (j = 0; j < idmy; j++){ if (fread(dt, sizeof(float), idmz, fp) != idmz){ printf("data read error "); exit(1); } if (sfg) swapper((char *)dt, idmz); for (k = 0; k < idmz; k++){ /*data = dt[k];*/ x = mbgx + mstx * i; y = mbgy + msty * j; z = mbgz + mstz * k; n = i + j * idmx + k * idmx * idmy; ids[n].x = x; ids[n].y = y; ids[n].z = z; ids[n].data = dt[k]; /* printf("%+f %+f %+f %+f\n",ids[i].x,ids[i].y,ids[i].z,ids[i].data); */ if (x > xmax) xmax = x; if (x < xmin) xmin = x; if (y > ymax) ymax = y; if (y < ymin) ymin = y; if (z > zmax) zmax = z; if (z < zmin) zmin = z; /*if (data > datamax) datamax = data;*/ /*if (data < datamin) datamin = data;*/ } } } printf("\n"); printf("x=%+e:%+e\n",xmin,xmax); printf("y=%+e:%+e\n",ymin,ymax); printf("z=%+e:%+e\n",zmin,zmax); printf("data=%+e:%+e\n",datamin,datamax); } void FlagSelect(f) int f; { int i, j, k, n; pgdct = 0; for(i = 0; i < idmx; i++){ for (j = 0; j < idmy; j++){ for (k = 0; k < idmz; k++){ n = i + j * idmx + k * idmx * idmy; ids[n].flag = 0; if (viso > 0){ if (ids[n].data > viso){ ids[n].flag = 1; incount[f]++; } } else { if (ids[n].data < viso){ ids[n].flag = 1; incount[f]++; } } if (ids[n].flag == 1){ if (i == 0 || i == idmx - 1) opencount++; if (j == 0 || j == idmy - 1) opencount++; if (k == 0 || k == idmz - 1) opencount++; } } } } } void NormalizeVector() { int i, j, k, l; int ct; int nb[100], xyz[100]; float length, x, y, z; for (i = 0; i < pgdct; i++){ if (i % 100 == 0){ printf("[%d]", i); fflush(stdout); } if (i % 1000 == 0){ printf("\n"); fflush(stdout); } for (j = 0; j < 3; j++){ if (pgd[i].vx[j] != 0.0 || pgd[i].vy[j] != 0.0 || pgd[i].vz[j] != 0.0) continue; ct = 1; x = pgd[i].vnx; y = pgd[i].vny; z = pgd[i].vnz; for (k = i + 1; k < pgdct; k++){ if (pgd[i].z[0] + mstz * 2 < pgd[k].z[0]) break; for (l = 0; l < 3; l++){ if (fabs(pgd[i].x[j] - pgd[k].x[l]) < mdvx && fabs(pgd[i].y[j] - pgd[k].y[l]) < mdvy && fabs(pgd[i].z[j] - pgd[k].z[l]) < mdvz && ct < 100){ nb[ct] = k; xyz[ct] = l; ct++; x += pgd[k].vnx; y += pgd[k].vny; z += pgd[k].vnz; break; } } } length = sqrt(x*x + y*y + z*z); pgd[i].vx[j] = (float)(x / length); pgd[i].vy[j] = (float)(y / length); pgd[i].vz[j] = (float)(z / length); for (k = 1; k < ct; k++){ pgd[nb[k]].vx[xyz[k]] = pgd[i].vx[j]; pgd[nb[k]].vy[xyz[k]] = pgd[i].vy[j]; pgd[nb[k]].vz[xyz[k]] = pgd[i].vz[j]; } } } printf("\n"); } XYZ CalcVertex(n1, n2) int n1; int n2; { float a; XYZ v; a = (float)(viso - ids[n1].data)/(ids[n2].data - ids[n1].data); v.x = ids[n1].x + a*(ids[n2].x - ids[n1].x); v.y = ids[n1].y + a*(ids[n2].y - ids[n1].y); v.z = ids[n1].z + a*(ids[n2].z - ids[n1].z); return v; } XYZ CalcNormal(v) XYZ v[3]; { int i; XYZ e[2], n; float length; float x, y, z; for (i = 0; i < 2; i++){ e[i].x = v[i + 1].x - v[0].x; e[i].y = v[i + 1].y - v[0].y; e[i].z = v[i + 1].z - v[0].z; } x = e[0].y * e[1].z - e[0].z * e[1].y; y = e[0].z * e[1].x - e[0].x * e[1].z; z = e[0].x * e[1].y - e[0].y * e[1].x; length = (float)sqrt(x*x + y*y + z*z); n.x = x / length; n.y = y / length; n.z = z / length; return n; } void SetTriangle(v) XYZ v[3]; { XYZ n; int i; float ck; n = CalcNormal(v); ck = n.x * n.x + n.y * n.y + n.z * n.z; if (0.95 < ck && ck < 1.05){ if (pgdct < pgdmx){ pgd[pgdct].vnx = n.x; pgd[pgdct].vny = n.y; pgd[pgdct].vnz = n.z; for (i = 0; i < 3; i++){ pgd[pgdct].x[i] = v[i].x; pgd[pgdct].y[i] = v[i].y; pgd[pgdct].z[i] = v[i].z; pgd[pgdct].vx[i] = 0.0; pgd[pgdct].vy[i] = 0.0; pgd[pgdct].vz[i] = 0.0; } pgdct++; } } /* else printf("E:%+f %+f %+f\n",n.x,n.y,n.z);*/ else printf("!"); } void DisplayTriangle(wnode) int *wnode; { XYZ v[3]; v[0] = CalcVertex(wnode[0],wnode[1]); v[2] = CalcVertex(wnode[2],wnode[3]); v[1] = CalcVertex(wnode[4],wnode[5]); SetTriangle(v); } void DisplayCorner(wnode) int wnode[4]; { XYZ v[3]; v[0] = CalcVertex(wnode[0],wnode[1]); v[2] = CalcVertex(wnode[0],wnode[2]); v[1] = CalcVertex(wnode[0],wnode[3]); SetTriangle(v); } void Display2Tri(wnode) int *wnode; { XYZ v[3]; v[0] = CalcVertex(wnode[0],wnode[1]); v[2] = CalcVertex(wnode[0],wnode[2]); v[1] = CalcVertex(wnode[3],wnode[5]); SetTriangle(v); v[0] = CalcVertex(wnode[0],wnode[2]); v[2] = CalcVertex(wnode[3],wnode[4]); v[1] = CalcVertex(wnode[3],wnode[5]); SetTriangle(v); } void SetCornerOn(cnode, i, j) int cnode[8][8]; int i,j; { int wnode[4]; if (j < 4){ wnode[0] = cnode[i][j]; wnode[1] = cnode[i][(j+1)%4]; wnode[2] = cnode[i][j+4]; wnode[3] = cnode[i][(j+3)%4]; } else { wnode[0] = cnode[i][j]; wnode[1] = cnode[i][4+(j+1)%4]; wnode[2] = cnode[i][4+(j+3)%4]; wnode[3] = cnode[i][j-4]; } DisplayCorner(wnode); } void SetCornerOff(cnode, i, j) int cnode[8][8]; int i,j; { int wnode[4]; if (j < 4){ wnode[0] = cnode[i][j]; wnode[1] = cnode[i][(j+1)%4]; wnode[2] = cnode[i][(j+3)%4]; wnode[3] = cnode[i][j+4]; } else { wnode[0] = cnode[i][j]; wnode[1] = cnode[i][4+(j+1)%4]; wnode[2] = cnode[i][j-4]; wnode[3] = cnode[i][4+(j+3)%4]; } DisplayCorner(wnode); } /* * １のパターン */ void Flag1(cnode) int cnode[8][8]; { int i; for (i = 0; i < 8; i++){ if(ids[cnode[i][0]].flag == 1){ SetCornerOn(cnode,i,0); break; } } } /* * １のパターンの反転パターン */ void Flag7(cnode) int cnode[8][8]; { int i; for (i = 0; i < 8; i++){ if(ids[cnode[i][0]].flag == 0){ SetCornerOff(cnode,i,0); break; } } } /* * 選ばれた２点が、格子の底面にある場合 */ void F2_Patern2_1(cnode,cyc,p0,p1) int cnode[8][8]; int cyc,p0,p1; { int wnode[6]; wnode[0] = cnode[cyc][p0]; wnode[1] = cnode[cyc][p0+4]; wnode[2] = cnode[cyc][(p0+3)%4]; wnode[3] = cnode[cyc][p1]; wnode[4] = cnode[cyc][(p1+1)%4]; wnode[5] = cnode[cyc][p1+4]; Display2Tri(wnode); } /* * 選ばれた２点が、格子の上面にある場合 */ void F2_Patern2_2(cnode, cyc, p0, p1) int cnode[8][8]; int cyc,p0,p1; { int wnode[6]; wnode[0] = cnode[cyc][p0]; wnode[1] = cnode[cyc][4+(p0+3)%4]; wnode[2] = cnode[cyc][p0-4]; wnode[3] = cnode[cyc][p1]; wnode[4] = cnode[cyc][p1-4]; wnode[5] = cnode[cyc][4+(p1+1)%4]; Display2Tri(wnode); } /* * 選ばれた２点が、格子の側面にある場合 */ void F2_Patern2_3(cnode, cyc, p0, p1) int cnode[8][8]; int cyc,p0,p1; { int wnode[6]; wnode[0] = cnode[cyc][p0]; wnode[1] = cnode[cyc][(p0+3)%4]; wnode[2] = cnode[cyc][(p0+1)%4]; wnode[3] = cnode[cyc][p1]; wnode[4] = cnode[cyc][4+(p1+1)%4]; wnode[5] = cnode[cyc][4+(p1+3)%4]; Display2Tri(wnode); } /* * ２、３、４のパターン */ void Flag2(cnode) int cnode[8][8]; { int i, j, cyc; int pnode[2]; for (i = 0; i < 8; i++){ if (ids[cnode[i][0]].flag == 1){ cyc = i; pnode[0] = 0; for(j=1;j<8;j++){ if (ids[cnode[i][j]].flag == 1){ pnode[1] = j; break; } } break; } } /* printf("Flag2:cyc=%d,0-%d \n",cyc,pnode[1]); */ switch(pnode[1]) { case 1: F2_Patern2_1(cnode,cyc,pnode[0],pnode[1]); break; case 3: F2_Patern2_1(cnode,cyc,pnode[1],pnode[0]); break; case 4: F2_Patern2_3(cnode,cyc,pnode[0],pnode[1]); break; default: SetCornerOn(cnode,cyc,pnode[0]); SetCornerOn(cnode,cyc,pnode[1]); printf("[WARNING PAT2]"); } } /* * 選ばれた２点が、格子の底面にある場合 */ void F6_Patern2_1(cnode, cyc, p0, p1) int cnode[8][8]; int cyc,p0,p1; { int wnode[6]; wnode[0] = cnode[cyc][p0]; wnode[2] = cnode[cyc][p0+4]; wnode[1] = cnode[cyc][(p0+3)%4]; wnode[3] = cnode[cyc][p1]; wnode[5] = cnode[cyc][(p1+1)%4]; wnode[4] = cnode[cyc][p1+4]; Display2Tri(wnode); } /* * 選ばれた２点が、格子の上面にある場合 */ void F6_Patern2_2(cnode, cyc, p0, p1) int cnode[8][8]; int cyc,p0,p1; { int wnode[6]; wnode[0] = cnode[cyc][p0]; wnode[2] = cnode[cyc][4+(p0+3)%4]; wnode[1] = cnode[cyc][p0-4]; wnode[3] = cnode[cyc][p1]; wnode[5] = cnode[cyc][p1-4]; wnode[4] = cnode[cyc][4+(p1+1)%4]; Display2Tri(wnode); } /* * 選ばれた２点が、格子の側面にある場合 */ void F6_Patern2_3(cnode, cyc, p0, p1) int cnode[8][8]; int cyc,p0,p1; { int wnode[6]; wnode[0] = cnode[cyc][p0]; wnode[2] = cnode[cyc][(p0+3)%4]; wnode[1] = cnode[cyc][(p0+1)%4]; wnode[3] = cnode[cyc][p1]; wnode[5] = cnode[cyc][4+(p1+1)%4]; wnode[4] = cnode[cyc][4+(p1+3)%4]; Display2Tri(wnode); } /* * ２、３、４のパターンの反転パターン */ void Flag6(cnode) int cnode[8][8]; { int i, j, cyc; int pnode[2]; for (i = 0; i < 8; i++){ if (ids[cnode[i][0]].flag == 0){ cyc = i; pnode[0] = 0; for(j=1;j<8;j++){ if (ids[cnode[i][j]].flag == 0){ pnode[1] = j; break; } } break; } } /* printf("Flag6:cyc=%d,0-%d \n",cyc,pnode[1]); */ switch(pnode[1]) { case 1: F6_Patern2_1(cnode,cyc,pnode[0],pnode[1]); break; case 3: F6_Patern2_1(cnode,cyc,pnode[1],pnode[0]); break; case 4: F6_Patern2_3(cnode,cyc,pnode[0],pnode[1]); break; default: SetCornerOff(cnode,cyc,pnode[0]); SetCornerOff(cnode,cyc,pnode[1]); printf("[WARNING PAT6]"); } } /* * 選ばれた6点で三角形を描く */ void F3_Patern5(cnode, cyc, p1, p2, p3, p4, p5, p6) int cnode[8][8]; int cyc,p1,p2,p3,p4,p5,p6; { int wnode[6]; wnode[0] = cnode[cyc][p1]; wnode[1] = cnode[cyc][p2]; wnode[2] = cnode[cyc][p3]; wnode[3] = cnode[cyc][p4]; wnode[4] = cnode[cyc][p5]; wnode[5] = cnode[cyc][p6]; DisplayTriangle(wnode); } /* * Patern 7 */ void F3_Patern7(cnode, cyc, pnode) int cnode[8][8]; int pnode[3]; { SetCornerOn(cnode,cyc,pnode[0]); SetCornerOn(cnode,cyc,pnode[1]); SetCornerOn(cnode,cyc,pnode[2]); } /* * For Case flag=3 * Patern 5,6,7 */ void Flag3(cnode) int cnode[8][8]; { int i, j, k, cyc; int pnode[3]; k = 1; for (i = 0; i < 8; i++){ if (ids[cnode[i][0]].flag == 1){ cyc = i; pnode[0] = 0; for(j = 1; j < 8; j++){ if (ids[cnode[i][j]].flag == 1){ pnode[k++] = j; } } break; } } /* printf("Flag3:cyc=%d,0-%d-%d \n",cyc,pnode[1],pnode[2]); */ switch(pnode[1]) { case 1: switch(pnode[2]) { case 2: /* Patern 5:0-1-2 */ F3_Patern5(cnode,cyc,0,4,2,6,1,5); F3_Patern5(cnode,cyc,0,4,0,3,2,6); F3_Patern5(cnode,cyc,2,6,0,3,2,3); break; case 3: /* Patern 5:0-1-3 */ F3_Patern5(cnode,cyc,1,5,0,4,3,7); F3_Patern5(cnode,cyc,1,5,3,7,1,2); F3_Patern5(cnode,cyc,1,2,3,7,2,3); break; case 4: /* Patern 5:0-1-4 */ F3_Patern5(cnode,cyc,4,7,0,3,1,2); F3_Patern5(cnode,cyc,4,7,1,2,4,5); F3_Patern5(cnode,cyc,4,5,1,2,1,5); break; case 5: /* Patern 5:0-1-5 */ F3_Patern5(cnode,cyc,0,3,1,2,5,6); F3_Patern5(cnode,cyc,0,3,5,6,0,4); F3_Patern5(cnode,cyc,0,4,5,6,4,5); break; case 6: /* Patern 6:0-1-6 */ case 7: /* Patern 6:0-1-7 */ F2_Patern2_1(cnode,cyc,pnode[0],pnode[1]); SetCornerOn(cnode,cyc,pnode[2]); break; } break; case 2: switch(pnode[2]) { case 3: /* Patern 5:0-2-3 */ F3_Patern5(cnode,cyc,0,4,3,7,2,6); F3_Patern5(cnode,cyc,0,4,2,6,0,1); F3_Patern5(cnode,cyc,0,1,2,6,1,2); break; case 4: /* Patern 6:0-2-4 */ F2_Patern2_3(cnode,cyc,pnode[0],pnode[2]); SetCornerOn(cnode,cyc,pnode[1]); break; case 5: /* Patern 7:0-2-5 */ F3_Patern7(cnode,cyc,pnode); break; case 6: /* Patern 6:0-2-6 */ F2_Patern2_3(cnode,cyc,pnode[1],pnode[2]); SetCornerOn(cnode,cyc,pnode[0]); break; case 7: /* Patern 7:0-2-7 */ F3_Patern7(cnode,cyc,pnode); break; } break; case 3: switch(pnode[2]) { case 4: /* Patern 5:0-3-4 */ F3_Patern5(cnode,cyc,2,3,0,1,4,5); F3_Patern5(cnode,cyc,2,3,4,5,3,7); F3_Patern5(cnode,cyc,3,7,4,5,4,7); break; case 5: /* Patern 6:0-3-5 */ case 6: /* Patern 6:0-3-6 */ F2_Patern2_1(cnode,cyc,pnode[1],pnode[0]); SetCornerOn(cnode,cyc,pnode[2]); break; case 7: /* Patern 5:0-3-7 */ F3_Patern5(cnode,cyc,6,7,2,3,0,1); F3_Patern5(cnode,cyc,6,7,0,1,4,7); F3_Patern5(cnode,cyc,4,7,0,1,0,4); break; } break; case 4: switch(pnode[2]) { case 5: /* Patern 5:0-4-5 */ F3_Patern5(cnode,cyc,5,6,4,7,0,3); F3_Patern5(cnode,cyc,5,6,0,3,1,5); F3_Patern5(cnode,cyc,1,5,0,3,0,1); break; case 6: /* Patern 6:0-4-6 */ F2_Patern2_3(cnode,cyc,pnode[0],pnode[1]); SetCornerOn(cnode,cyc,pnode[2]); break; case 7: /* Patern 5:0-4-7 */ F3_Patern5(cnode,cyc,0,1,4,5,6,7); F3_Patern5(cnode,cyc,0,1,6,7,0,3); F3_Patern5(cnode,cyc,0,3,6,7,3,7); break; } break; case 5: switch(pnode[2]) { case 6: /* Patern 6:0-5-6 */ F2_Patern2_2(cnode,cyc,pnode[1],pnode[2]); SetCornerOn(cnode,cyc,pnode[0]); break; case 7: /* Patern 7:0-5-7 */ F3_Patern7(cnode,cyc,pnode); break; } break; case 6: switch(pnode[2]) { case 7: /* Patern 6:0-6-7 */ F2_Patern2_2(cnode,cyc,pnode[1],pnode[2]); SetCornerOn(cnode,cyc,pnode[0]); break; } } /* end of switch(pnode[1]) */ } /* * Patern5 for flag=5 */ void F5_Patern5(cnode, cyc, p1, p2, p3, p4, p5, p6) int cnode[8][8]; int cyc,p1,p2,p3,p4,p5,p6; { int wnode[6]; wnode[0] = cnode[cyc][p1]; wnode[1] = cnode[cyc][p2]; wnode[4] = cnode[cyc][p3]; wnode[5] = cnode[cyc][p4]; wnode[2] = cnode[cyc][p5]; wnode[3] = cnode[cyc][p6]; DisplayTriangle(wnode); } /* * Patern 7 for flag=5 */ void F5_Patern7(cnode, cyc, pnode) int cnode[8][8]; int pnode[3]; { SetCornerOff(cnode,cyc,pnode[0]); SetCornerOff(cnode,cyc,pnode[1]); SetCornerOff(cnode,cyc,pnode[2]); } /* * For Case flag=5 * flag=3 の反転パターン * Patern 5,6,7 */ void Flag5(cnode) int cnode[8][8]; { int i, j, k, cyc; int pnode[3]; k = 1; for (i = 0; i < 8; i++){ if (ids[cnode[i][0]].flag == 0){ cyc = i; pnode[0] = 0; for(j = 1; j < 8; j++){ if (ids[cnode[i][j]].flag == 0){ pnode[k++] = j; } } break; } } /* printf("Flag5:cyc=%d,0-%d-%d \n",cyc,pnode[1],pnode[2]); */ switch(pnode[1]) { case 1: switch(pnode[2]) { case 2: /* Patern 5:0-1-2 */ F5_Patern5(cnode,cyc,0,4,2,6,1,5); F5_Patern5(cnode,cyc,0,4,0,3,2,6); F5_Patern5(cnode,cyc,2,6,0,3,2,3); break; case 3: /* Patern 5:0-1-3 */ F5_Patern5(cnode,cyc,1,5,0,4,3,7); F5_Patern5(cnode,cyc,1,5,3,7,1,2); F5_Patern5(cnode,cyc,1,2,3,7,2,3); break; case 4: /* Patern 5:0-1-4 */ F5_Patern5(cnode,cyc,4,7,0,3,1,2); F5_Patern5(cnode,cyc,4,7,1,2,4,5); F5_Patern5(cnode,cyc,4,5,1,2,1,5); break; case 5: /* Patern 5:0-1-5 */ F5_Patern5(cnode,cyc,0,3,1,2,5,6); F5_Patern5(cnode,cyc,0,3,5,6,0,4); F5_Patern5(cnode,cyc,0,4,5,6,4,5); break; case 6: /* Patern 6:0-1-6 */ case 7: /* Patern 6:0-1-7 */ F6_Patern2_1(cnode,cyc,pnode[0],pnode[1]); SetCornerOff(cnode,cyc,pnode[2]); break; } break; case 2: switch(pnode[2]) { case 3: /* Patern 5:0-2-3 */ F5_Patern5(cnode,cyc,0,4,3,7,2,6); F5_Patern5(cnode,cyc,0,4,2,6,0,1); F5_Patern5(cnode,cyc,0,1,2,6,1,2); break; case 4: /* Patern 6:0-2-4 */ F6_Patern2_3(cnode,cyc,pnode[0],pnode[2]); SetCornerOff(cnode,cyc,pnode[1]); break; case 5: /* Patern 7:0-2-5 */ F5_Patern7(cnode,cyc,pnode); break; case 6: /* Patern 6:0-2-6 */ F6_Patern2_3(cnode,cyc,pnode[1],pnode[2]); SetCornerOff(cnode,cyc,pnode[0]); break; case 7: /* Patern 7:0-2-7 */ F5_Patern7(cnode,cyc,pnode); break; } break; case 3: switch(pnode[2]) { case 4: /* Patern 5:0-3-4 */ F5_Patern5(cnode,cyc,2,3,0,1,4,5); F5_Patern5(cnode,cyc,2,3,4,5,3,7); F5_Patern5(cnode,cyc,3,7,4,5,4,7); break; case 5: /* Patern 6:0-3-5 */ case 6: /* Patern 6:0-3-6 */ F6_Patern2_1(cnode,cyc,pnode[1],pnode[0]); SetCornerOff(cnode,cyc,pnode[2]); break; case 7: /* Patern 5:0-3-7 */ F5_Patern5(cnode,cyc,6,7,2,3,0,1); F5_Patern5(cnode,cyc,6,7,0,1,4,7); F5_Patern5(cnode,cyc,4,7,0,1,0,4); break; } break; case 4: switch(pnode[2]) { case 5: /* Patern 5:0-4-5 */ F5_Patern5(cnode,cyc,5,6,4,7,0,3); F5_Patern5(cnode,cyc,5,6,0,3,1,5); F5_Patern5(cnode,cyc,1,5,0,3,0,1); break; case 6: /* Patern 6:0-4-6 */ F6_Patern2_3(cnode,cyc,pnode[0],pnode[1]); SetCornerOff(cnode,cyc,pnode[2]); break; case 7: /* Patern 5:0-4-7 */ F5_Patern5(cnode,cyc,0,1,4,5,6,7); F5_Patern5(cnode,cyc,0,1,6,7,0,3); F5_Patern5(cnode,cyc,0,3,6,7,3,7); break; } break; case 5: switch(pnode[2]) { case 6: /* Patern 6:0-5-6 */ F6_Patern2_2(cnode,cyc,pnode[1],pnode[2]); SetCornerOff(cnode,cyc,pnode[0]); break; case 7: /* Patern 7:0-5-7 */ F5_Patern7(cnode,cyc,pnode); break; } break; case 6: switch(pnode[2]) { case 7: /* Patern 6:0-6-7 */ F6_Patern2_2(cnode,cyc,pnode[1],pnode[2]); SetCornerOff(cnode,cyc,pnode[0]); break; } } /* end of switch(pnode[1]) */ } /* * Patern 9 * F4:0-1-3-4 のcaseを基準 */ void F4_Patern9(cnode, cyc, p1, p2, p3, p4, p5, p7) int cnode[8][8]; int cyc,p1,p2,p3,p4,p5,p7; { F3_Patern5(cnode,cyc,p1,p5,p4,p5,p4,p7); F3_Patern5(cnode,cyc,p1,p5,p4,p7,p1,p2); F3_Patern5(cnode,cyc,p1,p2,p4,p7,p3,p7); F3_Patern5(cnode,cyc,p1,p2,p3,p7,p2,p3); } /* * Patern 11 * F4:0-1-3-7 のcaseを基準 */ void F4_Patern11(cnode, cyc, p0, p1, p2, p3, p4, p5, p6, p7, p8) int cnode[8][8]; int cyc,p0,p1,p2,p3,p4,p5,p6,p7,p8; { F3_Patern5(cnode,cyc,p1,p5,p0,p4,p1,p2); F3_Patern5(cnode,cyc,p0,p4,p6,p7,p1,p2); F3_Patern5(cnode,cyc,p0,p4,p4,p7,p6,p7); F3_Patern5(cnode,cyc,p1,p2,p6,p7,p2,p3); } /* * Patern 13 */ void F4_Patern13(cnode, cyc, pnode) int cnode[8][8]; int pnode[4]; { SetCornerOn(cnode,cyc,pnode[0]); SetCornerOn(cnode,cyc,pnode[1]); SetCornerOn(cnode,cyc,pnode[2]); SetCornerOn(cnode,cyc,pnode[3]); } /* * Patern 14 * F4:0-2-3-4 のcaseを基準 */ void F4_Patern14(cnode, cyc, p0, p1, p2, p3, p4, p5, p6, p7, p8) int cnode[8][8]; int cyc,p0,p1,p2,p3,p4,p5,p6,p7,p8; { F3_Patern5(cnode,cyc,p0,p1,p4,p5,p1,p2); F3_Patern5(cnode,cyc,p4,p5,p3,p7,p1,p2); F3_Patern5(cnode,cyc,p1,p2,p3,p7,p2,p6); F3_Patern5(cnode,cyc,p4,p5,p4,p7,p3,p7); } /* * For Case flag=4 * Patern 8,9,10,11,12,13,14 */ void Flag4(cnode) int cnode[8][8]; { int i, j, k, cyc; int pnode[4]; k = 1; for (i = 0; i < 8; i++){ if (ids[cnode[i][0]].flag == 1){ cyc = i; pnode[0] = 0; for(j = 1; j < 8; j++){ if (ids[cnode[i][j]].flag == 1){ pnode[k++] = j; } } break; } } /* printf("Flag4:cyc=%d,0-%d-%d-%d \n",cyc,pnode[1],pnode[2],pnode[3]); */ switch(pnode[1]) { case 1: switch(pnode[2]) { case 2: switch(pnode[3]) { case 3: /* patern 8:0-1-2-3 */ F3_Patern5(cnode,cyc,0,4,3,7,1,5); F3_Patern5(cnode,cyc,1,5,3,7,2,6); break; case 4: /* Patern 11:0-1-2-4 */ F4_Patern11(cnode,cyc,1,2,3,0,5,6,7,4); break; case 5: /* patern 9:0-1-2-5 */ F4_Patern9(cnode,cyc,2,3,0,5,6,4); break; case 6: /* Patern 14:0-1-2-6 */ F4_Patern14(cnode,cyc,2,3,0,1,6,7,4,5); break; case 7: /* Patern 12:0-1-2-7 */ F3_Patern5(cnode,cyc,0,4,2,6,1,5); F3_Patern5(cnode,cyc,0,4,0,3,2,6); F3_Patern5(cnode,cyc,2,6,0,3,2,3); SetCornerOn(cnode,cyc,7); break; } break; case 3: switch(pnode[3]) { case 4: /* patern 9:0-1-3-4 */ F4_Patern9(cnode,cyc,1,2,3,4,5,7); break; case 5: /* Patern 14:0-1-3-5 */ F4_Patern14(cnode,cyc,1,2,3,0,5,6,7,4); break; case 6: /* Patern 12:0-1-3-6 */ F3_Patern5(cnode,cyc,1,5,0,4,3,7); F3_Patern5(cnode,cyc,1,5,3,7,1,2); F3_Patern5(cnode,cyc,1,2,3,7,2,3); SetCornerOn(cnode,cyc,6); break; case 7: /* Patern 11:0-1-3-7 */ F4_Patern11(cnode,cyc,0,1,2,3,4,5,6,7); break; } break; case 4: switch(pnode[3]) { case 5: /* patern 8 0-1-4-5 */ F3_Patern5(cnode,cyc,4,7,0,3,5,6); F3_Patern5(cnode,cyc,5,6,0,3,1,2); break; case 6: /* Patern 12:0-1-4-6 */ F3_Patern5(cnode,cyc,4,7,0,3,1,2); F3_Patern5(cnode,cyc,4,7,1,2,4,5); F3_Patern5(cnode,cyc,4,5,1,2,1,5); SetCornerOn(cnode,cyc,6); break; case 7: /* Patern 14:0-1-4-7 */ F4_Patern14(cnode,cyc,4,5,1,0,7,6,2,3); break; } break; case 5: switch(pnode[3]) { case 6: /* Patern 11:0-1-5-6 */ F4_Patern11(cnode,cyc,1,0,4,5,2,3,7,6); break; case 7: /* Patern 12:0-1-5-7 */ F3_Patern5(cnode,cyc,0,3,1,2,5,6); F3_Patern5(cnode,cyc,0,3,5,6,0,4); F3_Patern5(cnode,cyc,0,4,5,6,4,5); SetCornerOn(cnode,cyc,7); break; } break; case 6: switch(pnode[3]) { case 7: /* patern 10:0-1-6-7 */ F2_Patern2_1(cnode,cyc,0,1); F2_Patern2_2(cnode,cyc,6,7); break; } break; case 7: break; }/* end of switch(pnode[2]) */ break; case 2: switch(pnode[2]) { case 3: switch(pnode[3]) { case 4: /* Patern 14:0-2-3-4 */ F4_Patern14(cnode,cyc,0,1,2,3,4,5,6,7,8); break; case 5: /* Patern 12:0-2-3-5 */ F3_Patern5(cnode,cyc,0,4,3,7,2,6); F3_Patern5(cnode,cyc,0,4,2,6,0,1); F3_Patern5(cnode,cyc,0,1,2,6,1,2); SetCornerOn(cnode,cyc,5); break; case 6: /* Patern 11:0-2-3-6 */ F4_Patern11(cnode,cyc,3,0,1,2,7,4,5,6); break; case 7: /* patern 9:0-2-3-7 */ F4_Patern9(cnode,cyc,0,1,2,7,4,6); break; } break; case 4: switch(pnode[3]) { case 5: /* Patern 12:0-2-4-5 */ F3_Patern5(cnode,cyc,5,6,4,7,0,3); F3_Patern5(cnode,cyc,5,6,0,3,1,5); F3_Patern5(cnode,cyc,1,5,0,3,0,1); SetCornerOn(cnode,cyc,2); break; case 6: /* patern 10:0-2-4-6 */ F2_Patern2_3(cnode,cyc,0,4); F2_Patern2_3(cnode,cyc,2,6); break; case 7: /* Patern 12:0-2-4-7 */ F3_Patern5(cnode,cyc,0,1,4,5,6,7); F3_Patern5(cnode,cyc,0,1,6,7,0,3); F3_Patern5(cnode,cyc,0,3,6,7,3,7); SetCornerOn(cnode,cyc,2); break; } break; case 5: switch(pnode[3]) { case 7: /* patern 13:0-2-5-7 */ F4_Patern13(cnode,cyc,pnode); break; } break; }/* end of switch(pnode[2]) */ break; case 3: switch(pnode[2]) { case 4: switch(pnode[3]) { case 5: /* Patern 11:0-3-4-5 */ F4_Patern11(cnode,cyc,0,3,7,4,1,2,6,5); break; case 6: /* Patern 12:0-3-4-6 */ F3_Patern5(cnode,cyc,2,3,0,1,4,5); F3_Patern5(cnode,cyc,2,3,4,5,3,7); F3_Patern5(cnode,cyc,3,7,4,5,4,7); SetCornerOn(cnode,cyc,6); break; case 7: /* patern 8:0-3-4-7 */ F3_Patern5(cnode,cyc,4,5,6,7,0,1); F3_Patern5(cnode,cyc,0,1,6,7,2,3); break; } break; case 5: switch(pnode[3]) { case 6: /* patern 10:0-3-5-6 */ F2_Patern2_1(cnode,cyc,3,0); F2_Patern2_2(cnode,cyc,5,6); break; case 7: /* Patern 12:0-3-5-7 */ F3_Patern5(cnode,cyc,6,7,2,3,0,1); F3_Patern5(cnode,cyc,6,7,0,1,4,7); F3_Patern5(cnode,cyc,4,7,0,1,0,4); SetCornerOn(cnode,cyc,5); break; } break; case 6: switch(pnode[3]) { case 7: /* Patern 14:0-3-6-7 */ F4_Patern14(cnode,cyc,3,2,6,7,0,1,5,4); break; } break; }/* end of switch(pnode[2]) */ break; case 4: switch(pnode[2]) { case 5: switch(pnode[3]) { case 6: /* Patern 14:0-4-5-6 */ F4_Patern14(cnode,cyc,4,7,6,5,0,3,2,1); break; case 7: /* patern 9:0-4-5-7 */ F4_Patern9(cnode,cyc,7,6,5,0,3,1); break; } break; case 6: switch(pnode[3]) { case 7: /* Patern 11:0-4-6-7 */ F4_Patern11(cnode,cyc,7,6,5,4,3,2,1,0); break; } break; }/* end of switch(pnode[2]) */ break; case 5: switch(pnode[2]) { case 6: break; }/* end of switch(pnode[2]) */ break; } /* end of switch(pnode[1]) */ } void SelectCell(stpoint) int stpoint; { int node[8], cnode[8][8]; int idz,i,j; int flag; node[0] = stpoint; node[1] = node[0]+1; node[2] = node[0]+idmx+1; node[3] = node[0]+idmx; idz = idmx*idmy; node[4] = node[0]+idz; node[5] = node[1]+idz; node[6] = node[2]+idz; node[7] = node[3]+idz; flag = 0; for (i=0;i<8;i++){ flag += ids[node[i]].flag; } if (flag == 0) return; /* // nodeのindexをcyclicにcnodeにセット */ for (i=0;i<4;i++){ for (j=0;j<4;j++){ cnode[i][j] = node[(i+j)%4]; } for (j=4;j<8;j++){ cnode[i][j] = node[4+(i+j)%4]; } } for (i=4;i<8;i++){ for (j=0;j<4;j++){ cnode[i][j] = node[4+(i+8-j)%4]; } for (j=4;j<8;j++){ cnode[i][j] = node[(i+8-j)%4]; } } switch(flag) { case 1: Flag1(cnode); break; case 2: Flag2(cnode); break; case 3: Flag3(cnode); break; case 4: Flag4(cnode); break; case 5: Flag5(cnode); break; case 6: Flag6(cnode); break; case 7: Flag7(cnode); break; } } MakeOtherPolygon(v, f, c, fn) float v; int f, c; char *fn; { int i, j, k; FILE *fo; long pl[PMAX]; float pf[PMAX], v1,v2,vs,vt1,vt2; viso = v; printf("-----------------------\n"); printf("Make polygon data No.%d\n", c+1); printf("isovalue = %+f\n", v); FlagSelect(f + (c - 1) * 2); printf("Z loop "); for (k = 0; k < idmz - 1; k++){ printf("[%d]", k + 1); fflush(stdout); for (j = 0; j < idmy - 1; j++){ for (i = 0; i < idmx - 1; i++){ SelectCell(i + idmx * (j + idmy * k)); } } } printf("\nNo.%d polygon total = %d\n", c+1, pgdct); printf("Make normal vector No.%d\n",c+1); NormalizeVector(); printf("Write polygon data No.%d\n",c+1); if ((fo = fopen(fn, "rb+")) != NULL){ fread(pl, sizeof(long), PMAX, fo); fread(pf, sizeof(float), PMAX, fo); if (cmode) swapper((char *)pl, PMAX); if (cmode) swapper((char *)pf, PMAX); pl[5]++; /* object count */ pl[20+c] = pgdct; /* polygon count 1 */ pf[20+c] = viso; /* iso value 1 */ fseek(fo, 0L, SEEK_SET); if (cmode) swapper((char *)pl, PMAX); if (cmode) swapper((char *)pf, PMAX); fwrite(pl, sizeof(long), PMAX, fo); fwrite(pf, sizeof(float), PMAX, fo); fseek(fo, 0L, SEEK_END); if (cmode) swapper((char *)pgd, pgdct * (sizeof(PGD)/sizeof(long))); fwrite(pgd, sizeof(PGD), pgdct, fo); fclose(fo); } } /* // input argument // *argv[0] : program name // *argv[1] : input file name // *argv[2] : iso value */ void main(argc, argv) int argc; char **argv; { int winx=800, winy=800; int i, j, k; int db, ac; char filein[1024], fileout[1024]; FILE *fi, *fo; float dt[18]; long pl[PMAX]; float pf[PMAX], v1,v2,vs,vt1,vt2; float mbox; printf("\n"); printf("D2P Copyright(C) 1996-2001 RCCS\n\n"); db = 0; if (argc != 3 && argc != 4) { printf("USAGE: D2P infile iso_value\n"); exit(1); } strcpy(filein, argv[1]); sscanf(argv[2],"%f",&vvvv); viso = vvvv; ac = 3; while (1){ if (ac >= argc) break; if (argv[ac][0] == '-' || argv[ac][0] == '/'){ if (toupper(argv[ac][1]) == 'D') dble = 2; if (toupper(argv[ac][1]) == 'T') dble = 3; } ac++; } strcpy(fileout, filein); for (i = strlen(fileout) - 1; i >= 0; i--){ if (fileout[i] == '.') break; } if (i > 0){ fileout[i] = 0x00; if (dble == 1) strcat(fileout, ".pgv"); if (dble == 2) strcat(fileout, "d.pgv"); if (dble == 3) strcat(fileout, "t.pgv"); } printf("argc= %d infile = %s outfile = %s iso_value=%f\n", argc, filein, fileout, viso); if((fi = fopen(filein, "rb"))==NULL){ perror("fopen error "); exit(1); } InputData(fi); fclose(fi); printf("-----------------------\n"); printf("Make polygon data\n"); FlagSelect(0); printf("Z loop "); for (k = 0; k < idmz - 1; k++){ printf("[%d]", k + 1); fflush(stdout); for (j = 0; j < idmy - 1; j++){ for (i = 0; i < idmx - 1; i++){ SelectCell(i + idmx * (j + idmy * k)); } } } printf("\nPolygon total = %d\n", pgdct); printf("Make normal vector\n"); NormalizeVector(); pl[0] = (long)GCODE; /* program code */ pl[1] = dtype; /* data type */ pl[2] = idmx; /* x resolution */ pl[3] = idmy; /* y resolution */ pl[4] = idmz; /* z resolution */ pl[5] = 1; /* object count */ pl[6] = oct; /* orbital count */ pl[7] = molct; /* moleculer count */ pl[14] = obt[0]; pl[15] = obt[1]; pl[16] = obt[2]; pl[17] = obt[3]; pl[18] = obt[4]; pl[19] = obt[5]; pl[20] = pgdct; /* polygon count 1 */ pf[0] = mbgx; pf[1] = mbgy; pf[2] = mbgz; pf[3] = mstx; pf[4] = msty; pf[5] = mstz; pf[6] = xmax; pf[7] = xmin; pf[8] = ymax; pf[9] = ymin; pf[10] = zmax; pf[11] = zmin; pf[12] = datamax; pf[13] = datamin; pf[20] = viso; /* iso value 1 */ if (DCODE == DOSCOD) cmode = FALSE; else cmode = TRUE; if (db == 0 && dtype == 1) db = 1; if (db == 1) printf("Double data mode\n"); if (cmode) printf("Swap byte mode\n"); else printf("Not swap byte mode\n"); printf("Write polygon data\n"); if ((fo = fopen(fileout, "wb")) != NULL){ if (cmode) swapper((char *)pl, PMAX); if (cmode) swapper((char *)pf, PMAX); if (cmode) swapper((char *)mol, molct * (sizeof(MOL)/sizeof(long))); if (cmode) swapper((char *)pgd, pgdct * (sizeof(PGD)/sizeof(long))); fwrite(pl, sizeof(long), PMAX, fo); fwrite(pf, sizeof(float), PMAX, fo); fwrite(mol, sizeof(MOL), molct, fo); fwrite(pgd, sizeof(PGD), pgdct, fo); fclose(fo); } if (dble != 1){ v1 = log10(datamax); v2 = log10(viso); vs = (v1 - v2) / dble; } if (db == 1){ MakeOtherPolygon(-vvvv, 1, 1, fileout); if (dble == 2){ vt1 = (float)pow((double)10.0, (double)(v1 - vs)); MakeOtherPolygon( vt1, 0, 2, fileout); MakeOtherPolygon(-vt1, 1, 3, fileout); } if (dble == 3){ vt1 = (float)pow((double)10.0, (double)(v1 - vs * 2)); vt2 = (float)pow((double)10.0, (double)(v1 - vs)); MakeOtherPolygon( vt1, 0, 2, fileout); MakeOtherPolygon(-vt1, 1, 3, fileout); MakeOtherPolygon( vt2, 0, 4, fileout); MakeOtherPolygon(-vt2, 1, 5, fileout); } } else { if (dble == 2){ vt1 = (float)pow((double)10.0, (double)(v1 - vs)); MakeOtherPolygon(vt1, 0, 2, fileout); } if (dble == 3){ vt1 = (float)pow((double)10.0, (double)(v1 - vs * 2)); vt2 = (float)pow((double)10.0, (double)(v1 - vs)); MakeOtherPolygon(vt1, 0, 2, fileout); MakeOtherPolygon(vt2, 0, 3, fileout); } } printf("-----------------------\n"); printf("Analysis for volume\n"); mbox = mstx * msty * mstz; if (db == 1){ printf("Inside boxel count = %d ( %d : %d )\n", incount[0] + incount[1], incount[0], incount[1]); printf("Volume of boxel = %f\n", mbox); printf("Volume of total = %f ( %f : %f)\n", (incount[0] + incount[1]) * mbox, incount[0] * mbox, incount[1] * mbox); } else { printf("Inside boxel count = %d\n", incount[0]); printf("Volume of boxel = %f\n", mbox); printf("Volume of total = %f\n", incount[0] * mbox); } if (opencount > 0) printf(" There are open areas.\n"); printf("-----------------------\n"); printf("Terminate\n"); }