src/StdMeshers/StdMeshers_CartesianParameters3D.cxx

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