src/StdMeshers/StdMeshers_CartesianParameters3D.cxx

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