MAGiC  V5.0
Mailleurs Automatiques de Géometries intégrés à la Cao
ot_fonctions.cpp
Aller à la documentation de ce fichier.
1 //####//------------------------------------------------------------
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 //####// ot_fonctions.cpp
15 //####//
16 //####//------------------------------------------------------------
17 //####//------------------------------------------------------------
18 //####// COPYRIGHT 2000-2025
19 //####// Derniere modification par francois
20 //####// mer 14 mai 2025 17:54:47 EDT
21 //####//------------------------------------------------------------
22 //####//------------------------------------------------------------
23 
24 
25 #include "ot_fonctions.h"
26 
27 
28 double2 moment (double2 theta,std::vector<OT_VECTEUR_4DD>& ctrpts,int plan)
29 {
30 
31 
32  int nb_pts=ctrpts.size();
33  double ordre=4.;
34  double2 zer=0.;
35  double2 sum=0.;
36  for (int i=0;i<nb_pts;i++)
37  {
38  //thi[i]=thi[i]*PI/180;
39  //double2 nori=ctrpts[i].norme();
40  double2 x=ctrpts[i].x()^2;
41  double2 y=ctrpts[i].y()^2;
42  double2 z=ctrpts[i].z()^2;
43  double2 nori;
44 
45  // double2 nori=ctrpts[i].norme();
46  double2 cosi=0.;
47  double2 sini=0.;
48  if (plan==12)
49  {
50  nori=(x+y)^0.5;
51  if (nori!=zer) {
52  cosi=ctrpts[i].get_x()/nori;
53  sini=ctrpts[i].get_y()/nori;
54  }
55  }
56  if (plan==23)
57  {
58  nori=(y+z)^0.5;
59  if (nori!=zer) {
60  cosi=ctrpts[i].get_y()/nori;
61  sini=ctrpts[i].get_z()/nori;
62  }
63 
64  }
65  if (plan==13)
66  {
67  nori=(x+z)^0.5;
68  if (nori!=zer) {
69  cosi=ctrpts[i].get_x()/nori;
70  sini=ctrpts[i].get_z()/nori;
71  }
72  }
73  double2 costh=cos(theta);
74  double2 sinth=sin(theta);
75  double2 d_theta=sini*costh-cosi*sinth;
76  //sum=sum+ pow(li[i],4)*pow(d_theta,4);
77  double2 val1=nori^ordre;
78  double2 val2=d_theta^ordre;
79  double2 val=val1*val2;
80  sum=sum+ val;
81  }
82  return sum;
83 }
84 
85 
86 double2 Max(double2 a,double2 b)
87 {
88  return((a<b)? b : a);
89 }
90 
91 double2 Min(double2 a,double2 b)
92 {
93  return ((a>b)? b : a);
94 }
95 
96 double2 Sign(double2 a,double2 b)
97 {
98  double2 zer=0.;
99  double2 c=-1.;
100  if (b>=zer)
101  return a;
102  else
103  return c*a;
104 }
105 
106 void swap(double2& a,double2& b)
107 {
108  double2 tmp;
109  tmp = a;
110  a = b;
111  b = tmp;
112 }
113 void Shift(double2&a,double2&b,double2&c,double2 d)
114 {
115  a=b;
116  b=c;
117  c=d;
118 }
119 
120 void min_Brak(double2 &ax,double2 &bx,double2 &cx,double2 &fa,double2 &fb,double2 &fc,std::vector<OT_VECTEUR_4DD>& ctrpts,int plan)
121 {
122  static const double2 GOLD=1.618034;
123  static const double2 GLIMIT=100.0;
124  static const double2 TINY=1.0e-20;
125  double2 ulim,u,r,q,fu;
126  double2(*func)(double2,std::vector<OT_VECTEUR_4DD>&,int)=&moment;
127  double2 dtheta=0.0001;
128  ax=0.;
129  bx=ax+dtheta;
130  cx=bx+dtheta;
131  double2 fva=(*func)(ax,ctrpts,plan);
132  double2 fvb=(*func)(bx,ctrpts,plan);
133  double2 fvc=(*func)(cx,ctrpts,plan);
134 
135  while (fvb> fva)
136  {
137  ax=bx;
138  bx=ax+dtheta;
139  fva=(*func)(ax,ctrpts,plan);
140  fvb=(*func)(bx,ctrpts,plan);
141  }
142 
143  if (fvb<fva)
144  {
145  cx=bx+dtheta;
146  while (fvb> fvc)
147  {
148  bx=cx;
149  cx=bx+dtheta;
150  fvb=(*func)(bx,ctrpts,plan);
151  fvc=(*func)(cx,ctrpts,plan);
152  }
153  }
154 
155  fa=(*func)(ax,ctrpts,plan);
156  fb=(*func)(bx,ctrpts,plan);
157  if (fb>fa)
158  {
159  swap(ax,bx);
160  swap(fb,fa);
161  }
162  cx=bx+GOLD*(bx-ax);
163  fc=(*func)(cx,ctrpts,plan);
164  double2 d=2.;
165  double2 zer=0.;
166  while (fb>fc)
167  {
168  r=(bx-ax)*(fb-fc);
169  q=(bx-cx)*(fb-fa);
170  u=bx-((bx-cx)*q-(bx-ax)*r)/(d*Sign(Max((q-r).get_fabs(),TINY),q-r));//(d*Sign(Max(abs(q-r),TINY),q-r));
171  ulim=bx+GLIMIT*(cx-bx);
172  if ((bx-u)*(u-cx)>zer)
173  {
174  fu=(*func)(u,ctrpts,plan);
175  if (fu<fc)
176  {
177  ax=bx;
178  bx=u;
179  fa=fb;
180  fb=fu;
181  return;
182  }
183  else if (fu>fb)
184  {
185  cx=u;
186  fc=fu;
187  return;
188  }
189  u=cx+GOLD*(cx-bx);
190  fu=(*func)(u,ctrpts,plan);
191  }
192  else if ((cx-u)*(u-ulim)>zer)
193  {
194  fu=(*func)(u,ctrpts,plan);
195  if (fu<fc)
196  {
197  bx=cx;
198  cx=u;
199  u=cx+GOLD*(cx-bx);
200  Shift(fb,fc,fu,(*func)(u,ctrpts,plan));
201  }
202  }
203  else if ((u-ulim)*(ulim-cx)>=zer)
204  {
205  u=ulim;
206  fu=(*func)(u,ctrpts,plan);
207  }
208  else
209  {
210  u=cx+GOLD*(cx-bx);
211  fu=(*func)(u,ctrpts,plan);
212  }
213  Shift(ax,bx,cx,u);
214  Shift(fa,fb,fc,fu);
215  }
216 }
217 
218 
219 
220 
221 std::vector<OT_VECTEUR_4DD> system_axes(std::vector<OT_VECTEUR_4DD>& ctrpts)
222 {
223  std::vector<OT_VECTEUR_4DD> syst_axe;
224  int nb_points=ctrpts.size();
225  std::vector<OT_VECTEUR_4DD> ctrpts_proj;
226  OT_VECTEUR_4DD axe1,axe2,axe3,axe4;
227  double2 ax,bx,cx,fa,fb,fc;
228  int plan=12;
229  double2 z0=ctrpts[0].get_z() ;
230  double2 y0=ctrpts[0].get_z() ;
231  double2 x0=ctrpts[0].get_z() ;
232  int cmpt1=0;
233  int cmpt2=0;
234  int cmpt3=0;
235  for (int i=0;i<nb_points;i++)
236  {
237  if (ctrpts[i].get_x()==x0) cmpt1++;
238  if (ctrpts[i].get_y()==y0) cmpt2++;
239  if (ctrpts[i].get_z()==z0) cmpt3++;
240  }
241 
242  if (cmpt1==nb_points) plan=23;
243  if (cmpt2==nb_points) plan=13;
244  if (cmpt3==nb_points) plan=12;
245 
246 
247 
248  min_Brak(ax,bx,cx,fa,fb,fc,ctrpts,plan);
249  double theta= bx.get_x();
250  double c=cos(theta);
251  double s=sin(theta);
252  axe1.change_x(c);
253  axe1.change_y(s);
254  axe1.change_z(0.);
255  axe1.change_w(0.);
256  double e1_proj[3];
257  double e2_proj[3];
258  double normal[3];
259  normal[0] =(axe1.get_x()).get_x();
260  normal[1] =(axe1.get_y()).get_x();
261  normal[2] =(axe1.get_z()).get_x();
262  double root[3]={0.,0.,0.};
263  OT_ALGORITHME_GEOMETRIQUE algo;
264  int pris=0;
265  OT_MATRICE_3D P, P_1;
266 
267  for (int i=0;i<nb_points;i++)
268  {
269  double pnt[3];
270  double proj_pnt[3];
271  pnt[0] =(ctrpts[i].get_x()).get_x();
272  pnt[1] =(ctrpts[i].get_y()).get_x();
273  pnt[2] =(ctrpts[i].get_z()).get_x();
274 
275 
276  algo.Proj3D_Point_Plan (normal, root, pnt, proj_pnt);
277 
278  OT_VECTEUR_4DD pt_proj;
279 
280 
281  if ((proj_pnt[0]!=0.||proj_pnt[1]!=0.||proj_pnt[2]!=0)&&!pris)
282  {
283  e1_proj[0]=proj_pnt[0];
284  e1_proj[1]=proj_pnt[1];
285  e1_proj[2]=proj_pnt[2];
286  double norm=sqrt(pow(proj_pnt[0],2)+pow(proj_pnt[1],2)+pow(proj_pnt[2],2));
287 
288  e1_proj[0]=e1_proj[0]/norm;
289  e1_proj[1]=e1_proj[1]/norm;
290  e1_proj[2]=e1_proj[2]/norm;
291 
292  e2_proj[0]=normal[1]*e1_proj[2]-normal[2]*e1_proj[1];
293  e2_proj[1]=normal[2]*e1_proj[0]-normal[2]*e1_proj[2];
294  e2_proj[2]=normal[0]*e1_proj[1]-normal[2]*e1_proj[0];
295 
296  OT_VECTEUR_3D V1(e1_proj[0],e2_proj[0],normal[0]);
297  OT_VECTEUR_3D V2(e1_proj[1],e2_proj[1],normal[1]);
298  OT_VECTEUR_3D V3(e1_proj[2],e2_proj[2],normal[2]);
299 
300 
301 
302  P.change_vecteur1(V1);
303  P.change_vecteur2(V2);
304  P.change_vecteur3(V3);
305 
306  P_1=P.inverse();
307  pris=1;
308  }
309 
310  OT_VECTEUR_3D V_PT(proj_pnt[0],proj_pnt[1],proj_pnt[2]);
311 
312  OT_VECTEUR_3D V=P_1*V_PT;
313 
314 
315  pt_proj.change_x(V[0]);
316  pt_proj.change_y(V[1]);
317  pt_proj.change_z(V[2]);
318 
319  ctrpts_proj.insert(ctrpts_proj.end(),pt_proj);
320  }
321 
322  min_Brak(ax,bx,cx,fa,fb,fc,ctrpts_proj,plan);
323 
324  theta= bx.get_x();
325  c=cos(theta);
326  s=sin(theta);
327  OT_VECTEUR_3D Vp(c,s,0.);
328  double nnv_x=Vp*P.get_vecteur1();
329  double nnv_y=Vp*P.get_vecteur2();
330  double nnv_z=Vp*P.get_vecteur3();
331 
332  axe2.change_x(nnv_x);
333  axe2.change_y(nnv_y);
334  axe2.change_z(nnv_z);
335 
336  double2 val=axe1*axe2;
337 
338  axe3[0]=axe1[1]*axe2[2]-axe1[2]*axe2[1];
339  axe3[1]=axe1[2]*axe2[0]-axe1[2]*axe2[2];
340  axe3[2]=axe1[0]*axe2[1]-axe1[2]*axe2[0];
341 
342  axe3[3]=0.;
343 
344 
345  syst_axe.insert( syst_axe.end(),axe1);
346  syst_axe.insert( syst_axe.end(),axe2);
347  syst_axe.insert( syst_axe.end(),axe3);
348 
349  return syst_axe;
350 }
351 
352 
353 OT_FONCTION_SYMBOLIQUE::OT_FONCTION_SYMBOLIQUE(char* f, char *var):formule(f)
354 {
355 parseur.SetExpr(formule.c_str());
356 char *p=var;
357 bool debut=true;
358 for (int i=0;i<strlen(var)+1;i++)
359  {
360  if (debut)
361  if ((var[i]!=' ') && (var[i]!=0))
362  debut=false;
363  if (!debut)
364  if ((var[i]==' ') || (var[i]==0))
365  {
366  char v[200];
367  sscanf(p,"%s",v);
368  variable.push_back(v);
369  p=var+i;
370  debut=true;
371  }
372 
373  }
374  nb_variable=variable.size();
375 }
376 
377 
378 
379 OT_FONCTION_SYMBOLIQUE::OT_FONCTION_SYMBOLIQUE(OT_FONCTION_SYMBOLIQUE& mdd):nb_variable(mdd.nb_variable),formule(mdd.formule)
380 {
381 }
382 
383 
384 
385 OT_FONCTION_SYMBOLIQUE::~OT_FONCTION_SYMBOLIQUE()
386 {
387 }
388 
389 int OT_FONCTION_SYMBOLIQUE::get_nb_variable()
390 {
391 return nb_variable;
392 }
393 
394 std::string OT_FONCTION_SYMBOLIQUE::get_nom_variable(int num)
395 {
396 return variable[num];
397 }
398 
399 
400 void OT_FONCTION_SYMBOLIQUE::def_fonction ( char* nom, double f(double,double))
401 {
402 parseur.DefineFun(nom,f);
403 }
404 
405 void OT_FONCTION_SYMBOLIQUE::def_fonction ( char* nom, double f( double ) )
406 {
407 parseur.DefineFun(nom,f);
408 }
409 
410 void OT_FONCTION_SYMBOLIQUE::def_fonction ( char* nom, double f( double, double, double ))
411 {
412 parseur.DefineFun(nom,f);
413 }
414 
415 double OT_FONCTION_SYMBOLIQUE::evalue ( double* var )
416 {
417 for (int i=0;i<nb_variable;i++)
418  parseur.DefineVar(variable[i], var+i);
419 double result;
420 try
421  {
422  result = parseur.Eval();
423  }
424 catch (mu::Parser::exception_type &e)
425  {
426  result=0.;
427  }
428 return result;
429 }
430 
431 double OT_FONCTION_SYMBOLIQUE::integre ( double a, double b, int nb_pas)
432 {
433 if (nb_variable!=1) return 0.;
434 
435 double intg=0.;
436 for (int i=0;i<nb_pas;i++)
437  {
438  double ti=a+(b-a)/nb_pas*i;
439  double tii=a+(b-a)/nb_pas*(i+1.);
440  double xsi=-1./sqrt(3);
441  double t=(tii-ti)/2.*xsi+(tii+ti)/2.;
442  double res=evalue(&t);
443  intg=intg+(tii-ti)/2.*(1.*res);
444  xsi=1./sqrt(3);
445  t=(tii-ti)/2.*xsi+(tii+ti)/2.;
446  res=evalue(&t);
447  intg=intg+(tii-ti)/2.*(1.*res);
448  }
449 return intg;
450 }
451 
OT_VECTEUR_4DD::change_x
virtual void change_x(double2 x)
Definition: ot_mathematique.cpp:1456
double2::x
double x
Definition: ot_doubleprecision.h:84
Sign
double2 Sign(double2 a, double2 b)
Definition: ot_fonctions.cpp:96
OT_VECTEUR_4DD::change_y
virtual void change_y(double2 y)
Definition: ot_mathematique.cpp:1460
OT_FONCTION_SYMBOLIQUE::evalue
virtual double evalue(double *var)
Definition: ot_fonctions.cpp:415
OT_VECTEUR_4DD::change_w
virtual void change_w(double2 w)
Definition: ot_mathematique.cpp:1468
a
#define a(i, j)
OT_FONCTION_SYMBOLIQUE::~OT_FONCTION_SYMBOLIQUE
virtual ~OT_FONCTION_SYMBOLIQUE()
Definition: ot_fonctions.cpp:385
ot_fonctions.h
double2::get_x
double get_x()
Definition: ot_doubleprecision.cpp:367
OT_FONCTION_SYMBOLIQUE::get_nb_variable
virtual int get_nb_variable(void)
Definition: ot_fonctions.cpp:389
OT_FONCTION_SYMBOLIQUE::formule
std::string formule
Definition: ot_fonctions.h:65
OT_FONCTION_SYMBOLIQUE::get_nom_variable
virtual std::string get_nom_variable(int num)
Definition: ot_fonctions.cpp:394
OT_ALGORITHME_GEOMETRIQUE::Proj3D_Point_Plan
static void Proj3D_Point_Plan(double *norm, double *root, double *pnt, double *proj_pnt)
Definition: ot_algorithme_geometrique.cpp:309
system_axes
std::vector< OT_VECTEUR_4DD > system_axes(std::vector< OT_VECTEUR_4DD > &ctrpts)
Definition: ot_fonctions.cpp:221
swap
void swap(double2 &a, double2 &b)
Definition: ot_fonctions.cpp:106
V2
bool V2(MCBody *_mcBody, MG_ELEMENT_TOPOLOGIQUE *topo)
Definition: CAD4FE_MCBody.cpp:804
OT_FONCTION_SYMBOLIQUE
Definition: ot_fonctions.h:46
min_Brak
void min_Brak(double2 &ax, double2 &bx, double2 &cx, double2 &fa, double2 &fb, double2 &fc, std::vector< OT_VECTEUR_4DD > &ctrpts, int plan)
Definition: ot_fonctions.cpp:120
OT_FONCTION_SYMBOLIQUE::variable
std::vector< std::string > variable
Definition: ot_fonctions.h:66
double2
Definition: ot_doubleprecision.h:29
Max
double2 Max(double2 a, double2 b)
Definition: ot_fonctions.cpp:86
OT_FONCTION_SYMBOLIQUE::def_fonction
virtual void def_fonction(char *nom, double f(double))
Definition: ot_fonctions.cpp:405
f
double f(double x, long nb, double *xfonc, double *fonc, double eng, double eni, double lambda, double nor, double *fonc2)
Definition: fct_generateur_calibrage.cpp:96
OT_VECTEUR_4DD::get_y
virtual double2 get_y(void) const
Definition: ot_mathematique.cpp:1432
OT_VECTEUR_4DD::get_z
virtual double2 get_z(void) const
Definition: ot_mathematique.cpp:1440
OT_FONCTION_SYMBOLIQUE::nb_variable
int nb_variable
Definition: ot_fonctions.h:64
OT_MATRICE_3D
Definition: ot_mathematique.h:160
V
void V(MCAA *mcaa)
Definition: CAD4FE_MCAA.cpp:1794
OT_ALGORITHME_GEOMETRIQUE
Definition: ot_algorithme_geometrique.h:32
OT_VECTEUR_4DD::get_x
virtual double2 get_x(void) const
Definition: ot_mathematique.cpp:1424
OT_FONCTION_SYMBOLIQUE::integre
virtual double integre(double a, double b, int nb_pas=32)
Definition: ot_fonctions.cpp:431
OT_VECTEUR_3D
Definition: ot_mathematique.h:94
OT_VECTEUR_4DD::change_z
virtual void change_z(double2 z)
Definition: ot_mathematique.cpp:1464
sqrt
double2 sqrt(double2 &val)
Definition: ot_doubleprecision.cpp:345
OT_FONCTION_SYMBOLIQUE::parseur
mu::Parser parseur
Definition: ot_fonctions.h:63
moment
double2 moment(double2 theta, std::vector< OT_VECTEUR_4DD > &ctrpts, int plan)
Definition: ot_fonctions.cpp:28
res
#define res(i, j)
OT_FONCTION_SYMBOLIQUE::OT_FONCTION_SYMBOLIQUE
OT_FONCTION_SYMBOLIQUE(char *f, char *var)
Definition: ot_fonctions.cpp:353
Shift
void Shift(double2 &a, double2 &b, double2 &c, double2 d)
Definition: ot_fonctions.cpp:113
cos
double2 cos(double2 &val)
Definition: ot_doubleprecision.cpp:206
OT_VECTEUR_4DD
Definition: ot_mathematique.h:284
P
#define P(i, j)
Min
double2 Min(double2 a, double2 b)
Definition: ot_fonctions.cpp:91
sin
double2 sin(double2 &val)
Definition: ot_doubleprecision.cpp:250