ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/mailleur/src/mailleur1d.cpp
Revision: 61
Committed: Fri Nov 16 16:39:29 2007 UTC (17 years, 6 months ago) by francois
Original Path: magic/lib/mailleur/mailleur/src/mailleur1d.cpp
File size: 12830 byte(s)
Log Message:
Bub mailleur 3D + visualisation du front optimisée

File Contents

# User Rev Content
1 5 //------------------------------------------------------------
2     //------------------------------------------------------------
3     // MAGiC
4     // Jean Christophe Cuilli�e et Vincent FRANCOIS
5     // D�artement de G�ie M�anique - UQTR
6     //------------------------------------------------------------
7     // Le projet MAGIC est un projet de recherche du d�artement
8     // de g�ie m�anique de l'Universit�du Qu�ec �
9     // Trois Rivi�es
10     // Les librairies ne peuvent �re utilis�s sans l'accord
11     // des auteurs (contact : francois@uqtr.ca)
12     //------------------------------------------------------------
13     //------------------------------------------------------------
14     //
15     // mailleur1d.cpp
16     //
17     //------------------------------------------------------------
18     //------------------------------------------------------------
19     // COPYRIGHT 2000
20     // Version du 02/03/2006 �11H23
21     //------------------------------------------------------------
22     //------------------------------------------------------------
23    
24    
25     #include "gestionversion.h"
26     #include <math.h>
27     #include <vector>
28     #include "mailleur1d.h"
29 francois 35 #include "tpl_set.h"
30 5 //#include "affiche.h"
31     //#include "message.h"
32     #include <fstream>
33     struct integrale
34     {
35     double ti;
36     double li;
37     };
38    
39    
40     MAILLEUR1D::MAILLEUR1D(MG_MAILLAGE* mgmai,MG_GEOMETRIE *mggeo,MG_ARETE* mgarete,FCT_TAILLE* fct_taille):MAILLEUR(),mg_maillage(mgmai),mg_geometrie(mggeo),mg_arete(mgarete),metrique(fct_taille)
41     {
42     }
43    
44    
45    
46     MAILLEUR1D::~MAILLEUR1D()
47     {
48     }
49    
50    
51    
52    
53     void MAILLEUR1D::maille(void)
54     {
55     //afficheur << MAILLAGEARETE << endaff;
56     if (mg_arete!=NULL) maille(mg_arete);
57     else
58     {
59     int nb_arete=mg_geometrie->get_nb_mg_arete();
60     for (int i=0;i<nb_arete;i++)
61     {
62     MG_ARETE* mgarete=mg_geometrie->get_mg_arete(i);
63     maille(mgarete);
64     }
65     }
66     }
67    
68    
69     void MAILLEUR1D::maille(MG_ARETE* mgarete,double t1,MG_NOEUD* noeud_depart,double t2,MG_NOEUD* noeud_arrivee)
70     {
71     int creation_metrique=0;
72     if (metrique==NULL)
73     {
74     metrique=new FCT_TAILLE_ARETE(epsilon,distance_maximale,mgarete);
75     creation_metrique=1;
76     }
77    
78     /* calcul de la longueur de la arete dans la metrique */
79     if (fabs(t1)<1e-14)
80     if (fabs(t2)<1e-14)
81     {
82     t1=mgarete->get_tmin();
83     t2=mgarete->get_tmax();
84     }
85    
86     std::vector<integrale> tab_integrale;
87     double ti,tii,t;
88     double longueur_metrique=0.;
89     tii=t1;
90     double dt=0.001;
91     do
92     {
93     refresh();
94     double coefficient_metrique_ti[9];
95     double coefficient_metrique_tii[9];
96     double coefficient_metrique_derive_ti[9];
97     ti=tii;
98     tii=ti+(t2-t1)*dt;
99     if (tii>t2) tii=t2;
100    
101     #define ON_A_ABANDONNE_ET_ON_PREND_DT_EGALE_T2_MOINS_T1_DIVISE_PAR_1000
102     #ifndef ON_A_ABANDONNE_ET_ON_PREND_DT_EGALE_T2_MOINS_T1_DIVISE_PAR_1000
103     int bon_dt=0;
104     do
105     {
106     tii=ti+(t2-t1)*dt;
107     if (tii>t2) tii=t2;
108     double ddxyz[3];
109     double dxyz[3];
110     double xyz[3];
111     double dxyzii[3];
112     double xyzii[3];
113     mgarete->deriver_seconde(ti,ddxyz,dxyz,xyz);
114     mgarete->evaluer(tii,xyzii);
115     mgarete->deriver(tii,dxyzii);
116     if (creation_metrique)
117     {
118     metrique->evaluer(&ti,coefficient_metrique_ti);
119     metrique->evaluer(&tii,coefficient_metrique_tii);
120     metrique->deriver(&ti,coefficient_metrique_derive_ti);
121     }
122     else
123     {
124     metrique->evaluer(xyz,coefficient_metrique_ti);
125     metrique->evaluer(xyzii,coefficient_metrique_tii);
126     metrique->deriver(xyz,coefficient_metrique_derive_ti);
127     }
128     double facteur1=dxyz[0]*dxyz[0]*coefficient_metrique_ti[0]+dxyz[0]*dxyz[1]*coefficient_metrique_ti[3]+dxyz[0]*dxyz[2]*coefficient_metrique_ti[6]+dxyz[1]*dxyz[0]*coefficient_metrique_ti[1]+dxyz[1]*dxyz[1]*coefficient_metrique_ti[4]+dxyz[1]*dxyz[2]*coefficient_metrique_ti[7]+dxyz[2]*dxyz[0]*coefficient_metrique_ti[2]+dxyz[2]*dxyz[1]*coefficient_metrique_ti[4]+dxyz[2]*dxyz[2]*coefficient_metrique_ti[8];
129     double facteur2=ddxyz[0]*dxyz[0]*coefficient_metrique_ti[0]+dxyz[0]*ddxyz[0]*coefficient_metrique_ti[0]+dxyz[0]*dxyz[0]*coefficient_metrique_derive_ti[0]+
130     ddxyz[0]*dxyz[1]*coefficient_metrique_ti[3]+dxyz[0]*ddxyz[1]*coefficient_metrique_ti[3]+dxyz[0]*dxyz[1]*coefficient_metrique_derive_ti[3]+
131     ddxyz[0]*dxyz[2]*coefficient_metrique_ti[6]+dxyz[0]*ddxyz[2]*coefficient_metrique_ti[6]+dxyz[0]*dxyz[2]*coefficient_metrique_derive_ti[6]+
132     ddxyz[1]*dxyz[0]*coefficient_metrique_ti[1]+dxyz[1]*ddxyz[0]*coefficient_metrique_ti[1]+dxyz[1]*dxyz[0]*coefficient_metrique_derive_ti[1]+
133     ddxyz[1]*dxyz[1]*coefficient_metrique_ti[4]+dxyz[1]*ddxyz[1]*coefficient_metrique_ti[4]+dxyz[1]*dxyz[1]*coefficient_metrique_derive_ti[4]+
134     ddxyz[1]*dxyz[2]*coefficient_metrique_ti[7]+dxyz[1]*ddxyz[2]*coefficient_metrique_ti[7]+dxyz[1]*dxyz[2]*coefficient_metrique_derive_ti[7]+
135     ddxyz[2]*dxyz[0]*coefficient_metrique_ti[2]+dxyz[2]*ddxyz[0]*coefficient_metrique_ti[2]+dxyz[2]*dxyz[0]*coefficient_metrique_derive_ti[2]+
136     ddxyz[2]*dxyz[1]*coefficient_metrique_ti[4]+dxyz[2]*ddxyz[1]*coefficient_metrique_ti[4]+dxyz[2]*dxyz[1]*coefficient_metrique_derive_ti[4]+
137     ddxyz[2]*dxyz[2]*coefficient_metrique_ti[8]+dxyz[2]*ddxyz[2]*coefficient_metrique_ti[8]+dxyz[2]*dxyz[2]*coefficient_metrique_derive_ti[8];
138     double facteur3=dxyzii[0]*dxyzii[0]*coefficient_metrique_tii[0]+dxyzii[0]*dxyzii[1]*coefficient_metrique_tii[3]+dxyzii[0]*dxyzii[2]*coefficient_metrique_tii[6]+dxyzii[1]*dxyzii[0]*coefficient_metrique_tii[1]+dxyzii[1]*dxyzii[1]*coefficient_metrique_tii[4]+dxyzii[1]*dxyzii[2]*coefficient_metrique_tii[7]+dxyzii[2]*dxyzii[0]*coefficient_metrique_tii[2]+dxyzii[2]*dxyzii[1]*coefficient_metrique_tii[4]+dxyzii[2]*dxyzii[2]*coefficient_metrique_tii[8];
139     double residu=fabs(sqrt(facteur3)-sqrt(facteur1)-1./2./sqrt(facteur1)*facteur2*(tii-ti));
140     //if (residu>0.01) dt=dt/2.; else {dt=dt*2;bon_dt=1;}
141     //if (dt<1e-4) {dt=1e-4;bon_dt=1;}
142     //if (dt>0.25) {dt=0.25;bon_dt=1;}
143     bon_dt=1;
144    
145     }
146     while (bon_dt!=1);
147     #endif // ON_A_ABANDONNE_ET_ON_PREND_DT_EGALE_T2_MOINS_T1_DIVISE_PAR_1000
148    
149     t=0.7886751345*ti+0.2113248654*tii;
150     double coord[3];
151     double coefficient_metrique[9];
152     mgarete->deriver(t,coord);
153     if (creation_metrique) metrique->evaluer(&t,coefficient_metrique);
154     else
155     {
156     double xyz[3];
157     mgarete->evaluer(t,xyz);
158     metrique->evaluer(xyz,coefficient_metrique);
159     }
160     double facteur1=coord[0]*coord[0]*coefficient_metrique[0]+coord[0]*coord[1]*coefficient_metrique[3]+coord[0]*coord[2]*coefficient_metrique[6]+coord[1]*coord[0]*coefficient_metrique[1]+coord[1]*coord[1]*coefficient_metrique[4]+coord[1]*coord[2]*coefficient_metrique[7]+coord[2]*coord[0]*coefficient_metrique[2]+coord[2]*coord[1]*coefficient_metrique[4]+coord[2]*coord[2]*coefficient_metrique[8];
161     longueur_metrique=longueur_metrique+0.5*(tii-ti)*sqrt(facteur1);
162     t=0.7886751345*tii+0.2113248654*ti;
163     mgarete->deriver(t,coord);
164     if (creation_metrique) metrique->evaluer(&t,coefficient_metrique);
165     else
166     {
167     double xyz[3];
168     mgarete->evaluer(t,xyz);
169     metrique->evaluer(xyz,coefficient_metrique);
170     }
171     facteur1=coord[0]*coord[0]*coefficient_metrique[0]+coord[0]*coord[1]*coefficient_metrique[3]+coord[0]*coord[2]*coefficient_metrique[6]+coord[1]*coord[0]*coefficient_metrique[1]+coord[1]*coord[1]*coefficient_metrique[4]+coord[1]*coord[2]*coefficient_metrique[7]+coord[2]*coord[0]*coefficient_metrique[2]+coord[2]*coord[1]*coefficient_metrique[4]+coord[2]*coord[2]*coefficient_metrique[8];
172     longueur_metrique=longueur_metrique+0.5*(tii-ti)*sqrt(facteur1);
173     integrale pas;
174     pas.ti=ti;
175     pas.li=longueur_metrique;
176     tab_integrale.insert(tab_integrale.end(),pas);
177     }
178     while (tii<t2);
179     int nombre_de_segment=(int)floor(longueur_metrique);
180     if (longueur_metrique-floor(longueur_metrique)>0.5) nombre_de_segment++;
181     if (nombre_de_segment<1) nombre_de_segment=1;
182    
183     /* discretisation */
184    
185     double valeur_cible=longueur_metrique/nombre_de_segment;
186     int bon_noeud;
187     if (noeud_depart==NULL) bon_noeud=0; else bon_noeud=1;
188     int numnoeud=0;
189     while (bon_noeud==0)
190     {
191     noeud_depart=(MG_NOEUD*)mgarete->get_cosommet1()->get_sommet()->get_lien_maillage()->get(numnoeud);
192     MG_NOEUD* noeudtemp=mg_maillage->get_mg_noeudid(noeud_depart->get_id());
193     if (noeudtemp==NULL) numnoeud++; else bon_noeud=1;
194     }
195     if (noeud_arrivee==NULL) bon_noeud=0; else bon_noeud=1;
196     numnoeud=0;
197     while (bon_noeud==0)
198     {
199     noeud_arrivee=(MG_NOEUD*)mgarete->get_cosommet2()->get_sommet()->get_lien_maillage()->get(numnoeud);
200     MG_NOEUD* noeudtemp=mg_maillage->get_mg_noeudid(noeud_arrivee->get_id());
201     if (noeudtemp==NULL) numnoeud++; else bon_noeud=1;
202     }
203    
204    
205    
206     MG_NOEUD* noeud_precedent=noeud_depart;
207    
208     int nb_segment_cree=0;
209     int nb_pas_integrale=tab_integrale.size();
210     for (int i=0;i<nb_pas_integrale;i++)
211     {
212     while ((tab_integrale[i].li>(nb_segment_cree+1)*valeur_cible) && (nb_segment_cree<nombre_de_segment-1))
213     {
214     refresh();
215     double ti=tab_integrale[i].ti;
216     double tii=t2;
217     double li=0.;
218     if (i!=0) li=tab_integrale[i-1].li;
219     double lii=tab_integrale[i].li;
220     if (i!=nb_pas_integrale-1) tii=tab_integrale[i+1].ti;
221     t=ti+(tii-ti)/(lii-li)*((nb_segment_cree+1)*valeur_cible-li);
222     double coo[3];
223     mgarete->evaluer(t,coo);
224 francois 35 MG_NOEUD* nouveau_noeud=mg_maillage->ajouter_mg_noeud(mgarete,coo[0],coo[1],coo[2],MAILLEUR_AUTO);
225     mg_maillage->ajouter_mg_segment(mgarete,noeud_precedent,nouveau_noeud,MAILLEUR_AUTO);
226 5 noeud_precedent=nouveau_noeud;
227     nb_segment_cree++;
228     //ofstream o2("test_1D.mai",ios::out|ios::trunc);
229     //o2.precision(16);
230     //o2.setf(ios::showpoint);
231     //mg_maillage->enregistrer_sous_mesh_1D(o2);
232     }
233     }
234 francois 35 mg_maillage->ajouter_mg_segment(mgarete,noeud_precedent,noeud_arrivee,MAILLEUR_AUTO);
235 5
236    
237    
238     if (creation_metrique==1)
239     {
240     delete metrique;
241     metrique=NULL;
242     }
243     }
244    
245    
246    
247 francois 35 void MAILLEUR1D::adapte(void)
248     {
249     MG_NOEUD* noeud_dep;
250     MG_NOEUD* noeud_arr;
251     int nb=mg_arete->get_cosommet1()->get_sommet()->get_lien_maillage()->get_nb();
252     for (int i=0;i<nb;i++)
253     if (mg_maillage->get_mg_noeudid(mg_arete->get_cosommet1()->get_sommet()->get_lien_maillage()->get(i)->get_id())!=NULL)
254     noeud_dep=(MG_NOEUD*)mg_arete->get_cosommet1()->get_sommet()->get_lien_maillage()->get(i);
255     nb=mg_arete->get_cosommet2()->get_sommet()->get_lien_maillage()->get_nb();
256     for (int i=0;i<nb;i++)
257     if (mg_maillage->get_mg_noeudid(mg_arete->get_cosommet2()->get_sommet()->get_lien_maillage()->get(i)->get_id())!=NULL)
258     noeud_arr=(MG_NOEUD*)mg_arete->get_cosommet2()->get_sommet()->get_lien_maillage()->get(i);
259     int nb_segment=mg_arete->get_lien_maillage()->get_nb();
260     double tdeb=mg_arete->get_tmin();
261     double tfin=mg_arete->get_tmax();
262     std::map<double,MG_NOEUD*,less<double> > debut;
263     std::map<double,MG_NOEUD*,less<double> > fin;
264     std::pair<double,MG_NOEUD*> tmp(tfin,noeud_arr);
265     debut.insert(tmp);
266     std::pair<double,MG_NOEUD*> tmp2(tdeb,noeud_dep);
267     fin.insert(tmp2);
268     for (int i=0;i<nb_segment;i++)
269     {
270     MG_SEGMENT* seg=(MG_SEGMENT*)mg_arete->get_lien_maillage()->get(i);
271     if (mg_maillage->get_mg_segmentid(seg->get_id())==NULL) continue;
272     MG_NOEUD* noeud1=seg->get_noeud1();
273     MG_NOEUD* noeud2=seg->get_noeud2();
274     double *xyz1=noeud1->get_coord();
275     double *xyz2=noeud2->get_coord();
276     double t1,t2;
277     mg_arete->inverser(t1,xyz1);
278     mg_arete->inverser(t2,xyz2);
279 francois 61 if (mg_arete->get_courbe()->est_periodique())
280     {
281     if (t1<tdeb) t1=t1+mg_arete->get_courbe()->get_periode();
282     if (t2<tdeb) t2=t2+mg_arete->get_courbe()->get_periode();
283     }
284 francois 35 std::pair<double,MG_NOEUD*> tmp(t1,noeud1);
285     std::pair<double,MG_NOEUD*> tmp2(t2,noeud2);
286     debut.insert(tmp);
287     fin.insert(tmp2);
288     }
289     std::map<double,MG_NOEUD*,less<double> > :: iterator itdebut;
290     std::map<double,MG_NOEUD*,less<double> > :: iterator itfin;
291     itdebut=debut.begin();
292     itfin=fin.begin();
293     do
294     {
295     if (itdebut->second!=itfin->second) maille(mg_arete,(*itfin).first,(*itfin).second,(*itdebut).first,(*itdebut).second);
296     itdebut++;
297     itfin++;
298     }
299     while (itdebut!=debut.end());
300     }
301 5
302    
303    
304