ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/mailleur_auto/src/mailleur1d.cpp
Revision: 5
Committed: Tue Jun 12 20:26:34 2007 UTC (17 years, 11 months ago)
Original Path: magic/lib/mailleur/mailleur/src/mailleur1d.cpp
File size: 10368 byte(s)
Log Message:

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     //#include "affiche.h"
30     //#include "message.h"
31     #include <fstream>
32     struct integrale
33     {
34     double ti;
35     double li;
36     };
37    
38    
39     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)
40     {
41     }
42    
43    
44    
45     MAILLEUR1D::~MAILLEUR1D()
46     {
47     }
48    
49    
50    
51    
52     void MAILLEUR1D::maille(void)
53     {
54     //afficheur << MAILLAGEARETE << endaff;
55     if (mg_arete!=NULL) maille(mg_arete);
56     else
57     {
58     int nb_arete=mg_geometrie->get_nb_mg_arete();
59     for (int i=0;i<nb_arete;i++)
60     {
61     MG_ARETE* mgarete=mg_geometrie->get_mg_arete(i);
62     maille(mgarete);
63     }
64     }
65     }
66    
67    
68     void MAILLEUR1D::maille(MG_ARETE* mgarete,double t1,MG_NOEUD* noeud_depart,double t2,MG_NOEUD* noeud_arrivee)
69     {
70     int creation_metrique=0;
71     if (metrique==NULL)
72     {
73     metrique=new FCT_TAILLE_ARETE(epsilon,distance_maximale,mgarete);
74     creation_metrique=1;
75     }
76    
77     /* calcul de la longueur de la arete dans la metrique */
78     if (fabs(t1)<1e-14)
79     if (fabs(t2)<1e-14)
80     {
81     t1=mgarete->get_tmin();
82     t2=mgarete->get_tmax();
83     }
84    
85     std::vector<integrale> tab_integrale;
86     double ti,tii,t;
87     double longueur_metrique=0.;
88     tii=t1;
89     double dt=0.001;
90     do
91     {
92     refresh();
93     double coefficient_metrique_ti[9];
94     double coefficient_metrique_tii[9];
95     double coefficient_metrique_derive_ti[9];
96     ti=tii;
97     tii=ti+(t2-t1)*dt;
98     if (tii>t2) tii=t2;
99    
100     #define ON_A_ABANDONNE_ET_ON_PREND_DT_EGALE_T2_MOINS_T1_DIVISE_PAR_1000
101     #ifndef ON_A_ABANDONNE_ET_ON_PREND_DT_EGALE_T2_MOINS_T1_DIVISE_PAR_1000
102     int bon_dt=0;
103     do
104     {
105     tii=ti+(t2-t1)*dt;
106     if (tii>t2) tii=t2;
107     double ddxyz[3];
108     double dxyz[3];
109     double xyz[3];
110     double dxyzii[3];
111     double xyzii[3];
112     mgarete->deriver_seconde(ti,ddxyz,dxyz,xyz);
113     mgarete->evaluer(tii,xyzii);
114     mgarete->deriver(tii,dxyzii);
115     if (creation_metrique)
116     {
117     metrique->evaluer(&ti,coefficient_metrique_ti);
118     metrique->evaluer(&tii,coefficient_metrique_tii);
119     metrique->deriver(&ti,coefficient_metrique_derive_ti);
120     }
121     else
122     {
123     metrique->evaluer(xyz,coefficient_metrique_ti);
124     metrique->evaluer(xyzii,coefficient_metrique_tii);
125     metrique->deriver(xyz,coefficient_metrique_derive_ti);
126     }
127     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];
128     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]+
129     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]+
130     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]+
131     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]+
132     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]+
133     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]+
134     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]+
135     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]+
136     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];
137     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];
138     double residu=fabs(sqrt(facteur3)-sqrt(facteur1)-1./2./sqrt(facteur1)*facteur2*(tii-ti));
139     //if (residu>0.01) dt=dt/2.; else {dt=dt*2;bon_dt=1;}
140     //if (dt<1e-4) {dt=1e-4;bon_dt=1;}
141     //if (dt>0.25) {dt=0.25;bon_dt=1;}
142     bon_dt=1;
143    
144     }
145     while (bon_dt!=1);
146     #endif // ON_A_ABANDONNE_ET_ON_PREND_DT_EGALE_T2_MOINS_T1_DIVISE_PAR_1000
147    
148     t=0.7886751345*ti+0.2113248654*tii;
149     double coord[3];
150     double coefficient_metrique[9];
151     mgarete->deriver(t,coord);
152     if (creation_metrique) metrique->evaluer(&t,coefficient_metrique);
153     else
154     {
155     double xyz[3];
156     mgarete->evaluer(t,xyz);
157     metrique->evaluer(xyz,coefficient_metrique);
158     }
159     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];
160     longueur_metrique=longueur_metrique+0.5*(tii-ti)*sqrt(facteur1);
161     t=0.7886751345*tii+0.2113248654*ti;
162     mgarete->deriver(t,coord);
163     if (creation_metrique) metrique->evaluer(&t,coefficient_metrique);
164     else
165     {
166     double xyz[3];
167     mgarete->evaluer(t,xyz);
168     metrique->evaluer(xyz,coefficient_metrique);
169     }
170     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];
171     longueur_metrique=longueur_metrique+0.5*(tii-ti)*sqrt(facteur1);
172     integrale pas;
173     pas.ti=ti;
174     pas.li=longueur_metrique;
175     tab_integrale.insert(tab_integrale.end(),pas);
176     }
177     while (tii<t2);
178     int nombre_de_segment=(int)floor(longueur_metrique);
179     if (longueur_metrique-floor(longueur_metrique)>0.5) nombre_de_segment++;
180     if (nombre_de_segment<1) nombre_de_segment=1;
181    
182     /* discretisation */
183    
184     double valeur_cible=longueur_metrique/nombre_de_segment;
185     int bon_noeud;
186     if (noeud_depart==NULL) bon_noeud=0; else bon_noeud=1;
187     int numnoeud=0;
188     while (bon_noeud==0)
189     {
190     noeud_depart=(MG_NOEUD*)mgarete->get_cosommet1()->get_sommet()->get_lien_maillage()->get(numnoeud);
191     MG_NOEUD* noeudtemp=mg_maillage->get_mg_noeudid(noeud_depart->get_id());
192     if (noeudtemp==NULL) numnoeud++; else bon_noeud=1;
193     }
194     if (noeud_arrivee==NULL) bon_noeud=0; else bon_noeud=1;
195     numnoeud=0;
196     while (bon_noeud==0)
197     {
198     noeud_arrivee=(MG_NOEUD*)mgarete->get_cosommet2()->get_sommet()->get_lien_maillage()->get(numnoeud);
199     MG_NOEUD* noeudtemp=mg_maillage->get_mg_noeudid(noeud_arrivee->get_id());
200     if (noeudtemp==NULL) numnoeud++; else bon_noeud=1;
201     }
202    
203    
204    
205     MG_NOEUD* noeud_precedent=noeud_depart;
206    
207     int nb_segment_cree=0;
208     int nb_pas_integrale=tab_integrale.size();
209     for (int i=0;i<nb_pas_integrale;i++)
210     {
211     while ((tab_integrale[i].li>(nb_segment_cree+1)*valeur_cible) && (nb_segment_cree<nombre_de_segment-1))
212     {
213     refresh();
214     double ti=tab_integrale[i].ti;
215     double tii=t2;
216     double li=0.;
217     if (i!=0) li=tab_integrale[i-1].li;
218     double lii=tab_integrale[i].li;
219     if (i!=nb_pas_integrale-1) tii=tab_integrale[i+1].ti;
220     t=ti+(tii-ti)/(lii-li)*((nb_segment_cree+1)*valeur_cible-li);
221     double coo[3];
222     mgarete->evaluer(t,coo);
223     MG_NOEUD* nouveau_noeud=mg_maillage->ajouter_mg_noeud(mgarete,coo[0],coo[1],coo[2]);
224     mg_maillage->ajouter_mg_segment(mgarete,noeud_precedent,nouveau_noeud);
225     noeud_precedent=nouveau_noeud;
226     nb_segment_cree++;
227     //ofstream o2("test_1D.mai",ios::out|ios::trunc);
228     //o2.precision(16);
229     //o2.setf(ios::showpoint);
230     //mg_maillage->enregistrer_sous_mesh_1D(o2);
231     }
232     }
233     mg_maillage->ajouter_mg_segment(mgarete,noeud_precedent,noeud_arrivee);
234    
235    
236    
237     if (creation_metrique==1)
238     {
239     delete metrique;
240     metrique=NULL;
241     }
242     }
243    
244    
245    
246    
247    
248    
249