ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/mailleur/src/mailleur_fem.cpp
Revision: 1180
Committed: Fri Sep 20 20:10:22 2024 UTC (7 months, 3 weeks ago) by francois
File size: 50343 byte(s)
Log Message:
Correction dans la quadratisation sans geometrie

File Contents

# User Rev Content
1 francois 1158 //####//------------------------------------------------------------
2     //####//------------------------------------------------------------
3     //####// MAGiC
4     //####// Jean Christophe Cuilliere et Vincent FRANCOIS
5     //####// Departement de Genie Mecanique - UQTR
6     //####//------------------------------------------------------------
7     //####// MAGIC est un projet de recherche de l equipe ERICCA
8     //####// du departement de genie mecanique de l Universite du Quebec a Trois Rivieres
9     //####// http://www.uqtr.ca/ericca
10     //####// http://www.uqtr.ca/
11     //####//------------------------------------------------------------
12     //####//------------------------------------------------------------
13     //####//
14     //####// mailleur_fem.cpp
15     //####//
16     //####//------------------------------------------------------------
17     //####//------------------------------------------------------------
18     //####// COPYRIGHT 2000-2024
19 francois 1164 //####// Derniere modification par francois
20 francois 1168 //####// jeu 25 jui 2024 18:16:33 EDT
21 francois 1158 //####//------------------------------------------------------------
22     //####//------------------------------------------------------------
23 francois 883
24    
25     #include "gestionversion.h"
26     #include "mailleur_fem.h"
27     #include "mg_gestionnaire.h"
28     #include "fem_segment2.h"
29     #include "fem_mini_segment2.h"
30     #include "fem_segment3.h"
31     #include "fem_triangle3.h"
32     #include "fem_triangle6.h"
33     #include "fem_quadrangle4.h"
34     #include "fem_quadrangle8.h"
35     #include "fem_tetra4.h"
36     #include "fem_tetra10.h"
37     #include "fem_hexa8.h"
38     #include "fem_hexa20.h"
39     #include "fem_penta6.h"
40     #include "fem_penta15.h"
41     #include "fem_graph_noeud.h"
42     #include "fem_maillage_quadratique_outils.h"
43     #include "ot_decalage_parametre.h"
44 francois 1164 #include "mg_arete_element.h"
45     #include "mg_face_element.h"
46     #include <math.h>
47 francois 883
48    
49    
50     MAILLEUR_FEM::MAILLEUR_FEM(OT_CPU* comp):MAILLEUR(false,comp)
51     {
52    
53     }
54    
55     MAILLEUR_FEM::MAILLEUR_FEM(MAILLEUR_FEM &mdd):MAILLEUR(mdd)
56     {
57     ini_param();
58     }
59    
60    
61     MAILLEUR_FEM::~MAILLEUR_FEM()
62     {
63     }
64    
65    
66 francois 1164 int MAILLEUR_FEM::maille(FEM_MAILLAGE* fem,bool courbure_discrete,int num)
67 francois 883 {
68     int degre=fem->get_degre();
69     MG_MAILLAGE* maillage=fem->get_mg_maillage();
70     MG_GEOMETRIE* mg_geometrie=maillage->get_mg_geometrie();
71     TPL_MAP_ENTITE<MG_SEGMENT*> mini;
72     if (mg_geometrie!=NULL)
73     if (mg_geometrie->get_nb_mg_poutre()!=0)
74 francois 1164 if (mg_geometrie->get_nb_mg_volume()!=0)
75     recherche_connexion_multidimension(mini,maillage,mg_geometrie);
76 francois 883 if (degre==1) construire_lineaire(fem,mini,maillage,mg_geometrie);
77 francois 1164 if (degre==2) construire_quadratique(fem,mini,maillage,mg_geometrie,courbure_discrete);
78 francois 883 if (num>0) optimise_numerotation(fem);
79     if (num>1) copie_numerotation_opt(fem);
80     BOITE_3D b;
81     int nx,ny,nz;
82     bool structure=maillage->get_param_structure(b,nx,ny,nz);
83     if (structure==true)
84     fem->change_param_structure(b,nx,ny,nz);
85 couturad 966 return OK;
86 francois 883 }
87    
88    
89     void MAILLEUR_FEM::recherche_connexion_multidimension(TPL_MAP_ENTITE<MG_SEGMENT*> &mini,MG_MAILLAGE* maillage,MG_GEOMETRIE* mg_geometrie)
90     {
91     LISTE_MG_NOEUD::iterator it;
92     for (MG_NOEUD* no=maillage->get_premier_noeud(it);no!=NULL;no=maillage->get_suivant_noeud(it))
93     {
94     if ((no->get_lien_hexa()->get_nb()>0) || (no->get_lien_tetra()->get_nb()))
95     if (no->get_lien_segment()->get_nb()>0)
96     {
97     bool estjoint=false;
98     for (int i=0;i<no->get_lien_segment()->get_nb();i++)
99     {
100     LISTE_MG_POUTRE::iterator itp;
101     for (MG_POUTRE* pou=mg_geometrie->get_premier_poutre(itp);pou!=NULL;pou=mg_geometrie->get_suivant_poutre(itp))
102     if (pou->contient_element(no->get_lien_segment()->get(i))) {estjoint=true;}
103     }
104     if (estjoint)
105     {
106     for (int i=0;i<no->get_lien_segment()->get_nb();i++)
107     {
108     bool bord=false;
109     LISTE_MG_POUTRE::iterator itp;
110     int nbele2=no->get_lien_segment()->get(i)->get_lien_triangle()->get_nb();
111     for (int j=0;j<nbele2;j++)
112     if (no->get_lien_segment()->get(i)->get_lien_triangle()->get(j)->get_lien_tetra()->get_nb()>0)
113     if (no->get_lien_segment()->get(i)->get_lien_triangle()->get(j)->get_lien_tetra()->get_nb()<2) bord=true;
114     nbele2=no->get_lien_segment()->get(i)->get_lien_quadrangle()->get_nb();
115     for (int j=0;j<nbele2;j++)
116     if (no->get_lien_segment()->get(i)->get_lien_quadrangle()->get(j)->get_lien_hexa()->get_nb()>0)
117     if (no->get_lien_segment()->get(i)->get_lien_quadrangle()->get(j)->get_lien_hexa()->get_nb()<2) bord=true;
118     if (bord==true) mini.ajouter(no->get_lien_segment()->get(i));
119     }
120     //std::cout << no->get_id() << " est un joint 1D-3D : " << mini.get_nb() << " mini " << std::endl;
121    
122     }
123     }
124    
125     }
126     LISTE_MG_TRIANGLE::iterator itt;
127     for (MG_TRIANGLE* tri=maillage->get_premier_triangle(itt);tri!=NULL;tri=maillage->get_suivant_triangle(itt))
128     if (tri->get_origine()==MAGIC::ORIGINE::SECTION)
129     {
130     mini.ajouter(tri->get_segment1());
131     mini.ajouter(tri->get_segment2());
132     mini.ajouter(tri->get_segment3());
133     }
134    
135     LISTE_MG_QUADRANGLE::iterator itq;
136     for (MG_QUADRANGLE* qua=maillage->get_premier_quadrangle(itq);qua!=NULL;qua=maillage->get_suivant_quadrangle(itq))
137     if (qua->get_origine()==MAGIC::ORIGINE::SECTION)
138     {
139     mini.ajouter(qua->get_segment1());
140     mini.ajouter(qua->get_segment2());
141     mini.ajouter(qua->get_segment3());
142     mini.ajouter(qua->get_segment4());
143     }
144    
145    
146     }
147    
148    
149    
150    
151 couturad 968 int MAILLEUR_FEM::construire_lineaire(FEM_MAILLAGE* fem,TPL_MAP_ENTITE<MG_SEGMENT*> &mini,MG_MAILLAGE* maillage,MG_GEOMETRIE* mg_geometrie)
152 francois 883 {
153     LISTE_MG_NOEUD::iterator it_noeud;
154     MG_NOEUD* mgnoeud=maillage->get_premier_noeud(it_noeud);
155     std::vector<FEM_NOEUD*> lstnoeuddirect;
156     unsigned int i=0;
157     while (mgnoeud)
158     {
159     mgnoeud->change_nouveau_numero(i);
160     FEM_NOEUD* femnoeud=new FEM_NOEUD(mgnoeud);
161     femnoeud->change_numero(i+1);
162     i++;
163     fem->ajouter_fem_noeud(femnoeud);
164     lstnoeuddirect.insert(lstnoeuddirect.end(),femnoeud);
165     if (mgnoeud->get_lien_topologie()!=NULL)
166     if (mgnoeud->get_lien_topologie()->get_dimension()==0)
167     {
168     FEM_ELEMENT0* femele=new FEM_ELEMENT0(mgnoeud->get_lien_topologie(),mgnoeud,&femnoeud);
169     fem-> ajouter_fem_element0(femele);
170     }
171     mgnoeud=maillage->get_suivant_noeud(it_noeud);
172    
173     }
174     int dimsansgeo=0;
175     if (mg_geometrie==NULL)
176     {
177     if (maillage->get_nb_mg_tetra()+maillage->get_nb_mg_hexa()!=0) dimsansgeo=3;
178     else if (maillage->get_nb_mg_triangle()+maillage->get_nb_mg_quadrangle()!=0) dimsansgeo=2;
179     else if (maillage->get_nb_mg_segment()!=0) dimsansgeo=1;
180     }
181     if (mg_geometrie!=NULL)
182     if (strcmp(mg_geometrie->get_type_geometrie(),"VIRTUEL")==0)
183     {
184     if (maillage->get_nb_mg_tetra()+maillage->get_nb_mg_hexa()!=0) dimsansgeo=3;
185     else if (maillage->get_nb_mg_triangle()+maillage->get_nb_mg_quadrangle()!=0) dimsansgeo=2;
186     else if (maillage->get_nb_mg_segment()!=0) dimsansgeo=1;
187     }
188    
189     if ((mg_geometrie!=NULL) || ( (mg_geometrie==NULL) && (dimsansgeo==1)))
190     {
191     LISTE_MG_SEGMENT::iterator it_seg;
192     MG_SEGMENT * mgseg = maillage->get_premier_segment(it_seg);
193     while (mgseg)
194     {
195     if (mgseg->get_lien_topologie()!=NULL)
196     if (mgseg->get_lien_topologie()->get_dimension()!=1)
197     {
198     mgseg=maillage->get_suivant_segment(it_seg);
199     continue;
200     }
201     if (mg_geometrie!=NULL)
202     if (mgseg->get_lien_topologie()==NULL)
203     if (dimsansgeo!=1)
204     {
205     mgseg=maillage->get_suivant_segment(it_seg);
206     continue;
207     }
208     FEM_NOEUD *tabnoeud[2];
209     tabnoeud[0]=lstnoeuddirect[mgseg->get_noeud1()->get_nouveau_numero()];
210     tabnoeud[1]=lstnoeuddirect[mgseg->get_noeud2()->get_nouveau_numero()];
211     FEM_SEGMENT2* seg=new FEM_SEGMENT2(mgseg,tabnoeud);
212     fem->ajouter_fem_element1(seg);
213     mgseg=maillage->get_suivant_segment(it_seg);
214     }
215     TPL_MAP_ENTITE<MG_SEGMENT*>::ITERATEUR itm;
216     for (MG_SEGMENT* mgseg=mini.get_premier(itm);mgseg!=NULL;mgseg=mini.get_suivant(itm))
217     {
218     FEM_NOEUD *tabnoeud[2];
219     tabnoeud[0]=lstnoeuddirect[mgseg->get_noeud1()->get_nouveau_numero()];
220     tabnoeud[1]=lstnoeuddirect[mgseg->get_noeud2()->get_nouveau_numero()];
221     FEM_SEGMENT2* seg=new FEM_MINI_SEGMENT2(mgseg,tabnoeud);
222     fem->ajouter_fem_element1(seg);
223     }
224     }
225     if ((mg_geometrie!=NULL) || ( (mg_geometrie==NULL) && (dimsansgeo==2)))
226     {
227     LISTE_MG_TRIANGLE::iterator it_tri;
228     MG_TRIANGLE* mgtri=maillage->get_premier_triangle(it_tri);
229     while (mgtri)
230     {
231     if (mgtri->get_lien_topologie()!=NULL) if (mgtri->get_lien_topologie()->get_dimension()!=2)
232     {
233     mgtri=maillage->get_suivant_triangle(it_tri);
234     continue;
235     }
236     if (mg_geometrie!=NULL)
237     if (mgtri->get_lien_topologie()==NULL)
238     if (dimsansgeo!=2)
239     {
240     mgtri=maillage->get_suivant_triangle(it_tri);
241     continue;
242     }
243     FEM_NOEUD *tabnoeud[3];
244     tabnoeud[0]=lstnoeuddirect[mgtri->get_noeud1()->get_nouveau_numero()];
245     tabnoeud[1]=lstnoeuddirect[mgtri->get_noeud2()->get_nouveau_numero()];
246     tabnoeud[2]=lstnoeuddirect[mgtri->get_noeud3()->get_nouveau_numero()];
247     FEM_TRIANGLE3* tri=new FEM_TRIANGLE3(mgtri,tabnoeud);
248     fem->ajouter_fem_element2(tri);
249     mgtri=maillage->get_suivant_triangle(it_tri);
250     }
251     LISTE_MG_QUADRANGLE::iterator it_quad;
252     MG_QUADRANGLE* mgquad=maillage->get_premier_quadrangle(it_quad);
253     while (mgquad)
254     {
255     if (mgquad->get_lien_topologie()!=NULL) if (mgquad->get_lien_topologie()->get_dimension()!=2) {
256     mgquad=maillage->get_suivant_quadrangle(it_quad);
257     continue;
258     }
259     if (mg_geometrie!=NULL)
260     if (mgquad->get_lien_topologie()==NULL)
261     if (dimsansgeo!=2)
262     {
263     mgquad=maillage->get_suivant_quadrangle(it_quad);
264     continue;
265     }
266     FEM_NOEUD *tabnoeud[4];
267     tabnoeud[0]=lstnoeuddirect[mgquad->get_noeud1()->get_nouveau_numero()];
268     tabnoeud[1]=lstnoeuddirect[mgquad->get_noeud2()->get_nouveau_numero()];
269     tabnoeud[2]=lstnoeuddirect[mgquad->get_noeud3()->get_nouveau_numero()];
270     tabnoeud[3]=lstnoeuddirect[mgquad->get_noeud4()->get_nouveau_numero()];
271     FEM_QUADRANGLE4* quad=new FEM_QUADRANGLE4(mgquad,tabnoeud);
272     fem->ajouter_fem_element2(quad);
273     mgquad=maillage->get_suivant_quadrangle(it_quad);
274     }
275     }
276     if ((mg_geometrie!=NULL) || ( (mg_geometrie==NULL) && (dimsansgeo==3)))
277     {
278     LISTE_MG_TETRA::iterator it_tetra;
279     MG_TETRA* mgtetra=maillage->get_premier_tetra(it_tetra);
280     while (mgtetra)
281     {
282     FEM_NOEUD *tabnoeud[4];
283     tabnoeud[0]=lstnoeuddirect[mgtetra->get_noeud1()->get_nouveau_numero()];
284     tabnoeud[1]=lstnoeuddirect[mgtetra->get_noeud2()->get_nouveau_numero()];
285     tabnoeud[2]=lstnoeuddirect[mgtetra->get_noeud3()->get_nouveau_numero()];
286     tabnoeud[3]=lstnoeuddirect[mgtetra->get_noeud4()->get_nouveau_numero()];
287     FEM_TETRA4* tet=new FEM_TETRA4(mgtetra,tabnoeud);
288     fem->ajouter_fem_element3(tet);
289     mgtetra=maillage->get_suivant_tetra(it_tetra);
290     }
291     LISTE_MG_HEXA::iterator it_hex;
292     MG_HEXA* mghex=maillage->get_premier_hexa(it_hex);
293     while (mghex)
294     {
295     FEM_NOEUD *tabnoeud[8];
296     tabnoeud[0]=lstnoeuddirect[mghex->get_noeud1()->get_nouveau_numero()];
297     tabnoeud[1]=lstnoeuddirect[mghex->get_noeud2()->get_nouveau_numero()];
298     tabnoeud[2]=lstnoeuddirect[mghex->get_noeud3()->get_nouveau_numero()];
299     tabnoeud[3]=lstnoeuddirect[mghex->get_noeud4()->get_nouveau_numero()];
300     tabnoeud[4]=lstnoeuddirect[mghex->get_noeud5()->get_nouveau_numero()];
301     tabnoeud[5]=lstnoeuddirect[mghex->get_noeud6()->get_nouveau_numero()];
302     tabnoeud[6]=lstnoeuddirect[mghex->get_noeud7()->get_nouveau_numero()];
303     tabnoeud[7]=lstnoeuddirect[mghex->get_noeud8()->get_nouveau_numero()];
304     FEM_HEXA8* hex=new FEM_HEXA8(mghex,tabnoeud);
305     fem->ajouter_fem_element3(hex);
306     mghex=maillage->get_suivant_hexa(it_hex);
307     }
308     LISTE_MG_PENTA::iterator it_pen;
309     MG_PENTA* mgpen=maillage->get_premier_penta(it_pen);
310     while (mgpen)
311     {
312     FEM_NOEUD *tabnoeud[8];
313     tabnoeud[0]=lstnoeuddirect[mgpen->get_noeud1()->get_nouveau_numero()];
314     tabnoeud[1]=lstnoeuddirect[mgpen->get_noeud2()->get_nouveau_numero()];
315     tabnoeud[2]=lstnoeuddirect[mgpen->get_noeud3()->get_nouveau_numero()];
316     tabnoeud[3]=lstnoeuddirect[mgpen->get_noeud4()->get_nouveau_numero()];
317     tabnoeud[4]=lstnoeuddirect[mgpen->get_noeud5()->get_nouveau_numero()];
318     tabnoeud[5]=lstnoeuddirect[mgpen->get_noeud6()->get_nouveau_numero()];
319     FEM_PENTA6* pen=new FEM_PENTA6(mgpen,tabnoeud);
320     fem->ajouter_fem_element3(pen);
321     mgpen=maillage->get_suivant_penta(it_pen);
322     }
323     }
324 couturad 968 return OK;
325 francois 883 }
326    
327    
328 couturad 968 int MAILLEUR_FEM::recal_element_quadratique(FEM_MAILLAGE* fem)
329 francois 915 {
330 couturad 933 std::map<long,FEM_ELEMENT3*> map_tet_replace;
331 francois 1180 int size_map_tet_replace=map_tet_replace.size();
332 couturad 933 int nbpas=(int)param.get_valeur("Quadratisation_pas");
333     double dismin=param.get_valeur("Quadratisation_dis");
334     int nbpassemax=(int)param.get_valeur("Quadratisation_nbpassemax");
335     int nbpasse=1;
336     for (int passe=0;passe<nbpasse;passe++)
337     {
338     char mess[255];
339     sprintf(mess," Passe %d",passe+1);
340     affiche(mess);
341     LISTE_FEM_ELEMENT3::iterator itele3;
342 francois 1180 //bool recale=false;
343 couturad 933 for (FEM_ELEMENT3* tet=fem->get_premier_element3(itele3);tet;tet=fem->get_suivant_element3(itele3))
344     {
345     FEM_MAILLAGE_QUADRATIQUE_OUTILS ou;
346     double dis=ou.get_distorsion2(tet);
347     double jmin=ou.get_jmin(tet);;
348     if ((dis<dismin-1e-12) || (jmin<0.))
349 francois 1180 {
350     OT_VECTEUR_3D vec[10];
351     OT_VECTEUR_3D depart[10];
352     for (int i=1;i<9;i++)
353     {
354     if ((i==2) || (i==4)) continue;
355     double* xyzori=tet->get_fem_noeud(i)->get_coord_ori();
356     double* xyz=tet->get_fem_noeud(i)->get_coord();
357     vec[i].change_x(xyz[0]-xyzori[0]);
358     vec[i].change_y(xyz[1]-xyzori[1]);
359     vec[i].change_z(xyz[2]-xyzori[2]);
360     depart[i].change_x(xyz[0]);
361     depart[i].change_y(xyz[1]);
362     depart[i].change_z(xyz[2]);
363     }
364     int pas=0;
365     while ((dis<dismin-1e-12)||(jmin<0.))
366     {
367     for (int i=1;i<9;i++)
368     {
369     if ((i==2) || (i==4)) continue;
370     double* xyz=tet->get_fem_noeud(i)->get_coord();
371     tet->get_fem_noeud(i)->change_x(depart[i].get_x()-pas*1./1./nbpas*vec[i].get_x());
372     tet->get_fem_noeud(i)->change_y(depart[i].get_y()-pas*1./1./nbpas*vec[i].get_y());
373     tet->get_fem_noeud(i)->change_z(depart[i].get_z()-pas*1./1./nbpas*vec[i].get_z());
374     }
375     dis=ou.get_distorsion2(tet);
376     jmin=ou.get_jmin(tet);
377     pas++;
378     }
379     map_tet_replace.insert(std::pair<long,FEM_ELEMENT3*>(tet->get_id(),tet));
380     if (passe==0) sprintf(mess," Tetra %lu replacé à %.2lf%% de sa position optimale\t jmin=%le distorsion=%le",tet->get_id(),100.-(pas-1.)*1./nbpas*100.,jmin,dis);
381     else sprintf(mess," Tetra %lu replacé à %.2lf%% de sa position precedente\t jmin=%le distorsion=%le",tet->get_id(),100.-(pas-1.)*1./nbpas*100.,jmin,dis);
382     affiche(mess);
383     }
384    
385 couturad 933 }
386 francois 1180 if (size_map_tet_replace!=map_tet_replace.size()) nbpasse++;
387     size_map_tet_replace=map_tet_replace.size();
388 couturad 933 if (nbpasse>nbpassemax) nbpasse=nbpassemax;
389     }
390 francois 963 if (fem->get_nb_fem_element3()!=0)
391     {
392 couturad 933 long nb_tetra_replace = map_tet_replace.size();
393     char mess[255];
394     sprintf(mess," Nb tetra replace : %li (%lf%%)",nb_tetra_replace,((double)nb_tetra_replace/fem->get_nb_fem_element3())*(double)100.0);
395     affiche(mess);
396 francois 963 }
397 couturad 968 return OK;
398 francois 915 }
399 francois 883
400    
401 francois 1164 int MAILLEUR_FEM::construire_quadratique(FEM_MAILLAGE* fem,TPL_MAP_ENTITE<MG_SEGMENT*> &mini,MG_MAILLAGE* maillage,MG_GEOMETRIE* mg_geometrie,bool courbure_discrete)
402 francois 883 {
403 francois 915 int recal=(int)param.get_valeur("Quadratisation_jmin");
404 francois 1164 int dimsansgeo=0;
405     bool virtuel =false;
406     if (mg_geometrie!=NULL)
407     if (strcmp(mg_geometrie->get_type_geometrie(),"VIRTUEL")==0) virtuel=true;
408     if (virtuel)
409     {
410     if (maillage->get_nb_mg_tetra()+maillage->get_nb_mg_hexa()!=0) dimsansgeo=3;
411     else if (maillage->get_nb_mg_triangle()+maillage->get_nb_mg_quadrangle()!=0) dimsansgeo=2;
412     else if (maillage->get_nb_mg_segment()!=0) dimsansgeo=1;
413     }
414 francois 883 LISTE_MG_NOEUD::iterator it_nd;
415     MG_NOEUD* mgnoeud=maillage->get_premier_noeud(it_nd);
416 francois 1164 if (courbure_discrete) prepare_quad_discrete(fem->get_mg_maillage(),virtuel);
417 francois 883 std::vector<FEM_NOEUD*> lstnoeuddirect;
418     unsigned int i=0;
419     while (mgnoeud)
420     {
421     mgnoeud->change_nouveau_numero(i);
422     FEM_NOEUD* femnoeud=new FEM_NOEUD(mgnoeud);
423     femnoeud->change_numero(i+1);
424     i++ ;
425     fem->ajouter_fem_noeud(femnoeud);
426     //mgnoeud=maillage->get_suivant_noeud(it_nd);
427     lstnoeuddirect.insert(lstnoeuddirect.end(),femnoeud);
428 francois 1164 if (mgnoeud->get_lien_topologie()!=NULL)
429     if (mgnoeud->get_lien_topologie()->get_dimension()==0)
430     {
431     FEM_ELEMENT0* femele=new FEM_ELEMENT0(mgnoeud->get_lien_topologie(),mgnoeud,&femnoeud);
432     fem->ajouter_fem_element0(femele);
433     }
434     mgnoeud=maillage->get_suivant_noeud(it_nd);
435 francois 883
436     }
437     if (mg_geometrie==NULL)
438     {
439     if (maillage->get_nb_mg_tetra()!=0) dimsansgeo=3;
440     else if (maillage->get_nb_mg_triangle()!=0) dimsansgeo=2;
441     else if (maillage->get_nb_mg_segment()!=0) dimsansgeo=1;
442     }
443 francois 1164 int nbmgnoeud = maillage->get_nb_mg_noeud();
444    
445 francois 883 LISTE_MG_SEGMENT::iterator it_seg;
446     MG_SEGMENT* mgseg=maillage->get_premier_segment(it_seg);
447     FEM_MAILLAGE_QUADRATIQUE_OUTILS ou ;
448     i=0;
449     while (mgseg)
450     {
451     mgseg->change_nouveau_numero(i++);
452 francois 1180 MG_NOEUD* no1t=mgseg->get_noeud1();
453     MG_NOEUD* no2t=mgseg->get_noeud2();
454 francois 883 MG_NOEUD *noeud1=mgseg->get_noeud1();
455     MG_NOEUD *noeud2=mgseg->get_noeud2();
456     double *xyz1=noeud1->get_coord();
457     double *xyz2=noeud2->get_coord();
458     double xori=0.5*(xyz1[0]+xyz2[0]);
459     double yori=0.5*(xyz1[1]+xyz2[1]);
460     double zori=0.5*(xyz1[2]+xyz2[2]);
461 ghazal 1152 FEM_NOEUD* femnoeud;
462 francois 883 if (mgseg->get_lien_topologie()!=NULL)
463     if (mgseg->get_lien_topologie()->get_dimension()==1)
464     {
465 francois 1164 double L3=0,L4=0,j=0.0;
466 ghazal 1152 double xyz[3]={xori,yori,zori};
467     if (!virtuel)
468     {
469 francois 1164 MG_ARETE* arete=(MG_ARETE*)(mgseg->get_lien_topologie());
470     double t1,t2,tm;
471     arete->inverser(t1,xyz1);
472     arete->inverser(t2,xyz2);
473     if (arete->get_courbe()->est_periodique())
474     if (t1>t2) t2=t2+arete->get_courbe()->get_periode();
475     double t=0.5*(t1+t2);
476     double xyz3[3];
477     arete->evaluer(t,xyz);
478     double L1= ou.get_distance_curviligne(0,0.5,xyz1,xyz,xyz2);
479     double L2= ou.get_distance_curviligne(0.5,1,xyz1,xyz,xyz2);
480     if (L1/L2 >1.01)
481     {
482     do
483     {
484     j+=0.0001;
485     tm=(-t+t1)*j+t;
486     arete->evaluer(tm,xyz3);
487     L3= ou.get_distance_curviligne(0,0.5,xyz1,xyz3,xyz2);
488     L4= ou.get_distance_curviligne(0.5,1,xyz1,xyz3,xyz2);
489    
490     }while (0.5-(L4/(L3+L4))>=0.0001);
491    
492     femnoeud= new FEM_NOEUD(mgseg->get_lien_topologie(),xyz3[0],xyz3[1],xyz3[2],xori,yori,zori);
493     }
494    
495     if (L2/L1 >1.01)
496     {
497     do
498     {
499     j+=0.0001;
500     tm=(-t+t2)*j+t;
501     arete->evaluer(tm,xyz3);
502     L3= ou.get_distance_curviligne(0,0.5,xyz1,xyz3,xyz2);
503     L4= ou.get_distance_curviligne(0.5,1,xyz1,xyz3,xyz2);
504    
505     }while (0.5-(L3/(L3+L4))>=0.0001);
506    
507     femnoeud= new FEM_NOEUD(mgseg->get_lien_topologie(),xyz3[0],xyz3[1],xyz3[2],xori,yori,zori);
508     }
509     }
510     if(j == 0) // geo virtuel
511     {
512     if (courbure_discrete) get_quad_seg(mgseg,virtuel,xyz);
513     femnoeud= new FEM_NOEUD(mgseg->get_lien_topologie(),xyz[0],xyz[1],xyz[2],xori,yori,zori);
514     }
515 francois 1095
516 francois 883 }
517     if (mgseg->get_lien_topologie()!=NULL)
518 francois 1164 if (mgseg->get_lien_topologie()->get_dimension()==2)
519     {
520     double L3=0,L4=0,j=0;
521     double xyz[3]={xori,yori,zori};
522     if (!virtuel)
523     {
524     MG_FACE* face=(MG_FACE*)(mgseg->get_lien_topologie());
525     double uv1[2],uv2[2];
526     OT_DECALAGE_PARAMETRE decalage(face->get_surface()->get_periode_u(),face->get_surface()->get_periode_v());
527     face->inverser(uv1,xyz1);
528     face->inverser(uv2,xyz2);
529     if (face->get_surface()->est_periodique_u()==1)
530     {
531     double eps=1e-10*face->get_surface()->get_periode_u();
532     if (uv1[0]<-eps) uv1[0]=uv1[0]+face->get_surface()->get_periode_u();
533     if (uv1[0]>face->get_surface()->get_periode_u()-eps) uv1[0]=uv1[0]-face->get_surface()->get_periode_u();
534     if (uv2[0]<-eps) uv2[0]=uv2[0]+face->get_surface()->get_periode_u();
535     if (uv2[0]>face->get_surface()->get_periode_u()-eps) uv2[0]=uv2[0]-face->get_surface()->get_periode_u();
536     }
537     if (face->get_surface()->est_periodique_v()==1)
538     {
539     double eps=1e-10*face->get_surface()->get_periode_v();
540     if (uv1[1]<-eps) uv1[1]=uv1[1]+face->get_surface()->get_periode_v();
541     if (uv1[1]>face->get_surface()->get_periode_v()-eps) uv1[1]=uv1[1]-face->get_surface()->get_periode_v();
542     if (uv2[1]<-eps) uv2[1]=uv2[1]+face->get_surface()->get_periode_v();
543     if (uv2[1]>face->get_surface()->get_periode_v()-eps) uv2[1]=uv2[1]-face->get_surface()->get_periode_v();
544     }
545     double du=decalage.calcul_decalage_parametre_u(uv1[0]);
546     double dv=decalage.calcul_decalage_parametre_v(uv1[1]);
547     double u1=decalage.decalage_parametre_u(uv1[0],du);
548     double v1=decalage.decalage_parametre_v(uv1[1],dv);
549     double u2=decalage.decalage_parametre_u(uv2[0],du);
550     double v2=decalage.decalage_parametre_v(uv2[1],dv);
551     double uv[2],uv3[2],xyz3[3];
552     double milieu[3];
553     milieu[0]=0.5*(xyz1[0]+xyz2[0]);
554     milieu[1]=0.5*(xyz1[1]+xyz2[1]);
555     milieu[2]=0.5*(xyz1[2]+xyz2[2]);
556     OT_VECTEUR_3D vec(xyz1,xyz2);
557     face->inverser(uv,milieu,0.25*vec.get_longueur());
558     face->evaluer(uv,xyz);
559     double u=decalage.decalage_parametre_u(uv[0],du);
560     double v=decalage.decalage_parametre_v(uv[1],dv);
561    
562     double L1= ou.get_distance_curviligne(0,0.5,xyz1,xyz,xyz2);
563     double L2= ou.get_distance_curviligne(0.5,1,xyz1,xyz,xyz2);
564     if (L1/L2 >1.01)
565     {
566     do
567     {
568     j+=0.0001;
569     if ((u1<u) && (u2<u) && (u1<u2)) u1=u1+face->get_surface()->get_periode_u();
570     if ((v1<v) && (v2<v) && (v1<v2)) v1=v1+face->get_surface()->get_periode_v();
571    
572     double u3=(-u+u1)*j+u;
573     double v3=(-v+v1)*j+v;
574     uv3[0]=decalage.decalage_parametre_u(u3,-du);
575     uv3[1]=decalage.decalage_parametre_v(v3,-dv);
576     face->evaluer(uv3,xyz3);
577     L3= ou.get_distance_curviligne(0,0.5,xyz1,xyz3,xyz2);
578     L4= ou.get_distance_curviligne(0.5,1,xyz1,xyz3,xyz2);
579    
580     }while ( 0.5-(L4/(L3+L4)) >=0.0001);
581    
582     femnoeud= new FEM_NOEUD(mgseg->get_lien_topologie(),xyz3[0],xyz3[1],xyz3[2],xori,yori,zori);
583     }
584    
585     if (L2/L1 >1.01)
586     {
587     do
588     {
589     j+=0.0001;
590     if ((u1<u) && (u2<u) && (u2<u1)) u2=u2+face->get_surface()->get_periode_u();
591     if ((v1<v) && (v2<v) && (v2<v1)) v2=v2+face->get_surface()->get_periode_v();
592    
593     double u3=(-u+u2)*j+u;
594     double v3=(-v+v2)*j+v;
595     uv3[0]=decalage.decalage_parametre_u(u3,-du);
596     uv3[1]=decalage.decalage_parametre_v(v3,-dv);
597    
598     face->evaluer(uv3,xyz3);
599     L3= ou.get_distance_curviligne(0,0.5,xyz1,xyz3,xyz2);
600     L4= ou.get_distance_curviligne(0.5,1,xyz1,xyz3,xyz2);
601    
602     }while (0.5-(L3/(L3+L4)) >=0.0001);
603    
604     femnoeud= new FEM_NOEUD(mgseg->get_lien_topologie(),xyz3[0],xyz3[1],xyz3[2],xori,yori,zori);
605     }
606     }
607     if (j == 0.0) // geo virtuel
608     {
609     if (courbure_discrete) get_quad_seg(mgseg,virtuel,xyz);
610     femnoeud= new FEM_NOEUD(mgseg->get_lien_topologie(),xyz[0],xyz[1],xyz[2],xori,yori,zori);
611     }
612    
613     }
614     if (mgseg->get_lien_topologie()!=NULL)
615     {
616     if (mgseg->get_lien_topologie()->get_dimension()==3)
617     {
618     double x=0.5*(xyz1[0]+xyz2[0]);
619     double y=0.5*(xyz1[1]+xyz2[1]);
620     double z=0.5*(xyz1[2]+xyz2[2]);
621     femnoeud= new FEM_NOEUD(mgseg->get_lien_topologie(),x,y,z,x,y,z);
622     }
623     }
624     else // pas de geo
625     {
626     double xyz[3];
627     xyz[0]=xori;
628     xyz[1]=yori;
629     xyz[2]=zori;
630     if (courbure_discrete) get_quad_seg(mgseg,virtuel,xyz);
631     femnoeud= new FEM_NOEUD(mgseg->get_lien_topologie(),xyz[0],xyz[1],xyz[2],xori,yori,zori);
632     }
633     femnoeud->change_numero(i+nbmgnoeud+1);
634     fem->ajouter_fem_noeud(femnoeud);
635     mgseg=maillage->get_suivant_segment(it_seg);
636     lstnoeuddirect.insert(lstnoeuddirect.end(),femnoeud);
637 francois 883 }
638    
639     if ((mg_geometrie!=NULL) || ( (mg_geometrie==NULL) && (dimsansgeo==1)))
640     {
641 francois 1164
642 francois 883 mgseg=maillage->get_premier_segment(it_seg);
643     while (mgseg)
644     {
645     if (mgseg->get_lien_topologie()!=NULL)
646 francois 1164 if (mgseg->get_lien_topologie()->get_dimension()!=1)
647     {
648 francois 883 mgseg=maillage->get_suivant_segment(it_seg);
649     continue;
650 francois 1164 }
651     if (mg_geometrie!=NULL)
652     if (mgseg->get_lien_topologie()==NULL)
653     if (dimsansgeo!=1)
654     {
655     mgseg=maillage->get_suivant_segment(it_seg);
656     continue;
657     }
658     FEM_NOEUD *tabnoeud[3];
659     tabnoeud[0]=lstnoeuddirect[mgseg->get_noeud1()->get_nouveau_numero()];
660     tabnoeud[1]=lstnoeuddirect[mgseg->get_nouveau_numero()+nbmgnoeud];
661     tabnoeud[2]=lstnoeuddirect[mgseg->get_noeud2()->get_nouveau_numero()];
662     FEM_SEGMENT3* seg=new FEM_SEGMENT3(mgseg,tabnoeud);
663     fem->ajouter_fem_element1(seg);
664 francois 883 mgseg=maillage->get_suivant_segment(it_seg);
665     }
666     }
667     if ((mg_geometrie!=NULL) || ( (mg_geometrie==NULL) && (dimsansgeo==2)))
668     {
669     LISTE_MG_TRIANGLE::iterator it_tri;
670     MG_TRIANGLE* mgtri=maillage->get_premier_triangle(it_tri);
671     while (mgtri)
672     {
673     if (mgtri->get_lien_topologie()!=NULL) if (mgtri->get_lien_topologie()->get_dimension()!=2)
674 francois 1164 {
675     mgtri=maillage->get_suivant_triangle(it_tri);
676     continue;
677     }
678     if (mg_geometrie!=NULL)
679     if (mgtri->get_lien_topologie()==NULL)
680     if (dimsansgeo!=2)
681     {
682     mgtri=maillage->get_suivant_triangle(it_tri);
683     continue;
684 francois 883 }
685 francois 1164 FEM_NOEUD *tabnoeud[6];
686     tabnoeud[0]=lstnoeuddirect[mgtri->get_noeud1()->get_nouveau_numero()];
687     tabnoeud[1]=lstnoeuddirect[mgtri->get_segment1()->get_nouveau_numero()+nbmgnoeud];
688     tabnoeud[2]=lstnoeuddirect[mgtri->get_noeud2()->get_nouveau_numero()];
689     tabnoeud[3]=lstnoeuddirect[mgtri->get_segment2()->get_nouveau_numero()+nbmgnoeud];
690     tabnoeud[4]=lstnoeuddirect[mgtri->get_noeud3()->get_nouveau_numero()];
691     tabnoeud[5]=lstnoeuddirect[mgtri->get_segment3()->get_nouveau_numero()+nbmgnoeud];
692     FEM_TRIANGLE6* tri=new FEM_TRIANGLE6(mgtri,tabnoeud);
693     fem->ajouter_fem_element2(tri);
694 francois 883 mgtri=maillage->get_suivant_triangle(it_tri);
695     }
696     LISTE_MG_QUADRANGLE::iterator it_quad;
697     MG_QUADRANGLE* mgquad=maillage->get_premier_quadrangle(it_quad);
698     while (mgquad)
699     {
700     if (mgquad->get_lien_topologie()!=NULL) if (mgquad->get_lien_topologie()->get_dimension()!=2) {
701 francois 1164 mgquad=maillage->get_suivant_quadrangle(it_quad);
702     continue;
703     }
704 francois 883 if (mg_geometrie!=NULL)
705 francois 1164 if (mgquad->get_lien_topologie()==NULL)
706     if (dimsansgeo!=2)
707     {
708     mgquad=maillage->get_suivant_quadrangle(it_quad);
709     continue;
710     }
711     FEM_NOEUD *tabnoeud[8];
712 francois 883 tabnoeud[0]=lstnoeuddirect[mgquad->get_noeud1()->get_nouveau_numero()];
713 francois 1164 tabnoeud[1]=lstnoeuddirect[mgquad->get_segment1()->get_nouveau_numero()+nbmgnoeud];
714     tabnoeud[2]=lstnoeuddirect[mgquad->get_noeud2()->get_nouveau_numero()];
715     tabnoeud[3]=lstnoeuddirect[mgquad->get_segment2()->get_nouveau_numero()+nbmgnoeud];
716     tabnoeud[4]=lstnoeuddirect[mgquad->get_noeud3()->get_nouveau_numero()];
717     tabnoeud[5]=lstnoeuddirect[mgquad->get_segment3()->get_nouveau_numero()+nbmgnoeud];
718     tabnoeud[6]=lstnoeuddirect[mgquad->get_noeud4()->get_nouveau_numero()];
719     tabnoeud[7]=lstnoeuddirect[mgquad->get_segment4()->get_nouveau_numero()+nbmgnoeud];
720     FEM_QUADRANGLE8* quad=new FEM_QUADRANGLE8(mgquad,tabnoeud);
721     fem->ajouter_fem_element2(quad);
722 francois 883 mgquad=maillage->get_suivant_quadrangle(it_quad);
723     }
724     }
725     if ((mg_geometrie!=NULL) || ( (mg_geometrie==NULL) && (dimsansgeo==3)))
726     {
727 francois 1164
728     LISTE_MG_TETRA::iterator it_tetra;
729 francois 883 MG_TETRA* mgtetra=maillage->get_premier_tetra(it_tetra);
730     while (mgtetra)
731     {
732     FEM_NOEUD *tabnoeud[10];
733     tabnoeud[0]=lstnoeuddirect[mgtetra->get_noeud1()->get_nouveau_numero()];
734     tabnoeud[1]=lstnoeuddirect[maillage->get_mg_segment(mgtetra->get_noeud1()->get_id(),mgtetra->get_noeud2()->get_id())->get_nouveau_numero()+nbmgnoeud];
735     tabnoeud[2]=lstnoeuddirect[mgtetra->get_noeud2()->get_nouveau_numero()];
736     tabnoeud[3]=lstnoeuddirect[maillage->get_mg_segment(mgtetra->get_noeud2()->get_id(),mgtetra->get_noeud3()->get_id())->get_nouveau_numero()+nbmgnoeud];
737     tabnoeud[4]=lstnoeuddirect[mgtetra->get_noeud3()->get_nouveau_numero()];
738     tabnoeud[5]=lstnoeuddirect[maillage->get_mg_segment(mgtetra->get_noeud3()->get_id(),mgtetra->get_noeud1()->get_id())->get_nouveau_numero()+nbmgnoeud];
739     tabnoeud[6]=lstnoeuddirect[maillage->get_mg_segment(mgtetra->get_noeud1()->get_id(),mgtetra->get_noeud4()->get_id())->get_nouveau_numero()+nbmgnoeud];
740     tabnoeud[7]=lstnoeuddirect[maillage->get_mg_segment(mgtetra->get_noeud2()->get_id(),mgtetra->get_noeud4()->get_id())->get_nouveau_numero()+nbmgnoeud];
741     tabnoeud[8]=lstnoeuddirect[maillage->get_mg_segment(mgtetra->get_noeud3()->get_id(),mgtetra->get_noeud4()->get_id())->get_nouveau_numero()+nbmgnoeud];
742     tabnoeud[9]=lstnoeuddirect[mgtetra->get_noeud4()->get_nouveau_numero()];
743     FEM_TETRA10* tet=new FEM_TETRA10(mgtetra,tabnoeud);
744     fem->ajouter_fem_element3(tet);
745 francois 1164 mgtetra=maillage->get_suivant_tetra(it_tetra);
746     }
747     if (recal==1) if(recal_element_quadratique(fem)==FAIL) return FAIL;
748 francois 915
749 francois 1164
750    
751 francois 883 LISTE_MG_HEXA::iterator it_hexa;
752     MG_HEXA* mghex=maillage->get_premier_hexa(it_hexa);
753     while (mghex)
754     {
755     FEM_NOEUD *tabnoeud[20];
756     tabnoeud[0]=lstnoeuddirect[mghex->get_noeud1()->get_nouveau_numero()];
757     tabnoeud[1]=lstnoeuddirect[maillage->get_mg_segment(mghex->get_noeud1()->get_id(),mghex->get_noeud2()->get_id())->get_nouveau_numero()+nbmgnoeud];
758     tabnoeud[2]=lstnoeuddirect[mghex->get_noeud2()->get_nouveau_numero()];
759     tabnoeud[3]=lstnoeuddirect[maillage->get_mg_segment(mghex->get_noeud2()->get_id(),mghex->get_noeud3()->get_id())->get_nouveau_numero()+nbmgnoeud];
760     tabnoeud[4]=lstnoeuddirect[mghex->get_noeud3()->get_nouveau_numero()];
761     tabnoeud[5]=lstnoeuddirect[maillage->get_mg_segment(mghex->get_noeud3()->get_id(),mghex->get_noeud4()->get_id())->get_nouveau_numero()+nbmgnoeud];
762     tabnoeud[6]=lstnoeuddirect[mghex->get_noeud4()->get_nouveau_numero()];
763     tabnoeud[7]=lstnoeuddirect[maillage->get_mg_segment(mghex->get_noeud4()->get_id(),mghex->get_noeud1()->get_id())->get_nouveau_numero()+nbmgnoeud];
764     tabnoeud[8]=lstnoeuddirect[maillage->get_mg_segment(mghex->get_noeud1()->get_id(),mghex->get_noeud5()->get_id())->get_nouveau_numero()+nbmgnoeud];
765     tabnoeud[9]=lstnoeuddirect[maillage->get_mg_segment(mghex->get_noeud2()->get_id(),mghex->get_noeud6()->get_id())->get_nouveau_numero()+nbmgnoeud];
766     tabnoeud[10]=lstnoeuddirect[maillage->get_mg_segment(mghex->get_noeud3()->get_id(),mghex->get_noeud7()->get_id())->get_nouveau_numero()+nbmgnoeud];
767     tabnoeud[11]=lstnoeuddirect[maillage->get_mg_segment(mghex->get_noeud4()->get_id(),mghex->get_noeud8()->get_id())->get_nouveau_numero()+nbmgnoeud];
768     tabnoeud[12]=lstnoeuddirect[mghex->get_noeud5()->get_nouveau_numero()];
769     tabnoeud[13]=lstnoeuddirect[maillage->get_mg_segment(mghex->get_noeud5()->get_id(),mghex->get_noeud6()->get_id())->get_nouveau_numero()+nbmgnoeud];
770     tabnoeud[14]=lstnoeuddirect[mghex->get_noeud6()->get_nouveau_numero()];
771     tabnoeud[15]=lstnoeuddirect[maillage->get_mg_segment(mghex->get_noeud6()->get_id(),mghex->get_noeud7()->get_id())->get_nouveau_numero()+nbmgnoeud];
772     tabnoeud[16]=lstnoeuddirect[mghex->get_noeud7()->get_nouveau_numero()];
773     tabnoeud[17]=lstnoeuddirect[maillage->get_mg_segment(mghex->get_noeud7()->get_id(),mghex->get_noeud8()->get_id())->get_nouveau_numero()+nbmgnoeud];
774     tabnoeud[18]=lstnoeuddirect[mghex->get_noeud8()->get_nouveau_numero()];
775     tabnoeud[19]=lstnoeuddirect[maillage->get_mg_segment(mghex->get_noeud8()->get_id(),mghex->get_noeud5()->get_id())->get_nouveau_numero()+nbmgnoeud];
776     FEM_HEXA20* hex=new FEM_HEXA20(mghex,tabnoeud);
777     fem->ajouter_fem_element3(hex);
778     mghex=maillage->get_suivant_hexa(it_hexa);
779     }
780     LISTE_MG_PENTA::iterator it_pen;
781     MG_PENTA* mgpen=maillage->get_premier_penta(it_pen);
782     while (mgpen)
783     {
784     FEM_NOEUD *tabnoeud[15];
785     tabnoeud[0]=lstnoeuddirect[mgpen->get_noeud1()->get_nouveau_numero()];
786     tabnoeud[1]=lstnoeuddirect[maillage->get_mg_segment(mgpen->get_noeud1()->get_id(),mgpen->get_noeud2()->get_id())->get_nouveau_numero()+nbmgnoeud];
787     tabnoeud[2]=lstnoeuddirect[mgpen->get_noeud2()->get_nouveau_numero()];
788     tabnoeud[3]=lstnoeuddirect[maillage->get_mg_segment(mgpen->get_noeud2()->get_id(),mgpen->get_noeud3()->get_id())->get_nouveau_numero()+nbmgnoeud];
789     tabnoeud[4]=lstnoeuddirect[mgpen->get_noeud3()->get_nouveau_numero()];
790     tabnoeud[5]=lstnoeuddirect[maillage->get_mg_segment(mgpen->get_noeud3()->get_id(),mgpen->get_noeud1()->get_id())->get_nouveau_numero()+nbmgnoeud];
791     tabnoeud[6]=lstnoeuddirect[maillage->get_mg_segment(mgpen->get_noeud1()->get_id(),mgpen->get_noeud4()->get_id())->get_nouveau_numero()+nbmgnoeud];
792     tabnoeud[7]=lstnoeuddirect[maillage->get_mg_segment(mgpen->get_noeud2()->get_id(),mgpen->get_noeud5()->get_id())->get_nouveau_numero()+nbmgnoeud];
793     tabnoeud[8]=lstnoeuddirect[maillage->get_mg_segment(mgpen->get_noeud3()->get_id(),mgpen->get_noeud6()->get_id())->get_nouveau_numero()+nbmgnoeud];
794     tabnoeud[9]=lstnoeuddirect[mgpen->get_noeud4()->get_nouveau_numero()];
795     tabnoeud[10]=lstnoeuddirect[maillage->get_mg_segment(mgpen->get_noeud4()->get_id(),mgpen->get_noeud5()->get_id())->get_nouveau_numero()+nbmgnoeud];
796     tabnoeud[11]=lstnoeuddirect[mgpen->get_noeud5()->get_nouveau_numero()];
797     tabnoeud[12]=lstnoeuddirect[maillage->get_mg_segment(mgpen->get_noeud5()->get_id(),mgpen->get_noeud6()->get_id())->get_nouveau_numero()+nbmgnoeud];
798     tabnoeud[13]=lstnoeuddirect[mgpen->get_noeud6()->get_nouveau_numero()];
799     tabnoeud[14]=lstnoeuddirect[maillage->get_mg_segment(mgpen->get_noeud4()->get_id(),mgpen->get_noeud6()->get_id())->get_nouveau_numero()+nbmgnoeud];
800     FEM_PENTA15* pen=new FEM_PENTA15(mgpen,tabnoeud);
801     fem->ajouter_fem_element3(pen);
802     mgpen=maillage->get_suivant_penta(it_pen);
803     }
804     }
805 couturad 968 return OK;
806 francois 883 }
807    
808    
809    
810    
811     void MAILLEUR_FEM::optimise_numerotation(FEM_MAILLAGE* fem)
812     {
813     FEM_NOEUD* noeud=fem->get_fem_noeud(0);
814     FEM_GRAPH_NOEUD *graph;
815     graph=new FEM_GRAPH_NOEUD(noeud,fem);
816     int e=graph->get_excentricite();
817     TPL_MAP_ENTITE<FEM_NOEUD*> dernierniveau=graph->get_dernier_niveau();
818     int nbdernierniveau=dernierniveau.get_nb();
819     for (int i=0;i<nbdernierniveau;i++)
820     {
821     FEM_GRAPH_NOEUD* graphtmp=new FEM_GRAPH_NOEUD(dernierniveau.get(i),fem);
822     int etmp=graphtmp->get_excentricite();
823     if (etmp>e)
824     {
825     delete graph;
826     graph=graphtmp;
827     dernierniveau=graph->get_dernier_niveau();
828     nbdernierniveau=dernierniveau.get_nb();
829     i=-1;
830     e=etmp;
831     }
832     else delete graphtmp;
833     }
834     int numero=fem->get_nb_fem_noeud();
835     for (int i=0;i<e+1;i++)
836     {
837     TPL_MAP_ENTITE<FEM_NOEUD*> niveau=graph->get_niveau(i);
838     int nbnoeud=niveau.get_nb();
839     for (int j=0;j<nbnoeud;j++)
840     {
841     FEM_NOEUD* noeud=niveau.get(j);
842     noeud->change_numero_opt(numero);
843     numero--;
844     }
845     }
846    
847    
848     delete graph;
849     }
850    
851    
852     void MAILLEUR_FEM::copie_numerotation_opt(FEM_MAILLAGE* fem)
853     {
854     int nbnoeud=fem->get_nb_fem_noeud();
855     for (int i=0;i<nbnoeud;i++)
856     {
857     FEM_NOEUD* noeud=fem->get_fem_noeud(i);
858     noeud->change_numero(noeud->get_numero_opt());
859     }
860    
861     }
862    
863 francois 1164 void MAILLEUR_FEM::prepare_quad_discrete(MG_MAILLAGE* mai,bool virtuel)
864     {
865     MG_GEOMETRIE *geo=mai->get_mg_geometrie();
866 francois 1180 liste_normale.clear();
867     segaprojete.clear();
868     segnonprojete.clear();
869 francois 1164 LISTE_MG_TETRA::iterator it;
870     LISTE_MG_NOEUD::iterator itn;
871     LISTE_MG_TRIANGLE::iterator itt;
872     LISTE_MG_SEGMENT::iterator its;
873     for (MG_TETRA* tet=mai->get_premier_tetra(it);tet!=NULL;tet=mai->get_suivant_tetra(it))
874     {
875     tet->change_nouveau_numero(0);
876     tet->get_triangle1()->change_nouveau_numero(0);
877     tet->get_triangle2()->change_nouveau_numero(0);
878     tet->get_triangle3()->change_nouveau_numero(0);
879     tet->get_triangle4()->change_nouveau_numero(0);
880     }
881     for (MG_TETRA* tet=mai->get_premier_tetra(it);tet!=NULL;tet=mai->get_suivant_tetra(it))
882     {
883     if (tet->get_lien_topologie()!=NULL)
884     {
885     if (tet->get_triangle1()->get_nouveau_numero()==0)
886     tet->get_triangle1()->change_nouveau_numero(tet->get_lien_topologie()->get_id());
887     else if (tet->get_triangle1()->get_nouveau_numero()>2)
888     {
889     if (tet->get_triangle1()->get_nouveau_numero()==tet->get_lien_topologie()->get_id())
890     tet->get_triangle1()->change_nouveau_numero(2);
891     else
892     tet->get_triangle1()->change_nouveau_numero(1);
893     }
894    
895     }
896     else tet->get_triangle1()->change_nouveau_numero(tet->get_triangle1()->get_nouveau_numero()+1);
897    
898     if (tet->get_lien_topologie()!=NULL)
899     {
900     if (tet->get_triangle2()->get_nouveau_numero()==0)
901     tet->get_triangle2()->change_nouveau_numero(tet->get_lien_topologie()->get_id());
902     else if (tet->get_triangle2()->get_nouveau_numero()>2)
903     {
904     if (tet->get_triangle2()->get_nouveau_numero()==tet->get_lien_topologie()->get_id())
905     tet->get_triangle2()->change_nouveau_numero(2);
906     else
907     tet->get_triangle2()->change_nouveau_numero(1);
908     }
909    
910     }
911     else tet->get_triangle2()->change_nouveau_numero(tet->get_triangle2()->get_nouveau_numero()+1);
912    
913     if (tet->get_lien_topologie()!=NULL)
914     {
915     if (tet->get_triangle3()->get_nouveau_numero()==0)
916     tet->get_triangle3()->change_nouveau_numero(tet->get_lien_topologie()->get_id());
917     else if (tet->get_triangle3()->get_nouveau_numero()>2)
918     {
919     if (tet->get_triangle3()->get_nouveau_numero()==tet->get_lien_topologie()->get_id())
920     tet->get_triangle3()->change_nouveau_numero(2);
921     else
922     tet->get_triangle3()->change_nouveau_numero(1);
923     }
924    
925     }
926     else tet->get_triangle3()->change_nouveau_numero(tet->get_triangle3()->get_nouveau_numero()+1);
927    
928     if (tet->get_lien_topologie()!=NULL)
929     {
930     if (tet->get_triangle4()->get_nouveau_numero()==0)
931     tet->get_triangle4()->change_nouveau_numero(tet->get_lien_topologie()->get_id());
932     else if (tet->get_triangle4()->get_nouveau_numero()>2)
933     {
934     if (tet->get_triangle4()->get_nouveau_numero()==tet->get_lien_topologie()->get_id())
935     tet->get_triangle4()->change_nouveau_numero(2);
936     else
937     tet->get_triangle4()->change_nouveau_numero(1);
938     }
939    
940     }
941     else tet->get_triangle4()->change_nouveau_numero(tet->get_triangle4()->get_nouveau_numero()+1);
942    
943     }
944     if (geo!=NULL)
945     {
946     LISTE_MG_FACE::iterator itf;
947     for (MG_FACE* face=geo->get_premier_face(itf);face!=NULL;face=geo->get_suivant_face(itf))
948     {
949     int nbele=face->get_lien_maillage()->get_nb();
950     for (int i=0;i<nbele;i++)
951     {
952     MG_ELEMENT_MAILLAGE* ele=face->get_lien_maillage()->get(i);
953     ele->change_nouveau_numero(1);
954     }
955    
956     }
957     }
958     for (MG_SEGMENT* seg=mai->get_premier_segment(its);seg!=NULL;seg=mai->get_suivant_segment(its))
959     seg->change_nouveau_numero(0);
960     for (MG_TRIANGLE* tri=mai->get_premier_triangle(itt);tri!=NULL;tri=mai->get_suivant_triangle(itt))
961     {
962     if (tri->get_nouveau_numero()>2) tri->change_nouveau_numero(1);
963     if (tri->get_lien_topologie()!=NULL)
964     if (tri->get_lien_topologie())
965     if (virtuel)
966     if ((((MG_FACE_ELEMENT*)tri->get_lien_topologie())->get_nb_contrainte()>0))
967     {
968     segnonprojete[tri->get_segment1()]=tri->get_segment1();
969     segnonprojete[tri->get_segment2()]=tri->get_segment2();
970     segnonprojete[tri->get_segment3()]=tri->get_segment3();
971     }
972     if (tri->get_nouveau_numero()==1)
973     {
974     segaprojete[tri->get_segment1()]=tri->get_segment1();
975     segaprojete[tri->get_segment2()]=tri->get_segment2();
976     segaprojete[tri->get_segment3()]=tri->get_segment3();
977     }
978     }
979     }
980    
981     bool MAILLEUR_FEM::get_courbure_noeud(MG_NOEUD* no,MG_NOEUD* dirno,OT_VECTEUR_3D &n,double &r)
982     {
983     OT_VECTEUR_3D PP1(no->get_coord(),dirno->get_coord());
984 francois 1180 if (liste_normale.count(no)==0)
985 francois 1164 {
986 francois 1180 double deno=0.;
987     OT_VECTEUR_3D n0(0.,0.,0.);
988     n=n0;
989     int passe=0;
990     std::map<MG_NOEUD*,std::pair<int,int> > etoilesegment;
991     std::map<MG_NOEUD*,std::pair<MG_TRIANGLE*,MG_TRIANGLE*> > etoilesegmenttriangle;
992    
993     for (int i=0;i<no->get_lien_triangle()->get_nb();i++)
994     {
995     MG_TRIANGLE *tri=no->get_lien_triangle()->get(i);
996     if (tri->get_nouveau_numero()!=1) continue;
997     MG_NOEUD* not1=tri->get_noeud1();
998     MG_NOEUD* not2=tri->get_noeud2();
999     MG_NOEUD* not3=tri->get_noeud3();
1000     passe++;
1001     MG_NOEUD *nobase1,*nobase2;
1002     if (not1==no) {nobase1=not2;nobase2=not3;}
1003     if (not2==no) {nobase1=not3;nobase2=not1;}
1004     if (not3==no) {nobase1=not1;nobase2=not2;}
1005     if (etoilesegment.count(nobase1)==0) etoilesegment[nobase1].first=1; else etoilesegment[nobase1].second=1;
1006     if (etoilesegment.count(nobase2)==0) etoilesegment[nobase2].first=-1; else etoilesegment[nobase2].second=-1;
1007     if (etoilesegmenttriangle.count(nobase1)==0) etoilesegmenttriangle[nobase1].first=tri; else etoilesegmenttriangle[nobase1].second=tri;
1008     if (etoilesegmenttriangle.count(nobase2)==0) etoilesegmenttriangle[nobase2].first=tri; else etoilesegmenttriangle[nobase2].second=tri;
1009     }
1010     if (passe==0) return false;
1011    
1012     std::map<MG_NOEUD*,std::pair<int,int> > :: iterator itseg=etoilesegment.begin();
1013     MG_NOEUD* segdepart=itseg->first;
1014     MG_NOEUD* segcour=itseg->first;
1015     MG_TRIANGLE* tricour=NULL;
1016     int sens=1;
1017     do
1018     {
1019     if (etoilesegmenttriangle[segcour].first!=tricour) tricour=etoilesegmenttriangle[segcour].first; else tricour=etoilesegmenttriangle[segcour].second;
1020     if (etoilesegment[segcour].first*etoilesegment[segcour].second>0) sens=-sens;
1021     MG_NOEUD* not1=tricour->get_noeud1();
1022     MG_NOEUD* not2=tricour->get_noeud2();
1023     MG_NOEUD* not3=tricour->get_noeud3();
1024     OT_VECTEUR_3D vec1(not1->get_coord(),not2->get_coord());
1025     OT_VECTEUR_3D vec2(not1->get_coord(),not3->get_coord());
1026     OT_VECTEUR_3D pvec=vec1&vec2;
1027     double aire=pvec.get_longueur();
1028     pvec.norme();
1029     pvec=sens*pvec;
1030     n=n+aire*pvec;
1031     deno=deno+aire;
1032     if ((tricour->get_noeud1()!=no) && (tricour->get_noeud1()!=segcour)) segcour=tricour->get_noeud1();
1033     else if ((tricour->get_noeud2()!=no) && (tricour->get_noeud2()!=segcour)) segcour=tricour->get_noeud2();
1034     else if ((tricour->get_noeud3()!=no) && (tricour->get_noeud3()!=segcour)) segcour=tricour->get_noeud3();
1035     }
1036     while (segcour!=segdepart);
1037     n=(1./deno)*n;
1038     n.norme();
1039     liste_normale[no]=n;
1040 francois 1164 }
1041 francois 1180 else n=liste_normale[no];
1042 francois 1179 r=PP1*PP1/2./(n*PP1);
1043 francois 1164 return true;
1044     }
1045    
1046     void MAILLEUR_FEM::get_quad_seg(MG_SEGMENT *seg,bool virtuel,double *xyz)
1047     {
1048     if (virtuel)
1049     if (seg->get_lien_topologie()!=NULL)
1050     if (((MG_ARETE_ELEMENT*)seg->get_lien_topologie())->get_nb_contrainte()>0) return;
1051     if (segaprojete[seg]==NULL)
1052     return;
1053     if (segnonprojete[seg]!=NULL)
1054     return;
1055     MG_NOEUD* no1=seg->get_noeud1();
1056     MG_NOEUD* no2=seg->get_noeud2();
1057     OT_VECTEUR_3D n1,n2;
1058     double r1,r2;
1059     bool def1=get_courbure_noeud(no1,no2,n1,r1);
1060     bool def2=get_courbure_noeud(no2,no1,n2,r2);
1061     if (def1==false) return;
1062     if (def2==false) return;
1063     double r=0.5*(r1+r2);
1064     OT_VECTEUR_3D n=0.5*(n1+n2);
1065 francois 1168 n.norme();
1066 francois 1164 double dist=0.;
1067     if (fabs(r)<1000.*seg->get_longueur())
1068     {
1069     dist=r*r-0.25*seg->get_longueur()*seg->get_longueur();
1070 francois 1170 dist=fabs(r)-sqrt(dist);
1071 francois 1179 if (r>0) dist=-dist;
1072 francois 1164 }
1073     xyz[0]=xyz[0]+dist*n.get_x();
1074     xyz[1]=xyz[1]+dist*n.get_y();
1075     xyz[2]=xyz[2]+dist*n.get_z();
1076     }