MAGiC  V5.0
Mailleurs Automatiques de Géometries intégrés à la Cao
fct_taille_metrique.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 //####// fct_taille_metrique.cpp
15 //####//
16 //####//------------------------------------------------------------
17 //####//------------------------------------------------------------
18 //####// COPYRIGHT 2000-2024
19 //####// jeu 13 jun 2024 11:58:53 EDT
20 //####//------------------------------------------------------------
21 //####//------------------------------------------------------------
22 #include "gestionversion.h"
23 #include "fct_taille_metrique.h"
24 #include <math.h>
25 #include <map>
26 
27 
28 
29 
30 
31 void FCT_TAILLE_METRIQUE::evaluer_decompose(double* metrique_depart, double* valeur_propre, double* vecteur_propre)
32 {
33 double a00 = metrique_depart[0];
34 double a01 = metrique_depart[1];
35 double a02 = metrique_depart[2];
36 double a11 = metrique_depart[4];
37 double a12 = metrique_depart[5];
38 double a22 = metrique_depart[8];
39 double inv3=1./3.;
40 double c0 = a00*a11*a22 + 2.0*a01*a02*a12 - a00*a12*a12 - a11*a02*a02 - a22*a01*a01;
41 double c1 = a00*a11 - a01*a01 + a00*a22 - a02*a02 + a11*a22 - a12*a12;
42 double c2 = a00 + a11 + a22;
43 double c2Div3 = c2*inv3;
44 double aDiv3 = (c1 - c2*c2Div3)*inv3;
45 if (aDiv3 > 0.0) { aDiv3 = 0.0; }
46 double mbDiv2 = 0.5*(c0 + c2Div3*(2.0*c2Div3*c2Div3 - c1));
47 double q = mbDiv2*mbDiv2 + aDiv3*aDiv3*aDiv3;
48 if (q > 0.0) { q = 0.0; }
49 double magnitude = sqrt(-aDiv3);
50 double angle = atan2(sqrt(-q),mbDiv2)*inv3;
51 double cs = cos(angle);
52 double sn = sin(angle);
53 valeur_propre[0] = c2Div3 + 2.0*magnitude*cs;
54 valeur_propre[1] = c2Div3 - magnitude*(cs + sqrt(3)*sn);
55 valeur_propre[2] = c2Div3 - magnitude*(cs - sqrt(3)*sn);
56 std::multimap<double,double> tab;
57 tab.insert(std::pair<double,double>(valeur_propre[0],valeur_propre[0]));
58 tab.insert(std::pair<double,double>(valeur_propre[1],valeur_propre[1]));
59 tab.insert(std::pair<double,double>(valeur_propre[2],valeur_propre[2]));
60 std::map<double,double>::iterator it=tab.begin();
61 valeur_propre[0]=it->second;it++;
62 valeur_propre[1]=it->second;it++;
63 valeur_propre[2]=it->second;
64 double metrique2[9];
65 metrique2[0]=metrique_depart[0]-valeur_propre[0];
66 metrique2[4]=metrique_depart[4]-valeur_propre[0];
67 metrique2[8]=metrique_depart[8]-valeur_propre[0];
68 int rank0 = ComputeRank(metrique2);
69 if (rank0 == 0)
70  {
71  vecteur_propre[0] =1.;
72  vecteur_propre[1] =0.;
73  vecteur_propre[2] =0.;
74  vecteur_propre[3] =0.;
75  vecteur_propre[4] =1.;
76  vecteur_propre[5] =0.;
77  vecteur_propre[6] =0.;
78  vecteur_propre[7] =0.;
79  vecteur_propre[8] =1.;
80  return;
81  }
82 OT_VECTEUR_3D v1,v2,v3;
83 if (rank0 == 1)
84  {
85  OT_VECTEUR_3D row0(metrique2[0],metrique2[3],metrique2[6]);
86  GetComplement2(row0,v1,v2);
87  v3=v1&v2;
88  vecteur_propre[0] =v1.get_x();
89  vecteur_propre[1] =v2.get_x();
90  vecteur_propre[2] =v3.get_x();
91  vecteur_propre[3] =v1.get_y();
92  vecteur_propre[4] =v2.get_y();
93  vecteur_propre[5] =v3.get_z();
94  vecteur_propre[6] =v1.get_x();
95  vecteur_propre[7] =v2.get_y();
96  vecteur_propre[8] =v3.get_z();
97  return;
98  }
99 OT_VECTEUR_3D row10(metrique2[0],metrique2[3],metrique2[6]);
100 OT_VECTEUR_3D row11(metrique2[1],metrique2[4],metrique2[7]);
101 GetComplement1(row10,row11,v1);
102 metrique2[0]=metrique_depart[0]-valeur_propre[1];
103 metrique2[4]=metrique_depart[4]-valeur_propre[1];
104 metrique2[8]=metrique_depart[8]-valeur_propre[1];
105 int rank1 = ComputeRank(metrique2); // zero rank detected earlier, rank1 must be positive
106 if (rank1 == 1)
107  {
108  GetComplement2(v1,v2,v3);
109  vecteur_propre[0] =v1.get_x();
110  vecteur_propre[1] =v2.get_x();
111  vecteur_propre[2] =v3.get_x();
112  vecteur_propre[3] =v1.get_y();
113  vecteur_propre[4] =v2.get_y();
114  vecteur_propre[5] =v3.get_z();
115  vecteur_propre[6] =v1.get_x();
116  vecteur_propre[7] =v2.get_y();
117  vecteur_propre[8] =v3.get_z();
118  return;
119  }
120 OT_VECTEUR_3D row20(metrique2[0],metrique2[3],metrique2[6]);
121 OT_VECTEUR_3D row21(metrique2[1],metrique2[4],metrique2[7]);
122 GetComplement1(row20,row21,v2);
123 v3=v1&v2;
124 vecteur_propre[0] =v1.get_x();
125 vecteur_propre[1] =v2.get_x();
126 vecteur_propre[2] =v3.get_x();
127 vecteur_propre[3] =v1.get_y();
128 vecteur_propre[4] =v2.get_y();
129 vecteur_propre[5] =v3.get_z();
130 vecteur_propre[6] =v1.get_x();
131 vecteur_propre[7] =v2.get_y();
132 vecteur_propre[8] =v3.get_z();
133 }
134 
135 void FCT_TAILLE_METRIQUE::GetComplement1(OT_VECTEUR_3D U,OT_VECTEUR_3D V, OT_VECTEUR_3D &W)
136 {
137 W = U&V;
138 W.norme();
139 }
140 
141 void FCT_TAILLE_METRIQUE::GetComplement2(OT_VECTEUR_3D U, OT_VECTEUR_3D &V, OT_VECTEUR_3D &W)
142 {
143 U.norme();
144 
145 if (fabs(U.get_x()) >= fabs(U.get_y()))
146  {
147  double invLength = 1/sqrt(U.get_x()*U.get_x() + U.get_z()*U.get_z());
148  V.change_x(-U.get_z()*invLength);
149  V.change_y(0.);
150  V.change_z(-U.get_x()*invLength);
151  W.change_x(U.get_y()*V.get_z());
152  W.change_y(U.get_z()*V.get_x()-U.get_x()*V.get_z());
153  W.change_z(-U.get_y()*V.get_x());
154  }
155 else
156  {
157  double invLength = 1/sqrt(U.get_y()*U.get_y() + U.get_z()*U.get_z());
158  V.change_x(0.);
159  V.change_y(U.get_z()*invLength);
160  V.change_z(-U.get_y()*invLength);
161  W.change_x(U.get_y()*V.get_z()-U.get_z()*V.get_y());
162  W.change_y(-U.get_x()*V.get_z());
163  W.change_z(U.get_x()*V.get_y());
164  }
165 }
166 
167 
168 int FCT_TAILLE_METRIQUE::ComputeRank(double* metrique)
169 {
170 double M[3][3];
171 M[0][0]=metrique[0];
172 M[0][1]=metrique[1];
173 M[0][2]=metrique[2];
174 M[1][0]=metrique[3];
175 M[1][1]=metrique[4];
176 M[1][2]=metrique[5];
177 M[2][0]=metrique[6];
178 M[2][1]=metrique[7];
179 M[2][2]=metrique[8];
180 double abs, save, max = -1;
181 int maxRow = -1, maxCol = -1;
182 for (int row = 0; row < 3; row++)
183  {
184  for (int col = row; col < 3; col++)
185  {
186  abs = fabs(M[row][col]);
187  if (abs > max)
188  {
189  max = abs;
190  maxRow = row;
191  maxCol = col;
192  }
193  }
194  }
195 if (max < 1e-10)
196 {
197 return 0;
198 }
199 if (maxRow != 0)
200  {
201  for (int col = 0; col < 3; col++)
202  {
203  save = M[0][col];
204  M[0][col] = M[maxRow][col];
205  M[maxRow][col] = save;
206  }
207  }
208 double invMax = 1./M[0][maxCol];
209 M[0][0] *= invMax;
210 M[0][1] *= invMax;
211 M[0][2] *= invMax;
212 if (maxCol == 0)
213  {
214  M[1][1] -= M[1][0]*M[0][1];
215  M[1][2] -= M[1][0]*M[0][2];
216  M[2][1] -= M[2][0]*M[0][1];
217  M[2][2] -= M[2][0]*M[0][2];
218  M[1][0] = 0;
219  M[2][0] = 0;
220  }
221 else if (maxCol == 1)
222  {
223  M[1][0] -= M[1][1]*M[0][0];
224  M[1][2] -= M[1][1]*M[0][2];
225  M[2][0] -= M[2][1]*M[0][0];
226  M[2][2] -= M[2][1]*M[0][2];
227  M[1][1] = 0;
228  M[2][1] = 0;
229  }
230 else
231  {
232  M[1][0]-= M[1][2]*M[0][0];
233  M[1][1]-= M[1][2]*M[0][1];
234  M[2][0]-= M[2][2]*M[0][0];
235  M[2][1]-= M[2][2]*M[0][1];
236  M[1][2]= 0;
237  M[2][2]= 0;
238  }
239 max = -1;
240 maxRow = -1;
241 maxCol = -1;
242 for (int row = 1; row < 3; row++)
243  {
244  for (int col = 0; col < 3; col++)
245  {
246  abs = fabs(M[row][col]);
247  if (abs > max)
248  {
249  max = abs;
250  maxRow = row;
251  maxCol = col;
252  }
253  }
254  }
255 if (max< 1e-10)
256  {
257  // The rank is 1. The eigenvalue has multiplicity 2.
258  return 1;
259  }
260 if (maxRow == 2)
261  {
262  for (int col = 0; col < 3; col++)
263  {
264  save = M[1][col];
265  M[1][col] = M[2][col];
266  M[2][col] = save;
267  }
268  }
269 invMax = 1/M[1][maxCol];
270 M[1][0] *= invMax;
271 M[1][1] *= invMax;
272 M[1][2] *= invMax;
273 return 2;
274 }
OT_VECTEUR_3D::change_z
virtual void change_z(double z)
Definition: ot_mathematique.cpp:444
gestionversion.h
fct_taille_metrique.h
FCT_TAILLE_METRIQUE::GetComplement2
void GetComplement2(OT_VECTEUR_3D U, OT_VECTEUR_3D &V, OT_VECTEUR_3D &W)
Definition: fct_taille_metrique.cpp:141
OT_VECTEUR_3D::change_y
virtual void change_y(double y)
Definition: ot_mathematique.cpp:439
OT_VECTEUR_3D::get_x
virtual double get_x(void) const
Definition: ot_mathematique.cpp:417
FCT_TAILLE_METRIQUE::GetComplement1
void GetComplement1(OT_VECTEUR_3D U, OT_VECTEUR_3D V, OT_VECTEUR_3D &W)
Definition: fct_taille_metrique.cpp:135
V
void V(MCAA *mcaa)
Definition: CAD4FE_MCAA.cpp:1794
OT_VECTEUR_3D::norme
virtual void norme(void)
Definition: ot_mathematique.cpp:494
OT_VECTEUR_3D::get_y
virtual double get_y(void) const
Definition: ot_mathematique.cpp:423
OT_VECTEUR_3D
Definition: ot_mathematique.h:94
FCT_TAILLE_METRIQUE::evaluer_decompose
virtual void evaluer_decompose(double *metrique_depart, double *valeur_propre, double *vecteur_propre)
Definition: fct_taille_metrique.cpp:31
sqrt
double2 sqrt(double2 &val)
Definition: ot_doubleprecision.cpp:345
OT_VECTEUR_3D::get_z
virtual double get_z(void) const
Definition: ot_mathematique.cpp:429
FCT_TAILLE_METRIQUE::ComputeRank
int ComputeRank(double *metrique)
Definition: fct_taille_metrique.cpp:168
cos
double2 cos(double2 &val)
Definition: ot_doubleprecision.cpp:206
OT_VECTEUR_3D::change_x
virtual void change_x(double x)
Definition: ot_mathematique.cpp:434
sin
double2 sin(double2 &val)
Definition: ot_doubleprecision.cpp:250