src/StdMeshers/StdMeshers_Distribution.cxx

   1 // Copyright (C) 2007-2021  CEA/DEN, EDF R&D, OPEN CASCADE
   2 //
   3 // Copyright (C) 2003-2007  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, or (at your option) any later version.
  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.salome-platform.org/ or email : webmaster.salome@opencascade.com
  21 //
  22
  23 //  SMESH StdMeshers : implementation of point distribution algorithm
  24 //  File   : StdMeshers_Distribution.cxx
  25 //  Author : Alexandre SOLOVYOV
  26 //  Module : SMESH
  27 //  $Header$
  28 //
  29 #include "StdMeshers_Distribution.hxx"
  30
  31 #include <math_GaussSingleIntegration.hxx>
  32 #include <utilities.h>
  33
  34 #include <Standard_Failure.hxx>
  35 #include <Expr_NamedUnknown.hxx>
  36 #include <Standard_ErrorHandler.hxx>
  37
  38 using namespace std;
  39
  40 namespace StdMeshers {
  41
  42 Function::Function( const int conv )
  43 : myConv( conv )
  44 {
  45 }
  46
  47 Function::~Function()
  48 {
  49 }
  50
  51 bool Function::value( const double, double& f ) const
  52 {
  53   bool ok = true;
  54   if (myConv == 0) {
  55     try {
  56       OCC_CATCH_SIGNALS;
  57       f = pow( 10., f );
  58     } catch(Standard_Failure&) {
  59       f = 0.0;
  60       ok = false;
  61     }
  62   }
  63   else if( myConv==1 && f<0.0 )
  64     f = 0.0;
  65
  66   return ok;
  67 }
  68
  69 FunctionIntegral::FunctionIntegral( const Function* f, const double st )
  70 : Function( -1 ),
  71   myFunc( const_cast<Function*>( f ) ),
  72   myStart( st )
  73 {
  74 }
  75
  76 FunctionIntegral::~FunctionIntegral()
  77 {
  78 }
  79
  80 bool FunctionIntegral::value( const double t, double& f ) const
  81 {
  82   f = myFunc ? myFunc->integral( myStart, t ) : 0;
  83   return myFunc!=0 && Function::value( t, f );
  84 }
  85
  86 double FunctionIntegral::integral( const double, const double ) const
  87 {
  88   return 0;
  89 }
  90
  91 FunctionTable::FunctionTable( const std::vector<double>& data, const int conv )
  92 : Function( conv )
  93 {
  94   myData = data;
  95 }
  96
  97 FunctionTable::~FunctionTable()
  98 {
  99 }
 100
 101 bool FunctionTable::value( const double t, double& f ) const
 102 {
 103   int i1, i2;
 104   if( !findBounds( t, i1, i2 ) )
 105     return false;
 106
 107   if( i1==i2 ) {
 108     f = myData[ 2*i1+1 ];
 109     Function::value( t, f );
 110     return true;
 111   }
 112
 113   double
 114     x1 = myData[2*i1], y1 = myData[2*i1+1],
 115     x2 = myData[2*i2], y2 = myData[2*i2+1];
 116
 117   Function::value( x1, y1 );
 118   Function::value( x2, y2 );
 119
 120   f = y1 + ( y2-y1 ) * ( t-x1 ) / ( x2-x1 );
 121   return true;
 122 }
 123
 124 double FunctionTable::integral( const int i ) const
 125 {
 126   if ( i >= 0  &&  i < (int)myData.size()-1 )
 127     return integral( i, myData[2*(i+1)] - myData[2*i] );
 128   else
 129     return 0;
 130 }
 131
 132 double FunctionTable::integral( const int i, const double d ) const
 133 {
 134   double f1,f2, res = 0.0;
 135   if( value( myData[2*i]+d, f1 ) )
 136     if(!value(myData[2*i], f2)) {
 137       f2 = myData[2*i+1];
 138       Function::value( 1, f2 );
 139     }
 140   res = (f2+f1) * d / 2.0;
 141   return res;
 142 }
 143
 144 double FunctionTable::integral( const double a, const double b ) const
 145 {
 146   int x1s, x1f, x2s, x2f;
 147   findBounds( a, x1s, x1f );
 148   findBounds( b, x2s, x2f );
 149   double J = 0;
 150   for( int i=x1s; i<x2s; i++ )
 151     J+=integral( i );
 152   J-=integral( x1s, a-myData[2*x1s] );
 153   J+=integral( x2s, b-myData[2*x2s] );
 154   return J;
 155 }
 156
 157 bool FunctionTable::findBounds( const double x, int& x_ind_1, int& x_ind_2 ) const
 158 {
 159   int n = (int) myData.size() / 2;
 160   if( n==0 || x<myData[0] )
 161   {
 162     x_ind_1 = x_ind_2 = 0;
 163     return false;
 164   }
 165
 166   for( int i=0; i<n-1; i++ )
 167     if( myData[2*i]<=x && x<myData[2*(i+1)] )
 168     {
 169       x_ind_1 = i;
 170       x_ind_2 = i+1;
 171       return true;
 172     }
 173   x_ind_1 = n-1;
 174   x_ind_2 = n-1;
 175   return ( fabs( x - myData[2*x_ind_2] ) < 1.e-10 );
 176 }
 177
 178 FunctionExpr::FunctionExpr( const char* str, const int conv )
 179 : Function( conv ),
 180   myVars( 1, 1 ),
 181   myValues( 1, 1 )
 182 {
 183   bool ok = true;
 184   try {
 185     OCC_CATCH_SIGNALS;
 186     myExpr = ExprIntrp_GenExp::Create();
 187     myExpr->Process( ( Standard_CString )str );
 188   } catch(Standard_Failure&) {
 189     ok = false;
 190   }
 191
 192   if( !ok || !myExpr->IsDone() )
 193     myExpr.Nullify();
 194
 195   myVars.ChangeValue( 1 ) = new Expr_NamedUnknown( "t" );
 196 }
 197
 198 FunctionExpr::~FunctionExpr()
 199 {
 200 }
 201
 202 Standard_Boolean FunctionExpr::Value( const Standard_Real T, Standard_Real& F )
 203 {
 204   double f;
 205   Standard_Boolean res = value( T, f );
 206   F = f;
 207   return res;
 208 }
 209
 210 bool FunctionExpr::value( const double t, double& f ) const
 211 {
 212   if( myExpr.IsNull() )
 213     return false;
 214
 215   ( ( TColStd_Array1OfReal& )myValues ).ChangeValue( 1 ) = t;
 216   bool ok = true;
 217   try {
 218     OCC_CATCH_SIGNALS;
 219     f = myExpr->Expression()->Evaluate( myVars, myValues );
 220   } catch(Standard_Failure&) {
 221     f = 0.0;
 222     ok = false;
 223   }
 224
 225   ok = Function::value( t, f ) && ok;
 226   return ok;
 227 }
 228
 229 double FunctionExpr::integral( const double a, const double b ) const
 230 {
 231   double res = 0.0;
 232   try {
 233     OCC_CATCH_SIGNALS;
 234     math_GaussSingleIntegration _int
 235       ( *static_cast<math_Function*>( const_cast<FunctionExpr*> (this) ), a, b, 20 );
 236     if( _int.IsDone() )
 237       res = _int.Value();
 238   } catch(Standard_Failure&) {
 239     res = 0.0;
 240     MESSAGE( "Exception in integral calculating" );
 241   }
 242   return res;
 243 }
 244
 245 double dihotomySolve( Function& f, const double val, const double _start, const double _fin, const double eps, bool& ok )
 246 {
 247   double start = _start, fin = _fin, start_val, fin_val; bool ok1, ok2;
 248   ok1 = f.value( start, start_val );
 249   ok2 = f.value( fin, fin_val );
 250
 251   if( !ok1 || !ok2 )
 252   {
 253     ok = false;
 254     return 0.0;
 255   }
 256
 257   bool start_pos = start_val>=val, fin_pos = fin_val>=val;
 258   ok = true;
 259
 260   while( fin-start>eps )
 261   {
 262     double mid = ( start+fin )/2.0, mid_val;
 263     ok = f.value( mid, mid_val );
 264     if( !ok )
 265       return 0.0;
 266
 267     //char buf[1024];
 268     //sprintf( buf, "start=%f\nfin=%f\nmid_val=%f\n", float( start ), float( fin ), float( mid_val ) );
 269     //MESSAGE( buf );
 270
 271     bool mid_pos = mid_val>=val;
 272     if( start_pos!=mid_pos )
 273     {
 274       fin_pos = mid_pos;
 275       fin = mid;
 276     }
 277     else if( fin_pos!=mid_pos )
 278     {
 279       start_pos = mid_pos;
 280       start = mid;
 281     }
 282     else
 283     {
 284       ok = false;
 285       break;
 286     }
 287   }
 288   return (start+fin)/2.0;
 289 }
 290
 291 bool buildDistribution( const TCollection_AsciiString& f, const int conv, const double start, const double end,
 292                         const smIdType nbSeg, vector<double>& data, const double eps )
 293 {
 294   FunctionExpr F( f.ToCString(), conv );
 295   return buildDistribution( F, start, end, nbSeg, data, eps );
 296 }
 297
 298 bool buildDistribution( const std::vector<double>& f, const int conv, const double start, const double end,
 299                         const smIdType nbSeg, vector<double>& data, const double eps )
 300 {
 301   FunctionTable F( f, conv );
 302   return buildDistribution( F, start, end, nbSeg, data, eps );
 303 }
 304
 305 bool buildDistribution( const Function& func, const double start, const double end,
 306                         const smIdType nbSeg, vector<double>& data, const double eps )
 307 {
 308   if( nbSeg<=0 )
 309     return false;
 310
 311   data.resize( nbSeg+1 );
 312   data[0] = start;
 313   double J = func.integral( start, end ) / double( nbSeg );
 314   if( J<1E-10 )
 315     return false;
 316
 317   bool ok;
 318   //MESSAGE( "distribution:" );
 319   //char buf[1024];
 320   for( smIdType i = 1; i < nbSeg; i++ )
 321   {
 322     FunctionIntegral f_int( &func, data[i-1] );
 323     data[i] = dihotomySolve( f_int, J, data[i-1], end, eps, ok );
 324     //sprintf( buf, "%f\n", float( data[i] ) );
 325     //MESSAGE( buf );
 326     if( !ok )
 327       return false;
 328   }
 329
 330   data[nbSeg] = end;
 331   return true;
 332 }
 333 }