diff --git a/Common/Views.cpp b/Common/Views.cpp index d86998068e736684581362a54a61bcb7eaed2f08..ef5c50bbd229696b335a8dbad2604c7683e299c0 100644 --- a/Common/Views.cpp +++ b/Common/Views.cpp @@ -1,5 +1,6 @@ -// $Id: Views.cpp,v 1.23 2001-01-12 13:28:54 geuzaine Exp $ +// $Id: Views.cpp,v 1.24 2001-01-19 22:32:31 remacle Exp $ +#include <set> #include "Gmsh.h" #include "Views.h" #include "Context.h" @@ -265,6 +266,9 @@ void EndView(int AddInUI, int Number, char *FileName, char *Name, List_Replace(Post_ViewList,ActualView,fcmpPostViewNum); } + // it's a test, smoothing views in the volume !!! + // in the future, we'll have normals smoothed + ActualView->smooth(); ActualView = NULL; } @@ -605,3 +609,173 @@ void Read_View(FILE *file, char *filename){ } + +/* + A little util for smoothing a view. +*/ + +using namespace std; +struct xyzv +{ +private: +public: + double x,y,z,*vals; + int nbvals; + int nboccurences; + static double eps; + void update (int nbVals, double *); + xyzv(double x, double y, double z); + ~xyzv(); + xyzv & operator = ( const xyzv &); + xyzv ( const xyzv &); +}; + +double xyzv::eps = 0.0; + +xyzv::xyzv (double xx, double yy, double zz) : x(xx),y(yy),z(zz),vals(0),nbvals(0),nboccurences(0){} +xyzv::~xyzv() +{ + if(vals)delete [] vals; +} +xyzv::xyzv(const xyzv &other) +{ + x = other.x; + y = other.y; + z = other.z; + nbvals = other.nbvals; + nboccurences = other.nboccurences; + if(other.vals && other.nbvals) + { + vals = new double[other.nbvals]; + for(int i=0;i<nbvals;i++)vals[i] = other.vals[i]; + } +} +xyzv & xyzv::operator = (const xyzv &other) +{ + if(this != &other) + { + x = other.x; + y = other.y; + z = other.z; + nbvals = other.nbvals; + nboccurences = other.nboccurences; + if(other.vals && other.nbvals) + { + vals = new double[other.nbvals]; + for(int i=0;i<nbvals;i++)vals[i] = other.vals[i]; + } + } + return *this; +} + +void xyzv::update (int n, double *v) +{ + int i; + if(!vals) + { + vals = new double[n]; + for(i=0;i<nbvals;i++)vals[i] = 0.0; + nbvals = n; + nboccurences = 0; + } + else if (nbvals != n) + { + throw n; + } + double x1 = (double)(nboccurences)/ (double)(nboccurences + 1); + double x2 = 1./(double)(nboccurences + 1); + for(i=0;i<nbvals;i++)vals[i] = (x1 * vals[i] + x2 * v[i]); + nboccurences++; +} + + +struct lessthanxyzv +{ + bool operator () (const xyzv & p2, const xyzv &p1) const + { + if( p1.x - p2.x > xyzv::eps)return true; + if( p1.x - p2.x <-xyzv::eps)return false; + if( p1.y - p2.y > xyzv::eps)return true; + if( p1.y - p2.y <-xyzv::eps)return false; + if( p1.z - p2.z > xyzv::eps)return true; + return false; + } +}; + + +void Post_View :: smooth () +{ + int i,nb,j; + xyzv::eps = CTX.lc * 1.e-6; + + typedef set<xyzv,lessthanxyzv> mycont; + typedef mycont::const_iterator iter; + + mycont connectivities; + + double *x,*y,*z,*v; + if(NbSS){ + Msg(INFO,"Smoothing the view (%d) ...",NbTimeStep); + nb = List_Nbr(SS) / NbSS ; + + double *vals = new double[NbTimeStep]; + + for(i = 0 ; i < List_Nbr(SS) ; i+=nb) + { + x = (double*)List_Pointer_Fast(SS,i); + y = (double*)List_Pointer_Fast(SS,i+4); + z = (double*)List_Pointer_Fast(SS,i+8); + v = (double*)List_Pointer_Fast(SS,i+12); + + for(j=0;j<4;j++) + { + for(int k=0;k<NbTimeStep;k++)vals[k] = v[j+k*8]; + xyzv x(x[j],y[j],z[j]); + iter it = connectivities.find(x); + if(it == connectivities.end()) + { + x.update(NbTimeStep,vals); + connectivities.insert(x); + } + else + { + xyzv *xx = (xyzv*) &(*it); // a little weird ... becaus we know that + // this will not destroy the set ordering + xx->update(NbTimeStep,vals); + } + } + } + + int n1 = 0; + int n2 = 0; + int n3 = 0; + for(i = 0 ; i < List_Nbr(SS) ; i+=nb) + { + x = (double*)List_Pointer_Fast(SS,i); + y = (double*)List_Pointer_Fast(SS,i+4); + z = (double*)List_Pointer_Fast(SS,i+8); + v = (double*)List_Pointer_Fast(SS,i+12); + for(j=0;j<4;j++) + { + xyzv xyz(x[j],y[j],z[j]); + iter it = connectivities.find(xyz); + if(it == connectivities.end()) + { + n3++; + } + else + { + n1++; + n2 += (*it).nboccurences; + // test + //double val = sqrt (x[j] * x[j] + y[j] * y[j] + z[j] * z[j]); + //for(int k=0;k<NbTimeStep;k++)v[j+k*8] = val; + for(int k=0;k<NbTimeStep;k++)v[j+k*8] = (*it).vals[k]; + } + } + } + + Msg(INFO,"nbpoints = %d mean = %d size = %d miss %d\n",n1,n2/n1,connectivities.size(),n3); + delete [] vals; + } +} diff --git a/Common/Views.h b/Common/Views.h index 1f304f3b8ccdab22233b8af2c4c5dd783ea3e97a..afc260a847dee02c77615c91656f4fead72efa18 100644 --- a/Common/Views.h +++ b/Common/Views.h @@ -37,6 +37,8 @@ typedef struct{ // dynamic double (*GVFI) (double min, double max, int nb, int index); int (*GIFV) (double min, double max, int nb, double value); + // smooth the view + void smooth(); }Post_View; // The static list with pointers to all views diff --git a/Graphics/Entity.cpp b/Graphics/Entity.cpp index 731c4aa43b450c3dbbfff87aae8a7c181ad53768..7a3eb9faad2909e21a958e381f37f916d1fb5b8a 100644 --- a/Graphics/Entity.cpp +++ b/Graphics/Entity.cpp @@ -1,4 +1,4 @@ -// $Id: Entity.cpp,v 1.5 2001-01-10 08:50:30 geuzaine Exp $ +// $Id: Entity.cpp,v 1.6 2001-01-19 22:32:31 remacle Exp $ #include "Gmsh.h" #include "GmshUI.h" @@ -46,6 +46,7 @@ void Draw_Triangle (double *x, double *y, double *z, double x1x0, y1y0, z1z0, x2x0, y2y0, z2z0, nn[3]; + glBegin(GL_TRIANGLES); if (shade){ x1x0 = (x[1]+Raise[0][1]) - (x[0]+Raise[0][0]); y1y0 = (y[1]+Raise[1][1]) - (y[0]+Raise[1][0]); @@ -56,17 +57,9 @@ void Draw_Triangle (double *x, double *y, double *z, nn[0] = y1y0 * z2z0 - z1z0 * y2y0 ; nn[1] = z1z0 * x2x0 - x1x0 * z2z0 ; nn[2] = x1x0 * y2y0 - y1y0 * x2x0 ; - /* BOF BOF - if(nn[2] < -0.1){ - nn[0] = -nn[0]; - nn[1] = -nn[1]; - nn[2] = -nn[2]; - } - */ glNormal3dv(nn); } - glBegin(GL_TRIANGLES); glVertex3d(x[0]+Offset[0]+Raise[0][0], y[0]+Offset[1]+Raise[1][0], z[0]+Offset[2]+Raise[2][0]); @@ -87,46 +80,18 @@ void Draw_Triangle (double *x, double *y, double *z, void Draw_Quadrangle (double *x, double *y, double *z, double *Offset, double Raise[3][5], int shade){ - double x1x0, y1y0, z1z0, x2x0, y2y0, z2z0, nn[3]; + /* + I think this gives better results + */ - if (shade){ - x1x0 = (x[1]+Raise[0][1]) - (x[0]+Raise[0][0]); - y1y0 = (y[1]+Raise[1][1]) - (y[0]+Raise[1][0]); - z1z0 = (z[1]+Raise[2][1]) - (z[0]+Raise[2][0]); - x2x0 = (x[2]+Raise[0][2]) - (x[0]+Raise[0][0]); - y2y0 = (y[2]+Raise[1][2]) - (y[0]+Raise[1][0]); - z2z0 = (z[2]+Raise[2][2]) - (z[0]+Raise[2][0]); - nn[0] = y1y0 * z2z0 - z1z0 * y2y0 ; - nn[1] = z1z0 * x2x0 - x1x0 * z2z0 ; - nn[2] = x1x0 * y2y0 - y1y0 * x2x0 ; - /* BOF BOF - if(nn[2] < -0.1){ - nn[0] = -nn[0]; - nn[1] = -nn[1]; - nn[2] = -nn[2]; - } - */ - glNormal3dv(nn); - } - /* dessin de quandrangles "non convexes" */ + double x2[3]={x[2],x[3],x[0]}; + double y2[3]={y[2],y[3],y[0]}; + double z2[3]={z[2],z[3],z[0]}; - glBegin(GL_TRIANGLE_FAN); - glVertex3d(x[0]+Offset[0]+Raise[0][0], - y[0]+Offset[1]+Raise[1][0], - z[0]+Offset[2]+Raise[2][0]); - glVertex3d(x[1]+Offset[0]+Raise[0][1], - y[1]+Offset[1]+Raise[1][1], - z[1]+Offset[2]+Raise[2][1]); - glVertex3d(x[2]+Offset[0]+Raise[0][2], - y[2]+Offset[1]+Raise[1][2], - z[2]+Offset[2]+Raise[2][2]); - glVertex3d(x[3]+Offset[0]+Raise[0][3], - y[3]+Offset[1]+Raise[1][3], - z[3]+Offset[2]+Raise[2][3]); - glVertex3d(x[1]+Offset[0]+Raise[0][1], - y[1]+Offset[1]+Raise[1][1], - z[1]+Offset[2]+Raise[2][1]); - glEnd(); + Draw_Triangle(x,y,z,Offset,Raise,shade); + Draw_Triangle(x2,y2,z2,Offset,Raise,shade); + + return; } /* ------------------------------------------------------------------------ */ diff --git a/Graphics/Iso.cpp b/Graphics/Iso.cpp index 7c98e10320afc9f14d73eb021093a0751ba496d5..4fef32791cca5b3bda8dff6e9bc9cbc1beffdb5b 100644 --- a/Graphics/Iso.cpp +++ b/Graphics/Iso.cpp @@ -1,8 +1,9 @@ -// $Id: Iso.cpp,v 1.4 2001-01-08 08:05:43 geuzaine Exp $ +// $Id: Iso.cpp,v 1.5 2001-01-19 22:32:31 remacle Exp $ #include "Gmsh.h" #include "Mesh.h" #include "Draw.h" +#include "Numeric.h" void RaiseFill(int i, double Val, double ValMin, double Raise[3][5]); @@ -33,8 +34,9 @@ void Interpolate(double *X, double *Y, double *Z, void IsoSimplex(double *X, double *Y, double *Z, double *Val, double V, double Vmin, double Vmax, double *Offset, double Raise[3][5], int shade){ - int nb; + int nb,i; double Xp[6],Yp[6],Zp[6]; + double Xpi[6],Ypi[6],Zpi[6]; if(V != Vmax){ nb = 0; @@ -79,6 +81,57 @@ void IsoSimplex(double *X, double *Y, double *Z, double *Val, } } + /* + for having a nice isosurface, we should have n . grad v > 0 + n = normal to the polygon + v = unknown field we wanna draw + */ + + if(nb > 2) + { + double v1[3] = {Xp[2]-Xp[0],Yp[2]-Yp[0],Zp[2]-Zp[0]}; + double v2[3] = {Xp[1]-Xp[0],Yp[1]-Yp[0],Zp[1]-Zp[0]}; + double n[3]; + prodve(v1,v2,n); + //test + + // now get the gradient (simplified version of course) + + double xx = 0.0; + if(Val[2] != Val[1]) + { + double gr[3] = {X[2]-X[1],Y[2]-Y[1],Z[2]-Z[1]}; + double xx = gr[0] * n[0] + gr[1] * n[1] + gr[2] + n[2]; + if(Val[2] > Val[1]) xx = -xx; + } + if(Val[2] != Val[0]) + { + double gr[3] = {X[2]-X[0],Y[2]-Y[0],Z[2]-Z[0]}; + double xx = gr[0] * n[0] + gr[1] * n[1] + gr[2] + n[2]; + if(Val[2] > Val[0]) xx = -xx; + } + + + + if(xx > 0) + { + for(i=0;i<nb;i++) + { + Xpi[i] = Xp[i]; + Ypi[i] = Yp[i]; + Zpi[i] = Zp[i]; + } + for(i=0;i<nb;i++) + { + Xp[i] = Xpi[nb-i-1]; + Yp[i] = Ypi[nb-i-1]; + Zp[i] = Zpi[nb-i-1]; + } + } + } + + + if(nb == 3) Draw_Triangle(Xp,Yp,Zp,Offset,Raise,shade); else if(nb == 4) diff --git a/Makefile b/Makefile index 7862ad9cc487fcbf5361a57ce738c0f236f75571..a230a65102c28d067a09b9a88920b296c4a928fe 100644 --- a/Makefile +++ b/Makefile @@ -1,4 +1,4 @@ -# $Id: Makefile,v 1.41 2001-01-17 21:49:25 remacle Exp $ +# $Id: Makefile,v 1.42 2001-01-19 22:32:31 remacle Exp $ # ---------------------------------------------------------------------- # Makefile for Gmsh # ---------------------------------------------------------------------- @@ -7,7 +7,7 @@ MAKE = make CC = g++ - FLAGS = -g -Wall + FLAGS = -O RM = rm RMFLAGS = -f