#include <stdio.h>
#include <iostream.h>
#include <stdlib.h>
#include <fstream.h>
#include "mailleur.h"
typedef struct _vect{
float xx,yy,zz;
} vect;
//*******************************************************************
//*******************************************************************
float norme(vect u){
return(sqrt(u.xx*u.xx+u.yy*u.yy+u.zz*u.zz));
}
//*********definition du constructeur**********************************
maillage::maillage(float test){
//init des vecteurs x et y
for (int i=0;i<=N-1;i++){
x[i]=X*i/(N-1);
y[i]=Y*i/(N-1);
}
//init des triangles
for(int i=0;i<=N-2;i++){
for(int j=0;j<=N-2;j++){
// premier triangle
tabtri[2*i*(N-1)+2*(j+1)-2].a.x=&x[j];
tabtri[2*i*(N-1)+2*(j+1)-2].a.y=&y[i+1];
tabtri[2*i*(N-1)+2*(j+1)-2].a.z=&z[j][i+1];
tabtri[2*i*(N-1)+2*(j+1)-2].b.x=&x[j];
tabtri[2*i*(N-1)+2*(j+1)-2].b.y=&y[i];
tabtri[2*i*(N-1)+2*(j+1)-2].b.z=&z[j][i];
tabtri[2*i*(N-1)+2*(j+1)-2].c.x=&x[j+1];
tabtri[2*i*(N-1)+2*(j+1)-2].c.y=&y[i+1];
tabtri[2*i*(N-1)+2*(j+1)-2].c.z=&z[j+1][i+1];
// deuxieme triangle
tabtri[2*i*(N-1)+2*(j+1)-1].a.x=&x[j+1];
tabtri[2*i*(N-1)+2*(j+1)-1].a.y=&y[i];
tabtri[2*i*(N-1)+2*(j+1)-1].a.z=&z[j+1][i];
tabtri[2*i*(N-1)+2*(j+1)-1].b.x=&x[j+1];
tabtri[2*i*(N-1)+2*(j+1)-1].b.y=&y[i+1];
tabtri[2*i*(N-1)+2*(j+1)-1].b.z=&z[j+1][i+1];
tabtri[2*i*(N-1)+2*(j+1)-1].c.x=&x[j];
tabtri[2*i*(N-1)+2*(j+1)-1].c.y=&y[i];
tabtri[2*i*(N-1)+2*(j+1)-1].c.z=&z[j][i];
}
}
// cout << "nono\n";
}
//******Definition des fonctions menmbres de la classe point************
//***********************************************************************
float maillage::df(float x1,float y1,float z1,float x2,float y2,float z2
,float x3,float y3,float z3){
float u,u1,u2,u3,u1prim,u2prim;
u1=(y2-y1)*(z3-z1)-(z2-z1)*(y3-y1);
u2=(z2-z1)*(x3-x1)-(x2-x1)*(z3-z1);
u3=(x2-x1)*(y3-y1)-(y2-y1)*(x3-x1);
u=u1*u1+u2*u2+u3*u3;
u1prim=y3-y2;
u2prim=x2-x3;
return((u1prim*u1+u2prim*u2)/(2*sqrt(u)));
}
//*************************************************************************
float maillage::surface_triangle(int i){
vect ab,ac,u;
ab.xx=*tabtri[i].b.x-*tabtri[i].a.x;
ab.yy=*tabtri[i].b.y-*tabtri[i].a.y;
ab.zz=*tabtri[i].b.z-*tabtri[i].a.z;
ac.xx=*tabtri[i].c.x-*tabtri[i].a.x;
ac.yy=*tabtri[i].c.y-*tabtri[i].a.y;
ac.zz=*tabtri[i].c.z-*tabtri[i].a.z;
u.xx=ab.yy*ac.zz-ab.zz*ac.yy;
u.yy=ab.zz*ac.xx-ab.xx*ac.zz;
u.zz=ab.xx*ac.yy-ab.yy*ac.xx;
return(norme(u)/2.);
}
//**************************************************************************
//genere un fichier executable par scilab contenant
//x mat (1xN) y mat (1xN) z mat (NxN,XxY )
//****************************************************************************************
void maillage::results()
{
ofstream sortie("resultats",ios::out);
sortie << "z1=[";
for(int j=0;j<N-1;j++){
sortie << z[0][j] << ",";
}
sortie << z[0][N-1] << "];" << "\n";
sortie <<"z=z1;" << "\n";
for(int i=1;i<N;i++){
sortie << "z1=[";
for(int j=0;j<N;j++){
sortie << z[i][j];
if(j<N-1) sortie<<",";
}
sortie << "];" << "\n";
sortie << "z=[z;z1];" << "\n";
}
sortie << "x=[";
for(int i=0;i<N;i++){
sortie << X*i/(N-1);
if(i<N-1) sortie<<",";
}
sortie << "];" << "\n";
sortie << "y=[";
for(int i=0;i<N;i++){
sortie << Y*i/(N-1);
if(i<N-1) sortie << ",";
}
sortie << "];" << "\n";
sortie << "xbasc(0);" << "\n";
sortie << "plot3d(x,y,z);" << "\n";
sortie.close();
}
//****************************************************************************
void maillage::refresh_curv(){
float delta=X/(N-1);
float a=-4./(X*X*X);
float b=-3.*a*X/2.;
float temp;
for(int i=0;i<N;i++) for(int j=0;j<N;j++) z[i][j]=0.0; //tout le mode a zero
for(int i=0;i<N;i++){
temp=a*(i*delta)*(i*delta)*(i*delta)+b*(i*delta)*(i*delta);
z[0][i]=temp;
z[i][0]=temp;
z[N-1][N-i-1]=temp;
z[N-1-i][N-1]=temp;
}
}
//******************************************************************************
void maillage::refresh()
{
for(int i=0;i<N;i++)
for(int j=0;j<N;j++)
z[i][j]=0.0;
ifstream cond_init("cond_init",ios::in);
cond_init >> z[0][0] >> z[(N-1)/2][0] >> z[N-1][0]
>> z[0][(N-1)/2] >> z[N-1][(N-1)/2]
>> z[0][N-1] >> z[(N-1)/2][N-1] >> z[N-1][N-1];
cond_init.close();
float delta;
delta=(-z[0][0]+z[(N-1)/2][0])/(N-1)*2;
for(int i=1;i<N-1;i++)
if(i==(N-1)/2) delta=(-z[(N-1)/2][0]+z[N-1][0])/(N-1)*2;
else z[i][0]=z[i-1][0]+delta;
delta=(-z[0][N-1]+z[(N-1)/2][N-1])/(N-1)*2;
for(int i=1;i<N-1;i++)
if(i==(N-1)/2) delta=(-z[(N-1)/2][N-1]+z[N-1][N-1])/(N-1)*2;
else z[i][N-1]=z[i-1][N-1]+delta;
delta=(-z[0][0]+z[0][(N-1)/2])/(N-1)*2;
for(int i=1;i<N-1;i++)
if(i==(N-1)/2) delta=(-z[0][(N-1)/2]+z[0][N-1])/(N-1)*2;
else z[0][i]=z[0][i-1]+delta;
delta=(-z[N-1][0]+z[N-1][(N-1)/2])/(N-1)*2;
for(int i=1;i<N-1;i++)
if(i==(N-1)/2) delta=(-z[N-1][(N-1)/2]+z[N-1][N-1])/(N-1)*2;
else z[N-1][i]=z[N-1][i-1]+delta;
}
//*************************************************************************
float maillage::surface(){
int i;
float s;
s=0.;
for(i=0;i<2*(N-1)*(N-1);i++) s=s+surface_triangle(i);
return(s);
}
//**************************************************************************
void maillage::calculer_gradient(){
int i,j;
for(i=0;i<N-2;i++){
for(j=0;j<N-2;j++){
gradient[i][j]=df(x[i+1],y[j+1],z[i+1][j+1],x[i],y[j+1]
,z[i][j+1],x[i],y[j],z[i][j]);
gradient[i][j]+=df(x[i+1],y[j+1],z[i+1][j+1],x[i+1],y[j]
,z[i+1][j],x[i],y[j],z[i][j]);
gradient[i][j]+=df(x[i+1],y[j+1],z[i+1][j+1],x[i+1],y[j]
,z[i+1][j],x[i+2],y[j+1],z[i+2][j+1]);
gradient[i][j]+=df(x[i+1],y[j+1],z[i+1][j+1],x[i+2],y[j+2],
z[i+2][j+2],x[i+2],y[j+1],z[i+2][j+1]);
gradient[i][j]+=df(x[i+1],y[j+1],z[i+1][j+1],x[i+2],y[j+2]
,z[i+2][j+2],x[i+1],y[j+2],z[i+1][j+2]);
gradient[i][j]+=df(x[i+1],y[j+1],z[i+1][j+1],x[i],y[j+1],
z[i][j+1],x[i+1],y[j+2],z[i+1][j+2]);
}
}
}