1 // SMESH SMESHGUI : GUI for SMESH component
3 // Copyright (C) 2003 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
4 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
6 // This library is free software; you can redistribute it and/or
7 // modify it under the terms of the GNU Lesser General Public
8 // License as published by the Free Software Foundation; either
9 // version 2.1 of the License.
11 // This library is distributed in the hope that it will be useful,
12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14 // Lesser General Public License for more details.
16 // You should have received a copy of the GNU Lesser General Public
17 // License along with this library; if not, write to the Free Software
18 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
20 // See http://www.opencascade.org/SALOME/ or email : webmaster.salome@opencascade.org
24 // File : SMESHGUI_ComputeScalarValue.cxx
25 // Author : Nicolas REJNERI
30 #include "SMESHGUI_ComputeScalarValue.h"
31 #include "utilities.h"
34 //=============================================================================
38 //=============================================================================
39 static double ComputeLength(float* p1, float* p2) {
40 float a1,a2,a3,b1,b2,b3;
47 // MESSAGE( a1 << " "<< a2 << " "<< a3 << " " << b1 << " "<< b2 << " "<< b3 );
48 float X1,Y1,Z1,X2,Z2,Y2;
52 // MESSAGE( X1 << " "<< Y1 << " "<< Z1 );
54 e1 = sqrt( X1*X1 + Y1*Y1 + Z1*Z1 ) ;
55 // MESSAGE( "Length = " << e1 );
59 //=============================================================================
63 //=============================================================================
64 double SMESHGUI_ComputeScalarValue::LengthEdges(vtkCell* theCell) {
65 int num_points = theCell->GetNumberOfPoints ();
66 vtkPoints* points = theCell->GetPoints();
67 if (num_points != 2 ) return 0;
68 float* pnt1 = points->GetPoint(0);
69 float* pnt2 = points->GetPoint(1);
70 return ComputeLength(pnt1,pnt2);
73 //=============================================================================
77 //=============================================================================
78 static double ComputeAreaOfTriangle(float* p1, float* p2, float* p3) {
79 double a1,a2,a3,b1,b2,b3,c1,c2,c3;
91 e1 = sqrt( (a1-b1)*(a1-b1) + (a2-b2)*(a2-b2) + (a3-b3)*(a3-b3) ) ;
92 e2 = sqrt( (b1-c1)*(b1-c1) + (b2-c2)*(b2-c2) + (b3-c3)*(b3-c3) ) ;
93 e3 = sqrt( (c1-a1)*(c1-a1) + (c2-a2)*(c2-a2) + (c3-a3)*(c3-a3) ) ;
95 // MESSAGE( "e = " << e1 << " " << e2 <<" " << e3 );
96 float s = (e1+e2+e3)/2;
97 double area = sqrt(s*(s-e1)*(s-e2)*(s-e3));
98 // MESSAGE( "area = " << area );
102 //=============================================================================
106 //=============================================================================
107 double SMESHGUI_ComputeScalarValue::AreaElements(vtkCell* theCell){
108 // MESSAGE ( " SMESHGUI_ComputeScalarValue::AreaElements " )
109 int num_points = theCell->GetNumberOfPoints ();
110 vtkPoints* points = theCell->GetPoints();
111 // MESSAGE( "num_points = "<< num_points );
112 for (int j = 0; j < theCell->GetNumberOfPoints (); j++) {
113 float* pnt = points->GetPoint(j);
114 // MESSAGE( pnt[0] << " " << pnt[1] << " " << pnt[2] );
116 if (num_points < 3 ) return 0;
117 if (num_points == 3) {
118 float* p1 = points->GetPoint(0);
119 float* p2 = points->GetPoint(1);
120 float* p3 = points->GetPoint(2);
121 double area = ComputeAreaOfTriangle(p1,p2,p3);
123 } else if (num_points == 4) {
124 float* p1 = points->GetPoint(0);
125 float* p2 = points->GetPoint(1);
126 float* p3 = points->GetPoint(2);
127 float* p4 = points->GetPoint(3);
128 double area1 = ComputeAreaOfTriangle(p1,p2,p3);
129 double area2 = ComputeAreaOfTriangle(p3,p4,p1);
134 //=============================================================================
138 //=============================================================================
139 double SMESHGUI_ComputeScalarValue::Taper(vtkCell* theCell){
140 int num_points = theCell->GetNumberOfPoints ();
141 vtkPoints* points = theCell->GetPoints();
142 if (num_points != 4 ) return 0;
143 float* p1 = points->GetPoint(0);
144 float* p2 = points->GetPoint(1);
145 float* p3 = points->GetPoint(2);
146 float* p4 = points->GetPoint(3);
147 double A1 = ComputeAreaOfTriangle(p4,p1,p2);
148 double A2 = ComputeAreaOfTriangle(p3,p1,p2);
149 double A3 = ComputeAreaOfTriangle(p2,p3,p4);
150 double A4 = ComputeAreaOfTriangle(p3,p4,p1);
151 double JA = 0.25 * (A1 + A2 + A3 + A4);
152 double taper = fabs(A1/(JA-1));
153 if (fabs(A2/(JA-1)) > taper) taper = fabs(A2/(JA-1));
154 if (fabs(A3/(JA-1)) > taper) taper = fabs(A3/(JA-1));
155 if (fabs(A4/(JA-1)) > taper) taper = fabs(A4/(JA-1));
156 // MESSAGE( "Taper = " << taper);
160 //=============================================================================
164 //=============================================================================
165 double SMESHGUI_ComputeScalarValue::AspectRatio(vtkCell* theCell) {
166 int num_points = theCell->GetNumberOfPoints ();
167 vtkPoints* points = theCell->GetPoints();
168 if (num_points < 3 ) return 0;
169 if (num_points == 3) {
170 float a1,a2,a3,b1,b2,b3,c1,c2,c3;
171 float* pnt = points->GetPoint(0);
175 pnt = points->GetPoint(1);
179 pnt = points->GetPoint(2);
185 e1 = sqrt( (a1-b1)*(a1-b1) + (a2-b2)*(a2-b2) + (a3-b3)*(a3-b3) ) ;
186 e2 = sqrt( (b1-c1)*(b1-c1) + (b2-c2)*(b2-c2) + (b3-c3)*(b3-c3) ) ;
187 e3 = sqrt( (c1-a1)*(c1-a1) + (c2-a2)*(c2-a2) + (c3-a3)*(c3-a3) ) ;
192 if (e2>amax) amax=e2;
193 if (e3>amax) amax=e3;
196 s=AreaElements(theCell);
198 double aspectRatio=amax*p*sqrt(double(3))/(s*6);
199 // MESSAGE( "aspectRatio = " << aspectRatio );
202 else if (num_points == 4) {
203 float a1,a2,a3,b1,b2,b3,c1,c2,c3,d1,d2,d3;
204 float* pnt = points->GetPoint(0);
208 pnt = points->GetPoint(1);
212 pnt = points->GetPoint(2);
216 pnt = points->GetPoint(3);
221 float e1, e2, e3, e4;
222 float len_min, len_max;
223 e1 = sqrt( (a1-b1)*(a1-b1) + (a2-b2)*(a2-b2) + (a3-b3)*(a3-b3) ) ;
224 e2 = sqrt( (b1-c1)*(b1-c1) + (b2-c2)*(b2-c2) + (b3-c3)*(b3-c3) ) ;
225 e3 = sqrt( (c1-d1)*(c1-d1) + (c2-d2)*(c2-d2) + (c3-d3)*(c3-d3) ) ;
226 e4 = sqrt( (d1-a1)*(d1-a1) + (d2-a2)*(d2-a2) + (d3-a3)*(d3-a3) ) ;
228 len_min = e1; len_max = e1;
230 if (e2 >len_max ) len_max = e2;
231 if (e3 >len_max ) len_max = e3;
232 if (e4 >len_max ) len_max = e4;
233 if (e2 <len_min ) len_min = e2;
234 if (e3 <len_min ) len_min = e3;
235 if (e4 <len_min ) len_min = e4;
237 return (len_max/len_min);
241 //=============================================================================
245 //=============================================================================
246 static double ComputeAngle(float* p1, float* p2, float* p3) {
247 const double pi=4*atan(double(1));
248 float a1,a2,a3,b1,b2,b3,c1,c2,c3;
258 float X1,Y1,Z1,X2,Z2,Y2;
268 e1 = sqrt( X1*X1 + Y1*Y1 + Z1*Z1 ) ;
269 e2 = sqrt( X2*X2 + Y2*Y2 + Z2*Z2 ) ;
270 double dot=(X1*(X2)+Y1*(Y2)+Z1*(Z2));
272 // MESSAGE( dot/(e1*e2) );
273 double cosinus = dot/(e1*e2);
274 cosinus = fabs(cosinus);
275 return 180*acos (cosinus)/pi;
278 //=============================================================================
282 //=============================================================================
283 double SMESHGUI_ComputeScalarValue::MinimumAngle(vtkCell* theCell) {
284 int num_points = theCell->GetNumberOfPoints ();
285 vtkPoints* points = theCell->GetPoints();
286 if (num_points < 3 ) return 0;
287 float* pnt1 = points->GetPoint(0);
288 float* pnt2 = points->GetPoint(1);
289 float* pnt3 = points->GetPoint(2);
290 if (num_points == 3) {
291 double a1,a2,a3,amin;
292 a1=fabs(ComputeAngle(pnt1,pnt2,pnt3));
294 a2=fabs(ComputeAngle(pnt2,pnt3,pnt1));
295 if (a2<amin) amin=a2;
296 a3=fabs(ComputeAngle(pnt3,pnt1,pnt2));
297 if (a3<amin) amin=a3;
298 // MESSAGE( "Minimal angle " << amin );
301 else if (num_points == 4) {
302 float* pnt4 = points->GetPoint(3);
303 double a1,a2,a3,a4,amin;
304 a1=fabs(ComputeAngle(pnt1,pnt2,pnt3));
306 a2=fabs(ComputeAngle(pnt2,pnt3,pnt4));
307 if (a2<amin) amin=a2;
308 a3=fabs(ComputeAngle(pnt3,pnt4,pnt1));
309 if (a3<amin) amin=a3;
310 a4=fabs(ComputeAngle(pnt4,pnt1,pnt2));
311 if (a4<amin) amin=a4;
313 // MESSAGE( "Minimal angle " << amin );
318 //=============================================================================
322 //=============================================================================
323 double SMESHGUI_ComputeScalarValue::Skew(vtkCell* theCell) {
324 int num_points = theCell->GetNumberOfPoints ();
325 vtkPoints* points = theCell->GetPoints();
326 if (num_points < 3 ) return 0;
328 if (num_points == 3) {
329 float* pnt1 = points->GetPoint(0);
330 float* pnt2 = points->GetPoint(1);
331 float* pnt3 = points->GetPoint(2);
332 double a1,a2,a3,amax;
333 a1=fabs(60 - fabs(ComputeAngle(pnt1,pnt2,pnt3)));
335 a2=fabs(60 - fabs(ComputeAngle(pnt2,pnt3,pnt1)));
336 if (a2>amax) amax=a2;
337 a3=fabs(60 - fabs(ComputeAngle(pnt3,pnt1,pnt2)));
338 if (a3>amax) amax=a3;
339 // MESSAGE( "Skew = " << amax );
343 else if (num_points == 4) {
344 float* pnt1 = points->GetPoint(0);
345 float* pnt2 = points->GetPoint(1);
346 float* pnt3 = points->GetPoint(2);
347 float* pnt4 = points->GetPoint(3);
349 double a1,a2,a3,a4,amax;
350 a1=fabs(90 - fabs(ComputeAngle(pnt1,pnt2,pnt3)));
352 a2=fabs(90 - fabs(ComputeAngle(pnt2,pnt3,pnt4)));
353 if (a2>amax) amax=a2;
354 a3=fabs(90 - fabs(ComputeAngle(pnt3,pnt4,pnt1)));
355 if (a3>amax) amax=a3;
356 a4=fabs(90 - fabs(ComputeAngle(pnt4,pnt1,pnt2)));
357 if (a4>amax) amax=a4;
358 // MESSAGE( "Skew = " << amax );
363 //=============================================================================
367 //=============================================================================
368 static double ComputeA(float* p1, float* p2, float* p3, float* G) {
369 double e1 = sqrt(pow(p2[0]-p1[0], 2)+pow(p2[1]-p1[1], 2)+pow(p2[2]-p1[2], 2));
370 double e2 = sqrt(pow(p3[0]-p2[0], 2)+pow(p3[1]-p2[1], 2)+pow(p3[2]-p2[2], 2));
372 if (e1 < e2) l = 0.5*e1;
374 float GI[3], GJ[3], N[3];;
375 GI[0] = (p2[0]-p1[0])/2-G[0];
376 GI[1] = (p2[1]-p1[1])/2-G[1];
377 GI[2] = (p2[2]-p1[2])/2-G[2];
379 GJ[0] = (p3[0]-p2[0])/2-G[0];
380 GJ[1] = (p3[1]-p2[1])/2-G[1];
381 GJ[2] = (p3[2]-p2[2])/2-G[2];
383 N[0] = GI[1]*GJ[2] - GI[2]*GJ[1];
384 N[1] = GI[2]*GJ[0] - GI[0]*GJ[2];
385 N[2] = GI[0]*GJ[1] - GI[1]*GJ[0];
389 T[0] = (p1[0]-G[0])*N[0];
390 T[1] = (p1[1]-G[1])*N[1];
391 T[2] = (p1[2]-G[2])*N[2];
393 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));
398 //=============================================================================
402 //=============================================================================
403 double SMESHGUI_ComputeScalarValue::Warp(vtkCell* theCell) {
404 int num_points = theCell->GetNumberOfPoints ();
405 vtkPoints* points = theCell->GetPoints();
406 if (num_points != 4 ) return 0;
407 float* p1 = points->GetPoint(0);
408 float* p2 = points->GetPoint(1);
409 float* p3 = points->GetPoint(2);
410 float* p4 = points->GetPoint(3);
413 G[0] = (p1[0]+p2[0]+p3[0]+p4[0])/4;
414 G[1] = (p1[1]+p2[1]+p3[1]+p4[1])/4;
415 G[2] = (p1[2]+p2[2]+p3[2]+p4[2])/4;
416 double amax = ComputeA(p1, p2, p3, G);
417 double nextA = ComputeA(p2, p3, p4, G);
418 if (nextA > amax) amax = nextA;
419 nextA = ComputeA(p3, p4, p1, G);
420 if (nextA > amax) amax = nextA;
421 nextA = ComputeA(p4, p1, p2, G);
422 if (nextA > amax) amax = nextA;
423 // MESSAGE( "Warp = " << amax );