ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/mailleur_auto/src/mailleur1d.cpp
Revision: 906
Committed: Mon Nov 13 22:30:18 2017 UTC (7 years, 9 months ago) by couturad
File size: 13407 byte(s)
Log Message:
Nouveau opencascade commit 1

File Contents

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