]> SALOME platform Git repositories - modules/smesh.git/blob - src/SMESHUtils/SMESH_Delaunay.cxx
Salome HOME
#19232 [CEA] overal Mesh quality crash on imported MED
[modules/smesh.git] / src / SMESHUtils / SMESH_Delaunay.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 // File      : SMESH_Delaunay.cxx
23 // Created   : Wed Apr 19 15:41:15 2017
24 // Author    : Edward AGAPOV (eap)
25
26 #include "SMESH_Delaunay.hxx"
27
28 #include "SMESH_Comment.hxx"
29 #include "SMESH_File.hxx"
30 #include "SMESH_MeshAlgos.hxx"
31
32 #include <BRepAdaptor_Surface.hxx>
33 #include <BRepMesh_Delaun.hxx>
34
35 #include <Basics_OCCTVersion.hxx>
36
37 //================================================================================
38 /*!
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
43  */
44 //================================================================================
45
46 SMESH_Delaunay::SMESH_Delaunay(const std::vector< const UVPtStructVec* > & boundaryNodes,
47                                const TopoDS_Face&                          face,
48                                const int                                   faceID)
49   : _face( face ), _faceID( faceID ), _scale( 1., 1. )
50 {
51   // compute _scale
52   BRepAdaptor_Surface surf( face );
53   if ( surf.GetType() != GeomAbs_Plane )
54   {
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 )
66     {
67       gp_Pnt p1U = surf.Value( u, vFixed );
68       lenU += p1U.Distance( p0U );
69       p0U = p1U;
70       gp_Pnt p1V = surf.Value( uFixed, v );
71       lenV += p1V.Distance( p0V );
72       p0V = p1V;
73     }
74     _scale.SetCoord( lenU / uRange, lenV / vRange );
75   }
76
77   // count boundary points
78   int iP = 1, nbP = 0;
79   for ( size_t iW = 0; iW < boundaryNodes.size(); ++iW ) // loop on wires
80   {
81     nbP += boundaryNodes[iW]->size();
82     if ( boundaryNodes[iW]->front().node == boundaryNodes[iW]->back().node )
83       --nbP; // 1st and last points coincide
84   }
85   _bndNodes.resize( nbP );
86
87   // fill boundary points
88 #if OCC_VERSION_LARGE <= 0x07030000
89   BRepMesh::Array1OfVertexOfDelaun bndVert( 1, 1 + nbP );
90 #else
91   IMeshData::Array1OfVertexOfDelaun bndVert( 1, 1 + nbP );
92 #endif
93   BRepMesh_Vertex v( 0, 0, BRepMesh_Frontier );
94   for ( size_t iW = 0; iW < boundaryNodes.size(); ++iW )
95   {
96     const UVPtStructVec& bndPnt = *boundaryNodes[iW];
97     int i = 0, nb = bndPnt.size();
98     if ( bndPnt[0].node == bndPnt.back().node )
99       --nb;
100     for ( ; i < nb;  ++i, ++iP )
101     {
102       _bndNodes[ iP-1 ] = bndPnt[i].node;
103       bndPnt[i].node->setIsMarked( true );
104
105       v.ChangeCoord() = bndPnt[i].UV().Multiplied( _scale );
106       bndVert( iP )   = v;
107     }
108   }
109
110   // triangulate the srcFace in 2D
111   BRepMesh_Delaun Delaunay( bndVert );
112   _triaDS = Delaunay.Result();
113 }
114
115 //================================================================================
116 /*!
117  * \brief Prepare to the exploration of nodes
118  */
119 //================================================================================
120
121 void SMESH_Delaunay::InitTraversal(const int nbNodesToVisit)
122 {
123   _nbNodesToVisit = (size_t) nbNodesToVisit;
124   _nbVisitedNodes = _iBndNode = 0;
125   _noTriQueue.clear();
126 }
127
128 //================================================================================
129 /*!
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
135  */
136 //================================================================================
137
138 const SMDS_MeshNode* SMESH_Delaunay::NextNode( double bc[3], int triaNodes[3] )
139 {
140   while ( _nbVisitedNodes < _nbNodesToVisit )
141   {
142     while ( !_noTriQueue.empty() )
143     {
144       const SMDS_MeshNode*     node = _noTriQueue.front().first;
145       const BRepMesh_Triangle* tria = _noTriQueue.front().second;
146       _noTriQueue.pop_front();
147       if ( node->isMarked() )
148         continue;
149       ++_nbVisitedNodes;
150       node->setIsMarked( true );
151
152       // find a Delaunay triangle containing the src node
153       gp_XY uv = getNodeUV( _face, node );
154       tria = FindTriangle( uv, tria, bc, triaNodes );
155       if ( tria )
156       {
157         addCloseNodes( node, tria, _faceID, _noTriQueue );
158         return node;
159       }
160     }
161     for ( ; _iBndNode < _bndNodes.size() &&  _noTriQueue.empty();  ++_iBndNode )
162     {
163       if ( const BRepMesh_Triangle* tria = GetTriangleNear( _iBndNode ))
164         addCloseNodes( _bndNodes[ _iBndNode ], tria, _faceID, _noTriQueue );
165     }
166     if ( _noTriQueue.empty() )
167       break;
168   }
169
170   // if ( _nbVisitedNodes < _nbNodesToVisit )
171   //   _nbVisitedNodes = std::numeric_limits<int>::max();
172   return NULL;
173 }
174
175 //================================================================================
176 /*!
177  * \brief Find a Delaunay triangle containing a given 2D point and return
178  *        barycentric coordinates within the found triangle
179  */
180 //================================================================================
181
182 const BRepMesh_Triangle* SMESH_Delaunay::FindTriangle( const gp_XY&             UV,
183                                                        const BRepMesh_Triangle* tria,
184                                                        double                   bc[3],
185                                                        int                      triaNodes[3] )
186 {
187   int   nodeIDs[3];
188   gp_XY nodeUVs[3];
189   int   linkIDs[3];
190   Standard_Boolean ori[3];
191
192   gp_XY uv = UV.Multiplied( _scale );
193
194   while ( tria )
195   {
196     // check if the uv is in tria
197
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();
202
203     SMESH_MeshAlgos::GetBarycentricCoords( uv,
204                                            nodeUVs[0], nodeUVs[1], nodeUVs[2],
205                                            bc[0], bc[1] );
206     if ( bc[0] >= 0 && bc[1] >= 0 && bc[0] + bc[1] <= 1 )
207     {
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 )
211       {
212         return 0;
213       }
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;
218       return tria;
219     }
220
221     // look for a neighbor triangle, which is adjacent to a link intersected
222     // by a segment( triangle center -> uv )
223
224     gp_XY gc = ( nodeUVs[0] + nodeUVs[1] + nodeUVs[2] ) / 3.;
225     gp_XY seg = uv - gc;
226
227     tria->Edges( linkIDs, ori );
228
229     const BRepMesh_Triangle* prevTria = tria;
230     tria = 0;
231
232     for ( int i = 0; i < 3; ++i )
233     {
234       const BRepMesh_PairOfIndex & triIDs = _triaDS->ElementsConnectedTo( linkIDs[i] );
235       if ( triIDs.Extent() < 2 )
236         continue; // no neighbor triangle
237
238       // check if a link intersects gc2uv
239       const BRepMesh_Edge & link = _triaDS->GetLink( linkIDs[i] );
240       const BRepMesh_Vertex & n1 = _triaDS->GetNode( link.FirstNode() );
241       const BRepMesh_Vertex & n2 = _triaDS->GetNode( link.LastNode() );
242       gp_XY uv1 = n1.Coord();
243       gp_XY lin = n2.Coord() - uv1; // link direction
244
245       double crossSegLin = seg ^ lin;
246       if ( Abs( crossSegLin ) < std::numeric_limits<double>::min() )
247         continue; // parallel
248
249       double uSeg = ( uv1 - gc ) ^ lin / crossSegLin;
250       if ( 0. <= uSeg && uSeg <= 1. )
251       {
252         tria = & _triaDS->GetElement( triIDs.Index( 1 ));
253         if ( tria == prevTria )
254           tria = & _triaDS->GetElement( triIDs.Index( 2 ));
255         if ( tria->Movability() != BRepMesh_Deleted )
256           break;
257       }
258     }
259   }
260   return tria;
261 }
262
263 //================================================================================
264 /*!
265  * \brief Return a triangle sharing a given boundary node
266  *  \param [in] iBndNode - index of the boundary node
267  *  \return const BRepMesh_Triangle* - a found triangle
268  */
269 //================================================================================
270
271 const BRepMesh_Triangle* SMESH_Delaunay::GetTriangleNear( int iBndNode )
272 {
273   if ( iBndNode >= _triaDS->NbNodes() )
274     return 0;
275   int nodeIDs[3];
276   int nbNbNodes = _bndNodes.size();
277 #if OCC_VERSION_LARGE <= 0x07030000
278   typedef BRepMesh::ListOfInteger TLinkList;
279 #else
280   typedef IMeshData::ListOfInteger TLinkList;
281 #endif
282   const TLinkList &       linkIds = _triaDS->LinksConnectedTo( iBndNode + 1 );
283   TLinkList::const_iterator iLink = linkIds.cbegin();
284   for ( ; iLink != linkIds.cend(); ++iLink )
285   {
286     const BRepMesh_PairOfIndex & triaIds = _triaDS->ElementsConnectedTo( *iLink );
287     {
288       const BRepMesh_Triangle& tria = _triaDS->GetElement( triaIds.Index(1) );
289       if ( tria.Movability() != BRepMesh_Deleted )
290       {
291         _triaDS->ElementNodes( tria, nodeIDs );
292         if ( nodeIDs[0]-1 < nbNbNodes &&
293              nodeIDs[1]-1 < nbNbNodes &&
294              nodeIDs[2]-1 < nbNbNodes )
295           return &tria;
296       }
297     }
298     if ( triaIds.Extent() > 1 )
299     {
300       const BRepMesh_Triangle& tria = _triaDS->GetElement( triaIds.Index(2) );
301       if ( tria.Movability() != BRepMesh_Deleted )
302       {
303         _triaDS->ElementNodes( tria, nodeIDs );
304         if ( nodeIDs[0]-1 < nbNbNodes &&
305              nodeIDs[1]-1 < nbNbNodes &&
306              nodeIDs[2]-1 < nbNbNodes )
307           return &tria;
308       }
309     }
310   }
311   return 0;
312 }
313
314 //================================================================================
315 /*!
316  * \brief Return UV of the i-th source boundary node (zero based)
317  */
318 //================================================================================
319
320 gp_XY SMESH_Delaunay::GetBndUV(const int iNode) const
321 {
322   return _triaDS->GetNode( iNode+1 ).Coord();
323 }
324
325 //================================================================================
326 /*!
327  * \brief Add non-marked nodes surrounding a given one to a queue
328  */
329 //================================================================================
330
331 void SMESH_Delaunay::addCloseNodes( const SMDS_MeshNode*     node,
332                                     const BRepMesh_Triangle* tria,
333                                     const int                faceID,
334                                     TNodeTriaList &          _noTriQueue )
335 {
336   // find in-FACE nodes
337   SMDS_ElemIteratorPtr elems = node->GetInverseElementIterator(SMDSAbs_Face);
338   while ( elems->more() )
339   {
340     const SMDS_MeshElement* elem = elems->next();
341     if ( elem->getshapeId() == faceID )
342     {
343       for ( int i = 0, nb = elem->NbNodes(); i < nb; ++i )
344       {
345         const SMDS_MeshNode* n = elem->GetNode( i );
346         if ( !n->isMarked() /*&& n->getshapeId() == faceID*/ )
347           _noTriQueue.push_back( std::make_pair( n, tria ));
348       }
349     }
350   }
351 }
352
353 //================================================================================
354 /*!
355  * \brief Write a python script that creates an equal mesh in Mesh module
356  */
357 //================================================================================
358
359 void SMESH_Delaunay::ToPython() const
360 {
361   SMESH_Comment text;
362   text << "import salome, SMESH\n";
363   text << "salome.salome_init()\n";
364   text << "from salome.smesh import smeshBuilder\n";
365   text << "smesh = smeshBuilder.New()\n";
366   text << "mesh=smesh.Mesh()\n";
367   const char* endl = "\n";
368
369   for ( int i = 0; i < _triaDS->NbNodes(); ++i )
370   {
371     const BRepMesh_Vertex& v = _triaDS->GetNode( i+1 );
372     text << "mesh.AddNode( " << v.Coord().X() << ", " << v.Coord().Y() << ", 0 )" << endl;
373   }
374
375   int nodeIDs[3];
376   for ( int i = 0; i < _triaDS->NbElements(); ++i )
377   {
378     const BRepMesh_Triangle& t = _triaDS->GetElement( i+1 );
379     if ( t.Movability() == BRepMesh_Deleted )
380       continue;
381     _triaDS->ElementNodes( t, nodeIDs );
382     text << "mesh.AddFace([ " << nodeIDs[0] << ", " << nodeIDs[1] << ", " << nodeIDs[2] << " ])" << endl;
383   }
384
385   const char* fileName = "/tmp/Delaunay.py";
386   SMESH_File file( fileName, false );
387   file.remove();
388   file.openForWriting();
389   file.write( text.c_str(), text.size() );
390   std::cout << "exec(open('" << fileName << "', 'rb').read())";
391 }