optimisation/src/mgopt_simp.cpp

#include "gestionversion.h"
#include "mgopt_simp.h"
#include "mg_file.h"
#include "mg_export.h"
#include <string.h>
#include <math.h>


class SIMP_TETRA
{
public:
FEM_ELEMENT3* tet;
double densite;
double new_densite;
int design;
double new_denergie;
double energie;
double denergie;
double volume;
int maille_niveau;
int indice;
double distance_ref;
double distance_ref2;
double centre[3];
BOITE_3D get_boite_3D(void) 
        {
        return tet->get_boite_3D();
        };
vector<SIMP_TETRA*> voisin;
vector<SIMP_TETRA*> voisin2;
unsigned long get_id(void) {return tet->get_id();};
double distance(SIMP_TETRA* tet2) 
        {
        double dx=centre[0]-tet2->centre[0];
        double dy=centre[1]-tet2->centre[1];
        double dz=centre[2]-tet2->centre[2];
        return sqrt(dx*dx+dy*dy+dz*dz);
        };
};


MGOPT_SIMP::MGOPT_SIMP():MGOPT()
{
params.ajouter("f",0.3,"Fraction volumique de la methode SIMP");
params.ajouter("rminc",0.,"0. Valeur de la carte de taille (si pas de maillage le parametre fichiercarte permet d'utiliser une carte de taille) sinon valeur de rmin pour le lissage de la compliance");
params.ajouter("rmind",0.,"0. Valeur de la carte de taille (si pas de maillage le parametre fichiercarte permet d'utiliser une carte de taille) sinon valeur de rmin pour le lissage de la densite");
params.ajouter("coefvoisinc",1.0,"Coefficient de multiplication de rminc");
params.ajouter("coefvoisind",1.0,"Coefficient de multiplication de rmind");
params.ajouter("nbrniveau",150.,"Nombre de matériaux utilisés");
params.ajouter("p",3.,"Coefficient de penalite du module d'Young");
params.ajouter("ro_void",0.001,"Densite minimale consideree comme nulle");
params.ajouter("m",0.2,"Limite d'evolution de densité entre deux itérations");
params.ajouter("eta",0.5,"Coefficent d'affectation du beta");
params.ajouter("type_lissage",1.,"0. Pas de lissage 1. Lissage de la compliance 2. Lissage de la densité 3. Lissage de la compliance et de la densite");
params.ajouter("lissage_densite",1.,"0. Complet 1. Derniere iteration");
params.ajouter("type_lissage_densite",1.,"0. Distance 1. Distance pondéré 2. Gaussien 3. Gaussien pondéré");
params.ajouter("type_lissage_compliance",1.,"0. Sigmund 1997 1. Sigmund 2007");
params.ajouter("kc",1.0,"Ponderation de la distance dans le lissage de la compliance");
params.ajouter("kd",1.0,"Ponderation de la distance dans le lissage de la densité");
params.ajouter("convergence_lagrange",0.1,"Critere de convergence de la recherche du multiplicateur de Lagrange en \%");
params.ajouter("iter_max",50.,"Nombre d'iteration maximale dans le processus SIMP");
params.ajouter("critere_convergence",0.5,"Critere de convergence de la méthode SIMP en \%");
params.ajouter("ro_min",0.8,"Seuil de conservation des éléments");
params.ajouter("iter_vue",0.,"generation d'un fihcier gmsh tous les iter_vue. 0  pas de generation");
}

MGOPT_SIMP::MGOPT_SIMP(MGOPT_SIMP &mdd):MGOPT(mdd)
{
}


MGOPT_SIMP::~MGOPT_SIMP()
{
}


void MGOPT_SIMP::optimisation(FEM_MAILLAGE* fem)
{
MG_VOLUME* vol=fem->get_mg_geometrie()->get_mg_volume(0);
affiche((char*)"     Initialisation");
vector<SIMP_TETRA*> lsttet;
double f=params.get_valeur("f");
double rminc=params.get_valeur("rminc");
double rmind=params.get_valeur("rmind");
if (((rminc<1e-16)||(rmind<1e-16)) && (carte==NULL))
        {
        affiche((char*)"     Lecture de la carte de taille");
        carte=new FCT_GENERATEUR_3D<4>;
        std::string fichiercarte=params.get_nom("fichiercarte");
        carte->lire((char *)fichiercarte.c_str());
        }
double volume_tot;
double volume_design=0;
double volume_non_design=0;
double unite=fem->get_mg_maillage()->get_mg_geometrie()->get_valeur_unite();
LISTE_FEM_ELEMENT3::iterator it;
int ntet=0;
for (FEM_ELEMENT3* tet=fem->get_premier_element3(it);tet!=NULL;tet=fem->get_suivant_element3(it))
        {
        SIMP_TETRA* tet2=new SIMP_TETRA;
        tet2->tet=tet;
        tet2->tet->change_numero(ntet);
        ntet++;
        tet2->centre[0]=0.;
        tet2->centre[1]=0.;
        tet2->centre[2]=0.;
        int nbnoeud=tet->get_nb_fem_noeud();
        for (int i=0;i<nbnoeud;i++)
                {
                tet2->centre[0]=tet2->centre[0]+tet->get_fem_noeud(i)->get_x()*unite;
                tet2->centre[1]=tet2->centre[1]+tet->get_fem_noeud(i)->get_y()*unite;
                tet2->centre[2]=tet2->centre[2]+tet->get_fem_noeud(i)->get_z()*unite;
                }
        tet2->centre[0]=tet2->centre[0]/nbnoeud;                
        tet2->centre[1]=tet2->centre[1]/nbnoeud;                
        tet2->centre[2]=tet2->centre[2]/nbnoeud;                
        int lig,col;
        double jac[9],uv[3]={1./4.,1./4.,1./4.};
        tet2->volume=1./6.*tet->get_jacobien(jac,uv,lig,col,unite);
        tet2->distance_ref=rminc;
        tet2->distance_ref2=rmind;
        double xyz[3]={tet2->centre[0],tet2->centre[1],tet2->centre[2]};
        double val[9];
        if (tet2->distance_ref<1e-16) 
                {
                carte->evaluer(xyz,val);
                tet2->distance_ref=1./sqrt(val[0]);
                }
        if (tet2->distance_ref2<1e-16) 
                {
                carte->evaluer(xyz,val);
                tet2->distance_ref2=1./sqrt(val[0]);
                }
        tet2->distance_ref=params.get_valeur("coefvoisinc")*tet2->distance_ref*unite;
        tet2->distance_ref2=params.get_valeur("coefvoisind")*tet2->distance_ref2*unite;
        if (((MG_TETRA*)tet->get_mg_element_maillage())->get_origine()!=IMPOSE) 
                {
                tet2->design=1;
                tet2->densite=f;
                volume_design=volume_design+tet2->volume;
                }
         else 
                {
                tet2->design=0;
                tet2->densite=1.;
                volume_non_design=volume_non_design+tet2->volume;
                }
        lsttet.insert(lsttet.end(),tet2);
        }
volume_tot=volume_design+volume_non_design;
double volumemoyen=volume_tot/fem->get_nb_fem_element3();
affiche((char*)"     Recherche de voisins");
int nbsimptet=lsttet.size();
for (int i=0;i<nbsimptet;i++)
        {
        SIMP_TETRA* tet=lsttet[i];
        if (tet->design==0) continue;
        ajouter_voisin(i,tet,lsttet);
        }
int nbiteration=1;
int ok=0;
int niveaumax=(int)params.get_valeur("nbrniveau");
int itervue=(int)params.get_valeur("iter_vue");
double p=params.get_valeur("p");
vector<double> Ctotiter;
while (ok==0)
        {
        if (nbiteration>params.get_valeur("iter_max")) break;
        char message[255];
        sprintf(message,"     Iteration %d",nbiteration);
        affiche(message);
        vector<FEM_ELEMENT3*>  *lstniveau=new vector<FEM_ELEMENT3*>[niveaumax+1];
        double densmin = params.get_valeur("ro_void");
        for (int i=0;i<nbsimptet;i++)
                {
                SIMP_TETRA* tet=lsttet[i];
                tet->maille_niveau=(int)((tet->densite-densmin)*niveaumax/(1.-densmin))+1;
                if (tet->maille_niveau>niveaumax) tet->maille_niveau=niveaumax;
                lstniveau[tet->maille_niveau].insert(lstniveau[tet->maille_niveau].end(),tet->tet);
                }
        MG_EXPORT exp;
        int version_aster=params.get_valeur("version_aster");
        exp.aster(vol,fem,nometude,2,(char*)"00000001",p,niveaumax,lstniveau);
        delete [] lstniveau;
        affiche((char*)"        Calcul aster");
        char nomfichiertmp[255];
        sprintf(nomfichiertmp,"%s/as_run %s.export  1>aster.log 2>&1",getenv("PATHASTER"),nometude);
        int code=system(nomfichiertmp);
        if (code!=0)
                {
                sprintf(nomfichiertmp,"          Code de sortie aster : %d",code);
                affiche(nomfichiertmp);
                }
        affiche((char*)"        Analyse resultat");
        recupere_energie(lsttet);
        double Ctot=0;
        for (int i=0;i<nbsimptet;i++)
                {
                SIMP_TETRA* tet=lsttet[i];
                tet->denergie= -p*2.*tet->energie/tet->densite;
                Ctot=Ctot+2.*tet->energie;
                }
        Ctotiter.insert(Ctotiter.end(),Ctot);
        sprintf(message,"             Compliance %le",Ctot);
        affiche(message);
        if (nbiteration!=1)
                {
                int nb=Ctotiter.size();
                double c1=Ctotiter[nb-2];
                double c2=Ctotiter[nb-1];
                double ecart=fabs((c2-c1)/c1)*100.;
                if (ecart<params.get_valeur("critere_convergence")) ok=1;
                sprintf(message,"             Ecart %lf%%",ecart);
                affiche(message);
                }
        // lissage compliance
        double kc=params.get_valeur("kc");
        int type_lissage =(int)params.get_valeur("type_lissage");
        int lissage_densite =(int)params.get_valeur("lissage_densite");
        int type_lissage_densite =(int)params.get_valeur("type_lissage_densite");
        int type_lissage_compliance=(int)params.get_valeur("type_lissage_compliance");
        for (int i=0;i<nbsimptet;i++)
                {
                SIMP_TETRA* tet=lsttet[i];
                if ((type_lissage==0) || (type_lissage==2))
                        tet->new_denergie=tet->denergie;
                if ((type_lissage==1) || (type_lissage==3))
                        {
                        if (tet->design == 1)
                                {
                                double hi= pow(tet->distance_ref,kc);
                                double hisensi = hi*tet->densite*tet->denergie;
                                if (type_lissage_compliance==1)
                                        {
                                        hi=hi/tet->volume;
                                        hisensi=hisensi/max(volumemoyen*0.00001,tet->volume);
                                        }
                                int nbvoisin=tet->voisin.size();
                                for (int j = 0 ; j<nbvoisin ; j++)
                                        {
                                        SIMP_TETRA* tet2=tet->voisin[j];
                                        if (tet2->design == 1)
                                                {
                                                double dist=tet->distance(tet2);
                                                double hj =fabs(pow(tet->distance_ref - dist,kc));
                                                if (type_lissage_compliance==0)
                                                        {
                                                        hisensi=hisensi+tet2->densite*hj*tet2->denergie;
                                                        hi=hi+hj;
                                                        }
                                                if (type_lissage_compliance==1)
                                                        {
                                                        hisensi=hisensi+tet2->densite*hj*tet2->denergie/max(volumemoyen,tet2->volume);
                                                        hi=hi+hj/tet2->volume;
                                                        }
                                                hisensi=hisensi+tet2->densite*hj*tet2->denergie;
                                                hi=hi+hj;
                                                }
                                        }
                                tet->new_denergie = hisensi/hi/max(tet->densite,densmin);
                                }
                        else tet->new_denergie=tet->denergie;
                        }
                }
        // application formule pas encore connue
        double l1= 1e-8 ;
        double l2= 1e8;
        double m=params.get_valeur("m");
        double eta=params.get_valeur("eta");
        double convergence_lagrange=params.get_valeur("convergence_lagrange");
        double critere_densite = 0.0;
        int oklagrange=0;
        int compteurlagrange=0;
        while (oklagrange==0)
                {
                compteurlagrange++;
                double lmid = 0.5*(l2+l1);
                critere_densite = 0.0;
                        for (int i=0;i<nbsimptet;i++)
                        {
                        SIMP_TETRA* tet=lsttet[i];
                        if (tet->design==1)
                                {
                                double x1 = tet->densite+m;
                                double beta=-tet->new_denergie/lmid/tet->volume;
                                double x2 = tet->densite*pow(beta,eta);
                                double x3 = min(x1,x2);
                                double x4 = 1.;
                                double x5 = min(x3,x4);
                                double x6 = tet->densite-m;
                                double x7 = max(x5,x6);
                                double x8 = densmin;
                                tet->new_densite = max(x7,x8);
                                critere_densite = critere_densite + tet->new_densite*tet->volume;
                                }
                        else
                                {
                                tet->new_densite=1.;
                                }
                        }
                if (critere_densite - f*volume_design  > 0.)
                        l1=lmid;
                else
                        l2=lmid;
                if (100.*fabs(critere_densite- f*volume_design )/(f*volume_design)<convergence_lagrange) oklagrange=1;
                if (compteurlagrange>500) oklagrange=2;
                }
        if (oklagrange==1) sprintf(message,"             Convergence multiplicateur lagrange %le%%",100.*(critere_densite- f*volume_design )/(f*volume_design));
        if (oklagrange==2) sprintf(message,"             Divergence multiplicateur lagrange %le %le",l1,l2);
        affiche(message);
        // lissage de densite.
        double kd=params.get_valeur("kd");
        for (int i=0;i<nbsimptet;i++)
                {
                SIMP_TETRA* tet=lsttet[i];
                if ((type_lissage==2) || (type_lissage==3))
                        {
                        if (((lissage_densite==1)&&(ok==1)) || (lissage_densite==0))
                                {
                                if (tet->design == 1)
                                        {
                                        double widensite=0.;
                                        double wi= 0.;
                                        wi=poid_lissage(0.,tet->distance_ref2,kd,tet->volume,type_lissage_densite);
                                        widensite = wi*tet->new_densite;
                                        int nbvoisin=tet->voisin2.size();
                                        for (int j = 0 ; j<nbvoisin ; j++)
                                                {
                                                SIMP_TETRA* tet2=tet->voisin2[j];
                                                if (tet2->design == 1)
                                                        {
                                                        double dist=tet->distance(tet2);
                                                        double wj=poid_lissage(dist,tet->distance_ref2,kd,tet2->volume,type_lissage_densite);
                                                        wi = wi+wj;
                                                        widensite = widensite + tet2->new_densite*wj;
                                                        }
                                                }
                                        tet->densite = widensite/wi;
                                        }
                                else tet->densite=tet->new_densite;
                                }
                        else tet->densite=tet->new_densite;
                        }
                else tet->densite=tet->new_densite;
                }
        if (itervue!=0)
          if (nbiteration%itervue==0)
              {
              affiche((char*)"        Sauvegarde résultat iteration");
              sprintf(message,"%s_densite_iter%d.soltmp",nometudesortie,nbiteration);
              //int nbsolution=gestd->get_nb_fem_solution();
              //for (int i=nbsolution;i>0;i--)
                    //gestd->supprimer_fem_solution_du_gestionnaire(i-1);
              char message2[50];
              sprintf(message2,"Iteration %d",nbiteration);
              FEM_SOLUTION* solution=new FEM_SOLUTION(fem,1,message,fem->get_nb_fem_element3(),message2,ENTITE_ELEMENT3);
              gestd->ajouter_fem_solution(solution);
              sprintf(message,"Densite");
              solution->change_legende(0,message); 
              for (int i=0;i<nbsimptet;i++)
                  {
                  SIMP_TETRA* tet=lsttet[i];
                  solution->ecrire(i,0,tet->densite);
                  }
              }
        nbiteration++;
        
        
        }
if (itervue!=0)
  {     
  affiche((char*)"     Sauvegarde GMSH résultats iterations");
  MG_EXPORT exp;
  char nomfichier[500];
  sprintf(nomfichier,"%s_iterations",nometudesortie);
  exp.gmsh(fem,nomfichier);
  }
double seuil=params.get_valeur("ro_min");
affiche((char*)"     Ecriture des donnees finales");
char message[255];
sprintf(message,"%s.compliance",nometudesortie);
FILE *out=fopen(message,"wt");
time_t heurefin=time(NULL);
struct tm  tfin = *localtime(&heurefin);
fprintf(out,"*****************************************************\n");
fprintf(out,"           Optimisateur de topologie\n");
fprintf(out,"                   ERICCA - UQTR\n");
fprintf(out,"*****************************************************\n");
fprintf(out,"    Etude        : %s\n",nometude);
fprintf(out,"    Date         : %02u-%02u-%04u\n", tdebut.tm_mday, tdebut.tm_mon+1, 1900 + tdebut.tm_year);
fprintf(out,"    Heure debut  : %02u:%02u:%02u\n", tdebut.tm_hour, tdebut.tm_min, tdebut.tm_sec);
fprintf(out,"    Heure fin    : %02u:%02u:%02u\n", tfin.tm_hour, tfin.tm_min, tfin.tm_sec);
fprintf(out,"    Parametres   : \n");
int nbparam=params.get_nb();
for (int i=0;i<nbparam;i++)
  {
  fprintf(out,"                   %s=%lf\n",params.get_nom(i).c_str(),params.get_valeur(i));
  }
fprintf(out,"*****************************************************\n");
fprintf(out,"         Compliance après chaque itération\n");
fprintf(out,"*****************************************************\n");
for (int i=0;i<Ctotiter.size();i++)
        fprintf(out,"%le\n",Ctotiter[i]);
double crit1=0.;
double crit2=0.;
double crit3=0.;
double critv1[10]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
double critv2[10]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
for (int i=0;i<nbsimptet;i++)
        {
        SIMP_TETRA* tet=lsttet[i];
        if (tet->design==1)
                {
                crit1=crit1+tet->densite*tet->volume;
                if (tet->densite>seuil)
                        {
                        crit2=crit2+tet->densite*tet->volume;
                        crit3=crit3+tet->volume;
                        }
                for (int j=1;j<10;j++)
                        if (tet->densite>0.1*j)
                                {
                                critv1[j]=critv1[j]+tet->densite*tet->volume;
                                critv2[j]=critv2[j]+tet->volume;
                                }
                }
        }
fprintf(out,"*****************************************************\n");
fprintf(out,"         Résultat\n");
fprintf(out,"*****************************************************\n");
fprintf(out,"Volume du design : %le \n",volume_design);
fprintf(out,"Objectif du volume de design : %le \n",volume_design*f);
fprintf(out,"Volume de design obtenu : %le \n",crit1);
fprintf(out,"Volume de design obtenu avec le seuil: %le \n",crit2);
fprintf(out,"Volume reel de design obtenu avec le seuil: %le \n\n",crit3);
fprintf(out," : Volume : Volume reel : Objectif\n");
fprintf(out," : %le : %le : %le\n",volume_design,volume_design,volume_design);
for (int j=1;j<10;j++)
        fprintf(out,"Volume de design obtenu avec le seuil de %.1f : %le : %le : %le \n",j*0.1,critv1[j],critv2[j],volume_design*f);
fprintf(out,"*************************************************************\n");
fclose(out);
LISTE_FEM_NOEUD::iterator itnoeud;
int nbfemnoeud=fem->get_nb_fem_noeud();
/*double *nume=new double[nbfemnoeud];
double *deno=new double[nbfemnoeud];
int cpt=0;
for (FEM_NOEUD *nd=fem->get_premier_noeud(itnoeud);nd!=NULL;nd=fem->get_suivant_noeud(itnoeud))
        {
        nd->change_numero(cpt);
        nume[cpt]=0.;
        deno[cpt]=0.;
        cpt++;
        }*/
sprintf(message,"%s_densite1.sol",nometudesortie);
int nbsolution=gestd->get_nb_fem_solution();
for (int i=nbsolution;i>0;i--)
        {
          FEM_SOLUTION* sol=gestd->get_fem_solution(i-1);
          std::string nom=sol->get_nom_fichier();
          char message[500],extension[500];
          sprintf(message,"%s",nom.c_str());
          char *p=strrchr(message,'.');
          strncpy(extension,p,strlen(message)+message-p);
          if (strcmp(extension,".soltmp")==0)
               gestd->supprimer_fem_solution(i-1);
            else gestd->supprimer_fem_solution_du_gestionnaire(i-1);
        }


FEM_SOLUTION* solution=new FEM_SOLUTION(fem,1,message,fem->get_nb_fem_element3(),"Optimisation",ENTITE_ELEMENT3);
gestd->ajouter_fem_solution(solution);
solution->change_legende(0,"Densite"); 
for (int i=0;i<nbsimptet;i++)
        {
        SIMP_TETRA* tet=lsttet[i];
        if (tet->design==1) 
                if (tet->densite>seuil) 
                        ((MG_TETRA*)tet->tet->get_mg_element_maillage())->change_origine(OPTIMISE); 
                else
                        ((MG_TETRA*)tet->tet->get_mg_element_maillage())->change_origine(MAILLEUR_AUTO); 
        solution->ecrire(i,0,tet->densite);
        /*for (int j=0;j<tet->tet->get_nb_fem_noeud();j++)
                {
                FEM_NOEUD* noeud=tet->tet->get_fem_noeud(j);
                nume[noeud->get_numero()]=nume[noeud->get_numero()]+tet->volume*tet->densite;
                deno[noeud->get_numero()]=deno[noeud->get_numero()]+tet->volume;
                }       */
        }
/*sprintf(message,"%s_densite2.sol",nometudesortie);
FEM_SOLUTION* solution2=new FEM_SOLUTION(fem,1,message,fem->get_nb_fem_noeud(),"Optimisation",ENTITE_NOEUD);
gestd->ajouter_fem_solution(solution2);
solution2->change_legende(0,"Densite");
cpt=0;
for (FEM_NOEUD *nd=fem->get_premier_noeud(itnoeud);nd!=NULL;nd=fem->get_suivant_noeud(itnoeud))
        {
        solution2->ecrire(cpt,0,nume[cpt]/deno[cpt]);
        cpt++;
        }
delete [] deno;
delete [] nume;*/
int nb=lsttet.size();
for (int i=0;i<nb;i++) delete lsttet[i];
}

double MGOPT_SIMP::poid_lissage(double dist,double distref,double k,double volume,int type)
{
double wi;
if (type==0) wi=pow(distref-dist,k);
if (type==1) wi=pow(distref-dist,k)*volume;
if (type==2) wi=exp(-dist*dist/2./distref/distref/9.)/2./M_PI/(distref/3.);
if (type==3) wi=volume*exp(-dist*dist/2./distref/distref/9.)/2./M_PI/(distref/3.);
return fabs(wi);
        
}

void MGOPT_SIMP::adapte_resultat(char *nomgestd,char *nomparam)
{
if (nomparam!=NULL) lire_params(nomparam);
affiche((char*)"");
affiche((char*)"*************************");
affiche((char*)"Optimisation de topologie");
affiche((char*)"*************************");
affiche((char*)"");
affiche((char*)"");
affiche((char*)"Changement du seuil dans les resultats");
double seuil=params.get_valeur("seuil");
gestd=new MG_FILE(nomgestd);
FEM_MAILLAGE* fem=gestd->get_fem_maillage(0);
FEM_SOLUTION* solution=gestd->get_fem_solution(0);
solution->active_solution(0);
LISTE_FEM_ELEMENT3::iterator it;
for (FEM_ELEMENT3 *tet=fem->get_premier_element3(it);tet!=NULL;tet=fem->get_suivant_element3(it))
        {
        if (((MG_TETRA*)tet->get_mg_element_maillage())->get_origine()!=IMPOSE) 
                if (tet->get_solution()>seuil) 
                        ((MG_TETRA*)tet->get_mg_element_maillage())->change_origine(OPTIMISE); 
                else
                        ((MG_TETRA*)tet->get_mg_element_maillage())->change_origine(MAILLEUR_AUTO); 
        
        }
affiche((char*)"Enregistrement");
gestd->enregistrer(nomgestd);
affiche((char*)"Enregistrement sous GMSH");
char *p=strchr(nomgestd,'.');
strncpy(nometude,nomgestd,p-nomgestd);
nometude[p-nomgestd]=0;
MG_EXPORT exp;
char nomfichier[500];
sprintf(nomfichier,"%s_mg",nometude);
exp.gmsh(fem->get_mg_maillage(),nomfichier);
sprintf(nomfichier,"%s_fem",nometude);
exp.gmsh(fem,nomfichier);
affiche((char*)"Fin");
}


void MGOPT_SIMP::recupere_energie(vector<class SIMP_TETRA*> lsttet)
{
char message[750];
sprintf(message,"%s.resu",nometude);
FILE* in=fopen(message,"rt");
int fin=0;
do
        {
        char* res=fgets(message,750,in);
        if (feof(in)) fin=1;
        char mot1[100]; 
        char mot2[100]; 
        char mot3[100]; 
        char mot4[100]; 
        char mot5[100]; 
        char mot6[100]; 
        char mot7[100]; 
        char mot8[100]; 
        char mot9[100]; 
        char mot10[100]; 
        int numlu=sscanf(message,"%s %s %s %s %s %s %s %s %s %s",mot1,mot2,mot3,mot4,mot5,mot6,mot7,mot8,mot9,mot10);
        if (numlu>9)
        if (strcmp(mot1,"CHAMP")==0)
         if (strcmp(mot2,"PAR")==0) 
          if (strcmp(mot3,"ELEMENT")==0) 
           if (strcmp(mot4,"CONSTANT")==0) 
            if (strcmp(mot5,"SUR")==0) 
             if (strcmp(mot6,"ELEMENT")==0) 
              if (strcmp(mot7,"DE")==0) 
               if (strcmp(mot8,"NOM")==0) 
                if (strcmp(mot9,"SYMBOLIQUE")==0) 
                 if (strcmp(mot10,"ENERGIE")==0) 
                        {
                        int fin2=0;
                        int passe=0;
                        int nbelement=0;
                        do 
                                {
                                char message[750];
                                char* res=fgets(message,750,in);
                                char mot1[500];
                                char mot2[500];
                                int numlu=sscanf(message,"%s %s",mot1,mot2);
                                int decalage;
                                if ((numlu==2) && (strcmp(mot2,"TOTALE")==0))
                                        {
                                        char *p=strchr(mot1,'M')+1;
                                        int num=atoi(p);
                                        if (passe==0) {passe=1;decalage=num;}
                                        char* res=fgets(message,750,in);
                                        double val;
                                        sscanf(message,"%lf",&val);
                                        lsttet[num-decalage]->energie=val;
                                        nbelement++;
                                        }
                                if (nbelement == lsttet.size()) {fin2=1;fin=0;}
                                }
                                
                        while (fin2==0);
                        }
        }
while (fin==0);
fclose(in);
}


void MGOPT_SIMP::ajouter_voisin(int i,SIMP_TETRA* tet,vector<SIMP_TETRA*> &lst)
{
tet->indice=i;
int nbnoeud=tet->tet->get_nb_fem_noeud();
int correspondance[4];
if (nbnoeud==4) 
        {
        correspondance[0]=0;
        correspondance[1]=1;
        correspondance[2]=2;
        correspondance[3]=3;
        }
if (nbnoeud==10) 
        {
        correspondance[0]=0;
        correspondance[1]=2;
        correspondance[2]=4;
        correspondance[3]=9;
        }
SIMP_TETRA* tetcour=tet;
int ok=0;
int compteur=0;
vector<SIMP_TETRA*>& lstvoisin=tet->voisin;
if (tet->distance_ref<tet->distance_ref2)
        lstvoisin=tet->voisin2;
while (ok==0)
        {
        for (int j=0;j<4;j++)
                {
                int num=correspondance[j];
                FEM_NOEUD *noeud=tetcour->tet->get_fem_noeud(num);
                int nbtetra=noeud->get_lien_element3()->get_nb();
                for (int k=0;k<nbtetra;k++)
                        {
                        FEM_ELEMENT3* ftet=noeud->get_lien_element3()->get(k);
                        SIMP_TETRA* stet=lst[ftet->get_numero()];
                        if (stet->indice!=i)
                                {
                                stet->indice=i;
                                double dist=tet->distance(stet);
                                if (dist<tet->distance_ref) 
                                                tet->voisin.insert(tet->voisin.end(),stet);
                                if (dist<tet->distance_ref2) 
                                                tet->voisin2.insert(tet->voisin2.end(),stet);
                                }
                        }
                }
        if (compteur>=lstvoisin.size()) ok=1;
        else 
                {
                tetcour=lstvoisin[compteur];
                compteur++;
                }
                
        }

}
Revision:	342
Committed:	Thu May 31 20:19:01 2012 UTC (12 years, 11 months ago) by francois
File size:	22163 byte(s)
Log Message:	correction bug precedent mise a jour