ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/mailleur/src/mailleur_fem.cpp
Revision: 1182
Committed: Thu Nov 7 21:43:35 2024 UTC (6 months ago) by ghazal
File size: 50680 byte(s)
Log Message:
résolution du calcul de courbure pour un noeud entre fibre et face

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