fem_maillage_outils.cpp

Aller à la documentation de ce fichier.

//####//------------------------------------------------------------
//####//------------------------------------------------------------
//####//  MAGiC
//####//  Jean Christophe Cuilliere et Vincent FRANCOIS
//####//  Departement de Genie Mecanique - UQTR
//####//------------------------------------------------------------
//####//  MAGIC est un projet de recherche de l equipe ERICCA
//####//  du departement de genie mecanique de l Universite du Quebec a Trois Rivieres
//####//  http://www.uqtr.ca/ericca
//####//  http://www.uqtr.ca/
//####//------------------------------------------------------------
//####//------------------------------------------------------------
//####//
//####//       fem_maillage_outils.cpp
//####//
//####//------------------------------------------------------------
//####//------------------------------------------------------------
//####//    COPYRIGHT 2000-2024
//####//  jeu 13 jun 2024 11:58:55 EDT
//####//------------------------------------------------------------
//####//------------------------------------------------------------
#include "gestionversion.h"
#include "fem_maillage_outils.h"
#include "mg_element_maillage.h"
#include "fem_tetra4.h"
#include "fem_element2.h"
#include "fem_noeud.h"
#include "xfem_tetra4.h"
#include "ot_mathematique.h"
#include "tpl_fonctions_generiques.h"
#include "fem_solution.h"
#include "fem_maillage.h"
#include "mg_gestionnaire.h"
#include <tpl_map_entite.h>
#include <tpl_liste_entite.h>
#include <math.h>
#include <stdio.h>
#include <string.h>
#include "ot_tenseur.h"
#include "fem_triangle3.h"
#include "fem_template_element.h"
#include "tpl_octree.h"
 
 
FEM_MAILLAGE_OUTILS::FEM_MAILLAGE_OUTILS()
{
}
 
FEM_MAILLAGE_OUTILS::FEM_MAILLAGE_OUTILS(FEM_MAILLAGE_OUTILS &mdd)
{
}
 
FEM_MAILLAGE_OUTILS::~FEM_MAILLAGE_OUTILS()
{
}
 
 
 
 
 
 
int FEM_MAILLAGE_OUTILS::projeteestdansletriangle(FEM_ELEMENT2 *tri,double x,double y, double z)
{
double *xyz1,*xyz2,*xyz3;
if (tri->get_nb_fem_noeud()==3)
        {
        xyz1=tri->get_fem_noeud(0)->get_coord();
        xyz2=tri->get_fem_noeud(1)->get_coord();
        xyz3=tri->get_fem_noeud(2)->get_coord();
        }
if (tri->get_nb_fem_noeud()==6)
        {
        xyz1=tri->get_fem_noeud(0)->get_coord();
        xyz2=tri->get_fem_noeud(2)->get_coord();
        xyz3=tri->get_fem_noeud(4)->get_coord();
        }
return OPERATEUR::projeteestdansletriangle(xyz1,xyz2,xyz3,x,y,z); 
}
 
int FEM_MAILLAGE_OUTILS::estdansletriangle(FEM_ELEMENT2 *tri,double x,double y, double z)
{
double *xyz1,*xyz2,*xyz3;
if (tri->get_nb_fem_noeud()==3)
        {
        xyz1=tri->get_fem_noeud(0)->get_coord();
        xyz2=tri->get_fem_noeud(1)->get_coord();
        xyz3=tri->get_fem_noeud(2)->get_coord();
        }
if (tri->get_nb_fem_noeud()==6)
        {
        xyz1=tri->get_fem_noeud(0)->get_coord();
        xyz2=tri->get_fem_noeud(2)->get_coord();
        xyz3=tri->get_fem_noeud(4)->get_coord();
        }
return OPERATEUR::estdansletriangle(xyz1,xyz2,xyz3,x,y,z);
}
 
int FEM_MAILLAGE_OUTILS::estdansletetra(FEM_ELEMENT3 *tet,double x,double y, double z)
{
double *xyz1,*xyz2,*xyz3,*xyz4;
if (tet->get_nb_fem_noeud()==4)
        {
        xyz1=tet->get_fem_noeud(0)->get_coord();
        xyz2=tet->get_fem_noeud(1)->get_coord();
        xyz3=tet->get_fem_noeud(2)->get_coord();
        xyz4=tet->get_fem_noeud(3)->get_coord();
        }
if (tet->get_nb_fem_noeud()==10)
        {
        xyz1=tet->get_fem_noeud(0)->get_coord();
        xyz2=tet->get_fem_noeud(2)->get_coord();
        xyz3=tet->get_fem_noeud(4)->get_coord();
        xyz4=tet->get_fem_noeud(9)->get_coord();
        }
return OPERATEUR::estdansletetra(xyz1,xyz2,xyz3,xyz4,x,y,z);
}
 
int FEM_MAILLAGE_OUTILS::projeteestdansletriangle(XFEM_ELEMENT2 *tri,double x,double y, double z)
{
double *xyz1,*xyz2,*xyz3;
if (tri->get_nb_fem_noeud()==3)
        {
        xyz1=tri->get_fem_noeud(0)->get_coord();
        xyz2=tri->get_fem_noeud(1)->get_coord();
        xyz3=tri->get_fem_noeud(2)->get_coord();
        }
if (tri->get_nb_fem_noeud()==6)
        {
        xyz1=tri->get_fem_noeud(0)->get_coord();
        xyz2=tri->get_fem_noeud(2)->get_coord();
        xyz3=tri->get_fem_noeud(4)->get_coord();
        }
return OPERATEUR::projeteestdansletriangle(xyz1,xyz2,xyz3,x,y,z); 
}
 
 
int FEM_MAILLAGE_OUTILS::estdansletriangle(XFEM_ELEMENT2 *tri,double x,double y, double z)
{
double *xyz1,*xyz2,*xyz3;
if (tri->get_nb_fem_noeud()==3)
        {
        xyz1=tri->get_fem_noeud(0)->get_coord();
        xyz2=tri->get_fem_noeud(1)->get_coord();
        xyz3=tri->get_fem_noeud(2)->get_coord();
        }
if (tri->get_nb_fem_noeud()==6)
        {
        xyz1=tri->get_fem_noeud(0)->get_coord();
        xyz2=tri->get_fem_noeud(2)->get_coord();
        xyz3=tri->get_fem_noeud(4)->get_coord();
        }
return OPERATEUR::estdansletriangle(xyz1,xyz2,xyz3,x,y,z);
}
  
 
 
 
int FEM_MAILLAGE_OUTILS::estdansletetra(XFEM_ELEMENT3 *tet,double x,double y, double z)
{
double *xyz1,*xyz2,*xyz3,*xyz4;
if (tet->get_nb_fem_noeud()==4)
        {
        xyz1=tet->get_fem_noeud(0)->get_coord();
        xyz2=tet->get_fem_noeud(1)->get_coord();
        xyz3=tet->get_fem_noeud(2)->get_coord();
        xyz4=tet->get_fem_noeud(3)->get_coord();
        }
if (tet->get_nb_fem_noeud()==10)
        {
        xyz1=tet->get_fem_noeud(0)->get_coord();
        xyz2=tet->get_fem_noeud(2)->get_coord();
        xyz3=tet->get_fem_noeud(4)->get_coord();
        xyz4=tet->get_fem_noeud(9)->get_coord();
        }
return OPERATEUR::estdansletetra(xyz1,xyz2,xyz3,xyz4,x,y,z);
  
}
 
 
int FEM_MAILLAGE_OUTILS::compare_etat_tetra(int etat,int valeur)
{
return OPERATEUR::compare_etat_tetra(etat,valeur);  
}
 
 
int FEM_MAILLAGE_OUTILS::compare_etat_triangle(int etat,int valeur)
{
return OPERATEUR::compare_etat_triangle(etat,valeur);  
}
 
 
double FEM_MAILLAGE_OUTILS::get_angle(FEM_ELEMENT2* ele1,FEM_ELEMENT2* ele2)
{
    if (ele1->get_nb_fem_noeud()==3)
    if (ele2->get_nb_fem_noeud()==3)
      {
      FEM_NOEUD* noeud1=ele1->get_fem_noeud(0);
      FEM_NOEUD* noeud2=ele1->get_fem_noeud(1);
      FEM_NOEUD* noeud3=ele1->get_fem_noeud(2);
      FEM_NOEUD* noeud4=ele2->get_fem_noeud(0);
      FEM_NOEUD* noeud5=ele2->get_fem_noeud(1);
      FEM_NOEUD* noeud6=ele2->get_fem_noeud(2);
      return get_angle_par_noeud<FEM_NOEUD*>(noeud1,noeud2,noeud3,noeud4,noeud5,noeud6);
      }
    return 0.;
}
 
 
void FEM_MAILLAGE_OUTILS::calcul_courbure(FEM_MAILLAGE* femmai, MG_GESTIONNAIRE& gest,char *nom)
{
        FEM_SOLUTION* sol=new FEM_SOLUTION(femmai,5,nom,1,"Courbure",MAGIC::ENTITE_SOLUTION::ENTITE_NOEUD,MAGIC::TYPE_SOLUTION::SCALAIRE);
        gest.ajouter_fem_solution(sol);                     
 
        sol->change_legende(0,"Gauss");  
        sol->change_legende(1,"Moyenne");  
        sol->change_legende(2,"Maximale");    
        sol->change_legende(3,"Minimale");
        sol->change_legende(4,"frontiere");
        
        int h=0;      
        LISTE_FEM_NOEUD::iterator it_no;
        for (FEM_NOEUD* noeud=femmai->get_premier_noeud(it_no);noeud!=NULL;noeud=femmai->get_suivant_noeud(it_no))
                {          
                if (noeud->get_lien_topologie()->get_dimension()==2)
                        {
                        int nb_ele2=noeud->get_lien_element2()->get_nb();     
                        TPL_MAP_ENTITE<FEM_ELEMENT2*> liste_map_element2;                    
                        TPL_MAP_ENTITE<FEM_ELEMENT2*>::ITERATEUR it;                 
                        TPL_LISTE_ENTITE<FEM_ELEMENT2*> liste_tpl_element2;
 
                        for (int k=0;k<nb_ele2;k++)
                                {   
                                FEM_ELEMENT2* tri1=noeud->get_lien_element2()->get(k); 
                                liste_map_element2.ajouter(tri1);
                                }
 
                        FEM_ELEMENT2* tri2=liste_map_element2.get_premier(it);
                        FEM_NOEUD* na=tri2->get_fem_noeud(0);
                        FEM_NOEUD* nb=tri2->get_fem_noeud(1);
                        FEM_NOEUD* nc=tri2->get_fem_noeud(2);
                        FEM_NOEUD* n1=NULL;
                        FEM_NOEUD* n2=NULL;
                        FEM_NOEUD* n3=NULL;
 
                        if (na==noeud) {n1=na;n2=nb;n3=nc;}
                        if (nb==noeud) {n1=nb;n2=nc;n3=na;}
                        if (nc==noeud) {n1=nc;n2=na;n3=nb;}
 
                        liste_tpl_element2.ajouter(tri2);
 
                        FEM_ELEMENT2* tri=NULL;
                        while (liste_map_element2.get_nb()!=liste_tpl_element2.get_nb())
                        for (tri=liste_map_element2.get_premier(it);tri;tri=liste_map_element2.get_suivant(it)) 
                                {
                                if (tri!=tri2)
                                        {
                                        FEM_NOEUD* noa=tri->get_fem_noeud(0);
                                        FEM_NOEUD* nob=tri->get_fem_noeud(1);
                                        FEM_NOEUD* noc=tri->get_fem_noeud(2);
                                        FEM_NOEUD* n4=NULL;
                                        FEM_NOEUD* n5=NULL;
                                        FEM_NOEUD* n6=NULL;
 
                                        if (noa==noeud) {n4=noa;n5=nob;n6=noc;}
                                        if (nob==noeud) {n4=nob;n5=noc;n6=noa;}
                                        if (noc==noeud) {n4=noc;n5=noa;n6=nob;}
 
                                        if (n1==n4)
                                                if (n3==n5)
                                                        {
                                                        liste_tpl_element2.ajouter(tri);
                                                        n1=n4;
                                                        n2=n5;
                                                        n3=n6;
                                                        if (liste_map_element2.get_nb()==liste_tpl_element2.get_nb()) break;
                                                        }
                                        }
                                } 
 
                        double K=0.;
                        double H=0.;
                        double kM=0.;
                        double km=0.;     
 
                        double beta=0.;
                        double surface=0.;
                        double som1=0.;   
                        double som2=0.;   
                        double som3=0.;   
                        double som4=0.; 
                        
                        double def_ang=0.;
   
                        for (int i=0;i<liste_tpl_element2.get_nb();i++)
                                {
                                FEM_ELEMENT2* tri00=NULL;
                                FEM_ELEMENT2* tri01=NULL;
 
                                if (!(i==(liste_tpl_element2.get_nb()-1))) {tri00=liste_tpl_element2.get(i);tri01=liste_tpl_element2.get(i+1);}
                                if (i==(liste_tpl_element2.get_nb()-1))    {tri00=liste_tpl_element2.get(i);tri01=tri2;}
 
                                FEM_NOEUD* noeud1=NULL;
                                FEM_NOEUD* noeud2=NULL;
                                FEM_NOEUD* noeud3=NULL;
                                FEM_NOEUD* noeud4=NULL;
                                FEM_NOEUD* noeud5=NULL;
                                FEM_NOEUD* noeud6=NULL;
 
                                FEM_NOEUD* nd1=tri00->get_fem_noeud(0);
                                FEM_NOEUD* nd2=tri00->get_fem_noeud(1);
                                FEM_NOEUD* nd3=tri00->get_fem_noeud(2);   
                                FEM_NOEUD* nd4=tri01->get_fem_noeud(0);
                                FEM_NOEUD* nd5=tri01->get_fem_noeud(1);
                                FEM_NOEUD* nd6=tri01->get_fem_noeud(2); 
    
                                if (nd1==noeud) {noeud1=nd1;noeud2=nd2;noeud3=nd3;}
                                if (nd2==noeud) {noeud1=nd2;noeud2=nd3;noeud3=nd1;}
                                if (nd3==noeud) {noeud1=nd3;noeud2=nd1;noeud3=nd2;}
                                if (nd4==noeud) {noeud4=nd4;noeud5=nd5;noeud6=nd6;}
                                if (nd5==noeud) {noeud4=nd5;noeud5=nd6;noeud6=nd4;}
                                if (nd6==noeud) {noeud4=nd6;noeud5=nd4;noeud6=nd5;}
 
                                double *xyz1=noeud1->get_coord();
                                double *xyz2=noeud2->get_coord();
                                double *xyz3=noeud3->get_coord();
                                double *xyz4=noeud4->get_coord();
                                double *xyz5=noeud5->get_coord();
                                double *xyz6=noeud6->get_coord();
 
                                OT_VECTEUR_3D vec1(xyz1,xyz2);
                                OT_VECTEUR_3D vec2(xyz1,xyz3);
                                OT_VECTEUR_3D vec3(xyz2,xyz3);
                                OT_VECTEUR_3D vec4(xyz4,xyz5);
                                OT_VECTEUR_3D vec5(xyz4,xyz6);
                                OT_VECTEUR_3D vec6(xyz5,xyz6);
            
            
                                // COURBURE GAUSSIENNE
                                // Calcul de alpha (i) de chaque triangle lié au noeud
                                double cosalpha=((vec1*vec2)/(vec1.get_longueur()*vec2.get_longueur()));
                                double alpha=acos(cosalpha);
                                double cotalpha=(cosalpha/sin(alpha));
                                som1=som1+alpha;
 
                                // Calcul Aire triangle 1
                                OT_VECTEUR_3D pvec1=vec1&vec2;
                                surface=0.5*pvec1.get_longueur();
                                som2=som2+surface;
    
                                // Calcul de I = longueur de l'arête opposée à l'angle alpha
                                double I=vec3.get_longueur();
                                som3=som3+(I*I*cotalpha);
            
                                // COURBURE MOYENNE
                                // Calcul de la longueur de l'arête commune aux deux triangles adjacents tri00 et tri01
                                double L=vec2.get_longueur();
            
                                // Calcul de beta, l'angle entre les normales de 2 triangles adjacents
                                FEM_MAILLAGE_OUTILS ot;
                                beta=M_PI-ot.get_angle(tri00,tri01);
                                som4=som4+(beta*L);
                                }    
                        
                        // Courbures de Gauss, Moyenne, Min et Max
                        K=((2.*M_PI-som1)/((0.5*som2)-(0.125*som3)));
                        H=((0.25*som4)/((0.5*som2)-(0.125*som3)));
                        if (H*H-K>0.) kM=H+sqrt((H)*(H)-(K)); else kM=H;
                        if (kM<0.000001) kM=0.;
                        if (H*H-K>0.) km=H-sqrt((H)*(H)-(K)); else km=H;
                        if (km>-0.000001) km=0.;
 
                        // Solutions
                        sol->ecrire(K,h,0);
                        sol->ecrire(H,h,1);      
                        sol->ecrire(kM,h,2);
                        sol->ecrire(km,h,3); 
                        }
/*
                // Calcul de la courbure sur les frontières
                // Les courbures de Gauss et Moyenne (et donc les Min et Max aussi) étant inexistantes aux frontières, 
                // Calladine propose une approximation de la courbure de la courbe représentée par la frontière
                if (noeud->get_lien_topologie()->get_dimension()==1)
                        {
                        FEM_ELEMENT1* seg1=noeud->get_lien_element1()->get(0); 
                        FEM_ELEMENT1* seg2=noeud->get_lien_element1()->get(1);
        
                        FEM_NOEUD* ned2=NULL;
                        FEM_NOEUD* ned3=NULL;
        
                        FEM_NOEUD* nedX=NULL;
                        FEM_NOEUD* nedY=NULL;
                        FEM_NOEUD* nedZ=NULL;
                        FEM_NOEUD* nedT=NULL;
        
                        nedX=seg1->get_fem_noeud(0);
                        nedY=seg1->get_fem_noeud(1);
                        nedZ=seg2->get_fem_noeud(0);
                        nedT=seg2->get_fem_noeud(1);
        
                        if (nedX==noeud) ned2=nedY;
                        if (nedY==noeud) ned2=nedX;
                        if (nedZ==noeud) ned3=nedT;
                        if (nedT==noeud) ned3=nedZ;
        
                        double *coord1=noeud->get_coord();
                        double *coord2=ned2->get_coord();
                        double *coord3=ned3->get_coord();
        
                        OT_VECTEUR_3D vec1(coord1,coord2);
                        OT_VECTEUR_3D vec2(coord1,coord3);
                        double cosalpha=((vec1*vec2)/(vec1.get_longueur()*vec2.get_longueur()));
                        double alpha=acos(cosalpha);
        
                        double L1=vec1.get_longueur();
                        double L2=vec2.get_longueur();
                        double L=0.;
                        if (L1>L2) L=L1;
                        else L=L2;
                        double Cf=((M_PI-alpha)/L);
          
                        // Solutions
                        sol->ecrire(Cf,h,4); 
                        }
*/
 
                h++;
                }
}
 
 
 
double FEM_MAILLAGE_OUTILS::get_gradient_element(FEM_ELEMENT3* tet)
{ 
  double jac[9],j[9];
  if ((tet->get_nb_fem_noeud()==4))
        {
        FEM_NOEUD* n1=tet->get_fem_noeud(0);
        FEM_NOEUD* n2=tet->get_fem_noeud(1);
        FEM_NOEUD* n3=tet->get_fem_noeud(2);
        FEM_NOEUD* n4=tet->get_fem_noeud(3);
        double *xyz1=n1->get_coord();
        double *xyz2=n2->get_coord();
        double *xyz3=n3->get_coord();
        double *xyz4=n4->get_coord();
        jac[0]=(xyz2[0]-xyz1[0]);
        jac[1]=(xyz2[1]-xyz1[1]);
        jac[2]=(xyz2[2]-xyz1[2]);
        jac[3]=(xyz3[0]-xyz1[0]);
        jac[4]=(xyz3[1]-xyz1[1]);
        jac[5]=(xyz3[2]-xyz1[2]);
        jac[6]=(xyz4[0]-xyz1[0]);
        jac[7]=(xyz4[1]-xyz1[1]);
        jac[8]=(xyz4[2]-xyz1[2]);
        double detj=jac[0]*jac[4]*jac[8]+jac[3]*jac[7]*jac[2]+jac[6]*jac[1]*jac[5]-jac[2]*jac[4]*jac[6]-jac[5]*jac[7]*jac[0]-jac[8]*jac[1]*jac[3];
 
        j[0*3+0] =  (jac[1*3+1]*jac[2*3+2]-jac[1*3+2]*jac[2*3+1])/(detj);
        j[0*3+1] = -(jac[0*3+1]*jac[2*3+2]-jac[0*3+2]*jac[2*3+1])/(detj);
        j[0*3+2] =-(-jac[0*3+1]*jac[1*3+2]+jac[0*3+2]*jac[1*3+1])/(detj);
        j[1*3+0] = -(jac[1*3+0]*jac[2*3+2]-jac[1*3+2]*jac[2*3+0])/(detj);
        j[1*3+1] =  (jac[0*3+0]*jac[2*3+2]-jac[0*3+2]*jac[2*3+0])/(detj);
        j[1*3+2] = -(jac[0*3+0]*jac[1*3+2]-jac[0*3+2]*jac[1*3+0])/(detj);
        j[2*3+0] =-(-jac[1*3+0]*jac[2*3+1]+jac[1*3+1]*jac[2*3+0])/(detj);
        j[2*3+1] = -(jac[0*3+0]*jac[2*3+1]-jac[0*3+1]*jac[2*3+0])/(detj);
        j[2*3+2] =  (jac[0*3+0]*jac[1*3+1]-jac[0*3+1]*jac[1*3+0])/(detj);
 
        //Gradient par element Design
        double d0 = n1->get_solution();
        double d1 = n2->get_solution();
        double d2 = n3->get_solution();
        double d3 = n4->get_solution();
        
        double gd_x = j[0*3+0]*(d1-d0)+j[0*3+1]*(d2-d0)+j[0*3+2]*(d3-d0);
        double gd_y = j[1*3+0]*(d1-d0)+j[1*3+1]*(d2-d0)+j[1*3+2]*(d3-d0);
        double gd_z = j[2*3+0]*(d1-d0)+j[2*3+1]*(d2-d0)+j[2*3+2]*(d3-d0);
                
        OT_VECTEUR_3D gd(gd_x,gd_y,gd_z);
        return(gd.get_longueur());
        }
  return 0.;
}
 
 
 
 
 
void FEM_MAILLAGE_OUTILS::compare_champs_solution(class FEM_SOLUTION* sol1, int numchamp1,FEM_SOLUTION* sol2,int numchamp2)
{
FEM_MAILLAGE* fem1=sol1->get_maillage();  
FEM_MAILLAGE* fem2=sol2->get_maillage();  
MG_GESTIONNAIRE* gest=fem1->get_mg_maillage()->get_gestionnaire();
if (fem1!=fem2) return;
char nom[255];
sprintf(nom,"%s_%s-%s_%s",sol1->get_nom().c_str(),sol1->get_legende(numchamp1).c_str(),sol2->get_nom().c_str(),sol2->get_legende(numchamp1).c_str());
char nomfichier[255];
sprintf(nomfichier,"%s%d_diff_%d%s",sol1->get_nom_fichier().c_str(),numchamp1,numchamp2,sol2->get_nom_fichier().c_str());
FEM_SOLUTION *sol3=new FEM_SOLUTION(fem1,2,nomfichier,1,nom,MAGIC::ENTITE_SOLUTION::ENTITE_NOEUD,MAGIC::TYPE_SOLUTION::SCALAIRE);
sol3->change_legende(0,"Valeur");
sol3->change_legende(1,"|Valeur|");
gest->ajouter_fem_solution(sol3);
int numno=0;
LISTE_FEM_NOEUD::iterator it;
for (FEM_NOEUD* no=fem1->get_premier_noeud(it);no!=NULL;no=fem1->get_suivant_noeud(it))
  {
  double val1=sol1->lire(numno,numchamp1);
  double val2=sol2->lire(numno,numchamp2);
  double erreur=(val1-val2)/val2;
  sol3->ecrire(erreur,numno,0);
  sol3->ecrire(fabs(erreur),numno,1);
  numno++;
  }
}
 
/*void FEM_MAILLAGE_OUTILS::operation_champs_solution(class FEM_SOLUTION* sol1, int numchamp1,FEM_SOLUTION* sol3,int numchamp3,int operation)
{
FEM_MAILLAGE* fem=sol1->get_maillage();
sol1->active_solution(numchamp1);
int ligne=0;
int dim=sol1->get_dim_solution();
 
if (operation==MAGIC::OPERATION_FEM_SOLUTION::EXTRAPOLATION)
  {
  LISTE_FEM_NOEUD::iterator itn;
  for (FEM_NOEUD* no=fem->get_premier_noeud(itn);no!=NULL;no=fem->get_suivant_noeud(itn))
      for (int i=0;i<dim;i++)
        no->change_solution(0.,i);
  }
if ((operation==MAGIC::OPERATION_FEM_SOLUTION::EXTRAPOLATION) || (operation==MAGIC::OPERATION_FEM_SOLUTION::GRAD))
{
  LISTE_FEM_ELEMENT2::iterator it;
  for (FEM_ELEMENT2 *ele=fem->get_premier_element2(it);ele!=NULL;ele=fem->get_suivant_element2(it))
      {
      int nb=ele->get_nb_fem_noeud();
      double x=0,y=0,z=0;
      double uv[2];
      if ((operation==MAGIC::OPERATION_FEM_SOLUTION::GRAD) && (dim==1))
        {
        for (int i=0;i<nb;i++)
              {
              x=x+ele->get_fonction_derive_interpolation(i+1,1,uv)*ele->get_fem_noeud(i)->get_solution();
              y=y+ele->get_fonction_derive_interpolation(i+1,2,uv)*ele->get_fem_noeud(i)->get_solution();
              //z=z+ele->get_fonction_derive_interpolation(i+1,2,uv)*ele->get_fem_noeud(i)->get_solution();
              }
        double jacinv[9];
        double unite=1.;
        if (fem->get_mg_geometrie()!=NULL) unite=fem->get_mg_geometrie()->get_valeur_unite();
        ele->get_inverse_jacob(jacinv,uv,unite);
        double xx=jacinv[0]*x+jacinv[1]*y+jacinv[2]*z;
        double yy=jacinv[3]*x+jacinv[4]*y+jacinv[5]*z;
        double zz=jacinv[6]*x+jacinv[7]*y+jacinv[8]*z;
        sol3->ecrire(xx,ligne,numchamp3,0);
        sol3->ecrire(yy,ligne,numchamp3,1);
        sol3->ecrire(zz,ligne,numchamp3,2);
        ligne++;
        }
      if ((operation==MAGIC::OPERATION_FEM_SOLUTION::GRAD) && (dim==3))
        {
        for (int j=0;j<3;j++)
          {
          for (int i=0;i<nb;i++)
              {
              x=x+ele->get_fonction_derive_interpolation(i+1,1,uv)*ele->get_fem_noeud(i)->get_solution(j);
              y=y+ele->get_fonction_derive_interpolation(i+1,2,uv)*ele->get_fem_noeud(i)->get_solution(j);
              //z=z+ele->get_fonction_derive_interpolation(i+1,2,uv)*ele->get_fem_noeud(i)->get_solution();
              }
          double jacinv[9];
          double unite=1.;
          if (fem->get_mg_geometrie()!=NULL) unite=fem->get_mg_geometrie()->get_valeur_unite();
          ele->get_inverse_jacob(jacinv,uv,unite);
          double xx=jacinv[0]*x+jacinv[1]*y+jacinv[2]*z;
          double yy=jacinv[3]*x+jacinv[4]*y+jacinv[5]*z;
          double zz=jacinv[6]*x+jacinv[7]*y+jacinv[8]*z;
          sol3->ecrire(xx,ligne,numchamp3,0,3*j);
          sol3->ecrire(yy,ligne,numchamp3,1,3*j+1);
          sol3->ecrire(zz,ligne,numchamp3,2,3*j+2);
          }
        ligne++;
        }
      if (operation==MAGIC::OPERATION_FEM_SOLUTION::EXTRAPOLATION)
        {
        for (int i=0;i<nb;i++)
          for (int j=0;j<dim;j++)
            ele->get_fem_noeud(i)->change_solution(ele->get_fem_noeud(i)->get_solution(j)+ele->get_solution(j),j);
        }
      }
}    
if (operation==MAGIC::OPERATION_FEM_SOLUTION::EXTRAPOLATION)
  {
  LISTE_FEM_NOEUD::iterator itn;
  for (FEM_NOEUD* no=fem->get_premier_noeud(itn);no!=NULL;no=fem->get_suivant_noeud(itn))
    {
    for (int i=0;i<dim;i++)
        {
          no->change_solution(no->get_solution(i)/no->get_lien_element2()->get_nb(),i);
          sol3->ecrire(no->get_solution(i),ligne,numchamp3,i);
          
        }
    ligne++;
      
    }
  }
if (operation==MAGIC::OPERATION_FEM_SOLUTION::NORMESUP)
{
sol1->active_solution(numchamp1);
if (sol1->get_entite_solution()==MAGIC::ENTITE_SOLUTION::ENTITE_ELEMENT2)
    {
    LISTE_FEM_ELEMENT2::iterator it;
      for (FEM_ELEMENT2 *ele=fem->get_premier_element2(it);ele!=NULL;ele=fem->get_suivant_element2(it)) 
          {
          double max=fabs(ele->get_solution(0));
          for (int i=1;i<dim;i++)
            if (fabs(ele->get_solution(i))>max) 
              max=fabs(ele->get_solution(i));
          sol3->ecrire(max,ligne,numchamp3);
          ligne++;  
          }
    }
}
 
}
*/
 
void FEM_MAILLAGE_OUTILS::operation_champs_solution(class FEM_SOLUTION* sol1, int numchamp1,FEM_SOLUTION* sol2,int numchamp2,FEM_SOLUTION* sol3,int numchamp3,int operation)
{
FEM_MAILLAGE* fem1=sol1->get_maillage();  
FEM_MAILLAGE* fem2=sol2->get_maillage();  
MG_GESTIONNAIRE* gest=fem1->get_mg_maillage()->get_gestionnaire();
if (fem1!=fem2) return;
int numno=0;
LISTE_FEM_NOEUD::iterator it;
for (FEM_NOEUD* no=fem1->get_premier_noeud(it);no!=NULL;no=fem1->get_suivant_noeud(it))
  {
  double val1=sol1->lire(numno,numchamp1);
  double val2=sol2->lire(numno,numchamp2);
  double val;
  switch (operation)
      {
      case MAGIC::OPERATION_FEM_SOLUTION::ADDITION:
            val=val1+val2;
            break;
      case MAGIC::OPERATION_FEM_SOLUTION::SOUSTRACTION:
            val=val1-val2;
            break;
      case MAGIC::OPERATION_FEM_SOLUTION::MULTIPLICATION:
            val=val1*val2;
            break;
      case MAGIC::OPERATION_FEM_SOLUTION::DIVISION:
            val=val1/val2;
            break;
      case MAGIC::OPERATION_FEM_SOLUTION::COMPARAISON:
            val=fabs(val1-val2);
            break;
      case MAGIC::OPERATION_FEM_SOLUTION::COMPARAISON_RELATIVE:
            val=fabs(val1-val2)/val1;
            break;
      case MAGIC::OPERATION_FEM_SOLUTION::VALEUR_ABSOLUE:
            val=fabs(val1);
            break;
      }
  sol3->ecrire(val,numno,numchamp3);
  numno++;
  }
}
 
 
int FEM_MAILLAGE_OUTILS::projete_solution_maillage(FEM_SOLUTION* sol, int numchamp,FEM_MAILLAGE* fem)
{
if (sol->get_entite_solution()!=MAGIC::ENTITE_SOLUTION::ENTITE_NOEUD) return MAGIC::PROJECTION::TYPE_SOLUTION_INCORRECT;
int dimdest;
if (fem->get_nb_fem_element3()!=0) dimdest=3;
else if (fem->get_nb_fem_element2()!=0) dimdest=2;
else if (fem->get_nb_fem_element1()!=0) dimdest=1;
int dimori=0;
if (sol->get_maillage()->get_nb_fem_element3()!=0) dimori=3;
else if (sol->get_maillage()->get_nb_fem_element2()!=0) dimori=2;
else if (sol->get_maillage()->get_nb_fem_element1()!=0) dimori=1;
if (dimdest!=dimori) return MAGIC::PROJECTION::DIMENSION_MAILLAGE_INCOMPATIBLE;
if (dimdest<2) return MAGIC::PROJECTION::NON_IMPLEMENTER;
if (dimdest==2) return projete_solution_maillage_2D(sol,numchamp,fem);  
if (dimdest==3) return projete_solution_maillage_3D(sol,numchamp,fem);  
return 0;
}
 
 
int FEM_MAILLAGE_OUTILS::projete_solution_maillage_2D(FEM_SOLUTION* sol, int numchamp,FEM_MAILLAGE* fem)
{
sol->active_solution(numchamp);
FEM_MAILLAGE* fem1=sol->get_maillage();
TPL_LISTE_ENTITE<FEM_NOEUD*> lstnoeud;
LISTE_FEM_NOEUD::iterator it;
double xmin=1e308,ymin=1e308,zmin=1e308;
double xmax=-1e308,ymax=-1e308,zmax=-1e308;
for (FEM_NOEUD* no=fem1->get_premier_noeud(it);no!=NULL;no=fem1->get_suivant_noeud(it))
      {
        lstnoeud.ajouter(no);
        if (no->get_x()>xmax) xmax=no->get_x();
        if (no->get_y()>ymax) ymax=no->get_y();
        if (no->get_z()>zmax) zmax=no->get_z();
        if (no->get_x()<xmin) xmin=no->get_x();
        if (no->get_y()<ymin) ymin=no->get_y();
        if (no->get_z()<zmin) zmin=no->get_z();
      }
 
BOITE_3D boite(xmin,ymin,zmin,xmax,ymax,zmax);
boite.change_grosseur(1.1);
TPL_OCTREE<FEM_ELEMENT2*,FEM_NOEUD*> oct;
oct.initialiser2D(3,&lstnoeud,1,boite.get_xmin(),boite.get_ymin(),boite.get_zmin(),boite.get_xmax(),boite.get_ymax(),boite.get_zmax());
 
LISTE_FEM_ELEMENT2::iterator ite;
for (FEM_ELEMENT2* ele=fem1->get_premier_element2(ite);ele!=NULL;ele=fem1->get_suivant_element2(ite))
  oct.inserer(ele);
 
char nom[255];
sprintf(nom,"%sproj",sol->get_nom().c_str());
char nomfichier[255];
sprintf(nomfichier,"%sproj%lu",sol->get_nom_fichier().c_str(),fem->get_id());
 
int dimsol=sol->get_dim_solution();
FEM_SOLUTION *sol2=new FEM_SOLUTION(fem,1,nomfichier,1,nom,MAGIC::ENTITE_SOLUTION::ENTITE_NOEUD,dimsol);
fem->get_mg_maillage()->get_gestionnaire()->ajouter_fem_solution(sol2);
sol2->change_legende(0,sol->get_legende(numchamp).c_str());
int numno=0;
for (FEM_NOEUD* no=fem->get_premier_noeud(it);no!=NULL;no=fem->get_suivant_noeud(it))
    {
    TPL_MAP_ENTITE<FEM_ELEMENT2*> rechercheele;
    oct.rechercher(no->get_x(),no->get_y(),no->get_z(),0.,rechercheele);        
    TPL_MAP_ENTITE<FEM_NOEUD*>::ITERATEUR itno;
    TPL_MAP_ENTITE<FEM_ELEMENT2*>::ITERATEUR itele;
    bool passe=false;
    for (FEM_ELEMENT2* ele2=rechercheele.get_premier(itele);ele2!=NULL;ele2=rechercheele.get_suivant(itele))
          {
          double uv[2];
          bool res=ele2->get_param_element_fini_2D(no->get_coord(),uv);  
          if (res==true)
            if (ele2->valide_parametre_EF(uv)==true)
                {
                double *sol1val=new double[dimsol];
                for (int j=0;j<dimsol;j++)
                    sol1val[j]=0.;
                 for (int i=0;i<ele2->get_nb_fem_noeud();i++)
                  {
                  double interpol=ele2->get_fonction_interpolation(i+1,uv);
                  for (int j=0;j<dimsol;j++)
                    sol1val[j]=sol1val[j]+interpol*ele2->get_fem_noeud(i)->get_solution(j);
                  
                    
                  }
                for (int j=0;j<dimsol;j++)
                  sol2->ecrire(sol1val[j],numno,0,j);
                delete [] sol1val;
                passe=true;
                }
          }
   if (passe==false)
      { 
      std::cerr << "non projete  " <<std::endl;
      }
    numno++;
    }
return MAGIC::PROJECTION::SUCCES;
}
 
 
int FEM_MAILLAGE_OUTILS::projete_solution_maillage_3D(FEM_SOLUTION* sol, int numchamp,FEM_MAILLAGE* fem)
{
sol->active_solution(numchamp);
FEM_MAILLAGE* fem1=sol->get_maillage();
TPL_LISTE_ENTITE<FEM_NOEUD*> lstnoeud;
LISTE_FEM_NOEUD::iterator it;
double xmin=1e308,ymin=1e308,zmin=1e308;
double xmax=-1e308,ymax=-1e308,zmax=-1e308;
for (FEM_NOEUD* no=fem1->get_premier_noeud(it);no!=NULL;no=fem1->get_suivant_noeud(it))
      {
        lstnoeud.ajouter(no);
        if (no->get_x()>xmax) xmax=no->get_x();
        if (no->get_y()>ymax) ymax=no->get_y();
        if (no->get_z()>zmax) zmax=no->get_z();
        if (no->get_x()<xmin) xmin=no->get_x();
        if (no->get_y()<ymin) ymin=no->get_y();
        if (no->get_z()<zmin) zmin=no->get_z();
      }
 
BOITE_3D boite(xmin,ymin,zmin,xmax,ymax,zmax);
boite.change_grosseur(1.1);
TPL_OCTREE<FEM_ELEMENT3*,FEM_NOEUD*> oct;
oct.initialiser(&lstnoeud,1,boite.get_xmin(),boite.get_ymin(),boite.get_zmin(),boite.get_xmax(),boite.get_ymax(),boite.get_zmax());
LISTE_FEM_ELEMENT3::iterator ite;
for (FEM_ELEMENT3* ele=fem1->get_premier_element3(ite);ele!=NULL;ele=fem1->get_suivant_element3(ite))
  oct.inserer(ele);
 
char nom[255];
sprintf(nom,"%sproj",sol->get_nom().c_str());
char nomfichier[255];
sprintf(nomfichier,"%sproj%lu",sol->get_nom_fichier().c_str(),fem->get_id());
 
int dimsol=sol->get_dim_solution();
FEM_SOLUTION *sol2=new FEM_SOLUTION(fem,1,nomfichier,1,nom,MAGIC::ENTITE_SOLUTION::ENTITE_NOEUD,dimsol);
fem->get_mg_maillage()->get_gestionnaire()->ajouter_fem_solution(sol2);
sol2->change_legende(0,sol->get_legende(numchamp).c_str());
int numno=0;
for (FEM_NOEUD* no=fem->get_premier_noeud(it);no!=NULL;no=fem->get_suivant_noeud(it))
    {
    TPL_MAP_ENTITE<FEM_ELEMENT3*> rechercheele;
    oct.rechercher(no->get_x(),no->get_y(),no->get_z(),0.,rechercheele);        
    TPL_MAP_ENTITE<FEM_ELEMENT3*>::ITERATEUR itele;
    bool passe=false;
    int num=rechercheele.get_nb();
    for (FEM_ELEMENT3* ele3=rechercheele.get_premier(itele);ele3!=NULL;ele3=rechercheele.get_suivant(itele))
          {
          double uv[3];
          bool res=ele3->get_param_element_fini(no->get_coord(),uv);  
          if (res==false) std::cout << " non convergence " << std::endl;
          if (res==true)
            if (ele3->valide_parametre_EF(uv)==true)
                {
                double *sol1val=new double[dimsol];
                for (int j=0;j<dimsol;j++)
                    sol1val[j]=0.;
                 for (int i=0;i<ele3->get_nb_fem_noeud();i++)
                  {
                  double interpol=ele3->get_fonction_interpolation(i+1,uv);
                  for (int j=0;j<dimsol;j++)
                    sol1val[j]=sol1val[j]+interpol*ele3->get_fem_noeud(i)->get_solution(j);
                  
                    
                  }
                for (int j=0;j<dimsol;j++)
                  sol2->ecrire(sol1val[j],numno,0,j);
                delete [] sol1val;
                passe=true;
                }
          }
    if (passe==false)
      { 
      std::cerr << "non projete  " <<std::endl;
      }
    numno++;
    }
return MAGIC::PROJECTION::SUCCES;
}