src/SMESHGUI/SMESHGUI_ComputeScalarValue.cxx

   1 //  SMESH SMESHGUI : GUI for SMESH component
   2 //
   3 //  Copyright (C) 2003  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
   4 //  CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
   5 //
   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.
  10 //
  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.
  15 //
  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
  19 //
  20 //  See http://www.opencascade.org/SALOME/ or email : webmaster.salome@opencascade.org
  21 //
  22 //
  23 //
  24 //  File   : SMESHGUI_ComputeScalarValue.cxx
  25 //  Author : Nicolas REJNERI
  26 //  Module : SMESH
  27 //  $Header$
  28
  29 using namespace std;
  30 #include "SMESHGUI_ComputeScalarValue.h"
  31 #include "utilities.h"
  32 #include <math.h>
  33
  34 //=============================================================================
  35 /*!
  36  *
  37  */
  38 //=============================================================================
  39 static double ComputeLength(float* p1, float* p2) {
  40   float a1,a2,a3,b1,b2,b3;
  41   a1 =  p1[0];
  42   a2 =  p1[1];
  43   a3 =  p1[2];
  44   b1 =  p2[0];
  45   b2 =  p2[1];
  46   b3 =  p2[2];
  47   // MESSAGE( a1 << " "<< a2 << " "<< a3 << " " << b1 << " "<< b2 << " "<< b3 );
  48   float X1,Y1,Z1,X2,Z2,Y2;
  49   X1 = b1 - a1;
  50   Y1 = b2 - a2;
  51   Z1 = b3 - a3;
  52   // MESSAGE( X1 << " "<< Y1 << " "<< Z1 );
  53   float e1;
  54   e1 = sqrt( X1*X1 + Y1*Y1 + Z1*Z1 ) ;
  55   // MESSAGE( "Length = " << e1 );
  56   return e1;
  57 }
  58
  59 //=============================================================================
  60 /*!
  61  *
  62  */
  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);
  71 };
  72
  73 //=============================================================================
  74 /*!
  75  *
  76  */
  77 //=============================================================================
  78 static double ComputeAreaOfTriangle(float* p1, float* p2, float* p3) {
  79   double a1,a2,a3,b1,b2,b3,c1,c2,c3;
  80   a1 =  p1[0];
  81   a2 =  p1[1];
  82   a3 =  p1[2];
  83   b1 =  p2[0];
  84   b2 =  p2[1];
  85   b3 =  p2[2];
  86   c1 =  p3[0];
  87   c2 =  p3[1];
  88   c3 =  p3[2];
  89
  90   float e1, e2, e3;
  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) ) ;
  94
  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 );
  99   return area;
 100 }
 101
 102 //=============================================================================
 103 /*!
 104  *
 105  */
 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] );
 115   }
 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);
 122     return area;
 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);
 130     return area1+area2;
 131   }
 132 };
 133
 134 //=============================================================================
 135 /*!
 136  *
 137  */
 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);
 157   return taper;
 158 };
 159
 160 //=============================================================================
 161 /*!
 162  *
 163  */
 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);
 172     a1 =  pnt[0];
 173     a2 =  pnt[1];
 174     a3 =  pnt[2];
 175     pnt = points->GetPoint(1);
 176     b1 =  pnt[0];
 177     b2 =  pnt[1];
 178     b3 =  pnt[2];
 179     pnt = points->GetPoint(2);
 180     c1 =  pnt[0];
 181     c2 =  pnt[1];
 182     c3 =  pnt[2];
 183
 184     float e1, e2, e3;
 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) ) ;
 188
 189     float amax,p,s;
 190
 191     amax=e1;
 192     if (e2>amax) amax=e2;
 193     if (e3>amax) amax=e3;
 194
 195     p=(e1+e2+e3)/2;
 196     s=AreaElements(theCell);
 197
 198     double aspectRatio=amax*p*sqrt(double(3))/(s*6);
 199     // MESSAGE( "aspectRatio = " << aspectRatio );
 200     return(aspectRatio);
 201   }
 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);
 205     a1 =  pnt[0];
 206     a2 =  pnt[1];
 207     a3 =  pnt[2];
 208     pnt = points->GetPoint(1);
 209     b1 =  pnt[0];
 210     b2 =  pnt[1];
 211     b3 =  pnt[2];
 212     pnt = points->GetPoint(2);
 213     c1 =  pnt[0];
 214     c2 =  pnt[1];
 215     c3 =  pnt[2];
 216     pnt = points->GetPoint(3);
 217     d1 =  pnt[0];
 218     d2 =  pnt[1];
 219     d3 =  pnt[2];
 220
 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) ) ;
 227
 228     len_min = e1; len_max = e1;
 229
 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;
 236
 237     return (len_max/len_min);
 238   }
 239 };
 240
 241 //=============================================================================
 242 /*!
 243  *
 244  */
 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;
 249   a1 =  p1[0];
 250   a2 =  p1[1];
 251   a3 =  p1[2];
 252   b1 =  p2[0];
 253   b2 =  p2[1];
 254   b3 =  p2[2];
 255   c1 =  p3[0];
 256   c2 =  p3[1];
 257   c3 =  p3[2];
 258   float X1,Y1,Z1,X2,Z2,Y2;
 259   X1 = b1 - a1;
 260   X2 = c1 - b1;
 261   Y1 = b2 - a2;
 262   Y2 = c2 - b2;
 263   Z1 = b3 - a3;
 264   Z2 = c3 - b3;
 265
 266
 267   float e1, e2, e3;
 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));
 271   //  MESSAGE( pi );
 272   //  MESSAGE( dot/(e1*e2) );
 273   double cosinus = dot/(e1*e2);
 274   cosinus = fabs(cosinus);
 275   return 180*acos (cosinus)/pi;
 276 }
 277
 278 //=============================================================================
 279 /*!
 280  *
 281  */
 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));
 293     amin=a1;
 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 );
 299     return amin;
 300   }
 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));
 305     amin=a1;
 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;
 312
 313     // MESSAGE( "Minimal angle " << amin );
 314     return amin;
 315   }
 316 };
 317
 318 //=============================================================================
 319 /*!
 320  *
 321  */
 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;
 327   //triangle case
 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)));
 334     amax=a1;
 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 );
 340     return amax;
 341   }
 342   //quadrangle case
 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);
 348
 349     double a1,a2,a3,a4,amax;
 350     a1=fabs(90 - fabs(ComputeAngle(pnt1,pnt2,pnt3)));
 351     amax=a1;
 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 );
 359     return amax;
 360   }
 361 };
 362
 363 //=============================================================================
 364 /*!
 365  *
 366  */
 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));
 371   double l;
 372   if (e1 < e2) l = 0.5*e1;
 373   else l = 0.5*e2;
 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];
 378
 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];
 382
 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];
 386
 387   double H;
 388   float T[3];
 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];
 392
 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));
 394   double A = H/l;
 395   return A;
 396 }
 397
 398 //=============================================================================
 399 /*!
 400  *
 401  */
 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);
 411   double G1, G2, G3;
 412   float G[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 );
 424   return amax;
 425 }