src/StdMeshers/StdMeshers_CartesianParameters3D.cxx

   1 // Copyright (C) 2007-2013  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 //  File   : StdMeshers_CartesianParameters3D.cxx
  24 //  Author : Edward AGAPOV
  25 //  Module : SMESH
  26 //
  27 #include "StdMeshers_CartesianParameters3D.hxx"
  28
  29 #include "StdMeshers_NumberOfSegments.hxx"
  30 #include "StdMeshers_Distribution.hxx"
  31 #include "SMESH_Gen.hxx"
  32
  33 #include "utilities.h"
  34
  35 #include <map>
  36 #include <limits>
  37
  38 #include <BRepGProp.hxx>
  39 #include <BRep_Tool.hxx>
  40 #include <Bnd_Box.hxx>
  41 #include <GProp_GProps.hxx>
  42 #include <GeomLib_IsPlanarSurface.hxx>
  43 #include <Geom_Surface.hxx>
  44 #include <Precision.hxx>
  45 #include <TopExp_Explorer.hxx>
  46 #include <TopLoc_Location.hxx>
  47 #include <TopTools_MapIteratorOfMapOfShape.hxx>
  48 #include <TopTools_MapOfShape.hxx>
  49 #include <TopoDS.hxx>
  50 #include <TopoDS_Face.hxx>
  51 #include <gp_Dir.hxx>
  52 #include <gp_Pln.hxx>
  53 #include <gp_Vec.hxx>
  54
  55 using namespace std;
  56
  57 //=======================================================================
  58 //function : StdMeshers_CartesianParameters3D
  59 //purpose  : Constructor
  60 //=======================================================================
  61
  62 StdMeshers_CartesianParameters3D::StdMeshers_CartesianParameters3D(int         hypId,
  63                                                                    int         studyId,
  64                                                                    SMESH_Gen * gen)
  65   : SMESH_Hypothesis(hypId, studyId, gen),
  66     _sizeThreshold( 4.0 ), // default according to the customer specification
  67     _toAddEdges( false )
  68 {
  69   _name = "CartesianParameters3D"; // used by "Cartesian_3D"
  70   _param_algo_dim = 3; // 3D
  71
  72   _axisDirs[0] = 1.;
  73   _axisDirs[1] = 0.;
  74   _axisDirs[2] = 0.;
  75
  76   _axisDirs[3] = 0.;
  77   _axisDirs[4] = 1.;
  78   _axisDirs[5] = 0.;
  79
  80   _axisDirs[6] = 0.;
  81   _axisDirs[7] = 0.;
  82   _axisDirs[8] = 1.;
  83
  84   _fixedPoint[0] = 0.;
  85   _fixedPoint[1] = 0.;
  86   _fixedPoint[2] = 0.;
  87   SetFixedPoint( _fixedPoint, /*toUnset=*/true );
  88 }
  89
  90
  91 namespace
  92 {
  93   const char* axisName[3] = { "X", "Y", "Z" };
  94
  95   typedef std::pair< double, std::pair< double, double > > TCooTriple;
  96
  97 #define gpXYZ( cTriple ) gp_XYZ( (cTriple).first, (cTriple).second.first, (cTriple).second.second )
  98
  99   //================================================================================
 100   /*!
 101    * \brief Compare two normals
 102    */
 103   //================================================================================
 104
 105   bool sameDir( const TCooTriple& n1, const TCooTriple& n2 )
 106   {
 107     gp_XYZ xyz1 = gpXYZ( n1 ), xyz2 = gpXYZ( n2 );
 108     return ( xyz1 - xyz2 ).Modulus() < 0.01;
 109   }
 110
 111   //================================================================================
 112   /*!
 113    * \brief Checks validity of an axis index, throws in case of invalidity
 114    */
 115   //================================================================================
 116
 117   void checkAxis(const int axis)
 118   {
 119     if ( axis < 0 || axis > 2 )
 120       throw SALOME_Exception(SMESH_Comment("Invalid axis index ") << axis <<
 121                              ". Valid axis indices are 0, 1 and 2");
 122   }
 123
 124   //================================================================================
 125   /*!
 126    * \brief Checks validity of spacing data, throws in case of invalidity
 127    */
 128   //================================================================================
 129
 130   void checkGridSpacing(std::vector<std::string>& spaceFunctions,
 131                         std::vector<double>&      internalPoints,
 132                         const std::string&        axis)
 133     throw ( SALOME_Exception )
 134   {
 135     if ( spaceFunctions.empty() )
 136       throw SALOME_Exception(SMESH_Comment("Empty space function for ") << axis );
 137
 138     for ( size_t i = 1; i < internalPoints.size(); ++i )
 139       if ( internalPoints[i] - internalPoints[i-1] < 0 )
 140         throw SALOME_Exception(SMESH_Comment("Wrong order of internal points along ") << axis);
 141       else if ( internalPoints[i] - internalPoints[i-1] < 1e-3 )
 142         throw SALOME_Exception(SMESH_Comment("Too close internal points along ") << axis );
 143
 144     const double tol = Precision::Confusion();
 145     if ( !internalPoints.empty() &&
 146          ( internalPoints.front() < -tol || internalPoints.back() > 1 + tol ))
 147       throw SALOME_Exception(SMESH_Comment("Invalid internal points along ") << axis);
 148
 149     if ( internalPoints.empty() || internalPoints.front() > tol )
 150       internalPoints.insert( internalPoints.begin(), 0. );
 151     if ( internalPoints.size() < 2 || internalPoints.back() < 1 - tol )
 152       internalPoints.push_back( 1. );
 153
 154     if ( internalPoints.size() != spaceFunctions.size() + 1 )
 155       throw SALOME_Exception
 156         (SMESH_Comment("Numbre of internal points mismatch number of functions for ") << axis);
 157
 158     for ( size_t i = 0; i < spaceFunctions.size(); ++i )
 159       spaceFunctions[i] =
 160         StdMeshers_NumberOfSegments::CheckExpressionFunction( spaceFunctions[i], -1 );
 161   }
 162 }
 163
 164 //=======================================================================
 165 //function : SetGrid
 166 //purpose  : Sets coordinates of node positions along an axes
 167 //=======================================================================
 168
 169 void StdMeshers_CartesianParameters3D::SetGrid(std::vector<double>& coords, int axis)
 170   throw ( SALOME_Exception )
 171 {
 172   checkAxis( axis );
 173
 174   if ( coords.size() < 2 )
 175     throw SALOME_Exception(LOCALIZED("Wrong number of grid coordinates"));
 176
 177   std::sort( coords.begin(), coords.end() );
 178
 179   bool changed = ( _coords[axis] != coords );
 180   if ( changed )
 181   {
 182     _coords[axis] = coords;
 183     NotifySubMeshesHypothesisModification();
 184   }
 185
 186   _spaceFunctions[axis].clear();
 187   _internalPoints[axis].clear();
 188 }
 189
 190 //=======================================================================
 191 //function : SetGridSpacing
 192 //purpose  : Set grid spacing along the three axes
 193 //=======================================================================
 194
 195 void StdMeshers_CartesianParameters3D::SetGridSpacing(std::vector<string>& xSpaceFuns,
 196                                                       std::vector<double>& xInternalPoints,
 197                                                       const int            axis)
 198   throw ( SALOME_Exception )
 199 {
 200   checkAxis( axis );
 201
 202   checkGridSpacing( xSpaceFuns, xInternalPoints, axisName[axis] );
 203
 204   bool changed = ( xSpaceFuns      != _spaceFunctions[axis] ||
 205                    xInternalPoints != _internalPoints[axis] );
 206
 207   _spaceFunctions[axis] = xSpaceFuns;
 208   _internalPoints[axis] = xInternalPoints;
 209   _coords[axis].clear();
 210
 211   if ( changed )
 212     NotifySubMeshesHypothesisModification();
 213 }
 214
 215 //=======================================================================
 216 //function : SetFixedPoint
 217 //purpose  : * Set/unset a fixed point, at which a node will be created provided that grid
 218 //           * is defined by spacing in all directions
 219 //=======================================================================
 220
 221 void StdMeshers_CartesianParameters3D::SetFixedPoint(const double p[3], bool toUnset)
 222 {
 223   if ( toUnset != Precision::IsInfinite( _fixedPoint[0] ))
 224     NotifySubMeshesHypothesisModification();
 225
 226   if ( toUnset )
 227     _fixedPoint[0] = Precision::Infinite();
 228   else
 229     std::copy( &p[0], &p[0]+3, &_fixedPoint[0] );
 230 }
 231
 232 //=======================================================================
 233 //function : GetFixedPoint
 234 //purpose  : Returns either false or (true + point coordinates)
 235 //=======================================================================
 236
 237 bool StdMeshers_CartesianParameters3D::GetFixedPoint(double p[3])
 238 {
 239   if ( Precision::IsInfinite( _fixedPoint[0] ))
 240     return false;
 241   std::copy( &_fixedPoint[0], &_fixedPoint[0]+3, &p[0] );
 242 }
 243
 244
 245 //=======================================================================
 246 //function : SetSizeThreshold
 247 //purpose  : Set size threshold
 248 //=======================================================================
 249
 250 void StdMeshers_CartesianParameters3D::SetSizeThreshold(const double threshold)
 251   throw ( SALOME_Exception )
 252 {
 253   if ( threshold <= 1.0 )
 254     throw SALOME_Exception(LOCALIZED("threshold must be > 1.0"));
 255
 256   bool changed = fabs( _sizeThreshold - threshold ) > 1e-6;
 257   _sizeThreshold = threshold;
 258
 259   if ( changed )
 260     NotifySubMeshesHypothesisModification();
 261 }
 262
 263 //=======================================================================
 264 //function : GetGridSpacing
 265 //purpose  : return spacing
 266 //=======================================================================
 267
 268 void StdMeshers_CartesianParameters3D::GetGridSpacing(std::vector<std::string>& spaceFunctions,
 269                                                       std::vector<double>&      internalPoints,
 270                                                       const int                 axis) const
 271   throw ( SALOME_Exception )
 272 {
 273   if ( !IsGridBySpacing(axis) )
 274     throw SALOME_Exception(LOCALIZED("The grid is defined by coordinates and not by spacing"));
 275
 276   spaceFunctions = _spaceFunctions[axis];
 277   internalPoints = _internalPoints[axis];
 278 }
 279
 280 //=======================================================================
 281 //function : IsGridBySpacing
 282 //=======================================================================
 283
 284 bool StdMeshers_CartesianParameters3D::IsGridBySpacing(const int axis) const
 285   throw ( SALOME_Exception )
 286 {
 287   checkAxis(axis);
 288   return !_spaceFunctions[axis].empty();
 289 }
 290
 291
 292 //=======================================================================
 293 //function : ComputeCoordinates
 294 //purpose  : Computes node coordinates by spacing functions
 295 //=======================================================================
 296
 297 void StdMeshers_CartesianParameters3D::ComputeCoordinates(const double         x0,
 298                                                           const double         x1,
 299                                                           vector<std::string>& spaceFuns,
 300                                                           vector<double>&      points,
 301                                                           vector<double>&      coords,
 302                                                           const std::string&   axis )
 303   throw ( SALOME_Exception )
 304 {
 305   checkGridSpacing( spaceFuns, points, axis );
 306
 307   coords.clear();
 308   for ( size_t i = 0; i < spaceFuns.size(); ++i )
 309   {
 310     FunctionExpr fun( spaceFuns[i].c_str(), /*convMode=*/-1 );
 311
 312     const double p0 = x0 * ( 1. - points[i])   + x1 * points[i];
 313     const double p1 = x0 * ( 1. - points[i+1]) + x1 * points[i+1];
 314     const double length = p1 - p0;
 315
 316     const size_t nbSections = 1000;
 317     const double sectionLen = ( p1 - p0 ) / nbSections;
 318     vector< double > nbSegments( nbSections + 1 );
 319     nbSegments[ 0 ] = 0.;
 320
 321     double t, spacing = 0;
 322     for ( size_t i = 1; i <= nbSections; ++i )
 323     {
 324       t = double( i ) / nbSections;
 325       if ( !fun.value( t, spacing ) || spacing < std::numeric_limits<double>::min() )
 326         throw SALOME_Exception(LOCALIZED("Invalid spacing function"));
 327       nbSegments[ i ] = nbSegments[ i-1 ] + std::min( 1., sectionLen / spacing );
 328     }
 329
 330     const int nbCells = max (1, int(floor(nbSegments.back()+0.5)));
 331     const double corr = nbCells / nbSegments.back();
 332
 333     if ( coords.empty() ) coords.push_back( p0 );
 334
 335     for ( size_t iCell = 1, i = 1; i <= nbSections; ++i )
 336     {
 337       if ( nbSegments[i]*corr >= iCell )
 338       {
 339         t = (i - ( nbSegments[i] - iCell/corr )/( nbSegments[i] - nbSegments[i-1] )) / nbSections;
 340         coords.push_back( p0 + t * length );
 341         ++iCell;
 342       }
 343     }
 344     const double lastCellLen = coords.back() - coords[ coords.size() - 2 ];
 345     if ( fabs( coords.back() - p1 ) > 0.5 * lastCellLen )
 346       coords.push_back ( p1 );
 347   }
 348 }
 349
 350 //=======================================================================
 351 //function : GetCoordinates
 352 //purpose  : Return coordinates of node positions along the three axes.
 353 //           If the grid is defined by spacing functions, the coordinates are computed
 354 //=======================================================================
 355
 356 void StdMeshers_CartesianParameters3D::GetCoordinates(std::vector<double>& xNodes,
 357                                                       std::vector<double>& yNodes,
 358                                                       std::vector<double>& zNodes,
 359                                                       const Bnd_Box&       bndBox) const
 360   throw ( SALOME_Exception )
 361 {
 362   double x0,y0,z0, x1,y1,z1;
 363   if ( IsGridBySpacing(0) || IsGridBySpacing(1) || IsGridBySpacing(2))
 364   {
 365     if ( bndBox.IsVoid() ||
 366          bndBox.IsXThin( Precision::Confusion() ) ||
 367          bndBox.IsYThin( Precision::Confusion() ) ||
 368          bndBox.IsZThin( Precision::Confusion() ) )
 369       throw SALOME_Exception(LOCALIZED("Invalid bounding box"));
 370     bndBox.Get(x0,y0,z0, x1,y1,z1);
 371   }
 372
 373   StdMeshers_CartesianParameters3D* me = const_cast<StdMeshers_CartesianParameters3D*>(this);
 374   if ( IsGridBySpacing(0) )
 375     ComputeCoordinates( x0, x1, me->_spaceFunctions[0], me->_internalPoints[0], xNodes, "X" );
 376   else
 377     xNodes = _coords[0];
 378
 379   if ( IsGridBySpacing(1) )
 380     ComputeCoordinates( y0, y1, me->_spaceFunctions[1], me->_internalPoints[1], yNodes, "Y" );
 381   else
 382     yNodes = _coords[1];
 383
 384   if ( IsGridBySpacing(2) )
 385     ComputeCoordinates( z0, z1, me->_spaceFunctions[2], me->_internalPoints[2], zNodes, "Z" );
 386   else
 387     zNodes = _coords[2];
 388 }
 389
 390 //=======================================================================
 391 //function : ComputeOptimalAxesDirs
 392 //purpose  : Returns axes at which number of hexahedra is maximal
 393 //=======================================================================
 394
 395 void StdMeshers_CartesianParameters3D::
 396 ComputeOptimalAxesDirs(const TopoDS_Shape& shape,
 397                        const bool          isOrthogonal,
 398                        double              dirCoords[9])
 399 {
 400   for ( int i = 0; i < 9; ++i ) dirCoords[i] = 0.;
 401   dirCoords[0] = dirCoords[4] = dirCoords[8] = 1.;
 402
 403   if ( shape.IsNull() ) return;
 404
 405   TopLoc_Location loc;
 406   TopExp_Explorer exp;
 407
 408   // get external FACEs of the shape
 409   TopTools_MapOfShape faceMap;
 410   for ( exp.Init( shape, TopAbs_FACE ); exp.More(); exp.Next() )
 411     if ( !faceMap.Add( exp.Current() ))
 412       faceMap.Remove( exp.Current() );
 413
 414   // sort areas of planar faces by normal direction
 415
 416   std::multimap< TCooTriple, double > areasByNormal;
 417
 418   TopTools_MapIteratorOfMapOfShape fIt ( faceMap );
 419   for ( ; fIt.More(); fIt.Next() )
 420   {
 421     const TopoDS_Face&   face = TopoDS::Face( fIt.Key() );
 422     Handle(Geom_Surface) surf = BRep_Tool::Surface( face, loc );
 423     if ( surf.IsNull() ) continue;
 424
 425     GeomLib_IsPlanarSurface check( surf, 1e-5 );
 426     if ( !check.IsPlanar() ) continue;
 427
 428     GProp_GProps SProps;
 429     BRepGProp::SurfaceProperties( face, SProps );
 430     double area = SProps.Mass();
 431
 432     gp_Pln pln  = check.Plan();
 433     gp_Dir norm = pln.Axis().Direction().Transformed( loc );
 434     if ( norm.X() < -1e-3 ) { // negative X
 435       norm.Reverse();
 436     } else if ( norm.X() < 1e-3 ) { // zero X
 437       if ( norm.Y() < -1e-3 ) { // negative Y
 438         norm.Reverse();
 439       } else if ( norm.Y() < 1e-3 ) { // zero X && zero Y
 440         if ( norm.Y() < -1e-3 ) // negative Z
 441           norm.Reverse();
 442       }
 443     }
 444     TCooTriple coo3( norm.X(), make_pair( norm.Y(), norm.Z() ));
 445     areasByNormal.insert( make_pair( coo3, area ));
 446   }
 447
 448   // group coplanar normals and sort groups by sum area
 449
 450   std::multimap< double, vector< const TCooTriple* > > normsByArea;
 451   std::multimap< TCooTriple, double >::iterator norm2a = areasByNormal.begin();
 452   const TCooTriple*           norm1 = 0;
 453   double                      sumArea = 0;
 454   vector< const TCooTriple* > norms;
 455   for ( int iF = 1; norm2a != areasByNormal.end(); ++norm2a, ++iF )
 456   {
 457
 458     if ( !norm1 || !sameDir( *norm1, norm2a->first ))
 459     {
 460       if ( !norms.empty() )
 461       {
 462         normsByArea.insert( make_pair( sumArea, norms ));
 463         norms.clear();
 464       }
 465       norm1   = & norm2a->first;
 466       sumArea = norm2a->second;
 467       norms.push_back( norm1 );
 468     }
 469     else
 470     {
 471       sumArea += norm2a->second;
 472       norms.push_back( & norm2a->first );
 473     }
 474     if ( iF == areasByNormal.size() )
 475       normsByArea.insert( make_pair( sumArea, norms ));
 476   }
 477
 478   // try to set dirs by planar faces
 479
 480   gp_XYZ normDirs[3]; // normals to largest planes
 481
 482   if ( !normsByArea.empty() )
 483   {
 484     norm1 = normsByArea.rbegin()->second[0];
 485     normDirs[0] = gpXYZ( *norm1 );
 486
 487     if ( normsByArea.size() == 1 )
 488     {
 489       normDirs[1] = normDirs[0];
 490       if ( Abs( normDirs[0].Y() ) < 1e-100 &&
 491            Abs( normDirs[0].Z() ) < 1e-100 ) // normDirs[0] || OX
 492         normDirs[1].SetY( normDirs[0].Y() + 1. );
 493       else
 494         normDirs[1].SetX( normDirs[0].X() + 1. );
 495     }
 496     else
 497     {
 498       // look for 2 other directions
 499       gp_XYZ testDir = normDirs[0], minDir, maxDir;
 500       for ( int is2nd = 0; is2nd < 2; ++is2nd )
 501       {
 502         double maxMetric = 0, minMetric = 1e100;
 503         std::multimap< double, vector< const TCooTriple* > >::iterator a2n;
 504         for ( a2n = normsByArea.begin(); a2n != normsByArea.end(); ++a2n )
 505         {
 506           gp_XYZ n = gpXYZ( *( a2n->second[0]) );
 507           double dot = Abs( n * testDir );
 508           double metric = ( 1. - dot ) * ( isOrthogonal ? 1 : a2n->first );
 509           if ( metric > maxMetric )
 510           {
 511             maxDir = n;
 512             maxMetric = metric;
 513           }
 514           if ( metric < minMetric )
 515           {
 516             minDir = n;
 517             minMetric = metric;
 518           }
 519         }
 520         if ( is2nd )
 521         {
 522           normDirs[2] = minDir;
 523         }
 524         else
 525         {
 526           normDirs[1] = maxDir;
 527           normDirs[2] = normDirs[0] ^ normDirs[1];
 528           if ( isOrthogonal || normsByArea.size() < 3 )
 529             break;
 530           testDir = normDirs[2];
 531         }
 532       }
 533     }
 534     if ( isOrthogonal || normsByArea.size() == 1 )
 535     {
 536       normDirs[2] = normDirs[0] ^ normDirs[1];
 537       normDirs[1] = normDirs[2] ^ normDirs[0];
 538     }
 539   }
 540   else
 541   {
 542     return;
 543   }
 544
 545   gp_XYZ dirs[3];
 546   dirs[0] = normDirs[0] ^ normDirs[1];
 547   dirs[1] = normDirs[1] ^ normDirs[2];
 548   dirs[2] = normDirs[2] ^ normDirs[0];
 549
 550   dirs[0].Normalize();
 551   dirs[1].Normalize();
 552   dirs[2].Normalize();
 553
 554   // Select dirs for X, Y and Z axes
 555   int iX = ( Abs( dirs[0].X() ) > Abs( dirs[1].X() )) ? 0 : 1;
 556   if ( Abs( dirs[iX].X() ) < Abs( dirs[2].X() ))
 557     iX = 2;
 558   int iY = ( iX == 0 ) ? 1 : (( Abs( dirs[0].Y() ) > Abs( dirs[1].Y() )) ? 0 : 1 );
 559   if ( Abs( dirs[iY].Y() ) < Abs( dirs[2].Y() ) && iX != 2 )
 560     iY = 2;
 561   int iZ = 3 - iX - iY;
 562
 563   if ( dirs[iX].X() < 0 ) dirs[iX].Reverse();
 564   if ( dirs[iY].Y() < 0 ) dirs[iY].Reverse();
 565   gp_XYZ zDir = dirs[iX] ^ dirs[iY];
 566   if ( dirs[iZ] * zDir < 0 )
 567     dirs[iZ].Reverse();
 568
 569   dirCoords[0] = dirs[iX].X();
 570   dirCoords[1] = dirs[iX].Y();
 571   dirCoords[2] = dirs[iX].Z();
 572   dirCoords[3] = dirs[iY].X();
 573   dirCoords[4] = dirs[iY].Y();
 574   dirCoords[5] = dirs[iY].Z();
 575   dirCoords[6] = dirs[iZ].X();
 576   dirCoords[7] = dirs[iZ].Y();
 577   dirCoords[8] = dirs[iZ].Z();
 578 }
 579
 580 //=======================================================================
 581 //function : SetAxisDirs
 582 //purpose  : Sets custom direction of axes
 583 //=======================================================================
 584
 585 void StdMeshers_CartesianParameters3D::SetAxisDirs(const double* the9DirComps)
 586   throw ( SALOME_Exception )
 587 {
 588   gp_Vec x( the9DirComps[0],
 589             the9DirComps[1],
 590             the9DirComps[2] );
 591   gp_Vec y( the9DirComps[3],
 592             the9DirComps[4],
 593             the9DirComps[5] );
 594   gp_Vec z( the9DirComps[6],
 595             the9DirComps[7],
 596             the9DirComps[8] );
 597   if ( x.Magnitude() < RealSmall() ||
 598        y.Magnitude() < RealSmall() ||
 599        z.Magnitude() < RealSmall() )
 600     throw SALOME_Exception("Zero magnitude of axis direction");
 601
 602   if ( x.IsParallel( y, M_PI / 180. ) ||
 603        x.IsParallel( z, M_PI / 180. ) ||
 604        y.IsParallel( z, M_PI / 180. ))
 605     throw SALOME_Exception("Parallel axis directions");
 606
 607   gp_Vec normXY = x ^ y, normYZ = y ^ z;
 608   if ( normXY.IsParallel( normYZ, M_PI / 180. ))
 609     throw SALOME_Exception("Axes lie in one plane");
 610
 611   bool isChanged = false;
 612   for ( int i = 0; i < 9; ++i )
 613   {
 614     if ( Abs( _axisDirs[i] - the9DirComps[i] ) > 1e-7 )
 615       isChanged = true;
 616     _axisDirs[i] = the9DirComps[i];
 617   }
 618   if ( isChanged )
 619     NotifySubMeshesHypothesisModification();
 620 }
 621
 622 //=======================================================================
 623 //function : GetGrid
 624 //purpose  : Return coordinates of node positions along the three axes
 625 //=======================================================================
 626
 627 void StdMeshers_CartesianParameters3D::GetGrid(std::vector<double>& coords, int axis) const
 628   throw ( SALOME_Exception )
 629 {
 630   if ( IsGridBySpacing(axis) )
 631     throw SALOME_Exception(LOCALIZED("The grid is defined by spacing and not by coordinates"));
 632
 633   coords = _coords[axis];
 634 }
 635
 636 //=======================================================================
 637 //function : GetSizeThreshold
 638 //purpose  : Return size threshold
 639 //=======================================================================
 640
 641 double StdMeshers_CartesianParameters3D::GetSizeThreshold() const
 642 {
 643   return _sizeThreshold;
 644 }
 645
 646 //=======================================================================
 647 //function : SetToAddEdges
 648 //purpose  : Enables implementation of geometrical edges into the mesh. If this feature
 649 //           is disabled, sharp edges of the shape are lost ("smoothed") in the mesh if
 650 //           they don't coincide with the grid lines
 651 //=======================================================================
 652
 653 void StdMeshers_CartesianParameters3D::SetToAddEdges(bool toAdd)
 654 {
 655   if ( _toAddEdges != toAdd )
 656   {
 657     _toAddEdges = toAdd;
 658     NotifySubMeshesHypothesisModification();
 659   }
 660 }
 661
 662 //=======================================================================
 663 //function : GetToAddEdges
 664 //purpose  : Returns true if implementation of geometrical edges into the
 665 //           mesh is enabled
 666 //=======================================================================
 667
 668 bool StdMeshers_CartesianParameters3D::GetToAddEdges() const
 669 {
 670   return _toAddEdges;
 671 }
 672
 673 //=======================================================================
 674 //function : IsDefined
 675 //purpose  : Return true if parameters are well defined
 676 //=======================================================================
 677
 678 bool StdMeshers_CartesianParameters3D::IsDefined() const
 679 {
 680   for ( int i = 0; i < 3; ++i )
 681     if (_coords[i].empty() && (_spaceFunctions[i].empty() || _internalPoints[i].empty()))
 682       return false;
 683
 684   return ( _sizeThreshold > 1.0 );
 685 }
 686
 687 //=======================================================================
 688 //function : SaveTo
 689 //purpose  : store my parameters into a stream
 690 //=======================================================================
 691
 692 std::ostream & StdMeshers_CartesianParameters3D::SaveTo(std::ostream & save)
 693 {
 694   save << _sizeThreshold << " ";
 695
 696   for ( int i = 0; i < 3; ++i )
 697   {
 698     save << _coords[i].size() << " ";
 699     for ( size_t j = 0; j < _coords[i].size(); ++j )
 700       save << _coords[i][j] << " ";
 701
 702     save << _internalPoints[i].size() << " ";
 703     for ( size_t j = 0; j < _internalPoints[i].size(); ++j )
 704       save << _internalPoints[i][j] << " ";
 705
 706     save << _spaceFunctions[i].size() << " ";
 707     for ( size_t j = 0; j < _spaceFunctions[i].size(); ++j )
 708       save << _spaceFunctions[i][j] << " ";
 709   }
 710   save << _toAddEdges << " ";
 711
 712   save.setf( save.scientific );
 713   save.precision( 12 );
 714   for ( int i = 0; i < 9; ++i )
 715     save << _axisDirs[i] << " ";
 716
 717   for ( int i = 0; i < 3; ++i )
 718     save << _fixedPoint[i] << " ";
 719
 720   return save;
 721 }
 722
 723 //=======================================================================
 724 //function : LoadFrom
 725 //purpose  : resore my parameters from a stream
 726 //=======================================================================
 727
 728 std::istream & StdMeshers_CartesianParameters3D::LoadFrom(std::istream & load)
 729 {
 730   bool ok;
 731
 732   ok = ( load >> _sizeThreshold );
 733   for ( int ax = 0; ax < 3; ++ax )
 734   {
 735     if (ok)
 736     {
 737       size_t i = 0;
 738       ok = (load >> i  );
 739       if ( i > 0 && ok )
 740       {
 741         _coords[ax].resize( i );
 742         for ( i = 0; i < _coords[ax].size() && ok; ++i )
 743           ok = (load >> _coords[ax][i]  );
 744       }
 745     }
 746     if (ok)
 747     {
 748       size_t i = 0;
 749       ok = (load >> i  );
 750       if ( i > 0 && ok )
 751       {
 752         _internalPoints[ax].resize( i );
 753         for ( i = 0; i < _internalPoints[ax].size() && ok; ++i )
 754           ok = (load >> _internalPoints[ax][i]  );
 755       }
 756     }
 757     if (ok)
 758     {
 759       size_t i = 0;
 760       ok = (load >> i  );
 761       if ( i > 0 && ok )
 762       {
 763         _spaceFunctions[ax].resize( i );
 764         for ( i = 0; i < _spaceFunctions[ax].size() && ok; ++i )
 765           ok = (load >> _spaceFunctions[ax][i]  );
 766       }
 767     }
 768   }
 769
 770   ok = ( load >> _toAddEdges );
 771
 772   for ( int i = 0; i < 9 && ok; ++i )
 773     ok = ( load >> _axisDirs[i]);
 774
 775   for ( int i = 0; i < 3 && ok ; ++i )
 776     ok = ( load >> _fixedPoint[i]);
 777
 778   return load;
 779 }
 780
 781 //=======================================================================
 782 //function : SetParametersByMesh
 783 //=======================================================================
 784
 785 bool StdMeshers_CartesianParameters3D::SetParametersByMesh(const SMESH_Mesh*   ,
 786                                                            const TopoDS_Shape& )
 787 {
 788   return false;
 789 }
 790
 791 //=======================================================================
 792 //function : SetParametersByDefaults
 793 //=======================================================================
 794
 795 bool StdMeshers_CartesianParameters3D::SetParametersByDefaults(const TDefaults&  dflts,
 796                                                                const SMESH_Mesh* /*theMesh*/)
 797 {
 798   if ( dflts._elemLength > 1e-100 )
 799   {
 800     vector<string> spacing( 1, SMESH_Comment(dflts._elemLength));
 801     vector<double> intPnts;
 802     SetGridSpacing( spacing, intPnts, 0 );
 803     SetGridSpacing( spacing, intPnts, 1 );
 804     SetGridSpacing( spacing, intPnts, 2 );
 805     return true;
 806   }
 807   return false;
 808 }
 809