1 // SMESH SMESH : idl implementation based on 'SMESH' unit's classes
3 // Copyright (C) 2003 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.
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
24 // File : SMESH_MeshEditor.cxx
25 // Created : Mon Apr 12 16:10:22 2004
26 // Author : Edward AGAPOV (eap)
29 #include "SMESH_MeshEditor.hxx"
31 #include "SMDS_FaceOfNodes.hxx"
32 #include "SMDS_VolumeTool.hxx"
33 #include "SMDS_EdgePosition.hxx"
34 #include "SMDS_PolyhedralVolumeOfNodes.hxx"
35 #include "SMDS_FacePosition.hxx"
36 #include "SMDS_SpacePosition.hxx"
37 #include "SMDS_QuadraticFaceOfNodes.hxx"
38 #include "SMDS_MeshGroup.hxx"
40 #include "SMESHDS_Group.hxx"
41 #include "SMESHDS_Mesh.hxx"
43 #include "SMESH_subMesh.hxx"
44 #include "SMESH_ControlsDef.hxx"
45 #include "SMESH_MesherHelper.hxx"
46 #include "SMESH_OctreeNode.hxx"
47 #include "SMESH_Group.hxx"
49 #include "utilities.h"
51 #include <BRep_Tool.hxx>
53 #include <Extrema_GenExtPS.hxx>
54 #include <Extrema_POnSurf.hxx>
55 #include <Geom2d_Curve.hxx>
56 #include <GeomAdaptor_Surface.hxx>
57 #include <Geom_Curve.hxx>
58 #include <Geom_Surface.hxx>
59 #include <TColStd_ListOfInteger.hxx>
61 #include <TopExp_Explorer.hxx>
62 #include <TopTools_ListIteratorOfListOfShape.hxx>
63 #include <TopTools_ListOfShape.hxx>
65 #include <TopoDS_Face.hxx>
71 #include <gp_Trsf.hxx>
80 #define cast2Node(elem) static_cast<const SMDS_MeshNode*>( elem )
83 using namespace SMESH::Controls;
85 typedef map<const SMDS_MeshElement*, list<const SMDS_MeshNode*> > TElemOfNodeListMap;
86 typedef map<const SMDS_MeshElement*, list<const SMDS_MeshElement*> > TElemOfElemListMap;
87 //typedef map<const SMDS_MeshNode*, vector<const SMDS_MeshNode*> > TNodeOfNodeVecMap;
88 //typedef TNodeOfNodeVecMap::iterator TNodeOfNodeVecMapItr;
89 //typedef map<const SMDS_MeshElement*, vector<TNodeOfNodeVecMapItr> > TElemOfVecOfMapNodesMap;
91 struct TNodeXYZ : public gp_XYZ {
92 TNodeXYZ( const SMDS_MeshNode* n ):gp_XYZ( n->X(), n->Y(), n->Z() ) {}
95 //=======================================================================
96 //function : SMESH_MeshEditor
98 //=======================================================================
100 SMESH_MeshEditor::SMESH_MeshEditor( SMESH_Mesh* theMesh )
101 :myMesh( theMesh ) // theMesh may be NULL
105 //=======================================================================
109 //=======================================================================
112 SMESH_MeshEditor::AddElement(const vector<const SMDS_MeshNode*> & node,
113 const SMDSAbs_ElementType type,
117 SMDS_MeshElement* e = 0;
118 int nbnode = node.size();
119 SMESHDS_Mesh* mesh = GetMeshDS();
123 if ( ID ) e = mesh->AddEdgeWithID(node[0], node[1], ID);
124 else e = mesh->AddEdge (node[0], node[1] );
125 else if ( nbnode == 3 )
126 if ( ID ) e = mesh->AddEdgeWithID(node[0], node[1], node[2], ID);
127 else e = mesh->AddEdge (node[0], node[1], node[2] );
132 if ( ID ) e = mesh->AddFaceWithID(node[0], node[1], node[2], ID);
133 else e = mesh->AddFace (node[0], node[1], node[2] );
134 else if (nbnode == 4)
135 if ( ID ) e = mesh->AddFaceWithID(node[0], node[1], node[2], node[3], ID);
136 else e = mesh->AddFace (node[0], node[1], node[2], node[3] );
137 else if (nbnode == 6)
138 if ( ID ) e = mesh->AddFaceWithID(node[0], node[1], node[2], node[3],
139 node[4], node[5], ID);
140 else e = mesh->AddFace (node[0], node[1], node[2], node[3],
142 else if (nbnode == 8)
143 if ( ID ) e = mesh->AddFaceWithID(node[0], node[1], node[2], node[3],
144 node[4], node[5], node[6], node[7], ID);
145 else e = mesh->AddFace (node[0], node[1], node[2], node[3],
146 node[4], node[5], node[6], node[7] );
148 if ( ID ) e = mesh->AddPolygonalFaceWithID(node, ID);
149 else e = mesh->AddPolygonalFace (node );
155 if ( ID ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3], ID);
156 else e = mesh->AddVolume (node[0], node[1], node[2], node[3] );
157 else if (nbnode == 5)
158 if ( ID ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3],
160 else e = mesh->AddVolume (node[0], node[1], node[2], node[3],
162 else if (nbnode == 6)
163 if ( ID ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3],
164 node[4], node[5], ID);
165 else e = mesh->AddVolume (node[0], node[1], node[2], node[3],
167 else if (nbnode == 8)
168 if ( ID ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3],
169 node[4], node[5], node[6], node[7], ID);
170 else e = mesh->AddVolume (node[0], node[1], node[2], node[3],
171 node[4], node[5], node[6], node[7] );
172 else if (nbnode == 10)
173 if ( ID ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3],
174 node[4], node[5], node[6], node[7],
175 node[8], node[9], ID);
176 else e = mesh->AddVolume (node[0], node[1], node[2], node[3],
177 node[4], node[5], node[6], node[7],
179 else if (nbnode == 13)
180 if ( ID ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3],
181 node[4], node[5], node[6], node[7],
182 node[8], node[9], node[10],node[11],
184 else e = mesh->AddVolume (node[0], node[1], node[2], node[3],
185 node[4], node[5], node[6], node[7],
186 node[8], node[9], node[10],node[11],
188 else if (nbnode == 15)
189 if ( ID ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3],
190 node[4], node[5], node[6], node[7],
191 node[8], node[9], node[10],node[11],
192 node[12],node[13],node[14],ID);
193 else e = mesh->AddVolume (node[0], node[1], node[2], node[3],
194 node[4], node[5], node[6], node[7],
195 node[8], node[9], node[10],node[11],
196 node[12],node[13],node[14] );
197 else if (nbnode == 20)
198 if ( ID ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3],
199 node[4], node[5], node[6], node[7],
200 node[8], node[9], node[10],node[11],
201 node[12],node[13],node[14],node[15],
202 node[16],node[17],node[18],node[19],ID);
203 else e = mesh->AddVolume (node[0], node[1], node[2], node[3],
204 node[4], node[5], node[6], node[7],
205 node[8], node[9], node[10],node[11],
206 node[12],node[13],node[14],node[15],
207 node[16],node[17],node[18],node[19] );
213 //=======================================================================
217 //=======================================================================
219 SMDS_MeshElement* SMESH_MeshEditor::AddElement(const vector<int> & nodeIDs,
220 const SMDSAbs_ElementType type,
224 vector<const SMDS_MeshNode*> nodes;
225 nodes.reserve( nodeIDs.size() );
226 vector<int>::const_iterator id = nodeIDs.begin();
227 while ( id != nodeIDs.end() ) {
228 if ( const SMDS_MeshNode* node = GetMeshDS()->FindNode( *id++ ))
229 nodes.push_back( node );
233 return AddElement( nodes, type, isPoly, ID );
236 //=======================================================================
238 //purpose : Remove a node or an element.
239 // Modify a compute state of sub-meshes which become empty
240 //=======================================================================
242 bool SMESH_MeshEditor::Remove (const list< int >& theIDs,
245 myLastCreatedElems.Clear();
246 myLastCreatedNodes.Clear();
248 SMESHDS_Mesh* aMesh = GetMeshDS();
249 set< SMESH_subMesh *> smmap;
251 list<int>::const_iterator it = theIDs.begin();
252 for ( ; it != theIDs.end(); it++ ) {
253 const SMDS_MeshElement * elem;
255 elem = aMesh->FindNode( *it );
257 elem = aMesh->FindElement( *it );
261 // Notify VERTEX sub-meshes about modification
263 const SMDS_MeshNode* node = cast2Node( elem );
264 if ( node->GetPosition()->GetTypeOfPosition() == SMDS_TOP_VERTEX )
265 if ( int aShapeID = node->GetPosition()->GetShapeId() )
266 if ( SMESH_subMesh * sm = GetMesh()->GetSubMeshContaining( aShapeID ) )
269 // Find sub-meshes to notify about modification
270 // SMDS_ElemIteratorPtr nodeIt = elem->nodesIterator();
271 // while ( nodeIt->more() ) {
272 // const SMDS_MeshNode* node = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
273 // const SMDS_PositionPtr& aPosition = node->GetPosition();
274 // if ( aPosition.get() ) {
275 // if ( int aShapeID = aPosition->GetShapeId() ) {
276 // if ( SMESH_subMesh * sm = GetMesh()->GetSubMeshContaining( aShapeID ) )
277 // smmap.insert( sm );
284 aMesh->RemoveNode( static_cast< const SMDS_MeshNode* >( elem ));
286 aMesh->RemoveElement( elem );
289 // Notify sub-meshes about modification
290 if ( !smmap.empty() ) {
291 set< SMESH_subMesh *>::iterator smIt;
292 for ( smIt = smmap.begin(); smIt != smmap.end(); smIt++ )
293 (*smIt)->ComputeStateEngine( SMESH_subMesh::MESH_ENTITY_REMOVED );
296 // // Check if the whole mesh becomes empty
297 // if ( SMESH_subMesh * sm = GetMesh()->GetSubMeshContaining( 1 ) )
298 // sm->ComputeStateEngine( SMESH_subMesh::CHECK_COMPUTE_STATE );
303 //=======================================================================
304 //function : FindShape
305 //purpose : Return an index of the shape theElem is on
306 // or zero if a shape not found
307 //=======================================================================
309 int SMESH_MeshEditor::FindShape (const SMDS_MeshElement * theElem)
311 myLastCreatedElems.Clear();
312 myLastCreatedNodes.Clear();
314 SMESHDS_Mesh * aMesh = GetMeshDS();
315 if ( aMesh->ShapeToMesh().IsNull() )
318 if ( theElem->GetType() == SMDSAbs_Node ) {
319 const SMDS_PositionPtr& aPosition =
320 static_cast<const SMDS_MeshNode*>( theElem )->GetPosition();
321 if ( aPosition.get() )
322 return aPosition->GetShapeId();
327 TopoDS_Shape aShape; // the shape a node is on
328 SMDS_ElemIteratorPtr nodeIt = theElem->nodesIterator();
329 while ( nodeIt->more() ) {
330 const SMDS_MeshNode* node = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
331 const SMDS_PositionPtr& aPosition = node->GetPosition();
332 if ( aPosition.get() ) {
333 int aShapeID = aPosition->GetShapeId();
334 SMESHDS_SubMesh * sm = aMesh->MeshElements( aShapeID );
336 if ( sm->Contains( theElem ))
338 if ( aShape.IsNull() )
339 aShape = aMesh->IndexToShape( aShapeID );
342 //MESSAGE ( "::FindShape() No SubShape for aShapeID " << aShapeID );
347 // None of nodes is on a proper shape,
348 // find the shape among ancestors of aShape on which a node is
349 if ( aShape.IsNull() ) {
350 //MESSAGE ("::FindShape() - NONE node is on shape")
353 TopTools_ListIteratorOfListOfShape ancIt( GetMesh()->GetAncestors( aShape ));
354 for ( ; ancIt.More(); ancIt.Next() ) {
355 SMESHDS_SubMesh * sm = aMesh->MeshElements( ancIt.Value() );
356 if ( sm && sm->Contains( theElem ))
357 return aMesh->ShapeToIndex( ancIt.Value() );
360 //MESSAGE ("::FindShape() - SHAPE NOT FOUND")
364 //=======================================================================
365 //function : IsMedium
367 //=======================================================================
369 bool SMESH_MeshEditor::IsMedium(const SMDS_MeshNode* node,
370 const SMDSAbs_ElementType typeToCheck)
372 bool isMedium = false;
373 SMDS_ElemIteratorPtr it = node->GetInverseElementIterator(typeToCheck);
374 while (it->more() && !isMedium ) {
375 const SMDS_MeshElement* elem = it->next();
376 isMedium = elem->IsMediumNode(node);
381 //=======================================================================
382 //function : ShiftNodesQuadTria
384 // Shift nodes in the array corresponded to quadratic triangle
385 // example: (0,1,2,3,4,5) -> (1,2,0,4,5,3)
386 //=======================================================================
387 static void ShiftNodesQuadTria(const SMDS_MeshNode* aNodes[])
389 const SMDS_MeshNode* nd1 = aNodes[0];
390 aNodes[0] = aNodes[1];
391 aNodes[1] = aNodes[2];
393 const SMDS_MeshNode* nd2 = aNodes[3];
394 aNodes[3] = aNodes[4];
395 aNodes[4] = aNodes[5];
399 //=======================================================================
400 //function : GetNodesFromTwoTria
402 // Shift nodes in the array corresponded to quadratic triangle
403 // example: (0,1,2,3,4,5) -> (1,2,0,4,5,3)
404 //=======================================================================
405 static bool GetNodesFromTwoTria(const SMDS_MeshElement * theTria1,
406 const SMDS_MeshElement * theTria2,
407 const SMDS_MeshNode* N1[],
408 const SMDS_MeshNode* N2[])
410 SMDS_ElemIteratorPtr it = theTria1->nodesIterator();
413 N1[i] = static_cast<const SMDS_MeshNode*>( it->next() );
416 if(it->more()) return false;
417 it = theTria2->nodesIterator();
420 N2[i] = static_cast<const SMDS_MeshNode*>( it->next() );
423 if(it->more()) return false;
425 int sames[3] = {-1,-1,-1};
437 if(nbsames!=2) return false;
439 ShiftNodesQuadTria(N1);
441 ShiftNodesQuadTria(N1);
444 i = sames[0] + sames[1] + sames[2];
446 ShiftNodesQuadTria(N2);
448 // now we receive following N1 and N2 (using numeration as above image)
449 // tria1 : (1 2 4 5 9 7) and tria2 : (3 4 2 8 9 6)
450 // i.e. first nodes from both arrays determ new diagonal
454 //=======================================================================
455 //function : InverseDiag
456 //purpose : Replace two neighbour triangles with ones built on the same 4 nodes
457 // but having other common link.
458 // Return False if args are improper
459 //=======================================================================
461 bool SMESH_MeshEditor::InverseDiag (const SMDS_MeshElement * theTria1,
462 const SMDS_MeshElement * theTria2 )
464 myLastCreatedElems.Clear();
465 myLastCreatedNodes.Clear();
467 if (!theTria1 || !theTria2)
470 const SMDS_FaceOfNodes* F1 = dynamic_cast<const SMDS_FaceOfNodes*>( theTria1 );
471 const SMDS_FaceOfNodes* F2 = dynamic_cast<const SMDS_FaceOfNodes*>( theTria2 );
474 // 1 +--+ A theTria1: ( 1 A B ) A->2 ( 1 2 B ) 1 +--+ A
475 // | /| theTria2: ( B A 2 ) B->1 ( 1 A 2 ) |\ |
479 // put nodes in array and find out indices of the same ones
480 const SMDS_MeshNode* aNodes [6];
481 int sameInd [] = { 0, 0, 0, 0, 0, 0 };
483 SMDS_ElemIteratorPtr it = theTria1->nodesIterator();
484 while ( it->more() ) {
485 aNodes[ i ] = static_cast<const SMDS_MeshNode*>( it->next() );
487 if ( i > 2 ) // theTria2
488 // find same node of theTria1
489 for ( int j = 0; j < 3; j++ )
490 if ( aNodes[ i ] == aNodes[ j ]) {
499 return false; // theTria1 is not a triangle
500 it = theTria2->nodesIterator();
502 if ( i == 6 && it->more() )
503 return false; // theTria2 is not a triangle
506 // find indices of 1,2 and of A,B in theTria1
507 int iA = 0, iB = 0, i1 = 0, i2 = 0;
508 for ( i = 0; i < 6; i++ ) {
509 if ( sameInd [ i ] == 0 )
516 // nodes 1 and 2 should not be the same
517 if ( aNodes[ i1 ] == aNodes[ i2 ] )
521 aNodes[ iA ] = aNodes[ i2 ];
523 aNodes[ sameInd[ iB ]] = aNodes[ i1 ];
525 //MESSAGE( theTria1 << theTria2 );
527 GetMeshDS()->ChangeElementNodes( theTria1, aNodes, 3 );
528 GetMeshDS()->ChangeElementNodes( theTria2, &aNodes[ 3 ], 3 );
530 //MESSAGE( theTria1 << theTria2 );
534 } // end if(F1 && F2)
536 // check case of quadratic faces
537 const SMDS_QuadraticFaceOfNodes* QF1 =
538 dynamic_cast<const SMDS_QuadraticFaceOfNodes*> (theTria1);
539 if(!QF1) return false;
540 const SMDS_QuadraticFaceOfNodes* QF2 =
541 dynamic_cast<const SMDS_QuadraticFaceOfNodes*> (theTria2);
542 if(!QF2) return false;
545 // 1 +--+--+ 2 theTria1: (1 2 4 5 9 7) or (2 4 1 9 7 5) or (4 1 2 7 5 9)
546 // | /| theTria2: (2 3 4 6 8 9) or (3 4 2 8 9 6) or (4 2 3 9 6 8)
554 const SMDS_MeshNode* N1 [6];
555 const SMDS_MeshNode* N2 [6];
556 if(!GetNodesFromTwoTria(theTria1,theTria2,N1,N2))
558 // now we receive following N1 and N2 (using numeration as above image)
559 // tria1 : (1 2 4 5 9 7) and tria2 : (3 4 2 8 9 6)
560 // i.e. first nodes from both arrays determ new diagonal
562 const SMDS_MeshNode* N1new [6];
563 const SMDS_MeshNode* N2new [6];
576 // replaces nodes in faces
577 GetMeshDS()->ChangeElementNodes( theTria1, N1new, 6 );
578 GetMeshDS()->ChangeElementNodes( theTria2, N2new, 6 );
583 //=======================================================================
584 //function : findTriangles
585 //purpose : find triangles sharing theNode1-theNode2 link
586 //=======================================================================
588 static bool findTriangles(const SMDS_MeshNode * theNode1,
589 const SMDS_MeshNode * theNode2,
590 const SMDS_MeshElement*& theTria1,
591 const SMDS_MeshElement*& theTria2)
593 if ( !theNode1 || !theNode2 ) return false;
595 theTria1 = theTria2 = 0;
597 set< const SMDS_MeshElement* > emap;
598 SMDS_ElemIteratorPtr it = theNode1->GetInverseElementIterator(SMDSAbs_Face);
600 const SMDS_MeshElement* elem = it->next();
601 if ( elem->NbNodes() == 3 )
604 it = theNode2->GetInverseElementIterator(SMDSAbs_Face);
606 const SMDS_MeshElement* elem = it->next();
607 if ( emap.find( elem ) != emap.end() )
609 // theTria1 must be element with minimum ID
610 if( theTria1->GetID() < elem->GetID() ) {
623 return ( theTria1 && theTria2 );
626 //=======================================================================
627 //function : InverseDiag
628 //purpose : Replace two neighbour triangles sharing theNode1-theNode2 link
629 // with ones built on the same 4 nodes but having other common link.
630 // Return false if proper faces not found
631 //=======================================================================
633 bool SMESH_MeshEditor::InverseDiag (const SMDS_MeshNode * theNode1,
634 const SMDS_MeshNode * theNode2)
636 myLastCreatedElems.Clear();
637 myLastCreatedNodes.Clear();
639 MESSAGE( "::InverseDiag()" );
641 const SMDS_MeshElement *tr1, *tr2;
642 if ( !findTriangles( theNode1, theNode2, tr1, tr2 ))
645 const SMDS_FaceOfNodes* F1 = dynamic_cast<const SMDS_FaceOfNodes*>( tr1 );
646 //if (!F1) return false;
647 const SMDS_FaceOfNodes* F2 = dynamic_cast<const SMDS_FaceOfNodes*>( tr2 );
648 //if (!F2) return false;
651 // 1 +--+ A tr1: ( 1 A B ) A->2 ( 1 2 B ) 1 +--+ A
652 // | /| tr2: ( B A 2 ) B->1 ( 1 A 2 ) |\ |
656 // put nodes in array
657 // and find indices of 1,2 and of A in tr1 and of B in tr2
658 int i, iA1 = 0, i1 = 0;
659 const SMDS_MeshNode* aNodes1 [3];
660 SMDS_ElemIteratorPtr it;
661 for (i = 0, it = tr1->nodesIterator(); it->more(); i++ ) {
662 aNodes1[ i ] = static_cast<const SMDS_MeshNode*>( it->next() );
663 if ( aNodes1[ i ] == theNode1 )
664 iA1 = i; // node A in tr1
665 else if ( aNodes1[ i ] != theNode2 )
669 const SMDS_MeshNode* aNodes2 [3];
670 for (i = 0, it = tr2->nodesIterator(); it->more(); i++ ) {
671 aNodes2[ i ] = static_cast<const SMDS_MeshNode*>( it->next() );
672 if ( aNodes2[ i ] == theNode2 )
673 iB2 = i; // node B in tr2
674 else if ( aNodes2[ i ] != theNode1 )
678 // nodes 1 and 2 should not be the same
679 if ( aNodes1[ i1 ] == aNodes2[ i2 ] )
683 aNodes1[ iA1 ] = aNodes2[ i2 ];
685 aNodes2[ iB2 ] = aNodes1[ i1 ];
687 //MESSAGE( tr1 << tr2 );
689 GetMeshDS()->ChangeElementNodes( tr1, aNodes1, 3 );
690 GetMeshDS()->ChangeElementNodes( tr2, aNodes2, 3 );
692 //MESSAGE( tr1 << tr2 );
697 // check case of quadratic faces
698 const SMDS_QuadraticFaceOfNodes* QF1 =
699 dynamic_cast<const SMDS_QuadraticFaceOfNodes*> (tr1);
700 if(!QF1) return false;
701 const SMDS_QuadraticFaceOfNodes* QF2 =
702 dynamic_cast<const SMDS_QuadraticFaceOfNodes*> (tr2);
703 if(!QF2) return false;
704 return InverseDiag(tr1,tr2);
707 //=======================================================================
708 //function : getQuadrangleNodes
709 //purpose : fill theQuadNodes - nodes of a quadrangle resulting from
710 // fusion of triangles tr1 and tr2 having shared link on
711 // theNode1 and theNode2
712 //=======================================================================
714 bool getQuadrangleNodes(const SMDS_MeshNode * theQuadNodes [],
715 const SMDS_MeshNode * theNode1,
716 const SMDS_MeshNode * theNode2,
717 const SMDS_MeshElement * tr1,
718 const SMDS_MeshElement * tr2 )
720 if( tr1->NbNodes() != tr2->NbNodes() )
722 // find the 4-th node to insert into tr1
723 const SMDS_MeshNode* n4 = 0;
724 SMDS_ElemIteratorPtr it = tr2->nodesIterator();
726 while ( !n4 && i<3 ) {
727 const SMDS_MeshNode * n = cast2Node( it->next() );
729 bool isDiag = ( n == theNode1 || n == theNode2 );
733 // Make an array of nodes to be in a quadrangle
734 int iNode = 0, iFirstDiag = -1;
735 it = tr1->nodesIterator();
738 const SMDS_MeshNode * n = cast2Node( it->next() );
740 bool isDiag = ( n == theNode1 || n == theNode2 );
742 if ( iFirstDiag < 0 )
744 else if ( iNode - iFirstDiag == 1 )
745 theQuadNodes[ iNode++ ] = n4; // insert the 4-th node between diagonal nodes
747 else if ( n == n4 ) {
748 return false; // tr1 and tr2 should not have all the same nodes
750 theQuadNodes[ iNode++ ] = n;
752 if ( iNode == 3 ) // diagonal nodes have 0 and 2 indices
753 theQuadNodes[ iNode ] = n4;
758 //=======================================================================
759 //function : DeleteDiag
760 //purpose : Replace two neighbour triangles sharing theNode1-theNode2 link
761 // with a quadrangle built on the same 4 nodes.
762 // Return false if proper faces not found
763 //=======================================================================
765 bool SMESH_MeshEditor::DeleteDiag (const SMDS_MeshNode * theNode1,
766 const SMDS_MeshNode * theNode2)
768 myLastCreatedElems.Clear();
769 myLastCreatedNodes.Clear();
771 MESSAGE( "::DeleteDiag()" );
773 const SMDS_MeshElement *tr1, *tr2;
774 if ( !findTriangles( theNode1, theNode2, tr1, tr2 ))
777 const SMDS_FaceOfNodes* F1 = dynamic_cast<const SMDS_FaceOfNodes*>( tr1 );
778 //if (!F1) return false;
779 const SMDS_FaceOfNodes* F2 = dynamic_cast<const SMDS_FaceOfNodes*>( tr2 );
780 //if (!F2) return false;
783 const SMDS_MeshNode* aNodes [ 4 ];
784 if ( ! getQuadrangleNodes( aNodes, theNode1, theNode2, tr1, tr2 ))
787 //MESSAGE( endl << tr1 << tr2 );
789 GetMeshDS()->ChangeElementNodes( tr1, aNodes, 4 );
790 myLastCreatedElems.Append(tr1);
791 GetMeshDS()->RemoveElement( tr2 );
793 //MESSAGE( endl << tr1 );
798 // check case of quadratic faces
799 const SMDS_QuadraticFaceOfNodes* QF1 =
800 dynamic_cast<const SMDS_QuadraticFaceOfNodes*> (tr1);
801 if(!QF1) return false;
802 const SMDS_QuadraticFaceOfNodes* QF2 =
803 dynamic_cast<const SMDS_QuadraticFaceOfNodes*> (tr2);
804 if(!QF2) return false;
807 // 1 +--+--+ 2 tr1: (1 2 4 5 9 7) or (2 4 1 9 7 5) or (4 1 2 7 5 9)
808 // | /| tr2: (2 3 4 6 8 9) or (3 4 2 8 9 6) or (4 2 3 9 6 8)
816 const SMDS_MeshNode* N1 [6];
817 const SMDS_MeshNode* N2 [6];
818 if(!GetNodesFromTwoTria(tr1,tr2,N1,N2))
820 // now we receive following N1 and N2 (using numeration as above image)
821 // tria1 : (1 2 4 5 9 7) and tria2 : (3 4 2 8 9 6)
822 // i.e. first nodes from both arrays determ new diagonal
824 const SMDS_MeshNode* aNodes[8];
834 GetMeshDS()->ChangeElementNodes( tr1, aNodes, 8 );
835 myLastCreatedElems.Append(tr1);
836 GetMeshDS()->RemoveElement( tr2 );
838 // remove middle node (9)
839 GetMeshDS()->RemoveNode( N1[4] );
844 //=======================================================================
845 //function : Reorient
846 //purpose : Reverse theElement orientation
847 //=======================================================================
849 bool SMESH_MeshEditor::Reorient (const SMDS_MeshElement * theElem)
851 myLastCreatedElems.Clear();
852 myLastCreatedNodes.Clear();
856 SMDS_ElemIteratorPtr it = theElem->nodesIterator();
857 if ( !it || !it->more() )
860 switch ( theElem->GetType() ) {
864 if(!theElem->IsQuadratic()) {
865 int i = theElem->NbNodes();
866 vector<const SMDS_MeshNode*> aNodes( i );
868 aNodes[ --i ]= static_cast<const SMDS_MeshNode*>( it->next() );
869 return GetMeshDS()->ChangeElementNodes( theElem, &aNodes[0], theElem->NbNodes() );
872 // quadratic elements
873 if(theElem->GetType()==SMDSAbs_Edge) {
874 vector<const SMDS_MeshNode*> aNodes(3);
875 aNodes[1]= static_cast<const SMDS_MeshNode*>( it->next() );
876 aNodes[0]= static_cast<const SMDS_MeshNode*>( it->next() );
877 aNodes[2]= static_cast<const SMDS_MeshNode*>( it->next() );
878 return GetMeshDS()->ChangeElementNodes( theElem, &aNodes[0], 3 );
881 int nbn = theElem->NbNodes();
882 vector<const SMDS_MeshNode*> aNodes(nbn);
883 aNodes[0]= static_cast<const SMDS_MeshNode*>( it->next() );
885 for(; i<nbn/2; i++) {
886 aNodes[nbn/2-i]= static_cast<const SMDS_MeshNode*>( it->next() );
888 for(i=0; i<nbn/2; i++) {
889 aNodes[nbn-i-1]= static_cast<const SMDS_MeshNode*>( it->next() );
891 return GetMeshDS()->ChangeElementNodes( theElem, &aNodes[0], nbn );
895 case SMDSAbs_Volume: {
896 if (theElem->IsPoly()) {
897 const SMDS_PolyhedralVolumeOfNodes* aPolyedre =
898 static_cast<const SMDS_PolyhedralVolumeOfNodes*>( theElem );
900 MESSAGE("Warning: bad volumic element");
904 int nbFaces = aPolyedre->NbFaces();
905 vector<const SMDS_MeshNode *> poly_nodes;
906 vector<int> quantities (nbFaces);
908 // reverse each face of the polyedre
909 for (int iface = 1; iface <= nbFaces; iface++) {
910 int inode, nbFaceNodes = aPolyedre->NbFaceNodes(iface);
911 quantities[iface - 1] = nbFaceNodes;
913 for (inode = nbFaceNodes; inode >= 1; inode--) {
914 const SMDS_MeshNode* curNode = aPolyedre->GetFaceNode(iface, inode);
915 poly_nodes.push_back(curNode);
919 return GetMeshDS()->ChangePolyhedronNodes( theElem, poly_nodes, quantities );
923 SMDS_VolumeTool vTool;
924 if ( !vTool.Set( theElem ))
927 return GetMeshDS()->ChangeElementNodes( theElem, vTool.GetNodes(), vTool.NbNodes() );
936 //=======================================================================
937 //function : getBadRate
939 //=======================================================================
941 static double getBadRate (const SMDS_MeshElement* theElem,
942 SMESH::Controls::NumericalFunctorPtr& theCrit)
944 SMESH::Controls::TSequenceOfXYZ P;
945 if ( !theElem || !theCrit->GetPoints( theElem, P ))
947 return theCrit->GetBadRate( theCrit->GetValue( P ), theElem->NbNodes() );
948 //return theCrit->GetBadRate( theCrit->GetValue( theElem->GetID() ), theElem->NbNodes() );
951 //=======================================================================
952 //function : QuadToTri
953 //purpose : Cut quadrangles into triangles.
954 // theCrit is used to select a diagonal to cut
955 //=======================================================================
957 bool SMESH_MeshEditor::QuadToTri (TIDSortedElemSet & theElems,
958 SMESH::Controls::NumericalFunctorPtr theCrit)
960 myLastCreatedElems.Clear();
961 myLastCreatedNodes.Clear();
963 MESSAGE( "::QuadToTri()" );
965 if ( !theCrit.get() )
968 SMESHDS_Mesh * aMesh = GetMeshDS();
970 Handle(Geom_Surface) surface;
971 SMESH_MesherHelper helper( *GetMesh() );
973 TIDSortedElemSet::iterator itElem;
974 for ( itElem = theElems.begin(); itElem != theElems.end(); itElem++ ) {
975 const SMDS_MeshElement* elem = *itElem;
976 if ( !elem || elem->GetType() != SMDSAbs_Face )
978 if ( elem->NbNodes() != ( elem->IsQuadratic() ? 8 : 4 ))
981 // retrieve element nodes
982 const SMDS_MeshNode* aNodes [8];
983 SMDS_ElemIteratorPtr itN = elem->nodesIterator();
985 while ( itN->more() )
986 aNodes[ i++ ] = static_cast<const SMDS_MeshNode*>( itN->next() );
988 // compare two sets of possible triangles
989 double aBadRate1, aBadRate2; // to what extent a set is bad
990 SMDS_FaceOfNodes tr1 ( aNodes[0], aNodes[1], aNodes[2] );
991 SMDS_FaceOfNodes tr2 ( aNodes[2], aNodes[3], aNodes[0] );
992 aBadRate1 = getBadRate( &tr1, theCrit ) + getBadRate( &tr2, theCrit );
994 SMDS_FaceOfNodes tr3 ( aNodes[1], aNodes[2], aNodes[3] );
995 SMDS_FaceOfNodes tr4 ( aNodes[3], aNodes[0], aNodes[1] );
996 aBadRate2 = getBadRate( &tr3, theCrit ) + getBadRate( &tr4, theCrit );
998 int aShapeId = FindShape( elem );
999 const SMDS_MeshElement* newElem = 0;
1001 if( !elem->IsQuadratic() ) {
1003 // split liner quadrangle
1005 if ( aBadRate1 <= aBadRate2 ) {
1006 // tr1 + tr2 is better
1007 aMesh->ChangeElementNodes( elem, aNodes, 3 );
1008 newElem = aMesh->AddFace( aNodes[2], aNodes[3], aNodes[0] );
1011 // tr3 + tr4 is better
1012 aMesh->ChangeElementNodes( elem, &aNodes[1], 3 );
1013 newElem = aMesh->AddFace( aNodes[3], aNodes[0], aNodes[1] );
1018 // split quadratic quadrangle
1020 // get surface elem is on
1021 if ( aShapeId != helper.GetSubShapeID() ) {
1025 shape = aMesh->IndexToShape( aShapeId );
1026 if ( !shape.IsNull() && shape.ShapeType() == TopAbs_FACE ) {
1027 TopoDS_Face face = TopoDS::Face( shape );
1028 surface = BRep_Tool::Surface( face );
1029 if ( !surface.IsNull() )
1030 helper.SetSubShape( shape );
1034 const SMDS_MeshNode* aNodes [8];
1035 const SMDS_MeshNode* inFaceNode = 0;
1036 SMDS_ElemIteratorPtr itN = elem->nodesIterator();
1038 while ( itN->more() ) {
1039 aNodes[ i++ ] = static_cast<const SMDS_MeshNode*>( itN->next() );
1040 if ( !inFaceNode && helper.GetNodeUVneedInFaceNode() &&
1041 aNodes[ i-1 ]->GetPosition()->GetTypeOfPosition() == SMDS_TOP_FACE )
1043 inFaceNode = aNodes[ i-1 ];
1046 // find middle point for (0,1,2,3)
1047 // and create a node in this point;
1049 if ( surface.IsNull() ) {
1051 p += gp_XYZ(aNodes[i]->X(), aNodes[i]->Y(), aNodes[i]->Z() );
1055 TopoDS_Face face = TopoDS::Face( helper.GetSubShape() );
1058 uv += helper.GetNodeUV( face, aNodes[i], inFaceNode );
1060 p = surface->Value( uv.X(), uv.Y() ).XYZ();
1062 const SMDS_MeshNode* newN = aMesh->AddNode( p.X(), p.Y(), p.Z() );
1063 myLastCreatedNodes.Append(newN);
1065 // create a new element
1066 const SMDS_MeshNode* N[6];
1067 if ( aBadRate1 <= aBadRate2 ) {
1074 newElem = aMesh->AddFace(aNodes[2], aNodes[3], aNodes[0],
1075 aNodes[6], aNodes[7], newN );
1084 newElem = aMesh->AddFace(aNodes[3], aNodes[0], aNodes[1],
1085 aNodes[7], aNodes[4], newN );
1087 aMesh->ChangeElementNodes( elem, N, 6 );
1091 // care of a new element
1093 myLastCreatedElems.Append(newElem);
1094 AddToSameGroups( newElem, elem, aMesh );
1096 // put a new triangle on the same shape
1098 aMesh->SetMeshElementOnShape( newElem, aShapeId );
1103 //=======================================================================
1104 //function : BestSplit
1105 //purpose : Find better diagonal for cutting.
1106 //=======================================================================
1107 int SMESH_MeshEditor::BestSplit (const SMDS_MeshElement* theQuad,
1108 SMESH::Controls::NumericalFunctorPtr theCrit)
1110 myLastCreatedElems.Clear();
1111 myLastCreatedNodes.Clear();
1116 if (!theQuad || theQuad->GetType() != SMDSAbs_Face )
1119 if( theQuad->NbNodes()==4 ||
1120 (theQuad->NbNodes()==8 && theQuad->IsQuadratic()) ) {
1122 // retrieve element nodes
1123 const SMDS_MeshNode* aNodes [4];
1124 SMDS_ElemIteratorPtr itN = theQuad->nodesIterator();
1126 //while (itN->more())
1128 aNodes[ i++ ] = static_cast<const SMDS_MeshNode*>( itN->next() );
1130 // compare two sets of possible triangles
1131 double aBadRate1, aBadRate2; // to what extent a set is bad
1132 SMDS_FaceOfNodes tr1 ( aNodes[0], aNodes[1], aNodes[2] );
1133 SMDS_FaceOfNodes tr2 ( aNodes[2], aNodes[3], aNodes[0] );
1134 aBadRate1 = getBadRate( &tr1, theCrit ) + getBadRate( &tr2, theCrit );
1136 SMDS_FaceOfNodes tr3 ( aNodes[1], aNodes[2], aNodes[3] );
1137 SMDS_FaceOfNodes tr4 ( aNodes[3], aNodes[0], aNodes[1] );
1138 aBadRate2 = getBadRate( &tr3, theCrit ) + getBadRate( &tr4, theCrit );
1140 if (aBadRate1 <= aBadRate2) // tr1 + tr2 is better
1141 return 1; // diagonal 1-3
1143 return 2; // diagonal 2-4
1148 //=======================================================================
1149 //function : AddToSameGroups
1150 //purpose : add elemToAdd to the groups the elemInGroups belongs to
1151 //=======================================================================
1153 void SMESH_MeshEditor::AddToSameGroups (const SMDS_MeshElement* elemToAdd,
1154 const SMDS_MeshElement* elemInGroups,
1155 SMESHDS_Mesh * aMesh)
1157 const set<SMESHDS_GroupBase*>& groups = aMesh->GetGroups();
1158 set<SMESHDS_GroupBase*>::const_iterator grIt = groups.begin();
1159 for ( ; grIt != groups.end(); grIt++ ) {
1160 SMESHDS_Group* group = dynamic_cast<SMESHDS_Group*>( *grIt );
1161 if ( group && group->Contains( elemInGroups ))
1162 group->SMDSGroup().Add( elemToAdd );
1167 //=======================================================================
1168 //function : RemoveElemFromGroups
1169 //purpose : Remove removeelem to the groups the elemInGroups belongs to
1170 //=======================================================================
1171 void SMESH_MeshEditor::RemoveElemFromGroups (const SMDS_MeshElement* removeelem,
1172 SMESHDS_Mesh * aMesh)
1174 const set<SMESHDS_GroupBase*>& groups = aMesh->GetGroups();
1175 if (!groups.empty())
1177 set<SMESHDS_GroupBase*>::const_iterator GrIt = groups.begin();
1178 for (; GrIt != groups.end(); GrIt++)
1180 SMESHDS_Group* grp = dynamic_cast<SMESHDS_Group*>(*GrIt);
1181 if (!grp || grp->IsEmpty()) continue;
1182 grp->SMDSGroup().Remove(removeelem);
1188 //=======================================================================
1189 //function : QuadToTri
1190 //purpose : Cut quadrangles into triangles.
1191 // theCrit is used to select a diagonal to cut
1192 //=======================================================================
1194 bool SMESH_MeshEditor::QuadToTri (TIDSortedElemSet & theElems,
1195 const bool the13Diag)
1197 myLastCreatedElems.Clear();
1198 myLastCreatedNodes.Clear();
1200 MESSAGE( "::QuadToTri()" );
1202 SMESHDS_Mesh * aMesh = GetMeshDS();
1204 Handle(Geom_Surface) surface;
1205 SMESH_MesherHelper helper( *GetMesh() );
1207 TIDSortedElemSet::iterator itElem;
1208 for ( itElem = theElems.begin(); itElem != theElems.end(); itElem++ ) {
1209 const SMDS_MeshElement* elem = *itElem;
1210 if ( !elem || elem->GetType() != SMDSAbs_Face )
1212 bool isquad = elem->NbNodes()==4 || elem->NbNodes()==8;
1213 if(!isquad) continue;
1215 if(elem->NbNodes()==4) {
1216 // retrieve element nodes
1217 const SMDS_MeshNode* aNodes [4];
1218 SMDS_ElemIteratorPtr itN = elem->nodesIterator();
1220 while ( itN->more() )
1221 aNodes[ i++ ] = static_cast<const SMDS_MeshNode*>( itN->next() );
1223 int aShapeId = FindShape( elem );
1224 const SMDS_MeshElement* newElem = 0;
1226 aMesh->ChangeElementNodes( elem, aNodes, 3 );
1227 newElem = aMesh->AddFace( aNodes[2], aNodes[3], aNodes[0] );
1230 aMesh->ChangeElementNodes( elem, &aNodes[1], 3 );
1231 newElem = aMesh->AddFace( aNodes[3], aNodes[0], aNodes[1] );
1233 myLastCreatedElems.Append(newElem);
1234 // put a new triangle on the same shape and add to the same groups
1236 aMesh->SetMeshElementOnShape( newElem, aShapeId );
1237 AddToSameGroups( newElem, elem, aMesh );
1240 // Quadratic quadrangle
1242 if( elem->NbNodes()==8 && elem->IsQuadratic() ) {
1244 // get surface elem is on
1245 int aShapeId = FindShape( elem );
1246 if ( aShapeId != helper.GetSubShapeID() ) {
1250 shape = aMesh->IndexToShape( aShapeId );
1251 if ( !shape.IsNull() && shape.ShapeType() == TopAbs_FACE ) {
1252 TopoDS_Face face = TopoDS::Face( shape );
1253 surface = BRep_Tool::Surface( face );
1254 if ( !surface.IsNull() )
1255 helper.SetSubShape( shape );
1259 const SMDS_MeshNode* aNodes [8];
1260 const SMDS_MeshNode* inFaceNode = 0;
1261 SMDS_ElemIteratorPtr itN = elem->nodesIterator();
1263 while ( itN->more() ) {
1264 aNodes[ i++ ] = static_cast<const SMDS_MeshNode*>( itN->next() );
1265 if ( !inFaceNode && helper.GetNodeUVneedInFaceNode() &&
1266 aNodes[ i-1 ]->GetPosition()->GetTypeOfPosition() == SMDS_TOP_FACE )
1268 inFaceNode = aNodes[ i-1 ];
1272 // find middle point for (0,1,2,3)
1273 // and create a node in this point;
1275 if ( surface.IsNull() ) {
1277 p += gp_XYZ(aNodes[i]->X(), aNodes[i]->Y(), aNodes[i]->Z() );
1281 TopoDS_Face geomFace = TopoDS::Face( helper.GetSubShape() );
1284 uv += helper.GetNodeUV( geomFace, aNodes[i], inFaceNode );
1286 p = surface->Value( uv.X(), uv.Y() ).XYZ();
1288 const SMDS_MeshNode* newN = aMesh->AddNode( p.X(), p.Y(), p.Z() );
1289 myLastCreatedNodes.Append(newN);
1291 // create a new element
1292 const SMDS_MeshElement* newElem = 0;
1293 const SMDS_MeshNode* N[6];
1301 newElem = aMesh->AddFace(aNodes[2], aNodes[3], aNodes[0],
1302 aNodes[6], aNodes[7], newN );
1311 newElem = aMesh->AddFace(aNodes[3], aNodes[0], aNodes[1],
1312 aNodes[7], aNodes[4], newN );
1314 myLastCreatedElems.Append(newElem);
1315 aMesh->ChangeElementNodes( elem, N, 6 );
1316 // put a new triangle on the same shape and add to the same groups
1318 aMesh->SetMeshElementOnShape( newElem, aShapeId );
1319 AddToSameGroups( newElem, elem, aMesh );
1326 //=======================================================================
1327 //function : getAngle
1329 //=======================================================================
1331 double getAngle(const SMDS_MeshElement * tr1,
1332 const SMDS_MeshElement * tr2,
1333 const SMDS_MeshNode * n1,
1334 const SMDS_MeshNode * n2)
1336 double angle = 2*PI; // bad angle
1339 SMESH::Controls::TSequenceOfXYZ P1, P2;
1340 if ( !SMESH::Controls::NumericalFunctor::GetPoints( tr1, P1 ) ||
1341 !SMESH::Controls::NumericalFunctor::GetPoints( tr2, P2 ))
1344 if(!tr1->IsQuadratic())
1345 N1 = gp_Vec( P1(2) - P1(1) ) ^ gp_Vec( P1(3) - P1(1) );
1347 N1 = gp_Vec( P1(3) - P1(1) ) ^ gp_Vec( P1(5) - P1(1) );
1348 if ( N1.SquareMagnitude() <= gp::Resolution() )
1350 if(!tr2->IsQuadratic())
1351 N2 = gp_Vec( P2(2) - P2(1) ) ^ gp_Vec( P2(3) - P2(1) );
1353 N2 = gp_Vec( P2(3) - P2(1) ) ^ gp_Vec( P2(5) - P2(1) );
1354 if ( N2.SquareMagnitude() <= gp::Resolution() )
1357 // find the first diagonal node n1 in the triangles:
1358 // take in account a diagonal link orientation
1359 const SMDS_MeshElement *nFirst[2], *tr[] = { tr1, tr2 };
1360 for ( int t = 0; t < 2; t++ ) {
1361 SMDS_ElemIteratorPtr it = tr[ t ]->nodesIterator();
1362 int i = 0, iDiag = -1;
1363 while ( it->more()) {
1364 const SMDS_MeshElement *n = it->next();
1365 if ( n == n1 || n == n2 )
1369 if ( i - iDiag == 1 )
1370 nFirst[ t ] = ( n == n1 ? n2 : n1 );
1378 if ( nFirst[ 0 ] == nFirst[ 1 ] )
1381 angle = N1.Angle( N2 );
1386 // =================================================
1387 // class generating a unique ID for a pair of nodes
1388 // and able to return nodes by that ID
1389 // =================================================
1393 LinkID_Gen( const SMESHDS_Mesh* theMesh )
1394 :myMesh( theMesh ), myMaxID( theMesh->MaxNodeID() + 1)
1397 long GetLinkID (const SMDS_MeshNode * n1,
1398 const SMDS_MeshNode * n2) const
1400 return ( Min(n1->GetID(),n2->GetID()) * myMaxID + Max(n1->GetID(),n2->GetID()));
1403 bool GetNodes (const long theLinkID,
1404 const SMDS_MeshNode* & theNode1,
1405 const SMDS_MeshNode* & theNode2) const
1407 theNode1 = myMesh->FindNode( theLinkID / myMaxID );
1408 if ( !theNode1 ) return false;
1409 theNode2 = myMesh->FindNode( theLinkID % myMaxID );
1410 if ( !theNode2 ) return false;
1416 const SMESHDS_Mesh* myMesh;
1421 //=======================================================================
1422 //function : TriToQuad
1423 //purpose : Fuse neighbour triangles into quadrangles.
1424 // theCrit is used to select a neighbour to fuse with.
1425 // theMaxAngle is a max angle between element normals at which
1426 // fusion is still performed.
1427 //=======================================================================
1429 bool SMESH_MeshEditor::TriToQuad (TIDSortedElemSet & theElems,
1430 SMESH::Controls::NumericalFunctorPtr theCrit,
1431 const double theMaxAngle)
1433 myLastCreatedElems.Clear();
1434 myLastCreatedNodes.Clear();
1436 MESSAGE( "::TriToQuad()" );
1438 if ( !theCrit.get() )
1441 SMESHDS_Mesh * aMesh = GetMeshDS();
1443 // Prepare data for algo: build
1444 // 1. map of elements with their linkIDs
1445 // 2. map of linkIDs with their elements
1447 map< SMESH_TLink, list< const SMDS_MeshElement* > > mapLi_listEl;
1448 map< SMESH_TLink, list< const SMDS_MeshElement* > >::iterator itLE;
1449 map< const SMDS_MeshElement*, set< SMESH_TLink > > mapEl_setLi;
1450 map< const SMDS_MeshElement*, set< SMESH_TLink > >::iterator itEL;
1452 TIDSortedElemSet::iterator itElem;
1453 for ( itElem = theElems.begin(); itElem != theElems.end(); itElem++ ) {
1454 const SMDS_MeshElement* elem = *itElem;
1455 if(!elem || elem->GetType() != SMDSAbs_Face ) continue;
1456 bool IsTria = elem->NbNodes()==3 || (elem->NbNodes()==6 && elem->IsQuadratic());
1457 if(!IsTria) continue;
1459 // retrieve element nodes
1460 const SMDS_MeshNode* aNodes [4];
1461 SMDS_ElemIteratorPtr itN = elem->nodesIterator();
1464 aNodes[ i++ ] = cast2Node( itN->next() );
1465 aNodes[ 3 ] = aNodes[ 0 ];
1468 for ( i = 0; i < 3; i++ ) {
1469 SMESH_TLink link( aNodes[i], aNodes[i+1] );
1470 // check if elements sharing a link can be fused
1471 itLE = mapLi_listEl.find( link );
1472 if ( itLE != mapLi_listEl.end() ) {
1473 if ((*itLE).second.size() > 1 ) // consider only 2 elems adjacent by a link
1475 const SMDS_MeshElement* elem2 = (*itLE).second.front();
1476 //if ( FindShape( elem ) != FindShape( elem2 ))
1477 // continue; // do not fuse triangles laying on different shapes
1478 if ( getAngle( elem, elem2, aNodes[i], aNodes[i+1] ) > theMaxAngle )
1479 continue; // avoid making badly shaped quads
1480 (*itLE).second.push_back( elem );
1483 mapLi_listEl[ link ].push_back( elem );
1485 mapEl_setLi [ elem ].insert( link );
1488 // Clean the maps from the links shared by a sole element, ie
1489 // links to which only one element is bound in mapLi_listEl
1491 for ( itLE = mapLi_listEl.begin(); itLE != mapLi_listEl.end(); itLE++ ) {
1492 int nbElems = (*itLE).second.size();
1493 if ( nbElems < 2 ) {
1494 const SMDS_MeshElement* elem = (*itLE).second.front();
1495 SMESH_TLink link = (*itLE).first;
1496 mapEl_setLi[ elem ].erase( link );
1497 if ( mapEl_setLi[ elem ].empty() )
1498 mapEl_setLi.erase( elem );
1502 // Algo: fuse triangles into quadrangles
1504 while ( ! mapEl_setLi.empty() ) {
1505 // Look for the start element:
1506 // the element having the least nb of shared links
1507 const SMDS_MeshElement* startElem = 0;
1509 for ( itEL = mapEl_setLi.begin(); itEL != mapEl_setLi.end(); itEL++ ) {
1510 int nbLinks = (*itEL).second.size();
1511 if ( nbLinks < minNbLinks ) {
1512 startElem = (*itEL).first;
1513 minNbLinks = nbLinks;
1514 if ( minNbLinks == 1 )
1519 // search elements to fuse starting from startElem or links of elements
1520 // fused earlyer - startLinks
1521 list< SMESH_TLink > startLinks;
1522 while ( startElem || !startLinks.empty() ) {
1523 while ( !startElem && !startLinks.empty() ) {
1524 // Get an element to start, by a link
1525 SMESH_TLink linkId = startLinks.front();
1526 startLinks.pop_front();
1527 itLE = mapLi_listEl.find( linkId );
1528 if ( itLE != mapLi_listEl.end() ) {
1529 list< const SMDS_MeshElement* > & listElem = (*itLE).second;
1530 list< const SMDS_MeshElement* >::iterator itE = listElem.begin();
1531 for ( ; itE != listElem.end() ; itE++ )
1532 if ( mapEl_setLi.find( (*itE) ) != mapEl_setLi.end() )
1534 mapLi_listEl.erase( itLE );
1539 // Get candidates to be fused
1540 const SMDS_MeshElement *tr1 = startElem, *tr2 = 0, *tr3 = 0;
1541 const SMESH_TLink *link12, *link13;
1543 ASSERT( mapEl_setLi.find( tr1 ) != mapEl_setLi.end() );
1544 set< SMESH_TLink >& setLi = mapEl_setLi[ tr1 ];
1545 ASSERT( !setLi.empty() );
1546 set< SMESH_TLink >::iterator itLi;
1547 for ( itLi = setLi.begin(); itLi != setLi.end(); itLi++ )
1549 const SMESH_TLink & link = (*itLi);
1550 itLE = mapLi_listEl.find( link );
1551 if ( itLE == mapLi_listEl.end() )
1554 const SMDS_MeshElement* elem = (*itLE).second.front();
1556 elem = (*itLE).second.back();
1557 mapLi_listEl.erase( itLE );
1558 if ( mapEl_setLi.find( elem ) == mapEl_setLi.end())
1569 // add other links of elem to list of links to re-start from
1570 set< SMESH_TLink >& links = mapEl_setLi[ elem ];
1571 set< SMESH_TLink >::iterator it;
1572 for ( it = links.begin(); it != links.end(); it++ ) {
1573 const SMESH_TLink& link2 = (*it);
1574 if ( link2 != link )
1575 startLinks.push_back( link2 );
1579 // Get nodes of possible quadrangles
1580 const SMDS_MeshNode *n12 [4], *n13 [4];
1581 bool Ok12 = false, Ok13 = false;
1582 const SMDS_MeshNode *linkNode1, *linkNode2;
1584 linkNode1 = link12->first;
1585 linkNode2 = link12->second;
1586 if ( tr2 && getQuadrangleNodes( n12, linkNode1, linkNode2, tr1, tr2 ))
1590 linkNode1 = link13->first;
1591 linkNode2 = link13->second;
1592 if ( tr3 && getQuadrangleNodes( n13, linkNode1, linkNode2, tr1, tr3 ))
1596 // Choose a pair to fuse
1597 if ( Ok12 && Ok13 ) {
1598 SMDS_FaceOfNodes quad12 ( n12[ 0 ], n12[ 1 ], n12[ 2 ], n12[ 3 ] );
1599 SMDS_FaceOfNodes quad13 ( n13[ 0 ], n13[ 1 ], n13[ 2 ], n13[ 3 ] );
1600 double aBadRate12 = getBadRate( &quad12, theCrit );
1601 double aBadRate13 = getBadRate( &quad13, theCrit );
1602 if ( aBadRate13 < aBadRate12 )
1609 // and remove fused elems and removed links from the maps
1610 mapEl_setLi.erase( tr1 );
1612 mapEl_setLi.erase( tr2 );
1613 mapLi_listEl.erase( *link12 );
1614 if(tr1->NbNodes()==3) {
1615 if( tr1->GetID() < tr2->GetID() ) {
1616 aMesh->ChangeElementNodes( tr1, n12, 4 );
1617 myLastCreatedElems.Append(tr1);
1618 aMesh->RemoveElement( tr2 );
1621 aMesh->ChangeElementNodes( tr2, n12, 4 );
1622 myLastCreatedElems.Append(tr2);
1623 aMesh->RemoveElement( tr1);
1627 const SMDS_MeshNode* N1 [6];
1628 const SMDS_MeshNode* N2 [6];
1629 GetNodesFromTwoTria(tr1,tr2,N1,N2);
1630 // now we receive following N1 and N2 (using numeration as above image)
1631 // tria1 : (1 2 4 5 9 7) and tria2 : (3 4 2 8 9 6)
1632 // i.e. first nodes from both arrays determ new diagonal
1633 const SMDS_MeshNode* aNodes[8];
1642 if( tr1->GetID() < tr2->GetID() ) {
1643 GetMeshDS()->ChangeElementNodes( tr1, aNodes, 8 );
1644 myLastCreatedElems.Append(tr1);
1645 GetMeshDS()->RemoveElement( tr2 );
1648 GetMeshDS()->ChangeElementNodes( tr2, aNodes, 8 );
1649 myLastCreatedElems.Append(tr2);
1650 GetMeshDS()->RemoveElement( tr1 );
1652 // remove middle node (9)
1653 GetMeshDS()->RemoveNode( N1[4] );
1657 mapEl_setLi.erase( tr3 );
1658 mapLi_listEl.erase( *link13 );
1659 if(tr1->NbNodes()==3) {
1660 if( tr1->GetID() < tr2->GetID() ) {
1661 aMesh->ChangeElementNodes( tr1, n13, 4 );
1662 myLastCreatedElems.Append(tr1);
1663 aMesh->RemoveElement( tr3 );
1666 aMesh->ChangeElementNodes( tr3, n13, 4 );
1667 myLastCreatedElems.Append(tr3);
1668 aMesh->RemoveElement( tr1 );
1672 const SMDS_MeshNode* N1 [6];
1673 const SMDS_MeshNode* N2 [6];
1674 GetNodesFromTwoTria(tr1,tr3,N1,N2);
1675 // now we receive following N1 and N2 (using numeration as above image)
1676 // tria1 : (1 2 4 5 9 7) and tria2 : (3 4 2 8 9 6)
1677 // i.e. first nodes from both arrays determ new diagonal
1678 const SMDS_MeshNode* aNodes[8];
1687 if( tr1->GetID() < tr2->GetID() ) {
1688 GetMeshDS()->ChangeElementNodes( tr1, aNodes, 8 );
1689 myLastCreatedElems.Append(tr1);
1690 GetMeshDS()->RemoveElement( tr3 );
1693 GetMeshDS()->ChangeElementNodes( tr3, aNodes, 8 );
1694 myLastCreatedElems.Append(tr3);
1695 GetMeshDS()->RemoveElement( tr1 );
1697 // remove middle node (9)
1698 GetMeshDS()->RemoveNode( N1[4] );
1702 // Next element to fuse: the rejected one
1704 startElem = Ok12 ? tr3 : tr2;
1706 } // if ( startElem )
1707 } // while ( startElem || !startLinks.empty() )
1708 } // while ( ! mapEl_setLi.empty() )
1714 /*#define DUMPSO(txt) \
1715 // cout << txt << endl;
1716 //=============================================================================
1720 //=============================================================================
1721 static void swap( int i1, int i2, int idNodes[], gp_Pnt P[] )
1725 int tmp = idNodes[ i1 ];
1726 idNodes[ i1 ] = idNodes[ i2 ];
1727 idNodes[ i2 ] = tmp;
1728 gp_Pnt Ptmp = P[ i1 ];
1731 DUMPSO( i1 << "(" << idNodes[ i2 ] << ") <-> " << i2 << "(" << idNodes[ i1 ] << ")");
1734 //=======================================================================
1735 //function : SortQuadNodes
1736 //purpose : Set 4 nodes of a quadrangle face in a good order.
1737 // Swap 1<->2 or 2<->3 nodes and correspondingly return
1739 //=======================================================================
1741 int SMESH_MeshEditor::SortQuadNodes (const SMDS_Mesh * theMesh,
1746 for ( i = 0; i < 4; i++ ) {
1747 const SMDS_MeshNode *n = theMesh->FindNode( idNodes[i] );
1749 P[ i ].SetCoord( n->X(), n->Y(), n->Z() );
1752 gp_Vec V1(P[0], P[1]);
1753 gp_Vec V2(P[0], P[2]);
1754 gp_Vec V3(P[0], P[3]);
1756 gp_Vec Cross1 = V1 ^ V2;
1757 gp_Vec Cross2 = V2 ^ V3;
1760 if (Cross1.Dot(Cross2) < 0)
1765 if (Cross1.Dot(Cross2) < 0)
1769 swap ( i, i + 1, idNodes, P );
1771 // for ( int ii = 0; ii < 4; ii++ ) {
1772 // const SMDS_MeshNode *n = theMesh->FindNode( idNodes[ii] );
1773 // DUMPSO( ii << "(" << idNodes[ii] <<") : "<<n->X()<<" "<<n->Y()<<" "<<n->Z());
1779 //=======================================================================
1780 //function : SortHexaNodes
1781 //purpose : Set 8 nodes of a hexahedron in a good order.
1782 // Return success status
1783 //=======================================================================
1785 bool SMESH_MeshEditor::SortHexaNodes (const SMDS_Mesh * theMesh,
1790 DUMPSO( "INPUT: ========================================");
1791 for ( i = 0; i < 8; i++ ) {
1792 const SMDS_MeshNode *n = theMesh->FindNode( idNodes[i] );
1793 if ( !n ) return false;
1794 P[ i ].SetCoord( n->X(), n->Y(), n->Z() );
1795 DUMPSO( i << "(" << idNodes[i] <<") : "<<n->X()<<" "<<n->Y()<<" "<<n->Z());
1797 DUMPSO( "========================================");
1800 set<int> faceNodes; // ids of bottom face nodes, to be found
1801 set<int> checkedId1; // ids of tried 2-nd nodes
1802 Standard_Real leastDist = DBL_MAX; // dist of the 4-th node from 123 plane
1803 const Standard_Real tol = 1.e-6; // tolerance to find nodes in plane
1804 int iMin, iLoop1 = 0;
1806 // Loop to try the 2-nd nodes
1808 while ( leastDist > DBL_MIN && ++iLoop1 < 8 )
1810 // Find not checked 2-nd node
1811 for ( i = 1; i < 8; i++ )
1812 if ( checkedId1.find( idNodes[i] ) == checkedId1.end() ) {
1813 int id1 = idNodes[i];
1814 swap ( 1, i, idNodes, P );
1815 checkedId1.insert ( id1 );
1819 // Find the 3-d node so that 1-2-3 triangle to be on a hexa face,
1820 // ie that all but meybe one (id3 which is on the same face) nodes
1821 // lay on the same side from the triangle plane.
1823 bool manyInPlane = false; // more than 4 nodes lay in plane
1825 while ( ++iLoop2 < 6 ) {
1827 // get 1-2-3 plane coeffs
1828 Standard_Real A, B, C, D;
1829 gp_Vec N = gp_Vec (P[0], P[1]).Crossed( gp_Vec (P[0], P[2]) );
1830 if ( N.SquareMagnitude() > gp::Resolution() )
1832 gp_Pln pln ( P[0], N );
1833 pln.Coefficients( A, B, C, D );
1835 // find the node (iMin) closest to pln
1836 Standard_Real dist[ 8 ], minDist = DBL_MAX;
1838 for ( i = 3; i < 8; i++ ) {
1839 dist[i] = A * P[i].X() + B * P[i].Y() + C * P[i].Z() + D;
1840 if ( fabs( dist[i] ) < minDist ) {
1841 minDist = fabs( dist[i] );
1844 if ( fabs( dist[i] ) <= tol )
1845 idInPln.insert( idNodes[i] );
1848 // there should not be more than 4 nodes in bottom plane
1849 if ( idInPln.size() > 1 )
1851 DUMPSO( "### idInPln.size() = " << idInPln.size());
1852 // idInPlane does not contain the first 3 nodes
1853 if ( manyInPlane || idInPln.size() == 5)
1854 return false; // all nodes in one plane
1857 // set the 1-st node to be not in plane
1858 for ( i = 3; i < 8; i++ ) {
1859 if ( idInPln.find( idNodes[ i ] ) == idInPln.end() ) {
1860 DUMPSO( "### Reset 0-th node");
1861 swap( 0, i, idNodes, P );
1866 // reset to re-check second nodes
1867 leastDist = DBL_MAX;
1871 break; // from iLoop2;
1874 // check that the other 4 nodes are on the same side
1875 bool sameSide = true;
1876 bool isNeg = dist[ iMin == 3 ? 4 : 3 ] <= 0.;
1877 for ( i = 3; sameSide && i < 8; i++ ) {
1879 sameSide = ( isNeg == dist[i] <= 0.);
1882 // keep best solution
1883 if ( sameSide && minDist < leastDist ) {
1884 leastDist = minDist;
1886 faceNodes.insert( idNodes[ 1 ] );
1887 faceNodes.insert( idNodes[ 2 ] );
1888 faceNodes.insert( idNodes[ iMin ] );
1889 DUMPSO( "loop " << iLoop2 << " id2 " << idNodes[ 1 ] << " id3 " << idNodes[ 2 ]
1890 << " leastDist = " << leastDist);
1891 if ( leastDist <= DBL_MIN )
1896 // set next 3-d node to check
1897 int iNext = 2 + iLoop2;
1899 DUMPSO( "Try 2-nd");
1900 swap ( 2, iNext, idNodes, P );
1902 } // while ( iLoop2 < 6 )
1905 if ( faceNodes.empty() ) return false;
1907 // Put the faceNodes in proper places
1908 for ( i = 4; i < 8; i++ ) {
1909 if ( faceNodes.find( idNodes[ i ] ) != faceNodes.end() ) {
1910 // find a place to put
1912 while ( faceNodes.find( idNodes[ iTo ] ) != faceNodes.end() )
1914 DUMPSO( "Set faceNodes");
1915 swap ( iTo, i, idNodes, P );
1920 // Set nodes of the found bottom face in good order
1921 DUMPSO( " Found bottom face: ");
1922 i = SortQuadNodes( theMesh, idNodes );
1924 gp_Pnt Ptmp = P[ i ];
1929 // for ( int ii = 0; ii < 4; ii++ ) {
1930 // const SMDS_MeshNode *n = theMesh->FindNode( idNodes[ii] );
1931 // DUMPSO( ii << "(" << idNodes[ii] <<") : "<<n->X()<<" "<<n->Y()<<" "<<n->Z());
1934 // Gravity center of the top and bottom faces
1935 gp_Pnt aGCb = ( P[0].XYZ() + P[1].XYZ() + P[2].XYZ() + P[3].XYZ() ) / 4.;
1936 gp_Pnt aGCt = ( P[4].XYZ() + P[5].XYZ() + P[6].XYZ() + P[7].XYZ() ) / 4.;
1938 // Get direction from the bottom to the top face
1939 gp_Vec upDir ( aGCb, aGCt );
1940 Standard_Real upDirSize = upDir.Magnitude();
1941 if ( upDirSize <= gp::Resolution() ) return false;
1944 // Assure that the bottom face normal points up
1945 gp_Vec Nb = gp_Vec (P[0], P[1]).Crossed( gp_Vec (P[0], P[2]) );
1946 Nb += gp_Vec (P[0], P[2]).Crossed( gp_Vec (P[0], P[3]) );
1947 if ( Nb.Dot( upDir ) < 0 ) {
1948 DUMPSO( "Reverse bottom face");
1949 swap( 1, 3, idNodes, P );
1952 // Find 5-th node - the one closest to the 1-st among the last 4 nodes.
1953 Standard_Real minDist = DBL_MAX;
1954 for ( i = 4; i < 8; i++ ) {
1955 // projection of P[i] to the plane defined by P[0] and upDir
1956 gp_Pnt Pp = P[i].Translated( upDir * ( upDir.Dot( gp_Vec( P[i], P[0] ))));
1957 Standard_Real sqDist = P[0].SquareDistance( Pp );
1958 if ( sqDist < minDist ) {
1963 DUMPSO( "Set 4-th");
1964 swap ( 4, iMin, idNodes, P );
1966 // Set nodes of the top face in good order
1967 DUMPSO( "Sort top face");
1968 i = SortQuadNodes( theMesh, &idNodes[4] );
1971 gp_Pnt Ptmp = P[ i ];
1976 // Assure that direction of the top face normal is from the bottom face
1977 gp_Vec Nt = gp_Vec (P[4], P[5]).Crossed( gp_Vec (P[4], P[6]) );
1978 Nt += gp_Vec (P[4], P[6]).Crossed( gp_Vec (P[4], P[7]) );
1979 if ( Nt.Dot( upDir ) < 0 ) {
1980 DUMPSO( "Reverse top face");
1981 swap( 5, 7, idNodes, P );
1984 // DUMPSO( "OUTPUT: ========================================");
1985 // for ( i = 0; i < 8; i++ ) {
1986 // float *p = ugrid->GetPoint(idNodes[i]);
1987 // DUMPSO( i << "(" << idNodes[i] << ") : " << p[0] << " " << p[1] << " " << p[2]);
1993 //================================================================================
1995 * \brief Return nodes linked to the given one
1996 * \param theNode - the node
1997 * \param linkedNodes - the found nodes
1998 * \param type - the type of elements to check
2000 * Medium nodes are ignored
2002 //================================================================================
2004 void SMESH_MeshEditor::GetLinkedNodes( const SMDS_MeshNode* theNode,
2005 TIDSortedElemSet & linkedNodes,
2006 SMDSAbs_ElementType type )
2008 SMDS_ElemIteratorPtr elemIt = theNode->GetInverseElementIterator(type);
2009 while ( elemIt->more() )
2011 const SMDS_MeshElement* elem = elemIt->next();
2012 SMDS_ElemIteratorPtr nodeIt = elem->nodesIterator();
2013 if ( elem->GetType() == SMDSAbs_Volume )
2015 SMDS_VolumeTool vol( elem );
2016 while ( nodeIt->more() ) {
2017 const SMDS_MeshNode* n = cast2Node( nodeIt->next() );
2018 if ( theNode != n && vol.IsLinked( theNode, n ))
2019 linkedNodes.insert( n );
2024 for ( int i = 0; nodeIt->more(); ++i ) {
2025 const SMDS_MeshNode* n = cast2Node( nodeIt->next() );
2026 if ( n == theNode ) {
2027 int iBefore = i - 1;
2029 if ( elem->IsQuadratic() ) {
2030 int nb = elem->NbNodes() / 2;
2031 iAfter = SMESH_MesherHelper::WrapIndex( iAfter, nb );
2032 iBefore = SMESH_MesherHelper::WrapIndex( iBefore, nb );
2034 linkedNodes.insert( elem->GetNode( iAfter ));
2035 linkedNodes.insert( elem->GetNode( iBefore ));
2042 //=======================================================================
2043 //function : laplacianSmooth
2044 //purpose : pulls theNode toward the center of surrounding nodes directly
2045 // connected to that node along an element edge
2046 //=======================================================================
2048 void laplacianSmooth(const SMDS_MeshNode* theNode,
2049 const Handle(Geom_Surface)& theSurface,
2050 map< const SMDS_MeshNode*, gp_XY* >& theUVMap)
2052 // find surrounding nodes
2054 TIDSortedElemSet nodeSet;
2055 SMESH_MeshEditor::GetLinkedNodes( theNode, nodeSet, SMDSAbs_Face );
2057 // compute new coodrs
2059 double coord[] = { 0., 0., 0. };
2060 TIDSortedElemSet::iterator nodeSetIt = nodeSet.begin();
2061 for ( ; nodeSetIt != nodeSet.end(); nodeSetIt++ ) {
2062 const SMDS_MeshNode* node = cast2Node(*nodeSetIt);
2063 if ( theSurface.IsNull() ) { // smooth in 3D
2064 coord[0] += node->X();
2065 coord[1] += node->Y();
2066 coord[2] += node->Z();
2068 else { // smooth in 2D
2069 ASSERT( theUVMap.find( node ) != theUVMap.end() );
2070 gp_XY* uv = theUVMap[ node ];
2071 coord[0] += uv->X();
2072 coord[1] += uv->Y();
2075 int nbNodes = nodeSet.size();
2078 coord[0] /= nbNodes;
2079 coord[1] /= nbNodes;
2081 if ( !theSurface.IsNull() ) {
2082 ASSERT( theUVMap.find( theNode ) != theUVMap.end() );
2083 theUVMap[ theNode ]->SetCoord( coord[0], coord[1] );
2084 gp_Pnt p3d = theSurface->Value( coord[0], coord[1] );
2090 coord[2] /= nbNodes;
2094 const_cast< SMDS_MeshNode* >( theNode )->setXYZ(coord[0],coord[1],coord[2]);
2097 //=======================================================================
2098 //function : centroidalSmooth
2099 //purpose : pulls theNode toward the element-area-weighted centroid of the
2100 // surrounding elements
2101 //=======================================================================
2103 void centroidalSmooth(const SMDS_MeshNode* theNode,
2104 const Handle(Geom_Surface)& theSurface,
2105 map< const SMDS_MeshNode*, gp_XY* >& theUVMap)
2107 gp_XYZ aNewXYZ(0.,0.,0.);
2108 SMESH::Controls::Area anAreaFunc;
2109 double totalArea = 0.;
2114 SMDS_ElemIteratorPtr elemIt = theNode->GetInverseElementIterator(SMDSAbs_Face);
2115 while ( elemIt->more() )
2117 const SMDS_MeshElement* elem = elemIt->next();
2120 gp_XYZ elemCenter(0.,0.,0.);
2121 SMESH::Controls::TSequenceOfXYZ aNodePoints;
2122 SMDS_ElemIteratorPtr itN = elem->nodesIterator();
2123 int nn = elem->NbNodes();
2124 if(elem->IsQuadratic()) nn = nn/2;
2126 //while ( itN->more() ) {
2128 const SMDS_MeshNode* aNode = static_cast<const SMDS_MeshNode*>( itN->next() );
2130 gp_XYZ aP( aNode->X(), aNode->Y(), aNode->Z() );
2131 aNodePoints.push_back( aP );
2132 if ( !theSurface.IsNull() ) { // smooth in 2D
2133 ASSERT( theUVMap.find( aNode ) != theUVMap.end() );
2134 gp_XY* uv = theUVMap[ aNode ];
2135 aP.SetCoord( uv->X(), uv->Y(), 0. );
2139 double elemArea = anAreaFunc.GetValue( aNodePoints );
2140 totalArea += elemArea;
2142 aNewXYZ += elemCenter * elemArea;
2144 aNewXYZ /= totalArea;
2145 if ( !theSurface.IsNull() ) {
2146 theUVMap[ theNode ]->SetCoord( aNewXYZ.X(), aNewXYZ.Y() );
2147 aNewXYZ = theSurface->Value( aNewXYZ.X(), aNewXYZ.Y() ).XYZ();
2152 const_cast< SMDS_MeshNode* >( theNode )->setXYZ(aNewXYZ.X(),aNewXYZ.Y(),aNewXYZ.Z());
2155 //=======================================================================
2156 //function : getClosestUV
2157 //purpose : return UV of closest projection
2158 //=======================================================================
2160 static bool getClosestUV (Extrema_GenExtPS& projector,
2161 const gp_Pnt& point,
2164 projector.Perform( point );
2165 if ( projector.IsDone() ) {
2166 double u, v, minVal = DBL_MAX;
2167 for ( int i = projector.NbExt(); i > 0; i-- )
2168 if ( projector.Value( i ) < minVal ) {
2169 minVal = projector.Value( i );
2170 projector.Point( i ).Parameter( u, v );
2172 result.SetCoord( u, v );
2178 //=======================================================================
2180 //purpose : Smooth theElements during theNbIterations or until a worst
2181 // element has aspect ratio <= theTgtAspectRatio.
2182 // Aspect Ratio varies in range [1.0, inf].
2183 // If theElements is empty, the whole mesh is smoothed.
2184 // theFixedNodes contains additionally fixed nodes. Nodes built
2185 // on edges and boundary nodes are always fixed.
2186 //=======================================================================
2188 void SMESH_MeshEditor::Smooth (TIDSortedElemSet & theElems,
2189 set<const SMDS_MeshNode*> & theFixedNodes,
2190 const SmoothMethod theSmoothMethod,
2191 const int theNbIterations,
2192 double theTgtAspectRatio,
2195 myLastCreatedElems.Clear();
2196 myLastCreatedNodes.Clear();
2198 MESSAGE((theSmoothMethod==LAPLACIAN ? "LAPLACIAN" : "CENTROIDAL") << "--::Smooth()");
2200 if ( theTgtAspectRatio < 1.0 )
2201 theTgtAspectRatio = 1.0;
2203 const double disttol = 1.e-16;
2205 SMESH::Controls::AspectRatio aQualityFunc;
2207 SMESHDS_Mesh* aMesh = GetMeshDS();
2209 if ( theElems.empty() ) {
2210 // add all faces to theElems
2211 SMDS_FaceIteratorPtr fIt = aMesh->facesIterator();
2212 while ( fIt->more() ) {
2213 const SMDS_MeshElement* face = fIt->next();
2214 theElems.insert( face );
2217 // get all face ids theElems are on
2218 set< int > faceIdSet;
2219 TIDSortedElemSet::iterator itElem;
2221 for ( itElem = theElems.begin(); itElem != theElems.end(); itElem++ ) {
2222 int fId = FindShape( *itElem );
2223 // check that corresponding submesh exists and a shape is face
2225 faceIdSet.find( fId ) == faceIdSet.end() &&
2226 aMesh->MeshElements( fId )) {
2227 TopoDS_Shape F = aMesh->IndexToShape( fId );
2228 if ( !F.IsNull() && F.ShapeType() == TopAbs_FACE )
2229 faceIdSet.insert( fId );
2232 faceIdSet.insert( 0 ); // to smooth elements that are not on any TopoDS_Face
2234 // ===============================================
2235 // smooth elements on each TopoDS_Face separately
2236 // ===============================================
2238 set< int >::reverse_iterator fId = faceIdSet.rbegin(); // treate 0 fId at the end
2239 for ( ; fId != faceIdSet.rend(); ++fId ) {
2240 // get face surface and submesh
2241 Handle(Geom_Surface) surface;
2242 SMESHDS_SubMesh* faceSubMesh = 0;
2244 double fToler2 = 0, vPeriod = 0., uPeriod = 0., f,l;
2245 double u1 = 0, u2 = 0, v1 = 0, v2 = 0;
2246 bool isUPeriodic = false, isVPeriodic = false;
2248 face = TopoDS::Face( aMesh->IndexToShape( *fId ));
2249 surface = BRep_Tool::Surface( face );
2250 faceSubMesh = aMesh->MeshElements( *fId );
2251 fToler2 = BRep_Tool::Tolerance( face );
2252 fToler2 *= fToler2 * 10.;
2253 isUPeriodic = surface->IsUPeriodic();
2255 vPeriod = surface->UPeriod();
2256 isVPeriodic = surface->IsVPeriodic();
2258 uPeriod = surface->VPeriod();
2259 surface->Bounds( u1, u2, v1, v2 );
2261 // ---------------------------------------------------------
2262 // for elements on a face, find movable and fixed nodes and
2263 // compute UV for them
2264 // ---------------------------------------------------------
2265 bool checkBoundaryNodes = false;
2266 bool isQuadratic = false;
2267 set<const SMDS_MeshNode*> setMovableNodes;
2268 map< const SMDS_MeshNode*, gp_XY* > uvMap, uvMap2;
2269 list< gp_XY > listUV; // uvs the 2 uvMaps refer to
2270 list< const SMDS_MeshElement* > elemsOnFace;
2272 Extrema_GenExtPS projector;
2273 GeomAdaptor_Surface surfAdaptor;
2274 if ( !surface.IsNull() ) {
2275 surfAdaptor.Load( surface );
2276 projector.Initialize( surfAdaptor, 20,20, 1e-5,1e-5 );
2278 int nbElemOnFace = 0;
2279 itElem = theElems.begin();
2280 // loop on not yet smoothed elements: look for elems on a face
2281 while ( itElem != theElems.end() ) {
2282 if ( faceSubMesh && nbElemOnFace == faceSubMesh->NbElements() )
2283 break; // all elements found
2285 const SMDS_MeshElement* elem = *itElem;
2286 if ( !elem || elem->GetType() != SMDSAbs_Face || elem->NbNodes() < 3 ||
2287 ( faceSubMesh && !faceSubMesh->Contains( elem ))) {
2291 elemsOnFace.push_back( elem );
2292 theElems.erase( itElem++ );
2296 isQuadratic = elem->IsQuadratic();
2298 // get movable nodes of elem
2299 const SMDS_MeshNode* node;
2300 SMDS_TypeOfPosition posType;
2301 SMDS_ElemIteratorPtr itN = elem->nodesIterator();
2302 int nn = 0, nbn = elem->NbNodes();
2303 if(elem->IsQuadratic())
2305 while ( nn++ < nbn ) {
2306 node = static_cast<const SMDS_MeshNode*>( itN->next() );
2307 const SMDS_PositionPtr& pos = node->GetPosition();
2308 posType = pos.get() ? pos->GetTypeOfPosition() : SMDS_TOP_3DSPACE;
2309 if (posType != SMDS_TOP_EDGE &&
2310 posType != SMDS_TOP_VERTEX &&
2311 theFixedNodes.find( node ) == theFixedNodes.end())
2313 // check if all faces around the node are on faceSubMesh
2314 // because a node on edge may be bound to face
2315 SMDS_ElemIteratorPtr eIt = node->GetInverseElementIterator(SMDSAbs_Face);
2317 if ( faceSubMesh ) {
2318 while ( eIt->more() && all ) {
2319 const SMDS_MeshElement* e = eIt->next();
2320 all = faceSubMesh->Contains( e );
2324 setMovableNodes.insert( node );
2326 checkBoundaryNodes = true;
2328 if ( posType == SMDS_TOP_3DSPACE )
2329 checkBoundaryNodes = true;
2332 if ( surface.IsNull() )
2335 // get nodes to check UV
2336 list< const SMDS_MeshNode* > uvCheckNodes;
2337 itN = elem->nodesIterator();
2338 nn = 0; nbn = elem->NbNodes();
2339 if(elem->IsQuadratic())
2341 while ( nn++ < nbn ) {
2342 node = static_cast<const SMDS_MeshNode*>( itN->next() );
2343 if ( uvMap.find( node ) == uvMap.end() )
2344 uvCheckNodes.push_back( node );
2345 // add nodes of elems sharing node
2346 // SMDS_ElemIteratorPtr eIt = node->GetInverseElementIterator(SMDSAbs_Face);
2347 // while ( eIt->more() ) {
2348 // const SMDS_MeshElement* e = eIt->next();
2349 // if ( e != elem ) {
2350 // SMDS_ElemIteratorPtr nIt = e->nodesIterator();
2351 // while ( nIt->more() ) {
2352 // const SMDS_MeshNode* n =
2353 // static_cast<const SMDS_MeshNode*>( nIt->next() );
2354 // if ( uvMap.find( n ) == uvMap.end() )
2355 // uvCheckNodes.push_back( n );
2361 list< const SMDS_MeshNode* >::iterator n = uvCheckNodes.begin();
2362 for ( ; n != uvCheckNodes.end(); ++n ) {
2365 const SMDS_PositionPtr& pos = node->GetPosition();
2366 posType = pos.get() ? pos->GetTypeOfPosition() : SMDS_TOP_3DSPACE;
2368 switch ( posType ) {
2369 case SMDS_TOP_FACE: {
2370 SMDS_FacePosition* fPos = ( SMDS_FacePosition* ) pos.get();
2371 uv.SetCoord( fPos->GetUParameter(), fPos->GetVParameter() );
2374 case SMDS_TOP_EDGE: {
2375 TopoDS_Shape S = aMesh->IndexToShape( pos->GetShapeId() );
2376 Handle(Geom2d_Curve) pcurve;
2377 if ( !S.IsNull() && S.ShapeType() == TopAbs_EDGE )
2378 pcurve = BRep_Tool::CurveOnSurface( TopoDS::Edge( S ), face, f,l );
2379 if ( !pcurve.IsNull() ) {
2380 double u = (( SMDS_EdgePosition* ) pos.get() )->GetUParameter();
2381 uv = pcurve->Value( u ).XY();
2385 case SMDS_TOP_VERTEX: {
2386 TopoDS_Shape S = aMesh->IndexToShape( pos->GetShapeId() );
2387 if ( !S.IsNull() && S.ShapeType() == TopAbs_VERTEX )
2388 uv = BRep_Tool::Parameters( TopoDS::Vertex( S ), face ).XY();
2393 // check existing UV
2394 bool project = true;
2395 gp_Pnt pNode ( node->X(), node->Y(), node->Z() );
2396 double dist1 = DBL_MAX, dist2 = 0;
2397 if ( posType != SMDS_TOP_3DSPACE ) {
2398 dist1 = pNode.SquareDistance( surface->Value( uv.X(), uv.Y() ));
2399 project = dist1 > fToler2;
2401 if ( project ) { // compute new UV
2403 if ( !getClosestUV( projector, pNode, newUV )) {
2404 MESSAGE("Node Projection Failed " << node);
2408 newUV.SetX( ElCLib::InPeriod( newUV.X(), u1, u2 ));
2410 newUV.SetY( ElCLib::InPeriod( newUV.Y(), v1, v2 ));
2412 if ( posType != SMDS_TOP_3DSPACE )
2413 dist2 = pNode.SquareDistance( surface->Value( newUV.X(), newUV.Y() ));
2414 if ( dist2 < dist1 )
2418 // store UV in the map
2419 listUV.push_back( uv );
2420 uvMap.insert( make_pair( node, &listUV.back() ));
2422 } // loop on not yet smoothed elements
2424 if ( !faceSubMesh || nbElemOnFace != faceSubMesh->NbElements() )
2425 checkBoundaryNodes = true;
2427 // fix nodes on mesh boundary
2429 if ( checkBoundaryNodes ) {
2430 map< NLink, int > linkNbMap; // how many times a link encounters in elemsOnFace
2431 map< NLink, int >::iterator link_nb;
2432 // put all elements links to linkNbMap
2433 list< const SMDS_MeshElement* >::iterator elemIt = elemsOnFace.begin();
2434 for ( ; elemIt != elemsOnFace.end(); ++elemIt ) {
2435 const SMDS_MeshElement* elem = (*elemIt);
2436 int nbn = elem->NbNodes();
2437 if(elem->IsQuadratic())
2439 // loop on elem links: insert them in linkNbMap
2440 const SMDS_MeshNode* curNode, *prevNode = elem->GetNode( nbn );
2441 for ( int iN = 0; iN < nbn; ++iN ) {
2442 curNode = elem->GetNode( iN );
2444 if ( curNode < prevNode ) link = make_pair( curNode , prevNode );
2445 else link = make_pair( prevNode , curNode );
2447 link_nb = linkNbMap.find( link );
2448 if ( link_nb == linkNbMap.end() )
2449 linkNbMap.insert( make_pair ( link, 1 ));
2454 // remove nodes that are in links encountered only once from setMovableNodes
2455 for ( link_nb = linkNbMap.begin(); link_nb != linkNbMap.end(); ++link_nb ) {
2456 if ( link_nb->second == 1 ) {
2457 setMovableNodes.erase( link_nb->first.first );
2458 setMovableNodes.erase( link_nb->first.second );
2463 // -----------------------------------------------------
2464 // for nodes on seam edge, compute one more UV ( uvMap2 );
2465 // find movable nodes linked to nodes on seam and which
2466 // are to be smoothed using the second UV ( uvMap2 )
2467 // -----------------------------------------------------
2469 set<const SMDS_MeshNode*> nodesNearSeam; // to smooth using uvMap2
2470 if ( !surface.IsNull() ) {
2471 TopExp_Explorer eExp( face, TopAbs_EDGE );
2472 for ( ; eExp.More(); eExp.Next() ) {
2473 TopoDS_Edge edge = TopoDS::Edge( eExp.Current() );
2474 if ( !BRep_Tool::IsClosed( edge, face ))
2476 SMESHDS_SubMesh* sm = aMesh->MeshElements( edge );
2477 if ( !sm ) continue;
2478 // find out which parameter varies for a node on seam
2481 Handle(Geom2d_Curve) pcurve = BRep_Tool::CurveOnSurface( edge, face, f, l );
2482 if ( pcurve.IsNull() ) continue;
2483 uv1 = pcurve->Value( f );
2485 pcurve = BRep_Tool::CurveOnSurface( edge, face, f, l );
2486 if ( pcurve.IsNull() ) continue;
2487 uv2 = pcurve->Value( f );
2488 int iPar = Abs( uv1.X() - uv2.X() ) > Abs( uv1.Y() - uv2.Y() ) ? 1 : 2;
2490 if ( uv1.Coord( iPar ) > uv2.Coord( iPar )) {
2491 gp_Pnt2d tmp = uv1; uv1 = uv2; uv2 = tmp;
2493 // get nodes on seam and its vertices
2494 list< const SMDS_MeshNode* > seamNodes;
2495 SMDS_NodeIteratorPtr nSeamIt = sm->GetNodes();
2496 while ( nSeamIt->more() ) {
2497 const SMDS_MeshNode* node = nSeamIt->next();
2498 if ( !isQuadratic || !IsMedium( node ))
2499 seamNodes.push_back( node );
2501 TopExp_Explorer vExp( edge, TopAbs_VERTEX );
2502 for ( ; vExp.More(); vExp.Next() ) {
2503 sm = aMesh->MeshElements( vExp.Current() );
2505 nSeamIt = sm->GetNodes();
2506 while ( nSeamIt->more() )
2507 seamNodes.push_back( nSeamIt->next() );
2510 // loop on nodes on seam
2511 list< const SMDS_MeshNode* >::iterator noSeIt = seamNodes.begin();
2512 for ( ; noSeIt != seamNodes.end(); ++noSeIt ) {
2513 const SMDS_MeshNode* nSeam = *noSeIt;
2514 map< const SMDS_MeshNode*, gp_XY* >::iterator n_uv = uvMap.find( nSeam );
2515 if ( n_uv == uvMap.end() )
2518 n_uv->second->SetCoord( iPar, uv1.Coord( iPar ));
2519 // set the second UV
2520 listUV.push_back( *n_uv->second );
2521 listUV.back().SetCoord( iPar, uv2.Coord( iPar ));
2522 if ( uvMap2.empty() )
2523 uvMap2 = uvMap; // copy the uvMap contents
2524 uvMap2[ nSeam ] = &listUV.back();
2526 // collect movable nodes linked to ones on seam in nodesNearSeam
2527 SMDS_ElemIteratorPtr eIt = nSeam->GetInverseElementIterator(SMDSAbs_Face);
2528 while ( eIt->more() ) {
2529 const SMDS_MeshElement* e = eIt->next();
2530 int nbUseMap1 = 0, nbUseMap2 = 0;
2531 SMDS_ElemIteratorPtr nIt = e->nodesIterator();
2532 int nn = 0, nbn = e->NbNodes();
2533 if(e->IsQuadratic()) nbn = nbn/2;
2534 while ( nn++ < nbn )
2536 const SMDS_MeshNode* n =
2537 static_cast<const SMDS_MeshNode*>( nIt->next() );
2539 setMovableNodes.find( n ) == setMovableNodes.end() )
2541 // add only nodes being closer to uv2 than to uv1
2542 gp_Pnt pMid (0.5 * ( n->X() + nSeam->X() ),
2543 0.5 * ( n->Y() + nSeam->Y() ),
2544 0.5 * ( n->Z() + nSeam->Z() ));
2546 getClosestUV( projector, pMid, uv );
2547 if ( uv.Coord( iPar ) > uvMap[ n ]->Coord( iPar ) ) {
2548 nodesNearSeam.insert( n );
2554 // for centroidalSmooth all element nodes must
2555 // be on one side of a seam
2556 if ( theSmoothMethod == CENTROIDAL && nbUseMap1 && nbUseMap2 ) {
2557 SMDS_ElemIteratorPtr nIt = e->nodesIterator();
2559 while ( nn++ < nbn ) {
2560 const SMDS_MeshNode* n =
2561 static_cast<const SMDS_MeshNode*>( nIt->next() );
2562 setMovableNodes.erase( n );
2566 } // loop on nodes on seam
2567 } // loop on edge of a face
2568 } // if ( !face.IsNull() )
2570 if ( setMovableNodes.empty() ) {
2571 MESSAGE( "Face id : " << *fId << " - NO SMOOTHING: no nodes to move!!!");
2572 continue; // goto next face
2580 double maxRatio = -1., maxDisplacement = -1.;
2581 set<const SMDS_MeshNode*>::iterator nodeToMove;
2582 for ( it = 0; it < theNbIterations; it++ ) {
2583 maxDisplacement = 0.;
2584 nodeToMove = setMovableNodes.begin();
2585 for ( ; nodeToMove != setMovableNodes.end(); nodeToMove++ ) {
2586 const SMDS_MeshNode* node = (*nodeToMove);
2587 gp_XYZ aPrevPos ( node->X(), node->Y(), node->Z() );
2590 bool map2 = ( nodesNearSeam.find( node ) != nodesNearSeam.end() );
2591 if ( theSmoothMethod == LAPLACIAN )
2592 laplacianSmooth( node, surface, map2 ? uvMap2 : uvMap );
2594 centroidalSmooth( node, surface, map2 ? uvMap2 : uvMap );
2596 // node displacement
2597 gp_XYZ aNewPos ( node->X(), node->Y(), node->Z() );
2598 Standard_Real aDispl = (aPrevPos - aNewPos).SquareModulus();
2599 if ( aDispl > maxDisplacement )
2600 maxDisplacement = aDispl;
2602 // no node movement => exit
2603 //if ( maxDisplacement < 1.e-16 ) {
2604 if ( maxDisplacement < disttol ) {
2605 MESSAGE("-- no node movement --");
2609 // check elements quality
2611 list< const SMDS_MeshElement* >::iterator elemIt = elemsOnFace.begin();
2612 for ( ; elemIt != elemsOnFace.end(); ++elemIt ) {
2613 const SMDS_MeshElement* elem = (*elemIt);
2614 if ( !elem || elem->GetType() != SMDSAbs_Face )
2616 SMESH::Controls::TSequenceOfXYZ aPoints;
2617 if ( aQualityFunc.GetPoints( elem, aPoints )) {
2618 double aValue = aQualityFunc.GetValue( aPoints );
2619 if ( aValue > maxRatio )
2623 if ( maxRatio <= theTgtAspectRatio ) {
2624 MESSAGE("-- quality achived --");
2627 if (it+1 == theNbIterations) {
2628 MESSAGE("-- Iteration limit exceeded --");
2630 } // smoothing iterations
2632 MESSAGE(" Face id: " << *fId <<
2633 " Nb iterstions: " << it <<
2634 " Displacement: " << maxDisplacement <<
2635 " Aspect Ratio " << maxRatio);
2637 // ---------------------------------------
2638 // new nodes positions are computed,
2639 // record movement in DS and set new UV
2640 // ---------------------------------------
2641 nodeToMove = setMovableNodes.begin();
2642 for ( ; nodeToMove != setMovableNodes.end(); nodeToMove++ ) {
2643 SMDS_MeshNode* node = const_cast< SMDS_MeshNode* > (*nodeToMove);
2644 aMesh->MoveNode( node, node->X(), node->Y(), node->Z() );
2645 map< const SMDS_MeshNode*, gp_XY* >::iterator node_uv = uvMap.find( node );
2646 if ( node_uv != uvMap.end() ) {
2647 gp_XY* uv = node_uv->second;
2649 ( SMDS_PositionPtr( new SMDS_FacePosition( *fId, uv->X(), uv->Y() )));
2653 // move medium nodes of quadratic elements
2656 SMESH_MesherHelper helper( *GetMesh() );
2657 if ( !face.IsNull() )
2658 helper.SetSubShape( face );
2659 list< const SMDS_MeshElement* >::iterator elemIt = elemsOnFace.begin();
2660 for ( ; elemIt != elemsOnFace.end(); ++elemIt ) {
2661 const SMDS_QuadraticFaceOfNodes* QF =
2662 dynamic_cast<const SMDS_QuadraticFaceOfNodes*> (*elemIt);
2664 vector<const SMDS_MeshNode*> Ns;
2665 Ns.reserve(QF->NbNodes()+1);
2666 SMDS_NodeIteratorPtr anIter = QF->interlacedNodesIterator();
2667 while ( anIter->more() )
2668 Ns.push_back( anIter->next() );
2669 Ns.push_back( Ns[0] );
2671 for(int i=0; i<QF->NbNodes(); i=i+2) {
2672 if ( !surface.IsNull() ) {
2673 gp_XY uv1 = helper.GetNodeUV( face, Ns[i], Ns[i+2] );
2674 gp_XY uv2 = helper.GetNodeUV( face, Ns[i+2], Ns[i] );
2675 gp_XY uv = ( uv1 + uv2 ) / 2.;
2676 gp_Pnt xyz = surface->Value( uv.X(), uv.Y() );
2677 x = xyz.X(); y = xyz.Y(); z = xyz.Z();
2680 x = (Ns[i]->X() + Ns[i+2]->X())/2;
2681 y = (Ns[i]->Y() + Ns[i+2]->Y())/2;
2682 z = (Ns[i]->Z() + Ns[i+2]->Z())/2;
2684 if( fabs( Ns[i+1]->X() - x ) > disttol ||
2685 fabs( Ns[i+1]->Y() - y ) > disttol ||
2686 fabs( Ns[i+1]->Z() - z ) > disttol ) {
2687 // we have to move i+1 node
2688 aMesh->MoveNode( Ns[i+1], x, y, z );
2695 } // loop on face ids
2699 //=======================================================================
2700 //function : isReverse
2701 //purpose : Return true if normal of prevNodes is not co-directied with
2702 // gp_Vec(prevNodes[iNotSame],nextNodes[iNotSame]).
2703 // iNotSame is where prevNodes and nextNodes are different
2704 //=======================================================================
2706 static bool isReverse(vector<const SMDS_MeshNode*> prevNodes,
2707 vector<const SMDS_MeshNode*> nextNodes,
2711 int iBeforeNotSame = ( iNotSame == 0 ? nbNodes - 1 : iNotSame - 1 );
2712 int iAfterNotSame = ( iNotSame + 1 == nbNodes ? 0 : iNotSame + 1 );
2714 const SMDS_MeshNode* nB = prevNodes[ iBeforeNotSame ];
2715 const SMDS_MeshNode* nA = prevNodes[ iAfterNotSame ];
2716 const SMDS_MeshNode* nP = prevNodes[ iNotSame ];
2717 const SMDS_MeshNode* nN = nextNodes[ iNotSame ];
2719 gp_Pnt pB ( nB->X(), nB->Y(), nB->Z() );
2720 gp_Pnt pA ( nA->X(), nA->Y(), nA->Z() );
2721 gp_Pnt pP ( nP->X(), nP->Y(), nP->Z() );
2722 gp_Pnt pN ( nN->X(), nN->Y(), nN->Z() );
2724 gp_Vec vB ( pP, pB ), vA ( pP, pA ), vN ( pP, pN );
2726 return (vA ^ vB) * vN < 0.0;
2729 //=======================================================================
2731 * \brief Create elements by sweeping an element
2732 * \param elem - element to sweep
2733 * \param newNodesItVec - nodes generated from each node of the element
2734 * \param newElems - generated elements
2735 * \param nbSteps - number of sweeping steps
2736 * \param srcElements - to append elem for each generated element
2738 //=======================================================================
2740 void SMESH_MeshEditor::sweepElement(const SMDS_MeshElement* elem,
2741 const vector<TNodeOfNodeListMapItr> & newNodesItVec,
2742 list<const SMDS_MeshElement*>& newElems,
2744 SMESH_SequenceOfElemPtr& srcElements)
2746 SMESHDS_Mesh* aMesh = GetMeshDS();
2748 // Loop on elem nodes:
2749 // find new nodes and detect same nodes indices
2750 int nbNodes = elem->NbNodes();
2751 vector < list< const SMDS_MeshNode* >::const_iterator > itNN( nbNodes );
2752 vector<const SMDS_MeshNode*> prevNod( nbNodes );
2753 vector<const SMDS_MeshNode*> nextNod( nbNodes );
2754 vector<const SMDS_MeshNode*> midlNod( nbNodes );
2756 int iNode, nbSame = 0, iNotSameNode = 0, iSameNode = 0;
2757 vector<int> sames(nbNodes);
2759 //bool issimple[nbNodes];
2760 vector<bool> issimple(nbNodes);
2762 for ( iNode = 0; iNode < nbNodes; iNode++ ) {
2763 TNodeOfNodeListMapItr nnIt = newNodesItVec[ iNode ];
2764 const SMDS_MeshNode* node = nnIt->first;
2765 const list< const SMDS_MeshNode* > & listNewNodes = nnIt->second;
2766 if ( listNewNodes.empty() )
2769 if(listNewNodes.size()==nbSteps) {
2770 issimple[iNode] = true;
2773 issimple[iNode] = false;
2776 itNN[ iNode ] = listNewNodes.begin();
2777 prevNod[ iNode ] = node;
2778 nextNod[ iNode ] = listNewNodes.front();
2779 //cout<<"iNode="<<iNode<<endl;
2780 //cout<<" prevNod[iNode]="<< prevNod[iNode]<<" nextNod[iNode]="<< nextNod[iNode]<<endl;
2781 if ( prevNod[ iNode ] != nextNod [ iNode ])
2782 iNotSameNode = iNode;
2786 sames[nbSame++] = iNode;
2789 //cout<<"1 nbSame="<<nbSame<<endl;
2790 if ( nbSame == nbNodes || nbSame > 2) {
2791 MESSAGE( " Too many same nodes of element " << elem->GetID() );
2795 // if( elem->IsQuadratic() && nbSame>0 ) {
2796 // MESSAGE( "Can not rotate quadratic element " << elem->GetID() );
2800 int iBeforeSame = 0, iAfterSame = 0, iOpposSame = 0;
2802 iBeforeSame = ( iSameNode == 0 ? nbNodes - 1 : iSameNode - 1 );
2803 iAfterSame = ( iSameNode + 1 == nbNodes ? 0 : iSameNode + 1 );
2804 iOpposSame = ( iSameNode - 2 < 0 ? iSameNode + 2 : iSameNode - 2 );
2808 //cout<<" prevNod[0]="<< prevNod[0]<<" prevNod[1]="<< prevNod[1]
2809 // <<" prevNod[2]="<< prevNod[2]<<" prevNod[3]="<< prevNod[4]
2810 // <<" prevNod[4]="<< prevNod[4]<<" prevNod[5]="<< prevNod[5]
2811 // <<" prevNod[6]="<< prevNod[6]<<" prevNod[7]="<< prevNod[7]<<endl;
2813 // check element orientation
2815 if ( nbNodes > 2 && !isReverse( prevNod, nextNod, nbNodes, iNotSameNode )) {
2816 //MESSAGE("Reversed elem " << elem );
2820 int iAB = iAfterSame + iBeforeSame;
2821 iBeforeSame = iAB - iBeforeSame;
2822 iAfterSame = iAB - iAfterSame;
2826 // make new elements
2827 for (int iStep = 0; iStep < nbSteps; iStep++ ) {
2829 for ( iNode = 0; iNode < nbNodes; iNode++ ) {
2830 if(issimple[iNode]) {
2831 nextNod[ iNode ] = *itNN[ iNode ];
2835 if( elem->GetType()==SMDSAbs_Node ) {
2836 // we have to use two nodes
2837 midlNod[ iNode ] = *itNN[ iNode ];
2839 nextNod[ iNode ] = *itNN[ iNode ];
2842 else if(!elem->IsQuadratic() ||
2843 elem->IsQuadratic() && elem->IsMediumNode(prevNod[iNode]) ) {
2844 // we have to use each second node
2846 nextNod[ iNode ] = *itNN[ iNode ];
2850 // we have to use two nodes
2851 midlNod[ iNode ] = *itNN[ iNode ];
2853 nextNod[ iNode ] = *itNN[ iNode ];
2858 SMDS_MeshElement* aNewElem = 0;
2859 if(!elem->IsPoly()) {
2860 switch ( nbNodes ) {
2864 if ( nbSame == 0 ) {
2866 aNewElem = aMesh->AddEdge( prevNod[ 0 ], nextNod[ 0 ] );
2868 aNewElem = aMesh->AddEdge( prevNod[ 0 ], nextNod[ 0 ], midlNod[ 0 ] );
2874 aNewElem = aMesh->AddFace(prevNod[ 0 ], prevNod[ 1 ],
2875 nextNod[ 1 ], nextNod[ 0 ] );
2877 aNewElem = aMesh->AddFace(prevNod[ 0 ], prevNod[ 1 ],
2878 nextNod[ iNotSameNode ] );
2882 case 3: { // TRIANGLE or quadratic edge
2883 if(elem->GetType() == SMDSAbs_Face) { // TRIANGLE
2885 if ( nbSame == 0 ) // --- pentahedron
2886 aNewElem = aMesh->AddVolume (prevNod[ i0 ], prevNod[ 1 ], prevNod[ i2 ],
2887 nextNod[ i0 ], nextNod[ 1 ], nextNod[ i2 ] );
2889 else if ( nbSame == 1 ) // --- pyramid
2890 aNewElem = aMesh->AddVolume (prevNod[ iBeforeSame ], prevNod[ iAfterSame ],
2891 nextNod[ iAfterSame ], nextNod[ iBeforeSame ],
2892 nextNod[ iSameNode ]);
2894 else // 2 same nodes: --- tetrahedron
2895 aNewElem = aMesh->AddVolume (prevNod[ i0 ], prevNod[ 1 ], prevNod[ i2 ],
2896 nextNod[ iNotSameNode ]);
2898 else { // quadratic edge
2899 if(nbSame==0) { // quadratic quadrangle
2900 aNewElem = aMesh->AddFace(prevNod[0], nextNod[0], nextNod[1], prevNod[1],
2901 midlNod[0], nextNod[2], midlNod[1], prevNod[2]);
2903 else if(nbSame==1) { // quadratic triangle
2905 return; // medium node on axis
2906 else if(sames[0]==0) {
2907 aNewElem = aMesh->AddFace(prevNod[0], nextNod[1], prevNod[1],
2908 nextNod[2], midlNod[1], prevNod[2]);
2910 else { // sames[0]==1
2911 aNewElem = aMesh->AddFace(prevNod[0], nextNod[0], prevNod[1],
2912 midlNod[0], nextNod[2], prevNod[2]);
2920 case 4: { // QUADRANGLE
2922 if ( nbSame == 0 ) // --- hexahedron
2923 aNewElem = aMesh->AddVolume (prevNod[ i0 ], prevNod[ 1 ], prevNod[ i2 ], prevNod[ 3 ],
2924 nextNod[ i0 ], nextNod[ 1 ], nextNod[ i2 ], nextNod[ 3 ]);
2926 else if ( nbSame == 1 ) { // --- pyramid + pentahedron
2927 aNewElem = aMesh->AddVolume (prevNod[ iBeforeSame ], prevNod[ iAfterSame ],
2928 nextNod[ iAfterSame ], nextNod[ iBeforeSame ],
2929 nextNod[ iSameNode ]);
2930 newElems.push_back( aNewElem );
2931 aNewElem = aMesh->AddVolume (prevNod[ iAfterSame ], prevNod[ iOpposSame ],
2932 prevNod[ iBeforeSame ], nextNod[ iAfterSame ],
2933 nextNod[ iOpposSame ], nextNod[ iBeforeSame ] );
2935 else if ( nbSame == 2 ) { // pentahedron
2936 if ( prevNod[ iBeforeSame ] == nextNod[ iBeforeSame ] )
2937 // iBeforeSame is same too
2938 aNewElem = aMesh->AddVolume (prevNod[ iBeforeSame ], prevNod[ iOpposSame ],
2939 nextNod[ iOpposSame ], prevNod[ iSameNode ],
2940 prevNod[ iAfterSame ], nextNod[ iAfterSame ]);
2942 // iAfterSame is same too
2943 aNewElem = aMesh->AddVolume (prevNod[ iSameNode ], prevNod[ iBeforeSame ],
2944 nextNod[ iBeforeSame ], prevNod[ iAfterSame ],
2945 prevNod[ iOpposSame ], nextNod[ iOpposSame ]);
2949 case 6: { // quadratic triangle
2950 // create pentahedron with 15 nodes
2951 if(i0>0) { // reversed case
2952 aNewElem = aMesh->AddVolume (prevNod[0], prevNod[2], prevNod[1],
2953 nextNod[0], nextNod[2], nextNod[1],
2954 prevNod[5], prevNod[4], prevNod[3],
2955 nextNod[5], nextNod[4], nextNod[3],
2956 midlNod[0], midlNod[2], midlNod[1]);
2958 else { // not reversed case
2959 aNewElem = aMesh->AddVolume (prevNod[0], prevNod[1], prevNod[2],
2960 nextNod[0], nextNod[1], nextNod[2],
2961 prevNod[3], prevNod[4], prevNod[5],
2962 nextNod[3], nextNod[4], nextNod[5],
2963 midlNod[0], midlNod[1], midlNod[2]);
2967 case 8: { // quadratic quadrangle
2968 // create hexahedron with 20 nodes
2969 if(i0>0) { // reversed case
2970 aNewElem = aMesh->AddVolume (prevNod[0], prevNod[3], prevNod[2], prevNod[1],
2971 nextNod[0], nextNod[3], nextNod[2], nextNod[1],
2972 prevNod[7], prevNod[6], prevNod[5], prevNod[4],
2973 nextNod[7], nextNod[6], nextNod[5], nextNod[4],
2974 midlNod[0], midlNod[3], midlNod[2], midlNod[1]);
2976 else { // not reversed case
2977 aNewElem = aMesh->AddVolume (prevNod[0], prevNod[1], prevNod[2], prevNod[3],
2978 nextNod[0], nextNod[1], nextNod[2], nextNod[3],
2979 prevNod[4], prevNod[5], prevNod[6], prevNod[7],
2980 nextNod[4], nextNod[5], nextNod[6], nextNod[7],
2981 midlNod[0], midlNod[1], midlNod[2], midlNod[3]);
2986 // realized for extrusion only
2987 //vector<const SMDS_MeshNode*> polyedre_nodes (nbNodes*2 + 4*nbNodes);
2988 //vector<int> quantities (nbNodes + 2);
2990 //quantities[0] = nbNodes; // bottom of prism
2991 //for (int inode = 0; inode < nbNodes; inode++) {
2992 // polyedre_nodes[inode] = prevNod[inode];
2995 //quantities[1] = nbNodes; // top of prism
2996 //for (int inode = 0; inode < nbNodes; inode++) {
2997 // polyedre_nodes[nbNodes + inode] = nextNod[inode];
3000 //for (int iface = 0; iface < nbNodes; iface++) {
3001 // quantities[iface + 2] = 4;
3002 // int inextface = (iface == nbNodes - 1) ? 0 : iface + 1;
3003 // polyedre_nodes[2*nbNodes + 4*iface + 0] = prevNod[iface];
3004 // polyedre_nodes[2*nbNodes + 4*iface + 1] = prevNod[inextface];
3005 // polyedre_nodes[2*nbNodes + 4*iface + 2] = nextNod[inextface];
3006 // polyedre_nodes[2*nbNodes + 4*iface + 3] = nextNod[iface];
3008 //aNewElem = aMesh->AddPolyhedralVolume (polyedre_nodes, quantities);
3015 // realized for extrusion only
3016 vector<const SMDS_MeshNode*> polyedre_nodes (nbNodes*2 + 4*nbNodes);
3017 vector<int> quantities (nbNodes + 2);
3019 quantities[0] = nbNodes; // bottom of prism
3020 for (int inode = 0; inode < nbNodes; inode++) {
3021 polyedre_nodes[inode] = prevNod[inode];
3024 quantities[1] = nbNodes; // top of prism
3025 for (int inode = 0; inode < nbNodes; inode++) {
3026 polyedre_nodes[nbNodes + inode] = nextNod[inode];
3029 for (int iface = 0; iface < nbNodes; iface++) {
3030 quantities[iface + 2] = 4;
3031 int inextface = (iface == nbNodes - 1) ? 0 : iface + 1;
3032 polyedre_nodes[2*nbNodes + 4*iface + 0] = prevNod[iface];
3033 polyedre_nodes[2*nbNodes + 4*iface + 1] = prevNod[inextface];
3034 polyedre_nodes[2*nbNodes + 4*iface + 2] = nextNod[inextface];
3035 polyedre_nodes[2*nbNodes + 4*iface + 3] = nextNod[iface];
3037 aNewElem = aMesh->AddPolyhedralVolume (polyedre_nodes, quantities);
3041 newElems.push_back( aNewElem );
3042 myLastCreatedElems.Append(aNewElem);
3043 srcElements.Append( elem );
3046 // set new prev nodes
3047 for ( iNode = 0; iNode < nbNodes; iNode++ )
3048 prevNod[ iNode ] = nextNod[ iNode ];
3053 //=======================================================================
3055 * \brief Create 1D and 2D elements around swept elements
3056 * \param mapNewNodes - source nodes and ones generated from them
3057 * \param newElemsMap - source elements and ones generated from them
3058 * \param elemNewNodesMap - nodes generated from each node of each element
3059 * \param elemSet - all swept elements
3060 * \param nbSteps - number of sweeping steps
3061 * \param srcElements - to append elem for each generated element
3063 //=======================================================================
3065 void SMESH_MeshEditor::makeWalls (TNodeOfNodeListMap & mapNewNodes,
3066 TElemOfElemListMap & newElemsMap,
3067 TElemOfVecOfNnlmiMap & elemNewNodesMap,
3068 TIDSortedElemSet& elemSet,
3070 SMESH_SequenceOfElemPtr& srcElements)
3072 ASSERT( newElemsMap.size() == elemNewNodesMap.size() );
3073 SMESHDS_Mesh* aMesh = GetMeshDS();
3075 // Find nodes belonging to only one initial element - sweep them to get edges.
3077 TNodeOfNodeListMapItr nList = mapNewNodes.begin();
3078 for ( ; nList != mapNewNodes.end(); nList++ ) {
3079 const SMDS_MeshNode* node =
3080 static_cast<const SMDS_MeshNode*>( nList->first );
3081 SMDS_ElemIteratorPtr eIt = node->GetInverseElementIterator();
3082 int nbInitElems = 0;
3083 const SMDS_MeshElement* el = 0;
3084 SMDSAbs_ElementType highType = SMDSAbs_Edge; // count most complex elements only
3085 while ( eIt->more() && nbInitElems < 2 ) {
3087 SMDSAbs_ElementType type = el->GetType();
3088 if ( type == SMDSAbs_Volume || type < highType ) continue;
3089 if ( type > highType ) {
3093 if ( elemSet.find(el) != elemSet.end() )
3096 if ( nbInitElems < 2 ) {
3097 bool NotCreateEdge = el && el->IsQuadratic() && el->IsMediumNode(node);
3098 if(!NotCreateEdge) {
3099 vector<TNodeOfNodeListMapItr> newNodesItVec( 1, nList );
3100 list<const SMDS_MeshElement*> newEdges;
3101 sweepElement( node, newNodesItVec, newEdges, nbSteps, srcElements );
3106 // Make a ceiling for each element ie an equal element of last new nodes.
3107 // Find free links of faces - make edges and sweep them into faces.
3109 TElemOfElemListMap::iterator itElem = newElemsMap.begin();
3110 TElemOfVecOfNnlmiMap::iterator itElemNodes = elemNewNodesMap.begin();
3111 for ( ; itElem != newElemsMap.end(); itElem++, itElemNodes++ ) {
3112 const SMDS_MeshElement* elem = itElem->first;
3113 vector<TNodeOfNodeListMapItr>& vecNewNodes = itElemNodes->second;
3115 if ( elem->GetType() == SMDSAbs_Edge ) {
3116 // create a ceiling edge
3117 if (!elem->IsQuadratic()) {
3118 if ( !aMesh->FindEdge( vecNewNodes[ 0 ]->second.back(),
3119 vecNewNodes[ 1 ]->second.back())) {
3120 myLastCreatedElems.Append(aMesh->AddEdge(vecNewNodes[ 0 ]->second.back(),
3121 vecNewNodes[ 1 ]->second.back()));
3122 srcElements.Append( myLastCreatedElems.Last() );
3126 if ( !aMesh->FindEdge( vecNewNodes[ 0 ]->second.back(),
3127 vecNewNodes[ 1 ]->second.back(),
3128 vecNewNodes[ 2 ]->second.back())) {
3129 myLastCreatedElems.Append(aMesh->AddEdge(vecNewNodes[ 0 ]->second.back(),
3130 vecNewNodes[ 1 ]->second.back(),
3131 vecNewNodes[ 2 ]->second.back()));
3132 srcElements.Append( myLastCreatedElems.Last() );
3136 if ( elem->GetType() != SMDSAbs_Face )
3139 if(itElem->second.size()==0) continue;
3141 bool hasFreeLinks = false;
3143 TIDSortedElemSet avoidSet;
3144 avoidSet.insert( elem );
3146 set<const SMDS_MeshNode*> aFaceLastNodes;
3147 int iNode, nbNodes = vecNewNodes.size();
3148 if(!elem->IsQuadratic()) {
3149 // loop on the face nodes
3150 for ( iNode = 0; iNode < nbNodes; iNode++ ) {
3151 aFaceLastNodes.insert( vecNewNodes[ iNode ]->second.back() );
3152 // look for free links of the face
3153 int iNext = ( iNode + 1 == nbNodes ) ? 0 : iNode + 1;
3154 const SMDS_MeshNode* n1 = vecNewNodes[ iNode ]->first;
3155 const SMDS_MeshNode* n2 = vecNewNodes[ iNext ]->first;
3156 // check if a link is free
3157 if ( ! SMESH_MeshEditor::FindFaceInSet ( n1, n2, elemSet, avoidSet )) {
3158 hasFreeLinks = true;
3159 // make an edge and a ceiling for a new edge
3160 if ( !aMesh->FindEdge( n1, n2 )) {
3161 myLastCreatedElems.Append(aMesh->AddEdge( n1, n2 )); // free link edge
3162 srcElements.Append( myLastCreatedElems.Last() );
3164 n1 = vecNewNodes[ iNode ]->second.back();
3165 n2 = vecNewNodes[ iNext ]->second.back();
3166 if ( !aMesh->FindEdge( n1, n2 )) {
3167 myLastCreatedElems.Append(aMesh->AddEdge( n1, n2 )); // ceiling edge
3168 srcElements.Append( myLastCreatedElems.Last() );
3173 else { // elem is quadratic face
3174 int nbn = nbNodes/2;
3175 for ( iNode = 0; iNode < nbn; iNode++ ) {
3176 aFaceLastNodes.insert( vecNewNodes[ iNode ]->second.back() );
3177 int iNext = ( iNode + 1 == nbn ) ? 0 : iNode + 1;
3178 const SMDS_MeshNode* n1 = vecNewNodes[ iNode ]->first;
3179 const SMDS_MeshNode* n2 = vecNewNodes[ iNext ]->first;
3180 // check if a link is free
3181 if ( ! SMESH_MeshEditor::FindFaceInSet ( n1, n2, elemSet, avoidSet )) {
3182 hasFreeLinks = true;
3183 // make an edge and a ceiling for a new edge
3185 const SMDS_MeshNode* n3 = vecNewNodes[ iNode+nbn ]->first;
3186 if ( !aMesh->FindEdge( n1, n2, n3 )) {
3187 myLastCreatedElems.Append(aMesh->AddEdge( n1, n2, n3 )); // free link edge
3188 srcElements.Append( myLastCreatedElems.Last() );
3190 n1 = vecNewNodes[ iNode ]->second.back();
3191 n2 = vecNewNodes[ iNext ]->second.back();
3192 n3 = vecNewNodes[ iNode+nbn ]->second.back();
3193 if ( !aMesh->FindEdge( n1, n2, n3 )) {
3194 myLastCreatedElems.Append(aMesh->AddEdge( n1, n2, n3 )); // ceiling edge
3195 srcElements.Append( myLastCreatedElems.Last() );
3199 for ( iNode = nbn; iNode < 2*nbn; iNode++ ) {
3200 aFaceLastNodes.insert( vecNewNodes[ iNode ]->second.back() );
3204 // sweep free links into faces
3206 if ( hasFreeLinks ) {
3207 list<const SMDS_MeshElement*> & newVolumes = itElem->second;
3208 int iVol, volNb, nbVolumesByStep = newVolumes.size() / nbSteps;
3210 set<const SMDS_MeshNode*> initNodeSet, topNodeSet, faceNodeSet;
3211 for ( iNode = 0; iNode < nbNodes; iNode++ ) {
3212 initNodeSet.insert( vecNewNodes[ iNode ]->first );
3213 topNodeSet .insert( vecNewNodes[ iNode ]->second.back() );
3215 for ( volNb = 0; volNb < nbVolumesByStep; volNb++ ) {
3216 list<const SMDS_MeshElement*>::iterator v = newVolumes.begin();
3218 while ( iVol++ < volNb ) v++;
3219 // find indices of free faces of a volume and their source edges
3220 list< int > freeInd;
3221 list< const SMDS_MeshElement* > srcEdges; // source edges of free faces
3222 SMDS_VolumeTool vTool( *v );
3223 int iF, nbF = vTool.NbFaces();
3224 for ( iF = 0; iF < nbF; iF ++ ) {
3225 if (vTool.IsFreeFace( iF ) &&
3226 vTool.GetFaceNodes( iF, faceNodeSet ) &&
3227 initNodeSet != faceNodeSet) // except an initial face
3229 if ( nbSteps == 1 && faceNodeSet == topNodeSet )
3231 freeInd.push_back( iF );
3232 // find source edge of a free face iF
3233 vector<const SMDS_MeshNode*> commonNodes; // shared by the initial and free faces
3234 commonNodes.resize( initNodeSet.size(), NULL ); // avoid spoiling memory
3235 std::set_intersection( faceNodeSet.begin(), faceNodeSet.end(),
3236 initNodeSet.begin(), initNodeSet.end(),
3237 commonNodes.begin());
3238 if ( (*v)->IsQuadratic() )
3239 srcEdges.push_back(aMesh->FindEdge (commonNodes[0],commonNodes[1],commonNodes[2]));
3241 srcEdges.push_back(aMesh->FindEdge (commonNodes[0],commonNodes[1]));
3243 if ( !srcEdges.back() )
3245 cout << "SMESH_MeshEditor::makeWalls(), no source edge found for a free face #"
3246 << iF << " of volume #" << vTool.ID() << endl;
3251 if ( freeInd.empty() )
3254 // create faces for all steps;
3255 // if such a face has been already created by sweep of edge,
3256 // assure that its orientation is OK
3257 for ( int iStep = 0; iStep < nbSteps; iStep++ ) {
3259 vTool.SetExternalNormal();
3260 list< int >::iterator ind = freeInd.begin();
3261 list< const SMDS_MeshElement* >::iterator srcEdge = srcEdges.begin();
3262 for ( ; ind != freeInd.end(); ++ind, ++srcEdge ) // loop on free faces
3264 const SMDS_MeshNode** nodes = vTool.GetFaceNodes( *ind );
3265 int nbn = vTool.NbFaceNodes( *ind );
3267 case 3: { ///// triangle
3268 const SMDS_MeshFace * f = aMesh->FindFace( nodes[ 0 ], nodes[ 1 ], nodes[ 2 ]);
3270 myLastCreatedElems.Append(aMesh->AddFace( nodes[ 0 ], nodes[ 1 ], nodes[ 2 ] ));
3271 else if ( nodes[ 1 ] != f->GetNode( f->GetNodeIndex( nodes[ 0 ] ) + 1 ))
3272 aMesh->ChangeElementNodes( f, nodes, nbn );
3275 case 4: { ///// quadrangle
3276 const SMDS_MeshFace * f = aMesh->FindFace( nodes[ 0 ], nodes[ 1 ], nodes[ 2 ], nodes[ 3 ]);
3278 myLastCreatedElems.Append(aMesh->AddFace( nodes[ 0 ], nodes[ 1 ], nodes[ 2 ], nodes[ 3 ] ));
3279 else if ( nodes[ 1 ] != f->GetNode( f->GetNodeIndex( nodes[ 0 ] ) + 1 ))
3280 aMesh->ChangeElementNodes( f, nodes, nbn );
3284 if( (*v)->IsQuadratic() ) {
3285 if(nbn==6) { /////// quadratic triangle
3286 const SMDS_MeshFace * f = aMesh->FindFace( nodes[0], nodes[2], nodes[4],
3287 nodes[1], nodes[3], nodes[5] );
3289 myLastCreatedElems.Append(aMesh->AddFace(nodes[0], nodes[2], nodes[4],
3290 nodes[1], nodes[3], nodes[5]));
3291 else if ( nodes[ 2 ] != f->GetNode( f->GetNodeIndex( nodes[ 0 ] ) + 1 ))
3292 aMesh->ChangeElementNodes( f, nodes, nbn );
3294 else { /////// quadratic quadrangle
3295 const SMDS_MeshFace * f = aMesh->FindFace( nodes[0], nodes[2], nodes[4], nodes[6],
3296 nodes[1], nodes[3], nodes[5], nodes[7] );
3298 myLastCreatedElems.Append(aMesh->AddFace(nodes[0], nodes[2], nodes[4], nodes[6],
3299 nodes[1], nodes[3], nodes[5], nodes[7]));
3300 else if ( nodes[ 2 ] != f->GetNode( f->GetNodeIndex( nodes[ 0 ] ) + 1 ))
3301 aMesh->ChangeElementNodes( f, nodes, nbn );
3304 else { //////// polygon
3305 vector<const SMDS_MeshNode*> polygon_nodes ( nodes, &nodes[nbn] );
3306 const SMDS_MeshFace * f = aMesh->FindFace( polygon_nodes );
3308 myLastCreatedElems.Append(aMesh->AddPolygonalFace(polygon_nodes));
3309 else if ( nodes[ 1 ] != f->GetNode( f->GetNodeIndex( nodes[ 0 ] ) + 1 ))
3310 aMesh->ChangeElementNodes( f, nodes, nbn );
3313 while ( srcElements.Length() < myLastCreatedElems.Length() )
3314 srcElements.Append( *srcEdge );
3316 } // loop on free faces
3318 // go to the next volume
3320 while ( iVol++ < nbVolumesByStep ) v++;
3323 } // sweep free links into faces
3325 // Make a ceiling face with a normal external to a volume
3327 SMDS_VolumeTool lastVol( itElem->second.back() );
3329 int iF = lastVol.GetFaceIndex( aFaceLastNodes );
3331 lastVol.SetExternalNormal();
3332 const SMDS_MeshNode** nodes = lastVol.GetFaceNodes( iF );
3333 int nbn = lastVol.NbFaceNodes( iF );
3336 if (!hasFreeLinks ||
3337 !aMesh->FindFace( nodes[ 0 ], nodes[ 1 ], nodes[ 2 ]))
3338 myLastCreatedElems.Append(aMesh->AddFace( nodes[ 0 ], nodes[ 1 ], nodes[ 2 ] ));
3341 if (!hasFreeLinks ||
3342 !aMesh->FindFace( nodes[ 0 ], nodes[ 1 ], nodes[ 2 ], nodes[ 3 ]))
3343 myLastCreatedElems.Append(aMesh->AddFace( nodes[ 0 ], nodes[ 1 ], nodes[ 2 ], nodes[ 3 ] ));
3346 if(itElem->second.back()->IsQuadratic()) {
3348 if (!hasFreeLinks ||
3349 !aMesh->FindFace(nodes[0], nodes[2], nodes[4],
3350 nodes[1], nodes[3], nodes[5]) ) {
3351 myLastCreatedElems.Append(aMesh->AddFace(nodes[0], nodes[2], nodes[4],
3352 nodes[1], nodes[3], nodes[5]));
3356 if (!hasFreeLinks ||
3357 !aMesh->FindFace(nodes[0], nodes[2], nodes[4], nodes[6],
3358 nodes[1], nodes[3], nodes[5], nodes[7]) )
3359 myLastCreatedElems.Append(aMesh->AddFace(nodes[0], nodes[2], nodes[4], nodes[6],
3360 nodes[1], nodes[3], nodes[5], nodes[7]));
3364 vector<const SMDS_MeshNode*> polygon_nodes ( nodes, &nodes[nbn] );
3365 if (!hasFreeLinks || !aMesh->FindFace(polygon_nodes))
3366 myLastCreatedElems.Append(aMesh->AddPolygonalFace(polygon_nodes));
3370 while ( srcElements.Length() < myLastCreatedElems.Length() )
3371 srcElements.Append( myLastCreatedElems.Last() );
3373 } // loop on swept elements
3376 //=======================================================================
3377 //function : RotationSweep
3379 //=======================================================================
3381 SMESH_MeshEditor::PGroupIDs
3382 SMESH_MeshEditor::RotationSweep(TIDSortedElemSet & theElems,
3383 const gp_Ax1& theAxis,
3384 const double theAngle,
3385 const int theNbSteps,
3386 const double theTol,
3387 const bool theMakeGroups,
3388 const bool theMakeWalls)
3390 myLastCreatedElems.Clear();
3391 myLastCreatedNodes.Clear();
3393 // source elements for each generated one
3394 SMESH_SequenceOfElemPtr srcElems, srcNodes;
3396 MESSAGE( "RotationSweep()");
3398 aTrsf.SetRotation( theAxis, theAngle );
3400 aTrsf2.SetRotation( theAxis, theAngle/2. );
3402 gp_Lin aLine( theAxis );
3403 double aSqTol = theTol * theTol;
3405 SMESHDS_Mesh* aMesh = GetMeshDS();
3407 TNodeOfNodeListMap mapNewNodes;
3408 TElemOfVecOfNnlmiMap mapElemNewNodes;
3409 TElemOfElemListMap newElemsMap;
3412 TIDSortedElemSet::iterator itElem;
3413 for ( itElem = theElems.begin(); itElem != theElems.end(); itElem++ ) {
3414 const SMDS_MeshElement* elem = *itElem;
3415 if ( !elem || elem->GetType() == SMDSAbs_Volume )
3417 vector<TNodeOfNodeListMapItr> & newNodesItVec = mapElemNewNodes[ elem ];
3418 newNodesItVec.reserve( elem->NbNodes() );
3420 // loop on elem nodes
3421 SMDS_ElemIteratorPtr itN = elem->nodesIterator();
3422 while ( itN->more() )
3424 // check if a node has been already sweeped
3425 const SMDS_MeshNode* node = cast2Node( itN->next() );
3426 TNodeOfNodeListMapItr nIt = mapNewNodes.find( node );
3427 if ( nIt == mapNewNodes.end() ) {
3428 nIt = mapNewNodes.insert( make_pair( node, list<const SMDS_MeshNode*>() )).first;
3429 list<const SMDS_MeshNode*>& listNewNodes = nIt->second;
3432 gp_XYZ aXYZ( node->X(), node->Y(), node->Z() );
3434 aXYZ.Coord( coord[0], coord[1], coord[2] );
3435 bool isOnAxis = ( aLine.SquareDistance( aXYZ ) <= aSqTol );
3436 const SMDS_MeshNode * newNode = node;
3437 for ( int i = 0; i < theNbSteps; i++ ) {
3439 if( elem->IsQuadratic() && !elem->IsMediumNode(node) ) {
3441 aTrsf2.Transforms( coord[0], coord[1], coord[2] );
3442 //aTrsf.Transforms( coord[0], coord[1], coord[2] );
3443 newNode = aMesh->AddNode( coord[0], coord[1], coord[2] );
3444 myLastCreatedNodes.Append(newNode);
3445 srcNodes.Append( node );
3446 listNewNodes.push_back( newNode );
3447 aTrsf2.Transforms( coord[0], coord[1], coord[2] );
3448 //aTrsf.Transforms( coord[0], coord[1], coord[2] );
3451 aTrsf.Transforms( coord[0], coord[1], coord[2] );
3453 newNode = aMesh->AddNode( coord[0], coord[1], coord[2] );
3454 myLastCreatedNodes.Append(newNode);
3455 srcNodes.Append( node );
3457 listNewNodes.push_back( newNode );
3461 // if current elem is quadratic and current node is not medium
3462 // we have to check - may be it is needed to insert additional nodes
3463 if( elem->IsQuadratic() && !elem->IsMediumNode(node) ) {
3464 list< const SMDS_MeshNode* > & listNewNodes = nIt->second;
3465 if(listNewNodes.size()==theNbSteps) {
3466 listNewNodes.clear();
3468 gp_XYZ aXYZ( node->X(), node->Y(), node->Z() );
3470 aXYZ.Coord( coord[0], coord[1], coord[2] );
3471 const SMDS_MeshNode * newNode = node;
3472 for(int i = 0; i<theNbSteps; i++) {
3473 aTrsf2.Transforms( coord[0], coord[1], coord[2] );
3474 newNode = aMesh->AddNode( coord[0], coord[1], coord[2] );
3475 myLastCreatedNodes.Append(newNode);
3476 listNewNodes.push_back( newNode );
3477 srcNodes.Append( node );
3478 aTrsf2.Transforms( coord[0], coord[1], coord[2] );
3479 newNode = aMesh->AddNode( coord[0], coord[1], coord[2] );
3480 myLastCreatedNodes.Append(newNode);
3481 srcNodes.Append( node );
3482 listNewNodes.push_back( newNode );
3487 newNodesItVec.push_back( nIt );
3489 // make new elements
3490 sweepElement( elem, newNodesItVec, newElemsMap[elem], theNbSteps, srcElems );
3494 makeWalls( mapNewNodes, newElemsMap, mapElemNewNodes, theElems, theNbSteps, srcElems );
3496 PGroupIDs newGroupIDs;
3497 if ( theMakeGroups )
3498 newGroupIDs = generateGroups( srcNodes, srcElems, "rotated");
3504 //=======================================================================
3505 //function : CreateNode
3507 //=======================================================================
3508 const SMDS_MeshNode* SMESH_MeshEditor::CreateNode(const double x,
3511 const double tolnode,
3512 SMESH_SequenceOfNode& aNodes)
3514 myLastCreatedElems.Clear();
3515 myLastCreatedNodes.Clear();
3518 SMESHDS_Mesh * aMesh = myMesh->GetMeshDS();
3520 // try to search in sequence of existing nodes
3521 // if aNodes.Length()>0 we 'nave to use given sequence
3522 // else - use all nodes of mesh
3523 if(aNodes.Length()>0) {
3525 for(i=1; i<=aNodes.Length(); i++) {
3526 gp_Pnt P2(aNodes.Value(i)->X(),aNodes.Value(i)->Y(),aNodes.Value(i)->Z());
3527 if(P1.Distance(P2)<tolnode)
3528 return aNodes.Value(i);
3532 SMDS_NodeIteratorPtr itn = aMesh->nodesIterator();
3533 while(itn->more()) {
3534 const SMDS_MeshNode* aN = static_cast<const SMDS_MeshNode*> (itn->next());
3535 gp_Pnt P2(aN->X(),aN->Y(),aN->Z());
3536 if(P1.Distance(P2)<tolnode)
3541 // create new node and return it
3542 const SMDS_MeshNode* NewNode = aMesh->AddNode(x,y,z);
3543 myLastCreatedNodes.Append(NewNode);
3548 //=======================================================================
3549 //function : ExtrusionSweep
3551 //=======================================================================
3553 SMESH_MeshEditor::PGroupIDs
3554 SMESH_MeshEditor::ExtrusionSweep (TIDSortedElemSet & theElems,
3555 const gp_Vec& theStep,
3556 const int theNbSteps,
3557 TElemOfElemListMap& newElemsMap,
3558 const bool theMakeGroups,
3560 const double theTolerance)
3562 ExtrusParam aParams;
3563 aParams.myDir = gp_Dir(theStep);
3564 aParams.myNodes.Clear();
3565 aParams.mySteps = new TColStd_HSequenceOfReal;
3567 for(i=1; i<=theNbSteps; i++)
3568 aParams.mySteps->Append(theStep.Magnitude());
3571 ExtrusionSweep(theElems,aParams,newElemsMap,theMakeGroups,theFlags,theTolerance);
3575 //=======================================================================
3576 //function : ExtrusionSweep
3578 //=======================================================================
3580 SMESH_MeshEditor::PGroupIDs
3581 SMESH_MeshEditor::ExtrusionSweep (TIDSortedElemSet & theElems,
3582 ExtrusParam& theParams,
3583 TElemOfElemListMap& newElemsMap,
3584 const bool theMakeGroups,
3586 const double theTolerance)
3588 myLastCreatedElems.Clear();
3589 myLastCreatedNodes.Clear();
3591 // source elements for each generated one
3592 SMESH_SequenceOfElemPtr srcElems, srcNodes;
3594 SMESHDS_Mesh* aMesh = GetMeshDS();
3596 int nbsteps = theParams.mySteps->Length();
3598 TNodeOfNodeListMap mapNewNodes;
3599 //TNodeOfNodeVecMap mapNewNodes;
3600 TElemOfVecOfNnlmiMap mapElemNewNodes;
3601 //TElemOfVecOfMapNodesMap mapElemNewNodes;
3604 TIDSortedElemSet::iterator itElem;
3605 for ( itElem = theElems.begin(); itElem != theElems.end(); itElem++ ) {
3606 // check element type
3607 const SMDS_MeshElement* elem = *itElem;
3608 if ( !elem || elem->GetType() == SMDSAbs_Volume )
3611 vector<TNodeOfNodeListMapItr> & newNodesItVec = mapElemNewNodes[ elem ];
3612 //vector<TNodeOfNodeVecMapItr> & newNodesItVec = mapElemNewNodes[ elem ];
3613 newNodesItVec.reserve( elem->NbNodes() );
3615 // loop on elem nodes
3616 SMDS_ElemIteratorPtr itN = elem->nodesIterator();
3617 while ( itN->more() )
3619 // check if a node has been already sweeped
3620 const SMDS_MeshNode* node = cast2Node( itN->next() );
3621 TNodeOfNodeListMap::iterator nIt = mapNewNodes.find( node );
3622 //TNodeOfNodeVecMap::iterator nIt = mapNewNodes.find( node );
3623 if ( nIt == mapNewNodes.end() ) {
3624 nIt = mapNewNodes.insert( make_pair( node, list<const SMDS_MeshNode*>() )).first;
3625 //nIt = mapNewNodes.insert( make_pair( node, vector<const SMDS_MeshNode*>() )).first;
3626 list<const SMDS_MeshNode*>& listNewNodes = nIt->second;
3627 //vector<const SMDS_MeshNode*>& vecNewNodes = nIt->second;
3628 //vecNewNodes.reserve(nbsteps);
3631 double coord[] = { node->X(), node->Y(), node->Z() };
3632 //int nbsteps = theParams.mySteps->Length();
3633 for ( int i = 0; i < nbsteps; i++ ) {
3634 if( elem->IsQuadratic() && !elem->IsMediumNode(node) ) {
3635 // create additional node
3636 double x = coord[0] + theParams.myDir.X()*theParams.mySteps->Value(i+1)/2.;
3637 double y = coord[1] + theParams.myDir.Y()*theParams.mySteps->Value(i+1)/2.;
3638 double z = coord[2] + theParams.myDir.Z()*theParams.mySteps->Value(i+1)/2.;
3639 if( theFlags & EXTRUSION_FLAG_SEW ) {
3640 const SMDS_MeshNode * newNode = CreateNode(x, y, z,
3641 theTolerance, theParams.myNodes);
3642 listNewNodes.push_back( newNode );
3645 const SMDS_MeshNode * newNode = aMesh->AddNode(x, y, z);
3646 myLastCreatedNodes.Append(newNode);
3647 srcNodes.Append( node );
3648 listNewNodes.push_back( newNode );
3651 //aTrsf.Transforms( coord[0], coord[1], coord[2] );
3652 coord[0] = coord[0] + theParams.myDir.X()*theParams.mySteps->Value(i+1);
3653 coord[1] = coord[1] + theParams.myDir.Y()*theParams.mySteps->Value(i+1);
3654 coord[2] = coord[2] + theParams.myDir.Z()*theParams.mySteps->Value(i+1);
3655 if( theFlags & EXTRUSION_FLAG_SEW ) {
3656 const SMDS_MeshNode * newNode = CreateNode(coord[0], coord[1], coord[2],
3657 theTolerance, theParams.myNodes);
3658 listNewNodes.push_back( newNode );
3659 //vecNewNodes[i]=newNode;
3662 const SMDS_MeshNode * newNode = aMesh->AddNode( coord[0], coord[1], coord[2] );
3663 myLastCreatedNodes.Append(newNode);
3664 srcNodes.Append( node );
3665 listNewNodes.push_back( newNode );
3666 //vecNewNodes[i]=newNode;
3671 // if current elem is quadratic and current node is not medium
3672 // we have to check - may be it is needed to insert additional nodes
3673 if( elem->IsQuadratic() && !elem->IsMediumNode(node) ) {
3674 list< const SMDS_MeshNode* > & listNewNodes = nIt->second;
3675 if(listNewNodes.size()==nbsteps) {
3676 listNewNodes.clear();
3677 double coord[] = { node->X(), node->Y(), node->Z() };
3678 for ( int i = 0; i < nbsteps; i++ ) {
3679 double x = coord[0] + theParams.myDir.X()*theParams.mySteps->Value(i+1);
3680 double y = coord[1] + theParams.myDir.Y()*theParams.mySteps->Value(i+1);
3681 double z = coord[2] + theParams.myDir.Z()*theParams.mySteps->Value(i+1);
3682 if( theFlags & EXTRUSION_FLAG_SEW ) {
3683 const SMDS_MeshNode * newNode = CreateNode(x, y, z,
3684 theTolerance, theParams.myNodes);
3685 listNewNodes.push_back( newNode );
3688 const SMDS_MeshNode * newNode = aMesh->AddNode(x, y, z);
3689 myLastCreatedNodes.Append(newNode);
3690 srcNodes.Append( node );
3691 listNewNodes.push_back( newNode );
3693 coord[0] = coord[0] + theParams.myDir.X()*theParams.mySteps->Value(i+1);
3694 coord[1] = coord[1] + theParams.myDir.Y()*theParams.mySteps->Value(i+1);
3695 coord[2] = coord[2] + theParams.myDir.Z()*theParams.mySteps->Value(i+1);
3696 if( theFlags & EXTRUSION_FLAG_SEW ) {
3697 const SMDS_MeshNode * newNode = CreateNode(coord[0], coord[1], coord[2],
3698 theTolerance, theParams.myNodes);
3699 listNewNodes.push_back( newNode );
3702 const SMDS_MeshNode * newNode = aMesh->AddNode( coord[0], coord[1], coord[2] );
3703 myLastCreatedNodes.Append(newNode);
3704 srcNodes.Append( node );
3705 listNewNodes.push_back( newNode );
3711 newNodesItVec.push_back( nIt );
3713 // make new elements
3714 sweepElement( elem, newNodesItVec, newElemsMap[elem], nbsteps, srcElems );
3717 if( theFlags & EXTRUSION_FLAG_BOUNDARY ) {
3718 makeWalls( mapNewNodes, newElemsMap, mapElemNewNodes, theElems, nbsteps, srcElems );
3720 PGroupIDs newGroupIDs;
3721 if ( theMakeGroups )
3722 newGroupIDs = generateGroups( srcNodes, srcElems, "extruded");
3728 //=======================================================================
3729 //class : SMESH_MeshEditor_PathPoint
3730 //purpose : auxiliary class
3731 //=======================================================================
3732 class SMESH_MeshEditor_PathPoint {
3734 SMESH_MeshEditor_PathPoint() {
3735 myPnt.SetCoord(99., 99., 99.);
3736 myTgt.SetCoord(1.,0.,0.);
3740 void SetPnt(const gp_Pnt& aP3D){
3743 void SetTangent(const gp_Dir& aTgt){
3746 void SetAngle(const double& aBeta){
3749 void SetParameter(const double& aPrm){
3752 const gp_Pnt& Pnt()const{
3755 const gp_Dir& Tangent()const{
3758 double Angle()const{
3761 double Parameter()const{
3772 //=======================================================================
3773 //function : ExtrusionAlongTrack
3775 //=======================================================================
3776 SMESH_MeshEditor::Extrusion_Error
3777 SMESH_MeshEditor::ExtrusionAlongTrack (TIDSortedElemSet & theElements,
3778 SMESH_subMesh* theTrack,
3779 const SMDS_MeshNode* theN1,
3780 const bool theHasAngles,
3781 list<double>& theAngles,
3782 const bool theHasRefPoint,
3783 const gp_Pnt& theRefPoint,
3784 const bool theMakeGroups)
3786 myLastCreatedElems.Clear();
3787 myLastCreatedNodes.Clear();
3789 // source elements for each generated one
3790 SMESH_SequenceOfElemPtr srcElems, srcNodes;
3792 int j, aNbTP, aNbE, aNb;
3793 double aT1, aT2, aT, aAngle, aX, aY, aZ;
3794 std::list<double> aPrms;
3795 std::list<double>::iterator aItD;
3796 TIDSortedElemSet::iterator itElem;
3798 Standard_Real aTx1, aTx2, aL2, aTolVec, aTolVec2;
3802 Handle(Geom_Curve) aC3D;
3803 TopoDS_Edge aTrackEdge;
3804 TopoDS_Vertex aV1, aV2;
3806 SMDS_ElemIteratorPtr aItE;
3807 SMDS_NodeIteratorPtr aItN;
3808 SMDSAbs_ElementType aTypeE;
3810 TNodeOfNodeListMap mapNewNodes;
3811 TElemOfVecOfNnlmiMap mapElemNewNodes;
3812 TElemOfElemListMap newElemsMap;
3815 aTolVec2=aTolVec*aTolVec;
3818 aNbE = theElements.size();
3821 return EXTR_NO_ELEMENTS;
3823 // 1.1 Track Pattern
3826 SMESHDS_SubMesh* pSubMeshDS=theTrack->GetSubMeshDS();
3828 aItE = pSubMeshDS->GetElements();
3829 while ( aItE->more() ) {
3830 const SMDS_MeshElement* pE = aItE->next();
3831 aTypeE = pE->GetType();
3832 // Pattern must contain links only
3833 if ( aTypeE != SMDSAbs_Edge )
3834 return EXTR_PATH_NOT_EDGE;
3837 const TopoDS_Shape& aS = theTrack->GetSubShape();
3838 // Sub shape for the Pattern must be an Edge
3839 if ( aS.ShapeType() != TopAbs_EDGE )
3840 return EXTR_BAD_PATH_SHAPE;
3842 aTrackEdge = TopoDS::Edge( aS );
3843 // the Edge must not be degenerated
3844 if ( BRep_Tool::Degenerated( aTrackEdge ) )
3845 return EXTR_BAD_PATH_SHAPE;
3847 TopExp::Vertices( aTrackEdge, aV1, aV2 );
3848 aT1=BRep_Tool::Parameter( aV1, aTrackEdge );
3849 aT2=BRep_Tool::Parameter( aV2, aTrackEdge );
3851 aItN = theTrack->GetFather()->GetSubMesh( aV1 )->GetSubMeshDS()->GetNodes();
3852 const SMDS_MeshNode* aN1 = aItN->next();
3854 aItN = theTrack->GetFather()->GetSubMesh( aV2 )->GetSubMeshDS()->GetNodes();
3855 const SMDS_MeshNode* aN2 = aItN->next();
3857 // starting node must be aN1 or aN2
3858 if ( !( aN1 == theN1 || aN2 == theN1 ) )
3859 return EXTR_BAD_STARTING_NODE;
3861 aNbTP = pSubMeshDS->NbNodes() + 2;
3864 vector<double> aAngles( aNbTP );
3866 for ( j=0; j < aNbTP; ++j ) {
3870 if ( theHasAngles ) {
3871 aItD = theAngles.begin();
3872 for ( j=1; (aItD != theAngles.end()) && (j<aNbTP); ++aItD, ++j ) {
3874 aAngles[j] = aAngle;
3878 // 2. Collect parameters on the track edge
3879 aPrms.push_back( aT1 );
3880 aPrms.push_back( aT2 );
3882 aItN = pSubMeshDS->GetNodes();
3883 while ( aItN->more() ) {
3884 const SMDS_MeshNode* pNode = aItN->next();
3885 const SMDS_EdgePosition* pEPos =
3886 static_cast<const SMDS_EdgePosition*>( pNode->GetPosition().get() );
3887 aT = pEPos->GetUParameter();
3888 aPrms.push_back( aT );
3893 if ( aN1 == theN1 ) {
3905 SMESH_MeshEditor_PathPoint aPP;
3906 vector<SMESH_MeshEditor_PathPoint> aPPs( aNbTP );
3908 aC3D = BRep_Tool::Curve( aTrackEdge, aTx1, aTx2 );
3910 aItD = aPrms.begin();
3911 for ( j=0; aItD != aPrms.end(); ++aItD, ++j ) {
3913 aC3D->D1( aT, aP3D, aVec );
3914 aL2 = aVec.SquareMagnitude();
3915 if ( aL2 < aTolVec2 )
3916 return EXTR_CANT_GET_TANGENT;
3918 gp_Dir aTgt( aVec );
3919 aAngle = aAngles[j];
3922 aPP.SetTangent( aTgt );
3923 aPP.SetAngle( aAngle );
3924 aPP.SetParameter( aT );
3928 // 3. Center of rotation aV0
3930 if ( !theHasRefPoint ) {
3932 aGC.SetCoord( 0.,0.,0. );
3934 itElem = theElements.begin();
3935 for ( ; itElem != theElements.end(); itElem++ ) {
3936 const SMDS_MeshElement* elem = *itElem;
3938 SMDS_ElemIteratorPtr itN = elem->nodesIterator();
3939 while ( itN->more() ) {
3940 const SMDS_MeshNode* node = static_cast<const SMDS_MeshNode*>( itN->next() );
3945 if ( mapNewNodes.find( node ) == mapNewNodes.end() ) {
3946 list<const SMDS_MeshNode*> aLNx;
3947 mapNewNodes[node] = aLNx;
3949 gp_XYZ aXYZ( aX, aY, aZ );
3957 } // if (!theHasRefPoint) {
3958 mapNewNodes.clear();
3960 // 4. Processing the elements
3961 SMESHDS_Mesh* aMesh = GetMeshDS();
3963 for ( itElem = theElements.begin(); itElem != theElements.end(); itElem++ ) {
3964 // check element type
3965 const SMDS_MeshElement* elem = *itElem;
3966 aTypeE = elem->GetType();
3967 if ( !elem || ( aTypeE != SMDSAbs_Face && aTypeE != SMDSAbs_Edge ) )
3970 vector<TNodeOfNodeListMapItr> & newNodesItVec = mapElemNewNodes[ elem ];
3971 newNodesItVec.reserve( elem->NbNodes() );
3973 // loop on elem nodes
3975 SMDS_ElemIteratorPtr itN = elem->nodesIterator();
3976 while ( itN->more() )
3979 // check if a node has been already processed
3980 const SMDS_MeshNode* node =
3981 static_cast<const SMDS_MeshNode*>( itN->next() );
3982 TNodeOfNodeListMap::iterator nIt = mapNewNodes.find( node );
3983 if ( nIt == mapNewNodes.end() ) {
3984 nIt = mapNewNodes.insert( make_pair( node, list<const SMDS_MeshNode*>() )).first;
3985 list<const SMDS_MeshNode*>& listNewNodes = nIt->second;
3988 aX = node->X(); aY = node->Y(); aZ = node->Z();
3990 Standard_Real aAngle1x, aAngleT1T0, aTolAng;
3991 gp_Pnt aP0x, aP1x, aPN0, aPN1, aV0x, aV1x;
3992 gp_Ax1 anAx1, anAxT1T0;
3993 gp_Dir aDT1x, aDT0x, aDT1T0;
3998 aPN0.SetCoord(aX, aY, aZ);
4000 const SMESH_MeshEditor_PathPoint& aPP0 = aPPs[0];
4002 aDT0x= aPP0.Tangent();
4004 for ( j = 1; j < aNbTP; ++j ) {
4005 const SMESH_MeshEditor_PathPoint& aPP1 = aPPs[j];
4007 aDT1x = aPP1.Tangent();
4008 aAngle1x = aPP1.Angle();
4010 gp_Trsf aTrsf, aTrsfRot, aTrsfRotT1T0;
4012 gp_Vec aV01x( aP0x, aP1x );
4013 aTrsf.SetTranslation( aV01x );
4016 aV1x = aV0x.Transformed( aTrsf );
4017 aPN1 = aPN0.Transformed( aTrsf );
4019 // rotation 1 [ T1,T0 ]
4020 aAngleT1T0=-aDT1x.Angle( aDT0x );
4021 if (fabs(aAngleT1T0) > aTolAng) {
4023 anAxT1T0.SetLocation( aV1x );
4024 anAxT1T0.SetDirection( aDT1T0 );
4025 aTrsfRotT1T0.SetRotation( anAxT1T0, aAngleT1T0 );
4027 aPN1 = aPN1.Transformed( aTrsfRotT1T0 );
4031 if ( theHasAngles ) {
4032 anAx1.SetLocation( aV1x );
4033 anAx1.SetDirection( aDT1x );
4034 aTrsfRot.SetRotation( anAx1, aAngle1x );
4036 aPN1 = aPN1.Transformed( aTrsfRot );
4040 if( elem->IsQuadratic() && !elem->IsMediumNode(node) ) {
4041 // create additional node
4042 double x = ( aPN1.X() + aPN0.X() )/2.;
4043 double y = ( aPN1.Y() + aPN0.Y() )/2.;
4044 double z = ( aPN1.Z() + aPN0.Z() )/2.;
4045 const SMDS_MeshNode* newNode = aMesh->AddNode(x,y,z);
4046 myLastCreatedNodes.Append(newNode);
4047 srcNodes.Append( node );
4048 listNewNodes.push_back( newNode );
4053 const SMDS_MeshNode* newNode = aMesh->AddNode( aX, aY, aZ );
4054 myLastCreatedNodes.Append(newNode);
4055 srcNodes.Append( node );
4056 listNewNodes.push_back( newNode );
4066 // if current elem is quadratic and current node is not medium
4067 // we have to check - may be it is needed to insert additional nodes
4068 if( elem->IsQuadratic() && !elem->IsMediumNode(node) ) {
4069 list< const SMDS_MeshNode* > & listNewNodes = nIt->second;
4070 if(listNewNodes.size()==aNbTP-1) {
4071 vector<const SMDS_MeshNode*> aNodes(2*(aNbTP-1));
4072 gp_XYZ P(node->X(), node->Y(), node->Z());
4073 list< const SMDS_MeshNode* >::iterator it = listNewNodes.begin();
4075 for(i=0; i<aNbTP-1; i++) {
4076 const SMDS_MeshNode* N = *it;
4077 double x = ( N->X() + P.X() )/2.;
4078 double y = ( N->Y() + P.Y() )/2.;
4079 double z = ( N->Z() + P.Z() )/2.;
4080 const SMDS_MeshNode* newN = aMesh->AddNode(x,y,z);
4081 srcNodes.Append( node );
4082 myLastCreatedNodes.Append(newN);
4085 P = gp_XYZ(N->X(),N->Y(),N->Z());
4087 listNewNodes.clear();
4088 for(i=0; i<2*(aNbTP-1); i++) {
4089 listNewNodes.push_back(aNodes[i]);
4095 newNodesItVec.push_back( nIt );
4097 // make new elements
4098 //sweepElement( aMesh, elem, newNodesItVec, newElemsMap[elem],
4099 // newNodesItVec[0]->second.size(), myLastCreatedElems );
4100 sweepElement( elem, newNodesItVec, newElemsMap[elem], aNbTP-1, srcElems );
4103 makeWalls( mapNewNodes, newElemsMap, mapElemNewNodes, theElements, aNbTP-1, srcElems );
4105 if ( theMakeGroups )
4106 generateGroups( srcNodes, srcElems, "extruded");
4111 //=======================================================================
4112 //function : Transform
4114 //=======================================================================
4116 SMESH_MeshEditor::PGroupIDs
4117 SMESH_MeshEditor::Transform (TIDSortedElemSet & theElems,
4118 const gp_Trsf& theTrsf,
4120 const bool theMakeGroups,
4121 SMESH_Mesh* theTargetMesh)
4123 myLastCreatedElems.Clear();
4124 myLastCreatedNodes.Clear();
4126 bool needReverse = false;
4127 string groupPostfix;
4128 switch ( theTrsf.Form() ) {
4133 groupPostfix = "mirrored";
4136 groupPostfix = "rotated";
4138 case gp_Translation:
4139 groupPostfix = "translated";
4142 groupPostfix = "scaled";
4145 needReverse = false;
4146 groupPostfix = "transformed";
4149 SMESH_MeshEditor targetMeshEditor( theTargetMesh );
4150 SMESHDS_Mesh* aTgtMesh = theTargetMesh ? theTargetMesh->GetMeshDS() : 0;
4151 SMESHDS_Mesh* aMesh = GetMeshDS();
4154 // map old node to new one
4155 TNodeNodeMap nodeMap;
4157 // elements sharing moved nodes; those of them which have all
4158 // nodes mirrored but are not in theElems are to be reversed
4159 TIDSortedElemSet inverseElemSet;
4161 // source elements for each generated one
4162 SMESH_SequenceOfElemPtr srcElems, srcNodes;
4165 TIDSortedElemSet::iterator itElem;
4166 for ( itElem = theElems.begin(); itElem != theElems.end(); itElem++ ) {
4167 const SMDS_MeshElement* elem = *itElem;
4171 // loop on elem nodes
4172 SMDS_ElemIteratorPtr itN = elem->nodesIterator();
4173 while ( itN->more() ) {
4175 // check if a node has been already transformed
4176 const SMDS_MeshNode* node = cast2Node( itN->next() );
4177 pair<TNodeNodeMap::iterator,bool> n2n_isnew =
4178 nodeMap.insert( make_pair ( node, node ));
4179 if ( !n2n_isnew.second )
4183 coord[0] = node->X();
4184 coord[1] = node->Y();
4185 coord[2] = node->Z();
4186 theTrsf.Transforms( coord[0], coord[1], coord[2] );
4187 if ( theTargetMesh ) {
4188 const SMDS_MeshNode * newNode = aTgtMesh->AddNode( coord[0], coord[1], coord[2] );
4189 n2n_isnew.first->second = newNode;
4190 myLastCreatedNodes.Append(newNode);
4191 srcNodes.Append( node );
4193 else if ( theCopy ) {
4194 const SMDS_MeshNode * newNode = aMesh->AddNode( coord[0], coord[1], coord[2] );
4195 n2n_isnew.first->second = newNode;
4196 myLastCreatedNodes.Append(newNode);
4197 srcNodes.Append( node );
4200 aMesh->MoveNode( node, coord[0], coord[1], coord[2] );
4201 // node position on shape becomes invalid
4202 const_cast< SMDS_MeshNode* > ( node )->SetPosition
4203 ( SMDS_SpacePosition::originSpacePosition() );
4206 // keep inverse elements
4207 if ( !theCopy && !theTargetMesh && needReverse ) {
4208 SMDS_ElemIteratorPtr invElemIt = node->GetInverseElementIterator();
4209 while ( invElemIt->more() ) {
4210 const SMDS_MeshElement* iel = invElemIt->next();
4211 inverseElemSet.insert( iel );
4217 // either create new elements or reverse mirrored ones
4218 if ( !theCopy && !needReverse && !theTargetMesh )
4221 TIDSortedElemSet::iterator invElemIt = inverseElemSet.begin();
4222 for ( ; invElemIt != inverseElemSet.end(); invElemIt++ )
4223 theElems.insert( *invElemIt );
4225 // replicate or reverse elements
4228 REV_TETRA = 0, // = nbNodes - 4
4229 REV_PYRAMID = 1, // = nbNodes - 4
4230 REV_PENTA = 2, // = nbNodes - 4
4232 REV_HEXA = 4, // = nbNodes - 4
4236 { 2, 1, 0, 3, 4, 0, 0, 0 }, // REV_TETRA
4237 { 2, 1, 0, 3, 4, 0, 0, 0 }, // REV_PYRAMID
4238 { 2, 1, 0, 5, 4, 3, 0, 0 }, // REV_PENTA
4239 { 2, 1, 0, 3, 0, 0, 0, 0 }, // REV_FACE
4240 { 2, 1, 0, 3, 6, 5, 4, 7 }, // REV_HEXA
4241 { 0, 1, 2, 3, 4, 5, 6, 7 } // FORWARD
4244 for ( itElem = theElems.begin(); itElem != theElems.end(); itElem++ )
4246 const SMDS_MeshElement* elem = *itElem;
4247 if ( !elem || elem->GetType() == SMDSAbs_Node )
4250 int nbNodes = elem->NbNodes();
4251 int elemType = elem->GetType();
4253 if (elem->IsPoly()) {
4254 // Polygon or Polyhedral Volume
4255 switch ( elemType ) {
4258 vector<const SMDS_MeshNode*> poly_nodes (nbNodes);
4260 SMDS_ElemIteratorPtr itN = elem->nodesIterator();
4261 while (itN->more()) {
4262 const SMDS_MeshNode* node =
4263 static_cast<const SMDS_MeshNode*>(itN->next());
4264 TNodeNodeMap::iterator nodeMapIt = nodeMap.find(node);
4265 if (nodeMapIt == nodeMap.end())
4266 break; // not all nodes transformed
4268 // reverse mirrored faces and volumes
4269 poly_nodes[nbNodes - iNode - 1] = (*nodeMapIt).second;
4271 poly_nodes[iNode] = (*nodeMapIt).second;
4275 if ( iNode != nbNodes )
4276 continue; // not all nodes transformed
4278 if ( theTargetMesh ) {
4279 myLastCreatedElems.Append(aTgtMesh->AddPolygonalFace(poly_nodes));
4280 srcElems.Append( elem );
4282 else if ( theCopy ) {
4283 myLastCreatedElems.Append(aMesh->AddPolygonalFace(poly_nodes));
4284 srcElems.Append( elem );
4287 aMesh->ChangePolygonNodes(elem, poly_nodes);
4291 case SMDSAbs_Volume:
4293 // ATTENTION: Reversing is not yet done!!!
4294 const SMDS_PolyhedralVolumeOfNodes* aPolyedre =
4295 dynamic_cast<const SMDS_PolyhedralVolumeOfNodes*>( elem );
4297 MESSAGE("Warning: bad volumic element");
4301 vector<const SMDS_MeshNode*> poly_nodes;
4302 vector<int> quantities;
4304 bool allTransformed = true;
4305 int nbFaces = aPolyedre->NbFaces();
4306 for (int iface = 1; iface <= nbFaces && allTransformed; iface++) {
4307 int nbFaceNodes = aPolyedre->NbFaceNodes(iface);
4308 for (int inode = 1; inode <= nbFaceNodes && allTransformed; inode++) {
4309 const SMDS_MeshNode* node = aPolyedre->GetFaceNode(iface, inode);
4310 TNodeNodeMap::iterator nodeMapIt = nodeMap.find(node);
4311 if (nodeMapIt == nodeMap.end()) {
4312 allTransformed = false; // not all nodes transformed
4314 poly_nodes.push_back((*nodeMapIt).second);
4317 quantities.push_back(nbFaceNodes);
4319 if ( !allTransformed )
4320 continue; // not all nodes transformed
4322 if ( theTargetMesh ) {
4323 myLastCreatedElems.Append(aTgtMesh->AddPolyhedralVolume(poly_nodes, quantities));
4324 srcElems.Append( elem );
4326 else if ( theCopy ) {
4327 myLastCreatedElems.Append(aMesh->AddPolyhedralVolume(poly_nodes, quantities));
4328 srcElems.Append( elem );
4331 aMesh->ChangePolyhedronNodes(elem, poly_nodes, quantities);
4341 int* i = index[ FORWARD ];
4342 if ( needReverse && nbNodes > 2) // reverse mirrored faces and volumes
4343 if ( elemType == SMDSAbs_Face )
4344 i = index[ REV_FACE ];
4346 i = index[ nbNodes - 4 ];
4348 if(elem->IsQuadratic()) {
4349 static int anIds[] = {0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19};
4352 if(nbNodes==3) { // quadratic edge
4353 static int anIds[] = {1,0,2};
4356 else if(nbNodes==6) { // quadratic triangle
4357 static int anIds[] = {0,2,1,5,4,3};
4360 else if(nbNodes==8) { // quadratic quadrangle
4361 static int anIds[] = {0,3,2,1,7,6,5,4};
4364 else if(nbNodes==10) { // quadratic tetrahedron of 10 nodes
4365 static int anIds[] = {0,2,1,3,6,5,4,7,9,8};
4368 else if(nbNodes==13) { // quadratic pyramid of 13 nodes
4369 static int anIds[] = {0,3,2,1,4,8,7,6,5,9,12,11,10};
4372 else if(nbNodes==15) { // quadratic pentahedron with 15 nodes
4373 static int anIds[] = {0,2,1,3,5,4,8,7,6,11,10,9,12,14,13};
4376 else { // nbNodes==20 - quadratic hexahedron with 20 nodes
4377 static int anIds[] = {0,3,2,1,4,7,6,5,11,10,9,8,15,14,13,12,16,19,18,17};
4383 // find transformed nodes
4384 vector<const SMDS_MeshNode*> nodes(nbNodes);
4386 SMDS_ElemIteratorPtr itN = elem->nodesIterator();
4387 while ( itN->more() ) {
4388 const SMDS_MeshNode* node =
4389 static_cast<const SMDS_MeshNode*>( itN->next() );
4390 TNodeNodeMap::iterator nodeMapIt = nodeMap.find( node );
4391 if ( nodeMapIt == nodeMap.end() )
4392 break; // not all nodes transformed
4393 nodes[ i [ iNode++ ]] = (*nodeMapIt).second;
4395 if ( iNode != nbNodes )
4396 continue; // not all nodes transformed
4398 if ( theTargetMesh ) {
4399 if ( SMDS_MeshElement* copy =
4400 targetMeshEditor.AddElement( nodes, elem->GetType(), elem->IsPoly() )) {
4401 myLastCreatedElems.Append( copy );
4402 srcElems.Append( elem );
4405 else if ( theCopy ) {
4406 if ( SMDS_MeshElement* copy = AddElement( nodes, elem->GetType(), elem->IsPoly() )) {
4407 myLastCreatedElems.Append( copy );
4408 srcElems.Append( elem );
4412 // reverse element as it was reversed by transformation
4414 aMesh->ChangeElementNodes( elem, &nodes[0], nbNodes );
4418 PGroupIDs newGroupIDs;
4420 if ( theMakeGroups && theCopy ||
4421 theMakeGroups && theTargetMesh )
4422 newGroupIDs = generateGroups( srcNodes, srcElems, groupPostfix, theTargetMesh );
4427 //=======================================================================
4429 * \brief Create groups of elements made during transformation
4430 * \param nodeGens - nodes making corresponding myLastCreatedNodes
4431 * \param elemGens - elements making corresponding myLastCreatedElems
4432 * \param postfix - to append to names of new groups
4434 //=======================================================================
4436 SMESH_MeshEditor::PGroupIDs
4437 SMESH_MeshEditor::generateGroups(const SMESH_SequenceOfElemPtr& nodeGens,
4438 const SMESH_SequenceOfElemPtr& elemGens,
4439 const std::string& postfix,
4440 SMESH_Mesh* targetMesh)
4442 PGroupIDs newGroupIDs( new list<int> );
4443 SMESH_Mesh* mesh = targetMesh ? targetMesh : GetMesh();
4445 // Sort existing groups by types and collect their names
4447 // to store an old group and a generated new one
4448 typedef pair< SMESHDS_GroupBase*, SMDS_MeshGroup* > TOldNewGroup;
4449 vector< list< TOldNewGroup > > groupsByType( SMDSAbs_NbElementTypes );
4451 set< string > groupNames;
4453 SMDS_MeshGroup* nullNewGroup = (SMDS_MeshGroup*) 0;
4454 SMESH_Mesh::GroupIteratorPtr groupIt = GetMesh()->GetGroups();
4455 while ( groupIt->more() ) {
4456 SMESH_Group * group = groupIt->next();
4457 if ( !group ) continue;
4458 SMESHDS_GroupBase* groupDS = group->GetGroupDS();
4459 if ( !groupDS || groupDS->IsEmpty() ) continue;
4460 groupNames.insert( group->GetName() );
4461 groupDS->SetStoreName( group->GetName() );
4462 groupsByType[ groupDS->GetType() ].push_back( make_pair( groupDS, nullNewGroup ));
4467 // loop on nodes and elements
4468 for ( int isNodes = 0; isNodes < 2; ++isNodes )
4470 const SMESH_SequenceOfElemPtr& gens = isNodes ? nodeGens : elemGens;
4471 const SMESH_SequenceOfElemPtr& elems = isNodes ? myLastCreatedNodes : myLastCreatedElems;
4472 if ( gens.Length() != elems.Length() )
4473 throw SALOME_Exception(LOCALIZED("invalid args"));
4475 // loop on created elements
4476 for (int iElem = 1; iElem <= elems.Length(); ++iElem )
4478 const SMDS_MeshElement* sourceElem = gens( iElem );
4479 if ( !sourceElem ) {
4480 MESSAGE("generateGroups(): NULL source element");
4483 list< TOldNewGroup > & groupsOldNew = groupsByType[ sourceElem->GetType() ];
4484 if ( groupsOldNew.empty() ) {
4485 while ( iElem < gens.Length() && gens( iElem+1 ) == sourceElem )
4486 ++iElem; // skip all elements made by sourceElem
4489 // collect all elements made by sourceElem
4490 list< const SMDS_MeshElement* > resultElems;
4491 if ( const SMDS_MeshElement* resElem = elems( iElem ))
4492 if ( resElem != sourceElem )
4493 resultElems.push_back( resElem );
4494 while ( iElem < gens.Length() && gens( iElem+1 ) == sourceElem )
4495 if ( const SMDS_MeshElement* resElem = elems( ++iElem ))
4496 if ( resElem != sourceElem )
4497 resultElems.push_back( resElem );
4498 // do not generate element groups from node ones
4499 if ( sourceElem->GetType() == SMDSAbs_Node &&
4500 elems( iElem )->GetType() != SMDSAbs_Node )
4503 // add resultElems to groups made by ones the sourceElem belongs to
4504 list< TOldNewGroup >::iterator gOldNew, gLast = groupsOldNew.end();
4505 for ( gOldNew = groupsOldNew.begin(); gOldNew != gLast; ++gOldNew )
4507 SMESHDS_GroupBase* oldGroup = gOldNew->first;
4508 if ( oldGroup->Contains( sourceElem )) // sourceElem in oldGroup
4510 SMDS_MeshGroup* & newGroup = gOldNew->second;
4511 if ( !newGroup )// create a new group
4514 string name = oldGroup->GetStoreName();
4515 if ( !targetMesh ) {
4519 while ( !groupNames.insert( name ).second ) // name exists
4525 TCollection_AsciiString nbStr(nb+1);
4526 name.resize( name.rfind('_')+1 );
4527 name += nbStr.ToCString();
4534 SMESH_Group* group = mesh->AddGroup( resultElems.back()->GetType(),
4536 SMESHDS_Group* groupDS = static_cast<SMESHDS_Group*>(group->GetGroupDS());
4537 newGroup = & groupDS->SMDSGroup();
4538 newGroupIDs->push_back( id );
4541 // fill in a new group
4542 list< const SMDS_MeshElement* >::iterator resLast = resultElems.end(), resElemIt;
4543 for ( resElemIt = resultElems.begin(); resElemIt != resLast; ++resElemIt )
4544 newGroup->Add( *resElemIt );
4547 } // loop on created elements
4548 }// loop on nodes and elements
4553 //=======================================================================
4554 //function : FindCoincidentNodes
4555 //purpose : Return list of group of nodes close to each other within theTolerance
4556 // Search among theNodes or in the whole mesh if theNodes is empty using
4557 // an Octree algorithm
4558 //=======================================================================
4560 void SMESH_MeshEditor::FindCoincidentNodes (set<const SMDS_MeshNode*> & theNodes,
4561 const double theTolerance,
4562 TListOfListOfNodes & theGroupsOfNodes)
4564 myLastCreatedElems.Clear();
4565 myLastCreatedNodes.Clear();
4567 set<const SMDS_MeshNode*> nodes;
4568 if ( theNodes.empty() )
4569 { // get all nodes in the mesh
4570 SMDS_NodeIteratorPtr nIt = GetMeshDS()->nodesIterator();
4571 while ( nIt->more() )
4572 nodes.insert( nodes.end(),nIt->next());
4576 SMESH_OctreeNode::FindCoincidentNodes ( nodes, &theGroupsOfNodes, theTolerance);
4580 //=======================================================================
4582 * \brief Implementation of search for the node closest to point
4584 //=======================================================================
4586 struct SMESH_NodeSearcherImpl: public SMESH_NodeSearcher
4589 * \brief Constructor
4591 SMESH_NodeSearcherImpl( const SMESHDS_Mesh* theMesh )
4593 set<const SMDS_MeshNode*> nodes;
4595 SMDS_NodeIteratorPtr nIt = theMesh->nodesIterator();
4596 while ( nIt->more() )
4597 nodes.insert( nodes.end(), nIt->next() );
4599 myOctreeNode = new SMESH_OctreeNode(nodes) ;
4602 * \brief Do it's job
4604 const SMDS_MeshNode* FindClosestTo( const gp_Pnt& thePnt )
4606 SMDS_MeshNode tgtNode( thePnt.X(), thePnt.Y(), thePnt.Z() );
4607 list<const SMDS_MeshNode*> nodes;
4608 const double precision = 1e-6;
4609 myOctreeNode->NodesAround( &tgtNode, &nodes, precision );
4611 double minSqDist = DBL_MAX;
4613 if ( nodes.empty() ) // get all nodes of OctreeNode's closest to thePnt
4615 // sort leafs by their distance from thePnt
4616 typedef map< double, SMESH_OctreeNode* > TDistTreeMap;
4617 TDistTreeMap treeMap;
4618 list< SMESH_OctreeNode* > treeList;
4619 list< SMESH_OctreeNode* >::iterator trIt;
4620 treeList.push_back( myOctreeNode );
4621 for ( trIt = treeList.begin(); trIt != treeList.end(); ++trIt)
4623 SMESH_OctreeNode* tree = *trIt;
4624 if ( !tree->isLeaf() ) { // put children to the queue
4625 SMESH_OctreeNodeIteratorPtr cIt = tree->GetChildrenIterator();
4626 while ( cIt->more() )
4627 treeList.push_back( cIt->next() );
4629 else if ( tree->NbNodes() ) { // put tree to treeMap
4630 tree->getBox( box );
4631 double sqDist = thePnt.SquareDistance( 0.5 * ( box.CornerMin() + box.CornerMax() ));
4632 pair<TDistTreeMap::iterator,bool> it_in = treeMap.insert( make_pair( sqDist, tree ));
4633 if ( !it_in.second ) // not unique distance to box center
4634 treeMap.insert( it_in.first, make_pair( sqDist - 1e-13*treeMap.size(), tree ));
4637 // find distance after which there is no sense to check tree's
4638 double sqLimit = DBL_MAX;
4639 TDistTreeMap::iterator sqDist_tree = treeMap.begin();
4640 if ( treeMap.size() > 5 ) {
4641 SMESH_OctreeNode* closestTree = sqDist_tree->second;
4642 closestTree->getBox( box );
4643 double limit = sqrt( sqDist_tree->first ) + sqrt ( box.SquareExtent() );
4644 sqLimit = limit * limit;
4646 // get all nodes from trees
4647 for ( ; sqDist_tree != treeMap.end(); ++sqDist_tree) {
4648 if ( sqDist_tree->first > sqLimit )
4650 SMESH_OctreeNode* tree = sqDist_tree->second;
4651 tree->NodesAround( tree->GetNodeIterator()->next(), &nodes );
4654 // find closest among nodes
4655 minSqDist = DBL_MAX;
4656 const SMDS_MeshNode* closestNode = 0;
4657 list<const SMDS_MeshNode*>::iterator nIt = nodes.begin();
4658 for ( ; nIt != nodes.end(); ++nIt ) {
4659 double sqDist = thePnt.SquareDistance( TNodeXYZ( *nIt ) );
4660 if ( minSqDist > sqDist ) {
4670 ~SMESH_NodeSearcherImpl() { delete myOctreeNode; }
4672 SMESH_OctreeNode* myOctreeNode;
4675 //=======================================================================
4677 * \brief Return SMESH_NodeSearcher
4679 //=======================================================================
4681 SMESH_NodeSearcher* SMESH_MeshEditor::GetNodeSearcher()
4683 return new SMESH_NodeSearcherImpl( GetMeshDS() );
4686 //=======================================================================
4687 //function : SimplifyFace
4689 //=======================================================================
4690 int SMESH_MeshEditor::SimplifyFace (const vector<const SMDS_MeshNode *> faceNodes,
4691 vector<const SMDS_MeshNode *>& poly_nodes,
4692 vector<int>& quantities) const
4694 int nbNodes = faceNodes.size();
4699 set<const SMDS_MeshNode*> nodeSet;
4701 // get simple seq of nodes
4702 //const SMDS_MeshNode* simpleNodes[ nbNodes ];
4703 vector<const SMDS_MeshNode*> simpleNodes( nbNodes );
4704 int iSimple = 0, nbUnique = 0;
4706 simpleNodes[iSimple++] = faceNodes[0];
4708 for (int iCur = 1; iCur < nbNodes; iCur++) {
4709 if (faceNodes[iCur] != simpleNodes[iSimple - 1]) {
4710 simpleNodes[iSimple++] = faceNodes[iCur];
4711 if (nodeSet.insert( faceNodes[iCur] ).second)
4715 int nbSimple = iSimple;
4716 if (simpleNodes[nbSimple - 1] == simpleNodes[0]) {
4726 bool foundLoop = (nbSimple > nbUnique);
4729 set<const SMDS_MeshNode*> loopSet;
4730 for (iSimple = 0; iSimple < nbSimple && !foundLoop; iSimple++) {
4731 const SMDS_MeshNode* n = simpleNodes[iSimple];
4732 if (!loopSet.insert( n ).second) {
4736 int iC = 0, curLast = iSimple;
4737 for (; iC < curLast; iC++) {
4738 if (simpleNodes[iC] == n) break;
4740 int loopLen = curLast - iC;
4742 // create sub-element
4744 quantities.push_back(loopLen);
4745 for (; iC < curLast; iC++) {
4746 poly_nodes.push_back(simpleNodes[iC]);
4749 // shift the rest nodes (place from the first loop position)
4750 for (iC = curLast + 1; iC < nbSimple; iC++) {
4751 simpleNodes[iC - loopLen] = simpleNodes[iC];
4753 nbSimple -= loopLen;
4756 } // for (iSimple = 0; iSimple < nbSimple; iSimple++)
4757 } // while (foundLoop)
4761 quantities.push_back(iSimple);
4762 for (int i = 0; i < iSimple; i++)
4763 poly_nodes.push_back(simpleNodes[i]);
4769 //=======================================================================
4770 //function : MergeNodes
4771 //purpose : In each group, the cdr of nodes are substituted by the first one
4773 //=======================================================================
4775 void SMESH_MeshEditor::MergeNodes (TListOfListOfNodes & theGroupsOfNodes)
4777 myLastCreatedElems.Clear();
4778 myLastCreatedNodes.Clear();
4780 SMESHDS_Mesh* aMesh = GetMeshDS();
4782 TNodeNodeMap nodeNodeMap; // node to replace - new node
4783 set<const SMDS_MeshElement*> elems; // all elements with changed nodes
4784 list< int > rmElemIds, rmNodeIds;
4786 // Fill nodeNodeMap and elems
4788 TListOfListOfNodes::iterator grIt = theGroupsOfNodes.begin();
4789 for ( ; grIt != theGroupsOfNodes.end(); grIt++ ) {
4790 list<const SMDS_MeshNode*>& nodes = *grIt;
4791 list<const SMDS_MeshNode*>::iterator nIt = nodes.begin();
4792 const SMDS_MeshNode* nToKeep = *nIt;
4793 for ( ++nIt; nIt != nodes.end(); nIt++ ) {
4794 const SMDS_MeshNode* nToRemove = *nIt;
4795 nodeNodeMap.insert( TNodeNodeMap::value_type( nToRemove, nToKeep ));
4796 if ( nToRemove != nToKeep ) {
4797 rmNodeIds.push_back( nToRemove->GetID() );
4798 AddToSameGroups( nToKeep, nToRemove, aMesh );
4801 SMDS_ElemIteratorPtr invElemIt = nToRemove->GetInverseElementIterator();
4802 while ( invElemIt->more() ) {
4803 const SMDS_MeshElement* elem = invElemIt->next();
4808 // Change element nodes or remove an element
4810 set<const SMDS_MeshElement*>::iterator eIt = elems.begin();
4811 for ( ; eIt != elems.end(); eIt++ ) {
4812 const SMDS_MeshElement* elem = *eIt;
4813 int nbNodes = elem->NbNodes();
4814 int aShapeId = FindShape( elem );
4816 set<const SMDS_MeshNode*> nodeSet;
4817 vector< const SMDS_MeshNode*> curNodes( nbNodes ), uniqueNodes( nbNodes );
4818 int iUnique = 0, iCur = 0, nbRepl = 0;
4819 vector<int> iRepl( nbNodes );
4821 // get new seq of nodes
4822 SMDS_ElemIteratorPtr itN = elem->nodesIterator();
4823 while ( itN->more() ) {
4824 const SMDS_MeshNode* n =
4825 static_cast<const SMDS_MeshNode*>( itN->next() );
4827 TNodeNodeMap::iterator nnIt = nodeNodeMap.find( n );
4828 if ( nnIt != nodeNodeMap.end() ) { // n sticks
4830 iRepl[ nbRepl++ ] = iCur;
4832 curNodes[ iCur ] = n;
4833 bool isUnique = nodeSet.insert( n ).second;
4835 uniqueNodes[ iUnique++ ] = n;
4839 // Analyse element topology after replacement
4842 int nbUniqueNodes = nodeSet.size();
4843 if ( nbNodes != nbUniqueNodes ) { // some nodes stick
4844 // Polygons and Polyhedral volumes
4845 if (elem->IsPoly()) {
4847 if (elem->GetType() == SMDSAbs_Face) {
4849 vector<const SMDS_MeshNode *> face_nodes (nbNodes);
4851 for (; inode < nbNodes; inode++) {
4852 face_nodes[inode] = curNodes[inode];
4855 vector<const SMDS_MeshNode *> polygons_nodes;
4856 vector<int> quantities;
4857 int nbNew = SimplifyFace(face_nodes, polygons_nodes, quantities);
4861 for (int iface = 0; iface < nbNew - 1; iface++) {
4862 int nbNodes = quantities[iface];
4863 vector<const SMDS_MeshNode *> poly_nodes (nbNodes);
4864 for (int ii = 0; ii < nbNodes; ii++, inode++) {
4865 poly_nodes[ii] = polygons_nodes[inode];
4867 SMDS_MeshElement* newElem = aMesh->AddPolygonalFace(poly_nodes);
4868 myLastCreatedElems.Append(newElem);
4870 aMesh->SetMeshElementOnShape(newElem, aShapeId);
4872 aMesh->ChangeElementNodes(elem, &polygons_nodes[inode], quantities[nbNew - 1]);
4875 rmElemIds.push_back(elem->GetID());
4879 else if (elem->GetType() == SMDSAbs_Volume) {
4880 // Polyhedral volume
4881 if (nbUniqueNodes < 4) {
4882 rmElemIds.push_back(elem->GetID());
4885 // each face has to be analized in order to check volume validity
4886 const SMDS_PolyhedralVolumeOfNodes* aPolyedre =
4887 static_cast<const SMDS_PolyhedralVolumeOfNodes*>( elem );
4889 int nbFaces = aPolyedre->NbFaces();
4891 vector<const SMDS_MeshNode *> poly_nodes;
4892 vector<int> quantities;
4894 for (int iface = 1; iface <= nbFaces; iface++) {
4895 int nbFaceNodes = aPolyedre->NbFaceNodes(iface);
4896 vector<const SMDS_MeshNode *> faceNodes (nbFaceNodes);
4898 for (int inode = 1; inode <= nbFaceNodes; inode++) {
4899 const SMDS_MeshNode * faceNode = aPolyedre->GetFaceNode(iface, inode);
4900 TNodeNodeMap::iterator nnIt = nodeNodeMap.find(faceNode);
4901 if (nnIt != nodeNodeMap.end()) { // faceNode sticks
4902 faceNode = (*nnIt).second;
4904 faceNodes[inode - 1] = faceNode;
4907 SimplifyFace(faceNodes, poly_nodes, quantities);
4910 if (quantities.size() > 3) {
4911 // to be done: remove coincident faces
4914 if (quantities.size() > 3)
4915 aMesh->ChangePolyhedronNodes(elem, poly_nodes, quantities);
4917 rmElemIds.push_back(elem->GetID());
4921 rmElemIds.push_back(elem->GetID());
4932 switch ( nbNodes ) {
4933 case 2: ///////////////////////////////////// EDGE
4934 isOk = false; break;
4935 case 3: ///////////////////////////////////// TRIANGLE
4936 isOk = false; break;
4938 if ( elem->GetType() == SMDSAbs_Volume ) // TETRAHEDRON
4940 else { //////////////////////////////////// QUADRANGLE
4941 if ( nbUniqueNodes < 3 )
4943 else if ( nbRepl == 2 && iRepl[ 1 ] - iRepl[ 0 ] == 2 )
4944 isOk = false; // opposite nodes stick
4947 case 6: ///////////////////////////////////// PENTAHEDRON
4948 if ( nbUniqueNodes == 4 ) {
4949 // ---------------------------------> tetrahedron
4951 iRepl[ 0 ] > 2 && iRepl[ 1 ] > 2 && iRepl[ 2 ] > 2 ) {
4952 // all top nodes stick: reverse a bottom
4953 uniqueNodes[ 0 ] = curNodes [ 1 ];
4954 uniqueNodes[ 1 ] = curNodes [ 0 ];
4956 else if (nbRepl == 3 &&
4957 iRepl[ 0 ] < 3 && iRepl[ 1 ] < 3 && iRepl[ 2 ] < 3 ) {
4958 // all bottom nodes stick: set a top before
4959 uniqueNodes[ 3 ] = uniqueNodes [ 0 ];
4960 uniqueNodes[ 0 ] = curNodes [ 3 ];
4961 uniqueNodes[ 1 ] = curNodes [ 4 ];
4962 uniqueNodes[ 2 ] = curNodes [ 5 ];
4964 else if (nbRepl == 4 &&
4965 iRepl[ 2 ] - iRepl [ 0 ] == 3 && iRepl[ 3 ] - iRepl [ 1 ] == 3 ) {
4966 // a lateral face turns into a line: reverse a bottom
4967 uniqueNodes[ 0 ] = curNodes [ 1 ];
4968 uniqueNodes[ 1 ] = curNodes [ 0 ];
4973 else if ( nbUniqueNodes == 5 ) {
4974 // PENTAHEDRON --------------------> 2 tetrahedrons
4975 if ( nbRepl == 2 && iRepl[ 1 ] - iRepl [ 0 ] == 3 ) {
4976 // a bottom node sticks with a linked top one
4978 SMDS_MeshElement* newElem =
4979 aMesh->AddVolume(curNodes[ 3 ],
4982 curNodes[ iRepl[ 0 ] == 2 ? 1 : 2 ]);
4983 myLastCreatedElems.Append(newElem);
4985 aMesh->SetMeshElementOnShape( newElem, aShapeId );
4986 // 2. : reverse a bottom
4987 uniqueNodes[ 0 ] = curNodes [ 1 ];
4988 uniqueNodes[ 1 ] = curNodes [ 0 ];
4998 if(elem->IsQuadratic()) { // Quadratic quadrangle
5011 if( iRepl[0]==0 && iRepl[1]==1 && iRepl[2]==4 ) {
5012 uniqueNodes[0] = curNodes[0];
5013 uniqueNodes[1] = curNodes[2];
5014 uniqueNodes[2] = curNodes[3];
5015 uniqueNodes[3] = curNodes[5];
5016 uniqueNodes[4] = curNodes[6];
5017 uniqueNodes[5] = curNodes[7];
5020 if( iRepl[0]==0 && iRepl[1]==3 && iRepl[2]==7 ) {
5021 uniqueNodes[0] = curNodes[0];
5022 uniqueNodes[1] = curNodes[1];
5023 uniqueNodes[2] = curNodes[2];
5024 uniqueNodes[3] = curNodes[4];
5025 uniqueNodes[4] = curNodes[5];
5026 uniqueNodes[5] = curNodes[6];
5029 if( iRepl[0]==0 && iRepl[1]==4 && iRepl[2]==7 ) {
5030 uniqueNodes[0] = curNodes[1];
5031 uniqueNodes[1] = curNodes[2];
5032 uniqueNodes[2] = curNodes[3];
5033 uniqueNodes[3] = curNodes[5];
5034 uniqueNodes[4] = curNodes[6];
5035 uniqueNodes[5] = curNodes[0];
5038 if( iRepl[0]==1 && iRepl[1]==2 && iRepl[2]==5 ) {
5039 uniqueNodes[0] = curNodes[0];
5040 uniqueNodes[1] = curNodes[1];
5041 uniqueNodes[2] = curNodes[3];
5042 uniqueNodes[3] = curNodes[4];
5043 uniqueNodes[4] = curNodes[6];
5044 uniqueNodes[5] = curNodes[7];
5047 if( iRepl[0]==1 && iRepl[1]==4 && iRepl[2]==5 ) {
5048 uniqueNodes[0] = curNodes[0];
5049 uniqueNodes[1] = curNodes[2];
5050 uniqueNodes[2] = curNodes[3];
5051 uniqueNodes[3] = curNodes[1];
5052 uniqueNodes[4] = curNodes[6];
5053 uniqueNodes[5] = curNodes[7];
5056 if( iRepl[0]==2 && iRepl[1]==3 && iRepl[2]==6 ) {
5057 uniqueNodes[0] = curNodes[0];
5058 uniqueNodes[1] = curNodes[1];
5059 uniqueNodes[2] = curNodes[2];
5060 uniqueNodes[3] = curNodes[4];
5061 uniqueNodes[4] = curNodes[5];
5062 uniqueNodes[5] = curNodes[7];
5065 if( iRepl[0]==2 && iRepl[1]==5 && iRepl[2]==6 ) {
5066 uniqueNodes[0] = curNodes[0];
5067 uniqueNodes[1] = curNodes[1];
5068 uniqueNodes[2] = curNodes[3];
5069 uniqueNodes[3] = curNodes[4];
5070 uniqueNodes[4] = curNodes[2];
5071 uniqueNodes[5] = curNodes[7];
5074 if( iRepl[0]==3 && iRepl[1]==6 && iRepl[2]==7 ) {
5075 uniqueNodes[0] = curNodes[0];
5076 uniqueNodes[1] = curNodes[1];
5077 uniqueNodes[2] = curNodes[2];
5078 uniqueNodes[3] = curNodes[4];
5079 uniqueNodes[4] = curNodes[5];
5080 uniqueNodes[5] = curNodes[3];
5086 //////////////////////////////////// HEXAHEDRON
5088 SMDS_VolumeTool hexa (elem);
5089 hexa.SetExternalNormal();
5090 if ( nbUniqueNodes == 4 && nbRepl == 6 ) {
5091 //////////////////////// ---> tetrahedron
5092 for ( int iFace = 0; iFace < 6; iFace++ ) {
5093 const int *ind = hexa.GetFaceNodesIndices( iFace ); // indices of face nodes
5094 if (curNodes[ind[ 0 ]] == curNodes[ind[ 1 ]] &&
5095 curNodes[ind[ 0 ]] == curNodes[ind[ 2 ]] &&
5096 curNodes[ind[ 0 ]] == curNodes[ind[ 3 ]] ) {
5097 // one face turns into a point ...
5098 int iOppFace = hexa.GetOppFaceIndex( iFace );
5099 ind = hexa.GetFaceNodesIndices( iOppFace );
5101 iUnique = 2; // reverse a tetrahedron bottom
5102 for ( iCur = 0; iCur < 4 && nbStick < 2; iCur++ ) {
5103 if ( curNodes[ind[ iCur ]] == curNodes[ind[ iCur + 1 ]] )
5105 else if ( iUnique >= 0 )
5106 uniqueNodes[ iUnique-- ] = curNodes[ind[ iCur ]];
5108 if ( nbStick == 1 ) {
5109 // ... and the opposite one - into a triangle.
5111 ind = hexa.GetFaceNodesIndices( iFace );
5112 uniqueNodes[ 3 ] = curNodes[ind[ 0 ]];
5119 else if (nbUniqueNodes == 5 && nbRepl == 4 ) {
5120 //////////////////// HEXAHEDRON ---> 2 tetrahedrons
5121 for ( int iFace = 0; iFace < 6; iFace++ ) {
5122 const int *ind = hexa.GetFaceNodesIndices( iFace ); // indices of face nodes
5123 if (curNodes[ind[ 0 ]] == curNodes[ind[ 1 ]] &&
5124 curNodes[ind[ 0 ]] == curNodes[ind[ 2 ]] &&
5125 curNodes[ind[ 0 ]] == curNodes[ind[ 3 ]] ) {
5126 // one face turns into a point ...
5127 int iOppFace = hexa.GetOppFaceIndex( iFace );
5128 ind = hexa.GetFaceNodesIndices( iOppFace );
5130 iUnique = 2; // reverse a tetrahedron 1 bottom
5131 for ( iCur = 0; iCur < 4 && nbStick == 0; iCur++ ) {
5132 if ( curNodes[ind[ iCur ]] == curNodes[ind[ iCur + 1 ]] )
5134 else if ( iUnique >= 0 )
5135 uniqueNodes[ iUnique-- ] = curNodes[ind[ iCur ]];
5137 if ( nbStick == 0 ) {
5138 // ... and the opposite one is a quadrangle
5140 const int* indTop = hexa.GetFaceNodesIndices( iFace );
5141 uniqueNodes[ 3 ] = curNodes[indTop[ 0 ]];
5144 SMDS_MeshElement* newElem =
5145 aMesh->AddVolume(curNodes[ind[ 0 ]],
5148 curNodes[indTop[ 0 ]]);
5149 myLastCreatedElems.Append(newElem);
5151 aMesh->SetMeshElementOnShape( newElem, aShapeId );
5158 else if ( nbUniqueNodes == 6 && nbRepl == 4 ) {
5159 ////////////////// HEXAHEDRON ---> 2 tetrahedrons or 1 prism
5160 // find indices of quad and tri faces
5161 int iQuadFace[ 6 ], iTriFace[ 6 ], nbQuad = 0, nbTri = 0, iFace;
5162 for ( iFace = 0; iFace < 6; iFace++ ) {
5163 const int *ind = hexa.GetFaceNodesIndices( iFace ); // indices of face nodes
5165 for ( iCur = 0; iCur < 4; iCur++ )
5166 nodeSet.insert( curNodes[ind[ iCur ]] );
5167 nbUniqueNodes = nodeSet.size();
5168 if ( nbUniqueNodes == 3 )
5169 iTriFace[ nbTri++ ] = iFace;
5170 else if ( nbUniqueNodes == 4 )
5171 iQuadFace[ nbQuad++ ] = iFace;
5173 if (nbQuad == 2 && nbTri == 4 &&
5174 hexa.GetOppFaceIndex( iQuadFace[ 0 ] ) == iQuadFace[ 1 ]) {
5175 // 2 opposite quadrangles stuck with a diagonal;
5176 // sample groups of merged indices: (0-4)(2-6)
5177 // --------------------------------------------> 2 tetrahedrons
5178 const int *ind1 = hexa.GetFaceNodesIndices( iQuadFace[ 0 ]); // indices of quad1 nodes
5179 const int *ind2 = hexa.GetFaceNodesIndices( iQuadFace[ 1 ]);
5180 int i0, i1d, i2, i3d, i0t, i2t; // d-daigonal, t-top
5181 if (curNodes[ind1[ 0 ]] == curNodes[ind2[ 0 ]] &&
5182 curNodes[ind1[ 2 ]] == curNodes[ind2[ 2 ]]) {
5183 // stuck with 0-2 diagonal
5191 else if (curNodes[ind1[ 1 ]] == curNodes[ind2[ 3 ]] &&
5192 curNodes[ind1[ 3 ]] == curNodes[ind2[ 1 ]]) {
5193 // stuck with 1-3 diagonal
5205 uniqueNodes[ 0 ] = curNodes [ i0 ];
5206 uniqueNodes[ 1 ] = curNodes [ i1d ];
5207 uniqueNodes[ 2 ] = curNodes [ i3d ];
5208 uniqueNodes[ 3 ] = curNodes [ i0t ];
5211 SMDS_MeshElement* newElem = aMesh->AddVolume(curNodes[ i1d ],
5215 myLastCreatedElems.Append(newElem);
5217 aMesh->SetMeshElementOnShape( newElem, aShapeId );
5220 else if (( nbTri == 2 && nbQuad == 3 ) || // merged (0-4)(1-5)
5221 ( nbTri == 4 && nbQuad == 2 )) { // merged (7-4)(1-5)
5222 // --------------------------------------------> prism
5223 // find 2 opposite triangles
5225 for ( iFace = 0; iFace + 1 < nbTri; iFace++ ) {
5226 if ( hexa.GetOppFaceIndex( iTriFace[ iFace ] ) == iTriFace[ iFace + 1 ]) {
5227 // find indices of kept and replaced nodes
5228 // and fill unique nodes of 2 opposite triangles
5229 const int *ind1 = hexa.GetFaceNodesIndices( iTriFace[ iFace ]);
5230 const int *ind2 = hexa.GetFaceNodesIndices( iTriFace[ iFace + 1 ]);
5231 const SMDS_MeshNode** hexanodes = hexa.GetNodes();
5232 // fill unique nodes
5235 for ( iCur = 0; iCur < 4 && isOk; iCur++ ) {
5236 const SMDS_MeshNode* n = curNodes[ind1[ iCur ]];
5237 const SMDS_MeshNode* nInit = hexanodes[ind1[ iCur ]];
5239 // iCur of a linked node of the opposite face (make normals co-directed):
5240 int iCurOpp = ( iCur == 1 || iCur == 3 ) ? 4 - iCur : iCur;
5241 // check that correspondent corners of triangles are linked
5242 if ( !hexa.IsLinked( ind1[ iCur ], ind2[ iCurOpp ] ))
5245 uniqueNodes[ iUnique ] = n;
5246 uniqueNodes[ iUnique + 3 ] = curNodes[ind2[ iCurOpp ]];
5255 } // if ( nbUniqueNodes == 6 && nbRepl == 4 )
5261 } // switch ( nbNodes )
5263 } // if ( nbNodes != nbUniqueNodes ) // some nodes stick
5266 if (elem->IsPoly() && elem->GetType() == SMDSAbs_Volume) {
5267 // Change nodes of polyedre
5268 const SMDS_PolyhedralVolumeOfNodes* aPolyedre =
5269 static_cast<const SMDS_PolyhedralVolumeOfNodes*>( elem );
5271 int nbFaces = aPolyedre->NbFaces();
5273 vector<const SMDS_MeshNode *> poly_nodes;
5274 vector<int> quantities (nbFaces);
5276 for (int iface = 1; iface <= nbFaces; iface++) {
5277 int inode, nbFaceNodes = aPolyedre->NbFaceNodes(iface);
5278 quantities[iface - 1] = nbFaceNodes;
5280 for (inode = 1; inode <= nbFaceNodes; inode++) {
5281 const SMDS_MeshNode* curNode = aPolyedre->GetFaceNode(iface, inode);
5283 TNodeNodeMap::iterator nnIt = nodeNodeMap.find( curNode );
5284 if (nnIt != nodeNodeMap.end()) { // curNode sticks
5285 curNode = (*nnIt).second;
5287 poly_nodes.push_back(curNode);
5290 aMesh->ChangePolyhedronNodes( elem, poly_nodes, quantities );
5294 // Change regular element or polygon
5295 aMesh->ChangeElementNodes( elem, & uniqueNodes[0], nbUniqueNodes );
5299 // Remove invalid regular element or invalid polygon
5300 rmElemIds.push_back( elem->GetID() );
5303 } // loop on elements
5305 // Remove equal nodes and bad elements
5307 Remove( rmNodeIds, true );
5308 Remove( rmElemIds, false );
5313 // ========================================================
5314 // class : SortableElement
5315 // purpose : allow sorting elements basing on their nodes
5316 // ========================================================
5317 class SortableElement : public set <const SMDS_MeshElement*>
5321 SortableElement( const SMDS_MeshElement* theElem )
5324 SMDS_ElemIteratorPtr nodeIt = theElem->nodesIterator();
5325 while ( nodeIt->more() )
5326 this->insert( nodeIt->next() );
5329 const SMDS_MeshElement* Get() const
5332 void Set(const SMDS_MeshElement* e) const
5337 mutable const SMDS_MeshElement* myElem;
5340 //=======================================================================
5341 //function : FindEqualElements
5342 //purpose : Return list of group of elements built on the same nodes.
5343 // Search among theElements or in the whole mesh if theElements is empty
5344 //=======================================================================
5345 void SMESH_MeshEditor::FindEqualElements(set<const SMDS_MeshElement*> & theElements,
5346 TListOfListOfElementsID & theGroupsOfElementsID)
5348 myLastCreatedElems.Clear();
5349 myLastCreatedNodes.Clear();
5351 typedef set<const SMDS_MeshElement*> TElemsSet;
5352 typedef map< SortableElement, int > TMapOfNodeSet;
5353 typedef list<int> TGroupOfElems;
5356 if ( theElements.empty() )
5357 { // get all elements in the mesh
5358 SMDS_ElemIteratorPtr eIt = GetMeshDS()->elementsIterator();
5359 while ( eIt->more() )
5360 elems.insert( elems.end(), eIt->next());
5363 elems = theElements;
5365 vector< TGroupOfElems > arrayOfGroups;
5366 TGroupOfElems groupOfElems;
5367 TMapOfNodeSet mapOfNodeSet;
5369 TElemsSet::iterator elemIt = elems.begin();
5370 for ( int i = 0, j=0; elemIt != elems.end(); ++elemIt, ++j ) {
5371 const SMDS_MeshElement* curElem = *elemIt;
5372 SortableElement SE(curElem);
5375 pair< TMapOfNodeSet::iterator, bool> pp = mapOfNodeSet.insert(make_pair(SE, i));
5376 if( !(pp.second) ) {
5377 TMapOfNodeSet::iterator& itSE = pp.first;
5378 ind = (*itSE).second;
5379 arrayOfGroups[ind].push_back(curElem->GetID());
5382 groupOfElems.clear();
5383 groupOfElems.push_back(curElem->GetID());
5384 arrayOfGroups.push_back(groupOfElems);
5389 vector< TGroupOfElems >::iterator groupIt = arrayOfGroups.begin();
5390 for ( ; groupIt != arrayOfGroups.end(); ++groupIt ) {
5391 groupOfElems = *groupIt;
5392 if ( groupOfElems.size() > 1 ) {
5393 groupOfElems.sort();
5394 theGroupsOfElementsID.push_back(groupOfElems);
5399 //=======================================================================
5400 //function : MergeElements
5401 //purpose : In each given group, substitute all elements by the first one.
5402 //=======================================================================
5404 void SMESH_MeshEditor::MergeElements(TListOfListOfElementsID & theGroupsOfElementsID)
5406 myLastCreatedElems.Clear();
5407 myLastCreatedNodes.Clear();
5409 typedef list<int> TListOfIDs;
5410 TListOfIDs rmElemIds; // IDs of elems to remove
5412 SMESHDS_Mesh* aMesh = GetMeshDS();
5414 TListOfListOfElementsID::iterator groupsIt = theGroupsOfElementsID.begin();
5415 while ( groupsIt != theGroupsOfElementsID.end() ) {
5416 TListOfIDs& aGroupOfElemID = *groupsIt;
5417 aGroupOfElemID.sort();
5418 int elemIDToKeep = aGroupOfElemID.front();
5419 const SMDS_MeshElement* elemToKeep = aMesh->FindElement(elemIDToKeep);
5420 aGroupOfElemID.pop_front();
5421 TListOfIDs::iterator idIt = aGroupOfElemID.begin();
5422 while ( idIt != aGroupOfElemID.end() ) {
5423 int elemIDToRemove = *idIt;
5424 const SMDS_MeshElement* elemToRemove = aMesh->FindElement(elemIDToRemove);
5425 // add the kept element in groups of removed one (PAL15188)
5426 AddToSameGroups( elemToKeep, elemToRemove, aMesh );
5427 rmElemIds.push_back( elemIDToRemove );
5433 Remove( rmElemIds, false );
5436 //=======================================================================
5437 //function : MergeEqualElements
5438 //purpose : Remove all but one of elements built on the same nodes.
5439 //=======================================================================
5441 void SMESH_MeshEditor::MergeEqualElements()
5443 set<const SMDS_MeshElement*> aMeshElements; /* empty input -
5444 to merge equal elements in the whole mesh */
5445 TListOfListOfElementsID aGroupsOfElementsID;
5446 FindEqualElements(aMeshElements, aGroupsOfElementsID);
5447 MergeElements(aGroupsOfElementsID);
5450 //=======================================================================
5451 //function : FindFaceInSet
5452 //purpose : Return a face having linked nodes n1 and n2 and which is
5453 // - not in avoidSet,
5454 // - in elemSet provided that !elemSet.empty()
5455 //=======================================================================
5457 const SMDS_MeshElement*
5458 SMESH_MeshEditor::FindFaceInSet(const SMDS_MeshNode* n1,
5459 const SMDS_MeshNode* n2,
5460 const TIDSortedElemSet& elemSet,
5461 const TIDSortedElemSet& avoidSet)
5464 SMDS_ElemIteratorPtr invElemIt = n1->GetInverseElementIterator(SMDSAbs_Face);
5465 while ( invElemIt->more() ) { // loop on inverse elements of n1
5466 const SMDS_MeshElement* elem = invElemIt->next();
5467 if (avoidSet.find( elem ) != avoidSet.end() )
5469 if ( !elemSet.empty() && elemSet.find( elem ) == elemSet.end())
5471 // get face nodes and find index of n1
5472 int i1, nbN = elem->NbNodes(), iNode = 0;
5473 //const SMDS_MeshNode* faceNodes[ nbN ], *n;
5474 vector<const SMDS_MeshNode*> faceNodes( nbN );
5475 const SMDS_MeshNode* n;
5476 SMDS_ElemIteratorPtr nIt = elem->nodesIterator();
5477 while ( nIt->more() ) {
5478 faceNodes[ iNode ] = static_cast<const SMDS_MeshNode*>( nIt->next() );
5479 if ( faceNodes[ iNode++ ] == n1 )
5482 // find a n2 linked to n1
5483 if(!elem->IsQuadratic()) {
5484 for ( iNode = 0; iNode < 2; iNode++ ) {
5485 if ( iNode ) // node before n1
5486 n = faceNodes[ i1 == 0 ? nbN - 1 : i1 - 1 ];
5487 else // node after n1
5488 n = faceNodes[ i1 + 1 == nbN ? 0 : i1 + 1 ];
5493 else { // analysis for quadratic elements
5494 bool IsFind = false;
5495 // check using only corner nodes
5496 for ( iNode = 0; iNode < 2; iNode++ ) {
5497 if ( iNode ) // node before n1
5498 n = faceNodes[ i1 == 0 ? nbN/2 - 1 : i1 - 1 ];
5499 else // node after n1
5500 n = faceNodes[ i1 + 1 == nbN/2 ? 0 : i1 + 1 ];
5508 // check using all nodes
5509 const SMDS_QuadraticFaceOfNodes* F =
5510 static_cast<const SMDS_QuadraticFaceOfNodes*>(elem);
5511 // use special nodes iterator
5513 SMDS_NodeIteratorPtr anIter = F->interlacedNodesIterator();
5514 while ( anIter->more() ) {
5515 faceNodes[iNode] = static_cast<const SMDS_MeshNode*>(anIter->next());
5516 if ( faceNodes[ iNode++ ] == n1 )
5519 for ( iNode = 0; iNode < 2; iNode++ ) {
5520 if ( iNode ) // node before n1
5521 n = faceNodes[ i1 == 0 ? nbN - 1 : i1 - 1 ];
5522 else // node after n1
5523 n = faceNodes[ i1 + 1 == nbN ? 0 : i1 + 1 ];
5529 } // end analysis for quadratic elements
5534 //=======================================================================
5535 //function : findAdjacentFace
5537 //=======================================================================
5539 static const SMDS_MeshElement* findAdjacentFace(const SMDS_MeshNode* n1,
5540 const SMDS_MeshNode* n2,
5541 const SMDS_MeshElement* elem)
5543 TIDSortedElemSet elemSet, avoidSet;
5545 avoidSet.insert ( elem );
5546 return SMESH_MeshEditor::FindFaceInSet( n1, n2, elemSet, avoidSet );
5549 //=======================================================================
5550 //function : FindFreeBorder
5552 //=======================================================================
5554 #define ControlFreeBorder SMESH::Controls::FreeEdges::IsFreeEdge
5556 bool SMESH_MeshEditor::FindFreeBorder (const SMDS_MeshNode* theFirstNode,
5557 const SMDS_MeshNode* theSecondNode,
5558 const SMDS_MeshNode* theLastNode,
5559 list< const SMDS_MeshNode* > & theNodes,
5560 list< const SMDS_MeshElement* >& theFaces)
5562 if ( !theFirstNode || !theSecondNode )
5564 // find border face between theFirstNode and theSecondNode
5565 const SMDS_MeshElement* curElem = findAdjacentFace( theFirstNode, theSecondNode, 0 );
5569 theFaces.push_back( curElem );
5570 theNodes.push_back( theFirstNode );
5571 theNodes.push_back( theSecondNode );
5573 //vector<const SMDS_MeshNode*> nodes;
5574 const SMDS_MeshNode *nIgnore = theFirstNode, *nStart = theSecondNode;
5575 set < const SMDS_MeshElement* > foundElems;
5576 bool needTheLast = ( theLastNode != 0 );
5578 while ( nStart != theLastNode ) {
5579 if ( nStart == theFirstNode )
5580 return !needTheLast;
5582 // find all free border faces sharing form nStart
5584 list< const SMDS_MeshElement* > curElemList;
5585 list< const SMDS_MeshNode* > nStartList;
5586 SMDS_ElemIteratorPtr invElemIt = nStart->GetInverseElementIterator(SMDSAbs_Face);
5587 while ( invElemIt->more() ) {
5588 const SMDS_MeshElement* e = invElemIt->next();
5589 if ( e == curElem || foundElems.insert( e ).second ) {
5591 int iNode = 0, nbNodes = e->NbNodes();
5592 //const SMDS_MeshNode* nodes[nbNodes+1];
5593 vector<const SMDS_MeshNode*> nodes(nbNodes+1);
5595 if(e->IsQuadratic()) {
5596 const SMDS_QuadraticFaceOfNodes* F =
5597 static_cast<const SMDS_QuadraticFaceOfNodes*>(e);
5598 // use special nodes iterator
5599 SMDS_NodeIteratorPtr anIter = F->interlacedNodesIterator();
5600 while( anIter->more() ) {
5601 nodes[ iNode++ ] = anIter->next();
5605 SMDS_ElemIteratorPtr nIt = e->nodesIterator();
5606 while ( nIt->more() )
5607 nodes[ iNode++ ] = static_cast<const SMDS_MeshNode*>( nIt->next() );
5609 nodes[ iNode ] = nodes[ 0 ];
5611 for ( iNode = 0; iNode < nbNodes; iNode++ )
5612 if (((nodes[ iNode ] == nStart && nodes[ iNode + 1] != nIgnore ) ||
5613 (nodes[ iNode + 1] == nStart && nodes[ iNode ] != nIgnore )) &&
5614 ControlFreeBorder( &nodes[ iNode ], e->GetID() ))
5616 nStartList.push_back( nodes[ iNode + ( nodes[ iNode ] == nStart ? 1 : 0 )]);
5617 curElemList.push_back( e );
5621 // analyse the found
5623 int nbNewBorders = curElemList.size();
5624 if ( nbNewBorders == 0 ) {
5625 // no free border furthermore
5626 return !needTheLast;
5628 else if ( nbNewBorders == 1 ) {
5629 // one more element found
5631 nStart = nStartList.front();
5632 curElem = curElemList.front();
5633 theFaces.push_back( curElem );
5634 theNodes.push_back( nStart );
5637 // several continuations found
5638 list< const SMDS_MeshElement* >::iterator curElemIt;
5639 list< const SMDS_MeshNode* >::iterator nStartIt;
5640 // check if one of them reached the last node
5641 if ( needTheLast ) {
5642 for (curElemIt = curElemList.begin(), nStartIt = nStartList.begin();
5643 curElemIt!= curElemList.end();
5644 curElemIt++, nStartIt++ )
5645 if ( *nStartIt == theLastNode ) {
5646 theFaces.push_back( *curElemIt );
5647 theNodes.push_back( *nStartIt );
5651 // find the best free border by the continuations
5652 list<const SMDS_MeshNode*> contNodes[ 2 ], *cNL;
5653 list<const SMDS_MeshElement*> contFaces[ 2 ], *cFL;
5654 for (curElemIt = curElemList.begin(), nStartIt = nStartList.begin();
5655 curElemIt!= curElemList.end();
5656 curElemIt++, nStartIt++ )
5658 cNL = & contNodes[ contNodes[0].empty() ? 0 : 1 ];
5659 cFL = & contFaces[ contFaces[0].empty() ? 0 : 1 ];
5660 // find one more free border
5661 if ( ! FindFreeBorder( nStart, *nStartIt, theLastNode, *cNL, *cFL )) {
5665 else if ( !contNodes[0].empty() && !contNodes[1].empty() ) {
5666 // choice: clear a worse one
5667 int iLongest = ( contNodes[0].size() < contNodes[1].size() ? 1 : 0 );
5668 int iWorse = ( needTheLast ? 1 - iLongest : iLongest );
5669 contNodes[ iWorse ].clear();
5670 contFaces[ iWorse ].clear();
5673 if ( contNodes[0].empty() && contNodes[1].empty() )
5676 // append the best free border
5677 cNL = & contNodes[ contNodes[0].empty() ? 1 : 0 ];
5678 cFL = & contFaces[ contFaces[0].empty() ? 1 : 0 ];
5679 theNodes.pop_back(); // remove nIgnore
5680 theNodes.pop_back(); // remove nStart
5681 theFaces.pop_back(); // remove curElem
5682 list< const SMDS_MeshNode* >::iterator nIt = cNL->begin();
5683 list< const SMDS_MeshElement* >::iterator fIt = cFL->begin();
5684 for ( ; nIt != cNL->end(); nIt++ ) theNodes.push_back( *nIt );
5685 for ( ; fIt != cFL->end(); fIt++ ) theFaces.push_back( *fIt );
5688 } // several continuations found
5689 } // while ( nStart != theLastNode )
5694 //=======================================================================
5695 //function : CheckFreeBorderNodes
5696 //purpose : Return true if the tree nodes are on a free border
5697 //=======================================================================
5699 bool SMESH_MeshEditor::CheckFreeBorderNodes(const SMDS_MeshNode* theNode1,
5700 const SMDS_MeshNode* theNode2,
5701 const SMDS_MeshNode* theNode3)
5703 list< const SMDS_MeshNode* > nodes;
5704 list< const SMDS_MeshElement* > faces;
5705 return FindFreeBorder( theNode1, theNode2, theNode3, nodes, faces);
5708 //=======================================================================
5709 //function : SewFreeBorder
5711 //=======================================================================
5713 SMESH_MeshEditor::Sew_Error
5714 SMESH_MeshEditor::SewFreeBorder (const SMDS_MeshNode* theBordFirstNode,
5715 const SMDS_MeshNode* theBordSecondNode,
5716 const SMDS_MeshNode* theBordLastNode,
5717 const SMDS_MeshNode* theSideFirstNode,
5718 const SMDS_MeshNode* theSideSecondNode,
5719 const SMDS_MeshNode* theSideThirdNode,
5720 const bool theSideIsFreeBorder,
5721 const bool toCreatePolygons,
5722 const bool toCreatePolyedrs)
5724 myLastCreatedElems.Clear();
5725 myLastCreatedNodes.Clear();
5727 MESSAGE("::SewFreeBorder()");
5728 Sew_Error aResult = SEW_OK;
5730 // ====================================
5731 // find side nodes and elements
5732 // ====================================
5734 list< const SMDS_MeshNode* > nSide[ 2 ];
5735 list< const SMDS_MeshElement* > eSide[ 2 ];
5736 list< const SMDS_MeshNode* >::iterator nIt[ 2 ];
5737 list< const SMDS_MeshElement* >::iterator eIt[ 2 ];
5741 if (!FindFreeBorder(theBordFirstNode,theBordSecondNode,theBordLastNode,
5742 nSide[0], eSide[0])) {
5743 MESSAGE(" Free Border 1 not found " );
5744 aResult = SEW_BORDER1_NOT_FOUND;
5746 if (theSideIsFreeBorder) {
5749 if (!FindFreeBorder(theSideFirstNode, theSideSecondNode, theSideThirdNode,
5750 nSide[1], eSide[1])) {
5751 MESSAGE(" Free Border 2 not found " );
5752 aResult = ( aResult != SEW_OK ? SEW_BOTH_BORDERS_NOT_FOUND : SEW_BORDER2_NOT_FOUND );
5755 if ( aResult != SEW_OK )
5758 if (!theSideIsFreeBorder) {
5762 // -------------------------------------------------------------------------
5764 // 1. If nodes to merge are not coincident, move nodes of the free border
5765 // from the coord sys defined by the direction from the first to last
5766 // nodes of the border to the correspondent sys of the side 2
5767 // 2. On the side 2, find the links most co-directed with the correspondent
5768 // links of the free border
5769 // -------------------------------------------------------------------------
5771 // 1. Since sewing may brake if there are volumes to split on the side 2,
5772 // we wont move nodes but just compute new coordinates for them
5773 typedef map<const SMDS_MeshNode*, gp_XYZ> TNodeXYZMap;
5774 TNodeXYZMap nBordXYZ;
5775 list< const SMDS_MeshNode* >& bordNodes = nSide[ 0 ];
5776 list< const SMDS_MeshNode* >::iterator nBordIt;
5778 gp_XYZ Pb1( theBordFirstNode->X(), theBordFirstNode->Y(), theBordFirstNode->Z() );
5779 gp_XYZ Pb2( theBordLastNode->X(), theBordLastNode->Y(), theBordLastNode->Z() );
5780 gp_XYZ Ps1( theSideFirstNode->X(), theSideFirstNode->Y(), theSideFirstNode->Z() );
5781 gp_XYZ Ps2( theSideSecondNode->X(), theSideSecondNode->Y(), theSideSecondNode->Z() );
5782 double tol2 = 1.e-8;
5783 gp_Vec Vbs1( Pb1 - Ps1 ),Vbs2( Pb2 - Ps2 );
5784 if ( Vbs1.SquareMagnitude() > tol2 || Vbs2.SquareMagnitude() > tol2 ) {
5785 // Need node movement.
5787 // find X and Z axes to create trsf
5788 gp_Vec Zb( Pb1 - Pb2 ), Zs( Ps1 - Ps2 );
5790 if ( X.SquareMagnitude() <= gp::Resolution() * gp::Resolution() )
5792 X = gp_Ax2( gp::Origin(), Zb ).XDirection();
5795 gp_Ax3 toBordAx( Pb1, Zb, X );
5796 gp_Ax3 fromSideAx( Ps1, Zs, X );
5797 gp_Ax3 toGlobalAx( gp::Origin(), gp::DZ(), gp::DX() );
5799 gp_Trsf toBordSys, fromSide2Sys;
5800 toBordSys.SetTransformation( toBordAx );
5801 fromSide2Sys.SetTransformation( fromSideAx, toGlobalAx );
5802 fromSide2Sys.SetScaleFactor( Zs.Magnitude() / Zb.Magnitude() );
5805 for ( nBordIt = bordNodes.begin(); nBordIt != bordNodes.end(); nBordIt++ ) {
5806 const SMDS_MeshNode* n = *nBordIt;
5807 gp_XYZ xyz( n->X(),n->Y(),n->Z() );
5808 toBordSys.Transforms( xyz );
5809 fromSide2Sys.Transforms( xyz );
5810 nBordXYZ.insert( TNodeXYZMap::value_type( n, xyz ));
5814 // just insert nodes XYZ in the nBordXYZ map
5815 for ( nBordIt = bordNodes.begin(); nBordIt != bordNodes.end(); nBordIt++ ) {
5816 const SMDS_MeshNode* n = *nBordIt;
5817 nBordXYZ.insert( TNodeXYZMap::value_type( n, gp_XYZ( n->X(),n->Y(),n->Z() )));
5821 // 2. On the side 2, find the links most co-directed with the correspondent
5822 // links of the free border
5824 list< const SMDS_MeshElement* >& sideElems = eSide[ 1 ];
5825 list< const SMDS_MeshNode* >& sideNodes = nSide[ 1 ];
5826 sideNodes.push_back( theSideFirstNode );
5828 bool hasVolumes = false;
5829 LinkID_Gen aLinkID_Gen( GetMeshDS() );
5830 set<long> foundSideLinkIDs, checkedLinkIDs;
5831 SMDS_VolumeTool volume;
5832 //const SMDS_MeshNode* faceNodes[ 4 ];
5834 const SMDS_MeshNode* sideNode;
5835 const SMDS_MeshElement* sideElem;
5836 const SMDS_MeshNode* prevSideNode = theSideFirstNode;
5837 const SMDS_MeshNode* prevBordNode = theBordFirstNode;
5838 nBordIt = bordNodes.begin();
5840 // border node position and border link direction to compare with
5841 gp_XYZ bordPos = nBordXYZ[ *nBordIt ];
5842 gp_XYZ bordDir = bordPos - nBordXYZ[ prevBordNode ];
5843 // choose next side node by link direction or by closeness to
5844 // the current border node:
5845 bool searchByDir = ( *nBordIt != theBordLastNode );
5847 // find the next node on the Side 2
5849 double maxDot = -DBL_MAX, minDist = DBL_MAX;
5851 checkedLinkIDs.clear();
5852 gp_XYZ prevXYZ( prevSideNode->X(), prevSideNode->Y(), prevSideNode->Z() );
5854 // loop on inverse elements of current node (prevSideNode) on the Side 2
5855 SMDS_ElemIteratorPtr invElemIt = prevSideNode->GetInverseElementIterator();
5856 while ( invElemIt->more() )
5858 const SMDS_MeshElement* elem = invElemIt->next();
5859 // prepare data for a loop on links coming to prevSideNode, of a face or a volume
5860 int iPrevNode, iNode = 0, nbNodes = elem->NbNodes();
5861 vector< const SMDS_MeshNode* > faceNodes( nbNodes, (const SMDS_MeshNode*)0 );
5862 bool isVolume = volume.Set( elem );
5863 const SMDS_MeshNode** nodes = isVolume ? volume.GetNodes() : & faceNodes[0];
5864 if ( isVolume ) // --volume
5866 else if ( elem->GetType()==SMDSAbs_Face ) { // --face
5867 // retrieve all face nodes and find iPrevNode - an index of the prevSideNode
5868 if(elem->IsQuadratic()) {
5869 const SMDS_QuadraticFaceOfNodes* F =
5870 static_cast<const SMDS_QuadraticFaceOfNodes*>(elem);
5871 // use special nodes iterator
5872 SMDS_NodeIteratorPtr anIter = F->interlacedNodesIterator();
5873 while( anIter->more() ) {
5874 nodes[ iNode ] = anIter->next();
5875 if ( nodes[ iNode++ ] == prevSideNode )
5876 iPrevNode = iNode - 1;
5880 SMDS_ElemIteratorPtr nIt = elem->nodesIterator();
5881 while ( nIt->more() ) {
5882 nodes[ iNode ] = cast2Node( nIt->next() );
5883 if ( nodes[ iNode++ ] == prevSideNode )
5884 iPrevNode = iNode - 1;
5887 // there are 2 links to check
5892 // loop on links, to be precise, on the second node of links
5893 for ( iNode = 0; iNode < nbNodes; iNode++ ) {
5894 const SMDS_MeshNode* n = nodes[ iNode ];
5896 if ( !volume.IsLinked( n, prevSideNode ))
5900 if ( iNode ) // a node before prevSideNode
5901 n = nodes[ iPrevNode == 0 ? elem->NbNodes() - 1 : iPrevNode - 1 ];
5902 else // a node after prevSideNode
5903 n = nodes[ iPrevNode + 1 == elem->NbNodes() ? 0 : iPrevNode + 1 ];
5905 // check if this link was already used
5906 long iLink = aLinkID_Gen.GetLinkID( prevSideNode, n );
5907 bool isJustChecked = !checkedLinkIDs.insert( iLink ).second;
5908 if (!isJustChecked &&
5909 foundSideLinkIDs.find( iLink ) == foundSideLinkIDs.end() )
5911 // test a link geometrically
5912 gp_XYZ nextXYZ ( n->X(), n->Y(), n->Z() );
5913 bool linkIsBetter = false;
5914 double dot = 0.0, dist = 0.0;
5915 if ( searchByDir ) { // choose most co-directed link
5916 dot = bordDir * ( nextXYZ - prevXYZ ).Normalized();
5917 linkIsBetter = ( dot > maxDot );
5919 else { // choose link with the node closest to bordPos
5920 dist = ( nextXYZ - bordPos ).SquareModulus();
5921 linkIsBetter = ( dist < minDist );
5923 if ( linkIsBetter ) {
5932 } // loop on inverse elements of prevSideNode
5935 MESSAGE(" Cant find path by links of the Side 2 ");
5936 return SEW_BAD_SIDE_NODES;
5938 sideNodes.push_back( sideNode );
5939 sideElems.push_back( sideElem );
5940 foundSideLinkIDs.insert ( linkID );
5941 prevSideNode = sideNode;
5943 if ( *nBordIt == theBordLastNode )
5944 searchByDir = false;
5946 // find the next border link to compare with
5947 gp_XYZ sidePos( sideNode->X(), sideNode->Y(), sideNode->Z() );
5948 searchByDir = ( bordDir * ( sidePos - bordPos ) <= 0 );
5949 // move to next border node if sideNode is before forward border node (bordPos)
5950 while ( *nBordIt != theBordLastNode && !searchByDir ) {
5951 prevBordNode = *nBordIt;
5953 bordPos = nBordXYZ[ *nBordIt ];
5954 bordDir = bordPos - nBordXYZ[ prevBordNode ];
5955 searchByDir = ( bordDir * ( sidePos - bordPos ) <= 0 );
5959 while ( sideNode != theSideSecondNode );
5961 if ( hasVolumes && sideNodes.size () != bordNodes.size() && !toCreatePolyedrs) {
5962 MESSAGE("VOLUME SPLITTING IS FORBIDDEN");
5963 return SEW_VOLUMES_TO_SPLIT; // volume splitting is forbidden
5965 } // end nodes search on the side 2
5967 // ============================
5968 // sew the border to the side 2
5969 // ============================
5971 int nbNodes[] = { nSide[0].size(), nSide[1].size() };
5972 int maxNbNodes = Max( nbNodes[0], nbNodes[1] );
5974 TListOfListOfNodes nodeGroupsToMerge;
5975 if ( nbNodes[0] == nbNodes[1] ||
5976 ( theSideIsFreeBorder && !theSideThirdNode)) {
5978 // all nodes are to be merged
5980 for (nIt[0] = nSide[0].begin(), nIt[1] = nSide[1].begin();
5981 nIt[0] != nSide[0].end() && nIt[1] != nSide[1].end();
5982 nIt[0]++, nIt[1]++ )
5984 nodeGroupsToMerge.push_back( list<const SMDS_MeshNode*>() );
5985 nodeGroupsToMerge.back().push_back( *nIt[1] ); // to keep
5986 nodeGroupsToMerge.back().push_back( *nIt[0] ); // to remove
5991 // insert new nodes into the border and the side to get equal nb of segments
5993 // get normalized parameters of nodes on the borders
5994 //double param[ 2 ][ maxNbNodes ];
5996 param[0] = new double [ maxNbNodes ];
5997 param[1] = new double [ maxNbNodes ];
5999 for ( iBord = 0; iBord < 2; iBord++ ) { // loop on 2 borders
6000 list< const SMDS_MeshNode* >& nodes = nSide[ iBord ];
6001 list< const SMDS_MeshNode* >::iterator nIt = nodes.begin();
6002 const SMDS_MeshNode* nPrev = *nIt;
6003 double bordLength = 0;
6004 for ( iNode = 0; nIt != nodes.end(); nIt++, iNode++ ) { // loop on border nodes
6005 const SMDS_MeshNode* nCur = *nIt;
6006 gp_XYZ segment (nCur->X() - nPrev->X(),
6007 nCur->Y() - nPrev->Y(),
6008 nCur->Z() - nPrev->Z());
6009 double segmentLen = segment.Modulus();
6010 bordLength += segmentLen;
6011 param[ iBord ][ iNode ] = bordLength;
6014 // normalize within [0,1]
6015 for ( iNode = 0; iNode < nbNodes[ iBord ]; iNode++ ) {
6016 param[ iBord ][ iNode ] /= bordLength;
6020 // loop on border segments
6021 const SMDS_MeshNode *nPrev[ 2 ] = { 0, 0 };
6022 int i[ 2 ] = { 0, 0 };
6023 nIt[0] = nSide[0].begin(); eIt[0] = eSide[0].begin();
6024 nIt[1] = nSide[1].begin(); eIt[1] = eSide[1].begin();
6026 TElemOfNodeListMap insertMap;
6027 TElemOfNodeListMap::iterator insertMapIt;
6029 // key: elem to insert nodes into
6030 // value: 2 nodes to insert between + nodes to be inserted
6032 bool next[ 2 ] = { false, false };
6034 // find min adjacent segment length after sewing
6035 double nextParam = 10., prevParam = 0;
6036 for ( iBord = 0; iBord < 2; iBord++ ) { // loop on 2 borders
6037 if ( i[ iBord ] + 1 < nbNodes[ iBord ])
6038 nextParam = Min( nextParam, param[iBord][ i[iBord] + 1 ]);
6039 if ( i[ iBord ] > 0 )
6040 prevParam = Max( prevParam, param[iBord][ i[iBord] - 1 ]);
6042 double minParam = Min( param[ 0 ][ i[0] ], param[ 1 ][ i[1] ]);
6043 double maxParam = Max( param[ 0 ][ i[0] ], param[ 1 ][ i[1] ]);
6044 double minSegLen = Min( nextParam - minParam, maxParam - prevParam );
6046 // choose to insert or to merge nodes
6047 double du = param[ 1 ][ i[1] ] - param[ 0 ][ i[0] ];
6048 if ( Abs( du ) <= minSegLen * 0.2 ) {
6051 nodeGroupsToMerge.push_back( list<const SMDS_MeshNode*>() );
6052 const SMDS_MeshNode* n0 = *nIt[0];
6053 const SMDS_MeshNode* n1 = *nIt[1];
6054 nodeGroupsToMerge.back().push_back( n1 );
6055 nodeGroupsToMerge.back().push_back( n0 );
6056 // position of node of the border changes due to merge
6057 param[ 0 ][ i[0] ] += du;
6058 // move n1 for the sake of elem shape evaluation during insertion.
6059 // n1 will be removed by MergeNodes() anyway
6060 const_cast<SMDS_MeshNode*>( n0 )->setXYZ( n1->X(), n1->Y(), n1->Z() );
6061 next[0] = next[1] = true;
6066 int intoBord = ( du < 0 ) ? 0 : 1;
6067 const SMDS_MeshElement* elem = *eIt[ intoBord ];
6068 const SMDS_MeshNode* n1 = nPrev[ intoBord ];
6069 const SMDS_MeshNode* n2 = *nIt[ intoBord ];
6070 const SMDS_MeshNode* nIns = *nIt[ 1 - intoBord ];
6071 if ( intoBord == 1 ) {
6072 // move node of the border to be on a link of elem of the side
6073 gp_XYZ p1 (n1->X(), n1->Y(), n1->Z());
6074 gp_XYZ p2 (n2->X(), n2->Y(), n2->Z());
6075 double ratio = du / ( param[ 1 ][ i[1] ] - param[ 1 ][ i[1]-1 ]);
6076 gp_XYZ p = p2 * ( 1 - ratio ) + p1 * ratio;
6077 GetMeshDS()->MoveNode( nIns, p.X(), p.Y(), p.Z() );
6079 insertMapIt = insertMap.find( elem );
6080 bool notFound = ( insertMapIt == insertMap.end() );
6081 bool otherLink = ( !notFound && (*insertMapIt).second.front() != n1 );
6083 // insert into another link of the same element:
6084 // 1. perform insertion into the other link of the elem
6085 list<const SMDS_MeshNode*> & nodeList = (*insertMapIt).second;
6086 const SMDS_MeshNode* n12 = nodeList.front(); nodeList.pop_front();
6087 const SMDS_MeshNode* n22 = nodeList.front(); nodeList.pop_front();
6088 InsertNodesIntoLink( elem, n12, n22, nodeList, toCreatePolygons );
6089 // 2. perform insertion into the link of adjacent faces
6091 const SMDS_MeshElement* adjElem = findAdjacentFace( n12, n22, elem );
6093 InsertNodesIntoLink( adjElem, n12, n22, nodeList, toCreatePolygons );
6097 if (toCreatePolyedrs) {
6098 // perform insertion into the links of adjacent volumes
6099 UpdateVolumes(n12, n22, nodeList);
6101 // 3. find an element appeared on n1 and n2 after the insertion
6102 insertMap.erase( elem );
6103 elem = findAdjacentFace( n1, n2, 0 );
6105 if ( notFound || otherLink ) {
6106 // add element and nodes of the side into the insertMap
6107 insertMapIt = insertMap.insert
6108 ( TElemOfNodeListMap::value_type( elem, list<const SMDS_MeshNode*>() )).first;
6109 (*insertMapIt).second.push_back( n1 );
6110 (*insertMapIt).second.push_back( n2 );
6112 // add node to be inserted into elem
6113 (*insertMapIt).second.push_back( nIns );
6114 next[ 1 - intoBord ] = true;
6117 // go to the next segment
6118 for ( iBord = 0; iBord < 2; iBord++ ) { // loop on 2 borders
6119 if ( next[ iBord ] ) {
6120 if ( i[ iBord ] != 0 && eIt[ iBord ] != eSide[ iBord ].end())
6122 nPrev[ iBord ] = *nIt[ iBord ];
6123 nIt[ iBord ]++; i[ iBord ]++;
6127 while ( nIt[0] != nSide[0].end() && nIt[1] != nSide[1].end());
6129 // perform insertion of nodes into elements
6131 for (insertMapIt = insertMap.begin();
6132 insertMapIt != insertMap.end();
6135 const SMDS_MeshElement* elem = (*insertMapIt).first;
6136 list<const SMDS_MeshNode*> & nodeList = (*insertMapIt).second;
6137 const SMDS_MeshNode* n1 = nodeList.front(); nodeList.pop_front();
6138 const SMDS_MeshNode* n2 = nodeList.front(); nodeList.pop_front();
6140 InsertNodesIntoLink( elem, n1, n2, nodeList, toCreatePolygons );
6142 if ( !theSideIsFreeBorder ) {
6143 // look for and insert nodes into the faces adjacent to elem
6145 const SMDS_MeshElement* adjElem = findAdjacentFace( n1, n2, elem );
6147 InsertNodesIntoLink( adjElem, n1, n2, nodeList, toCreatePolygons );
6152 if (toCreatePolyedrs) {
6153 // perform insertion into the links of adjacent volumes
6154 UpdateVolumes(n1, n2, nodeList);
6160 } // end: insert new nodes
6162 MergeNodes ( nodeGroupsToMerge );
6167 //=======================================================================
6168 //function : InsertNodesIntoLink
6169 //purpose : insert theNodesToInsert into theFace between theBetweenNode1
6170 // and theBetweenNode2 and split theElement
6171 //=======================================================================
6173 void SMESH_MeshEditor::InsertNodesIntoLink(const SMDS_MeshElement* theFace,
6174 const SMDS_MeshNode* theBetweenNode1,
6175 const SMDS_MeshNode* theBetweenNode2,
6176 list<const SMDS_MeshNode*>& theNodesToInsert,
6177 const bool toCreatePoly)
6179 if ( theFace->GetType() != SMDSAbs_Face ) return;
6181 // find indices of 2 link nodes and of the rest nodes
6182 int iNode = 0, il1, il2, i3, i4;
6183 il1 = il2 = i3 = i4 = -1;
6184 //const SMDS_MeshNode* nodes[ theFace->NbNodes() ];
6185 vector<const SMDS_MeshNode*> nodes( theFace->NbNodes() );
6187 if(theFace->IsQuadratic()) {
6188 const SMDS_QuadraticFaceOfNodes* F =
6189 static_cast<const SMDS_QuadraticFaceOfNodes*>(theFace);
6190 // use special nodes iterator
6191 SMDS_NodeIteratorPtr anIter = F->interlacedNodesIterator();
6192 while( anIter->more() ) {
6193 const SMDS_MeshNode* n = anIter->next();
6194 if ( n == theBetweenNode1 )
6196 else if ( n == theBetweenNode2 )
6202 nodes[ iNode++ ] = n;
6206 SMDS_ElemIteratorPtr nodeIt = theFace->nodesIterator();
6207 while ( nodeIt->more() ) {
6208 const SMDS_MeshNode* n = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
6209 if ( n == theBetweenNode1 )
6211 else if ( n == theBetweenNode2 )
6217 nodes[ iNode++ ] = n;
6220 if ( il1 < 0 || il2 < 0 || i3 < 0 )
6223 // arrange link nodes to go one after another regarding the face orientation
6224 bool reverse = ( Abs( il2 - il1 ) == 1 ? il2 < il1 : il1 < il2 );
6225 list<const SMDS_MeshNode *> aNodesToInsert = theNodesToInsert;
6230 aNodesToInsert.reverse();
6232 // check that not link nodes of a quadrangles are in good order
6233 int nbFaceNodes = theFace->NbNodes();
6234 if ( nbFaceNodes == 4 && i4 - i3 != 1 ) {
6240 if (toCreatePoly || theFace->IsPoly()) {
6243 vector<const SMDS_MeshNode *> poly_nodes (nbFaceNodes + aNodesToInsert.size());
6245 // add nodes of face up to first node of link
6248 if(theFace->IsQuadratic()) {
6249 const SMDS_QuadraticFaceOfNodes* F =
6250 static_cast<const SMDS_QuadraticFaceOfNodes*>(theFace);
6251 // use special nodes iterator
6252 SMDS_NodeIteratorPtr anIter = F->interlacedNodesIterator();
6253 while( anIter->more() && !isFLN ) {
6254 const SMDS_MeshNode* n = anIter->next();
6255 poly_nodes[iNode++] = n;
6256 if (n == nodes[il1]) {
6260 // add nodes to insert
6261 list<const SMDS_MeshNode*>::iterator nIt = aNodesToInsert.begin();
6262 for (; nIt != aNodesToInsert.end(); nIt++) {
6263 poly_nodes[iNode++] = *nIt;
6265 // add nodes of face starting from last node of link
6266 while ( anIter->more() ) {
6267 poly_nodes[iNode++] = anIter->next();
6271 SMDS_ElemIteratorPtr nodeIt = theFace->nodesIterator();
6272 while ( nodeIt->more() && !isFLN ) {
6273 const SMDS_MeshNode* n = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
6274 poly_nodes[iNode++] = n;
6275 if (n == nodes[il1]) {
6279 // add nodes to insert
6280 list<const SMDS_MeshNode*>::iterator nIt = aNodesToInsert.begin();
6281 for (; nIt != aNodesToInsert.end(); nIt++) {
6282 poly_nodes[iNode++] = *nIt;
6284 // add nodes of face starting from last node of link
6285 while ( nodeIt->more() ) {
6286 const SMDS_MeshNode* n = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
6287 poly_nodes[iNode++] = n;
6291 // edit or replace the face
6292 SMESHDS_Mesh *aMesh = GetMeshDS();
6294 if (theFace->IsPoly()) {
6295 aMesh->ChangePolygonNodes(theFace, poly_nodes);
6298 int aShapeId = FindShape( theFace );
6300 SMDS_MeshElement* newElem = aMesh->AddPolygonalFace(poly_nodes);
6301 myLastCreatedElems.Append(newElem);
6302 if ( aShapeId && newElem )
6303 aMesh->SetMeshElementOnShape( newElem, aShapeId );
6305 aMesh->RemoveElement(theFace);
6310 if( !theFace->IsQuadratic() ) {
6312 // put aNodesToInsert between theBetweenNode1 and theBetweenNode2
6313 int nbLinkNodes = 2 + aNodesToInsert.size();
6314 //const SMDS_MeshNode* linkNodes[ nbLinkNodes ];
6315 vector<const SMDS_MeshNode*> linkNodes( nbLinkNodes );
6316 linkNodes[ 0 ] = nodes[ il1 ];
6317 linkNodes[ nbLinkNodes - 1 ] = nodes[ il2 ];
6318 list<const SMDS_MeshNode*>::iterator nIt = aNodesToInsert.begin();
6319 for ( iNode = 1; nIt != aNodesToInsert.end(); nIt++ ) {
6320 linkNodes[ iNode++ ] = *nIt;
6322 // decide how to split a quadrangle: compare possible variants
6323 // and choose which of splits to be a quadrangle
6324 int i1, i2, iSplit, nbSplits = nbLinkNodes - 1, iBestQuad;
6325 if ( nbFaceNodes == 3 ) {
6326 iBestQuad = nbSplits;
6329 else if ( nbFaceNodes == 4 ) {
6330 SMESH::Controls::NumericalFunctorPtr aCrit( new SMESH::Controls::AspectRatio);
6331 double aBestRate = DBL_MAX;
6332 for ( int iQuad = 0; iQuad < nbSplits; iQuad++ ) {
6334 double aBadRate = 0;
6335 // evaluate elements quality
6336 for ( iSplit = 0; iSplit < nbSplits; iSplit++ ) {
6337 if ( iSplit == iQuad ) {
6338 SMDS_FaceOfNodes quad (linkNodes[ i1++ ],
6342 aBadRate += getBadRate( &quad, aCrit );
6345 SMDS_FaceOfNodes tria (linkNodes[ i1++ ],
6347 nodes[ iSplit < iQuad ? i4 : i3 ]);
6348 aBadRate += getBadRate( &tria, aCrit );
6352 if ( aBadRate < aBestRate ) {
6354 aBestRate = aBadRate;
6359 // create new elements
6360 SMESHDS_Mesh *aMesh = GetMeshDS();
6361 int aShapeId = FindShape( theFace );
6364 for ( iSplit = 0; iSplit < nbSplits - 1; iSplit++ ) {
6365 SMDS_MeshElement* newElem = 0;
6366 if ( iSplit == iBestQuad )
6367 newElem = aMesh->AddFace (linkNodes[ i1++ ],
6372 newElem = aMesh->AddFace (linkNodes[ i1++ ],
6374 nodes[ iSplit < iBestQuad ? i4 : i3 ]);
6375 myLastCreatedElems.Append(newElem);
6376 if ( aShapeId && newElem )
6377 aMesh->SetMeshElementOnShape( newElem, aShapeId );
6380 // change nodes of theFace
6381 const SMDS_MeshNode* newNodes[ 4 ];
6382 newNodes[ 0 ] = linkNodes[ i1 ];
6383 newNodes[ 1 ] = linkNodes[ i2 ];
6384 newNodes[ 2 ] = nodes[ iSplit >= iBestQuad ? i3 : i4 ];
6385 newNodes[ 3 ] = nodes[ i4 ];
6386 aMesh->ChangeElementNodes( theFace, newNodes, iSplit == iBestQuad ? 4 : 3 );
6387 } // end if(!theFace->IsQuadratic())
6388 else { // theFace is quadratic
6389 // we have to split theFace on simple triangles and one simple quadrangle
6391 int nbshift = tmp*2;
6392 // shift nodes in nodes[] by nbshift
6394 for(i=0; i<nbshift; i++) {
6395 const SMDS_MeshNode* n = nodes[0];
6396 for(j=0; j<nbFaceNodes-1; j++) {
6397 nodes[j] = nodes[j+1];
6399 nodes[nbFaceNodes-1] = n;
6401 il1 = il1 - nbshift;
6402 // now have to insert nodes between n0 and n1 or n1 and n2 (see below)
6403 // n0 n1 n2 n0 n1 n2
6404 // +-----+-----+ +-----+-----+
6413 // create new elements
6414 SMESHDS_Mesh *aMesh = GetMeshDS();
6415 int aShapeId = FindShape( theFace );
6418 if(nbFaceNodes==6) { // quadratic triangle
6419 SMDS_MeshElement* newElem =
6420 aMesh->AddFace(nodes[3],nodes[4],nodes[5]);
6421 myLastCreatedElems.Append(newElem);
6422 if ( aShapeId && newElem )
6423 aMesh->SetMeshElementOnShape( newElem, aShapeId );
6424 if(theFace->IsMediumNode(nodes[il1])) {
6425 // create quadrangle
6426 newElem = aMesh->AddFace(nodes[0],nodes[1],nodes[3],nodes[5]);
6427 myLastCreatedElems.Append(newElem);
6428 if ( aShapeId && newElem )
6429 aMesh->SetMeshElementOnShape( newElem, aShapeId );
6435 // create quadrangle
6436 newElem = aMesh->AddFace(nodes[1],nodes[2],nodes[3],nodes[5]);
6437 myLastCreatedElems.Append(newElem);
6438 if ( aShapeId && newElem )
6439 aMesh->SetMeshElementOnShape( newElem, aShapeId );
6445 else { // nbFaceNodes==8 - quadratic quadrangle
6446 SMDS_MeshElement* newElem =
6447 aMesh->AddFace(nodes[3],nodes[4],nodes[5]);
6448 myLastCreatedElems.Append(newElem);
6449 if ( aShapeId && newElem )
6450 aMesh->SetMeshElementOnShape( newElem, aShapeId );
6451 newElem = aMesh->AddFace(nodes[5],nodes[6],nodes[7]);
6452 myLastCreatedElems.Append(newElem);
6453 if ( aShapeId && newElem )
6454 aMesh->SetMeshElementOnShape( newElem, aShapeId );
6455 newElem = aMesh->AddFace(nodes[5],nodes[7],nodes[3]);
6456 myLastCreatedElems.Append(newElem);
6457 if ( aShapeId && newElem )
6458 aMesh->SetMeshElementOnShape( newElem, aShapeId );
6459 if(theFace->IsMediumNode(nodes[il1])) {
6460 // create quadrangle
6461 newElem = aMesh->AddFace(nodes[0],nodes[1],nodes[3],nodes[7]);
6462 myLastCreatedElems.Append(newElem);
6463 if ( aShapeId && newElem )
6464 aMesh->SetMeshElementOnShape( newElem, aShapeId );
6470 // create quadrangle
6471 newElem = aMesh->AddFace(nodes[1],nodes[2],nodes[3],nodes[7]);
6472 myLastCreatedElems.Append(newElem);
6473 if ( aShapeId && newElem )
6474 aMesh->SetMeshElementOnShape( newElem, aShapeId );
6480 // create needed triangles using n1,n2,n3 and inserted nodes
6481 int nbn = 2 + aNodesToInsert.size();
6482 //const SMDS_MeshNode* aNodes[nbn];
6483 vector<const SMDS_MeshNode*> aNodes(nbn);
6484 aNodes[0] = nodes[n1];
6485 aNodes[nbn-1] = nodes[n2];
6486 list<const SMDS_MeshNode*>::iterator nIt = aNodesToInsert.begin();
6487 for ( iNode = 1; nIt != aNodesToInsert.end(); nIt++ ) {
6488 aNodes[iNode++] = *nIt;
6490 for(i=1; i<nbn; i++) {
6491 SMDS_MeshElement* newElem =
6492 aMesh->AddFace(aNodes[i-1],aNodes[i],nodes[n3]);
6493 myLastCreatedElems.Append(newElem);
6494 if ( aShapeId && newElem )
6495 aMesh->SetMeshElementOnShape( newElem, aShapeId );
6497 // remove old quadratic face
6498 aMesh->RemoveElement(theFace);
6502 //=======================================================================
6503 //function : UpdateVolumes
6505 //=======================================================================
6506 void SMESH_MeshEditor::UpdateVolumes (const SMDS_MeshNode* theBetweenNode1,
6507 const SMDS_MeshNode* theBetweenNode2,
6508 list<const SMDS_MeshNode*>& theNodesToInsert)
6510 myLastCreatedElems.Clear();
6511 myLastCreatedNodes.Clear();
6513 SMDS_ElemIteratorPtr invElemIt = theBetweenNode1->GetInverseElementIterator(SMDSAbs_Volume);
6514 while (invElemIt->more()) { // loop on inverse elements of theBetweenNode1
6515 const SMDS_MeshElement* elem = invElemIt->next();
6517 // check, if current volume has link theBetweenNode1 - theBetweenNode2
6518 SMDS_VolumeTool aVolume (elem);
6519 if (!aVolume.IsLinked(theBetweenNode1, theBetweenNode2))
6522 // insert new nodes in all faces of the volume, sharing link theBetweenNode1 - theBetweenNode2
6523 int iface, nbFaces = aVolume.NbFaces();
6524 vector<const SMDS_MeshNode *> poly_nodes;
6525 vector<int> quantities (nbFaces);
6527 for (iface = 0; iface < nbFaces; iface++) {
6528 int nbFaceNodes = aVolume.NbFaceNodes(iface), nbInserted = 0;
6529 // faceNodes will contain (nbFaceNodes + 1) nodes, last = first
6530 const SMDS_MeshNode** faceNodes = aVolume.GetFaceNodes(iface);
6532 for (int inode = 0; inode < nbFaceNodes; inode++) {
6533 poly_nodes.push_back(faceNodes[inode]);
6535 if (nbInserted == 0) {
6536 if (faceNodes[inode] == theBetweenNode1) {
6537 if (faceNodes[inode + 1] == theBetweenNode2) {
6538 nbInserted = theNodesToInsert.size();
6540 // add nodes to insert
6541 list<const SMDS_MeshNode*>::iterator nIt = theNodesToInsert.begin();
6542 for (; nIt != theNodesToInsert.end(); nIt++) {
6543 poly_nodes.push_back(*nIt);
6547 else if (faceNodes[inode] == theBetweenNode2) {
6548 if (faceNodes[inode + 1] == theBetweenNode1) {
6549 nbInserted = theNodesToInsert.size();
6551 // add nodes to insert in reversed order
6552 list<const SMDS_MeshNode*>::iterator nIt = theNodesToInsert.end();
6554 for (; nIt != theNodesToInsert.begin(); nIt--) {
6555 poly_nodes.push_back(*nIt);
6557 poly_nodes.push_back(*nIt);
6564 quantities[iface] = nbFaceNodes + nbInserted;
6567 // Replace or update the volume
6568 SMESHDS_Mesh *aMesh = GetMeshDS();
6570 if (elem->IsPoly()) {
6571 aMesh->ChangePolyhedronNodes(elem, poly_nodes, quantities);
6575 int aShapeId = FindShape( elem );
6577 SMDS_MeshElement* newElem =
6578 aMesh->AddPolyhedralVolume(poly_nodes, quantities);
6579 myLastCreatedElems.Append(newElem);
6580 if (aShapeId && newElem)
6581 aMesh->SetMeshElementOnShape(newElem, aShapeId);
6583 aMesh->RemoveElement(elem);
6588 //=======================================================================
6590 * \brief Convert elements contained in a submesh to quadratic
6591 * \retval int - nb of checked elements
6593 //=======================================================================
6595 int SMESH_MeshEditor::convertElemToQuadratic(SMESHDS_SubMesh * theSm,
6596 SMESH_MesherHelper& theHelper,
6597 const bool theForce3d)
6600 if( !theSm ) return nbElem;
6601 SMDS_ElemIteratorPtr ElemItr = theSm->GetElements();
6602 while(ElemItr->more())
6605 const SMDS_MeshElement* elem = ElemItr->next();
6606 if( !elem || elem->IsQuadratic() ) continue;
6608 int id = elem->GetID();
6609 int nbNodes = elem->NbNodes();
6610 vector<const SMDS_MeshNode *> aNds (nbNodes);
6612 for(int i = 0; i < nbNodes; i++)
6614 aNds[i] = elem->GetNode(i);
6616 SMDSAbs_ElementType aType = elem->GetType();
6618 theSm->RemoveElement(elem);
6619 GetMeshDS()->SMDS_Mesh::RemoveFreeElement(elem);
6621 const SMDS_MeshElement* NewElem = 0;
6627 NewElem = theHelper.AddEdge(aNds[0], aNds[1], id, theForce3d);
6635 NewElem = theHelper.AddFace(aNds[0], aNds[1], aNds[2], id, theForce3d);
6638 NewElem = theHelper.AddFace(aNds[0], aNds[1], aNds[2], aNds[3], id, theForce3d);
6645 case SMDSAbs_Volume :
6650 NewElem = theHelper.AddVolume(aNds[0], aNds[1], aNds[2], aNds[3], id, true);
6653 NewElem = theHelper.AddVolume(aNds[0], aNds[1], aNds[2], aNds[3], aNds[4], aNds[5], id, true);
6656 NewElem = theHelper.AddVolume(aNds[0], aNds[1], aNds[2], aNds[3],
6657 aNds[4], aNds[5], aNds[6], aNds[7], id, true);
6669 AddToSameGroups( NewElem, elem, GetMeshDS());
6670 theSm->AddElement( NewElem );
6672 if ( NewElem != elem )
6673 RemoveElemFromGroups (elem, GetMeshDS());
6678 //=======================================================================
6679 //function : ConvertToQuadratic
6681 //=======================================================================
6682 void SMESH_MeshEditor::ConvertToQuadratic(const bool theForce3d)
6684 SMESHDS_Mesh* meshDS = GetMeshDS();
6686 SMESH_MesherHelper aHelper(*myMesh);
6687 aHelper.SetIsQuadratic( true );
6689 int nbCheckedElems = 0;
6690 if ( myMesh->HasShapeToMesh() )
6692 if ( SMESH_subMesh *aSubMesh = myMesh->GetSubMeshContaining(myMesh->GetShapeToMesh()))
6694 SMESH_subMeshIteratorPtr smIt = aSubMesh->getDependsOnIterator(true,false);
6695 while ( smIt->more() ) {
6696 SMESH_subMesh* sm = smIt->next();
6697 if ( SMESHDS_SubMesh *smDS = sm->GetSubMeshDS() ) {
6698 aHelper.SetSubShape( sm->GetSubShape() );
6699 nbCheckedElems += convertElemToQuadratic(smDS, aHelper, theForce3d);
6704 int totalNbElems = meshDS->NbEdges() + meshDS->NbFaces() + meshDS->NbVolumes();
6705 if ( nbCheckedElems < totalNbElems ) // not all elements in submeshes
6707 SMDS_EdgeIteratorPtr aEdgeItr = meshDS->edgesIterator();
6708 while(aEdgeItr->more())
6710 const SMDS_MeshEdge* edge = aEdgeItr->next();
6711 if(edge && !edge->IsQuadratic())
6713 int id = edge->GetID();
6714 const SMDS_MeshNode* n1 = edge->GetNode(0);
6715 const SMDS_MeshNode* n2 = edge->GetNode(1);
6717 meshDS->SMDS_Mesh::RemoveFreeElement(edge);
6719 const SMDS_MeshEdge* NewEdge = aHelper.AddEdge(n1, n2, id, theForce3d);
6721 AddToSameGroups(NewEdge, edge, meshDS);
6722 if ( NewEdge != edge )
6723 RemoveElemFromGroups (edge, meshDS);
6726 SMDS_FaceIteratorPtr aFaceItr = meshDS->facesIterator();
6727 while(aFaceItr->more())
6729 const SMDS_MeshFace* face = aFaceItr->next();
6730 if(!face || face->IsQuadratic() ) continue;
6732 int id = face->GetID();
6733 int nbNodes = face->NbNodes();
6734 vector<const SMDS_MeshNode *> aNds (nbNodes);
6736 for(int i = 0; i < nbNodes; i++)
6738 aNds[i] = face->GetNode(i);
6741 meshDS->SMDS_Mesh::RemoveFreeElement(face);
6743 SMDS_MeshFace * NewFace = 0;
6747 NewFace = aHelper.AddFace(aNds[0], aNds[1], aNds[2], id, theForce3d);
6750 NewFace = aHelper.AddFace(aNds[0], aNds[1], aNds[2], aNds[3], id, theForce3d);
6756 AddToSameGroups(NewFace, face, meshDS);
6757 if ( NewFace != face )
6758 RemoveElemFromGroups (face, meshDS);
6760 SMDS_VolumeIteratorPtr aVolumeItr = meshDS->volumesIterator();
6761 while(aVolumeItr->more())
6763 const SMDS_MeshVolume* volume = aVolumeItr->next();
6764 if(!volume || volume->IsQuadratic() ) continue;
6766 int id = volume->GetID();
6767 int nbNodes = volume->NbNodes();
6768 vector<const SMDS_MeshNode *> aNds (nbNodes);
6770 for(int i = 0; i < nbNodes; i++)
6772 aNds[i] = volume->GetNode(i);
6775 meshDS->SMDS_Mesh::RemoveFreeElement(volume);
6777 SMDS_MeshVolume * NewVolume = 0;
6781 NewVolume = aHelper.AddVolume(aNds[0], aNds[1], aNds[2],
6782 aNds[3], id, true );
6785 NewVolume = aHelper.AddVolume(aNds[0], aNds[1], aNds[2],
6786 aNds[3], aNds[4], aNds[5], id, true);
6789 NewVolume = aHelper.AddVolume(aNds[0], aNds[1], aNds[2], aNds[3],
6790 aNds[4], aNds[5], aNds[6], aNds[7], id, true);
6796 AddToSameGroups(NewVolume, volume, meshDS);
6797 if ( NewVolume != volume )
6798 RemoveElemFromGroups (volume, meshDS);
6803 //=======================================================================
6805 * \brief Convert quadratic elements to linear ones and remove quadratic nodes
6806 * \retval int - nb of checked elements
6808 //=======================================================================
6810 int SMESH_MeshEditor::removeQuadElem(SMESHDS_SubMesh * theSm,
6811 SMDS_ElemIteratorPtr theItr,
6812 const int theShapeID)
6815 SMESHDS_Mesh* meshDS = GetMeshDS();
6816 while( theItr->more() )
6818 const SMDS_MeshElement* elem = theItr->next();
6820 if( elem && elem->IsQuadratic())
6822 int id = elem->GetID();
6823 int nbNodes = elem->NbNodes();
6824 vector<const SMDS_MeshNode *> aNds, mediumNodes;
6825 aNds.reserve( nbNodes );
6826 mediumNodes.reserve( nbNodes );
6828 for(int i = 0; i < nbNodes; i++)
6830 const SMDS_MeshNode* n = elem->GetNode(i);
6832 if( elem->IsMediumNode( n ) )
6833 mediumNodes.push_back( n );
6835 aNds.push_back( n );
6837 if( aNds.empty() ) continue;
6838 SMDSAbs_ElementType aType = elem->GetType();
6840 //remove old quadratic element
6841 meshDS->SMDS_Mesh::RemoveFreeElement( elem );
6843 theSm->RemoveElement( elem );
6845 SMDS_MeshElement * NewElem = AddElement( aNds, aType, false, id );
6847 AddToSameGroups(NewElem, elem, meshDS);
6848 if ( NewElem != elem )
6849 RemoveElemFromGroups (elem, meshDS);
6850 if( theSm && NewElem )
6851 theSm->AddElement( NewElem );
6853 // remove medium nodes
6854 vector<const SMDS_MeshNode*>::iterator nIt = mediumNodes.begin();
6855 for ( ; nIt != mediumNodes.end(); ++nIt ) {
6856 const SMDS_MeshNode* n = *nIt;
6857 if ( n->NbInverseElements() == 0 ) {
6858 if ( n->GetPosition()->GetShapeId() != theShapeID )
6859 meshDS->RemoveFreeNode( n, meshDS->MeshElements
6860 ( n->GetPosition()->GetShapeId() ));
6862 meshDS->RemoveFreeNode( n, theSm );
6870 //=======================================================================
6871 //function : ConvertFromQuadratic
6873 //=======================================================================
6874 bool SMESH_MeshEditor::ConvertFromQuadratic()
6876 int nbCheckedElems = 0;
6877 if ( myMesh->HasShapeToMesh() )
6879 if ( SMESH_subMesh *aSubMesh = myMesh->GetSubMeshContaining(myMesh->GetShapeToMesh()))
6881 SMESH_subMeshIteratorPtr smIt = aSubMesh->getDependsOnIterator(true,false);
6882 while ( smIt->more() ) {
6883 SMESH_subMesh* sm = smIt->next();
6884 if ( SMESHDS_SubMesh *smDS = sm->GetSubMeshDS() )
6885 nbCheckedElems += removeQuadElem( smDS, smDS->GetElements(), sm->GetId() );
6891 GetMeshDS()->NbEdges() + GetMeshDS()->NbFaces() + GetMeshDS()->NbVolumes();
6892 if ( nbCheckedElems < totalNbElems ) // not all elements in submeshes
6894 SMESHDS_SubMesh *aSM = 0;
6895 removeQuadElem( aSM, GetMeshDS()->elementsIterator(), 0 );
6901 //=======================================================================
6902 //function : SewSideElements
6904 //=======================================================================
6906 SMESH_MeshEditor::Sew_Error
6907 SMESH_MeshEditor::SewSideElements (TIDSortedElemSet& theSide1,
6908 TIDSortedElemSet& theSide2,
6909 const SMDS_MeshNode* theFirstNode1,
6910 const SMDS_MeshNode* theFirstNode2,
6911 const SMDS_MeshNode* theSecondNode1,
6912 const SMDS_MeshNode* theSecondNode2)
6914 myLastCreatedElems.Clear();
6915 myLastCreatedNodes.Clear();
6917 MESSAGE ("::::SewSideElements()");
6918 if ( theSide1.size() != theSide2.size() )
6919 return SEW_DIFF_NB_OF_ELEMENTS;
6921 Sew_Error aResult = SEW_OK;
6923 // 1. Build set of faces representing each side
6924 // 2. Find which nodes of the side 1 to merge with ones on the side 2
6925 // 3. Replace nodes in elements of the side 1 and remove replaced nodes
6927 // =======================================================================
6928 // 1. Build set of faces representing each side:
6929 // =======================================================================
6930 // a. build set of nodes belonging to faces
6931 // b. complete set of faces: find missing fices whose nodes are in set of nodes
6932 // c. create temporary faces representing side of volumes if correspondent
6933 // face does not exist
6935 SMESHDS_Mesh* aMesh = GetMeshDS();
6936 SMDS_Mesh aTmpFacesMesh;
6937 set<const SMDS_MeshElement*> faceSet1, faceSet2;
6938 set<const SMDS_MeshElement*> volSet1, volSet2;
6939 set<const SMDS_MeshNode*> nodeSet1, nodeSet2;
6940 set<const SMDS_MeshElement*> * faceSetPtr[] = { &faceSet1, &faceSet2 };
6941 set<const SMDS_MeshElement*> * volSetPtr[] = { &volSet1, &volSet2 };
6942 set<const SMDS_MeshNode*> * nodeSetPtr[] = { &nodeSet1, &nodeSet2 };
6943 TIDSortedElemSet * elemSetPtr[] = { &theSide1, &theSide2 };
6944 int iSide, iFace, iNode;
6946 for ( iSide = 0; iSide < 2; iSide++ ) {
6947 set<const SMDS_MeshNode*> * nodeSet = nodeSetPtr[ iSide ];
6948 TIDSortedElemSet * elemSet = elemSetPtr[ iSide ];
6949 set<const SMDS_MeshElement*> * faceSet = faceSetPtr[ iSide ];
6950 set<const SMDS_MeshElement*> * volSet = volSetPtr [ iSide ];
6951 set<const SMDS_MeshElement*>::iterator vIt;
6952 TIDSortedElemSet::iterator eIt;
6953 set<const SMDS_MeshNode*>::iterator nIt;
6955 // check that given nodes belong to given elements
6956 const SMDS_MeshNode* n1 = ( iSide == 0 ) ? theFirstNode1 : theFirstNode2;
6957 const SMDS_MeshNode* n2 = ( iSide == 0 ) ? theSecondNode1 : theSecondNode2;
6958 int firstIndex = -1, secondIndex = -1;
6959 for (eIt = elemSet->begin(); eIt != elemSet->end(); eIt++ ) {
6960 const SMDS_MeshElement* elem = *eIt;
6961 if ( firstIndex < 0 ) firstIndex = elem->GetNodeIndex( n1 );
6962 if ( secondIndex < 0 ) secondIndex = elem->GetNodeIndex( n2 );
6963 if ( firstIndex > -1 && secondIndex > -1 ) break;
6965 if ( firstIndex < 0 || secondIndex < 0 ) {
6966 // we can simply return until temporary faces created
6967 return (iSide == 0 ) ? SEW_BAD_SIDE1_NODES : SEW_BAD_SIDE2_NODES;
6970 // -----------------------------------------------------------
6971 // 1a. Collect nodes of existing faces
6972 // and build set of face nodes in order to detect missing
6973 // faces corresponing to sides of volumes
6974 // -----------------------------------------------------------
6976 set< set <const SMDS_MeshNode*> > setOfFaceNodeSet;
6978 // loop on the given element of a side
6979 for (eIt = elemSet->begin(); eIt != elemSet->end(); eIt++ ) {
6980 //const SMDS_MeshElement* elem = *eIt;
6981 const SMDS_MeshElement* elem = *eIt;
6982 if ( elem->GetType() == SMDSAbs_Face ) {
6983 faceSet->insert( elem );
6984 set <const SMDS_MeshNode*> faceNodeSet;
6985 SMDS_ElemIteratorPtr nodeIt = elem->nodesIterator();
6986 while ( nodeIt->more() ) {
6987 const SMDS_MeshNode* n = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
6988 nodeSet->insert( n );
6989 faceNodeSet.insert( n );
6991 setOfFaceNodeSet.insert( faceNodeSet );
6993 else if ( elem->GetType() == SMDSAbs_Volume )
6994 volSet->insert( elem );
6996 // ------------------------------------------------------------------------------
6997 // 1b. Complete set of faces: find missing fices whose nodes are in set of nodes
6998 // ------------------------------------------------------------------------------
7000 for ( nIt = nodeSet->begin(); nIt != nodeSet->end(); nIt++ ) { // loop on nodes of iSide
7001 SMDS_ElemIteratorPtr fIt = (*nIt)->GetInverseElementIterator(SMDSAbs_Face);
7002 while ( fIt->more() ) { // loop on faces sharing a node
7003 const SMDS_MeshElement* f = fIt->next();
7004 if ( faceSet->find( f ) == faceSet->end() ) {
7005 // check if all nodes are in nodeSet and
7006 // complete setOfFaceNodeSet if they are
7007 set <const SMDS_MeshNode*> faceNodeSet;
7008 SMDS_ElemIteratorPtr nodeIt = f->nodesIterator();
7009 bool allInSet = true;
7010 while ( nodeIt->more() && allInSet ) { // loop on nodes of a face
7011 const SMDS_MeshNode* n = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
7012 if ( nodeSet->find( n ) == nodeSet->end() )
7015 faceNodeSet.insert( n );
7018 faceSet->insert( f );
7019 setOfFaceNodeSet.insert( faceNodeSet );
7025 // -------------------------------------------------------------------------
7026 // 1c. Create temporary faces representing sides of volumes if correspondent
7027 // face does not exist
7028 // -------------------------------------------------------------------------
7030 if ( !volSet->empty() ) {
7031 //int nodeSetSize = nodeSet->size();
7033 // loop on given volumes
7034 for ( vIt = volSet->begin(); vIt != volSet->end(); vIt++ ) {
7035 SMDS_VolumeTool vol (*vIt);
7036 // loop on volume faces: find free faces
7037 // --------------------------------------
7038 list<const SMDS_MeshElement* > freeFaceList;
7039 for ( iFace = 0; iFace < vol.NbFaces(); iFace++ ) {
7040 if ( !vol.IsFreeFace( iFace ))
7042 // check if there is already a face with same nodes in a face set
7043 const SMDS_MeshElement* aFreeFace = 0;
7044 const SMDS_MeshNode** fNodes = vol.GetFaceNodes( iFace );
7045 int nbNodes = vol.NbFaceNodes( iFace );
7046 set <const SMDS_MeshNode*> faceNodeSet;
7047 vol.GetFaceNodes( iFace, faceNodeSet );
7048 bool isNewFace = setOfFaceNodeSet.insert( faceNodeSet ).second;
7050 // no such a face is given but it still can exist, check it
7051 if ( nbNodes == 3 ) {
7052 aFreeFace = aMesh->FindFace( fNodes[0],fNodes[1],fNodes[2] );
7054 else if ( nbNodes == 4 ) {
7055 aFreeFace = aMesh->FindFace( fNodes[0],fNodes[1],fNodes[2],fNodes[3] );
7058 vector<const SMDS_MeshNode *> poly_nodes ( fNodes, & fNodes[nbNodes]);
7059 aFreeFace = aMesh->FindFace(poly_nodes);
7063 // create a temporary face
7064 if ( nbNodes == 3 ) {
7065 aFreeFace = aTmpFacesMesh.AddFace( fNodes[0],fNodes[1],fNodes[2] );
7067 else if ( nbNodes == 4 ) {
7068 aFreeFace = aTmpFacesMesh.AddFace( fNodes[0],fNodes[1],fNodes[2],fNodes[3] );
7071 vector<const SMDS_MeshNode *> poly_nodes ( fNodes, & fNodes[nbNodes]);
7072 aFreeFace = aTmpFacesMesh.AddPolygonalFace(poly_nodes);
7076 freeFaceList.push_back( aFreeFace );
7078 } // loop on faces of a volume
7080 // choose one of several free faces
7081 // --------------------------------------
7082 if ( freeFaceList.size() > 1 ) {
7083 // choose a face having max nb of nodes shared by other elems of a side
7084 int maxNbNodes = -1/*, nbExcludedFaces = 0*/;
7085 list<const SMDS_MeshElement* >::iterator fIt = freeFaceList.begin();
7086 while ( fIt != freeFaceList.end() ) { // loop on free faces
7087 int nbSharedNodes = 0;
7088 SMDS_ElemIteratorPtr nodeIt = (*fIt)->nodesIterator();
7089 while ( nodeIt->more() ) { // loop on free face nodes
7090 const SMDS_MeshNode* n =
7091 static_cast<const SMDS_MeshNode*>( nodeIt->next() );
7092 SMDS_ElemIteratorPtr invElemIt = n->GetInverseElementIterator();
7093 while ( invElemIt->more() ) {
7094 const SMDS_MeshElement* e = invElemIt->next();
7095 if ( faceSet->find( e ) != faceSet->end() )
7097 if ( elemSet->find( e ) != elemSet->end() )
7101 if ( nbSharedNodes >= maxNbNodes ) {
7102 maxNbNodes = nbSharedNodes;
7106 freeFaceList.erase( fIt++ ); // here fIt++ occures before erase
7108 if ( freeFaceList.size() > 1 )
7110 // could not choose one face, use another way
7111 // choose a face most close to the bary center of the opposite side
7112 gp_XYZ aBC( 0., 0., 0. );
7113 set <const SMDS_MeshNode*> addedNodes;
7114 TIDSortedElemSet * elemSet2 = elemSetPtr[ 1 - iSide ];
7115 eIt = elemSet2->begin();
7116 for ( eIt = elemSet2->begin(); eIt != elemSet2->end(); eIt++ ) {
7117 SMDS_ElemIteratorPtr nodeIt = (*eIt)->nodesIterator();
7118 while ( nodeIt->more() ) { // loop on free face nodes
7119 const SMDS_MeshNode* n =
7120 static_cast<const SMDS_MeshNode*>( nodeIt->next() );
7121 if ( addedNodes.insert( n ).second )
7122 aBC += gp_XYZ( n->X(),n->Y(),n->Z() );
7125 aBC /= addedNodes.size();
7126 double minDist = DBL_MAX;
7127 fIt = freeFaceList.begin();
7128 while ( fIt != freeFaceList.end() ) { // loop on free faces
7130 SMDS_ElemIteratorPtr nodeIt = (*fIt)->nodesIterator();
7131 while ( nodeIt->more() ) { // loop on free face nodes
7132 const SMDS_MeshNode* n =
7133 static_cast<const SMDS_MeshNode*>( nodeIt->next() );
7134 gp_XYZ p( n->X(),n->Y(),n->Z() );
7135 dist += ( aBC - p ).SquareModulus();
7137 if ( dist < minDist ) {
7139 freeFaceList.erase( freeFaceList.begin(), fIt++ );
7142 fIt = freeFaceList.erase( fIt++ );
7145 } // choose one of several free faces of a volume
7147 if ( freeFaceList.size() == 1 ) {
7148 const SMDS_MeshElement* aFreeFace = freeFaceList.front();
7149 faceSet->insert( aFreeFace );
7150 // complete a node set with nodes of a found free face
7151 // for ( iNode = 0; iNode < ; iNode++ )
7152 // nodeSet->insert( fNodes[ iNode ] );
7155 } // loop on volumes of a side
7157 // // complete a set of faces if new nodes in a nodeSet appeared
7158 // // ----------------------------------------------------------
7159 // if ( nodeSetSize != nodeSet->size() ) {
7160 // for ( ; nIt != nodeSet->end(); nIt++ ) { // loop on nodes of iSide
7161 // SMDS_ElemIteratorPtr fIt = (*nIt)->GetInverseElementIterator(SMDSAbs_Face);
7162 // while ( fIt->more() ) { // loop on faces sharing a node
7163 // const SMDS_MeshElement* f = fIt->next();
7164 // if ( faceSet->find( f ) == faceSet->end() ) {
7165 // // check if all nodes are in nodeSet and
7166 // // complete setOfFaceNodeSet if they are
7167 // set <const SMDS_MeshNode*> faceNodeSet;
7168 // SMDS_ElemIteratorPtr nodeIt = f->nodesIterator();
7169 // bool allInSet = true;
7170 // while ( nodeIt->more() && allInSet ) { // loop on nodes of a face
7171 // const SMDS_MeshNode* n = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
7172 // if ( nodeSet->find( n ) == nodeSet->end() )
7173 // allInSet = false;
7175 // faceNodeSet.insert( n );
7177 // if ( allInSet ) {
7178 // faceSet->insert( f );
7179 // setOfFaceNodeSet.insert( faceNodeSet );
7185 } // Create temporary faces, if there are volumes given
7188 if ( faceSet1.size() != faceSet2.size() ) {
7189 // delete temporary faces: they are in reverseElements of actual nodes
7190 SMDS_FaceIteratorPtr tmpFaceIt = aTmpFacesMesh.facesIterator();
7191 while ( tmpFaceIt->more() )
7192 aTmpFacesMesh.RemoveElement( tmpFaceIt->next() );
7193 MESSAGE("Diff nb of faces");
7194 return SEW_TOPO_DIFF_SETS_OF_ELEMENTS;
7197 // ============================================================
7198 // 2. Find nodes to merge:
7199 // bind a node to remove to a node to put instead
7200 // ============================================================
7202 TNodeNodeMap nReplaceMap; // bind a node to remove to a node to put instead
7203 if ( theFirstNode1 != theFirstNode2 )
7204 nReplaceMap.insert( TNodeNodeMap::value_type( theFirstNode1, theFirstNode2 ));
7205 if ( theSecondNode1 != theSecondNode2 )
7206 nReplaceMap.insert( TNodeNodeMap::value_type( theSecondNode1, theSecondNode2 ));
7208 LinkID_Gen aLinkID_Gen( GetMeshDS() );
7209 set< long > linkIdSet; // links to process
7210 linkIdSet.insert( aLinkID_Gen.GetLinkID( theFirstNode1, theSecondNode1 ));
7212 typedef pair< const SMDS_MeshNode*, const SMDS_MeshNode* > NLink;
7213 list< NLink > linkList[2];
7214 linkList[0].push_back( NLink( theFirstNode1, theSecondNode1 ));
7215 linkList[1].push_back( NLink( theFirstNode2, theSecondNode2 ));
7216 // loop on links in linkList; find faces by links and append links
7217 // of the found faces to linkList
7218 list< NLink >::iterator linkIt[] = { linkList[0].begin(), linkList[1].begin() } ;
7219 for ( ; linkIt[0] != linkList[0].end(); linkIt[0]++, linkIt[1]++ ) {
7220 NLink link[] = { *linkIt[0], *linkIt[1] };
7221 long linkID = aLinkID_Gen.GetLinkID( link[0].first, link[0].second );
7222 if ( linkIdSet.find( linkID ) == linkIdSet.end() )
7225 // by links, find faces in the face sets,
7226 // and find indices of link nodes in the found faces;
7227 // in a face set, there is only one or no face sharing a link
7228 // ---------------------------------------------------------------
7230 const SMDS_MeshElement* face[] = { 0, 0 };
7231 //const SMDS_MeshNode* faceNodes[ 2 ][ 5 ];
7232 vector<const SMDS_MeshNode*> fnodes1(9);
7233 vector<const SMDS_MeshNode*> fnodes2(9);
7234 //const SMDS_MeshNode* notLinkNodes[ 2 ][ 2 ] = {{ 0, 0 },{ 0, 0 }} ;
7235 vector<const SMDS_MeshNode*> notLinkNodes1(6);
7236 vector<const SMDS_MeshNode*> notLinkNodes2(6);
7237 int iLinkNode[2][2];
7238 for ( iSide = 0; iSide < 2; iSide++ ) { // loop on 2 sides
7239 const SMDS_MeshNode* n1 = link[iSide].first;
7240 const SMDS_MeshNode* n2 = link[iSide].second;
7241 set<const SMDS_MeshElement*> * faceSet = faceSetPtr[ iSide ];
7242 set< const SMDS_MeshElement* > fMap;
7243 for ( int i = 0; i < 2; i++ ) { // loop on 2 nodes of a link
7244 const SMDS_MeshNode* n = i ? n1 : n2; // a node of a link
7245 SMDS_ElemIteratorPtr fIt = n->GetInverseElementIterator(SMDSAbs_Face);
7246 while ( fIt->more() ) { // loop on faces sharing a node
7247 const SMDS_MeshElement* f = fIt->next();
7248 if (faceSet->find( f ) != faceSet->end() && // f is in face set
7249 ! fMap.insert( f ).second ) // f encounters twice
7251 if ( face[ iSide ] ) {
7252 MESSAGE( "2 faces per link " );
7253 aResult = iSide ? SEW_BAD_SIDE2_NODES : SEW_BAD_SIDE1_NODES;
7257 faceSet->erase( f );
7258 // get face nodes and find ones of a link
7263 fnodes1.resize(f->NbNodes()+1);
7264 notLinkNodes1.resize(f->NbNodes()-2);
7267 fnodes2.resize(f->NbNodes()+1);
7268 notLinkNodes2.resize(f->NbNodes()-2);
7271 if(!f->IsQuadratic()) {
7272 SMDS_ElemIteratorPtr nIt = f->nodesIterator();
7273 while ( nIt->more() ) {
7274 const SMDS_MeshNode* n =
7275 static_cast<const SMDS_MeshNode*>( nIt->next() );
7277 iLinkNode[ iSide ][ 0 ] = iNode;
7279 else if ( n == n2 ) {
7280 iLinkNode[ iSide ][ 1 ] = iNode;
7282 //else if ( notLinkNodes[ iSide ][ 0 ] )
7283 // notLinkNodes[ iSide ][ 1 ] = n;
7285 // notLinkNodes[ iSide ][ 0 ] = n;
7289 notLinkNodes1[nbl] = n;
7290 //notLinkNodes1.push_back(n);
7292 notLinkNodes2[nbl] = n;
7293 //notLinkNodes2.push_back(n);
7295 //faceNodes[ iSide ][ iNode++ ] = n;
7297 fnodes1[iNode++] = n;
7300 fnodes2[iNode++] = n;
7304 else { // f->IsQuadratic()
7305 const SMDS_QuadraticFaceOfNodes* F =
7306 static_cast<const SMDS_QuadraticFaceOfNodes*>(f);
7307 // use special nodes iterator
7308 SMDS_NodeIteratorPtr anIter = F->interlacedNodesIterator();
7309 while ( anIter->more() ) {
7310 const SMDS_MeshNode* n =
7311 static_cast<const SMDS_MeshNode*>( anIter->next() );
7313 iLinkNode[ iSide ][ 0 ] = iNode;
7315 else if ( n == n2 ) {
7316 iLinkNode[ iSide ][ 1 ] = iNode;
7321 notLinkNodes1[nbl] = n;
7324 notLinkNodes2[nbl] = n;
7328 fnodes1[iNode++] = n;
7331 fnodes2[iNode++] = n;
7335 //faceNodes[ iSide ][ iNode ] = faceNodes[ iSide ][ 0 ];
7337 fnodes1[iNode] = fnodes1[0];
7340 fnodes2[iNode] = fnodes1[0];
7347 // check similarity of elements of the sides
7348 if (aResult == SEW_OK && ( face[0] && !face[1] ) || ( !face[0] && face[1] )) {
7349 MESSAGE("Correspondent face not found on side " << ( face[0] ? 1 : 0 ));
7350 if ( nReplaceMap.size() == 2 ) { // faces on input nodes not found
7351 aResult = ( face[0] ? SEW_BAD_SIDE2_NODES : SEW_BAD_SIDE1_NODES );
7354 aResult = SEW_TOPO_DIFF_SETS_OF_ELEMENTS;
7356 break; // do not return because it s necessary to remove tmp faces
7359 // set nodes to merge
7360 // -------------------
7362 if ( face[0] && face[1] ) {
7363 int nbNodes = face[0]->NbNodes();
7364 if ( nbNodes != face[1]->NbNodes() ) {
7365 MESSAGE("Diff nb of face nodes");
7366 aResult = SEW_TOPO_DIFF_SETS_OF_ELEMENTS;
7367 break; // do not return because it s necessary to remove tmp faces
7369 bool reverse[] = { false, false }; // order of notLinkNodes of quadrangle
7370 if ( nbNodes == 3 ) {
7371 //nReplaceMap.insert( TNodeNodeMap::value_type
7372 // ( notLinkNodes[0][0], notLinkNodes[1][0] ));
7373 nReplaceMap.insert( TNodeNodeMap::value_type
7374 ( notLinkNodes1[0], notLinkNodes2[0] ));
7377 for ( iSide = 0; iSide < 2; iSide++ ) { // loop on 2 sides
7378 // analyse link orientation in faces
7379 int i1 = iLinkNode[ iSide ][ 0 ];
7380 int i2 = iLinkNode[ iSide ][ 1 ];
7381 reverse[ iSide ] = Abs( i1 - i2 ) == 1 ? i1 > i2 : i2 > i1;
7382 // if notLinkNodes are the first and the last ones, then
7383 // their order does not correspond to the link orientation
7384 if (( i1 == 1 && i2 == 2 ) ||
7385 ( i1 == 2 && i2 == 1 ))
7386 reverse[ iSide ] = !reverse[ iSide ];
7388 if ( reverse[0] == reverse[1] ) {
7389 //nReplaceMap.insert( TNodeNodeMap::value_type
7390 // ( notLinkNodes[0][0], notLinkNodes[1][0] ));
7391 //nReplaceMap.insert( TNodeNodeMap::value_type
7392 // ( notLinkNodes[0][1], notLinkNodes[1][1] ));
7393 for(int nn=0; nn<nbNodes-2; nn++) {
7394 nReplaceMap.insert( TNodeNodeMap::value_type
7395 ( notLinkNodes1[nn], notLinkNodes2[nn] ));
7399 //nReplaceMap.insert( TNodeNodeMap::value_type
7400 // ( notLinkNodes[0][0], notLinkNodes[1][1] ));
7401 //nReplaceMap.insert( TNodeNodeMap::value_type
7402 // ( notLinkNodes[0][1], notLinkNodes[1][0] ));
7403 for(int nn=0; nn<nbNodes-2; nn++) {
7404 nReplaceMap.insert( TNodeNodeMap::value_type
7405 ( notLinkNodes1[nn], notLinkNodes2[nbNodes-3-nn] ));
7410 // add other links of the faces to linkList
7411 // -----------------------------------------
7413 //const SMDS_MeshNode** nodes = faceNodes[ 0 ];
7414 for ( iNode = 0; iNode < nbNodes; iNode++ ) {
7415 //linkID = aLinkID_Gen.GetLinkID( nodes[iNode], nodes[iNode+1] );
7416 linkID = aLinkID_Gen.GetLinkID( fnodes1[iNode], fnodes1[iNode+1] );
7417 pair< set<long>::iterator, bool > iter_isnew = linkIdSet.insert( linkID );
7418 if ( !iter_isnew.second ) { // already in a set: no need to process
7419 linkIdSet.erase( iter_isnew.first );
7421 else // new in set == encountered for the first time: add
7423 //const SMDS_MeshNode* n1 = nodes[ iNode ];
7424 //const SMDS_MeshNode* n2 = nodes[ iNode + 1];
7425 const SMDS_MeshNode* n1 = fnodes1[ iNode ];
7426 const SMDS_MeshNode* n2 = fnodes1[ iNode + 1];
7427 linkList[0].push_back ( NLink( n1, n2 ));
7428 linkList[1].push_back ( NLink( nReplaceMap[n1], nReplaceMap[n2] ));
7432 } // loop on link lists
7434 if ( aResult == SEW_OK &&
7435 ( linkIt[0] != linkList[0].end() ||
7436 !faceSetPtr[0]->empty() || !faceSetPtr[1]->empty() )) {
7437 MESSAGE( (linkIt[0] != linkList[0].end()) <<" "<< (faceSetPtr[0]->empty()) <<
7438 " " << (faceSetPtr[1]->empty()));
7439 aResult = SEW_TOPO_DIFF_SETS_OF_ELEMENTS;
7442 // ====================================================================
7443 // 3. Replace nodes in elements of the side 1 and remove replaced nodes
7444 // ====================================================================
7446 // delete temporary faces: they are in reverseElements of actual nodes
7447 SMDS_FaceIteratorPtr tmpFaceIt = aTmpFacesMesh.facesIterator();
7448 while ( tmpFaceIt->more() )
7449 aTmpFacesMesh.RemoveElement( tmpFaceIt->next() );
7451 if ( aResult != SEW_OK)
7454 list< int > nodeIDsToRemove/*, elemIDsToRemove*/;
7455 // loop on nodes replacement map
7456 TNodeNodeMap::iterator nReplaceMapIt = nReplaceMap.begin(), nnIt;
7457 for ( ; nReplaceMapIt != nReplaceMap.end(); nReplaceMapIt++ )
7458 if ( (*nReplaceMapIt).first != (*nReplaceMapIt).second ) {
7459 const SMDS_MeshNode* nToRemove = (*nReplaceMapIt).first;
7460 nodeIDsToRemove.push_back( nToRemove->GetID() );
7461 // loop on elements sharing nToRemove
7462 SMDS_ElemIteratorPtr invElemIt = nToRemove->GetInverseElementIterator();
7463 while ( invElemIt->more() ) {
7464 const SMDS_MeshElement* e = invElemIt->next();
7465 // get a new suite of nodes: make replacement
7466 int nbReplaced = 0, i = 0, nbNodes = e->NbNodes();
7467 vector< const SMDS_MeshNode*> nodes( nbNodes );
7468 SMDS_ElemIteratorPtr nIt = e->nodesIterator();
7469 while ( nIt->more() ) {
7470 const SMDS_MeshNode* n =
7471 static_cast<const SMDS_MeshNode*>( nIt->next() );
7472 nnIt = nReplaceMap.find( n );
7473 if ( nnIt != nReplaceMap.end() ) {
7479 // if ( nbReplaced == nbNodes && e->GetType() == SMDSAbs_Face )
7480 // elemIDsToRemove.push_back( e->GetID() );
7483 aMesh->ChangeElementNodes( e, & nodes[0], nbNodes );
7487 Remove( nodeIDsToRemove, true );
7492 //================================================================================
7494 * \brief Find corresponding nodes in two sets of faces
7495 * \param theSide1 - first face set
7496 * \param theSide2 - second first face
7497 * \param theFirstNode1 - a boundary node of set 1
7498 * \param theFirstNode2 - a node of set 2 corresponding to theFirstNode1
7499 * \param theSecondNode1 - a boundary node of set 1 linked with theFirstNode1
7500 * \param theSecondNode2 - a node of set 2 corresponding to theSecondNode1
7501 * \param nReplaceMap - output map of corresponding nodes
7502 * \retval bool - is a success or not
7504 //================================================================================
7507 //#define DEBUG_MATCHING_NODES
7510 SMESH_MeshEditor::Sew_Error
7511 SMESH_MeshEditor::FindMatchingNodes(set<const SMDS_MeshElement*>& theSide1,
7512 set<const SMDS_MeshElement*>& theSide2,
7513 const SMDS_MeshNode* theFirstNode1,
7514 const SMDS_MeshNode* theFirstNode2,
7515 const SMDS_MeshNode* theSecondNode1,
7516 const SMDS_MeshNode* theSecondNode2,
7517 TNodeNodeMap & nReplaceMap)
7519 set<const SMDS_MeshElement*> * faceSetPtr[] = { &theSide1, &theSide2 };
7521 nReplaceMap.clear();
7522 if ( theFirstNode1 != theFirstNode2 )
7523 nReplaceMap.insert( make_pair( theFirstNode1, theFirstNode2 ));
7524 if ( theSecondNode1 != theSecondNode2 )
7525 nReplaceMap.insert( make_pair( theSecondNode1, theSecondNode2 ));
7527 set< SMESH_TLink > linkSet; // set of nodes where order of nodes is ignored
7528 linkSet.insert( SMESH_TLink( theFirstNode1, theSecondNode1 ));
7530 list< NLink > linkList[2];
7531 linkList[0].push_back( NLink( theFirstNode1, theSecondNode1 ));
7532 linkList[1].push_back( NLink( theFirstNode2, theSecondNode2 ));
7534 // loop on links in linkList; find faces by links and append links
7535 // of the found faces to linkList
7536 list< NLink >::iterator linkIt[] = { linkList[0].begin(), linkList[1].begin() } ;
7537 for ( ; linkIt[0] != linkList[0].end(); linkIt[0]++, linkIt[1]++ ) {
7538 NLink link[] = { *linkIt[0], *linkIt[1] };
7539 if ( linkSet.find( link[0] ) == linkSet.end() )
7542 // by links, find faces in the face sets,
7543 // and find indices of link nodes in the found faces;
7544 // in a face set, there is only one or no face sharing a link
7545 // ---------------------------------------------------------------
7547 const SMDS_MeshElement* face[] = { 0, 0 };
7548 list<const SMDS_MeshNode*> notLinkNodes[2];
7549 //bool reverse[] = { false, false }; // order of notLinkNodes
7551 for ( int iSide = 0; iSide < 2; iSide++ ) // loop on 2 sides
7553 const SMDS_MeshNode* n1 = link[iSide].first;
7554 const SMDS_MeshNode* n2 = link[iSide].second;
7555 set<const SMDS_MeshElement*> * faceSet = faceSetPtr[ iSide ];
7556 set< const SMDS_MeshElement* > facesOfNode1;
7557 for ( int iNode = 0; iNode < 2; iNode++ ) // loop on 2 nodes of a link
7559 // during a loop of the first node, we find all faces around n1,
7560 // during a loop of the second node, we find one face sharing both n1 and n2
7561 const SMDS_MeshNode* n = iNode ? n1 : n2; // a node of a link
7562 SMDS_ElemIteratorPtr fIt = n->GetInverseElementIterator(SMDSAbs_Face);
7563 while ( fIt->more() ) { // loop on faces sharing a node
7564 const SMDS_MeshElement* f = fIt->next();
7565 if (faceSet->find( f ) != faceSet->end() && // f is in face set
7566 ! facesOfNode1.insert( f ).second ) // f encounters twice
7568 if ( face[ iSide ] ) {
7569 MESSAGE( "2 faces per link " );
7570 return ( iSide ? SEW_BAD_SIDE2_NODES : SEW_BAD_SIDE1_NODES );
7573 faceSet->erase( f );
7575 // get not link nodes
7576 int nbN = f->NbNodes();
7577 if ( f->IsQuadratic() )
7579 nbNodes[ iSide ] = nbN;
7580 list< const SMDS_MeshNode* > & nodes = notLinkNodes[ iSide ];
7581 int i1 = f->GetNodeIndex( n1 );
7582 int i2 = f->GetNodeIndex( n2 );
7583 int iEnd = nbN, iBeg = -1, iDelta = 1;
7584 bool reverse = ( Abs( i1 - i2 ) == 1 ? i1 > i2 : i2 > i1 );
7586 std::swap( iEnd, iBeg ); iDelta = -1;
7591 if ( i == iEnd ) i = iBeg + iDelta;
7592 if ( i == i1 ) break;
7593 nodes.push_back ( f->GetNode( i ) );
7599 // check similarity of elements of the sides
7600 if (( face[0] && !face[1] ) || ( !face[0] && face[1] )) {
7601 MESSAGE("Correspondent face not found on side " << ( face[0] ? 1 : 0 ));
7602 if ( nReplaceMap.size() == 2 ) { // faces on input nodes not found
7603 return ( face[0] ? SEW_BAD_SIDE2_NODES : SEW_BAD_SIDE1_NODES );
7606 return SEW_TOPO_DIFF_SETS_OF_ELEMENTS;
7610 // set nodes to merge
7611 // -------------------
7613 if ( face[0] && face[1] ) {
7614 if ( nbNodes[0] != nbNodes[1] ) {
7615 MESSAGE("Diff nb of face nodes");
7616 return SEW_TOPO_DIFF_SETS_OF_ELEMENTS;
7618 #ifdef DEBUG_MATCHING_NODES
7619 MESSAGE ( " Link 1: " << link[0].first->GetID() <<" "<< link[0].second->GetID()
7620 << " F 1: " << face[0] << "| Link 2: " << link[1].first->GetID() <<" "
7621 << link[1].second->GetID() << " F 2: " << face[1] << " | Bind: " ) ;
7623 int nbN = nbNodes[0];
7625 list<const SMDS_MeshNode*>::iterator n1 = notLinkNodes[0].begin();
7626 list<const SMDS_MeshNode*>::iterator n2 = notLinkNodes[1].begin();
7627 for ( int i = 0 ; i < nbN - 2; ++i ) {
7628 #ifdef DEBUG_MATCHING_NODES
7629 MESSAGE ( (*n1)->GetID() << " to " << (*n2)->GetID() );
7631 nReplaceMap.insert( make_pair( *(n1++), *(n2++) ));
7635 // add other links of the face 1 to linkList
7636 // -----------------------------------------
7638 const SMDS_MeshElement* f0 = face[0];
7639 const SMDS_MeshNode* n1 = f0->GetNode( nbN - 1 );
7640 for ( int i = 0; i < nbN; i++ )
7642 const SMDS_MeshNode* n2 = f0->GetNode( i );
7643 pair< set< SMESH_TLink >::iterator, bool > iter_isnew =
7644 linkSet.insert( SMESH_TLink( n1, n2 ));
7645 if ( !iter_isnew.second ) { // already in a set: no need to process
7646 linkSet.erase( iter_isnew.first );
7648 else // new in set == encountered for the first time: add
7650 #ifdef DEBUG_MATCHING_NODES
7651 MESSAGE ( "Add link 1: " << n1->GetID() << " " << n2->GetID() << " "
7652 << " | link 2: " << nReplaceMap[n1]->GetID() << " " << nReplaceMap[n2]->GetID() << " " );
7654 linkList[0].push_back ( NLink( n1, n2 ));
7655 linkList[1].push_back ( NLink( nReplaceMap[n1], nReplaceMap[n2] ));
7660 } // loop on link lists
7666 \brief Creates a hole in a mesh by doubling the nodes of some particular elements
7667 \param theNodes - identifiers of nodes to be doubled
7668 \param theModifiedElems - identifiers of elements to be updated by the new (doubled)
7669 nodes. If list of element identifiers is empty then nodes are doubled but
7670 they not assigned to elements
7671 \return TRUE if operation has been completed successfully, FALSE otherwise
7673 bool SMESH_MeshEditor::DoubleNodes( const std::list< int >& theListOfNodes,
7674 const std::list< int >& theListOfModifiedElems )
7676 myLastCreatedElems.Clear();
7677 myLastCreatedNodes.Clear();
7679 if ( theListOfNodes.size() == 0 )
7682 SMESHDS_Mesh* aMeshDS = GetMeshDS();
7686 // iterate through nodes and duplicate them
7688 std::map< const SMDS_MeshNode*, const SMDS_MeshNode* > anOldNodeToNewNode;
7690 std::list< int >::const_iterator aNodeIter;
7691 for ( aNodeIter = theListOfNodes.begin(); aNodeIter != theListOfNodes.end(); ++aNodeIter )
7693 int aCurr = *aNodeIter;
7694 SMDS_MeshNode* aNode = (SMDS_MeshNode*)aMeshDS->FindNode( aCurr );
7700 const SMDS_MeshNode* aNewNode = aMeshDS->AddNode( aNode->X(), aNode->Y(), aNode->Z() );
7703 anOldNodeToNewNode[ aNode ] = aNewNode;
7704 myLastCreatedNodes.Append( aNewNode );
7708 // Create map of new nodes for modified elements
7710 std::map< SMDS_MeshElement*, vector<const SMDS_MeshNode*> > anElemToNodes;
7712 std::list< int >::const_iterator anElemIter;
7713 for ( anElemIter = theListOfModifiedElems.begin();
7714 anElemIter != theListOfModifiedElems.end(); ++anElemIter )
7716 int aCurr = *anElemIter;
7717 SMDS_MeshElement* anElem = (SMDS_MeshElement*)aMeshDS->FindElement( aCurr );
7721 vector<const SMDS_MeshNode*> aNodeArr( anElem->NbNodes() );
7723 SMDS_ElemIteratorPtr anIter = anElem->nodesIterator();
7725 while ( anIter->more() )
7727 SMDS_MeshNode* aCurrNode = (SMDS_MeshNode*)anIter->next();
7728 if ( aCurr && anOldNodeToNewNode.find( aCurrNode ) != anOldNodeToNewNode.end() )
7730 const SMDS_MeshNode* aNewNode = anOldNodeToNewNode[ aCurrNode ];
7731 aNodeArr[ ind++ ] = aNewNode;
7734 aNodeArr[ ind++ ] = aCurrNode;
7736 anElemToNodes[ anElem ] = aNodeArr;
7739 // Change nodes of elements
7741 std::map< SMDS_MeshElement*, vector<const SMDS_MeshNode*> >::iterator
7742 anElemToNodesIter = anElemToNodes.begin();
7743 for ( ; anElemToNodesIter != anElemToNodes.end(); ++anElemToNodesIter )
7745 const SMDS_MeshElement* anElem = anElemToNodesIter->first;
7746 vector<const SMDS_MeshNode*> aNodeArr = anElemToNodesIter->second;
7748 aMeshDS->ChangeElementNodes( anElem, &aNodeArr[ 0 ], anElem->NbNodes() );