ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/mailleur/src/mailleur_fem.cpp
Revision: 1168
Committed: Thu Jul 25 22:21:32 2024 UTC (9 months, 2 weeks ago) by francois
File size: 48009 byte(s)
Log Message:
bur dans le maillage quadratique discret et mailleur affichage de -examine pour les mg_cg_assemblage

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     int nbpas=(int)param.get_valeur("Quadratisation_pas");
332     double dismin=param.get_valeur("Quadratisation_dis");
333     int nbpassemax=(int)param.get_valeur("Quadratisation_nbpassemax");
334     int nbpasse=1;
335     for (int passe=0;passe<nbpasse;passe++)
336     {
337     char mess[255];
338     sprintf(mess," Passe %d",passe+1);
339     affiche(mess);
340     LISTE_FEM_ELEMENT3::iterator itele3;
341     bool recale=false;
342     for (FEM_ELEMENT3* tet=fem->get_premier_element3(itele3);tet;tet=fem->get_suivant_element3(itele3))
343     {
344     FEM_MAILLAGE_QUADRATIQUE_OUTILS ou;
345     double dis=ou.get_distorsion2(tet);
346     double jmin=ou.get_jmin(tet);;
347     if ((dis<dismin-1e-12) || (jmin<0.))
348     {
349     OT_VECTEUR_3D vec[10];
350     OT_VECTEUR_3D depart[10];
351 francois 1006 for (int i=1;i<9;i++)
352 francois 915 {
353 francois 1006 if ((i==2) || (i==4)) continue;
354 couturad 933 double* xyzori=tet->get_fem_noeud(i)->get_coord_ori();
355     double* xyz=tet->get_fem_noeud(i)->get_coord();
356     vec[i].change_x(xyz[0]-xyzori[0]);
357     vec[i].change_y(xyz[1]-xyzori[1]);
358     vec[i].change_z(xyz[2]-xyzori[2]);
359     depart[i].change_x(xyz[0]);
360     depart[i].change_y(xyz[1]);
361     depart[i].change_z(xyz[2]);
362     }
363     int pas=0;
364     while ((dis<dismin-1e-12)||(jmin<0.))
365     {
366 francois 1006 for (int i=1;i<9;i++)
367 couturad 933 {
368 francois 1006 if ((i==2) || (i==4)) continue;
369 couturad 933 double* xyz=tet->get_fem_noeud(i)->get_coord();
370     tet->get_fem_noeud(i)->change_x(depart[i].get_x()-pas*1./1./nbpas*vec[i].get_x());
371     tet->get_fem_noeud(i)->change_y(depart[i].get_y()-pas*1./1./nbpas*vec[i].get_y());
372     tet->get_fem_noeud(i)->change_z(depart[i].get_z()-pas*1./1./nbpas*vec[i].get_z());
373     }
374     dis=ou.get_distorsion2(tet);
375     jmin=ou.get_jmin(tet);
376     pas++;
377 francois 1006 /* if(pas>nbpas)
378 couturad 968 {
379     std::cerr << "*** MAILLEUR_FEM::recal_element_quadratique : pas>nbpas ! ***" << std::endl;
380     fem->get_mg_maillage()->get_gestionnaire()->enregistrer("void.magic");
381     return FAIL;
382 francois 1006 }*/
383 couturad 933 }
384     map_tet_replace.insert(std::pair<long,FEM_ELEMENT3*>(tet->get_id(),tet));
385     if (passe==0) sprintf(mess," Tetra %lu replacé à %.2lf%% de sa position optimale\t jmin=%le distorsion=%le",tet->get_id(),(pas-1.)*1./nbpas*100.,jmin,dis);
386     else sprintf(mess," Tetra %lu replacé à %.2lf%% de sa position precedente\t jmin=%le distorsion=%le",tet->get_id(),(pas-1.)*1./nbpas*100.,jmin,dis);
387     affiche(mess);
388     recale=true;
389     }
390     }
391     nbpasse++;
392     if (nbpasse>nbpassemax) nbpasse=nbpassemax;
393     }
394 francois 963 if (fem->get_nb_fem_element3()!=0)
395     {
396 couturad 933 long nb_tetra_replace = map_tet_replace.size();
397     char mess[255];
398     sprintf(mess," Nb tetra replace : %li (%lf%%)",nb_tetra_replace,((double)nb_tetra_replace/fem->get_nb_fem_element3())*(double)100.0);
399     affiche(mess);
400 francois 963 }
401 couturad 968 return OK;
402 francois 915 }
403 francois 883
404    
405 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)
406 francois 883 {
407 francois 915 int recal=(int)param.get_valeur("Quadratisation_jmin");
408 francois 1164 int dimsansgeo=0;
409     bool virtuel =false;
410     if (mg_geometrie!=NULL)
411     if (strcmp(mg_geometrie->get_type_geometrie(),"VIRTUEL")==0) virtuel=true;
412     if (virtuel)
413     {
414     if (maillage->get_nb_mg_tetra()+maillage->get_nb_mg_hexa()!=0) dimsansgeo=3;
415     else if (maillage->get_nb_mg_triangle()+maillage->get_nb_mg_quadrangle()!=0) dimsansgeo=2;
416     else if (maillage->get_nb_mg_segment()!=0) dimsansgeo=1;
417     }
418 francois 883 LISTE_MG_NOEUD::iterator it_nd;
419     MG_NOEUD* mgnoeud=maillage->get_premier_noeud(it_nd);
420 francois 1164 if (courbure_discrete) prepare_quad_discrete(fem->get_mg_maillage(),virtuel);
421 francois 883 std::vector<FEM_NOEUD*> lstnoeuddirect;
422     unsigned int i=0;
423     while (mgnoeud)
424     {
425     mgnoeud->change_nouveau_numero(i);
426     FEM_NOEUD* femnoeud=new FEM_NOEUD(mgnoeud);
427     femnoeud->change_numero(i+1);
428     i++ ;
429     fem->ajouter_fem_noeud(femnoeud);
430     //mgnoeud=maillage->get_suivant_noeud(it_nd);
431     lstnoeuddirect.insert(lstnoeuddirect.end(),femnoeud);
432 francois 1164 if (mgnoeud->get_lien_topologie()!=NULL)
433     if (mgnoeud->get_lien_topologie()->get_dimension()==0)
434     {
435     FEM_ELEMENT0* femele=new FEM_ELEMENT0(mgnoeud->get_lien_topologie(),mgnoeud,&femnoeud);
436     fem->ajouter_fem_element0(femele);
437     }
438     mgnoeud=maillage->get_suivant_noeud(it_nd);
439 francois 883
440     }
441     if (mg_geometrie==NULL)
442     {
443     if (maillage->get_nb_mg_tetra()!=0) dimsansgeo=3;
444     else if (maillage->get_nb_mg_triangle()!=0) dimsansgeo=2;
445     else if (maillage->get_nb_mg_segment()!=0) dimsansgeo=1;
446     }
447 francois 1164 int nbmgnoeud = maillage->get_nb_mg_noeud();
448    
449 francois 883 LISTE_MG_SEGMENT::iterator it_seg;
450     MG_SEGMENT* mgseg=maillage->get_premier_segment(it_seg);
451     FEM_MAILLAGE_QUADRATIQUE_OUTILS ou ;
452     i=0;
453     while (mgseg)
454     {
455     mgseg->change_nouveau_numero(i++);
456     MG_NOEUD *noeud1=mgseg->get_noeud1();
457     MG_NOEUD *noeud2=mgseg->get_noeud2();
458     double *xyz1=noeud1->get_coord();
459     double *xyz2=noeud2->get_coord();
460     double xori=0.5*(xyz1[0]+xyz2[0]);
461     double yori=0.5*(xyz1[1]+xyz2[1]);
462     double zori=0.5*(xyz1[2]+xyz2[2]);
463 ghazal 1152 FEM_NOEUD* femnoeud;
464 francois 883 if (mgseg->get_lien_topologie()!=NULL)
465     if (mgseg->get_lien_topologie()->get_dimension()==1)
466     {
467 francois 1164 double L3=0,L4=0,j=0.0;
468 ghazal 1152 double xyz[3]={xori,yori,zori};
469     if (!virtuel)
470     {
471 francois 1164 MG_ARETE* arete=(MG_ARETE*)(mgseg->get_lien_topologie());
472     double t1,t2,tm;
473     arete->inverser(t1,xyz1);
474     arete->inverser(t2,xyz2);
475     if (arete->get_courbe()->est_periodique())
476     if (t1>t2) t2=t2+arete->get_courbe()->get_periode();
477     double t=0.5*(t1+t2);
478     double xyz3[3];
479     arete->evaluer(t,xyz);
480     double L1= ou.get_distance_curviligne(0,0.5,xyz1,xyz,xyz2);
481     double L2= ou.get_distance_curviligne(0.5,1,xyz1,xyz,xyz2);
482     if (L1/L2 >1.01)
483     {
484     do
485     {
486     j+=0.0001;
487     tm=(-t+t1)*j+t;
488     arete->evaluer(tm,xyz3);
489     L3= ou.get_distance_curviligne(0,0.5,xyz1,xyz3,xyz2);
490     L4= ou.get_distance_curviligne(0.5,1,xyz1,xyz3,xyz2);
491    
492     }while (0.5-(L4/(L3+L4))>=0.0001);
493    
494     femnoeud= new FEM_NOEUD(mgseg->get_lien_topologie(),xyz3[0],xyz3[1],xyz3[2],xori,yori,zori);
495     }
496    
497     if (L2/L1 >1.01)
498     {
499     do
500     {
501     j+=0.0001;
502     tm=(-t+t2)*j+t;
503     arete->evaluer(tm,xyz3);
504     L3= ou.get_distance_curviligne(0,0.5,xyz1,xyz3,xyz2);
505     L4= ou.get_distance_curviligne(0.5,1,xyz1,xyz3,xyz2);
506    
507     }while (0.5-(L3/(L3+L4))>=0.0001);
508    
509     femnoeud= new FEM_NOEUD(mgseg->get_lien_topologie(),xyz3[0],xyz3[1],xyz3[2],xori,yori,zori);
510     }
511     }
512     if(j == 0) // geo virtuel
513     {
514     if (courbure_discrete) get_quad_seg(mgseg,virtuel,xyz);
515     femnoeud= new FEM_NOEUD(mgseg->get_lien_topologie(),xyz[0],xyz[1],xyz[2],xori,yori,zori);
516     }
517 francois 1095
518 francois 883 }
519     if (mgseg->get_lien_topologie()!=NULL)
520 francois 1164 if (mgseg->get_lien_topologie()->get_dimension()==2)
521     {
522     double L3=0,L4=0,j=0;
523     double xyz[3]={xori,yori,zori};
524     if (!virtuel)
525     {
526     MG_FACE* face=(MG_FACE*)(mgseg->get_lien_topologie());
527     double uv1[2],uv2[2];
528     OT_DECALAGE_PARAMETRE decalage(face->get_surface()->get_periode_u(),face->get_surface()->get_periode_v());
529     face->inverser(uv1,xyz1);
530     face->inverser(uv2,xyz2);
531     if (face->get_surface()->est_periodique_u()==1)
532     {
533     double eps=1e-10*face->get_surface()->get_periode_u();
534     if (uv1[0]<-eps) uv1[0]=uv1[0]+face->get_surface()->get_periode_u();
535     if (uv1[0]>face->get_surface()->get_periode_u()-eps) uv1[0]=uv1[0]-face->get_surface()->get_periode_u();
536     if (uv2[0]<-eps) uv2[0]=uv2[0]+face->get_surface()->get_periode_u();
537     if (uv2[0]>face->get_surface()->get_periode_u()-eps) uv2[0]=uv2[0]-face->get_surface()->get_periode_u();
538     }
539     if (face->get_surface()->est_periodique_v()==1)
540     {
541     double eps=1e-10*face->get_surface()->get_periode_v();
542     if (uv1[1]<-eps) uv1[1]=uv1[1]+face->get_surface()->get_periode_v();
543     if (uv1[1]>face->get_surface()->get_periode_v()-eps) uv1[1]=uv1[1]-face->get_surface()->get_periode_v();
544     if (uv2[1]<-eps) uv2[1]=uv2[1]+face->get_surface()->get_periode_v();
545     if (uv2[1]>face->get_surface()->get_periode_v()-eps) uv2[1]=uv2[1]-face->get_surface()->get_periode_v();
546     }
547     double du=decalage.calcul_decalage_parametre_u(uv1[0]);
548     double dv=decalage.calcul_decalage_parametre_v(uv1[1]);
549     double u1=decalage.decalage_parametre_u(uv1[0],du);
550     double v1=decalage.decalage_parametre_v(uv1[1],dv);
551     double u2=decalage.decalage_parametre_u(uv2[0],du);
552     double v2=decalage.decalage_parametre_v(uv2[1],dv);
553     double uv[2],uv3[2],xyz3[3];
554     double milieu[3];
555     milieu[0]=0.5*(xyz1[0]+xyz2[0]);
556     milieu[1]=0.5*(xyz1[1]+xyz2[1]);
557     milieu[2]=0.5*(xyz1[2]+xyz2[2]);
558     OT_VECTEUR_3D vec(xyz1,xyz2);
559     face->inverser(uv,milieu,0.25*vec.get_longueur());
560     face->evaluer(uv,xyz);
561     double u=decalage.decalage_parametre_u(uv[0],du);
562     double v=decalage.decalage_parametre_v(uv[1],dv);
563    
564     double L1= ou.get_distance_curviligne(0,0.5,xyz1,xyz,xyz2);
565     double L2= ou.get_distance_curviligne(0.5,1,xyz1,xyz,xyz2);
566     if (L1/L2 >1.01)
567     {
568     do
569     {
570     j+=0.0001;
571     if ((u1<u) && (u2<u) && (u1<u2)) u1=u1+face->get_surface()->get_periode_u();
572     if ((v1<v) && (v2<v) && (v1<v2)) v1=v1+face->get_surface()->get_periode_v();
573    
574     double u3=(-u+u1)*j+u;
575     double v3=(-v+v1)*j+v;
576     uv3[0]=decalage.decalage_parametre_u(u3,-du);
577     uv3[1]=decalage.decalage_parametre_v(v3,-dv);
578     face->evaluer(uv3,xyz3);
579     L3= ou.get_distance_curviligne(0,0.5,xyz1,xyz3,xyz2);
580     L4= ou.get_distance_curviligne(0.5,1,xyz1,xyz3,xyz2);
581    
582     }while ( 0.5-(L4/(L3+L4)) >=0.0001);
583    
584     femnoeud= new FEM_NOEUD(mgseg->get_lien_topologie(),xyz3[0],xyz3[1],xyz3[2],xori,yori,zori);
585     }
586    
587     if (L2/L1 >1.01)
588     {
589     do
590     {
591     j+=0.0001;
592     if ((u1<u) && (u2<u) && (u2<u1)) u2=u2+face->get_surface()->get_periode_u();
593     if ((v1<v) && (v2<v) && (v2<v1)) v2=v2+face->get_surface()->get_periode_v();
594    
595     double u3=(-u+u2)*j+u;
596     double v3=(-v+v2)*j+v;
597     uv3[0]=decalage.decalage_parametre_u(u3,-du);
598     uv3[1]=decalage.decalage_parametre_v(v3,-dv);
599    
600     face->evaluer(uv3,xyz3);
601     L3= ou.get_distance_curviligne(0,0.5,xyz1,xyz3,xyz2);
602     L4= ou.get_distance_curviligne(0.5,1,xyz1,xyz3,xyz2);
603    
604     }while (0.5-(L3/(L3+L4)) >=0.0001);
605    
606     femnoeud= new FEM_NOEUD(mgseg->get_lien_topologie(),xyz3[0],xyz3[1],xyz3[2],xori,yori,zori);
607     }
608     }
609     if (j == 0.0) // geo virtuel
610     {
611     if (courbure_discrete) get_quad_seg(mgseg,virtuel,xyz);
612     femnoeud= new FEM_NOEUD(mgseg->get_lien_topologie(),xyz[0],xyz[1],xyz[2],xori,yori,zori);
613     }
614    
615     }
616     if (mgseg->get_lien_topologie()!=NULL)
617     {
618     if (mgseg->get_lien_topologie()->get_dimension()==3)
619     {
620     double x=0.5*(xyz1[0]+xyz2[0]);
621     double y=0.5*(xyz1[1]+xyz2[1]);
622     double z=0.5*(xyz1[2]+xyz2[2]);
623     femnoeud= new FEM_NOEUD(mgseg->get_lien_topologie(),x,y,z,x,y,z);
624     }
625     }
626     else // pas de geo
627     {
628     double xyz[3];
629     xyz[0]=xori;
630     xyz[1]=yori;
631     xyz[2]=zori;
632     if (courbure_discrete) get_quad_seg(mgseg,virtuel,xyz);
633     femnoeud= new FEM_NOEUD(mgseg->get_lien_topologie(),xyz[0],xyz[1],xyz[2],xori,yori,zori);
634     }
635     femnoeud->change_numero(i+nbmgnoeud+1);
636     fem->ajouter_fem_noeud(femnoeud);
637     mgseg=maillage->get_suivant_segment(it_seg);
638     lstnoeuddirect.insert(lstnoeuddirect.end(),femnoeud);
639 francois 883 }
640    
641     if ((mg_geometrie!=NULL) || ( (mg_geometrie==NULL) && (dimsansgeo==1)))
642     {
643 francois 1164
644 francois 883 mgseg=maillage->get_premier_segment(it_seg);
645     while (mgseg)
646     {
647     if (mgseg->get_lien_topologie()!=NULL)
648 francois 1164 if (mgseg->get_lien_topologie()->get_dimension()!=1)
649     {
650 francois 883 mgseg=maillage->get_suivant_segment(it_seg);
651     continue;
652 francois 1164 }
653     if (mg_geometrie!=NULL)
654     if (mgseg->get_lien_topologie()==NULL)
655     if (dimsansgeo!=1)
656     {
657     mgseg=maillage->get_suivant_segment(it_seg);
658     continue;
659     }
660     FEM_NOEUD *tabnoeud[3];
661     tabnoeud[0]=lstnoeuddirect[mgseg->get_noeud1()->get_nouveau_numero()];
662     tabnoeud[1]=lstnoeuddirect[mgseg->get_nouveau_numero()+nbmgnoeud];
663     tabnoeud[2]=lstnoeuddirect[mgseg->get_noeud2()->get_nouveau_numero()];
664     FEM_SEGMENT3* seg=new FEM_SEGMENT3(mgseg,tabnoeud);
665     fem->ajouter_fem_element1(seg);
666 francois 883 mgseg=maillage->get_suivant_segment(it_seg);
667     }
668     }
669     if ((mg_geometrie!=NULL) || ( (mg_geometrie==NULL) && (dimsansgeo==2)))
670     {
671     LISTE_MG_TRIANGLE::iterator it_tri;
672     MG_TRIANGLE* mgtri=maillage->get_premier_triangle(it_tri);
673     while (mgtri)
674     {
675     if (mgtri->get_lien_topologie()!=NULL) if (mgtri->get_lien_topologie()->get_dimension()!=2)
676 francois 1164 {
677     mgtri=maillage->get_suivant_triangle(it_tri);
678     continue;
679     }
680     if (mg_geometrie!=NULL)
681     if (mgtri->get_lien_topologie()==NULL)
682     if (dimsansgeo!=2)
683     {
684     mgtri=maillage->get_suivant_triangle(it_tri);
685     continue;
686 francois 883 }
687 francois 1164 FEM_NOEUD *tabnoeud[6];
688     tabnoeud[0]=lstnoeuddirect[mgtri->get_noeud1()->get_nouveau_numero()];
689     tabnoeud[1]=lstnoeuddirect[mgtri->get_segment1()->get_nouveau_numero()+nbmgnoeud];
690     tabnoeud[2]=lstnoeuddirect[mgtri->get_noeud2()->get_nouveau_numero()];
691     tabnoeud[3]=lstnoeuddirect[mgtri->get_segment2()->get_nouveau_numero()+nbmgnoeud];
692     tabnoeud[4]=lstnoeuddirect[mgtri->get_noeud3()->get_nouveau_numero()];
693     tabnoeud[5]=lstnoeuddirect[mgtri->get_segment3()->get_nouveau_numero()+nbmgnoeud];
694     FEM_TRIANGLE6* tri=new FEM_TRIANGLE6(mgtri,tabnoeud);
695     fem->ajouter_fem_element2(tri);
696 francois 883 mgtri=maillage->get_suivant_triangle(it_tri);
697     }
698     LISTE_MG_QUADRANGLE::iterator it_quad;
699     MG_QUADRANGLE* mgquad=maillage->get_premier_quadrangle(it_quad);
700     while (mgquad)
701     {
702     if (mgquad->get_lien_topologie()!=NULL) if (mgquad->get_lien_topologie()->get_dimension()!=2) {
703 francois 1164 mgquad=maillage->get_suivant_quadrangle(it_quad);
704     continue;
705     }
706 francois 883 if (mg_geometrie!=NULL)
707 francois 1164 if (mgquad->get_lien_topologie()==NULL)
708     if (dimsansgeo!=2)
709     {
710     mgquad=maillage->get_suivant_quadrangle(it_quad);
711     continue;
712     }
713     FEM_NOEUD *tabnoeud[8];
714 francois 883 tabnoeud[0]=lstnoeuddirect[mgquad->get_noeud1()->get_nouveau_numero()];
715 francois 1164 tabnoeud[1]=lstnoeuddirect[mgquad->get_segment1()->get_nouveau_numero()+nbmgnoeud];
716     tabnoeud[2]=lstnoeuddirect[mgquad->get_noeud2()->get_nouveau_numero()];
717     tabnoeud[3]=lstnoeuddirect[mgquad->get_segment2()->get_nouveau_numero()+nbmgnoeud];
718     tabnoeud[4]=lstnoeuddirect[mgquad->get_noeud3()->get_nouveau_numero()];
719     tabnoeud[5]=lstnoeuddirect[mgquad->get_segment3()->get_nouveau_numero()+nbmgnoeud];
720     tabnoeud[6]=lstnoeuddirect[mgquad->get_noeud4()->get_nouveau_numero()];
721     tabnoeud[7]=lstnoeuddirect[mgquad->get_segment4()->get_nouveau_numero()+nbmgnoeud];
722     FEM_QUADRANGLE8* quad=new FEM_QUADRANGLE8(mgquad,tabnoeud);
723     fem->ajouter_fem_element2(quad);
724 francois 883 mgquad=maillage->get_suivant_quadrangle(it_quad);
725     }
726     }
727     if ((mg_geometrie!=NULL) || ( (mg_geometrie==NULL) && (dimsansgeo==3)))
728     {
729 francois 1164
730     LISTE_MG_TETRA::iterator it_tetra;
731 francois 883 MG_TETRA* mgtetra=maillage->get_premier_tetra(it_tetra);
732     while (mgtetra)
733     {
734     FEM_NOEUD *tabnoeud[10];
735     tabnoeud[0]=lstnoeuddirect[mgtetra->get_noeud1()->get_nouveau_numero()];
736     tabnoeud[1]=lstnoeuddirect[maillage->get_mg_segment(mgtetra->get_noeud1()->get_id(),mgtetra->get_noeud2()->get_id())->get_nouveau_numero()+nbmgnoeud];
737     tabnoeud[2]=lstnoeuddirect[mgtetra->get_noeud2()->get_nouveau_numero()];
738     tabnoeud[3]=lstnoeuddirect[maillage->get_mg_segment(mgtetra->get_noeud2()->get_id(),mgtetra->get_noeud3()->get_id())->get_nouveau_numero()+nbmgnoeud];
739     tabnoeud[4]=lstnoeuddirect[mgtetra->get_noeud3()->get_nouveau_numero()];
740     tabnoeud[5]=lstnoeuddirect[maillage->get_mg_segment(mgtetra->get_noeud3()->get_id(),mgtetra->get_noeud1()->get_id())->get_nouveau_numero()+nbmgnoeud];
741     tabnoeud[6]=lstnoeuddirect[maillage->get_mg_segment(mgtetra->get_noeud1()->get_id(),mgtetra->get_noeud4()->get_id())->get_nouveau_numero()+nbmgnoeud];
742     tabnoeud[7]=lstnoeuddirect[maillage->get_mg_segment(mgtetra->get_noeud2()->get_id(),mgtetra->get_noeud4()->get_id())->get_nouveau_numero()+nbmgnoeud];
743     tabnoeud[8]=lstnoeuddirect[maillage->get_mg_segment(mgtetra->get_noeud3()->get_id(),mgtetra->get_noeud4()->get_id())->get_nouveau_numero()+nbmgnoeud];
744     tabnoeud[9]=lstnoeuddirect[mgtetra->get_noeud4()->get_nouveau_numero()];
745     FEM_TETRA10* tet=new FEM_TETRA10(mgtetra,tabnoeud);
746     fem->ajouter_fem_element3(tet);
747 francois 1164 mgtetra=maillage->get_suivant_tetra(it_tetra);
748     }
749     if (recal==1) if(recal_element_quadratique(fem)==FAIL) return FAIL;
750 francois 915
751 francois 1164
752    
753 francois 883 LISTE_MG_HEXA::iterator it_hexa;
754     MG_HEXA* mghex=maillage->get_premier_hexa(it_hexa);
755     while (mghex)
756     {
757     FEM_NOEUD *tabnoeud[20];
758     tabnoeud[0]=lstnoeuddirect[mghex->get_noeud1()->get_nouveau_numero()];
759     tabnoeud[1]=lstnoeuddirect[maillage->get_mg_segment(mghex->get_noeud1()->get_id(),mghex->get_noeud2()->get_id())->get_nouveau_numero()+nbmgnoeud];
760     tabnoeud[2]=lstnoeuddirect[mghex->get_noeud2()->get_nouveau_numero()];
761     tabnoeud[3]=lstnoeuddirect[maillage->get_mg_segment(mghex->get_noeud2()->get_id(),mghex->get_noeud3()->get_id())->get_nouveau_numero()+nbmgnoeud];
762     tabnoeud[4]=lstnoeuddirect[mghex->get_noeud3()->get_nouveau_numero()];
763     tabnoeud[5]=lstnoeuddirect[maillage->get_mg_segment(mghex->get_noeud3()->get_id(),mghex->get_noeud4()->get_id())->get_nouveau_numero()+nbmgnoeud];
764     tabnoeud[6]=lstnoeuddirect[mghex->get_noeud4()->get_nouveau_numero()];
765     tabnoeud[7]=lstnoeuddirect[maillage->get_mg_segment(mghex->get_noeud4()->get_id(),mghex->get_noeud1()->get_id())->get_nouveau_numero()+nbmgnoeud];
766     tabnoeud[8]=lstnoeuddirect[maillage->get_mg_segment(mghex->get_noeud1()->get_id(),mghex->get_noeud5()->get_id())->get_nouveau_numero()+nbmgnoeud];
767     tabnoeud[9]=lstnoeuddirect[maillage->get_mg_segment(mghex->get_noeud2()->get_id(),mghex->get_noeud6()->get_id())->get_nouveau_numero()+nbmgnoeud];
768     tabnoeud[10]=lstnoeuddirect[maillage->get_mg_segment(mghex->get_noeud3()->get_id(),mghex->get_noeud7()->get_id())->get_nouveau_numero()+nbmgnoeud];
769     tabnoeud[11]=lstnoeuddirect[maillage->get_mg_segment(mghex->get_noeud4()->get_id(),mghex->get_noeud8()->get_id())->get_nouveau_numero()+nbmgnoeud];
770     tabnoeud[12]=lstnoeuddirect[mghex->get_noeud5()->get_nouveau_numero()];
771     tabnoeud[13]=lstnoeuddirect[maillage->get_mg_segment(mghex->get_noeud5()->get_id(),mghex->get_noeud6()->get_id())->get_nouveau_numero()+nbmgnoeud];
772     tabnoeud[14]=lstnoeuddirect[mghex->get_noeud6()->get_nouveau_numero()];
773     tabnoeud[15]=lstnoeuddirect[maillage->get_mg_segment(mghex->get_noeud6()->get_id(),mghex->get_noeud7()->get_id())->get_nouveau_numero()+nbmgnoeud];
774     tabnoeud[16]=lstnoeuddirect[mghex->get_noeud7()->get_nouveau_numero()];
775     tabnoeud[17]=lstnoeuddirect[maillage->get_mg_segment(mghex->get_noeud7()->get_id(),mghex->get_noeud8()->get_id())->get_nouveau_numero()+nbmgnoeud];
776     tabnoeud[18]=lstnoeuddirect[mghex->get_noeud8()->get_nouveau_numero()];
777     tabnoeud[19]=lstnoeuddirect[maillage->get_mg_segment(mghex->get_noeud8()->get_id(),mghex->get_noeud5()->get_id())->get_nouveau_numero()+nbmgnoeud];
778     FEM_HEXA20* hex=new FEM_HEXA20(mghex,tabnoeud);
779     fem->ajouter_fem_element3(hex);
780     mghex=maillage->get_suivant_hexa(it_hexa);
781     }
782     LISTE_MG_PENTA::iterator it_pen;
783     MG_PENTA* mgpen=maillage->get_premier_penta(it_pen);
784     while (mgpen)
785     {
786     FEM_NOEUD *tabnoeud[15];
787     tabnoeud[0]=lstnoeuddirect[mgpen->get_noeud1()->get_nouveau_numero()];
788     tabnoeud[1]=lstnoeuddirect[maillage->get_mg_segment(mgpen->get_noeud1()->get_id(),mgpen->get_noeud2()->get_id())->get_nouveau_numero()+nbmgnoeud];
789     tabnoeud[2]=lstnoeuddirect[mgpen->get_noeud2()->get_nouveau_numero()];
790     tabnoeud[3]=lstnoeuddirect[maillage->get_mg_segment(mgpen->get_noeud2()->get_id(),mgpen->get_noeud3()->get_id())->get_nouveau_numero()+nbmgnoeud];
791     tabnoeud[4]=lstnoeuddirect[mgpen->get_noeud3()->get_nouveau_numero()];
792     tabnoeud[5]=lstnoeuddirect[maillage->get_mg_segment(mgpen->get_noeud3()->get_id(),mgpen->get_noeud1()->get_id())->get_nouveau_numero()+nbmgnoeud];
793     tabnoeud[6]=lstnoeuddirect[maillage->get_mg_segment(mgpen->get_noeud1()->get_id(),mgpen->get_noeud4()->get_id())->get_nouveau_numero()+nbmgnoeud];
794     tabnoeud[7]=lstnoeuddirect[maillage->get_mg_segment(mgpen->get_noeud2()->get_id(),mgpen->get_noeud5()->get_id())->get_nouveau_numero()+nbmgnoeud];
795     tabnoeud[8]=lstnoeuddirect[maillage->get_mg_segment(mgpen->get_noeud3()->get_id(),mgpen->get_noeud6()->get_id())->get_nouveau_numero()+nbmgnoeud];
796     tabnoeud[9]=lstnoeuddirect[mgpen->get_noeud4()->get_nouveau_numero()];
797     tabnoeud[10]=lstnoeuddirect[maillage->get_mg_segment(mgpen->get_noeud4()->get_id(),mgpen->get_noeud5()->get_id())->get_nouveau_numero()+nbmgnoeud];
798     tabnoeud[11]=lstnoeuddirect[mgpen->get_noeud5()->get_nouveau_numero()];
799     tabnoeud[12]=lstnoeuddirect[maillage->get_mg_segment(mgpen->get_noeud5()->get_id(),mgpen->get_noeud6()->get_id())->get_nouveau_numero()+nbmgnoeud];
800     tabnoeud[13]=lstnoeuddirect[mgpen->get_noeud6()->get_nouveau_numero()];
801     tabnoeud[14]=lstnoeuddirect[maillage->get_mg_segment(mgpen->get_noeud4()->get_id(),mgpen->get_noeud6()->get_id())->get_nouveau_numero()+nbmgnoeud];
802     FEM_PENTA15* pen=new FEM_PENTA15(mgpen,tabnoeud);
803     fem->ajouter_fem_element3(pen);
804     mgpen=maillage->get_suivant_penta(it_pen);
805     }
806     }
807 couturad 968 return OK;
808 francois 883 }
809    
810    
811    
812    
813     void MAILLEUR_FEM::optimise_numerotation(FEM_MAILLAGE* fem)
814     {
815     FEM_NOEUD* noeud=fem->get_fem_noeud(0);
816     FEM_GRAPH_NOEUD *graph;
817     graph=new FEM_GRAPH_NOEUD(noeud,fem);
818     int e=graph->get_excentricite();
819     TPL_MAP_ENTITE<FEM_NOEUD*> dernierniveau=graph->get_dernier_niveau();
820     int nbdernierniveau=dernierniveau.get_nb();
821     for (int i=0;i<nbdernierniveau;i++)
822     {
823     FEM_GRAPH_NOEUD* graphtmp=new FEM_GRAPH_NOEUD(dernierniveau.get(i),fem);
824     int etmp=graphtmp->get_excentricite();
825     if (etmp>e)
826     {
827     delete graph;
828     graph=graphtmp;
829     dernierniveau=graph->get_dernier_niveau();
830     nbdernierniveau=dernierniveau.get_nb();
831     i=-1;
832     e=etmp;
833     }
834     else delete graphtmp;
835     }
836     int numero=fem->get_nb_fem_noeud();
837     for (int i=0;i<e+1;i++)
838     {
839     TPL_MAP_ENTITE<FEM_NOEUD*> niveau=graph->get_niveau(i);
840     int nbnoeud=niveau.get_nb();
841     for (int j=0;j<nbnoeud;j++)
842     {
843     FEM_NOEUD* noeud=niveau.get(j);
844     noeud->change_numero_opt(numero);
845     numero--;
846     }
847     }
848    
849    
850     delete graph;
851     }
852    
853    
854     void MAILLEUR_FEM::copie_numerotation_opt(FEM_MAILLAGE* fem)
855     {
856     int nbnoeud=fem->get_nb_fem_noeud();
857     for (int i=0;i<nbnoeud;i++)
858     {
859     FEM_NOEUD* noeud=fem->get_fem_noeud(i);
860     noeud->change_numero(noeud->get_numero_opt());
861     }
862    
863     }
864    
865 francois 1164 void MAILLEUR_FEM::prepare_quad_discrete(MG_MAILLAGE* mai,bool virtuel)
866     {
867     MG_GEOMETRIE *geo=mai->get_mg_geometrie();
868     LISTE_MG_TETRA::iterator it;
869     LISTE_MG_NOEUD::iterator itn;
870     LISTE_MG_TRIANGLE::iterator itt;
871     LISTE_MG_SEGMENT::iterator its;
872     for (MG_TETRA* tet=mai->get_premier_tetra(it);tet!=NULL;tet=mai->get_suivant_tetra(it))
873     {
874     tet->change_nouveau_numero(0);
875     tet->get_triangle1()->change_nouveau_numero(0);
876     tet->get_triangle2()->change_nouveau_numero(0);
877     tet->get_triangle3()->change_nouveau_numero(0);
878     tet->get_triangle4()->change_nouveau_numero(0);
879     }
880     for (MG_TETRA* tet=mai->get_premier_tetra(it);tet!=NULL;tet=mai->get_suivant_tetra(it))
881     {
882     if (tet->get_lien_topologie()!=NULL)
883     {
884     if (tet->get_triangle1()->get_nouveau_numero()==0)
885     tet->get_triangle1()->change_nouveau_numero(tet->get_lien_topologie()->get_id());
886     else if (tet->get_triangle1()->get_nouveau_numero()>2)
887     {
888     if (tet->get_triangle1()->get_nouveau_numero()==tet->get_lien_topologie()->get_id())
889     tet->get_triangle1()->change_nouveau_numero(2);
890     else
891     tet->get_triangle1()->change_nouveau_numero(1);
892     }
893    
894     }
895     else tet->get_triangle1()->change_nouveau_numero(tet->get_triangle1()->get_nouveau_numero()+1);
896    
897     if (tet->get_lien_topologie()!=NULL)
898     {
899     if (tet->get_triangle2()->get_nouveau_numero()==0)
900     tet->get_triangle2()->change_nouveau_numero(tet->get_lien_topologie()->get_id());
901     else if (tet->get_triangle2()->get_nouveau_numero()>2)
902     {
903     if (tet->get_triangle2()->get_nouveau_numero()==tet->get_lien_topologie()->get_id())
904     tet->get_triangle2()->change_nouveau_numero(2);
905     else
906     tet->get_triangle2()->change_nouveau_numero(1);
907     }
908    
909     }
910     else tet->get_triangle2()->change_nouveau_numero(tet->get_triangle2()->get_nouveau_numero()+1);
911    
912     if (tet->get_lien_topologie()!=NULL)
913     {
914     if (tet->get_triangle3()->get_nouveau_numero()==0)
915     tet->get_triangle3()->change_nouveau_numero(tet->get_lien_topologie()->get_id());
916     else if (tet->get_triangle3()->get_nouveau_numero()>2)
917     {
918     if (tet->get_triangle3()->get_nouveau_numero()==tet->get_lien_topologie()->get_id())
919     tet->get_triangle3()->change_nouveau_numero(2);
920     else
921     tet->get_triangle3()->change_nouveau_numero(1);
922     }
923    
924     }
925     else tet->get_triangle3()->change_nouveau_numero(tet->get_triangle3()->get_nouveau_numero()+1);
926    
927     if (tet->get_lien_topologie()!=NULL)
928     {
929     if (tet->get_triangle4()->get_nouveau_numero()==0)
930     tet->get_triangle4()->change_nouveau_numero(tet->get_lien_topologie()->get_id());
931     else if (tet->get_triangle4()->get_nouveau_numero()>2)
932     {
933     if (tet->get_triangle4()->get_nouveau_numero()==tet->get_lien_topologie()->get_id())
934     tet->get_triangle4()->change_nouveau_numero(2);
935     else
936     tet->get_triangle4()->change_nouveau_numero(1);
937     }
938    
939     }
940     else tet->get_triangle4()->change_nouveau_numero(tet->get_triangle4()->get_nouveau_numero()+1);
941    
942     }
943     if (geo!=NULL)
944     {
945     LISTE_MG_FACE::iterator itf;
946     for (MG_FACE* face=geo->get_premier_face(itf);face!=NULL;face=geo->get_suivant_face(itf))
947     {
948     int nbele=face->get_lien_maillage()->get_nb();
949     for (int i=0;i<nbele;i++)
950     {
951     MG_ELEMENT_MAILLAGE* ele=face->get_lien_maillage()->get(i);
952     ele->change_nouveau_numero(1);
953     }
954    
955     }
956     }
957     for (MG_SEGMENT* seg=mai->get_premier_segment(its);seg!=NULL;seg=mai->get_suivant_segment(its))
958     seg->change_nouveau_numero(0);
959     for (MG_TRIANGLE* tri=mai->get_premier_triangle(itt);tri!=NULL;tri=mai->get_suivant_triangle(itt))
960     {
961     if (tri->get_nouveau_numero()>2) tri->change_nouveau_numero(1);
962     if (tri->get_lien_topologie()!=NULL)
963     if (tri->get_lien_topologie())
964     if (virtuel)
965     if ((((MG_FACE_ELEMENT*)tri->get_lien_topologie())->get_nb_contrainte()>0))
966     {
967     segnonprojete[tri->get_segment1()]=tri->get_segment1();
968     segnonprojete[tri->get_segment2()]=tri->get_segment2();
969     segnonprojete[tri->get_segment3()]=tri->get_segment3();
970     }
971     if (tri->get_nouveau_numero()==1)
972     {
973     segaprojete[tri->get_segment1()]=tri->get_segment1();
974     segaprojete[tri->get_segment2()]=tri->get_segment2();
975     segaprojete[tri->get_segment3()]=tri->get_segment3();
976     }
977     }
978     }
979    
980     bool MAILLEUR_FEM::get_courbure_noeud(MG_NOEUD* no,MG_NOEUD* dirno,OT_VECTEUR_3D &n,double &r)
981     {
982     double deno=0.;
983     OT_VECTEUR_3D PP1(no->get_coord(),dirno->get_coord());
984     OT_VECTEUR_3D n0(0.,0.,0.);
985     n=n0;
986     int passe=0;
987     for (int i=0;i<no->get_lien_triangle()->get_nb();i++)
988     {
989     MG_TRIANGLE *tri=no->get_lien_triangle()->get(i);
990     if (tri->get_nouveau_numero()!=1) continue;
991     MG_NOEUD* not1=tri->get_noeud1();
992     MG_NOEUD* not2=tri->get_noeud2();
993     MG_NOEUD* not3=tri->get_noeud3();
994     OT_VECTEUR_3D vec1(not1->get_coord(),not2->get_coord());
995     OT_VECTEUR_3D vec2(not1->get_coord(),not3->get_coord());
996     OT_VECTEUR_3D pvec=vec1&vec2;
997     //if (PP1*pvec>0) pvec=-pvec;
998     double aire=pvec.get_longueur();
999     pvec.norme();
1000     n=n+aire*pvec;
1001     deno=deno+aire;
1002     passe++;
1003     }
1004     if (passe==0) return false;
1005     n=(1./deno)*n;
1006     n.norme();
1007     double signe=1.;
1008     if (PP1*n>0) signe=-1.;
1009     r=PP1*PP1/2./signe/(n*PP1);
1010     return true;
1011     }
1012    
1013     void MAILLEUR_FEM::get_quad_seg(MG_SEGMENT *seg,bool virtuel,double *xyz)
1014     {
1015     if (virtuel)
1016     if (seg->get_lien_topologie()!=NULL)
1017     if (((MG_ARETE_ELEMENT*)seg->get_lien_topologie())->get_nb_contrainte()>0) return;
1018     if (segaprojete[seg]==NULL)
1019     return;
1020     if (segnonprojete[seg]!=NULL)
1021     return;
1022     MG_NOEUD* no1=seg->get_noeud1();
1023     MG_NOEUD* no2=seg->get_noeud2();
1024     OT_VECTEUR_3D n1,n2;
1025     double r1,r2;
1026     bool def1=get_courbure_noeud(no1,no2,n1,r1);
1027     bool def2=get_courbure_noeud(no2,no1,n2,r2);
1028     if (def1==false) return;
1029     if (def2==false) return;
1030     double r=0.5*(r1+r2);
1031     OT_VECTEUR_3D n=0.5*(n1+n2);
1032 francois 1168 n.norme();
1033 francois 1164 double dist=0.;
1034     if (fabs(r)<1000.*seg->get_longueur())
1035     {
1036     dist=r*r-0.25*seg->get_longueur()*seg->get_longueur();
1037     dist=fabs(r)-fsqrt(dist);
1038     if (r<0) dist=-dist;
1039     }
1040     xyz[0]=xyz[0]+dist*n.get_x();
1041     xyz[1]=xyz[1]+dist*n.get_y();
1042     xyz[2]=xyz[2]+dist*n.get_z();
1043     }