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 (13 months 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

# Content
1 //####//------------------------------------------------------------
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 //####// Derniere modification par francois
20 //####// jeu 25 jui 2024 18:16:33 EDT
21 //####//------------------------------------------------------------
22 //####//------------------------------------------------------------
23
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 #include "mg_arete_element.h"
45 #include "mg_face_element.h"
46 #include <math.h>
47
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 int MAILLEUR_FEM::maille(FEM_MAILLAGE* fem,bool courbure_discrete,int num)
67 {
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 if (mg_geometrie->get_nb_mg_volume()!=0)
75 recherche_connexion_multidimension(mini,maillage,mg_geometrie);
76 if (degre==1) construire_lineaire(fem,mini,maillage,mg_geometrie);
77 if (degre==2) construire_quadratique(fem,mini,maillage,mg_geometrie,courbure_discrete);
78 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 return OK;
86 }
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 int MAILLEUR_FEM::construire_lineaire(FEM_MAILLAGE* fem,TPL_MAP_ENTITE<MG_SEGMENT*> &mini,MG_MAILLAGE* maillage,MG_GEOMETRIE* mg_geometrie)
152 {
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 return OK;
325 }
326
327
328 int MAILLEUR_FEM::recal_element_quadratique(FEM_MAILLAGE* fem)
329 {
330 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 for (int i=1;i<9;i++)
352 {
353 if ((i==2) || (i==4)) continue;
354 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 for (int i=1;i<9;i++)
367 {
368 if ((i==2) || (i==4)) continue;
369 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 /* if(pas>nbpas)
378 {
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 }*/
383 }
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 if (fem->get_nb_fem_element3()!=0)
395 {
396 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 }
401 return OK;
402 }
403
404
405 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 {
407 int recal=(int)param.get_valeur("Quadratisation_jmin");
408 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 LISTE_MG_NOEUD::iterator it_nd;
419 MG_NOEUD* mgnoeud=maillage->get_premier_noeud(it_nd);
420 if (courbure_discrete) prepare_quad_discrete(fem->get_mg_maillage(),virtuel);
421 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 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
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 int nbmgnoeud = maillage->get_nb_mg_noeud();
448
449 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 FEM_NOEUD* femnoeud;
464 if (mgseg->get_lien_topologie()!=NULL)
465 if (mgseg->get_lien_topologie()->get_dimension()==1)
466 {
467 double L3=0,L4=0,j=0.0;
468 double xyz[3]={xori,yori,zori};
469 if (!virtuel)
470 {
471 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
518 }
519 if (mgseg->get_lien_topologie()!=NULL)
520 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 }
640
641 if ((mg_geometrie!=NULL) || ( (mg_geometrie==NULL) && (dimsansgeo==1)))
642 {
643
644 mgseg=maillage->get_premier_segment(it_seg);
645 while (mgseg)
646 {
647 if (mgseg->get_lien_topologie()!=NULL)
648 if (mgseg->get_lien_topologie()->get_dimension()!=1)
649 {
650 mgseg=maillage->get_suivant_segment(it_seg);
651 continue;
652 }
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 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 {
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 }
687 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 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 mgquad=maillage->get_suivant_quadrangle(it_quad);
704 continue;
705 }
706 if (mg_geometrie!=NULL)
707 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 tabnoeud[0]=lstnoeuddirect[mgquad->get_noeud1()->get_nouveau_numero()];
715 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 mgquad=maillage->get_suivant_quadrangle(it_quad);
725 }
726 }
727 if ((mg_geometrie!=NULL) || ( (mg_geometrie==NULL) && (dimsansgeo==3)))
728 {
729
730 LISTE_MG_TETRA::iterator it_tetra;
731 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 mgtetra=maillage->get_suivant_tetra(it_tetra);
748 }
749 if (recal==1) if(recal_element_quadratique(fem)==FAIL) return FAIL;
750
751
752
753 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 return OK;
808 }
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 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 n.norme();
1033 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 }