X-Git-Url: http://git.salome-platform.org/gitweb/?p=modules%2Fsmesh.git;a=blobdiff_plain;f=src%2FStdMeshers%2FStdMeshers_CartesianParameters3D.cxx;h=b2e15df8c1a8ec1f5d5b9b782f800150c3047a97;hp=369d4e424151d5fc98bc5a89379911ca15926c25;hb=a0f09b9f1b8f5eac0e1c9277f76d65eb643cac94;hpb=f7aba4830d53719b963fdb7fccc98b760fdef2d1 diff --git a/src/StdMeshers/StdMeshers_CartesianParameters3D.cxx b/src/StdMeshers/StdMeshers_CartesianParameters3D.cxx index 369d4e424..b2e15df8c 100644 --- a/src/StdMeshers/StdMeshers_CartesianParameters3D.cxx +++ b/src/StdMeshers/StdMeshers_CartesianParameters3D.cxx @@ -1,4 +1,4 @@ -// Copyright (C) 2007-2013 CEA/DEN, EDF R&D, OPEN CASCADE +// Copyright (C) 2007-2014 CEA/DEN, EDF R&D, OPEN CASCADE // // Copyright (C) 2003-2007 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN, // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS @@ -6,7 +6,7 @@ // This library is free software; you can redistribute it and/or // modify it under the terms of the GNU Lesser General Public // License as published by the Free Software Foundation; either -// version 2.1 of the License. +// version 2.1 of the License, or (at your option) any later version. // // This library is distributed in the hope that it will be useful, // but WITHOUT ANY WARRANTY; without even the implied warranty of @@ -32,10 +32,25 @@ #include "utilities.h" +#include #include +#include +#include #include +#include +#include +#include #include +#include +#include +#include +#include +#include +#include +#include +#include +#include #include using namespace std; @@ -66,6 +81,11 @@ StdMeshers_CartesianParameters3D::StdMeshers_CartesianParameters3D(int h _axisDirs[6] = 0.; _axisDirs[7] = 0.; _axisDirs[8] = 1.; + + _fixedPoint[0] = 0.; + _fixedPoint[1] = 0.; + _fixedPoint[2] = 0.; + SetFixedPoint( _fixedPoint, /*toUnset=*/true ); } @@ -73,6 +93,22 @@ namespace { const char* axisName[3] = { "X", "Y", "Z" }; + typedef std::pair< double, std::pair< double, double > > TCooTriple; + +#define gpXYZ( cTriple ) gp_XYZ( (cTriple).first, (cTriple).second.first, (cTriple).second.second ) + + //================================================================================ + /*! + * \brief Compare two normals + */ + //================================================================================ + + bool sameDir( const TCooTriple& n1, const TCooTriple& n2 ) + { + gp_XYZ xyz1 = gpXYZ( n1 ), xyz2 = gpXYZ( n2 ); + return ( xyz1 - xyz2 ).Modulus() < 0.01; + } + //================================================================================ /*! * \brief Checks validity of an axis index, throws in case of invalidity @@ -177,6 +213,36 @@ void StdMeshers_CartesianParameters3D::SetGridSpacing(std::vector& xSpac NotifySubMeshesHypothesisModification(); } +//======================================================================= +//function : SetFixedPoint +//purpose : * Set/unset a fixed point, at which a node will be created provided that grid +// * is defined by spacing in all directions +//======================================================================= + +void StdMeshers_CartesianParameters3D::SetFixedPoint(const double p[3], bool toUnset) +{ + if ( toUnset != Precision::IsInfinite( _fixedPoint[0] )) + NotifySubMeshesHypothesisModification(); + + if ( toUnset ) + _fixedPoint[0] = Precision::Infinite(); + else + std::copy( &p[0], &p[0]+3, &_fixedPoint[0] ); +} + +//======================================================================= +//function : GetFixedPoint +//purpose : Returns either false or (true + point coordinates) +//======================================================================= + +bool StdMeshers_CartesianParameters3D::GetFixedPoint(double p[3]) const +{ + if ( Precision::IsInfinite( _fixedPoint[0] )) + return false; + std::copy( &_fixedPoint[0], &_fixedPoint[0]+3, &p[0] ); +} + + //======================================================================= //function : SetSizeThreshold //purpose : Set size threshold @@ -229,20 +295,41 @@ bool StdMeshers_CartesianParameters3D::IsGridBySpacing(const int axis) const //purpose : Computes node coordinates by spacing functions //======================================================================= -void StdMeshers_CartesianParameters3D::ComputeCoordinates(const double x0, - const double x1, - vector& spaceFuns, - vector& points, - vector& coords, - const std::string& axis ) +void StdMeshers_CartesianParameters3D::ComputeCoordinates(const double x0, + const double x1, + vector& theSpaceFuns, + vector& thePoints, + vector& coords, + const string& axis, + const double* xForced ) throw ( SALOME_Exception ) { - checkGridSpacing( spaceFuns, points, axis ); + checkGridSpacing( theSpaceFuns, thePoints, axis ); + + vector spaceFuns = theSpaceFuns; + vector points = thePoints; + + bool forced = false; + if (( forced = ( xForced && ( x0 < *xForced ) && ( *xForced < x1 )))) + { + // divide a range at xForced + + // find a range to insert xForced + double pos = ( *xForced - x0 ) / ( x1 - x0 ); + int iR = 1; + while ( pos > points[ iR ] ) ++iR; + + // insert xForced + vector::iterator pntIt = points.begin() + iR; + points.insert( pntIt, pos ); + vector::iterator funIt = spaceFuns.begin() + iR; + spaceFuns.insert( funIt, spaceFuns[ iR-1 ]); + } coords.clear(); for ( size_t i = 0; i < spaceFuns.size(); ++i ) { - FunctionExpr fun( spaceFuns[i].c_str(), /*convMode=*/-1 ); + StdMeshers::FunctionExpr fun( spaceFuns[i].c_str(), /*convMode=*/-1 ); const double p0 = x0 * ( 1. - points[i]) + x1 * points[i]; const double p1 = x0 * ( 1. - points[i+1]) + x1 * points[i+1]; @@ -280,6 +367,28 @@ void StdMeshers_CartesianParameters3D::ComputeCoordinates(const double x if ( fabs( coords.back() - p1 ) > 0.5 * lastCellLen ) coords.push_back ( p1 ); } + + // correct coords if a forced point is too close to a neighbor node + if ( forced ) + { + int iF = 0; + double minLen = ( x1 - x0 ); + for ( size_t i = 1; i < coords.size(); ++i ) + { + if ( !iF && Abs( coords[i] - *xForced ) < 1e-20 ) + iF = i++; // xForced found + else + minLen = Min( minLen, coords[i] - coords[i-1] ); + } + const double tol = minLen * 1e-3; + int iRem = -1; + if (( iF > 1 ) && ( coords[iF] - coords[iF-1] < tol )) + iRem = iF-1; + else if (( iF < coords.size()-2 ) && ( coords[iF+1] - coords[iF] < tol )) + iRem = iF+1; + if ( iRem > 0 ) + coords.erase( coords.begin() + iRem ); + } } //======================================================================= @@ -305,26 +414,242 @@ void StdMeshers_CartesianParameters3D::GetCoordinates(std::vector& xNode bndBox.Get(x0,y0,z0, x1,y1,z1); } + double fp[3], *pfp[3] = { NULL, NULL, NULL }; + if ( GetFixedPoint( fp )) + { + // convert fp into a basis defined by _axisDirs + gp_XYZ axis[3] = { gp_XYZ( _axisDirs[0], _axisDirs[1], _axisDirs[2] ), + gp_XYZ( _axisDirs[3], _axisDirs[4], _axisDirs[5] ), + gp_XYZ( _axisDirs[6], _axisDirs[7], _axisDirs[8] ) }; + axis[0].Normalize(); + axis[1].Normalize(); + axis[2].Normalize(); + + gp_Mat basis( axis[0], axis[1], axis[2] ); + gp_Mat bi = basis.Inverted(); + + gp_XYZ p( fp[0], fp[1], fp[2] ); + p *= bi; + p.Coord( fp[0], fp[1], fp[2] ); + + pfp[0] = & fp[0]; + pfp[1] = & fp[1]; + pfp[2] = & fp[2]; + } + StdMeshers_CartesianParameters3D* me = const_cast(this); if ( IsGridBySpacing(0) ) - ComputeCoordinates( x0, x1, me->_spaceFunctions[0], me->_internalPoints[0], xNodes, "X" ); + ComputeCoordinates + ( x0, x1, me->_spaceFunctions[0], me->_internalPoints[0], xNodes, "X", pfp[0] ); else xNodes = _coords[0]; if ( IsGridBySpacing(1) ) - ComputeCoordinates( y0, y1, me->_spaceFunctions[1], me->_internalPoints[1], yNodes, "Y" ); + ComputeCoordinates + ( y0, y1, me->_spaceFunctions[1], me->_internalPoints[1], yNodes, "Y", pfp[1] ); else yNodes = _coords[1]; if ( IsGridBySpacing(2) ) - ComputeCoordinates( z0, z1, me->_spaceFunctions[2], me->_internalPoints[2], zNodes, "Z" ); + ComputeCoordinates + ( z0, z1, me->_spaceFunctions[2], me->_internalPoints[2], zNodes, "Z", pfp[2] ); else zNodes = _coords[2]; } +//======================================================================= +//function : ComputeOptimalAxesDirs +//purpose : Returns axes at which number of hexahedra is maximal +//======================================================================= + +void StdMeshers_CartesianParameters3D:: +ComputeOptimalAxesDirs(const TopoDS_Shape& shape, + const bool isOrthogonal, + double dirCoords[9]) +{ + for ( int i = 0; i < 9; ++i ) dirCoords[i] = 0.; + dirCoords[0] = dirCoords[4] = dirCoords[8] = 1.; + + if ( shape.IsNull() ) return; + + TopLoc_Location loc; + TopExp_Explorer exp; + + // get external FACEs of the shape + TopTools_MapOfShape faceMap; + for ( exp.Init( shape, TopAbs_FACE ); exp.More(); exp.Next() ) + if ( !faceMap.Add( exp.Current() )) + faceMap.Remove( exp.Current() ); + + // sort areas of planar faces by normal direction + + std::multimap< TCooTriple, double > areasByNormal; + + TopTools_MapIteratorOfMapOfShape fIt ( faceMap ); + for ( ; fIt.More(); fIt.Next() ) + { + const TopoDS_Face& face = TopoDS::Face( fIt.Key() ); + Handle(Geom_Surface) surf = BRep_Tool::Surface( face, loc ); + if ( surf.IsNull() ) continue; + + GeomLib_IsPlanarSurface check( surf, 1e-5 ); + if ( !check.IsPlanar() ) continue; + + GProp_GProps SProps; + BRepGProp::SurfaceProperties( face, SProps ); + double area = SProps.Mass(); + + gp_Pln pln = check.Plan(); + gp_Dir norm = pln.Axis().Direction().Transformed( loc ); + if ( norm.X() < -1e-3 ) { // negative X + norm.Reverse(); + } else if ( norm.X() < 1e-3 ) { // zero X + if ( norm.Y() < -1e-3 ) { // negative Y + norm.Reverse(); + } else if ( norm.Y() < 1e-3 ) { // zero X && zero Y + if ( norm.Y() < -1e-3 ) // negative Z + norm.Reverse(); + } + } + TCooTriple coo3( norm.X(), make_pair( norm.Y(), norm.Z() )); + areasByNormal.insert( make_pair( coo3, area )); + } + + // group coplanar normals and sort groups by sum area + + std::multimap< double, vector< const TCooTriple* > > normsByArea; + std::multimap< TCooTriple, double >::iterator norm2a = areasByNormal.begin(); + const TCooTriple* norm1 = 0; + double sumArea = 0; + vector< const TCooTriple* > norms; + for ( int iF = 1; norm2a != areasByNormal.end(); ++norm2a, ++iF ) + { + + if ( !norm1 || !sameDir( *norm1, norm2a->first )) + { + if ( !norms.empty() ) + { + normsByArea.insert( make_pair( sumArea, norms )); + norms.clear(); + } + norm1 = & norm2a->first; + sumArea = norm2a->second; + norms.push_back( norm1 ); + } + else + { + sumArea += norm2a->second; + norms.push_back( & norm2a->first ); + } + if ( iF == areasByNormal.size() ) + normsByArea.insert( make_pair( sumArea, norms )); + } + + // try to set dirs by planar faces + + gp_XYZ normDirs[3]; // normals to largest planes + + if ( !normsByArea.empty() ) + { + norm1 = normsByArea.rbegin()->second[0]; + normDirs[0] = gpXYZ( *norm1 ); + + if ( normsByArea.size() == 1 ) + { + normDirs[1] = normDirs[0]; + if ( Abs( normDirs[0].Y() ) < 1e-100 && + Abs( normDirs[0].Z() ) < 1e-100 ) // normDirs[0] || OX + normDirs[1].SetY( normDirs[0].Y() + 1. ); + else + normDirs[1].SetX( normDirs[0].X() + 1. ); + } + else + { + // look for 2 other directions + gp_XYZ testDir = normDirs[0], minDir, maxDir; + for ( int is2nd = 0; is2nd < 2; ++is2nd ) + { + double maxMetric = 0, minMetric = 1e100; + std::multimap< double, vector< const TCooTriple* > >::iterator a2n; + for ( a2n = normsByArea.begin(); a2n != normsByArea.end(); ++a2n ) + { + gp_XYZ n = gpXYZ( *( a2n->second[0]) ); + double dot = Abs( n * testDir ); + double metric = ( 1. - dot ) * ( isOrthogonal ? 1 : a2n->first ); + if ( metric > maxMetric ) + { + maxDir = n; + maxMetric = metric; + } + if ( metric < minMetric ) + { + minDir = n; + minMetric = metric; + } + } + if ( is2nd ) + { + normDirs[2] = minDir; + } + else + { + normDirs[1] = maxDir; + normDirs[2] = normDirs[0] ^ normDirs[1]; + if ( isOrthogonal || normsByArea.size() < 3 ) + break; + testDir = normDirs[2]; + } + } + } + if ( isOrthogonal || normsByArea.size() == 1 ) + { + normDirs[2] = normDirs[0] ^ normDirs[1]; + normDirs[1] = normDirs[2] ^ normDirs[0]; + } + } + else + { + return; + } + + gp_XYZ dirs[3]; + dirs[0] = normDirs[0] ^ normDirs[1]; + dirs[1] = normDirs[1] ^ normDirs[2]; + dirs[2] = normDirs[2] ^ normDirs[0]; + + dirs[0].Normalize(); + dirs[1].Normalize(); + dirs[2].Normalize(); + + // Select dirs for X, Y and Z axes + int iX = ( Abs( dirs[0].X() ) > Abs( dirs[1].X() )) ? 0 : 1; + if ( Abs( dirs[iX].X() ) < Abs( dirs[2].X() )) + iX = 2; + int iY = ( iX == 0 ) ? 1 : (( Abs( dirs[0].Y() ) > Abs( dirs[1].Y() )) ? 0 : 1 ); + if ( Abs( dirs[iY].Y() ) < Abs( dirs[2].Y() ) && iX != 2 ) + iY = 2; + int iZ = 3 - iX - iY; + + if ( dirs[iX].X() < 0 ) dirs[iX].Reverse(); + if ( dirs[iY].Y() < 0 ) dirs[iY].Reverse(); + gp_XYZ zDir = dirs[iX] ^ dirs[iY]; + if ( dirs[iZ] * zDir < 0 ) + dirs[iZ].Reverse(); + + dirCoords[0] = dirs[iX].X(); + dirCoords[1] = dirs[iX].Y(); + dirCoords[2] = dirs[iX].Z(); + dirCoords[3] = dirs[iY].X(); + dirCoords[4] = dirs[iY].Y(); + dirCoords[5] = dirs[iY].Z(); + dirCoords[6] = dirs[iZ].X(); + dirCoords[7] = dirs[iZ].Y(); + dirCoords[8] = dirs[iZ].Z(); +} + //======================================================================= //function : SetAxisDirs -//purpose : Sets directions of axes +//purpose : Sets custom direction of axes //======================================================================= void StdMeshers_CartesianParameters3D::SetAxisDirs(const double* the9DirComps) @@ -349,6 +674,10 @@ void StdMeshers_CartesianParameters3D::SetAxisDirs(const double* the9DirComps) y.IsParallel( z, M_PI / 180. )) throw SALOME_Exception("Parallel axis directions"); + gp_Vec normXY = x ^ y, normYZ = y ^ z; + if ( normXY.IsParallel( normYZ, M_PI / 180. )) + throw SALOME_Exception("Axes lie in one plane"); + bool isChanged = false; for ( int i = 0; i < 9; ++i ) { @@ -450,6 +779,14 @@ std::ostream & StdMeshers_CartesianParameters3D::SaveTo(std::ostream & save) } save << _toAddEdges << " "; + save.setf( save.scientific ); + save.precision( 12 ); + for ( int i = 0; i < 9; ++i ) + save << _axisDirs[i] << " "; + + for ( int i = 0; i < 3; ++i ) + save << _fixedPoint[i] << " "; + return save; } @@ -500,7 +837,13 @@ std::istream & StdMeshers_CartesianParameters3D::LoadFrom(std::istream & load) } } - load >> _toAddEdges; + ok = ( load >> _toAddEdges ); + + for ( int i = 0; i < 9 && ok; ++i ) + ok = ( load >> _axisDirs[i]); + + for ( int i = 0; i < 3 && ok ; ++i ) + ok = ( load >> _fixedPoint[i]); return load; }