1 line
26 KiB
C++
Executable file
1 line
26 KiB
C++
Executable file
// FICHIER : Met_Sfe1.cp
|
|
// CLASSE : Met_Sfe1
|
|
|
|
|
|
#include <iostream>
|
|
#include <math.h>
|
|
#include <stdlib.h>
|
|
#include "Sortie.h"
|
|
#include "Util.h"
|
|
#include "MathUtil.h"
|
|
|
|
|
|
#include "Met_Sfe1.h"
|
|
|
|
|
|
// =========================================================================================
|
|
// vu la taille des executables le fichier est decompose en deux
|
|
// le premier : Met_Sfe1s1.cp concerne les constructeurs, destructeur et
|
|
// la gestion des variables
|
|
// le second : Met_Sfe1s2.cp concerne le calcul
|
|
// =========================================================================================
|
|
|
|
// METHODES PUBLIQUES :
|
|
// ------------------------ calculs ----------------------------------------
|
|
// cas explicite, toutes les grandeurs sont a 0 ou t
|
|
/* Met_Sfe1::Expli Met_Sfe1::Cal_explicit
|
|
( Tableau<Noeud *>& tab_noeud, Mat_pleine& dphi,int nombre_noeud)
|
|
{ Met_abstraite::Calcul_giB_0(tab_noeud,dphi,nombre_noeud); // calcul de la base a t=0
|
|
Calcul_giB_t(tab_noeud,dphi,nombre_noeud); // calcul de la base a t
|
|
Calcul_gijBB_0 (); // metrique base naturelle
|
|
Calcul_gijBB_t (); // "
|
|
Calcul_gijHH_t (); // composantes controvariantes
|
|
Jacobien_t(); // calcul du jacobien a t
|
|
|
|
D_giB_t(dphi); //
|
|
D_gijBB_t(); //variation de la metrique en BB
|
|
// retour des infos
|
|
Expli ex;
|
|
ex.gijBB_0 = gijBB_0;
|
|
ex.gijBB_t = gijBB_t;
|
|
ex.gijHH_t = gijHH_t;
|
|
ex.d_gijBB_t = &d_gijBB_t;
|
|
ex.jacobien = jacobien_t;
|
|
// liberation des tenseurs intermediaires
|
|
LibereTenseur();
|
|
return ex;
|
|
}; */
|
|
|
|
// cas implicite, toutes les grandeurs sont a 0 ou t+dt
|
|
/*Met_abstraite::Impli Met_Sfe1::Cal_implicit
|
|
( Tableau<Noeud *>& tab_noeud, Mat_pleine& dphi, int nombre_noeud)
|
|
{ // tout d'abord on calcul les elements relatifs a la facette plane
|
|
Met_abstraite::Calcul_giB_0(tab_noeud,dphi,nombre_noeud); // calcul de la base a t=0
|
|
*aiB_0 = *giB_0;
|
|
Met_abstraite::Calcul_giB_t(tab_noeud,dphi,nombre_noeud); // calcul de la base a t
|
|
*aiB_t = *giB_t;
|
|
Met_abstraite::Calcul_giB_tdt(tab_noeud,dphi,nombre_noeud); // calcul de la base a tdt
|
|
*aiB_tdt = *giB_tdt;
|
|
Met_abstraite::Calcul_gijBB_0 (); // metrique base naturelle
|
|
Met_abstraite::Calcul_gijBB_t (); // "
|
|
Met_abstraite::Calcul_gijBB_tdt (); // "
|
|
Met_abstraite::Calcul_gijHH_tdt (); // composantes controvariantes
|
|
*aijBB_0 = *gijBB_0;
|
|
*aijBB_t = *gijBB_t;
|
|
*aijBB_tdt = *gijBB_tdt;
|
|
*aijHH_tdt = *gijHH_tdt;
|
|
Met_abstraite::Jacobien_tdt(); // calcul du jacobien a t
|
|
ajacobien_tdt = jacobien_tdt;
|
|
*/
|
|
|
|
/* D_giB_tdt(dphi); //
|
|
D_gijBB_tdt(); //variation de la metrique en BB
|
|
Calcul_giH_tdt(); // base duale
|
|
D_giH_tdt(); // variation de la base duale
|
|
D_gijHH_tdt(); //variation de la metrique en HH
|
|
Djacobien_tdt(); // variation du jacobien */
|
|
// retour des infos
|
|
/* Impli ex;
|
|
ex.gijBB_0 = gijBB_0;
|
|
ex.gijBB_t = gijBB_t; ex.gijBB_tdt = gijBB_tdt;
|
|
ex.gijHH_tdt = gijHH_tdt;
|
|
ex.d_gijBB_tdt = &d_gijBB_tdt; ex.d_gijHH_tdt = &d_gijHH_tdt;
|
|
ex.jacobien = jacobien_tdt; ex.d_jacobien_tdt = &d_jacobien_tdt;
|
|
// liberation des tenseurs intermediaires
|
|
LibereTenseur();
|
|
return ex;
|
|
|
|
};*/
|
|
/* // pour la remontee au infos duaux
|
|
Met_abstraite::InfoImp Met_Sfe1::Cal_InfoImp
|
|
( Tableau<Noeud *>& tab_noeud, Mat_pleine& dphi,Vecteur& phi, int nombre_noeud)
|
|
{ Calcul_M0(tab_noeud,phi,nombre_noeud);
|
|
Calcul_Mtdt(tab_noeud,phi,nombre_noeud);
|
|
Calcul_giB_0(tab_noeud,dphi,nombre_noeud); // calcul de la base a t=0
|
|
Calcul_giB_tdt(tab_noeud,dphi,nombre_noeud); // calcul de la base a tdt
|
|
Calcul_gijBB_0 (); // "
|
|
Calcul_gijHH_0 (); // composantes controvariantes
|
|
Calcul_gijBB_tdt (); // "
|
|
Calcul_gijHH_tdt (); // composantes controvariantes
|
|
// bases duales
|
|
Calcul_giH_0();
|
|
Calcul_giH_tdt();
|
|
|
|
// retour des infos
|
|
InfoImp ex;
|
|
ex.M0 = M0;
|
|
ex.Mtdt = Mtdt;
|
|
ex.giB_0 = giB_0;
|
|
ex.giB_tdt = giB_tdt;
|
|
ex.gijBB_tdt = gijBB_tdt;
|
|
ex.gijHH_tdt = gijHH_tdt;
|
|
ex.giH_0 = giH_0;
|
|
ex.giH_tdt = giH_tdt;
|
|
// liberation des tenseurs intermediaires
|
|
LibereTenseur();
|
|
return ex;
|
|
};
|
|
|
|
Met_abstraite::InfoExp Met_Sfe1::Cal_InfoExp
|
|
(Tableau<Noeud *>& tab_noeud, Mat_pleine& dphi,Vecteur& phi, int nombre_noeud)
|
|
{ Calcul_M0(tab_noeud,phi,nombre_noeud);
|
|
Calcul_Mt(tab_noeud,phi,nombre_noeud);
|
|
Calcul_giB_0(tab_noeud,dphi,nombre_noeud); // calcul de la base a t=0
|
|
Calcul_giB_t(tab_noeud,dphi,nombre_noeud); // calcul de la base a t
|
|
Calcul_gijBB_0 (); // "
|
|
Calcul_gijHH_0 (); // composantes controvariantes
|
|
Calcul_gijBB_t (); // "
|
|
Calcul_gijHH_t (); // composantes controvariantes
|
|
// bases duales
|
|
Calcul_giH_0();
|
|
Calcul_giH_t();
|
|
// retour des infos
|
|
InfoExp ex;
|
|
ex.M0 = M0;
|
|
ex.Mt = Mt;
|
|
ex.giB_0 = giB_0;
|
|
ex.giB_t = giB_t;
|
|
ex.gijBB_t = gijBB_t;
|
|
ex.gijHH_t = gijHH_t;
|
|
ex.giH_0 = giH_0;
|
|
ex.giH_t = giH_t;
|
|
// liberation des tenseurs intermediaires
|
|
LibereTenseur();
|
|
return ex;
|
|
};
|
|
|
|
// ========== utilise par le contact =========================
|
|
|
|
// calcul d'un point M0 en fonction de phi et des coordonnees a 0
|
|
Coordonnee & Met_Sfe1::PointM_0 // ( stockage a t=0)
|
|
(Tableau<Noeud *>& tab_noeud,Vecteur& Phi)
|
|
{ Calcul_M0(tab_noeud,Phi,tab_noeud.Taille());
|
|
return *M0;
|
|
};
|
|
|
|
// calcul d'un point M0 en fonction de phi et des coordonnees a t
|
|
Coordonnee & Met_Sfe1::PointM_t // ( stockage a t)
|
|
(Tableau<Noeud *>& tab_noeud,Vecteur& Phi)
|
|
{ Calcul_Mt(tab_noeud,Phi,tab_noeud.Taille());
|
|
return *Mt;
|
|
};
|
|
|
|
// calcul d'un point Mt en fonction de phi et des coordonnees a tdt
|
|
Coordonnee & Met_Sfe1::PointM_tdt // ( stockage a tdt)
|
|
(Tableau<Noeud *>& tab_noeud,Vecteur& Phi)
|
|
{ Calcul_Mtdt(tab_noeud,Phi,tab_noeud.Taille());
|
|
return *Mtdt;
|
|
};
|
|
|
|
// calcul de la base naturel ( stockage a t=0) au point correspondant au dphi
|
|
BaseB& Met_Sfe1::BaseNat_0 // en fonction des coordonnees a t=0
|
|
(Tableau<Noeud *>& tab_noeud,Mat_pleine& dphi)
|
|
{ Calcul_giB_0 (tab_noeud,dphi,tab_noeud.Taille());
|
|
return *giB_0;
|
|
};
|
|
|
|
// calcul de la base naturel ( stockage a t) au point correspondant au dphi
|
|
BaseB& Met_Sfe1::BaseNat_t // en fonction des coordonnees a t
|
|
(Tableau<Noeud *>& tab_noeud,Mat_pleine& dphi)
|
|
{ Calcul_giB_t (tab_noeud,dphi,tab_noeud.Taille());
|
|
return *giB_t;
|
|
};
|
|
|
|
// calcul de la base naturel ( stockage a t=tdt) au point correspondant au dphi
|
|
BaseB& Met_Sfe1::BaseNat_tdt // en fonction des coordonnees a tdt
|
|
(Tableau<Noeud *>& tab_noeud,Mat_pleine& dphi)
|
|
{ Calcul_giB_tdt (tab_noeud,dphi,tab_noeud.Taille());
|
|
return *giB_tdt;
|
|
};
|
|
|
|
// calcul de la base naturelle et duale ( stockage a t=0) en fonction des coord a 0
|
|
void Met_Sfe1::BaseND_0(Tableau<Noeud *>& tab_noeud,Mat_pleine& dphi,BaseB& bB,BaseH& bH)
|
|
{ Calcul_giB_0 (tab_noeud,dphi,tab_noeud.Taille());
|
|
Calcul_gijBB_0 (); // "
|
|
Calcul_gijHH_0 (); // composantes controvariantes
|
|
Calcul_giH_0 ();
|
|
bB = *giB_0;
|
|
bH = *giH_0;
|
|
};
|
|
// calcul de la base naturelle et duale ( stockage a t) en fonction des coord a t
|
|
void Met_Sfe1::BaseND_t(Tableau<Noeud *>& tab_noeud,Mat_pleine& dphi,BaseB& bB,BaseH& bH)
|
|
{ Calcul_giB_t (tab_noeud,dphi,tab_noeud.Taille());
|
|
Calcul_gijBB_t (); // "
|
|
Calcul_gijHH_t (); // composantes controvariantes
|
|
Calcul_giH_t ();
|
|
bB = *giB_t;
|
|
bH = *giH_t;
|
|
};
|
|
// calcul de la base naturelle et duale ( stockage a tdt) en fonction des coord a tdt
|
|
void Met_Sfe1::BaseND_tdt(Tableau<Noeud *>& tab_noeud,Mat_pleine& dphi,BaseB& bB,BaseH& bH)
|
|
{ Calcul_giB_tdt (tab_noeud,dphi,tab_noeud.Taille());
|
|
Calcul_gijBB_tdt (); // "
|
|
Calcul_gijHH_tdt (); // composantes controvariantes
|
|
Calcul_giH_tdt ();
|
|
bB = *giB_tdt;
|
|
bH = *giH_tdt;
|
|
}; */
|
|
|
|
//==================== METHODES PROTEGEES===============================
|
|
|
|
/* // calcul de la base naturel a t0 // a changer
|
|
void Met_Sfe1::Calcul_giB_0
|
|
{}; */
|
|
//-----------// calcul du tenseur de courbure dans la base naturelle
|
|
// plusieurs cas sont etudies suivant l'instant considere
|
|
// a l'instant t = 0
|
|
Vecteur& Met_Sfe1::courbure_0 (Tableau<Noeud *>& tab_noeud)
|
|
{ Vecteur & aiB1 = (*aiB_0)(1); Vecteur & aiB2 = (*aiB_0)(2); // pour simplifier
|
|
Vecteur & aiH1 = (*aiH_0)(1); Vecteur & aiH2 = (*aiH_0)(2); // pour simplifier
|
|
Tableau<Coordonnee> tab_coor(6);
|
|
for (int i=1;i<=6;i++)
|
|
tab_coor(i) = tab_noeud(i)->Coord0();
|
|
return courbure1(tab_coor,aiB1,aiB2,aiH1,aiH2,tab_noeud);
|
|
};
|
|
// a l'instant t
|
|
Vecteur& Met_Sfe1::courbure_t (Tableau<Noeud *>& tab_noeud)
|
|
{ Vecteur & aiB1 = (*aiB_t)(1); Vecteur & aiB2 = (*aiB_t)(2); // pour simplifier
|
|
Vecteur & aiH1 = (*aiH_t)(1); Vecteur & aiH2 = (*aiH_t)(2); // pour simplifier
|
|
Tableau<Coordonnee> tab_coor(6);
|
|
for (int i=1;i<=6;i++)
|
|
tab_coor(i) = tab_noeud(i)->Coord1();
|
|
return courbure1(tab_coor,aiB1,aiB2,aiH1,aiH2,tab_noeud);
|
|
};
|
|
// a l'instant t+dt
|
|
Vecteur& Met_Sfe1::courbure_tdt (Tableau<Noeud *>& tab_noeud)
|
|
{ Vecteur & aiB1 = (*aiB_tdt)(1); Vecteur & aiB2 = (*aiB_tdt)(2); // pour simplifier
|
|
Vecteur & aiH1 = (*aiH_tdt)(1); Vecteur & aiH2 = (*aiH_tdt)(2); // pour simplifier
|
|
Tableau<Coordonnee> tab_coor(6);
|
|
for (int i=1;i<=6;i++)
|
|
tab_coor(i) = tab_noeud(i)->Coord2();
|
|
return courbure1(tab_coor,aiB1,aiB2,aiH1,aiH2,tab_noeud);
|
|
};
|
|
|
|
// routine generale de calcul de la courbure
|
|
Vecteur& Met_Sfe1::courbure
|
|
(Tableau<Coordonnee>& tab_coor,Vecteur & aiB1,Vecteur & aiB2,
|
|
Vecteur & aiH1,Vecteur & aiH2,Tableau<Noeud *>& tab_noeud)
|
|
{ Vecteur curbp(3),curb(3); // : lestenseurs de courbures d'abord dans les axes
|
|
// locaux des arretes du triangle puis dans le repere naturel: b11,b12,b22
|
|
Mat_pleine mat(3,3); // matrice pour calculer le tenseur de courbure
|
|
// calcul de la normale centrale
|
|
Vecteur N = (ProdVec(aiB1,aiB2)).Normer();
|
|
// calcul du centre de gravite
|
|
Coordonnee G = (tab_coor(1) + tab_coor(2) + tab_coor(3)) / 3.;
|
|
// le meme calcul est effectue pour les trois cotes du triangle principale
|
|
//1- tout d'abord un tableau d'indice pour ne pas etre embete par les modulos 3
|
|
Vecteur indi(5);indi(1) = 1; indi(2) = 2; indi(3) = 3; indi(4) = 1; indi(5) = 2;
|
|
// 2- on boucle sur le cotes du triangle
|
|
Vecteur Ni(3); // normale du triangle exterieur
|
|
for (int ncot=1;ncot<=3;ncot++)
|
|
{ // on utilise les notations : triangle principal DBC et secondaire BAC
|
|
// on cherche a calculer la longueur H2G2
|
|
// simplification de l'ecriture
|
|
Coordonnee& D = tab_noeud(indi(ncot))->Coord0();
|
|
Coordonnee& B = tab_noeud(indi(ncot+1))->Coord0();
|
|
Coordonnee& C = tab_noeud(indi(ncot+2))->Coord0();
|
|
Coordonnee& A = tab_noeud(ncot+1)->Coord0();
|
|
// calcul de G2, centre de gravite du triangle exterieur
|
|
Coordonnee G2 = (B + C + A) / 3.;
|
|
// calcul de U1
|
|
Vecteur CB = B - C;
|
|
Vecteur U1 = CB / CB.Norme();
|
|
// vecteur CG2
|
|
Vecteur CG2 = G2 - C;
|
|
// longueur H2G2 et GH1
|
|
Vecteur H2G2 = CG2 - (CG2 * U1) * U1;
|
|
double h2G2 = CG2.Norme();
|
|
Vecteur CG = G - C;
|
|
double GH1 = sqrt ( CG * CG - SQR( CG * U1)) ;
|
|
// calcul de la normale du triangle exterieur
|
|
Vecteur CA = A - C;
|
|
Ni = (ProdVec(CB,CA)).Normer();
|
|
// calcul de la derivee de la normale
|
|
Vecteur dNdV = ( Ni - N) / ( GH1 + h2G2);
|
|
// courbure dans la direction u2 c-a-d H2G2
|
|
curbp(ncot) = dNdV * H2G2;
|
|
// calcul de la matrice de passage de Ualpha a aalpha
|
|
double Beta21 = H2G2 * aiH1;
|
|
double Beta22 = H2G2 * aiH2;
|
|
// remplissage de la matrice mat
|
|
mat(ncot,1) = Beta21*Beta21;
|
|
mat(ncot,2) = Beta21*Beta22;
|
|
mat(ncot,3) = Beta22*Beta22;
|
|
};
|
|
|
|
// calcul du tenseur de courbure dans la base ai
|
|
Mat_pleine mati = mat.Inverse();
|
|
curb = mati * curbp;
|
|
// retour de la courbure
|
|
return curb;
|
|
};
|
|
|
|
//--------------// calcul du tenseur de courbure et de sa variation
|
|
// plusieurs cas sont etudies suivant l'instant considere
|
|
// a l'instant t
|
|
void Met_Sfe1::Dcourbure_t
|
|
(Tableau<Noeud *>& tab_noeud,Vecteur& curb,TabOper<Vecteur>& dcurb)
|
|
{ Vecteur & aiB1 = (*aiB_t)(1); Vecteur & aiB2 = (*aiB_t)(2); // pour simplifier
|
|
Vecteur & aiH1 = (*aiH_t)(1); Vecteur & aiH2 = (*aiH_t)(2); // pour simplifier
|
|
Tableau<Coordonnee> tab_coor(6);
|
|
for (int i=1;i<=6;i++)
|
|
tab_coor(i) = tab_noeud(i)->Coord1();
|
|
Dcourbure1 (tab_coor,aiB1,aiB2,aiH1,aiH2,*d_aiB_t,*d_aiH_t,tab_noeud,curb,dcurb);
|
|
};
|
|
// a l'instant t+dt
|
|
void Met_Sfe1::Dcourbure_tdt
|
|
(Tableau<Noeud *>& tab_noeud,Vecteur& curb,TabOper<Vecteur>& dcurb)
|
|
{ Vecteur & aiB1 = (*aiB_tdt)(1); Vecteur & aiB2 = (*aiB_tdt)(2); // pour simplifier
|
|
Vecteur & aiH1 = (*aiH_tdt)(1); Vecteur & aiH2 = (*aiH_tdt)(2); // pour simplifier
|
|
Tableau<Coordonnee> tab_coor(6);
|
|
for (int i=1;i<=6;i++)
|
|
tab_coor(i) = tab_noeud(i)->Coord2();
|
|
Dcourbure1 (tab_coor,aiB1,aiB2,aiH1,aiH2,*d_aiB_t,*d_aiH_t,tab_noeud,curb,dcurb);
|
|
};
|
|
|
|
// routine generale de calcul de la courbure et de sa variation
|
|
// en sortie : curb , la courbure b11,b12,b22
|
|
// dcurb , la variation de courbure.
|
|
void Met_Sfe1::Dcourbure
|
|
(Tableau<Coordonnee>& tab_coor,Vecteur & aiB1,Vecteur & aiB2,
|
|
Vecteur & aiH1,Vecteur & aiH2,Tableau <BaseB>& DaiB,Tableau <BaseH>& DaiH,
|
|
Tableau<Noeud *>& tab_noeud,Vecteur& curb,TabOper <Vecteur>& dcurb)
|
|
{ int nbddl = 18; // nombre de ddl total
|
|
int nbddlint = 9; // nombre de ddl interne
|
|
Vecteur Nul(3); // le vecteur nul de dim 3, pour les initialisations
|
|
|
|
Vecteur curbp(3); // : les tenseurs de courbures dans les axes
|
|
// locaux des arretes du triangle
|
|
Tableau <Vecteur> dcurbp(nbddl); // variation des courbures
|
|
Mat_pleine mat(3,3); // matrice pour calculer le tenseur de courbure
|
|
Tableau <Mat_pleine> dmat(nbddl,mat);
|
|
// calcul de la normale centrale
|
|
Vecteur NN = ProdVec(aiB1,aiB2); // normale non normee
|
|
double nor = NN.Norme(); // la norme
|
|
Vecteur N = NN / nor;
|
|
Vecteur DNN; // vecteur intermediair
|
|
Tableau <Vecteur> dN(nbddl,Nul);
|
|
for (int i1=1;i1<= nbddlint;i1++)
|
|
{ DNN = ProdVec(DaiB(i1)(1),aiB2) + ProdVec(aiB1,DaiB(i1)(2));
|
|
dN(i1) = VarUnVect(NN,DNN,nor);
|
|
}
|
|
// calcul du centre de gravite
|
|
Coordonnee G = (tab_coor(1) + tab_coor(2) + tab_coor(3)) / 3.;
|
|
Tableau <Coordonnee> dG(nbddl);
|
|
for (int ia=1;ia<= 3;ia++)
|
|
{dG(ia) = Ia(ia)/3. ;dG(ia+3) = Ia(ia)/3. ;dG(ia+6) = Ia(ia)/3. ;}
|
|
|
|
// le meme calcul est effectue pour les trois cotes du triangle principale
|
|
//1- tout d'abord un tableau d'indice pour ne pas etre embete par les modulos 3
|
|
Tableau <int> indi(5);indi(1) = 1; indi(2) = 2; indi(3) = 3; indi(4) = 1; indi(5) = 2;
|
|
// 2- on boucle sur le cotes du triangle
|
|
Vecteur Ni(3); // normale du triangle exterieur
|
|
// -------- declaration pour la boucle sur les 3 cotes ------------
|
|
Coordonnee G2(3); // init a zero et a 3 composantes
|
|
Tableau <Coordonnee> dD(nbddl,Nul),dB(nbddl,Nul),dC(nbddl,Nul),dA(nbddl,Nul) ; // variation des points
|
|
Tableau <Coordonnee> dG2(nbddl,Nul);
|
|
Vecteur CB(3),U1(3),CG2(3),H2G2(3);
|
|
Tableau <Vecteur> dCB(nbddl,Nul),dU1(nbddl,Nul),dH2G2(nbddl,Nul);
|
|
Tableau <double> dh2G2(nbddl);
|
|
Vecteur CG(3),CA(3),dNdV(3);
|
|
Tableau <Vecteur> dCA(nbddl,Nul),dNi(nbddl,Nul);
|
|
double under3 = 1./3.;
|
|
Vecteur dCG(3),dGH1(3),ddNdv(3),dCG2_jnB,dCG2_jnC,dCG2_jnA;
|
|
for (int ncot=1;ncot<=3;ncot++)
|
|
{ // on utilise les notations : triangle principal DBC et secondaire BAC
|
|
// on cherche a calculer la longueur H2G2
|
|
// simplification de l'ecriture
|
|
Coordonnee& D = tab_noeud(indi(ncot))->Coord0();
|
|
Coordonnee& B = tab_noeud(indi(ncot+1))->Coord0();
|
|
Coordonnee& C = tab_noeud(indi(ncot+2))->Coord0();
|
|
Coordonnee& A = tab_noeud(ncot+1)->Coord0();
|
|
// calcul de G2, centre de gravite du triangle exterieur
|
|
G2 = (B + C + A) * under3;
|
|
// calcul de U1
|
|
CB = B - C;
|
|
double norCB = CB.Norme();
|
|
U1 = CB / norCB;
|
|
// vecteur CG2
|
|
CG2 = G2 - C;
|
|
// longueur H2G2 et GH1
|
|
H2G2 = CG2 - (CG2 * U1) * U1;
|
|
double h2G2 = CG2.Norme();
|
|
double unSurh2G2= 1./h2G2;
|
|
CG = G - C;
|
|
double GH1 = sqrt ( CG * CG - SQR( CG * U1)) ;
|
|
// calcul de la normale du triangle exterieur
|
|
CA = A - C;
|
|
Ni = ProdVec(CB,CA);
|
|
double norNi = Ni.Norme();
|
|
Ni /= norNi;
|
|
// calcul de la derivee de la normale
|
|
dNdV = ( Ni - N) / ( GH1 + h2G2);
|
|
// courbure dans la direction u2 c-a-d H2G2
|
|
curbp(ncot) = dNdV * H2G2;
|
|
// calcul de la matrice de passage de Ualpha a aalpha
|
|
double Beta21 = H2G2 * aiH1;
|
|
double Beta22 = H2G2 * aiH2;
|
|
// remplissage de la matrice mat
|
|
mat(ncot,1) = Beta21*Beta21;
|
|
mat(ncot,2) = Beta21*Beta22;
|
|
mat(ncot,3) = Beta22*Beta22;
|
|
|
|
// calcul des variations
|
|
for (int ia=1;ia<= 3;ia++)
|
|
{ // on a 4 indices pour chaque ia, pour lesquels la variation est non nulle
|
|
// dans cette boucle on n'utilise pas systematiquement les 4 indices
|
|
int jnD = (indi(ncot)-1)*3+ia; int jnB = (indi(ncot+1)-1)*3+ia;
|
|
int jnC = (indi(ncot+2)-1)*3+ia; int jnA = 9+(ncot-1)*3+ia;
|
|
// tout d'abord les variations des points
|
|
dD(jnD) = Ia(ia);dB(jnB) = Ia(ia);
|
|
dC(jnC) = Ia(ia);dA(jnA) = Ia(ia);
|
|
// centre de gravite
|
|
dG2(jnB) = Ia(ia)*under3;
|
|
dG2(jnC) = Ia(ia)*under3;
|
|
dG2(jnA) = Ia(ia)*under3;
|
|
// vecteur CB
|
|
dCB(jnB) = dB(jnB); dCB(jnC) = - dC(jnC);
|
|
// vecteur U1
|
|
dU1(jnB) = VarUnVect(CB,dCB(jnB),norCB);
|
|
dU1(jnC) = VarUnVect(CB,dCB(jnC),norCB);
|
|
// vecteur CG2
|
|
dCG2_jnB = dG2(jnB) - dC(jnB);
|
|
dCG2_jnC = dG2(jnC) - dC(jnC);
|
|
dCG2_jnA = dG2(jnA);
|
|
// vecteur H2G2= CG2 - (CG2 * U1) * U1;
|
|
dH2G2(jnB) = dCG2_jnB-(dCG2_jnB*U1 + CG2*dU1(jnB))* U1- (CG2 * U1)* dU1(jnB);
|
|
dH2G2(jnC) = dCG2_jnC-(dCG2_jnC*U1 + CG2*dU1(jnC))* U1- (CG2 * U1)* dU1(jnC);
|
|
dH2G2(jnA) = dCG2_jnA-(dCG2_jnA*U1 + CG2*dU1(jnA))* U1- (CG2 * U1)* dU1(jnA);
|
|
// scalaire h2G2 = CG2.Norme();
|
|
dh2G2(jnB) = dH2G2(jnB) * H2G2 /unSurh2G2;
|
|
dh2G2(jnC) = dH2G2(jnC) * H2G2 /unSurh2G2;
|
|
dh2G2(jnA) = dH2G2(jnA) * H2G2 /unSurh2G2;
|
|
// Vecteur CA = A - C
|
|
dCA(jnC) = -dC(jnC); dCA(jnA) = dA(jnA);
|
|
// Vecteur Ni = (ProdVec(CB,CA)).Normer()
|
|
dNi(jnC) = VarUnVect(Ni,(ProdVec(dCB(jnC),CA)+ProdVec(CA,dCA(jnC))),norNi);
|
|
dNi(jnB) = VarUnVect(Ni,(ProdVec(dCB(jnB),CA)+ProdVec(CA,dCA(jnB))),norNi);
|
|
dNi(jnA) = VarUnVect(Ni,(ProdVec(dCB(jnA),CA)+ProdVec(CA,dCA(jnA))),norNi);
|
|
}
|
|
for (int ib=1;ib<= 3;ib++)
|
|
{ // on a 4 indices pour chaque ia, pour lesquels la variation est non nulle
|
|
// mais ici on utilise systematiquement les 4 indices d'ou une seconde boucle
|
|
// egalement moins de stockage intermediaire
|
|
int jnb[4] = {(indi(ncot)-1)*3+ib,(indi(ncot+1)-1)*3+ib,
|
|
(indi(ncot+2)-1)*3+ib,9+(ncot-1)*3+ib};
|
|
for (int jb=0;jb<=3;jb++)
|
|
{ // vecteur CG = G - C
|
|
dCG = dG(jnb[jb]) - dC(jnb[jb]);
|
|
// scalaire GH1 = sqrt ( CG * CG - SQR( CG * U1))
|
|
dGH1 = (dCG * CG - (CG * U1) * (dCG * U1 + CG * dU1(jnb[jb])))/ GH1 ;
|
|
// vecteur dNdV = ( Ni - N) / ( GH1 + h2G2)
|
|
ddNdv = (dNi(jnb[jb]) - dN(jnb[jb])) / ( GH1 + h2G2)
|
|
- ( dGH1 + dh2G2(jnb[jb])) * (Ni - N) / ( SQR ( GH1 + h2G2));
|
|
// courbure curbp(ncot) = dNdV * H2G2;
|
|
dcurbp(jnb[jb])(ncot) = ddNdv * H2G2 + dNdV * dH2G2(jnb[jb]);
|
|
// matrice de passage : Beta21 = H2G2 * aiH1 et Beta22 = H2G2 * aiH2
|
|
double dBeta21 = dH2G2(jnb[jb]) * aiH1 + H2G2 * DaiH(1)(jnb[jb]);
|
|
double dBeta22 = dH2G2(jnb[jb]) * aiH2 + H2G2 * DaiH(2)(jnb[jb]);
|
|
// matrice mat : mat(ncot,1) = Beta21*Beta21 ; mat(ncot,2) = Beta21*Beta22;
|
|
// mat(ncot,3) = Beta22*Beta22;
|
|
dmat(jnb[jb])(ncot,1) = 2. * dBeta21 * Beta21;
|
|
dmat(jnb[jb])(ncot,2) = dBeta21 * Beta22 + Beta21 * dBeta22;
|
|
dmat(jnb[jb])(ncot,3) = dBeta22 * Beta22 + Beta22 * dBeta22;
|
|
}
|
|
}
|
|
};
|
|
|
|
// calcul du tenseur de courbure dans la base ai
|
|
Mat_pleine mati = mat.Inverse();
|
|
curb = mati * curbp;
|
|
// variations
|
|
Mat_pleine dmati(mati);
|
|
for (int ddl=1;ddl<=nbddl;ddl++)
|
|
{ dmati = - mat * ( dmat(ddl) * mati);
|
|
dcurb(ddl) = dmati * curbp + mati * dcurbp(ddl);
|
|
}
|
|
};
|
|
|
|
|
|
|
|
/*//
|
|
void Met_Sfe1::courbure2
|
|
(Tableau<Noeud *>& tab_noeud, Mat_pleine& dphi, int nombre_noeud)
|
|
{ // tout d'abord on calcul la base naturelle du triangle central
|
|
if (giB_0 =! NULL)
|
|
Calcul_giB_0 (tab_noeud,dphi,3); // 3 pour les 3 noeuds centrals
|
|
Vecteur & aiB1 = (*giB_0)(1); Vecteur & aiB2 = (*giB_0)(2); // pour simplicite
|
|
Calcul_gijBB_0 (); // puis la metrique "
|
|
Calcul_gijHH_0 (); // composantes controvariantes
|
|
// les bases duales
|
|
Calcul_giH_0();
|
|
Vecteur & aiH1 = (*giH_0)(1); Vecteur & aiH2 = (*giH_0)(2); // pour simplicite
|
|
|
|
Vecteur curbp(3),curb(3); // lestenseurs de courbures d'abord dans les axes
|
|
// locaux des arretes du triangle puis dans le repere naturel: b11,b12,b22
|
|
Mat_pleine mat(3,3); // matrice pour calculer le tenseur de courbure
|
|
// calcul de la normale centrale
|
|
Vecteur N = (ProdVec(aiB1,aiB2)).Normer();
|
|
// calcul du centre de gravite
|
|
Coordonnee G = (tab_noeud(1)->Coord0() + tab_noeud(2)->Coord0()
|
|
+ tab_noeud(3)->Coord0()) / 3.;
|
|
// le meme calcul est effectue pour les trois cotes du triangle principale
|
|
//1- tout d'abord un tableau d'indice pour ne pas etre embete par les modulos 3
|
|
Vecteur indi(5);indi(1) = 1; indi(2) = 2; indi(3) = 3; indi(4) = 1; indi(5) = 2;
|
|
// 2- on boucle sur le cotes du triangle
|
|
Vecteur Ni(3); // normale du triangle exterieur
|
|
for (int ncot=1;ncot<=3;ncot++)
|
|
{ // on utilise les notations : triangle principal DBC et secondaire BAC
|
|
// on cherche l'image de par rotation du plan BAC dans le plan DBC
|
|
// le point A devient alors Ap
|
|
// calcul des vecteurs cote
|
|
Vecteur CA = tab_noeud(ncot+3)->Coord0() - tab_noeud(indi(ncot)+2)->Coord0();
|
|
Vecteur CB = tab_noeud(indi(ncot)+1)->Coord0() - tab_noeud(indi(ncot)+2)->Coord0();
|
|
Vecteur BA = tab_noeud(ncot+3)->Coord0() - tab_noeud(indi(ncot)+1)->Coord0();
|
|
// calcul du terme (i) appele xi
|
|
double CA2 = CA * CA;
|
|
double xi = CA2 + CB * CB - BA * BA;
|
|
// calcul des coefficients a et b
|
|
double g1CB = aiB1 * CB;
|
|
double a = -(aiB2 * CB) / g1CB ;
|
|
double b = xi * 0.5 / g1CB;
|
|
// calcul des coefficients de l'equation du second degre
|
|
Vecteur ag1g2 = a * aiB1 + aiB2;
|
|
double aa = ag1g2 * ag1g2;
|
|
double bb = (2 * b) * (aiB1 * ag1g2);
|
|
double cc = (b*b) * (aiB1 * aiB1) - CA2;
|
|
// resolution de l'equation du second degre
|
|
double x1,x2; int isol;
|
|
// ResoEqua2(aa,bb,cc,x1,x2,isol);
|
|
// calcul de la position du point externe sur le plan de la facette centrale
|
|
double teta1 = x1 * a + b; // choix de x1 a priori pour teta2
|
|
Coordonnee Ap = teta1 * aiB1 + x1 * aiB2;
|
|
// verif du point
|
|
Vecteur DC = tab_noeud(indi(ncot)+2)->Coord0() - tab_noeud(indi(ncot))->Coord0();
|
|
Vecteur CAp = Ap - tab_noeud(indi(ncot)+2)->Coord0();
|
|
if (DC * CAp < 0)
|
|
{ // cas ou c'etait x2 le teta2
|
|
teta1 = x2 * a + b;
|
|
Ap = teta1 * aiB1 + x2 * aiB2;
|
|
CAp = Ap - tab_noeud(indi(ncot)+2)->Coord0();
|
|
}
|
|
// calcul de la normale du triangle exterieur
|
|
Ni = (ProdVec(CB,CAp)).Normer();
|
|
// calcul de GH1 et GpH2
|
|
Vecteur CG = G - tab_noeud(indi(ncot)+2)->Coord0();
|
|
Vecteur U1 = CB / CB.Norme();
|
|
double GH1 = sqrt ( CG * CG - ( CG * U1) * ( CG * U2)) ; // a modifier !!!!
|
|
double H2Gp = sqrt ( CG * CGp - ( CGp * U2) * ( CGp * U2)) ; // a modifier !!!!
|
|
// calcul de la derivee de la normale
|
|
Vecteur dNdV = ( Ni - N) / ( GH1 + H2Gp);
|
|
// calcul du vecteur normal a U1, c'a dire U2, mais qui n'est pas norme
|
|
Vecteur U2 = CAp - (CAp * U1) * U1;
|
|
// courbure dans la direction u22
|
|
curbp(ncot) = dNdV * U2;
|
|
// calcul de la matrice de passage de Ualpha a aalpha
|
|
double Beta21 = U2 * giH(1);
|
|
double Beta22 = U2 * giH(2);
|
|
// remplissage de la matrice mat
|
|
mat(ncot,1) = Beta21*Beta21;
|
|
mat(ncot,2) = Beta21*Beta22;
|
|
mat(ncot,3) = Beta22*Beta22;
|
|
};
|
|
|
|
// calcul du tenseur de courbure dans la base gi
|
|
Mat_pleine mati = mat.Inverse();
|
|
curb = mati * curbp;
|
|
};
|
|
|
|
|
|
// calcul de la base naturel a t
|
|
void Met_Sfe1::Calcul_giB_t
|
|
( Tableau<Noeud *>& tab_noeud, Mat_pleine& dphi, int nombre_noeud)
|
|
{
|
|
for (int a=1;a<= dim_base;a++)
|
|
{ (*(giB_t))(1)(a)= 0.;
|
|
for (int r=1;r<=nombre_noeud;r++)
|
|
{ (*(giB_t))(1)(a)=(*(giB_t))(1)(a) + tab_noeud(r)->Coord1()(a) * dphi(1,r);
|
|
};
|
|
};
|
|
};
|
|
// calcul de la base naturel a tdt
|
|
void Met_Sfe1::Calcul_giB_tdt
|
|
( Tableau<Noeud *>& tab_noeud, Mat_pleine& dphi, int nombre_noeud)
|
|
{
|
|
for (int a=1;a<= dim_base;a++)
|
|
{ (*(giB_tdt))(1)(a)= 0.;
|
|
for (int r=1;r<=nombre_noeud;r++)
|
|
{ (*(giB_tdt))(1)(a)=(*(giB_tdt))(1)(a) + tab_noeud(r)->Coord2()(a) * dphi(1,r);
|
|
};
|
|
};
|
|
}; */
|
|
|
|
|