ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/mailleur/src/mailleur1d.cpp
Revision: 966
Committed: Thu Sep 6 16:46:34 2018 UTC (6 years, 8 months ago) by couturad
Original Path: magic/lib/mailleur_auto/src/mailleur1d.cpp
File size: 13502 byte(s)
Log Message:
Ajout de l'histogramme a MAGIC_PLOT
Ajout d'une sortie OK ou FAIL (int) au MAILLEUR afin de gerer certaines exceptions
Ajout d'une phase RSA a la fin du generateur DCR

File Contents

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