src/StdMeshers/StdMeshers_CartesianParameters3D.cxx

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