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