1 // Copyright (C) 2007-2019 CEA/DEN, EDF R&D, OPEN CASCADE
3 // Copyright (C) 2003-2007 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
4 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
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.
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.
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
20 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
22 // File : SMESH_Delaunay.cxx
23 // Created : Wed Apr 19 15:41:15 2017
24 // Author : Edward AGAPOV (eap)
26 #include "SMESH_Delaunay.hxx"
28 #include "SMESH_Comment.hxx"
29 #include "SMESH_File.hxx"
30 #include "SMESH_MeshAlgos.hxx"
32 #include <BRepAdaptor_Surface.hxx>
33 #include <BRepMesh_Delaun.hxx>
35 #include <Basics_OCCTVersion.hxx>
37 //================================================================================
39 * \brief Construct a Delaunay triangulation of given boundary nodes
40 * \param [in] boundaryNodes - vector of nodes of a wire
41 * \param [in] face - the face
42 * \param [in] faceID - the face ID
44 //================================================================================
46 SMESH_Delaunay::SMESH_Delaunay(const std::vector< const UVPtStructVec* > & boundaryNodes,
47 const TopoDS_Face& face,
49 : _face( face ), _faceID( faceID ), _scale( 1., 1. )
52 BRepAdaptor_Surface surf( face );
53 if ( surf.GetType() != GeomAbs_Plane )
55 const int nbDiv = 100;
56 const double uRange = surf.LastUParameter() - surf.FirstUParameter();
57 const double vRange = surf.LastVParameter() - surf.FirstVParameter();
58 const double uFixed = surf.FirstUParameter() + 0.5 * uRange;
59 const double vFixed = surf.FirstVParameter() + 0.5 * vRange;
60 const double dU = uRange / nbDiv;
61 const double dV = vRange / nbDiv;
62 double u = surf.FirstUParameter(), v = surf.FirstVParameter();
63 gp_Pnt p0U = surf.Value( u, v ), p0V = p0U;
64 double lenU = 0, lenV = 0;
65 for ( ; u < surf.LastUParameter(); u += dU, v += dV )
67 gp_Pnt p1U = surf.Value( u, vFixed );
68 lenU += p1U.Distance( p0U );
70 gp_Pnt p1V = surf.Value( uFixed, v );
71 lenV += p1V.Distance( p0V );
74 _scale.SetCoord( lenU / uRange, lenV / vRange );
77 // count boundary points
79 for ( size_t iW = 0; iW < boundaryNodes.size(); ++iW ) // loop on wires
81 nbP += boundaryNodes[iW]->size();
82 if ( boundaryNodes[iW]->front().node == boundaryNodes[iW]->back().node )
83 --nbP; // 1st and last points coincide
85 _bndNodes.resize( nbP );
87 // fill boundary points
88 #if OCC_VERSION_LARGE <= 0x07030000
89 BRepMesh::Array1OfVertexOfDelaun bndVert( 1, 1 + nbP );
91 IMeshData::Array1OfVertexOfDelaun bndVert( 1, 1 + nbP );
93 BRepMesh_Vertex v( 0, 0, BRepMesh_Frontier );
94 for ( size_t iW = 0; iW < boundaryNodes.size(); ++iW )
96 const UVPtStructVec& bndPnt = *boundaryNodes[iW];
97 int i = 0, nb = bndPnt.size();
98 if ( bndPnt[0].node == bndPnt.back().node )
100 for ( ; i < nb; ++i, ++iP )
102 _bndNodes[ iP-1 ] = bndPnt[i].node;
103 bndPnt[i].node->setIsMarked( true );
105 v.ChangeCoord() = bndPnt[i].UV().Multiplied( _scale );
110 // triangulate the srcFace in 2D
111 BRepMesh_Delaun Delaunay( bndVert );
112 _triaDS = Delaunay.Result();
115 //================================================================================
117 * \brief Prepare to the exploration of nodes
119 //================================================================================
121 void SMESH_Delaunay::InitTraversal(const int nbNodesToVisit)
123 _nbNodesToVisit = (size_t) nbNodesToVisit;
124 _nbVisitedNodes = _iBndNode = 0;
128 //================================================================================
130 * \brief Return a node with its Barycentric Coordinates within the triangle
131 * defined by its node indices (zero based)
132 * \param [out] bc - Barycentric Coordinates of the returned node
133 * \param [out] triaNodes - indices of triangle nodes
134 * \return const SMDS_MeshNode* - the next node or NULL
136 //================================================================================
138 const SMDS_MeshNode* SMESH_Delaunay::NextNode( double bc[3], int triaNodes[3] )
140 while ( _nbVisitedNodes < _nbNodesToVisit )
142 while ( !_noTriQueue.empty() )
144 const SMDS_MeshNode* node = _noTriQueue.front().first;
145 const BRepMesh_Triangle* tria = _noTriQueue.front().second;
146 _noTriQueue.pop_front();
147 if ( node->isMarked() )
150 node->setIsMarked( true );
152 // find a Delaunay triangle containing the src node
153 gp_XY uv = getNodeUV( _face, node );
154 tria = FindTriangle( uv, tria, bc, triaNodes );
157 addCloseNodes( node, tria, _faceID, _noTriQueue );
161 for ( ; _iBndNode < _bndNodes.size() && _noTriQueue.empty(); ++_iBndNode )
163 if ( const BRepMesh_Triangle* tria = GetTriangleNear( _iBndNode ))
164 addCloseNodes( _bndNodes[ _iBndNode ], tria, _faceID, _noTriQueue );
166 if ( _noTriQueue.empty() )
170 // if ( _nbVisitedNodes < _nbNodesToVisit )
171 // _nbVisitedNodes = std::numeric_limits<int>::max();
175 //================================================================================
177 * \brief Find a Delaunay triangle containing a given 2D point and return
178 * barycentric coordinates within the found triangle
180 //================================================================================
182 const BRepMesh_Triangle* SMESH_Delaunay::FindTriangle( const gp_XY& UV,
183 const BRepMesh_Triangle* tria,
190 Standard_Boolean ori[3];
192 gp_XY uv = UV.Multiplied( _scale );
196 // check if the uv is in tria
198 _triaDS->ElementNodes( *tria, nodeIDs );
199 nodeUVs[0] = _triaDS->GetNode( nodeIDs[0] ).Coord();
200 nodeUVs[1] = _triaDS->GetNode( nodeIDs[1] ).Coord();
201 nodeUVs[2] = _triaDS->GetNode( nodeIDs[2] ).Coord();
203 SMESH_MeshAlgos::GetBarycentricCoords( uv,
204 nodeUVs[0], nodeUVs[1], nodeUVs[2],
206 if ( bc[0] >= 0 && bc[1] >= 0 && bc[0] + bc[1] <= 1 )
208 if ( _triaDS->GetNode( nodeIDs[0] ).Movability() != BRepMesh_Frontier ||
209 _triaDS->GetNode( nodeIDs[1] ).Movability() != BRepMesh_Frontier ||
210 _triaDS->GetNode( nodeIDs[2] ).Movability() != BRepMesh_Frontier )
214 bc[2] = 1 - bc[0] - bc[1];
215 triaNodes[0] = nodeIDs[0] - 1;
216 triaNodes[1] = nodeIDs[1] - 1;
217 triaNodes[2] = nodeIDs[2] - 1;
221 // look for a neighbor triangle, which is adjacent to a link intersected
222 // by a segment( triangle center -> uv )
224 gp_XY gc = ( nodeUVs[0] + nodeUVs[1] + nodeUVs[2] ) / 3.;
227 tria->Edges( linkIDs, ori );
228 #if OCC_VERSION_LARGE <= 0x07030000
229 int triaID = _triaDS->IndexOf( *tria );
231 int triaID = tria - & ( _triaDS->GetElement( 0 ));
235 for ( int i = 0; i < 3; ++i )
237 const BRepMesh_PairOfIndex & triIDs = _triaDS->ElementsConnectedTo( linkIDs[i] );
238 if ( triIDs.Extent() < 2 )
239 continue; // no neighbor triangle
241 // check if a link intersects gc2uv
242 const BRepMesh_Edge & link = _triaDS->GetLink( linkIDs[i] );
243 const BRepMesh_Vertex & n1 = _triaDS->GetNode( link.FirstNode() );
244 const BRepMesh_Vertex & n2 = _triaDS->GetNode( link.LastNode() );
245 gp_XY uv1 = n1.Coord();
246 gp_XY lin = n2.Coord() - uv1; // link direction
248 double crossSegLin = seg ^ lin;
249 if ( Abs( crossSegLin ) < std::numeric_limits<double>::min() )
250 continue; // parallel
252 double uSeg = ( uv1 - gc ) ^ lin / crossSegLin;
253 if ( 0. <= uSeg && uSeg <= 1. )
255 tria = & _triaDS->GetElement( triIDs.Index( 1 + ( triIDs.Index(1) == triaID )));
256 if ( tria->Movability() != BRepMesh_Deleted )
264 //================================================================================
266 * \brief Return a triangle sharing a given boundary node
267 * \param [in] iBndNode - index of the boundary node
268 * \return const BRepMesh_Triangle* - a found triangle
270 //================================================================================
272 const BRepMesh_Triangle* SMESH_Delaunay::GetTriangleNear( int iBndNode )
274 if ( iBndNode >= _triaDS->NbNodes() )
277 int nbNbNodes = _bndNodes.size();
278 #if OCC_VERSION_LARGE <= 0x07030000
279 typedef BRepMesh::ListOfInteger TLinkList;
281 typedef IMeshData::ListOfInteger TLinkList;
283 const TLinkList & linkIds = _triaDS->LinksConnectedTo( iBndNode + 1 );
284 TLinkList::const_iterator iLink = linkIds.cbegin();
285 for ( ; iLink != linkIds.cend(); ++iLink )
287 const BRepMesh_PairOfIndex & triaIds = _triaDS->ElementsConnectedTo( *iLink );
289 const BRepMesh_Triangle& tria = _triaDS->GetElement( triaIds.Index(1) );
290 if ( tria.Movability() != BRepMesh_Deleted )
292 _triaDS->ElementNodes( tria, nodeIDs );
293 if ( nodeIDs[0]-1 < nbNbNodes &&
294 nodeIDs[1]-1 < nbNbNodes &&
295 nodeIDs[2]-1 < nbNbNodes )
299 if ( triaIds.Extent() > 1 )
301 const BRepMesh_Triangle& tria = _triaDS->GetElement( triaIds.Index(2) );
302 if ( tria.Movability() != BRepMesh_Deleted )
304 _triaDS->ElementNodes( tria, nodeIDs );
305 if ( nodeIDs[0]-1 < nbNbNodes &&
306 nodeIDs[1]-1 < nbNbNodes &&
307 nodeIDs[2]-1 < nbNbNodes )
315 //================================================================================
317 * \brief Return UV of the i-th source boundary node (zero based)
319 //================================================================================
321 gp_XY SMESH_Delaunay::GetBndUV(const int iNode) const
323 return _triaDS->GetNode( iNode+1 ).Coord();
326 //================================================================================
328 * \brief Add non-marked nodes surrounding a given one to a queue
330 //================================================================================
332 void SMESH_Delaunay::addCloseNodes( const SMDS_MeshNode* node,
333 const BRepMesh_Triangle* tria,
335 TNodeTriaList & _noTriQueue )
337 // find in-FACE nodes
338 SMDS_ElemIteratorPtr elems = node->GetInverseElementIterator(SMDSAbs_Face);
339 while ( elems->more() )
341 const SMDS_MeshElement* elem = elems->next();
342 if ( elem->getshapeId() == faceID )
344 for ( int i = 0, nb = elem->NbNodes(); i < nb; ++i )
346 const SMDS_MeshNode* n = elem->GetNode( i );
347 if ( !n->isMarked() /*&& n->getshapeId() == faceID*/ )
348 _noTriQueue.push_back( std::make_pair( n, tria ));
354 //================================================================================
356 * \brief Write a python script that creates an equal mesh in Mesh module
358 //================================================================================
360 void SMESH_Delaunay::ToPython() const
363 text << "import salome, SMESH\n";
364 text << "salome.salome_init()\n";
365 text << "from salome.smesh import smeshBuilder\n";
366 text << "smesh = smeshBuilder.New()\n";
367 text << "mesh=smesh.Mesh()\n";
368 const char* endl = "\n";
370 for ( int i = 0; i < _triaDS->NbNodes(); ++i )
372 const BRepMesh_Vertex& v = _triaDS->GetNode( i+1 );
373 text << "mesh.AddNode( " << v.Coord().X() << ", " << v.Coord().Y() << ", 0 )" << endl;
377 for ( int i = 0; i < _triaDS->NbElements(); ++i )
379 const BRepMesh_Triangle& t = _triaDS->GetElement( i+1 );
380 if ( t.Movability() == BRepMesh_Deleted )
382 _triaDS->ElementNodes( t, nodeIDs );
383 text << "mesh.AddFace([ " << nodeIDs[0] << ", " << nodeIDs[1] << ", " << nodeIDs[2] << " ])" << endl;
386 const char* fileName = "/tmp/Delaunay.py";
387 SMESH_File file( fileName, false );
389 file.openForWriting();
390 file.write( text.c_str(), text.size() );
391 std::cout << "exec(open('" << fileName << "', 'rb').read())";