ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/mailleur/src/mailleur1d.cpp
Revision: 253
Committed: Tue Jul 13 19:40:46 2010 UTC (14 years, 10 months ago) by francois
File size: 13371 byte(s)
Log Message:
changement de hiearchie et utilisation de ccmake + mise a jour

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 francois 210 void MAILLEUR1D::maille(MG_GROUPE_TOPOLOGIQUE* mggt)
54 5 {
55     //afficheur << MAILLAGEARETE << endaff;
56     if (mg_arete!=NULL) maille(mg_arete);
57     else
58     {
59 francois 210 TPL_MAP_ENTITE<MG_ELEMENT_TOPOLOGIQUE*> lst;
60     if (mggt!=NULL)
61     {
62     int nb=mggt->get_nb();
63     for (int i=0;i<nb;i++)
64     {
65     lst.ajouter(mggt->get(i));
66     mggt->get(i)->get_topologie_sousjacente(&lst);
67     }
68     }
69 5 int nb_arete=mg_geometrie->get_nb_mg_arete();
70     for (int i=0;i<nb_arete;i++)
71     {
72     MG_ARETE* mgarete=mg_geometrie->get_mg_arete(i);
73 francois 210 if (mggt!=NULL)
74     if (lst.existe(mgarete)==0) continue;
75     maille(mgarete);
76 5 }
77     }
78     }
79    
80    
81     void MAILLEUR1D::maille(MG_ARETE* mgarete,double t1,MG_NOEUD* noeud_depart,double t2,MG_NOEUD* noeud_arrivee)
82     {
83     int creation_metrique=0;
84     if (metrique==NULL)
85     {
86     metrique=new FCT_TAILLE_ARETE(epsilon,distance_maximale,mgarete);
87     creation_metrique=1;
88     }
89    
90     /* calcul de la longueur de la arete dans la metrique */
91     if (fabs(t1)<1e-14)
92     if (fabs(t2)<1e-14)
93     {
94     t1=mgarete->get_tmin();
95     t2=mgarete->get_tmax();
96     }
97    
98     std::vector<integrale> tab_integrale;
99     double ti,tii,t;
100     double longueur_metrique=0.;
101     tii=t1;
102     double dt=0.001;
103     do
104     {
105     refresh();
106     double coefficient_metrique_ti[9];
107     double coefficient_metrique_tii[9];
108     double coefficient_metrique_derive_ti[9];
109     ti=tii;
110     tii=ti+(t2-t1)*dt;
111     if (tii>t2) tii=t2;
112    
113     #define ON_A_ABANDONNE_ET_ON_PREND_DT_EGALE_T2_MOINS_T1_DIVISE_PAR_1000
114     #ifndef ON_A_ABANDONNE_ET_ON_PREND_DT_EGALE_T2_MOINS_T1_DIVISE_PAR_1000
115     int bon_dt=0;
116     do
117     {
118     tii=ti+(t2-t1)*dt;
119     if (tii>t2) tii=t2;
120     double ddxyz[3];
121     double dxyz[3];
122     double xyz[3];
123     double dxyzii[3];
124     double xyzii[3];
125     mgarete->deriver_seconde(ti,ddxyz,dxyz,xyz);
126     mgarete->evaluer(tii,xyzii);
127     mgarete->deriver(tii,dxyzii);
128     if (creation_metrique)
129     {
130     metrique->evaluer(&ti,coefficient_metrique_ti);
131     metrique->evaluer(&tii,coefficient_metrique_tii);
132     metrique->deriver(&ti,coefficient_metrique_derive_ti);
133     }
134     else
135     {
136     metrique->evaluer(xyz,coefficient_metrique_ti);
137     metrique->evaluer(xyzii,coefficient_metrique_tii);
138     metrique->deriver(xyz,coefficient_metrique_derive_ti);
139     }
140     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];
141     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]+
142     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]+
143     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]+
144     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]+
145     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]+
146     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]+
147     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]+
148     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]+
149     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];
150     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];
151     double residu=fabs(sqrt(facteur3)-sqrt(facteur1)-1./2./sqrt(facteur1)*facteur2*(tii-ti));
152     //if (residu>0.01) dt=dt/2.; else {dt=dt*2;bon_dt=1;}
153     //if (dt<1e-4) {dt=1e-4;bon_dt=1;}
154     //if (dt>0.25) {dt=0.25;bon_dt=1;}
155     bon_dt=1;
156    
157     }
158     while (bon_dt!=1);
159     #endif // ON_A_ABANDONNE_ET_ON_PREND_DT_EGALE_T2_MOINS_T1_DIVISE_PAR_1000
160    
161     t=0.7886751345*ti+0.2113248654*tii;
162     double coord[3];
163     double coefficient_metrique[9];
164     mgarete->deriver(t,coord);
165     if (creation_metrique) metrique->evaluer(&t,coefficient_metrique);
166     else
167     {
168     double xyz[3];
169     mgarete->evaluer(t,xyz);
170     metrique->evaluer(xyz,coefficient_metrique);
171     }
172     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];
173     longueur_metrique=longueur_metrique+0.5*(tii-ti)*sqrt(facteur1);
174     t=0.7886751345*tii+0.2113248654*ti;
175     mgarete->deriver(t,coord);
176     if (creation_metrique) metrique->evaluer(&t,coefficient_metrique);
177     else
178     {
179     double xyz[3];
180     mgarete->evaluer(t,xyz);
181     metrique->evaluer(xyz,coefficient_metrique);
182     }
183     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];
184     longueur_metrique=longueur_metrique+0.5*(tii-ti)*sqrt(facteur1);
185     integrale pas;
186     pas.ti=ti;
187     pas.li=longueur_metrique;
188     tab_integrale.insert(tab_integrale.end(),pas);
189     }
190     while (tii<t2);
191     int nombre_de_segment=(int)floor(longueur_metrique);
192     if (longueur_metrique-floor(longueur_metrique)>0.5) nombre_de_segment++;
193     if (nombre_de_segment<1) nombre_de_segment=1;
194    
195     /* discretisation */
196    
197     double valeur_cible=longueur_metrique/nombre_de_segment;
198     int bon_noeud;
199     if (noeud_depart==NULL) bon_noeud=0; else bon_noeud=1;
200     int numnoeud=0;
201     while (bon_noeud==0)
202     {
203     noeud_depart=(MG_NOEUD*)mgarete->get_cosommet1()->get_sommet()->get_lien_maillage()->get(numnoeud);
204     MG_NOEUD* noeudtemp=mg_maillage->get_mg_noeudid(noeud_depart->get_id());
205     if (noeudtemp==NULL) numnoeud++; else bon_noeud=1;
206     }
207     if (noeud_arrivee==NULL) bon_noeud=0; else bon_noeud=1;
208     numnoeud=0;
209     while (bon_noeud==0)
210     {
211     noeud_arrivee=(MG_NOEUD*)mgarete->get_cosommet2()->get_sommet()->get_lien_maillage()->get(numnoeud);
212     MG_NOEUD* noeudtemp=mg_maillage->get_mg_noeudid(noeud_arrivee->get_id());
213     if (noeudtemp==NULL) numnoeud++; else bon_noeud=1;
214     }
215    
216    
217    
218     MG_NOEUD* noeud_precedent=noeud_depart;
219    
220     int nb_segment_cree=0;
221     int nb_pas_integrale=tab_integrale.size();
222     for (int i=0;i<nb_pas_integrale;i++)
223     {
224     while ((tab_integrale[i].li>(nb_segment_cree+1)*valeur_cible) && (nb_segment_cree<nombre_de_segment-1))
225     {
226     refresh();
227     double ti=tab_integrale[i].ti;
228     double tii=t2;
229     double li=0.;
230     if (i!=0) li=tab_integrale[i-1].li;
231     double lii=tab_integrale[i].li;
232     if (i!=nb_pas_integrale-1) tii=tab_integrale[i+1].ti;
233     t=ti+(tii-ti)/(lii-li)*((nb_segment_cree+1)*valeur_cible-li);
234     double coo[3];
235     mgarete->evaluer(t,coo);
236 francois 35 MG_NOEUD* nouveau_noeud=mg_maillage->ajouter_mg_noeud(mgarete,coo[0],coo[1],coo[2],MAILLEUR_AUTO);
237     mg_maillage->ajouter_mg_segment(mgarete,noeud_precedent,nouveau_noeud,MAILLEUR_AUTO);
238 5 noeud_precedent=nouveau_noeud;
239     nb_segment_cree++;
240     //ofstream o2("test_1D.mai",ios::out|ios::trunc);
241     //o2.precision(16);
242     //o2.setf(ios::showpoint);
243     //mg_maillage->enregistrer_sous_mesh_1D(o2);
244     }
245     }
246 francois 35 mg_maillage->ajouter_mg_segment(mgarete,noeud_precedent,noeud_arrivee,MAILLEUR_AUTO);
247 5
248    
249    
250     if (creation_metrique==1)
251     {
252     delete metrique;
253     metrique=NULL;
254     }
255     }
256    
257    
258    
259 francois 35 void MAILLEUR1D::adapte(void)
260     {
261     MG_NOEUD* noeud_dep;
262     MG_NOEUD* noeud_arr;
263     int nb=mg_arete->get_cosommet1()->get_sommet()->get_lien_maillage()->get_nb();
264     for (int i=0;i<nb;i++)
265     if (mg_maillage->get_mg_noeudid(mg_arete->get_cosommet1()->get_sommet()->get_lien_maillage()->get(i)->get_id())!=NULL)
266     noeud_dep=(MG_NOEUD*)mg_arete->get_cosommet1()->get_sommet()->get_lien_maillage()->get(i);
267     nb=mg_arete->get_cosommet2()->get_sommet()->get_lien_maillage()->get_nb();
268     for (int i=0;i<nb;i++)
269     if (mg_maillage->get_mg_noeudid(mg_arete->get_cosommet2()->get_sommet()->get_lien_maillage()->get(i)->get_id())!=NULL)
270     noeud_arr=(MG_NOEUD*)mg_arete->get_cosommet2()->get_sommet()->get_lien_maillage()->get(i);
271     int nb_segment=mg_arete->get_lien_maillage()->get_nb();
272     double tdeb=mg_arete->get_tmin();
273     double tfin=mg_arete->get_tmax();
274 francois 102 std::map<double,MG_NOEUD*,std::less<double> > debut;
275     std::map<double,MG_NOEUD*,std::less<double> > fin;
276 francois 35 std::pair<double,MG_NOEUD*> tmp(tfin,noeud_arr);
277     debut.insert(tmp);
278     std::pair<double,MG_NOEUD*> tmp2(tdeb,noeud_dep);
279     fin.insert(tmp2);
280     for (int i=0;i<nb_segment;i++)
281     {
282     MG_SEGMENT* seg=(MG_SEGMENT*)mg_arete->get_lien_maillage()->get(i);
283     if (mg_maillage->get_mg_segmentid(seg->get_id())==NULL) continue;
284     MG_NOEUD* noeud1=seg->get_noeud1();
285     MG_NOEUD* noeud2=seg->get_noeud2();
286     double t1,t2;
287 francois 164 if (noeud1==noeud_dep) t1=tdeb;
288     else if (noeud1==noeud_arr) t1=tfin;
289     else
290     {
291     double *xyz1=noeud1->get_coord();
292     mg_arete->inverser(t1,xyz1);
293     }
294     if (noeud2==noeud_dep) t2=tdeb;
295     else if (noeud2==noeud_arr) t2=tfin;
296     else
297     {
298     double *xyz2=noeud2->get_coord();
299     mg_arete->inverser(t2,xyz2);
300     }
301 francois 61 if (mg_arete->get_courbe()->est_periodique())
302     {
303     if (t1<tdeb) t1=t1+mg_arete->get_courbe()->get_periode();
304 francois 164 if (t2<t1) t2=t2+mg_arete->get_courbe()->get_periode();
305 francois 61 }
306 francois 35 std::pair<double,MG_NOEUD*> tmp(t1,noeud1);
307     std::pair<double,MG_NOEUD*> tmp2(t2,noeud2);
308     debut.insert(tmp);
309     fin.insert(tmp2);
310     }
311 francois 102 std::map<double,MG_NOEUD*,std::less<double> > :: iterator itdebut;
312     std::map<double,MG_NOEUD*,std::less<double> > :: iterator itfin;
313 francois 35 itdebut=debut.begin();
314     itfin=fin.begin();
315     do
316     {
317 francois 164 if ((itdebut->second!=itfin->second) || (nb_segment==0) ) maille(mg_arete,(*itfin).first,(*itfin).second,(*itdebut).first,(*itdebut).second);
318 francois 35 itdebut++;
319     itfin++;
320     }
321     while (itdebut!=debut.end());
322     }
323 5
324    
325    
326