#include <limits>
-#undef WITH_TBB
+// #undef WITH_TBB
#ifdef WITH_TBB
#include <tbb/parallel_for.h>
//#include <tbb/enumerable_thread_specific.h>
*/
struct GridPlanes
{
- gp_XYZ _uNorm, _vNorm, _zNorm;
+ gp_XYZ _zNorm;
vector< gp_XYZ > _origins; // origin points of all planes in one direction
vector< double > _zProjs; // projections of origins to _zNorm
};
double _tol, _minCellSize;
gp_XYZ _origin;
gp_Mat _invB; // inverted basis of _axes
- //bool _isOrthogonalAxes;
vector< const SMDS_MeshNode* > _nodes; // mesh nodes at grid nodes
vector< const F_IntersectPoint* > _gridIntP; // grid node intersection with geometry
FaceGridIntersector(): _grid(0), _surfaceInt(0) {}
void Intersect();
- bool IsInGrid(const Bnd_Box& gridBox);
void StoreIntersections()
{
}
else if ( _link->_faces[0] == f )
{
- _link->_faces[0];
+ _link->_faces[0] = 0;
if ( _link->_faces[1] )
{
_link->_faces[0] = _link->_faces[1];
_invB.SetCols( _axes[0], _axes[1], _axes[2] );
_invB.Invert();
- // _isOrthogonalAxes = ( Abs( _axes[0] * _axes[1] ) < 1e-20 &&
- // Abs( _axes[1] * _axes[2] ) < 1e-20 &&
- // Abs( _axes[2] * _axes[0] ) < 1e-20 );
-
// compute tolerance
_minCellSize = Precision::Infinite();
for ( int iDir = 0; iDir < 3; ++iDir ) // loop on 3 line directions
*/
void Grid::ComputeUVW(const gp_XYZ& P, double UVW[3])
{
- // gp_XYZ p = P - _origin;
- // UVW[ 0 ] = p.X() * _invB( 1, 1 ) + p.Y() * _invB( 1, 2 ) + p.Z() * _invB( 1, 3 );
- // UVW[ 1 ] = p.X() * _invB( 2, 1 ) + p.Y() * _invB( 2, 2 ) + p.Z() * _invB( 2, 3 );
- // UVW[ 2 ] = p.X() * _invB( 3, 1 ) + p.Y() * _invB( 3, 2 ) + p.Z() * _invB( 3, 3 );
- // UVW[ 0 ] += _coords[0][0];
- // UVW[ 1 ] += _coords[1][0];
- // UVW[ 2 ] += _coords[2][0];
gp_XYZ p = P * _invB;
p.Coord( UVW[0], UVW[1], UVW[2] );
}
const gp_XYZ lineLoc = line._line.Location().XYZ();
const gp_XYZ lineDir = line._line.Direction().XYZ();
line.RemoveExcessIntPoints( _tol );
- multiset< F_IntersectPoint >& intPnts = _lines[ iDir ][ li.LineIndex() ]._intPoints;
+ multiset< F_IntersectPoint >& intPnts = line._intPoints;
multiset< F_IntersectPoint >::iterator ip = intPnts.begin();
bool isOut = true;
#endif
}
- //=============================================================================
- /*
- * Checks if the face is encosed by the grid
- */
- bool FaceGridIntersector::IsInGrid(const Bnd_Box& gridBox)
- {
- // double x0,y0,z0, x1,y1,z1;
- // const Bnd_Box& faceBox = GetFaceBndBox();
- // faceBox.Get(x0,y0,z0, x1,y1,z1);
-
- // if ( !gridBox.IsOut( gp_Pnt( x0,y0,z0 )) &&
- // !gridBox.IsOut( gp_Pnt( x1,y1,z1 )))
- // return true;
-
- // double X0,Y0,Z0, X1,Y1,Z1;
- // gridBox.Get(X0,Y0,Z0, X1,Y1,Z1);
- // double faceP[6] = { x0,y0,z0, x1,y1,z1 };
- // double gridP[6] = { X0,Y0,Z0, X1,Y1,Z1 };
- // gp_Dir axes[3] = { gp::DX(), gp::DY(), gp::DZ() };
- // for ( int iDir = 0; iDir < 6; ++iDir )
- // {
- // if ( iDir < 3 && gridP[ iDir ] <= faceP[ iDir ] ) continue;
- // if ( iDir >= 3 && gridP[ iDir ] >= faceP[ iDir ] ) continue;
-
- // // check if the face intersects a side of a gridBox
-
- // gp_Pnt p = iDir < 3 ? gp_Pnt( X0,Y0,Z0 ) : gp_Pnt( X1,Y1,Z1 );
- // gp_Ax1 norm( p, axes[ iDir % 3 ] );
- // if ( iDir < 3 ) norm.Reverse();
-
- // gp_XYZ O = norm.Location().XYZ(), N = norm.Direction().XYZ();
-
- // TopLoc_Location loc = _face.Location();
- // Handle(Poly_Triangulation) aPoly = BRep_Tool::Triangulation(_face,loc);
- // if ( !aPoly.IsNull() )
- // {
- // if ( !loc.IsIdentity() )
- // {
- // norm.Transform( loc.Transformation().Inverted() );
- // O = norm.Location().XYZ(), N = norm.Direction().XYZ();
- // }
- // const double deflection = aPoly->Deflection();
-
- // const TColgp_Array1OfPnt& nodes = aPoly->Nodes();
- // for ( int i = nodes.Lower(); i <= nodes.Upper(); ++i )
- // if (( nodes( i ).XYZ() - O ) * N > _grid->_tol + deflection )
- // return false;
- // }
- // else
- // {
- // BRepAdaptor_Surface surf( _face );
- // double u0, u1, v0, v1, du, dv, u, v;
- // BRepTools::UVBounds( _face, u0, u1, v0, v1);
- // if ( surf.GetType() == GeomAbs_Plane ) {
- // du = u1 - u0, dv = v1 - v0;
- // }
- // else {
- // du = surf.UResolution( _grid->_minCellSize / 10. );
- // dv = surf.VResolution( _grid->_minCellSize / 10. );
- // }
- // for ( u = u0, v = v0; u <= u1 && v <= v1; u += du, v += dv )
- // {
- // gp_Pnt p = surf.Value( u, v );
- // if (( p.XYZ() - O ) * N > _grid->_tol )
- // {
- // TopAbs_State state = GetCurveFaceIntersector()->ClassifyUVPoint(gp_Pnt2d( u, v ));
- // if ( state == TopAbs_IN || state == TopAbs_ON )
- // return false;
- // }
- // }
- // }
- // }
- return true;
- }
//=============================================================================
/*
* Intersects TopoDS_Face with all GridLine's
typedef void (FaceLineIntersector::* PIntFun )(const GridLine& gridLine);
PIntFun interFunction;
+ bool isDirect = true;
BRepAdaptor_Surface surf( _face );
switch ( surf.GetType() ) {
case GeomAbs_Plane:
intersector._plane = surf.Plane();
interFunction = &FaceLineIntersector::IntersectWithPlane;
+ isDirect = intersector._plane.Direct();
break;
case GeomAbs_Cylinder:
intersector._cylinder = surf.Cylinder();
interFunction = &FaceLineIntersector::IntersectWithCylinder;
+ isDirect = intersector._cylinder.Direct();
break;
case GeomAbs_Cone:
intersector._cone = surf.Cone();
interFunction = &FaceLineIntersector::IntersectWithCone;
+ //isDirect = intersector._cone.Direct();
break;
case GeomAbs_Sphere:
intersector._sphere = surf.Sphere();
interFunction = &FaceLineIntersector::IntersectWithSphere;
+ isDirect = intersector._sphere.Direct();
break;
case GeomAbs_Torus:
intersector._torus = surf.Torus();
interFunction = &FaceLineIntersector::IntersectWithTorus;
+ //isDirect = intersector._torus.Direct();
break;
default:
interFunction = &FaceLineIntersector::IntersectWithSurface;
}
+ if ( !isDirect )
+ std::swap( intersector._transOut, intersector._transIn );
_intersections.clear();
for ( int iDir = 0; iDir < 3; ++iDir ) // loop on 3 line directions
}
else // there are intersections with EDGEs
{
- // get a remaining link to start from, one lying on minimal
- // nb of FACEs
+ // get a remaining link to start from, one lying on minimal nb of FACEs
{
vector< pair< TGeomID, int > > facesOfLink[3];
pair< TGeomID, int > faceOfLink( -1, -1 );
polygon._links.push_back( *curLink );
--nbFreeLinks;
- // find all links bounding a FACE of curLink
+ // find all links lying on a curFace
do
{
// go forward from curLink
if ( polygon._links.back().NbFaces() > 0 )
polygon._links.front().AddFace( polygon._links.back()._link->_faces[0] );
- _polygons.pop_back();
+ if ( iPolygon == _polygons.size()-1 )
+ _polygons.pop_back();
}
else // polygon._links.size() >= 2
{
if ( iL2 == polygon._links.size() )
coplanarPolyg = 0;
}
- if ( 0 /*coplanarPolyg*/ ) // coplanar polygon found
+ if ( coplanarPolyg ) // coplanar polygon found
{
freeLinks.resize( freeLinks.size() - polygon._polyLinks.size() );
nbFreeLinks -= polygon._polyLinks.size();
+ // an artificial E_IntersectPoint used to mark nodes of coplanarPolyg
+ // as lying on curFace while they are not at intersection with geometry
+ E_IntersectPoint* ip = & _grid->_edgeIntP.back();
+ ip->_faceIDs.resize(1);
+ ip->_faceIDs[0] = curFace;
+
// fill freeLinks with links not shared by coplanarPolyg and polygon
for ( size_t iL = 0; iL < polygon._links.size(); ++iL )
if ( polygon._links[ iL ]._link->_faces[1] &&
for ( size_t iL2 = 0; iL2 < p->_links.size(); ++iL2 )
if ( p->_links[ iL2 ]._link == coplanarPolyg->_links[ iL ]._link )
{
- freeLinks.push_back( & p->_links[ iL2 ] );
+ // set links of coplanarPolyg in place of used freeLinks
+ // to re-create coplanarPolyg next
+ size_t iL3 = 0;
+ for ( ; iL3 < freeLinks.size() && freeLinks[ iL3 ]; ++iL3 );
+ if ( iL3 < freeLinks.size() )
+ freeLinks[ iL3 ] = ( & p->_links[ iL2 ] );
+ else
+ freeLinks.push_back( & p->_links[ iL2 ] );
++nbFreeLinks;
- freeLinks.back()->RemoveFace( coplanarPolyg );
+ freeLinks[ iL3 ]->RemoveFace( coplanarPolyg );
+ // mark nodes of coplanarPolyg as lying on curFace
+ for ( int iN = 0; iN < 2; ++iN )
+ {
+ _Node* n = freeLinks[ iL3 ]->_link->_nodes[ iN ];
+ if ( n->_intPoint ) n->_intPoint->Add( ip->_faceIDs );
+ else n->_intPoint = ip;
+ }
break;
}
}
usedFaceIDs.erase( curFace );
continue;
} // if ( coplanarPolyg )
- } // if ( hasEdgeIntersections )
+ } // if ( hasEdgeIntersections ) - search for coplanarPolyg
iPolygon = _polygons.size();
}
const int iLink = iL + iDir * 4;
hex->_hexLinks[iLink]._fIntPoints.push_back( &(*ip) );
- //hex->_hexLinks[iLink]._fIntNodes.push_back( _Node( 0, &(*ip) ));
hex->_nbFaceIntNodes += bool( ip->_node );
}
}
{
// all intersection of hex with geometry are at grid nodes
hex = new Hexahedron( *this );
- //hex->init( i );
hex->_i = _i;
hex->_j = _j;
hex->_k = _k;
GridPlanes& planes = pln[ iDirZ ];
int iDirX = ( iDirZ + 1 ) % 3;
int iDirY = ( iDirZ + 2 ) % 3;
- // planes._uNorm = ( _grid->_axes[ iDirY ] ^ _grid->_axes[ iDirZ ] ).Normalized();
- // planes._vNorm = ( _grid->_axes[ iDirZ ] ^ _grid->_axes[ iDirX ] ).Normalized();
planes._zNorm = ( _grid->_axes[ iDirX ] ^ _grid->_axes[ iDirY ] ).Normalized();
planes._zProjs.resize ( _grid->_coords[ iDirZ ].size() );
planes._zProjs [0] = 0;
double xLen = _grid->_coords[ iDirX ].back() - _grid->_coords[ iDirX ][0];
double yLen = _grid->_coords[ iDirY ].back() - _grid->_coords[ iDirY ][0];
double zLen = _grid->_coords[ iDirZ ].back() - _grid->_coords[ iDirZ ][0];
- //double zFactor = _grid->_axes[ iDirZ ] * planes._zNorm;
int dIJK[3], d000[3] = { 0,0,0 };
double o[3] = { _grid->_coords[0][0],
_grid->_coords[1][0],
} // loop on 3 grid directions
} // loop on EDGEs
- // Create nodes at found intersections
- // const E_IntersectPoint* eip;
- // for ( size_t i = 0; i < hexes.size(); ++i )
- // {
- // Hexahedron* h = hexes[i];
- // if ( !h ) continue;
- // for ( int iF = 0; iF < 6; ++iF )
- // {
- // _Face& quad = h->_hexQuads[ iF ];
- // for ( size_t iP = 0; iP < quad._eIntNodes.size(); ++iP )
- // if ( !quad._eIntNodes[ iP ]._node )
- // if (( eip = quad._eIntNodes[ iP ].EdgeIntPnt() ))
- // quad._eIntNodes[ iP ]._intPoint->_node = helper.AddNode( eip->_point.X(),
- // eip->_point.Y(),
- // eip->_point.Z() );
- // }
- // for ( size_t iP = 0; iP < hexes[i]->_vIntNodes.size(); ++iP )
- // if (( eip = h->_vIntNodes[ iP ].EdgeIntPnt() ))
- // h->_vIntNodes[ iP ]._intPoint->_node = helper.AddNode( eip->_point.X(),
- // eip->_point.Y(),
- // eip->_point.Z() );
- // }
+ // add an artificial E_IntersectPoint used in Hexahedron::ComputeElements()
+ ip._shapeID = -1;
+ ip._faceIDs.clear();
+ ip._node = 0;
+ _grid->_edgeIntP.push_back( ip );
}
//================================================================================
default: ; // tangent transition
}
}
+ if ( !link._nodes[1]->Node() )
+ return true;
+
gp_Pnt p1 = link._nodes[0]->Point();
gp_Pnt p2 = link._nodes[1]->Point();
gp_Pnt testPnt = 0.8 * p1.XYZ() + 0.2 * p2.XYZ();
return v * normal > 0;
}
}
- // if ( !_hexLinks[ iLink ]._nodes[0]->Node() ) // no node
- // return true;
- // if ( !_hexLinks[ iLink ]._nodes[0]->_intPoint ) // no intersection with geometry
- // return false;
- // switch ( _hexLinks[ iLink ]._nodes[0]->FaceIntPnt()->_transition ) {
- // case Trans_OUT: return true;
- // case Trans_IN : return false;
- // default: ; // tangent transition
- // }
-
-// // ijk of a GridLine corresponding to the link
-// int iDir = iLink / 4;
-// int indSub = iLink % 4;
-// LineIndexer li = _grid->GetLineIndexer( iDir );
-// li.SetIJK( _i,_j,_k );
-// size_t lineIndex[4] = { li.LineIndex (),
-// li.LineIndex10(),
-// li.LineIndex01(),
-// li.LineIndex11() };
-// GridLine& line = _grid->_lines[ iDir ][ lineIndex[ indSub ]];
-
-// // analyze transition of previous ip
-// bool isOut = true;
-// multiset< F_IntersectPoint >::const_iterator ip = line._intPoints.begin();
-// for ( ; ip != line._intPoints.end(); ++ip )
-// {
-// if ( &(*ip) == _hexLinks[ iLink ]._nodes[0]->_intPoint )
-// break;
-// switch ( ip->_transition ) {
-// case Trans_OUT: isOut = true;
-// case Trans_IN : isOut = false;
-// default:;
-// }
-// }
-// #ifdef _DEBUG_
-// if ( ip == line._intPoints.end() )
-// cout << "BUG: Wrong GridLine. IKJ = ( "<< _i << " " << _j << " " << _k << " )" << endl;
-// #endif
return isOut;
}
//================================================================================
// find a top node above the base node
_Link* link = _polygons[0]._links[iL]._link;
- //ASSERT( link->_faces.size() > 1 );
if ( !link->_faces[0] || !link->_faces[1] )
return debugDumpLink( link );
// a quadrangle sharing <link> with _polygons[0]
nodes[2] = _polygons[0]._links[2].FirstNode();
_Link* link = _polygons[0]._links[0]._link;
- //ASSERT( link->_faces.size() > 1 );
- if ( !link->_faces[0] || !link->_faces[1] )
+ if ( !link->_faces[0] || !link->_faces[1] )
return debugDumpLink( link );
// a triangle sharing <link> with _polygons[0]
// find a top node above the base node
_Link* link = _polygons[ iTri ]._links[iL]._link;
- //ASSERT( link->_faces.size() > 1 );
if ( !link->_faces[0] || !link->_faces[1] )
return debugDumpLink( link );
// a quadrangle sharing <link> with a base triangle
nodes[3] = _polygons[iQuad]._links[3].FirstNode();
_Link* link = _polygons[iQuad]._links[0]._link;
- //ASSERT( link->_faces.size() > 1 );
if ( !link->_faces[0] || !link->_faces[1] )
return debugDumpLink( link );
