Skip to content
Snippets Groups Projects
Commit 9708e82a authored by Jean-François Remacle's avatar Jean-François Remacle
Browse files

Enhencing isosurfaces

parent 634e6594
No related branches found
No related tags found
No related merge requests found
// $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;
}
}
......@@ -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
......
// $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];
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];
}
/*
I think this gives better results
*/
glNormal3dv(nn);
}
/* dessin de quandrangles "non convexes" */
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();
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]};
Draw_Triangle(x,y,z,Offset,Raise,shade);
Draw_Triangle(x2,y2,z2,Offset,Raise,shade);
return;
}
/* ------------------------------------------------------------------------ */
......
// $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)
......
# $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
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment