ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/CAD4FE/src/CAD4FE_mailleur1d.cpp
Revision: 1158
Committed: Thu Jun 13 22:18:49 2024 UTC (11 months ago) by francois
File size: 16271 byte(s)
Log Message: 
compatibilité Ubuntu 22.04
Suppression des refeences à Windows
Ajout d'une banière

File Contents

# User Rev Content
1 francois 1158 //####//------------------------------------------------------------
2     //####//------------------------------------------------------------
3     //####// MAGiC
4     //####// Jean Christophe Cuilliere et Vincent FRANCOIS
5     //####// Departement de Genie Mecanique - UQTR
6     //####//------------------------------------------------------------
7     //####// MAGIC est un projet de recherche de l equipe ERICCA
8     //####// du departement de genie mecanique de l Universite du Quebec a Trois Rivieres
9     //####// http://www.uqtr.ca/ericca
10     //####// http://www.uqtr.ca/
11     //####//------------------------------------------------------------
12     //####//------------------------------------------------------------
13     //####//
14     //####// CAD4FE_mailleur1d.cpp
15     //####//
16     //####//------------------------------------------------------------
17     //####//------------------------------------------------------------
18     //####// COPYRIGHT 2000-2024
19     //####// jeu 13 jun 2024 11:58:56 EDT
20     //####//------------------------------------------------------------
21     //####//------------------------------------------------------------
22 foucault 27
23    
24     #include "gestionversion.h"
25     #include <math.h>
26     #include <vector>
27     #include "CAD4FE_mailleur1d.h"
28    
29     #include "CAD4FE_MCNode.h"
30     #include "CAD4FE_MCSegment.h"
31     #include "CAD4FE_MCEdge.h"
32     #include "CAD4FE_MCVertex.h"
33     #include "CAD4FE_PolyCurve.h"
34     #include <fstream>
35    
36     using namespace CAD4FE;
37    
38 foucault 569 struct integrale
39 foucault 27 {
40     double ti;
41     double li;
42     };
43    
44    
45     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)
46     {
47     }
48    
49    
50    
51     MAILLEUR1D::~MAILLEUR1D()
52     {
53     }
54    
55     MG_SOMMET * MAILLEUR1D::IsPointInRefVertex(PolyCurve * __polycurve, double __nodeCoord[3], double __relativeTolerance)
56     {
57     int nbRefV=__polycurve->GetRefEdgeCount();
58     double sqAbsoluteTolerance = __relativeTolerance*__polycurve->get_longueur();
59     sqAbsoluteTolerance *= sqAbsoluteTolerance;
60     OT_VECTEUR_3D xNode(__nodeCoord);
61     int i;
62     for (i=0;i<nbRefV;i++)
63     {
64     OT_VECTEUR_3D xRefVertex(0,0,0);
65     MG_SOMMET * vertex = __polycurve->GetRefVertex(i);
66     vertex->get_point()->evaluer(xRefVertex);
67     OT_VECTEUR_3D dx = xRefVertex-xNode;
68     double sqdistNodeToVertex = dx.get_longueur2();
69     if (sqdistNodeToVertex < sqAbsoluteTolerance)
70     {
71     return vertex;
72     }
73     }
74     return NULL;
75     }
76    
77     void MAILLEUR1D::maille()
78     {
79 foucault 64 if (mcEdge->GetPolyCurve()->get_longueur() == 0)
80     return;
81    
82 foucault 27 int creation_metrique=0;
83     if (metrique==NULL)
84     {
85 foucault 569 printf("Erreur : metrique = NULL\n");
86     return;
87     /* metrique=new FCT_TAILLE_ARETE(epsilon,distance_maximale,mcEdge);
88     creation_metrique=1;*/
89 foucault 27 }
90    
91     /* calcul de la longueur de la arete dans la metrique */
92     double t1=mcEdge->get_cosommet1()->get_t();
93     double t2=mcEdge->get_cosommet2()->get_t();
94     if ((t2<t1+0.0001)&&(mcEdge->get_courbe()->est_periodique())) t2=t2+mcEdge->get_courbe()->get_periode();
95    
96 foucault 569 std::vector<double> tab_integrale_ti, tab_integrale_li;
97 foucault 27 double ti,tii,t;
98     double longueur_metrique=0.;
99     tii=t1;
100     double dt=0.001;
101     do
102     {
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     mcEdge->deriver_seconde(ti,ddxyz,dxyz,xyz);
123     mcEdge->evaluer(tii,xyzii);
124     mcEdge->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 foucault 569
148    
149 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];
150     double residu=fabs(sqrt(facteur3)-sqrt(facteur1)-1./2./sqrt(facteur1)*facteur2*(tii-ti));
151     //if (residu>0.01) dt=dt/2.; else {dt=dt*2;bon_dt=1;}
152     //if (dt<1e-4) {dt=1e-4;bon_dt=1;}
153     //if (dt>0.25) {dt=0.25;bon_dt=1;}
154     bon_dt=1;
155    
156     }
157     while (bon_dt!=1);
158     #endif // ON_A_ABANDONNE_ET_ON_PREND_DT_EGALE_T2_MOINS_T1_DIVISE_PAR_1000
159    
160     t=0.7886751345*ti+0.2113248654*tii;
161     double coord[3];
162     double coefficient_metrique[9];
163     mcEdge->deriver(t,coord);
164 foucault 569
165 foucault 27 if (creation_metrique) metrique->evaluer(&t,coefficient_metrique);
166     else
167     {
168     double xyz[3];
169     mcEdge->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 foucault 569 mcEdge->deriver(t,coord);
176    
177 foucault 27 if (creation_metrique) metrique->evaluer(&t,coefficient_metrique);
178     else
179     {
180     double xyz[3];
181     mcEdge->evaluer(t,xyz);
182     metrique->evaluer(xyz,coefficient_metrique);
183 foucault 569 }
184 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];
185 foucault 569
186 foucault 27 longueur_metrique=longueur_metrique+0.5*(tii-ti)*sqrt(facteur1);
187 foucault 569
188     tab_integrale_ti.push_back(ti);
189     tab_integrale_li.push_back(longueur_metrique);
190 foucault 27 }
191     while (tii<t2);
192     int nombre_de_segment=(int)floor(longueur_metrique);
193     if (longueur_metrique-floor(longueur_metrique)>0.5) nombre_de_segment++;
194     if (nombre_de_segment<1) nombre_de_segment=1;
195    
196     /* discretisation */
197    
198     double valeur_cible=longueur_metrique/nombre_de_segment;
199     int bon_noeud=0;
200     int numnoeud=0;
201 foucault 569 MG_NOEUD* mg_noeud_depart;
202 foucault 27 while (bon_noeud==0)
203     {
204 foucault 569 mg_noeud_depart=(MG_NOEUD*)mcEdge->get_cosommet1()->get_sommet()->get_lien_maillage()->get(numnoeud);
205     MG_NOEUD * mg_noeudtemp=mg_maillage->get_mg_noeudid(mg_noeud_depart->get_id());
206     if (mg_noeudtemp==NULL) numnoeud++; else bon_noeud=1;
207 foucault 27 }
208     bon_noeud=0;
209     numnoeud=0;
210 foucault 569 MG_NOEUD* mg_noeud_arrivee;
211 foucault 27 while (bon_noeud==0)
212     {
213 foucault 569 mg_noeud_arrivee=(MG_NOEUD*)mcEdge->get_cosommet2()->get_sommet()->get_lien_maillage()->get(numnoeud);
214     MG_NOEUD* mg_noeudtemp=mg_maillage->get_mg_noeudid(mg_noeud_arrivee->get_id());
215     if (mg_noeudtemp==NULL) numnoeud++; else bon_noeud=1;
216 foucault 27 }
217 foucault 569 MCNode * noeud_depart = (MCNode*) mg_noeud_depart;
218     MCNode * noeud_arrivee = (MCNode*) mg_noeud_arrivee;
219 foucault 27
220     MCNode* noeud_precedent=noeud_depart;
221    
222     int nb_segment_cree=0;
223 foucault 569 int nb_pas_integrale=tab_integrale_ti.size();
224 foucault 27 for (int i=0;i<nb_pas_integrale;i++)
225     {
226 foucault 569 while ((tab_integrale_li[i]>(nb_segment_cree+1)*valeur_cible) && (nb_segment_cree<nombre_de_segment-1))
227 foucault 27 {
228 foucault 569 double ti=tab_integrale_ti[i];
229 foucault 27 double tii=t2;
230     double li=0.;
231 foucault 569 if (i!=0) li=tab_integrale_li[i-1];
232     double lii=tab_integrale_li[i];
233     if (i!=nb_pas_integrale-1) tii=tab_integrale_ti[i+1];
234 foucault 27 t=ti+(tii-ti)/(lii-li)*((nb_segment_cree+1)*valeur_cible-li);
235     double coo[3];
236     mcEdge->evaluer(t,coo);
237    
238     MG_ARETE * nouveau_noeud_refEdge = NULL;
239     // parametre de l'arete de reference
240     double nouveau_noeud_refEdgeT = 0, Parameter_SToRefEdgeDerT;
241     mcEdge->GetPolyCurve()->Parameter_SToRefEdgeT(t,&nouveau_noeud_refEdge, &nouveau_noeud_refEdgeT, &Parameter_SToRefEdgeDerT, false);
242     // reference topology mapping
243     MCNode * nouveau_noeud = NULL;
244     MG_SOMMET * refVertex = IsPointInRefVertex(mcEdge->GetPolyCurve(),coo,1E-6);
245     if (refVertex == NULL)
246     {
247     nouveau_noeud = new MCNode(mcEdge,nouveau_noeud_refEdge,nouveau_noeud_refEdgeT,coo);
248     }
249     else
250     {
251     nouveau_noeud = new MCNode(mcEdge,refVertex,coo[0],coo[1],coo[2]);
252     }
253 foucault 569 mg_maillage->ajouter_mg_noeud(nouveau_noeud);
254     nouveau_noeud->Creator = MCNode::mailleur_1d;
255 foucault 27 MCSegment * nouveau_segment = new MCSegment (mcEdge,noeud_precedent,nouveau_noeud);
256    
257     mg_maillage->ajouter_mg_segment(nouveau_segment);
258     noeud_precedent=nouveau_noeud;
259     nb_segment_cree++;
260     }
261     }
262    
263     MCSegment * nouveau_segment = new MCSegment (mcEdge,noeud_precedent,noeud_arrivee);
264    
265     mg_maillage->ajouter_mg_segment(nouveau_segment);
266    
267     if (debug)
268     {
269     char fichierProgression1D [1024];
270     sprintf(fichierProgression1D, "%s%s", "c:\\temp\\void", "_1D.mai");
271     ofstream o3(fichierProgression1D,ios::out|ios::trunc);
272     o3.precision(16);
273     o3.setf(ios::showpoint);
274     mg_maillage->enregistrer_sous_mesh_1D(o3);
275     }
276    
277     if (creation_metrique==1)
278     {
279     delete metrique;
280     metrique=NULL;
281     }
282    
283     }
284    
285 foucault 569 void MAILLEUR1D::maille_points(std::vector <OT_VECTEUR_3D> & __points)
286     {
287     if (mcEdge->GetPolyCurve()->get_longueur() == 0)
288     return;
289 foucault 27
290 foucault 569 /* calcul de la longueur de la arete dans la metrique */
291     double t1=mcEdge->get_cosommet1()->get_t();
292     double t2=mcEdge->get_cosommet2()->get_t();
293     if ((t2<t1+0.0001)&&(mcEdge->get_courbe()->est_periodique())) t2=t2+mcEdge->get_courbe()->get_periode();
294 foucault 27
295 foucault 569 std::set <double> tab_t;
296     std::set <double>::iterator tab_t_it1, tab_t_it2, tab_t_it;
297 foucault 27
298 foucault 569 tab_t.insert(mcEdge->get_cosommet1()->get_t());
299     tab_t.insert(mcEdge->get_cosommet2()->get_t());
300 foucault 27
301 foucault 569 for (int i=0; i<__points.size(); i++)
302     {
303     double t;
304     mcEdge->inverser(t, __points[i]);
305     OT_VECTEUR_3D xyz;
306     mcEdge->evaluer(t, xyz);
307     double distance = (xyz-__points[i]).get_longueur();
308     tab_t_it1 = tab_t.lower_bound(t);
309     tab_t_it2 = tab_t.upper_bound(t);
310     double t_dist = 1E99;
311     for (tab_t_it = tab_t_it1; tab_t_it != tab_t.end(); tab_t_it++)
312     {
313     t_dist = std::min(t_dist, fabs(t-*tab_t_it));
314     if (tab_t_it == tab_t_it2) break;
315     }
316     if (distance > 0.001 * mcEdge->GetPolyCurve()->get_longueur())
317     {
318     printf("node #%d/%lu is too far from edge %lu\n", i, __points.size(), mcEdge->get_id());
319     continue;
320     }
321     if (t_dist < 0.001 * t2-t1)
322     {
323     printf("node #%d/%lu already exists in edge %lu\n", i, __points.size(), mcEdge->get_id());
324     continue;
325     }
326 foucault 27
327 foucault 569 tab_t.insert(t);
328     }
329 foucault 27
330 foucault 569 /* discretisation */
331     int bon_noeud=0;
332     int numnoeud=0;
333     MG_NOEUD* mg_noeud_depart;
334     while (bon_noeud==0)
335     {
336     mg_noeud_depart=(MG_NOEUD*)mcEdge->get_cosommet1()->get_sommet()->get_lien_maillage()->get(numnoeud);
337     MG_NOEUD * mg_noeudtemp=mg_maillage->get_mg_noeudid(mg_noeud_depart->get_id());
338     if (mg_noeudtemp==NULL) numnoeud++; else bon_noeud=1;
339     }
340     bon_noeud=0;
341     numnoeud=0;
342     MG_NOEUD* mg_noeud_arrivee;
343     while (bon_noeud==0)
344     {
345     mg_noeud_arrivee=(MG_NOEUD*)mcEdge->get_cosommet2()->get_sommet()->get_lien_maillage()->get(numnoeud);
346     MG_NOEUD* mg_noeudtemp=mg_maillage->get_mg_noeudid(mg_noeud_arrivee->get_id());
347     if (mg_noeudtemp==NULL) numnoeud++; else bon_noeud=1;
348     }
349     MCNode * noeud_depart = (MCNode*) mg_noeud_depart;
350     MCNode * noeud_arrivee = (MCNode*) mg_noeud_arrivee;
351    
352     MCNode* noeud_precedent=noeud_depart;
353    
354     int nb_segment_cree=0;
355     tab_t_it2 = tab_t.begin();
356     for (tab_t_it1 = tab_t.begin(); nb_segment_cree+2 < tab_t.size(); tab_t_it1++)
357     {
358     tab_t_it2++;
359     if (tab_t_it2 == tab_t.end()) break;
360    
361     double ti=*tab_t_it1;
362     double tii=*tab_t_it2;
363    
364     MG_ARETE * nouveau_noeud_refEdge = NULL;
365     // parametre de l'arete de reference
366     double nouveau_noeud_refEdgeT = 0, Parameter_SToRefEdgeDerT;
367     mcEdge->GetPolyCurve()->Parameter_SToRefEdgeT(tii,&nouveau_noeud_refEdge, &nouveau_noeud_refEdgeT, &Parameter_SToRefEdgeDerT, false);
368     // reference topology mapping
369     MCNode * nouveau_noeud = NULL;
370     double coo[3]; mcEdge->evaluer(tii, coo);
371     MG_SOMMET * refVertex = IsPointInRefVertex(mcEdge->GetPolyCurve(),coo,1E-6);
372     if (refVertex == NULL)
373     {
374     nouveau_noeud = new MCNode(mcEdge,nouveau_noeud_refEdge,nouveau_noeud_refEdgeT,coo);
375     }
376     else
377     {
378     nouveau_noeud = new MCNode(mcEdge,refVertex,coo[0],coo[1],coo[2]);
379     }
380     mg_maillage->ajouter_mg_noeud(nouveau_noeud);
381     nouveau_noeud->Creator = MCNode::mailleur_1d;
382     MCSegment * nouveau_segment = new MCSegment (mcEdge,noeud_precedent,nouveau_noeud);
383    
384     mg_maillage->ajouter_mg_segment(nouveau_segment);
385     noeud_precedent=nouveau_noeud;
386     nb_segment_cree++;
387     }
388    
389     MCSegment * nouveau_segment = new MCSegment (mcEdge,noeud_precedent,noeud_arrivee);
390    
391     mg_maillage->ajouter_mg_segment(nouveau_segment);
392    
393     if (debug)
394     {
395     char fichierProgression1D [1024];
396     sprintf(fichierProgression1D, "%s%s", "c:\\temp\\void", "_1D.mai");
397     ofstream o3(fichierProgression1D,ios::out|ios::trunc);
398     o3.precision(16);
399     o3.setf(ios::showpoint);
400     mg_maillage->enregistrer_sous_mesh_1D(o3);
401     }
402     }