geometrie/src/fem_triangle3.cpp

//------------------------------------------------------------
//------------------------------------------------------------
//  MAGiC
//  Jean Christophe Cuilli�re et Vincent FRANCOIS
//  D�partement de G�nie M�canique - UQTR
//------------------------------------------------------------
//  Le projet MAGIC est un projet de recherche du d�partement
//  de g�nie m�canique de l'Universit� du Qu�bec �
//  Trois Rivi�res
//  Les librairies ne peuvent �tre utilis�es sans l'accord
//  des auteurs (contact : francois@uqtr.ca)
//------------------------------------------------------------
//------------------------------------------------------------
//
//       fem_triangle3.cpp
//
//------------------------------------------------------------
//------------------------------------------------------------
//  COPYRIGHT 2000
//  Version du 02/03/2006 � 11H22
//------------------------------------------------------------
//------------------------------------------------------------


#include "gestionversion.h"
#include "fem_triangle3.h"
#include "fem_maillage.h"
#include "fem_noeud.h"
#include "mg_element_maillage.h"
#include "math.h"


FEM_TRIANGLE3::FEM_TRIANGLE3(unsigned long num,class MG_ELEMENT_MAILLAGE* mai,class FEM_NOEUD** tabnoeud):FEM_ELEMENT2(num,mai),FEM_TEMPLATE_ELEMENT<3>(tabnoeud)
{
    if (liaison_topologique!=NULL) if (liaison_topologique->get_dimension()==2) liaison_topologique->get_lien_fem_maillage()->ajouter(this);
    tab[0]->get_lien_element2()->ajouter(this);
    tab[1]->get_lien_element2()->ajouter(this);
    tab[2]->get_lien_element2()->ajouter(this);
    get_fem_noeudpetitid()->get_lien_petit_element2()->ajouter(this);
}

FEM_TRIANGLE3::FEM_TRIANGLE3(class MG_ELEMENT_MAILLAGE* mai,FEM_NOEUD** tabnoeud):FEM_ELEMENT2(mai),FEM_TEMPLATE_ELEMENT<3>(tabnoeud)
{
    if (liaison_topologique!=NULL) if (liaison_topologique->get_dimension()==2) liaison_topologique->get_lien_fem_maillage()->ajouter(this);
    tab[0]->get_lien_element2()->ajouter(this);
    tab[1]->get_lien_element2()->ajouter(this);
    tab[2]->get_lien_element2()->ajouter(this);
    get_fem_noeudpetitid()->get_lien_petit_element2()->ajouter(this);
}
FEM_TRIANGLE3::FEM_TRIANGLE3(unsigned long num,class MG_ELEMENT_TOPOLOGIQUE* topo,class FEM_NOEUD** tabnoeud):FEM_ELEMENT2(num,topo),FEM_TEMPLATE_ELEMENT<3>(tabnoeud)
{
    if (liaison_topologique!=NULL) if (liaison_topologique->get_dimension()==2) liaison_topologique->get_lien_fem_maillage()->ajouter(this);
    tab[0]->get_lien_element2()->ajouter(this);
    tab[1]->get_lien_element2()->ajouter(this);
    tab[2]->get_lien_element2()->ajouter(this);
    get_fem_noeudpetitid()->get_lien_petit_element2()->ajouter(this);
}

FEM_TRIANGLE3::FEM_TRIANGLE3(class MG_ELEMENT_TOPOLOGIQUE* topo,FEM_NOEUD** tabnoeud):FEM_ELEMENT2(topo),FEM_TEMPLATE_ELEMENT<3>(tabnoeud)
{
    if (liaison_topologique!=NULL) if (liaison_topologique->get_dimension()==2) liaison_topologique->get_lien_fem_maillage()->ajouter(this);
    tab[0]->get_lien_element2()->ajouter(this);
    tab[1]->get_lien_element2()->ajouter(this);
    tab[2]->get_lien_element2()->ajouter(this);
    get_fem_noeudpetitid()->get_lien_petit_element2()->ajouter(this);
}
FEM_TRIANGLE3::FEM_TRIANGLE3(unsigned long num,class MG_ELEMENT_TOPOLOGIQUE* topo,class MG_ELEMENT_MAILLAGE* mai,class FEM_NOEUD** tabnoeud):FEM_ELEMENT2(num,topo,mai),FEM_TEMPLATE_ELEMENT<3>(tabnoeud)
{
    if (liaison_topologique!=NULL) if (liaison_topologique->get_dimension()==2) liaison_topologique->get_lien_fem_maillage()->ajouter(this);
    tab[0]->get_lien_element2()->ajouter(this);
    tab[1]->get_lien_element2()->ajouter(this);
    tab[2]->get_lien_element2()->ajouter(this);
    get_fem_noeudpetitid()->get_lien_petit_element2()->ajouter(this);
}

FEM_TRIANGLE3::FEM_TRIANGLE3(class MG_ELEMENT_TOPOLOGIQUE* topo,class MG_ELEMENT_MAILLAGE* mai,FEM_NOEUD** tabnoeud):FEM_ELEMENT2(topo,mai),FEM_TEMPLATE_ELEMENT<3>(tabnoeud)
{
    if (liaison_topologique!=NULL) if (liaison_topologique->get_dimension()==2) liaison_topologique->get_lien_fem_maillage()->ajouter(this);
    tab[0]->get_lien_element2()->ajouter(this);
    tab[1]->get_lien_element2()->ajouter(this);
    tab[2]->get_lien_element2()->ajouter(this);
    get_fem_noeudpetitid()->get_lien_petit_element2()->ajouter(this);
}
FEM_TRIANGLE3::FEM_TRIANGLE3(FEM_TRIANGLE3& mdd):FEM_ELEMENT2(mdd),FEM_TEMPLATE_ELEMENT<3>(mdd)
{
    if (liaison_topologique!=NULL) if (liaison_topologique->get_dimension()==2) liaison_topologique->get_lien_fem_maillage()->ajouter(this);
    tab[0]->get_lien_element2()->ajouter(this);
    tab[1]->get_lien_element2()->ajouter(this);
    tab[2]->get_lien_element2()->ajouter(this);
    get_fem_noeudpetitid()->get_lien_petit_element2()->ajouter(this);
}
FEM_TRIANGLE3::~FEM_TRIANGLE3()
{
    if (liaison_topologique==NULL) return;
    if (liaison_topologique->get_dimension()==0) liaison_topologique->get_lien_fem_maillage()->supprimer(this);
    tab[0]->get_lien_element2()->supprimer(this);
    tab[1]->get_lien_element2()->supprimer(this);
    tab[2]->get_lien_element2()->supprimer(this);
    get_fem_noeudpetitid()->get_lien_petit_element2()->supprimer(this);
}


FEM_ELEMENT_MAILLAGE* FEM_TRIANGLE3::dupliquer(FEM_MAILLAGE *femmai,long decalage)
{
    FEM_NOEUD* tabnoeud[3];
    tabnoeud[0]=femmai->get_fem_noeudid(tab[0]->get_id()+decalage);
    tabnoeud[1]=femmai->get_fem_noeudid(tab[1]->get_id()+decalage);
    tabnoeud[2]=femmai->get_fem_noeudid(tab[2]->get_id()+decalage);
    FEM_TRIANGLE3* tri=new FEM_TRIANGLE3(get_id()+decalage,maillage,tabnoeud);
    femmai->ajouter_fem_element2(tri);
    return tri;
}


int FEM_TRIANGLE3::get_type_entite(void)
{
    return   IDFEM_TRIANGLE3;
}

int FEM_TRIANGLE3::get_dimension(void)
{
    return 2;
}


int FEM_TRIANGLE3::get_nb_fem_noeud(void)
{
    return FEM_TEMPLATE_ELEMENT<3>::get_nb_fem_noeud();
}

FEM_NOEUD* FEM_TRIANGLE3::get_fem_noeud(int num)
{
    return FEM_TEMPLATE_ELEMENT<3>::get_fem_noeud(num);
}

void FEM_TRIANGLE3::change_noeud(int num,FEM_NOEUD* noeud)
{
    FEM_TEMPLATE_ELEMENT<3>::change_noeud(num,noeud);
}

BOITE_3D& FEM_TRIANGLE3::get_boite_3D(void)
{
    return FEM_TEMPLATE_ELEMENT<3>::get_boite_3D();
}


void FEM_TRIANGLE3::enregistrer(std::ostream& o)
{
    if (maillage!=NULL)
      if (get_lien_topologie()!=NULL) o << "%" << get_id() << "=FEM_TRIANGLE3($"<< get_lien_topologie()->get_id() << ",$" << maillage->get_id() << ",$" << tab[0]->get_id() << ",$" << tab[1]->get_id() << ",$" << tab[2]->get_id()<<    ");" << std::endl;
      else o << "%" << get_id() << "=FEM_TRIANGLE3(NULL,$" << maillage->get_id() << ",$" << tab[0]->get_id() << ",$" << tab[1]->get_id() << ",$" << tab[2]->get_id()<< ");" << std::endl;
    else
      if (get_lien_topologie()!=NULL) o << "%" << get_id() << "=FEM_TRIANGLE3($"<< get_lien_topologie()->get_id() << ",NULL,$" << tab[0]->get_id() << ",$" << tab[1]->get_id() << ",$" << tab[2]->get_id()<<    ");" << std::endl;
      else o << "%" << get_id() << "=FEM_TRIANGLE3(NULL,NULL,$" << tab[0]->get_id() << ",$" << tab[1]->get_id() << ",$" << tab[2]->get_id()<< ");" << std::endl;
    
}       

int FEM_TRIANGLE3::nb_fonction_interpolation(void)
{
    return 3;
}

double FEM_TRIANGLE3::get_fonction_interpolation(int num,double *uv)
{
    double val;
    switch (num)
    {
    case 1:
        val=1-uv[0]-uv[1];
        break;
    case 2:
        val=uv[0];
        break;
    case 3:
        val=uv[1];
        break;
    
    }
    return val;
}

double FEM_TRIANGLE3::get_fonction_derive_interpolation(int num,int num_variable,double *uv)
{
    double val;

    switch (num)
    {
    case 1:
        switch (num_variable)
        {
        case 1:
            val=-1;
            break;
        case 2:
            val=-1;
            break;
        } break;
    case 2:
        switch (num_variable)
        {
        case 1:
            val=1;
            break;
        case 2:
            val=0.;
            break;
        }break;
    case 3:
        switch (num_variable)
        {
        case 1:
            val=0.;
            break;
        case 2:
            val=1.;
            break;
        }break;
    
    }
    return val;
}

double FEM_TRIANGLE3::get_jacobien(double* jac,double *uv,int& li,int& col,double unite)
{   
    // Il faut exprimer les coordonnées x,y,z du repère 3D initial Ri dans un repère 2D final Rf de coordonnées
    // x',y',z' où z' est nul afin de pouvoir utiliser les fonctions d'interpolation de l'élément de référence 2D
    FEM_NOEUD* n1=get_fem_noeud(0);
    FEM_NOEUD* n2=get_fem_noeud(1);
    FEM_NOEUD* n3=get_fem_noeud(2);
    
    double xyzn1[3];
    double xyzn2[3];
    double xyzn3[3];
    
    xyzn1[0]=n1->get_x();
    xyzn1[1]=n1->get_y();
    xyzn1[2]=n1->get_z();
    
    xyzn2[0]=n2->get_x();
    xyzn2[1]=n2->get_y();
    xyzn2[2]=n2->get_z();
    
    xyzn3[0]=n3->get_x();
    xyzn3[1]=n3->get_y();
    xyzn3[2]=n3->get_z();
    
   // Position des noeuds dans le repère initial (i)
   OT_VECTEUR_3D n1_i(xyzn1);
   OT_VECTEUR_3D n2_i(xyzn2);
   OT_VECTEUR_3D n3_i(xyzn3);
   OT_MATRICE_3D pos_i(n1_i,n2_i,n3_i);
    
    // Le nouveau repère est défini par les vecteurs vec1 et vec2 et son origine est le noeud n1
    OT_VECTEUR_3D vec1(xyzn1,xyzn2);
    OT_VECTEUR_3D vec3(xyzn1,xyzn3);
    vec1.norme();
    OT_VECTEUR_3D n=vec1&vec3;          // Normale au triangle au noeud n1
    n.norme();
    OT_VECTEUR_3D vec2=n&vec1;          // Vecteur normal à vec1 dans le plan du triangle
    vec2.norme();
    OT_MATRICE_3D pif(vec1,vec2,n);     // Matrice de passage du repère Ri au repère Rf
    
    // Position des noeuds dans le repère final (f)
    OT_MATRICE_3D pif_inverse=pif.inverse();
    OT_MATRICE_3D pos_f=pif_inverse*pos_i;
    OT_VECTEUR_3D vecn1=pos_f.get_vecteur1();
    pos_f.change_vecteur1(pos_f.get_vecteur1()-vecn1);  // Étant donné que le noeud 1 est l'origine du repère final, on soustrait ses coordonnées à chaque noeud
    pos_f.change_vecteur2(pos_f.get_vecteur2()-vecn1);
    pos_f.change_vecteur3(pos_f.get_vecteur3()-vecn1);
    
    // Jacobien dans le repère final
    OT_VECTEUR_3D dNi_de(get_fonction_derive_interpolation(1,1,uv),get_fonction_derive_interpolation(2,1,uv),get_fonction_derive_interpolation(3,1,uv));
    OT_VECTEUR_3D dNi_dn(get_fonction_derive_interpolation(1,2,uv),get_fonction_derive_interpolation(2,2,uv),get_fonction_derive_interpolation(3,2,uv));
    OT_VECTEUR_3D dNi_dZ(0.,0.,0.);
    
    OT_MATRICE_3D dNi(dNi_de,dNi_dn,dNi_dZ);
    OT_MATRICE_3D dNi_transpose=dNi.transpose();
    OT_MATRICE_3D pos_f_transpose=pos_f.transpose();
    OT_MATRICE_3D J_f=dNi_transpose*pos_f_transpose;
    
    // Jacobien dans le repère initial
    OT_MATRICE_3D J_i=pif*J_f;
    
    jac[0]=J_i(0,0);
    jac[1]=J_i(1,0);
    jac[2]=J_i(2,0);
    
    jac[3]=J_i(0,1);
    jac[4]=J_i(1,1);
    jac[5]=J_i(2,1);
    
    jac[6]=0.;
    jac[7]=0.;
    jac[8]=0.;
}


bool FEM_TRIANGLE3::valide_parametre_EF(double* uvw)
{
if (uvw[0]>=-1e-10)
   if (uvw[1]>=-1e-10)
     if (uvw[0]+uvw[1]<=1.+1e-10)  
          return true;
return false;
}

void FEM_TRIANGLE3::get_inverse_jacob(double* j,double *uv,double unite)
{
    // Il faut exprimer les coordonnées x,y,z du repère 3D initial Ri dans un repère 2D final Rf de coordonnées
    // x',y',z' où z' est nul afin de pouvoir utiliser les fonctions d'interpolation de l'élément de référence
    FEM_NOEUD* n1=get_fem_noeud(0);
    FEM_NOEUD* n2=get_fem_noeud(1);
    FEM_NOEUD* n3=get_fem_noeud(2);
    
    double xyzn1[3];
    double xyzn2[3];
    double xyzn3[3];
    
    xyzn1[0]=n1->get_x();
    xyzn1[1]=n1->get_y();
    xyzn1[2]=n1->get_z();
    
    xyzn2[0]=n2->get_x();
    xyzn2[1]=n2->get_y();
    xyzn2[2]=n2->get_z();
    
    xyzn3[0]=n3->get_x();
    xyzn3[1]=n3->get_y();
    xyzn3[2]=n3->get_z();
        
    // Position des noeuds dans le repère initial (i)
    OT_VECTEUR_3D n1_i(xyzn1);
    OT_VECTEUR_3D n2_i(xyzn2);
    OT_VECTEUR_3D n3_i(xyzn3);
    OT_MATRICE_3D pos_i(n1_i,n2_i,n3_i);
    
    // Le nouveau repère est défini par les vecteurs vec1 et vec2 et son origine est le noeud n1
    OT_VECTEUR_3D vec1(xyzn1,xyzn2);
    OT_VECTEUR_3D vec3(xyzn1,xyzn3);
    vec1.norme();
    OT_VECTEUR_3D n=vec1&vec3;          // Normale au triangle au noeud n1
    n.norme();
    OT_VECTEUR_3D vec2=n&vec1;          // Vecteur normal à vec1 dans le plan du triangle
    vec2.norme();
    OT_MATRICE_3D pif(vec1,vec2,n);     // Matrice de passage du repère Ri au repère Rf
    
    // Position des noeuds dans le repère final (f)
    OT_MATRICE_3D pif_inverse=pif.inverse();
    OT_MATRICE_3D pos_f=pif_inverse*pos_i;
    OT_VECTEUR_3D vecn1=pos_f.get_vecteur1();
    pos_f.change_vecteur1(pos_f.get_vecteur1()-vecn1);  // Étant donné que le noeud 1 est l'origine du repère final, on soustrait ses coordonnées à chaque noeud
    pos_f.change_vecteur2(pos_f.get_vecteur2()-vecn1);
    pos_f.change_vecteur3(pos_f.get_vecteur3()-vecn1);
    
    // Déterminant du jacobien 2D dans le repère final
    OT_VECTEUR_3D dNi_de(get_fonction_derive_interpolation(1,1,uv),get_fonction_derive_interpolation(2,1,uv),get_fonction_derive_interpolation(3,1,uv));
    OT_VECTEUR_3D dNi_dn(get_fonction_derive_interpolation(1,2,uv),get_fonction_derive_interpolation(2,2,uv),get_fonction_derive_interpolation(3,2,uv));
    OT_VECTEUR_3D dNi_dZ(0.,0.,0.);
    
    OT_MATRICE_3D dNi(dNi_de,dNi_dn,dNi_dZ);
    OT_MATRICE_3D dNi_transpose=dNi.transpose();
    OT_MATRICE_3D pos_f_transpose=pos_f.transpose();
    OT_MATRICE_3D J=dNi_transpose*pos_f_transpose;
    
    double J11=J(0,0);
    double J21=J(1,0);
    double J12=J(0,1);
    double J22=J(1,1);
    double det_J=J11*J22-J12*J21;
    
    // Jacobien inverse 3D dans le repère final
    OT_VECTEUR_3D col1_jf_3d(J22/det_J,-J21/det_J,0.);
    OT_VECTEUR_3D col2_jf_3d(-J12/det_J,J11/det_J,0.);
    OT_VECTEUR_3D col3_jf_3d(0.,0.,0.);
    OT_MATRICE_3D j_f_3d(col1_jf_3d,col2_jf_3d,col3_jf_3d);
    
    // Jacobien inverse 3D dans le repère initial
    OT_MATRICE_3D j_i=pif*j_f_3d;
    
    j[0]=j_i(0,0);
    j[1]=j_i(1,0);
    j[2]=j_i(2,0);
    
    j[3]=j_i(0,1);
    j[4]=j_i(1,1);
    j[5]=j_i(2,1);
    
    j[6]=j_i(0,2);
    j[7]=j_i(1,2);
    j[8]=j_i(2,2);
}
Revision:	635
Committed:	Fri Jan 16 22:48:47 2015 UTC (10 years, 3 months ago) by francois
File size:	14223 byte(s)
Log Message:	Generalisation de la projection d'une solution sur un maillage avec gestion d'erreur
#	User	Rev	Content
1	francois	283	//------------------------------------------------------------
2			//------------------------------------------------------------
3			// MAGiC
4			// Jean Christophe Cuilli�re et Vincent FRANCOIS
5			// D�partement de G�nie M�canique - UQTR
6			//------------------------------------------------------------
7			// Le projet MAGIC est un projet de recherche du d�partement
8			// de g�nie m�canique de l'Universit� du Qu�bec �
9			// Trois Rivi�res
10			// Les librairies ne peuvent �tre utilis�es sans l'accord
11			// des auteurs (contact : francois@uqtr.ca)
12			//------------------------------------------------------------
13			//------------------------------------------------------------
14			//
15			// fem_triangle3.cpp
16			//
17			//------------------------------------------------------------
18			//------------------------------------------------------------
19			// COPYRIGHT 2000
20			// Version du 02/03/2006 � 11H22
21			//------------------------------------------------------------
22			//------------------------------------------------------------
23
24
25			#include "gestionversion.h"
26			#include "fem_triangle3.h"
27			#include "fem_maillage.h"
28			#include "fem_noeud.h"
29			#include "mg_element_maillage.h"
30	gervaislavoie	382	#include "math.h"
31	francois	283
32
33	francois	309	FEM_TRIANGLE3::FEM_TRIANGLE3(unsigned long num,class MG_ELEMENT_MAILLAGE* mai,class FEM_NOEUD** tabnoeud):FEM_ELEMENT2(num,mai),FEM_TEMPLATE_ELEMENT<3>(tabnoeud)
34	francois	283	{
35			if (liaison_topologique!=NULL) if (liaison_topologique->get_dimension()==2) liaison_topologique->get_lien_fem_maillage()->ajouter(this);
36	francois	309	tab[0]->get_lien_element2()->ajouter(this);
37			tab[1]->get_lien_element2()->ajouter(this);
38			tab[2]->get_lien_element2()->ajouter(this);
39			get_fem_noeudpetitid()->get_lien_petit_element2()->ajouter(this);
40	francois	283	}
41
42	francois	309	FEM_TRIANGLE3::FEM_TRIANGLE3(class MG_ELEMENT_MAILLAGE* mai,FEM_NOEUD** tabnoeud):FEM_ELEMENT2(mai),FEM_TEMPLATE_ELEMENT<3>(tabnoeud)
43	francois	283	{
44			if (liaison_topologique!=NULL) if (liaison_topologique->get_dimension()==2) liaison_topologique->get_lien_fem_maillage()->ajouter(this);
45	francois	309	tab[0]->get_lien_element2()->ajouter(this);
46			tab[1]->get_lien_element2()->ajouter(this);
47			tab[2]->get_lien_element2()->ajouter(this);
48			get_fem_noeudpetitid()->get_lien_petit_element2()->ajouter(this);
49	francois	283	}
50	francois	378	FEM_TRIANGLE3::FEM_TRIANGLE3(unsigned long num,class MG_ELEMENT_TOPOLOGIQUE* topo,class FEM_NOEUD** tabnoeud):FEM_ELEMENT2(num,topo),FEM_TEMPLATE_ELEMENT<3>(tabnoeud)
51			{
52			if (liaison_topologique!=NULL) if (liaison_topologique->get_dimension()==2) liaison_topologique->get_lien_fem_maillage()->ajouter(this);
53			tab[0]->get_lien_element2()->ajouter(this);
54			tab[1]->get_lien_element2()->ajouter(this);
55			tab[2]->get_lien_element2()->ajouter(this);
56			get_fem_noeudpetitid()->get_lien_petit_element2()->ajouter(this);
57			}
58	francois	283
59	francois	378	FEM_TRIANGLE3::FEM_TRIANGLE3(class MG_ELEMENT_TOPOLOGIQUE* topo,FEM_NOEUD** tabnoeud):FEM_ELEMENT2(topo),FEM_TEMPLATE_ELEMENT<3>(tabnoeud)
60			{
61			if (liaison_topologique!=NULL) if (liaison_topologique->get_dimension()==2) liaison_topologique->get_lien_fem_maillage()->ajouter(this);
62			tab[0]->get_lien_element2()->ajouter(this);
63			tab[1]->get_lien_element2()->ajouter(this);
64			tab[2]->get_lien_element2()->ajouter(this);
65			get_fem_noeudpetitid()->get_lien_petit_element2()->ajouter(this);
66			}
67			FEM_TRIANGLE3::FEM_TRIANGLE3(unsigned long num,class MG_ELEMENT_TOPOLOGIQUE* topo,class MG_ELEMENT_MAILLAGE* mai,class FEM_NOEUD** tabnoeud):FEM_ELEMENT2(num,topo,mai),FEM_TEMPLATE_ELEMENT<3>(tabnoeud)
68			{
69			if (liaison_topologique!=NULL) if (liaison_topologique->get_dimension()==2) liaison_topologique->get_lien_fem_maillage()->ajouter(this);
70			tab[0]->get_lien_element2()->ajouter(this);
71			tab[1]->get_lien_element2()->ajouter(this);
72			tab[2]->get_lien_element2()->ajouter(this);
73			get_fem_noeudpetitid()->get_lien_petit_element2()->ajouter(this);
74			}
75
76			FEM_TRIANGLE3::FEM_TRIANGLE3(class MG_ELEMENT_TOPOLOGIQUE* topo,class MG_ELEMENT_MAILLAGE* mai,FEM_NOEUD** tabnoeud):FEM_ELEMENT2(topo,mai),FEM_TEMPLATE_ELEMENT<3>(tabnoeud)
77			{
78			if (liaison_topologique!=NULL) if (liaison_topologique->get_dimension()==2) liaison_topologique->get_lien_fem_maillage()->ajouter(this);
79			tab[0]->get_lien_element2()->ajouter(this);
80			tab[1]->get_lien_element2()->ajouter(this);
81			tab[2]->get_lien_element2()->ajouter(this);
82			get_fem_noeudpetitid()->get_lien_petit_element2()->ajouter(this);
83			}
84	francois	309	FEM_TRIANGLE3::FEM_TRIANGLE3(FEM_TRIANGLE3& mdd):FEM_ELEMENT2(mdd),FEM_TEMPLATE_ELEMENT<3>(mdd)
85	francois	283	{
86			if (liaison_topologique!=NULL) if (liaison_topologique->get_dimension()==2) liaison_topologique->get_lien_fem_maillage()->ajouter(this);
87	francois	309	tab[0]->get_lien_element2()->ajouter(this);
88			tab[1]->get_lien_element2()->ajouter(this);
89			tab[2]->get_lien_element2()->ajouter(this);
90			get_fem_noeudpetitid()->get_lien_petit_element2()->ajouter(this);
91	francois	283	}
92			FEM_TRIANGLE3::~FEM_TRIANGLE3()
93			{
94			if (liaison_topologique==NULL) return;
95			if (liaison_topologique->get_dimension()==0) liaison_topologique->get_lien_fem_maillage()->supprimer(this);
96	francois	309	tab[0]->get_lien_element2()->supprimer(this);
97			tab[1]->get_lien_element2()->supprimer(this);
98			tab[2]->get_lien_element2()->supprimer(this);
99			get_fem_noeudpetitid()->get_lien_petit_element2()->supprimer(this);
100	francois	283	}
101
102
103			FEM_ELEMENT_MAILLAGE* FEM_TRIANGLE3::dupliquer(FEM_MAILLAGE *femmai,long decalage)
104			{
105			FEM_NOEUD* tabnoeud[3];
106			tabnoeud[0]=femmai->get_fem_noeudid(tab[0]->get_id()+decalage);
107			tabnoeud[1]=femmai->get_fem_noeudid(tab[1]->get_id()+decalage);
108			tabnoeud[2]=femmai->get_fem_noeudid(tab[2]->get_id()+decalage);
109			FEM_TRIANGLE3* tri=new FEM_TRIANGLE3(get_id()+decalage,maillage,tabnoeud);
110	francois	309	femmai->ajouter_fem_element2(tri);
111	francois	283	return tri;
112			}
113
114
115
116
117
118			int FEM_TRIANGLE3::get_type_entite(void)
119			{
120			return IDFEM_TRIANGLE3;
121			}
122
123			int FEM_TRIANGLE3::get_dimension(void)
124			{
125			return 2;
126			}
127
128
129
130			int FEM_TRIANGLE3::get_nb_fem_noeud(void)
131			{
132			return FEM_TEMPLATE_ELEMENT<3>::get_nb_fem_noeud();
133			}
134
135			FEM_NOEUD* FEM_TRIANGLE3::get_fem_noeud(int num)
136			{
137			return FEM_TEMPLATE_ELEMENT<3>::get_fem_noeud(num);
138			}
139
140			void FEM_TRIANGLE3::change_noeud(int num,FEM_NOEUD* noeud)
141			{
142			FEM_TEMPLATE_ELEMENT<3>::change_noeud(num,noeud);
143			}
144
145			BOITE_3D& FEM_TRIANGLE3::get_boite_3D(void)
146			{
147			return FEM_TEMPLATE_ELEMENT<3>::get_boite_3D();
148			}
149
150
151			void FEM_TRIANGLE3::enregistrer(std::ostream& o)
152			{
153	francois	378	if (maillage!=NULL)
154			if (get_lien_topologie()!=NULL) o << "%" << get_id() << "=FEM_TRIANGLE3($"<< get_lien_topologie()->get_id() << ",$" << maillage->get_id() << ",$" << tab[0]->get_id() << ",$" << tab[1]->get_id() << ",$" << tab[2]->get_id()<< ");" << std::endl;
155			else o << "%" << get_id() << "=FEM_TRIANGLE3(NULL,$" << maillage->get_id() << ",$" << tab[0]->get_id() << ",$" << tab[1]->get_id() << ",$" << tab[2]->get_id()<< ");" << std::endl;
156			else
157			if (get_lien_topologie()!=NULL) o << "%" << get_id() << "=FEM_TRIANGLE3($"<< get_lien_topologie()->get_id() << ",NULL,$" << tab[0]->get_id() << ",$" << tab[1]->get_id() << ",$" << tab[2]->get_id()<< ");" << std::endl;
158			else o << "%" << get_id() << "=FEM_TRIANGLE3(NULL,NULL,$" << tab[0]->get_id() << ",$" << tab[1]->get_id() << ",$" << tab[2]->get_id()<< ");" << std::endl;
159
160			}
161	francois	310
162			int FEM_TRIANGLE3::nb_fonction_interpolation(void)
163			{
164			return 3;
165			}
166
167			double FEM_TRIANGLE3::get_fonction_interpolation(int num,double *uv)
168			{
169			double val;
170			switch (num)
171			{
172			case 1:
173			val=1-uv[0]-uv[1];
174			break;
175			case 2:
176			val=uv[0];
177			break;
178			case 3:
179			val=uv[1];
180			break;
181
182			}
183			return val;
184			}
185
186			double FEM_TRIANGLE3::get_fonction_derive_interpolation(int num,int num_variable,double *uv)
187			{
188			double val;
189
190			switch (num)
191			{
192			case 1:
193			switch (num_variable)
194			{
195			case 1:
196			val=-1;
197			break;
198			case 2:
199			val=-1;
200			break;
201			} break;
202			case 2:
203			switch (num_variable)
204			{
205			case 1:
206			val=1;
207			break;
208			case 2:
209			val=0.;
210			break;
211			}break;
212			case 3:
213			switch (num_variable)
214			{
215			case 1:
216			val=0.;
217			break;
218			case 2:
219			val=1.;
220			break;
221			}break;
222
223			}
224			return val;
225			}
226
227			double FEM_TRIANGLE3::get_jacobien(double* jac,double *uv,int& li,int& col,double unite)
228	gervaislavoie	382	{
229			// Il faut exprimer les coordonnées x,y,z du repère 3D initial Ri dans un repère 2D final Rf de coordonnées
230	gervaislavoie	389	// x',y',z' où z' est nul afin de pouvoir utiliser les fonctions d'interpolation de l'élément de référence 2D
231	gervaislavoie	382	FEM_NOEUD* n1=get_fem_noeud(0);
232			FEM_NOEUD* n2=get_fem_noeud(1);
233			FEM_NOEUD* n3=get_fem_noeud(2);
234
235	gervaislavoie	389	double xyzn1[3];
236			double xyzn2[3];
237			double xyzn3[3];
238
239	gervaislavoie	390	xyzn1[0]=n1->get_x();
240			xyzn1[1]=n1->get_y();
241			xyzn1[2]=n1->get_z();
242	gervaislavoie	389
243	gervaislavoie	390	xyzn2[0]=n2->get_x();
244			xyzn2[1]=n2->get_y();
245			xyzn2[2]=n2->get_z();
246	gervaislavoie	389
247	gervaislavoie	390	xyzn3[0]=n3->get_x();
248			xyzn3[1]=n3->get_y();
249			xyzn3[2]=n3->get_z();
250	gervaislavoie	389
251	gervaislavoie	392	// Position des noeuds dans le repère initial (i)
252			OT_VECTEUR_3D n1_i(xyzn1);
253			OT_VECTEUR_3D n2_i(xyzn2);
254			OT_VECTEUR_3D n3_i(xyzn3);
255			OT_MATRICE_3D pos_i(n1_i,n2_i,n3_i);
256
257	gervaislavoie	382	// Le nouveau repère est défini par les vecteurs vec1 et vec2 et son origine est le noeud n1
258			OT_VECTEUR_3D vec1(xyzn1,xyzn2);
259			OT_VECTEUR_3D vec3(xyzn1,xyzn3);
260			vec1.norme();
261			OT_VECTEUR_3D n=vec1&vec3; // Normale au triangle au noeud n1
262			n.norme();
263			OT_VECTEUR_3D vec2=n&vec1; // Vecteur normal à vec1 dans le plan du triangle
264			vec2.norme();
265			OT_MATRICE_3D pif(vec1,vec2,n); // Matrice de passage du repère Ri au repère Rf
266
267	gervaislavoie	392	// Position des noeuds dans le repère final (f)
268	gervaislavoie	382	OT_MATRICE_3D pif_inverse=pif.inverse();
269	gervaislavoie	392	OT_MATRICE_3D pos_f=pif_inverse*pos_i;
270			OT_VECTEUR_3D vecn1=pos_f.get_vecteur1();
271			pos_f.change_vecteur1(pos_f.get_vecteur1()-vecn1); // Étant donné que le noeud 1 est l'origine du repère final, on soustrait ses coordonnées à chaque noeud
272			pos_f.change_vecteur2(pos_f.get_vecteur2()-vecn1);
273			pos_f.change_vecteur3(pos_f.get_vecteur3()-vecn1);
274	gervaislavoie	382
275	gervaislavoie	389	// Jacobien dans le repère final
276			OT_VECTEUR_3D dNi_de(get_fonction_derive_interpolation(1,1,uv),get_fonction_derive_interpolation(2,1,uv),get_fonction_derive_interpolation(3,1,uv));
277			OT_VECTEUR_3D dNi_dn(get_fonction_derive_interpolation(1,2,uv),get_fonction_derive_interpolation(2,2,uv),get_fonction_derive_interpolation(3,2,uv));
278			OT_VECTEUR_3D dNi_dZ(0.,0.,0.);
279	gervaislavoie	382
280	gervaislavoie	389	OT_MATRICE_3D dNi(dNi_de,dNi_dn,dNi_dZ);
281			OT_MATRICE_3D dNi_transpose=dNi.transpose();
282	gervaislavoie	392	OT_MATRICE_3D pos_f_transpose=pos_f.transpose();
283			OT_MATRICE_3D J_f=dNi_transpose*pos_f_transpose;
284	gervaislavoie	382
285	gervaislavoie	389	// Jacobien dans le repère initial
286			OT_MATRICE_3D J_i=pif*J_f;
287	gervaislavoie	382
288			jac[0]=J_i(0,0);
289			jac[1]=J_i(1,0);
290			jac[2]=J_i(2,0);
291
292			jac[3]=J_i(0,1);
293			jac[4]=J_i(1,1);
294			jac[5]=J_i(2,1);
295
296			jac[6]=0.;
297			jac[7]=0.;
298			jac[8]=0.;
299	francois	310	}
300
301	francois	635
302			bool FEM_TRIANGLE3::valide_parametre_EF(double* uvw)
303			{
304			if (uvw[0]>=-1e-10)
305			if (uvw[1]>=-1e-10)
306			if (uvw[0]+uvw[1]<=1.+1e-10)
307			return true;
308			return false;
309			}
310
311	francois	310	void FEM_TRIANGLE3::get_inverse_jacob(double* j,double *uv,double unite)
312			{
313	gervaislavoie	382	// Il faut exprimer les coordonnées x,y,z du repère 3D initial Ri dans un repère 2D final Rf de coordonnées
314	gervaislavoie	389	// x',y',z' où z' est nul afin de pouvoir utiliser les fonctions d'interpolation de l'élément de référence
315	gervaislavoie	382	FEM_NOEUD* n1=get_fem_noeud(0);
316			FEM_NOEUD* n2=get_fem_noeud(1);
317			FEM_NOEUD* n3=get_fem_noeud(2);
318
319	gervaislavoie	389	double xyzn1[3];
320			double xyzn2[3];
321			double xyzn3[3];
322	gervaislavoie	382
323	gervaislavoie	390	xyzn1[0]=n1->get_x();
324			xyzn1[1]=n1->get_y();
325			xyzn1[2]=n1->get_z();
326	gervaislavoie	382
327	gervaislavoie	390	xyzn2[0]=n2->get_x();
328			xyzn2[1]=n2->get_y();
329			xyzn2[2]=n2->get_z();
330	gervaislavoie	382
331	gervaislavoie	390	xyzn3[0]=n3->get_x();
332			xyzn3[1]=n3->get_y();
333			xyzn3[2]=n3->get_z();
334	gervaislavoie	392
335			// Position des noeuds dans le repère initial (i)
336			OT_VECTEUR_3D n1_i(xyzn1);
337			OT_VECTEUR_3D n2_i(xyzn2);
338			OT_VECTEUR_3D n3_i(xyzn3);
339			OT_MATRICE_3D pos_i(n1_i,n2_i,n3_i);
340	gervaislavoie	382
341	gervaislavoie	392	// Le nouveau repère est défini par les vecteurs vec1 et vec2 et son origine est le noeud n1
342	gervaislavoie	382	OT_VECTEUR_3D vec1(xyzn1,xyzn2);
343			OT_VECTEUR_3D vec3(xyzn1,xyzn3);
344			vec1.norme();
345			OT_VECTEUR_3D n=vec1&vec3; // Normale au triangle au noeud n1
346			n.norme();
347			OT_VECTEUR_3D vec2=n&vec1; // Vecteur normal à vec1 dans le plan du triangle
348			vec2.norme();
349			OT_MATRICE_3D pif(vec1,vec2,n); // Matrice de passage du repère Ri au repère Rf
350
351	gervaislavoie	392	// Position des noeuds dans le repère final (f)
352	gervaislavoie	382	OT_MATRICE_3D pif_inverse=pif.inverse();
353	gervaislavoie	392	OT_MATRICE_3D pos_f=pif_inverse*pos_i;
354			OT_VECTEUR_3D vecn1=pos_f.get_vecteur1();
355			pos_f.change_vecteur1(pos_f.get_vecteur1()-vecn1); // Étant donné que le noeud 1 est l'origine du repère final, on soustrait ses coordonnées à chaque noeud
356			pos_f.change_vecteur2(pos_f.get_vecteur2()-vecn1);
357			pos_f.change_vecteur3(pos_f.get_vecteur3()-vecn1);
358	gervaislavoie	382
359	gervaislavoie	389	// Déterminant du jacobien 2D dans le repère final
360			OT_VECTEUR_3D dNi_de(get_fonction_derive_interpolation(1,1,uv),get_fonction_derive_interpolation(2,1,uv),get_fonction_derive_interpolation(3,1,uv));
361			OT_VECTEUR_3D dNi_dn(get_fonction_derive_interpolation(1,2,uv),get_fonction_derive_interpolation(2,2,uv),get_fonction_derive_interpolation(3,2,uv));
362			OT_VECTEUR_3D dNi_dZ(0.,0.,0.);
363	gervaislavoie	382
364	gervaislavoie	389	OT_MATRICE_3D dNi(dNi_de,dNi_dn,dNi_dZ);
365			OT_MATRICE_3D dNi_transpose=dNi.transpose();
366	gervaislavoie	392	OT_MATRICE_3D pos_f_transpose=pos_f.transpose();
367			OT_MATRICE_3D J=dNi_transpose*pos_f_transpose;
368	gervaislavoie	382
369	gervaislavoie	389	double J11=J(0,0);
370			double J21=J(1,0);
371			double J12=J(0,1);
372			double J22=J(1,1);
373	gervaislavoie	382	double det_J=J11J22-J12J21;
374
375	gervaislavoie	389	// Jacobien inverse 3D dans le repère final
376	gervaislavoie	382	OT_VECTEUR_3D col1_jf_3d(J22/det_J,-J21/det_J,0.);
377			OT_VECTEUR_3D col2_jf_3d(-J12/det_J,J11/det_J,0.);
378			OT_VECTEUR_3D col3_jf_3d(0.,0.,0.);
379	gervaislavoie	389	OT_MATRICE_3D j_f_3d(col1_jf_3d,col2_jf_3d,col3_jf_3d);
380	gervaislavoie	382
381	gervaislavoie	389	// Jacobien inverse 3D dans le repère initial
382			OT_MATRICE_3D j_i=pif*j_f_3d;
383	gervaislavoie	382
384			j[0]=j_i(0,0);
385			j[1]=j_i(1,0);
386			j[2]=j_i(2,0);
387
388			j[3]=j_i(0,1);
389			j[4]=j_i(1,1);
390			j[5]=j_i(2,1);
391
392			j[6]=j_i(0,2);
393			j[7]=j_i(1,2);
394			j[8]=j_i(2,2);
395	francois	310	}