Salome HOME
21bc0bbb180e147ae57cd93170045406d5d8247e
[modules/smesh.git] / 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 {
69   _name = "CartesianParameters3D"; // used by "Cartesian_3D"
70   _param_algo_dim = 3; // 3D
71
72   _axisDirs[0] = 1.;
73   _axisDirs[1] = 0.;
74   _axisDirs[2] = 0.;
75
76   _axisDirs[3] = 0.;
77   _axisDirs[4] = 1.;
78   _axisDirs[5] = 0.;
79
80   _axisDirs[6] = 0.;
81   _axisDirs[7] = 0.;
82   _axisDirs[8] = 1.;
83
84   _fixedPoint[0] = 0.;
85   _fixedPoint[1] = 0.;
86   _fixedPoint[2] = 0.;
87   SetFixedPoint( _fixedPoint, /*toUnset=*/true );
88 }
89
90
91 namespace
92 {
93   const char* axisName[3] = { "X", "Y", "Z" };
94
95   typedef std::pair< double, std::pair< double, double > > TCooTriple;
96
97 #define gpXYZ( cTriple ) gp_XYZ( (cTriple).first, (cTriple).second.first, (cTriple).second.second )
98
99   //================================================================================
100   /*!
101    * \brief Compare two normals
102    */
103   //================================================================================
104
105   bool sameDir( const TCooTriple& n1, const TCooTriple& n2 )
106   {
107     gp_XYZ xyz1 = gpXYZ( n1 ), xyz2 = gpXYZ( n2 );
108     return ( xyz1 - xyz2 ).Modulus() < 0.01;
109   }
110
111   //================================================================================
112   /*!
113    * \brief Checks validity of an axis index, throws in case of invalidity
114    */
115   //================================================================================
116
117   void checkAxis(const int axis)
118   {
119     if ( axis < 0 || axis > 2 )
120       throw SALOME_Exception(SMESH_Comment("Invalid axis index ") << axis <<
121                              ". Valid axis indices are 0, 1 and 2");
122   }
123
124   //================================================================================
125   /*!
126    * \brief Checks validity of spacing data, throws in case of invalidity
127    */
128   //================================================================================
129
130   void checkGridSpacing(std::vector<std::string>& spaceFunctions,
131                         std::vector<double>&      internalPoints,
132                         const std::string&        axis)
133     throw ( SALOME_Exception )
134   {
135     if ( spaceFunctions.empty() )
136       throw SALOME_Exception(SMESH_Comment("Empty space function for ") << axis );
137
138     for ( size_t i = 1; i < internalPoints.size(); ++i )
139       if ( internalPoints[i] - internalPoints[i-1] < 0 )
140         throw SALOME_Exception(SMESH_Comment("Wrong order of internal points along ") << axis);
141       else if ( internalPoints[i] - internalPoints[i-1] < 1e-3 )
142         throw SALOME_Exception(SMESH_Comment("Too close internal points along ") << axis );
143
144     const double tol = Precision::Confusion();
145     if ( !internalPoints.empty() &&
146          ( internalPoints.front() < -tol || internalPoints.back() > 1 + tol ))
147       throw SALOME_Exception(SMESH_Comment("Invalid internal points along ") << axis);
148
149     if ( internalPoints.empty() || internalPoints.front() > tol )
150       internalPoints.insert( internalPoints.begin(), 0. );
151     if ( internalPoints.size() < 2 || internalPoints.back() < 1 - tol )
152       internalPoints.push_back( 1. );
153
154     if ( internalPoints.size() != spaceFunctions.size() + 1 )
155       throw SALOME_Exception
156         (SMESH_Comment("Numbre of internal points mismatch number of functions for ") << axis);
157
158     for ( size_t i = 0; i < spaceFunctions.size(); ++i )
159       spaceFunctions[i] =
160         StdMeshers_NumberOfSegments::CheckExpressionFunction( spaceFunctions[i], -1 );
161   }
162 }
163
164 //=======================================================================
165 //function : SetGrid
166 //purpose  : Sets coordinates of node positions along an axes
167 //=======================================================================
168
169 void StdMeshers_CartesianParameters3D::SetGrid(std::vector<double>& coords, int axis)
170   throw ( SALOME_Exception )
171 {
172   checkAxis( axis );
173
174   if ( coords.size() < 2 )
175     throw SALOME_Exception(LOCALIZED("Wrong number of grid coordinates"));
176
177   std::sort( coords.begin(), coords.end() );
178
179   bool changed = ( _coords[axis] != coords );
180   if ( changed )
181   {
182     _coords[axis] = coords;
183     NotifySubMeshesHypothesisModification();
184   }
185
186   _spaceFunctions[axis].clear();
187   _internalPoints[axis].clear();
188 }
189
190 //=======================================================================
191 //function : SetGridSpacing
192 //purpose  : Set grid spacing along the three axes
193 //=======================================================================
194
195 void StdMeshers_CartesianParameters3D::SetGridSpacing(std::vector<string>& xSpaceFuns,
196                                                       std::vector<double>& xInternalPoints,
197                                                       const int            axis)
198   throw ( SALOME_Exception )
199 {
200   checkAxis( axis );
201
202   checkGridSpacing( xSpaceFuns, xInternalPoints, axisName[axis] );
203
204   bool changed = ( xSpaceFuns      != _spaceFunctions[axis] ||
205                    xInternalPoints != _internalPoints[axis] );
206
207   _spaceFunctions[axis] = xSpaceFuns;
208   _internalPoints[axis] = xInternalPoints;
209   _coords[axis].clear();
210
211   if ( changed )
212     NotifySubMeshesHypothesisModification();
213 }
214
215 //=======================================================================
216 //function : SetFixedPoint
217 //purpose  : * Set/unset a fixed point, at which a node will be created provided that grid
218 //           * is defined by spacing in all directions
219 //=======================================================================
220
221 void StdMeshers_CartesianParameters3D::SetFixedPoint(const double p[3], bool toUnset)
222 {
223   if ( toUnset != Precision::IsInfinite( _fixedPoint[0] ))
224     NotifySubMeshesHypothesisModification();
225
226   if ( toUnset )
227     _fixedPoint[0] = Precision::Infinite();
228   else
229     std::copy( &p[0], &p[0]+3, &_fixedPoint[0] );
230 }
231
232 //=======================================================================
233 //function : GetFixedPoint
234 //purpose  : Returns either false or (true + point coordinates)
235 //=======================================================================
236
237 bool StdMeshers_CartesianParameters3D::GetFixedPoint(double p[3]) 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   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     StdMeshers::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     size_t 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 ( size_t iF = 1; norm2a != areasByNormal.end(); ++norm2a, ++iF )
526   {
527     if ( !norm1 || !sameDir( *norm1, norm2a->first ))
528     {
529       if ( !norms.empty() )
530       {
531         normsByArea.insert( make_pair( sumArea, norms ));
532         norms.clear();
533       }
534       norm1   = & norm2a->first;
535       sumArea = norm2a->second;
536       norms.push_back( norm1 );
537     }
538     else
539     {
540       sumArea += norm2a->second;
541       norms.push_back( & norm2a->first );
542     }
543     if ( iF == areasByNormal.size() )
544       normsByArea.insert( make_pair( sumArea, norms ));
545   }
546
547   // try to set dirs by planar faces
548
549   gp_XYZ normDirs[3]; // normals to largest planes
550
551   if ( !normsByArea.empty() )
552   {
553     norm1 = normsByArea.rbegin()->second[0];
554     normDirs[0] = gpXYZ( *norm1 );
555
556     if ( normsByArea.size() == 1 )
557     {
558       normDirs[1] = normDirs[0];
559       if ( Abs( normDirs[0].Y() ) < 1e-100 &&
560            Abs( normDirs[0].Z() ) < 1e-100 ) // normDirs[0] || OX
561         normDirs[1].SetY( normDirs[0].Y() + 1. );
562       else
563         normDirs[1].SetX( normDirs[0].X() + 1. );
564     }
565     else
566     {
567       // look for 2 other directions
568       gp_XYZ testDir = normDirs[0], minDir, maxDir;
569       for ( int is2nd = 0; is2nd < 2; ++is2nd )
570       {
571         double maxMetric = 0, minMetric = 1e100;
572         std::multimap< double, vector< const TCooTriple* > >::iterator a2n;
573         for ( a2n = normsByArea.begin(); a2n != normsByArea.end(); ++a2n )
574         {
575           gp_XYZ n = gpXYZ( *( a2n->second[0]) );
576           double dot = Abs( n * testDir );
577           double metric = ( 1. - dot ) * ( isOrthogonal ? 1 : a2n->first );
578           if ( metric > maxMetric )
579           {
580             maxDir = n;
581             maxMetric = metric;
582           }
583           if ( metric < minMetric )
584           {
585             minDir = n;
586             minMetric = metric;
587           }
588         }
589         if ( is2nd )
590         {
591           normDirs[2] = minDir;
592         }
593         else
594         {
595           normDirs[1] = maxDir;
596           normDirs[2] = normDirs[0] ^ normDirs[1];
597           if ( isOrthogonal || normsByArea.size() < 3 )
598             break;
599           testDir = normDirs[2];
600         }
601       }
602     }
603     if ( isOrthogonal || normsByArea.size() == 1 )
604     {
605       normDirs[2] = normDirs[0] ^ normDirs[1];
606       normDirs[1] = normDirs[2] ^ normDirs[0];
607     }
608   }
609   else
610   {
611     return;
612   }
613
614   gp_XYZ dirs[3];
615   dirs[0] = normDirs[0] ^ normDirs[1];
616   dirs[1] = normDirs[1] ^ normDirs[2];
617   dirs[2] = normDirs[2] ^ normDirs[0];
618
619   dirs[0].Normalize();
620   dirs[1].Normalize();
621   dirs[2].Normalize();
622
623   // Select dirs for X, Y and Z axes
624   int iX = ( Abs( dirs[0].X() ) > Abs( dirs[1].X() )) ? 0 : 1;
625   if ( Abs( dirs[iX].X() ) < Abs( dirs[2].X() ))
626     iX = 2;
627   int iY = ( iX == 0 ) ? 1 : (( Abs( dirs[0].Y() ) > Abs( dirs[1].Y() )) ? 0 : 1 );
628   if ( Abs( dirs[iY].Y() ) < Abs( dirs[2].Y() ) && iX != 2 )
629     iY = 2;
630   int iZ = 3 - iX - iY;
631
632   if ( dirs[iX].X() < 0 ) dirs[iX].Reverse();
633   if ( dirs[iY].Y() < 0 ) dirs[iY].Reverse();
634   gp_XYZ zDir = dirs[iX] ^ dirs[iY];
635   if ( dirs[iZ] * zDir < 0 )
636     dirs[iZ].Reverse();
637
638   dirCoords[0] = dirs[iX].X();
639   dirCoords[1] = dirs[iX].Y();
640   dirCoords[2] = dirs[iX].Z();
641   dirCoords[3] = dirs[iY].X();
642   dirCoords[4] = dirs[iY].Y();
643   dirCoords[5] = dirs[iY].Z();
644   dirCoords[6] = dirs[iZ].X();
645   dirCoords[7] = dirs[iZ].Y();
646   dirCoords[8] = dirs[iZ].Z();
647 }
648
649 //=======================================================================
650 //function : SetAxisDirs
651 //purpose  : Sets custom direction of axes
652 //=======================================================================
653
654 void StdMeshers_CartesianParameters3D::SetAxisDirs(const double* the9DirComps)
655   throw ( SALOME_Exception )
656 {
657   gp_Vec x( the9DirComps[0],
658             the9DirComps[1],
659             the9DirComps[2] );
660   gp_Vec y( the9DirComps[3],
661             the9DirComps[4],
662             the9DirComps[5] );
663   gp_Vec z( the9DirComps[6],
664             the9DirComps[7],
665             the9DirComps[8] );
666   if ( x.Magnitude() < RealSmall() ||
667        y.Magnitude() < RealSmall() ||
668        z.Magnitude() < RealSmall() )
669     throw SALOME_Exception("Zero magnitude of axis direction");
670
671   if ( x.IsParallel( y, M_PI / 180. ) ||
672        x.IsParallel( z, M_PI / 180. ) ||
673        y.IsParallel( z, M_PI / 180. ))
674     throw SALOME_Exception("Parallel axis directions");
675
676   gp_Vec normXY = x ^ y, normYZ = y ^ z;
677   if ( normXY.IsParallel( normYZ, M_PI / 180. ))
678     throw SALOME_Exception("Axes lie in one plane");
679
680   bool isChanged = false;
681   for ( int i = 0; i < 9; ++i )
682   {
683     if ( Abs( _axisDirs[i] - the9DirComps[i] ) > 1e-7 )
684       isChanged = true;
685     _axisDirs[i] = the9DirComps[i];
686   }
687   if ( isChanged )
688     NotifySubMeshesHypothesisModification();
689 }
690
691 //=======================================================================
692 //function : GetGrid
693 //purpose  : Return coordinates of node positions along the three axes
694 //=======================================================================
695
696 void StdMeshers_CartesianParameters3D::GetGrid(std::vector<double>& coords, int axis) const
697   throw ( SALOME_Exception )
698 {
699   if ( IsGridBySpacing(axis) )
700     throw SALOME_Exception(LOCALIZED("The grid is defined by spacing and not by coordinates"));
701
702   coords = _coords[axis];
703 }
704
705 //=======================================================================
706 //function : GetSizeThreshold
707 //purpose  : Return size threshold
708 //=======================================================================
709
710 double StdMeshers_CartesianParameters3D::GetSizeThreshold() const
711 {
712   return _sizeThreshold;
713 }
714
715 //=======================================================================
716 //function : SetToAddEdges
717 //purpose  : Enables implementation of geometrical edges into the mesh. If this feature
718 //           is disabled, sharp edges of the shape are lost ("smoothed") in the mesh if
719 //           they don't coincide with the grid lines
720 //=======================================================================
721
722 void StdMeshers_CartesianParameters3D::SetToAddEdges(bool toAdd)
723 {
724   if ( _toAddEdges != toAdd )
725   {
726     _toAddEdges = toAdd;
727     NotifySubMeshesHypothesisModification();
728   }
729 }
730
731 //=======================================================================
732 //function : GetToAddEdges
733 //purpose  : Returns true if implementation of geometrical edges into the
734 //           mesh is enabled
735 //=======================================================================
736
737 bool StdMeshers_CartesianParameters3D::GetToAddEdges() const
738 {
739   return _toAddEdges;
740 }
741
742 //=======================================================================
743 //function : IsDefined
744 //purpose  : Return true if parameters are well defined
745 //=======================================================================
746
747 bool StdMeshers_CartesianParameters3D::IsDefined() const
748 {
749   for ( int i = 0; i < 3; ++i )
750     if (_coords[i].empty() && (_spaceFunctions[i].empty() || _internalPoints[i].empty()))
751       return false;
752
753   return ( _sizeThreshold > 1.0 );
754 }
755
756 //=======================================================================
757 //function : SaveTo
758 //purpose  : store my parameters into a stream
759 //=======================================================================
760
761 std::ostream & StdMeshers_CartesianParameters3D::SaveTo(std::ostream & save)
762 {
763   save << _sizeThreshold << " ";
764
765   for ( int i = 0; i < 3; ++i )
766   {
767     save << _coords[i].size() << " ";
768     for ( size_t j = 0; j < _coords[i].size(); ++j )
769       save << _coords[i][j] << " ";
770
771     save << _internalPoints[i].size() << " ";
772     for ( size_t j = 0; j < _internalPoints[i].size(); ++j )
773       save << _internalPoints[i][j] << " ";
774
775     save << _spaceFunctions[i].size() << " ";
776     for ( size_t j = 0; j < _spaceFunctions[i].size(); ++j )
777       save << _spaceFunctions[i][j] << " ";
778   }
779   save << _toAddEdges << " ";
780
781   save.setf( save.scientific );
782   save.precision( 12 );
783   for ( int i = 0; i < 9; ++i )
784     save << _axisDirs[i] << " ";
785
786   for ( int i = 0; i < 3; ++i )
787     save << _fixedPoint[i] << " ";
788
789   return save;
790 }
791
792 //=======================================================================
793 //function : LoadFrom
794 //purpose  : restore my parameters from a stream
795 //=======================================================================
796
797 std::istream & StdMeshers_CartesianParameters3D::LoadFrom(std::istream & load)
798 {
799   bool ok;
800
801   ok = static_cast<bool>( load >> _sizeThreshold );
802   for ( int ax = 0; ax < 3; ++ax )
803   {
804     if (ok)
805     {
806       size_t i = 0;
807       ok = static_cast<bool>(load >> i  );
808       if ( i > 0 && ok )
809       {
810         _coords[ax].resize( i );
811         for ( i = 0; i < _coords[ax].size() && ok; ++i )
812           ok = static_cast<bool>(load >> _coords[ax][i]  );
813       }
814     }
815     if (ok)
816     {
817       size_t i = 0;
818       ok = static_cast<bool>(load >> i  );
819       if ( i > 0 && ok )
820       {
821         _internalPoints[ax].resize( i );
822         for ( i = 0; i < _internalPoints[ax].size() && ok; ++i )
823           ok = static_cast<bool>(load >> _internalPoints[ax][i]  );
824       }
825     }
826     if (ok)
827     {
828       size_t i = 0;
829       ok = static_cast<bool>(load >> i  );
830       if ( i > 0 && ok )
831       {
832         _spaceFunctions[ax].resize( i );
833         for ( i = 0; i < _spaceFunctions[ax].size() && ok; ++i )
834           ok = static_cast<bool>(load >> _spaceFunctions[ax][i]  );
835       }
836     }
837   }
838
839   ok = static_cast<bool>( load >> _toAddEdges );
840
841   for ( int i = 0; i < 9 && ok; ++i )
842     ok = static_cast<bool>( load >> _axisDirs[i]);
843
844   for ( int i = 0; i < 3 && ok ; ++i )
845     ok = static_cast<bool>( load >> _fixedPoint[i]);
846
847   return load;
848 }
849
850 //=======================================================================
851 //function : SetParametersByMesh
852 //=======================================================================
853
854 bool StdMeshers_CartesianParameters3D::SetParametersByMesh(const SMESH_Mesh*   ,
855                                                            const TopoDS_Shape& )
856 {
857   return false;
858 }
859
860 //=======================================================================
861 //function : SetParametersByDefaults
862 //=======================================================================
863
864 bool StdMeshers_CartesianParameters3D::SetParametersByDefaults(const TDefaults&  dflts,
865                                                                const SMESH_Mesh* /*theMesh*/)
866 {
867   if ( dflts._elemLength > 1e-100 )
868   {
869     vector<string> spacing( 1, SMESH_Comment(dflts._elemLength));
870     vector<double> intPnts;
871     SetGridSpacing( spacing, intPnts, 0 );
872     SetGridSpacing( spacing, intPnts, 1 );
873     SetGridSpacing( spacing, intPnts, 2 );
874     return true;
875   }
876   return false;
877 }
878