]> SALOME platform Git repositories - modules/smesh.git/blob - SMESHGUI/SMESHGUI_ComputeScalarValue.cxx
Salome HOME
bos #43434 [CEA][Windows] auto in function declaration
[modules/smesh.git] / SMESHGUI / SMESHGUI_ComputeScalarValue.cxx
1 using namespace std;
2 //  File      : SMESHGUI_ComputeScalarValue.cxx
3 //  Created   : Mon Jun 24 14:06:00 2002
4 //  Author    : Nicolas REJNERI
5
6 //  Project   : SALOME
7 //  Module    : SMESH
8 //  Copyright : Open CASCADE 2002
9 //  $Header$
10
11 #include "SMESHGUI_ComputeScalarValue.h"
12 #include "utilities.h"
13 #include <math.h>
14
15
16
17
18 //=============================================================================
19 /*!
20  *
21  */
22 //=============================================================================
23 static double ComputeLength(float* p1, float* p2) {
24   float a1,a2,a3,b1,b2,b3;
25   a1 =  p1[0];
26   a2 =  p1[1];
27   a3 =  p1[2];
28   b1 =  p2[0];
29   b2 =  p2[1];
30   b3 =  p2[2];
31   // MESSAGE( a1 << " "<< a2 << " "<< a3 << " " << b1 << " "<< b2 << " "<< b3 );
32   float X1,Y1,Z1,X2,Z2,Y2;
33   X1 = b1 - a1;
34   Y1 = b2 - a2;
35   Z1 = b3 - a3;
36   // MESSAGE( X1 << " "<< Y1 << " "<< Z1 );
37   float e1;
38   e1 = sqrt( X1*X1 + Y1*Y1 + Z1*Z1 ) ;
39   // MESSAGE( "Length = " << e1 );
40   return e1;
41 }
42
43 //=============================================================================
44 /*!
45  *
46  */
47 //=============================================================================
48 double SMESHGUI_ComputeScalarValue::LengthEdges(vtkCell* theCell) {
49   int num_points = theCell->GetNumberOfPoints ();
50   vtkPoints* points = theCell->GetPoints();
51   if (num_points != 2 ) return 0;
52   float* pnt1 = points->GetPoint(0);
53   float* pnt2 = points->GetPoint(1);
54   return ComputeLength(pnt1,pnt2);
55 };
56
57 //=============================================================================
58 /*!
59  *
60  */
61 //=============================================================================
62 static double ComputeAreaOfTriangle(float* p1, float* p2, float* p3) {
63   double a1,a2,a3,b1,b2,b3,c1,c2,c3;
64   a1 =  p1[0];
65   a2 =  p1[1];
66   a3 =  p1[2];
67   b1 =  p2[0];
68   b2 =  p2[1];
69   b3 =  p2[2];
70   c1 =  p3[0];
71   c2 =  p3[1];
72   c3 =  p3[2];
73   
74   float e1, e2, e3;
75   e1 = sqrt( (a1-b1)*(a1-b1) + (a2-b2)*(a2-b2) + (a3-b3)*(a3-b3) ) ;
76   e2 = sqrt( (b1-c1)*(b1-c1) + (b2-c2)*(b2-c2) + (b3-c3)*(b3-c3) ) ;
77   e3 = sqrt( (c1-a1)*(c1-a1) + (c2-a2)*(c2-a2) + (c3-a3)*(c3-a3) ) ;
78   
79   // MESSAGE( "e = " << e1 << " " << e2 <<"  " << e3 );
80   float s = (e1+e2+e3)/2;
81   double area = sqrt(s*(s-e1)*(s-e2)*(s-e3));
82   // MESSAGE( "area = " << area );
83   return area;
84 }
85
86 //=============================================================================
87 /*!
88  *
89  */
90 //=============================================================================
91 double SMESHGUI_ComputeScalarValue::AreaElements(vtkCell* theCell){
92   // MESSAGE ( " SMESHGUI_ComputeScalarValue::AreaElements " )
93   int num_points = theCell->GetNumberOfPoints ();
94   vtkPoints* points = theCell->GetPoints();
95   // MESSAGE( "num_points =  "<< num_points );
96   for (int j = 0; j < theCell->GetNumberOfPoints (); j++) {
97     float* pnt = points->GetPoint(j);
98     // MESSAGE( pnt[0] << " " << pnt[1] << " " <<  pnt[2] );
99   }
100   if (num_points < 3 ) return 0;
101   if (num_points == 3) {
102     float* p1 = points->GetPoint(0);
103     float* p2 = points->GetPoint(1);
104     float* p3 = points->GetPoint(2);
105     double area = ComputeAreaOfTriangle(p1,p2,p3);
106     return area;
107   } else if (num_points == 4) {
108     float* p1 = points->GetPoint(0);
109     float* p2 = points->GetPoint(1);
110     float* p3 = points->GetPoint(2);
111     float* p4 = points->GetPoint(3);
112     double area1 = ComputeAreaOfTriangle(p1,p2,p3);
113     double area2 = ComputeAreaOfTriangle(p3,p4,p1);
114     return area1+area2;
115   }
116 };
117
118 //=============================================================================
119 /*!
120  *
121  */
122 //=============================================================================
123 double  SMESHGUI_ComputeScalarValue::Taper(vtkCell* theCell){
124   int num_points = theCell->GetNumberOfPoints ();
125   vtkPoints* points = theCell->GetPoints();
126   if (num_points != 4 ) return 0;
127   float* p1 = points->GetPoint(0);
128   float* p2 = points->GetPoint(1);
129   float* p3 = points->GetPoint(2);
130   float* p4 = points->GetPoint(3);
131   double A1 = ComputeAreaOfTriangle(p4,p1,p2);
132   double A2 = ComputeAreaOfTriangle(p3,p1,p2);
133   double A3 = ComputeAreaOfTriangle(p2,p3,p4);
134   double A4 = ComputeAreaOfTriangle(p3,p4,p1);
135   double JA = 0.25 * (A1 + A2 + A3 + A4);
136   double taper = fabs(A1/(JA-1));
137   if (fabs(A2/(JA-1)) > taper) taper = fabs(A2/(JA-1));
138   if (fabs(A3/(JA-1)) > taper) taper = fabs(A3/(JA-1));
139   if (fabs(A4/(JA-1)) > taper) taper = fabs(A4/(JA-1));
140   // MESSAGE( "Taper = " << taper);
141   return taper;
142 };
143
144 //=============================================================================
145 /*!
146  *
147  */
148 //=============================================================================
149 double  SMESHGUI_ComputeScalarValue::AspectRatio(vtkCell* theCell) {
150   int num_points = theCell->GetNumberOfPoints ();
151   vtkPoints* points = theCell->GetPoints();
152   if (num_points < 3 ) return 0;
153   if (num_points == 3) {
154     float a1,a2,a3,b1,b2,b3,c1,c2,c3;
155     float* pnt = points->GetPoint(0);
156     a1 =  pnt[0];
157     a2 =  pnt[1];
158     a3 =  pnt[2];
159     pnt = points->GetPoint(1);
160     b1 =  pnt[0];
161     b2 =  pnt[1];
162     b3 =  pnt[2];
163     pnt = points->GetPoint(2);
164     c1 =  pnt[0];
165     c2 =  pnt[1];
166     c3 =  pnt[2];  
167     
168     float e1, e2, e3;    
169     e1 = sqrt( (a1-b1)*(a1-b1) + (a2-b2)*(a2-b2) + (a3-b3)*(a3-b3) ) ;
170     e2 = sqrt( (b1-c1)*(b1-c1) + (b2-c2)*(b2-c2) + (b3-c3)*(b3-c3) ) ;
171     e3 = sqrt( (c1-a1)*(c1-a1) + (c2-a2)*(c2-a2) + (c3-a3)*(c3-a3) ) ;
172     
173     float amax,p,s;
174
175     amax=e1;
176     if (e2>amax) amax=e2;
177     if (e3>amax) amax=e3;
178     
179     p=(e1+e2+e3)/2;
180     s=AreaElements(theCell);
181   
182     double aspectRatio=amax*p*sqrt(double(3))/(s*6);
183     // MESSAGE( "aspectRatio = " << aspectRatio );
184     return(aspectRatio);
185   }
186   else if (num_points == 4) {
187     float a1,a2,a3,b1,b2,b3,c1,c2,c3,d1,d2,d3;
188     float* pnt = points->GetPoint(0);
189     a1 =  pnt[0];
190     a2 =  pnt[1];
191     a3 =  pnt[2];
192     pnt = points->GetPoint(1);
193     b1 =  pnt[0];
194     b2 =  pnt[1];
195     b3 =  pnt[2];
196     pnt = points->GetPoint(2);
197     c1 =  pnt[0];
198     c2 =  pnt[1];
199     c3 =  pnt[2];
200     pnt = points->GetPoint(3);
201     d1 =  pnt[0];
202     d2 =  pnt[1];
203     d3 =  pnt[2];
204     
205     float e1, e2, e3, e4;
206     float len_min, len_max;
207     e1 = sqrt( (a1-b1)*(a1-b1) + (a2-b2)*(a2-b2) + (a3-b3)*(a3-b3) ) ;
208     e2 = sqrt( (b1-c1)*(b1-c1) + (b2-c2)*(b2-c2) + (b3-c3)*(b3-c3) ) ;
209     e3 = sqrt( (c1-d1)*(c1-d1) + (c2-d2)*(c2-d2) + (c3-d3)*(c3-d3) ) ;
210     e4 = sqrt( (d1-a1)*(d1-a1) + (d2-a2)*(d2-a2) + (d3-a3)*(d3-a3) ) ;
211
212     len_min = e1; len_max = e1;
213
214     if (e2 >len_max ) len_max = e2;
215     if (e3 >len_max ) len_max = e3;
216     if (e4 >len_max ) len_max = e4;
217     if (e2 <len_min ) len_min = e2;
218     if (e3 <len_min ) len_min = e3;
219     if (e4 <len_min ) len_min = e4;
220     
221     return (len_max/len_min);
222   }
223 };
224
225 //=============================================================================
226 /*!
227  *
228  */
229 //=============================================================================
230 static double ComputeAngle(float* p1, float* p2, float* p3) {
231   const double pi=4*atan(double(1));
232   float a1,a2,a3,b1,b2,b3,c1,c2,c3;
233   a1 =  p1[0];
234   a2 =  p1[1];
235   a3 =  p1[2];
236   b1 =  p2[0];
237   b2 =  p2[1];
238   b3 =  p2[2];
239   c1 =  p3[0];
240   c2 =  p3[1];
241   c3 =  p3[2];
242   float X1,Y1,Z1,X2,Z2,Y2;
243   X1 = b1 - a1;
244   X2 = c1 - b1;
245   Y1 = b2 - a2;
246   Y2 = c2 - b2;
247   Z1 = b3 - a3;
248   Z2 = c3 - b3;
249   
250
251   float e1, e2, e3;
252   e1 = sqrt( X1*X1 + Y1*Y1 + Z1*Z1 ) ;
253   e2 = sqrt( X2*X2 + Y2*Y2 + Z2*Z2 ) ;
254   double dot=(X1*(X2)+Y1*(Y2)+Z1*(Z2));
255   //  MESSAGE( pi );
256   //  MESSAGE( dot/(e1*e2) );
257   double cosinus = dot/(e1*e2);
258   cosinus = fabs(cosinus);
259   return 180*acos (cosinus)/pi;
260 }
261
262 //=============================================================================
263 /*!
264  *
265  */
266 //=============================================================================
267 double SMESHGUI_ComputeScalarValue::MinimumAngle(vtkCell* theCell) {
268   int num_points = theCell->GetNumberOfPoints ();
269   vtkPoints* points = theCell->GetPoints();
270   if (num_points < 3 ) return 0;
271   float* pnt1 = points->GetPoint(0);
272   float* pnt2 = points->GetPoint(1);
273   float* pnt3 = points->GetPoint(2);
274   if (num_points == 3) {
275     double a1,a2,a3,amin;
276     a1=fabs(ComputeAngle(pnt1,pnt2,pnt3));
277     amin=a1;
278     a2=fabs(ComputeAngle(pnt2,pnt3,pnt1));
279     if (a2<amin) amin=a2;
280     a3=fabs(ComputeAngle(pnt3,pnt1,pnt2));
281     if (a3<amin) amin=a3;
282     // MESSAGE( "Minimal angle " << amin );
283     return amin;
284   }
285   else if (num_points == 4) {
286     float* pnt4 = points->GetPoint(3);
287     double a1,a2,a3,a4,amin;
288     a1=fabs(ComputeAngle(pnt1,pnt2,pnt3));
289     amin=a1;
290     a2=fabs(ComputeAngle(pnt2,pnt3,pnt4));
291     if (a2<amin) amin=a2;
292     a3=fabs(ComputeAngle(pnt3,pnt4,pnt1));
293     if (a3<amin) amin=a3;
294     a4=fabs(ComputeAngle(pnt4,pnt1,pnt2));
295     if (a4<amin) amin=a4;
296     
297     // MESSAGE( "Minimal angle " << amin );
298     return amin;
299   }
300 };
301
302 //=============================================================================
303 /*!
304  *
305  */
306 //=============================================================================
307 double  SMESHGUI_ComputeScalarValue::Skew(vtkCell* theCell) {
308   int num_points = theCell->GetNumberOfPoints ();
309   vtkPoints* points = theCell->GetPoints();
310   if (num_points < 3 ) return 0;
311   //triangle case
312   if (num_points == 3) {
313     float* pnt1 = points->GetPoint(0);
314     float* pnt2 = points->GetPoint(1);
315     float* pnt3 = points->GetPoint(2);
316     double a1,a2,a3,amax;
317     a1=fabs(60 - fabs(ComputeAngle(pnt1,pnt2,pnt3)));
318     amax=a1;
319     a2=fabs(60 - fabs(ComputeAngle(pnt2,pnt3,pnt1)));
320     if (a2>amax) amax=a2;
321     a3=fabs(60 - fabs(ComputeAngle(pnt3,pnt1,pnt2)));
322     if (a3>amax) amax=a3;
323     // MESSAGE( "Skew = " << amax );
324     return amax;
325   } 
326   //quadrangle case
327   else if (num_points == 4) {
328     float* pnt1 = points->GetPoint(0);
329     float* pnt2 = points->GetPoint(1);
330     float* pnt3 = points->GetPoint(2);
331     float* pnt4 = points->GetPoint(3);
332     
333     double a1,a2,a3,a4,amax;
334     a1=fabs(90 - fabs(ComputeAngle(pnt1,pnt2,pnt3)));
335     amax=a1;
336     a2=fabs(90 - fabs(ComputeAngle(pnt2,pnt3,pnt4)));
337     if (a2>amax) amax=a2;
338     a3=fabs(90 - fabs(ComputeAngle(pnt3,pnt4,pnt1)));
339     if (a3>amax) amax=a3;
340     a4=fabs(90 - fabs(ComputeAngle(pnt4,pnt1,pnt2)));
341     if (a4>amax) amax=a4;
342     // MESSAGE( "Skew = " << amax );
343     return amax;
344   }
345 };
346
347 //=============================================================================
348 /*!
349  *
350  */
351 //=============================================================================
352 static double ComputeA(float* p1, float* p2, float* p3, float* G) {
353   double e1 = sqrt(pow(p2[0]-p1[0], 2)+pow(p2[1]-p1[1], 2)+pow(p2[2]-p1[2], 2));
354   double e2 = sqrt(pow(p3[0]-p2[0], 2)+pow(p3[1]-p2[1], 2)+pow(p3[2]-p2[2], 2));
355   double l;
356   if (e1 < e2) l = 0.5*e1;
357   else l = 0.5*e2;
358   float GI[3], GJ[3], N[3];;
359   GI[0] = (p2[0]-p1[0])/2-G[0];
360   GI[1] = (p2[1]-p1[1])/2-G[1];
361   GI[2] = (p2[2]-p1[2])/2-G[2];
362   
363   GJ[0] = (p3[0]-p2[0])/2-G[0];
364   GJ[1] = (p3[1]-p2[1])/2-G[1];
365   GJ[2] = (p3[2]-p2[2])/2-G[2];
366   
367   N[0] = GI[1]*GJ[2] - GI[2]*GJ[1];
368   N[1] = GI[2]*GJ[0] - GI[0]*GJ[2];
369   N[2] = GI[0]*GJ[1] - GI[1]*GJ[0];
370   
371   double H;
372   float T[3];
373   T[0] = (p1[0]-G[0])*N[0];
374   T[1] = (p1[1]-G[1])*N[1];
375   T[2] = (p1[2]-G[2])*N[2];
376   
377   H = sqrt(pow(T[0],2)+pow(T[1],2)+pow(T[2],2))/sqrt(pow(N[0],2)+pow(N[1],2)+pow(N[2],2));
378   double A = H/l;
379   return A;
380 }
381
382 //=============================================================================
383 /*!
384  *
385  */
386 //=============================================================================
387 double  SMESHGUI_ComputeScalarValue::Warp(vtkCell* theCell) {
388   int num_points = theCell->GetNumberOfPoints ();
389   vtkPoints* points = theCell->GetPoints();
390   if (num_points != 4 ) return 0;
391   float* p1 = points->GetPoint(0);
392   float* p2 = points->GetPoint(1);
393   float* p3 = points->GetPoint(2);
394   float* p4 = points->GetPoint(3);
395   double G1, G2, G3;
396   float G[3];
397   G[0] = (p1[0]+p2[0]+p3[0]+p4[0])/4;
398   G[1] = (p1[1]+p2[1]+p3[1]+p4[1])/4;
399   G[2] = (p1[2]+p2[2]+p3[2]+p4[2])/4;
400   double amax = ComputeA(p1, p2, p3, G);
401   double nextA = ComputeA(p2, p3, p4, G);
402   if (nextA > amax) amax = nextA;
403   nextA = ComputeA(p3, p4, p1, G);
404   if (nextA > amax) amax = nextA;
405   nextA = ComputeA(p4, p1, p2, G);
406   if (nextA > amax) amax = nextA;
407   // MESSAGE( "Warp = " << amax );
408   return amax;
409 }