2 // File : SMESHGUI_ComputeScalarValue.cxx
3 // Created : Mon Jun 24 14:06:00 2002
4 // Author : Nicolas REJNERI
8 // Copyright : Open CASCADE 2002
11 #include "SMESHGUI_ComputeScalarValue.h"
12 #include "utilities.h"
18 //=============================================================================
22 //=============================================================================
23 static double ComputeLength(float* p1, float* p2) {
24 float a1,a2,a3,b1,b2,b3;
31 // MESSAGE( a1 << " "<< a2 << " "<< a3 << " " << b1 << " "<< b2 << " "<< b3 );
32 float X1,Y1,Z1,X2,Z2,Y2;
36 // MESSAGE( X1 << " "<< Y1 << " "<< Z1 );
38 e1 = sqrt( X1*X1 + Y1*Y1 + Z1*Z1 ) ;
39 // MESSAGE( "Length = " << e1 );
43 //=============================================================================
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);
57 //=============================================================================
61 //=============================================================================
62 static double ComputeAreaOfTriangle(float* p1, float* p2, float* p3) {
63 double a1,a2,a3,b1,b2,b3,c1,c2,c3;
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) ) ;
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 );
86 //=============================================================================
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] );
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);
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);
118 //=============================================================================
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);
144 //=============================================================================
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);
159 pnt = points->GetPoint(1);
163 pnt = points->GetPoint(2);
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) ) ;
176 if (e2>amax) amax=e2;
177 if (e3>amax) amax=e3;
180 s=AreaElements(theCell);
182 double aspectRatio=amax*p*sqrt(double(3))/(s*6);
183 // MESSAGE( "aspectRatio = " << aspectRatio );
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);
192 pnt = points->GetPoint(1);
196 pnt = points->GetPoint(2);
200 pnt = points->GetPoint(3);
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) ) ;
212 len_min = e1; len_max = e1;
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;
221 return (len_max/len_min);
225 //=============================================================================
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;
242 float X1,Y1,Z1,X2,Z2,Y2;
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));
256 // MESSAGE( dot/(e1*e2) );
257 double cosinus = dot/(e1*e2);
258 cosinus = fabs(cosinus);
259 return 180*acos (cosinus)/pi;
262 //=============================================================================
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));
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 );
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));
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;
297 // MESSAGE( "Minimal angle " << amin );
302 //=============================================================================
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;
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)));
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 );
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);
333 double a1,a2,a3,a4,amax;
334 a1=fabs(90 - fabs(ComputeAngle(pnt1,pnt2,pnt3)));
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 );
347 //=============================================================================
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));
356 if (e1 < e2) l = 0.5*e1;
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];
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];
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];
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];
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));
382 //=============================================================================
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);
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 );