ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/CAD4FE/src/CAD4FE_mailleur1d.cpp
Revision: 569
Committed: Thu Oct 16 14:36:31 2014 UTC (10 years, 6 months ago) by foucault
File size: 16558 byte(s)
Log Message:
Mise à jour pour CAD4FE (Gilles) : operation 1 (tentative)

File Contents

# User Rev Content
1 foucault 27 //------------------------------------------------------------
2     //------------------------------------------------------------
3     // MAGiC
4     // Jean Christophe Cuillière et Vincent FRANCOIS
5     // Département de Génie Mécanique - UQTR
6     //------------------------------------------------------------
7     // Le projet MAGIC est un projet de recherche du département
8     // de génie mécanique de l'Université du Québec à
9     // Trois Rivières
10     // Les librairies ne peuvent être utilisées 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 "CAD4FE_mailleur1d.h"
29    
30     #include "CAD4FE_MCNode.h"
31     #include "CAD4FE_MCSegment.h"
32     #include "CAD4FE_MCEdge.h"
33     #include "CAD4FE_MCVertex.h"
34     #include "CAD4FE_PolyCurve.h"
35     //#include "affiche.h"
36     //#include "message.h"
37     #include <fstream>
38    
39     using namespace CAD4FE;
40    
41 foucault 569 struct integrale
42 foucault 27 {
43     double ti;
44     double li;
45     };
46    
47    
48     MAILLEUR1D::MAILLEUR1D(MG_MAILLAGE* mgmai,MG_GEOMETRIE *mggeo,MCEdge* __mcEdge,FCT_TAILLE* fct_taille):MAILLEUR(),mg_maillage(mgmai),mg_geometrie(mggeo),mcEdge(__mcEdge),metrique(fct_taille),debug(0)
49     {
50     }
51    
52    
53    
54     MAILLEUR1D::~MAILLEUR1D()
55     {
56     }
57    
58     MG_SOMMET * MAILLEUR1D::IsPointInRefVertex(PolyCurve * __polycurve, double __nodeCoord[3], double __relativeTolerance)
59     {
60     int nbRefV=__polycurve->GetRefEdgeCount();
61     double sqAbsoluteTolerance = __relativeTolerance*__polycurve->get_longueur();
62     sqAbsoluteTolerance *= sqAbsoluteTolerance;
63     OT_VECTEUR_3D xNode(__nodeCoord);
64     int i;
65     for (i=0;i<nbRefV;i++)
66     {
67     OT_VECTEUR_3D xRefVertex(0,0,0);
68     MG_SOMMET * vertex = __polycurve->GetRefVertex(i);
69     vertex->get_point()->evaluer(xRefVertex);
70     OT_VECTEUR_3D dx = xRefVertex-xNode;
71     double sqdistNodeToVertex = dx.get_longueur2();
72     if (sqdistNodeToVertex < sqAbsoluteTolerance)
73     {
74     return vertex;
75     }
76     }
77     return NULL;
78     }
79    
80     void MAILLEUR1D::maille()
81     {
82 foucault 64 // If this edge is a point isolated in a face, then quit
83     if (mcEdge->GetPolyCurve()->get_longueur() == 0)
84     return;
85    
86 foucault 27 int creation_metrique=0;
87     if (metrique==NULL)
88     {
89 foucault 569 printf("Erreur : metrique = NULL\n");
90     return;
91     /* metrique=new FCT_TAILLE_ARETE(epsilon,distance_maximale,mcEdge);
92     creation_metrique=1;*/
93 foucault 27 }
94    
95     /* calcul de la longueur de la arete dans la metrique */
96     double t1=mcEdge->get_cosommet1()->get_t();
97     double t2=mcEdge->get_cosommet2()->get_t();
98     if ((t2<t1+0.0001)&&(mcEdge->get_courbe()->est_periodique())) t2=t2+mcEdge->get_courbe()->get_periode();
99    
100 foucault 569 std::vector<double> tab_integrale_ti, tab_integrale_li;
101 foucault 27 double ti,tii,t;
102     double longueur_metrique=0.;
103     tii=t1;
104     double dt=0.001;
105     do
106     {
107     double coefficient_metrique_ti[9];
108     double coefficient_metrique_tii[9];
109     double coefficient_metrique_derive_ti[9];
110     ti=tii;
111     tii=ti+(t2-t1)*dt;
112     if (tii>t2) tii=t2;
113    
114     #define ON_A_ABANDONNE_ET_ON_PREND_DT_EGALE_T2_MOINS_T1_DIVISE_PAR_1000
115     #ifndef ON_A_ABANDONNE_ET_ON_PREND_DT_EGALE_T2_MOINS_T1_DIVISE_PAR_1000
116     int bon_dt=0;
117     do
118     {
119     tii=ti+(t2-t1)*dt;
120     if (tii>t2) tii=t2;
121     double ddxyz[3];
122     double dxyz[3];
123     double xyz[3];
124     double dxyzii[3];
125     double xyzii[3];
126     mcEdge->deriver_seconde(ti,ddxyz,dxyz,xyz);
127     mcEdge->evaluer(tii,xyzii);
128     mcEdge->deriver(tii,dxyzii);
129     if (creation_metrique)
130     {
131     metrique->evaluer(&ti,coefficient_metrique_ti);
132     metrique->evaluer(&tii,coefficient_metrique_tii);
133     metrique->deriver(&ti,coefficient_metrique_derive_ti);
134     }
135     else
136     {
137     metrique->evaluer(xyz,coefficient_metrique_ti);
138     metrique->evaluer(xyzii,coefficient_metrique_tii);
139     metrique->deriver(xyz,coefficient_metrique_derive_ti);
140     }
141     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];
142     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]+
143     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]+
144     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]+
145     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]+
146     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]+
147     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]+
148     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]+
149     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]+
150     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];
151 foucault 569
152    
153 foucault 27 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];
154     double residu=fabs(sqrt(facteur3)-sqrt(facteur1)-1./2./sqrt(facteur1)*facteur2*(tii-ti));
155     //if (residu>0.01) dt=dt/2.; else {dt=dt*2;bon_dt=1;}
156     //if (dt<1e-4) {dt=1e-4;bon_dt=1;}
157     //if (dt>0.25) {dt=0.25;bon_dt=1;}
158     bon_dt=1;
159    
160     }
161     while (bon_dt!=1);
162     #endif // ON_A_ABANDONNE_ET_ON_PREND_DT_EGALE_T2_MOINS_T1_DIVISE_PAR_1000
163    
164     t=0.7886751345*ti+0.2113248654*tii;
165     double coord[3];
166     double coefficient_metrique[9];
167     mcEdge->deriver(t,coord);
168 foucault 569
169 foucault 27 if (creation_metrique) metrique->evaluer(&t,coefficient_metrique);
170     else
171     {
172     double xyz[3];
173     mcEdge->evaluer(t,xyz);
174     metrique->evaluer(xyz,coefficient_metrique);
175     }
176     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];
177     longueur_metrique=longueur_metrique+0.5*(tii-ti)*sqrt(facteur1);
178     t=0.7886751345*tii+0.2113248654*ti;
179 foucault 569 mcEdge->deriver(t,coord);
180    
181 foucault 27 if (creation_metrique) metrique->evaluer(&t,coefficient_metrique);
182     else
183     {
184     double xyz[3];
185     mcEdge->evaluer(t,xyz);
186     metrique->evaluer(xyz,coefficient_metrique);
187 foucault 569 }
188 foucault 27 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];
189 foucault 569
190 foucault 27 longueur_metrique=longueur_metrique+0.5*(tii-ti)*sqrt(facteur1);
191 foucault 569
192     tab_integrale_ti.push_back(ti);
193     tab_integrale_li.push_back(longueur_metrique);
194 foucault 27 }
195     while (tii<t2);
196     int nombre_de_segment=(int)floor(longueur_metrique);
197     if (longueur_metrique-floor(longueur_metrique)>0.5) nombre_de_segment++;
198     if (nombre_de_segment<1) nombre_de_segment=1;
199    
200     /* discretisation */
201    
202     double valeur_cible=longueur_metrique/nombre_de_segment;
203     int bon_noeud=0;
204     int numnoeud=0;
205 foucault 569 MG_NOEUD* mg_noeud_depart;
206 foucault 27 while (bon_noeud==0)
207     {
208 foucault 569 mg_noeud_depart=(MG_NOEUD*)mcEdge->get_cosommet1()->get_sommet()->get_lien_maillage()->get(numnoeud);
209     MG_NOEUD * mg_noeudtemp=mg_maillage->get_mg_noeudid(mg_noeud_depart->get_id());
210     if (mg_noeudtemp==NULL) numnoeud++; else bon_noeud=1;
211 foucault 27 }
212     bon_noeud=0;
213     numnoeud=0;
214 foucault 569 MG_NOEUD* mg_noeud_arrivee;
215 foucault 27 while (bon_noeud==0)
216     {
217 foucault 569 mg_noeud_arrivee=(MG_NOEUD*)mcEdge->get_cosommet2()->get_sommet()->get_lien_maillage()->get(numnoeud);
218     MG_NOEUD* mg_noeudtemp=mg_maillage->get_mg_noeudid(mg_noeud_arrivee->get_id());
219     if (mg_noeudtemp==NULL) numnoeud++; else bon_noeud=1;
220 foucault 27 }
221 foucault 569 MCNode * noeud_depart = (MCNode*) mg_noeud_depart;
222     MCNode * noeud_arrivee = (MCNode*) mg_noeud_arrivee;
223 foucault 27
224     MCNode* noeud_precedent=noeud_depart;
225    
226     int nb_segment_cree=0;
227 foucault 569 int nb_pas_integrale=tab_integrale_ti.size();
228 foucault 27 for (int i=0;i<nb_pas_integrale;i++)
229     {
230 foucault 569 while ((tab_integrale_li[i]>(nb_segment_cree+1)*valeur_cible) && (nb_segment_cree<nombre_de_segment-1))
231 foucault 27 {
232 foucault 569 // refresh();
233     double ti=tab_integrale_ti[i];
234 foucault 27 double tii=t2;
235     double li=0.;
236 foucault 569 if (i!=0) li=tab_integrale_li[i-1];
237     double lii=tab_integrale_li[i];
238     if (i!=nb_pas_integrale-1) tii=tab_integrale_ti[i+1];
239 foucault 27 t=ti+(tii-ti)/(lii-li)*((nb_segment_cree+1)*valeur_cible-li);
240     double coo[3];
241     mcEdge->evaluer(t,coo);
242    
243     // arete de référence
244     MG_ARETE * nouveau_noeud_refEdge = NULL;
245     // parametre de l'arete de reference
246     double nouveau_noeud_refEdgeT = 0, Parameter_SToRefEdgeDerT;
247     mcEdge->GetPolyCurve()->Parameter_SToRefEdgeT(t,&nouveau_noeud_refEdge, &nouveau_noeud_refEdgeT, &Parameter_SToRefEdgeDerT, false);
248     // reference topology mapping
249     MCNode * nouveau_noeud = NULL;
250     MG_SOMMET * refVertex = IsPointInRefVertex(mcEdge->GetPolyCurve(),coo,1E-6);
251     if (refVertex == NULL)
252     {
253     nouveau_noeud = new MCNode(mcEdge,nouveau_noeud_refEdge,nouveau_noeud_refEdgeT,coo);
254     }
255     else
256     {
257     nouveau_noeud = new MCNode(mcEdge,refVertex,coo[0],coo[1],coo[2]);
258     }
259 foucault 569 mg_maillage->ajouter_mg_noeud(nouveau_noeud);
260     nouveau_noeud->Creator = MCNode::mailleur_1d;
261 foucault 27 // création d'un segment de MC arête (arête composite)
262     MCSegment * nouveau_segment = new MCSegment (mcEdge,noeud_precedent,nouveau_noeud);
263    
264     mg_maillage->ajouter_mg_segment(nouveau_segment);
265     noeud_precedent=nouveau_noeud;
266     nb_segment_cree++;
267     }
268     }
269    
270     MCSegment * nouveau_segment = new MCSegment (mcEdge,noeud_precedent,noeud_arrivee);
271    
272     mg_maillage->ajouter_mg_segment(nouveau_segment);
273    
274     if (debug)
275     {
276     char fichierProgression1D [1024];
277     sprintf(fichierProgression1D, "%s%s", "c:\\temp\\void", "_1D.mai");
278     ofstream o3(fichierProgression1D,ios::out|ios::trunc);
279     o3.precision(16);
280     o3.setf(ios::showpoint);
281     mg_maillage->enregistrer_sous_mesh_1D(o3);
282     }
283    
284     if (creation_metrique==1)
285     {
286     delete metrique;
287     metrique=NULL;
288     }
289    
290     }
291    
292 foucault 569 void MAILLEUR1D::maille_points(std::vector <OT_VECTEUR_3D> & __points)
293     {
294     // If this edge is a point isolated in a face, then quit
295     if (mcEdge->GetPolyCurve()->get_longueur() == 0)
296     return;
297 foucault 27
298 foucault 569 /* calcul de la longueur de la arete dans la metrique */
299     double t1=mcEdge->get_cosommet1()->get_t();
300     double t2=mcEdge->get_cosommet2()->get_t();
301     if ((t2<t1+0.0001)&&(mcEdge->get_courbe()->est_periodique())) t2=t2+mcEdge->get_courbe()->get_periode();
302 foucault 27
303 foucault 569 std::set <double> tab_t;
304     std::set <double>::iterator tab_t_it1, tab_t_it2, tab_t_it;
305 foucault 27
306 foucault 569 tab_t.insert(mcEdge->get_cosommet1()->get_t());
307     tab_t.insert(mcEdge->get_cosommet2()->get_t());
308 foucault 27
309 foucault 569 for (int i=0; i<__points.size(); i++)
310     {
311     double t;
312     mcEdge->inverser(t, __points[i]);
313     OT_VECTEUR_3D xyz;
314     mcEdge->evaluer(t, xyz);
315     double distance = (xyz-__points[i]).get_longueur();
316     tab_t_it1 = tab_t.lower_bound(t);
317     tab_t_it2 = tab_t.upper_bound(t);
318     double t_dist = 1E99;
319     for (tab_t_it = tab_t_it1; tab_t_it != tab_t.end(); tab_t_it++)
320     {
321     t_dist = std::min(t_dist, fabs(t-*tab_t_it));
322     if (tab_t_it == tab_t_it2) break;
323     }
324     if (distance > 0.001 * mcEdge->GetPolyCurve()->get_longueur())
325     {
326     printf("node #%d/%lu is too far from edge %lu\n", i, __points.size(), mcEdge->get_id());
327     continue;
328     }
329     if (t_dist < 0.001 * t2-t1)
330     {
331     printf("node #%d/%lu already exists in edge %lu\n", i, __points.size(), mcEdge->get_id());
332     continue;
333     }
334 foucault 27
335 foucault 569 tab_t.insert(t);
336     }
337 foucault 27
338 foucault 569 /* discretisation */
339     int bon_noeud=0;
340     int numnoeud=0;
341     MG_NOEUD* mg_noeud_depart;
342     while (bon_noeud==0)
343     {
344     mg_noeud_depart=(MG_NOEUD*)mcEdge->get_cosommet1()->get_sommet()->get_lien_maillage()->get(numnoeud);
345     MG_NOEUD * mg_noeudtemp=mg_maillage->get_mg_noeudid(mg_noeud_depart->get_id());
346     if (mg_noeudtemp==NULL) numnoeud++; else bon_noeud=1;
347     }
348     bon_noeud=0;
349     numnoeud=0;
350     MG_NOEUD* mg_noeud_arrivee;
351     while (bon_noeud==0)
352     {
353     mg_noeud_arrivee=(MG_NOEUD*)mcEdge->get_cosommet2()->get_sommet()->get_lien_maillage()->get(numnoeud);
354     MG_NOEUD* mg_noeudtemp=mg_maillage->get_mg_noeudid(mg_noeud_arrivee->get_id());
355     if (mg_noeudtemp==NULL) numnoeud++; else bon_noeud=1;
356     }
357     MCNode * noeud_depart = (MCNode*) mg_noeud_depart;
358     MCNode * noeud_arrivee = (MCNode*) mg_noeud_arrivee;
359    
360     MCNode* noeud_precedent=noeud_depart;
361    
362     int nb_segment_cree=0;
363     tab_t_it2 = tab_t.begin();
364     for (tab_t_it1 = tab_t.begin(); nb_segment_cree+2 < tab_t.size(); tab_t_it1++)
365     {
366     tab_t_it2++;
367     if (tab_t_it2 == tab_t.end()) break;
368    
369     double ti=*tab_t_it1;
370     double tii=*tab_t_it2;
371    
372     // arete de référence
373     MG_ARETE * nouveau_noeud_refEdge = NULL;
374     // parametre de l'arete de reference
375     double nouveau_noeud_refEdgeT = 0, Parameter_SToRefEdgeDerT;
376     mcEdge->GetPolyCurve()->Parameter_SToRefEdgeT(tii,&nouveau_noeud_refEdge, &nouveau_noeud_refEdgeT, &Parameter_SToRefEdgeDerT, false);
377     // reference topology mapping
378     MCNode * nouveau_noeud = NULL;
379     double coo[3]; mcEdge->evaluer(tii, coo);
380     MG_SOMMET * refVertex = IsPointInRefVertex(mcEdge->GetPolyCurve(),coo,1E-6);
381     if (refVertex == NULL)
382     {
383     nouveau_noeud = new MCNode(mcEdge,nouveau_noeud_refEdge,nouveau_noeud_refEdgeT,coo);
384     }
385     else
386     {
387     nouveau_noeud = new MCNode(mcEdge,refVertex,coo[0],coo[1],coo[2]);
388     }
389     mg_maillage->ajouter_mg_noeud(nouveau_noeud);
390     nouveau_noeud->Creator = MCNode::mailleur_1d;
391     // création d'un segment de MC arête (arête composite)
392     MCSegment * nouveau_segment = new MCSegment (mcEdge,noeud_precedent,nouveau_noeud);
393    
394     mg_maillage->ajouter_mg_segment(nouveau_segment);
395     noeud_precedent=nouveau_noeud;
396     nb_segment_cree++;
397     }
398    
399     MCSegment * nouveau_segment = new MCSegment (mcEdge,noeud_precedent,noeud_arrivee);
400    
401     mg_maillage->ajouter_mg_segment(nouveau_segment);
402    
403     if (debug)
404     {
405     char fichierProgression1D [1024];
406     sprintf(fichierProgression1D, "%s%s", "c:\\temp\\void", "_1D.mai");
407     ofstream o3(fichierProgression1D,ios::out|ios::trunc);
408     o3.precision(16);
409     o3.setf(ios::showpoint);
410     mg_maillage->enregistrer_sous_mesh_1D(o3);
411     }
412     }