ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/optimisation/src/mgopt_simp.cpp
Revision: 267
Committed: Fri Oct 29 20:14:04 2010 UTC (14 years, 6 months ago) by francois
File size: 19039 byte(s)
Log Message:
Ajout d'un filtre dans le lissage de la derivee de la compliance en optimisation

File Contents

# User Rev Content
1 francois 239 #include "gestionversion.h"
2     #include "mgopt_simp.h"
3     #include "mg_file.h"
4     #include "mg_export.h"
5     #include <string.h>
6 francois 258 #include <math.h>
7 francois 239
8    
9    
10    
11     class SIMP_TETRA
12     {
13     public:
14     FEM_TETRA* tet;
15     double densite;
16 francois 242 double new_densite;
17 francois 239 int design;
18 francois 242 double new_denergie;
19 francois 239 double energie;
20     double denergie;
21     double volume;
22     int maille_niveau;
23     int indice;
24     double distance_ref;
25 francois 247 double distance_ref2;
26 francois 239 double centre[3];
27     BOITE_3D get_boite_3D(void)
28     {
29     return tet->get_boite_3D();
30     };
31     vector<SIMP_TETRA*> voisin;
32 francois 247 vector<SIMP_TETRA*> voisin2;
33 francois 239 unsigned long get_id(void) {return tet->get_id();};
34     double distance(SIMP_TETRA* tet2)
35     {
36     double dx=centre[0]-tet2->centre[0];
37     double dy=centre[1]-tet2->centre[1];
38     double dz=centre[2]-tet2->centre[2];
39     return sqrt(dx*dx+dy*dy+dz*dz);
40     };
41     };
42    
43    
44    
45    
46    
47    
48    
49    
50    
51    
52     MGOPT_SIMP::MGOPT_SIMP():MGOPT()
53     {
54 francois 247 params.ajouter("f",0.3,"Fraction volumique de la methode SIMP");
55     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");
56     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");
57     params.ajouter("coefvoisinc",1.0,"Coefficient de multiplication de rminc");
58     params.ajouter("coefvoisind",1.0,"Coefficient de multiplication de rmind");
59 francois 243 params.ajouter("nbrniveau",150.,"Nombre de matériaux utilisés");
60 francois 247 params.ajouter("p",3.,"Coefficient de penalite du module d'Young");
61     params.ajouter("ro_void",0.001,"Densite minimale consideree comme nulle");
62 francois 243 params.ajouter("m",0.2,"Limite d'evolution de densité entre deux itérations");
63     params.ajouter("eta",0.5,"Coefficent d'affectation du beta");
64 francois 247 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");
65     params.ajouter("lissage_densite",1.,"0. Complet 1. Derniere iteration");
66     params.ajouter("type_lissage_densite",1.,"0. Distance 1. Distance pondéré 2. Gaussien 3. Gaussien pondéré");
67 francois 267 params.ajouter("type_lissage_compliance",1.,"0. Sigmund 1997 1. Sigmund 2007");
68 francois 247 params.ajouter("kc",1.0,"Ponderation de la distance dans le lissage de la compliance");
69     params.ajouter("kd",1.0,"Ponderation de la distance dans le lissage de la densité");
70 francois 243 params.ajouter("convergence_lagrange",0.1,"Critere de convergence de la recherche du multiplicateur de Lagrange en \%");
71     params.ajouter("iter_max",50.,"Nombre d'iteration maximale dans le processus SIMP");
72     params.ajouter("critere_convergence",0.5,"Critere de convergence de la méthode SIMP en \%");
73 francois 247 params.ajouter("ro_min",0.8,"Seuil de conservation des éléments");
74 francois 239 }
75    
76     MGOPT_SIMP::MGOPT_SIMP(MGOPT_SIMP &mdd):MGOPT(mdd)
77     {
78     }
79    
80    
81     MGOPT_SIMP::~MGOPT_SIMP()
82     {
83     }
84    
85    
86     void MGOPT_SIMP::optimisation(FEM_MAILLAGE* fem)
87     {
88     affiche(" Initialisation");
89     vector<SIMP_TETRA*> lsttet;
90 francois 247 double f=params.get_valeur("f");
91     double rminc=params.get_valeur("rminc");
92     double rmind=params.get_valeur("rmind");
93     if (((rminc<1e-16)||(rmind<1e-16)) && (carte==NULL))
94 francois 246 {
95     affiche(" Lecture de la carte de taille");
96     carte=new FCT_GENERATEUR_3D<4>;
97     std::string fichiercarte=params.get_nom("fichiercarte");
98     carte->lire((char *)fichiercarte.c_str());
99     }
100 francois 239 double volume_tot;
101     double volume_design=0;
102     double volume_non_design=0;
103     double unite=fem->get_mg_maillage()->get_mg_geometrie()->get_valeur_unite();
104     LISTE_FEM_TETRA::iterator it;
105     int ntet=0;
106     for (FEM_TETRA* tet=fem->get_premier_tetra(it);tet!=NULL;tet=fem->get_suivant_tetra(it))
107     {
108     SIMP_TETRA* tet2=new SIMP_TETRA;
109     tet2->tet=tet;
110     tet2->tet->change_numero(ntet);
111     ntet++;
112     tet2->centre[0]=0.;
113     tet2->centre[1]=0.;
114     tet2->centre[2]=0.;
115     int nbnoeud=tet->get_nb_fem_noeud();
116     for (int i=0;i<nbnoeud;i++)
117     {
118 francois 242 tet2->centre[0]=tet2->centre[0]+tet->get_fem_noeud(i)->get_x()*unite;
119     tet2->centre[1]=tet2->centre[1]+tet->get_fem_noeud(i)->get_y()*unite;
120     tet2->centre[2]=tet2->centre[2]+tet->get_fem_noeud(i)->get_z()*unite;
121 francois 239 }
122     tet2->centre[0]=tet2->centre[0]/nbnoeud;
123     tet2->centre[1]=tet2->centre[1]/nbnoeud;
124     tet2->centre[2]=tet2->centre[2]/nbnoeud;
125     int lig,col;
126     double jac[9],uv[3]={1./4.,1./4.,1./4.};
127 francois 242 tet2->volume=1./6.*tet->get_jacobien(jac,uv,lig,col,unite);
128 francois 247 tet2->distance_ref=rminc;
129     tet2->distance_ref2=rmind;
130 francois 246 double xyz[3]={tet2->centre[0],tet2->centre[1],tet2->centre[2]};
131     double val[9];
132     if (tet2->distance_ref<1e-16)
133 francois 242 {
134 francois 246 carte->evaluer(xyz,val);
135     tet2->distance_ref=1./sqrt(val[0]);
136 francois 242 }
137 francois 247 if (tet2->distance_ref2<1e-16)
138     {
139     carte->evaluer(xyz,val);
140     tet2->distance_ref2=1./sqrt(val[0]);
141     }
142     tet2->distance_ref=params.get_valeur("coefvoisinc")*tet2->distance_ref*unite;
143     tet2->distance_ref2=params.get_valeur("coefvoisind")*tet2->distance_ref2*unite;
144 francois 239 if (((MG_TETRA*)tet->get_mg_element_maillage())->get_origine()!=IMPOSE)
145     {
146     tet2->design=1;
147 francois 247 tet2->densite=f;
148 francois 239 volume_design=volume_design+tet2->volume;
149     }
150     else
151     {
152     tet2->design=0;
153     tet2->densite=1.;
154     volume_non_design=volume_non_design+tet2->volume;
155     }
156     lsttet.insert(lsttet.end(),tet2);
157     }
158     affiche(" Recherche de voisins");
159     int nbsimptet=lsttet.size();
160     for (int i=0;i<nbsimptet;i++)
161     {
162     SIMP_TETRA* tet=lsttet[i];
163     if (tet->design==0) continue;
164     ajouter_voisin(i,tet,lsttet);
165     }
166     int nbiteration=1;
167     int ok=0;
168     int niveaumax=(int)params.get_valeur("nbrniveau");
169 francois 247 double p=params.get_valeur("p");
170 francois 239 vector<double> Ctotiter;
171     while (ok==0)
172     {
173 francois 247 if (nbiteration>params.get_valeur("iter_max")) break;
174 francois 239 char message[255];
175     sprintf(message," Iteration %d",nbiteration);
176     affiche(message);
177     vector<FEM_TETRA*> *lstniveau=new vector<FEM_TETRA*>[niveaumax+1];
178 francois 247 double densmin = params.get_valeur("ro_void");
179 francois 239 for (int i=0;i<nbsimptet;i++)
180     {
181     SIMP_TETRA* tet=lsttet[i];
182     tet->maille_niveau=(int)((tet->densite-densmin)*niveaumax/(1.-densmin))+1;
183     if (tet->maille_niveau>niveaumax) tet->maille_niveau=niveaumax;
184     lstniveau[tet->maille_niveau].insert(lstniveau[tet->maille_niveau].end(),tet->tet);
185     }
186     MG_EXPORT exp;
187 francois 247 exp.aster(fem,nometude,2,"00000001",p,niveaumax,lstniveau);
188 francois 239 delete [] lstniveau;
189     affiche(" Calcul aster");
190     char nomfichiertmp[255];
191     sprintf(nomfichiertmp,"%s/as_run %s.export 1>aster.log 2>&1",getenv("PATHASTER"),nometude);
192 francois 242 int code=system(nomfichiertmp);
193     if (code!=0)
194     {
195     sprintf(nomfichiertmp," Code de sortie aster : %d",code);
196     affiche(nomfichiertmp);
197     }
198 francois 239 affiche(" Analyse resultat");
199     recupere_energie(lsttet);
200     double Ctot=0;
201     for (int i=0;i<nbsimptet;i++)
202     {
203     SIMP_TETRA* tet=lsttet[i];
204 francois 247 tet->denergie= -p*2.*tet->energie/tet->densite;
205 francois 239 Ctot=Ctot+2.*tet->energie;
206     }
207     Ctotiter.insert(Ctotiter.end(),Ctot);
208     sprintf(message," Compliance %le",Ctot);
209     affiche(message);
210     if (nbiteration!=1)
211     {
212     int nb=Ctotiter.size();
213     double c1=Ctotiter[nb-2];
214     double c2=Ctotiter[nb-1];
215     double ecart=fabs((c2-c1)/c1)*100.;
216     if (ecart<params.get_valeur("critere_convergence")) ok=1;
217     sprintf(message," Ecart %lf%%",ecart);
218     affiche(message);
219     }
220 francois 247 // lissage compliance
221     double kc=params.get_valeur("kc");
222 francois 243 int type_lissage =(int)params.get_valeur("type_lissage");
223 francois 247 int lissage_densite =(int)params.get_valeur("lissage_densite");
224     int type_lissage_densite =(int)params.get_valeur("type_lissage_densite");
225 francois 267 int type_lissage_compliance=(int)params.get_valeur("type_lissage_compliance");
226 francois 242 for (int i=0;i<nbsimptet;i++)
227     {
228     SIMP_TETRA* tet=lsttet[i];
229 francois 247 if ((type_lissage==0) || (type_lissage==2))
230 francois 243 tet->new_denergie=tet->denergie;
231 francois 247 if ((type_lissage==1) || (type_lissage==3))
232 francois 242 {
233 francois 243 if (tet->design == 1)
234 francois 242 {
235 francois 247 double hi= pow(tet->distance_ref,kc);
236     double hisensi = hi*tet->densite*tet->denergie;
237 francois 267 if (type_lissage_compliance==1)
238     {
239     hi=hi/tet->volume;
240     hisensi=hisensi/tet->volume;
241     }
242     int nbvoisin=tet->voisin.size();
243 francois 243 for (int j = 0 ; j<nbvoisin ; j++)
244     {
245     SIMP_TETRA* tet2=tet->voisin[j];
246     if (tet2->design == 1)
247     {
248 francois 247 double dist=tet->distance(tet2);
249     double hj =fabs(pow(tet->distance_ref - dist,kc));
250 francois 267 if (type_lissage_compliance==0)
251     {
252     hisensi=hisensi+tet2->densite*hj*tet2->denergie;
253     hi=hi+hj;
254     }
255     if (type_lissage_compliance==1)
256     {
257     hisensi=hisensi+tet2->densite*hj*tet2->denergie/tet2->volume;
258     hi=hi+hj/tet2->volume;
259     }
260     hisensi=hisensi+tet2->densite*hj*tet2->denergie;
261 francois 247 hi=hi+hj;
262 francois 243 }
263 francois 242 }
264 francois 247 tet->new_denergie = hisensi/hi/tet->densite;
265 francois 243 }
266     else tet->new_denergie=tet->denergie;
267 francois 242 }
268     }
269 francois 239 // application formule pas encore connue
270     double l1= 1e-8 ;
271     double l2= 1e8;
272     double m=params.get_valeur("m");
273 francois 241 double eta=params.get_valeur("eta");
274 francois 242 double convergence_lagrange=params.get_valeur("convergence_lagrange");
275     double critere_densite = 0.0;
276     int oklagrange=0;
277     int compteurlagrange=0;
278     while (oklagrange==0)
279 francois 239 {
280 francois 242 compteurlagrange++;
281     double lmid = 0.5*(l2+l1);
282     critere_densite = 0.0;
283     for (int i=0;i<nbsimptet;i++)
284 francois 239 {
285     SIMP_TETRA* tet=lsttet[i];
286     if (tet->design==1)
287     {
288     double x1 = tet->densite+m;
289 francois 242 double beta=-tet->new_denergie/lmid/tet->volume;
290 francois 241 double x2 = tet->densite*pow(beta,eta);
291     double x3 = min(x1,x2);
292     double x4 = 1.;
293     double x5 = min(x3,x4);
294     double x6 = tet->densite-m;
295     double x7 = max(x5,x6);
296     double x8 = densmin;
297 francois 239 tet->new_densite = max(x7,x8);
298 francois 242 critere_densite = critere_densite + tet->new_densite*tet->volume;
299     }
300 francois 239 else
301     {
302     tet->new_densite=1.;
303     }
304 francois 242 }
305 francois 247 if (critere_densite - f*volume_design > 0.)
306 francois 239 l1=lmid;
307     else
308     l2=lmid;
309 francois 247 if (100.*fabs(critere_densite- f*volume_design )/(f*volume_design)<convergence_lagrange) oklagrange=1;
310 francois 243 if (compteurlagrange>500) oklagrange=2;
311 francois 242 }
312 francois 247 if (oklagrange==1) sprintf(message," Convergence multiplicateur lagrange %le%%",100.*(critere_densite- f*volume_design )/(f*volume_design));
313 francois 242 if (oklagrange==2) sprintf(message," Divergence multiplicateur lagrange %le %le",l1,l2);
314     affiche(message);
315 francois 247 // lissage de densite.
316     double kd=params.get_valeur("kd");
317 francois 239 for (int i=0;i<nbsimptet;i++)
318 francois 242 {
319 francois 239 SIMP_TETRA* tet=lsttet[i];
320 francois 247 if ((type_lissage==2) || (type_lissage==3))
321     {
322     if (((lissage_densite==1)&&(ok==1)) || (lissage_densite==0))
323     {
324     if (tet->design == 1)
325     {
326     double widensite=0.;
327     double wi= 0.;
328     wi=poid_lissage(0.,tet->distance_ref2,kd,tet->volume,type_lissage_densite);
329     widensite = wi*tet->new_densite;
330     int nbvoisin=tet->voisin2.size();
331     for (int j = 0 ; j<nbvoisin ; j++)
332     {
333     SIMP_TETRA* tet2=tet->voisin2[j];
334     if (tet2->design == 1)
335     {
336     double dist=tet->distance(tet2);
337     double wj=poid_lissage(dist,tet->distance_ref2,kd,tet2->volume,type_lissage_densite);
338     wi = wi+wj;
339     widensite = widensite + tet2->new_densite*wj;
340     }
341     }
342     tet->densite = widensite/wi;
343     }
344     else tet->densite=tet->new_densite;
345     }
346     else tet->densite=tet->new_densite;
347     }
348     else tet->densite=tet->new_densite;
349 francois 242 }
350 francois 239 nbiteration++;
351     }
352 francois 247 double seuil=params.get_valeur("ro_min");
353 francois 245 affiche(" Ecriture des donnees finales");
354 francois 239 char message[255];
355 francois 244 sprintf(message,"%s.compliance",nometudesortie);
356 francois 239 FILE *out=fopen(message,"wt");
357     for (int i=0;i<Ctotiter.size();i++)
358     fprintf(out,"%le\n",Ctotiter[i]);
359 francois 245 double crit1=0.;
360     double crit2=0.;
361     double crit3=0.;
362     double critv1[10]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
363     double critv2[10]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
364     for (int i=0;i<nbsimptet;i++)
365     {
366     SIMP_TETRA* tet=lsttet[i];
367     if (tet->design==1)
368     {
369     crit1=crit1+tet->densite*tet->volume;
370     if (tet->densite>seuil)
371     {
372     crit2=crit2+tet->densite*tet->volume;
373     crit3=crit3+tet->volume;
374     }
375     for (int j=1;j<10;j++)
376     if (tet->densite>0.1*j)
377     {
378     critv1[j]=critv1[j]+tet->densite*tet->volume;
379     critv2[j]=critv2[j]+tet->volume;
380     }
381     }
382     }
383     fprintf(out,"\n\n\n**************************************************************\n");
384     fprintf(out,"Volume du design : %le \n",volume_design);
385 francois 247 fprintf(out,"Objectif du volume de design : %le \n",volume_design*f);
386 francois 245 fprintf(out,"Volume de design obtenu : %le \n",crit1);
387     fprintf(out,"Volume de design obtenu avec le seuil: %le \n",crit2);
388     fprintf(out,"Volume reel de design obtenu avec le seuil: %le \n\n",crit3);
389     fprintf(out," : Volume : Volume reel : Objectif\n");
390     fprintf(out," : %le : %le : %le\n",volume_design,volume_design,volume_design);
391     for (int j=1;j<10;j++)
392 francois 247 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);
393 francois 245 fprintf(out,"*************************************************************\n");
394     fclose(out);
395 francois 239 LISTE_FEM_NOEUD::iterator itnoeud;
396     int nbfemnoeud=fem->get_nb_fem_noeud();
397 francois 242 /*double *nume=new double[nbfemnoeud];
398 francois 239 double *deno=new double[nbfemnoeud];
399     int cpt=0;
400     for (FEM_NOEUD *nd=fem->get_premier_noeud(itnoeud);nd!=NULL;nd=fem->get_suivant_noeud(itnoeud))
401     {
402     nd->change_numero(cpt);
403     nume[cpt]=0.;
404     deno[cpt]=0.;
405     cpt++;
406 francois 242 }*/
407 francois 243 sprintf(message,"%s_densite1.sol",nometudesortie);
408 francois 241 int nbsolution=gestd->get_nb_fem_solution();
409 francois 242 for (int i=nbsolution;i>0;i--)
410 francois 243 gestd->supprimer_fem_solution_du_gestionnaire(i-1);
411 francois 239 FEM_SOLUTION* solution=new FEM_SOLUTION(fem,1,message,fem->get_nb_fem_tetra(),"Optimisation",ENTITE_TETRA);
412     gestd->ajouter_fem_solution(solution);
413     solution->change_legende(0,"Densite");
414     for (int i=0;i<nbsimptet;i++)
415     {
416     SIMP_TETRA* tet=lsttet[i];
417     if (tet->design==1)
418     if (tet->densite>seuil)
419     ((MG_TETRA*)tet->tet->get_mg_element_maillage())->change_origine(OPTIMISE);
420     else
421     ((MG_TETRA*)tet->tet->get_mg_element_maillage())->change_origine(MAILLEUR_AUTO);
422     solution->ecrire(i,0,tet->densite);
423 francois 242 /*for (int j=0;j<tet->tet->get_nb_fem_noeud();j++)
424 francois 239 {
425     FEM_NOEUD* noeud=tet->tet->get_fem_noeud(j);
426     nume[noeud->get_numero()]=nume[noeud->get_numero()]+tet->volume*tet->densite;
427     deno[noeud->get_numero()]=deno[noeud->get_numero()]+tet->volume;
428 francois 242 } */
429 francois 239 }
430 francois 243 /*sprintf(message,"%s_densite2.sol",nometudesortie);
431 francois 239 FEM_SOLUTION* solution2=new FEM_SOLUTION(fem,1,message,fem->get_nb_fem_noeud(),"Optimisation",ENTITE_NOEUD);
432     gestd->ajouter_fem_solution(solution2);
433     solution2->change_legende(0,"Densite");
434     cpt=0;
435     for (FEM_NOEUD *nd=fem->get_premier_noeud(itnoeud);nd!=NULL;nd=fem->get_suivant_noeud(itnoeud))
436     {
437     solution2->ecrire(cpt,0,nume[cpt]/deno[cpt]);
438     cpt++;
439     }
440     delete [] deno;
441 francois 242 delete [] nume;*/
442 francois 239 int nb=lsttet.size();
443     for (int i=0;i<nb;i++) delete lsttet[i];
444     }
445    
446 francois 247 double MGOPT_SIMP::poid_lissage(double dist,double distref,double k,double volume,int type)
447     {
448     double wi;
449     if (type==0) wi=pow(distref-dist,k);
450     if (type==1) wi=pow(distref-dist,k)*volume;
451     if (type==2) wi=exp(-dist*dist/2./distref/distref/9.)/2./M_PI/(distref/3.);
452     if (type==3) wi=volume*exp(-dist*dist/2./distref/distref/9.)/2./M_PI/(distref/3.);
453     return fabs(wi);
454    
455     }
456 francois 239
457 francois 240 void MGOPT_SIMP::adapte_resultat(char *nomgestd,char *nomparam)
458     {
459     if (nomparam!=NULL) lire_params(nomparam);
460     affiche("");
461     affiche("*************************");
462     affiche("Optimisation de topologie");
463     affiche("*************************");
464     affiche("");
465     affiche("");
466     affiche("Changement du seuil dans les resultats");
467     double seuil=params.get_valeur("seuil");
468     gestd=new MG_FILE(nomgestd);
469     FEM_MAILLAGE* fem=gestd->get_fem_maillage(0);
470     FEM_SOLUTION* solution=gestd->get_fem_solution(0);
471     solution->active_solution(0);
472     LISTE_FEM_TETRA::iterator it;
473     for (FEM_TETRA *tet=fem->get_premier_tetra(it);tet!=NULL;tet=fem->get_suivant_tetra(it))
474     {
475     if (((MG_TETRA*)tet->get_mg_element_maillage())->get_origine()!=IMPOSE)
476     if (tet->get_solution()>seuil)
477     ((MG_TETRA*)tet->get_mg_element_maillage())->change_origine(OPTIMISE);
478     else
479     ((MG_TETRA*)tet->get_mg_element_maillage())->change_origine(MAILLEUR_AUTO);
480    
481     }
482     affiche("Enregistrement");
483     gestd->enregistrer(nomgestd);
484 francois 245 affiche("Enregistrement sous GMSH");
485     char *p=strchr(nomgestd,'.');
486     strncpy(nometude,nomgestd,p-nomgestd);
487     nometude[p-nomgestd]=0;
488     MG_EXPORT exp;
489     char nomfichier[500];
490     sprintf(nomfichier,"%s_mg",nometude);
491     exp.gmsh(fem->get_mg_maillage(),nomfichier);
492     sprintf(nomfichier,"%s_fem",nometude);
493     exp.gmsh(fem,nomfichier);
494 francois 240 affiche("Fin");
495     }
496 francois 239
497    
498     void MGOPT_SIMP::recupere_energie(vector<class SIMP_TETRA*> lsttet)
499     {
500     char message[750];
501     sprintf(message,"%s.resu",nometude);
502     FILE* in=fopen(message,"rt");
503     int fin=0;
504     do
505     {
506     fgets(message,750,in);
507     if (feof(in)) fin=1;
508     char mot1[100];
509     char mot2[100];
510     char mot3[100];
511     char mot4[100];
512     char mot5[100];
513     char mot6[100];
514     char mot7[100];
515     char mot8[100];
516     char mot9[100];
517     char mot10[100];
518     int numlu=sscanf(message,"%s %s %s %s %s %s %s %s %s %s",mot1,mot2,mot3,mot4,mot5,mot6,mot7,mot8,mot9,mot10);
519     if (numlu>9)
520     if (strcmp(mot1,"CHAMP")==0)
521     if (strcmp(mot2,"PAR")==0)
522     if (strcmp(mot3,"ELEMENT")==0)
523     if (strcmp(mot4,"CONSTANT")==0)
524     if (strcmp(mot5,"SUR")==0)
525     if (strcmp(mot6,"L'ELEMENT")==0)
526     if (strcmp(mot7,"DE")==0)
527     if (strcmp(mot8,"NOM")==0)
528     if (strcmp(mot9,"SYMBOLIQUE")==0)
529     if (strcmp(mot10,"EPOT_ELEM_DEPL")==0)
530     {
531     int fin2=0;
532     int passe=0;
533     int nbelement=0;
534     do
535     {
536     char message[750];
537     fgets(message,750,in);
538     char mot1[500];
539     char mot2[500];
540     int numlu=sscanf(message,"%s %s",mot1,mot2);
541     int decalage;
542     if ((numlu==2) && (strcmp(mot2,"TOTALE")==0))
543     {
544     char *p=strchr(mot1,'M')+1;
545     int num=atoi(p);
546     if (passe==0) {passe=1;decalage=num;}
547     fgets(message,750,in);
548     double val;
549     sscanf(message,"%lf",&val);
550     lsttet[num-decalage]->energie=val;
551     nbelement++;
552     }
553     if (nbelement == lsttet.size()) {fin2=1;fin=0;}
554     }
555    
556     while (fin2==0);
557     }
558     }
559     while (fin==0);
560     fclose(in);
561     }
562    
563    
564    
565     void MGOPT_SIMP::ajouter_voisin(int i,SIMP_TETRA* tet,vector<SIMP_TETRA*> &lst)
566     {
567     tet->indice=i;
568     int nbnoeud=tet->tet->get_nb_fem_noeud();
569     int correspondance[4];
570     if (nbnoeud==4)
571     {
572     correspondance[0]=0;
573     correspondance[1]=1;
574     correspondance[2]=2;
575     correspondance[3]=3;
576     }
577     if (nbnoeud==10)
578     {
579     correspondance[0]=0;
580     correspondance[1]=2;
581     correspondance[2]=4;
582     correspondance[3]=9;
583     }
584     SIMP_TETRA* tetcour=tet;
585     int ok=0;
586     int compteur=0;
587 francois 247 vector<SIMP_TETRA*>& lstvoisin=tet->voisin;
588     if (tet->distance_ref<tet->distance_ref2)
589     lstvoisin=tet->voisin2;
590 francois 239 while (ok==0)
591     {
592     for (int j=0;j<4;j++)
593     {
594     int num=correspondance[j];
595     FEM_NOEUD *noeud=tetcour->tet->get_fem_noeud(num);
596     int nbtetra=noeud->get_lien_tetra()->get_nb();
597     for (int k=0;k<nbtetra;k++)
598     {
599     FEM_TETRA* ftet=noeud->get_lien_tetra()->get(k);
600     SIMP_TETRA* stet=lst[ftet->get_numero()];
601     if (stet->indice!=i)
602     {
603     stet->indice=i;
604 francois 247 double dist=tet->distance(stet);
605     if (dist<tet->distance_ref)
606 francois 239 tet->voisin.insert(tet->voisin.end(),stet);
607 francois 247 if (dist<tet->distance_ref2)
608     tet->voisin2.insert(tet->voisin2.end(),stet);
609 francois 239 }
610     }
611     }
612 francois 247 if (compteur>=lstvoisin.size()) ok=1;
613 francois 239 else
614     {
615 francois 247 tetcour=lstvoisin[compteur];
616 francois 239 compteur++;
617     }
618 francois 247
619 francois 239 }
620    
621     }