src/StdMeshers/StdMeshers_AutomaticLength.cxx

   1 // Copyright (C) 2007-2012  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.
  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 SMESH : implementaion of SMESH idl descriptions
  24 //  File   : StdMeshers_AutomaticLength.cxx
  25 //  Author : Edward AGAPOV, OCC
  26 //  Module : SMESH
  27
  28 #include "StdMeshers_AutomaticLength.hxx"
  29
  30 #include "SMESH_Mesh.hxx"
  31 #include "SMESHDS_Mesh.hxx"
  32 #include "SMESH_Algo.hxx"
  33 #include "SMESHDS_SubMesh.hxx"
  34
  35 #include "utilities.h"
  36
  37 #include <TopTools_IndexedMapOfShape.hxx>
  38 #include <TopExp.hxx>
  39 #include <TopoDS.hxx>
  40 #include <TopoDS_Edge.hxx>
  41
  42 using namespace std;
  43
  44 //=============================================================================
  45 /*!
  46  *
  47  */
  48 //=============================================================================
  49
  50 StdMeshers_AutomaticLength::StdMeshers_AutomaticLength(int hypId, int studyId, SMESH_Gen * gen)
  51   :SMESH_Hypothesis(hypId, studyId, gen)
  52 {
  53   _name = "AutomaticLength";
  54   _param_algo_dim = 1; // is used by SMESH_Regular_1D
  55
  56   _mesh = 0;
  57   _fineness = 0;
  58 }
  59
  60 //=============================================================================
  61 /*!
  62  *
  63  */
  64 //=============================================================================
  65
  66 StdMeshers_AutomaticLength::~StdMeshers_AutomaticLength()
  67 {
  68 }
  69
  70 //================================================================================
  71 /*!
  72  * \brief Set Fineness
  73  * \param theFineness - The Fineness value [0.0-1.0],
  74  *                        0 - coarse mesh
  75  *                        1 - fine mesh
  76  *
  77  * Raise if theFineness is out of range
  78  * The "Initial Number of Elements on the Shortest Edge" (S0)
  79  * is divided by (0.5 + 4.5 x theFineness)
  80  */
  81 //================================================================================
  82
  83 const double theCoarseConst = 0.5;
  84 const double theFineConst   = 4.5;
  85
  86 void StdMeshers_AutomaticLength::SetFineness(double theFineness)
  87   throw(SALOME_Exception)
  88 {
  89   if ( theFineness < 0.0 || theFineness > 1.0 )
  90     throw SALOME_Exception(LOCALIZED("theFineness is out of range [0.0-1.0]"));
  91
  92   if ( _fineness != theFineness )
  93   {
  94     NotifySubMeshesHypothesisModification();
  95     _fineness = theFineness;
  96   }
  97 }
  98
  99 namespace {
 100
 101   //================================================================================
 102   /*!
 103    * \brief Return pointer to TopoDS_TShape
 104    * \param theShape - The TopoDS_Shape
 105    * \retval inline const TopoDS_TShape* - result
 106    */
 107   //================================================================================
 108
 109   inline const TopoDS_TShape* getTShape(const TopoDS_Shape& theShape)
 110   {
 111     return theShape.TShape().operator->();
 112   }
 113
 114   //================================================================================
 115   /*!
 116    * \brief computes segment length by S0 and edge length
 117    */
 118   //================================================================================
 119
 120   //const double a14divPI = 14. / M_PI;
 121   const double a2div7divPI = 2. / 7. / M_PI;
 122
 123   inline double segLength(double S0, double edgeLen, double minLen )
 124   {
 125     // PAL10237
 126     // S = S0 * f(L/Lmin) where
 127     // f(x) = 1 + (7 * 2/Pi * atan(x/5))
 128     // =>
 129     // S = S0 * ( 1 + 14/PI * atan( L / ( 5 * Lmin )))
 130     //
 131     // return S0 * ( 1. + a14divPI * atan( edgeLen / ( 5 * minLen )));
 132
 133     // The above formular gives too short segments when Lmax/Lmin is too high
 134     // because by this formular the largest segment is only 8 times longer than the
 135     // shortest one ( 2/Pi * atan(x/5) varies within [0,1] ). So a new formular is:
 136     //
 137     // f(x) = 1 + (x/7 * 2/Pi * atan(x/5))
 138     // =>
 139     // S = S0 * ( 1 + 2/7/PI * L/Lmin * atan( 5 * L/Lmin ))
 140     //
 141     const double Lratio = edgeLen / minLen;
 142     return S0 * ( 1. + a2div7divPI * Lratio * atan( 5 * Lratio ));
 143   }
 144
 145   //================================================================================
 146   /*!
 147    * \brief Compute segment length for all edges
 148    * \param theMesh - The mesh
 149    * \param theTShapeToLengthMap - The map of edge to segment length
 150    */
 151   //================================================================================
 152
 153   void computeLengths( SMESHDS_Mesh*                       aMesh,
 154                        map<const TopoDS_TShape*, double> & theTShapeToLengthMap,
 155                        double &                            theS0,
 156                        double &                            theMinLen)
 157   {
 158     theTShapeToLengthMap.clear();
 159
 160     TopoDS_Shape aMainShape = aMesh->ShapeToMesh();
 161
 162     // Find length of longest and shortest edge
 163     double Lmin = DBL_MAX, Lmax = -DBL_MAX;
 164     TopTools_IndexedMapOfShape edgeMap;
 165     TopExp::MapShapes( aMainShape, TopAbs_EDGE, edgeMap);
 166     for ( int i = 1; i <= edgeMap.Extent(); ++i )
 167     {
 168       TopoDS_Edge edge = TopoDS::Edge( edgeMap(i) );
 169       //if ( BRep_Tool::Degenerated( edge )) continue;
 170
 171       Standard_Real L = SMESH_Algo::EdgeLength( edge );
 172       if ( L < DBL_MIN ) continue;
 173
 174       if ( L > Lmax ) Lmax = L;
 175       if ( L < Lmin ) Lmin = L;
 176
 177       // remember i-th edge length
 178       theTShapeToLengthMap.insert( make_pair( getTShape( edge ), L ));
 179     }
 180
 181     // Compute S0 - minimal segement length, is computed by the shortest EDGE
 182
 183     /* image attached to PAL10237
 184
 185        NbSeg (on the shortest EDGE)
 186          ^
 187          |
 188        10|\
 189          | \
 190          |  \
 191          |   \
 192         5|    --------
 193          |
 194          +------------>
 195          1    10       Lmax/Lmin
 196     */
 197     const int NbSegMin = 5, NbSegMax = 10; //  on axis NbSeg
 198     const double Lrat1 = 1., Lrat2 = 10.;  //  on axis Lmax/Lmin
 199
 200     double Lratio = Lmax/Lmin;
 201     double NbSeg = NbSegMin;
 202     if ( Lratio < Lrat2 )
 203       NbSeg += ( Lrat2 - Lratio ) / ( Lrat2 - Lrat1 )  * ( NbSegMax - NbSegMin );
 204
 205     double S0 = Lmin / (int) NbSeg;
 206     MESSAGE( "S0 = " << S0 << ", Lmin = " << Lmin << ", Nbseg = " << (int) NbSeg);
 207
 208     // Compute segments length for all edges
 209
 210     map<const TopoDS_TShape*, double>::iterator tshape_length = theTShapeToLengthMap.begin();
 211     for ( ; tshape_length != theTShapeToLengthMap.end(); ++tshape_length )
 212     {
 213       double & L = tshape_length->second;
 214       L = segLength( S0, L, Lmin );
 215     }
 216     theS0 = S0;
 217     theMinLen = Lmin;
 218   }
 219 }
 220
 221 //=============================================================================
 222 /*!
 223  * \brief Computes segment length for an edge of given length
 224  */
 225 //=============================================================================
 226
 227 double StdMeshers_AutomaticLength::GetLength(const SMESH_Mesh* theMesh,
 228                                              const double      theEdgeLength)
 229   throw(SALOME_Exception)
 230 {
 231   if ( !theMesh ) throw SALOME_Exception(LOCALIZED("NULL Mesh"));
 232
 233   SMESHDS_Mesh* aMeshDS = const_cast< SMESH_Mesh* > ( theMesh )->GetMeshDS();
 234   if ( theMesh != _mesh )
 235   {
 236     computeLengths( aMeshDS, _TShapeToLength, _S0, _minLen );
 237     _mesh = theMesh;
 238   }
 239   double L = segLength( _S0, theEdgeLength, _minLen );
 240   return L / (theCoarseConst + theFineConst * _fineness);
 241 }
 242
 243 //=============================================================================
 244 /*!
 245  *
 246  */
 247 //=============================================================================
 248
 249 double StdMeshers_AutomaticLength::GetLength(const SMESH_Mesh*   theMesh,
 250                                              const TopoDS_Shape& anEdge)
 251   throw(SALOME_Exception)
 252 {
 253   if ( !theMesh ) throw SALOME_Exception(LOCALIZED("NULL Mesh"));
 254
 255   if ( anEdge.IsNull() || anEdge.ShapeType() != TopAbs_EDGE )
 256     throw SALOME_Exception(LOCALIZED("Bad edge shape"));
 257
 258   if ( theMesh != _mesh )
 259   {
 260     SMESHDS_Mesh* aMeshDS = const_cast< SMESH_Mesh* > ( theMesh )->GetMeshDS();
 261     computeLengths( aMeshDS, _TShapeToLength, _S0, _minLen );
 262     _mesh = theMesh;
 263   }
 264
 265   map<const TopoDS_TShape*, double>::iterator tshape_length =
 266     _TShapeToLength.find( getTShape( anEdge ));
 267
 268   if ( tshape_length == _TShapeToLength.end() )
 269     return 1; // it is a dgenerated edge
 270
 271   return tshape_length->second / (theCoarseConst + theFineConst * _fineness);
 272 }
 273
 274 //=============================================================================
 275 /*!
 276  *
 277  */
 278 //=============================================================================
 279
 280 ostream & StdMeshers_AutomaticLength::SaveTo(ostream & save)
 281 {
 282   save << _fineness;
 283   return save;
 284 }
 285
 286 //=============================================================================
 287 /*!
 288  *
 289  */
 290 //=============================================================================
 291
 292 istream & StdMeshers_AutomaticLength::LoadFrom(istream & load)
 293 {
 294   if ( ! ( load >> _fineness ))
 295     load.clear(ios::badbit | load.rdstate());
 296   return load;
 297 }
 298
 299 //=============================================================================
 300 /*!
 301  *
 302  */
 303 //=============================================================================
 304
 305 ostream & operator <<(ostream & save, StdMeshers_AutomaticLength & hyp)
 306 {
 307   return hyp.SaveTo( save );
 308 }
 309
 310 //=============================================================================
 311 /*!
 312  *
 313  */
 314 //=============================================================================
 315
 316 istream & operator >>(istream & load, StdMeshers_AutomaticLength & hyp)
 317 {
 318   return hyp.LoadFrom( load );
 319 }
 320
 321 //================================================================================
 322 /*!
 323  * \brief Initialize Fineness by the mesh built on the geometry
 324  * \param theMesh - the built mesh
 325  * \param theShape - the geometry of interest
 326  * \retval bool - true if parameter values have been successfully defined
 327  */
 328 //================================================================================
 329
 330 bool StdMeshers_AutomaticLength::SetParametersByMesh(const SMESH_Mesh*   theMesh,
 331                                                      const TopoDS_Shape& theShape)
 332 {
 333   if ( !theMesh || theShape.IsNull() )
 334     return false;
 335
 336   _fineness = 0;
 337
 338   SMESHDS_Mesh* aMeshDS = const_cast< SMESH_Mesh* >( theMesh )->GetMeshDS();
 339
 340   int nbEdges = 0;
 341   TopTools_IndexedMapOfShape edgeMap;
 342   TopExp::MapShapes( theShape, TopAbs_EDGE, edgeMap );
 343   for ( int i = 1; i <= edgeMap.Extent(); ++i )
 344   {
 345     const TopoDS_Edge& edge = TopoDS::Edge( edgeMap( i ));
 346
 347     // assure the base automatic length is stored in _TShapeToLength
 348     if ( i == 1 )
 349       GetLength( theMesh, edge );
 350
 351     // get current segment length
 352     double L = SMESH_Algo::EdgeLength( edge );
 353     if ( L <= DBL_MIN )
 354       continue;
 355     SMESHDS_SubMesh * eSubMesh = aMeshDS->MeshElements( edge );
 356     if ( !eSubMesh )
 357       return false;
 358     int nbSeg = eSubMesh->NbElements();
 359     if ( nbSeg < 1 )
 360       continue;
 361     double segLen = L / nbSeg;
 362
 363     // get segment length from _TShapeToLength
 364     map<const TopoDS_TShape*, double>::iterator tshape_length =
 365       _TShapeToLength.find( getTShape( edge ));
 366     if ( tshape_length == _TShapeToLength.end() )
 367       continue;
 368     double autoLen = tshape_length->second;
 369
 370     // segLen = autoLen / (theCoarseConst + theFineConst * _fineness) -->
 371     _fineness += ( autoLen / segLen - theCoarseConst ) / theFineConst;
 372
 373     ++nbEdges;
 374   }
 375   if ( nbEdges )
 376     _fineness /= nbEdges;
 377
 378   if (_fineness > 1.0)
 379     _fineness = 1.0;
 380   else if (_fineness < 0.0)
 381     _fineness = 0.0;
 382
 383   return nbEdges;
 384 }
 385
 386 //================================================================================
 387 /*!
 388  * \brief Initialize my parameter values by default parameters.
 389  *  \retval bool - true if parameter values have been successfully defined
 390  */
 391 //================================================================================
 392
 393 bool StdMeshers_AutomaticLength::SetParametersByDefaults(const TDefaults&  /*theDflts*/,
 394                                                          const SMESH_Mesh* /*theMesh*/)
 395 {
 396   return false;
 397
 398   // assure the base automatic length is stored in _TShapeToLength
 399 //   GetLength( theMesh, elemLenght );
 400
 401 //   // find maximal edge length
 402 //   double maxLen = 0;
 403 //   map<const TopoDS_TShape*, double>::iterator
 404 //     tshape_length = _TShapeToLength.begin(), slEnd = _TShapeToLength.end();
 405 //   for ( ; tshape_length != slEnd; ++tshape_length )
 406 //     if ( tshape_length->second > maxLen )
 407 //       maxLen = tshape_length->second;
 408
 409 //   // automatic length for longest element
 410 //   double autoLen = GetLength( theMesh, maxLen );
 411
 412 //   // elemLenght = autoLen / (theCoarseConst + theFineConst * _fineness) -->
 413 //   _fineness = ( autoLen / elemLenght - theCoarseConst ) / theFineConst;
 414
 415 //   return true;
 416 }