#include #include #include #include "isofile.h" #include "mymalloc.h" enum boolean {FALSE, TRUE}; /* calcule le gradient du champ level au point x, y, z second ordre aux limites, 4eme ordre a l'interieur */ static void gradient(float ***level, int x, int y, int z, int lx, int ly, int lz, double *nx, double *ny, double *nz) { if (x == 0 || x == 1) *nx = 0.5 * (-3.*level[x][y][z]+4.*level[x+1][y][z]-level[x+2][y][z]); else if (x == lx - 1 || x == lx - 2) *nx = 0.5 * (3.*level[x][y][z]-4.*level[x-1][y][z]+level[x-2][y][z]); else *nx = (-level[x+2][y][z]+8.*(level[x+1][y][z]-level[x-1][y][z]) +level[x-2][y][z])/12.; if (y == 0 || y == 1) *ny = 0.5 * (-3.*level[x][y][z]+4.*level[x][y+1][z]-level[x][y+2][z]); else if (y == ly - 1 || y == ly - 2) *ny = 0.5 * (3.*level[x][y][z]-4.*level[x][y-1][z]+level[x][y-2][z]); else *ny = (-level[x][y+2][z]+8.*(level[x][y+1][z]-level[x][y-1][z]) +level[x][y-2][z])/12.; if (z == 0 || z == 1) *nz = 0.5 * (-3.*level[x][y][z]+4.*level[x][y][z+1]-level[x][y][z+2]); else if (z == lz - 1 || z == lz - 2) *nz = 0.5 * (3.*level[x][y][z]-4.*level[x][y][z-1]+level[x][y][z-2]); else *nz = (-level[x][y][z+2]+8.*(level[x][y][z+1]-level[x][y][z-1]) +level[x][y][z-2])/12.; } #include "mirror.h" int isomeshfile(float ***level, int lx, int ly, int lz, double cut_level, FILE *fptr, int order) { static char xs[12] = {0, 0, 0, 0, 0, 0, 1, 1, 0, 1, 1, 0}; static char ys[12] = {0, 0, 1, 1, 1, 0, 0, 1, 0, 0, 1, 1}; static char zs[12] = {0, 1, 1, 0, 0, 1, 1, 0, 0, 0, 1, 1}; static char xe[12] = {1, 1, 1, 1, 0, 0, 1, 1, 0, 1, 1, 0}; static char ye[12] = {0, 0, 1, 1, 0, 1, 1, 0, 0, 0, 1, 1}; static char ze[12] = {0, 1, 1, 0, 0, 1, 1, 0, 1, 1, 0, 0}; double xpt[12], ypt[12], zpt[12], nx[12], ny[12], nz[12]; static char fi[12][2] = { {0, 3}, {0, 1}, {1, 2}, {2, 3}, {3, 4}, {1, 4}, {1, 5}, {3, 5}, {0, 4}, {0, 5}, {2, 5}, {2, 4} }; static int face[6][4] = { {0, 9, 1, 8}, {1, 6, 2, 5}, {2, 10, 3, 11}, {3, 7, 0, 4}, {5, 11, 4, 8}, {7, 10, 6, 9} }; static char fx[6][5] = { {0, 1, 1, 0, 0}, {0, 1, 1, 0, 0}, {0, 1, 1, 0, 0}, {0, 1, 1, 0, 0}, {0, 0, 0, 0, 0}, {1, 1, 1, 1, 1} }; static char fy[6][5] = { {0, 0, 0, 0, 0}, {0, 0, 1, 1, 0}, {1, 1, 1, 1, 1}, {1, 1, 0, 0, 1}, {0, 1, 1, 0, 0}, {0, 1, 1, 0, 0} }; static char fz[6][5] = { {0, 0, 1, 1, 0}, {1, 1, 1, 1, 1}, {1, 1, 0, 0, 1}, {0, 0, 0, 0, 0}, {1, 1, 0, 0, 1}, {0, 0, 1, 1, 0} }; int seg[6][12]; int seg1[12]; int cut, count, ntriangle = 0; double levels, levele, a; int x, y, z, i, j, k, nb, f, ar, max; double poly[12][3], norm[12][3]; double nxe, nye, nze, nxs, nys, nzs; double scale; /* double xc = (double)(lx - 1)/2., yc = (double)(ly - 1)/2., zc = (double)(lz - 1)/2.; */ double xc = 0., yc=0., zc=0.; int kv=0; int meshtype=2; int maxtriangles = 4000000; int maxvertex = 4000000; double **vertex_table, **norm_table; int **vertex_index; fprintf(stderr, "Generating triangles, 1 dot per 10000 vertices\n"); if (!( vertex_table = (double **)mymalloc(sizeof(double), maxvertex,3))) { fprintf(stderr, "ZM3D: can not allocate array\n"); exit(1); } if (!( norm_table = (double **)mymalloc(sizeof(double), maxvertex,3))) { fprintf(stderr, "ZM3D: can not allocate array\n"); exit(1); } if (!( vertex_index = (int **)mymalloc(sizeof(int), maxtriangles,3))) { fprintf(stderr, "ZM3D: can not allocate array\n"); exit(1); } max = lx; if (ly > max) max = ly; if (lz > max) max = lz; scale = 1./(double)(max - 1); if(meshtype == 1) fprintf(fptr, "#declare Goutte = mesh { \n"); /* calcul pour chaque point */ /* x, y, z sont des entiers de 0 a lx -1 */ for (x = 0; x < lx - 1; x++) { for (y = 0; y < ly - 1; y++) for (z = 0; z < lz - 1; z++) { /* calcul des points d'intersection avec chaque arrete */ cut = FALSE; /* xs,xe,ys, ... sont des coordonnées entières de sommets relatives au cube considéré */ for (i = 0; i < 12; i++) { levels = level[x + xs[i]][y + ys[i]][z + zs[i]]; levele = level[x + xe[i]][y + ye[i]][z + ze[i]]; if (((levels >= cut_level) && (levele < cut_level)) || ((levels < cut_level) && (levele >= cut_level))) { a = (cut_level - levels) / (levele - levels); /* xpt sont des coordonnées flottantes relatives relatives au cube considéré */ xpt[i] = a * ((double)xe[i] - (double)xs[i]) + (double)xs[i]; ypt[i] = a * ((double)ye[i] - (double)ys[i]) + (double)ys[i]; zpt[i] = a * ((double)ze[i] - (double)zs[i]) + (double)zs[i]; gradient(level, x + xs[i], y + ys[i], z + zs[i], lx, ly, lz, &nxs, &nys, &nzs); gradient(level, x + xe[i], y + ye[i], z + ze[i], lx, ly, lz, &nxe, &nye, &nze); nx[i] = a * (nxe - nxs) + nxs; ny[i] = a * (nye - nys) + nys; nz[i] = a * (nze - nzs) + nzs; a = sqrt(nx[i]*nx[i] + ny[i]*ny[i] + nz[i]*nz[i]); #ifdef DEBUG if(a==0) { printf(" x y z i %d %d %d %d \n",x,y,z,i); printf(" nx ny nz %g %g %g \n",nx[i],ny[i],nz[i]); printf(" nxe nye nze %g %g %g \n",nxe,nye,nze); printf(" nxs nys nzs %g %g %g \n",nxs,nys,nzs); printf(" levels levele %g %g \n",levels,levele); nrerror("zero normal"); } #endif cut = TRUE; } else xpt[i] = -10.0; } /* aucune intersection */ if (cut == FALSE) continue; /* relie les points pour chaque face */ for (i = 0; i < 6; i++) { for (j = 0; j < 12; j++) seg[i][j] = -1; k = -1; count = 0; do { for (j = 0; j < 4; j++) { if (xpt[face[i][j]] == -10.0) continue; if (k == -1) k = face[i][j]; else { if (level[x+fx[i][j]][y+fy[i][j]][z+fz[i][j]] >= cut_level) { seg[i][k] = face[i][j]; k = -1; } else k = face[i][j]; } if (count == 1 && k == -1) break; } count++; } while(k != -1); } /* cree les facettes et les ecrit dans le fichier */ for (i = 0; i < 6; i++) for (j = 0; j < 12; j++) { f = i; ar = j; seg1[0] = j; nb = 1; while (seg[f][ar] != -1) { seg1[nb] = seg[f][ar]; seg[f][ar] = -1; ar = seg1[nb++]; f = (fi[ar][0] == f) ? fi[ar][1] : fi[ar][0]; } /* nb = nombre de cotés du polygone + 1 */ /* ecriture */ if (nb > 3) { for (k = 0; k < nb; k++) { /* poly sont des coordonnées flottantes correspondant au coordonnées entières x,y,z */ poly[k][0] = ((double)x + xpt[seg1[k]] - xc)*scale; poly[k][1] = ((double)y + ypt[seg1[k]] - yc)*scale; poly[k][2] = ((double)z + zpt[seg1[k]] - zc)*scale; a = sqrt(nx[seg1[k]]*nx[seg1[k]] + ny[seg1[k]]*ny[seg1[k]] + nz[seg1[k]]*nz[seg1[k]]); if(a==0) nrerror("zero normal"); norm[k][0] = nx[seg1[k]]/a; norm[k][1] = ny[seg1[k]]/a; norm[k][2] = nz[seg1[k]]/a; } /* cree des triangles a partir du polygone */ /* ecriture du tableau complet des sommets et des normes */ if((kv + nb - 1) > maxvertex) { printf(" x y z i %d %d %d %d \n",x,y,z,i); nrerror("max number of vertices exceeded"); } if((ntriangle + nb - 3) > maxtriangles) nrerror("max number of triangles exceeded"); for (k = 0; k < nb - 1 ; k++) for (i=0; i<3; i++) { vertex_table[kv + k][i] = poly[k][i]; norm_table[kv + k][i] = norm[k][i]; } /* ecriture du tableau complet des indices de face */ for (k = 2; k < nb - 1; k++) { vertex_index[ntriangle][0] = kv; vertex_index[ntriangle][1] = kv + k - 1; vertex_index[ntriangle][2] = kv + k; ntriangle++; } if(kv%10000 == 0) { fprintf(stderr,"."); fflush(stderr); } kv += nb-1; } else if (nb > 1) return -nb; } } /* We are at y=ly-1 and z=lz-1 thus at the corner we print the height of the local vertices to get an approximation of the hieght of the plane at infinity */ if(cut==TRUE) fprintf(stderr,"\n corner height %e %e %e \n",poly[0][0],poly[1][0],poly[2][0]); } if(ntriangle <= 0) return ntriangle; fprintf(fptr, "#declare Goutte = mesh2 { \n vertex_vectors { %d, \n",order*kv); printf("\nWriting vertices to file\n"); fflush(stderr); fprintf(stderr,"."); fflush(stderr); for(k=0 ; k < kv; k++) fprintf(fptr, " < %g, %g, %g > \n",vertex_table[k][0], vertex_table[k][1],vertex_table[k][2]); if(order> 1) { fprintf(stderr,"."); fflush(stderr); for(k=0 ; k < kv; k++) fprintf(fptr, " < %g, %g, %g > \n",MIRROR_P_X_0(vertex_table[k][0]), MIRROR_P_X_1(vertex_table[k][1]),MIRROR_P_X_2(vertex_table[k][2])); if(order > 2) { fprintf(stderr,"."); fflush(stderr); for(k=0 ; k < kv; k++) fprintf(fptr, " < %g, %g, %g > \n",MIRROR_P_Y_0(vertex_table[k][0]), MIRROR_P_Y_1(vertex_table[k][1]),MIRROR_P_Y_2(vertex_table[k][2])); fprintf(stderr,"."); fflush(stderr); for(k=0 ; k < kv; k++) fprintf(fptr, " < %g, %g, %g > \n",MIRROR_P_XY_0(vertex_table[k][0]), MIRROR_P_XY_1(vertex_table[k][1]),MIRROR_P_XY_2(vertex_table[k][2])); } } fprintf(fptr, "}\n"); fprintf(fptr, "normal_vectors { %d, \n",order*kv); printf("\nWriting normals to file\n"); for(k=0 ; k < kv; k++) fprintf(fptr, " < %g, %g, %g > \n",norm_table[k][0], norm_table[k][1],norm_table[k][2]); if(order > 1) { fprintf(stderr,"."); fflush(stderr); for(k=0 ; k < kv; k++) fprintf(fptr, " < %g, %g, %g > \n",MIRROR_V_X_0(norm_table[k][0]), MIRROR_V_X_1(norm_table[k][1]),MIRROR_V_X_2(norm_table[k][2])); if(order > 2) { fprintf(stderr,"."); fflush(stderr); for(k=0 ; k < kv; k++) fprintf(fptr, " < %g, %g, %g > \n",MIRROR_V_Y_0(norm_table[k][0]), MIRROR_V_Y_1(norm_table[k][1]),MIRROR_V_Y_2(norm_table[k][2])); fprintf(stderr,"."); fflush(stderr); for(k=0 ; k < kv; k++) fprintf(fptr, " < %g, %g, %g > \n",MIRROR_V_XY_0(norm_table[k][0]), MIRROR_V_XY_1(norm_table[k][1]),MIRROR_V_XY_2(norm_table[k][2])); } } fprintf(stderr,"."); fflush(stderr); fprintf(fptr, "}\n"); fprintf(fptr, "face_indices { %d, \n",order*ntriangle); printf("\nWriting face indices to file\n"); fflush(stderr); for(i=0;i \n",i*kv+ vertex_index[k][0], i*kv+vertex_index[k][1],i*kv+vertex_index[k][2]); } fprintf(stderr,"."); fflush(stderr); } fprintf(fptr, "} \n} \n"); /* myfree(vertex_table); myfree(norm_table); myfree(vertex_index); */ fprintf(stderr,"\nvertices %d \n",order*kv); fprintf(stderr,"triangles %d \n",order*ntriangle); return ntriangle; }