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 typedef pair< const SMDS_MeshNode*, const SMDS_MeshNode* > NLink;
97 //=======================================================================
99 * \brief A sorted pair of nodes
101 //=======================================================================
103 struct TLink: public NLink
105 TLink(const SMDS_MeshNode* n1, const SMDS_MeshNode* n2 ):NLink( n1, n2 )
106 { if ( n1->GetID() < n2->GetID() ) std::swap( first, second ); }
107 TLink(const NLink& link ):NLink( link )
108 { if ( first->GetID() < second->GetID() ) std::swap( first, second ); }
111 //=======================================================================
112 //function : SMESH_MeshEditor
114 //=======================================================================
116 SMESH_MeshEditor::SMESH_MeshEditor( SMESH_Mesh* theMesh )
117 :myMesh( theMesh ) // theMesh may be NULL
121 //=======================================================================
125 //=======================================================================
128 SMESH_MeshEditor::AddElement(const vector<const SMDS_MeshNode*> & node,
129 const SMDSAbs_ElementType type,
133 SMDS_MeshElement* e = 0;
134 int nbnode = node.size();
135 SMESHDS_Mesh* mesh = GetMeshDS();
139 if ( ID ) e = mesh->AddEdgeWithID(node[0], node[1], ID);
140 else e = mesh->AddEdge (node[0], node[1] );
141 else if ( nbnode == 3 )
142 if ( ID ) e = mesh->AddEdgeWithID(node[0], node[1], node[2], ID);
143 else e = mesh->AddEdge (node[0], node[1], node[2] );
148 if ( ID ) e = mesh->AddFaceWithID(node[0], node[1], node[2], ID);
149 else e = mesh->AddFace (node[0], node[1], node[2] );
150 else if (nbnode == 4)
151 if ( ID ) e = mesh->AddFaceWithID(node[0], node[1], node[2], node[3], ID);
152 else e = mesh->AddFace (node[0], node[1], node[2], node[3] );
153 else if (nbnode == 6)
154 if ( ID ) e = mesh->AddFaceWithID(node[0], node[1], node[2], node[3],
155 node[4], node[5], ID);
156 else e = mesh->AddFace (node[0], node[1], node[2], node[3],
158 else if (nbnode == 8)
159 if ( ID ) e = mesh->AddFaceWithID(node[0], node[1], node[2], node[3],
160 node[4], node[5], node[6], node[7], ID);
161 else e = mesh->AddFace (node[0], node[1], node[2], node[3],
162 node[4], node[5], node[6], node[7] );
164 if ( ID ) e = mesh->AddPolygonalFaceWithID(node, ID);
165 else e = mesh->AddPolygonalFace (node );
171 if ( ID ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3], ID);
172 else e = mesh->AddVolume (node[0], node[1], node[2], node[3] );
173 else if (nbnode == 5)
174 if ( ID ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3],
176 else e = mesh->AddVolume (node[0], node[1], node[2], node[3],
178 else if (nbnode == 6)
179 if ( ID ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3],
180 node[4], node[5], ID);
181 else e = mesh->AddVolume (node[0], node[1], node[2], node[3],
183 else if (nbnode == 8)
184 if ( ID ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3],
185 node[4], node[5], node[6], node[7], ID);
186 else e = mesh->AddVolume (node[0], node[1], node[2], node[3],
187 node[4], node[5], node[6], node[7] );
188 else if (nbnode == 10)
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], ID);
192 else e = mesh->AddVolume (node[0], node[1], node[2], node[3],
193 node[4], node[5], node[6], node[7],
195 else if (nbnode == 13)
196 if ( ID ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3],
197 node[4], node[5], node[6], node[7],
198 node[8], node[9], node[10],node[11],
200 else e = mesh->AddVolume (node[0], node[1], node[2], node[3],
201 node[4], node[5], node[6], node[7],
202 node[8], node[9], node[10],node[11],
204 else if (nbnode == 15)
205 if ( ID ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3],
206 node[4], node[5], node[6], node[7],
207 node[8], node[9], node[10],node[11],
208 node[12],node[13],node[14],ID);
209 else e = mesh->AddVolume (node[0], node[1], node[2], node[3],
210 node[4], node[5], node[6], node[7],
211 node[8], node[9], node[10],node[11],
212 node[12],node[13],node[14] );
213 else if (nbnode == 20)
214 if ( ID ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3],
215 node[4], node[5], node[6], node[7],
216 node[8], node[9], node[10],node[11],
217 node[12],node[13],node[14],node[15],
218 node[16],node[17],node[18],node[19],ID);
219 else e = mesh->AddVolume (node[0], node[1], node[2], node[3],
220 node[4], node[5], node[6], node[7],
221 node[8], node[9], node[10],node[11],
222 node[12],node[13],node[14],node[15],
223 node[16],node[17],node[18],node[19] );
229 //=======================================================================
233 //=======================================================================
235 SMDS_MeshElement* SMESH_MeshEditor::AddElement(const vector<int> & nodeIDs,
236 const SMDSAbs_ElementType type,
240 vector<const SMDS_MeshNode*> nodes;
241 nodes.reserve( nodeIDs.size() );
242 vector<int>::const_iterator id = nodeIDs.begin();
243 while ( id != nodeIDs.end() ) {
244 if ( const SMDS_MeshNode* node = GetMeshDS()->FindNode( *id++ ))
245 nodes.push_back( node );
249 return AddElement( nodes, type, isPoly, ID );
252 //=======================================================================
254 //purpose : Remove a node or an element.
255 // Modify a compute state of sub-meshes which become empty
256 //=======================================================================
258 bool SMESH_MeshEditor::Remove (const list< int >& theIDs,
261 myLastCreatedElems.Clear();
262 myLastCreatedNodes.Clear();
264 SMESHDS_Mesh* aMesh = GetMeshDS();
265 set< SMESH_subMesh *> smmap;
267 list<int>::const_iterator it = theIDs.begin();
268 for ( ; it != theIDs.end(); it++ ) {
269 const SMDS_MeshElement * elem;
271 elem = aMesh->FindNode( *it );
273 elem = aMesh->FindElement( *it );
277 // Notify VERTEX sub-meshes about modification
279 const SMDS_MeshNode* node = cast2Node( elem );
280 if ( node->GetPosition()->GetTypeOfPosition() == SMDS_TOP_VERTEX )
281 if ( int aShapeID = node->GetPosition()->GetShapeId() )
282 if ( SMESH_subMesh * sm = GetMesh()->GetSubMeshContaining( aShapeID ) )
285 // Find sub-meshes to notify about modification
286 // SMDS_ElemIteratorPtr nodeIt = elem->nodesIterator();
287 // while ( nodeIt->more() ) {
288 // const SMDS_MeshNode* node = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
289 // const SMDS_PositionPtr& aPosition = node->GetPosition();
290 // if ( aPosition.get() ) {
291 // if ( int aShapeID = aPosition->GetShapeId() ) {
292 // if ( SMESH_subMesh * sm = GetMesh()->GetSubMeshContaining( aShapeID ) )
293 // smmap.insert( sm );
300 aMesh->RemoveNode( static_cast< const SMDS_MeshNode* >( elem ));
302 aMesh->RemoveElement( elem );
305 // Notify sub-meshes about modification
306 if ( !smmap.empty() ) {
307 set< SMESH_subMesh *>::iterator smIt;
308 for ( smIt = smmap.begin(); smIt != smmap.end(); smIt++ )
309 (*smIt)->ComputeStateEngine( SMESH_subMesh::MESH_ENTITY_REMOVED );
312 // // Check if the whole mesh becomes empty
313 // if ( SMESH_subMesh * sm = GetMesh()->GetSubMeshContaining( 1 ) )
314 // sm->ComputeStateEngine( SMESH_subMesh::CHECK_COMPUTE_STATE );
319 //=======================================================================
320 //function : FindShape
321 //purpose : Return an index of the shape theElem is on
322 // or zero if a shape not found
323 //=======================================================================
325 int SMESH_MeshEditor::FindShape (const SMDS_MeshElement * theElem)
327 myLastCreatedElems.Clear();
328 myLastCreatedNodes.Clear();
330 SMESHDS_Mesh * aMesh = GetMeshDS();
331 if ( aMesh->ShapeToMesh().IsNull() )
334 if ( theElem->GetType() == SMDSAbs_Node ) {
335 const SMDS_PositionPtr& aPosition =
336 static_cast<const SMDS_MeshNode*>( theElem )->GetPosition();
337 if ( aPosition.get() )
338 return aPosition->GetShapeId();
343 TopoDS_Shape aShape; // the shape a node is on
344 SMDS_ElemIteratorPtr nodeIt = theElem->nodesIterator();
345 while ( nodeIt->more() ) {
346 const SMDS_MeshNode* node = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
347 const SMDS_PositionPtr& aPosition = node->GetPosition();
348 if ( aPosition.get() ) {
349 int aShapeID = aPosition->GetShapeId();
350 SMESHDS_SubMesh * sm = aMesh->MeshElements( aShapeID );
352 if ( sm->Contains( theElem ))
354 if ( aShape.IsNull() )
355 aShape = aMesh->IndexToShape( aShapeID );
358 //MESSAGE ( "::FindShape() No SubShape for aShapeID " << aShapeID );
363 // None of nodes is on a proper shape,
364 // find the shape among ancestors of aShape on which a node is
365 if ( aShape.IsNull() ) {
366 //MESSAGE ("::FindShape() - NONE node is on shape")
369 TopTools_ListIteratorOfListOfShape ancIt( GetMesh()->GetAncestors( aShape ));
370 for ( ; ancIt.More(); ancIt.Next() ) {
371 SMESHDS_SubMesh * sm = aMesh->MeshElements( ancIt.Value() );
372 if ( sm && sm->Contains( theElem ))
373 return aMesh->ShapeToIndex( ancIt.Value() );
376 //MESSAGE ("::FindShape() - SHAPE NOT FOUND")
380 //=======================================================================
381 //function : IsMedium
383 //=======================================================================
385 bool SMESH_MeshEditor::IsMedium(const SMDS_MeshNode* node,
386 const SMDSAbs_ElementType typeToCheck)
388 bool isMedium = false;
389 SMDS_ElemIteratorPtr it = node->GetInverseElementIterator(typeToCheck);
390 while (it->more() && !isMedium ) {
391 const SMDS_MeshElement* elem = it->next();
392 isMedium = elem->IsMediumNode(node);
397 //=======================================================================
398 //function : ShiftNodesQuadTria
400 // Shift nodes in the array corresponded to quadratic triangle
401 // example: (0,1,2,3,4,5) -> (1,2,0,4,5,3)
402 //=======================================================================
403 static void ShiftNodesQuadTria(const SMDS_MeshNode* aNodes[])
405 const SMDS_MeshNode* nd1 = aNodes[0];
406 aNodes[0] = aNodes[1];
407 aNodes[1] = aNodes[2];
409 const SMDS_MeshNode* nd2 = aNodes[3];
410 aNodes[3] = aNodes[4];
411 aNodes[4] = aNodes[5];
415 //=======================================================================
416 //function : GetNodesFromTwoTria
418 // Shift nodes in the array corresponded to quadratic triangle
419 // example: (0,1,2,3,4,5) -> (1,2,0,4,5,3)
420 //=======================================================================
421 static bool GetNodesFromTwoTria(const SMDS_MeshElement * theTria1,
422 const SMDS_MeshElement * theTria2,
423 const SMDS_MeshNode* N1[],
424 const SMDS_MeshNode* N2[])
426 SMDS_ElemIteratorPtr it = theTria1->nodesIterator();
429 N1[i] = static_cast<const SMDS_MeshNode*>( it->next() );
432 if(it->more()) return false;
433 it = theTria2->nodesIterator();
436 N2[i] = static_cast<const SMDS_MeshNode*>( it->next() );
439 if(it->more()) return false;
441 int sames[3] = {-1,-1,-1};
453 if(nbsames!=2) return false;
455 ShiftNodesQuadTria(N1);
457 ShiftNodesQuadTria(N1);
460 i = sames[0] + sames[1] + sames[2];
462 ShiftNodesQuadTria(N2);
464 // now we receive following N1 and N2 (using numeration as above image)
465 // tria1 : (1 2 4 5 9 7) and tria2 : (3 4 2 8 9 6)
466 // i.e. first nodes from both arrays determ new diagonal
470 //=======================================================================
471 //function : InverseDiag
472 //purpose : Replace two neighbour triangles with ones built on the same 4 nodes
473 // but having other common link.
474 // Return False if args are improper
475 //=======================================================================
477 bool SMESH_MeshEditor::InverseDiag (const SMDS_MeshElement * theTria1,
478 const SMDS_MeshElement * theTria2 )
480 myLastCreatedElems.Clear();
481 myLastCreatedNodes.Clear();
483 if (!theTria1 || !theTria2)
486 const SMDS_FaceOfNodes* F1 = dynamic_cast<const SMDS_FaceOfNodes*>( theTria1 );
487 const SMDS_FaceOfNodes* F2 = dynamic_cast<const SMDS_FaceOfNodes*>( theTria2 );
490 // 1 +--+ A theTria1: ( 1 A B ) A->2 ( 1 2 B ) 1 +--+ A
491 // | /| theTria2: ( B A 2 ) B->1 ( 1 A 2 ) |\ |
495 // put nodes in array and find out indices of the same ones
496 const SMDS_MeshNode* aNodes [6];
497 int sameInd [] = { 0, 0, 0, 0, 0, 0 };
499 SMDS_ElemIteratorPtr it = theTria1->nodesIterator();
500 while ( it->more() ) {
501 aNodes[ i ] = static_cast<const SMDS_MeshNode*>( it->next() );
503 if ( i > 2 ) // theTria2
504 // find same node of theTria1
505 for ( int j = 0; j < 3; j++ )
506 if ( aNodes[ i ] == aNodes[ j ]) {
515 return false; // theTria1 is not a triangle
516 it = theTria2->nodesIterator();
518 if ( i == 6 && it->more() )
519 return false; // theTria2 is not a triangle
522 // find indices of 1,2 and of A,B in theTria1
523 int iA = 0, iB = 0, i1 = 0, i2 = 0;
524 for ( i = 0; i < 6; i++ ) {
525 if ( sameInd [ i ] == 0 )
532 // nodes 1 and 2 should not be the same
533 if ( aNodes[ i1 ] == aNodes[ i2 ] )
537 aNodes[ iA ] = aNodes[ i2 ];
539 aNodes[ sameInd[ iB ]] = aNodes[ i1 ];
541 //MESSAGE( theTria1 << theTria2 );
543 GetMeshDS()->ChangeElementNodes( theTria1, aNodes, 3 );
544 GetMeshDS()->ChangeElementNodes( theTria2, &aNodes[ 3 ], 3 );
546 //MESSAGE( theTria1 << theTria2 );
550 } // end if(F1 && F2)
552 // check case of quadratic faces
553 const SMDS_QuadraticFaceOfNodes* QF1 =
554 dynamic_cast<const SMDS_QuadraticFaceOfNodes*> (theTria1);
555 if(!QF1) return false;
556 const SMDS_QuadraticFaceOfNodes* QF2 =
557 dynamic_cast<const SMDS_QuadraticFaceOfNodes*> (theTria2);
558 if(!QF2) return false;
561 // 1 +--+--+ 2 theTria1: (1 2 4 5 9 7) or (2 4 1 9 7 5) or (4 1 2 7 5 9)
562 // | /| theTria2: (2 3 4 6 8 9) or (3 4 2 8 9 6) or (4 2 3 9 6 8)
570 const SMDS_MeshNode* N1 [6];
571 const SMDS_MeshNode* N2 [6];
572 if(!GetNodesFromTwoTria(theTria1,theTria2,N1,N2))
574 // now we receive following N1 and N2 (using numeration as above image)
575 // tria1 : (1 2 4 5 9 7) and tria2 : (3 4 2 8 9 6)
576 // i.e. first nodes from both arrays determ new diagonal
578 const SMDS_MeshNode* N1new [6];
579 const SMDS_MeshNode* N2new [6];
592 // replaces nodes in faces
593 GetMeshDS()->ChangeElementNodes( theTria1, N1new, 6 );
594 GetMeshDS()->ChangeElementNodes( theTria2, N2new, 6 );
599 //=======================================================================
600 //function : findTriangles
601 //purpose : find triangles sharing theNode1-theNode2 link
602 //=======================================================================
604 static bool findTriangles(const SMDS_MeshNode * theNode1,
605 const SMDS_MeshNode * theNode2,
606 const SMDS_MeshElement*& theTria1,
607 const SMDS_MeshElement*& theTria2)
609 if ( !theNode1 || !theNode2 ) return false;
611 theTria1 = theTria2 = 0;
613 set< const SMDS_MeshElement* > emap;
614 SMDS_ElemIteratorPtr it = theNode1->GetInverseElementIterator(SMDSAbs_Face);
616 const SMDS_MeshElement* elem = it->next();
617 if ( elem->NbNodes() == 3 )
620 it = theNode2->GetInverseElementIterator(SMDSAbs_Face);
622 const SMDS_MeshElement* elem = it->next();
623 if ( emap.find( elem ) != emap.end() )
625 // theTria1 must be element with minimum ID
626 if( theTria1->GetID() < elem->GetID() ) {
639 return ( theTria1 && theTria2 );
642 //=======================================================================
643 //function : InverseDiag
644 //purpose : Replace two neighbour triangles sharing theNode1-theNode2 link
645 // with ones built on the same 4 nodes but having other common link.
646 // Return false if proper faces not found
647 //=======================================================================
649 bool SMESH_MeshEditor::InverseDiag (const SMDS_MeshNode * theNode1,
650 const SMDS_MeshNode * theNode2)
652 myLastCreatedElems.Clear();
653 myLastCreatedNodes.Clear();
655 MESSAGE( "::InverseDiag()" );
657 const SMDS_MeshElement *tr1, *tr2;
658 if ( !findTriangles( theNode1, theNode2, tr1, tr2 ))
661 const SMDS_FaceOfNodes* F1 = dynamic_cast<const SMDS_FaceOfNodes*>( tr1 );
662 //if (!F1) return false;
663 const SMDS_FaceOfNodes* F2 = dynamic_cast<const SMDS_FaceOfNodes*>( tr2 );
664 //if (!F2) return false;
667 // 1 +--+ A tr1: ( 1 A B ) A->2 ( 1 2 B ) 1 +--+ A
668 // | /| tr2: ( B A 2 ) B->1 ( 1 A 2 ) |\ |
672 // put nodes in array
673 // and find indices of 1,2 and of A in tr1 and of B in tr2
674 int i, iA1 = 0, i1 = 0;
675 const SMDS_MeshNode* aNodes1 [3];
676 SMDS_ElemIteratorPtr it;
677 for (i = 0, it = tr1->nodesIterator(); it->more(); i++ ) {
678 aNodes1[ i ] = static_cast<const SMDS_MeshNode*>( it->next() );
679 if ( aNodes1[ i ] == theNode1 )
680 iA1 = i; // node A in tr1
681 else if ( aNodes1[ i ] != theNode2 )
685 const SMDS_MeshNode* aNodes2 [3];
686 for (i = 0, it = tr2->nodesIterator(); it->more(); i++ ) {
687 aNodes2[ i ] = static_cast<const SMDS_MeshNode*>( it->next() );
688 if ( aNodes2[ i ] == theNode2 )
689 iB2 = i; // node B in tr2
690 else if ( aNodes2[ i ] != theNode1 )
694 // nodes 1 and 2 should not be the same
695 if ( aNodes1[ i1 ] == aNodes2[ i2 ] )
699 aNodes1[ iA1 ] = aNodes2[ i2 ];
701 aNodes2[ iB2 ] = aNodes1[ i1 ];
703 //MESSAGE( tr1 << tr2 );
705 GetMeshDS()->ChangeElementNodes( tr1, aNodes1, 3 );
706 GetMeshDS()->ChangeElementNodes( tr2, aNodes2, 3 );
708 //MESSAGE( tr1 << tr2 );
713 // check case of quadratic faces
714 const SMDS_QuadraticFaceOfNodes* QF1 =
715 dynamic_cast<const SMDS_QuadraticFaceOfNodes*> (tr1);
716 if(!QF1) return false;
717 const SMDS_QuadraticFaceOfNodes* QF2 =
718 dynamic_cast<const SMDS_QuadraticFaceOfNodes*> (tr2);
719 if(!QF2) return false;
720 return InverseDiag(tr1,tr2);
723 //=======================================================================
724 //function : getQuadrangleNodes
725 //purpose : fill theQuadNodes - nodes of a quadrangle resulting from
726 // fusion of triangles tr1 and tr2 having shared link on
727 // theNode1 and theNode2
728 //=======================================================================
730 bool getQuadrangleNodes(const SMDS_MeshNode * theQuadNodes [],
731 const SMDS_MeshNode * theNode1,
732 const SMDS_MeshNode * theNode2,
733 const SMDS_MeshElement * tr1,
734 const SMDS_MeshElement * tr2 )
736 if( tr1->NbNodes() != tr2->NbNodes() )
738 // find the 4-th node to insert into tr1
739 const SMDS_MeshNode* n4 = 0;
740 SMDS_ElemIteratorPtr it = tr2->nodesIterator();
742 while ( !n4 && i<3 ) {
743 const SMDS_MeshNode * n = cast2Node( it->next() );
745 bool isDiag = ( n == theNode1 || n == theNode2 );
749 // Make an array of nodes to be in a quadrangle
750 int iNode = 0, iFirstDiag = -1;
751 it = tr1->nodesIterator();
754 const SMDS_MeshNode * n = cast2Node( it->next() );
756 bool isDiag = ( n == theNode1 || n == theNode2 );
758 if ( iFirstDiag < 0 )
760 else if ( iNode - iFirstDiag == 1 )
761 theQuadNodes[ iNode++ ] = n4; // insert the 4-th node between diagonal nodes
763 else if ( n == n4 ) {
764 return false; // tr1 and tr2 should not have all the same nodes
766 theQuadNodes[ iNode++ ] = n;
768 if ( iNode == 3 ) // diagonal nodes have 0 and 2 indices
769 theQuadNodes[ iNode ] = n4;
774 //=======================================================================
775 //function : DeleteDiag
776 //purpose : Replace two neighbour triangles sharing theNode1-theNode2 link
777 // with a quadrangle built on the same 4 nodes.
778 // Return false if proper faces not found
779 //=======================================================================
781 bool SMESH_MeshEditor::DeleteDiag (const SMDS_MeshNode * theNode1,
782 const SMDS_MeshNode * theNode2)
784 myLastCreatedElems.Clear();
785 myLastCreatedNodes.Clear();
787 MESSAGE( "::DeleteDiag()" );
789 const SMDS_MeshElement *tr1, *tr2;
790 if ( !findTriangles( theNode1, theNode2, tr1, tr2 ))
793 const SMDS_FaceOfNodes* F1 = dynamic_cast<const SMDS_FaceOfNodes*>( tr1 );
794 //if (!F1) return false;
795 const SMDS_FaceOfNodes* F2 = dynamic_cast<const SMDS_FaceOfNodes*>( tr2 );
796 //if (!F2) return false;
799 const SMDS_MeshNode* aNodes [ 4 ];
800 if ( ! getQuadrangleNodes( aNodes, theNode1, theNode2, tr1, tr2 ))
803 //MESSAGE( endl << tr1 << tr2 );
805 GetMeshDS()->ChangeElementNodes( tr1, aNodes, 4 );
806 myLastCreatedElems.Append(tr1);
807 GetMeshDS()->RemoveElement( tr2 );
809 //MESSAGE( endl << tr1 );
814 // check case of quadratic faces
815 const SMDS_QuadraticFaceOfNodes* QF1 =
816 dynamic_cast<const SMDS_QuadraticFaceOfNodes*> (tr1);
817 if(!QF1) return false;
818 const SMDS_QuadraticFaceOfNodes* QF2 =
819 dynamic_cast<const SMDS_QuadraticFaceOfNodes*> (tr2);
820 if(!QF2) return false;
823 // 1 +--+--+ 2 tr1: (1 2 4 5 9 7) or (2 4 1 9 7 5) or (4 1 2 7 5 9)
824 // | /| tr2: (2 3 4 6 8 9) or (3 4 2 8 9 6) or (4 2 3 9 6 8)
832 const SMDS_MeshNode* N1 [6];
833 const SMDS_MeshNode* N2 [6];
834 if(!GetNodesFromTwoTria(tr1,tr2,N1,N2))
836 // now we receive following N1 and N2 (using numeration as above image)
837 // tria1 : (1 2 4 5 9 7) and tria2 : (3 4 2 8 9 6)
838 // i.e. first nodes from both arrays determ new diagonal
840 const SMDS_MeshNode* aNodes[8];
850 GetMeshDS()->ChangeElementNodes( tr1, aNodes, 8 );
851 myLastCreatedElems.Append(tr1);
852 GetMeshDS()->RemoveElement( tr2 );
854 // remove middle node (9)
855 GetMeshDS()->RemoveNode( N1[4] );
860 //=======================================================================
861 //function : Reorient
862 //purpose : Reverse theElement orientation
863 //=======================================================================
865 bool SMESH_MeshEditor::Reorient (const SMDS_MeshElement * theElem)
867 myLastCreatedElems.Clear();
868 myLastCreatedNodes.Clear();
872 SMDS_ElemIteratorPtr it = theElem->nodesIterator();
873 if ( !it || !it->more() )
876 switch ( theElem->GetType() ) {
880 if(!theElem->IsQuadratic()) {
881 int i = theElem->NbNodes();
882 vector<const SMDS_MeshNode*> aNodes( i );
884 aNodes[ --i ]= static_cast<const SMDS_MeshNode*>( it->next() );
885 return GetMeshDS()->ChangeElementNodes( theElem, &aNodes[0], theElem->NbNodes() );
888 // quadratic elements
889 if(theElem->GetType()==SMDSAbs_Edge) {
890 vector<const SMDS_MeshNode*> aNodes(3);
891 aNodes[1]= static_cast<const SMDS_MeshNode*>( it->next() );
892 aNodes[0]= static_cast<const SMDS_MeshNode*>( it->next() );
893 aNodes[2]= static_cast<const SMDS_MeshNode*>( it->next() );
894 return GetMeshDS()->ChangeElementNodes( theElem, &aNodes[0], 3 );
897 int nbn = theElem->NbNodes();
898 vector<const SMDS_MeshNode*> aNodes(nbn);
899 aNodes[0]= static_cast<const SMDS_MeshNode*>( it->next() );
901 for(; i<nbn/2; i++) {
902 aNodes[nbn/2-i]= static_cast<const SMDS_MeshNode*>( it->next() );
904 for(i=0; i<nbn/2; i++) {
905 aNodes[nbn-i-1]= static_cast<const SMDS_MeshNode*>( it->next() );
907 return GetMeshDS()->ChangeElementNodes( theElem, &aNodes[0], nbn );
911 case SMDSAbs_Volume: {
912 if (theElem->IsPoly()) {
913 const SMDS_PolyhedralVolumeOfNodes* aPolyedre =
914 static_cast<const SMDS_PolyhedralVolumeOfNodes*>( theElem );
916 MESSAGE("Warning: bad volumic element");
920 int nbFaces = aPolyedre->NbFaces();
921 vector<const SMDS_MeshNode *> poly_nodes;
922 vector<int> quantities (nbFaces);
924 // reverse each face of the polyedre
925 for (int iface = 1; iface <= nbFaces; iface++) {
926 int inode, nbFaceNodes = aPolyedre->NbFaceNodes(iface);
927 quantities[iface - 1] = nbFaceNodes;
929 for (inode = nbFaceNodes; inode >= 1; inode--) {
930 const SMDS_MeshNode* curNode = aPolyedre->GetFaceNode(iface, inode);
931 poly_nodes.push_back(curNode);
935 return GetMeshDS()->ChangePolyhedronNodes( theElem, poly_nodes, quantities );
939 SMDS_VolumeTool vTool;
940 if ( !vTool.Set( theElem ))
943 return GetMeshDS()->ChangeElementNodes( theElem, vTool.GetNodes(), vTool.NbNodes() );
952 //=======================================================================
953 //function : getBadRate
955 //=======================================================================
957 static double getBadRate (const SMDS_MeshElement* theElem,
958 SMESH::Controls::NumericalFunctorPtr& theCrit)
960 SMESH::Controls::TSequenceOfXYZ P;
961 if ( !theElem || !theCrit->GetPoints( theElem, P ))
963 return theCrit->GetBadRate( theCrit->GetValue( P ), theElem->NbNodes() );
964 //return theCrit->GetBadRate( theCrit->GetValue( theElem->GetID() ), theElem->NbNodes() );
967 //=======================================================================
968 //function : QuadToTri
969 //purpose : Cut quadrangles into triangles.
970 // theCrit is used to select a diagonal to cut
971 //=======================================================================
973 bool SMESH_MeshEditor::QuadToTri (TIDSortedElemSet & theElems,
974 SMESH::Controls::NumericalFunctorPtr theCrit)
976 myLastCreatedElems.Clear();
977 myLastCreatedNodes.Clear();
979 MESSAGE( "::QuadToTri()" );
981 if ( !theCrit.get() )
984 SMESHDS_Mesh * aMesh = GetMeshDS();
986 Handle(Geom_Surface) surface;
987 SMESH_MesherHelper helper( *GetMesh() );
989 TIDSortedElemSet::iterator itElem;
990 for ( itElem = theElems.begin(); itElem != theElems.end(); itElem++ ) {
991 const SMDS_MeshElement* elem = *itElem;
992 if ( !elem || elem->GetType() != SMDSAbs_Face )
994 if ( elem->NbNodes() != ( elem->IsQuadratic() ? 8 : 4 ))
997 // retrieve element nodes
998 const SMDS_MeshNode* aNodes [8];
999 SMDS_ElemIteratorPtr itN = elem->nodesIterator();
1001 while ( itN->more() )
1002 aNodes[ i++ ] = static_cast<const SMDS_MeshNode*>( itN->next() );
1004 // compare two sets of possible triangles
1005 double aBadRate1, aBadRate2; // to what extent a set is bad
1006 SMDS_FaceOfNodes tr1 ( aNodes[0], aNodes[1], aNodes[2] );
1007 SMDS_FaceOfNodes tr2 ( aNodes[2], aNodes[3], aNodes[0] );
1008 aBadRate1 = getBadRate( &tr1, theCrit ) + getBadRate( &tr2, theCrit );
1010 SMDS_FaceOfNodes tr3 ( aNodes[1], aNodes[2], aNodes[3] );
1011 SMDS_FaceOfNodes tr4 ( aNodes[3], aNodes[0], aNodes[1] );
1012 aBadRate2 = getBadRate( &tr3, theCrit ) + getBadRate( &tr4, theCrit );
1014 int aShapeId = FindShape( elem );
1015 const SMDS_MeshElement* newElem = 0;
1017 if( !elem->IsQuadratic() ) {
1019 // split liner quadrangle
1021 if ( aBadRate1 <= aBadRate2 ) {
1022 // tr1 + tr2 is better
1023 aMesh->ChangeElementNodes( elem, aNodes, 3 );
1024 newElem = aMesh->AddFace( aNodes[2], aNodes[3], aNodes[0] );
1027 // tr3 + tr4 is better
1028 aMesh->ChangeElementNodes( elem, &aNodes[1], 3 );
1029 newElem = aMesh->AddFace( aNodes[3], aNodes[0], aNodes[1] );
1034 // split quadratic quadrangle
1036 // get surface elem is on
1037 if ( aShapeId != helper.GetSubShapeID() ) {
1041 shape = aMesh->IndexToShape( aShapeId );
1042 if ( !shape.IsNull() && shape.ShapeType() == TopAbs_FACE ) {
1043 TopoDS_Face face = TopoDS::Face( shape );
1044 surface = BRep_Tool::Surface( face );
1045 if ( !surface.IsNull() )
1046 helper.SetSubShape( shape );
1050 const SMDS_MeshNode* aNodes [8];
1051 const SMDS_MeshNode* inFaceNode = 0;
1052 SMDS_ElemIteratorPtr itN = elem->nodesIterator();
1054 while ( itN->more() ) {
1055 aNodes[ i++ ] = static_cast<const SMDS_MeshNode*>( itN->next() );
1056 if ( !inFaceNode && helper.GetNodeUVneedInFaceNode() &&
1057 aNodes[ i-1 ]->GetPosition()->GetTypeOfPosition() == SMDS_TOP_FACE )
1059 inFaceNode = aNodes[ i-1 ];
1062 // find middle point for (0,1,2,3)
1063 // and create a node in this point;
1065 if ( surface.IsNull() ) {
1067 p += gp_XYZ(aNodes[i]->X(), aNodes[i]->Y(), aNodes[i]->Z() );
1071 TopoDS_Face face = TopoDS::Face( helper.GetSubShape() );
1074 uv += helper.GetNodeUV( face, aNodes[i], inFaceNode );
1076 p = surface->Value( uv.X(), uv.Y() ).XYZ();
1078 const SMDS_MeshNode* newN = aMesh->AddNode( p.X(), p.Y(), p.Z() );
1079 myLastCreatedNodes.Append(newN);
1081 // create a new element
1082 const SMDS_MeshNode* N[6];
1083 if ( aBadRate1 <= aBadRate2 ) {
1090 newElem = aMesh->AddFace(aNodes[2], aNodes[3], aNodes[0],
1091 aNodes[6], aNodes[7], newN );
1100 newElem = aMesh->AddFace(aNodes[3], aNodes[0], aNodes[1],
1101 aNodes[7], aNodes[4], newN );
1103 aMesh->ChangeElementNodes( elem, N, 6 );
1107 // care of a new element
1109 myLastCreatedElems.Append(newElem);
1110 AddToSameGroups( newElem, elem, aMesh );
1112 // put a new triangle on the same shape
1114 aMesh->SetMeshElementOnShape( newElem, aShapeId );
1119 //=======================================================================
1120 //function : BestSplit
1121 //purpose : Find better diagonal for cutting.
1122 //=======================================================================
1123 int SMESH_MeshEditor::BestSplit (const SMDS_MeshElement* theQuad,
1124 SMESH::Controls::NumericalFunctorPtr theCrit)
1126 myLastCreatedElems.Clear();
1127 myLastCreatedNodes.Clear();
1132 if (!theQuad || theQuad->GetType() != SMDSAbs_Face )
1135 if( theQuad->NbNodes()==4 ||
1136 (theQuad->NbNodes()==8 && theQuad->IsQuadratic()) ) {
1138 // retrieve element nodes
1139 const SMDS_MeshNode* aNodes [4];
1140 SMDS_ElemIteratorPtr itN = theQuad->nodesIterator();
1142 //while (itN->more())
1144 aNodes[ i++ ] = static_cast<const SMDS_MeshNode*>( itN->next() );
1146 // compare two sets of possible triangles
1147 double aBadRate1, aBadRate2; // to what extent a set is bad
1148 SMDS_FaceOfNodes tr1 ( aNodes[0], aNodes[1], aNodes[2] );
1149 SMDS_FaceOfNodes tr2 ( aNodes[2], aNodes[3], aNodes[0] );
1150 aBadRate1 = getBadRate( &tr1, theCrit ) + getBadRate( &tr2, theCrit );
1152 SMDS_FaceOfNodes tr3 ( aNodes[1], aNodes[2], aNodes[3] );
1153 SMDS_FaceOfNodes tr4 ( aNodes[3], aNodes[0], aNodes[1] );
1154 aBadRate2 = getBadRate( &tr3, theCrit ) + getBadRate( &tr4, theCrit );
1156 if (aBadRate1 <= aBadRate2) // tr1 + tr2 is better
1157 return 1; // diagonal 1-3
1159 return 2; // diagonal 2-4
1164 //=======================================================================
1165 //function : AddToSameGroups
1166 //purpose : add elemToAdd to the groups the elemInGroups belongs to
1167 //=======================================================================
1169 void SMESH_MeshEditor::AddToSameGroups (const SMDS_MeshElement* elemToAdd,
1170 const SMDS_MeshElement* elemInGroups,
1171 SMESHDS_Mesh * aMesh)
1173 const set<SMESHDS_GroupBase*>& groups = aMesh->GetGroups();
1174 set<SMESHDS_GroupBase*>::const_iterator grIt = groups.begin();
1175 for ( ; grIt != groups.end(); grIt++ ) {
1176 SMESHDS_Group* group = dynamic_cast<SMESHDS_Group*>( *grIt );
1177 if ( group && group->Contains( elemInGroups ))
1178 group->SMDSGroup().Add( elemToAdd );
1183 //=======================================================================
1184 //function : RemoveElemFromGroups
1185 //purpose : Remove removeelem to the groups the elemInGroups belongs to
1186 //=======================================================================
1187 void SMESH_MeshEditor::RemoveElemFromGroups (const SMDS_MeshElement* removeelem,
1188 SMESHDS_Mesh * aMesh)
1190 const set<SMESHDS_GroupBase*>& groups = aMesh->GetGroups();
1191 if (!groups.empty())
1193 set<SMESHDS_GroupBase*>::const_iterator GrIt = groups.begin();
1194 for (; GrIt != groups.end(); GrIt++)
1196 SMESHDS_Group* grp = dynamic_cast<SMESHDS_Group*>(*GrIt);
1197 if (!grp || grp->IsEmpty()) continue;
1198 grp->SMDSGroup().Remove(removeelem);
1204 //=======================================================================
1205 //function : QuadToTri
1206 //purpose : Cut quadrangles into triangles.
1207 // theCrit is used to select a diagonal to cut
1208 //=======================================================================
1210 bool SMESH_MeshEditor::QuadToTri (TIDSortedElemSet & theElems,
1211 const bool the13Diag)
1213 myLastCreatedElems.Clear();
1214 myLastCreatedNodes.Clear();
1216 MESSAGE( "::QuadToTri()" );
1218 SMESHDS_Mesh * aMesh = GetMeshDS();
1220 Handle(Geom_Surface) surface;
1221 SMESH_MesherHelper helper( *GetMesh() );
1223 TIDSortedElemSet::iterator itElem;
1224 for ( itElem = theElems.begin(); itElem != theElems.end(); itElem++ ) {
1225 const SMDS_MeshElement* elem = *itElem;
1226 if ( !elem || elem->GetType() != SMDSAbs_Face )
1228 bool isquad = elem->NbNodes()==4 || elem->NbNodes()==8;
1229 if(!isquad) continue;
1231 if(elem->NbNodes()==4) {
1232 // retrieve element nodes
1233 const SMDS_MeshNode* aNodes [4];
1234 SMDS_ElemIteratorPtr itN = elem->nodesIterator();
1236 while ( itN->more() )
1237 aNodes[ i++ ] = static_cast<const SMDS_MeshNode*>( itN->next() );
1239 int aShapeId = FindShape( elem );
1240 const SMDS_MeshElement* newElem = 0;
1242 aMesh->ChangeElementNodes( elem, aNodes, 3 );
1243 newElem = aMesh->AddFace( aNodes[2], aNodes[3], aNodes[0] );
1246 aMesh->ChangeElementNodes( elem, &aNodes[1], 3 );
1247 newElem = aMesh->AddFace( aNodes[3], aNodes[0], aNodes[1] );
1249 myLastCreatedElems.Append(newElem);
1250 // put a new triangle on the same shape and add to the same groups
1252 aMesh->SetMeshElementOnShape( newElem, aShapeId );
1253 AddToSameGroups( newElem, elem, aMesh );
1256 // Quadratic quadrangle
1258 if( elem->NbNodes()==8 && elem->IsQuadratic() ) {
1260 // get surface elem is on
1261 int aShapeId = FindShape( elem );
1262 if ( aShapeId != helper.GetSubShapeID() ) {
1266 shape = aMesh->IndexToShape( aShapeId );
1267 if ( !shape.IsNull() && shape.ShapeType() == TopAbs_FACE ) {
1268 TopoDS_Face face = TopoDS::Face( shape );
1269 surface = BRep_Tool::Surface( face );
1270 if ( !surface.IsNull() )
1271 helper.SetSubShape( shape );
1275 const SMDS_MeshNode* aNodes [8];
1276 const SMDS_MeshNode* inFaceNode = 0;
1277 SMDS_ElemIteratorPtr itN = elem->nodesIterator();
1279 while ( itN->more() ) {
1280 aNodes[ i++ ] = static_cast<const SMDS_MeshNode*>( itN->next() );
1281 if ( !inFaceNode && helper.GetNodeUVneedInFaceNode() &&
1282 aNodes[ i-1 ]->GetPosition()->GetTypeOfPosition() == SMDS_TOP_FACE )
1284 inFaceNode = aNodes[ i-1 ];
1288 // find middle point for (0,1,2,3)
1289 // and create a node in this point;
1291 if ( surface.IsNull() ) {
1293 p += gp_XYZ(aNodes[i]->X(), aNodes[i]->Y(), aNodes[i]->Z() );
1297 TopoDS_Face geomFace = TopoDS::Face( helper.GetSubShape() );
1300 uv += helper.GetNodeUV( geomFace, aNodes[i], inFaceNode );
1302 p = surface->Value( uv.X(), uv.Y() ).XYZ();
1304 const SMDS_MeshNode* newN = aMesh->AddNode( p.X(), p.Y(), p.Z() );
1305 myLastCreatedNodes.Append(newN);
1307 // create a new element
1308 const SMDS_MeshElement* newElem = 0;
1309 const SMDS_MeshNode* N[6];
1317 newElem = aMesh->AddFace(aNodes[2], aNodes[3], aNodes[0],
1318 aNodes[6], aNodes[7], newN );
1327 newElem = aMesh->AddFace(aNodes[3], aNodes[0], aNodes[1],
1328 aNodes[7], aNodes[4], newN );
1330 myLastCreatedElems.Append(newElem);
1331 aMesh->ChangeElementNodes( elem, N, 6 );
1332 // put a new triangle on the same shape and add to the same groups
1334 aMesh->SetMeshElementOnShape( newElem, aShapeId );
1335 AddToSameGroups( newElem, elem, aMesh );
1342 //=======================================================================
1343 //function : getAngle
1345 //=======================================================================
1347 double getAngle(const SMDS_MeshElement * tr1,
1348 const SMDS_MeshElement * tr2,
1349 const SMDS_MeshNode * n1,
1350 const SMDS_MeshNode * n2)
1352 double angle = 2*PI; // bad angle
1355 SMESH::Controls::TSequenceOfXYZ P1, P2;
1356 if ( !SMESH::Controls::NumericalFunctor::GetPoints( tr1, P1 ) ||
1357 !SMESH::Controls::NumericalFunctor::GetPoints( tr2, P2 ))
1360 if(!tr1->IsQuadratic())
1361 N1 = gp_Vec( P1(2) - P1(1) ) ^ gp_Vec( P1(3) - P1(1) );
1363 N1 = gp_Vec( P1(3) - P1(1) ) ^ gp_Vec( P1(5) - P1(1) );
1364 if ( N1.SquareMagnitude() <= gp::Resolution() )
1366 if(!tr2->IsQuadratic())
1367 N2 = gp_Vec( P2(2) - P2(1) ) ^ gp_Vec( P2(3) - P2(1) );
1369 N2 = gp_Vec( P2(3) - P2(1) ) ^ gp_Vec( P2(5) - P2(1) );
1370 if ( N2.SquareMagnitude() <= gp::Resolution() )
1373 // find the first diagonal node n1 in the triangles:
1374 // take in account a diagonal link orientation
1375 const SMDS_MeshElement *nFirst[2], *tr[] = { tr1, tr2 };
1376 for ( int t = 0; t < 2; t++ ) {
1377 SMDS_ElemIteratorPtr it = tr[ t ]->nodesIterator();
1378 int i = 0, iDiag = -1;
1379 while ( it->more()) {
1380 const SMDS_MeshElement *n = it->next();
1381 if ( n == n1 || n == n2 )
1385 if ( i - iDiag == 1 )
1386 nFirst[ t ] = ( n == n1 ? n2 : n1 );
1394 if ( nFirst[ 0 ] == nFirst[ 1 ] )
1397 angle = N1.Angle( N2 );
1402 // =================================================
1403 // class generating a unique ID for a pair of nodes
1404 // and able to return nodes by that ID
1405 // =================================================
1409 LinkID_Gen( const SMESHDS_Mesh* theMesh )
1410 :myMesh( theMesh ), myMaxID( theMesh->MaxNodeID() + 1)
1413 long GetLinkID (const SMDS_MeshNode * n1,
1414 const SMDS_MeshNode * n2) const
1416 return ( Min(n1->GetID(),n2->GetID()) * myMaxID + Max(n1->GetID(),n2->GetID()));
1419 bool GetNodes (const long theLinkID,
1420 const SMDS_MeshNode* & theNode1,
1421 const SMDS_MeshNode* & theNode2) const
1423 theNode1 = myMesh->FindNode( theLinkID / myMaxID );
1424 if ( !theNode1 ) return false;
1425 theNode2 = myMesh->FindNode( theLinkID % myMaxID );
1426 if ( !theNode2 ) return false;
1432 const SMESHDS_Mesh* myMesh;
1437 //=======================================================================
1438 //function : TriToQuad
1439 //purpose : Fuse neighbour triangles into quadrangles.
1440 // theCrit is used to select a neighbour to fuse with.
1441 // theMaxAngle is a max angle between element normals at which
1442 // fusion is still performed.
1443 //=======================================================================
1445 bool SMESH_MeshEditor::TriToQuad (TIDSortedElemSet & theElems,
1446 SMESH::Controls::NumericalFunctorPtr theCrit,
1447 const double theMaxAngle)
1449 myLastCreatedElems.Clear();
1450 myLastCreatedNodes.Clear();
1452 MESSAGE( "::TriToQuad()" );
1454 if ( !theCrit.get() )
1457 SMESHDS_Mesh * aMesh = GetMeshDS();
1459 // Prepare data for algo: build
1460 // 1. map of elements with their linkIDs
1461 // 2. map of linkIDs with their elements
1463 map< TLink, list< const SMDS_MeshElement* > > mapLi_listEl;
1464 map< TLink, list< const SMDS_MeshElement* > >::iterator itLE;
1465 map< const SMDS_MeshElement*, set< TLink > > mapEl_setLi;
1466 map< const SMDS_MeshElement*, set< TLink > >::iterator itEL;
1468 TIDSortedElemSet::iterator itElem;
1469 for ( itElem = theElems.begin(); itElem != theElems.end(); itElem++ ) {
1470 const SMDS_MeshElement* elem = *itElem;
1471 if(!elem || elem->GetType() != SMDSAbs_Face ) continue;
1472 bool IsTria = elem->NbNodes()==3 || (elem->NbNodes()==6 && elem->IsQuadratic());
1473 if(!IsTria) continue;
1475 // retrieve element nodes
1476 const SMDS_MeshNode* aNodes [4];
1477 SMDS_ElemIteratorPtr itN = elem->nodesIterator();
1480 aNodes[ i++ ] = cast2Node( itN->next() );
1481 aNodes[ 3 ] = aNodes[ 0 ];
1484 for ( i = 0; i < 3; i++ ) {
1485 TLink link( aNodes[i], aNodes[i+1] );
1486 // check if elements sharing a link can be fused
1487 itLE = mapLi_listEl.find( link );
1488 if ( itLE != mapLi_listEl.end() ) {
1489 if ((*itLE).second.size() > 1 ) // consider only 2 elems adjacent by a link
1491 const SMDS_MeshElement* elem2 = (*itLE).second.front();
1492 //if ( FindShape( elem ) != FindShape( elem2 ))
1493 // continue; // do not fuse triangles laying on different shapes
1494 if ( getAngle( elem, elem2, aNodes[i], aNodes[i+1] ) > theMaxAngle )
1495 continue; // avoid making badly shaped quads
1496 (*itLE).second.push_back( elem );
1499 mapLi_listEl[ link ].push_back( elem );
1501 mapEl_setLi [ elem ].insert( link );
1504 // Clean the maps from the links shared by a sole element, ie
1505 // links to which only one element is bound in mapLi_listEl
1507 for ( itLE = mapLi_listEl.begin(); itLE != mapLi_listEl.end(); itLE++ ) {
1508 int nbElems = (*itLE).second.size();
1509 if ( nbElems < 2 ) {
1510 const SMDS_MeshElement* elem = (*itLE).second.front();
1511 TLink link = (*itLE).first;
1512 mapEl_setLi[ elem ].erase( link );
1513 if ( mapEl_setLi[ elem ].empty() )
1514 mapEl_setLi.erase( elem );
1518 // Algo: fuse triangles into quadrangles
1520 while ( ! mapEl_setLi.empty() ) {
1521 // Look for the start element:
1522 // the element having the least nb of shared links
1523 const SMDS_MeshElement* startElem = 0;
1525 for ( itEL = mapEl_setLi.begin(); itEL != mapEl_setLi.end(); itEL++ ) {
1526 int nbLinks = (*itEL).second.size();
1527 if ( nbLinks < minNbLinks ) {
1528 startElem = (*itEL).first;
1529 minNbLinks = nbLinks;
1530 if ( minNbLinks == 1 )
1535 // search elements to fuse starting from startElem or links of elements
1536 // fused earlyer - startLinks
1537 list< TLink > startLinks;
1538 while ( startElem || !startLinks.empty() ) {
1539 while ( !startElem && !startLinks.empty() ) {
1540 // Get an element to start, by a link
1541 TLink linkId = startLinks.front();
1542 startLinks.pop_front();
1543 itLE = mapLi_listEl.find( linkId );
1544 if ( itLE != mapLi_listEl.end() ) {
1545 list< const SMDS_MeshElement* > & listElem = (*itLE).second;
1546 list< const SMDS_MeshElement* >::iterator itE = listElem.begin();
1547 for ( ; itE != listElem.end() ; itE++ )
1548 if ( mapEl_setLi.find( (*itE) ) != mapEl_setLi.end() )
1550 mapLi_listEl.erase( itLE );
1555 // Get candidates to be fused
1556 const SMDS_MeshElement *tr1 = startElem, *tr2 = 0, *tr3 = 0;
1557 const TLink *link12, *link13;
1559 ASSERT( mapEl_setLi.find( tr1 ) != mapEl_setLi.end() );
1560 set< TLink >& setLi = mapEl_setLi[ tr1 ];
1561 ASSERT( !setLi.empty() );
1562 set< TLink >::iterator itLi;
1563 for ( itLi = setLi.begin(); itLi != setLi.end(); itLi++ )
1565 const TLink & link = (*itLi);
1566 itLE = mapLi_listEl.find( link );
1567 if ( itLE == mapLi_listEl.end() )
1570 const SMDS_MeshElement* elem = (*itLE).second.front();
1572 elem = (*itLE).second.back();
1573 mapLi_listEl.erase( itLE );
1574 if ( mapEl_setLi.find( elem ) == mapEl_setLi.end())
1585 // add other links of elem to list of links to re-start from
1586 set< TLink >& links = mapEl_setLi[ elem ];
1587 set< TLink >::iterator it;
1588 for ( it = links.begin(); it != links.end(); it++ ) {
1589 const TLink& link2 = (*it);
1590 if ( link2 != link )
1591 startLinks.push_back( link2 );
1595 // Get nodes of possible quadrangles
1596 const SMDS_MeshNode *n12 [4], *n13 [4];
1597 bool Ok12 = false, Ok13 = false;
1598 const SMDS_MeshNode *linkNode1, *linkNode2;
1600 linkNode1 = link12->first;
1601 linkNode2 = link12->second;
1602 if ( tr2 && getQuadrangleNodes( n12, linkNode1, linkNode2, tr1, tr2 ))
1606 linkNode1 = link13->first;
1607 linkNode2 = link13->second;
1608 if ( tr3 && getQuadrangleNodes( n13, linkNode1, linkNode2, tr1, tr3 ))
1612 // Choose a pair to fuse
1613 if ( Ok12 && Ok13 ) {
1614 SMDS_FaceOfNodes quad12 ( n12[ 0 ], n12[ 1 ], n12[ 2 ], n12[ 3 ] );
1615 SMDS_FaceOfNodes quad13 ( n13[ 0 ], n13[ 1 ], n13[ 2 ], n13[ 3 ] );
1616 double aBadRate12 = getBadRate( &quad12, theCrit );
1617 double aBadRate13 = getBadRate( &quad13, theCrit );
1618 if ( aBadRate13 < aBadRate12 )
1625 // and remove fused elems and removed links from the maps
1626 mapEl_setLi.erase( tr1 );
1628 mapEl_setLi.erase( tr2 );
1629 mapLi_listEl.erase( *link12 );
1630 if(tr1->NbNodes()==3) {
1631 if( tr1->GetID() < tr2->GetID() ) {
1632 aMesh->ChangeElementNodes( tr1, n12, 4 );
1633 myLastCreatedElems.Append(tr1);
1634 aMesh->RemoveElement( tr2 );
1637 aMesh->ChangeElementNodes( tr2, n12, 4 );
1638 myLastCreatedElems.Append(tr2);
1639 aMesh->RemoveElement( tr1);
1643 const SMDS_MeshNode* N1 [6];
1644 const SMDS_MeshNode* N2 [6];
1645 GetNodesFromTwoTria(tr1,tr2,N1,N2);
1646 // now we receive following N1 and N2 (using numeration as above image)
1647 // tria1 : (1 2 4 5 9 7) and tria2 : (3 4 2 8 9 6)
1648 // i.e. first nodes from both arrays determ new diagonal
1649 const SMDS_MeshNode* aNodes[8];
1658 if( tr1->GetID() < tr2->GetID() ) {
1659 GetMeshDS()->ChangeElementNodes( tr1, aNodes, 8 );
1660 myLastCreatedElems.Append(tr1);
1661 GetMeshDS()->RemoveElement( tr2 );
1664 GetMeshDS()->ChangeElementNodes( tr2, aNodes, 8 );
1665 myLastCreatedElems.Append(tr2);
1666 GetMeshDS()->RemoveElement( tr1 );
1668 // remove middle node (9)
1669 GetMeshDS()->RemoveNode( N1[4] );
1673 mapEl_setLi.erase( tr3 );
1674 mapLi_listEl.erase( *link13 );
1675 if(tr1->NbNodes()==3) {
1676 if( tr1->GetID() < tr2->GetID() ) {
1677 aMesh->ChangeElementNodes( tr1, n13, 4 );
1678 myLastCreatedElems.Append(tr1);
1679 aMesh->RemoveElement( tr3 );
1682 aMesh->ChangeElementNodes( tr3, n13, 4 );
1683 myLastCreatedElems.Append(tr3);
1684 aMesh->RemoveElement( tr1 );
1688 const SMDS_MeshNode* N1 [6];
1689 const SMDS_MeshNode* N2 [6];
1690 GetNodesFromTwoTria(tr1,tr3,N1,N2);
1691 // now we receive following N1 and N2 (using numeration as above image)
1692 // tria1 : (1 2 4 5 9 7) and tria2 : (3 4 2 8 9 6)
1693 // i.e. first nodes from both arrays determ new diagonal
1694 const SMDS_MeshNode* aNodes[8];
1703 if( tr1->GetID() < tr2->GetID() ) {
1704 GetMeshDS()->ChangeElementNodes( tr1, aNodes, 8 );
1705 myLastCreatedElems.Append(tr1);
1706 GetMeshDS()->RemoveElement( tr3 );
1709 GetMeshDS()->ChangeElementNodes( tr3, aNodes, 8 );
1710 myLastCreatedElems.Append(tr3);
1711 GetMeshDS()->RemoveElement( tr1 );
1713 // remove middle node (9)
1714 GetMeshDS()->RemoveNode( N1[4] );
1718 // Next element to fuse: the rejected one
1720 startElem = Ok12 ? tr3 : tr2;
1722 } // if ( startElem )
1723 } // while ( startElem || !startLinks.empty() )
1724 } // while ( ! mapEl_setLi.empty() )
1730 /*#define DUMPSO(txt) \
1731 // cout << txt << endl;
1732 //=============================================================================
1736 //=============================================================================
1737 static void swap( int i1, int i2, int idNodes[], gp_Pnt P[] )
1741 int tmp = idNodes[ i1 ];
1742 idNodes[ i1 ] = idNodes[ i2 ];
1743 idNodes[ i2 ] = tmp;
1744 gp_Pnt Ptmp = P[ i1 ];
1747 DUMPSO( i1 << "(" << idNodes[ i2 ] << ") <-> " << i2 << "(" << idNodes[ i1 ] << ")");
1750 //=======================================================================
1751 //function : SortQuadNodes
1752 //purpose : Set 4 nodes of a quadrangle face in a good order.
1753 // Swap 1<->2 or 2<->3 nodes and correspondingly return
1755 //=======================================================================
1757 int SMESH_MeshEditor::SortQuadNodes (const SMDS_Mesh * theMesh,
1762 for ( i = 0; i < 4; i++ ) {
1763 const SMDS_MeshNode *n = theMesh->FindNode( idNodes[i] );
1765 P[ i ].SetCoord( n->X(), n->Y(), n->Z() );
1768 gp_Vec V1(P[0], P[1]);
1769 gp_Vec V2(P[0], P[2]);
1770 gp_Vec V3(P[0], P[3]);
1772 gp_Vec Cross1 = V1 ^ V2;
1773 gp_Vec Cross2 = V2 ^ V3;
1776 if (Cross1.Dot(Cross2) < 0)
1781 if (Cross1.Dot(Cross2) < 0)
1785 swap ( i, i + 1, idNodes, P );
1787 // for ( int ii = 0; ii < 4; ii++ ) {
1788 // const SMDS_MeshNode *n = theMesh->FindNode( idNodes[ii] );
1789 // DUMPSO( ii << "(" << idNodes[ii] <<") : "<<n->X()<<" "<<n->Y()<<" "<<n->Z());
1795 //=======================================================================
1796 //function : SortHexaNodes
1797 //purpose : Set 8 nodes of a hexahedron in a good order.
1798 // Return success status
1799 //=======================================================================
1801 bool SMESH_MeshEditor::SortHexaNodes (const SMDS_Mesh * theMesh,
1806 DUMPSO( "INPUT: ========================================");
1807 for ( i = 0; i < 8; i++ ) {
1808 const SMDS_MeshNode *n = theMesh->FindNode( idNodes[i] );
1809 if ( !n ) return false;
1810 P[ i ].SetCoord( n->X(), n->Y(), n->Z() );
1811 DUMPSO( i << "(" << idNodes[i] <<") : "<<n->X()<<" "<<n->Y()<<" "<<n->Z());
1813 DUMPSO( "========================================");
1816 set<int> faceNodes; // ids of bottom face nodes, to be found
1817 set<int> checkedId1; // ids of tried 2-nd nodes
1818 Standard_Real leastDist = DBL_MAX; // dist of the 4-th node from 123 plane
1819 const Standard_Real tol = 1.e-6; // tolerance to find nodes in plane
1820 int iMin, iLoop1 = 0;
1822 // Loop to try the 2-nd nodes
1824 while ( leastDist > DBL_MIN && ++iLoop1 < 8 )
1826 // Find not checked 2-nd node
1827 for ( i = 1; i < 8; i++ )
1828 if ( checkedId1.find( idNodes[i] ) == checkedId1.end() ) {
1829 int id1 = idNodes[i];
1830 swap ( 1, i, idNodes, P );
1831 checkedId1.insert ( id1 );
1835 // Find the 3-d node so that 1-2-3 triangle to be on a hexa face,
1836 // ie that all but meybe one (id3 which is on the same face) nodes
1837 // lay on the same side from the triangle plane.
1839 bool manyInPlane = false; // more than 4 nodes lay in plane
1841 while ( ++iLoop2 < 6 ) {
1843 // get 1-2-3 plane coeffs
1844 Standard_Real A, B, C, D;
1845 gp_Vec N = gp_Vec (P[0], P[1]).Crossed( gp_Vec (P[0], P[2]) );
1846 if ( N.SquareMagnitude() > gp::Resolution() )
1848 gp_Pln pln ( P[0], N );
1849 pln.Coefficients( A, B, C, D );
1851 // find the node (iMin) closest to pln
1852 Standard_Real dist[ 8 ], minDist = DBL_MAX;
1854 for ( i = 3; i < 8; i++ ) {
1855 dist[i] = A * P[i].X() + B * P[i].Y() + C * P[i].Z() + D;
1856 if ( fabs( dist[i] ) < minDist ) {
1857 minDist = fabs( dist[i] );
1860 if ( fabs( dist[i] ) <= tol )
1861 idInPln.insert( idNodes[i] );
1864 // there should not be more than 4 nodes in bottom plane
1865 if ( idInPln.size() > 1 )
1867 DUMPSO( "### idInPln.size() = " << idInPln.size());
1868 // idInPlane does not contain the first 3 nodes
1869 if ( manyInPlane || idInPln.size() == 5)
1870 return false; // all nodes in one plane
1873 // set the 1-st node to be not in plane
1874 for ( i = 3; i < 8; i++ ) {
1875 if ( idInPln.find( idNodes[ i ] ) == idInPln.end() ) {
1876 DUMPSO( "### Reset 0-th node");
1877 swap( 0, i, idNodes, P );
1882 // reset to re-check second nodes
1883 leastDist = DBL_MAX;
1887 break; // from iLoop2;
1890 // check that the other 4 nodes are on the same side
1891 bool sameSide = true;
1892 bool isNeg = dist[ iMin == 3 ? 4 : 3 ] <= 0.;
1893 for ( i = 3; sameSide && i < 8; i++ ) {
1895 sameSide = ( isNeg == dist[i] <= 0.);
1898 // keep best solution
1899 if ( sameSide && minDist < leastDist ) {
1900 leastDist = minDist;
1902 faceNodes.insert( idNodes[ 1 ] );
1903 faceNodes.insert( idNodes[ 2 ] );
1904 faceNodes.insert( idNodes[ iMin ] );
1905 DUMPSO( "loop " << iLoop2 << " id2 " << idNodes[ 1 ] << " id3 " << idNodes[ 2 ]
1906 << " leastDist = " << leastDist);
1907 if ( leastDist <= DBL_MIN )
1912 // set next 3-d node to check
1913 int iNext = 2 + iLoop2;
1915 DUMPSO( "Try 2-nd");
1916 swap ( 2, iNext, idNodes, P );
1918 } // while ( iLoop2 < 6 )
1921 if ( faceNodes.empty() ) return false;
1923 // Put the faceNodes in proper places
1924 for ( i = 4; i < 8; i++ ) {
1925 if ( faceNodes.find( idNodes[ i ] ) != faceNodes.end() ) {
1926 // find a place to put
1928 while ( faceNodes.find( idNodes[ iTo ] ) != faceNodes.end() )
1930 DUMPSO( "Set faceNodes");
1931 swap ( iTo, i, idNodes, P );
1936 // Set nodes of the found bottom face in good order
1937 DUMPSO( " Found bottom face: ");
1938 i = SortQuadNodes( theMesh, idNodes );
1940 gp_Pnt Ptmp = P[ i ];
1945 // for ( int ii = 0; ii < 4; ii++ ) {
1946 // const SMDS_MeshNode *n = theMesh->FindNode( idNodes[ii] );
1947 // DUMPSO( ii << "(" << idNodes[ii] <<") : "<<n->X()<<" "<<n->Y()<<" "<<n->Z());
1950 // Gravity center of the top and bottom faces
1951 gp_Pnt aGCb = ( P[0].XYZ() + P[1].XYZ() + P[2].XYZ() + P[3].XYZ() ) / 4.;
1952 gp_Pnt aGCt = ( P[4].XYZ() + P[5].XYZ() + P[6].XYZ() + P[7].XYZ() ) / 4.;
1954 // Get direction from the bottom to the top face
1955 gp_Vec upDir ( aGCb, aGCt );
1956 Standard_Real upDirSize = upDir.Magnitude();
1957 if ( upDirSize <= gp::Resolution() ) return false;
1960 // Assure that the bottom face normal points up
1961 gp_Vec Nb = gp_Vec (P[0], P[1]).Crossed( gp_Vec (P[0], P[2]) );
1962 Nb += gp_Vec (P[0], P[2]).Crossed( gp_Vec (P[0], P[3]) );
1963 if ( Nb.Dot( upDir ) < 0 ) {
1964 DUMPSO( "Reverse bottom face");
1965 swap( 1, 3, idNodes, P );
1968 // Find 5-th node - the one closest to the 1-st among the last 4 nodes.
1969 Standard_Real minDist = DBL_MAX;
1970 for ( i = 4; i < 8; i++ ) {
1971 // projection of P[i] to the plane defined by P[0] and upDir
1972 gp_Pnt Pp = P[i].Translated( upDir * ( upDir.Dot( gp_Vec( P[i], P[0] ))));
1973 Standard_Real sqDist = P[0].SquareDistance( Pp );
1974 if ( sqDist < minDist ) {
1979 DUMPSO( "Set 4-th");
1980 swap ( 4, iMin, idNodes, P );
1982 // Set nodes of the top face in good order
1983 DUMPSO( "Sort top face");
1984 i = SortQuadNodes( theMesh, &idNodes[4] );
1987 gp_Pnt Ptmp = P[ i ];
1992 // Assure that direction of the top face normal is from the bottom face
1993 gp_Vec Nt = gp_Vec (P[4], P[5]).Crossed( gp_Vec (P[4], P[6]) );
1994 Nt += gp_Vec (P[4], P[6]).Crossed( gp_Vec (P[4], P[7]) );
1995 if ( Nt.Dot( upDir ) < 0 ) {
1996 DUMPSO( "Reverse top face");
1997 swap( 5, 7, idNodes, P );
2000 // DUMPSO( "OUTPUT: ========================================");
2001 // for ( i = 0; i < 8; i++ ) {
2002 // float *p = ugrid->GetPoint(idNodes[i]);
2003 // DUMPSO( i << "(" << idNodes[i] << ") : " << p[0] << " " << p[1] << " " << p[2]);
2009 //================================================================================
2011 * \brief Return nodes linked to the given one
2012 * \param theNode - the node
2013 * \param linkedNodes - the found nodes
2014 * \param type - the type of elements to check
2016 * Medium nodes are ignored
2018 //================================================================================
2020 void SMESH_MeshEditor::GetLinkedNodes( const SMDS_MeshNode* theNode,
2021 TIDSortedElemSet & linkedNodes,
2022 SMDSAbs_ElementType type )
2024 SMDS_ElemIteratorPtr elemIt = theNode->GetInverseElementIterator(type);
2025 while ( elemIt->more() )
2027 const SMDS_MeshElement* elem = elemIt->next();
2028 SMDS_ElemIteratorPtr nodeIt = elem->nodesIterator();
2029 if ( elem->GetType() == SMDSAbs_Volume )
2031 SMDS_VolumeTool vol( elem );
2032 while ( nodeIt->more() ) {
2033 const SMDS_MeshNode* n = cast2Node( nodeIt->next() );
2034 if ( theNode != n && vol.IsLinked( theNode, n ))
2035 linkedNodes.insert( n );
2040 for ( int i = 0; nodeIt->more(); ++i ) {
2041 const SMDS_MeshNode* n = cast2Node( nodeIt->next() );
2042 if ( n == theNode ) {
2043 int iBefore = i - 1;
2045 if ( elem->IsQuadratic() ) {
2046 int nb = elem->NbNodes() / 2;
2047 iAfter = SMESH_MesherHelper::WrapIndex( iAfter, nb );
2048 iBefore = SMESH_MesherHelper::WrapIndex( iBefore, nb );
2050 linkedNodes.insert( elem->GetNode( iAfter ));
2051 linkedNodes.insert( elem->GetNode( iBefore ));
2058 //=======================================================================
2059 //function : laplacianSmooth
2060 //purpose : pulls theNode toward the center of surrounding nodes directly
2061 // connected to that node along an element edge
2062 //=======================================================================
2064 void laplacianSmooth(const SMDS_MeshNode* theNode,
2065 const Handle(Geom_Surface)& theSurface,
2066 map< const SMDS_MeshNode*, gp_XY* >& theUVMap)
2068 // find surrounding nodes
2070 TIDSortedElemSet nodeSet;
2071 SMESH_MeshEditor::GetLinkedNodes( theNode, nodeSet, SMDSAbs_Face );
2073 // compute new coodrs
2075 double coord[] = { 0., 0., 0. };
2076 TIDSortedElemSet::iterator nodeSetIt = nodeSet.begin();
2077 for ( ; nodeSetIt != nodeSet.end(); nodeSetIt++ ) {
2078 const SMDS_MeshNode* node = cast2Node(*nodeSetIt);
2079 if ( theSurface.IsNull() ) { // smooth in 3D
2080 coord[0] += node->X();
2081 coord[1] += node->Y();
2082 coord[2] += node->Z();
2084 else { // smooth in 2D
2085 ASSERT( theUVMap.find( node ) != theUVMap.end() );
2086 gp_XY* uv = theUVMap[ node ];
2087 coord[0] += uv->X();
2088 coord[1] += uv->Y();
2091 int nbNodes = nodeSet.size();
2094 coord[0] /= nbNodes;
2095 coord[1] /= nbNodes;
2097 if ( !theSurface.IsNull() ) {
2098 ASSERT( theUVMap.find( theNode ) != theUVMap.end() );
2099 theUVMap[ theNode ]->SetCoord( coord[0], coord[1] );
2100 gp_Pnt p3d = theSurface->Value( coord[0], coord[1] );
2106 coord[2] /= nbNodes;
2110 const_cast< SMDS_MeshNode* >( theNode )->setXYZ(coord[0],coord[1],coord[2]);
2113 //=======================================================================
2114 //function : centroidalSmooth
2115 //purpose : pulls theNode toward the element-area-weighted centroid of the
2116 // surrounding elements
2117 //=======================================================================
2119 void centroidalSmooth(const SMDS_MeshNode* theNode,
2120 const Handle(Geom_Surface)& theSurface,
2121 map< const SMDS_MeshNode*, gp_XY* >& theUVMap)
2123 gp_XYZ aNewXYZ(0.,0.,0.);
2124 SMESH::Controls::Area anAreaFunc;
2125 double totalArea = 0.;
2130 SMDS_ElemIteratorPtr elemIt = theNode->GetInverseElementIterator(SMDSAbs_Face);
2131 while ( elemIt->more() )
2133 const SMDS_MeshElement* elem = elemIt->next();
2136 gp_XYZ elemCenter(0.,0.,0.);
2137 SMESH::Controls::TSequenceOfXYZ aNodePoints;
2138 SMDS_ElemIteratorPtr itN = elem->nodesIterator();
2139 int nn = elem->NbNodes();
2140 if(elem->IsQuadratic()) nn = nn/2;
2142 //while ( itN->more() ) {
2144 const SMDS_MeshNode* aNode = static_cast<const SMDS_MeshNode*>( itN->next() );
2146 gp_XYZ aP( aNode->X(), aNode->Y(), aNode->Z() );
2147 aNodePoints.push_back( aP );
2148 if ( !theSurface.IsNull() ) { // smooth in 2D
2149 ASSERT( theUVMap.find( aNode ) != theUVMap.end() );
2150 gp_XY* uv = theUVMap[ aNode ];
2151 aP.SetCoord( uv->X(), uv->Y(), 0. );
2155 double elemArea = anAreaFunc.GetValue( aNodePoints );
2156 totalArea += elemArea;
2158 aNewXYZ += elemCenter * elemArea;
2160 aNewXYZ /= totalArea;
2161 if ( !theSurface.IsNull() ) {
2162 theUVMap[ theNode ]->SetCoord( aNewXYZ.X(), aNewXYZ.Y() );
2163 aNewXYZ = theSurface->Value( aNewXYZ.X(), aNewXYZ.Y() ).XYZ();
2168 const_cast< SMDS_MeshNode* >( theNode )->setXYZ(aNewXYZ.X(),aNewXYZ.Y(),aNewXYZ.Z());
2171 //=======================================================================
2172 //function : getClosestUV
2173 //purpose : return UV of closest projection
2174 //=======================================================================
2176 static bool getClosestUV (Extrema_GenExtPS& projector,
2177 const gp_Pnt& point,
2180 projector.Perform( point );
2181 if ( projector.IsDone() ) {
2182 double u, v, minVal = DBL_MAX;
2183 for ( int i = projector.NbExt(); i > 0; i-- )
2184 if ( projector.Value( i ) < minVal ) {
2185 minVal = projector.Value( i );
2186 projector.Point( i ).Parameter( u, v );
2188 result.SetCoord( u, v );
2194 //=======================================================================
2196 //purpose : Smooth theElements during theNbIterations or until a worst
2197 // element has aspect ratio <= theTgtAspectRatio.
2198 // Aspect Ratio varies in range [1.0, inf].
2199 // If theElements is empty, the whole mesh is smoothed.
2200 // theFixedNodes contains additionally fixed nodes. Nodes built
2201 // on edges and boundary nodes are always fixed.
2202 //=======================================================================
2204 void SMESH_MeshEditor::Smooth (TIDSortedElemSet & theElems,
2205 set<const SMDS_MeshNode*> & theFixedNodes,
2206 const SmoothMethod theSmoothMethod,
2207 const int theNbIterations,
2208 double theTgtAspectRatio,
2211 myLastCreatedElems.Clear();
2212 myLastCreatedNodes.Clear();
2214 MESSAGE((theSmoothMethod==LAPLACIAN ? "LAPLACIAN" : "CENTROIDAL") << "--::Smooth()");
2216 if ( theTgtAspectRatio < 1.0 )
2217 theTgtAspectRatio = 1.0;
2219 const double disttol = 1.e-16;
2221 SMESH::Controls::AspectRatio aQualityFunc;
2223 SMESHDS_Mesh* aMesh = GetMeshDS();
2225 if ( theElems.empty() ) {
2226 // add all faces to theElems
2227 SMDS_FaceIteratorPtr fIt = aMesh->facesIterator();
2228 while ( fIt->more() ) {
2229 const SMDS_MeshElement* face = fIt->next();
2230 theElems.insert( face );
2233 // get all face ids theElems are on
2234 set< int > faceIdSet;
2235 TIDSortedElemSet::iterator itElem;
2237 for ( itElem = theElems.begin(); itElem != theElems.end(); itElem++ ) {
2238 int fId = FindShape( *itElem );
2239 // check that corresponding submesh exists and a shape is face
2241 faceIdSet.find( fId ) == faceIdSet.end() &&
2242 aMesh->MeshElements( fId )) {
2243 TopoDS_Shape F = aMesh->IndexToShape( fId );
2244 if ( !F.IsNull() && F.ShapeType() == TopAbs_FACE )
2245 faceIdSet.insert( fId );
2248 faceIdSet.insert( 0 ); // to smooth elements that are not on any TopoDS_Face
2250 // ===============================================
2251 // smooth elements on each TopoDS_Face separately
2252 // ===============================================
2254 set< int >::reverse_iterator fId = faceIdSet.rbegin(); // treate 0 fId at the end
2255 for ( ; fId != faceIdSet.rend(); ++fId ) {
2256 // get face surface and submesh
2257 Handle(Geom_Surface) surface;
2258 SMESHDS_SubMesh* faceSubMesh = 0;
2260 double fToler2 = 0, vPeriod = 0., uPeriod = 0., f,l;
2261 double u1 = 0, u2 = 0, v1 = 0, v2 = 0;
2262 bool isUPeriodic = false, isVPeriodic = false;
2264 face = TopoDS::Face( aMesh->IndexToShape( *fId ));
2265 surface = BRep_Tool::Surface( face );
2266 faceSubMesh = aMesh->MeshElements( *fId );
2267 fToler2 = BRep_Tool::Tolerance( face );
2268 fToler2 *= fToler2 * 10.;
2269 isUPeriodic = surface->IsUPeriodic();
2271 vPeriod = surface->UPeriod();
2272 isVPeriodic = surface->IsVPeriodic();
2274 uPeriod = surface->VPeriod();
2275 surface->Bounds( u1, u2, v1, v2 );
2277 // ---------------------------------------------------------
2278 // for elements on a face, find movable and fixed nodes and
2279 // compute UV for them
2280 // ---------------------------------------------------------
2281 bool checkBoundaryNodes = false;
2282 bool isQuadratic = false;
2283 set<const SMDS_MeshNode*> setMovableNodes;
2284 map< const SMDS_MeshNode*, gp_XY* > uvMap, uvMap2;
2285 list< gp_XY > listUV; // uvs the 2 uvMaps refer to
2286 list< const SMDS_MeshElement* > elemsOnFace;
2288 Extrema_GenExtPS projector;
2289 GeomAdaptor_Surface surfAdaptor;
2290 if ( !surface.IsNull() ) {
2291 surfAdaptor.Load( surface );
2292 projector.Initialize( surfAdaptor, 20,20, 1e-5,1e-5 );
2294 int nbElemOnFace = 0;
2295 itElem = theElems.begin();
2296 // loop on not yet smoothed elements: look for elems on a face
2297 while ( itElem != theElems.end() ) {
2298 if ( faceSubMesh && nbElemOnFace == faceSubMesh->NbElements() )
2299 break; // all elements found
2301 const SMDS_MeshElement* elem = *itElem;
2302 if ( !elem || elem->GetType() != SMDSAbs_Face || elem->NbNodes() < 3 ||
2303 ( faceSubMesh && !faceSubMesh->Contains( elem ))) {
2307 elemsOnFace.push_back( elem );
2308 theElems.erase( itElem++ );
2312 isQuadratic = elem->IsQuadratic();
2314 // get movable nodes of elem
2315 const SMDS_MeshNode* node;
2316 SMDS_TypeOfPosition posType;
2317 SMDS_ElemIteratorPtr itN = elem->nodesIterator();
2318 int nn = 0, nbn = elem->NbNodes();
2319 if(elem->IsQuadratic())
2321 while ( nn++ < nbn ) {
2322 node = static_cast<const SMDS_MeshNode*>( itN->next() );
2323 const SMDS_PositionPtr& pos = node->GetPosition();
2324 posType = pos.get() ? pos->GetTypeOfPosition() : SMDS_TOP_3DSPACE;
2325 if (posType != SMDS_TOP_EDGE &&
2326 posType != SMDS_TOP_VERTEX &&
2327 theFixedNodes.find( node ) == theFixedNodes.end())
2329 // check if all faces around the node are on faceSubMesh
2330 // because a node on edge may be bound to face
2331 SMDS_ElemIteratorPtr eIt = node->GetInverseElementIterator(SMDSAbs_Face);
2333 if ( faceSubMesh ) {
2334 while ( eIt->more() && all ) {
2335 const SMDS_MeshElement* e = eIt->next();
2336 all = faceSubMesh->Contains( e );
2340 setMovableNodes.insert( node );
2342 checkBoundaryNodes = true;
2344 if ( posType == SMDS_TOP_3DSPACE )
2345 checkBoundaryNodes = true;
2348 if ( surface.IsNull() )
2351 // get nodes to check UV
2352 list< const SMDS_MeshNode* > uvCheckNodes;
2353 itN = elem->nodesIterator();
2354 nn = 0; nbn = elem->NbNodes();
2355 if(elem->IsQuadratic())
2357 while ( nn++ < nbn ) {
2358 node = static_cast<const SMDS_MeshNode*>( itN->next() );
2359 if ( uvMap.find( node ) == uvMap.end() )
2360 uvCheckNodes.push_back( node );
2361 // add nodes of elems sharing node
2362 // SMDS_ElemIteratorPtr eIt = node->GetInverseElementIterator(SMDSAbs_Face);
2363 // while ( eIt->more() ) {
2364 // const SMDS_MeshElement* e = eIt->next();
2365 // if ( e != elem ) {
2366 // SMDS_ElemIteratorPtr nIt = e->nodesIterator();
2367 // while ( nIt->more() ) {
2368 // const SMDS_MeshNode* n =
2369 // static_cast<const SMDS_MeshNode*>( nIt->next() );
2370 // if ( uvMap.find( n ) == uvMap.end() )
2371 // uvCheckNodes.push_back( n );
2377 list< const SMDS_MeshNode* >::iterator n = uvCheckNodes.begin();
2378 for ( ; n != uvCheckNodes.end(); ++n ) {
2381 const SMDS_PositionPtr& pos = node->GetPosition();
2382 posType = pos.get() ? pos->GetTypeOfPosition() : SMDS_TOP_3DSPACE;
2384 switch ( posType ) {
2385 case SMDS_TOP_FACE: {
2386 SMDS_FacePosition* fPos = ( SMDS_FacePosition* ) pos.get();
2387 uv.SetCoord( fPos->GetUParameter(), fPos->GetVParameter() );
2390 case SMDS_TOP_EDGE: {
2391 TopoDS_Shape S = aMesh->IndexToShape( pos->GetShapeId() );
2392 Handle(Geom2d_Curve) pcurve;
2393 if ( !S.IsNull() && S.ShapeType() == TopAbs_EDGE )
2394 pcurve = BRep_Tool::CurveOnSurface( TopoDS::Edge( S ), face, f,l );
2395 if ( !pcurve.IsNull() ) {
2396 double u = (( SMDS_EdgePosition* ) pos.get() )->GetUParameter();
2397 uv = pcurve->Value( u ).XY();
2401 case SMDS_TOP_VERTEX: {
2402 TopoDS_Shape S = aMesh->IndexToShape( pos->GetShapeId() );
2403 if ( !S.IsNull() && S.ShapeType() == TopAbs_VERTEX )
2404 uv = BRep_Tool::Parameters( TopoDS::Vertex( S ), face ).XY();
2409 // check existing UV
2410 bool project = true;
2411 gp_Pnt pNode ( node->X(), node->Y(), node->Z() );
2412 double dist1 = DBL_MAX, dist2 = 0;
2413 if ( posType != SMDS_TOP_3DSPACE ) {
2414 dist1 = pNode.SquareDistance( surface->Value( uv.X(), uv.Y() ));
2415 project = dist1 > fToler2;
2417 if ( project ) { // compute new UV
2419 if ( !getClosestUV( projector, pNode, newUV )) {
2420 MESSAGE("Node Projection Failed " << node);
2424 newUV.SetX( ElCLib::InPeriod( newUV.X(), u1, u2 ));
2426 newUV.SetY( ElCLib::InPeriod( newUV.Y(), v1, v2 ));
2428 if ( posType != SMDS_TOP_3DSPACE )
2429 dist2 = pNode.SquareDistance( surface->Value( newUV.X(), newUV.Y() ));
2430 if ( dist2 < dist1 )
2434 // store UV in the map
2435 listUV.push_back( uv );
2436 uvMap.insert( make_pair( node, &listUV.back() ));
2438 } // loop on not yet smoothed elements
2440 if ( !faceSubMesh || nbElemOnFace != faceSubMesh->NbElements() )
2441 checkBoundaryNodes = true;
2443 // fix nodes on mesh boundary
2445 if ( checkBoundaryNodes ) {
2446 typedef pair<const SMDS_MeshNode*, const SMDS_MeshNode*> TLink;
2447 map< TLink, int > linkNbMap; // how many times a link encounters in elemsOnFace
2448 map< TLink, int >::iterator link_nb;
2449 // put all elements links to linkNbMap
2450 list< const SMDS_MeshElement* >::iterator elemIt = elemsOnFace.begin();
2451 for ( ; elemIt != elemsOnFace.end(); ++elemIt ) {
2452 const SMDS_MeshElement* elem = (*elemIt);
2453 int nbn = elem->NbNodes();
2454 if(elem->IsQuadratic())
2456 // loop on elem links: insert them in linkNbMap
2457 const SMDS_MeshNode* curNode, *prevNode = elem->GetNode( nbn );
2458 for ( int iN = 0; iN < nbn; ++iN ) {
2459 curNode = elem->GetNode( iN );
2461 if ( curNode < prevNode ) link = make_pair( curNode , prevNode );
2462 else link = make_pair( prevNode , curNode );
2464 link_nb = linkNbMap.find( link );
2465 if ( link_nb == linkNbMap.end() )
2466 linkNbMap.insert( make_pair ( link, 1 ));
2471 // remove nodes that are in links encountered only once from setMovableNodes
2472 for ( link_nb = linkNbMap.begin(); link_nb != linkNbMap.end(); ++link_nb ) {
2473 if ( link_nb->second == 1 ) {
2474 setMovableNodes.erase( link_nb->first.first );
2475 setMovableNodes.erase( link_nb->first.second );
2480 // -----------------------------------------------------
2481 // for nodes on seam edge, compute one more UV ( uvMap2 );
2482 // find movable nodes linked to nodes on seam and which
2483 // are to be smoothed using the second UV ( uvMap2 )
2484 // -----------------------------------------------------
2486 set<const SMDS_MeshNode*> nodesNearSeam; // to smooth using uvMap2
2487 if ( !surface.IsNull() ) {
2488 TopExp_Explorer eExp( face, TopAbs_EDGE );
2489 for ( ; eExp.More(); eExp.Next() ) {
2490 TopoDS_Edge edge = TopoDS::Edge( eExp.Current() );
2491 if ( !BRep_Tool::IsClosed( edge, face ))
2493 SMESHDS_SubMesh* sm = aMesh->MeshElements( edge );
2494 if ( !sm ) continue;
2495 // find out which parameter varies for a node on seam
2498 Handle(Geom2d_Curve) pcurve = BRep_Tool::CurveOnSurface( edge, face, f, l );
2499 if ( pcurve.IsNull() ) continue;
2500 uv1 = pcurve->Value( f );
2502 pcurve = BRep_Tool::CurveOnSurface( edge, face, f, l );
2503 if ( pcurve.IsNull() ) continue;
2504 uv2 = pcurve->Value( f );
2505 int iPar = Abs( uv1.X() - uv2.X() ) > Abs( uv1.Y() - uv2.Y() ) ? 1 : 2;
2507 if ( uv1.Coord( iPar ) > uv2.Coord( iPar )) {
2508 gp_Pnt2d tmp = uv1; uv1 = uv2; uv2 = tmp;
2510 // get nodes on seam and its vertices
2511 list< const SMDS_MeshNode* > seamNodes;
2512 SMDS_NodeIteratorPtr nSeamIt = sm->GetNodes();
2513 while ( nSeamIt->more() ) {
2514 const SMDS_MeshNode* node = nSeamIt->next();
2515 if ( !isQuadratic || !IsMedium( node ))
2516 seamNodes.push_back( node );
2518 TopExp_Explorer vExp( edge, TopAbs_VERTEX );
2519 for ( ; vExp.More(); vExp.Next() ) {
2520 sm = aMesh->MeshElements( vExp.Current() );
2522 nSeamIt = sm->GetNodes();
2523 while ( nSeamIt->more() )
2524 seamNodes.push_back( nSeamIt->next() );
2527 // loop on nodes on seam
2528 list< const SMDS_MeshNode* >::iterator noSeIt = seamNodes.begin();
2529 for ( ; noSeIt != seamNodes.end(); ++noSeIt ) {
2530 const SMDS_MeshNode* nSeam = *noSeIt;
2531 map< const SMDS_MeshNode*, gp_XY* >::iterator n_uv = uvMap.find( nSeam );
2532 if ( n_uv == uvMap.end() )
2535 n_uv->second->SetCoord( iPar, uv1.Coord( iPar ));
2536 // set the second UV
2537 listUV.push_back( *n_uv->second );
2538 listUV.back().SetCoord( iPar, uv2.Coord( iPar ));
2539 if ( uvMap2.empty() )
2540 uvMap2 = uvMap; // copy the uvMap contents
2541 uvMap2[ nSeam ] = &listUV.back();
2543 // collect movable nodes linked to ones on seam in nodesNearSeam
2544 SMDS_ElemIteratorPtr eIt = nSeam->GetInverseElementIterator(SMDSAbs_Face);
2545 while ( eIt->more() ) {
2546 const SMDS_MeshElement* e = eIt->next();
2547 int nbUseMap1 = 0, nbUseMap2 = 0;
2548 SMDS_ElemIteratorPtr nIt = e->nodesIterator();
2549 int nn = 0, nbn = e->NbNodes();
2550 if(e->IsQuadratic()) nbn = nbn/2;
2551 while ( nn++ < nbn )
2553 const SMDS_MeshNode* n =
2554 static_cast<const SMDS_MeshNode*>( nIt->next() );
2556 setMovableNodes.find( n ) == setMovableNodes.end() )
2558 // add only nodes being closer to uv2 than to uv1
2559 gp_Pnt pMid (0.5 * ( n->X() + nSeam->X() ),
2560 0.5 * ( n->Y() + nSeam->Y() ),
2561 0.5 * ( n->Z() + nSeam->Z() ));
2563 getClosestUV( projector, pMid, uv );
2564 if ( uv.Coord( iPar ) > uvMap[ n ]->Coord( iPar ) ) {
2565 nodesNearSeam.insert( n );
2571 // for centroidalSmooth all element nodes must
2572 // be on one side of a seam
2573 if ( theSmoothMethod == CENTROIDAL && nbUseMap1 && nbUseMap2 ) {
2574 SMDS_ElemIteratorPtr nIt = e->nodesIterator();
2576 while ( nn++ < nbn ) {
2577 const SMDS_MeshNode* n =
2578 static_cast<const SMDS_MeshNode*>( nIt->next() );
2579 setMovableNodes.erase( n );
2583 } // loop on nodes on seam
2584 } // loop on edge of a face
2585 } // if ( !face.IsNull() )
2587 if ( setMovableNodes.empty() ) {
2588 MESSAGE( "Face id : " << *fId << " - NO SMOOTHING: no nodes to move!!!");
2589 continue; // goto next face
2597 double maxRatio = -1., maxDisplacement = -1.;
2598 set<const SMDS_MeshNode*>::iterator nodeToMove;
2599 for ( it = 0; it < theNbIterations; it++ ) {
2600 maxDisplacement = 0.;
2601 nodeToMove = setMovableNodes.begin();
2602 for ( ; nodeToMove != setMovableNodes.end(); nodeToMove++ ) {
2603 const SMDS_MeshNode* node = (*nodeToMove);
2604 gp_XYZ aPrevPos ( node->X(), node->Y(), node->Z() );
2607 bool map2 = ( nodesNearSeam.find( node ) != nodesNearSeam.end() );
2608 if ( theSmoothMethod == LAPLACIAN )
2609 laplacianSmooth( node, surface, map2 ? uvMap2 : uvMap );
2611 centroidalSmooth( node, surface, map2 ? uvMap2 : uvMap );
2613 // node displacement
2614 gp_XYZ aNewPos ( node->X(), node->Y(), node->Z() );
2615 Standard_Real aDispl = (aPrevPos - aNewPos).SquareModulus();
2616 if ( aDispl > maxDisplacement )
2617 maxDisplacement = aDispl;
2619 // no node movement => exit
2620 //if ( maxDisplacement < 1.e-16 ) {
2621 if ( maxDisplacement < disttol ) {
2622 MESSAGE("-- no node movement --");
2626 // check elements quality
2628 list< const SMDS_MeshElement* >::iterator elemIt = elemsOnFace.begin();
2629 for ( ; elemIt != elemsOnFace.end(); ++elemIt ) {
2630 const SMDS_MeshElement* elem = (*elemIt);
2631 if ( !elem || elem->GetType() != SMDSAbs_Face )
2633 SMESH::Controls::TSequenceOfXYZ aPoints;
2634 if ( aQualityFunc.GetPoints( elem, aPoints )) {
2635 double aValue = aQualityFunc.GetValue( aPoints );
2636 if ( aValue > maxRatio )
2640 if ( maxRatio <= theTgtAspectRatio ) {
2641 MESSAGE("-- quality achived --");
2644 if (it+1 == theNbIterations) {
2645 MESSAGE("-- Iteration limit exceeded --");
2647 } // smoothing iterations
2649 MESSAGE(" Face id: " << *fId <<
2650 " Nb iterstions: " << it <<
2651 " Displacement: " << maxDisplacement <<
2652 " Aspect Ratio " << maxRatio);
2654 // ---------------------------------------
2655 // new nodes positions are computed,
2656 // record movement in DS and set new UV
2657 // ---------------------------------------
2658 nodeToMove = setMovableNodes.begin();
2659 for ( ; nodeToMove != setMovableNodes.end(); nodeToMove++ ) {
2660 SMDS_MeshNode* node = const_cast< SMDS_MeshNode* > (*nodeToMove);
2661 aMesh->MoveNode( node, node->X(), node->Y(), node->Z() );
2662 map< const SMDS_MeshNode*, gp_XY* >::iterator node_uv = uvMap.find( node );
2663 if ( node_uv != uvMap.end() ) {
2664 gp_XY* uv = node_uv->second;
2666 ( SMDS_PositionPtr( new SMDS_FacePosition( *fId, uv->X(), uv->Y() )));
2670 // move medium nodes of quadratic elements
2673 SMESH_MesherHelper helper( *GetMesh() );
2674 if ( !face.IsNull() )
2675 helper.SetSubShape( face );
2676 list< const SMDS_MeshElement* >::iterator elemIt = elemsOnFace.begin();
2677 for ( ; elemIt != elemsOnFace.end(); ++elemIt ) {
2678 const SMDS_QuadraticFaceOfNodes* QF =
2679 dynamic_cast<const SMDS_QuadraticFaceOfNodes*> (*elemIt);
2681 vector<const SMDS_MeshNode*> Ns;
2682 Ns.reserve(QF->NbNodes()+1);
2683 SMDS_NodeIteratorPtr anIter = QF->interlacedNodesIterator();
2684 while ( anIter->more() )
2685 Ns.push_back( anIter->next() );
2686 Ns.push_back( Ns[0] );
2688 for(int i=0; i<QF->NbNodes(); i=i+2) {
2689 if ( !surface.IsNull() ) {
2690 gp_XY uv1 = helper.GetNodeUV( face, Ns[i], Ns[i+2] );
2691 gp_XY uv2 = helper.GetNodeUV( face, Ns[i+2], Ns[i] );
2692 gp_XY uv = ( uv1 + uv2 ) / 2.;
2693 gp_Pnt xyz = surface->Value( uv.X(), uv.Y() );
2694 x = xyz.X(); y = xyz.Y(); z = xyz.Z();
2697 x = (Ns[i]->X() + Ns[i+2]->X())/2;
2698 y = (Ns[i]->Y() + Ns[i+2]->Y())/2;
2699 z = (Ns[i]->Z() + Ns[i+2]->Z())/2;
2701 if( fabs( Ns[i+1]->X() - x ) > disttol ||
2702 fabs( Ns[i+1]->Y() - y ) > disttol ||
2703 fabs( Ns[i+1]->Z() - z ) > disttol ) {
2704 // we have to move i+1 node
2705 aMesh->MoveNode( Ns[i+1], x, y, z );
2712 } // loop on face ids
2716 //=======================================================================
2717 //function : isReverse
2718 //purpose : Return true if normal of prevNodes is not co-directied with
2719 // gp_Vec(prevNodes[iNotSame],nextNodes[iNotSame]).
2720 // iNotSame is where prevNodes and nextNodes are different
2721 //=======================================================================
2723 static bool isReverse(vector<const SMDS_MeshNode*> prevNodes,
2724 vector<const SMDS_MeshNode*> nextNodes,
2728 int iBeforeNotSame = ( iNotSame == 0 ? nbNodes - 1 : iNotSame - 1 );
2729 int iAfterNotSame = ( iNotSame + 1 == nbNodes ? 0 : iNotSame + 1 );
2731 const SMDS_MeshNode* nB = prevNodes[ iBeforeNotSame ];
2732 const SMDS_MeshNode* nA = prevNodes[ iAfterNotSame ];
2733 const SMDS_MeshNode* nP = prevNodes[ iNotSame ];
2734 const SMDS_MeshNode* nN = nextNodes[ iNotSame ];
2736 gp_Pnt pB ( nB->X(), nB->Y(), nB->Z() );
2737 gp_Pnt pA ( nA->X(), nA->Y(), nA->Z() );
2738 gp_Pnt pP ( nP->X(), nP->Y(), nP->Z() );
2739 gp_Pnt pN ( nN->X(), nN->Y(), nN->Z() );
2741 gp_Vec vB ( pP, pB ), vA ( pP, pA ), vN ( pP, pN );
2743 return (vA ^ vB) * vN < 0.0;
2746 //=======================================================================
2748 * \brief Create elements by sweeping an element
2749 * \param elem - element to sweep
2750 * \param newNodesItVec - nodes generated from each node of the element
2751 * \param newElems - generated elements
2752 * \param nbSteps - number of sweeping steps
2753 * \param srcElements - to append elem for each generated element
2755 //=======================================================================
2757 void SMESH_MeshEditor::sweepElement(const SMDS_MeshElement* elem,
2758 const vector<TNodeOfNodeListMapItr> & newNodesItVec,
2759 list<const SMDS_MeshElement*>& newElems,
2761 SMESH_SequenceOfElemPtr& srcElements)
2763 SMESHDS_Mesh* aMesh = GetMeshDS();
2765 // Loop on elem nodes:
2766 // find new nodes and detect same nodes indices
2767 int nbNodes = elem->NbNodes();
2768 vector < list< const SMDS_MeshNode* >::const_iterator > itNN( nbNodes );
2769 vector<const SMDS_MeshNode*> prevNod( nbNodes );
2770 vector<const SMDS_MeshNode*> nextNod( nbNodes );
2771 vector<const SMDS_MeshNode*> midlNod( nbNodes );
2773 int iNode, nbSame = 0, iNotSameNode = 0, iSameNode = 0;
2774 vector<int> sames(nbNodes);
2776 //bool issimple[nbNodes];
2777 vector<bool> issimple(nbNodes);
2779 for ( iNode = 0; iNode < nbNodes; iNode++ ) {
2780 TNodeOfNodeListMapItr nnIt = newNodesItVec[ iNode ];
2781 const SMDS_MeshNode* node = nnIt->first;
2782 const list< const SMDS_MeshNode* > & listNewNodes = nnIt->second;
2783 if ( listNewNodes.empty() )
2786 if(listNewNodes.size()==nbSteps) {
2787 issimple[iNode] = true;
2790 issimple[iNode] = false;
2793 itNN[ iNode ] = listNewNodes.begin();
2794 prevNod[ iNode ] = node;
2795 nextNod[ iNode ] = listNewNodes.front();
2796 //cout<<"iNode="<<iNode<<endl;
2797 //cout<<" prevNod[iNode]="<< prevNod[iNode]<<" nextNod[iNode]="<< nextNod[iNode]<<endl;
2798 if ( prevNod[ iNode ] != nextNod [ iNode ])
2799 iNotSameNode = iNode;
2803 sames[nbSame++] = iNode;
2806 //cout<<"1 nbSame="<<nbSame<<endl;
2807 if ( nbSame == nbNodes || nbSame > 2) {
2808 MESSAGE( " Too many same nodes of element " << elem->GetID() );
2812 // if( elem->IsQuadratic() && nbSame>0 ) {
2813 // MESSAGE( "Can not rotate quadratic element " << elem->GetID() );
2817 int iBeforeSame = 0, iAfterSame = 0, iOpposSame = 0;
2819 iBeforeSame = ( iSameNode == 0 ? nbNodes - 1 : iSameNode - 1 );
2820 iAfterSame = ( iSameNode + 1 == nbNodes ? 0 : iSameNode + 1 );
2821 iOpposSame = ( iSameNode - 2 < 0 ? iSameNode + 2 : iSameNode - 2 );
2825 //cout<<" prevNod[0]="<< prevNod[0]<<" prevNod[1]="<< prevNod[1]
2826 // <<" prevNod[2]="<< prevNod[2]<<" prevNod[3]="<< prevNod[4]
2827 // <<" prevNod[4]="<< prevNod[4]<<" prevNod[5]="<< prevNod[5]
2828 // <<" prevNod[6]="<< prevNod[6]<<" prevNod[7]="<< prevNod[7]<<endl;
2830 // check element orientation
2832 if ( nbNodes > 2 && !isReverse( prevNod, nextNod, nbNodes, iNotSameNode )) {
2833 //MESSAGE("Reversed elem " << elem );
2837 int iAB = iAfterSame + iBeforeSame;
2838 iBeforeSame = iAB - iBeforeSame;
2839 iAfterSame = iAB - iAfterSame;
2843 // make new elements
2844 for (int iStep = 0; iStep < nbSteps; iStep++ ) {
2846 for ( iNode = 0; iNode < nbNodes; iNode++ ) {
2847 if(issimple[iNode]) {
2848 nextNod[ iNode ] = *itNN[ iNode ];
2852 if( elem->GetType()==SMDSAbs_Node ) {
2853 // we have to use two nodes
2854 midlNod[ iNode ] = *itNN[ iNode ];
2856 nextNod[ iNode ] = *itNN[ iNode ];
2859 else if(!elem->IsQuadratic() ||
2860 elem->IsQuadratic() && elem->IsMediumNode(prevNod[iNode]) ) {
2861 // we have to use each second node
2863 nextNod[ iNode ] = *itNN[ iNode ];
2867 // we have to use two nodes
2868 midlNod[ iNode ] = *itNN[ iNode ];
2870 nextNod[ iNode ] = *itNN[ iNode ];
2875 SMDS_MeshElement* aNewElem = 0;
2876 if(!elem->IsPoly()) {
2877 switch ( nbNodes ) {
2881 if ( nbSame == 0 ) {
2883 aNewElem = aMesh->AddEdge( prevNod[ 0 ], nextNod[ 0 ] );
2885 aNewElem = aMesh->AddEdge( prevNod[ 0 ], nextNod[ 0 ], midlNod[ 0 ] );
2891 aNewElem = aMesh->AddFace(prevNod[ 0 ], prevNod[ 1 ],
2892 nextNod[ 1 ], nextNod[ 0 ] );
2894 aNewElem = aMesh->AddFace(prevNod[ 0 ], prevNod[ 1 ],
2895 nextNod[ iNotSameNode ] );
2899 case 3: { // TRIANGLE or quadratic edge
2900 if(elem->GetType() == SMDSAbs_Face) { // TRIANGLE
2902 if ( nbSame == 0 ) // --- pentahedron
2903 aNewElem = aMesh->AddVolume (prevNod[ i0 ], prevNod[ 1 ], prevNod[ i2 ],
2904 nextNod[ i0 ], nextNod[ 1 ], nextNod[ i2 ] );
2906 else if ( nbSame == 1 ) // --- pyramid
2907 aNewElem = aMesh->AddVolume (prevNod[ iBeforeSame ], prevNod[ iAfterSame ],
2908 nextNod[ iAfterSame ], nextNod[ iBeforeSame ],
2909 nextNod[ iSameNode ]);
2911 else // 2 same nodes: --- tetrahedron
2912 aNewElem = aMesh->AddVolume (prevNod[ i0 ], prevNod[ 1 ], prevNod[ i2 ],
2913 nextNod[ iNotSameNode ]);
2915 else { // quadratic edge
2916 if(nbSame==0) { // quadratic quadrangle
2917 aNewElem = aMesh->AddFace(prevNod[0], nextNod[0], nextNod[1], prevNod[1],
2918 midlNod[0], nextNod[2], midlNod[1], prevNod[2]);
2920 else if(nbSame==1) { // quadratic triangle
2922 return; // medium node on axis
2923 else if(sames[0]==0) {
2924 aNewElem = aMesh->AddFace(prevNod[0], nextNod[1], prevNod[1],
2925 nextNod[2], midlNod[1], prevNod[2]);
2927 else { // sames[0]==1
2928 aNewElem = aMesh->AddFace(prevNod[0], nextNod[0], prevNod[1],
2929 midlNod[0], nextNod[2], prevNod[2]);
2937 case 4: { // QUADRANGLE
2939 if ( nbSame == 0 ) // --- hexahedron
2940 aNewElem = aMesh->AddVolume (prevNod[ i0 ], prevNod[ 1 ], prevNod[ i2 ], prevNod[ 3 ],
2941 nextNod[ i0 ], nextNod[ 1 ], nextNod[ i2 ], nextNod[ 3 ]);
2943 else if ( nbSame == 1 ) { // --- pyramid + pentahedron
2944 aNewElem = aMesh->AddVolume (prevNod[ iBeforeSame ], prevNod[ iAfterSame ],
2945 nextNod[ iAfterSame ], nextNod[ iBeforeSame ],
2946 nextNod[ iSameNode ]);
2947 newElems.push_back( aNewElem );
2948 aNewElem = aMesh->AddVolume (prevNod[ iAfterSame ], prevNod[ iOpposSame ],
2949 prevNod[ iBeforeSame ], nextNod[ iAfterSame ],
2950 nextNod[ iOpposSame ], nextNod[ iBeforeSame ] );
2952 else if ( nbSame == 2 ) { // pentahedron
2953 if ( prevNod[ iBeforeSame ] == nextNod[ iBeforeSame ] )
2954 // iBeforeSame is same too
2955 aNewElem = aMesh->AddVolume (prevNod[ iBeforeSame ], prevNod[ iOpposSame ],
2956 nextNod[ iOpposSame ], prevNod[ iSameNode ],
2957 prevNod[ iAfterSame ], nextNod[ iAfterSame ]);
2959 // iAfterSame is same too
2960 aNewElem = aMesh->AddVolume (prevNod[ iSameNode ], prevNod[ iBeforeSame ],
2961 nextNod[ iBeforeSame ], prevNod[ iAfterSame ],
2962 prevNod[ iOpposSame ], nextNod[ iOpposSame ]);
2966 case 6: { // quadratic triangle
2967 // create pentahedron with 15 nodes
2968 if(i0>0) { // reversed case
2969 aNewElem = aMesh->AddVolume (prevNod[0], prevNod[2], prevNod[1],
2970 nextNod[0], nextNod[2], nextNod[1],
2971 prevNod[5], prevNod[4], prevNod[3],
2972 nextNod[5], nextNod[4], nextNod[3],
2973 midlNod[0], midlNod[2], midlNod[1]);
2975 else { // not reversed case
2976 aNewElem = aMesh->AddVolume (prevNod[0], prevNod[1], prevNod[2],
2977 nextNod[0], nextNod[1], nextNod[2],
2978 prevNod[3], prevNod[4], prevNod[5],
2979 nextNod[3], nextNod[4], nextNod[5],
2980 midlNod[0], midlNod[1], midlNod[2]);
2984 case 8: { // quadratic quadrangle
2985 // create hexahedron with 20 nodes
2986 if(i0>0) { // reversed case
2987 aNewElem = aMesh->AddVolume (prevNod[0], prevNod[3], prevNod[2], prevNod[1],
2988 nextNod[0], nextNod[3], nextNod[2], nextNod[1],
2989 prevNod[7], prevNod[6], prevNod[5], prevNod[4],
2990 nextNod[7], nextNod[6], nextNod[5], nextNod[4],
2991 midlNod[0], midlNod[3], midlNod[2], midlNod[1]);
2993 else { // not reversed case
2994 aNewElem = aMesh->AddVolume (prevNod[0], prevNod[1], prevNod[2], prevNod[3],
2995 nextNod[0], nextNod[1], nextNod[2], nextNod[3],
2996 prevNod[4], prevNod[5], prevNod[6], prevNod[7],
2997 nextNod[4], nextNod[5], nextNod[6], nextNod[7],
2998 midlNod[0], midlNod[1], midlNod[2], midlNod[3]);
3003 // realized for extrusion only
3004 //vector<const SMDS_MeshNode*> polyedre_nodes (nbNodes*2 + 4*nbNodes);
3005 //vector<int> quantities (nbNodes + 2);
3007 //quantities[0] = nbNodes; // bottom of prism
3008 //for (int inode = 0; inode < nbNodes; inode++) {
3009 // polyedre_nodes[inode] = prevNod[inode];
3012 //quantities[1] = nbNodes; // top of prism
3013 //for (int inode = 0; inode < nbNodes; inode++) {
3014 // polyedre_nodes[nbNodes + inode] = nextNod[inode];
3017 //for (int iface = 0; iface < nbNodes; iface++) {
3018 // quantities[iface + 2] = 4;
3019 // int inextface = (iface == nbNodes - 1) ? 0 : iface + 1;
3020 // polyedre_nodes[2*nbNodes + 4*iface + 0] = prevNod[iface];
3021 // polyedre_nodes[2*nbNodes + 4*iface + 1] = prevNod[inextface];
3022 // polyedre_nodes[2*nbNodes + 4*iface + 2] = nextNod[inextface];
3023 // polyedre_nodes[2*nbNodes + 4*iface + 3] = nextNod[iface];
3025 //aNewElem = aMesh->AddPolyhedralVolume (polyedre_nodes, quantities);
3032 // realized for extrusion only
3033 vector<const SMDS_MeshNode*> polyedre_nodes (nbNodes*2 + 4*nbNodes);
3034 vector<int> quantities (nbNodes + 2);
3036 quantities[0] = nbNodes; // bottom of prism
3037 for (int inode = 0; inode < nbNodes; inode++) {
3038 polyedre_nodes[inode] = prevNod[inode];
3041 quantities[1] = nbNodes; // top of prism
3042 for (int inode = 0; inode < nbNodes; inode++) {
3043 polyedre_nodes[nbNodes + inode] = nextNod[inode];
3046 for (int iface = 0; iface < nbNodes; iface++) {
3047 quantities[iface + 2] = 4;
3048 int inextface = (iface == nbNodes - 1) ? 0 : iface + 1;
3049 polyedre_nodes[2*nbNodes + 4*iface + 0] = prevNod[iface];
3050 polyedre_nodes[2*nbNodes + 4*iface + 1] = prevNod[inextface];
3051 polyedre_nodes[2*nbNodes + 4*iface + 2] = nextNod[inextface];
3052 polyedre_nodes[2*nbNodes + 4*iface + 3] = nextNod[iface];
3054 aNewElem = aMesh->AddPolyhedralVolume (polyedre_nodes, quantities);
3058 newElems.push_back( aNewElem );
3059 myLastCreatedElems.Append(aNewElem);
3060 srcElements.Append( elem );
3063 // set new prev nodes
3064 for ( iNode = 0; iNode < nbNodes; iNode++ )
3065 prevNod[ iNode ] = nextNod[ iNode ];
3070 //=======================================================================
3072 * \brief Create 1D and 2D elements around swept elements
3073 * \param mapNewNodes - source nodes and ones generated from them
3074 * \param newElemsMap - source elements and ones generated from them
3075 * \param elemNewNodesMap - nodes generated from each node of each element
3076 * \param elemSet - all swept elements
3077 * \param nbSteps - number of sweeping steps
3078 * \param srcElements - to append elem for each generated element
3080 //=======================================================================
3082 void SMESH_MeshEditor::makeWalls (TNodeOfNodeListMap & mapNewNodes,
3083 TElemOfElemListMap & newElemsMap,
3084 TElemOfVecOfNnlmiMap & elemNewNodesMap,
3085 TIDSortedElemSet& elemSet,
3087 SMESH_SequenceOfElemPtr& srcElements)
3089 ASSERT( newElemsMap.size() == elemNewNodesMap.size() );
3090 SMESHDS_Mesh* aMesh = GetMeshDS();
3092 // Find nodes belonging to only one initial element - sweep them to get edges.
3094 TNodeOfNodeListMapItr nList = mapNewNodes.begin();
3095 for ( ; nList != mapNewNodes.end(); nList++ ) {
3096 const SMDS_MeshNode* node =
3097 static_cast<const SMDS_MeshNode*>( nList->first );
3098 SMDS_ElemIteratorPtr eIt = node->GetInverseElementIterator();
3099 int nbInitElems = 0;
3100 const SMDS_MeshElement* el = 0;
3101 SMDSAbs_ElementType highType = SMDSAbs_Edge; // count most complex elements only
3102 while ( eIt->more() && nbInitElems < 2 ) {
3104 SMDSAbs_ElementType type = el->GetType();
3105 if ( type == SMDSAbs_Volume || type < highType ) continue;
3106 if ( type > highType ) {
3110 if ( elemSet.find(el) != elemSet.end() )
3113 if ( nbInitElems < 2 ) {
3114 bool NotCreateEdge = el && el->IsQuadratic() && el->IsMediumNode(node);
3115 if(!NotCreateEdge) {
3116 vector<TNodeOfNodeListMapItr> newNodesItVec( 1, nList );
3117 list<const SMDS_MeshElement*> newEdges;
3118 sweepElement( node, newNodesItVec, newEdges, nbSteps, srcElements );
3123 // Make a ceiling for each element ie an equal element of last new nodes.
3124 // Find free links of faces - make edges and sweep them into faces.
3126 TElemOfElemListMap::iterator itElem = newElemsMap.begin();
3127 TElemOfVecOfNnlmiMap::iterator itElemNodes = elemNewNodesMap.begin();
3128 for ( ; itElem != newElemsMap.end(); itElem++, itElemNodes++ ) {
3129 const SMDS_MeshElement* elem = itElem->first;
3130 vector<TNodeOfNodeListMapItr>& vecNewNodes = itElemNodes->second;
3132 if ( elem->GetType() == SMDSAbs_Edge ) {
3133 // create a ceiling edge
3134 if (!elem->IsQuadratic()) {
3135 if ( !aMesh->FindEdge( vecNewNodes[ 0 ]->second.back(),
3136 vecNewNodes[ 1 ]->second.back())) {
3137 myLastCreatedElems.Append(aMesh->AddEdge(vecNewNodes[ 0 ]->second.back(),
3138 vecNewNodes[ 1 ]->second.back()));
3139 srcElements.Append( myLastCreatedElems.Last() );
3143 if ( !aMesh->FindEdge( vecNewNodes[ 0 ]->second.back(),
3144 vecNewNodes[ 1 ]->second.back(),
3145 vecNewNodes[ 2 ]->second.back())) {
3146 myLastCreatedElems.Append(aMesh->AddEdge(vecNewNodes[ 0 ]->second.back(),
3147 vecNewNodes[ 1 ]->second.back(),
3148 vecNewNodes[ 2 ]->second.back()));
3149 srcElements.Append( myLastCreatedElems.Last() );
3153 if ( elem->GetType() != SMDSAbs_Face )
3156 if(itElem->second.size()==0) continue;
3158 bool hasFreeLinks = false;
3160 TIDSortedElemSet avoidSet;
3161 avoidSet.insert( elem );
3163 set<const SMDS_MeshNode*> aFaceLastNodes;
3164 int iNode, nbNodes = vecNewNodes.size();
3165 if(!elem->IsQuadratic()) {
3166 // loop on the face nodes
3167 for ( iNode = 0; iNode < nbNodes; iNode++ ) {
3168 aFaceLastNodes.insert( vecNewNodes[ iNode ]->second.back() );
3169 // look for free links of the face
3170 int iNext = ( iNode + 1 == nbNodes ) ? 0 : iNode + 1;
3171 const SMDS_MeshNode* n1 = vecNewNodes[ iNode ]->first;
3172 const SMDS_MeshNode* n2 = vecNewNodes[ iNext ]->first;
3173 // check if a link is free
3174 if ( ! SMESH_MeshEditor::FindFaceInSet ( n1, n2, elemSet, avoidSet )) {
3175 hasFreeLinks = true;
3176 // make an edge and a ceiling for a new edge
3177 if ( !aMesh->FindEdge( n1, n2 )) {
3178 myLastCreatedElems.Append(aMesh->AddEdge( n1, n2 )); // free link edge
3179 srcElements.Append( myLastCreatedElems.Last() );
3181 n1 = vecNewNodes[ iNode ]->second.back();
3182 n2 = vecNewNodes[ iNext ]->second.back();
3183 if ( !aMesh->FindEdge( n1, n2 )) {
3184 myLastCreatedElems.Append(aMesh->AddEdge( n1, n2 )); // ceiling edge
3185 srcElements.Append( myLastCreatedElems.Last() );
3190 else { // elem is quadratic face
3191 int nbn = nbNodes/2;
3192 for ( iNode = 0; iNode < nbn; iNode++ ) {
3193 aFaceLastNodes.insert( vecNewNodes[ iNode ]->second.back() );
3194 int iNext = ( iNode + 1 == nbn ) ? 0 : iNode + 1;
3195 const SMDS_MeshNode* n1 = vecNewNodes[ iNode ]->first;
3196 const SMDS_MeshNode* n2 = vecNewNodes[ iNext ]->first;
3197 // check if a link is free
3198 if ( ! SMESH_MeshEditor::FindFaceInSet ( n1, n2, elemSet, avoidSet )) {
3199 hasFreeLinks = true;
3200 // make an edge and a ceiling for a new edge
3202 const SMDS_MeshNode* n3 = vecNewNodes[ iNode+nbn ]->first;
3203 if ( !aMesh->FindEdge( n1, n2, n3 )) {
3204 myLastCreatedElems.Append(aMesh->AddEdge( n1, n2, n3 )); // free link edge
3205 srcElements.Append( myLastCreatedElems.Last() );
3207 n1 = vecNewNodes[ iNode ]->second.back();
3208 n2 = vecNewNodes[ iNext ]->second.back();
3209 n3 = vecNewNodes[ iNode+nbn ]->second.back();
3210 if ( !aMesh->FindEdge( n1, n2, n3 )) {
3211 myLastCreatedElems.Append(aMesh->AddEdge( n1, n2, n3 )); // ceiling edge
3212 srcElements.Append( myLastCreatedElems.Last() );
3216 for ( iNode = nbn; iNode < 2*nbn; iNode++ ) {
3217 aFaceLastNodes.insert( vecNewNodes[ iNode ]->second.back() );
3221 // sweep free links into faces
3223 if ( hasFreeLinks ) {
3224 list<const SMDS_MeshElement*> & newVolumes = itElem->second;
3225 int iVol, volNb, nbVolumesByStep = newVolumes.size() / nbSteps;
3227 set<const SMDS_MeshNode*> initNodeSet, topNodeSet, faceNodeSet;
3228 for ( iNode = 0; iNode < nbNodes; iNode++ ) {
3229 initNodeSet.insert( vecNewNodes[ iNode ]->first );
3230 topNodeSet .insert( vecNewNodes[ iNode ]->second.back() );
3232 for ( volNb = 0; volNb < nbVolumesByStep; volNb++ ) {
3233 list<const SMDS_MeshElement*>::iterator v = newVolumes.begin();
3235 while ( iVol++ < volNb ) v++;
3236 // find indices of free faces of a volume and their source edges
3237 list< int > freeInd;
3238 list< const SMDS_MeshElement* > srcEdges; // source edges of free faces
3239 SMDS_VolumeTool vTool( *v );
3240 int iF, nbF = vTool.NbFaces();
3241 for ( iF = 0; iF < nbF; iF ++ ) {
3242 if (vTool.IsFreeFace( iF ) &&
3243 vTool.GetFaceNodes( iF, faceNodeSet ) &&
3244 initNodeSet != faceNodeSet) // except an initial face
3246 if ( nbSteps == 1 && faceNodeSet == topNodeSet )
3248 freeInd.push_back( iF );
3249 // find source edge of a free face iF
3250 vector<const SMDS_MeshNode*> commonNodes; // shared by the initial and free faces
3251 commonNodes.resize( initNodeSet.size(), NULL ); // avoid spoiling memory
3252 std::set_intersection( faceNodeSet.begin(), faceNodeSet.end(),
3253 initNodeSet.begin(), initNodeSet.end(),
3254 commonNodes.begin());
3255 if ( (*v)->IsQuadratic() )
3256 srcEdges.push_back(aMesh->FindEdge (commonNodes[0],commonNodes[1],commonNodes[2]));
3258 srcEdges.push_back(aMesh->FindEdge (commonNodes[0],commonNodes[1]));
3260 if ( !srcEdges.back() )
3262 cout << "SMESH_MeshEditor::makeWalls(), no source edge found for a free face #"
3263 << iF << " of volume #" << vTool.ID() << endl;
3268 if ( freeInd.empty() )
3271 // create faces for all steps;
3272 // if such a face has been already created by sweep of edge,
3273 // assure that its orientation is OK
3274 for ( int iStep = 0; iStep < nbSteps; iStep++ ) {
3276 vTool.SetExternalNormal();
3277 list< int >::iterator ind = freeInd.begin();
3278 list< const SMDS_MeshElement* >::iterator srcEdge = srcEdges.begin();
3279 for ( ; ind != freeInd.end(); ++ind, ++srcEdge ) // loop on free faces
3281 const SMDS_MeshNode** nodes = vTool.GetFaceNodes( *ind );
3282 int nbn = vTool.NbFaceNodes( *ind );
3284 case 3: { ///// triangle
3285 const SMDS_MeshFace * f = aMesh->FindFace( nodes[ 0 ], nodes[ 1 ], nodes[ 2 ]);
3287 myLastCreatedElems.Append(aMesh->AddFace( nodes[ 0 ], nodes[ 1 ], nodes[ 2 ] ));
3288 else if ( nodes[ 1 ] != f->GetNode( f->GetNodeIndex( nodes[ 0 ] ) + 1 ))
3289 aMesh->ChangeElementNodes( f, nodes, nbn );
3292 case 4: { ///// quadrangle
3293 const SMDS_MeshFace * f = aMesh->FindFace( nodes[ 0 ], nodes[ 1 ], nodes[ 2 ], nodes[ 3 ]);
3295 myLastCreatedElems.Append(aMesh->AddFace( nodes[ 0 ], nodes[ 1 ], nodes[ 2 ], nodes[ 3 ] ));
3296 else if ( nodes[ 1 ] != f->GetNode( f->GetNodeIndex( nodes[ 0 ] ) + 1 ))
3297 aMesh->ChangeElementNodes( f, nodes, nbn );
3301 if( (*v)->IsQuadratic() ) {
3302 if(nbn==6) { /////// quadratic triangle
3303 const SMDS_MeshFace * f = aMesh->FindFace( nodes[0], nodes[2], nodes[4],
3304 nodes[1], nodes[3], nodes[5] );
3306 myLastCreatedElems.Append(aMesh->AddFace(nodes[0], nodes[2], nodes[4],
3307 nodes[1], nodes[3], nodes[5]));
3308 else if ( nodes[ 2 ] != f->GetNode( f->GetNodeIndex( nodes[ 0 ] ) + 1 ))
3309 aMesh->ChangeElementNodes( f, nodes, nbn );
3311 else { /////// quadratic quadrangle
3312 const SMDS_MeshFace * f = aMesh->FindFace( nodes[0], nodes[2], nodes[4], nodes[6],
3313 nodes[1], nodes[3], nodes[5], nodes[7] );
3315 myLastCreatedElems.Append(aMesh->AddFace(nodes[0], nodes[2], nodes[4], nodes[6],
3316 nodes[1], nodes[3], nodes[5], nodes[7]));
3317 else if ( nodes[ 2 ] != f->GetNode( f->GetNodeIndex( nodes[ 0 ] ) + 1 ))
3318 aMesh->ChangeElementNodes( f, nodes, nbn );
3321 else { //////// polygon
3322 vector<const SMDS_MeshNode*> polygon_nodes ( nodes, &nodes[nbn] );
3323 const SMDS_MeshFace * f = aMesh->FindFace( polygon_nodes );
3325 myLastCreatedElems.Append(aMesh->AddPolygonalFace(polygon_nodes));
3326 else if ( nodes[ 1 ] != f->GetNode( f->GetNodeIndex( nodes[ 0 ] ) + 1 ))
3327 aMesh->ChangeElementNodes( f, nodes, nbn );
3330 while ( srcElements.Length() < myLastCreatedElems.Length() )
3331 srcElements.Append( *srcEdge );
3333 } // loop on free faces
3335 // go to the next volume
3337 while ( iVol++ < nbVolumesByStep ) v++;
3340 } // sweep free links into faces
3342 // Make a ceiling face with a normal external to a volume
3344 SMDS_VolumeTool lastVol( itElem->second.back() );
3346 int iF = lastVol.GetFaceIndex( aFaceLastNodes );
3348 lastVol.SetExternalNormal();
3349 const SMDS_MeshNode** nodes = lastVol.GetFaceNodes( iF );
3350 int nbn = lastVol.NbFaceNodes( iF );
3353 if (!hasFreeLinks ||
3354 !aMesh->FindFace( nodes[ 0 ], nodes[ 1 ], nodes[ 2 ]))
3355 myLastCreatedElems.Append(aMesh->AddFace( nodes[ 0 ], nodes[ 1 ], nodes[ 2 ] ));
3358 if (!hasFreeLinks ||
3359 !aMesh->FindFace( nodes[ 0 ], nodes[ 1 ], nodes[ 2 ], nodes[ 3 ]))
3360 myLastCreatedElems.Append(aMesh->AddFace( nodes[ 0 ], nodes[ 1 ], nodes[ 2 ], nodes[ 3 ] ));
3363 if(itElem->second.back()->IsQuadratic()) {
3365 if (!hasFreeLinks ||
3366 !aMesh->FindFace(nodes[0], nodes[2], nodes[4],
3367 nodes[1], nodes[3], nodes[5]) ) {
3368 myLastCreatedElems.Append(aMesh->AddFace(nodes[0], nodes[2], nodes[4],
3369 nodes[1], nodes[3], nodes[5]));
3373 if (!hasFreeLinks ||
3374 !aMesh->FindFace(nodes[0], nodes[2], nodes[4], nodes[6],
3375 nodes[1], nodes[3], nodes[5], nodes[7]) )
3376 myLastCreatedElems.Append(aMesh->AddFace(nodes[0], nodes[2], nodes[4], nodes[6],
3377 nodes[1], nodes[3], nodes[5], nodes[7]));
3381 vector<const SMDS_MeshNode*> polygon_nodes ( nodes, &nodes[nbn] );
3382 if (!hasFreeLinks || !aMesh->FindFace(polygon_nodes))
3383 myLastCreatedElems.Append(aMesh->AddPolygonalFace(polygon_nodes));
3387 while ( srcElements.Length() < myLastCreatedElems.Length() )
3388 srcElements.Append( myLastCreatedElems.Last() );
3390 } // loop on swept elements
3393 //=======================================================================
3394 //function : RotationSweep
3396 //=======================================================================
3398 SMESH_MeshEditor::PGroupIDs
3399 SMESH_MeshEditor::RotationSweep(TIDSortedElemSet & theElems,
3400 const gp_Ax1& theAxis,
3401 const double theAngle,
3402 const int theNbSteps,
3403 const double theTol,
3404 const bool theMakeGroups,
3405 const bool theMakeWalls)
3407 myLastCreatedElems.Clear();
3408 myLastCreatedNodes.Clear();
3410 // source elements for each generated one
3411 SMESH_SequenceOfElemPtr srcElems, srcNodes;
3413 MESSAGE( "RotationSweep()");
3415 aTrsf.SetRotation( theAxis, theAngle );
3417 aTrsf2.SetRotation( theAxis, theAngle/2. );
3419 gp_Lin aLine( theAxis );
3420 double aSqTol = theTol * theTol;
3422 SMESHDS_Mesh* aMesh = GetMeshDS();
3424 TNodeOfNodeListMap mapNewNodes;
3425 TElemOfVecOfNnlmiMap mapElemNewNodes;
3426 TElemOfElemListMap newElemsMap;
3429 TIDSortedElemSet::iterator itElem;
3430 for ( itElem = theElems.begin(); itElem != theElems.end(); itElem++ ) {
3431 const SMDS_MeshElement* elem = *itElem;
3432 if ( !elem || elem->GetType() == SMDSAbs_Volume )
3434 vector<TNodeOfNodeListMapItr> & newNodesItVec = mapElemNewNodes[ elem ];
3435 newNodesItVec.reserve( elem->NbNodes() );
3437 // loop on elem nodes
3438 SMDS_ElemIteratorPtr itN = elem->nodesIterator();
3439 while ( itN->more() )
3441 // check if a node has been already sweeped
3442 const SMDS_MeshNode* node = cast2Node( itN->next() );
3443 TNodeOfNodeListMapItr nIt = mapNewNodes.find( node );
3444 if ( nIt == mapNewNodes.end() ) {
3445 nIt = mapNewNodes.insert( make_pair( node, list<const SMDS_MeshNode*>() )).first;
3446 list<const SMDS_MeshNode*>& listNewNodes = nIt->second;
3449 gp_XYZ aXYZ( node->X(), node->Y(), node->Z() );
3451 aXYZ.Coord( coord[0], coord[1], coord[2] );
3452 bool isOnAxis = ( aLine.SquareDistance( aXYZ ) <= aSqTol );
3453 const SMDS_MeshNode * newNode = node;
3454 for ( int i = 0; i < theNbSteps; i++ ) {
3456 if( elem->IsQuadratic() && !elem->IsMediumNode(node) ) {
3458 aTrsf2.Transforms( coord[0], coord[1], coord[2] );
3459 //aTrsf.Transforms( coord[0], coord[1], coord[2] );
3460 newNode = aMesh->AddNode( coord[0], coord[1], coord[2] );
3461 myLastCreatedNodes.Append(newNode);
3462 srcNodes.Append( node );
3463 listNewNodes.push_back( newNode );
3464 aTrsf2.Transforms( coord[0], coord[1], coord[2] );
3465 //aTrsf.Transforms( coord[0], coord[1], coord[2] );
3468 aTrsf.Transforms( coord[0], coord[1], coord[2] );
3470 newNode = aMesh->AddNode( coord[0], coord[1], coord[2] );
3471 myLastCreatedNodes.Append(newNode);
3472 srcNodes.Append( node );
3474 listNewNodes.push_back( newNode );
3478 // if current elem is quadratic and current node is not medium
3479 // we have to check - may be it is needed to insert additional nodes
3480 if( elem->IsQuadratic() && !elem->IsMediumNode(node) ) {
3481 list< const SMDS_MeshNode* > & listNewNodes = nIt->second;
3482 if(listNewNodes.size()==theNbSteps) {
3483 listNewNodes.clear();
3485 gp_XYZ aXYZ( node->X(), node->Y(), node->Z() );
3487 aXYZ.Coord( coord[0], coord[1], coord[2] );
3488 const SMDS_MeshNode * newNode = node;
3489 for(int i = 0; i<theNbSteps; i++) {
3490 aTrsf2.Transforms( coord[0], coord[1], coord[2] );
3491 newNode = aMesh->AddNode( coord[0], coord[1], coord[2] );
3492 myLastCreatedNodes.Append(newNode);
3493 listNewNodes.push_back( newNode );
3494 srcNodes.Append( node );
3495 aTrsf2.Transforms( coord[0], coord[1], coord[2] );
3496 newNode = aMesh->AddNode( coord[0], coord[1], coord[2] );
3497 myLastCreatedNodes.Append(newNode);
3498 srcNodes.Append( node );
3499 listNewNodes.push_back( newNode );
3504 newNodesItVec.push_back( nIt );
3506 // make new elements
3507 sweepElement( elem, newNodesItVec, newElemsMap[elem], theNbSteps, srcElems );
3511 makeWalls( mapNewNodes, newElemsMap, mapElemNewNodes, theElems, theNbSteps, srcElems );
3513 PGroupIDs newGroupIDs;
3514 if ( theMakeGroups )
3515 newGroupIDs = generateGroups( srcNodes, srcElems, "rotated");
3521 //=======================================================================
3522 //function : CreateNode
3524 //=======================================================================
3525 const SMDS_MeshNode* SMESH_MeshEditor::CreateNode(const double x,
3528 const double tolnode,
3529 SMESH_SequenceOfNode& aNodes)
3531 myLastCreatedElems.Clear();
3532 myLastCreatedNodes.Clear();
3535 SMESHDS_Mesh * aMesh = myMesh->GetMeshDS();
3537 // try to search in sequence of existing nodes
3538 // if aNodes.Length()>0 we 'nave to use given sequence
3539 // else - use all nodes of mesh
3540 if(aNodes.Length()>0) {
3542 for(i=1; i<=aNodes.Length(); i++) {
3543 gp_Pnt P2(aNodes.Value(i)->X(),aNodes.Value(i)->Y(),aNodes.Value(i)->Z());
3544 if(P1.Distance(P2)<tolnode)
3545 return aNodes.Value(i);
3549 SMDS_NodeIteratorPtr itn = aMesh->nodesIterator();
3550 while(itn->more()) {
3551 const SMDS_MeshNode* aN = static_cast<const SMDS_MeshNode*> (itn->next());
3552 gp_Pnt P2(aN->X(),aN->Y(),aN->Z());
3553 if(P1.Distance(P2)<tolnode)
3558 // create new node and return it
3559 const SMDS_MeshNode* NewNode = aMesh->AddNode(x,y,z);
3560 myLastCreatedNodes.Append(NewNode);
3565 //=======================================================================
3566 //function : ExtrusionSweep
3568 //=======================================================================
3570 SMESH_MeshEditor::PGroupIDs
3571 SMESH_MeshEditor::ExtrusionSweep (TIDSortedElemSet & theElems,
3572 const gp_Vec& theStep,
3573 const int theNbSteps,
3574 TElemOfElemListMap& newElemsMap,
3575 const bool theMakeGroups,
3577 const double theTolerance)
3579 ExtrusParam aParams;
3580 aParams.myDir = gp_Dir(theStep);
3581 aParams.myNodes.Clear();
3582 aParams.mySteps = new TColStd_HSequenceOfReal;
3584 for(i=1; i<=theNbSteps; i++)
3585 aParams.mySteps->Append(theStep.Magnitude());
3588 ExtrusionSweep(theElems,aParams,newElemsMap,theMakeGroups,theFlags,theTolerance);
3592 //=======================================================================
3593 //function : ExtrusionSweep
3595 //=======================================================================
3597 SMESH_MeshEditor::PGroupIDs
3598 SMESH_MeshEditor::ExtrusionSweep (TIDSortedElemSet & theElems,
3599 ExtrusParam& theParams,
3600 TElemOfElemListMap& newElemsMap,
3601 const bool theMakeGroups,
3603 const double theTolerance)
3605 myLastCreatedElems.Clear();
3606 myLastCreatedNodes.Clear();
3608 // source elements for each generated one
3609 SMESH_SequenceOfElemPtr srcElems, srcNodes;
3611 SMESHDS_Mesh* aMesh = GetMeshDS();
3613 int nbsteps = theParams.mySteps->Length();
3615 TNodeOfNodeListMap mapNewNodes;
3616 //TNodeOfNodeVecMap mapNewNodes;
3617 TElemOfVecOfNnlmiMap mapElemNewNodes;
3618 //TElemOfVecOfMapNodesMap mapElemNewNodes;
3621 TIDSortedElemSet::iterator itElem;
3622 for ( itElem = theElems.begin(); itElem != theElems.end(); itElem++ ) {
3623 // check element type
3624 const SMDS_MeshElement* elem = *itElem;
3625 if ( !elem || elem->GetType() == SMDSAbs_Volume )
3628 vector<TNodeOfNodeListMapItr> & newNodesItVec = mapElemNewNodes[ elem ];
3629 //vector<TNodeOfNodeVecMapItr> & newNodesItVec = mapElemNewNodes[ elem ];
3630 newNodesItVec.reserve( elem->NbNodes() );
3632 // loop on elem nodes
3633 SMDS_ElemIteratorPtr itN = elem->nodesIterator();
3634 while ( itN->more() )
3636 // check if a node has been already sweeped
3637 const SMDS_MeshNode* node = cast2Node( itN->next() );
3638 TNodeOfNodeListMap::iterator nIt = mapNewNodes.find( node );
3639 //TNodeOfNodeVecMap::iterator nIt = mapNewNodes.find( node );
3640 if ( nIt == mapNewNodes.end() ) {
3641 nIt = mapNewNodes.insert( make_pair( node, list<const SMDS_MeshNode*>() )).first;
3642 //nIt = mapNewNodes.insert( make_pair( node, vector<const SMDS_MeshNode*>() )).first;
3643 list<const SMDS_MeshNode*>& listNewNodes = nIt->second;
3644 //vector<const SMDS_MeshNode*>& vecNewNodes = nIt->second;
3645 //vecNewNodes.reserve(nbsteps);
3648 double coord[] = { node->X(), node->Y(), node->Z() };
3649 //int nbsteps = theParams.mySteps->Length();
3650 for ( int i = 0; i < nbsteps; i++ ) {
3651 if( elem->IsQuadratic() && !elem->IsMediumNode(node) ) {
3652 // create additional node
3653 double x = coord[0] + theParams.myDir.X()*theParams.mySteps->Value(i+1)/2.;
3654 double y = coord[1] + theParams.myDir.Y()*theParams.mySteps->Value(i+1)/2.;
3655 double z = coord[2] + theParams.myDir.Z()*theParams.mySteps->Value(i+1)/2.;
3656 if( theFlags & EXTRUSION_FLAG_SEW ) {
3657 const SMDS_MeshNode * newNode = CreateNode(x, y, z,
3658 theTolerance, theParams.myNodes);
3659 listNewNodes.push_back( newNode );
3662 const SMDS_MeshNode * newNode = aMesh->AddNode(x, y, z);
3663 myLastCreatedNodes.Append(newNode);
3664 srcNodes.Append( node );
3665 listNewNodes.push_back( newNode );
3668 //aTrsf.Transforms( coord[0], coord[1], coord[2] );
3669 coord[0] = coord[0] + theParams.myDir.X()*theParams.mySteps->Value(i+1);
3670 coord[1] = coord[1] + theParams.myDir.Y()*theParams.mySteps->Value(i+1);
3671 coord[2] = coord[2] + theParams.myDir.Z()*theParams.mySteps->Value(i+1);
3672 if( theFlags & EXTRUSION_FLAG_SEW ) {
3673 const SMDS_MeshNode * newNode = CreateNode(coord[0], coord[1], coord[2],
3674 theTolerance, theParams.myNodes);
3675 listNewNodes.push_back( newNode );
3676 //vecNewNodes[i]=newNode;
3679 const SMDS_MeshNode * newNode = aMesh->AddNode( coord[0], coord[1], coord[2] );
3680 myLastCreatedNodes.Append(newNode);
3681 srcNodes.Append( node );
3682 listNewNodes.push_back( newNode );
3683 //vecNewNodes[i]=newNode;
3688 // if current elem is quadratic and current node is not medium
3689 // we have to check - may be it is needed to insert additional nodes
3690 if( elem->IsQuadratic() && !elem->IsMediumNode(node) ) {
3691 list< const SMDS_MeshNode* > & listNewNodes = nIt->second;
3692 if(listNewNodes.size()==nbsteps) {
3693 listNewNodes.clear();
3694 double coord[] = { node->X(), node->Y(), node->Z() };
3695 for ( int i = 0; i < nbsteps; i++ ) {
3696 double x = coord[0] + theParams.myDir.X()*theParams.mySteps->Value(i+1);
3697 double y = coord[1] + theParams.myDir.Y()*theParams.mySteps->Value(i+1);
3698 double z = coord[2] + theParams.myDir.Z()*theParams.mySteps->Value(i+1);
3699 if( theFlags & EXTRUSION_FLAG_SEW ) {
3700 const SMDS_MeshNode * newNode = CreateNode(x, y, z,
3701 theTolerance, theParams.myNodes);
3702 listNewNodes.push_back( newNode );
3705 const SMDS_MeshNode * newNode = aMesh->AddNode(x, y, z);
3706 myLastCreatedNodes.Append(newNode);
3707 srcNodes.Append( node );
3708 listNewNodes.push_back( newNode );
3710 coord[0] = coord[0] + theParams.myDir.X()*theParams.mySteps->Value(i+1);
3711 coord[1] = coord[1] + theParams.myDir.Y()*theParams.mySteps->Value(i+1);
3712 coord[2] = coord[2] + theParams.myDir.Z()*theParams.mySteps->Value(i+1);
3713 if( theFlags & EXTRUSION_FLAG_SEW ) {
3714 const SMDS_MeshNode * newNode = CreateNode(coord[0], coord[1], coord[2],
3715 theTolerance, theParams.myNodes);
3716 listNewNodes.push_back( newNode );
3719 const SMDS_MeshNode * newNode = aMesh->AddNode( coord[0], coord[1], coord[2] );
3720 myLastCreatedNodes.Append(newNode);
3721 srcNodes.Append( node );
3722 listNewNodes.push_back( newNode );
3728 newNodesItVec.push_back( nIt );
3730 // make new elements
3731 sweepElement( elem, newNodesItVec, newElemsMap[elem], nbsteps, srcElems );
3734 if( theFlags & EXTRUSION_FLAG_BOUNDARY ) {
3735 makeWalls( mapNewNodes, newElemsMap, mapElemNewNodes, theElems, nbsteps, srcElems );
3737 PGroupIDs newGroupIDs;
3738 if ( theMakeGroups )
3739 newGroupIDs = generateGroups( srcNodes, srcElems, "extruded");
3745 //=======================================================================
3746 //class : SMESH_MeshEditor_PathPoint
3747 //purpose : auxiliary class
3748 //=======================================================================
3749 class SMESH_MeshEditor_PathPoint {
3751 SMESH_MeshEditor_PathPoint() {
3752 myPnt.SetCoord(99., 99., 99.);
3753 myTgt.SetCoord(1.,0.,0.);
3757 void SetPnt(const gp_Pnt& aP3D){
3760 void SetTangent(const gp_Dir& aTgt){
3763 void SetAngle(const double& aBeta){
3766 void SetParameter(const double& aPrm){
3769 const gp_Pnt& Pnt()const{
3772 const gp_Dir& Tangent()const{
3775 double Angle()const{
3778 double Parameter()const{
3789 //=======================================================================
3790 //function : ExtrusionAlongTrack
3792 //=======================================================================
3793 SMESH_MeshEditor::Extrusion_Error
3794 SMESH_MeshEditor::ExtrusionAlongTrack (TIDSortedElemSet & theElements,
3795 SMESH_subMesh* theTrack,
3796 const SMDS_MeshNode* theN1,
3797 const bool theHasAngles,
3798 list<double>& theAngles,
3799 const bool theHasRefPoint,
3800 const gp_Pnt& theRefPoint,
3801 const bool theMakeGroups)
3803 myLastCreatedElems.Clear();
3804 myLastCreatedNodes.Clear();
3806 // source elements for each generated one
3807 SMESH_SequenceOfElemPtr srcElems, srcNodes;
3809 int j, aNbTP, aNbE, aNb;
3810 double aT1, aT2, aT, aAngle, aX, aY, aZ;
3811 std::list<double> aPrms;
3812 std::list<double>::iterator aItD;
3813 TIDSortedElemSet::iterator itElem;
3815 Standard_Real aTx1, aTx2, aL2, aTolVec, aTolVec2;
3819 Handle(Geom_Curve) aC3D;
3820 TopoDS_Edge aTrackEdge;
3821 TopoDS_Vertex aV1, aV2;
3823 SMDS_ElemIteratorPtr aItE;
3824 SMDS_NodeIteratorPtr aItN;
3825 SMDSAbs_ElementType aTypeE;
3827 TNodeOfNodeListMap mapNewNodes;
3828 TElemOfVecOfNnlmiMap mapElemNewNodes;
3829 TElemOfElemListMap newElemsMap;
3832 aTolVec2=aTolVec*aTolVec;
3835 aNbE = theElements.size();
3838 return EXTR_NO_ELEMENTS;
3840 // 1.1 Track Pattern
3843 SMESHDS_SubMesh* pSubMeshDS=theTrack->GetSubMeshDS();
3845 aItE = pSubMeshDS->GetElements();
3846 while ( aItE->more() ) {
3847 const SMDS_MeshElement* pE = aItE->next();
3848 aTypeE = pE->GetType();
3849 // Pattern must contain links only
3850 if ( aTypeE != SMDSAbs_Edge )
3851 return EXTR_PATH_NOT_EDGE;
3854 const TopoDS_Shape& aS = theTrack->GetSubShape();
3855 // Sub shape for the Pattern must be an Edge
3856 if ( aS.ShapeType() != TopAbs_EDGE )
3857 return EXTR_BAD_PATH_SHAPE;
3859 aTrackEdge = TopoDS::Edge( aS );
3860 // the Edge must not be degenerated
3861 if ( BRep_Tool::Degenerated( aTrackEdge ) )
3862 return EXTR_BAD_PATH_SHAPE;
3864 TopExp::Vertices( aTrackEdge, aV1, aV2 );
3865 aT1=BRep_Tool::Parameter( aV1, aTrackEdge );
3866 aT2=BRep_Tool::Parameter( aV2, aTrackEdge );
3868 aItN = theTrack->GetFather()->GetSubMesh( aV1 )->GetSubMeshDS()->GetNodes();
3869 const SMDS_MeshNode* aN1 = aItN->next();
3871 aItN = theTrack->GetFather()->GetSubMesh( aV2 )->GetSubMeshDS()->GetNodes();
3872 const SMDS_MeshNode* aN2 = aItN->next();
3874 // starting node must be aN1 or aN2
3875 if ( !( aN1 == theN1 || aN2 == theN1 ) )
3876 return EXTR_BAD_STARTING_NODE;
3878 aNbTP = pSubMeshDS->NbNodes() + 2;
3881 vector<double> aAngles( aNbTP );
3883 for ( j=0; j < aNbTP; ++j ) {
3887 if ( theHasAngles ) {
3888 aItD = theAngles.begin();
3889 for ( j=1; (aItD != theAngles.end()) && (j<aNbTP); ++aItD, ++j ) {
3891 aAngles[j] = aAngle;
3895 // 2. Collect parameters on the track edge
3896 aPrms.push_back( aT1 );
3897 aPrms.push_back( aT2 );
3899 aItN = pSubMeshDS->GetNodes();
3900 while ( aItN->more() ) {
3901 const SMDS_MeshNode* pNode = aItN->next();
3902 const SMDS_EdgePosition* pEPos =
3903 static_cast<const SMDS_EdgePosition*>( pNode->GetPosition().get() );
3904 aT = pEPos->GetUParameter();
3905 aPrms.push_back( aT );
3910 if ( aN1 == theN1 ) {
3922 SMESH_MeshEditor_PathPoint aPP;
3923 vector<SMESH_MeshEditor_PathPoint> aPPs( aNbTP );
3925 aC3D = BRep_Tool::Curve( aTrackEdge, aTx1, aTx2 );
3927 aItD = aPrms.begin();
3928 for ( j=0; aItD != aPrms.end(); ++aItD, ++j ) {
3930 aC3D->D1( aT, aP3D, aVec );
3931 aL2 = aVec.SquareMagnitude();
3932 if ( aL2 < aTolVec2 )
3933 return EXTR_CANT_GET_TANGENT;
3935 gp_Dir aTgt( aVec );
3936 aAngle = aAngles[j];
3939 aPP.SetTangent( aTgt );
3940 aPP.SetAngle( aAngle );
3941 aPP.SetParameter( aT );
3945 // 3. Center of rotation aV0
3947 if ( !theHasRefPoint ) {
3949 aGC.SetCoord( 0.,0.,0. );
3951 itElem = theElements.begin();
3952 for ( ; itElem != theElements.end(); itElem++ ) {
3953 const SMDS_MeshElement* elem = *itElem;
3955 SMDS_ElemIteratorPtr itN = elem->nodesIterator();
3956 while ( itN->more() ) {
3957 const SMDS_MeshNode* node = static_cast<const SMDS_MeshNode*>( itN->next() );
3962 if ( mapNewNodes.find( node ) == mapNewNodes.end() ) {
3963 list<const SMDS_MeshNode*> aLNx;
3964 mapNewNodes[node] = aLNx;
3966 gp_XYZ aXYZ( aX, aY, aZ );
3974 } // if (!theHasRefPoint) {
3975 mapNewNodes.clear();
3977 // 4. Processing the elements
3978 SMESHDS_Mesh* aMesh = GetMeshDS();
3980 for ( itElem = theElements.begin(); itElem != theElements.end(); itElem++ ) {
3981 // check element type
3982 const SMDS_MeshElement* elem = *itElem;
3983 aTypeE = elem->GetType();
3984 if ( !elem || ( aTypeE != SMDSAbs_Face && aTypeE != SMDSAbs_Edge ) )
3987 vector<TNodeOfNodeListMapItr> & newNodesItVec = mapElemNewNodes[ elem ];
3988 newNodesItVec.reserve( elem->NbNodes() );
3990 // loop on elem nodes
3992 SMDS_ElemIteratorPtr itN = elem->nodesIterator();
3993 while ( itN->more() )
3996 // check if a node has been already processed
3997 const SMDS_MeshNode* node =
3998 static_cast<const SMDS_MeshNode*>( itN->next() );
3999 TNodeOfNodeListMap::iterator nIt = mapNewNodes.find( node );
4000 if ( nIt == mapNewNodes.end() ) {
4001 nIt = mapNewNodes.insert( make_pair( node, list<const SMDS_MeshNode*>() )).first;
4002 list<const SMDS_MeshNode*>& listNewNodes = nIt->second;
4005 aX = node->X(); aY = node->Y(); aZ = node->Z();
4007 Standard_Real aAngle1x, aAngleT1T0, aTolAng;
4008 gp_Pnt aP0x, aP1x, aPN0, aPN1, aV0x, aV1x;
4009 gp_Ax1 anAx1, anAxT1T0;
4010 gp_Dir aDT1x, aDT0x, aDT1T0;
4015 aPN0.SetCoord(aX, aY, aZ);
4017 const SMESH_MeshEditor_PathPoint& aPP0 = aPPs[0];
4019 aDT0x= aPP0.Tangent();
4021 for ( j = 1; j < aNbTP; ++j ) {
4022 const SMESH_MeshEditor_PathPoint& aPP1 = aPPs[j];
4024 aDT1x = aPP1.Tangent();
4025 aAngle1x = aPP1.Angle();
4027 gp_Trsf aTrsf, aTrsfRot, aTrsfRotT1T0;
4029 gp_Vec aV01x( aP0x, aP1x );
4030 aTrsf.SetTranslation( aV01x );
4033 aV1x = aV0x.Transformed( aTrsf );
4034 aPN1 = aPN0.Transformed( aTrsf );
4036 // rotation 1 [ T1,T0 ]
4037 aAngleT1T0=-aDT1x.Angle( aDT0x );
4038 if (fabs(aAngleT1T0) > aTolAng) {
4040 anAxT1T0.SetLocation( aV1x );
4041 anAxT1T0.SetDirection( aDT1T0 );
4042 aTrsfRotT1T0.SetRotation( anAxT1T0, aAngleT1T0 );
4044 aPN1 = aPN1.Transformed( aTrsfRotT1T0 );
4048 if ( theHasAngles ) {
4049 anAx1.SetLocation( aV1x );
4050 anAx1.SetDirection( aDT1x );
4051 aTrsfRot.SetRotation( anAx1, aAngle1x );
4053 aPN1 = aPN1.Transformed( aTrsfRot );
4057 if( elem->IsQuadratic() && !elem->IsMediumNode(node) ) {
4058 // create additional node
4059 double x = ( aPN1.X() + aPN0.X() )/2.;
4060 double y = ( aPN1.Y() + aPN0.Y() )/2.;
4061 double z = ( aPN1.Z() + aPN0.Z() )/2.;
4062 const SMDS_MeshNode* newNode = aMesh->AddNode(x,y,z);
4063 myLastCreatedNodes.Append(newNode);
4064 srcNodes.Append( node );
4065 listNewNodes.push_back( newNode );
4070 const SMDS_MeshNode* newNode = aMesh->AddNode( aX, aY, aZ );
4071 myLastCreatedNodes.Append(newNode);
4072 srcNodes.Append( node );
4073 listNewNodes.push_back( newNode );
4083 // if current elem is quadratic and current node is not medium
4084 // we have to check - may be it is needed to insert additional nodes
4085 if( elem->IsQuadratic() && !elem->IsMediumNode(node) ) {
4086 list< const SMDS_MeshNode* > & listNewNodes = nIt->second;
4087 if(listNewNodes.size()==aNbTP-1) {
4088 vector<const SMDS_MeshNode*> aNodes(2*(aNbTP-1));
4089 gp_XYZ P(node->X(), node->Y(), node->Z());
4090 list< const SMDS_MeshNode* >::iterator it = listNewNodes.begin();
4092 for(i=0; i<aNbTP-1; i++) {
4093 const SMDS_MeshNode* N = *it;
4094 double x = ( N->X() + P.X() )/2.;
4095 double y = ( N->Y() + P.Y() )/2.;
4096 double z = ( N->Z() + P.Z() )/2.;
4097 const SMDS_MeshNode* newN = aMesh->AddNode(x,y,z);
4098 srcNodes.Append( node );
4099 myLastCreatedNodes.Append(newN);
4102 P = gp_XYZ(N->X(),N->Y(),N->Z());
4104 listNewNodes.clear();
4105 for(i=0; i<2*(aNbTP-1); i++) {
4106 listNewNodes.push_back(aNodes[i]);
4112 newNodesItVec.push_back( nIt );
4114 // make new elements
4115 //sweepElement( aMesh, elem, newNodesItVec, newElemsMap[elem],
4116 // newNodesItVec[0]->second.size(), myLastCreatedElems );
4117 sweepElement( elem, newNodesItVec, newElemsMap[elem], aNbTP-1, srcElems );
4120 makeWalls( mapNewNodes, newElemsMap, mapElemNewNodes, theElements, aNbTP-1, srcElems );
4122 if ( theMakeGroups )
4123 generateGroups( srcNodes, srcElems, "extruded");
4128 //=======================================================================
4129 //function : Transform
4131 //=======================================================================
4133 SMESH_MeshEditor::PGroupIDs
4134 SMESH_MeshEditor::Transform (TIDSortedElemSet & theElems,
4135 const gp_Trsf& theTrsf,
4137 const bool theMakeGroups,
4138 SMESH_Mesh* theTargetMesh)
4140 myLastCreatedElems.Clear();
4141 myLastCreatedNodes.Clear();
4143 bool needReverse = false;
4144 string groupPostfix;
4145 switch ( theTrsf.Form() ) {
4150 groupPostfix = "mirrored";
4153 groupPostfix = "rotated";
4155 case gp_Translation:
4156 groupPostfix = "translated";
4159 groupPostfix = "scaled";
4162 needReverse = false;
4163 groupPostfix = "transformed";
4166 SMESH_MeshEditor targetMeshEditor( theTargetMesh );
4167 SMESHDS_Mesh* aTgtMesh = theTargetMesh ? theTargetMesh->GetMeshDS() : 0;
4168 SMESHDS_Mesh* aMesh = GetMeshDS();
4171 // map old node to new one
4172 TNodeNodeMap nodeMap;
4174 // elements sharing moved nodes; those of them which have all
4175 // nodes mirrored but are not in theElems are to be reversed
4176 TIDSortedElemSet inverseElemSet;
4178 // source elements for each generated one
4179 SMESH_SequenceOfElemPtr srcElems, srcNodes;
4182 TIDSortedElemSet::iterator itElem;
4183 for ( itElem = theElems.begin(); itElem != theElems.end(); itElem++ ) {
4184 const SMDS_MeshElement* elem = *itElem;
4188 // loop on elem nodes
4189 SMDS_ElemIteratorPtr itN = elem->nodesIterator();
4190 while ( itN->more() ) {
4192 // check if a node has been already transformed
4193 const SMDS_MeshNode* node = cast2Node( itN->next() );
4194 pair<TNodeNodeMap::iterator,bool> n2n_isnew =
4195 nodeMap.insert( make_pair ( node, node ));
4196 if ( !n2n_isnew.second )
4200 coord[0] = node->X();
4201 coord[1] = node->Y();
4202 coord[2] = node->Z();
4203 theTrsf.Transforms( coord[0], coord[1], coord[2] );
4204 if ( theTargetMesh ) {
4205 const SMDS_MeshNode * newNode = aTgtMesh->AddNode( coord[0], coord[1], coord[2] );
4206 n2n_isnew.first->second = newNode;
4207 myLastCreatedNodes.Append(newNode);
4208 srcNodes.Append( node );
4210 else if ( theCopy ) {
4211 const SMDS_MeshNode * newNode = aMesh->AddNode( coord[0], coord[1], coord[2] );
4212 n2n_isnew.first->second = newNode;
4213 myLastCreatedNodes.Append(newNode);
4214 srcNodes.Append( node );
4217 aMesh->MoveNode( node, coord[0], coord[1], coord[2] );
4218 // node position on shape becomes invalid
4219 const_cast< SMDS_MeshNode* > ( node )->SetPosition
4220 ( SMDS_SpacePosition::originSpacePosition() );
4223 // keep inverse elements
4224 if ( !theCopy && !theTargetMesh && needReverse ) {
4225 SMDS_ElemIteratorPtr invElemIt = node->GetInverseElementIterator();
4226 while ( invElemIt->more() ) {
4227 const SMDS_MeshElement* iel = invElemIt->next();
4228 inverseElemSet.insert( iel );
4234 // either create new elements or reverse mirrored ones
4235 if ( !theCopy && !needReverse && !theTargetMesh )
4238 TIDSortedElemSet::iterator invElemIt = inverseElemSet.begin();
4239 for ( ; invElemIt != inverseElemSet.end(); invElemIt++ )
4240 theElems.insert( *invElemIt );
4242 // replicate or reverse elements
4245 REV_TETRA = 0, // = nbNodes - 4
4246 REV_PYRAMID = 1, // = nbNodes - 4
4247 REV_PENTA = 2, // = nbNodes - 4
4249 REV_HEXA = 4, // = nbNodes - 4
4253 { 2, 1, 0, 3, 4, 0, 0, 0 }, // REV_TETRA
4254 { 2, 1, 0, 3, 4, 0, 0, 0 }, // REV_PYRAMID
4255 { 2, 1, 0, 5, 4, 3, 0, 0 }, // REV_PENTA
4256 { 2, 1, 0, 3, 0, 0, 0, 0 }, // REV_FACE
4257 { 2, 1, 0, 3, 6, 5, 4, 7 }, // REV_HEXA
4258 { 0, 1, 2, 3, 4, 5, 6, 7 } // FORWARD
4261 for ( itElem = theElems.begin(); itElem != theElems.end(); itElem++ )
4263 const SMDS_MeshElement* elem = *itElem;
4264 if ( !elem || elem->GetType() == SMDSAbs_Node )
4267 int nbNodes = elem->NbNodes();
4268 int elemType = elem->GetType();
4270 if (elem->IsPoly()) {
4271 // Polygon or Polyhedral Volume
4272 switch ( elemType ) {
4275 vector<const SMDS_MeshNode*> poly_nodes (nbNodes);
4277 SMDS_ElemIteratorPtr itN = elem->nodesIterator();
4278 while (itN->more()) {
4279 const SMDS_MeshNode* node =
4280 static_cast<const SMDS_MeshNode*>(itN->next());
4281 TNodeNodeMap::iterator nodeMapIt = nodeMap.find(node);
4282 if (nodeMapIt == nodeMap.end())
4283 break; // not all nodes transformed
4285 // reverse mirrored faces and volumes
4286 poly_nodes[nbNodes - iNode - 1] = (*nodeMapIt).second;
4288 poly_nodes[iNode] = (*nodeMapIt).second;
4292 if ( iNode != nbNodes )
4293 continue; // not all nodes transformed
4295 if ( theTargetMesh ) {
4296 myLastCreatedElems.Append(aTgtMesh->AddPolygonalFace(poly_nodes));
4297 srcElems.Append( elem );
4299 else if ( theCopy ) {
4300 myLastCreatedElems.Append(aMesh->AddPolygonalFace(poly_nodes));
4301 srcElems.Append( elem );
4304 aMesh->ChangePolygonNodes(elem, poly_nodes);
4308 case SMDSAbs_Volume:
4310 // ATTENTION: Reversing is not yet done!!!
4311 const SMDS_PolyhedralVolumeOfNodes* aPolyedre =
4312 dynamic_cast<const SMDS_PolyhedralVolumeOfNodes*>( elem );
4314 MESSAGE("Warning: bad volumic element");
4318 vector<const SMDS_MeshNode*> poly_nodes;
4319 vector<int> quantities;
4321 bool allTransformed = true;
4322 int nbFaces = aPolyedre->NbFaces();
4323 for (int iface = 1; iface <= nbFaces && allTransformed; iface++) {
4324 int nbFaceNodes = aPolyedre->NbFaceNodes(iface);
4325 for (int inode = 1; inode <= nbFaceNodes && allTransformed; inode++) {
4326 const SMDS_MeshNode* node = aPolyedre->GetFaceNode(iface, inode);
4327 TNodeNodeMap::iterator nodeMapIt = nodeMap.find(node);
4328 if (nodeMapIt == nodeMap.end()) {
4329 allTransformed = false; // not all nodes transformed
4331 poly_nodes.push_back((*nodeMapIt).second);
4334 quantities.push_back(nbFaceNodes);
4336 if ( !allTransformed )
4337 continue; // not all nodes transformed
4339 if ( theTargetMesh ) {
4340 myLastCreatedElems.Append(aTgtMesh->AddPolyhedralVolume(poly_nodes, quantities));
4341 srcElems.Append( elem );
4343 else if ( theCopy ) {
4344 myLastCreatedElems.Append(aMesh->AddPolyhedralVolume(poly_nodes, quantities));
4345 srcElems.Append( elem );
4348 aMesh->ChangePolyhedronNodes(elem, poly_nodes, quantities);
4358 int* i = index[ FORWARD ];
4359 if ( needReverse && nbNodes > 2) // reverse mirrored faces and volumes
4360 if ( elemType == SMDSAbs_Face )
4361 i = index[ REV_FACE ];
4363 i = index[ nbNodes - 4 ];
4365 if(elem->IsQuadratic()) {
4366 static int anIds[] = {0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19};
4369 if(nbNodes==3) { // quadratic edge
4370 static int anIds[] = {1,0,2};
4373 else if(nbNodes==6) { // quadratic triangle
4374 static int anIds[] = {0,2,1,5,4,3};
4377 else if(nbNodes==8) { // quadratic quadrangle
4378 static int anIds[] = {0,3,2,1,7,6,5,4};
4381 else if(nbNodes==10) { // quadratic tetrahedron of 10 nodes
4382 static int anIds[] = {0,2,1,3,6,5,4,7,9,8};
4385 else if(nbNodes==13) { // quadratic pyramid of 13 nodes
4386 static int anIds[] = {0,3,2,1,4,8,7,6,5,9,12,11,10};
4389 else if(nbNodes==15) { // quadratic pentahedron with 15 nodes
4390 static int anIds[] = {0,2,1,3,5,4,8,7,6,11,10,9,12,14,13};
4393 else { // nbNodes==20 - quadratic hexahedron with 20 nodes
4394 static int anIds[] = {0,3,2,1,4,7,6,5,11,10,9,8,15,14,13,12,16,19,18,17};
4400 // find transformed nodes
4401 vector<const SMDS_MeshNode*> nodes(nbNodes);
4403 SMDS_ElemIteratorPtr itN = elem->nodesIterator();
4404 while ( itN->more() ) {
4405 const SMDS_MeshNode* node =
4406 static_cast<const SMDS_MeshNode*>( itN->next() );
4407 TNodeNodeMap::iterator nodeMapIt = nodeMap.find( node );
4408 if ( nodeMapIt == nodeMap.end() )
4409 break; // not all nodes transformed
4410 nodes[ i [ iNode++ ]] = (*nodeMapIt).second;
4412 if ( iNode != nbNodes )
4413 continue; // not all nodes transformed
4415 if ( theTargetMesh ) {
4416 if ( SMDS_MeshElement* copy =
4417 targetMeshEditor.AddElement( nodes, elem->GetType(), elem->IsPoly() )) {
4418 myLastCreatedElems.Append( copy );
4419 srcElems.Append( elem );
4422 else if ( theCopy ) {
4423 if ( SMDS_MeshElement* copy = AddElement( nodes, elem->GetType(), elem->IsPoly() )) {
4424 myLastCreatedElems.Append( copy );
4425 srcElems.Append( elem );
4429 // reverse element as it was reversed by transformation
4431 aMesh->ChangeElementNodes( elem, &nodes[0], nbNodes );
4435 PGroupIDs newGroupIDs;
4437 if ( theMakeGroups && theCopy ||
4438 theMakeGroups && theTargetMesh )
4439 newGroupIDs = generateGroups( srcNodes, srcElems, groupPostfix, theTargetMesh );
4444 //=======================================================================
4446 * \brief Create groups of elements made during transformation
4447 * \param nodeGens - nodes making corresponding myLastCreatedNodes
4448 * \param elemGens - elements making corresponding myLastCreatedElems
4449 * \param postfix - to append to names of new groups
4451 //=======================================================================
4453 SMESH_MeshEditor::PGroupIDs
4454 SMESH_MeshEditor::generateGroups(const SMESH_SequenceOfElemPtr& nodeGens,
4455 const SMESH_SequenceOfElemPtr& elemGens,
4456 const std::string& postfix,
4457 SMESH_Mesh* targetMesh)
4459 PGroupIDs newGroupIDs( new list<int> );
4460 SMESH_Mesh* mesh = targetMesh ? targetMesh : GetMesh();
4462 // Sort existing groups by types and collect their names
4464 // to store an old group and a generated new one
4465 typedef pair< SMESHDS_GroupBase*, SMDS_MeshGroup* > TOldNewGroup;
4466 vector< list< TOldNewGroup > > groupsByType( SMDSAbs_NbElementTypes );
4468 set< string > groupNames;
4470 SMDS_MeshGroup* nullNewGroup = (SMDS_MeshGroup*) 0;
4471 SMESH_Mesh::GroupIteratorPtr groupIt = GetMesh()->GetGroups();
4472 while ( groupIt->more() ) {
4473 SMESH_Group * group = groupIt->next();
4474 if ( !group ) continue;
4475 SMESHDS_GroupBase* groupDS = group->GetGroupDS();
4476 if ( !groupDS || groupDS->IsEmpty() ) continue;
4477 groupNames.insert( group->GetName() );
4478 groupDS->SetStoreName( group->GetName() );
4479 groupsByType[ groupDS->GetType() ].push_back( make_pair( groupDS, nullNewGroup ));
4484 // loop on nodes and elements
4485 for ( int isNodes = 0; isNodes < 2; ++isNodes )
4487 const SMESH_SequenceOfElemPtr& gens = isNodes ? nodeGens : elemGens;
4488 const SMESH_SequenceOfElemPtr& elems = isNodes ? myLastCreatedNodes : myLastCreatedElems;
4489 if ( gens.Length() != elems.Length() )
4490 throw SALOME_Exception(LOCALIZED("invalid args"));
4492 // loop on created elements
4493 for (int iElem = 1; iElem <= elems.Length(); ++iElem )
4495 const SMDS_MeshElement* sourceElem = gens( iElem );
4496 if ( !sourceElem ) {
4497 MESSAGE("generateGroups(): NULL source element");
4500 list< TOldNewGroup > & groupsOldNew = groupsByType[ sourceElem->GetType() ];
4501 if ( groupsOldNew.empty() ) {
4502 while ( iElem < gens.Length() && gens( iElem+1 ) == sourceElem )
4503 ++iElem; // skip all elements made by sourceElem
4506 // collect all elements made by sourceElem
4507 list< const SMDS_MeshElement* > resultElems;
4508 if ( const SMDS_MeshElement* resElem = elems( iElem ))
4509 if ( resElem != sourceElem )
4510 resultElems.push_back( resElem );
4511 while ( iElem < gens.Length() && gens( iElem+1 ) == sourceElem )
4512 if ( const SMDS_MeshElement* resElem = elems( ++iElem ))
4513 if ( resElem != sourceElem )
4514 resultElems.push_back( resElem );
4515 // do not generate element groups from node ones
4516 if ( sourceElem->GetType() == SMDSAbs_Node &&
4517 elems( iElem )->GetType() != SMDSAbs_Node )
4520 // add resultElems to groups made by ones the sourceElem belongs to
4521 list< TOldNewGroup >::iterator gOldNew, gLast = groupsOldNew.end();
4522 for ( gOldNew = groupsOldNew.begin(); gOldNew != gLast; ++gOldNew )
4524 SMESHDS_GroupBase* oldGroup = gOldNew->first;
4525 if ( oldGroup->Contains( sourceElem )) // sourceElem in oldGroup
4527 SMDS_MeshGroup* & newGroup = gOldNew->second;
4528 if ( !newGroup )// create a new group
4531 string name = oldGroup->GetStoreName();
4532 if ( !targetMesh ) {
4536 while ( !groupNames.insert( name ).second ) // name exists
4542 TCollection_AsciiString nbStr(nb+1);
4543 name.resize( name.rfind('_')+1 );
4544 name += nbStr.ToCString();
4551 SMESH_Group* group = mesh->AddGroup( resultElems.back()->GetType(),
4553 SMESHDS_Group* groupDS = static_cast<SMESHDS_Group*>(group->GetGroupDS());
4554 newGroup = & groupDS->SMDSGroup();
4555 newGroupIDs->push_back( id );
4558 // fill in a new group
4559 list< const SMDS_MeshElement* >::iterator resLast = resultElems.end(), resElemIt;
4560 for ( resElemIt = resultElems.begin(); resElemIt != resLast; ++resElemIt )
4561 newGroup->Add( *resElemIt );
4564 } // loop on created elements
4565 }// loop on nodes and elements
4570 //=======================================================================
4571 //function : FindCoincidentNodes
4572 //purpose : Return list of group of nodes close to each other within theTolerance
4573 // Search among theNodes or in the whole mesh if theNodes is empty using
4574 // an Octree algorithm
4575 //=======================================================================
4577 void SMESH_MeshEditor::FindCoincidentNodes (set<const SMDS_MeshNode*> & theNodes,
4578 const double theTolerance,
4579 TListOfListOfNodes & theGroupsOfNodes)
4581 myLastCreatedElems.Clear();
4582 myLastCreatedNodes.Clear();
4584 set<const SMDS_MeshNode*> nodes;
4585 if ( theNodes.empty() )
4586 { // get all nodes in the mesh
4587 SMDS_NodeIteratorPtr nIt = GetMeshDS()->nodesIterator();
4588 while ( nIt->more() )
4589 nodes.insert( nodes.end(),nIt->next());
4593 SMESH_OctreeNode::FindCoincidentNodes ( nodes, &theGroupsOfNodes, theTolerance);
4597 //=======================================================================
4599 * \brief Implementation of search for the node closest to point
4601 //=======================================================================
4603 struct SMESH_NodeSearcherImpl: public SMESH_NodeSearcher
4606 * \brief Constructor
4608 SMESH_NodeSearcherImpl( const SMESHDS_Mesh* theMesh )
4610 set<const SMDS_MeshNode*> nodes;
4612 SMDS_NodeIteratorPtr nIt = theMesh->nodesIterator();
4613 while ( nIt->more() )
4614 nodes.insert( nodes.end(), nIt->next() );
4616 myOctreeNode = new SMESH_OctreeNode(nodes) ;
4619 * \brief Do it's job
4621 const SMDS_MeshNode* FindClosestTo( const gp_Pnt& thePnt )
4623 SMDS_MeshNode tgtNode( thePnt.X(), thePnt.Y(), thePnt.Z() );
4624 list<const SMDS_MeshNode*> nodes;
4625 const double precision = 1e-6;
4626 myOctreeNode->NodesAround( &tgtNode, &nodes, precision );
4628 double minSqDist = DBL_MAX;
4630 if ( nodes.empty() ) // get all nodes of OctreeNode's closest to thePnt
4632 // sort leafs by their distance from thePnt
4633 typedef map< double, SMESH_OctreeNode* > TDistTreeMap;
4634 TDistTreeMap treeMap;
4635 list< SMESH_OctreeNode* > treeList;
4636 list< SMESH_OctreeNode* >::iterator trIt;
4637 treeList.push_back( myOctreeNode );
4638 for ( trIt = treeList.begin(); trIt != treeList.end(); ++trIt)
4640 SMESH_OctreeNode* tree = *trIt;
4641 if ( !tree->isLeaf() ) { // put children to the queue
4642 SMESH_OctreeNodeIteratorPtr cIt = tree->GetChildrenIterator();
4643 while ( cIt->more() )
4644 treeList.push_back( cIt->next() );
4646 else if ( tree->NbNodes() ) { // put tree to treeMap
4647 tree->getBox( box );
4648 double sqDist = thePnt.SquareDistance( 0.5 * ( box.CornerMin() + box.CornerMax() ));
4649 pair<TDistTreeMap::iterator,bool> it_in = treeMap.insert( make_pair( sqDist, tree ));
4650 if ( !it_in.second ) // not unique distance to box center
4651 treeMap.insert( it_in.first, make_pair( sqDist - 1e-13*treeMap.size(), tree ));
4654 // find distance after which there is no sense to check tree's
4655 double sqLimit = DBL_MAX;
4656 TDistTreeMap::iterator sqDist_tree = treeMap.begin();
4657 if ( treeMap.size() > 5 ) {
4658 SMESH_OctreeNode* closestTree = sqDist_tree->second;
4659 closestTree->getBox( box );
4660 double limit = sqrt( sqDist_tree->first ) + sqrt ( box.SquareExtent() );
4661 sqLimit = limit * limit;
4663 // get all nodes from trees
4664 for ( ; sqDist_tree != treeMap.end(); ++sqDist_tree) {
4665 if ( sqDist_tree->first > sqLimit )
4667 SMESH_OctreeNode* tree = sqDist_tree->second;
4668 tree->NodesAround( tree->GetNodeIterator()->next(), &nodes );
4671 // find closest among nodes
4672 minSqDist = DBL_MAX;
4673 const SMDS_MeshNode* closestNode = 0;
4674 list<const SMDS_MeshNode*>::iterator nIt = nodes.begin();
4675 for ( ; nIt != nodes.end(); ++nIt ) {
4676 double sqDist = thePnt.SquareDistance( TNodeXYZ( *nIt ) );
4677 if ( minSqDist > sqDist ) {
4687 ~SMESH_NodeSearcherImpl() { delete myOctreeNode; }
4689 SMESH_OctreeNode* myOctreeNode;
4692 //=======================================================================
4694 * \brief Return SMESH_NodeSearcher
4696 //=======================================================================
4698 SMESH_NodeSearcher* SMESH_MeshEditor::GetNodeSearcher()
4700 return new SMESH_NodeSearcherImpl( GetMeshDS() );
4703 //=======================================================================
4704 //function : SimplifyFace
4706 //=======================================================================
4707 int SMESH_MeshEditor::SimplifyFace (const vector<const SMDS_MeshNode *> faceNodes,
4708 vector<const SMDS_MeshNode *>& poly_nodes,
4709 vector<int>& quantities) const
4711 int nbNodes = faceNodes.size();
4716 set<const SMDS_MeshNode*> nodeSet;
4718 // get simple seq of nodes
4719 //const SMDS_MeshNode* simpleNodes[ nbNodes ];
4720 vector<const SMDS_MeshNode*> simpleNodes( nbNodes );
4721 int iSimple = 0, nbUnique = 0;
4723 simpleNodes[iSimple++] = faceNodes[0];
4725 for (int iCur = 1; iCur < nbNodes; iCur++) {
4726 if (faceNodes[iCur] != simpleNodes[iSimple - 1]) {
4727 simpleNodes[iSimple++] = faceNodes[iCur];
4728 if (nodeSet.insert( faceNodes[iCur] ).second)
4732 int nbSimple = iSimple;
4733 if (simpleNodes[nbSimple - 1] == simpleNodes[0]) {
4743 bool foundLoop = (nbSimple > nbUnique);
4746 set<const SMDS_MeshNode*> loopSet;
4747 for (iSimple = 0; iSimple < nbSimple && !foundLoop; iSimple++) {
4748 const SMDS_MeshNode* n = simpleNodes[iSimple];
4749 if (!loopSet.insert( n ).second) {
4753 int iC = 0, curLast = iSimple;
4754 for (; iC < curLast; iC++) {
4755 if (simpleNodes[iC] == n) break;
4757 int loopLen = curLast - iC;
4759 // create sub-element
4761 quantities.push_back(loopLen);
4762 for (; iC < curLast; iC++) {
4763 poly_nodes.push_back(simpleNodes[iC]);
4766 // shift the rest nodes (place from the first loop position)
4767 for (iC = curLast + 1; iC < nbSimple; iC++) {
4768 simpleNodes[iC - loopLen] = simpleNodes[iC];
4770 nbSimple -= loopLen;
4773 } // for (iSimple = 0; iSimple < nbSimple; iSimple++)
4774 } // while (foundLoop)
4778 quantities.push_back(iSimple);
4779 for (int i = 0; i < iSimple; i++)
4780 poly_nodes.push_back(simpleNodes[i]);
4786 //=======================================================================
4787 //function : MergeNodes
4788 //purpose : In each group, the cdr of nodes are substituted by the first one
4790 //=======================================================================
4792 void SMESH_MeshEditor::MergeNodes (TListOfListOfNodes & theGroupsOfNodes)
4794 myLastCreatedElems.Clear();
4795 myLastCreatedNodes.Clear();
4797 SMESHDS_Mesh* aMesh = GetMeshDS();
4799 TNodeNodeMap nodeNodeMap; // node to replace - new node
4800 set<const SMDS_MeshElement*> elems; // all elements with changed nodes
4801 list< int > rmElemIds, rmNodeIds;
4803 // Fill nodeNodeMap and elems
4805 TListOfListOfNodes::iterator grIt = theGroupsOfNodes.begin();
4806 for ( ; grIt != theGroupsOfNodes.end(); grIt++ ) {
4807 list<const SMDS_MeshNode*>& nodes = *grIt;
4808 list<const SMDS_MeshNode*>::iterator nIt = nodes.begin();
4809 const SMDS_MeshNode* nToKeep = *nIt;
4810 for ( ++nIt; nIt != nodes.end(); nIt++ ) {
4811 const SMDS_MeshNode* nToRemove = *nIt;
4812 nodeNodeMap.insert( TNodeNodeMap::value_type( nToRemove, nToKeep ));
4813 if ( nToRemove != nToKeep ) {
4814 rmNodeIds.push_back( nToRemove->GetID() );
4815 AddToSameGroups( nToKeep, nToRemove, aMesh );
4818 SMDS_ElemIteratorPtr invElemIt = nToRemove->GetInverseElementIterator();
4819 while ( invElemIt->more() ) {
4820 const SMDS_MeshElement* elem = invElemIt->next();
4825 // Change element nodes or remove an element
4827 set<const SMDS_MeshElement*>::iterator eIt = elems.begin();
4828 for ( ; eIt != elems.end(); eIt++ ) {
4829 const SMDS_MeshElement* elem = *eIt;
4830 int nbNodes = elem->NbNodes();
4831 int aShapeId = FindShape( elem );
4833 set<const SMDS_MeshNode*> nodeSet;
4834 vector< const SMDS_MeshNode*> curNodes( nbNodes ), uniqueNodes( nbNodes );
4835 int iUnique = 0, iCur = 0, nbRepl = 0;
4836 vector<int> iRepl( nbNodes );
4838 // get new seq of nodes
4839 SMDS_ElemIteratorPtr itN = elem->nodesIterator();
4840 while ( itN->more() ) {
4841 const SMDS_MeshNode* n =
4842 static_cast<const SMDS_MeshNode*>( itN->next() );
4844 TNodeNodeMap::iterator nnIt = nodeNodeMap.find( n );
4845 if ( nnIt != nodeNodeMap.end() ) { // n sticks
4847 iRepl[ nbRepl++ ] = iCur;
4849 curNodes[ iCur ] = n;
4850 bool isUnique = nodeSet.insert( n ).second;
4852 uniqueNodes[ iUnique++ ] = n;
4856 // Analyse element topology after replacement
4859 int nbUniqueNodes = nodeSet.size();
4860 if ( nbNodes != nbUniqueNodes ) { // some nodes stick
4861 // Polygons and Polyhedral volumes
4862 if (elem->IsPoly()) {
4864 if (elem->GetType() == SMDSAbs_Face) {
4866 vector<const SMDS_MeshNode *> face_nodes (nbNodes);
4868 for (; inode < nbNodes; inode++) {
4869 face_nodes[inode] = curNodes[inode];
4872 vector<const SMDS_MeshNode *> polygons_nodes;
4873 vector<int> quantities;
4874 int nbNew = SimplifyFace(face_nodes, polygons_nodes, quantities);
4878 for (int iface = 0; iface < nbNew - 1; iface++) {
4879 int nbNodes = quantities[iface];
4880 vector<const SMDS_MeshNode *> poly_nodes (nbNodes);
4881 for (int ii = 0; ii < nbNodes; ii++, inode++) {
4882 poly_nodes[ii] = polygons_nodes[inode];
4884 SMDS_MeshElement* newElem = aMesh->AddPolygonalFace(poly_nodes);
4885 myLastCreatedElems.Append(newElem);
4887 aMesh->SetMeshElementOnShape(newElem, aShapeId);
4889 aMesh->ChangeElementNodes(elem, &polygons_nodes[inode], quantities[nbNew - 1]);
4892 rmElemIds.push_back(elem->GetID());
4896 else if (elem->GetType() == SMDSAbs_Volume) {
4897 // Polyhedral volume
4898 if (nbUniqueNodes < 4) {
4899 rmElemIds.push_back(elem->GetID());
4902 // each face has to be analized in order to check volume validity
4903 const SMDS_PolyhedralVolumeOfNodes* aPolyedre =
4904 static_cast<const SMDS_PolyhedralVolumeOfNodes*>( elem );
4906 int nbFaces = aPolyedre->NbFaces();
4908 vector<const SMDS_MeshNode *> poly_nodes;
4909 vector<int> quantities;
4911 for (int iface = 1; iface <= nbFaces; iface++) {
4912 int nbFaceNodes = aPolyedre->NbFaceNodes(iface);
4913 vector<const SMDS_MeshNode *> faceNodes (nbFaceNodes);
4915 for (int inode = 1; inode <= nbFaceNodes; inode++) {
4916 const SMDS_MeshNode * faceNode = aPolyedre->GetFaceNode(iface, inode);
4917 TNodeNodeMap::iterator nnIt = nodeNodeMap.find(faceNode);
4918 if (nnIt != nodeNodeMap.end()) { // faceNode sticks
4919 faceNode = (*nnIt).second;
4921 faceNodes[inode - 1] = faceNode;
4924 SimplifyFace(faceNodes, poly_nodes, quantities);
4927 if (quantities.size() > 3) {
4928 // to be done: remove coincident faces
4931 if (quantities.size() > 3)
4932 aMesh->ChangePolyhedronNodes(elem, poly_nodes, quantities);
4934 rmElemIds.push_back(elem->GetID());
4938 rmElemIds.push_back(elem->GetID());
4949 switch ( nbNodes ) {
4950 case 2: ///////////////////////////////////// EDGE
4951 isOk = false; break;
4952 case 3: ///////////////////////////////////// TRIANGLE
4953 isOk = false; break;
4955 if ( elem->GetType() == SMDSAbs_Volume ) // TETRAHEDRON
4957 else { //////////////////////////////////// QUADRANGLE
4958 if ( nbUniqueNodes < 3 )
4960 else if ( nbRepl == 2 && iRepl[ 1 ] - iRepl[ 0 ] == 2 )
4961 isOk = false; // opposite nodes stick
4964 case 6: ///////////////////////////////////// PENTAHEDRON
4965 if ( nbUniqueNodes == 4 ) {
4966 // ---------------------------------> tetrahedron
4968 iRepl[ 0 ] > 2 && iRepl[ 1 ] > 2 && iRepl[ 2 ] > 2 ) {
4969 // all top nodes stick: reverse a bottom
4970 uniqueNodes[ 0 ] = curNodes [ 1 ];
4971 uniqueNodes[ 1 ] = curNodes [ 0 ];
4973 else if (nbRepl == 3 &&
4974 iRepl[ 0 ] < 3 && iRepl[ 1 ] < 3 && iRepl[ 2 ] < 3 ) {
4975 // all bottom nodes stick: set a top before
4976 uniqueNodes[ 3 ] = uniqueNodes [ 0 ];
4977 uniqueNodes[ 0 ] = curNodes [ 3 ];
4978 uniqueNodes[ 1 ] = curNodes [ 4 ];
4979 uniqueNodes[ 2 ] = curNodes [ 5 ];
4981 else if (nbRepl == 4 &&
4982 iRepl[ 2 ] - iRepl [ 0 ] == 3 && iRepl[ 3 ] - iRepl [ 1 ] == 3 ) {
4983 // a lateral face turns into a line: reverse a bottom
4984 uniqueNodes[ 0 ] = curNodes [ 1 ];
4985 uniqueNodes[ 1 ] = curNodes [ 0 ];
4990 else if ( nbUniqueNodes == 5 ) {
4991 // PENTAHEDRON --------------------> 2 tetrahedrons
4992 if ( nbRepl == 2 && iRepl[ 1 ] - iRepl [ 0 ] == 3 ) {
4993 // a bottom node sticks with a linked top one
4995 SMDS_MeshElement* newElem =
4996 aMesh->AddVolume(curNodes[ 3 ],
4999 curNodes[ iRepl[ 0 ] == 2 ? 1 : 2 ]);
5000 myLastCreatedElems.Append(newElem);
5002 aMesh->SetMeshElementOnShape( newElem, aShapeId );
5003 // 2. : reverse a bottom
5004 uniqueNodes[ 0 ] = curNodes [ 1 ];
5005 uniqueNodes[ 1 ] = curNodes [ 0 ];
5015 if(elem->IsQuadratic()) { // Quadratic quadrangle
5028 if( iRepl[0]==0 && iRepl[1]==1 && iRepl[2]==4 ) {
5029 uniqueNodes[0] = curNodes[0];
5030 uniqueNodes[1] = curNodes[2];
5031 uniqueNodes[2] = curNodes[3];
5032 uniqueNodes[3] = curNodes[5];
5033 uniqueNodes[4] = curNodes[6];
5034 uniqueNodes[5] = curNodes[7];
5037 if( iRepl[0]==0 && iRepl[1]==3 && iRepl[2]==7 ) {
5038 uniqueNodes[0] = curNodes[0];
5039 uniqueNodes[1] = curNodes[1];
5040 uniqueNodes[2] = curNodes[2];
5041 uniqueNodes[3] = curNodes[4];
5042 uniqueNodes[4] = curNodes[5];
5043 uniqueNodes[5] = curNodes[6];
5046 if( iRepl[0]==0 && iRepl[1]==4 && iRepl[2]==7 ) {
5047 uniqueNodes[0] = curNodes[1];
5048 uniqueNodes[1] = curNodes[2];
5049 uniqueNodes[2] = curNodes[3];
5050 uniqueNodes[3] = curNodes[5];
5051 uniqueNodes[4] = curNodes[6];
5052 uniqueNodes[5] = curNodes[0];
5055 if( iRepl[0]==1 && iRepl[1]==2 && iRepl[2]==5 ) {
5056 uniqueNodes[0] = curNodes[0];
5057 uniqueNodes[1] = curNodes[1];
5058 uniqueNodes[2] = curNodes[3];
5059 uniqueNodes[3] = curNodes[4];
5060 uniqueNodes[4] = curNodes[6];
5061 uniqueNodes[5] = curNodes[7];
5064 if( iRepl[0]==1 && iRepl[1]==4 && iRepl[2]==5 ) {
5065 uniqueNodes[0] = curNodes[0];
5066 uniqueNodes[1] = curNodes[2];
5067 uniqueNodes[2] = curNodes[3];
5068 uniqueNodes[3] = curNodes[1];
5069 uniqueNodes[4] = curNodes[6];
5070 uniqueNodes[5] = curNodes[7];
5073 if( iRepl[0]==2 && iRepl[1]==3 && iRepl[2]==6 ) {
5074 uniqueNodes[0] = curNodes[0];
5075 uniqueNodes[1] = curNodes[1];
5076 uniqueNodes[2] = curNodes[2];
5077 uniqueNodes[3] = curNodes[4];
5078 uniqueNodes[4] = curNodes[5];
5079 uniqueNodes[5] = curNodes[7];
5082 if( iRepl[0]==2 && iRepl[1]==5 && iRepl[2]==6 ) {
5083 uniqueNodes[0] = curNodes[0];
5084 uniqueNodes[1] = curNodes[1];
5085 uniqueNodes[2] = curNodes[3];
5086 uniqueNodes[3] = curNodes[4];
5087 uniqueNodes[4] = curNodes[2];
5088 uniqueNodes[5] = curNodes[7];
5091 if( iRepl[0]==3 && iRepl[1]==6 && iRepl[2]==7 ) {
5092 uniqueNodes[0] = curNodes[0];
5093 uniqueNodes[1] = curNodes[1];
5094 uniqueNodes[2] = curNodes[2];
5095 uniqueNodes[3] = curNodes[4];
5096 uniqueNodes[4] = curNodes[5];
5097 uniqueNodes[5] = curNodes[3];
5103 //////////////////////////////////// HEXAHEDRON
5105 SMDS_VolumeTool hexa (elem);
5106 hexa.SetExternalNormal();
5107 if ( nbUniqueNodes == 4 && nbRepl == 6 ) {
5108 //////////////////////// ---> tetrahedron
5109 for ( int iFace = 0; iFace < 6; iFace++ ) {
5110 const int *ind = hexa.GetFaceNodesIndices( iFace ); // indices of face nodes
5111 if (curNodes[ind[ 0 ]] == curNodes[ind[ 1 ]] &&
5112 curNodes[ind[ 0 ]] == curNodes[ind[ 2 ]] &&
5113 curNodes[ind[ 0 ]] == curNodes[ind[ 3 ]] ) {
5114 // one face turns into a point ...
5115 int iOppFace = hexa.GetOppFaceIndex( iFace );
5116 ind = hexa.GetFaceNodesIndices( iOppFace );
5118 iUnique = 2; // reverse a tetrahedron bottom
5119 for ( iCur = 0; iCur < 4 && nbStick < 2; iCur++ ) {
5120 if ( curNodes[ind[ iCur ]] == curNodes[ind[ iCur + 1 ]] )
5122 else if ( iUnique >= 0 )
5123 uniqueNodes[ iUnique-- ] = curNodes[ind[ iCur ]];
5125 if ( nbStick == 1 ) {
5126 // ... and the opposite one - into a triangle.
5128 ind = hexa.GetFaceNodesIndices( iFace );
5129 uniqueNodes[ 3 ] = curNodes[ind[ 0 ]];
5136 else if (nbUniqueNodes == 5 && nbRepl == 4 ) {
5137 //////////////////// HEXAHEDRON ---> 2 tetrahedrons
5138 for ( int iFace = 0; iFace < 6; iFace++ ) {
5139 const int *ind = hexa.GetFaceNodesIndices( iFace ); // indices of face nodes
5140 if (curNodes[ind[ 0 ]] == curNodes[ind[ 1 ]] &&
5141 curNodes[ind[ 0 ]] == curNodes[ind[ 2 ]] &&
5142 curNodes[ind[ 0 ]] == curNodes[ind[ 3 ]] ) {
5143 // one face turns into a point ...
5144 int iOppFace = hexa.GetOppFaceIndex( iFace );
5145 ind = hexa.GetFaceNodesIndices( iOppFace );
5147 iUnique = 2; // reverse a tetrahedron 1 bottom
5148 for ( iCur = 0; iCur < 4 && nbStick == 0; iCur++ ) {
5149 if ( curNodes[ind[ iCur ]] == curNodes[ind[ iCur + 1 ]] )
5151 else if ( iUnique >= 0 )
5152 uniqueNodes[ iUnique-- ] = curNodes[ind[ iCur ]];
5154 if ( nbStick == 0 ) {
5155 // ... and the opposite one is a quadrangle
5157 const int* indTop = hexa.GetFaceNodesIndices( iFace );
5158 uniqueNodes[ 3 ] = curNodes[indTop[ 0 ]];
5161 SMDS_MeshElement* newElem =
5162 aMesh->AddVolume(curNodes[ind[ 0 ]],
5165 curNodes[indTop[ 0 ]]);
5166 myLastCreatedElems.Append(newElem);
5168 aMesh->SetMeshElementOnShape( newElem, aShapeId );
5175 else if ( nbUniqueNodes == 6 && nbRepl == 4 ) {
5176 ////////////////// HEXAHEDRON ---> 2 tetrahedrons or 1 prism
5177 // find indices of quad and tri faces
5178 int iQuadFace[ 6 ], iTriFace[ 6 ], nbQuad = 0, nbTri = 0, iFace;
5179 for ( iFace = 0; iFace < 6; iFace++ ) {
5180 const int *ind = hexa.GetFaceNodesIndices( iFace ); // indices of face nodes
5182 for ( iCur = 0; iCur < 4; iCur++ )
5183 nodeSet.insert( curNodes[ind[ iCur ]] );
5184 nbUniqueNodes = nodeSet.size();
5185 if ( nbUniqueNodes == 3 )
5186 iTriFace[ nbTri++ ] = iFace;
5187 else if ( nbUniqueNodes == 4 )
5188 iQuadFace[ nbQuad++ ] = iFace;
5190 if (nbQuad == 2 && nbTri == 4 &&
5191 hexa.GetOppFaceIndex( iQuadFace[ 0 ] ) == iQuadFace[ 1 ]) {
5192 // 2 opposite quadrangles stuck with a diagonal;
5193 // sample groups of merged indices: (0-4)(2-6)
5194 // --------------------------------------------> 2 tetrahedrons
5195 const int *ind1 = hexa.GetFaceNodesIndices( iQuadFace[ 0 ]); // indices of quad1 nodes
5196 const int *ind2 = hexa.GetFaceNodesIndices( iQuadFace[ 1 ]);
5197 int i0, i1d, i2, i3d, i0t, i2t; // d-daigonal, t-top
5198 if (curNodes[ind1[ 0 ]] == curNodes[ind2[ 0 ]] &&
5199 curNodes[ind1[ 2 ]] == curNodes[ind2[ 2 ]]) {
5200 // stuck with 0-2 diagonal
5208 else if (curNodes[ind1[ 1 ]] == curNodes[ind2[ 3 ]] &&
5209 curNodes[ind1[ 3 ]] == curNodes[ind2[ 1 ]]) {
5210 // stuck with 1-3 diagonal
5222 uniqueNodes[ 0 ] = curNodes [ i0 ];
5223 uniqueNodes[ 1 ] = curNodes [ i1d ];
5224 uniqueNodes[ 2 ] = curNodes [ i3d ];
5225 uniqueNodes[ 3 ] = curNodes [ i0t ];
5228 SMDS_MeshElement* newElem = aMesh->AddVolume(curNodes[ i1d ],
5232 myLastCreatedElems.Append(newElem);
5234 aMesh->SetMeshElementOnShape( newElem, aShapeId );
5237 else if (( nbTri == 2 && nbQuad == 3 ) || // merged (0-4)(1-5)
5238 ( nbTri == 4 && nbQuad == 2 )) { // merged (7-4)(1-5)
5239 // --------------------------------------------> prism
5240 // find 2 opposite triangles
5242 for ( iFace = 0; iFace + 1 < nbTri; iFace++ ) {
5243 if ( hexa.GetOppFaceIndex( iTriFace[ iFace ] ) == iTriFace[ iFace + 1 ]) {
5244 // find indices of kept and replaced nodes
5245 // and fill unique nodes of 2 opposite triangles
5246 const int *ind1 = hexa.GetFaceNodesIndices( iTriFace[ iFace ]);
5247 const int *ind2 = hexa.GetFaceNodesIndices( iTriFace[ iFace + 1 ]);
5248 const SMDS_MeshNode** hexanodes = hexa.GetNodes();
5249 // fill unique nodes
5252 for ( iCur = 0; iCur < 4 && isOk; iCur++ ) {
5253 const SMDS_MeshNode* n = curNodes[ind1[ iCur ]];
5254 const SMDS_MeshNode* nInit = hexanodes[ind1[ iCur ]];
5256 // iCur of a linked node of the opposite face (make normals co-directed):
5257 int iCurOpp = ( iCur == 1 || iCur == 3 ) ? 4 - iCur : iCur;
5258 // check that correspondent corners of triangles are linked
5259 if ( !hexa.IsLinked( ind1[ iCur ], ind2[ iCurOpp ] ))
5262 uniqueNodes[ iUnique ] = n;
5263 uniqueNodes[ iUnique + 3 ] = curNodes[ind2[ iCurOpp ]];
5272 } // if ( nbUniqueNodes == 6 && nbRepl == 4 )
5278 } // switch ( nbNodes )
5280 } // if ( nbNodes != nbUniqueNodes ) // some nodes stick
5283 if (elem->IsPoly() && elem->GetType() == SMDSAbs_Volume) {
5284 // Change nodes of polyedre
5285 const SMDS_PolyhedralVolumeOfNodes* aPolyedre =
5286 static_cast<const SMDS_PolyhedralVolumeOfNodes*>( elem );
5288 int nbFaces = aPolyedre->NbFaces();
5290 vector<const SMDS_MeshNode *> poly_nodes;
5291 vector<int> quantities (nbFaces);
5293 for (int iface = 1; iface <= nbFaces; iface++) {
5294 int inode, nbFaceNodes = aPolyedre->NbFaceNodes(iface);
5295 quantities[iface - 1] = nbFaceNodes;
5297 for (inode = 1; inode <= nbFaceNodes; inode++) {
5298 const SMDS_MeshNode* curNode = aPolyedre->GetFaceNode(iface, inode);
5300 TNodeNodeMap::iterator nnIt = nodeNodeMap.find( curNode );
5301 if (nnIt != nodeNodeMap.end()) { // curNode sticks
5302 curNode = (*nnIt).second;
5304 poly_nodes.push_back(curNode);
5307 aMesh->ChangePolyhedronNodes( elem, poly_nodes, quantities );
5311 // Change regular element or polygon
5312 aMesh->ChangeElementNodes( elem, & uniqueNodes[0], nbUniqueNodes );
5316 // Remove invalid regular element or invalid polygon
5317 rmElemIds.push_back( elem->GetID() );
5320 } // loop on elements
5322 // Remove equal nodes and bad elements
5324 Remove( rmNodeIds, true );
5325 Remove( rmElemIds, false );
5330 // ========================================================
5331 // class : SortableElement
5332 // purpose : allow sorting elements basing on their nodes
5333 // ========================================================
5334 class SortableElement : public set <const SMDS_MeshElement*>
5338 SortableElement( const SMDS_MeshElement* theElem )
5341 SMDS_ElemIteratorPtr nodeIt = theElem->nodesIterator();
5342 while ( nodeIt->more() )
5343 this->insert( nodeIt->next() );
5346 const SMDS_MeshElement* Get() const
5349 void Set(const SMDS_MeshElement* e) const
5354 mutable const SMDS_MeshElement* myElem;
5357 //=======================================================================
5358 //function : FindEqualElements
5359 //purpose : Return list of group of elements built on the same nodes.
5360 // Search among theElements or in the whole mesh if theElements is empty
5361 //=======================================================================
5362 void SMESH_MeshEditor::FindEqualElements(set<const SMDS_MeshElement*> & theElements,
5363 TListOfListOfElementsID & theGroupsOfElementsID)
5365 myLastCreatedElems.Clear();
5366 myLastCreatedNodes.Clear();
5368 typedef set<const SMDS_MeshElement*> TElemsSet;
5369 typedef map< SortableElement, int > TMapOfNodeSet;
5370 typedef list<int> TGroupOfElems;
5373 if ( theElements.empty() )
5374 { // get all elements in the mesh
5375 SMDS_ElemIteratorPtr eIt = GetMeshDS()->elementsIterator();
5376 while ( eIt->more() )
5377 elems.insert( elems.end(), eIt->next());
5380 elems = theElements;
5382 vector< TGroupOfElems > arrayOfGroups;
5383 TGroupOfElems groupOfElems;
5384 TMapOfNodeSet mapOfNodeSet;
5386 TElemsSet::iterator elemIt = elems.begin();
5387 for ( int i = 0, j=0; elemIt != elems.end(); ++elemIt, ++j ) {
5388 const SMDS_MeshElement* curElem = *elemIt;
5389 SortableElement SE(curElem);
5392 pair< TMapOfNodeSet::iterator, bool> pp = mapOfNodeSet.insert(make_pair(SE, i));
5393 if( !(pp.second) ) {
5394 TMapOfNodeSet::iterator& itSE = pp.first;
5395 ind = (*itSE).second;
5396 arrayOfGroups[ind].push_back(curElem->GetID());
5399 groupOfElems.clear();
5400 groupOfElems.push_back(curElem->GetID());
5401 arrayOfGroups.push_back(groupOfElems);
5406 vector< TGroupOfElems >::iterator groupIt = arrayOfGroups.begin();
5407 for ( ; groupIt != arrayOfGroups.end(); ++groupIt ) {
5408 groupOfElems = *groupIt;
5409 if ( groupOfElems.size() > 1 ) {
5410 groupOfElems.sort();
5411 theGroupsOfElementsID.push_back(groupOfElems);
5416 //=======================================================================
5417 //function : MergeElements
5418 //purpose : In each given group, substitute all elements by the first one.
5419 //=======================================================================
5421 void SMESH_MeshEditor::MergeElements(TListOfListOfElementsID & theGroupsOfElementsID)
5423 myLastCreatedElems.Clear();
5424 myLastCreatedNodes.Clear();
5426 typedef list<int> TListOfIDs;
5427 TListOfIDs rmElemIds; // IDs of elems to remove
5429 SMESHDS_Mesh* aMesh = GetMeshDS();
5431 TListOfListOfElementsID::iterator groupsIt = theGroupsOfElementsID.begin();
5432 while ( groupsIt != theGroupsOfElementsID.end() ) {
5433 TListOfIDs& aGroupOfElemID = *groupsIt;
5434 aGroupOfElemID.sort();
5435 int elemIDToKeep = aGroupOfElemID.front();
5436 const SMDS_MeshElement* elemToKeep = aMesh->FindElement(elemIDToKeep);
5437 aGroupOfElemID.pop_front();
5438 TListOfIDs::iterator idIt = aGroupOfElemID.begin();
5439 while ( idIt != aGroupOfElemID.end() ) {
5440 int elemIDToRemove = *idIt;
5441 const SMDS_MeshElement* elemToRemove = aMesh->FindElement(elemIDToRemove);
5442 // add the kept element in groups of removed one (PAL15188)
5443 AddToSameGroups( elemToKeep, elemToRemove, aMesh );
5444 rmElemIds.push_back( elemIDToRemove );
5450 Remove( rmElemIds, false );
5453 //=======================================================================
5454 //function : MergeEqualElements
5455 //purpose : Remove all but one of elements built on the same nodes.
5456 //=======================================================================
5458 void SMESH_MeshEditor::MergeEqualElements()
5460 set<const SMDS_MeshElement*> aMeshElements; /* empty input -
5461 to merge equal elements in the whole mesh */
5462 TListOfListOfElementsID aGroupsOfElementsID;
5463 FindEqualElements(aMeshElements, aGroupsOfElementsID);
5464 MergeElements(aGroupsOfElementsID);
5467 //=======================================================================
5468 //function : FindFaceInSet
5469 //purpose : Return a face having linked nodes n1 and n2 and which is
5470 // - not in avoidSet,
5471 // - in elemSet provided that !elemSet.empty()
5472 //=======================================================================
5474 const SMDS_MeshElement*
5475 SMESH_MeshEditor::FindFaceInSet(const SMDS_MeshNode* n1,
5476 const SMDS_MeshNode* n2,
5477 const TIDSortedElemSet& elemSet,
5478 const TIDSortedElemSet& avoidSet)
5481 SMDS_ElemIteratorPtr invElemIt = n1->GetInverseElementIterator(SMDSAbs_Face);
5482 while ( invElemIt->more() ) { // loop on inverse elements of n1
5483 const SMDS_MeshElement* elem = invElemIt->next();
5484 if (avoidSet.find( elem ) != avoidSet.end() )
5486 if ( !elemSet.empty() && elemSet.find( elem ) == elemSet.end())
5488 // get face nodes and find index of n1
5489 int i1, nbN = elem->NbNodes(), iNode = 0;
5490 //const SMDS_MeshNode* faceNodes[ nbN ], *n;
5491 vector<const SMDS_MeshNode*> faceNodes( nbN );
5492 const SMDS_MeshNode* n;
5493 SMDS_ElemIteratorPtr nIt = elem->nodesIterator();
5494 while ( nIt->more() ) {
5495 faceNodes[ iNode ] = static_cast<const SMDS_MeshNode*>( nIt->next() );
5496 if ( faceNodes[ iNode++ ] == n1 )
5499 // find a n2 linked to n1
5500 if(!elem->IsQuadratic()) {
5501 for ( iNode = 0; iNode < 2; iNode++ ) {
5502 if ( iNode ) // node before n1
5503 n = faceNodes[ i1 == 0 ? nbN - 1 : i1 - 1 ];
5504 else // node after n1
5505 n = faceNodes[ i1 + 1 == nbN ? 0 : i1 + 1 ];
5510 else { // analysis for quadratic elements
5511 bool IsFind = false;
5512 // check using only corner nodes
5513 for ( iNode = 0; iNode < 2; iNode++ ) {
5514 if ( iNode ) // node before n1
5515 n = faceNodes[ i1 == 0 ? nbN/2 - 1 : i1 - 1 ];
5516 else // node after n1
5517 n = faceNodes[ i1 + 1 == nbN/2 ? 0 : i1 + 1 ];
5525 // check using all nodes
5526 const SMDS_QuadraticFaceOfNodes* F =
5527 static_cast<const SMDS_QuadraticFaceOfNodes*>(elem);
5528 // use special nodes iterator
5530 SMDS_NodeIteratorPtr anIter = F->interlacedNodesIterator();
5531 while ( anIter->more() ) {
5532 faceNodes[iNode] = static_cast<const SMDS_MeshNode*>(anIter->next());
5533 if ( faceNodes[ iNode++ ] == n1 )
5536 for ( iNode = 0; iNode < 2; iNode++ ) {
5537 if ( iNode ) // node before n1
5538 n = faceNodes[ i1 == 0 ? nbN - 1 : i1 - 1 ];
5539 else // node after n1
5540 n = faceNodes[ i1 + 1 == nbN ? 0 : i1 + 1 ];
5546 } // end analysis for quadratic elements
5551 //=======================================================================
5552 //function : findAdjacentFace
5554 //=======================================================================
5556 static const SMDS_MeshElement* findAdjacentFace(const SMDS_MeshNode* n1,
5557 const SMDS_MeshNode* n2,
5558 const SMDS_MeshElement* elem)
5560 TIDSortedElemSet elemSet, avoidSet;
5562 avoidSet.insert ( elem );
5563 return SMESH_MeshEditor::FindFaceInSet( n1, n2, elemSet, avoidSet );
5566 //=======================================================================
5567 //function : FindFreeBorder
5569 //=======================================================================
5571 #define ControlFreeBorder SMESH::Controls::FreeEdges::IsFreeEdge
5573 bool SMESH_MeshEditor::FindFreeBorder (const SMDS_MeshNode* theFirstNode,
5574 const SMDS_MeshNode* theSecondNode,
5575 const SMDS_MeshNode* theLastNode,
5576 list< const SMDS_MeshNode* > & theNodes,
5577 list< const SMDS_MeshElement* >& theFaces)
5579 if ( !theFirstNode || !theSecondNode )
5581 // find border face between theFirstNode and theSecondNode
5582 const SMDS_MeshElement* curElem = findAdjacentFace( theFirstNode, theSecondNode, 0 );
5586 theFaces.push_back( curElem );
5587 theNodes.push_back( theFirstNode );
5588 theNodes.push_back( theSecondNode );
5590 //vector<const SMDS_MeshNode*> nodes;
5591 const SMDS_MeshNode *nIgnore = theFirstNode, *nStart = theSecondNode;
5592 set < const SMDS_MeshElement* > foundElems;
5593 bool needTheLast = ( theLastNode != 0 );
5595 while ( nStart != theLastNode ) {
5596 if ( nStart == theFirstNode )
5597 return !needTheLast;
5599 // find all free border faces sharing form nStart
5601 list< const SMDS_MeshElement* > curElemList;
5602 list< const SMDS_MeshNode* > nStartList;
5603 SMDS_ElemIteratorPtr invElemIt = nStart->GetInverseElementIterator(SMDSAbs_Face);
5604 while ( invElemIt->more() ) {
5605 const SMDS_MeshElement* e = invElemIt->next();
5606 if ( e == curElem || foundElems.insert( e ).second ) {
5608 int iNode = 0, nbNodes = e->NbNodes();
5609 //const SMDS_MeshNode* nodes[nbNodes+1];
5610 vector<const SMDS_MeshNode*> nodes(nbNodes+1);
5612 if(e->IsQuadratic()) {
5613 const SMDS_QuadraticFaceOfNodes* F =
5614 static_cast<const SMDS_QuadraticFaceOfNodes*>(e);
5615 // use special nodes iterator
5616 SMDS_NodeIteratorPtr anIter = F->interlacedNodesIterator();
5617 while( anIter->more() ) {
5618 nodes[ iNode++ ] = anIter->next();
5622 SMDS_ElemIteratorPtr nIt = e->nodesIterator();
5623 while ( nIt->more() )
5624 nodes[ iNode++ ] = static_cast<const SMDS_MeshNode*>( nIt->next() );
5626 nodes[ iNode ] = nodes[ 0 ];
5628 for ( iNode = 0; iNode < nbNodes; iNode++ )
5629 if (((nodes[ iNode ] == nStart && nodes[ iNode + 1] != nIgnore ) ||
5630 (nodes[ iNode + 1] == nStart && nodes[ iNode ] != nIgnore )) &&
5631 ControlFreeBorder( &nodes[ iNode ], e->GetID() ))
5633 nStartList.push_back( nodes[ iNode + ( nodes[ iNode ] == nStart ? 1 : 0 )]);
5634 curElemList.push_back( e );
5638 // analyse the found
5640 int nbNewBorders = curElemList.size();
5641 if ( nbNewBorders == 0 ) {
5642 // no free border furthermore
5643 return !needTheLast;
5645 else if ( nbNewBorders == 1 ) {
5646 // one more element found
5648 nStart = nStartList.front();
5649 curElem = curElemList.front();
5650 theFaces.push_back( curElem );
5651 theNodes.push_back( nStart );
5654 // several continuations found
5655 list< const SMDS_MeshElement* >::iterator curElemIt;
5656 list< const SMDS_MeshNode* >::iterator nStartIt;
5657 // check if one of them reached the last node
5658 if ( needTheLast ) {
5659 for (curElemIt = curElemList.begin(), nStartIt = nStartList.begin();
5660 curElemIt!= curElemList.end();
5661 curElemIt++, nStartIt++ )
5662 if ( *nStartIt == theLastNode ) {
5663 theFaces.push_back( *curElemIt );
5664 theNodes.push_back( *nStartIt );
5668 // find the best free border by the continuations
5669 list<const SMDS_MeshNode*> contNodes[ 2 ], *cNL;
5670 list<const SMDS_MeshElement*> contFaces[ 2 ], *cFL;
5671 for (curElemIt = curElemList.begin(), nStartIt = nStartList.begin();
5672 curElemIt!= curElemList.end();
5673 curElemIt++, nStartIt++ )
5675 cNL = & contNodes[ contNodes[0].empty() ? 0 : 1 ];
5676 cFL = & contFaces[ contFaces[0].empty() ? 0 : 1 ];
5677 // find one more free border
5678 if ( ! FindFreeBorder( nStart, *nStartIt, theLastNode, *cNL, *cFL )) {
5682 else if ( !contNodes[0].empty() && !contNodes[1].empty() ) {
5683 // choice: clear a worse one
5684 int iLongest = ( contNodes[0].size() < contNodes[1].size() ? 1 : 0 );
5685 int iWorse = ( needTheLast ? 1 - iLongest : iLongest );
5686 contNodes[ iWorse ].clear();
5687 contFaces[ iWorse ].clear();
5690 if ( contNodes[0].empty() && contNodes[1].empty() )
5693 // append the best free border
5694 cNL = & contNodes[ contNodes[0].empty() ? 1 : 0 ];
5695 cFL = & contFaces[ contFaces[0].empty() ? 1 : 0 ];
5696 theNodes.pop_back(); // remove nIgnore
5697 theNodes.pop_back(); // remove nStart
5698 theFaces.pop_back(); // remove curElem
5699 list< const SMDS_MeshNode* >::iterator nIt = cNL->begin();
5700 list< const SMDS_MeshElement* >::iterator fIt = cFL->begin();
5701 for ( ; nIt != cNL->end(); nIt++ ) theNodes.push_back( *nIt );
5702 for ( ; fIt != cFL->end(); fIt++ ) theFaces.push_back( *fIt );
5705 } // several continuations found
5706 } // while ( nStart != theLastNode )
5711 //=======================================================================
5712 //function : CheckFreeBorderNodes
5713 //purpose : Return true if the tree nodes are on a free border
5714 //=======================================================================
5716 bool SMESH_MeshEditor::CheckFreeBorderNodes(const SMDS_MeshNode* theNode1,
5717 const SMDS_MeshNode* theNode2,
5718 const SMDS_MeshNode* theNode3)
5720 list< const SMDS_MeshNode* > nodes;
5721 list< const SMDS_MeshElement* > faces;
5722 return FindFreeBorder( theNode1, theNode2, theNode3, nodes, faces);
5725 //=======================================================================
5726 //function : SewFreeBorder
5728 //=======================================================================
5730 SMESH_MeshEditor::Sew_Error
5731 SMESH_MeshEditor::SewFreeBorder (const SMDS_MeshNode* theBordFirstNode,
5732 const SMDS_MeshNode* theBordSecondNode,
5733 const SMDS_MeshNode* theBordLastNode,
5734 const SMDS_MeshNode* theSideFirstNode,
5735 const SMDS_MeshNode* theSideSecondNode,
5736 const SMDS_MeshNode* theSideThirdNode,
5737 const bool theSideIsFreeBorder,
5738 const bool toCreatePolygons,
5739 const bool toCreatePolyedrs)
5741 myLastCreatedElems.Clear();
5742 myLastCreatedNodes.Clear();
5744 MESSAGE("::SewFreeBorder()");
5745 Sew_Error aResult = SEW_OK;
5747 // ====================================
5748 // find side nodes and elements
5749 // ====================================
5751 list< const SMDS_MeshNode* > nSide[ 2 ];
5752 list< const SMDS_MeshElement* > eSide[ 2 ];
5753 list< const SMDS_MeshNode* >::iterator nIt[ 2 ];
5754 list< const SMDS_MeshElement* >::iterator eIt[ 2 ];
5758 if (!FindFreeBorder(theBordFirstNode,theBordSecondNode,theBordLastNode,
5759 nSide[0], eSide[0])) {
5760 MESSAGE(" Free Border 1 not found " );
5761 aResult = SEW_BORDER1_NOT_FOUND;
5763 if (theSideIsFreeBorder) {
5766 if (!FindFreeBorder(theSideFirstNode, theSideSecondNode, theSideThirdNode,
5767 nSide[1], eSide[1])) {
5768 MESSAGE(" Free Border 2 not found " );
5769 aResult = ( aResult != SEW_OK ? SEW_BOTH_BORDERS_NOT_FOUND : SEW_BORDER2_NOT_FOUND );
5772 if ( aResult != SEW_OK )
5775 if (!theSideIsFreeBorder) {
5779 // -------------------------------------------------------------------------
5781 // 1. If nodes to merge are not coincident, move nodes of the free border
5782 // from the coord sys defined by the direction from the first to last
5783 // nodes of the border to the correspondent sys of the side 2
5784 // 2. On the side 2, find the links most co-directed with the correspondent
5785 // links of the free border
5786 // -------------------------------------------------------------------------
5788 // 1. Since sewing may brake if there are volumes to split on the side 2,
5789 // we wont move nodes but just compute new coordinates for them
5790 typedef map<const SMDS_MeshNode*, gp_XYZ> TNodeXYZMap;
5791 TNodeXYZMap nBordXYZ;
5792 list< const SMDS_MeshNode* >& bordNodes = nSide[ 0 ];
5793 list< const SMDS_MeshNode* >::iterator nBordIt;
5795 gp_XYZ Pb1( theBordFirstNode->X(), theBordFirstNode->Y(), theBordFirstNode->Z() );
5796 gp_XYZ Pb2( theBordLastNode->X(), theBordLastNode->Y(), theBordLastNode->Z() );
5797 gp_XYZ Ps1( theSideFirstNode->X(), theSideFirstNode->Y(), theSideFirstNode->Z() );
5798 gp_XYZ Ps2( theSideSecondNode->X(), theSideSecondNode->Y(), theSideSecondNode->Z() );
5799 double tol2 = 1.e-8;
5800 gp_Vec Vbs1( Pb1 - Ps1 ),Vbs2( Pb2 - Ps2 );
5801 if ( Vbs1.SquareMagnitude() > tol2 || Vbs2.SquareMagnitude() > tol2 ) {
5802 // Need node movement.
5804 // find X and Z axes to create trsf
5805 gp_Vec Zb( Pb1 - Pb2 ), Zs( Ps1 - Ps2 );
5807 if ( X.SquareMagnitude() <= gp::Resolution() * gp::Resolution() )
5809 X = gp_Ax2( gp::Origin(), Zb ).XDirection();
5812 gp_Ax3 toBordAx( Pb1, Zb, X );
5813 gp_Ax3 fromSideAx( Ps1, Zs, X );
5814 gp_Ax3 toGlobalAx( gp::Origin(), gp::DZ(), gp::DX() );
5816 gp_Trsf toBordSys, fromSide2Sys;
5817 toBordSys.SetTransformation( toBordAx );
5818 fromSide2Sys.SetTransformation( fromSideAx, toGlobalAx );
5819 fromSide2Sys.SetScaleFactor( Zs.Magnitude() / Zb.Magnitude() );
5822 for ( nBordIt = bordNodes.begin(); nBordIt != bordNodes.end(); nBordIt++ ) {
5823 const SMDS_MeshNode* n = *nBordIt;
5824 gp_XYZ xyz( n->X(),n->Y(),n->Z() );
5825 toBordSys.Transforms( xyz );
5826 fromSide2Sys.Transforms( xyz );
5827 nBordXYZ.insert( TNodeXYZMap::value_type( n, xyz ));
5831 // just insert nodes XYZ in the nBordXYZ map
5832 for ( nBordIt = bordNodes.begin(); nBordIt != bordNodes.end(); nBordIt++ ) {
5833 const SMDS_MeshNode* n = *nBordIt;
5834 nBordXYZ.insert( TNodeXYZMap::value_type( n, gp_XYZ( n->X(),n->Y(),n->Z() )));
5838 // 2. On the side 2, find the links most co-directed with the correspondent
5839 // links of the free border
5841 list< const SMDS_MeshElement* >& sideElems = eSide[ 1 ];
5842 list< const SMDS_MeshNode* >& sideNodes = nSide[ 1 ];
5843 sideNodes.push_back( theSideFirstNode );
5845 bool hasVolumes = false;
5846 LinkID_Gen aLinkID_Gen( GetMeshDS() );
5847 set<long> foundSideLinkIDs, checkedLinkIDs;
5848 SMDS_VolumeTool volume;
5849 //const SMDS_MeshNode* faceNodes[ 4 ];
5851 const SMDS_MeshNode* sideNode;
5852 const SMDS_MeshElement* sideElem;
5853 const SMDS_MeshNode* prevSideNode = theSideFirstNode;
5854 const SMDS_MeshNode* prevBordNode = theBordFirstNode;
5855 nBordIt = bordNodes.begin();
5857 // border node position and border link direction to compare with
5858 gp_XYZ bordPos = nBordXYZ[ *nBordIt ];
5859 gp_XYZ bordDir = bordPos - nBordXYZ[ prevBordNode ];
5860 // choose next side node by link direction or by closeness to
5861 // the current border node:
5862 bool searchByDir = ( *nBordIt != theBordLastNode );
5864 // find the next node on the Side 2
5866 double maxDot = -DBL_MAX, minDist = DBL_MAX;
5868 checkedLinkIDs.clear();
5869 gp_XYZ prevXYZ( prevSideNode->X(), prevSideNode->Y(), prevSideNode->Z() );
5871 // loop on inverse elements of current node (prevSideNode) on the Side 2
5872 SMDS_ElemIteratorPtr invElemIt = prevSideNode->GetInverseElementIterator();
5873 while ( invElemIt->more() )
5875 const SMDS_MeshElement* elem = invElemIt->next();
5876 // prepare data for a loop on links coming to prevSideNode, of a face or a volume
5877 int iPrevNode, iNode = 0, nbNodes = elem->NbNodes();
5878 vector< const SMDS_MeshNode* > faceNodes( nbNodes, (const SMDS_MeshNode*)0 );
5879 bool isVolume = volume.Set( elem );
5880 const SMDS_MeshNode** nodes = isVolume ? volume.GetNodes() : & faceNodes[0];
5881 if ( isVolume ) // --volume
5883 else if ( elem->GetType()==SMDSAbs_Face ) { // --face
5884 // retrieve all face nodes and find iPrevNode - an index of the prevSideNode
5885 if(elem->IsQuadratic()) {
5886 const SMDS_QuadraticFaceOfNodes* F =
5887 static_cast<const SMDS_QuadraticFaceOfNodes*>(elem);
5888 // use special nodes iterator
5889 SMDS_NodeIteratorPtr anIter = F->interlacedNodesIterator();
5890 while( anIter->more() ) {
5891 nodes[ iNode ] = anIter->next();
5892 if ( nodes[ iNode++ ] == prevSideNode )
5893 iPrevNode = iNode - 1;
5897 SMDS_ElemIteratorPtr nIt = elem->nodesIterator();
5898 while ( nIt->more() ) {
5899 nodes[ iNode ] = cast2Node( nIt->next() );
5900 if ( nodes[ iNode++ ] == prevSideNode )
5901 iPrevNode = iNode - 1;
5904 // there are 2 links to check
5909 // loop on links, to be precise, on the second node of links
5910 for ( iNode = 0; iNode < nbNodes; iNode++ ) {
5911 const SMDS_MeshNode* n = nodes[ iNode ];
5913 if ( !volume.IsLinked( n, prevSideNode ))
5917 if ( iNode ) // a node before prevSideNode
5918 n = nodes[ iPrevNode == 0 ? elem->NbNodes() - 1 : iPrevNode - 1 ];
5919 else // a node after prevSideNode
5920 n = nodes[ iPrevNode + 1 == elem->NbNodes() ? 0 : iPrevNode + 1 ];
5922 // check if this link was already used
5923 long iLink = aLinkID_Gen.GetLinkID( prevSideNode, n );
5924 bool isJustChecked = !checkedLinkIDs.insert( iLink ).second;
5925 if (!isJustChecked &&
5926 foundSideLinkIDs.find( iLink ) == foundSideLinkIDs.end() )
5928 // test a link geometrically
5929 gp_XYZ nextXYZ ( n->X(), n->Y(), n->Z() );
5930 bool linkIsBetter = false;
5931 double dot = 0.0, dist = 0.0;
5932 if ( searchByDir ) { // choose most co-directed link
5933 dot = bordDir * ( nextXYZ - prevXYZ ).Normalized();
5934 linkIsBetter = ( dot > maxDot );
5936 else { // choose link with the node closest to bordPos
5937 dist = ( nextXYZ - bordPos ).SquareModulus();
5938 linkIsBetter = ( dist < minDist );
5940 if ( linkIsBetter ) {
5949 } // loop on inverse elements of prevSideNode
5952 MESSAGE(" Cant find path by links of the Side 2 ");
5953 return SEW_BAD_SIDE_NODES;
5955 sideNodes.push_back( sideNode );
5956 sideElems.push_back( sideElem );
5957 foundSideLinkIDs.insert ( linkID );
5958 prevSideNode = sideNode;
5960 if ( *nBordIt == theBordLastNode )
5961 searchByDir = false;
5963 // find the next border link to compare with
5964 gp_XYZ sidePos( sideNode->X(), sideNode->Y(), sideNode->Z() );
5965 searchByDir = ( bordDir * ( sidePos - bordPos ) <= 0 );
5966 // move to next border node if sideNode is before forward border node (bordPos)
5967 while ( *nBordIt != theBordLastNode && !searchByDir ) {
5968 prevBordNode = *nBordIt;
5970 bordPos = nBordXYZ[ *nBordIt ];
5971 bordDir = bordPos - nBordXYZ[ prevBordNode ];
5972 searchByDir = ( bordDir * ( sidePos - bordPos ) <= 0 );
5976 while ( sideNode != theSideSecondNode );
5978 if ( hasVolumes && sideNodes.size () != bordNodes.size() && !toCreatePolyedrs) {
5979 MESSAGE("VOLUME SPLITTING IS FORBIDDEN");
5980 return SEW_VOLUMES_TO_SPLIT; // volume splitting is forbidden
5982 } // end nodes search on the side 2
5984 // ============================
5985 // sew the border to the side 2
5986 // ============================
5988 int nbNodes[] = { nSide[0].size(), nSide[1].size() };
5989 int maxNbNodes = Max( nbNodes[0], nbNodes[1] );
5991 TListOfListOfNodes nodeGroupsToMerge;
5992 if ( nbNodes[0] == nbNodes[1] ||
5993 ( theSideIsFreeBorder && !theSideThirdNode)) {
5995 // all nodes are to be merged
5997 for (nIt[0] = nSide[0].begin(), nIt[1] = nSide[1].begin();
5998 nIt[0] != nSide[0].end() && nIt[1] != nSide[1].end();
5999 nIt[0]++, nIt[1]++ )
6001 nodeGroupsToMerge.push_back( list<const SMDS_MeshNode*>() );
6002 nodeGroupsToMerge.back().push_back( *nIt[1] ); // to keep
6003 nodeGroupsToMerge.back().push_back( *nIt[0] ); // to remove
6008 // insert new nodes into the border and the side to get equal nb of segments
6010 // get normalized parameters of nodes on the borders
6011 //double param[ 2 ][ maxNbNodes ];
6013 param[0] = new double [ maxNbNodes ];
6014 param[1] = new double [ maxNbNodes ];
6016 for ( iBord = 0; iBord < 2; iBord++ ) { // loop on 2 borders
6017 list< const SMDS_MeshNode* >& nodes = nSide[ iBord ];
6018 list< const SMDS_MeshNode* >::iterator nIt = nodes.begin();
6019 const SMDS_MeshNode* nPrev = *nIt;
6020 double bordLength = 0;
6021 for ( iNode = 0; nIt != nodes.end(); nIt++, iNode++ ) { // loop on border nodes
6022 const SMDS_MeshNode* nCur = *nIt;
6023 gp_XYZ segment (nCur->X() - nPrev->X(),
6024 nCur->Y() - nPrev->Y(),
6025 nCur->Z() - nPrev->Z());
6026 double segmentLen = segment.Modulus();
6027 bordLength += segmentLen;
6028 param[ iBord ][ iNode ] = bordLength;
6031 // normalize within [0,1]
6032 for ( iNode = 0; iNode < nbNodes[ iBord ]; iNode++ ) {
6033 param[ iBord ][ iNode ] /= bordLength;
6037 // loop on border segments
6038 const SMDS_MeshNode *nPrev[ 2 ] = { 0, 0 };
6039 int i[ 2 ] = { 0, 0 };
6040 nIt[0] = nSide[0].begin(); eIt[0] = eSide[0].begin();
6041 nIt[1] = nSide[1].begin(); eIt[1] = eSide[1].begin();
6043 TElemOfNodeListMap insertMap;
6044 TElemOfNodeListMap::iterator insertMapIt;
6046 // key: elem to insert nodes into
6047 // value: 2 nodes to insert between + nodes to be inserted
6049 bool next[ 2 ] = { false, false };
6051 // find min adjacent segment length after sewing
6052 double nextParam = 10., prevParam = 0;
6053 for ( iBord = 0; iBord < 2; iBord++ ) { // loop on 2 borders
6054 if ( i[ iBord ] + 1 < nbNodes[ iBord ])
6055 nextParam = Min( nextParam, param[iBord][ i[iBord] + 1 ]);
6056 if ( i[ iBord ] > 0 )
6057 prevParam = Max( prevParam, param[iBord][ i[iBord] - 1 ]);
6059 double minParam = Min( param[ 0 ][ i[0] ], param[ 1 ][ i[1] ]);
6060 double maxParam = Max( param[ 0 ][ i[0] ], param[ 1 ][ i[1] ]);
6061 double minSegLen = Min( nextParam - minParam, maxParam - prevParam );
6063 // choose to insert or to merge nodes
6064 double du = param[ 1 ][ i[1] ] - param[ 0 ][ i[0] ];
6065 if ( Abs( du ) <= minSegLen * 0.2 ) {
6068 nodeGroupsToMerge.push_back( list<const SMDS_MeshNode*>() );
6069 const SMDS_MeshNode* n0 = *nIt[0];
6070 const SMDS_MeshNode* n1 = *nIt[1];
6071 nodeGroupsToMerge.back().push_back( n1 );
6072 nodeGroupsToMerge.back().push_back( n0 );
6073 // position of node of the border changes due to merge
6074 param[ 0 ][ i[0] ] += du;
6075 // move n1 for the sake of elem shape evaluation during insertion.
6076 // n1 will be removed by MergeNodes() anyway
6077 const_cast<SMDS_MeshNode*>( n0 )->setXYZ( n1->X(), n1->Y(), n1->Z() );
6078 next[0] = next[1] = true;
6083 int intoBord = ( du < 0 ) ? 0 : 1;
6084 const SMDS_MeshElement* elem = *eIt[ intoBord ];
6085 const SMDS_MeshNode* n1 = nPrev[ intoBord ];
6086 const SMDS_MeshNode* n2 = *nIt[ intoBord ];
6087 const SMDS_MeshNode* nIns = *nIt[ 1 - intoBord ];
6088 if ( intoBord == 1 ) {
6089 // move node of the border to be on a link of elem of the side
6090 gp_XYZ p1 (n1->X(), n1->Y(), n1->Z());
6091 gp_XYZ p2 (n2->X(), n2->Y(), n2->Z());
6092 double ratio = du / ( param[ 1 ][ i[1] ] - param[ 1 ][ i[1]-1 ]);
6093 gp_XYZ p = p2 * ( 1 - ratio ) + p1 * ratio;
6094 GetMeshDS()->MoveNode( nIns, p.X(), p.Y(), p.Z() );
6096 insertMapIt = insertMap.find( elem );
6097 bool notFound = ( insertMapIt == insertMap.end() );
6098 bool otherLink = ( !notFound && (*insertMapIt).second.front() != n1 );
6100 // insert into another link of the same element:
6101 // 1. perform insertion into the other link of the elem
6102 list<const SMDS_MeshNode*> & nodeList = (*insertMapIt).second;
6103 const SMDS_MeshNode* n12 = nodeList.front(); nodeList.pop_front();
6104 const SMDS_MeshNode* n22 = nodeList.front(); nodeList.pop_front();
6105 InsertNodesIntoLink( elem, n12, n22, nodeList, toCreatePolygons );
6106 // 2. perform insertion into the link of adjacent faces
6108 const SMDS_MeshElement* adjElem = findAdjacentFace( n12, n22, elem );
6110 InsertNodesIntoLink( adjElem, n12, n22, nodeList, toCreatePolygons );
6114 if (toCreatePolyedrs) {
6115 // perform insertion into the links of adjacent volumes
6116 UpdateVolumes(n12, n22, nodeList);
6118 // 3. find an element appeared on n1 and n2 after the insertion
6119 insertMap.erase( elem );
6120 elem = findAdjacentFace( n1, n2, 0 );
6122 if ( notFound || otherLink ) {
6123 // add element and nodes of the side into the insertMap
6124 insertMapIt = insertMap.insert
6125 ( TElemOfNodeListMap::value_type( elem, list<const SMDS_MeshNode*>() )).first;
6126 (*insertMapIt).second.push_back( n1 );
6127 (*insertMapIt).second.push_back( n2 );
6129 // add node to be inserted into elem
6130 (*insertMapIt).second.push_back( nIns );
6131 next[ 1 - intoBord ] = true;
6134 // go to the next segment
6135 for ( iBord = 0; iBord < 2; iBord++ ) { // loop on 2 borders
6136 if ( next[ iBord ] ) {
6137 if ( i[ iBord ] != 0 && eIt[ iBord ] != eSide[ iBord ].end())
6139 nPrev[ iBord ] = *nIt[ iBord ];
6140 nIt[ iBord ]++; i[ iBord ]++;
6144 while ( nIt[0] != nSide[0].end() && nIt[1] != nSide[1].end());
6146 // perform insertion of nodes into elements
6148 for (insertMapIt = insertMap.begin();
6149 insertMapIt != insertMap.end();
6152 const SMDS_MeshElement* elem = (*insertMapIt).first;
6153 list<const SMDS_MeshNode*> & nodeList = (*insertMapIt).second;
6154 const SMDS_MeshNode* n1 = nodeList.front(); nodeList.pop_front();
6155 const SMDS_MeshNode* n2 = nodeList.front(); nodeList.pop_front();
6157 InsertNodesIntoLink( elem, n1, n2, nodeList, toCreatePolygons );
6159 if ( !theSideIsFreeBorder ) {
6160 // look for and insert nodes into the faces adjacent to elem
6162 const SMDS_MeshElement* adjElem = findAdjacentFace( n1, n2, elem );
6164 InsertNodesIntoLink( adjElem, n1, n2, nodeList, toCreatePolygons );
6169 if (toCreatePolyedrs) {
6170 // perform insertion into the links of adjacent volumes
6171 UpdateVolumes(n1, n2, nodeList);
6177 } // end: insert new nodes
6179 MergeNodes ( nodeGroupsToMerge );
6184 //=======================================================================
6185 //function : InsertNodesIntoLink
6186 //purpose : insert theNodesToInsert into theFace between theBetweenNode1
6187 // and theBetweenNode2 and split theElement
6188 //=======================================================================
6190 void SMESH_MeshEditor::InsertNodesIntoLink(const SMDS_MeshElement* theFace,
6191 const SMDS_MeshNode* theBetweenNode1,
6192 const SMDS_MeshNode* theBetweenNode2,
6193 list<const SMDS_MeshNode*>& theNodesToInsert,
6194 const bool toCreatePoly)
6196 if ( theFace->GetType() != SMDSAbs_Face ) return;
6198 // find indices of 2 link nodes and of the rest nodes
6199 int iNode = 0, il1, il2, i3, i4;
6200 il1 = il2 = i3 = i4 = -1;
6201 //const SMDS_MeshNode* nodes[ theFace->NbNodes() ];
6202 vector<const SMDS_MeshNode*> nodes( theFace->NbNodes() );
6204 if(theFace->IsQuadratic()) {
6205 const SMDS_QuadraticFaceOfNodes* F =
6206 static_cast<const SMDS_QuadraticFaceOfNodes*>(theFace);
6207 // use special nodes iterator
6208 SMDS_NodeIteratorPtr anIter = F->interlacedNodesIterator();
6209 while( anIter->more() ) {
6210 const SMDS_MeshNode* n = anIter->next();
6211 if ( n == theBetweenNode1 )
6213 else if ( n == theBetweenNode2 )
6219 nodes[ iNode++ ] = n;
6223 SMDS_ElemIteratorPtr nodeIt = theFace->nodesIterator();
6224 while ( nodeIt->more() ) {
6225 const SMDS_MeshNode* n = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
6226 if ( n == theBetweenNode1 )
6228 else if ( n == theBetweenNode2 )
6234 nodes[ iNode++ ] = n;
6237 if ( il1 < 0 || il2 < 0 || i3 < 0 )
6240 // arrange link nodes to go one after another regarding the face orientation
6241 bool reverse = ( Abs( il2 - il1 ) == 1 ? il2 < il1 : il1 < il2 );
6242 list<const SMDS_MeshNode *> aNodesToInsert = theNodesToInsert;
6247 aNodesToInsert.reverse();
6249 // check that not link nodes of a quadrangles are in good order
6250 int nbFaceNodes = theFace->NbNodes();
6251 if ( nbFaceNodes == 4 && i4 - i3 != 1 ) {
6257 if (toCreatePoly || theFace->IsPoly()) {
6260 vector<const SMDS_MeshNode *> poly_nodes (nbFaceNodes + aNodesToInsert.size());
6262 // add nodes of face up to first node of link
6265 if(theFace->IsQuadratic()) {
6266 const SMDS_QuadraticFaceOfNodes* F =
6267 static_cast<const SMDS_QuadraticFaceOfNodes*>(theFace);
6268 // use special nodes iterator
6269 SMDS_NodeIteratorPtr anIter = F->interlacedNodesIterator();
6270 while( anIter->more() && !isFLN ) {
6271 const SMDS_MeshNode* n = anIter->next();
6272 poly_nodes[iNode++] = n;
6273 if (n == nodes[il1]) {
6277 // add nodes to insert
6278 list<const SMDS_MeshNode*>::iterator nIt = aNodesToInsert.begin();
6279 for (; nIt != aNodesToInsert.end(); nIt++) {
6280 poly_nodes[iNode++] = *nIt;
6282 // add nodes of face starting from last node of link
6283 while ( anIter->more() ) {
6284 poly_nodes[iNode++] = anIter->next();
6288 SMDS_ElemIteratorPtr nodeIt = theFace->nodesIterator();
6289 while ( nodeIt->more() && !isFLN ) {
6290 const SMDS_MeshNode* n = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
6291 poly_nodes[iNode++] = n;
6292 if (n == nodes[il1]) {
6296 // add nodes to insert
6297 list<const SMDS_MeshNode*>::iterator nIt = aNodesToInsert.begin();
6298 for (; nIt != aNodesToInsert.end(); nIt++) {
6299 poly_nodes[iNode++] = *nIt;
6301 // add nodes of face starting from last node of link
6302 while ( nodeIt->more() ) {
6303 const SMDS_MeshNode* n = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
6304 poly_nodes[iNode++] = n;
6308 // edit or replace the face
6309 SMESHDS_Mesh *aMesh = GetMeshDS();
6311 if (theFace->IsPoly()) {
6312 aMesh->ChangePolygonNodes(theFace, poly_nodes);
6315 int aShapeId = FindShape( theFace );
6317 SMDS_MeshElement* newElem = aMesh->AddPolygonalFace(poly_nodes);
6318 myLastCreatedElems.Append(newElem);
6319 if ( aShapeId && newElem )
6320 aMesh->SetMeshElementOnShape( newElem, aShapeId );
6322 aMesh->RemoveElement(theFace);
6327 if( !theFace->IsQuadratic() ) {
6329 // put aNodesToInsert between theBetweenNode1 and theBetweenNode2
6330 int nbLinkNodes = 2 + aNodesToInsert.size();
6331 //const SMDS_MeshNode* linkNodes[ nbLinkNodes ];
6332 vector<const SMDS_MeshNode*> linkNodes( nbLinkNodes );
6333 linkNodes[ 0 ] = nodes[ il1 ];
6334 linkNodes[ nbLinkNodes - 1 ] = nodes[ il2 ];
6335 list<const SMDS_MeshNode*>::iterator nIt = aNodesToInsert.begin();
6336 for ( iNode = 1; nIt != aNodesToInsert.end(); nIt++ ) {
6337 linkNodes[ iNode++ ] = *nIt;
6339 // decide how to split a quadrangle: compare possible variants
6340 // and choose which of splits to be a quadrangle
6341 int i1, i2, iSplit, nbSplits = nbLinkNodes - 1, iBestQuad;
6342 if ( nbFaceNodes == 3 ) {
6343 iBestQuad = nbSplits;
6346 else if ( nbFaceNodes == 4 ) {
6347 SMESH::Controls::NumericalFunctorPtr aCrit( new SMESH::Controls::AspectRatio);
6348 double aBestRate = DBL_MAX;
6349 for ( int iQuad = 0; iQuad < nbSplits; iQuad++ ) {
6351 double aBadRate = 0;
6352 // evaluate elements quality
6353 for ( iSplit = 0; iSplit < nbSplits; iSplit++ ) {
6354 if ( iSplit == iQuad ) {
6355 SMDS_FaceOfNodes quad (linkNodes[ i1++ ],
6359 aBadRate += getBadRate( &quad, aCrit );
6362 SMDS_FaceOfNodes tria (linkNodes[ i1++ ],
6364 nodes[ iSplit < iQuad ? i4 : i3 ]);
6365 aBadRate += getBadRate( &tria, aCrit );
6369 if ( aBadRate < aBestRate ) {
6371 aBestRate = aBadRate;
6376 // create new elements
6377 SMESHDS_Mesh *aMesh = GetMeshDS();
6378 int aShapeId = FindShape( theFace );
6381 for ( iSplit = 0; iSplit < nbSplits - 1; iSplit++ ) {
6382 SMDS_MeshElement* newElem = 0;
6383 if ( iSplit == iBestQuad )
6384 newElem = aMesh->AddFace (linkNodes[ i1++ ],
6389 newElem = aMesh->AddFace (linkNodes[ i1++ ],
6391 nodes[ iSplit < iBestQuad ? i4 : i3 ]);
6392 myLastCreatedElems.Append(newElem);
6393 if ( aShapeId && newElem )
6394 aMesh->SetMeshElementOnShape( newElem, aShapeId );
6397 // change nodes of theFace
6398 const SMDS_MeshNode* newNodes[ 4 ];
6399 newNodes[ 0 ] = linkNodes[ i1 ];
6400 newNodes[ 1 ] = linkNodes[ i2 ];
6401 newNodes[ 2 ] = nodes[ iSplit >= iBestQuad ? i3 : i4 ];
6402 newNodes[ 3 ] = nodes[ i4 ];
6403 aMesh->ChangeElementNodes( theFace, newNodes, iSplit == iBestQuad ? 4 : 3 );
6404 } // end if(!theFace->IsQuadratic())
6405 else { // theFace is quadratic
6406 // we have to split theFace on simple triangles and one simple quadrangle
6408 int nbshift = tmp*2;
6409 // shift nodes in nodes[] by nbshift
6411 for(i=0; i<nbshift; i++) {
6412 const SMDS_MeshNode* n = nodes[0];
6413 for(j=0; j<nbFaceNodes-1; j++) {
6414 nodes[j] = nodes[j+1];
6416 nodes[nbFaceNodes-1] = n;
6418 il1 = il1 - nbshift;
6419 // now have to insert nodes between n0 and n1 or n1 and n2 (see below)
6420 // n0 n1 n2 n0 n1 n2
6421 // +-----+-----+ +-----+-----+
6430 // create new elements
6431 SMESHDS_Mesh *aMesh = GetMeshDS();
6432 int aShapeId = FindShape( theFace );
6435 if(nbFaceNodes==6) { // quadratic triangle
6436 SMDS_MeshElement* newElem =
6437 aMesh->AddFace(nodes[3],nodes[4],nodes[5]);
6438 myLastCreatedElems.Append(newElem);
6439 if ( aShapeId && newElem )
6440 aMesh->SetMeshElementOnShape( newElem, aShapeId );
6441 if(theFace->IsMediumNode(nodes[il1])) {
6442 // create quadrangle
6443 newElem = aMesh->AddFace(nodes[0],nodes[1],nodes[3],nodes[5]);
6444 myLastCreatedElems.Append(newElem);
6445 if ( aShapeId && newElem )
6446 aMesh->SetMeshElementOnShape( newElem, aShapeId );
6452 // create quadrangle
6453 newElem = aMesh->AddFace(nodes[1],nodes[2],nodes[3],nodes[5]);
6454 myLastCreatedElems.Append(newElem);
6455 if ( aShapeId && newElem )
6456 aMesh->SetMeshElementOnShape( newElem, aShapeId );
6462 else { // nbFaceNodes==8 - quadratic quadrangle
6463 SMDS_MeshElement* newElem =
6464 aMesh->AddFace(nodes[3],nodes[4],nodes[5]);
6465 myLastCreatedElems.Append(newElem);
6466 if ( aShapeId && newElem )
6467 aMesh->SetMeshElementOnShape( newElem, aShapeId );
6468 newElem = aMesh->AddFace(nodes[5],nodes[6],nodes[7]);
6469 myLastCreatedElems.Append(newElem);
6470 if ( aShapeId && newElem )
6471 aMesh->SetMeshElementOnShape( newElem, aShapeId );
6472 newElem = aMesh->AddFace(nodes[5],nodes[7],nodes[3]);
6473 myLastCreatedElems.Append(newElem);
6474 if ( aShapeId && newElem )
6475 aMesh->SetMeshElementOnShape( newElem, aShapeId );
6476 if(theFace->IsMediumNode(nodes[il1])) {
6477 // create quadrangle
6478 newElem = aMesh->AddFace(nodes[0],nodes[1],nodes[3],nodes[7]);
6479 myLastCreatedElems.Append(newElem);
6480 if ( aShapeId && newElem )
6481 aMesh->SetMeshElementOnShape( newElem, aShapeId );
6487 // create quadrangle
6488 newElem = aMesh->AddFace(nodes[1],nodes[2],nodes[3],nodes[7]);
6489 myLastCreatedElems.Append(newElem);
6490 if ( aShapeId && newElem )
6491 aMesh->SetMeshElementOnShape( newElem, aShapeId );
6497 // create needed triangles using n1,n2,n3 and inserted nodes
6498 int nbn = 2 + aNodesToInsert.size();
6499 //const SMDS_MeshNode* aNodes[nbn];
6500 vector<const SMDS_MeshNode*> aNodes(nbn);
6501 aNodes[0] = nodes[n1];
6502 aNodes[nbn-1] = nodes[n2];
6503 list<const SMDS_MeshNode*>::iterator nIt = aNodesToInsert.begin();
6504 for ( iNode = 1; nIt != aNodesToInsert.end(); nIt++ ) {
6505 aNodes[iNode++] = *nIt;
6507 for(i=1; i<nbn; i++) {
6508 SMDS_MeshElement* newElem =
6509 aMesh->AddFace(aNodes[i-1],aNodes[i],nodes[n3]);
6510 myLastCreatedElems.Append(newElem);
6511 if ( aShapeId && newElem )
6512 aMesh->SetMeshElementOnShape( newElem, aShapeId );
6514 // remove old quadratic face
6515 aMesh->RemoveElement(theFace);
6519 //=======================================================================
6520 //function : UpdateVolumes
6522 //=======================================================================
6523 void SMESH_MeshEditor::UpdateVolumes (const SMDS_MeshNode* theBetweenNode1,
6524 const SMDS_MeshNode* theBetweenNode2,
6525 list<const SMDS_MeshNode*>& theNodesToInsert)
6527 myLastCreatedElems.Clear();
6528 myLastCreatedNodes.Clear();
6530 SMDS_ElemIteratorPtr invElemIt = theBetweenNode1->GetInverseElementIterator(SMDSAbs_Volume);
6531 while (invElemIt->more()) { // loop on inverse elements of theBetweenNode1
6532 const SMDS_MeshElement* elem = invElemIt->next();
6534 // check, if current volume has link theBetweenNode1 - theBetweenNode2
6535 SMDS_VolumeTool aVolume (elem);
6536 if (!aVolume.IsLinked(theBetweenNode1, theBetweenNode2))
6539 // insert new nodes in all faces of the volume, sharing link theBetweenNode1 - theBetweenNode2
6540 int iface, nbFaces = aVolume.NbFaces();
6541 vector<const SMDS_MeshNode *> poly_nodes;
6542 vector<int> quantities (nbFaces);
6544 for (iface = 0; iface < nbFaces; iface++) {
6545 int nbFaceNodes = aVolume.NbFaceNodes(iface), nbInserted = 0;
6546 // faceNodes will contain (nbFaceNodes + 1) nodes, last = first
6547 const SMDS_MeshNode** faceNodes = aVolume.GetFaceNodes(iface);
6549 for (int inode = 0; inode < nbFaceNodes; inode++) {
6550 poly_nodes.push_back(faceNodes[inode]);
6552 if (nbInserted == 0) {
6553 if (faceNodes[inode] == theBetweenNode1) {
6554 if (faceNodes[inode + 1] == theBetweenNode2) {
6555 nbInserted = theNodesToInsert.size();
6557 // add nodes to insert
6558 list<const SMDS_MeshNode*>::iterator nIt = theNodesToInsert.begin();
6559 for (; nIt != theNodesToInsert.end(); nIt++) {
6560 poly_nodes.push_back(*nIt);
6564 else if (faceNodes[inode] == theBetweenNode2) {
6565 if (faceNodes[inode + 1] == theBetweenNode1) {
6566 nbInserted = theNodesToInsert.size();
6568 // add nodes to insert in reversed order
6569 list<const SMDS_MeshNode*>::iterator nIt = theNodesToInsert.end();
6571 for (; nIt != theNodesToInsert.begin(); nIt--) {
6572 poly_nodes.push_back(*nIt);
6574 poly_nodes.push_back(*nIt);
6581 quantities[iface] = nbFaceNodes + nbInserted;
6584 // Replace or update the volume
6585 SMESHDS_Mesh *aMesh = GetMeshDS();
6587 if (elem->IsPoly()) {
6588 aMesh->ChangePolyhedronNodes(elem, poly_nodes, quantities);
6592 int aShapeId = FindShape( elem );
6594 SMDS_MeshElement* newElem =
6595 aMesh->AddPolyhedralVolume(poly_nodes, quantities);
6596 myLastCreatedElems.Append(newElem);
6597 if (aShapeId && newElem)
6598 aMesh->SetMeshElementOnShape(newElem, aShapeId);
6600 aMesh->RemoveElement(elem);
6605 //=======================================================================
6607 * \brief Convert elements contained in a submesh to quadratic
6608 * \retval int - nb of checked elements
6610 //=======================================================================
6612 int SMESH_MeshEditor::convertElemToQuadratic(SMESHDS_SubMesh * theSm,
6613 SMESH_MesherHelper& theHelper,
6614 const bool theForce3d)
6617 if( !theSm ) return nbElem;
6618 SMDS_ElemIteratorPtr ElemItr = theSm->GetElements();
6619 while(ElemItr->more())
6622 const SMDS_MeshElement* elem = ElemItr->next();
6623 if( !elem || elem->IsQuadratic() ) continue;
6625 int id = elem->GetID();
6626 int nbNodes = elem->NbNodes();
6627 vector<const SMDS_MeshNode *> aNds (nbNodes);
6629 for(int i = 0; i < nbNodes; i++)
6631 aNds[i] = elem->GetNode(i);
6633 SMDSAbs_ElementType aType = elem->GetType();
6635 theSm->RemoveElement(elem);
6636 GetMeshDS()->SMDS_Mesh::RemoveFreeElement(elem);
6638 const SMDS_MeshElement* NewElem = 0;
6644 NewElem = theHelper.AddEdge(aNds[0], aNds[1], id, theForce3d);
6652 NewElem = theHelper.AddFace(aNds[0], aNds[1], aNds[2], id, theForce3d);
6655 NewElem = theHelper.AddFace(aNds[0], aNds[1], aNds[2], aNds[3], id, theForce3d);
6662 case SMDSAbs_Volume :
6667 NewElem = theHelper.AddVolume(aNds[0], aNds[1], aNds[2], aNds[3], id, true);
6670 NewElem = theHelper.AddVolume(aNds[0], aNds[1], aNds[2], aNds[3], aNds[4], aNds[5], id, true);
6673 NewElem = theHelper.AddVolume(aNds[0], aNds[1], aNds[2], aNds[3],
6674 aNds[4], aNds[5], aNds[6], aNds[7], id, true);
6686 AddToSameGroups( NewElem, elem, GetMeshDS());
6687 theSm->AddElement( NewElem );
6689 if ( NewElem != elem )
6690 RemoveElemFromGroups (elem, GetMeshDS());
6695 //=======================================================================
6696 //function : ConvertToQuadratic
6698 //=======================================================================
6699 void SMESH_MeshEditor::ConvertToQuadratic(const bool theForce3d)
6701 SMESHDS_Mesh* meshDS = GetMeshDS();
6703 SMESH_MesherHelper aHelper(*myMesh);
6704 aHelper.SetIsQuadratic( true );
6706 int nbCheckedElems = 0;
6707 if ( myMesh->HasShapeToMesh() )
6709 if ( SMESH_subMesh *aSubMesh = myMesh->GetSubMeshContaining(myMesh->GetShapeToMesh()))
6711 SMESH_subMeshIteratorPtr smIt = aSubMesh->getDependsOnIterator(true,false);
6712 while ( smIt->more() ) {
6713 SMESH_subMesh* sm = smIt->next();
6714 if ( SMESHDS_SubMesh *smDS = sm->GetSubMeshDS() ) {
6715 aHelper.SetSubShape( sm->GetSubShape() );
6716 nbCheckedElems += convertElemToQuadratic(smDS, aHelper, theForce3d);
6721 int totalNbElems = meshDS->NbEdges() + meshDS->NbFaces() + meshDS->NbVolumes();
6722 if ( nbCheckedElems < totalNbElems ) // not all elements in submeshes
6724 SMDS_EdgeIteratorPtr aEdgeItr = meshDS->edgesIterator();
6725 while(aEdgeItr->more())
6727 const SMDS_MeshEdge* edge = aEdgeItr->next();
6728 if(edge && !edge->IsQuadratic())
6730 int id = edge->GetID();
6731 const SMDS_MeshNode* n1 = edge->GetNode(0);
6732 const SMDS_MeshNode* n2 = edge->GetNode(1);
6734 meshDS->SMDS_Mesh::RemoveFreeElement(edge);
6736 const SMDS_MeshEdge* NewEdge = aHelper.AddEdge(n1, n2, id, theForce3d);
6738 AddToSameGroups(NewEdge, edge, meshDS);
6739 if ( NewEdge != edge )
6740 RemoveElemFromGroups (edge, meshDS);
6743 SMDS_FaceIteratorPtr aFaceItr = meshDS->facesIterator();
6744 while(aFaceItr->more())
6746 const SMDS_MeshFace* face = aFaceItr->next();
6747 if(!face || face->IsQuadratic() ) continue;
6749 int id = face->GetID();
6750 int nbNodes = face->NbNodes();
6751 vector<const SMDS_MeshNode *> aNds (nbNodes);
6753 for(int i = 0; i < nbNodes; i++)
6755 aNds[i] = face->GetNode(i);
6758 meshDS->SMDS_Mesh::RemoveFreeElement(face);
6760 SMDS_MeshFace * NewFace = 0;
6764 NewFace = aHelper.AddFace(aNds[0], aNds[1], aNds[2], id, theForce3d);
6767 NewFace = aHelper.AddFace(aNds[0], aNds[1], aNds[2], aNds[3], id, theForce3d);
6773 AddToSameGroups(NewFace, face, meshDS);
6774 if ( NewFace != face )
6775 RemoveElemFromGroups (face, meshDS);
6777 SMDS_VolumeIteratorPtr aVolumeItr = meshDS->volumesIterator();
6778 while(aVolumeItr->more())
6780 const SMDS_MeshVolume* volume = aVolumeItr->next();
6781 if(!volume || volume->IsQuadratic() ) continue;
6783 int id = volume->GetID();
6784 int nbNodes = volume->NbNodes();
6785 vector<const SMDS_MeshNode *> aNds (nbNodes);
6787 for(int i = 0; i < nbNodes; i++)
6789 aNds[i] = volume->GetNode(i);
6792 meshDS->SMDS_Mesh::RemoveFreeElement(volume);
6794 SMDS_MeshVolume * NewVolume = 0;
6798 NewVolume = aHelper.AddVolume(aNds[0], aNds[1], aNds[2],
6799 aNds[3], id, true );
6802 NewVolume = aHelper.AddVolume(aNds[0], aNds[1], aNds[2],
6803 aNds[3], aNds[4], aNds[5], id, true);
6806 NewVolume = aHelper.AddVolume(aNds[0], aNds[1], aNds[2], aNds[3],
6807 aNds[4], aNds[5], aNds[6], aNds[7], id, true);
6813 AddToSameGroups(NewVolume, volume, meshDS);
6814 if ( NewVolume != volume )
6815 RemoveElemFromGroups (volume, meshDS);
6820 //=======================================================================
6822 * \brief Convert quadratic elements to linear ones and remove quadratic nodes
6823 * \retval int - nb of checked elements
6825 //=======================================================================
6827 int SMESH_MeshEditor::removeQuadElem(SMESHDS_SubMesh * theSm,
6828 SMDS_ElemIteratorPtr theItr,
6829 const int theShapeID)
6832 SMESHDS_Mesh* meshDS = GetMeshDS();
6833 while( theItr->more() )
6835 const SMDS_MeshElement* elem = theItr->next();
6837 if( elem && elem->IsQuadratic())
6839 int id = elem->GetID();
6840 int nbNodes = elem->NbNodes();
6841 vector<const SMDS_MeshNode *> aNds, mediumNodes;
6842 aNds.reserve( nbNodes );
6843 mediumNodes.reserve( nbNodes );
6845 for(int i = 0; i < nbNodes; i++)
6847 const SMDS_MeshNode* n = elem->GetNode(i);
6849 if( elem->IsMediumNode( n ) )
6850 mediumNodes.push_back( n );
6852 aNds.push_back( n );
6854 if( aNds.empty() ) continue;
6855 SMDSAbs_ElementType aType = elem->GetType();
6857 //remove old quadratic element
6858 meshDS->SMDS_Mesh::RemoveFreeElement( elem );
6860 theSm->RemoveElement( elem );
6862 SMDS_MeshElement * NewElem = AddElement( aNds, aType, false, id );
6864 AddToSameGroups(NewElem, elem, meshDS);
6865 if ( NewElem != elem )
6866 RemoveElemFromGroups (elem, meshDS);
6867 if( theSm && NewElem )
6868 theSm->AddElement( NewElem );
6870 // remove medium nodes
6871 vector<const SMDS_MeshNode*>::iterator nIt = mediumNodes.begin();
6872 for ( ; nIt != mediumNodes.end(); ++nIt ) {
6873 const SMDS_MeshNode* n = *nIt;
6874 if ( n->NbInverseElements() == 0 ) {
6875 if ( n->GetPosition()->GetShapeId() != theShapeID )
6876 meshDS->RemoveFreeNode( n, meshDS->MeshElements
6877 ( n->GetPosition()->GetShapeId() ));
6879 meshDS->RemoveFreeNode( n, theSm );
6887 //=======================================================================
6888 //function : ConvertFromQuadratic
6890 //=======================================================================
6891 bool SMESH_MeshEditor::ConvertFromQuadratic()
6893 int nbCheckedElems = 0;
6894 if ( myMesh->HasShapeToMesh() )
6896 if ( SMESH_subMesh *aSubMesh = myMesh->GetSubMeshContaining(myMesh->GetShapeToMesh()))
6898 SMESH_subMeshIteratorPtr smIt = aSubMesh->getDependsOnIterator(true,false);
6899 while ( smIt->more() ) {
6900 SMESH_subMesh* sm = smIt->next();
6901 if ( SMESHDS_SubMesh *smDS = sm->GetSubMeshDS() )
6902 nbCheckedElems += removeQuadElem( smDS, smDS->GetElements(), sm->GetId() );
6908 GetMeshDS()->NbEdges() + GetMeshDS()->NbFaces() + GetMeshDS()->NbVolumes();
6909 if ( nbCheckedElems < totalNbElems ) // not all elements in submeshes
6911 SMESHDS_SubMesh *aSM = 0;
6912 removeQuadElem( aSM, GetMeshDS()->elementsIterator(), 0 );
6918 //=======================================================================
6919 //function : SewSideElements
6921 //=======================================================================
6923 SMESH_MeshEditor::Sew_Error
6924 SMESH_MeshEditor::SewSideElements (TIDSortedElemSet& theSide1,
6925 TIDSortedElemSet& theSide2,
6926 const SMDS_MeshNode* theFirstNode1,
6927 const SMDS_MeshNode* theFirstNode2,
6928 const SMDS_MeshNode* theSecondNode1,
6929 const SMDS_MeshNode* theSecondNode2)
6931 myLastCreatedElems.Clear();
6932 myLastCreatedNodes.Clear();
6934 MESSAGE ("::::SewSideElements()");
6935 if ( theSide1.size() != theSide2.size() )
6936 return SEW_DIFF_NB_OF_ELEMENTS;
6938 Sew_Error aResult = SEW_OK;
6940 // 1. Build set of faces representing each side
6941 // 2. Find which nodes of the side 1 to merge with ones on the side 2
6942 // 3. Replace nodes in elements of the side 1 and remove replaced nodes
6944 // =======================================================================
6945 // 1. Build set of faces representing each side:
6946 // =======================================================================
6947 // a. build set of nodes belonging to faces
6948 // b. complete set of faces: find missing fices whose nodes are in set of nodes
6949 // c. create temporary faces representing side of volumes if correspondent
6950 // face does not exist
6952 SMESHDS_Mesh* aMesh = GetMeshDS();
6953 SMDS_Mesh aTmpFacesMesh;
6954 set<const SMDS_MeshElement*> faceSet1, faceSet2;
6955 set<const SMDS_MeshElement*> volSet1, volSet2;
6956 set<const SMDS_MeshNode*> nodeSet1, nodeSet2;
6957 set<const SMDS_MeshElement*> * faceSetPtr[] = { &faceSet1, &faceSet2 };
6958 set<const SMDS_MeshElement*> * volSetPtr[] = { &volSet1, &volSet2 };
6959 set<const SMDS_MeshNode*> * nodeSetPtr[] = { &nodeSet1, &nodeSet2 };
6960 TIDSortedElemSet * elemSetPtr[] = { &theSide1, &theSide2 };
6961 int iSide, iFace, iNode;
6963 for ( iSide = 0; iSide < 2; iSide++ ) {
6964 set<const SMDS_MeshNode*> * nodeSet = nodeSetPtr[ iSide ];
6965 TIDSortedElemSet * elemSet = elemSetPtr[ iSide ];
6966 set<const SMDS_MeshElement*> * faceSet = faceSetPtr[ iSide ];
6967 set<const SMDS_MeshElement*> * volSet = volSetPtr [ iSide ];
6968 set<const SMDS_MeshElement*>::iterator vIt;
6969 TIDSortedElemSet::iterator eIt;
6970 set<const SMDS_MeshNode*>::iterator nIt;
6972 // check that given nodes belong to given elements
6973 const SMDS_MeshNode* n1 = ( iSide == 0 ) ? theFirstNode1 : theFirstNode2;
6974 const SMDS_MeshNode* n2 = ( iSide == 0 ) ? theSecondNode1 : theSecondNode2;
6975 int firstIndex = -1, secondIndex = -1;
6976 for (eIt = elemSet->begin(); eIt != elemSet->end(); eIt++ ) {
6977 const SMDS_MeshElement* elem = *eIt;
6978 if ( firstIndex < 0 ) firstIndex = elem->GetNodeIndex( n1 );
6979 if ( secondIndex < 0 ) secondIndex = elem->GetNodeIndex( n2 );
6980 if ( firstIndex > -1 && secondIndex > -1 ) break;
6982 if ( firstIndex < 0 || secondIndex < 0 ) {
6983 // we can simply return until temporary faces created
6984 return (iSide == 0 ) ? SEW_BAD_SIDE1_NODES : SEW_BAD_SIDE2_NODES;
6987 // -----------------------------------------------------------
6988 // 1a. Collect nodes of existing faces
6989 // and build set of face nodes in order to detect missing
6990 // faces corresponing to sides of volumes
6991 // -----------------------------------------------------------
6993 set< set <const SMDS_MeshNode*> > setOfFaceNodeSet;
6995 // loop on the given element of a side
6996 for (eIt = elemSet->begin(); eIt != elemSet->end(); eIt++ ) {
6997 //const SMDS_MeshElement* elem = *eIt;
6998 const SMDS_MeshElement* elem = *eIt;
6999 if ( elem->GetType() == SMDSAbs_Face ) {
7000 faceSet->insert( elem );
7001 set <const SMDS_MeshNode*> faceNodeSet;
7002 SMDS_ElemIteratorPtr nodeIt = elem->nodesIterator();
7003 while ( nodeIt->more() ) {
7004 const SMDS_MeshNode* n = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
7005 nodeSet->insert( n );
7006 faceNodeSet.insert( n );
7008 setOfFaceNodeSet.insert( faceNodeSet );
7010 else if ( elem->GetType() == SMDSAbs_Volume )
7011 volSet->insert( elem );
7013 // ------------------------------------------------------------------------------
7014 // 1b. Complete set of faces: find missing fices whose nodes are in set of nodes
7015 // ------------------------------------------------------------------------------
7017 for ( nIt = nodeSet->begin(); nIt != nodeSet->end(); nIt++ ) { // loop on nodes of iSide
7018 SMDS_ElemIteratorPtr fIt = (*nIt)->GetInverseElementIterator(SMDSAbs_Face);
7019 while ( fIt->more() ) { // loop on faces sharing a node
7020 const SMDS_MeshElement* f = fIt->next();
7021 if ( faceSet->find( f ) == faceSet->end() ) {
7022 // check if all nodes are in nodeSet and
7023 // complete setOfFaceNodeSet if they are
7024 set <const SMDS_MeshNode*> faceNodeSet;
7025 SMDS_ElemIteratorPtr nodeIt = f->nodesIterator();
7026 bool allInSet = true;
7027 while ( nodeIt->more() && allInSet ) { // loop on nodes of a face
7028 const SMDS_MeshNode* n = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
7029 if ( nodeSet->find( n ) == nodeSet->end() )
7032 faceNodeSet.insert( n );
7035 faceSet->insert( f );
7036 setOfFaceNodeSet.insert( faceNodeSet );
7042 // -------------------------------------------------------------------------
7043 // 1c. Create temporary faces representing sides of volumes if correspondent
7044 // face does not exist
7045 // -------------------------------------------------------------------------
7047 if ( !volSet->empty() ) {
7048 //int nodeSetSize = nodeSet->size();
7050 // loop on given volumes
7051 for ( vIt = volSet->begin(); vIt != volSet->end(); vIt++ ) {
7052 SMDS_VolumeTool vol (*vIt);
7053 // loop on volume faces: find free faces
7054 // --------------------------------------
7055 list<const SMDS_MeshElement* > freeFaceList;
7056 for ( iFace = 0; iFace < vol.NbFaces(); iFace++ ) {
7057 if ( !vol.IsFreeFace( iFace ))
7059 // check if there is already a face with same nodes in a face set
7060 const SMDS_MeshElement* aFreeFace = 0;
7061 const SMDS_MeshNode** fNodes = vol.GetFaceNodes( iFace );
7062 int nbNodes = vol.NbFaceNodes( iFace );
7063 set <const SMDS_MeshNode*> faceNodeSet;
7064 vol.GetFaceNodes( iFace, faceNodeSet );
7065 bool isNewFace = setOfFaceNodeSet.insert( faceNodeSet ).second;
7067 // no such a face is given but it still can exist, check it
7068 if ( nbNodes == 3 ) {
7069 aFreeFace = aMesh->FindFace( fNodes[0],fNodes[1],fNodes[2] );
7071 else if ( nbNodes == 4 ) {
7072 aFreeFace = aMesh->FindFace( fNodes[0],fNodes[1],fNodes[2],fNodes[3] );
7075 vector<const SMDS_MeshNode *> poly_nodes ( fNodes, & fNodes[nbNodes]);
7076 aFreeFace = aMesh->FindFace(poly_nodes);
7080 // create a temporary face
7081 if ( nbNodes == 3 ) {
7082 aFreeFace = aTmpFacesMesh.AddFace( fNodes[0],fNodes[1],fNodes[2] );
7084 else if ( nbNodes == 4 ) {
7085 aFreeFace = aTmpFacesMesh.AddFace( fNodes[0],fNodes[1],fNodes[2],fNodes[3] );
7088 vector<const SMDS_MeshNode *> poly_nodes ( fNodes, & fNodes[nbNodes]);
7089 aFreeFace = aTmpFacesMesh.AddPolygonalFace(poly_nodes);
7093 freeFaceList.push_back( aFreeFace );
7095 } // loop on faces of a volume
7097 // choose one of several free faces
7098 // --------------------------------------
7099 if ( freeFaceList.size() > 1 ) {
7100 // choose a face having max nb of nodes shared by other elems of a side
7101 int maxNbNodes = -1/*, nbExcludedFaces = 0*/;
7102 list<const SMDS_MeshElement* >::iterator fIt = freeFaceList.begin();
7103 while ( fIt != freeFaceList.end() ) { // loop on free faces
7104 int nbSharedNodes = 0;
7105 SMDS_ElemIteratorPtr nodeIt = (*fIt)->nodesIterator();
7106 while ( nodeIt->more() ) { // loop on free face nodes
7107 const SMDS_MeshNode* n =
7108 static_cast<const SMDS_MeshNode*>( nodeIt->next() );
7109 SMDS_ElemIteratorPtr invElemIt = n->GetInverseElementIterator();
7110 while ( invElemIt->more() ) {
7111 const SMDS_MeshElement* e = invElemIt->next();
7112 if ( faceSet->find( e ) != faceSet->end() )
7114 if ( elemSet->find( e ) != elemSet->end() )
7118 if ( nbSharedNodes >= maxNbNodes ) {
7119 maxNbNodes = nbSharedNodes;
7123 freeFaceList.erase( fIt++ ); // here fIt++ occures before erase
7125 if ( freeFaceList.size() > 1 )
7127 // could not choose one face, use another way
7128 // choose a face most close to the bary center of the opposite side
7129 gp_XYZ aBC( 0., 0., 0. );
7130 set <const SMDS_MeshNode*> addedNodes;
7131 TIDSortedElemSet * elemSet2 = elemSetPtr[ 1 - iSide ];
7132 eIt = elemSet2->begin();
7133 for ( eIt = elemSet2->begin(); eIt != elemSet2->end(); eIt++ ) {
7134 SMDS_ElemIteratorPtr nodeIt = (*eIt)->nodesIterator();
7135 while ( nodeIt->more() ) { // loop on free face nodes
7136 const SMDS_MeshNode* n =
7137 static_cast<const SMDS_MeshNode*>( nodeIt->next() );
7138 if ( addedNodes.insert( n ).second )
7139 aBC += gp_XYZ( n->X(),n->Y(),n->Z() );
7142 aBC /= addedNodes.size();
7143 double minDist = DBL_MAX;
7144 fIt = freeFaceList.begin();
7145 while ( fIt != freeFaceList.end() ) { // loop on free faces
7147 SMDS_ElemIteratorPtr nodeIt = (*fIt)->nodesIterator();
7148 while ( nodeIt->more() ) { // loop on free face nodes
7149 const SMDS_MeshNode* n =
7150 static_cast<const SMDS_MeshNode*>( nodeIt->next() );
7151 gp_XYZ p( n->X(),n->Y(),n->Z() );
7152 dist += ( aBC - p ).SquareModulus();
7154 if ( dist < minDist ) {
7156 freeFaceList.erase( freeFaceList.begin(), fIt++ );
7159 fIt = freeFaceList.erase( fIt++ );
7162 } // choose one of several free faces of a volume
7164 if ( freeFaceList.size() == 1 ) {
7165 const SMDS_MeshElement* aFreeFace = freeFaceList.front();
7166 faceSet->insert( aFreeFace );
7167 // complete a node set with nodes of a found free face
7168 // for ( iNode = 0; iNode < ; iNode++ )
7169 // nodeSet->insert( fNodes[ iNode ] );
7172 } // loop on volumes of a side
7174 // // complete a set of faces if new nodes in a nodeSet appeared
7175 // // ----------------------------------------------------------
7176 // if ( nodeSetSize != nodeSet->size() ) {
7177 // for ( ; nIt != nodeSet->end(); nIt++ ) { // loop on nodes of iSide
7178 // SMDS_ElemIteratorPtr fIt = (*nIt)->GetInverseElementIterator(SMDSAbs_Face);
7179 // while ( fIt->more() ) { // loop on faces sharing a node
7180 // const SMDS_MeshElement* f = fIt->next();
7181 // if ( faceSet->find( f ) == faceSet->end() ) {
7182 // // check if all nodes are in nodeSet and
7183 // // complete setOfFaceNodeSet if they are
7184 // set <const SMDS_MeshNode*> faceNodeSet;
7185 // SMDS_ElemIteratorPtr nodeIt = f->nodesIterator();
7186 // bool allInSet = true;
7187 // while ( nodeIt->more() && allInSet ) { // loop on nodes of a face
7188 // const SMDS_MeshNode* n = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
7189 // if ( nodeSet->find( n ) == nodeSet->end() )
7190 // allInSet = false;
7192 // faceNodeSet.insert( n );
7194 // if ( allInSet ) {
7195 // faceSet->insert( f );
7196 // setOfFaceNodeSet.insert( faceNodeSet );
7202 } // Create temporary faces, if there are volumes given
7205 if ( faceSet1.size() != faceSet2.size() ) {
7206 // delete temporary faces: they are in reverseElements of actual nodes
7207 SMDS_FaceIteratorPtr tmpFaceIt = aTmpFacesMesh.facesIterator();
7208 while ( tmpFaceIt->more() )
7209 aTmpFacesMesh.RemoveElement( tmpFaceIt->next() );
7210 MESSAGE("Diff nb of faces");
7211 return SEW_TOPO_DIFF_SETS_OF_ELEMENTS;
7214 // ============================================================
7215 // 2. Find nodes to merge:
7216 // bind a node to remove to a node to put instead
7217 // ============================================================
7219 TNodeNodeMap nReplaceMap; // bind a node to remove to a node to put instead
7220 if ( theFirstNode1 != theFirstNode2 )
7221 nReplaceMap.insert( TNodeNodeMap::value_type( theFirstNode1, theFirstNode2 ));
7222 if ( theSecondNode1 != theSecondNode2 )
7223 nReplaceMap.insert( TNodeNodeMap::value_type( theSecondNode1, theSecondNode2 ));
7225 LinkID_Gen aLinkID_Gen( GetMeshDS() );
7226 set< long > linkIdSet; // links to process
7227 linkIdSet.insert( aLinkID_Gen.GetLinkID( theFirstNode1, theSecondNode1 ));
7229 typedef pair< const SMDS_MeshNode*, const SMDS_MeshNode* > NLink;
7230 list< NLink > linkList[2];
7231 linkList[0].push_back( NLink( theFirstNode1, theSecondNode1 ));
7232 linkList[1].push_back( NLink( theFirstNode2, theSecondNode2 ));
7233 // loop on links in linkList; find faces by links and append links
7234 // of the found faces to linkList
7235 list< NLink >::iterator linkIt[] = { linkList[0].begin(), linkList[1].begin() } ;
7236 for ( ; linkIt[0] != linkList[0].end(); linkIt[0]++, linkIt[1]++ ) {
7237 NLink link[] = { *linkIt[0], *linkIt[1] };
7238 long linkID = aLinkID_Gen.GetLinkID( link[0].first, link[0].second );
7239 if ( linkIdSet.find( linkID ) == linkIdSet.end() )
7242 // by links, find faces in the face sets,
7243 // and find indices of link nodes in the found faces;
7244 // in a face set, there is only one or no face sharing a link
7245 // ---------------------------------------------------------------
7247 const SMDS_MeshElement* face[] = { 0, 0 };
7248 //const SMDS_MeshNode* faceNodes[ 2 ][ 5 ];
7249 vector<const SMDS_MeshNode*> fnodes1(9);
7250 vector<const SMDS_MeshNode*> fnodes2(9);
7251 //const SMDS_MeshNode* notLinkNodes[ 2 ][ 2 ] = {{ 0, 0 },{ 0, 0 }} ;
7252 vector<const SMDS_MeshNode*> notLinkNodes1(6);
7253 vector<const SMDS_MeshNode*> notLinkNodes2(6);
7254 int iLinkNode[2][2];
7255 for ( iSide = 0; iSide < 2; iSide++ ) { // loop on 2 sides
7256 const SMDS_MeshNode* n1 = link[iSide].first;
7257 const SMDS_MeshNode* n2 = link[iSide].second;
7258 set<const SMDS_MeshElement*> * faceSet = faceSetPtr[ iSide ];
7259 set< const SMDS_MeshElement* > fMap;
7260 for ( int i = 0; i < 2; i++ ) { // loop on 2 nodes of a link
7261 const SMDS_MeshNode* n = i ? n1 : n2; // a node of a link
7262 SMDS_ElemIteratorPtr fIt = n->GetInverseElementIterator(SMDSAbs_Face);
7263 while ( fIt->more() ) { // loop on faces sharing a node
7264 const SMDS_MeshElement* f = fIt->next();
7265 if (faceSet->find( f ) != faceSet->end() && // f is in face set
7266 ! fMap.insert( f ).second ) // f encounters twice
7268 if ( face[ iSide ] ) {
7269 MESSAGE( "2 faces per link " );
7270 aResult = iSide ? SEW_BAD_SIDE2_NODES : SEW_BAD_SIDE1_NODES;
7274 faceSet->erase( f );
7275 // get face nodes and find ones of a link
7280 fnodes1.resize(f->NbNodes()+1);
7281 notLinkNodes1.resize(f->NbNodes()-2);
7284 fnodes2.resize(f->NbNodes()+1);
7285 notLinkNodes2.resize(f->NbNodes()-2);
7288 if(!f->IsQuadratic()) {
7289 SMDS_ElemIteratorPtr nIt = f->nodesIterator();
7290 while ( nIt->more() ) {
7291 const SMDS_MeshNode* n =
7292 static_cast<const SMDS_MeshNode*>( nIt->next() );
7294 iLinkNode[ iSide ][ 0 ] = iNode;
7296 else if ( n == n2 ) {
7297 iLinkNode[ iSide ][ 1 ] = iNode;
7299 //else if ( notLinkNodes[ iSide ][ 0 ] )
7300 // notLinkNodes[ iSide ][ 1 ] = n;
7302 // notLinkNodes[ iSide ][ 0 ] = n;
7306 notLinkNodes1[nbl] = n;
7307 //notLinkNodes1.push_back(n);
7309 notLinkNodes2[nbl] = n;
7310 //notLinkNodes2.push_back(n);
7312 //faceNodes[ iSide ][ iNode++ ] = n;
7314 fnodes1[iNode++] = n;
7317 fnodes2[iNode++] = n;
7321 else { // f->IsQuadratic()
7322 const SMDS_QuadraticFaceOfNodes* F =
7323 static_cast<const SMDS_QuadraticFaceOfNodes*>(f);
7324 // use special nodes iterator
7325 SMDS_NodeIteratorPtr anIter = F->interlacedNodesIterator();
7326 while ( anIter->more() ) {
7327 const SMDS_MeshNode* n =
7328 static_cast<const SMDS_MeshNode*>( anIter->next() );
7330 iLinkNode[ iSide ][ 0 ] = iNode;
7332 else if ( n == n2 ) {
7333 iLinkNode[ iSide ][ 1 ] = iNode;
7338 notLinkNodes1[nbl] = n;
7341 notLinkNodes2[nbl] = n;
7345 fnodes1[iNode++] = n;
7348 fnodes2[iNode++] = n;
7352 //faceNodes[ iSide ][ iNode ] = faceNodes[ iSide ][ 0 ];
7354 fnodes1[iNode] = fnodes1[0];
7357 fnodes2[iNode] = fnodes1[0];
7364 // check similarity of elements of the sides
7365 if (aResult == SEW_OK && ( face[0] && !face[1] ) || ( !face[0] && face[1] )) {
7366 MESSAGE("Correspondent face not found on side " << ( face[0] ? 1 : 0 ));
7367 if ( nReplaceMap.size() == 2 ) { // faces on input nodes not found
7368 aResult = ( face[0] ? SEW_BAD_SIDE2_NODES : SEW_BAD_SIDE1_NODES );
7371 aResult = SEW_TOPO_DIFF_SETS_OF_ELEMENTS;
7373 break; // do not return because it s necessary to remove tmp faces
7376 // set nodes to merge
7377 // -------------------
7379 if ( face[0] && face[1] ) {
7380 int nbNodes = face[0]->NbNodes();
7381 if ( nbNodes != face[1]->NbNodes() ) {
7382 MESSAGE("Diff nb of face nodes");
7383 aResult = SEW_TOPO_DIFF_SETS_OF_ELEMENTS;
7384 break; // do not return because it s necessary to remove tmp faces
7386 bool reverse[] = { false, false }; // order of notLinkNodes of quadrangle
7387 if ( nbNodes == 3 ) {
7388 //nReplaceMap.insert( TNodeNodeMap::value_type
7389 // ( notLinkNodes[0][0], notLinkNodes[1][0] ));
7390 nReplaceMap.insert( TNodeNodeMap::value_type
7391 ( notLinkNodes1[0], notLinkNodes2[0] ));
7394 for ( iSide = 0; iSide < 2; iSide++ ) { // loop on 2 sides
7395 // analyse link orientation in faces
7396 int i1 = iLinkNode[ iSide ][ 0 ];
7397 int i2 = iLinkNode[ iSide ][ 1 ];
7398 reverse[ iSide ] = Abs( i1 - i2 ) == 1 ? i1 > i2 : i2 > i1;
7399 // if notLinkNodes are the first and the last ones, then
7400 // their order does not correspond to the link orientation
7401 if (( i1 == 1 && i2 == 2 ) ||
7402 ( i1 == 2 && i2 == 1 ))
7403 reverse[ iSide ] = !reverse[ iSide ];
7405 if ( reverse[0] == reverse[1] ) {
7406 //nReplaceMap.insert( TNodeNodeMap::value_type
7407 // ( notLinkNodes[0][0], notLinkNodes[1][0] ));
7408 //nReplaceMap.insert( TNodeNodeMap::value_type
7409 // ( notLinkNodes[0][1], notLinkNodes[1][1] ));
7410 for(int nn=0; nn<nbNodes-2; nn++) {
7411 nReplaceMap.insert( TNodeNodeMap::value_type
7412 ( notLinkNodes1[nn], notLinkNodes2[nn] ));
7416 //nReplaceMap.insert( TNodeNodeMap::value_type
7417 // ( notLinkNodes[0][0], notLinkNodes[1][1] ));
7418 //nReplaceMap.insert( TNodeNodeMap::value_type
7419 // ( notLinkNodes[0][1], notLinkNodes[1][0] ));
7420 for(int nn=0; nn<nbNodes-2; nn++) {
7421 nReplaceMap.insert( TNodeNodeMap::value_type
7422 ( notLinkNodes1[nn], notLinkNodes2[nbNodes-3-nn] ));
7427 // add other links of the faces to linkList
7428 // -----------------------------------------
7430 //const SMDS_MeshNode** nodes = faceNodes[ 0 ];
7431 for ( iNode = 0; iNode < nbNodes; iNode++ ) {
7432 //linkID = aLinkID_Gen.GetLinkID( nodes[iNode], nodes[iNode+1] );
7433 linkID = aLinkID_Gen.GetLinkID( fnodes1[iNode], fnodes1[iNode+1] );
7434 pair< set<long>::iterator, bool > iter_isnew = linkIdSet.insert( linkID );
7435 if ( !iter_isnew.second ) { // already in a set: no need to process
7436 linkIdSet.erase( iter_isnew.first );
7438 else // new in set == encountered for the first time: add
7440 //const SMDS_MeshNode* n1 = nodes[ iNode ];
7441 //const SMDS_MeshNode* n2 = nodes[ iNode + 1];
7442 const SMDS_MeshNode* n1 = fnodes1[ iNode ];
7443 const SMDS_MeshNode* n2 = fnodes1[ iNode + 1];
7444 linkList[0].push_back ( NLink( n1, n2 ));
7445 linkList[1].push_back ( NLink( nReplaceMap[n1], nReplaceMap[n2] ));
7449 } // loop on link lists
7451 if ( aResult == SEW_OK &&
7452 ( linkIt[0] != linkList[0].end() ||
7453 !faceSetPtr[0]->empty() || !faceSetPtr[1]->empty() )) {
7454 MESSAGE( (linkIt[0] != linkList[0].end()) <<" "<< (faceSetPtr[0]->empty()) <<
7455 " " << (faceSetPtr[1]->empty()));
7456 aResult = SEW_TOPO_DIFF_SETS_OF_ELEMENTS;
7459 // ====================================================================
7460 // 3. Replace nodes in elements of the side 1 and remove replaced nodes
7461 // ====================================================================
7463 // delete temporary faces: they are in reverseElements of actual nodes
7464 SMDS_FaceIteratorPtr tmpFaceIt = aTmpFacesMesh.facesIterator();
7465 while ( tmpFaceIt->more() )
7466 aTmpFacesMesh.RemoveElement( tmpFaceIt->next() );
7468 if ( aResult != SEW_OK)
7471 list< int > nodeIDsToRemove/*, elemIDsToRemove*/;
7472 // loop on nodes replacement map
7473 TNodeNodeMap::iterator nReplaceMapIt = nReplaceMap.begin(), nnIt;
7474 for ( ; nReplaceMapIt != nReplaceMap.end(); nReplaceMapIt++ )
7475 if ( (*nReplaceMapIt).first != (*nReplaceMapIt).second ) {
7476 const SMDS_MeshNode* nToRemove = (*nReplaceMapIt).first;
7477 nodeIDsToRemove.push_back( nToRemove->GetID() );
7478 // loop on elements sharing nToRemove
7479 SMDS_ElemIteratorPtr invElemIt = nToRemove->GetInverseElementIterator();
7480 while ( invElemIt->more() ) {
7481 const SMDS_MeshElement* e = invElemIt->next();
7482 // get a new suite of nodes: make replacement
7483 int nbReplaced = 0, i = 0, nbNodes = e->NbNodes();
7484 vector< const SMDS_MeshNode*> nodes( nbNodes );
7485 SMDS_ElemIteratorPtr nIt = e->nodesIterator();
7486 while ( nIt->more() ) {
7487 const SMDS_MeshNode* n =
7488 static_cast<const SMDS_MeshNode*>( nIt->next() );
7489 nnIt = nReplaceMap.find( n );
7490 if ( nnIt != nReplaceMap.end() ) {
7496 // if ( nbReplaced == nbNodes && e->GetType() == SMDSAbs_Face )
7497 // elemIDsToRemove.push_back( e->GetID() );
7500 aMesh->ChangeElementNodes( e, & nodes[0], nbNodes );
7504 Remove( nodeIDsToRemove, true );
7509 //================================================================================
7511 * \brief Find corresponding nodes in two sets of faces
7512 * \param theSide1 - first face set
7513 * \param theSide2 - second first face
7514 * \param theFirstNode1 - a boundary node of set 1
7515 * \param theFirstNode2 - a node of set 2 corresponding to theFirstNode1
7516 * \param theSecondNode1 - a boundary node of set 1 linked with theFirstNode1
7517 * \param theSecondNode2 - a node of set 2 corresponding to theSecondNode1
7518 * \param nReplaceMap - output map of corresponding nodes
7519 * \retval bool - is a success or not
7521 //================================================================================
7524 //#define DEBUG_MATCHING_NODES
7527 SMESH_MeshEditor::Sew_Error
7528 SMESH_MeshEditor::FindMatchingNodes(set<const SMDS_MeshElement*>& theSide1,
7529 set<const SMDS_MeshElement*>& theSide2,
7530 const SMDS_MeshNode* theFirstNode1,
7531 const SMDS_MeshNode* theFirstNode2,
7532 const SMDS_MeshNode* theSecondNode1,
7533 const SMDS_MeshNode* theSecondNode2,
7534 TNodeNodeMap & nReplaceMap)
7536 set<const SMDS_MeshElement*> * faceSetPtr[] = { &theSide1, &theSide2 };
7538 nReplaceMap.clear();
7539 if ( theFirstNode1 != theFirstNode2 )
7540 nReplaceMap.insert( make_pair( theFirstNode1, theFirstNode2 ));
7541 if ( theSecondNode1 != theSecondNode2 )
7542 nReplaceMap.insert( make_pair( theSecondNode1, theSecondNode2 ));
7544 set< TLink > linkSet; // set of nodes where order of nodes is ignored
7545 linkSet.insert( TLink( theFirstNode1, theSecondNode1 ));
7547 list< NLink > linkList[2];
7548 linkList[0].push_back( NLink( theFirstNode1, theSecondNode1 ));
7549 linkList[1].push_back( NLink( theFirstNode2, theSecondNode2 ));
7551 // loop on links in linkList; find faces by links and append links
7552 // of the found faces to linkList
7553 list< NLink >::iterator linkIt[] = { linkList[0].begin(), linkList[1].begin() } ;
7554 for ( ; linkIt[0] != linkList[0].end(); linkIt[0]++, linkIt[1]++ ) {
7555 NLink link[] = { *linkIt[0], *linkIt[1] };
7556 if ( linkSet.find( link[0] ) == linkSet.end() )
7559 // by links, find faces in the face sets,
7560 // and find indices of link nodes in the found faces;
7561 // in a face set, there is only one or no face sharing a link
7562 // ---------------------------------------------------------------
7564 const SMDS_MeshElement* face[] = { 0, 0 };
7565 list<const SMDS_MeshNode*> notLinkNodes[2];
7566 //bool reverse[] = { false, false }; // order of notLinkNodes
7568 for ( int iSide = 0; iSide < 2; iSide++ ) // loop on 2 sides
7570 const SMDS_MeshNode* n1 = link[iSide].first;
7571 const SMDS_MeshNode* n2 = link[iSide].second;
7572 set<const SMDS_MeshElement*> * faceSet = faceSetPtr[ iSide ];
7573 set< const SMDS_MeshElement* > facesOfNode1;
7574 for ( int iNode = 0; iNode < 2; iNode++ ) // loop on 2 nodes of a link
7576 // during a loop of the first node, we find all faces around n1,
7577 // during a loop of the second node, we find one face sharing both n1 and n2
7578 const SMDS_MeshNode* n = iNode ? n1 : n2; // a node of a link
7579 SMDS_ElemIteratorPtr fIt = n->GetInverseElementIterator(SMDSAbs_Face);
7580 while ( fIt->more() ) { // loop on faces sharing a node
7581 const SMDS_MeshElement* f = fIt->next();
7582 if (faceSet->find( f ) != faceSet->end() && // f is in face set
7583 ! facesOfNode1.insert( f ).second ) // f encounters twice
7585 if ( face[ iSide ] ) {
7586 MESSAGE( "2 faces per link " );
7587 return ( iSide ? SEW_BAD_SIDE2_NODES : SEW_BAD_SIDE1_NODES );
7590 faceSet->erase( f );
7592 // get not link nodes
7593 int nbN = f->NbNodes();
7594 if ( f->IsQuadratic() )
7596 nbNodes[ iSide ] = nbN;
7597 list< const SMDS_MeshNode* > & nodes = notLinkNodes[ iSide ];
7598 int i1 = f->GetNodeIndex( n1 );
7599 int i2 = f->GetNodeIndex( n2 );
7600 int iEnd = nbN, iBeg = -1, iDelta = 1;
7601 bool reverse = ( Abs( i1 - i2 ) == 1 ? i1 > i2 : i2 > i1 );
7603 std::swap( iEnd, iBeg ); iDelta = -1;
7608 if ( i == iEnd ) i = iBeg + iDelta;
7609 if ( i == i1 ) break;
7610 nodes.push_back ( f->GetNode( i ) );
7616 // check similarity of elements of the sides
7617 if (( face[0] && !face[1] ) || ( !face[0] && face[1] )) {
7618 MESSAGE("Correspondent face not found on side " << ( face[0] ? 1 : 0 ));
7619 if ( nReplaceMap.size() == 2 ) { // faces on input nodes not found
7620 return ( face[0] ? SEW_BAD_SIDE2_NODES : SEW_BAD_SIDE1_NODES );
7623 return SEW_TOPO_DIFF_SETS_OF_ELEMENTS;
7627 // set nodes to merge
7628 // -------------------
7630 if ( face[0] && face[1] ) {
7631 if ( nbNodes[0] != nbNodes[1] ) {
7632 MESSAGE("Diff nb of face nodes");
7633 return SEW_TOPO_DIFF_SETS_OF_ELEMENTS;
7635 #ifdef DEBUG_MATCHING_NODES
7636 MESSAGE ( " Link 1: " << link[0].first->GetID() <<" "<< link[0].second->GetID()
7637 << " F 1: " << face[0] << "| Link 2: " << link[1].first->GetID() <<" "
7638 << link[1].second->GetID() << " F 2: " << face[1] << " | Bind: " ) ;
7640 int nbN = nbNodes[0];
7642 list<const SMDS_MeshNode*>::iterator n1 = notLinkNodes[0].begin();
7643 list<const SMDS_MeshNode*>::iterator n2 = notLinkNodes[1].begin();
7644 for ( int i = 0 ; i < nbN - 2; ++i ) {
7645 #ifdef DEBUG_MATCHING_NODES
7646 MESSAGE ( (*n1)->GetID() << " to " << (*n2)->GetID() );
7648 nReplaceMap.insert( make_pair( *(n1++), *(n2++) ));
7652 // add other links of the face 1 to linkList
7653 // -----------------------------------------
7655 const SMDS_MeshElement* f0 = face[0];
7656 const SMDS_MeshNode* n1 = f0->GetNode( nbN - 1 );
7657 for ( int i = 0; i < nbN; i++ )
7659 const SMDS_MeshNode* n2 = f0->GetNode( i );
7660 pair< set< TLink >::iterator, bool > iter_isnew =
7661 linkSet.insert( TLink( n1, n2 ));
7662 if ( !iter_isnew.second ) { // already in a set: no need to process
7663 linkSet.erase( iter_isnew.first );
7665 else // new in set == encountered for the first time: add
7667 #ifdef DEBUG_MATCHING_NODES
7668 MESSAGE ( "Add link 1: " << n1->GetID() << " " << n2->GetID() << " "
7669 << " | link 2: " << nReplaceMap[n1]->GetID() << " " << nReplaceMap[n2]->GetID() << " " );
7671 linkList[0].push_back ( NLink( n1, n2 ));
7672 linkList[1].push_back ( NLink( nReplaceMap[n1], nReplaceMap[n2] ));
7677 } // loop on link lists
7683 \brief Creates a hole in a mesh by doubling the nodes of some particular elements
7684 \param theNodes - identifiers of nodes to be doubled
7685 \param theModifiedElems - identifiers of elements to be updated by the new (doubled)
7686 nodes. If list of element identifiers is empty then nodes are doubled but
7687 they not assigned to elements
7688 \return TRUE if operation has been completed successfully, FALSE otherwise
7690 bool SMESH_MeshEditor::DoubleNodes( const std::list< int >& theListOfNodes,
7691 const std::list< int >& theListOfModifiedElems )
7693 myLastCreatedElems.Clear();
7694 myLastCreatedNodes.Clear();
7696 if ( theListOfNodes.size() == 0 )
7699 SMESHDS_Mesh* aMeshDS = GetMeshDS();
7703 // iterate through nodes and duplicate them
7705 std::map< const SMDS_MeshNode*, const SMDS_MeshNode* > anOldNodeToNewNode;
7707 std::list< int >::const_iterator aNodeIter;
7708 for ( aNodeIter = theListOfNodes.begin(); aNodeIter != theListOfNodes.end(); ++aNodeIter )
7710 int aCurr = *aNodeIter;
7711 SMDS_MeshNode* aNode = (SMDS_MeshNode*)aMeshDS->FindNode( aCurr );
7717 const SMDS_MeshNode* aNewNode = aMeshDS->AddNode( aNode->X(), aNode->Y(), aNode->Z() );
7720 anOldNodeToNewNode[ aNode ] = aNewNode;
7721 myLastCreatedNodes.Append( aNewNode );
7725 // Create map of new nodes for modified elements
7727 std::map< SMDS_MeshElement*, vector<const SMDS_MeshNode*> > anElemToNodes;
7729 std::list< int >::const_iterator anElemIter;
7730 for ( anElemIter = theListOfModifiedElems.begin();
7731 anElemIter != theListOfModifiedElems.end(); ++anElemIter )
7733 int aCurr = *anElemIter;
7734 SMDS_MeshElement* anElem = (SMDS_MeshElement*)aMeshDS->FindElement( aCurr );
7738 vector<const SMDS_MeshNode*> aNodeArr( anElem->NbNodes() );
7740 SMDS_ElemIteratorPtr anIter = anElem->nodesIterator();
7742 while ( anIter->more() )
7744 SMDS_MeshNode* aCurrNode = (SMDS_MeshNode*)anIter->next();
7745 if ( aCurr && anOldNodeToNewNode.find( aCurrNode ) != anOldNodeToNewNode.end() )
7747 const SMDS_MeshNode* aNewNode = anOldNodeToNewNode[ aCurrNode ];
7748 aNodeArr[ ind++ ] = aNewNode;
7751 aNodeArr[ ind++ ] = aCurrNode;
7753 anElemToNodes[ anElem ] = aNodeArr;
7756 // Change nodes of elements
7758 std::map< SMDS_MeshElement*, vector<const SMDS_MeshNode*> >::iterator
7759 anElemToNodesIter = anElemToNodes.begin();
7760 for ( ; anElemToNodesIter != anElemToNodes.end(); ++anElemToNodesIter )
7762 const SMDS_MeshElement* anElem = anElemToNodesIter->first;
7763 vector<const SMDS_MeshNode*> aNodeArr = anElemToNodesIter->second;
7765 aMeshDS->ChangeElementNodes( anElem, &aNodeArr[ 0 ], anElem->NbNodes() );