#include <stdlib.h>

R TailleArete=0.1; // >0 => Taille souhaitee des aretes du maillage
// =0 => Emploi de la fonction Taille_Ideale() (A FAIRE...)
// <0 => Emploi d'une carte de taille des aretes(A FAIRE...)
// sur un maillage de fond

int NoPoint=-1;
int NoLigne=-1;

R3 PtMIN(0.0, 0.0, 0.0); //Pt MIN par defaut Points de coordonnees extremales
R3 PtMAX(1.5, 1.1, 0.0); //Pt MAX par defaut

XPoint pxy0={0,0};
XPoint pxy;

Drawable fenoupix; //fenetre ou pixmap variable auxiliaire

int TraceTg=0; //Trace (1) ou non (0) des tangentes aux noeuds des courbes NURBS
int TraceNormales=0; //Trace (1) ou non (0) des normales et bi-normales aux noeuds des courbes NURBS

int NoSurface = -1;
int Ombrage = 2; //Ombrage Gouraud si 1, Filaire si 0, Gouraud-Filaire si 2

void NommerPoint( Nom & nom )
{
int nbc;

NoPoint++;
strcpy( nom, "P" );
itoa( NoPoint, nbc, nom+1 );
}

void NommerLigne( Nom & nom )
{
int nbc;

NoLigne++;
strcpy( nom, "L" );
itoa( NoLigne, nbc, nom+1 );
}

void NommerSurface( Nom & nom )
{
int nbc;
NoSurface++;
strcpy( nom, "S" );
itoa( NoSurface, nbc, nom+1 );
}

Surface::Surface() : PLSV ()
{
NommerSurface( nom );
PSuivant = NULL;
}

Surface::Surface( Nom NomSurface ) : PLSV()
{
strcpy( nom, NomSurface );
PSuivant = NULL;
}

Surface::~Surface()
{
if ( NoStEF ) delete [] NoStEF;
if ( XyzSommet ) delete [] XyzSommet;
if ( U ) delete [] U;
if ( V ) delete [] V;
if ( PC ) delete [] PC;
if ( Poids ) delete [] Poids;
if ( Courbes ) delete [] Courbes;

if ( PSuivant )
delete PSuivant;
}

R3 Surface::PointMilieuSurface(void)
{
if ( (XyzSommet == NULL) || (NoStEF == NULL) )
return R3(0,0,0);
int pos = (NbStEF*NbEF) / 2 ;
return ( (XyzSommet[ NoStEF[pos] ] + XyzSommet[ NoStEF[pos-1] ]) / 2 );
}

void Surface::PointSurface( R u, R v, R3 & P )
{
int i,j;
int k,l, m, n;
R Nu[1+64], Nv[1+64];
R denominateur = 0;
R Np2;
R denominateur2;

i = IntervalleNURBS( u, U, nbnu, degreu );
j = IntervalleNURBS( v, V, nbnv, degrev );

ValeurBaseNURBS( i, degreu, u, U, Nu );
ValeurBaseNURBS( j, degrev, v, V, Nv );

P.x = P.y = P.z = 0.0;

if ( Poids == NULL )
{
for (k=0; k<=degreu; k++)
{
l = i - degreu + k;
R3 P2(0,0,0);
for (m=0; m<=degrev; m++)
{
n = j - degrev + m;
P2 += Nv[m] * PC[l][n]->XyzSommet[0];
}
P += Nu[k] * P2;
}
}
else
{
for (k=0; k<=degreu; k++)
{
l = i - degreu + k;
R3 P2(0,0,0);
denominateur2 = 0;
for (m=0; m<=degrev; m++)
{
n = j - degrev + m;
Np2 = Nv[m] * Poids[l][n];
denominateur2 += Np2;
P2 += Np2 * PC[l][n]->XyzSommet[0];
}
P += Nu[k] * P2;
denominateur += Nu[k] * denominateur2;
}
P = P / denominateur;
}
}

int Surface::Mailler( Points * LesPoints )
{
R u, v, hu, hv;
int nbefu, nbefv; /* Nombre d'elements finis selon u et v */
bool renverse = false;

if ( TailleArete <= 0 )
return 2;

CalculNbEF( nbefu, nbefv ); /* fonction de la classe derivee */

NbEF = nbefu * nbefv;
if ( NbEF < 1 ) NbEF = 1;

NbSommet = (nbefu+1) * (nbefv+1);
XyzSommet = new R3[NbSommet];

hu = ( U[nbnu+1] - U[0] ) / nbefu;
hv = ( V[nbnv+1] - V[0] ) / nbefv;

u = U[0];
v = V[0];

for (int i=0; i<=nbefv; i++) {
u = U[0];
for (int j=0; j<=nbefu; j++) {
int pos = i*(nbefu+1)+j;
PointSurface( u, v, XyzSommet[pos] );
PtMIN = Min( PtMIN, XyzSommet[pos] );
PtMAX = Max( PtMAX, XyzSommet[pos] );
u += hu;
}
v += hv;
}

NbStEF = 4;
NoStEF = new Z[NbEF*NbStEF];

R3 vect1;
R3 vect2;
vect1 = XyzSommet[1]-XyzSommet[0];
vect2 = XyzSommet[nbefu+2]-XyzSommet[0];
if ( ( vect1.x*vect2.y - vect1.y*vect2.x ) < 0 ) /* calcul de l orientation d une surface*/
renverse = true;

int m = 0;
int no=0;

if ( !renverse )
{
for (int i=0; i<nbefv; i++)
{
for (int j=0; j<nbefu; j++)
{
int pos = j+i*(nbefu+1);
no++;

NoStEF[m++] = pos;
NoStEF[m++] = pos + 1; /* quadrangle parcouru dans le sens de la normale si elle est exterieure */
NoStEF[m++] = pos + nbefu + 2;
NoStEF[m++] = pos + nbefu + 1;
}
}
}
else
{
for (int i=0; i<nbefv; i++)
{
for (int j=0; j<nbefu; j++)
{
int pos = j+i*(nbefu+1);
no++;
NoStEF[m++] = pos; /* quadrangle parcouru dans l autre sens */
NoStEF[m++] = pos + nbefu + 1;
NoStEF[m++] = pos + nbefu + 2;
NoStEF[m++] = pos + 1;
}
}
}
return 0;
}

void Surface::TracerMaillage( MemoirePx3d & mpx )
{
R3 xyz;
R3 coords[4];
Z m = 0;
char texte[8];
R3 normale;
R luminosite;
R3 oeil = mpx.axoeil;
R3 ptv = mpx.axoptv;
R3 barycentre[NbEF];
Z noanc[NbEF];
R d[NbEF];
R3 tmp(0.0,0.0,0.0);

TracerPolyedre(mpx);

VuChoisirFonte ( hauttexte );
VuEpaisseur( 2 );

for (int i=0; i<NbEF ; i++)
{
barycentre[i]=tmp;
for (int j=0; j<NbStEF; j++)
{
int index = NoStEF[m++];
coords[j] = XyzSommet[index];
barycentre[i] += coords[j];
}
barycentre[i] = 0.25*barycentre[i];
d[i]=Norme2(oeil-barycentre[i]);
noanc[i]=i;
}
tritas(NbEF,d,noanc); /* suivant la distance du barycentre des quadrangles a la pos de loeil la fonction tritas ordonne les faces */

Z pos;

for (int i=NbEF-1; i>=0 ; i--)
{
pos=noanc[i];
pos *= NbStEF;
for (int j=0; j<NbStEF; j++)
{
int index = NoStEF[pos++];
coords[j] = XyzSommet[index];
}

if ( Ombrage )
{
normale = ( coords[1] - coords[0] ) ^ ( coords[3] - coords[0] ); /* calcul des quatres normales afin d affecter une couleur a chaque point du maillage */
luminosite = ( normale, oeil - ptv );
luminosite /= ( Norme2( normale ) * Norme2( oeil - ptv ) );

normale = ( coords[1] - coords[0] ) ^ ( coords[3] - coords[0] );
R Ca = ( normale, oeil - ptv );
Ca /= ( Norme2( normale ) * Norme2( oeil - ptv ) );
Ca = 127+255*fabs(Ca/2);

normale = ( coords[0] - coords[1] ) ^ ( coords[2] - coords[1] );
R Cb = ( normale, oeil - ptv );
Cb /= ( Norme2( normale ) * Norme2( oeil - ptv ) );
Cb = 127+255*fabs(Cb/2);

normale = ( coords[3] - coords[2] ) ^ ( coords[1] - coords[2] );
R Cc = ( normale, oeil - ptv );
Cc /= ( Norme2( normale ) * Norme2( oeil - ptv ) );
Cc = 127+255*fabs(Cc/2);

normale = ( coords[0] - coords[3] ) ^ ( coords[2] - coords[3] );
R Cd = ( normale, oeil - ptv );
Cd /= ( Norme2( normale ) * Norme2( oeil - ptv ) );
Cd = 127+255*fabs(Cd/2);

R grad = fabs(luminosite/2);
VuPaletteCouleurs( 11 );
VuCouleur( (int) (127+255*grad) );
// mpx.VuFace( 4, coords ); affichage par quadrangle
mpx.VuTriangleCouleur (coords[0],coords[1],coords[2],Ca,Cb,Cc); /* trace le triangle en tenant compte des couleurs intermediaires entre les sommets */
mpx.VuTriangleCouleur (coords[2],coords[3],coords[0],Cc,Cd,Ca);

if (Ombrage==2)
{
VuCouleur( Orange );
VuEpaisseur( 1 );
// mpx.VuTraits( 4, coords ); affichage par quadrangle
//mpx.VuBordTriangle(coords[0],coords[1],coords[2]);
//mpx.VuBordTriangle(coords[2],coords[3],coords[0]);

}

}
else
{
VuEpaisseur( 2 );
VuCouleur( Bleu );
mpx.VuTraits( 4, coords );
}
}
xyz = PointMilieuSurface();
sprintf( texte, "*%s", nom);
VuCouleur( Noir );
mpx.VuSymbole( xyz, texte );
}

void Surface::VerifierMaillage( Points * LesPoints )
{
for (int i=0; i<=nbnu; i++)
{
for (int j=0; j<=nbnv; j++)
{
if ( (PC[i][j]->date) > this->date ) /* controle de nonmodification posterieur au dernier calcul */
{
Mailler(LesPoints);
return;
}
}
}
}

void Surface::TracerPolyedre( MemoirePx3d & mpx )
{
VuEpaisseur( 2 );
VuCouleur( Rose );
for (int i=0; i<nbnu; i++ )
for (int j=0; j<nbnv; j++)
{
mpx.VuTrait( PC[i][j]->XyzSommet[0], PC[i+1][j]->XyzSommet[0] );
mpx.VuTrait( PC[i][j]->XyzSommet[0], PC[i][j+1]->XyzSommet[0] ); /* trace du maillage de pts de controle */
}
for (int i=0; i<nbnu; i++ )
mpx.VuTrait( PC[i][nbnv]->XyzSommet[0], PC[i+1][nbnv]->XyzSommet[0] );
for (int j=0; j<nbnv; j++)
mpx.VuTrait( PC[nbnu][j]->XyzSommet[0], PC[nbnu][j+1]->XyzSommet[0] );

for (int i=0;i<=nbnu;i++)
for (int j=0;j<=nbnv;j++)
PC[i][j]->TracerMaillage(mpx);
}

//=============================================================================
// Classe SurfaceNURBS
//=============================================================================

SurfaceNURBS::SurfaceNURBS(Lignes *LesLignes,Points *LesPoints):Surface()
{
Nom nomP;
Ligne *Base1,*Base2;
int ncvals;

mpxInvite.VuInvite( "Selectionnez la courbe NURBS n°1 de base ou tapez son nom : ");
if ( (Base1 = LesLignes->Lignes::CliquerLigne()) == NULL )
return;
mpxInvite.VuEffacerMemoirePx();

mpxInvite.VuInvite( "Selectionnez la courbe NURBS n°2 de base ou tapez son nom : ");
if ( (Base2 = LesLignes->Lignes::CliquerLigne()) == NULL )
return;
mpxInvite.VuEffacerMemoirePx();

nbnu=Base1->nbn;
nbnv=Base2->nbn;
degreu=Base1->degre;
degrev=Base2->degre;

//Initialisation des tableaux des points de controle et des poids

PC = new Point **[1+nbnu];
for (int i=0; i<=nbnu; i++)
PC[i] = new Point *[1+nbnv];

for (int i=0; i<=nbnu; i++)
{
NommerPoint(nomP);
PC[i][0]=new Pointxyz(nomP,
Base1->PC[i]->XyzSommet[0].x,
Base1->PC[i]->XyzSommet[0].y,
Base1->PC[i]->XyzSommet[0].z); /* affectation post rentree des infos par l utilisateur */
LesPoints->Ajouter1Point(PC[i][0]);
}
for (int j=1; j<=nbnv; j++)
for (int i=0; i<=nbnu; i++)
{
NommerPoint(nomP);
PC[i][j]=new Pointxyz(nomP,
PC[i][j-1]->XyzSommet[0].x+Base2->PC[j]->XyzSommet[0].x-Base2->PC[j-1]->XyzSommet[0].x,
PC[i][j-1]->XyzSommet[0].y+Base2->PC[j]->XyzSommet[0].y-Base2->PC[j-1]->XyzSommet[0].y,
PC[i][j-1]->XyzSommet[0].z+Base2->PC[j]->XyzSommet[0].z-Base2->PC[j-1]->XyzSommet[0].z);
LesPoints->Ajouter1Point(PC[i][j]);
}

if ( (Base1->Poids) && (Base2->Poids) )
{
Poids = new R *[1+nbnu];
for (int i=0; i<=nbnu; i++)
Poids[i] = new R[1+nbnv]; /* creationdu tableau des noeuds s ils existent */

for (int i=0; i<=nbnu; i++)
for (int j=0; j<=nbnv; j++)
Poids[i][j]=((Base1->Poids[i])+(Base2->Poids[j]))/2;
}

U = new R [1+nbnu+1+degreu];
for (int i=0; i<=1+nbnu+degreu; i++)
U[i] = Base1->U[i];

V = new R [1+nbnv+1+degrev];
for (int j=0; j<=1+nbnv+degrev; j++)
V[j] = Base2->U[j];

Courbes = new Ligne *[2];
Courbes[0] = Base1;
Courbes[1] = Base2;
Mailler( LesPoints );

}

SurfaceNURBS::SurfaceNURBS(Points *LesPoints,Lignes *LesLignes):Surface()
{
int nbc, ncvals;
Nom NomPt;
Point * PtAvant;
char texte[48];
char carlu;

mpxInvite.VuInvite( "Frapper le DEGRE en u de la surface NURBS" ); /* interface utilisateur habituelle */
do
{
mpxInvite.lireZ( ncvals, degreu );
if(ncvals==-1) {nbnu=0; return;}
}
while( degreu<=0 );

mpxInvite.VuInvite( "Frapper le DEGRE en v de la surface NURBS" );
do
{
mpxInvite.lireZ( ncvals, degrev );
if(ncvals==-1) {nbnv=0; return;}
}
while( degrev<=0 );

mpxInvite.VuInvite( "Frapper le NOMBRE de POINTS en u de la surface NURBS" );
do
{
mpxInvite.lireZ( ncvals, nbnu );
if(ncvals==-1) {nbnu=0; return;}
}
while( nbnu<=0 );
nbnu--;

mpxInvite.VuInvite( "Frapper le NOMBRE de POINTS en v de la surface NURBS" );
do
{
mpxInvite.lireZ( ncvals, nbnv );
if(ncvals==-1) {nbnv=0; return;}
}
while( nbnv<=0 );
nbnv--;

PC = new Point **[1+nbnu];
for (int i=0; i<=nbnu; i++)
PC[i] = new Point *[1+nbnv];

strcpy( texte, "Cliquer le POINT de CONTROLE[" );
for (int i=0; i<=nbnu; i++)
for (int j=0; j<=nbnv; j++)
{
itoa( i, nbc, (texte+29) );
strcat( texte, "][");
itoa( j, nbc, (texte+32) );
strcat(texte, "] ou frapper son NOM" );
mpxInvite.VuInvite( texte );
PC[i][j] = LesPoints->Points::CliquerPoint();
if( PC[i][j] == NULL ) { nbnu=0; nbnv=0; return; }
mpxInvite.VuEffacerMemoirePx();
}

U = new R[1+nbnu+1+degreu]; // declaration du tableau des noeuds
mpxInvite.VuInvite( "Frapper la valeur des NOEUDS u (o/N)" );
carlu = mpxInvite.lire1Caractere();
mpxInvite.VuEffacerMemoirePx();
if( carlu == '@' || int(carlu)== 27 ) {nbnu=0; return;} //ABANDON
if( carlu == 'o' || carlu == 'O' || carlu == '0' )
{ //saisie des valeurs des noeuds
strcpy( texte, "Frapper la valeur du noeud U[" ); //comme pour la saisie des lignes nurbs dans plsv
for (int i=degreu; i<=nbnu+1; i++)
{
itoa( i, nbc, (texte+29) ); //texte de l'invite
strcat( texte, "]" );
mpxInvite.VuInvite( texte ); //trace de l'invite
do {
mpxInvite.lireR( ncvals, U[i] );
if(ncvals==-1) {nbnu=0; return;} //ABANDON demande
if(i==0) break;
} while (U[i]<U[i-1]);
mpxInvite.VuEffacerMemoirePx();
}
}
else
{ // creation automatique des noeuds simples
U[degreu] = 0;
for (int i=degreu+1; i<=nbnu+1; i++)
{ U[i] = U[i-1] + Norme2( PC[i-degreu ][0]->XyzSommet[0]
- PC[i-degreu-1][0]->XyzSommet[0] ); }
}

for (int i=0; i<=degreu; i++ ) //noeud initial
{ U[i] = U[degreu]; }

for (int i=nbnu+2; i<=nbnu+1+degreu; i++) //noeud final
{ U[i]=U[nbnu+1]; }

V = new R[1+nbnv+1+degrev]; //declaration du tableau des noeuds
mpxInvite.VuInvite( "Frapper la valeur des NOEUDS v (o/N)" );
carlu = mpxInvite.lire1Caractere();
mpxInvite.VuEffacerMemoirePx();
if( carlu == '@' || int(carlu)== 27 ) {nbnv=0; return;} //ABANDON
if( carlu == 'o' || carlu == 'O' || carlu == '0' )
{ //saisie interactive des valeurs des noeuds
strcpy( texte, "Frapper la valeur du noeud V[" );
for (int j=degrev; j<=nbnv+1; j++)
{
itoa( j, nbc, (texte+29) ); //texte de l'invite
strcat( texte, "]" );
mpxInvite.VuInvite( texte ); //trace de l'invite
do {
mpxInvite.lireR( ncvals, V[j] );
if(ncvals==-1) {nbnv=0; return;} //ABANDON demande
if(j==0) break;
} while (V[j]<V[j-1]);
mpxInvite.VuEffacerMemoirePx();
}
}
else
{ // creation automatique des noeuds simples
V[degrev] = 0;
for (int j=degrev+1; j<=nbnv+1; j++)
{ V[j] = V[j-1] + Norme2( PC[0][j-degrev ]->XyzSommet[0]
- PC[0][j-degrev-1]->XyzSommet[0] ); }
}

for (int j=0; j<=degrev; j++ ) //noeud initial
{ V[j] = V[degrev]; }

for (int j=nbnv+2; j<=nbnv+1+degrev; j++) //noeud final
{ V[j]=V[nbnv+1]; }

mpxInvite.VuInvite( "Frapper la valeur des Poids (o/N)" );
carlu = mpxInvite.lire1Caractere();
if( carlu == '@' || int(carlu)== 27 ) {nbnv=0; return;} //ABANDON
if( carlu == 'o' || carlu == 'O' || carlu == '0' )
{
//saisie interactive des valeurs des Poids => courbe NURBS

Poids = new R *[1+nbnu];
for (int i=0; i<=nbnu; i++)
Poids[i] = new R[1+nbnv];

strcpy( texte, "Frapper le Poids[" );
for (int i=0; i<=nbnu; i++)
for (int j=0; j<=nbnv; j++)
{
itoa( i, nbc, (texte+17) ); //texte de l'invite
strcat( texte, "][" );
itoa( j, nbc, (texte+17) );
strcat( texte, "]" );
mpxInvite.VuInvite( texte );
mpxInvite.lireR( ncvals, Poids[i][j] );
if(ncvals==-1) {nbnu=0; nbnv=0; return;} //ABANDON demande
}
}
else
{ // pas de creation des Poids => courbe polynomiale
// Poids = NULL; cela signifie non pas que les poids soient nuls mais que le tableau des poids est vide

Poids = new R *[1+nbnu];
for (int i=0; i<=nbnu; i++)
Poids[i] = new R[1+nbnv];

for (int i=0; i<=nbnu; i++)
for (int j=0; j<=nbnv; j++)
Poids[i][j]=1; //valeur par defaut
}

mpxInvite.VuEffacerMemoirePx(); //on libere la varible qui stocke les donnees d interface

Mailler(LesPoints);
}

void SurfaceNURBS::CalculNbEF( int & nbefu, int & nbefv )
{
nbefu = (int) ( Courbes[0]->LongueurNURBS() / TailleArete + 0.5 );
nbefv = (int) ( Courbes[1]->LongueurNURBS() / TailleArete + 0.5 );
}

void SurfaceNURBS::AfficherDefinition()
{
cout << "Nom de la Surface NURBS= " << nom <<endl;
}

void SurfaceNURBS::VerifierMaillage(Points *LesPoints)
{
if( (Courbes[0]->date > this->date) || (Courbes[1]->date > this->date))
Mailler(LesPoints);
}

//=============================================================================
// Classe SurfaceBezier
//=============================================================================

SurfaceBezier::SurfaceBezier(Lignes *LesLignes,Points *LesPoints):Surface() //de maniere tres similaire a la classe precedante
{
Nom nomP;
Ligne *Base1,*Base2;
int ncvals;

mpxInvite.VuInvite( "Selectionnez la courbe NURBS n°1 de base ou tapez son nom : ");
if ( (Base1 = LesLignes->Lignes::CliquerLigne()) == NULL )
return;
mpxInvite.VuEffacerMemoirePx();

mpxInvite.VuInvite( "Selectionnez la courbe NURBS n°2 de base ou tapez son nom : ");
if ( (Base2 = LesLignes->Lignes::CliquerLigne()) == NULL )
return;
mpxInvite.VuEffacerMemoirePx();

nbnu=Base1->nbn;
nbnv=Base2->nbn;
degreu=Base1->degre;
degrev=Base2->degre;

//Initialisation des points de controle et des poids

PC = new Point **[1+nbnu];
for (int i=0; i<=nbnu; i++)
PC[i] = new Point *[1+nbnv];

for (int i=0; i<=nbnu; i++)
{
NommerPoint(nomP);
PC[i][0]=new Pointxyz(nomP,
Base1->PC[i]->XyzSommet[0].x,
Base1->PC[i]->XyzSommet[0].y,
Base1->PC[i]->XyzSommet[0].z);
LesPoints->Ajouter1Point(PC[i][0]);
}
for (int j=1; j<=nbnv; j++)
for (int i=0; i<=nbnu; i++)
{
NommerPoint(nomP);
PC[i][j]=new Pointxyz(nomP,
PC[i][j-1]->XyzSommet[0].x+Base2->PC[j]->XyzSommet[0].x-Base2->PC[j-1]->XyzSommet[0].x,
PC[i][j-1]->XyzSommet[0].y+Base2->PC[j]->XyzSommet[0].y-Base2->PC[j-1]->XyzSommet[0].y,
PC[i][j-1]->XyzSommet[0].z+Base2->PC[j]->XyzSommet[0].z-Base2->PC[j-1]->XyzSommet[0].z);
LesPoints->Ajouter1Point(PC[i][j]);
}

if ( (Base1->Poids) && (Base2->Poids) )
{
Poids = new R *[1+nbnu];
for (int i=0; i<=nbnu; i++)
Poids[i] = new R[1+nbnv];

for (int i=0; i<=nbnu; i++)
for (int j=0; j<=nbnv; j++)
Poids[i][j]=((Base1->Poids[i])+(Base2->Poids[j]))/2;
}

U = new R [1+nbnu+1+degreu];
for (int i=0; i<=(1+nbnu+degreu)/2; i++)
U[i] = 0;

for (int i=1+(1+nbnu+degreu)/2; i<=1+nbnu+degreu; i++)
U[i] = 1;

V = new R [1+nbnv+1+degrev];
for (int j=0; j<=(1+nbnu+degreu)/2; j++)
V[j] = 0;

for (int j=1+(1+nbnu+degreu)/2; j<=1+nbnu+degreu; j++)
V[j] = 1;

Courbes = new Ligne *[2];
Courbes[0] = Base1;
Courbes[1] = Base2;
Mailler( LesPoints );

}

void SurfaceBezier::CalculNbEF( int & nbefu, int & nbefv )
{
nbefu = (int) ( Courbes[0]->LongueurNURBS() / TailleArete + 0.5 );
nbefv = (int) ( Courbes[1]->LongueurNURBS() / TailleArete + 0.5 );
}

void SurfaceBezier::AfficherDefinition()
{
cout << "Nom de la Surface NURBS= " << nom <<endl;
}

void SurfaceBezier::VerifierMaillage(Points *LesPoints)
{
if( (Courbes[0]->date > this->date) || (Courbes[0]->date > this->date))
Mailler(LesPoints);
}

//=============================================================================
// Classe SurfaceReglee
//=============================================================================

SurfaceReglee::SurfaceReglee(Lignes *LesLignes,Points *LesPoints):Surface()
{
Nom nomP;
Ligne *Base1,*Base2;
int ncvals;

mpxInvite.VuInvite( "Selectionnez la premiere courbe de base ou tapez son nom : ");
if ( (Base1 = LesLignes->Lignes::CliquerLigne()) == NULL )
return;
mpxInvite.VuEffacerMemoirePx();

mpxInvite.VuInvite( "Selectionnez la deuxieme courbe de base ou tapez son nom : ");
if ( (Base2 = LesLignes->Lignes::CliquerLigne()) == NULL )
return;
mpxInvite.VuEffacerMemoirePx();

Z deg1,deg2;
deg1=Base1->degre;
deg2=Base2->degre;
while (Base1->degre!=Base2->degre)
{
if (Base1->degre < Base2->degre)
Base1->Elevation1Degre(LesPoints);
else
Base2->Elevation1Degre(LesPoints);

if ( (deg1==Base1->degre) && (deg2==Base2->degre))
break;
}

if (Base1->degre!=Base2->degre)
{
mpxInvite.VuInvite( "Degré non élevable");
return;
}

nbnu=Base1->nbn; //les chemins entre les deux courbes qui se font face sont tendus de degre 1
nbnv=1;
degreu=Base1->degre;
degrev=1;

PC = new Point **[1+nbnu];
for (int i=0; i<=nbnu; i++)
PC[i] = new Point *[1+nbnv];

for (int i=0; i<=nbnu; i++)
{
NommerPoint(nomP);
PC[i][0]=new Pointxyz(nomP,
Base1->PC[i]->XyzSommet[0].x,
Base1->PC[i]->XyzSommet[0].y,
Base1->PC[i]->XyzSommet[0].z);
LesPoints->Ajouter1Point(PC[i][0]);

NommerPoint(nomP);
PC[i][1]=new Pointxyz(nomP,
Base2->PC[i]->XyzSommet[0].x,
Base2->PC[i]->XyzSommet[0].y,
Base2->PC[i]->XyzSommet[0].z);
LesPoints->Ajouter1Point(PC[i][1]);
}

if ( (Base1->Poids) && (Base2->Poids) )
{
Poids = new R *[1+nbnu];
for (int i=0; i<=nbnu; i++)
Poids[i] = new R[1+nbnv];

for (int i=0; i<=nbnu; i++)
for (int j=0; j<=nbnv; j++)
Poids[i][j]=((Base1->Poids[i])+(Base2->Poids[j]))/2;
}

U = new R [1+nbnu+1+degreu];
for (int i=0; i<=1+nbnu+degreu; i++)
U[i] = Base1->U[i];

V = new R [1+nbnv+1+degrev];
V[0] = V[1] = 0;
V[2] = V[3] = 1;

Courbes = new Ligne *[1];
Courbes[0] = Base1;
Mailler( LesPoints );

}

void SurfaceReglee::CalculNbEF( int & nbefu, int & nbefv )
{
nbefu = (int) ( Courbes[0]->LongueurNURBS() / TailleArete + 0.5 );
nbefv=nbefu;
}

void SurfaceReglee::AfficherDefinition()
{
cout << "Nom de la Surface REGLEE= " << nom <<endl;
}

void SurfaceReglee::VerifierMaillage(Points *LesPoints)
{
if( (Courbes[0]->date > this->date) || (Courbes[0]->date > this->date))
Mailler(LesPoints);
}

//=============================================================================
// Classe SurfaceBilineaire
//=============================================================================
SurfaceBilineaire::SurfaceBilineaire(Lignes *LesLignes,Points *LesPoints):Surface()
{
Nom nomP;
int ncvals;
Ligne *Base1,*Base2;

mpxInvite.VuInvite( "Selectionnez la premiere ligne ou tapez son nom : ");
if ( (Base1 = LesLignes->Lignes::CliquerLigne()) == NULL )
return;
mpxInvite.VuEffacerMemoirePx();

mpxInvite.VuInvite( "Selectionnez la seconde ligne ou tapez son nom : ");
if ( (Base2 = LesLignes->Lignes::CliquerLigne()) == NULL )
return;
mpxInvite.VuEffacerMemoirePx();

degreu = 1; //les deux courbes sont ici tendues
nbnu = 1;
degrev = 1;
nbnv = 1;

PC = new Point **[1+nbnu];
for (int i=0; i<=nbnu; i++)
PC[i] = new Point *[1+nbnv];

PC[0][0]=Base1->PC[0];
PC[0][1]=Base1->PC[1];
PC[1][0]=Base2->PC[0];
PC[1][1]=Base2->PC[1];

U = new R [1+nbnu+1+degreu]; //pas de noeud intermediaire
U[0] = U[1] = 0;
U[2] = U[3] = 1;

V = new R [1+nbnv+1+degrev];
V[0] = V[1] = 0;
V[2] = V[3] = 1;

Poids = new R *[1+nbnu];
for (int i=0; i<=nbnu; i++)
Poids[i] = new R[1+nbnv];

for (int i=0; i<=nbnu; i++)
Poids[i][0] = Poids[i][1] = 1; //les poids sont egaux la valeur du noeud ne compte pas ici

Courbes = new Ligne *[2];
Courbes[0] = Base1;
Courbes[1] = Base2;

Mailler( LesPoints );

}

void SurfaceBilineaire::CalculNbEF( int & nbefu, int & nbefv )
{
R long1 = Courbes[0]->LongueurNURBS();
R long2 = Norme2( (Courbes[1]->PC[0]->XyzSommet[0]) - (Courbes[0]->PC[0]->XyzSommet[0]) );

nbefu = (int) ( long1 / TailleArete + 0.5 );
nbefv = (int) ( long2 / TailleArete + 0.5 );
}

void SurfaceBilineaire::AfficherDefinition()
{
cout << "Nom de la Surface BILINEAIRE= " << nom <<endl;
}

void SurfaceBilineaire::VerifierMaillage(Points *LesPoints)
{
if( (Courbes[0]->date > this->date) || (Courbes[1]->date > this->date) )
Mailler(LesPoints);
}

//=============================================================================
// Classe SurfaceExtrudee
//=============================================================================

SurfaceExtrudee::SurfaceExtrudee(Lignes *LesLignes,Points *LesPoints):Surface()
{
Nom nomP;
Ligne * Base;
int ncvals;

mpxInvite.VuInvite( "Selectionnez la courbe à extruder ou tapez son nom : ");
if ( (Base = LesLignes->Lignes::CliquerLigne()) == NULL )
return;
mpxInvite.VuEffacerMemoirePx();

mpxInvite.VuInvite( "Rentrez la hauteur d'extrusion : ");
mpxInvite.lireR( ncvals, Hauteur );
mpxInvite.VuEffacerMemoirePx();

degreu = Base->degre;
nbnu = Base->nbn;
degrev = 1;
nbnv = 1;

R3 tmpSCG;
R3 tmpSCU(0,0,Hauteur);
VecteurScuScg( SCU, tmpSCU, tmpSCG ); //extrude verticallement dans le SCU

PC = new Point **[1+nbnu];
for (int i=0; i<=nbnu; i++)
PC[i] = new Point *[1+nbnv];


for (int i=0; i<=nbnu; i++)
{
NommerPoint(nomP);
PC[i][0] = new Pointxyz( nomP,
Base->PC[i]->XyzSommet[0].x,
Base->PC[i]->XyzSommet[0].y,
Base->PC[i]->XyzSommet[0].z );
LesPoints->Ajouter1Point( PC[i][0] );

NommerPoint( nomP );
PC[i][1] = new Pointxyz( nomP,
Base->PC[i]->XyzSommet[0].x + tmpSCG.x,
Base->PC[i]->XyzSommet[0].y + tmpSCG.y,
Base->PC[i]->XyzSommet[0].z + tmpSCG.z ) ;
LesPoints->Ajouter1Point( PC[i][1] );
}

U = new R [1+nbnu+1+degreu]; //ceci est la courbe libre a extruder
for (int i=0; i<=1+nbnu+degreu; i++)
U[i] = Base->U[i];

V = new R [1+nbnv+1+degrev]; //la deuxieme courbe est un chemin elementaire
V[0] = V[1] = 0;
V[2] = V[3] = 1;

if ( Base->Poids )
{
Poids = new R *[1+nbnu];
for (int i=0; i<=nbnu; i++)
Poids[i] = new R[1+nbnv];

for (int i=0; i<=nbnu; i++)
Poids[i][0] = Poids[i][1] = Base->Poids[i];
}

Courbes = new Ligne *[1];
Courbes[0] = Base;
Mailler( LesPoints );

}

void SurfaceExtrudee::CalculNbEF( int & nbefu, int & nbefv )
{
nbefu = (int) ( Courbes[0]->LongueurNURBS() / TailleArete + 0.5 );
nbefv = 1 + (int) fabs( (Hauteur / TailleArete + 0.5) );
}

void SurfaceExtrudee::AfficherDefinition()
{
cout << "Nom de la Surface EXTRUDEE= " << nom
<< " Hauteur de la Surface EXTRUDEE= " << Hauteur <<endl;
}

void SurfaceExtrudee::VerifierMaillage(Points *LesPoints)
{
if( (Courbes[0]->date > this->date) ) Mailler(LesPoints);
}

//=============================================================================
// Classe SurfaceRevolution
//=============================================================================

SurfaceRevolution::SurfaceRevolution(Lignes *LesLignes,Points *LesPoints):Surface()
{

Nom nomP;
Ligne *Genere;
Point *PtI,*PtII;
int ncvals;

mpxInvite.VuInvite( "Selectionnez la courbe generatrice ou tapez son nom : ");
if ( (Genere = LesLignes->Lignes::CliquerLigne()) == NULL )
return;
mpxInvite.VuEffacerMemoirePx();

mpxInvite.VuInvite( "Premier point definissant l'axe de revolution : ");
if ( (PtI = LesPoints->Points::CliquerPoint()) == NULL )
return;
mpxInvite.VuEffacerMemoirePx();

PtII = PtI;
do
{
mpxInvite.VuInvite( "Second point definissant l'axe de revolution : ");
if ( (PtII = LesPoints->Points::CliquerPoint()) == NULL )
return;
mpxInvite.VuEffacerMemoirePx();
}
while ( PtII == PtI );

degreu = Genere->degre;
nbnu = Genere->nbn;
degrev = 2;
nbnv = 8;

R3 A,B,P,Q,V1,V2,W1,W2,Aux; //Quelques points et vecteurs auxiliaires
R3 Liste [nbnv+1]; //Une liste auxiliaire de R3
R lambda; //Une constante auxiliaire

//Calculs hors toute boucle

A = PtI->XyzSommet[0];
B = PtII->XyzSommet[0];
V2 = B - A;
Rayon = 0;

PC = new Point **[1+nbnu];
for (int i=0; i<=nbnu; i++)
PC[i] = new Point *[1+nbnv];

for (int i=0; i<=nbnu; i++)
{
//calculs intermediaires
P = Genere->PC[i]->XyzSommet[0];
V1 = P - A; //un R3 est manipulable indifferemment comme un point ou un vecteur
lambda = (V1,V2)/(V2,V2);
Q = A + lambda*V2;

W1 = P - Q;
if ((W1,W1)) //Point hors l'axe
{
Aux = W1^V2;
W2 = (Aux/Norme2(Aux))*Norme2(W1);
Rayon = Max(Rayon,Norme2(W1));

//Les 8 points du carre de controle du cercle

Liste[0] = Q+W1;
Liste[1] = Q+W1+W2;
Liste[2] = Q+W2;
Liste[3] = Q-W1+W2;
Liste[4] = Q-W1;
Liste[5] = Q-W1-W2;
Liste[6] = Q-W2;
Liste[7] = Q+W1-W2;
Liste[8] = Liste[0];
}
else //Point sur l'axe=>tous identiques
{
for (int k=0; k<=nbnv; k++)
Liste[k]=Q;
}

for (int j=0; j<=nbnv; j++)
{
NommerPoint(nomP);
PC[i][j] = new Pointxyz (nomP,Liste[j].x,Liste[j].y,Liste[j].z);
LesPoints->Ajouter1Point(PC[i][j]);
}
}

U = new R [1+nbnu+1+degreu];
for (int i=0; i<=1+nbnu+degreu; i++)
U[i] = Genere->U[i];

V = new R [1+nbnv+1+degrev]; //valeur des noeuds du cercle est un resultat connu
V[0] = V[1] = V[2] = 0.0;
V[3] = V[4] = 0.25;
V[5] = V[6] = 0.5;
V[7] = V[8] = 0.75;
V[9] = V[10] = V[11] = 1.0;

//Existence des poids donc creation dynamique d un tableau poitant sur les poids

Poids = new R *[1+nbnu];
for (int i=0; i<=nbnu; i++)
Poids[i] = new R[1+nbnv];

if (Genere->Poids)
{
for (int i=0; i<=nbnu; i++)
for (int j=0; j<=nbnv; j++)
Poids[i][j]=((j%2)*RAC2_2+(!(j%2)))*Genere->Poids[i]; //ici les valeurs des poids suivent une loi modulo
}
else
{
for (int i=0; i<=nbnu; i++)
for (int j=0; j<=nbnv; j++)
Poids[i][j]=(j%2)*RAC2_2+(!(j%2));
}

Courbes = new Ligne *[1];
Courbes[0] = Genere;
Mailler( LesPoints );

}

void SurfaceRevolution::CalculNbEF( int & nbefu, int & nbefv )
{
nbefu= (int) ( Courbes[0]->LongueurNURBS() / TailleArete + 0.5 );
nbefv=(int)( 1 + fabs( ( (2*PI*Rayon) / TailleArete + 0.5) ) );
}

void SurfaceRevolution::AfficherDefinition()
{
cout << "Nom de la surface de REVOLUTION= " << nom <<endl;
}

void SurfaceRevolution::VerifierMaillage(Points *LesPoints)
{
if( (Courbes[0]->date > this->date) )
Mailler(LesPoints);
}

//=============================================================================
// Classe SurfaceGauche
//=============================================================================

SurfaceGauche::SurfaceGauche(Lignes *LesLignes,Points *LesPoints):Surface()
{
Nom nomP;
Ligne *C1, *C2, *D1, *D2;
int ncvals;

mpxInvite.VuInvite( "Selectionnez la premiere courbe ou tapez son nom : "); //interface utilisateur habituelle
if ( (C1 = LesLignes->Lignes::CliquerLigne()) == NULL )
return;
mpxInvite.VuEffacerMemoirePx();

mpxInvite.VuInvite( "Selectionnez sa courbe opposée ou tapez son nom : ");
if ( (C2 = LesLignes->Lignes::CliquerLigne()) == NULL )
return;
mpxInvite.VuEffacerMemoirePx();

mpxInvite.VuInvite( "Selectionnez la deuxieme courbe ou tapez son nom : ");
if ( (D1 = LesLignes->Lignes::CliquerLigne()) == NULL )
return;
mpxInvite.VuEffacerMemoirePx();

mpxInvite.VuInvite( "Selectionnez sa courbe opposee ou tapez son nom : ");
if ( (D2 = LesLignes->Lignes::CliquerLigne()) == NULL )
return;
mpxInvite.VuEffacerMemoirePx();

if ( (C1->nbn != C2->nbn) || (D1->nbn != D2->nbn) ||
(C1->degre != C2->degre) || (D1->degre != D2->degre) )
return;

degreu = C1->degre; //4 courbes 2 a 2 en face les unes des autres
nbnu = C1->nbn;
degrev = D1->degre;
nbnv = D1->nbn;

R3 Aux;

PC = new Point **[1+nbnu];
for (int i=0; i<=nbnu; i++)
PC[i] = new Point *[1+nbnv];

for (int i=0; i<=nbnu; i++)
{
NommerPoint(nomP);
PC[i][0] = new Pointxyz( nomP,
C1->PC[i]->XyzSommet[0].x,
C1->PC[i]->XyzSommet[0].y,
C1->PC[i]->XyzSommet[0].z );
LesPoints->Ajouter1Point( PC[i][0] );

NommerPoint( nomP );
PC[i][nbnv] = new Pointxyz( nomP,
C2->PC[i]->XyzSommet[0].x,
C2->PC[i]->XyzSommet[0].y,
C2->PC[i]->XyzSommet[0].z ) ;
LesPoints->Ajouter1Point( PC[i][nbnv] );
}

for (int j=0; j<=nbnv; j++)
{
NommerPoint(nomP);
PC[0][j] = new Pointxyz( nomP,
D1->PC[j]->XyzSommet[0].x,
D1->PC[j]->XyzSommet[0].y,
D1->PC[j]->XyzSommet[0].z );
LesPoints->Ajouter1Point( PC[0][j] );

NommerPoint( nomP );
PC[nbnu][j] = new Pointxyz( nomP,
D2->PC[j]->XyzSommet[0].x,
D2->PC[j]->XyzSommet[0].y,
D2->PC[j]->XyzSommet[0].z ) ;
LesPoints->Ajouter1Point( PC[nbnu][j] );
}

R mu1, mu2; //Coefficients

for (int i=1; i<nbnu; i++)
{
mu1 = i/nbnu; //Coefficient multiplicatif des points de C1 et C2
for (int j=1; j<nbnv; j++)
{

mu2 = j/nbnv;//Coefficient multiplicatif des points de D1 et D2

Aux.x=((1-mu1)*(C1->PC[i]->XyzSommet[0].x)+
(mu1)*(C2->PC[i]->XyzSommet[0].x)
+(1-mu2)*(D1->PC[j]->XyzSommet[0].x)+
(mu2)*(D2->PC[j]->XyzSommet[0].x))/2; //c est la formule du cours a prpos ds surfaces reglees
Aux.y=((1-mu1)*(C1->PC[i]->XyzSommet[0].y)+
(mu1)*(C2->PC[i]->XyzSommet[0].y)
+(1-mu2)*(D1->PC[j]->XyzSommet[0].y)+
+(mu2)*(D2->PC[j]->XyzSommet[0].y))/2;
Aux.z=((1-mu1)*(C1->PC[i]->XyzSommet[0].z)+
(mu1)*(C2->PC[i]->XyzSommet[0].z)
+(1-mu2)*(D1->PC[j]->XyzSommet[0].z)+
(mu2)*(D2->PC[j]->XyzSommet[0].z))/2;
NommerPoint( nomP );
PC[i][j] = new Pointxyz( nomP, Aux.x, Aux.y, Aux.z ) ;
LesPoints->Ajouter1Point( PC[i][j] );
}
}

//on unie les jeux de liens qui sont en vis a vis

U = new R [1+nbnu+1+degreu];
for (int i=0; i<=1+nbnu+degreu; i++)
U[i] = C1->U[i];

V = new R [1+nbnv+1+degrev];
for (int j=0; j<=1+nbnv+degrev; j++)
V[j] = D1->U[j];

if ( C1->Poids && C2->Poids && D1->Poids && D2->Poids)
{
Poids = new R *[1+nbnu];
for (int i=0; i<=nbnu; i++)
Poids[i] = new R[1+nbnv];

for (int i=0; i<=nbnu; i++)
for (int j=0; j<=nbnv; j++)
Poids[i][j] =
(C1->Poids[i]+C2->Poids[i]+D1->Poids[j]+D2->Poids[j])/4;
}

Courbes = new Ligne *[4];
Courbes[0] = C1;
Courbes[1] = C2;
Courbes[2] = D1;
Courbes[3] = D2;

Mailler( LesPoints );

}

void SurfaceGauche::CalculNbEF( int & nbefu, int & nbefv )
{
nbefu = (int) ( Courbes[0]->LongueurNURBS() / TailleArete + 0.5 );
nbefv = (int) ( Courbes[2]->LongueurNURBS() / TailleArete + 0.5 );
}

void SurfaceGauche::AfficherDefinition()
{
cout << "Nom de la Surface GAUCHE= " << nom <<endl;
}

void SurfaceGauche::VerifierMaillage(Points *LesPoints)
{
if( (Courbes[0]->date > this->date) ) Mailler(LesPoints);
}

//=============================================================================
// Classe SurfaceTore
//=============================================================================

SurfaceTore::SurfaceTore(Lignes *LesLignes,Points *LesPoints):Surface()
{

Nom nomP;
Ligne *Genere;
Point *PtI,*PtII;
int ncvals;

do
{
mpxInvite.VuInvite( "Selectionnez le cercle ou tapez son nom : "); //interface utilisateur
if ( (Genere = LesLignes->Lignes::CliquerLigne()) == NULL )
return;
mpxInvite.VuEffacerMemoirePx();
}
while ( (Genere->nbn!=8) || (Genere->degre!=2) );

mpxInvite.VuInvite( "Premier point definissant l'axe de revolution : ");
if ( (PtI = LesPoints->Points::CliquerPoint()) == NULL )
return;
mpxInvite.VuEffacerMemoirePx();

PtII = PtI;
do
{
mpxInvite.VuInvite( "Second point definissant l'axe de revolution : ");
if ( (PtII = LesPoints->Points::CliquerPoint()) == NULL )
return;
mpxInvite.VuEffacerMemoirePx();
}
while ( PtII == PtI );

degreu = Genere->degre;
nbnu = Genere->nbn;
degrev = 2;
nbnv = 8; //car ici on aurait pu faire deriver cette classe de la classe des surfaces de revolution

R3 A,B,P,Q,V1,V2,W1,W2,Aux; //variables annexes
R3 Liste [nbnv+1]; //points de controle generes pour un point de la premiere courbe nurbs
R lambda;

A = PtI->XyzSommet[0];
B = PtII->XyzSommet[0];
V2 = B - A;
Rayon = 0;

//Une boucle exterieure pour calculer le carre de controle
//pour chaque point de controle de la courbe generatrice

PC = new Point **[1+nbnu];
for (int i=0; i<=nbnu; i++)
PC[i] = new Point *[1+nbnv];

for (int i=0; i<=nbnu; i++)
{
//Quelques points et vecteurs necessaires
P = Genere->PC[i]->XyzSommet[0];
V1 = P - A;
lambda = (V1,V2)/(V2,V2);
Q = A + lambda*V2;

W1 = P - Q;
if ((W1,W1)) //Point hors l'axe
{
Aux = W1^V2;
W2 = (Aux/Norme2(Aux))*Norme2(W1);
Rayon = Max(Rayon,Norme2(W1));

//Les 8 points du carre de controle du cercle

Liste[0] = Q+W1;
Liste[1] = Q+W1+W2;
Liste[2] = Q+W2;
Liste[3] = Q-W1+W2;
Liste[4] = Q-W1;
Liste[5] = Q-W1-W2;
Liste[6] = Q-W2;
Liste[7] = Q+W1-W2;
Liste[8] = Liste[0];
}
else //Point sur l'axe
{
for (int k=0; k<=nbnv; k++)
Liste[k]=Q;
}

for (int j=0; j<=nbnv; j++)
{
NommerPoint(nomP);
PC[i][j] = new Pointxyz (nomP,Liste[j].x,Liste[j].y,Liste[j].z);
LesPoints->Ajouter1Point(PC[i][j]);
}
}

U = new R [1+nbnu+1+degreu];
for (int i=0; i<=1+nbnu+degreu; i++)
U[i] = Genere->U[i];

V = new R [1+nbnv+1+degrev];
V[0] = V[1] = V[2] = 0.0;
V[3] = V[4] = 0.25;
V[5] = V[6] = 0.5;
V[7] = V[8] = 0.75;
V[9] = V[10] = V[11] = 1.0;

//Existence donc creation des poids

Poids = new R *[1+nbnu];
for (int i=0; i<=nbnu; i++)
Poids[i] = new R[1+nbnv];

if (Genere->Poids)
{
for (int i=0; i<=nbnu; i++)
for (int j=0; j<=nbnv; j++)
Poids[i][j]=((j%2)*RAC2_2+(!(j%2)))*Genere->Poids[i]; //regle en modulo
}
else
{
for (int i=0; i<=nbnu; i++)
for (int j=0; j<=nbnv; j++)
Poids[i][j]=(j%2)*RAC2_2+(!(j%2));
}

Courbes = new Ligne *[1];
Courbes[0] = Genere;
Mailler( LesPoints );
}

void SurfaceTore::CalculNbEF( int & nbefu, int & nbefv )
{
nbefu= (int) ( Courbes[0]->LongueurNURBS() / TailleArete + 0.5 );
nbefv=(int) (2*PI*nbefu);
}

void SurfaceTore::AfficherDefinition()
{
cout << "Nom de la Surface TORE= " << nom <<endl;
}
void SurfaceTore::VerifierMaillage(Points *LesPoints)
{
if( (Courbes[0]->date > this->date) )
Mailler(LesPoints);
}

//=============================================================================
// Classe SurfaceCylindre
//=============================================================================

SurfaceCylindre::SurfaceCylindre(Lignes * LesLignes, Points *LesPoints):Surface()
{
Nom nomP;
Ligne * BaseC;
int ncvals;

BaseC = new LigneCercle(LesPoints);
if( BaseC->NbSommet == 0 )
return ;
LesLignes->Ajouter1Ligne( BaseC );
mpxInvite.VuEffacerMemoirePx();

mpxInvite.VuInvite( "Rentrez la hauteur du CYLINDRE : ");
mpxInvite.lireR( ncvals, Hauteur );
mpxInvite.VuEffacerMemoirePx();

degreu = BaseC->degre;
nbnu = BaseC->nbn;
degrev = 1;
nbnv = 1;

R3 tmpSCG;
R3 tmpSCU(0,0,Hauteur);
VecteurScuScg( SCU, tmpSCU, tmpSCG );

PC = new Point **[1+nbnu];
for (int i=0; i<=nbnu; i++)
PC[i] = new Point *[1+nbnv];


for (int i=0; i<=nbnu; i++)
{
NommerPoint(nomP);
PC[i][0] = new Pointxyz( nomP,
BaseC->PC[i]->XyzSommet[0].x,
BaseC->PC[i]->XyzSommet[0].y,
BaseC->PC[i]->XyzSommet[0].z );
LesPoints->Ajouter1Point( PC[i][0] );

NommerPoint( nomP );
PC[i][1] = new Pointxyz( nomP,
BaseC->PC[i]->XyzSommet[0].x + tmpSCG.x,
BaseC->PC[i]->XyzSommet[0].y + tmpSCG.y,
BaseC->PC[i]->XyzSommet[0].z + tmpSCG.z ) ;
LesPoints->Ajouter1Point( PC[i][1] );
}

U = new R [1+nbnu+1+degreu];
for (int i=0; i<=1+nbnu+degreu; i++)
U[i] = BaseC->U[i];

V = new R [1+nbnv+1+degrev];
V[0] = V[1] = 0;
V[2] = V[3] = 1;

if ( BaseC->Poids )
{
Poids = new R *[1+nbnu];
for (int i=0; i<=nbnu; i++)
Poids[i] = new R[1+nbnv];

for (int i=0; i<=nbnu; i++)
Poids[i][0] = Poids[i][1] = BaseC->Poids[i];
}

Courbes = new Ligne *[1];
Courbes[0] = BaseC;
Mailler( LesPoints );

}

void SurfaceCylindre::CalculNbEF( int & nbefu, int & nbefv )
{
nbefu = (int) ( Courbes[0]->LongueurNURBS() / TailleArete + 0.5 );
nbefv = 1 + (int) fabs( (Hauteur / TailleArete + 0.5) );
}

void SurfaceCylindre::AfficherDefinition()
{
cout << "Nom de la Surface CYLINDRE= " << nom
<< " Hauteur du CYLINDRE= " << Hauteur <<endl;
}

void SurfaceCylindre::VerifierMaillage(Points *LesPoints)
{
if( (Courbes[0]->date > this->date) ) Mailler(LesPoints);
}

//=============================================================================
// Classe SurfaceCone
//=============================================================================

SurfaceCone::SurfaceCone(Lignes *LesLignes,Points *LesPoints):Surface()
{

Nom nomP;
Ligne *Genere;
Point *PtI,*PtII;
int ncvals;

do
{
mpxInvite.VuInvite( "Selectionnez la droite generatrice ou tapez son nom : ");
if ( (Genere = LesLignes->Lignes::CliquerLigne()) == NULL )
return;
mpxInvite.VuEffacerMemoirePx();
}
while ( (Genere->nbn!=1) || (Genere->degre!=1) );

mpxInvite.VuInvite( "Premier point definissant l'axe de revolution : ");
if ( (PtI = LesPoints->Points::CliquerPoint()) == NULL )
return;
mpxInvite.VuEffacerMemoirePx();

PtII = PtI;
do
{
mpxInvite.VuInvite( "Second point definissant l'axe de revolution : ");
if ( (PtII = LesPoints->Points::CliquerPoint()) == NULL )
return;
mpxInvite.VuEffacerMemoirePx();
}
while ( PtII == PtI );

degreu = Genere->degre;
nbnu = Genere->nbn;
degrev = 2;
nbnv = 8;

R3 A,B,P,Q,V1,V2,W1,W2,Aux;
R3 Liste [nbnv+1];
R lambda;

A = PtI->XyzSommet[0];
B = PtII->XyzSommet[0];
V2 = B - A;

PC = new Point **[1+nbnu];
for (int i=0; i<=nbnu; i++)
PC[i] = new Point *[1+nbnv];

for (int i=0; i<=nbnu; i++)
{
//calculs intermediaires
P = Genere->PC[i]->XyzSommet[0];
V1 = P - A;
lambda = (V1,V2)/(V2,V2);
Q = A + lambda*V2;

W1 = P - Q;
if ((W1,W1)) //Point hors axe
{
Aux = W1^V2;
W2 = (Aux/Norme2(Aux))*Norme2(W1);

//Les 8 points du carre de controle du cercle

Liste[0] = Q+W1;
Liste[1] = Q+W1+W2;
Liste[2] = Q+W2;
Liste[3] = Q-W1+W2;
Liste[4] = Q-W1;
Liste[5] = Q-W1-W2;
Liste[6] = Q-W2;
Liste[7] = Q+W1-W2;
Liste[8] = Liste[0];
}
else
{
for (int k=0; k<=nbnv; k++)
Liste[k]=Q;
}

for (int j=0; j<=nbnv; j++)
{
NommerPoint(nomP);
PC[i][j] = new Pointxyz (nomP,Liste[j].x,Liste[j].y,Liste[j].z);
LesPoints->Ajouter1Point(PC[i][j]);
}
}

U = new R [1+nbnu+1+degreu];
for (int i=0; i<=1+nbnu+degreu; i++)
U[i] = Genere->U[i];

V = new R [1+nbnv+1+degrev];
V[0] = V[1] = V[2] = 0.0;
V[3] = V[4] = 0.25;
V[5] = V[6] = 0.5;
V[7] = V[8] = 0.75;
V[9] = V[10] = V[11] = 1.0;

//Existence donc creation des poids

Poids = new R *[1+nbnu];
for (int i=0; i<=nbnu; i++)
Poids[i] = new R[1+nbnv];

if (Genere->Poids)
{
for (int i=0; i<=nbnu; i++)
for (int j=0; j<=nbnv; j++)
Poids[i][j]=((j%2)*RAC2_2+(!(j%2)))*Genere->Poids[i]; //regle du modulo pour les poids du cercle
}
else
{
for (int i=0; i<=nbnu; i++)
for (int j=0; j<=nbnv; j++)
Poids[i][j]=(j%2)*RAC2_2+(!(j%2));
}

Courbes = new Ligne *[1];
Courbes[0] = Genere;
Mailler( LesPoints );
}

void SurfaceCone::CalculNbEF( int & nbefu, int & nbefv )
{
nbefu= (int) ( Courbes[0]->LongueurNURBS() / TailleArete + 0.5 );
nbefv= nbefu;
}

void SurfaceCone::AfficherDefinition()
{
cout << "Nom de la Surface CONE= " << nom <<endl;
}
void SurfaceCone::VerifierMaillage(Points *LesPoints)
{
if( (Courbes[0]->date > this->date) )
Mailler(LesPoints);
}

//=============================================================================
// Classe Sphere
//=============================================================================

SurfaceSphere::SurfaceSphere(Lignes *LesLignes,Points *LesPoints):Surface()
{
Nom nomP,nomL,nomC;
Ligne *Genere;
Point *PtI,*PtII;
int ncvals;

mpxInvite.VuInvite( "Centre de la SPHERE : "); //interface utilisateur
if ( (PtCentre = LesPoints->Points::CliquerPoint()) == NULL )
return;
mpxInvite.VuEffacerMemoirePx();
mpxInvite.VuInvite( "Rentrez le rayon de la SPHERE : ");
mpxInvite.lireR( ncvals, Rayon );
mpxInvite.VuEffacerMemoirePx();

R aux[6]; //noeuds du demi-cercle
aux[0]=aux[1]=aux[2]=0;
aux[3]=aux[4]=0.5;
aux[5]=1;
Nom nompc[5]; //Noms des points du demi-cercle
R poidscercle[5];//Poids du demi-cercle=poids du cercle entier
Point ** PP;
PP= new Point *[5]; //Points de controle du demi-cercle=points de controle du cercle

NommerPoint(nomP);
PP[0]=new Pointxyz(nomP,PtCentre->XyzSommet[0].x-Rayon,
PtCentre->XyzSommet[0].y,
PtCentre->XyzSommet[0].z);
LesPoints->Ajouter1Point(PP[0]);

NommerPoint(nomP);
PP[1]=new Pointxyz(nomP,PtCentre->XyzSommet[0].x-Rayon,
PtCentre->XyzSommet[0].y+Rayon,
PtCentre->XyzSommet[0].z);
LesPoints->Ajouter1Point(PP[1]);

NommerPoint(nomP);
PP[2]=new Pointxyz(nomP,PtCentre->XyzSommet[0].x,
PtCentre->XyzSommet[0].y+Rayon,
PtCentre->XyzSommet[0].z);
LesPoints->Ajouter1Point(PP[2]);

NommerPoint(nomP);
PP[3]=new Pointxyz(nomP,PtCentre->XyzSommet[0].x+Rayon,
PtCentre->XyzSommet[0].y+Rayon,
PtCentre->XyzSommet[0].z);
LesPoints->Ajouter1Point(PP[3]);

NommerPoint(nomP);
PP[4]=new Pointxyz(nomP,PtCentre->XyzSommet[0].x+Rayon,
PtCentre->XyzSommet[0].y,
PtCentre->XyzSommet[0].z);
LesPoints->Ajouter1Point(PP[4]);

for (int k=0; k<=4; k++)
{
nompc[k]=PP[k]->nom;
poidscercle[k]=(k%2)*RAC2_2+(1-(k%2));
}

NommerLigne(nomL);
Genere = new LigneNURBS(LesPoints,nomL,2,5,aux,nompc,poidscercle);
LesLignes->Ajouter1Ligne(Genere);

NommerPoint(nomP);
PtI=new Pointxyz(nomP,PtCentre->XyzSommet[0].x-Rayon,
PtCentre->XyzSommet[0].y,
PtCentre->XyzSommet[0].z);
LesPoints->Ajouter1Point(PtI);

NommerPoint(nomP);
PtII=new Pointxyz(nomP,PtCentre->XyzSommet[0].x+Rayon,
PtCentre->XyzSommet[0].y,
PtCentre->XyzSommet[0].z);
LesPoints->Ajouter1Point(PtII);

degreu = 2;
nbnu = 4;
degrev = 2;
nbnv = 8;

R3 A,B,P,Q,V1,V2,W1,W2,Aux;
R3 Liste [nbnv+1];
R lambda;

A = PtI->XyzSommet[0];
B = PtII->XyzSommet[0];
V2 = B - A;
Rayon = 0;

PC = new Point **[1+nbnu];
for (int i=0; i<=nbnu; i++)
PC[i] = new Point *[1+nbnv];

for (int i=0; i<=nbnu; i++)
{
//calculs intermediaires
P = Genere->PC[i]->XyzSommet[0];
V1 = P - A;
lambda = (V1,V2)/(V2,V2);
Q = A + lambda*V2;
W1 = P - Q;

if ((W1,W1)) //Point hors axe
{
Aux = W1^V2;
W2 = (Aux/Norme2(Aux))*Norme2(W1);
Rayon = Max(Rayon,Norme2(W1));
//Les 8 points du carre de controle du cercle

Liste[0] = Q+W1;
Liste[1] = Q+W1+W2;
Liste[2] = Q+W2;
Liste[3] = Q-W1+W2;
Liste[4] = Q-W1;
Liste[5] = Q-W1-W2;
Liste[6] = Q-W2;
Liste[7] = Q+W1-W2;
Liste[8] = Liste[0];
}
else
{
for (int k=0; k<=nbnv; k++)
Liste[k]=Q;
}

for (int j=0; j<=nbnv; j++)
{
NommerPoint(nomP);
PC[i][j] = new Pointxyz (nomP,Liste[j].x,Liste[j].y,Liste[j].z);
LesPoints->Ajouter1Point(PC[i][j]);
}
}

U = new R [1+nbnu+1+degreu];
for (int i=0; i<=1+nbnu+degreu; i++)
U[i] = Genere->U[i];

V = new R [1+nbnv+1+degrev];
V[0] = V[1] = V[2] = 0.0;
V[3] = V[4] = 0.25;
V[5] = V[6] = 0.5;
V[7] = V[8] = 0.75;
V[9] = V[10] = V[11] = 1.0;

//Existence donc creation des poids

Poids = new R *[1+nbnu];
for (int i=0; i<=nbnu; i++)
Poids[i] = new R[1+nbnv];

if (Genere->Poids)
{
for (int i=0; i<=nbnu; i++)
for (int j=0; j<=nbnv; j++)
Poids[i][j]=((j%2)*RAC2_2+(1-(j%2)))*Genere->Poids[i];
}
else
{
for (int i=0; i<=nbnu; i++)
for (int j=0; j<=nbnv; j++)
Poids[i][j]=(j%2)*RAC2_2+(1-(j%2));
}

Courbes = new Ligne *[1];
Courbes[0] = Genere;
Mailler( LesPoints );

}

void SurfaceSphere::CalculNbEF( int & nbefu, int & nbefv )
{
nbefu= (int) ( Courbes[0]->LongueurNURBS() / TailleArete + 0.5 );
nbefv=(int)( 1 + fabs( ( (2*PI*Rayon) / TailleArete + 0.5) ) );
}

void SurfaceSphere::AfficherDefinition()
{
cout << "Nom de la Surface SPHERE= " << nom
<< " Nom du Point CENTRE= " << PtCentre->nom
<< " Rayon de la SPHERE= " << Rayon <<endl;
}
void SurfaceSphere::VerifierMaillage(Points *LesPoints)
{
if( (Courbes[0]->date > this->date) )
Mailler(LesPoints);
}

//============================================================================
// Classe Surfaces
//============================================================================

Surfaces::Surfaces()
{
P1Surface=NULL;
}

Surfaces::~Surfaces()
{
#ifdef VERBOSE
cout <<"Entrée dans Surfaces" <<endl;
#endif
if (Surfaces::P1Surface != NULL )
delete P1Surface;
#ifdef VERBOSE
cout << "Sortie de ~Surfaces" << endl;
#endif
}

Surface * Surfaces::RetrouverSurface( Nom NomSf, Surface * & SfAvant ) // retourne NULL si pas retrouve
{
Surface * SfCourant;

SfAvant=NULL;
SfCourant = P1Surface;
while ( SfCourant != NULL )
{
if ( strcmp( SfCourant->nom, NomSf )==0 )
return SfCourant;
else
{
SfAvant=SfCourant;
SfCourant=SfCourant->PSuivant;
}
}
return NULL;
}

Surface * Surfaces::CliquerSurface( void )
// retourne le pointeur sur la surface la plus proche de la souris
// NULL si la surface pas retrouvee ou abandon demande
{
int notyev, nbc;
XPoint pxy, pxymin0={-1,-1}, pxymin;
Surface *SfCourant, *SfMin;
R3 xyz, xyzm, axyz;
R d, dmin;
Nom nomSurface;

//saisie a la souris de la poignee de la Surface
SaisieT:
mpxInvite.VuSourisTexte( notyev, nbc, pxy, nomSurface );
if( notyev == 0 ) { return NULL; } //Abandon demande

if( notyev < 0 ) //souris deplacee ou bouton enfonce non relache
{
//PB ici en 3d
xyz.x = mpxPLSV.Vuxob2px( pxy.x ); //construction du point de R2 correspondant
xyz.y = mpxPLSV.Vuyob2px( pxy.y );
xyz.z = 0;

//recherche de la Surface actuelle la plus proche du clic souris
//============================================================
dmin = 1e100;
SfMin = NULL;
SfCourant = P1Surface;
while ( SfCourant != NULL )
{
//PB PML pour ligne et non pour surface
xyzm = SfCourant->PointMilieuSurface(); //la poignee de la ligne
if( mpxPLSV.visee == 2 ) //Visee 2d
{
d = Norme2_2( xyz - SfCourant->PointMilieuSurface() );
}
else
{ //Visee 3d Un point pxy clique est une droite dans R3
//recherche du point le plus proche de cette droite
// Passage en coordonnees axonometriques du point
mpxPLSV.Vuxyzaxo( xyzm, axyz );
//distance en XY axonometriques seulement
d = sqrt( pow( xyz.x-axyz.x, 2 )
+ pow( xyz.y-axyz.y, 2 ) );
}
if ( d < dmin ) { dmin = d; SfMin = SfCourant; }
SfCourant = SfCourant->PSuivant;
}

//Verification que SfMin n'est pas trop loin du point clique
if( SfMin != NULL )
{
//PB PML pour ligne et non pour surface
xyzm = SfMin->PointMilieuSurface();
if( mpxPLSV.visee == 2 ) //Visee 2d
{
pxymin.x = mpxPLSV.Vunupxex( xyzm.x );
pxymin.y = mpxPLSV.Vunupxey( xyzm.y );
if( abs( pxymin.x - pxy.x ) > 123
|| abs( pxymin.y - pxy.y ) > 123 ) {SfMin = NULL;}
}
else
{ //Visee 3d
mpxPLSV.Vuxyzaxo( xyzm, axyz );
pxymin.x = mpxPLSV.Vunupxex( axyz.x );
pxymin.y = mpxPLSV.Vunupxey( axyz.y );
if( abs( pxymin.x - pxy.x ) > 123
|| abs( pxymin.y - pxy.y ) > 123 )
{SfMin = NULL;}
}
if ( SfMin == NULL ) //Surface min trop loin => Accrochage a effacer
{
if( (pxymin0.x!=-1) && (pxymin0.y!=-1) )
{
//suppression du mpxAccroche de fenetre_vu
//GXxor 2 fois de suite redonne la fenetre_vu initiale
mpxAccroche.VuCombinerMemoirePxDansFenetre( GXxor,
pxy0, mpxAccroche.largeurpx, mpxAccroche.hauteurpx,
pxymin0 );
pxymin0.x = -1;
pxymin0.y = -1; //reinitialisation a vide
}
}
}
if( SfMin == NULL ) { goto SaisieT; }

//Une Surface est suffisamment proche => son accrochage est materialise
//si ce n'est deja fait
pxymin.x -= (int) mpxAccroche.largeurpx/2;
pxymin.y -= (int) mpxAccroche.hauteurpx/2;
if( (pxymin.x!=pxymin0.x) || (pxymin.y!=pxymin0.y) )
{
if( (pxymin0.x!=-1) && (pxymin0.y!=-1) )
{
//suppression du mpxAccroche de fenetre_vu
//GXxor 2 fois de suite redonne la fenetre_vu initiale
mpxAccroche.VuCombinerMemoirePxDansFenetre( GXxor,
pxy0, mpxAccroche.largeurpx, mpxAccroche.hauteurpx,
pxymin0 );
}
//ajout du mpxAccroche a fenetre_vu
pxymin0 = pxymin;
mpxAccroche.VuCombinerMemoirePxDansFenetre( GXxor,
pxy0, mpxAccroche.largeurpx, mpxAccroche.hauteurpx,
pxymin0 );
}
goto SaisieT;
}
else if( notyev == 1 ) //clic enfonce et relache d'un bouton de la souris
{
if( (pxymin0.x!=-1) && (pxymin0.y!=-1) )
{
//suppression du mpxAccroche de fenetre_vu
//GXxor 2 fois de suite redonne la fenetre_vu initiale
mpxAccroche.VuCombinerMemoirePxDansFenetre( GXxor,
pxy0, mpxAccroche.largeurpx, mpxAccroche.hauteurpx,
pxymin0 );
}
xyz.x = mpxPLSV.Vuxob2px( pxy.x ); //Construction du point de R2 correspondant
xyz.y = mpxPLSV.Vuyob2px( pxy.y );
xyz.z = 0;

//recherche de la Surface actuelle la plus proche du clic souris
dmin = 1e100;
SfMin = NULL;
SfCourant = P1Surface;
while ( SfCourant != NULL )
{
if( mpxPLSV.visee == 2 ) //Visee 2d
{
//PB PML pour ligne et non pour surface
d = Norme2_2( xyz - SfCourant->PointMilieuSurface() );
}
else
{ //Visee 3d Un point pxy clique est une droite dans R3
//recherche du point le plus proche de cette droite
// Passage en coordonnees axonometriques du point
mpxPLSV.Vuxyzaxo( SfCourant->PointMilieuSurface(), axyz );
//distance en XY axonometriques seulement
d = sqrt( pow( xyz.x-axyz.x , 2 )
+ pow( xyz.y-axyz.y , 2 ) );
}
if ( d < dmin ) { dmin = d; SfMin = SfCourant; }
SfCourant = SfCourant->PSuivant;
}
//verification que SfMin n'est pas trop loin du point clique
if( SfMin != NULL )
{
if( mpxPLSV.visee == 2 ) //Visee 2d
{
if( abs( mpxPLSV.Vunupxex( SfMin->PointMilieuSurface().x ) - pxy.x ) > 123
|| abs( mpxPLSV.Vunupxey( SfMin->PointMilieuSurface().y ) - pxy.y ) > 123 )
{SfMin = NULL;}
}
else
{ //Visee 3d
mpxPLSV.Vuxyzaxo( SfMin->PointMilieuSurface(), axyz );
if( abs( mpxPLSV.Vunupxex( axyz.x ) - pxy.x ) > 123
|| abs( mpxPLSV.Vunupxey( axyz.y ) - pxy.y ) > 123 )
{SfMin = NULL;}
}
}
}
else if ( notyev == 2 )
{
Surface *SfAvant;
SfMin = RetrouverSurface( nomSurface, SfAvant ); //par son nom
}

return SfMin;
}

void Surfaces::Ajouter1Surface( Surface *PSurface )
{
Surface * SfAvant;
Surface * SfNom;

SfNom=RetrouverSurface( PSurface->nom, SfAvant );
if( SfNom != NULL )
{ // La Surface etait deja dans la liste => Destruction
if (SfAvant==NULL) // Pas de Surface avant
P1Surface=SfNom->PSuivant; // La suivante de SfNom devient la premiere
else
SfAvant->PSuivant=SfNom->PSuivant;// Le suivant d'avant est le suivant de SfNom

SfNom->PSuivant=NULL;
delete [] SfNom->NoStEF; // La Surface SfNom est detruite
delete [] SfNom->XyzSommet;
}
// La Surface n'est pas dans la liste
PSurface->PSuivant = P1Surface; // La Surface devient la premiere
P1Surface=PSurface; // Surface de la liste
}

void Surfaces::AfficherSurfaces()
{
Surface *PCourant;
cout << endl;
cout << "Affichage des SURFACES:" << endl;
PCourant = P1Surface;
while( PCourant != NULL )
{
PCourant->AfficherDefinition();
PCourant = PCourant->PSuivant;
}
}

void Surfaces::TracerSurfaces( MemoirePx3d & mpx, Points *LesPoints )
{
Surface *PCourant;
PCourant = P1Surface;
while( PCourant != NULL )
{
PCourant->VerifierMaillage( LesPoints );
PCourant->TracerMaillage(mpx);
PCourant = PCourant->PSuivant;
}
}

void Surfaces::TracerPolyedre( MemoirePx3d & mpx )
{
Surface * SfCourant = P1Surface;
while ( SfCourant )
{
SfCourant->TracerPolyedre( mpx );
SfCourant = SfCourant->PSuivant;
}
}

//============================================================================
// Identique a Ligne::IntervalleNURBS, mais avec l'object explicite U

int IntervalleNURBS( R u, R *U, int nbn, int degre) {
int min, max, milieu;

if( u >= U[1+nbn] ) return nbn; //Cas particulier!

//recherche par dichotomie
min = degre;
max = nbn+1;

Iteration:
milieu = (min+max)/2;

if( u < U[milieu] ) {
max = milieu;
goto Iteration;
} else if( u >= U[milieu+1] ) {
min = milieu;
goto Iteration;
} else {
return milieu;
}
}

//============================================================================
// Identique a Ligne::ValeurBaseNURBS, mais avec l'objet explicite U
void ValeurBaseNURBS( int i, Z NbEtape, R u, R U[], R N[] )
{
int j,k;
R g[1+64], d[1+64];//numerateur gauche et droite de N[i,j]
R s, t; //denominateur=numerateur gauche + droite!

N[0] = 1.0;
for (j=1; j<=NbEtape; j++)
{
g[j] = u - U[i+1-j]; //numerateur gauche j
d[j] = U[i+j] - u; //numerateur droite j
s=0.0;
for (k=0; k<j; k++)
{
t = N[k] / ( d[k+1] + g[j-k] ); //la somme des indices vaut j+1
N[k]= s + d[k+1] * t;
s = g[j-k] * t;
}
N[j]=s;
}
}

void Surfaces::TracerPolyedre( MemoirePx3d & mpx )
{
Surface * SfCourant = P1Surface;
while ( SfCourant )
{
SfCourant->TracerPolyedre( mpx );
SfCourant = SfCourant->PSuivant;
}
}

//============================================================================
// Identique a Ligne::IntervalleNURBS, mais avec l'object explicite U

int IntervalleNURBS( R u, R *U, int nbn, int degre) {
int min, max, milieu;

if( u >= U[1+nbn] ) return nbn; //Cas particulier!

//recherche par dichotomie
min = degre;
max = nbn+1;

Iteration:
milieu = (min+max)/2;

if( u < U[milieu] ) {
max = milieu;
goto Iteration;
} else if( u >= U[milieu+1] ) {
min = milieu;
goto Iteration;
} else {
return milieu;
}
}

//============================================================================
// Identique a Ligne::ValeurBaseNURBS, mais avec l'objet explicite U
void ValeurBaseNURBS( int i, Z NbEtape, R u, R U[], R N[] )
{
int j,k;
R g[1+64], d[1+64];//numerateur gauche et droite de N[i,j]
R s, t; //denominateur=numerateur gauche + droite!

N[0] = 1.0;
for (j=1; j<=NbEtape; j++)
{
g[j] = u - U[i+1-j]; //numerateur gauche j
d[j] = U[i+j] - u; //numerateur droite j
s=0.0;
for (k=0; k<j; k++)
{
t = N[k] / ( d[k+1] + g[j-k] ); //la somme des indices vaut j+1
N[k]= s + d[k+1] * t;
s = g[j-k] * t;
}
N[j]=s;
}
}