1 // Copyright (C) 2007-2008 CEA/DEN, EDF R&D, OPEN CASCADE
3 // Copyright (C) 2003-2007 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
4 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
6 // This library is free software; you can redistribute it and/or
7 // modify it under the terms of the GNU Lesser General Public
8 // License as published by the Free Software Foundation; either
9 // version 2.1 of the License.
11 // This library is distributed in the hope that it will be useful,
12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14 // Lesser General Public License for more details.
16 // You should have received a copy of the GNU Lesser General Public
17 // License along with this library; if not, write to the Free Software
18 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
20 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
22 // SMESH SMESH : idl implementation based on 'SMESH' unit's classes
23 // File : SMESH_MeshEditor.cxx
24 // Created : Mon Apr 12 16:10:22 2004
25 // Author : Edward AGAPOV (eap)
27 #include "SMESH_MeshEditor.hxx"
29 #include "SMDS_FaceOfNodes.hxx"
30 #include "SMDS_VolumeTool.hxx"
31 #include "SMDS_EdgePosition.hxx"
32 #include "SMDS_PolyhedralVolumeOfNodes.hxx"
33 #include "SMDS_FacePosition.hxx"
34 #include "SMDS_SpacePosition.hxx"
35 #include "SMDS_QuadraticFaceOfNodes.hxx"
36 #include "SMDS_MeshGroup.hxx"
38 #include "SMESHDS_Group.hxx"
39 #include "SMESHDS_Mesh.hxx"
41 #include "SMESH_subMesh.hxx"
42 #include "SMESH_ControlsDef.hxx"
43 #include "SMESH_MesherHelper.hxx"
44 #include "SMESH_OctreeNode.hxx"
45 #include "SMESH_Group.hxx"
47 #include "utilities.h"
49 #include <BRep_Tool.hxx>
51 #include <Extrema_GenExtPS.hxx>
52 #include <Extrema_POnSurf.hxx>
53 #include <Geom2d_Curve.hxx>
54 #include <GeomAdaptor_Surface.hxx>
55 #include <Geom_Curve.hxx>
56 #include <Geom_Surface.hxx>
57 #include <TColStd_ListOfInteger.hxx>
59 #include <TopExp_Explorer.hxx>
60 #include <TopTools_ListIteratorOfListOfShape.hxx>
61 #include <TopTools_ListOfShape.hxx>
63 #include <TopoDS_Face.hxx>
69 #include <gp_Trsf.hxx>
78 #define cast2Node(elem) static_cast<const SMDS_MeshNode*>( elem )
81 using namespace SMESH::Controls;
83 typedef map<const SMDS_MeshElement*, list<const SMDS_MeshNode*> > TElemOfNodeListMap;
84 typedef map<const SMDS_MeshElement*, list<const SMDS_MeshElement*> > TElemOfElemListMap;
85 //typedef map<const SMDS_MeshNode*, vector<const SMDS_MeshNode*> > TNodeOfNodeVecMap;
86 //typedef TNodeOfNodeVecMap::iterator TNodeOfNodeVecMapItr;
87 //typedef map<const SMDS_MeshElement*, vector<TNodeOfNodeVecMapItr> > TElemOfVecOfMapNodesMap;
89 struct TNodeXYZ : public gp_XYZ {
90 TNodeXYZ( const SMDS_MeshNode* n ):gp_XYZ( n->X(), n->Y(), n->Z() ) {}
93 //=======================================================================
94 //function : SMESH_MeshEditor
96 //=======================================================================
98 SMESH_MeshEditor::SMESH_MeshEditor( SMESH_Mesh* theMesh )
99 :myMesh( theMesh ) // theMesh may be NULL
103 //=======================================================================
107 //=======================================================================
110 SMESH_MeshEditor::AddElement(const vector<const SMDS_MeshNode*> & node,
111 const SMDSAbs_ElementType type,
115 SMDS_MeshElement* e = 0;
116 int nbnode = node.size();
117 SMESHDS_Mesh* mesh = GetMeshDS();
121 if ( ID ) e = mesh->AddEdgeWithID(node[0], node[1], ID);
122 else e = mesh->AddEdge (node[0], node[1] );
123 else if ( nbnode == 3 )
124 if ( ID ) e = mesh->AddEdgeWithID(node[0], node[1], node[2], ID);
125 else e = mesh->AddEdge (node[0], node[1], node[2] );
130 if ( ID ) e = mesh->AddFaceWithID(node[0], node[1], node[2], ID);
131 else e = mesh->AddFace (node[0], node[1], node[2] );
132 else if (nbnode == 4)
133 if ( ID ) e = mesh->AddFaceWithID(node[0], node[1], node[2], node[3], ID);
134 else e = mesh->AddFace (node[0], node[1], node[2], node[3] );
135 else if (nbnode == 6)
136 if ( ID ) e = mesh->AddFaceWithID(node[0], node[1], node[2], node[3],
137 node[4], node[5], ID);
138 else e = mesh->AddFace (node[0], node[1], node[2], node[3],
140 else if (nbnode == 8)
141 if ( ID ) e = mesh->AddFaceWithID(node[0], node[1], node[2], node[3],
142 node[4], node[5], node[6], node[7], ID);
143 else e = mesh->AddFace (node[0], node[1], node[2], node[3],
144 node[4], node[5], node[6], node[7] );
146 if ( ID ) e = mesh->AddPolygonalFaceWithID(node, ID);
147 else e = mesh->AddPolygonalFace (node );
153 if ( ID ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3], ID);
154 else e = mesh->AddVolume (node[0], node[1], node[2], node[3] );
155 else if (nbnode == 5)
156 if ( ID ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3],
158 else e = mesh->AddVolume (node[0], node[1], node[2], node[3],
160 else if (nbnode == 6)
161 if ( ID ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3],
162 node[4], node[5], ID);
163 else e = mesh->AddVolume (node[0], node[1], node[2], node[3],
165 else if (nbnode == 8)
166 if ( ID ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3],
167 node[4], node[5], node[6], node[7], ID);
168 else e = mesh->AddVolume (node[0], node[1], node[2], node[3],
169 node[4], node[5], node[6], node[7] );
170 else if (nbnode == 10)
171 if ( ID ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3],
172 node[4], node[5], node[6], node[7],
173 node[8], node[9], ID);
174 else e = mesh->AddVolume (node[0], node[1], node[2], node[3],
175 node[4], node[5], node[6], node[7],
177 else if (nbnode == 13)
178 if ( ID ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3],
179 node[4], node[5], node[6], node[7],
180 node[8], node[9], node[10],node[11],
182 else e = mesh->AddVolume (node[0], node[1], node[2], node[3],
183 node[4], node[5], node[6], node[7],
184 node[8], node[9], node[10],node[11],
186 else if (nbnode == 15)
187 if ( ID ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3],
188 node[4], node[5], node[6], node[7],
189 node[8], node[9], node[10],node[11],
190 node[12],node[13],node[14],ID);
191 else e = mesh->AddVolume (node[0], node[1], node[2], node[3],
192 node[4], node[5], node[6], node[7],
193 node[8], node[9], node[10],node[11],
194 node[12],node[13],node[14] );
195 else if (nbnode == 20)
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],
199 node[12],node[13],node[14],node[15],
200 node[16],node[17],node[18],node[19],ID);
201 else e = mesh->AddVolume (node[0], node[1], node[2], node[3],
202 node[4], node[5], node[6], node[7],
203 node[8], node[9], node[10],node[11],
204 node[12],node[13],node[14],node[15],
205 node[16],node[17],node[18],node[19] );
211 //=======================================================================
215 //=======================================================================
217 SMDS_MeshElement* SMESH_MeshEditor::AddElement(const vector<int> & nodeIDs,
218 const SMDSAbs_ElementType type,
222 vector<const SMDS_MeshNode*> nodes;
223 nodes.reserve( nodeIDs.size() );
224 vector<int>::const_iterator id = nodeIDs.begin();
225 while ( id != nodeIDs.end() ) {
226 if ( const SMDS_MeshNode* node = GetMeshDS()->FindNode( *id++ ))
227 nodes.push_back( node );
231 return AddElement( nodes, type, isPoly, ID );
234 //=======================================================================
236 //purpose : Remove a node or an element.
237 // Modify a compute state of sub-meshes which become empty
238 //=======================================================================
240 bool SMESH_MeshEditor::Remove (const list< int >& theIDs,
243 myLastCreatedElems.Clear();
244 myLastCreatedNodes.Clear();
246 SMESHDS_Mesh* aMesh = GetMeshDS();
247 set< SMESH_subMesh *> smmap;
249 list<int>::const_iterator it = theIDs.begin();
250 for ( ; it != theIDs.end(); it++ ) {
251 const SMDS_MeshElement * elem;
253 elem = aMesh->FindNode( *it );
255 elem = aMesh->FindElement( *it );
259 // Notify VERTEX sub-meshes about modification
261 const SMDS_MeshNode* node = cast2Node( elem );
262 if ( node->GetPosition()->GetTypeOfPosition() == SMDS_TOP_VERTEX )
263 if ( int aShapeID = node->GetPosition()->GetShapeId() )
264 if ( SMESH_subMesh * sm = GetMesh()->GetSubMeshContaining( aShapeID ) )
267 // Find sub-meshes to notify about modification
268 // SMDS_ElemIteratorPtr nodeIt = elem->nodesIterator();
269 // while ( nodeIt->more() ) {
270 // const SMDS_MeshNode* node = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
271 // const SMDS_PositionPtr& aPosition = node->GetPosition();
272 // if ( aPosition.get() ) {
273 // if ( int aShapeID = aPosition->GetShapeId() ) {
274 // if ( SMESH_subMesh * sm = GetMesh()->GetSubMeshContaining( aShapeID ) )
275 // smmap.insert( sm );
282 aMesh->RemoveNode( static_cast< const SMDS_MeshNode* >( elem ));
284 aMesh->RemoveElement( elem );
287 // Notify sub-meshes about modification
288 if ( !smmap.empty() ) {
289 set< SMESH_subMesh *>::iterator smIt;
290 for ( smIt = smmap.begin(); smIt != smmap.end(); smIt++ )
291 (*smIt)->ComputeStateEngine( SMESH_subMesh::MESH_ENTITY_REMOVED );
294 // // Check if the whole mesh becomes empty
295 // if ( SMESH_subMesh * sm = GetMesh()->GetSubMeshContaining( 1 ) )
296 // sm->ComputeStateEngine( SMESH_subMesh::CHECK_COMPUTE_STATE );
301 //=======================================================================
302 //function : FindShape
303 //purpose : Return an index of the shape theElem is on
304 // or zero if a shape not found
305 //=======================================================================
307 int SMESH_MeshEditor::FindShape (const SMDS_MeshElement * theElem)
309 myLastCreatedElems.Clear();
310 myLastCreatedNodes.Clear();
312 SMESHDS_Mesh * aMesh = GetMeshDS();
313 if ( aMesh->ShapeToMesh().IsNull() )
316 if ( theElem->GetType() == SMDSAbs_Node ) {
317 const SMDS_PositionPtr& aPosition =
318 static_cast<const SMDS_MeshNode*>( theElem )->GetPosition();
319 if ( aPosition.get() )
320 return aPosition->GetShapeId();
325 TopoDS_Shape aShape; // the shape a node is on
326 SMDS_ElemIteratorPtr nodeIt = theElem->nodesIterator();
327 while ( nodeIt->more() ) {
328 const SMDS_MeshNode* node = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
329 const SMDS_PositionPtr& aPosition = node->GetPosition();
330 if ( aPosition.get() ) {
331 int aShapeID = aPosition->GetShapeId();
332 SMESHDS_SubMesh * sm = aMesh->MeshElements( aShapeID );
334 if ( sm->Contains( theElem ))
336 if ( aShape.IsNull() )
337 aShape = aMesh->IndexToShape( aShapeID );
340 //MESSAGE ( "::FindShape() No SubShape for aShapeID " << aShapeID );
345 // None of nodes is on a proper shape,
346 // find the shape among ancestors of aShape on which a node is
347 if ( aShape.IsNull() ) {
348 //MESSAGE ("::FindShape() - NONE node is on shape")
351 TopTools_ListIteratorOfListOfShape ancIt( GetMesh()->GetAncestors( aShape ));
352 for ( ; ancIt.More(); ancIt.Next() ) {
353 SMESHDS_SubMesh * sm = aMesh->MeshElements( ancIt.Value() );
354 if ( sm && sm->Contains( theElem ))
355 return aMesh->ShapeToIndex( ancIt.Value() );
358 //MESSAGE ("::FindShape() - SHAPE NOT FOUND")
362 //=======================================================================
363 //function : IsMedium
365 //=======================================================================
367 bool SMESH_MeshEditor::IsMedium(const SMDS_MeshNode* node,
368 const SMDSAbs_ElementType typeToCheck)
370 bool isMedium = false;
371 SMDS_ElemIteratorPtr it = node->GetInverseElementIterator(typeToCheck);
372 while (it->more() && !isMedium ) {
373 const SMDS_MeshElement* elem = it->next();
374 isMedium = elem->IsMediumNode(node);
379 //=======================================================================
380 //function : ShiftNodesQuadTria
382 // Shift nodes in the array corresponded to quadratic triangle
383 // example: (0,1,2,3,4,5) -> (1,2,0,4,5,3)
384 //=======================================================================
385 static void ShiftNodesQuadTria(const SMDS_MeshNode* aNodes[])
387 const SMDS_MeshNode* nd1 = aNodes[0];
388 aNodes[0] = aNodes[1];
389 aNodes[1] = aNodes[2];
391 const SMDS_MeshNode* nd2 = aNodes[3];
392 aNodes[3] = aNodes[4];
393 aNodes[4] = aNodes[5];
397 //=======================================================================
398 //function : GetNodesFromTwoTria
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 bool GetNodesFromTwoTria(const SMDS_MeshElement * theTria1,
404 const SMDS_MeshElement * theTria2,
405 const SMDS_MeshNode* N1[],
406 const SMDS_MeshNode* N2[])
408 SMDS_ElemIteratorPtr it = theTria1->nodesIterator();
411 N1[i] = static_cast<const SMDS_MeshNode*>( it->next() );
414 if(it->more()) return false;
415 it = theTria2->nodesIterator();
418 N2[i] = static_cast<const SMDS_MeshNode*>( it->next() );
421 if(it->more()) return false;
423 int sames[3] = {-1,-1,-1};
435 if(nbsames!=2) return false;
437 ShiftNodesQuadTria(N1);
439 ShiftNodesQuadTria(N1);
442 i = sames[0] + sames[1] + sames[2];
444 ShiftNodesQuadTria(N2);
446 // now we receive following N1 and N2 (using numeration as above image)
447 // tria1 : (1 2 4 5 9 7) and tria2 : (3 4 2 8 9 6)
448 // i.e. first nodes from both arrays determ new diagonal
452 //=======================================================================
453 //function : InverseDiag
454 //purpose : Replace two neighbour triangles with ones built on the same 4 nodes
455 // but having other common link.
456 // Return False if args are improper
457 //=======================================================================
459 bool SMESH_MeshEditor::InverseDiag (const SMDS_MeshElement * theTria1,
460 const SMDS_MeshElement * theTria2 )
462 myLastCreatedElems.Clear();
463 myLastCreatedNodes.Clear();
465 if (!theTria1 || !theTria2)
468 const SMDS_FaceOfNodes* F1 = dynamic_cast<const SMDS_FaceOfNodes*>( theTria1 );
469 const SMDS_FaceOfNodes* F2 = dynamic_cast<const SMDS_FaceOfNodes*>( theTria2 );
472 // 1 +--+ A theTria1: ( 1 A B ) A->2 ( 1 2 B ) 1 +--+ A
473 // | /| theTria2: ( B A 2 ) B->1 ( 1 A 2 ) |\ |
477 // put nodes in array and find out indices of the same ones
478 const SMDS_MeshNode* aNodes [6];
479 int sameInd [] = { 0, 0, 0, 0, 0, 0 };
481 SMDS_ElemIteratorPtr it = theTria1->nodesIterator();
482 while ( it->more() ) {
483 aNodes[ i ] = static_cast<const SMDS_MeshNode*>( it->next() );
485 if ( i > 2 ) // theTria2
486 // find same node of theTria1
487 for ( int j = 0; j < 3; j++ )
488 if ( aNodes[ i ] == aNodes[ j ]) {
497 return false; // theTria1 is not a triangle
498 it = theTria2->nodesIterator();
500 if ( i == 6 && it->more() )
501 return false; // theTria2 is not a triangle
504 // find indices of 1,2 and of A,B in theTria1
505 int iA = 0, iB = 0, i1 = 0, i2 = 0;
506 for ( i = 0; i < 6; i++ ) {
507 if ( sameInd [ i ] == 0 )
514 // nodes 1 and 2 should not be the same
515 if ( aNodes[ i1 ] == aNodes[ i2 ] )
519 aNodes[ iA ] = aNodes[ i2 ];
521 aNodes[ sameInd[ iB ]] = aNodes[ i1 ];
523 //MESSAGE( theTria1 << theTria2 );
525 GetMeshDS()->ChangeElementNodes( theTria1, aNodes, 3 );
526 GetMeshDS()->ChangeElementNodes( theTria2, &aNodes[ 3 ], 3 );
528 //MESSAGE( theTria1 << theTria2 );
532 } // end if(F1 && F2)
534 // check case of quadratic faces
535 const SMDS_QuadraticFaceOfNodes* QF1 =
536 dynamic_cast<const SMDS_QuadraticFaceOfNodes*> (theTria1);
537 if(!QF1) return false;
538 const SMDS_QuadraticFaceOfNodes* QF2 =
539 dynamic_cast<const SMDS_QuadraticFaceOfNodes*> (theTria2);
540 if(!QF2) return false;
543 // 1 +--+--+ 2 theTria1: (1 2 4 5 9 7) or (2 4 1 9 7 5) or (4 1 2 7 5 9)
544 // | /| theTria2: (2 3 4 6 8 9) or (3 4 2 8 9 6) or (4 2 3 9 6 8)
552 const SMDS_MeshNode* N1 [6];
553 const SMDS_MeshNode* N2 [6];
554 if(!GetNodesFromTwoTria(theTria1,theTria2,N1,N2))
556 // now we receive following N1 and N2 (using numeration as above image)
557 // tria1 : (1 2 4 5 9 7) and tria2 : (3 4 2 8 9 6)
558 // i.e. first nodes from both arrays determ new diagonal
560 const SMDS_MeshNode* N1new [6];
561 const SMDS_MeshNode* N2new [6];
574 // replaces nodes in faces
575 GetMeshDS()->ChangeElementNodes( theTria1, N1new, 6 );
576 GetMeshDS()->ChangeElementNodes( theTria2, N2new, 6 );
581 //=======================================================================
582 //function : findTriangles
583 //purpose : find triangles sharing theNode1-theNode2 link
584 //=======================================================================
586 static bool findTriangles(const SMDS_MeshNode * theNode1,
587 const SMDS_MeshNode * theNode2,
588 const SMDS_MeshElement*& theTria1,
589 const SMDS_MeshElement*& theTria2)
591 if ( !theNode1 || !theNode2 ) return false;
593 theTria1 = theTria2 = 0;
595 set< const SMDS_MeshElement* > emap;
596 SMDS_ElemIteratorPtr it = theNode1->GetInverseElementIterator(SMDSAbs_Face);
598 const SMDS_MeshElement* elem = it->next();
599 if ( elem->NbNodes() == 3 )
602 it = theNode2->GetInverseElementIterator(SMDSAbs_Face);
604 const SMDS_MeshElement* elem = it->next();
605 if ( emap.find( elem ) != emap.end() )
607 // theTria1 must be element with minimum ID
608 if( theTria1->GetID() < elem->GetID() ) {
621 return ( theTria1 && theTria2 );
624 //=======================================================================
625 //function : InverseDiag
626 //purpose : Replace two neighbour triangles sharing theNode1-theNode2 link
627 // with ones built on the same 4 nodes but having other common link.
628 // Return false if proper faces not found
629 //=======================================================================
631 bool SMESH_MeshEditor::InverseDiag (const SMDS_MeshNode * theNode1,
632 const SMDS_MeshNode * theNode2)
634 myLastCreatedElems.Clear();
635 myLastCreatedNodes.Clear();
637 MESSAGE( "::InverseDiag()" );
639 const SMDS_MeshElement *tr1, *tr2;
640 if ( !findTriangles( theNode1, theNode2, tr1, tr2 ))
643 const SMDS_FaceOfNodes* F1 = dynamic_cast<const SMDS_FaceOfNodes*>( tr1 );
644 //if (!F1) return false;
645 const SMDS_FaceOfNodes* F2 = dynamic_cast<const SMDS_FaceOfNodes*>( tr2 );
646 //if (!F2) return false;
649 // 1 +--+ A tr1: ( 1 A B ) A->2 ( 1 2 B ) 1 +--+ A
650 // | /| tr2: ( B A 2 ) B->1 ( 1 A 2 ) |\ |
654 // put nodes in array
655 // and find indices of 1,2 and of A in tr1 and of B in tr2
656 int i, iA1 = 0, i1 = 0;
657 const SMDS_MeshNode* aNodes1 [3];
658 SMDS_ElemIteratorPtr it;
659 for (i = 0, it = tr1->nodesIterator(); it->more(); i++ ) {
660 aNodes1[ i ] = static_cast<const SMDS_MeshNode*>( it->next() );
661 if ( aNodes1[ i ] == theNode1 )
662 iA1 = i; // node A in tr1
663 else if ( aNodes1[ i ] != theNode2 )
667 const SMDS_MeshNode* aNodes2 [3];
668 for (i = 0, it = tr2->nodesIterator(); it->more(); i++ ) {
669 aNodes2[ i ] = static_cast<const SMDS_MeshNode*>( it->next() );
670 if ( aNodes2[ i ] == theNode2 )
671 iB2 = i; // node B in tr2
672 else if ( aNodes2[ i ] != theNode1 )
676 // nodes 1 and 2 should not be the same
677 if ( aNodes1[ i1 ] == aNodes2[ i2 ] )
681 aNodes1[ iA1 ] = aNodes2[ i2 ];
683 aNodes2[ iB2 ] = aNodes1[ i1 ];
685 //MESSAGE( tr1 << tr2 );
687 GetMeshDS()->ChangeElementNodes( tr1, aNodes1, 3 );
688 GetMeshDS()->ChangeElementNodes( tr2, aNodes2, 3 );
690 //MESSAGE( tr1 << tr2 );
695 // check case of quadratic faces
696 const SMDS_QuadraticFaceOfNodes* QF1 =
697 dynamic_cast<const SMDS_QuadraticFaceOfNodes*> (tr1);
698 if(!QF1) return false;
699 const SMDS_QuadraticFaceOfNodes* QF2 =
700 dynamic_cast<const SMDS_QuadraticFaceOfNodes*> (tr2);
701 if(!QF2) return false;
702 return InverseDiag(tr1,tr2);
705 //=======================================================================
706 //function : getQuadrangleNodes
707 //purpose : fill theQuadNodes - nodes of a quadrangle resulting from
708 // fusion of triangles tr1 and tr2 having shared link on
709 // theNode1 and theNode2
710 //=======================================================================
712 bool getQuadrangleNodes(const SMDS_MeshNode * theQuadNodes [],
713 const SMDS_MeshNode * theNode1,
714 const SMDS_MeshNode * theNode2,
715 const SMDS_MeshElement * tr1,
716 const SMDS_MeshElement * tr2 )
718 if( tr1->NbNodes() != tr2->NbNodes() )
720 // find the 4-th node to insert into tr1
721 const SMDS_MeshNode* n4 = 0;
722 SMDS_ElemIteratorPtr it = tr2->nodesIterator();
724 while ( !n4 && i<3 ) {
725 const SMDS_MeshNode * n = cast2Node( it->next() );
727 bool isDiag = ( n == theNode1 || n == theNode2 );
731 // Make an array of nodes to be in a quadrangle
732 int iNode = 0, iFirstDiag = -1;
733 it = tr1->nodesIterator();
736 const SMDS_MeshNode * n = cast2Node( it->next() );
738 bool isDiag = ( n == theNode1 || n == theNode2 );
740 if ( iFirstDiag < 0 )
742 else if ( iNode - iFirstDiag == 1 )
743 theQuadNodes[ iNode++ ] = n4; // insert the 4-th node between diagonal nodes
745 else if ( n == n4 ) {
746 return false; // tr1 and tr2 should not have all the same nodes
748 theQuadNodes[ iNode++ ] = n;
750 if ( iNode == 3 ) // diagonal nodes have 0 and 2 indices
751 theQuadNodes[ iNode ] = n4;
756 //=======================================================================
757 //function : DeleteDiag
758 //purpose : Replace two neighbour triangles sharing theNode1-theNode2 link
759 // with a quadrangle built on the same 4 nodes.
760 // Return false if proper faces not found
761 //=======================================================================
763 bool SMESH_MeshEditor::DeleteDiag (const SMDS_MeshNode * theNode1,
764 const SMDS_MeshNode * theNode2)
766 myLastCreatedElems.Clear();
767 myLastCreatedNodes.Clear();
769 MESSAGE( "::DeleteDiag()" );
771 const SMDS_MeshElement *tr1, *tr2;
772 if ( !findTriangles( theNode1, theNode2, tr1, tr2 ))
775 const SMDS_FaceOfNodes* F1 = dynamic_cast<const SMDS_FaceOfNodes*>( tr1 );
776 //if (!F1) return false;
777 const SMDS_FaceOfNodes* F2 = dynamic_cast<const SMDS_FaceOfNodes*>( tr2 );
778 //if (!F2) return false;
781 const SMDS_MeshNode* aNodes [ 4 ];
782 if ( ! getQuadrangleNodes( aNodes, theNode1, theNode2, tr1, tr2 ))
785 //MESSAGE( endl << tr1 << tr2 );
787 GetMeshDS()->ChangeElementNodes( tr1, aNodes, 4 );
788 myLastCreatedElems.Append(tr1);
789 GetMeshDS()->RemoveElement( tr2 );
791 //MESSAGE( endl << tr1 );
796 // check case of quadratic faces
797 const SMDS_QuadraticFaceOfNodes* QF1 =
798 dynamic_cast<const SMDS_QuadraticFaceOfNodes*> (tr1);
799 if(!QF1) return false;
800 const SMDS_QuadraticFaceOfNodes* QF2 =
801 dynamic_cast<const SMDS_QuadraticFaceOfNodes*> (tr2);
802 if(!QF2) return false;
805 // 1 +--+--+ 2 tr1: (1 2 4 5 9 7) or (2 4 1 9 7 5) or (4 1 2 7 5 9)
806 // | /| tr2: (2 3 4 6 8 9) or (3 4 2 8 9 6) or (4 2 3 9 6 8)
814 const SMDS_MeshNode* N1 [6];
815 const SMDS_MeshNode* N2 [6];
816 if(!GetNodesFromTwoTria(tr1,tr2,N1,N2))
818 // now we receive following N1 and N2 (using numeration as above image)
819 // tria1 : (1 2 4 5 9 7) and tria2 : (3 4 2 8 9 6)
820 // i.e. first nodes from both arrays determ new diagonal
822 const SMDS_MeshNode* aNodes[8];
832 GetMeshDS()->ChangeElementNodes( tr1, aNodes, 8 );
833 myLastCreatedElems.Append(tr1);
834 GetMeshDS()->RemoveElement( tr2 );
836 // remove middle node (9)
837 GetMeshDS()->RemoveNode( N1[4] );
842 //=======================================================================
843 //function : Reorient
844 //purpose : Reverse theElement orientation
845 //=======================================================================
847 bool SMESH_MeshEditor::Reorient (const SMDS_MeshElement * theElem)
849 myLastCreatedElems.Clear();
850 myLastCreatedNodes.Clear();
854 SMDS_ElemIteratorPtr it = theElem->nodesIterator();
855 if ( !it || !it->more() )
858 switch ( theElem->GetType() ) {
862 if(!theElem->IsQuadratic()) {
863 int i = theElem->NbNodes();
864 vector<const SMDS_MeshNode*> aNodes( i );
866 aNodes[ --i ]= static_cast<const SMDS_MeshNode*>( it->next() );
867 return GetMeshDS()->ChangeElementNodes( theElem, &aNodes[0], theElem->NbNodes() );
870 // quadratic elements
871 if(theElem->GetType()==SMDSAbs_Edge) {
872 vector<const SMDS_MeshNode*> aNodes(3);
873 aNodes[1]= static_cast<const SMDS_MeshNode*>( it->next() );
874 aNodes[0]= static_cast<const SMDS_MeshNode*>( it->next() );
875 aNodes[2]= static_cast<const SMDS_MeshNode*>( it->next() );
876 return GetMeshDS()->ChangeElementNodes( theElem, &aNodes[0], 3 );
879 int nbn = theElem->NbNodes();
880 vector<const SMDS_MeshNode*> aNodes(nbn);
881 aNodes[0]= static_cast<const SMDS_MeshNode*>( it->next() );
883 for(; i<nbn/2; i++) {
884 aNodes[nbn/2-i]= static_cast<const SMDS_MeshNode*>( it->next() );
886 for(i=0; i<nbn/2; i++) {
887 aNodes[nbn-i-1]= static_cast<const SMDS_MeshNode*>( it->next() );
889 return GetMeshDS()->ChangeElementNodes( theElem, &aNodes[0], nbn );
893 case SMDSAbs_Volume: {
894 if (theElem->IsPoly()) {
895 const SMDS_PolyhedralVolumeOfNodes* aPolyedre =
896 static_cast<const SMDS_PolyhedralVolumeOfNodes*>( theElem );
898 MESSAGE("Warning: bad volumic element");
902 int nbFaces = aPolyedre->NbFaces();
903 vector<const SMDS_MeshNode *> poly_nodes;
904 vector<int> quantities (nbFaces);
906 // reverse each face of the polyedre
907 for (int iface = 1; iface <= nbFaces; iface++) {
908 int inode, nbFaceNodes = aPolyedre->NbFaceNodes(iface);
909 quantities[iface - 1] = nbFaceNodes;
911 for (inode = nbFaceNodes; inode >= 1; inode--) {
912 const SMDS_MeshNode* curNode = aPolyedre->GetFaceNode(iface, inode);
913 poly_nodes.push_back(curNode);
917 return GetMeshDS()->ChangePolyhedronNodes( theElem, poly_nodes, quantities );
921 SMDS_VolumeTool vTool;
922 if ( !vTool.Set( theElem ))
925 return GetMeshDS()->ChangeElementNodes( theElem, vTool.GetNodes(), vTool.NbNodes() );
934 //=======================================================================
935 //function : getBadRate
937 //=======================================================================
939 static double getBadRate (const SMDS_MeshElement* theElem,
940 SMESH::Controls::NumericalFunctorPtr& theCrit)
942 SMESH::Controls::TSequenceOfXYZ P;
943 if ( !theElem || !theCrit->GetPoints( theElem, P ))
945 return theCrit->GetBadRate( theCrit->GetValue( P ), theElem->NbNodes() );
946 //return theCrit->GetBadRate( theCrit->GetValue( theElem->GetID() ), theElem->NbNodes() );
949 //=======================================================================
950 //function : QuadToTri
951 //purpose : Cut quadrangles into triangles.
952 // theCrit is used to select a diagonal to cut
953 //=======================================================================
955 bool SMESH_MeshEditor::QuadToTri (TIDSortedElemSet & theElems,
956 SMESH::Controls::NumericalFunctorPtr theCrit)
958 myLastCreatedElems.Clear();
959 myLastCreatedNodes.Clear();
961 MESSAGE( "::QuadToTri()" );
963 if ( !theCrit.get() )
966 SMESHDS_Mesh * aMesh = GetMeshDS();
968 Handle(Geom_Surface) surface;
969 SMESH_MesherHelper helper( *GetMesh() );
971 TIDSortedElemSet::iterator itElem;
972 for ( itElem = theElems.begin(); itElem != theElems.end(); itElem++ ) {
973 const SMDS_MeshElement* elem = *itElem;
974 if ( !elem || elem->GetType() != SMDSAbs_Face )
976 if ( elem->NbNodes() != ( elem->IsQuadratic() ? 8 : 4 ))
979 // retrieve element nodes
980 const SMDS_MeshNode* aNodes [8];
981 SMDS_ElemIteratorPtr itN = elem->nodesIterator();
983 while ( itN->more() )
984 aNodes[ i++ ] = static_cast<const SMDS_MeshNode*>( itN->next() );
986 // compare two sets of possible triangles
987 double aBadRate1, aBadRate2; // to what extent a set is bad
988 SMDS_FaceOfNodes tr1 ( aNodes[0], aNodes[1], aNodes[2] );
989 SMDS_FaceOfNodes tr2 ( aNodes[2], aNodes[3], aNodes[0] );
990 aBadRate1 = getBadRate( &tr1, theCrit ) + getBadRate( &tr2, theCrit );
992 SMDS_FaceOfNodes tr3 ( aNodes[1], aNodes[2], aNodes[3] );
993 SMDS_FaceOfNodes tr4 ( aNodes[3], aNodes[0], aNodes[1] );
994 aBadRate2 = getBadRate( &tr3, theCrit ) + getBadRate( &tr4, theCrit );
996 int aShapeId = FindShape( elem );
997 const SMDS_MeshElement* newElem = 0;
999 if( !elem->IsQuadratic() ) {
1001 // split liner quadrangle
1003 if ( aBadRate1 <= aBadRate2 ) {
1004 // tr1 + tr2 is better
1005 aMesh->ChangeElementNodes( elem, aNodes, 3 );
1006 newElem = aMesh->AddFace( aNodes[2], aNodes[3], aNodes[0] );
1009 // tr3 + tr4 is better
1010 aMesh->ChangeElementNodes( elem, &aNodes[1], 3 );
1011 newElem = aMesh->AddFace( aNodes[3], aNodes[0], aNodes[1] );
1016 // split quadratic quadrangle
1018 // get surface elem is on
1019 if ( aShapeId != helper.GetSubShapeID() ) {
1023 shape = aMesh->IndexToShape( aShapeId );
1024 if ( !shape.IsNull() && shape.ShapeType() == TopAbs_FACE ) {
1025 TopoDS_Face face = TopoDS::Face( shape );
1026 surface = BRep_Tool::Surface( face );
1027 if ( !surface.IsNull() )
1028 helper.SetSubShape( shape );
1032 const SMDS_MeshNode* aNodes [8];
1033 const SMDS_MeshNode* inFaceNode = 0;
1034 SMDS_ElemIteratorPtr itN = elem->nodesIterator();
1036 while ( itN->more() ) {
1037 aNodes[ i++ ] = static_cast<const SMDS_MeshNode*>( itN->next() );
1038 if ( !inFaceNode && helper.GetNodeUVneedInFaceNode() &&
1039 aNodes[ i-1 ]->GetPosition()->GetTypeOfPosition() == SMDS_TOP_FACE )
1041 inFaceNode = aNodes[ i-1 ];
1044 // find middle point for (0,1,2,3)
1045 // and create a node in this point;
1047 if ( surface.IsNull() ) {
1049 p += gp_XYZ(aNodes[i]->X(), aNodes[i]->Y(), aNodes[i]->Z() );
1053 TopoDS_Face face = TopoDS::Face( helper.GetSubShape() );
1056 uv += helper.GetNodeUV( face, aNodes[i], inFaceNode );
1058 p = surface->Value( uv.X(), uv.Y() ).XYZ();
1060 const SMDS_MeshNode* newN = aMesh->AddNode( p.X(), p.Y(), p.Z() );
1061 myLastCreatedNodes.Append(newN);
1063 // create a new element
1064 const SMDS_MeshNode* N[6];
1065 if ( aBadRate1 <= aBadRate2 ) {
1072 newElem = aMesh->AddFace(aNodes[2], aNodes[3], aNodes[0],
1073 aNodes[6], aNodes[7], newN );
1082 newElem = aMesh->AddFace(aNodes[3], aNodes[0], aNodes[1],
1083 aNodes[7], aNodes[4], newN );
1085 aMesh->ChangeElementNodes( elem, N, 6 );
1089 // care of a new element
1091 myLastCreatedElems.Append(newElem);
1092 AddToSameGroups( newElem, elem, aMesh );
1094 // put a new triangle on the same shape
1096 aMesh->SetMeshElementOnShape( newElem, aShapeId );
1101 //=======================================================================
1102 //function : BestSplit
1103 //purpose : Find better diagonal for cutting.
1104 //=======================================================================
1105 int SMESH_MeshEditor::BestSplit (const SMDS_MeshElement* theQuad,
1106 SMESH::Controls::NumericalFunctorPtr theCrit)
1108 myLastCreatedElems.Clear();
1109 myLastCreatedNodes.Clear();
1114 if (!theQuad || theQuad->GetType() != SMDSAbs_Face )
1117 if( theQuad->NbNodes()==4 ||
1118 (theQuad->NbNodes()==8 && theQuad->IsQuadratic()) ) {
1120 // retrieve element nodes
1121 const SMDS_MeshNode* aNodes [4];
1122 SMDS_ElemIteratorPtr itN = theQuad->nodesIterator();
1124 //while (itN->more())
1126 aNodes[ i++ ] = static_cast<const SMDS_MeshNode*>( itN->next() );
1128 // compare two sets of possible triangles
1129 double aBadRate1, aBadRate2; // to what extent a set is bad
1130 SMDS_FaceOfNodes tr1 ( aNodes[0], aNodes[1], aNodes[2] );
1131 SMDS_FaceOfNodes tr2 ( aNodes[2], aNodes[3], aNodes[0] );
1132 aBadRate1 = getBadRate( &tr1, theCrit ) + getBadRate( &tr2, theCrit );
1134 SMDS_FaceOfNodes tr3 ( aNodes[1], aNodes[2], aNodes[3] );
1135 SMDS_FaceOfNodes tr4 ( aNodes[3], aNodes[0], aNodes[1] );
1136 aBadRate2 = getBadRate( &tr3, theCrit ) + getBadRate( &tr4, theCrit );
1138 if (aBadRate1 <= aBadRate2) // tr1 + tr2 is better
1139 return 1; // diagonal 1-3
1141 return 2; // diagonal 2-4
1146 //=======================================================================
1147 //function : AddToSameGroups
1148 //purpose : add elemToAdd to the groups the elemInGroups belongs to
1149 //=======================================================================
1151 void SMESH_MeshEditor::AddToSameGroups (const SMDS_MeshElement* elemToAdd,
1152 const SMDS_MeshElement* elemInGroups,
1153 SMESHDS_Mesh * aMesh)
1155 const set<SMESHDS_GroupBase*>& groups = aMesh->GetGroups();
1156 if (!groups.empty()) {
1157 set<SMESHDS_GroupBase*>::const_iterator grIt = groups.begin();
1158 for ( ; grIt != groups.end(); grIt++ ) {
1159 SMESHDS_Group* group = dynamic_cast<SMESHDS_Group*>( *grIt );
1160 if ( group && group->Contains( elemInGroups ))
1161 group->SMDSGroup().Add( elemToAdd );
1167 //=======================================================================
1168 //function : RemoveElemFromGroups
1169 //purpose : Remove removeelem to the groups the elemInGroups belongs to
1170 //=======================================================================
1171 void SMESH_MeshEditor::RemoveElemFromGroups (const SMDS_MeshElement* removeelem,
1172 SMESHDS_Mesh * aMesh)
1174 const set<SMESHDS_GroupBase*>& groups = aMesh->GetGroups();
1175 if (!groups.empty())
1177 set<SMESHDS_GroupBase*>::const_iterator GrIt = groups.begin();
1178 for (; GrIt != groups.end(); GrIt++)
1180 SMESHDS_Group* grp = dynamic_cast<SMESHDS_Group*>(*GrIt);
1181 if (!grp || grp->IsEmpty()) continue;
1182 grp->SMDSGroup().Remove(removeelem);
1187 //=======================================================================
1188 //function : ReplaceElemInGroups
1189 //purpose : replace elemToRm by elemToAdd in the all groups
1190 //=======================================================================
1192 void SMESH_MeshEditor::ReplaceElemInGroups (const SMDS_MeshElement* elemToRm,
1193 const SMDS_MeshElement* elemToAdd,
1194 SMESHDS_Mesh * aMesh)
1196 const set<SMESHDS_GroupBase*>& groups = aMesh->GetGroups();
1197 if (!groups.empty()) {
1198 set<SMESHDS_GroupBase*>::const_iterator grIt = groups.begin();
1199 for ( ; grIt != groups.end(); grIt++ ) {
1200 SMESHDS_Group* group = dynamic_cast<SMESHDS_Group*>( *grIt );
1201 if ( group && group->SMDSGroup().Remove( elemToRm ) && elemToAdd )
1202 group->SMDSGroup().Add( elemToAdd );
1207 //=======================================================================
1208 //function : QuadToTri
1209 //purpose : Cut quadrangles into triangles.
1210 // theCrit is used to select a diagonal to cut
1211 //=======================================================================
1213 bool SMESH_MeshEditor::QuadToTri (TIDSortedElemSet & theElems,
1214 const bool the13Diag)
1216 myLastCreatedElems.Clear();
1217 myLastCreatedNodes.Clear();
1219 MESSAGE( "::QuadToTri()" );
1221 SMESHDS_Mesh * aMesh = GetMeshDS();
1223 Handle(Geom_Surface) surface;
1224 SMESH_MesherHelper helper( *GetMesh() );
1226 TIDSortedElemSet::iterator itElem;
1227 for ( itElem = theElems.begin(); itElem != theElems.end(); itElem++ ) {
1228 const SMDS_MeshElement* elem = *itElem;
1229 if ( !elem || elem->GetType() != SMDSAbs_Face )
1231 bool isquad = elem->NbNodes()==4 || elem->NbNodes()==8;
1232 if(!isquad) continue;
1234 if(elem->NbNodes()==4) {
1235 // retrieve element nodes
1236 const SMDS_MeshNode* aNodes [4];
1237 SMDS_ElemIteratorPtr itN = elem->nodesIterator();
1239 while ( itN->more() )
1240 aNodes[ i++ ] = static_cast<const SMDS_MeshNode*>( itN->next() );
1242 int aShapeId = FindShape( elem );
1243 const SMDS_MeshElement* newElem = 0;
1245 aMesh->ChangeElementNodes( elem, aNodes, 3 );
1246 newElem = aMesh->AddFace( aNodes[2], aNodes[3], aNodes[0] );
1249 aMesh->ChangeElementNodes( elem, &aNodes[1], 3 );
1250 newElem = aMesh->AddFace( aNodes[3], aNodes[0], aNodes[1] );
1252 myLastCreatedElems.Append(newElem);
1253 // put a new triangle on the same shape and add to the same groups
1255 aMesh->SetMeshElementOnShape( newElem, aShapeId );
1256 AddToSameGroups( newElem, elem, aMesh );
1259 // Quadratic quadrangle
1261 if( elem->NbNodes()==8 && elem->IsQuadratic() ) {
1263 // get surface elem is on
1264 int aShapeId = FindShape( elem );
1265 if ( aShapeId != helper.GetSubShapeID() ) {
1269 shape = aMesh->IndexToShape( aShapeId );
1270 if ( !shape.IsNull() && shape.ShapeType() == TopAbs_FACE ) {
1271 TopoDS_Face face = TopoDS::Face( shape );
1272 surface = BRep_Tool::Surface( face );
1273 if ( !surface.IsNull() )
1274 helper.SetSubShape( shape );
1278 const SMDS_MeshNode* aNodes [8];
1279 const SMDS_MeshNode* inFaceNode = 0;
1280 SMDS_ElemIteratorPtr itN = elem->nodesIterator();
1282 while ( itN->more() ) {
1283 aNodes[ i++ ] = static_cast<const SMDS_MeshNode*>( itN->next() );
1284 if ( !inFaceNode && helper.GetNodeUVneedInFaceNode() &&
1285 aNodes[ i-1 ]->GetPosition()->GetTypeOfPosition() == SMDS_TOP_FACE )
1287 inFaceNode = aNodes[ i-1 ];
1291 // find middle point for (0,1,2,3)
1292 // and create a node in this point;
1294 if ( surface.IsNull() ) {
1296 p += gp_XYZ(aNodes[i]->X(), aNodes[i]->Y(), aNodes[i]->Z() );
1300 TopoDS_Face geomFace = TopoDS::Face( helper.GetSubShape() );
1303 uv += helper.GetNodeUV( geomFace, aNodes[i], inFaceNode );
1305 p = surface->Value( uv.X(), uv.Y() ).XYZ();
1307 const SMDS_MeshNode* newN = aMesh->AddNode( p.X(), p.Y(), p.Z() );
1308 myLastCreatedNodes.Append(newN);
1310 // create a new element
1311 const SMDS_MeshElement* newElem = 0;
1312 const SMDS_MeshNode* N[6];
1320 newElem = aMesh->AddFace(aNodes[2], aNodes[3], aNodes[0],
1321 aNodes[6], aNodes[7], newN );
1330 newElem = aMesh->AddFace(aNodes[3], aNodes[0], aNodes[1],
1331 aNodes[7], aNodes[4], newN );
1333 myLastCreatedElems.Append(newElem);
1334 aMesh->ChangeElementNodes( elem, N, 6 );
1335 // put a new triangle on the same shape and add to the same groups
1337 aMesh->SetMeshElementOnShape( newElem, aShapeId );
1338 AddToSameGroups( newElem, elem, aMesh );
1345 //=======================================================================
1346 //function : getAngle
1348 //=======================================================================
1350 double getAngle(const SMDS_MeshElement * tr1,
1351 const SMDS_MeshElement * tr2,
1352 const SMDS_MeshNode * n1,
1353 const SMDS_MeshNode * n2)
1355 double angle = 2*PI; // bad angle
1358 SMESH::Controls::TSequenceOfXYZ P1, P2;
1359 if ( !SMESH::Controls::NumericalFunctor::GetPoints( tr1, P1 ) ||
1360 !SMESH::Controls::NumericalFunctor::GetPoints( tr2, P2 ))
1363 if(!tr1->IsQuadratic())
1364 N1 = gp_Vec( P1(2) - P1(1) ) ^ gp_Vec( P1(3) - P1(1) );
1366 N1 = gp_Vec( P1(3) - P1(1) ) ^ gp_Vec( P1(5) - P1(1) );
1367 if ( N1.SquareMagnitude() <= gp::Resolution() )
1369 if(!tr2->IsQuadratic())
1370 N2 = gp_Vec( P2(2) - P2(1) ) ^ gp_Vec( P2(3) - P2(1) );
1372 N2 = gp_Vec( P2(3) - P2(1) ) ^ gp_Vec( P2(5) - P2(1) );
1373 if ( N2.SquareMagnitude() <= gp::Resolution() )
1376 // find the first diagonal node n1 in the triangles:
1377 // take in account a diagonal link orientation
1378 const SMDS_MeshElement *nFirst[2], *tr[] = { tr1, tr2 };
1379 for ( int t = 0; t < 2; t++ ) {
1380 SMDS_ElemIteratorPtr it = tr[ t ]->nodesIterator();
1381 int i = 0, iDiag = -1;
1382 while ( it->more()) {
1383 const SMDS_MeshElement *n = it->next();
1384 if ( n == n1 || n == n2 )
1388 if ( i - iDiag == 1 )
1389 nFirst[ t ] = ( n == n1 ? n2 : n1 );
1397 if ( nFirst[ 0 ] == nFirst[ 1 ] )
1400 angle = N1.Angle( N2 );
1405 // =================================================
1406 // class generating a unique ID for a pair of nodes
1407 // and able to return nodes by that ID
1408 // =================================================
1412 LinkID_Gen( const SMESHDS_Mesh* theMesh )
1413 :myMesh( theMesh ), myMaxID( theMesh->MaxNodeID() + 1)
1416 long GetLinkID (const SMDS_MeshNode * n1,
1417 const SMDS_MeshNode * n2) const
1419 return ( Min(n1->GetID(),n2->GetID()) * myMaxID + Max(n1->GetID(),n2->GetID()));
1422 bool GetNodes (const long theLinkID,
1423 const SMDS_MeshNode* & theNode1,
1424 const SMDS_MeshNode* & theNode2) const
1426 theNode1 = myMesh->FindNode( theLinkID / myMaxID );
1427 if ( !theNode1 ) return false;
1428 theNode2 = myMesh->FindNode( theLinkID % myMaxID );
1429 if ( !theNode2 ) return false;
1435 const SMESHDS_Mesh* myMesh;
1440 //=======================================================================
1441 //function : TriToQuad
1442 //purpose : Fuse neighbour triangles into quadrangles.
1443 // theCrit is used to select a neighbour to fuse with.
1444 // theMaxAngle is a max angle between element normals at which
1445 // fusion is still performed.
1446 //=======================================================================
1448 bool SMESH_MeshEditor::TriToQuad (TIDSortedElemSet & theElems,
1449 SMESH::Controls::NumericalFunctorPtr theCrit,
1450 const double theMaxAngle)
1452 myLastCreatedElems.Clear();
1453 myLastCreatedNodes.Clear();
1455 MESSAGE( "::TriToQuad()" );
1457 if ( !theCrit.get() )
1460 SMESHDS_Mesh * aMesh = GetMeshDS();
1462 // Prepare data for algo: build
1463 // 1. map of elements with their linkIDs
1464 // 2. map of linkIDs with their elements
1466 map< SMESH_TLink, list< const SMDS_MeshElement* > > mapLi_listEl;
1467 map< SMESH_TLink, list< const SMDS_MeshElement* > >::iterator itLE;
1468 map< const SMDS_MeshElement*, set< SMESH_TLink > > mapEl_setLi;
1469 map< const SMDS_MeshElement*, set< SMESH_TLink > >::iterator itEL;
1471 TIDSortedElemSet::iterator itElem;
1472 for ( itElem = theElems.begin(); itElem != theElems.end(); itElem++ ) {
1473 const SMDS_MeshElement* elem = *itElem;
1474 if(!elem || elem->GetType() != SMDSAbs_Face ) continue;
1475 bool IsTria = elem->NbNodes()==3 || (elem->NbNodes()==6 && elem->IsQuadratic());
1476 if(!IsTria) continue;
1478 // retrieve element nodes
1479 const SMDS_MeshNode* aNodes [4];
1480 SMDS_ElemIteratorPtr itN = elem->nodesIterator();
1483 aNodes[ i++ ] = cast2Node( itN->next() );
1484 aNodes[ 3 ] = aNodes[ 0 ];
1487 for ( i = 0; i < 3; i++ ) {
1488 SMESH_TLink link( aNodes[i], aNodes[i+1] );
1489 // check if elements sharing a link can be fused
1490 itLE = mapLi_listEl.find( link );
1491 if ( itLE != mapLi_listEl.end() ) {
1492 if ((*itLE).second.size() > 1 ) // consider only 2 elems adjacent by a link
1494 const SMDS_MeshElement* elem2 = (*itLE).second.front();
1495 //if ( FindShape( elem ) != FindShape( elem2 ))
1496 // continue; // do not fuse triangles laying on different shapes
1497 if ( getAngle( elem, elem2, aNodes[i], aNodes[i+1] ) > theMaxAngle )
1498 continue; // avoid making badly shaped quads
1499 (*itLE).second.push_back( elem );
1502 mapLi_listEl[ link ].push_back( elem );
1504 mapEl_setLi [ elem ].insert( link );
1507 // Clean the maps from the links shared by a sole element, ie
1508 // links to which only one element is bound in mapLi_listEl
1510 for ( itLE = mapLi_listEl.begin(); itLE != mapLi_listEl.end(); itLE++ ) {
1511 int nbElems = (*itLE).second.size();
1512 if ( nbElems < 2 ) {
1513 const SMDS_MeshElement* elem = (*itLE).second.front();
1514 SMESH_TLink link = (*itLE).first;
1515 mapEl_setLi[ elem ].erase( link );
1516 if ( mapEl_setLi[ elem ].empty() )
1517 mapEl_setLi.erase( elem );
1521 // Algo: fuse triangles into quadrangles
1523 while ( ! mapEl_setLi.empty() ) {
1524 // Look for the start element:
1525 // the element having the least nb of shared links
1526 const SMDS_MeshElement* startElem = 0;
1528 for ( itEL = mapEl_setLi.begin(); itEL != mapEl_setLi.end(); itEL++ ) {
1529 int nbLinks = (*itEL).second.size();
1530 if ( nbLinks < minNbLinks ) {
1531 startElem = (*itEL).first;
1532 minNbLinks = nbLinks;
1533 if ( minNbLinks == 1 )
1538 // search elements to fuse starting from startElem or links of elements
1539 // fused earlyer - startLinks
1540 list< SMESH_TLink > startLinks;
1541 while ( startElem || !startLinks.empty() ) {
1542 while ( !startElem && !startLinks.empty() ) {
1543 // Get an element to start, by a link
1544 SMESH_TLink linkId = startLinks.front();
1545 startLinks.pop_front();
1546 itLE = mapLi_listEl.find( linkId );
1547 if ( itLE != mapLi_listEl.end() ) {
1548 list< const SMDS_MeshElement* > & listElem = (*itLE).second;
1549 list< const SMDS_MeshElement* >::iterator itE = listElem.begin();
1550 for ( ; itE != listElem.end() ; itE++ )
1551 if ( mapEl_setLi.find( (*itE) ) != mapEl_setLi.end() )
1553 mapLi_listEl.erase( itLE );
1558 // Get candidates to be fused
1559 const SMDS_MeshElement *tr1 = startElem, *tr2 = 0, *tr3 = 0;
1560 const SMESH_TLink *link12, *link13;
1562 ASSERT( mapEl_setLi.find( tr1 ) != mapEl_setLi.end() );
1563 set< SMESH_TLink >& setLi = mapEl_setLi[ tr1 ];
1564 ASSERT( !setLi.empty() );
1565 set< SMESH_TLink >::iterator itLi;
1566 for ( itLi = setLi.begin(); itLi != setLi.end(); itLi++ )
1568 const SMESH_TLink & link = (*itLi);
1569 itLE = mapLi_listEl.find( link );
1570 if ( itLE == mapLi_listEl.end() )
1573 const SMDS_MeshElement* elem = (*itLE).second.front();
1575 elem = (*itLE).second.back();
1576 mapLi_listEl.erase( itLE );
1577 if ( mapEl_setLi.find( elem ) == mapEl_setLi.end())
1588 // add other links of elem to list of links to re-start from
1589 set< SMESH_TLink >& links = mapEl_setLi[ elem ];
1590 set< SMESH_TLink >::iterator it;
1591 for ( it = links.begin(); it != links.end(); it++ ) {
1592 const SMESH_TLink& link2 = (*it);
1593 if ( link2 != link )
1594 startLinks.push_back( link2 );
1598 // Get nodes of possible quadrangles
1599 const SMDS_MeshNode *n12 [4], *n13 [4];
1600 bool Ok12 = false, Ok13 = false;
1601 const SMDS_MeshNode *linkNode1, *linkNode2;
1603 linkNode1 = link12->first;
1604 linkNode2 = link12->second;
1605 if ( tr2 && getQuadrangleNodes( n12, linkNode1, linkNode2, tr1, tr2 ))
1609 linkNode1 = link13->first;
1610 linkNode2 = link13->second;
1611 if ( tr3 && getQuadrangleNodes( n13, linkNode1, linkNode2, tr1, tr3 ))
1615 // Choose a pair to fuse
1616 if ( Ok12 && Ok13 ) {
1617 SMDS_FaceOfNodes quad12 ( n12[ 0 ], n12[ 1 ], n12[ 2 ], n12[ 3 ] );
1618 SMDS_FaceOfNodes quad13 ( n13[ 0 ], n13[ 1 ], n13[ 2 ], n13[ 3 ] );
1619 double aBadRate12 = getBadRate( &quad12, theCrit );
1620 double aBadRate13 = getBadRate( &quad13, theCrit );
1621 if ( aBadRate13 < aBadRate12 )
1628 // and remove fused elems and removed links from the maps
1629 mapEl_setLi.erase( tr1 );
1631 mapEl_setLi.erase( tr2 );
1632 mapLi_listEl.erase( *link12 );
1633 if(tr1->NbNodes()==3) {
1634 if( tr1->GetID() < tr2->GetID() ) {
1635 aMesh->ChangeElementNodes( tr1, n12, 4 );
1636 myLastCreatedElems.Append(tr1);
1637 aMesh->RemoveElement( tr2 );
1640 aMesh->ChangeElementNodes( tr2, n12, 4 );
1641 myLastCreatedElems.Append(tr2);
1642 aMesh->RemoveElement( tr1);
1646 const SMDS_MeshNode* N1 [6];
1647 const SMDS_MeshNode* N2 [6];
1648 GetNodesFromTwoTria(tr1,tr2,N1,N2);
1649 // now we receive following N1 and N2 (using numeration as above image)
1650 // tria1 : (1 2 4 5 9 7) and tria2 : (3 4 2 8 9 6)
1651 // i.e. first nodes from both arrays determ new diagonal
1652 const SMDS_MeshNode* aNodes[8];
1661 if( tr1->GetID() < tr2->GetID() ) {
1662 GetMeshDS()->ChangeElementNodes( tr1, aNodes, 8 );
1663 myLastCreatedElems.Append(tr1);
1664 GetMeshDS()->RemoveElement( tr2 );
1667 GetMeshDS()->ChangeElementNodes( tr2, aNodes, 8 );
1668 myLastCreatedElems.Append(tr2);
1669 GetMeshDS()->RemoveElement( tr1 );
1671 // remove middle node (9)
1672 GetMeshDS()->RemoveNode( N1[4] );
1676 mapEl_setLi.erase( tr3 );
1677 mapLi_listEl.erase( *link13 );
1678 if(tr1->NbNodes()==3) {
1679 if( tr1->GetID() < tr2->GetID() ) {
1680 aMesh->ChangeElementNodes( tr1, n13, 4 );
1681 myLastCreatedElems.Append(tr1);
1682 aMesh->RemoveElement( tr3 );
1685 aMesh->ChangeElementNodes( tr3, n13, 4 );
1686 myLastCreatedElems.Append(tr3);
1687 aMesh->RemoveElement( tr1 );
1691 const SMDS_MeshNode* N1 [6];
1692 const SMDS_MeshNode* N2 [6];
1693 GetNodesFromTwoTria(tr1,tr3,N1,N2);
1694 // now we receive following N1 and N2 (using numeration as above image)
1695 // tria1 : (1 2 4 5 9 7) and tria2 : (3 4 2 8 9 6)
1696 // i.e. first nodes from both arrays determ new diagonal
1697 const SMDS_MeshNode* aNodes[8];
1706 if( tr1->GetID() < tr2->GetID() ) {
1707 GetMeshDS()->ChangeElementNodes( tr1, aNodes, 8 );
1708 myLastCreatedElems.Append(tr1);
1709 GetMeshDS()->RemoveElement( tr3 );
1712 GetMeshDS()->ChangeElementNodes( tr3, aNodes, 8 );
1713 myLastCreatedElems.Append(tr3);
1714 GetMeshDS()->RemoveElement( tr1 );
1716 // remove middle node (9)
1717 GetMeshDS()->RemoveNode( N1[4] );
1721 // Next element to fuse: the rejected one
1723 startElem = Ok12 ? tr3 : tr2;
1725 } // if ( startElem )
1726 } // while ( startElem || !startLinks.empty() )
1727 } // while ( ! mapEl_setLi.empty() )
1733 /*#define DUMPSO(txt) \
1734 // cout << txt << endl;
1735 //=============================================================================
1739 //=============================================================================
1740 static void swap( int i1, int i2, int idNodes[], gp_Pnt P[] )
1744 int tmp = idNodes[ i1 ];
1745 idNodes[ i1 ] = idNodes[ i2 ];
1746 idNodes[ i2 ] = tmp;
1747 gp_Pnt Ptmp = P[ i1 ];
1750 DUMPSO( i1 << "(" << idNodes[ i2 ] << ") <-> " << i2 << "(" << idNodes[ i1 ] << ")");
1753 //=======================================================================
1754 //function : SortQuadNodes
1755 //purpose : Set 4 nodes of a quadrangle face in a good order.
1756 // Swap 1<->2 or 2<->3 nodes and correspondingly return
1758 //=======================================================================
1760 int SMESH_MeshEditor::SortQuadNodes (const SMDS_Mesh * theMesh,
1765 for ( i = 0; i < 4; i++ ) {
1766 const SMDS_MeshNode *n = theMesh->FindNode( idNodes[i] );
1768 P[ i ].SetCoord( n->X(), n->Y(), n->Z() );
1771 gp_Vec V1(P[0], P[1]);
1772 gp_Vec V2(P[0], P[2]);
1773 gp_Vec V3(P[0], P[3]);
1775 gp_Vec Cross1 = V1 ^ V2;
1776 gp_Vec Cross2 = V2 ^ V3;
1779 if (Cross1.Dot(Cross2) < 0)
1784 if (Cross1.Dot(Cross2) < 0)
1788 swap ( i, i + 1, idNodes, P );
1790 // for ( int ii = 0; ii < 4; ii++ ) {
1791 // const SMDS_MeshNode *n = theMesh->FindNode( idNodes[ii] );
1792 // DUMPSO( ii << "(" << idNodes[ii] <<") : "<<n->X()<<" "<<n->Y()<<" "<<n->Z());
1798 //=======================================================================
1799 //function : SortHexaNodes
1800 //purpose : Set 8 nodes of a hexahedron in a good order.
1801 // Return success status
1802 //=======================================================================
1804 bool SMESH_MeshEditor::SortHexaNodes (const SMDS_Mesh * theMesh,
1809 DUMPSO( "INPUT: ========================================");
1810 for ( i = 0; i < 8; i++ ) {
1811 const SMDS_MeshNode *n = theMesh->FindNode( idNodes[i] );
1812 if ( !n ) return false;
1813 P[ i ].SetCoord( n->X(), n->Y(), n->Z() );
1814 DUMPSO( i << "(" << idNodes[i] <<") : "<<n->X()<<" "<<n->Y()<<" "<<n->Z());
1816 DUMPSO( "========================================");
1819 set<int> faceNodes; // ids of bottom face nodes, to be found
1820 set<int> checkedId1; // ids of tried 2-nd nodes
1821 Standard_Real leastDist = DBL_MAX; // dist of the 4-th node from 123 plane
1822 const Standard_Real tol = 1.e-6; // tolerance to find nodes in plane
1823 int iMin, iLoop1 = 0;
1825 // Loop to try the 2-nd nodes
1827 while ( leastDist > DBL_MIN && ++iLoop1 < 8 )
1829 // Find not checked 2-nd node
1830 for ( i = 1; i < 8; i++ )
1831 if ( checkedId1.find( idNodes[i] ) == checkedId1.end() ) {
1832 int id1 = idNodes[i];
1833 swap ( 1, i, idNodes, P );
1834 checkedId1.insert ( id1 );
1838 // Find the 3-d node so that 1-2-3 triangle to be on a hexa face,
1839 // ie that all but meybe one (id3 which is on the same face) nodes
1840 // lay on the same side from the triangle plane.
1842 bool manyInPlane = false; // more than 4 nodes lay in plane
1844 while ( ++iLoop2 < 6 ) {
1846 // get 1-2-3 plane coeffs
1847 Standard_Real A, B, C, D;
1848 gp_Vec N = gp_Vec (P[0], P[1]).Crossed( gp_Vec (P[0], P[2]) );
1849 if ( N.SquareMagnitude() > gp::Resolution() )
1851 gp_Pln pln ( P[0], N );
1852 pln.Coefficients( A, B, C, D );
1854 // find the node (iMin) closest to pln
1855 Standard_Real dist[ 8 ], minDist = DBL_MAX;
1857 for ( i = 3; i < 8; i++ ) {
1858 dist[i] = A * P[i].X() + B * P[i].Y() + C * P[i].Z() + D;
1859 if ( fabs( dist[i] ) < minDist ) {
1860 minDist = fabs( dist[i] );
1863 if ( fabs( dist[i] ) <= tol )
1864 idInPln.insert( idNodes[i] );
1867 // there should not be more than 4 nodes in bottom plane
1868 if ( idInPln.size() > 1 )
1870 DUMPSO( "### idInPln.size() = " << idInPln.size());
1871 // idInPlane does not contain the first 3 nodes
1872 if ( manyInPlane || idInPln.size() == 5)
1873 return false; // all nodes in one plane
1876 // set the 1-st node to be not in plane
1877 for ( i = 3; i < 8; i++ ) {
1878 if ( idInPln.find( idNodes[ i ] ) == idInPln.end() ) {
1879 DUMPSO( "### Reset 0-th node");
1880 swap( 0, i, idNodes, P );
1885 // reset to re-check second nodes
1886 leastDist = DBL_MAX;
1890 break; // from iLoop2;
1893 // check that the other 4 nodes are on the same side
1894 bool sameSide = true;
1895 bool isNeg = dist[ iMin == 3 ? 4 : 3 ] <= 0.;
1896 for ( i = 3; sameSide && i < 8; i++ ) {
1898 sameSide = ( isNeg == dist[i] <= 0.);
1901 // keep best solution
1902 if ( sameSide && minDist < leastDist ) {
1903 leastDist = minDist;
1905 faceNodes.insert( idNodes[ 1 ] );
1906 faceNodes.insert( idNodes[ 2 ] );
1907 faceNodes.insert( idNodes[ iMin ] );
1908 DUMPSO( "loop " << iLoop2 << " id2 " << idNodes[ 1 ] << " id3 " << idNodes[ 2 ]
1909 << " leastDist = " << leastDist);
1910 if ( leastDist <= DBL_MIN )
1915 // set next 3-d node to check
1916 int iNext = 2 + iLoop2;
1918 DUMPSO( "Try 2-nd");
1919 swap ( 2, iNext, idNodes, P );
1921 } // while ( iLoop2 < 6 )
1924 if ( faceNodes.empty() ) return false;
1926 // Put the faceNodes in proper places
1927 for ( i = 4; i < 8; i++ ) {
1928 if ( faceNodes.find( idNodes[ i ] ) != faceNodes.end() ) {
1929 // find a place to put
1931 while ( faceNodes.find( idNodes[ iTo ] ) != faceNodes.end() )
1933 DUMPSO( "Set faceNodes");
1934 swap ( iTo, i, idNodes, P );
1939 // Set nodes of the found bottom face in good order
1940 DUMPSO( " Found bottom face: ");
1941 i = SortQuadNodes( theMesh, idNodes );
1943 gp_Pnt Ptmp = P[ i ];
1948 // for ( int ii = 0; ii < 4; ii++ ) {
1949 // const SMDS_MeshNode *n = theMesh->FindNode( idNodes[ii] );
1950 // DUMPSO( ii << "(" << idNodes[ii] <<") : "<<n->X()<<" "<<n->Y()<<" "<<n->Z());
1953 // Gravity center of the top and bottom faces
1954 gp_Pnt aGCb = ( P[0].XYZ() + P[1].XYZ() + P[2].XYZ() + P[3].XYZ() ) / 4.;
1955 gp_Pnt aGCt = ( P[4].XYZ() + P[5].XYZ() + P[6].XYZ() + P[7].XYZ() ) / 4.;
1957 // Get direction from the bottom to the top face
1958 gp_Vec upDir ( aGCb, aGCt );
1959 Standard_Real upDirSize = upDir.Magnitude();
1960 if ( upDirSize <= gp::Resolution() ) return false;
1963 // Assure that the bottom face normal points up
1964 gp_Vec Nb = gp_Vec (P[0], P[1]).Crossed( gp_Vec (P[0], P[2]) );
1965 Nb += gp_Vec (P[0], P[2]).Crossed( gp_Vec (P[0], P[3]) );
1966 if ( Nb.Dot( upDir ) < 0 ) {
1967 DUMPSO( "Reverse bottom face");
1968 swap( 1, 3, idNodes, P );
1971 // Find 5-th node - the one closest to the 1-st among the last 4 nodes.
1972 Standard_Real minDist = DBL_MAX;
1973 for ( i = 4; i < 8; i++ ) {
1974 // projection of P[i] to the plane defined by P[0] and upDir
1975 gp_Pnt Pp = P[i].Translated( upDir * ( upDir.Dot( gp_Vec( P[i], P[0] ))));
1976 Standard_Real sqDist = P[0].SquareDistance( Pp );
1977 if ( sqDist < minDist ) {
1982 DUMPSO( "Set 4-th");
1983 swap ( 4, iMin, idNodes, P );
1985 // Set nodes of the top face in good order
1986 DUMPSO( "Sort top face");
1987 i = SortQuadNodes( theMesh, &idNodes[4] );
1990 gp_Pnt Ptmp = P[ i ];
1995 // Assure that direction of the top face normal is from the bottom face
1996 gp_Vec Nt = gp_Vec (P[4], P[5]).Crossed( gp_Vec (P[4], P[6]) );
1997 Nt += gp_Vec (P[4], P[6]).Crossed( gp_Vec (P[4], P[7]) );
1998 if ( Nt.Dot( upDir ) < 0 ) {
1999 DUMPSO( "Reverse top face");
2000 swap( 5, 7, idNodes, P );
2003 // DUMPSO( "OUTPUT: ========================================");
2004 // for ( i = 0; i < 8; i++ ) {
2005 // float *p = ugrid->GetPoint(idNodes[i]);
2006 // DUMPSO( i << "(" << idNodes[i] << ") : " << p[0] << " " << p[1] << " " << p[2]);
2012 //================================================================================
2014 * \brief Return nodes linked to the given one
2015 * \param theNode - the node
2016 * \param linkedNodes - the found nodes
2017 * \param type - the type of elements to check
2019 * Medium nodes are ignored
2021 //================================================================================
2023 void SMESH_MeshEditor::GetLinkedNodes( const SMDS_MeshNode* theNode,
2024 TIDSortedElemSet & linkedNodes,
2025 SMDSAbs_ElementType type )
2027 SMDS_ElemIteratorPtr elemIt = theNode->GetInverseElementIterator(type);
2028 while ( elemIt->more() )
2030 const SMDS_MeshElement* elem = elemIt->next();
2031 SMDS_ElemIteratorPtr nodeIt = elem->nodesIterator();
2032 if ( elem->GetType() == SMDSAbs_Volume )
2034 SMDS_VolumeTool vol( elem );
2035 while ( nodeIt->more() ) {
2036 const SMDS_MeshNode* n = cast2Node( nodeIt->next() );
2037 if ( theNode != n && vol.IsLinked( theNode, n ))
2038 linkedNodes.insert( n );
2043 for ( int i = 0; nodeIt->more(); ++i ) {
2044 const SMDS_MeshNode* n = cast2Node( nodeIt->next() );
2045 if ( n == theNode ) {
2046 int iBefore = i - 1;
2048 if ( elem->IsQuadratic() ) {
2049 int nb = elem->NbNodes() / 2;
2050 iAfter = SMESH_MesherHelper::WrapIndex( iAfter, nb );
2051 iBefore = SMESH_MesherHelper::WrapIndex( iBefore, nb );
2053 linkedNodes.insert( elem->GetNode( iAfter ));
2054 linkedNodes.insert( elem->GetNode( iBefore ));
2061 //=======================================================================
2062 //function : laplacianSmooth
2063 //purpose : pulls theNode toward the center of surrounding nodes directly
2064 // connected to that node along an element edge
2065 //=======================================================================
2067 void laplacianSmooth(const SMDS_MeshNode* theNode,
2068 const Handle(Geom_Surface)& theSurface,
2069 map< const SMDS_MeshNode*, gp_XY* >& theUVMap)
2071 // find surrounding nodes
2073 TIDSortedElemSet nodeSet;
2074 SMESH_MeshEditor::GetLinkedNodes( theNode, nodeSet, SMDSAbs_Face );
2076 // compute new coodrs
2078 double coord[] = { 0., 0., 0. };
2079 TIDSortedElemSet::iterator nodeSetIt = nodeSet.begin();
2080 for ( ; nodeSetIt != nodeSet.end(); nodeSetIt++ ) {
2081 const SMDS_MeshNode* node = cast2Node(*nodeSetIt);
2082 if ( theSurface.IsNull() ) { // smooth in 3D
2083 coord[0] += node->X();
2084 coord[1] += node->Y();
2085 coord[2] += node->Z();
2087 else { // smooth in 2D
2088 ASSERT( theUVMap.find( node ) != theUVMap.end() );
2089 gp_XY* uv = theUVMap[ node ];
2090 coord[0] += uv->X();
2091 coord[1] += uv->Y();
2094 int nbNodes = nodeSet.size();
2097 coord[0] /= nbNodes;
2098 coord[1] /= nbNodes;
2100 if ( !theSurface.IsNull() ) {
2101 ASSERT( theUVMap.find( theNode ) != theUVMap.end() );
2102 theUVMap[ theNode ]->SetCoord( coord[0], coord[1] );
2103 gp_Pnt p3d = theSurface->Value( coord[0], coord[1] );
2109 coord[2] /= nbNodes;
2113 const_cast< SMDS_MeshNode* >( theNode )->setXYZ(coord[0],coord[1],coord[2]);
2116 //=======================================================================
2117 //function : centroidalSmooth
2118 //purpose : pulls theNode toward the element-area-weighted centroid of the
2119 // surrounding elements
2120 //=======================================================================
2122 void centroidalSmooth(const SMDS_MeshNode* theNode,
2123 const Handle(Geom_Surface)& theSurface,
2124 map< const SMDS_MeshNode*, gp_XY* >& theUVMap)
2126 gp_XYZ aNewXYZ(0.,0.,0.);
2127 SMESH::Controls::Area anAreaFunc;
2128 double totalArea = 0.;
2133 SMDS_ElemIteratorPtr elemIt = theNode->GetInverseElementIterator(SMDSAbs_Face);
2134 while ( elemIt->more() )
2136 const SMDS_MeshElement* elem = elemIt->next();
2139 gp_XYZ elemCenter(0.,0.,0.);
2140 SMESH::Controls::TSequenceOfXYZ aNodePoints;
2141 SMDS_ElemIteratorPtr itN = elem->nodesIterator();
2142 int nn = elem->NbNodes();
2143 if(elem->IsQuadratic()) nn = nn/2;
2145 //while ( itN->more() ) {
2147 const SMDS_MeshNode* aNode = static_cast<const SMDS_MeshNode*>( itN->next() );
2149 gp_XYZ aP( aNode->X(), aNode->Y(), aNode->Z() );
2150 aNodePoints.push_back( aP );
2151 if ( !theSurface.IsNull() ) { // smooth in 2D
2152 ASSERT( theUVMap.find( aNode ) != theUVMap.end() );
2153 gp_XY* uv = theUVMap[ aNode ];
2154 aP.SetCoord( uv->X(), uv->Y(), 0. );
2158 double elemArea = anAreaFunc.GetValue( aNodePoints );
2159 totalArea += elemArea;
2161 aNewXYZ += elemCenter * elemArea;
2163 aNewXYZ /= totalArea;
2164 if ( !theSurface.IsNull() ) {
2165 theUVMap[ theNode ]->SetCoord( aNewXYZ.X(), aNewXYZ.Y() );
2166 aNewXYZ = theSurface->Value( aNewXYZ.X(), aNewXYZ.Y() ).XYZ();
2171 const_cast< SMDS_MeshNode* >( theNode )->setXYZ(aNewXYZ.X(),aNewXYZ.Y(),aNewXYZ.Z());
2174 //=======================================================================
2175 //function : getClosestUV
2176 //purpose : return UV of closest projection
2177 //=======================================================================
2179 static bool getClosestUV (Extrema_GenExtPS& projector,
2180 const gp_Pnt& point,
2183 projector.Perform( point );
2184 if ( projector.IsDone() ) {
2185 double u, v, minVal = DBL_MAX;
2186 for ( int i = projector.NbExt(); i > 0; i-- )
2187 if ( projector.Value( i ) < minVal ) {
2188 minVal = projector.Value( i );
2189 projector.Point( i ).Parameter( u, v );
2191 result.SetCoord( u, v );
2197 //=======================================================================
2199 //purpose : Smooth theElements during theNbIterations or until a worst
2200 // element has aspect ratio <= theTgtAspectRatio.
2201 // Aspect Ratio varies in range [1.0, inf].
2202 // If theElements is empty, the whole mesh is smoothed.
2203 // theFixedNodes contains additionally fixed nodes. Nodes built
2204 // on edges and boundary nodes are always fixed.
2205 //=======================================================================
2207 void SMESH_MeshEditor::Smooth (TIDSortedElemSet & theElems,
2208 set<const SMDS_MeshNode*> & theFixedNodes,
2209 const SmoothMethod theSmoothMethod,
2210 const int theNbIterations,
2211 double theTgtAspectRatio,
2214 myLastCreatedElems.Clear();
2215 myLastCreatedNodes.Clear();
2217 MESSAGE((theSmoothMethod==LAPLACIAN ? "LAPLACIAN" : "CENTROIDAL") << "--::Smooth()");
2219 if ( theTgtAspectRatio < 1.0 )
2220 theTgtAspectRatio = 1.0;
2222 const double disttol = 1.e-16;
2224 SMESH::Controls::AspectRatio aQualityFunc;
2226 SMESHDS_Mesh* aMesh = GetMeshDS();
2228 if ( theElems.empty() ) {
2229 // add all faces to theElems
2230 SMDS_FaceIteratorPtr fIt = aMesh->facesIterator();
2231 while ( fIt->more() ) {
2232 const SMDS_MeshElement* face = fIt->next();
2233 theElems.insert( face );
2236 // get all face ids theElems are on
2237 set< int > faceIdSet;
2238 TIDSortedElemSet::iterator itElem;
2240 for ( itElem = theElems.begin(); itElem != theElems.end(); itElem++ ) {
2241 int fId = FindShape( *itElem );
2242 // check that corresponding submesh exists and a shape is face
2244 faceIdSet.find( fId ) == faceIdSet.end() &&
2245 aMesh->MeshElements( fId )) {
2246 TopoDS_Shape F = aMesh->IndexToShape( fId );
2247 if ( !F.IsNull() && F.ShapeType() == TopAbs_FACE )
2248 faceIdSet.insert( fId );
2251 faceIdSet.insert( 0 ); // to smooth elements that are not on any TopoDS_Face
2253 // ===============================================
2254 // smooth elements on each TopoDS_Face separately
2255 // ===============================================
2257 set< int >::reverse_iterator fId = faceIdSet.rbegin(); // treate 0 fId at the end
2258 for ( ; fId != faceIdSet.rend(); ++fId ) {
2259 // get face surface and submesh
2260 Handle(Geom_Surface) surface;
2261 SMESHDS_SubMesh* faceSubMesh = 0;
2263 double fToler2 = 0, vPeriod = 0., uPeriod = 0., f,l;
2264 double u1 = 0, u2 = 0, v1 = 0, v2 = 0;
2265 bool isUPeriodic = false, isVPeriodic = false;
2267 face = TopoDS::Face( aMesh->IndexToShape( *fId ));
2268 surface = BRep_Tool::Surface( face );
2269 faceSubMesh = aMesh->MeshElements( *fId );
2270 fToler2 = BRep_Tool::Tolerance( face );
2271 fToler2 *= fToler2 * 10.;
2272 isUPeriodic = surface->IsUPeriodic();
2274 vPeriod = surface->UPeriod();
2275 isVPeriodic = surface->IsVPeriodic();
2277 uPeriod = surface->VPeriod();
2278 surface->Bounds( u1, u2, v1, v2 );
2280 // ---------------------------------------------------------
2281 // for elements on a face, find movable and fixed nodes and
2282 // compute UV for them
2283 // ---------------------------------------------------------
2284 bool checkBoundaryNodes = false;
2285 bool isQuadratic = false;
2286 set<const SMDS_MeshNode*> setMovableNodes;
2287 map< const SMDS_MeshNode*, gp_XY* > uvMap, uvMap2;
2288 list< gp_XY > listUV; // uvs the 2 uvMaps refer to
2289 list< const SMDS_MeshElement* > elemsOnFace;
2291 Extrema_GenExtPS projector;
2292 GeomAdaptor_Surface surfAdaptor;
2293 if ( !surface.IsNull() ) {
2294 surfAdaptor.Load( surface );
2295 projector.Initialize( surfAdaptor, 20,20, 1e-5,1e-5 );
2297 int nbElemOnFace = 0;
2298 itElem = theElems.begin();
2299 // loop on not yet smoothed elements: look for elems on a face
2300 while ( itElem != theElems.end() ) {
2301 if ( faceSubMesh && nbElemOnFace == faceSubMesh->NbElements() )
2302 break; // all elements found
2304 const SMDS_MeshElement* elem = *itElem;
2305 if ( !elem || elem->GetType() != SMDSAbs_Face || elem->NbNodes() < 3 ||
2306 ( faceSubMesh && !faceSubMesh->Contains( elem ))) {
2310 elemsOnFace.push_back( elem );
2311 theElems.erase( itElem++ );
2315 isQuadratic = elem->IsQuadratic();
2317 // get movable nodes of elem
2318 const SMDS_MeshNode* node;
2319 SMDS_TypeOfPosition posType;
2320 SMDS_ElemIteratorPtr itN = elem->nodesIterator();
2321 int nn = 0, nbn = elem->NbNodes();
2322 if(elem->IsQuadratic())
2324 while ( nn++ < nbn ) {
2325 node = static_cast<const SMDS_MeshNode*>( itN->next() );
2326 const SMDS_PositionPtr& pos = node->GetPosition();
2327 posType = pos.get() ? pos->GetTypeOfPosition() : SMDS_TOP_3DSPACE;
2328 if (posType != SMDS_TOP_EDGE &&
2329 posType != SMDS_TOP_VERTEX &&
2330 theFixedNodes.find( node ) == theFixedNodes.end())
2332 // check if all faces around the node are on faceSubMesh
2333 // because a node on edge may be bound to face
2334 SMDS_ElemIteratorPtr eIt = node->GetInverseElementIterator(SMDSAbs_Face);
2336 if ( faceSubMesh ) {
2337 while ( eIt->more() && all ) {
2338 const SMDS_MeshElement* e = eIt->next();
2339 all = faceSubMesh->Contains( e );
2343 setMovableNodes.insert( node );
2345 checkBoundaryNodes = true;
2347 if ( posType == SMDS_TOP_3DSPACE )
2348 checkBoundaryNodes = true;
2351 if ( surface.IsNull() )
2354 // get nodes to check UV
2355 list< const SMDS_MeshNode* > uvCheckNodes;
2356 itN = elem->nodesIterator();
2357 nn = 0; nbn = elem->NbNodes();
2358 if(elem->IsQuadratic())
2360 while ( nn++ < nbn ) {
2361 node = static_cast<const SMDS_MeshNode*>( itN->next() );
2362 if ( uvMap.find( node ) == uvMap.end() )
2363 uvCheckNodes.push_back( node );
2364 // add nodes of elems sharing node
2365 // SMDS_ElemIteratorPtr eIt = node->GetInverseElementIterator(SMDSAbs_Face);
2366 // while ( eIt->more() ) {
2367 // const SMDS_MeshElement* e = eIt->next();
2368 // if ( e != elem ) {
2369 // SMDS_ElemIteratorPtr nIt = e->nodesIterator();
2370 // while ( nIt->more() ) {
2371 // const SMDS_MeshNode* n =
2372 // static_cast<const SMDS_MeshNode*>( nIt->next() );
2373 // if ( uvMap.find( n ) == uvMap.end() )
2374 // uvCheckNodes.push_back( n );
2380 list< const SMDS_MeshNode* >::iterator n = uvCheckNodes.begin();
2381 for ( ; n != uvCheckNodes.end(); ++n ) {
2384 const SMDS_PositionPtr& pos = node->GetPosition();
2385 posType = pos.get() ? pos->GetTypeOfPosition() : SMDS_TOP_3DSPACE;
2387 switch ( posType ) {
2388 case SMDS_TOP_FACE: {
2389 SMDS_FacePosition* fPos = ( SMDS_FacePosition* ) pos.get();
2390 uv.SetCoord( fPos->GetUParameter(), fPos->GetVParameter() );
2393 case SMDS_TOP_EDGE: {
2394 TopoDS_Shape S = aMesh->IndexToShape( pos->GetShapeId() );
2395 Handle(Geom2d_Curve) pcurve;
2396 if ( !S.IsNull() && S.ShapeType() == TopAbs_EDGE )
2397 pcurve = BRep_Tool::CurveOnSurface( TopoDS::Edge( S ), face, f,l );
2398 if ( !pcurve.IsNull() ) {
2399 double u = (( SMDS_EdgePosition* ) pos.get() )->GetUParameter();
2400 uv = pcurve->Value( u ).XY();
2404 case SMDS_TOP_VERTEX: {
2405 TopoDS_Shape S = aMesh->IndexToShape( pos->GetShapeId() );
2406 if ( !S.IsNull() && S.ShapeType() == TopAbs_VERTEX )
2407 uv = BRep_Tool::Parameters( TopoDS::Vertex( S ), face ).XY();
2412 // check existing UV
2413 bool project = true;
2414 gp_Pnt pNode ( node->X(), node->Y(), node->Z() );
2415 double dist1 = DBL_MAX, dist2 = 0;
2416 if ( posType != SMDS_TOP_3DSPACE ) {
2417 dist1 = pNode.SquareDistance( surface->Value( uv.X(), uv.Y() ));
2418 project = dist1 > fToler2;
2420 if ( project ) { // compute new UV
2422 if ( !getClosestUV( projector, pNode, newUV )) {
2423 MESSAGE("Node Projection Failed " << node);
2427 newUV.SetX( ElCLib::InPeriod( newUV.X(), u1, u2 ));
2429 newUV.SetY( ElCLib::InPeriod( newUV.Y(), v1, v2 ));
2431 if ( posType != SMDS_TOP_3DSPACE )
2432 dist2 = pNode.SquareDistance( surface->Value( newUV.X(), newUV.Y() ));
2433 if ( dist2 < dist1 )
2437 // store UV in the map
2438 listUV.push_back( uv );
2439 uvMap.insert( make_pair( node, &listUV.back() ));
2441 } // loop on not yet smoothed elements
2443 if ( !faceSubMesh || nbElemOnFace != faceSubMesh->NbElements() )
2444 checkBoundaryNodes = true;
2446 // fix nodes on mesh boundary
2448 if ( checkBoundaryNodes ) {
2449 map< NLink, int > linkNbMap; // how many times a link encounters in elemsOnFace
2450 map< NLink, int >::iterator link_nb;
2451 // put all elements links to linkNbMap
2452 list< const SMDS_MeshElement* >::iterator elemIt = elemsOnFace.begin();
2453 for ( ; elemIt != elemsOnFace.end(); ++elemIt ) {
2454 const SMDS_MeshElement* elem = (*elemIt);
2455 int nbn = elem->NbNodes();
2456 if(elem->IsQuadratic())
2458 // loop on elem links: insert them in linkNbMap
2459 const SMDS_MeshNode* curNode, *prevNode = elem->GetNode( nbn );
2460 for ( int iN = 0; iN < nbn; ++iN ) {
2461 curNode = elem->GetNode( iN );
2463 if ( curNode < prevNode ) link = make_pair( curNode , prevNode );
2464 else link = make_pair( prevNode , curNode );
2466 link_nb = linkNbMap.find( link );
2467 if ( link_nb == linkNbMap.end() )
2468 linkNbMap.insert( make_pair ( link, 1 ));
2473 // remove nodes that are in links encountered only once from setMovableNodes
2474 for ( link_nb = linkNbMap.begin(); link_nb != linkNbMap.end(); ++link_nb ) {
2475 if ( link_nb->second == 1 ) {
2476 setMovableNodes.erase( link_nb->first.first );
2477 setMovableNodes.erase( link_nb->first.second );
2482 // -----------------------------------------------------
2483 // for nodes on seam edge, compute one more UV ( uvMap2 );
2484 // find movable nodes linked to nodes on seam and which
2485 // are to be smoothed using the second UV ( uvMap2 )
2486 // -----------------------------------------------------
2488 set<const SMDS_MeshNode*> nodesNearSeam; // to smooth using uvMap2
2489 if ( !surface.IsNull() ) {
2490 TopExp_Explorer eExp( face, TopAbs_EDGE );
2491 for ( ; eExp.More(); eExp.Next() ) {
2492 TopoDS_Edge edge = TopoDS::Edge( eExp.Current() );
2493 if ( !BRep_Tool::IsClosed( edge, face ))
2495 SMESHDS_SubMesh* sm = aMesh->MeshElements( edge );
2496 if ( !sm ) continue;
2497 // find out which parameter varies for a node on seam
2500 Handle(Geom2d_Curve) pcurve = BRep_Tool::CurveOnSurface( edge, face, f, l );
2501 if ( pcurve.IsNull() ) continue;
2502 uv1 = pcurve->Value( f );
2504 pcurve = BRep_Tool::CurveOnSurface( edge, face, f, l );
2505 if ( pcurve.IsNull() ) continue;
2506 uv2 = pcurve->Value( f );
2507 int iPar = Abs( uv1.X() - uv2.X() ) > Abs( uv1.Y() - uv2.Y() ) ? 1 : 2;
2509 if ( uv1.Coord( iPar ) > uv2.Coord( iPar )) {
2510 gp_Pnt2d tmp = uv1; uv1 = uv2; uv2 = tmp;
2512 // get nodes on seam and its vertices
2513 list< const SMDS_MeshNode* > seamNodes;
2514 SMDS_NodeIteratorPtr nSeamIt = sm->GetNodes();
2515 while ( nSeamIt->more() ) {
2516 const SMDS_MeshNode* node = nSeamIt->next();
2517 if ( !isQuadratic || !IsMedium( node ))
2518 seamNodes.push_back( node );
2520 TopExp_Explorer vExp( edge, TopAbs_VERTEX );
2521 for ( ; vExp.More(); vExp.Next() ) {
2522 sm = aMesh->MeshElements( vExp.Current() );
2524 nSeamIt = sm->GetNodes();
2525 while ( nSeamIt->more() )
2526 seamNodes.push_back( nSeamIt->next() );
2529 // loop on nodes on seam
2530 list< const SMDS_MeshNode* >::iterator noSeIt = seamNodes.begin();
2531 for ( ; noSeIt != seamNodes.end(); ++noSeIt ) {
2532 const SMDS_MeshNode* nSeam = *noSeIt;
2533 map< const SMDS_MeshNode*, gp_XY* >::iterator n_uv = uvMap.find( nSeam );
2534 if ( n_uv == uvMap.end() )
2537 n_uv->second->SetCoord( iPar, uv1.Coord( iPar ));
2538 // set the second UV
2539 listUV.push_back( *n_uv->second );
2540 listUV.back().SetCoord( iPar, uv2.Coord( iPar ));
2541 if ( uvMap2.empty() )
2542 uvMap2 = uvMap; // copy the uvMap contents
2543 uvMap2[ nSeam ] = &listUV.back();
2545 // collect movable nodes linked to ones on seam in nodesNearSeam
2546 SMDS_ElemIteratorPtr eIt = nSeam->GetInverseElementIterator(SMDSAbs_Face);
2547 while ( eIt->more() ) {
2548 const SMDS_MeshElement* e = eIt->next();
2549 int nbUseMap1 = 0, nbUseMap2 = 0;
2550 SMDS_ElemIteratorPtr nIt = e->nodesIterator();
2551 int nn = 0, nbn = e->NbNodes();
2552 if(e->IsQuadratic()) nbn = nbn/2;
2553 while ( nn++ < nbn )
2555 const SMDS_MeshNode* n =
2556 static_cast<const SMDS_MeshNode*>( nIt->next() );
2558 setMovableNodes.find( n ) == setMovableNodes.end() )
2560 // add only nodes being closer to uv2 than to uv1
2561 gp_Pnt pMid (0.5 * ( n->X() + nSeam->X() ),
2562 0.5 * ( n->Y() + nSeam->Y() ),
2563 0.5 * ( n->Z() + nSeam->Z() ));
2565 getClosestUV( projector, pMid, uv );
2566 if ( uv.Coord( iPar ) > uvMap[ n ]->Coord( iPar ) ) {
2567 nodesNearSeam.insert( n );
2573 // for centroidalSmooth all element nodes must
2574 // be on one side of a seam
2575 if ( theSmoothMethod == CENTROIDAL && nbUseMap1 && nbUseMap2 ) {
2576 SMDS_ElemIteratorPtr nIt = e->nodesIterator();
2578 while ( nn++ < nbn ) {
2579 const SMDS_MeshNode* n =
2580 static_cast<const SMDS_MeshNode*>( nIt->next() );
2581 setMovableNodes.erase( n );
2585 } // loop on nodes on seam
2586 } // loop on edge of a face
2587 } // if ( !face.IsNull() )
2589 if ( setMovableNodes.empty() ) {
2590 MESSAGE( "Face id : " << *fId << " - NO SMOOTHING: no nodes to move!!!");
2591 continue; // goto next face
2599 double maxRatio = -1., maxDisplacement = -1.;
2600 set<const SMDS_MeshNode*>::iterator nodeToMove;
2601 for ( it = 0; it < theNbIterations; it++ ) {
2602 maxDisplacement = 0.;
2603 nodeToMove = setMovableNodes.begin();
2604 for ( ; nodeToMove != setMovableNodes.end(); nodeToMove++ ) {
2605 const SMDS_MeshNode* node = (*nodeToMove);
2606 gp_XYZ aPrevPos ( node->X(), node->Y(), node->Z() );
2609 bool map2 = ( nodesNearSeam.find( node ) != nodesNearSeam.end() );
2610 if ( theSmoothMethod == LAPLACIAN )
2611 laplacianSmooth( node, surface, map2 ? uvMap2 : uvMap );
2613 centroidalSmooth( node, surface, map2 ? uvMap2 : uvMap );
2615 // node displacement
2616 gp_XYZ aNewPos ( node->X(), node->Y(), node->Z() );
2617 Standard_Real aDispl = (aPrevPos - aNewPos).SquareModulus();
2618 if ( aDispl > maxDisplacement )
2619 maxDisplacement = aDispl;
2621 // no node movement => exit
2622 //if ( maxDisplacement < 1.e-16 ) {
2623 if ( maxDisplacement < disttol ) {
2624 MESSAGE("-- no node movement --");
2628 // check elements quality
2630 list< const SMDS_MeshElement* >::iterator elemIt = elemsOnFace.begin();
2631 for ( ; elemIt != elemsOnFace.end(); ++elemIt ) {
2632 const SMDS_MeshElement* elem = (*elemIt);
2633 if ( !elem || elem->GetType() != SMDSAbs_Face )
2635 SMESH::Controls::TSequenceOfXYZ aPoints;
2636 if ( aQualityFunc.GetPoints( elem, aPoints )) {
2637 double aValue = aQualityFunc.GetValue( aPoints );
2638 if ( aValue > maxRatio )
2642 if ( maxRatio <= theTgtAspectRatio ) {
2643 MESSAGE("-- quality achived --");
2646 if (it+1 == theNbIterations) {
2647 MESSAGE("-- Iteration limit exceeded --");
2649 } // smoothing iterations
2651 MESSAGE(" Face id: " << *fId <<
2652 " Nb iterstions: " << it <<
2653 " Displacement: " << maxDisplacement <<
2654 " Aspect Ratio " << maxRatio);
2656 // ---------------------------------------
2657 // new nodes positions are computed,
2658 // record movement in DS and set new UV
2659 // ---------------------------------------
2660 nodeToMove = setMovableNodes.begin();
2661 for ( ; nodeToMove != setMovableNodes.end(); nodeToMove++ ) {
2662 SMDS_MeshNode* node = const_cast< SMDS_MeshNode* > (*nodeToMove);
2663 aMesh->MoveNode( node, node->X(), node->Y(), node->Z() );
2664 map< const SMDS_MeshNode*, gp_XY* >::iterator node_uv = uvMap.find( node );
2665 if ( node_uv != uvMap.end() ) {
2666 gp_XY* uv = node_uv->second;
2668 ( SMDS_PositionPtr( new SMDS_FacePosition( *fId, uv->X(), uv->Y() )));
2672 // move medium nodes of quadratic elements
2675 SMESH_MesherHelper helper( *GetMesh() );
2676 if ( !face.IsNull() )
2677 helper.SetSubShape( face );
2678 list< const SMDS_MeshElement* >::iterator elemIt = elemsOnFace.begin();
2679 for ( ; elemIt != elemsOnFace.end(); ++elemIt ) {
2680 const SMDS_QuadraticFaceOfNodes* QF =
2681 dynamic_cast<const SMDS_QuadraticFaceOfNodes*> (*elemIt);
2683 vector<const SMDS_MeshNode*> Ns;
2684 Ns.reserve(QF->NbNodes()+1);
2685 SMDS_NodeIteratorPtr anIter = QF->interlacedNodesIterator();
2686 while ( anIter->more() )
2687 Ns.push_back( anIter->next() );
2688 Ns.push_back( Ns[0] );
2690 for(int i=0; i<QF->NbNodes(); i=i+2) {
2691 if ( !surface.IsNull() ) {
2692 gp_XY uv1 = helper.GetNodeUV( face, Ns[i], Ns[i+2] );
2693 gp_XY uv2 = helper.GetNodeUV( face, Ns[i+2], Ns[i] );
2694 gp_XY uv = ( uv1 + uv2 ) / 2.;
2695 gp_Pnt xyz = surface->Value( uv.X(), uv.Y() );
2696 x = xyz.X(); y = xyz.Y(); z = xyz.Z();
2699 x = (Ns[i]->X() + Ns[i+2]->X())/2;
2700 y = (Ns[i]->Y() + Ns[i+2]->Y())/2;
2701 z = (Ns[i]->Z() + Ns[i+2]->Z())/2;
2703 if( fabs( Ns[i+1]->X() - x ) > disttol ||
2704 fabs( Ns[i+1]->Y() - y ) > disttol ||
2705 fabs( Ns[i+1]->Z() - z ) > disttol ) {
2706 // we have to move i+1 node
2707 aMesh->MoveNode( Ns[i+1], x, y, z );
2714 } // loop on face ids
2718 //=======================================================================
2719 //function : isReverse
2720 //purpose : Return true if normal of prevNodes is not co-directied with
2721 // gp_Vec(prevNodes[iNotSame],nextNodes[iNotSame]).
2722 // iNotSame is where prevNodes and nextNodes are different
2723 //=======================================================================
2725 static bool isReverse(vector<const SMDS_MeshNode*> prevNodes,
2726 vector<const SMDS_MeshNode*> nextNodes,
2730 int iBeforeNotSame = ( iNotSame == 0 ? nbNodes - 1 : iNotSame - 1 );
2731 int iAfterNotSame = ( iNotSame + 1 == nbNodes ? 0 : iNotSame + 1 );
2733 const SMDS_MeshNode* nB = prevNodes[ iBeforeNotSame ];
2734 const SMDS_MeshNode* nA = prevNodes[ iAfterNotSame ];
2735 const SMDS_MeshNode* nP = prevNodes[ iNotSame ];
2736 const SMDS_MeshNode* nN = nextNodes[ iNotSame ];
2738 gp_Pnt pB ( nB->X(), nB->Y(), nB->Z() );
2739 gp_Pnt pA ( nA->X(), nA->Y(), nA->Z() );
2740 gp_Pnt pP ( nP->X(), nP->Y(), nP->Z() );
2741 gp_Pnt pN ( nN->X(), nN->Y(), nN->Z() );
2743 gp_Vec vB ( pP, pB ), vA ( pP, pA ), vN ( pP, pN );
2745 return (vA ^ vB) * vN < 0.0;
2748 //=======================================================================
2750 * \brief Create elements by sweeping an element
2751 * \param elem - element to sweep
2752 * \param newNodesItVec - nodes generated from each node of the element
2753 * \param newElems - generated elements
2754 * \param nbSteps - number of sweeping steps
2755 * \param srcElements - to append elem for each generated element
2757 //=======================================================================
2759 void SMESH_MeshEditor::sweepElement(const SMDS_MeshElement* elem,
2760 const vector<TNodeOfNodeListMapItr> & newNodesItVec,
2761 list<const SMDS_MeshElement*>& newElems,
2763 SMESH_SequenceOfElemPtr& srcElements)
2765 SMESHDS_Mesh* aMesh = GetMeshDS();
2767 // Loop on elem nodes:
2768 // find new nodes and detect same nodes indices
2769 int nbNodes = elem->NbNodes();
2770 vector < list< const SMDS_MeshNode* >::const_iterator > itNN( nbNodes );
2771 vector<const SMDS_MeshNode*> prevNod( nbNodes );
2772 vector<const SMDS_MeshNode*> nextNod( nbNodes );
2773 vector<const SMDS_MeshNode*> midlNod( nbNodes );
2775 int iNode, nbSame = 0, iNotSameNode = 0, iSameNode = 0;
2776 vector<int> sames(nbNodes);
2778 //bool issimple[nbNodes];
2779 vector<bool> issimple(nbNodes);
2781 for ( iNode = 0; iNode < nbNodes; iNode++ ) {
2782 TNodeOfNodeListMapItr nnIt = newNodesItVec[ iNode ];
2783 const SMDS_MeshNode* node = nnIt->first;
2784 const list< const SMDS_MeshNode* > & listNewNodes = nnIt->second;
2785 if ( listNewNodes.empty() )
2788 if(listNewNodes.size()==nbSteps) {
2789 issimple[iNode] = true;
2792 issimple[iNode] = false;
2795 itNN[ iNode ] = listNewNodes.begin();
2796 prevNod[ iNode ] = node;
2797 nextNod[ iNode ] = listNewNodes.front();
2798 //cout<<"iNode="<<iNode<<endl;
2799 //cout<<" prevNod[iNode]="<< prevNod[iNode]<<" nextNod[iNode]="<< nextNod[iNode]<<endl;
2800 if ( prevNod[ iNode ] != nextNod [ iNode ])
2801 iNotSameNode = iNode;
2805 sames[nbSame++] = iNode;
2808 //cout<<"1 nbSame="<<nbSame<<endl;
2809 if ( nbSame == nbNodes || nbSame > 2) {
2810 MESSAGE( " Too many same nodes of element " << elem->GetID() );
2814 // if( elem->IsQuadratic() && nbSame>0 ) {
2815 // MESSAGE( "Can not rotate quadratic element " << elem->GetID() );
2819 int iBeforeSame = 0, iAfterSame = 0, iOpposSame = 0;
2821 iBeforeSame = ( iSameNode == 0 ? nbNodes - 1 : iSameNode - 1 );
2822 iAfterSame = ( iSameNode + 1 == nbNodes ? 0 : iSameNode + 1 );
2823 iOpposSame = ( iSameNode - 2 < 0 ? iSameNode + 2 : iSameNode - 2 );
2827 //cout<<" prevNod[0]="<< prevNod[0]<<" prevNod[1]="<< prevNod[1]
2828 // <<" prevNod[2]="<< prevNod[2]<<" prevNod[3]="<< prevNod[4]
2829 // <<" prevNod[4]="<< prevNod[4]<<" prevNod[5]="<< prevNod[5]
2830 // <<" prevNod[6]="<< prevNod[6]<<" prevNod[7]="<< prevNod[7]<<endl;
2832 // check element orientation
2834 if ( nbNodes > 2 && !isReverse( prevNod, nextNod, nbNodes, iNotSameNode )) {
2835 //MESSAGE("Reversed elem " << elem );
2839 int iAB = iAfterSame + iBeforeSame;
2840 iBeforeSame = iAB - iBeforeSame;
2841 iAfterSame = iAB - iAfterSame;
2845 // make new elements
2846 for (int iStep = 0; iStep < nbSteps; iStep++ ) {
2848 for ( iNode = 0; iNode < nbNodes; iNode++ ) {
2849 if(issimple[iNode]) {
2850 nextNod[ iNode ] = *itNN[ iNode ];
2854 if( elem->GetType()==SMDSAbs_Node ) {
2855 // we have to use two nodes
2856 midlNod[ iNode ] = *itNN[ iNode ];
2858 nextNod[ iNode ] = *itNN[ iNode ];
2861 else if(!elem->IsQuadratic() ||
2862 elem->IsQuadratic() && elem->IsMediumNode(prevNod[iNode]) ) {
2863 // we have to use each second node
2865 nextNod[ iNode ] = *itNN[ iNode ];
2869 // we have to use two nodes
2870 midlNod[ iNode ] = *itNN[ iNode ];
2872 nextNod[ iNode ] = *itNN[ iNode ];
2877 SMDS_MeshElement* aNewElem = 0;
2878 if(!elem->IsPoly()) {
2879 switch ( nbNodes ) {
2883 if ( nbSame == 0 ) {
2885 aNewElem = aMesh->AddEdge( prevNod[ 0 ], nextNod[ 0 ] );
2887 aNewElem = aMesh->AddEdge( prevNod[ 0 ], nextNod[ 0 ], midlNod[ 0 ] );
2893 aNewElem = aMesh->AddFace(prevNod[ 0 ], prevNod[ 1 ],
2894 nextNod[ 1 ], nextNod[ 0 ] );
2896 aNewElem = aMesh->AddFace(prevNod[ 0 ], prevNod[ 1 ],
2897 nextNod[ iNotSameNode ] );
2901 case 3: { // TRIANGLE or quadratic edge
2902 if(elem->GetType() == SMDSAbs_Face) { // TRIANGLE
2904 if ( nbSame == 0 ) // --- pentahedron
2905 aNewElem = aMesh->AddVolume (prevNod[ i0 ], prevNod[ 1 ], prevNod[ i2 ],
2906 nextNod[ i0 ], nextNod[ 1 ], nextNod[ i2 ] );
2908 else if ( nbSame == 1 ) // --- pyramid
2909 aNewElem = aMesh->AddVolume (prevNod[ iBeforeSame ], prevNod[ iAfterSame ],
2910 nextNod[ iAfterSame ], nextNod[ iBeforeSame ],
2911 nextNod[ iSameNode ]);
2913 else // 2 same nodes: --- tetrahedron
2914 aNewElem = aMesh->AddVolume (prevNod[ i0 ], prevNod[ 1 ], prevNod[ i2 ],
2915 nextNod[ iNotSameNode ]);
2917 else { // quadratic edge
2918 if(nbSame==0) { // quadratic quadrangle
2919 aNewElem = aMesh->AddFace(prevNod[0], nextNod[0], nextNod[1], prevNod[1],
2920 midlNod[0], nextNod[2], midlNod[1], prevNod[2]);
2922 else if(nbSame==1) { // quadratic triangle
2924 return; // medium node on axis
2925 else if(sames[0]==0) {
2926 aNewElem = aMesh->AddFace(prevNod[0], nextNod[1], prevNod[1],
2927 nextNod[2], midlNod[1], prevNod[2]);
2929 else { // sames[0]==1
2930 aNewElem = aMesh->AddFace(prevNod[0], nextNod[0], prevNod[1],
2931 midlNod[0], nextNod[2], prevNod[2]);
2939 case 4: { // QUADRANGLE
2941 if ( nbSame == 0 ) // --- hexahedron
2942 aNewElem = aMesh->AddVolume (prevNod[ i0 ], prevNod[ 1 ], prevNod[ i2 ], prevNod[ 3 ],
2943 nextNod[ i0 ], nextNod[ 1 ], nextNod[ i2 ], nextNod[ 3 ]);
2945 else if ( nbSame == 1 ) { // --- pyramid + pentahedron
2946 aNewElem = aMesh->AddVolume (prevNod[ iBeforeSame ], prevNod[ iAfterSame ],
2947 nextNod[ iAfterSame ], nextNod[ iBeforeSame ],
2948 nextNod[ iSameNode ]);
2949 newElems.push_back( aNewElem );
2950 aNewElem = aMesh->AddVolume (prevNod[ iAfterSame ], prevNod[ iOpposSame ],
2951 prevNod[ iBeforeSame ], nextNod[ iAfterSame ],
2952 nextNod[ iOpposSame ], nextNod[ iBeforeSame ] );
2954 else if ( nbSame == 2 ) { // pentahedron
2955 if ( prevNod[ iBeforeSame ] == nextNod[ iBeforeSame ] )
2956 // iBeforeSame is same too
2957 aNewElem = aMesh->AddVolume (prevNod[ iBeforeSame ], prevNod[ iOpposSame ],
2958 nextNod[ iOpposSame ], prevNod[ iSameNode ],
2959 prevNod[ iAfterSame ], nextNod[ iAfterSame ]);
2961 // iAfterSame is same too
2962 aNewElem = aMesh->AddVolume (prevNod[ iSameNode ], prevNod[ iBeforeSame ],
2963 nextNod[ iBeforeSame ], prevNod[ iAfterSame ],
2964 prevNod[ iOpposSame ], nextNod[ iOpposSame ]);
2968 case 6: { // quadratic triangle
2969 // create pentahedron with 15 nodes
2970 if(i0>0) { // reversed case
2971 aNewElem = aMesh->AddVolume (prevNod[0], prevNod[2], prevNod[1],
2972 nextNod[0], nextNod[2], nextNod[1],
2973 prevNod[5], prevNod[4], prevNod[3],
2974 nextNod[5], nextNod[4], nextNod[3],
2975 midlNod[0], midlNod[2], midlNod[1]);
2977 else { // not reversed case
2978 aNewElem = aMesh->AddVolume (prevNod[0], prevNod[1], prevNod[2],
2979 nextNod[0], nextNod[1], nextNod[2],
2980 prevNod[3], prevNod[4], prevNod[5],
2981 nextNod[3], nextNod[4], nextNod[5],
2982 midlNod[0], midlNod[1], midlNod[2]);
2986 case 8: { // quadratic quadrangle
2987 // create hexahedron with 20 nodes
2988 if(i0>0) { // reversed case
2989 aNewElem = aMesh->AddVolume (prevNod[0], prevNod[3], prevNod[2], prevNod[1],
2990 nextNod[0], nextNod[3], nextNod[2], nextNod[1],
2991 prevNod[7], prevNod[6], prevNod[5], prevNod[4],
2992 nextNod[7], nextNod[6], nextNod[5], nextNod[4],
2993 midlNod[0], midlNod[3], midlNod[2], midlNod[1]);
2995 else { // not reversed case
2996 aNewElem = aMesh->AddVolume (prevNod[0], prevNod[1], prevNod[2], prevNod[3],
2997 nextNod[0], nextNod[1], nextNod[2], nextNod[3],
2998 prevNod[4], prevNod[5], prevNod[6], prevNod[7],
2999 nextNod[4], nextNod[5], nextNod[6], nextNod[7],
3000 midlNod[0], midlNod[1], midlNod[2], midlNod[3]);
3005 // realized for extrusion only
3006 //vector<const SMDS_MeshNode*> polyedre_nodes (nbNodes*2 + 4*nbNodes);
3007 //vector<int> quantities (nbNodes + 2);
3009 //quantities[0] = nbNodes; // bottom of prism
3010 //for (int inode = 0; inode < nbNodes; inode++) {
3011 // polyedre_nodes[inode] = prevNod[inode];
3014 //quantities[1] = nbNodes; // top of prism
3015 //for (int inode = 0; inode < nbNodes; inode++) {
3016 // polyedre_nodes[nbNodes + inode] = nextNod[inode];
3019 //for (int iface = 0; iface < nbNodes; iface++) {
3020 // quantities[iface + 2] = 4;
3021 // int inextface = (iface == nbNodes - 1) ? 0 : iface + 1;
3022 // polyedre_nodes[2*nbNodes + 4*iface + 0] = prevNod[iface];
3023 // polyedre_nodes[2*nbNodes + 4*iface + 1] = prevNod[inextface];
3024 // polyedre_nodes[2*nbNodes + 4*iface + 2] = nextNod[inextface];
3025 // polyedre_nodes[2*nbNodes + 4*iface + 3] = nextNod[iface];
3027 //aNewElem = aMesh->AddPolyhedralVolume (polyedre_nodes, quantities);
3034 // realized for extrusion only
3035 vector<const SMDS_MeshNode*> polyedre_nodes (nbNodes*2 + 4*nbNodes);
3036 vector<int> quantities (nbNodes + 2);
3038 quantities[0] = nbNodes; // bottom of prism
3039 for (int inode = 0; inode < nbNodes; inode++) {
3040 polyedre_nodes[inode] = prevNod[inode];
3043 quantities[1] = nbNodes; // top of prism
3044 for (int inode = 0; inode < nbNodes; inode++) {
3045 polyedre_nodes[nbNodes + inode] = nextNod[inode];
3048 for (int iface = 0; iface < nbNodes; iface++) {
3049 quantities[iface + 2] = 4;
3050 int inextface = (iface == nbNodes - 1) ? 0 : iface + 1;
3051 polyedre_nodes[2*nbNodes + 4*iface + 0] = prevNod[iface];
3052 polyedre_nodes[2*nbNodes + 4*iface + 1] = prevNod[inextface];
3053 polyedre_nodes[2*nbNodes + 4*iface + 2] = nextNod[inextface];
3054 polyedre_nodes[2*nbNodes + 4*iface + 3] = nextNod[iface];
3056 aNewElem = aMesh->AddPolyhedralVolume (polyedre_nodes, quantities);
3060 newElems.push_back( aNewElem );
3061 myLastCreatedElems.Append(aNewElem);
3062 srcElements.Append( elem );
3065 // set new prev nodes
3066 for ( iNode = 0; iNode < nbNodes; iNode++ )
3067 prevNod[ iNode ] = nextNod[ iNode ];
3072 //=======================================================================
3074 * \brief Create 1D and 2D elements around swept elements
3075 * \param mapNewNodes - source nodes and ones generated from them
3076 * \param newElemsMap - source elements and ones generated from them
3077 * \param elemNewNodesMap - nodes generated from each node of each element
3078 * \param elemSet - all swept elements
3079 * \param nbSteps - number of sweeping steps
3080 * \param srcElements - to append elem for each generated element
3082 //=======================================================================
3084 void SMESH_MeshEditor::makeWalls (TNodeOfNodeListMap & mapNewNodes,
3085 TElemOfElemListMap & newElemsMap,
3086 TElemOfVecOfNnlmiMap & elemNewNodesMap,
3087 TIDSortedElemSet& elemSet,
3089 SMESH_SequenceOfElemPtr& srcElements)
3091 ASSERT( newElemsMap.size() == elemNewNodesMap.size() );
3092 SMESHDS_Mesh* aMesh = GetMeshDS();
3094 // Find nodes belonging to only one initial element - sweep them to get edges.
3096 TNodeOfNodeListMapItr nList = mapNewNodes.begin();
3097 for ( ; nList != mapNewNodes.end(); nList++ ) {
3098 const SMDS_MeshNode* node =
3099 static_cast<const SMDS_MeshNode*>( nList->first );
3100 SMDS_ElemIteratorPtr eIt = node->GetInverseElementIterator();
3101 int nbInitElems = 0;
3102 const SMDS_MeshElement* el = 0;
3103 SMDSAbs_ElementType highType = SMDSAbs_Edge; // count most complex elements only
3104 while ( eIt->more() && nbInitElems < 2 ) {
3106 SMDSAbs_ElementType type = el->GetType();
3107 if ( type == SMDSAbs_Volume || type < highType ) continue;
3108 if ( type > highType ) {
3112 if ( elemSet.find(el) != elemSet.end() )
3115 if ( nbInitElems < 2 ) {
3116 bool NotCreateEdge = el && el->IsQuadratic() && el->IsMediumNode(node);
3117 if(!NotCreateEdge) {
3118 vector<TNodeOfNodeListMapItr> newNodesItVec( 1, nList );
3119 list<const SMDS_MeshElement*> newEdges;
3120 sweepElement( node, newNodesItVec, newEdges, nbSteps, srcElements );
3125 // Make a ceiling for each element ie an equal element of last new nodes.
3126 // Find free links of faces - make edges and sweep them into faces.
3128 TElemOfElemListMap::iterator itElem = newElemsMap.begin();
3129 TElemOfVecOfNnlmiMap::iterator itElemNodes = elemNewNodesMap.begin();
3130 for ( ; itElem != newElemsMap.end(); itElem++, itElemNodes++ ) {
3131 const SMDS_MeshElement* elem = itElem->first;
3132 vector<TNodeOfNodeListMapItr>& vecNewNodes = itElemNodes->second;
3134 if ( elem->GetType() == SMDSAbs_Edge ) {
3135 // create a ceiling edge
3136 if (!elem->IsQuadratic()) {
3137 if ( !aMesh->FindEdge( vecNewNodes[ 0 ]->second.back(),
3138 vecNewNodes[ 1 ]->second.back())) {
3139 myLastCreatedElems.Append(aMesh->AddEdge(vecNewNodes[ 0 ]->second.back(),
3140 vecNewNodes[ 1 ]->second.back()));
3141 srcElements.Append( myLastCreatedElems.Last() );
3145 if ( !aMesh->FindEdge( vecNewNodes[ 0 ]->second.back(),
3146 vecNewNodes[ 1 ]->second.back(),
3147 vecNewNodes[ 2 ]->second.back())) {
3148 myLastCreatedElems.Append(aMesh->AddEdge(vecNewNodes[ 0 ]->second.back(),
3149 vecNewNodes[ 1 ]->second.back(),
3150 vecNewNodes[ 2 ]->second.back()));
3151 srcElements.Append( myLastCreatedElems.Last() );
3155 if ( elem->GetType() != SMDSAbs_Face )
3158 if(itElem->second.size()==0) continue;
3160 bool hasFreeLinks = false;
3162 TIDSortedElemSet avoidSet;
3163 avoidSet.insert( elem );
3165 set<const SMDS_MeshNode*> aFaceLastNodes;
3166 int iNode, nbNodes = vecNewNodes.size();
3167 if(!elem->IsQuadratic()) {
3168 // loop on the face nodes
3169 for ( iNode = 0; iNode < nbNodes; iNode++ ) {
3170 aFaceLastNodes.insert( vecNewNodes[ iNode ]->second.back() );
3171 // look for free links of the face
3172 int iNext = ( iNode + 1 == nbNodes ) ? 0 : iNode + 1;
3173 const SMDS_MeshNode* n1 = vecNewNodes[ iNode ]->first;
3174 const SMDS_MeshNode* n2 = vecNewNodes[ iNext ]->first;
3175 // check if a link is free
3176 if ( ! SMESH_MeshEditor::FindFaceInSet ( n1, n2, elemSet, avoidSet )) {
3177 hasFreeLinks = true;
3178 // make an edge and a ceiling for a new edge
3179 if ( !aMesh->FindEdge( n1, n2 )) {
3180 myLastCreatedElems.Append(aMesh->AddEdge( n1, n2 )); // free link edge
3181 srcElements.Append( myLastCreatedElems.Last() );
3183 n1 = vecNewNodes[ iNode ]->second.back();
3184 n2 = vecNewNodes[ iNext ]->second.back();
3185 if ( !aMesh->FindEdge( n1, n2 )) {
3186 myLastCreatedElems.Append(aMesh->AddEdge( n1, n2 )); // ceiling edge
3187 srcElements.Append( myLastCreatedElems.Last() );
3192 else { // elem is quadratic face
3193 int nbn = nbNodes/2;
3194 for ( iNode = 0; iNode < nbn; iNode++ ) {
3195 aFaceLastNodes.insert( vecNewNodes[ iNode ]->second.back() );
3196 int iNext = ( iNode + 1 == nbn ) ? 0 : iNode + 1;
3197 const SMDS_MeshNode* n1 = vecNewNodes[ iNode ]->first;
3198 const SMDS_MeshNode* n2 = vecNewNodes[ iNext ]->first;
3199 // check if a link is free
3200 if ( ! SMESH_MeshEditor::FindFaceInSet ( n1, n2, elemSet, avoidSet )) {
3201 hasFreeLinks = true;
3202 // make an edge and a ceiling for a new edge
3204 const SMDS_MeshNode* n3 = vecNewNodes[ iNode+nbn ]->first;
3205 if ( !aMesh->FindEdge( n1, n2, n3 )) {
3206 myLastCreatedElems.Append(aMesh->AddEdge( n1, n2, n3 )); // free link edge
3207 srcElements.Append( myLastCreatedElems.Last() );
3209 n1 = vecNewNodes[ iNode ]->second.back();
3210 n2 = vecNewNodes[ iNext ]->second.back();
3211 n3 = vecNewNodes[ iNode+nbn ]->second.back();
3212 if ( !aMesh->FindEdge( n1, n2, n3 )) {
3213 myLastCreatedElems.Append(aMesh->AddEdge( n1, n2, n3 )); // ceiling edge
3214 srcElements.Append( myLastCreatedElems.Last() );
3218 for ( iNode = nbn; iNode < 2*nbn; iNode++ ) {
3219 aFaceLastNodes.insert( vecNewNodes[ iNode ]->second.back() );
3223 // sweep free links into faces
3225 if ( hasFreeLinks ) {
3226 list<const SMDS_MeshElement*> & newVolumes = itElem->second;
3227 int iVol, volNb, nbVolumesByStep = newVolumes.size() / nbSteps;
3229 set<const SMDS_MeshNode*> initNodeSet, topNodeSet, faceNodeSet;
3230 for ( iNode = 0; iNode < nbNodes; iNode++ ) {
3231 initNodeSet.insert( vecNewNodes[ iNode ]->first );
3232 topNodeSet .insert( vecNewNodes[ iNode ]->second.back() );
3234 for ( volNb = 0; volNb < nbVolumesByStep; volNb++ ) {
3235 list<const SMDS_MeshElement*>::iterator v = newVolumes.begin();
3237 while ( iVol++ < volNb ) v++;
3238 // find indices of free faces of a volume and their source edges
3239 list< int > freeInd;
3240 list< const SMDS_MeshElement* > srcEdges; // source edges of free faces
3241 SMDS_VolumeTool vTool( *v );
3242 int iF, nbF = vTool.NbFaces();
3243 for ( iF = 0; iF < nbF; iF ++ ) {
3244 if (vTool.IsFreeFace( iF ) &&
3245 vTool.GetFaceNodes( iF, faceNodeSet ) &&
3246 initNodeSet != faceNodeSet) // except an initial face
3248 if ( nbSteps == 1 && faceNodeSet == topNodeSet )
3250 freeInd.push_back( iF );
3251 // find source edge of a free face iF
3252 vector<const SMDS_MeshNode*> commonNodes; // shared by the initial and free faces
3253 commonNodes.resize( initNodeSet.size(), NULL ); // avoid spoiling memory
3254 std::set_intersection( faceNodeSet.begin(), faceNodeSet.end(),
3255 initNodeSet.begin(), initNodeSet.end(),
3256 commonNodes.begin());
3257 if ( (*v)->IsQuadratic() )
3258 srcEdges.push_back(aMesh->FindEdge (commonNodes[0],commonNodes[1],commonNodes[2]));
3260 srcEdges.push_back(aMesh->FindEdge (commonNodes[0],commonNodes[1]));
3262 if ( !srcEdges.back() )
3264 cout << "SMESH_MeshEditor::makeWalls(), no source edge found for a free face #"
3265 << iF << " of volume #" << vTool.ID() << endl;
3270 if ( freeInd.empty() )
3273 // create faces for all steps;
3274 // if such a face has been already created by sweep of edge,
3275 // assure that its orientation is OK
3276 for ( int iStep = 0; iStep < nbSteps; iStep++ ) {
3278 vTool.SetExternalNormal();
3279 list< int >::iterator ind = freeInd.begin();
3280 list< const SMDS_MeshElement* >::iterator srcEdge = srcEdges.begin();
3281 for ( ; ind != freeInd.end(); ++ind, ++srcEdge ) // loop on free faces
3283 const SMDS_MeshNode** nodes = vTool.GetFaceNodes( *ind );
3284 int nbn = vTool.NbFaceNodes( *ind );
3286 case 3: { ///// triangle
3287 const SMDS_MeshFace * f = aMesh->FindFace( nodes[ 0 ], nodes[ 1 ], nodes[ 2 ]);
3289 myLastCreatedElems.Append(aMesh->AddFace( nodes[ 0 ], nodes[ 1 ], nodes[ 2 ] ));
3290 else if ( nodes[ 1 ] != f->GetNode( f->GetNodeIndex( nodes[ 0 ] ) + 1 ))
3291 aMesh->ChangeElementNodes( f, nodes, nbn );
3294 case 4: { ///// quadrangle
3295 const SMDS_MeshFace * f = aMesh->FindFace( nodes[ 0 ], nodes[ 1 ], nodes[ 2 ], nodes[ 3 ]);
3297 myLastCreatedElems.Append(aMesh->AddFace( nodes[ 0 ], nodes[ 1 ], nodes[ 2 ], nodes[ 3 ] ));
3298 else if ( nodes[ 1 ] != f->GetNode( f->GetNodeIndex( nodes[ 0 ] ) + 1 ))
3299 aMesh->ChangeElementNodes( f, nodes, nbn );
3303 if( (*v)->IsQuadratic() ) {
3304 if(nbn==6) { /////// quadratic triangle
3305 const SMDS_MeshFace * f = aMesh->FindFace( nodes[0], nodes[2], nodes[4],
3306 nodes[1], nodes[3], nodes[5] );
3308 myLastCreatedElems.Append(aMesh->AddFace(nodes[0], nodes[2], nodes[4],
3309 nodes[1], nodes[3], nodes[5]));
3310 else if ( nodes[ 2 ] != f->GetNode( f->GetNodeIndex( nodes[ 0 ] ) + 1 ))
3311 aMesh->ChangeElementNodes( f, nodes, nbn );
3313 else { /////// quadratic quadrangle
3314 const SMDS_MeshFace * f = aMesh->FindFace( nodes[0], nodes[2], nodes[4], nodes[6],
3315 nodes[1], nodes[3], nodes[5], nodes[7] );
3317 myLastCreatedElems.Append(aMesh->AddFace(nodes[0], nodes[2], nodes[4], nodes[6],
3318 nodes[1], nodes[3], nodes[5], nodes[7]));
3319 else if ( nodes[ 2 ] != f->GetNode( f->GetNodeIndex( nodes[ 0 ] ) + 1 ))
3320 aMesh->ChangeElementNodes( f, nodes, nbn );
3323 else { //////// polygon
3324 vector<const SMDS_MeshNode*> polygon_nodes ( nodes, &nodes[nbn] );
3325 const SMDS_MeshFace * f = aMesh->FindFace( polygon_nodes );
3327 myLastCreatedElems.Append(aMesh->AddPolygonalFace(polygon_nodes));
3328 else if ( nodes[ 1 ] != f->GetNode( f->GetNodeIndex( nodes[ 0 ] ) + 1 ))
3329 aMesh->ChangeElementNodes( f, nodes, nbn );
3332 while ( srcElements.Length() < myLastCreatedElems.Length() )
3333 srcElements.Append( *srcEdge );
3335 } // loop on free faces
3337 // go to the next volume
3339 while ( iVol++ < nbVolumesByStep ) v++;
3342 } // sweep free links into faces
3344 // Make a ceiling face with a normal external to a volume
3346 SMDS_VolumeTool lastVol( itElem->second.back() );
3348 int iF = lastVol.GetFaceIndex( aFaceLastNodes );
3350 lastVol.SetExternalNormal();
3351 const SMDS_MeshNode** nodes = lastVol.GetFaceNodes( iF );
3352 int nbn = lastVol.NbFaceNodes( iF );
3355 if (!hasFreeLinks ||
3356 !aMesh->FindFace( nodes[ 0 ], nodes[ 1 ], nodes[ 2 ]))
3357 myLastCreatedElems.Append(aMesh->AddFace( nodes[ 0 ], nodes[ 1 ], nodes[ 2 ] ));
3360 if (!hasFreeLinks ||
3361 !aMesh->FindFace( nodes[ 0 ], nodes[ 1 ], nodes[ 2 ], nodes[ 3 ]))
3362 myLastCreatedElems.Append(aMesh->AddFace( nodes[ 0 ], nodes[ 1 ], nodes[ 2 ], nodes[ 3 ] ));
3365 if(itElem->second.back()->IsQuadratic()) {
3367 if (!hasFreeLinks ||
3368 !aMesh->FindFace(nodes[0], nodes[2], nodes[4],
3369 nodes[1], nodes[3], nodes[5]) ) {
3370 myLastCreatedElems.Append(aMesh->AddFace(nodes[0], nodes[2], nodes[4],
3371 nodes[1], nodes[3], nodes[5]));
3375 if (!hasFreeLinks ||
3376 !aMesh->FindFace(nodes[0], nodes[2], nodes[4], nodes[6],
3377 nodes[1], nodes[3], nodes[5], nodes[7]) )
3378 myLastCreatedElems.Append(aMesh->AddFace(nodes[0], nodes[2], nodes[4], nodes[6],
3379 nodes[1], nodes[3], nodes[5], nodes[7]));
3383 vector<const SMDS_MeshNode*> polygon_nodes ( nodes, &nodes[nbn] );
3384 if (!hasFreeLinks || !aMesh->FindFace(polygon_nodes))
3385 myLastCreatedElems.Append(aMesh->AddPolygonalFace(polygon_nodes));
3389 while ( srcElements.Length() < myLastCreatedElems.Length() )
3390 srcElements.Append( myLastCreatedElems.Last() );
3392 } // loop on swept elements
3395 //=======================================================================
3396 //function : RotationSweep
3398 //=======================================================================
3400 SMESH_MeshEditor::PGroupIDs
3401 SMESH_MeshEditor::RotationSweep(TIDSortedElemSet & theElems,
3402 const gp_Ax1& theAxis,
3403 const double theAngle,
3404 const int theNbSteps,
3405 const double theTol,
3406 const bool theMakeGroups,
3407 const bool theMakeWalls)
3409 myLastCreatedElems.Clear();
3410 myLastCreatedNodes.Clear();
3412 // source elements for each generated one
3413 SMESH_SequenceOfElemPtr srcElems, srcNodes;
3415 MESSAGE( "RotationSweep()");
3417 aTrsf.SetRotation( theAxis, theAngle );
3419 aTrsf2.SetRotation( theAxis, theAngle/2. );
3421 gp_Lin aLine( theAxis );
3422 double aSqTol = theTol * theTol;
3424 SMESHDS_Mesh* aMesh = GetMeshDS();
3426 TNodeOfNodeListMap mapNewNodes;
3427 TElemOfVecOfNnlmiMap mapElemNewNodes;
3428 TElemOfElemListMap newElemsMap;
3431 TIDSortedElemSet::iterator itElem;
3432 for ( itElem = theElems.begin(); itElem != theElems.end(); itElem++ ) {
3433 const SMDS_MeshElement* elem = *itElem;
3434 if ( !elem || elem->GetType() == SMDSAbs_Volume )
3436 vector<TNodeOfNodeListMapItr> & newNodesItVec = mapElemNewNodes[ elem ];
3437 newNodesItVec.reserve( elem->NbNodes() );
3439 // loop on elem nodes
3440 SMDS_ElemIteratorPtr itN = elem->nodesIterator();
3441 while ( itN->more() )
3443 // check if a node has been already sweeped
3444 const SMDS_MeshNode* node = cast2Node( itN->next() );
3445 TNodeOfNodeListMapItr nIt = mapNewNodes.find( node );
3446 if ( nIt == mapNewNodes.end() ) {
3447 nIt = mapNewNodes.insert( make_pair( node, list<const SMDS_MeshNode*>() )).first;
3448 list<const SMDS_MeshNode*>& listNewNodes = nIt->second;
3451 gp_XYZ aXYZ( node->X(), node->Y(), node->Z() );
3453 aXYZ.Coord( coord[0], coord[1], coord[2] );
3454 bool isOnAxis = ( aLine.SquareDistance( aXYZ ) <= aSqTol );
3455 const SMDS_MeshNode * newNode = node;
3456 for ( int i = 0; i < theNbSteps; i++ ) {
3458 if( elem->IsQuadratic() && !elem->IsMediumNode(node) ) {
3460 aTrsf2.Transforms( coord[0], coord[1], coord[2] );
3461 //aTrsf.Transforms( coord[0], coord[1], coord[2] );
3462 newNode = aMesh->AddNode( coord[0], coord[1], coord[2] );
3463 myLastCreatedNodes.Append(newNode);
3464 srcNodes.Append( node );
3465 listNewNodes.push_back( newNode );
3466 aTrsf2.Transforms( coord[0], coord[1], coord[2] );
3467 //aTrsf.Transforms( coord[0], coord[1], coord[2] );
3470 aTrsf.Transforms( coord[0], coord[1], coord[2] );
3472 newNode = aMesh->AddNode( coord[0], coord[1], coord[2] );
3473 myLastCreatedNodes.Append(newNode);
3474 srcNodes.Append( node );
3476 listNewNodes.push_back( newNode );
3480 // if current elem is quadratic and current node is not medium
3481 // we have to check - may be it is needed to insert additional nodes
3482 if( elem->IsQuadratic() && !elem->IsMediumNode(node) ) {
3483 list< const SMDS_MeshNode* > & listNewNodes = nIt->second;
3484 if(listNewNodes.size()==theNbSteps) {
3485 listNewNodes.clear();
3487 gp_XYZ aXYZ( node->X(), node->Y(), node->Z() );
3489 aXYZ.Coord( coord[0], coord[1], coord[2] );
3490 const SMDS_MeshNode * newNode = node;
3491 for(int i = 0; i<theNbSteps; i++) {
3492 aTrsf2.Transforms( coord[0], coord[1], coord[2] );
3493 newNode = aMesh->AddNode( coord[0], coord[1], coord[2] );
3494 myLastCreatedNodes.Append(newNode);
3495 listNewNodes.push_back( newNode );
3496 srcNodes.Append( node );
3497 aTrsf2.Transforms( coord[0], coord[1], coord[2] );
3498 newNode = aMesh->AddNode( coord[0], coord[1], coord[2] );
3499 myLastCreatedNodes.Append(newNode);
3500 srcNodes.Append( node );
3501 listNewNodes.push_back( newNode );
3506 newNodesItVec.push_back( nIt );
3508 // make new elements
3509 sweepElement( elem, newNodesItVec, newElemsMap[elem], theNbSteps, srcElems );
3513 makeWalls( mapNewNodes, newElemsMap, mapElemNewNodes, theElems, theNbSteps, srcElems );
3515 PGroupIDs newGroupIDs;
3516 if ( theMakeGroups )
3517 newGroupIDs = generateGroups( srcNodes, srcElems, "rotated");
3523 //=======================================================================
3524 //function : CreateNode
3526 //=======================================================================
3527 const SMDS_MeshNode* SMESH_MeshEditor::CreateNode(const double x,
3530 const double tolnode,
3531 SMESH_SequenceOfNode& aNodes)
3533 myLastCreatedElems.Clear();
3534 myLastCreatedNodes.Clear();
3537 SMESHDS_Mesh * aMesh = myMesh->GetMeshDS();
3539 // try to search in sequence of existing nodes
3540 // if aNodes.Length()>0 we 'nave to use given sequence
3541 // else - use all nodes of mesh
3542 if(aNodes.Length()>0) {
3544 for(i=1; i<=aNodes.Length(); i++) {
3545 gp_Pnt P2(aNodes.Value(i)->X(),aNodes.Value(i)->Y(),aNodes.Value(i)->Z());
3546 if(P1.Distance(P2)<tolnode)
3547 return aNodes.Value(i);
3551 SMDS_NodeIteratorPtr itn = aMesh->nodesIterator();
3552 while(itn->more()) {
3553 const SMDS_MeshNode* aN = static_cast<const SMDS_MeshNode*> (itn->next());
3554 gp_Pnt P2(aN->X(),aN->Y(),aN->Z());
3555 if(P1.Distance(P2)<tolnode)
3560 // create new node and return it
3561 const SMDS_MeshNode* NewNode = aMesh->AddNode(x,y,z);
3562 myLastCreatedNodes.Append(NewNode);
3567 //=======================================================================
3568 //function : ExtrusionSweep
3570 //=======================================================================
3572 SMESH_MeshEditor::PGroupIDs
3573 SMESH_MeshEditor::ExtrusionSweep (TIDSortedElemSet & theElems,
3574 const gp_Vec& theStep,
3575 const int theNbSteps,
3576 TElemOfElemListMap& newElemsMap,
3577 const bool theMakeGroups,
3579 const double theTolerance)
3581 ExtrusParam aParams;
3582 aParams.myDir = gp_Dir(theStep);
3583 aParams.myNodes.Clear();
3584 aParams.mySteps = new TColStd_HSequenceOfReal;
3586 for(i=1; i<=theNbSteps; i++)
3587 aParams.mySteps->Append(theStep.Magnitude());
3590 ExtrusionSweep(theElems,aParams,newElemsMap,theMakeGroups,theFlags,theTolerance);
3594 //=======================================================================
3595 //function : ExtrusionSweep
3597 //=======================================================================
3599 SMESH_MeshEditor::PGroupIDs
3600 SMESH_MeshEditor::ExtrusionSweep (TIDSortedElemSet & theElems,
3601 ExtrusParam& theParams,
3602 TElemOfElemListMap& newElemsMap,
3603 const bool theMakeGroups,
3605 const double theTolerance)
3607 myLastCreatedElems.Clear();
3608 myLastCreatedNodes.Clear();
3610 // source elements for each generated one
3611 SMESH_SequenceOfElemPtr srcElems, srcNodes;
3613 SMESHDS_Mesh* aMesh = GetMeshDS();
3615 int nbsteps = theParams.mySteps->Length();
3617 TNodeOfNodeListMap mapNewNodes;
3618 //TNodeOfNodeVecMap mapNewNodes;
3619 TElemOfVecOfNnlmiMap mapElemNewNodes;
3620 //TElemOfVecOfMapNodesMap mapElemNewNodes;
3623 TIDSortedElemSet::iterator itElem;
3624 for ( itElem = theElems.begin(); itElem != theElems.end(); itElem++ ) {
3625 // check element type
3626 const SMDS_MeshElement* elem = *itElem;
3627 if ( !elem || elem->GetType() == SMDSAbs_Volume )
3630 vector<TNodeOfNodeListMapItr> & newNodesItVec = mapElemNewNodes[ elem ];
3631 //vector<TNodeOfNodeVecMapItr> & newNodesItVec = mapElemNewNodes[ elem ];
3632 newNodesItVec.reserve( elem->NbNodes() );
3634 // loop on elem nodes
3635 SMDS_ElemIteratorPtr itN = elem->nodesIterator();
3636 while ( itN->more() )
3638 // check if a node has been already sweeped
3639 const SMDS_MeshNode* node = cast2Node( itN->next() );
3640 TNodeOfNodeListMap::iterator nIt = mapNewNodes.find( node );
3641 //TNodeOfNodeVecMap::iterator nIt = mapNewNodes.find( node );
3642 if ( nIt == mapNewNodes.end() ) {
3643 nIt = mapNewNodes.insert( make_pair( node, list<const SMDS_MeshNode*>() )).first;
3644 //nIt = mapNewNodes.insert( make_pair( node, vector<const SMDS_MeshNode*>() )).first;
3645 list<const SMDS_MeshNode*>& listNewNodes = nIt->second;
3646 //vector<const SMDS_MeshNode*>& vecNewNodes = nIt->second;
3647 //vecNewNodes.reserve(nbsteps);
3650 double coord[] = { node->X(), node->Y(), node->Z() };
3651 //int nbsteps = theParams.mySteps->Length();
3652 for ( int i = 0; i < nbsteps; i++ ) {
3653 if( elem->IsQuadratic() && !elem->IsMediumNode(node) ) {
3654 // create additional node
3655 double x = coord[0] + theParams.myDir.X()*theParams.mySteps->Value(i+1)/2.;
3656 double y = coord[1] + theParams.myDir.Y()*theParams.mySteps->Value(i+1)/2.;
3657 double z = coord[2] + theParams.myDir.Z()*theParams.mySteps->Value(i+1)/2.;
3658 if( theFlags & EXTRUSION_FLAG_SEW ) {
3659 const SMDS_MeshNode * newNode = CreateNode(x, y, z,
3660 theTolerance, theParams.myNodes);
3661 listNewNodes.push_back( newNode );
3664 const SMDS_MeshNode * newNode = aMesh->AddNode(x, y, z);
3665 myLastCreatedNodes.Append(newNode);
3666 srcNodes.Append( node );
3667 listNewNodes.push_back( newNode );
3670 //aTrsf.Transforms( coord[0], coord[1], coord[2] );
3671 coord[0] = coord[0] + theParams.myDir.X()*theParams.mySteps->Value(i+1);
3672 coord[1] = coord[1] + theParams.myDir.Y()*theParams.mySteps->Value(i+1);
3673 coord[2] = coord[2] + theParams.myDir.Z()*theParams.mySteps->Value(i+1);
3674 if( theFlags & EXTRUSION_FLAG_SEW ) {
3675 const SMDS_MeshNode * newNode = CreateNode(coord[0], coord[1], coord[2],
3676 theTolerance, theParams.myNodes);
3677 listNewNodes.push_back( newNode );
3678 //vecNewNodes[i]=newNode;
3681 const SMDS_MeshNode * newNode = aMesh->AddNode( coord[0], coord[1], coord[2] );
3682 myLastCreatedNodes.Append(newNode);
3683 srcNodes.Append( node );
3684 listNewNodes.push_back( newNode );
3685 //vecNewNodes[i]=newNode;
3690 // if current elem is quadratic and current node is not medium
3691 // we have to check - may be it is needed to insert additional nodes
3692 if( elem->IsQuadratic() && !elem->IsMediumNode(node) ) {
3693 list< const SMDS_MeshNode* > & listNewNodes = nIt->second;
3694 if(listNewNodes.size()==nbsteps) {
3695 listNewNodes.clear();
3696 double coord[] = { node->X(), node->Y(), node->Z() };
3697 for ( int i = 0; i < nbsteps; i++ ) {
3698 double x = coord[0] + theParams.myDir.X()*theParams.mySteps->Value(i+1);
3699 double y = coord[1] + theParams.myDir.Y()*theParams.mySteps->Value(i+1);
3700 double z = coord[2] + theParams.myDir.Z()*theParams.mySteps->Value(i+1);
3701 if( theFlags & EXTRUSION_FLAG_SEW ) {
3702 const SMDS_MeshNode * newNode = CreateNode(x, y, z,
3703 theTolerance, theParams.myNodes);
3704 listNewNodes.push_back( newNode );
3707 const SMDS_MeshNode * newNode = aMesh->AddNode(x, y, z);
3708 myLastCreatedNodes.Append(newNode);
3709 srcNodes.Append( node );
3710 listNewNodes.push_back( newNode );
3712 coord[0] = coord[0] + theParams.myDir.X()*theParams.mySteps->Value(i+1);
3713 coord[1] = coord[1] + theParams.myDir.Y()*theParams.mySteps->Value(i+1);
3714 coord[2] = coord[2] + theParams.myDir.Z()*theParams.mySteps->Value(i+1);
3715 if( theFlags & EXTRUSION_FLAG_SEW ) {
3716 const SMDS_MeshNode * newNode = CreateNode(coord[0], coord[1], coord[2],
3717 theTolerance, theParams.myNodes);
3718 listNewNodes.push_back( newNode );
3721 const SMDS_MeshNode * newNode = aMesh->AddNode( coord[0], coord[1], coord[2] );
3722 myLastCreatedNodes.Append(newNode);
3723 srcNodes.Append( node );
3724 listNewNodes.push_back( newNode );
3730 newNodesItVec.push_back( nIt );
3732 // make new elements
3733 sweepElement( elem, newNodesItVec, newElemsMap[elem], nbsteps, srcElems );
3736 if( theFlags & EXTRUSION_FLAG_BOUNDARY ) {
3737 makeWalls( mapNewNodes, newElemsMap, mapElemNewNodes, theElems, nbsteps, srcElems );
3739 PGroupIDs newGroupIDs;
3740 if ( theMakeGroups )
3741 newGroupIDs = generateGroups( srcNodes, srcElems, "extruded");
3747 //=======================================================================
3748 //class : SMESH_MeshEditor_PathPoint
3749 //purpose : auxiliary class
3750 //=======================================================================
3751 class SMESH_MeshEditor_PathPoint {
3753 SMESH_MeshEditor_PathPoint() {
3754 myPnt.SetCoord(99., 99., 99.);
3755 myTgt.SetCoord(1.,0.,0.);
3759 void SetPnt(const gp_Pnt& aP3D){
3762 void SetTangent(const gp_Dir& aTgt){
3765 void SetAngle(const double& aBeta){
3768 void SetParameter(const double& aPrm){
3771 const gp_Pnt& Pnt()const{
3774 const gp_Dir& Tangent()const{
3777 double Angle()const{
3780 double Parameter()const{
3791 //=======================================================================
3792 //function : ExtrusionAlongTrack
3794 //=======================================================================
3795 SMESH_MeshEditor::Extrusion_Error
3796 SMESH_MeshEditor::ExtrusionAlongTrack (TIDSortedElemSet & theElements,
3797 SMESH_subMesh* theTrack,
3798 const SMDS_MeshNode* theN1,
3799 const bool theHasAngles,
3800 list<double>& theAngles,
3801 const bool theHasRefPoint,
3802 const gp_Pnt& theRefPoint,
3803 const bool theMakeGroups)
3805 myLastCreatedElems.Clear();
3806 myLastCreatedNodes.Clear();
3808 // source elements for each generated one
3809 SMESH_SequenceOfElemPtr srcElems, srcNodes;
3811 int j, aNbTP, aNbE, aNb;
3812 double aT1, aT2, aT, aAngle, aX, aY, aZ;
3813 std::list<double> aPrms;
3814 std::list<double>::iterator aItD;
3815 TIDSortedElemSet::iterator itElem;
3817 Standard_Real aTx1, aTx2, aL2, aTolVec, aTolVec2;
3821 Handle(Geom_Curve) aC3D;
3822 TopoDS_Edge aTrackEdge;
3823 TopoDS_Vertex aV1, aV2;
3825 SMDS_ElemIteratorPtr aItE;
3826 SMDS_NodeIteratorPtr aItN;
3827 SMDSAbs_ElementType aTypeE;
3829 TNodeOfNodeListMap mapNewNodes;
3830 TElemOfVecOfNnlmiMap mapElemNewNodes;
3831 TElemOfElemListMap newElemsMap;
3834 aTolVec2=aTolVec*aTolVec;
3837 aNbE = theElements.size();
3840 return EXTR_NO_ELEMENTS;
3842 // 1.1 Track Pattern
3845 SMESHDS_SubMesh* pSubMeshDS=theTrack->GetSubMeshDS();
3847 aItE = pSubMeshDS->GetElements();
3848 while ( aItE->more() ) {
3849 const SMDS_MeshElement* pE = aItE->next();
3850 aTypeE = pE->GetType();
3851 // Pattern must contain links only
3852 if ( aTypeE != SMDSAbs_Edge )
3853 return EXTR_PATH_NOT_EDGE;
3856 const TopoDS_Shape& aS = theTrack->GetSubShape();
3857 // Sub shape for the Pattern must be an Edge
3858 if ( aS.ShapeType() != TopAbs_EDGE )
3859 return EXTR_BAD_PATH_SHAPE;
3861 aTrackEdge = TopoDS::Edge( aS );
3862 // the Edge must not be degenerated
3863 if ( BRep_Tool::Degenerated( aTrackEdge ) )
3864 return EXTR_BAD_PATH_SHAPE;
3866 TopExp::Vertices( aTrackEdge, aV1, aV2 );
3867 aT1=BRep_Tool::Parameter( aV1, aTrackEdge );
3868 aT2=BRep_Tool::Parameter( aV2, aTrackEdge );
3870 aItN = theTrack->GetFather()->GetSubMesh( aV1 )->GetSubMeshDS()->GetNodes();
3871 const SMDS_MeshNode* aN1 = aItN->next();
3873 aItN = theTrack->GetFather()->GetSubMesh( aV2 )->GetSubMeshDS()->GetNodes();
3874 const SMDS_MeshNode* aN2 = aItN->next();
3876 // starting node must be aN1 or aN2
3877 if ( !( aN1 == theN1 || aN2 == theN1 ) )
3878 return EXTR_BAD_STARTING_NODE;
3880 aNbTP = pSubMeshDS->NbNodes() + 2;
3883 vector<double> aAngles( aNbTP );
3885 for ( j=0; j < aNbTP; ++j ) {
3889 if ( theHasAngles ) {
3890 aItD = theAngles.begin();
3891 for ( j=1; (aItD != theAngles.end()) && (j<aNbTP); ++aItD, ++j ) {
3893 aAngles[j] = aAngle;
3897 // 2. Collect parameters on the track edge
3898 aPrms.push_back( aT1 );
3899 aPrms.push_back( aT2 );
3901 aItN = pSubMeshDS->GetNodes();
3902 while ( aItN->more() ) {
3903 const SMDS_MeshNode* pNode = aItN->next();
3904 const SMDS_EdgePosition* pEPos =
3905 static_cast<const SMDS_EdgePosition*>( pNode->GetPosition().get() );
3906 aT = pEPos->GetUParameter();
3907 aPrms.push_back( aT );
3912 if ( aN1 == theN1 ) {
3924 SMESH_MeshEditor_PathPoint aPP;
3925 vector<SMESH_MeshEditor_PathPoint> aPPs( aNbTP );
3927 aC3D = BRep_Tool::Curve( aTrackEdge, aTx1, aTx2 );
3929 aItD = aPrms.begin();
3930 for ( j=0; aItD != aPrms.end(); ++aItD, ++j ) {
3932 aC3D->D1( aT, aP3D, aVec );
3933 aL2 = aVec.SquareMagnitude();
3934 if ( aL2 < aTolVec2 )
3935 return EXTR_CANT_GET_TANGENT;
3937 gp_Dir aTgt( aVec );
3938 aAngle = aAngles[j];
3941 aPP.SetTangent( aTgt );
3942 aPP.SetAngle( aAngle );
3943 aPP.SetParameter( aT );
3947 // 3. Center of rotation aV0
3949 if ( !theHasRefPoint ) {
3951 aGC.SetCoord( 0.,0.,0. );
3953 itElem = theElements.begin();
3954 for ( ; itElem != theElements.end(); itElem++ ) {
3955 const SMDS_MeshElement* elem = *itElem;
3957 SMDS_ElemIteratorPtr itN = elem->nodesIterator();
3958 while ( itN->more() ) {
3959 const SMDS_MeshNode* node = static_cast<const SMDS_MeshNode*>( itN->next() );
3964 if ( mapNewNodes.find( node ) == mapNewNodes.end() ) {
3965 list<const SMDS_MeshNode*> aLNx;
3966 mapNewNodes[node] = aLNx;
3968 gp_XYZ aXYZ( aX, aY, aZ );
3976 } // if (!theHasRefPoint) {
3977 mapNewNodes.clear();
3979 // 4. Processing the elements
3980 SMESHDS_Mesh* aMesh = GetMeshDS();
3982 for ( itElem = theElements.begin(); itElem != theElements.end(); itElem++ ) {
3983 // check element type
3984 const SMDS_MeshElement* elem = *itElem;
3985 aTypeE = elem->GetType();
3986 if ( !elem || ( aTypeE != SMDSAbs_Face && aTypeE != SMDSAbs_Edge ) )
3989 vector<TNodeOfNodeListMapItr> & newNodesItVec = mapElemNewNodes[ elem ];
3990 newNodesItVec.reserve( elem->NbNodes() );
3992 // loop on elem nodes
3994 SMDS_ElemIteratorPtr itN = elem->nodesIterator();
3995 while ( itN->more() )
3998 // check if a node has been already processed
3999 const SMDS_MeshNode* node =
4000 static_cast<const SMDS_MeshNode*>( itN->next() );
4001 TNodeOfNodeListMap::iterator nIt = mapNewNodes.find( node );
4002 if ( nIt == mapNewNodes.end() ) {
4003 nIt = mapNewNodes.insert( make_pair( node, list<const SMDS_MeshNode*>() )).first;
4004 list<const SMDS_MeshNode*>& listNewNodes = nIt->second;
4007 aX = node->X(); aY = node->Y(); aZ = node->Z();
4009 Standard_Real aAngle1x, aAngleT1T0, aTolAng;
4010 gp_Pnt aP0x, aP1x, aPN0, aPN1, aV0x, aV1x;
4011 gp_Ax1 anAx1, anAxT1T0;
4012 gp_Dir aDT1x, aDT0x, aDT1T0;
4017 aPN0.SetCoord(aX, aY, aZ);
4019 const SMESH_MeshEditor_PathPoint& aPP0 = aPPs[0];
4021 aDT0x= aPP0.Tangent();
4023 for ( j = 1; j < aNbTP; ++j ) {
4024 const SMESH_MeshEditor_PathPoint& aPP1 = aPPs[j];
4026 aDT1x = aPP1.Tangent();
4027 aAngle1x = aPP1.Angle();
4029 gp_Trsf aTrsf, aTrsfRot, aTrsfRotT1T0;
4031 gp_Vec aV01x( aP0x, aP1x );
4032 aTrsf.SetTranslation( aV01x );
4035 aV1x = aV0x.Transformed( aTrsf );
4036 aPN1 = aPN0.Transformed( aTrsf );
4038 // rotation 1 [ T1,T0 ]
4039 aAngleT1T0=-aDT1x.Angle( aDT0x );
4040 if (fabs(aAngleT1T0) > aTolAng) {
4042 anAxT1T0.SetLocation( aV1x );
4043 anAxT1T0.SetDirection( aDT1T0 );
4044 aTrsfRotT1T0.SetRotation( anAxT1T0, aAngleT1T0 );
4046 aPN1 = aPN1.Transformed( aTrsfRotT1T0 );
4050 if ( theHasAngles ) {
4051 anAx1.SetLocation( aV1x );
4052 anAx1.SetDirection( aDT1x );
4053 aTrsfRot.SetRotation( anAx1, aAngle1x );
4055 aPN1 = aPN1.Transformed( aTrsfRot );
4059 if( elem->IsQuadratic() && !elem->IsMediumNode(node) ) {
4060 // create additional node
4061 double x = ( aPN1.X() + aPN0.X() )/2.;
4062 double y = ( aPN1.Y() + aPN0.Y() )/2.;
4063 double z = ( aPN1.Z() + aPN0.Z() )/2.;
4064 const SMDS_MeshNode* newNode = aMesh->AddNode(x,y,z);
4065 myLastCreatedNodes.Append(newNode);
4066 srcNodes.Append( node );
4067 listNewNodes.push_back( newNode );
4072 const SMDS_MeshNode* newNode = aMesh->AddNode( aX, aY, aZ );
4073 myLastCreatedNodes.Append(newNode);
4074 srcNodes.Append( node );
4075 listNewNodes.push_back( newNode );
4085 // if current elem is quadratic and current node is not medium
4086 // we have to check - may be it is needed to insert additional nodes
4087 if( elem->IsQuadratic() && !elem->IsMediumNode(node) ) {
4088 list< const SMDS_MeshNode* > & listNewNodes = nIt->second;
4089 if(listNewNodes.size()==aNbTP-1) {
4090 vector<const SMDS_MeshNode*> aNodes(2*(aNbTP-1));
4091 gp_XYZ P(node->X(), node->Y(), node->Z());
4092 list< const SMDS_MeshNode* >::iterator it = listNewNodes.begin();
4094 for(i=0; i<aNbTP-1; i++) {
4095 const SMDS_MeshNode* N = *it;
4096 double x = ( N->X() + P.X() )/2.;
4097 double y = ( N->Y() + P.Y() )/2.;
4098 double z = ( N->Z() + P.Z() )/2.;
4099 const SMDS_MeshNode* newN = aMesh->AddNode(x,y,z);
4100 srcNodes.Append( node );
4101 myLastCreatedNodes.Append(newN);
4104 P = gp_XYZ(N->X(),N->Y(),N->Z());
4106 listNewNodes.clear();
4107 for(i=0; i<2*(aNbTP-1); i++) {
4108 listNewNodes.push_back(aNodes[i]);
4114 newNodesItVec.push_back( nIt );
4116 // make new elements
4117 //sweepElement( aMesh, elem, newNodesItVec, newElemsMap[elem],
4118 // newNodesItVec[0]->second.size(), myLastCreatedElems );
4119 sweepElement( elem, newNodesItVec, newElemsMap[elem], aNbTP-1, srcElems );
4122 makeWalls( mapNewNodes, newElemsMap, mapElemNewNodes, theElements, aNbTP-1, srcElems );
4124 if ( theMakeGroups )
4125 generateGroups( srcNodes, srcElems, "extruded");
4130 //=======================================================================
4131 //function : Transform
4133 //=======================================================================
4135 SMESH_MeshEditor::PGroupIDs
4136 SMESH_MeshEditor::Transform (TIDSortedElemSet & theElems,
4137 const gp_Trsf& theTrsf,
4139 const bool theMakeGroups,
4140 SMESH_Mesh* theTargetMesh)
4142 myLastCreatedElems.Clear();
4143 myLastCreatedNodes.Clear();
4145 bool needReverse = false;
4146 string groupPostfix;
4147 switch ( theTrsf.Form() ) {
4152 groupPostfix = "mirrored";
4155 groupPostfix = "rotated";
4157 case gp_Translation:
4158 groupPostfix = "translated";
4161 groupPostfix = "scaled";
4164 needReverse = false;
4165 groupPostfix = "transformed";
4168 SMESH_MeshEditor targetMeshEditor( theTargetMesh );
4169 SMESHDS_Mesh* aTgtMesh = theTargetMesh ? theTargetMesh->GetMeshDS() : 0;
4170 SMESHDS_Mesh* aMesh = GetMeshDS();
4173 // map old node to new one
4174 TNodeNodeMap nodeMap;
4176 // elements sharing moved nodes; those of them which have all
4177 // nodes mirrored but are not in theElems are to be reversed
4178 TIDSortedElemSet inverseElemSet;
4180 // source elements for each generated one
4181 SMESH_SequenceOfElemPtr srcElems, srcNodes;
4184 TIDSortedElemSet::iterator itElem;
4185 for ( itElem = theElems.begin(); itElem != theElems.end(); itElem++ ) {
4186 const SMDS_MeshElement* elem = *itElem;
4190 // loop on elem nodes
4191 SMDS_ElemIteratorPtr itN = elem->nodesIterator();
4192 while ( itN->more() ) {
4194 // check if a node has been already transformed
4195 const SMDS_MeshNode* node = cast2Node( itN->next() );
4196 pair<TNodeNodeMap::iterator,bool> n2n_isnew =
4197 nodeMap.insert( make_pair ( node, node ));
4198 if ( !n2n_isnew.second )
4202 coord[0] = node->X();
4203 coord[1] = node->Y();
4204 coord[2] = node->Z();
4205 theTrsf.Transforms( coord[0], coord[1], coord[2] );
4206 if ( theTargetMesh ) {
4207 const SMDS_MeshNode * newNode = aTgtMesh->AddNode( coord[0], coord[1], coord[2] );
4208 n2n_isnew.first->second = newNode;
4209 myLastCreatedNodes.Append(newNode);
4210 srcNodes.Append( node );
4212 else if ( theCopy ) {
4213 const SMDS_MeshNode * newNode = aMesh->AddNode( coord[0], coord[1], coord[2] );
4214 n2n_isnew.first->second = newNode;
4215 myLastCreatedNodes.Append(newNode);
4216 srcNodes.Append( node );
4219 aMesh->MoveNode( node, coord[0], coord[1], coord[2] );
4220 // node position on shape becomes invalid
4221 const_cast< SMDS_MeshNode* > ( node )->SetPosition
4222 ( SMDS_SpacePosition::originSpacePosition() );
4225 // keep inverse elements
4226 if ( !theCopy && !theTargetMesh && needReverse ) {
4227 SMDS_ElemIteratorPtr invElemIt = node->GetInverseElementIterator();
4228 while ( invElemIt->more() ) {
4229 const SMDS_MeshElement* iel = invElemIt->next();
4230 inverseElemSet.insert( iel );
4236 // either create new elements or reverse mirrored ones
4237 if ( !theCopy && !needReverse && !theTargetMesh )
4240 TIDSortedElemSet::iterator invElemIt = inverseElemSet.begin();
4241 for ( ; invElemIt != inverseElemSet.end(); invElemIt++ )
4242 theElems.insert( *invElemIt );
4244 // replicate or reverse elements
4247 REV_TETRA = 0, // = nbNodes - 4
4248 REV_PYRAMID = 1, // = nbNodes - 4
4249 REV_PENTA = 2, // = nbNodes - 4
4251 REV_HEXA = 4, // = nbNodes - 4
4255 { 2, 1, 0, 3, 4, 0, 0, 0 }, // REV_TETRA
4256 { 2, 1, 0, 3, 4, 0, 0, 0 }, // REV_PYRAMID
4257 { 2, 1, 0, 5, 4, 3, 0, 0 }, // REV_PENTA
4258 { 2, 1, 0, 3, 0, 0, 0, 0 }, // REV_FACE
4259 { 2, 1, 0, 3, 6, 5, 4, 7 }, // REV_HEXA
4260 { 0, 1, 2, 3, 4, 5, 6, 7 } // FORWARD
4263 for ( itElem = theElems.begin(); itElem != theElems.end(); itElem++ )
4265 const SMDS_MeshElement* elem = *itElem;
4266 if ( !elem || elem->GetType() == SMDSAbs_Node )
4269 int nbNodes = elem->NbNodes();
4270 int elemType = elem->GetType();
4272 if (elem->IsPoly()) {
4273 // Polygon or Polyhedral Volume
4274 switch ( elemType ) {
4277 vector<const SMDS_MeshNode*> poly_nodes (nbNodes);
4279 SMDS_ElemIteratorPtr itN = elem->nodesIterator();
4280 while (itN->more()) {
4281 const SMDS_MeshNode* node =
4282 static_cast<const SMDS_MeshNode*>(itN->next());
4283 TNodeNodeMap::iterator nodeMapIt = nodeMap.find(node);
4284 if (nodeMapIt == nodeMap.end())
4285 break; // not all nodes transformed
4287 // reverse mirrored faces and volumes
4288 poly_nodes[nbNodes - iNode - 1] = (*nodeMapIt).second;
4290 poly_nodes[iNode] = (*nodeMapIt).second;
4294 if ( iNode != nbNodes )
4295 continue; // not all nodes transformed
4297 if ( theTargetMesh ) {
4298 myLastCreatedElems.Append(aTgtMesh->AddPolygonalFace(poly_nodes));
4299 srcElems.Append( elem );
4301 else if ( theCopy ) {
4302 myLastCreatedElems.Append(aMesh->AddPolygonalFace(poly_nodes));
4303 srcElems.Append( elem );
4306 aMesh->ChangePolygonNodes(elem, poly_nodes);
4310 case SMDSAbs_Volume:
4312 // ATTENTION: Reversing is not yet done!!!
4313 const SMDS_PolyhedralVolumeOfNodes* aPolyedre =
4314 dynamic_cast<const SMDS_PolyhedralVolumeOfNodes*>( elem );
4316 MESSAGE("Warning: bad volumic element");
4320 vector<const SMDS_MeshNode*> poly_nodes;
4321 vector<int> quantities;
4323 bool allTransformed = true;
4324 int nbFaces = aPolyedre->NbFaces();
4325 for (int iface = 1; iface <= nbFaces && allTransformed; iface++) {
4326 int nbFaceNodes = aPolyedre->NbFaceNodes(iface);
4327 for (int inode = 1; inode <= nbFaceNodes && allTransformed; inode++) {
4328 const SMDS_MeshNode* node = aPolyedre->GetFaceNode(iface, inode);
4329 TNodeNodeMap::iterator nodeMapIt = nodeMap.find(node);
4330 if (nodeMapIt == nodeMap.end()) {
4331 allTransformed = false; // not all nodes transformed
4333 poly_nodes.push_back((*nodeMapIt).second);
4336 quantities.push_back(nbFaceNodes);
4338 if ( !allTransformed )
4339 continue; // not all nodes transformed
4341 if ( theTargetMesh ) {
4342 myLastCreatedElems.Append(aTgtMesh->AddPolyhedralVolume(poly_nodes, quantities));
4343 srcElems.Append( elem );
4345 else if ( theCopy ) {
4346 myLastCreatedElems.Append(aMesh->AddPolyhedralVolume(poly_nodes, quantities));
4347 srcElems.Append( elem );
4350 aMesh->ChangePolyhedronNodes(elem, poly_nodes, quantities);
4360 int* i = index[ FORWARD ];
4361 if ( needReverse && nbNodes > 2) // reverse mirrored faces and volumes
4362 if ( elemType == SMDSAbs_Face )
4363 i = index[ REV_FACE ];
4365 i = index[ nbNodes - 4 ];
4367 if(elem->IsQuadratic()) {
4368 static int anIds[] = {0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19};
4371 if(nbNodes==3) { // quadratic edge
4372 static int anIds[] = {1,0,2};
4375 else if(nbNodes==6) { // quadratic triangle
4376 static int anIds[] = {0,2,1,5,4,3};
4379 else if(nbNodes==8) { // quadratic quadrangle
4380 static int anIds[] = {0,3,2,1,7,6,5,4};
4383 else if(nbNodes==10) { // quadratic tetrahedron of 10 nodes
4384 static int anIds[] = {0,2,1,3,6,5,4,7,9,8};
4387 else if(nbNodes==13) { // quadratic pyramid of 13 nodes
4388 static int anIds[] = {0,3,2,1,4,8,7,6,5,9,12,11,10};
4391 else if(nbNodes==15) { // quadratic pentahedron with 15 nodes
4392 static int anIds[] = {0,2,1,3,5,4,8,7,6,11,10,9,12,14,13};
4395 else { // nbNodes==20 - quadratic hexahedron with 20 nodes
4396 static int anIds[] = {0,3,2,1,4,7,6,5,11,10,9,8,15,14,13,12,16,19,18,17};
4402 // find transformed nodes
4403 vector<const SMDS_MeshNode*> nodes(nbNodes);
4405 SMDS_ElemIteratorPtr itN = elem->nodesIterator();
4406 while ( itN->more() ) {
4407 const SMDS_MeshNode* node =
4408 static_cast<const SMDS_MeshNode*>( itN->next() );
4409 TNodeNodeMap::iterator nodeMapIt = nodeMap.find( node );
4410 if ( nodeMapIt == nodeMap.end() )
4411 break; // not all nodes transformed
4412 nodes[ i [ iNode++ ]] = (*nodeMapIt).second;
4414 if ( iNode != nbNodes )
4415 continue; // not all nodes transformed
4417 if ( theTargetMesh ) {
4418 if ( SMDS_MeshElement* copy =
4419 targetMeshEditor.AddElement( nodes, elem->GetType(), elem->IsPoly() )) {
4420 myLastCreatedElems.Append( copy );
4421 srcElems.Append( elem );
4424 else if ( theCopy ) {
4425 if ( SMDS_MeshElement* copy = AddElement( nodes, elem->GetType(), elem->IsPoly() )) {
4426 myLastCreatedElems.Append( copy );
4427 srcElems.Append( elem );
4431 // reverse element as it was reversed by transformation
4433 aMesh->ChangeElementNodes( elem, &nodes[0], nbNodes );
4437 PGroupIDs newGroupIDs;
4439 if ( theMakeGroups && theCopy ||
4440 theMakeGroups && theTargetMesh )
4441 newGroupIDs = generateGroups( srcNodes, srcElems, groupPostfix, theTargetMesh );
4446 //=======================================================================
4448 * \brief Create groups of elements made during transformation
4449 * \param nodeGens - nodes making corresponding myLastCreatedNodes
4450 * \param elemGens - elements making corresponding myLastCreatedElems
4451 * \param postfix - to append to names of new groups
4453 //=======================================================================
4455 SMESH_MeshEditor::PGroupIDs
4456 SMESH_MeshEditor::generateGroups(const SMESH_SequenceOfElemPtr& nodeGens,
4457 const SMESH_SequenceOfElemPtr& elemGens,
4458 const std::string& postfix,
4459 SMESH_Mesh* targetMesh)
4461 PGroupIDs newGroupIDs( new list<int> );
4462 SMESH_Mesh* mesh = targetMesh ? targetMesh : GetMesh();
4464 // Sort existing groups by types and collect their names
4466 // to store an old group and a generated new one
4467 typedef pair< SMESHDS_GroupBase*, SMDS_MeshGroup* > TOldNewGroup;
4468 vector< list< TOldNewGroup > > groupsByType( SMDSAbs_NbElementTypes );
4470 set< string > groupNames;
4472 SMDS_MeshGroup* nullNewGroup = (SMDS_MeshGroup*) 0;
4473 SMESH_Mesh::GroupIteratorPtr groupIt = GetMesh()->GetGroups();
4474 while ( groupIt->more() ) {
4475 SMESH_Group * group = groupIt->next();
4476 if ( !group ) continue;
4477 SMESHDS_GroupBase* groupDS = group->GetGroupDS();
4478 if ( !groupDS || groupDS->IsEmpty() ) continue;
4479 groupNames.insert( group->GetName() );
4480 groupDS->SetStoreName( group->GetName() );
4481 groupsByType[ groupDS->GetType() ].push_back( make_pair( groupDS, nullNewGroup ));
4486 // loop on nodes and elements
4487 for ( int isNodes = 0; isNodes < 2; ++isNodes )
4489 const SMESH_SequenceOfElemPtr& gens = isNodes ? nodeGens : elemGens;
4490 const SMESH_SequenceOfElemPtr& elems = isNodes ? myLastCreatedNodes : myLastCreatedElems;
4491 if ( gens.Length() != elems.Length() )
4492 throw SALOME_Exception(LOCALIZED("invalid args"));
4494 // loop on created elements
4495 for (int iElem = 1; iElem <= elems.Length(); ++iElem )
4497 const SMDS_MeshElement* sourceElem = gens( iElem );
4498 if ( !sourceElem ) {
4499 MESSAGE("generateGroups(): NULL source element");
4502 list< TOldNewGroup > & groupsOldNew = groupsByType[ sourceElem->GetType() ];
4503 if ( groupsOldNew.empty() ) {
4504 while ( iElem < gens.Length() && gens( iElem+1 ) == sourceElem )
4505 ++iElem; // skip all elements made by sourceElem
4508 // collect all elements made by sourceElem
4509 list< const SMDS_MeshElement* > resultElems;
4510 if ( const SMDS_MeshElement* resElem = elems( iElem ))
4511 if ( resElem != sourceElem )
4512 resultElems.push_back( resElem );
4513 while ( iElem < gens.Length() && gens( iElem+1 ) == sourceElem )
4514 if ( const SMDS_MeshElement* resElem = elems( ++iElem ))
4515 if ( resElem != sourceElem )
4516 resultElems.push_back( resElem );
4517 // do not generate element groups from node ones
4518 if ( sourceElem->GetType() == SMDSAbs_Node &&
4519 elems( iElem )->GetType() != SMDSAbs_Node )
4522 // add resultElems to groups made by ones the sourceElem belongs to
4523 list< TOldNewGroup >::iterator gOldNew, gLast = groupsOldNew.end();
4524 for ( gOldNew = groupsOldNew.begin(); gOldNew != gLast; ++gOldNew )
4526 SMESHDS_GroupBase* oldGroup = gOldNew->first;
4527 if ( oldGroup->Contains( sourceElem )) // sourceElem in oldGroup
4529 SMDS_MeshGroup* & newGroup = gOldNew->second;
4530 if ( !newGroup )// create a new group
4533 string name = oldGroup->GetStoreName();
4534 if ( !targetMesh ) {
4538 while ( !groupNames.insert( name ).second ) // name exists
4544 TCollection_AsciiString nbStr(nb+1);
4545 name.resize( name.rfind('_')+1 );
4546 name += nbStr.ToCString();
4553 SMESH_Group* group = mesh->AddGroup( resultElems.back()->GetType(),
4555 SMESHDS_Group* groupDS = static_cast<SMESHDS_Group*>(group->GetGroupDS());
4556 newGroup = & groupDS->SMDSGroup();
4557 newGroupIDs->push_back( id );
4560 // fill in a new group
4561 list< const SMDS_MeshElement* >::iterator resLast = resultElems.end(), resElemIt;
4562 for ( resElemIt = resultElems.begin(); resElemIt != resLast; ++resElemIt )
4563 newGroup->Add( *resElemIt );
4566 } // loop on created elements
4567 }// loop on nodes and elements
4572 //=======================================================================
4573 //function : FindCoincidentNodes
4574 //purpose : Return list of group of nodes close to each other within theTolerance
4575 // Search among theNodes or in the whole mesh if theNodes is empty using
4576 // an Octree algorithm
4577 //=======================================================================
4579 void SMESH_MeshEditor::FindCoincidentNodes (set<const SMDS_MeshNode*> & theNodes,
4580 const double theTolerance,
4581 TListOfListOfNodes & theGroupsOfNodes)
4583 myLastCreatedElems.Clear();
4584 myLastCreatedNodes.Clear();
4586 set<const SMDS_MeshNode*> nodes;
4587 if ( theNodes.empty() )
4588 { // get all nodes in the mesh
4589 SMDS_NodeIteratorPtr nIt = GetMeshDS()->nodesIterator();
4590 while ( nIt->more() )
4591 nodes.insert( nodes.end(),nIt->next());
4595 SMESH_OctreeNode::FindCoincidentNodes ( nodes, &theGroupsOfNodes, theTolerance);
4599 //=======================================================================
4601 * \brief Implementation of search for the node closest to point
4603 //=======================================================================
4605 struct SMESH_NodeSearcherImpl: public SMESH_NodeSearcher
4608 * \brief Constructor
4610 SMESH_NodeSearcherImpl( const SMESHDS_Mesh* theMesh )
4612 set<const SMDS_MeshNode*> nodes;
4614 SMDS_NodeIteratorPtr nIt = theMesh->nodesIterator();
4615 while ( nIt->more() )
4616 nodes.insert( nodes.end(), nIt->next() );
4618 myOctreeNode = new SMESH_OctreeNode(nodes) ;
4621 * \brief Do it's job
4623 const SMDS_MeshNode* FindClosestTo( const gp_Pnt& thePnt )
4625 SMDS_MeshNode tgtNode( thePnt.X(), thePnt.Y(), thePnt.Z() );
4626 list<const SMDS_MeshNode*> nodes;
4627 const double precision = 1e-6;
4628 myOctreeNode->NodesAround( &tgtNode, &nodes, precision );
4630 double minSqDist = DBL_MAX;
4632 if ( nodes.empty() ) // get all nodes of OctreeNode's closest to thePnt
4634 // sort leafs by their distance from thePnt
4635 typedef map< double, SMESH_OctreeNode* > TDistTreeMap;
4636 TDistTreeMap treeMap;
4637 list< SMESH_OctreeNode* > treeList;
4638 list< SMESH_OctreeNode* >::iterator trIt;
4639 treeList.push_back( myOctreeNode );
4640 for ( trIt = treeList.begin(); trIt != treeList.end(); ++trIt)
4642 SMESH_OctreeNode* tree = *trIt;
4643 if ( !tree->isLeaf() ) { // put children to the queue
4644 SMESH_OctreeNodeIteratorPtr cIt = tree->GetChildrenIterator();
4645 while ( cIt->more() )
4646 treeList.push_back( cIt->next() );
4648 else if ( tree->NbNodes() ) { // put tree to treeMap
4649 tree->getBox( box );
4650 double sqDist = thePnt.SquareDistance( 0.5 * ( box.CornerMin() + box.CornerMax() ));
4651 pair<TDistTreeMap::iterator,bool> it_in = treeMap.insert( make_pair( sqDist, tree ));
4652 if ( !it_in.second ) // not unique distance to box center
4653 treeMap.insert( it_in.first, make_pair( sqDist - 1e-13*treeMap.size(), tree ));
4656 // find distance after which there is no sense to check tree's
4657 double sqLimit = DBL_MAX;
4658 TDistTreeMap::iterator sqDist_tree = treeMap.begin();
4659 if ( treeMap.size() > 5 ) {
4660 SMESH_OctreeNode* closestTree = sqDist_tree->second;
4661 closestTree->getBox( box );
4662 double limit = sqrt( sqDist_tree->first ) + sqrt ( box.SquareExtent() );
4663 sqLimit = limit * limit;
4665 // get all nodes from trees
4666 for ( ; sqDist_tree != treeMap.end(); ++sqDist_tree) {
4667 if ( sqDist_tree->first > sqLimit )
4669 SMESH_OctreeNode* tree = sqDist_tree->second;
4670 tree->NodesAround( tree->GetNodeIterator()->next(), &nodes );
4673 // find closest among nodes
4674 minSqDist = DBL_MAX;
4675 const SMDS_MeshNode* closestNode = 0;
4676 list<const SMDS_MeshNode*>::iterator nIt = nodes.begin();
4677 for ( ; nIt != nodes.end(); ++nIt ) {
4678 double sqDist = thePnt.SquareDistance( TNodeXYZ( *nIt ) );
4679 if ( minSqDist > sqDist ) {
4689 ~SMESH_NodeSearcherImpl() { delete myOctreeNode; }
4691 SMESH_OctreeNode* myOctreeNode;
4694 //=======================================================================
4696 * \brief Return SMESH_NodeSearcher
4698 //=======================================================================
4700 SMESH_NodeSearcher* SMESH_MeshEditor::GetNodeSearcher()
4702 return new SMESH_NodeSearcherImpl( GetMeshDS() );
4705 //=======================================================================
4706 //function : SimplifyFace
4708 //=======================================================================
4709 int SMESH_MeshEditor::SimplifyFace (const vector<const SMDS_MeshNode *> faceNodes,
4710 vector<const SMDS_MeshNode *>& poly_nodes,
4711 vector<int>& quantities) const
4713 int nbNodes = faceNodes.size();
4718 set<const SMDS_MeshNode*> nodeSet;
4720 // get simple seq of nodes
4721 //const SMDS_MeshNode* simpleNodes[ nbNodes ];
4722 vector<const SMDS_MeshNode*> simpleNodes( nbNodes );
4723 int iSimple = 0, nbUnique = 0;
4725 simpleNodes[iSimple++] = faceNodes[0];
4727 for (int iCur = 1; iCur < nbNodes; iCur++) {
4728 if (faceNodes[iCur] != simpleNodes[iSimple - 1]) {
4729 simpleNodes[iSimple++] = faceNodes[iCur];
4730 if (nodeSet.insert( faceNodes[iCur] ).second)
4734 int nbSimple = iSimple;
4735 if (simpleNodes[nbSimple - 1] == simpleNodes[0]) {
4745 bool foundLoop = (nbSimple > nbUnique);
4748 set<const SMDS_MeshNode*> loopSet;
4749 for (iSimple = 0; iSimple < nbSimple && !foundLoop; iSimple++) {
4750 const SMDS_MeshNode* n = simpleNodes[iSimple];
4751 if (!loopSet.insert( n ).second) {
4755 int iC = 0, curLast = iSimple;
4756 for (; iC < curLast; iC++) {
4757 if (simpleNodes[iC] == n) break;
4759 int loopLen = curLast - iC;
4761 // create sub-element
4763 quantities.push_back(loopLen);
4764 for (; iC < curLast; iC++) {
4765 poly_nodes.push_back(simpleNodes[iC]);
4768 // shift the rest nodes (place from the first loop position)
4769 for (iC = curLast + 1; iC < nbSimple; iC++) {
4770 simpleNodes[iC - loopLen] = simpleNodes[iC];
4772 nbSimple -= loopLen;
4775 } // for (iSimple = 0; iSimple < nbSimple; iSimple++)
4776 } // while (foundLoop)
4780 quantities.push_back(iSimple);
4781 for (int i = 0; i < iSimple; i++)
4782 poly_nodes.push_back(simpleNodes[i]);
4788 //=======================================================================
4789 //function : MergeNodes
4790 //purpose : In each group, the cdr of nodes are substituted by the first one
4792 //=======================================================================
4794 void SMESH_MeshEditor::MergeNodes (TListOfListOfNodes & theGroupsOfNodes)
4796 myLastCreatedElems.Clear();
4797 myLastCreatedNodes.Clear();
4799 SMESHDS_Mesh* aMesh = GetMeshDS();
4801 TNodeNodeMap nodeNodeMap; // node to replace - new node
4802 set<const SMDS_MeshElement*> elems; // all elements with changed nodes
4803 list< int > rmElemIds, rmNodeIds;
4805 // Fill nodeNodeMap and elems
4807 TListOfListOfNodes::iterator grIt = theGroupsOfNodes.begin();
4808 for ( ; grIt != theGroupsOfNodes.end(); grIt++ ) {
4809 list<const SMDS_MeshNode*>& nodes = *grIt;
4810 list<const SMDS_MeshNode*>::iterator nIt = nodes.begin();
4811 const SMDS_MeshNode* nToKeep = *nIt;
4812 for ( ++nIt; nIt != nodes.end(); nIt++ ) {
4813 const SMDS_MeshNode* nToRemove = *nIt;
4814 nodeNodeMap.insert( TNodeNodeMap::value_type( nToRemove, nToKeep ));
4815 if ( nToRemove != nToKeep ) {
4816 rmNodeIds.push_back( nToRemove->GetID() );
4817 AddToSameGroups( nToKeep, nToRemove, aMesh );
4820 SMDS_ElemIteratorPtr invElemIt = nToRemove->GetInverseElementIterator();
4821 while ( invElemIt->more() ) {
4822 const SMDS_MeshElement* elem = invElemIt->next();
4827 // Change element nodes or remove an element
4829 set<const SMDS_MeshElement*>::iterator eIt = elems.begin();
4830 for ( ; eIt != elems.end(); eIt++ ) {
4831 const SMDS_MeshElement* elem = *eIt;
4832 int nbNodes = elem->NbNodes();
4833 int aShapeId = FindShape( elem );
4835 set<const SMDS_MeshNode*> nodeSet;
4836 vector< const SMDS_MeshNode*> curNodes( nbNodes ), uniqueNodes( nbNodes );
4837 int iUnique = 0, iCur = 0, nbRepl = 0;
4838 vector<int> iRepl( nbNodes );
4840 // get new seq of nodes
4841 SMDS_ElemIteratorPtr itN = elem->nodesIterator();
4842 while ( itN->more() ) {
4843 const SMDS_MeshNode* n =
4844 static_cast<const SMDS_MeshNode*>( itN->next() );
4846 TNodeNodeMap::iterator nnIt = nodeNodeMap.find( n );
4847 if ( nnIt != nodeNodeMap.end() ) { // n sticks
4849 // BUG 0020185: begin
4851 bool stopRecur = false;
4852 set<const SMDS_MeshNode*> nodesRecur;
4853 nodesRecur.insert(n);
4854 while (!stopRecur) {
4855 TNodeNodeMap::iterator nnIt_i = nodeNodeMap.find( n );
4856 if ( nnIt_i != nodeNodeMap.end() ) { // n sticks
4857 n = (*nnIt_i).second;
4858 if (!nodesRecur.insert(n).second) {
4859 // error: recursive dependancy
4868 iRepl[ nbRepl++ ] = iCur;
4870 curNodes[ iCur ] = n;
4871 bool isUnique = nodeSet.insert( n ).second;
4873 uniqueNodes[ iUnique++ ] = n;
4877 // Analyse element topology after replacement
4880 int nbUniqueNodes = nodeSet.size();
4881 if ( nbNodes != nbUniqueNodes ) { // some nodes stick
4882 // Polygons and Polyhedral volumes
4883 if (elem->IsPoly()) {
4885 if (elem->GetType() == SMDSAbs_Face) {
4887 vector<const SMDS_MeshNode *> face_nodes (nbNodes);
4889 for (; inode < nbNodes; inode++) {
4890 face_nodes[inode] = curNodes[inode];
4893 vector<const SMDS_MeshNode *> polygons_nodes;
4894 vector<int> quantities;
4895 int nbNew = SimplifyFace(face_nodes, polygons_nodes, quantities);
4899 for (int iface = 0; iface < nbNew - 1; iface++) {
4900 int nbNodes = quantities[iface];
4901 vector<const SMDS_MeshNode *> poly_nodes (nbNodes);
4902 for (int ii = 0; ii < nbNodes; ii++, inode++) {
4903 poly_nodes[ii] = polygons_nodes[inode];
4905 SMDS_MeshElement* newElem = aMesh->AddPolygonalFace(poly_nodes);
4906 myLastCreatedElems.Append(newElem);
4908 aMesh->SetMeshElementOnShape(newElem, aShapeId);
4910 aMesh->ChangeElementNodes(elem, &polygons_nodes[inode], quantities[nbNew - 1]);
4913 rmElemIds.push_back(elem->GetID());
4917 else if (elem->GetType() == SMDSAbs_Volume) {
4918 // Polyhedral volume
4919 if (nbUniqueNodes < 4) {
4920 rmElemIds.push_back(elem->GetID());
4923 // each face has to be analized in order to check volume validity
4924 const SMDS_PolyhedralVolumeOfNodes* aPolyedre =
4925 static_cast<const SMDS_PolyhedralVolumeOfNodes*>( elem );
4927 int nbFaces = aPolyedre->NbFaces();
4929 vector<const SMDS_MeshNode *> poly_nodes;
4930 vector<int> quantities;
4932 for (int iface = 1; iface <= nbFaces; iface++) {
4933 int nbFaceNodes = aPolyedre->NbFaceNodes(iface);
4934 vector<const SMDS_MeshNode *> faceNodes (nbFaceNodes);
4936 for (int inode = 1; inode <= nbFaceNodes; inode++) {
4937 const SMDS_MeshNode * faceNode = aPolyedre->GetFaceNode(iface, inode);
4938 TNodeNodeMap::iterator nnIt = nodeNodeMap.find(faceNode);
4939 if (nnIt != nodeNodeMap.end()) { // faceNode sticks
4940 faceNode = (*nnIt).second;
4942 faceNodes[inode - 1] = faceNode;
4945 SimplifyFace(faceNodes, poly_nodes, quantities);
4948 if (quantities.size() > 3) {
4949 // to be done: remove coincident faces
4952 if (quantities.size() > 3)
4953 aMesh->ChangePolyhedronNodes(elem, poly_nodes, quantities);
4955 rmElemIds.push_back(elem->GetID());
4959 rmElemIds.push_back(elem->GetID());
4970 switch ( nbNodes ) {
4971 case 2: ///////////////////////////////////// EDGE
4972 isOk = false; break;
4973 case 3: ///////////////////////////////////// TRIANGLE
4974 isOk = false; break;
4976 if ( elem->GetType() == SMDSAbs_Volume ) // TETRAHEDRON
4978 else { //////////////////////////////////// QUADRANGLE
4979 if ( nbUniqueNodes < 3 )
4981 else if ( nbRepl == 2 && iRepl[ 1 ] - iRepl[ 0 ] == 2 )
4982 isOk = false; // opposite nodes stick
4985 case 6: ///////////////////////////////////// PENTAHEDRON
4986 if ( nbUniqueNodes == 4 ) {
4987 // ---------------------------------> tetrahedron
4989 iRepl[ 0 ] > 2 && iRepl[ 1 ] > 2 && iRepl[ 2 ] > 2 ) {
4990 // all top nodes stick: reverse a bottom
4991 uniqueNodes[ 0 ] = curNodes [ 1 ];
4992 uniqueNodes[ 1 ] = curNodes [ 0 ];
4994 else if (nbRepl == 3 &&
4995 iRepl[ 0 ] < 3 && iRepl[ 1 ] < 3 && iRepl[ 2 ] < 3 ) {
4996 // all bottom nodes stick: set a top before
4997 uniqueNodes[ 3 ] = uniqueNodes [ 0 ];
4998 uniqueNodes[ 0 ] = curNodes [ 3 ];
4999 uniqueNodes[ 1 ] = curNodes [ 4 ];
5000 uniqueNodes[ 2 ] = curNodes [ 5 ];
5002 else if (nbRepl == 4 &&
5003 iRepl[ 2 ] - iRepl [ 0 ] == 3 && iRepl[ 3 ] - iRepl [ 1 ] == 3 ) {
5004 // a lateral face turns into a line: reverse a bottom
5005 uniqueNodes[ 0 ] = curNodes [ 1 ];
5006 uniqueNodes[ 1 ] = curNodes [ 0 ];
5011 else if ( nbUniqueNodes == 5 ) {
5012 // PENTAHEDRON --------------------> 2 tetrahedrons
5013 if ( nbRepl == 2 && iRepl[ 1 ] - iRepl [ 0 ] == 3 ) {
5014 // a bottom node sticks with a linked top one
5016 SMDS_MeshElement* newElem =
5017 aMesh->AddVolume(curNodes[ 3 ],
5020 curNodes[ iRepl[ 0 ] == 2 ? 1 : 2 ]);
5021 myLastCreatedElems.Append(newElem);
5023 aMesh->SetMeshElementOnShape( newElem, aShapeId );
5024 // 2. : reverse a bottom
5025 uniqueNodes[ 0 ] = curNodes [ 1 ];
5026 uniqueNodes[ 1 ] = curNodes [ 0 ];
5036 if(elem->IsQuadratic()) { // Quadratic quadrangle
5049 if( iRepl[0]==0 && iRepl[1]==1 && iRepl[2]==4 ) {
5050 uniqueNodes[0] = curNodes[0];
5051 uniqueNodes[1] = curNodes[2];
5052 uniqueNodes[2] = curNodes[3];
5053 uniqueNodes[3] = curNodes[5];
5054 uniqueNodes[4] = curNodes[6];
5055 uniqueNodes[5] = curNodes[7];
5058 if( iRepl[0]==0 && iRepl[1]==3 && iRepl[2]==7 ) {
5059 uniqueNodes[0] = curNodes[0];
5060 uniqueNodes[1] = curNodes[1];
5061 uniqueNodes[2] = curNodes[2];
5062 uniqueNodes[3] = curNodes[4];
5063 uniqueNodes[4] = curNodes[5];
5064 uniqueNodes[5] = curNodes[6];
5067 if( iRepl[0]==0 && iRepl[1]==4 && iRepl[2]==7 ) {
5068 uniqueNodes[0] = curNodes[1];
5069 uniqueNodes[1] = curNodes[2];
5070 uniqueNodes[2] = curNodes[3];
5071 uniqueNodes[3] = curNodes[5];
5072 uniqueNodes[4] = curNodes[6];
5073 uniqueNodes[5] = curNodes[0];
5076 if( iRepl[0]==1 && iRepl[1]==2 && iRepl[2]==5 ) {
5077 uniqueNodes[0] = curNodes[0];
5078 uniqueNodes[1] = curNodes[1];
5079 uniqueNodes[2] = curNodes[3];
5080 uniqueNodes[3] = curNodes[4];
5081 uniqueNodes[4] = curNodes[6];
5082 uniqueNodes[5] = curNodes[7];
5085 if( iRepl[0]==1 && iRepl[1]==4 && iRepl[2]==5 ) {
5086 uniqueNodes[0] = curNodes[0];
5087 uniqueNodes[1] = curNodes[2];
5088 uniqueNodes[2] = curNodes[3];
5089 uniqueNodes[3] = curNodes[1];
5090 uniqueNodes[4] = curNodes[6];
5091 uniqueNodes[5] = curNodes[7];
5094 if( iRepl[0]==2 && iRepl[1]==3 && iRepl[2]==6 ) {
5095 uniqueNodes[0] = curNodes[0];
5096 uniqueNodes[1] = curNodes[1];
5097 uniqueNodes[2] = curNodes[2];
5098 uniqueNodes[3] = curNodes[4];
5099 uniqueNodes[4] = curNodes[5];
5100 uniqueNodes[5] = curNodes[7];
5103 if( iRepl[0]==2 && iRepl[1]==5 && iRepl[2]==6 ) {
5104 uniqueNodes[0] = curNodes[0];
5105 uniqueNodes[1] = curNodes[1];
5106 uniqueNodes[2] = curNodes[3];
5107 uniqueNodes[3] = curNodes[4];
5108 uniqueNodes[4] = curNodes[2];
5109 uniqueNodes[5] = curNodes[7];
5112 if( iRepl[0]==3 && iRepl[1]==6 && iRepl[2]==7 ) {
5113 uniqueNodes[0] = curNodes[0];
5114 uniqueNodes[1] = curNodes[1];
5115 uniqueNodes[2] = curNodes[2];
5116 uniqueNodes[3] = curNodes[4];
5117 uniqueNodes[4] = curNodes[5];
5118 uniqueNodes[5] = curNodes[3];
5124 //////////////////////////////////// HEXAHEDRON
5126 SMDS_VolumeTool hexa (elem);
5127 hexa.SetExternalNormal();
5128 if ( nbUniqueNodes == 4 && nbRepl == 6 ) {
5129 //////////////////////// ---> tetrahedron
5130 for ( int iFace = 0; iFace < 6; iFace++ ) {
5131 const int *ind = hexa.GetFaceNodesIndices( iFace ); // indices of face nodes
5132 if (curNodes[ind[ 0 ]] == curNodes[ind[ 1 ]] &&
5133 curNodes[ind[ 0 ]] == curNodes[ind[ 2 ]] &&
5134 curNodes[ind[ 0 ]] == curNodes[ind[ 3 ]] ) {
5135 // one face turns into a point ...
5136 int iOppFace = hexa.GetOppFaceIndex( iFace );
5137 ind = hexa.GetFaceNodesIndices( iOppFace );
5139 iUnique = 2; // reverse a tetrahedron bottom
5140 for ( iCur = 0; iCur < 4 && nbStick < 2; iCur++ ) {
5141 if ( curNodes[ind[ iCur ]] == curNodes[ind[ iCur + 1 ]] )
5143 else if ( iUnique >= 0 )
5144 uniqueNodes[ iUnique-- ] = curNodes[ind[ iCur ]];
5146 if ( nbStick == 1 ) {
5147 // ... and the opposite one - into a triangle.
5149 ind = hexa.GetFaceNodesIndices( iFace );
5150 uniqueNodes[ 3 ] = curNodes[ind[ 0 ]];
5157 else if (nbUniqueNodes == 5 && nbRepl == 4 ) {
5158 //////////////////// HEXAHEDRON ---> 2 tetrahedrons
5159 for ( int iFace = 0; iFace < 6; iFace++ ) {
5160 const int *ind = hexa.GetFaceNodesIndices( iFace ); // indices of face nodes
5161 if (curNodes[ind[ 0 ]] == curNodes[ind[ 1 ]] &&
5162 curNodes[ind[ 0 ]] == curNodes[ind[ 2 ]] &&
5163 curNodes[ind[ 0 ]] == curNodes[ind[ 3 ]] ) {
5164 // one face turns into a point ...
5165 int iOppFace = hexa.GetOppFaceIndex( iFace );
5166 ind = hexa.GetFaceNodesIndices( iOppFace );
5168 iUnique = 2; // reverse a tetrahedron 1 bottom
5169 for ( iCur = 0; iCur < 4 && nbStick == 0; iCur++ ) {
5170 if ( curNodes[ind[ iCur ]] == curNodes[ind[ iCur + 1 ]] )
5172 else if ( iUnique >= 0 )
5173 uniqueNodes[ iUnique-- ] = curNodes[ind[ iCur ]];
5175 if ( nbStick == 0 ) {
5176 // ... and the opposite one is a quadrangle
5178 const int* indTop = hexa.GetFaceNodesIndices( iFace );
5179 uniqueNodes[ 3 ] = curNodes[indTop[ 0 ]];
5182 SMDS_MeshElement* newElem =
5183 aMesh->AddVolume(curNodes[ind[ 0 ]],
5186 curNodes[indTop[ 0 ]]);
5187 myLastCreatedElems.Append(newElem);
5189 aMesh->SetMeshElementOnShape( newElem, aShapeId );
5196 else if ( nbUniqueNodes == 6 && nbRepl == 4 ) {
5197 ////////////////// HEXAHEDRON ---> 2 tetrahedrons or 1 prism
5198 // find indices of quad and tri faces
5199 int iQuadFace[ 6 ], iTriFace[ 6 ], nbQuad = 0, nbTri = 0, iFace;
5200 for ( iFace = 0; iFace < 6; iFace++ ) {
5201 const int *ind = hexa.GetFaceNodesIndices( iFace ); // indices of face nodes
5203 for ( iCur = 0; iCur < 4; iCur++ )
5204 nodeSet.insert( curNodes[ind[ iCur ]] );
5205 nbUniqueNodes = nodeSet.size();
5206 if ( nbUniqueNodes == 3 )
5207 iTriFace[ nbTri++ ] = iFace;
5208 else if ( nbUniqueNodes == 4 )
5209 iQuadFace[ nbQuad++ ] = iFace;
5211 if (nbQuad == 2 && nbTri == 4 &&
5212 hexa.GetOppFaceIndex( iQuadFace[ 0 ] ) == iQuadFace[ 1 ]) {
5213 // 2 opposite quadrangles stuck with a diagonal;
5214 // sample groups of merged indices: (0-4)(2-6)
5215 // --------------------------------------------> 2 tetrahedrons
5216 const int *ind1 = hexa.GetFaceNodesIndices( iQuadFace[ 0 ]); // indices of quad1 nodes
5217 const int *ind2 = hexa.GetFaceNodesIndices( iQuadFace[ 1 ]);
5218 int i0, i1d, i2, i3d, i0t, i2t; // d-daigonal, t-top
5219 if (curNodes[ind1[ 0 ]] == curNodes[ind2[ 0 ]] &&
5220 curNodes[ind1[ 2 ]] == curNodes[ind2[ 2 ]]) {
5221 // stuck with 0-2 diagonal
5229 else if (curNodes[ind1[ 1 ]] == curNodes[ind2[ 3 ]] &&
5230 curNodes[ind1[ 3 ]] == curNodes[ind2[ 1 ]]) {
5231 // stuck with 1-3 diagonal
5243 uniqueNodes[ 0 ] = curNodes [ i0 ];
5244 uniqueNodes[ 1 ] = curNodes [ i1d ];
5245 uniqueNodes[ 2 ] = curNodes [ i3d ];
5246 uniqueNodes[ 3 ] = curNodes [ i0t ];
5249 SMDS_MeshElement* newElem = aMesh->AddVolume(curNodes[ i1d ],
5253 myLastCreatedElems.Append(newElem);
5255 aMesh->SetMeshElementOnShape( newElem, aShapeId );
5258 else if (( nbTri == 2 && nbQuad == 3 ) || // merged (0-4)(1-5)
5259 ( nbTri == 4 && nbQuad == 2 )) { // merged (7-4)(1-5)
5260 // --------------------------------------------> prism
5261 // find 2 opposite triangles
5263 for ( iFace = 0; iFace + 1 < nbTri; iFace++ ) {
5264 if ( hexa.GetOppFaceIndex( iTriFace[ iFace ] ) == iTriFace[ iFace + 1 ]) {
5265 // find indices of kept and replaced nodes
5266 // and fill unique nodes of 2 opposite triangles
5267 const int *ind1 = hexa.GetFaceNodesIndices( iTriFace[ iFace ]);
5268 const int *ind2 = hexa.GetFaceNodesIndices( iTriFace[ iFace + 1 ]);
5269 const SMDS_MeshNode** hexanodes = hexa.GetNodes();
5270 // fill unique nodes
5273 for ( iCur = 0; iCur < 4 && isOk; iCur++ ) {
5274 const SMDS_MeshNode* n = curNodes[ind1[ iCur ]];
5275 const SMDS_MeshNode* nInit = hexanodes[ind1[ iCur ]];
5277 // iCur of a linked node of the opposite face (make normals co-directed):
5278 int iCurOpp = ( iCur == 1 || iCur == 3 ) ? 4 - iCur : iCur;
5279 // check that correspondent corners of triangles are linked
5280 if ( !hexa.IsLinked( ind1[ iCur ], ind2[ iCurOpp ] ))
5283 uniqueNodes[ iUnique ] = n;
5284 uniqueNodes[ iUnique + 3 ] = curNodes[ind2[ iCurOpp ]];
5293 } // if ( nbUniqueNodes == 6 && nbRepl == 4 )
5299 } // switch ( nbNodes )
5301 } // if ( nbNodes != nbUniqueNodes ) // some nodes stick
5304 if (elem->IsPoly() && elem->GetType() == SMDSAbs_Volume) {
5305 // Change nodes of polyedre
5306 const SMDS_PolyhedralVolumeOfNodes* aPolyedre =
5307 static_cast<const SMDS_PolyhedralVolumeOfNodes*>( elem );
5309 int nbFaces = aPolyedre->NbFaces();
5311 vector<const SMDS_MeshNode *> poly_nodes;
5312 vector<int> quantities (nbFaces);
5314 for (int iface = 1; iface <= nbFaces; iface++) {
5315 int inode, nbFaceNodes = aPolyedre->NbFaceNodes(iface);
5316 quantities[iface - 1] = nbFaceNodes;
5318 for (inode = 1; inode <= nbFaceNodes; inode++) {
5319 const SMDS_MeshNode* curNode = aPolyedre->GetFaceNode(iface, inode);
5321 TNodeNodeMap::iterator nnIt = nodeNodeMap.find( curNode );
5322 if (nnIt != nodeNodeMap.end()) { // curNode sticks
5323 curNode = (*nnIt).second;
5325 poly_nodes.push_back(curNode);
5328 aMesh->ChangePolyhedronNodes( elem, poly_nodes, quantities );
5332 // Change regular element or polygon
5333 aMesh->ChangeElementNodes( elem, & uniqueNodes[0], nbUniqueNodes );
5337 // Remove invalid regular element or invalid polygon
5338 rmElemIds.push_back( elem->GetID() );
5341 } // loop on elements
5343 // Remove equal nodes and bad elements
5345 Remove( rmNodeIds, true );
5346 Remove( rmElemIds, false );
5351 // ========================================================
5352 // class : SortableElement
5353 // purpose : allow sorting elements basing on their nodes
5354 // ========================================================
5355 class SortableElement : public set <const SMDS_MeshElement*>
5359 SortableElement( const SMDS_MeshElement* theElem )
5362 SMDS_ElemIteratorPtr nodeIt = theElem->nodesIterator();
5363 while ( nodeIt->more() )
5364 this->insert( nodeIt->next() );
5367 const SMDS_MeshElement* Get() const
5370 void Set(const SMDS_MeshElement* e) const
5375 mutable const SMDS_MeshElement* myElem;
5378 //=======================================================================
5379 //function : FindEqualElements
5380 //purpose : Return list of group of elements built on the same nodes.
5381 // Search among theElements or in the whole mesh if theElements is empty
5382 //=======================================================================
5383 void SMESH_MeshEditor::FindEqualElements(set<const SMDS_MeshElement*> & theElements,
5384 TListOfListOfElementsID & theGroupsOfElementsID)
5386 myLastCreatedElems.Clear();
5387 myLastCreatedNodes.Clear();
5389 typedef set<const SMDS_MeshElement*> TElemsSet;
5390 typedef map< SortableElement, int > TMapOfNodeSet;
5391 typedef list<int> TGroupOfElems;
5394 if ( theElements.empty() )
5395 { // get all elements in the mesh
5396 SMDS_ElemIteratorPtr eIt = GetMeshDS()->elementsIterator();
5397 while ( eIt->more() )
5398 elems.insert( elems.end(), eIt->next());
5401 elems = theElements;
5403 vector< TGroupOfElems > arrayOfGroups;
5404 TGroupOfElems groupOfElems;
5405 TMapOfNodeSet mapOfNodeSet;
5407 TElemsSet::iterator elemIt = elems.begin();
5408 for ( int i = 0, j=0; elemIt != elems.end(); ++elemIt, ++j ) {
5409 const SMDS_MeshElement* curElem = *elemIt;
5410 SortableElement SE(curElem);
5413 pair< TMapOfNodeSet::iterator, bool> pp = mapOfNodeSet.insert(make_pair(SE, i));
5414 if( !(pp.second) ) {
5415 TMapOfNodeSet::iterator& itSE = pp.first;
5416 ind = (*itSE).second;
5417 arrayOfGroups[ind].push_back(curElem->GetID());
5420 groupOfElems.clear();
5421 groupOfElems.push_back(curElem->GetID());
5422 arrayOfGroups.push_back(groupOfElems);
5427 vector< TGroupOfElems >::iterator groupIt = arrayOfGroups.begin();
5428 for ( ; groupIt != arrayOfGroups.end(); ++groupIt ) {
5429 groupOfElems = *groupIt;
5430 if ( groupOfElems.size() > 1 ) {
5431 groupOfElems.sort();
5432 theGroupsOfElementsID.push_back(groupOfElems);
5437 //=======================================================================
5438 //function : MergeElements
5439 //purpose : In each given group, substitute all elements by the first one.
5440 //=======================================================================
5442 void SMESH_MeshEditor::MergeElements(TListOfListOfElementsID & theGroupsOfElementsID)
5444 myLastCreatedElems.Clear();
5445 myLastCreatedNodes.Clear();
5447 typedef list<int> TListOfIDs;
5448 TListOfIDs rmElemIds; // IDs of elems to remove
5450 SMESHDS_Mesh* aMesh = GetMeshDS();
5452 TListOfListOfElementsID::iterator groupsIt = theGroupsOfElementsID.begin();
5453 while ( groupsIt != theGroupsOfElementsID.end() ) {
5454 TListOfIDs& aGroupOfElemID = *groupsIt;
5455 aGroupOfElemID.sort();
5456 int elemIDToKeep = aGroupOfElemID.front();
5457 const SMDS_MeshElement* elemToKeep = aMesh->FindElement(elemIDToKeep);
5458 aGroupOfElemID.pop_front();
5459 TListOfIDs::iterator idIt = aGroupOfElemID.begin();
5460 while ( idIt != aGroupOfElemID.end() ) {
5461 int elemIDToRemove = *idIt;
5462 const SMDS_MeshElement* elemToRemove = aMesh->FindElement(elemIDToRemove);
5463 // add the kept element in groups of removed one (PAL15188)
5464 AddToSameGroups( elemToKeep, elemToRemove, aMesh );
5465 rmElemIds.push_back( elemIDToRemove );
5471 Remove( rmElemIds, false );
5474 //=======================================================================
5475 //function : MergeEqualElements
5476 //purpose : Remove all but one of elements built on the same nodes.
5477 //=======================================================================
5479 void SMESH_MeshEditor::MergeEqualElements()
5481 set<const SMDS_MeshElement*> aMeshElements; /* empty input -
5482 to merge equal elements in the whole mesh */
5483 TListOfListOfElementsID aGroupsOfElementsID;
5484 FindEqualElements(aMeshElements, aGroupsOfElementsID);
5485 MergeElements(aGroupsOfElementsID);
5488 //=======================================================================
5489 //function : FindFaceInSet
5490 //purpose : Return a face having linked nodes n1 and n2 and which is
5491 // - not in avoidSet,
5492 // - in elemSet provided that !elemSet.empty()
5493 //=======================================================================
5495 const SMDS_MeshElement*
5496 SMESH_MeshEditor::FindFaceInSet(const SMDS_MeshNode* n1,
5497 const SMDS_MeshNode* n2,
5498 const TIDSortedElemSet& elemSet,
5499 const TIDSortedElemSet& avoidSet)
5502 SMDS_ElemIteratorPtr invElemIt = n1->GetInverseElementIterator(SMDSAbs_Face);
5503 while ( invElemIt->more() ) { // loop on inverse elements of n1
5504 const SMDS_MeshElement* elem = invElemIt->next();
5505 if (avoidSet.find( elem ) != avoidSet.end() )
5507 if ( !elemSet.empty() && elemSet.find( elem ) == elemSet.end())
5509 // get face nodes and find index of n1
5510 int i1, nbN = elem->NbNodes(), iNode = 0;
5511 //const SMDS_MeshNode* faceNodes[ nbN ], *n;
5512 vector<const SMDS_MeshNode*> faceNodes( nbN );
5513 const SMDS_MeshNode* n;
5514 SMDS_ElemIteratorPtr nIt = elem->nodesIterator();
5515 while ( nIt->more() ) {
5516 faceNodes[ iNode ] = static_cast<const SMDS_MeshNode*>( nIt->next() );
5517 if ( faceNodes[ iNode++ ] == n1 )
5520 // find a n2 linked to n1
5521 if(!elem->IsQuadratic()) {
5522 for ( iNode = 0; iNode < 2; iNode++ ) {
5523 if ( iNode ) // node before n1
5524 n = faceNodes[ i1 == 0 ? nbN - 1 : i1 - 1 ];
5525 else // node after n1
5526 n = faceNodes[ i1 + 1 == nbN ? 0 : i1 + 1 ];
5531 else { // analysis for quadratic elements
5532 bool IsFind = false;
5533 // check using only corner nodes
5534 for ( iNode = 0; iNode < 2; iNode++ ) {
5535 if ( iNode ) // node before n1
5536 n = faceNodes[ i1 == 0 ? nbN/2 - 1 : i1 - 1 ];
5537 else // node after n1
5538 n = faceNodes[ i1 + 1 == nbN/2 ? 0 : i1 + 1 ];
5546 // check using all nodes
5547 const SMDS_QuadraticFaceOfNodes* F =
5548 static_cast<const SMDS_QuadraticFaceOfNodes*>(elem);
5549 // use special nodes iterator
5551 SMDS_NodeIteratorPtr anIter = F->interlacedNodesIterator();
5552 while ( anIter->more() ) {
5553 faceNodes[iNode] = static_cast<const SMDS_MeshNode*>(anIter->next());
5554 if ( faceNodes[ iNode++ ] == n1 )
5557 for ( iNode = 0; iNode < 2; iNode++ ) {
5558 if ( iNode ) // node before n1
5559 n = faceNodes[ i1 == 0 ? nbN - 1 : i1 - 1 ];
5560 else // node after n1
5561 n = faceNodes[ i1 + 1 == nbN ? 0 : i1 + 1 ];
5567 } // end analysis for quadratic elements
5572 //=======================================================================
5573 //function : findAdjacentFace
5575 //=======================================================================
5577 static const SMDS_MeshElement* findAdjacentFace(const SMDS_MeshNode* n1,
5578 const SMDS_MeshNode* n2,
5579 const SMDS_MeshElement* elem)
5581 TIDSortedElemSet elemSet, avoidSet;
5583 avoidSet.insert ( elem );
5584 return SMESH_MeshEditor::FindFaceInSet( n1, n2, elemSet, avoidSet );
5587 //=======================================================================
5588 //function : FindFreeBorder
5590 //=======================================================================
5592 #define ControlFreeBorder SMESH::Controls::FreeEdges::IsFreeEdge
5594 bool SMESH_MeshEditor::FindFreeBorder (const SMDS_MeshNode* theFirstNode,
5595 const SMDS_MeshNode* theSecondNode,
5596 const SMDS_MeshNode* theLastNode,
5597 list< const SMDS_MeshNode* > & theNodes,
5598 list< const SMDS_MeshElement* >& theFaces)
5600 if ( !theFirstNode || !theSecondNode )
5602 // find border face between theFirstNode and theSecondNode
5603 const SMDS_MeshElement* curElem = findAdjacentFace( theFirstNode, theSecondNode, 0 );
5607 theFaces.push_back( curElem );
5608 theNodes.push_back( theFirstNode );
5609 theNodes.push_back( theSecondNode );
5611 //vector<const SMDS_MeshNode*> nodes;
5612 const SMDS_MeshNode *nIgnore = theFirstNode, *nStart = theSecondNode;
5613 set < const SMDS_MeshElement* > foundElems;
5614 bool needTheLast = ( theLastNode != 0 );
5616 while ( nStart != theLastNode ) {
5617 if ( nStart == theFirstNode )
5618 return !needTheLast;
5620 // find all free border faces sharing form nStart
5622 list< const SMDS_MeshElement* > curElemList;
5623 list< const SMDS_MeshNode* > nStartList;
5624 SMDS_ElemIteratorPtr invElemIt = nStart->GetInverseElementIterator(SMDSAbs_Face);
5625 while ( invElemIt->more() ) {
5626 const SMDS_MeshElement* e = invElemIt->next();
5627 if ( e == curElem || foundElems.insert( e ).second ) {
5629 int iNode = 0, nbNodes = e->NbNodes();
5630 //const SMDS_MeshNode* nodes[nbNodes+1];
5631 vector<const SMDS_MeshNode*> nodes(nbNodes+1);
5633 if(e->IsQuadratic()) {
5634 const SMDS_QuadraticFaceOfNodes* F =
5635 static_cast<const SMDS_QuadraticFaceOfNodes*>(e);
5636 // use special nodes iterator
5637 SMDS_NodeIteratorPtr anIter = F->interlacedNodesIterator();
5638 while( anIter->more() ) {
5639 nodes[ iNode++ ] = anIter->next();
5643 SMDS_ElemIteratorPtr nIt = e->nodesIterator();
5644 while ( nIt->more() )
5645 nodes[ iNode++ ] = static_cast<const SMDS_MeshNode*>( nIt->next() );
5647 nodes[ iNode ] = nodes[ 0 ];
5649 for ( iNode = 0; iNode < nbNodes; iNode++ )
5650 if (((nodes[ iNode ] == nStart && nodes[ iNode + 1] != nIgnore ) ||
5651 (nodes[ iNode + 1] == nStart && nodes[ iNode ] != nIgnore )) &&
5652 ControlFreeBorder( &nodes[ iNode ], e->GetID() ))
5654 nStartList.push_back( nodes[ iNode + ( nodes[ iNode ] == nStart ? 1 : 0 )]);
5655 curElemList.push_back( e );
5659 // analyse the found
5661 int nbNewBorders = curElemList.size();
5662 if ( nbNewBorders == 0 ) {
5663 // no free border furthermore
5664 return !needTheLast;
5666 else if ( nbNewBorders == 1 ) {
5667 // one more element found
5669 nStart = nStartList.front();
5670 curElem = curElemList.front();
5671 theFaces.push_back( curElem );
5672 theNodes.push_back( nStart );
5675 // several continuations found
5676 list< const SMDS_MeshElement* >::iterator curElemIt;
5677 list< const SMDS_MeshNode* >::iterator nStartIt;
5678 // check if one of them reached the last node
5679 if ( needTheLast ) {
5680 for (curElemIt = curElemList.begin(), nStartIt = nStartList.begin();
5681 curElemIt!= curElemList.end();
5682 curElemIt++, nStartIt++ )
5683 if ( *nStartIt == theLastNode ) {
5684 theFaces.push_back( *curElemIt );
5685 theNodes.push_back( *nStartIt );
5689 // find the best free border by the continuations
5690 list<const SMDS_MeshNode*> contNodes[ 2 ], *cNL;
5691 list<const SMDS_MeshElement*> contFaces[ 2 ], *cFL;
5692 for (curElemIt = curElemList.begin(), nStartIt = nStartList.begin();
5693 curElemIt!= curElemList.end();
5694 curElemIt++, nStartIt++ )
5696 cNL = & contNodes[ contNodes[0].empty() ? 0 : 1 ];
5697 cFL = & contFaces[ contFaces[0].empty() ? 0 : 1 ];
5698 // find one more free border
5699 if ( ! FindFreeBorder( nStart, *nStartIt, theLastNode, *cNL, *cFL )) {
5703 else if ( !contNodes[0].empty() && !contNodes[1].empty() ) {
5704 // choice: clear a worse one
5705 int iLongest = ( contNodes[0].size() < contNodes[1].size() ? 1 : 0 );
5706 int iWorse = ( needTheLast ? 1 - iLongest : iLongest );
5707 contNodes[ iWorse ].clear();
5708 contFaces[ iWorse ].clear();
5711 if ( contNodes[0].empty() && contNodes[1].empty() )
5714 // append the best free border
5715 cNL = & contNodes[ contNodes[0].empty() ? 1 : 0 ];
5716 cFL = & contFaces[ contFaces[0].empty() ? 1 : 0 ];
5717 theNodes.pop_back(); // remove nIgnore
5718 theNodes.pop_back(); // remove nStart
5719 theFaces.pop_back(); // remove curElem
5720 list< const SMDS_MeshNode* >::iterator nIt = cNL->begin();
5721 list< const SMDS_MeshElement* >::iterator fIt = cFL->begin();
5722 for ( ; nIt != cNL->end(); nIt++ ) theNodes.push_back( *nIt );
5723 for ( ; fIt != cFL->end(); fIt++ ) theFaces.push_back( *fIt );
5726 } // several continuations found
5727 } // while ( nStart != theLastNode )
5732 //=======================================================================
5733 //function : CheckFreeBorderNodes
5734 //purpose : Return true if the tree nodes are on a free border
5735 //=======================================================================
5737 bool SMESH_MeshEditor::CheckFreeBorderNodes(const SMDS_MeshNode* theNode1,
5738 const SMDS_MeshNode* theNode2,
5739 const SMDS_MeshNode* theNode3)
5741 list< const SMDS_MeshNode* > nodes;
5742 list< const SMDS_MeshElement* > faces;
5743 return FindFreeBorder( theNode1, theNode2, theNode3, nodes, faces);
5746 //=======================================================================
5747 //function : SewFreeBorder
5749 //=======================================================================
5751 SMESH_MeshEditor::Sew_Error
5752 SMESH_MeshEditor::SewFreeBorder (const SMDS_MeshNode* theBordFirstNode,
5753 const SMDS_MeshNode* theBordSecondNode,
5754 const SMDS_MeshNode* theBordLastNode,
5755 const SMDS_MeshNode* theSideFirstNode,
5756 const SMDS_MeshNode* theSideSecondNode,
5757 const SMDS_MeshNode* theSideThirdNode,
5758 const bool theSideIsFreeBorder,
5759 const bool toCreatePolygons,
5760 const bool toCreatePolyedrs)
5762 myLastCreatedElems.Clear();
5763 myLastCreatedNodes.Clear();
5765 MESSAGE("::SewFreeBorder()");
5766 Sew_Error aResult = SEW_OK;
5768 // ====================================
5769 // find side nodes and elements
5770 // ====================================
5772 list< const SMDS_MeshNode* > nSide[ 2 ];
5773 list< const SMDS_MeshElement* > eSide[ 2 ];
5774 list< const SMDS_MeshNode* >::iterator nIt[ 2 ];
5775 list< const SMDS_MeshElement* >::iterator eIt[ 2 ];
5779 if (!FindFreeBorder(theBordFirstNode,theBordSecondNode,theBordLastNode,
5780 nSide[0], eSide[0])) {
5781 MESSAGE(" Free Border 1 not found " );
5782 aResult = SEW_BORDER1_NOT_FOUND;
5784 if (theSideIsFreeBorder) {
5787 if (!FindFreeBorder(theSideFirstNode, theSideSecondNode, theSideThirdNode,
5788 nSide[1], eSide[1])) {
5789 MESSAGE(" Free Border 2 not found " );
5790 aResult = ( aResult != SEW_OK ? SEW_BOTH_BORDERS_NOT_FOUND : SEW_BORDER2_NOT_FOUND );
5793 if ( aResult != SEW_OK )
5796 if (!theSideIsFreeBorder) {
5800 // -------------------------------------------------------------------------
5802 // 1. If nodes to merge are not coincident, move nodes of the free border
5803 // from the coord sys defined by the direction from the first to last
5804 // nodes of the border to the correspondent sys of the side 2
5805 // 2. On the side 2, find the links most co-directed with the correspondent
5806 // links of the free border
5807 // -------------------------------------------------------------------------
5809 // 1. Since sewing may brake if there are volumes to split on the side 2,
5810 // we wont move nodes but just compute new coordinates for them
5811 typedef map<const SMDS_MeshNode*, gp_XYZ> TNodeXYZMap;
5812 TNodeXYZMap nBordXYZ;
5813 list< const SMDS_MeshNode* >& bordNodes = nSide[ 0 ];
5814 list< const SMDS_MeshNode* >::iterator nBordIt;
5816 gp_XYZ Pb1( theBordFirstNode->X(), theBordFirstNode->Y(), theBordFirstNode->Z() );
5817 gp_XYZ Pb2( theBordLastNode->X(), theBordLastNode->Y(), theBordLastNode->Z() );
5818 gp_XYZ Ps1( theSideFirstNode->X(), theSideFirstNode->Y(), theSideFirstNode->Z() );
5819 gp_XYZ Ps2( theSideSecondNode->X(), theSideSecondNode->Y(), theSideSecondNode->Z() );
5820 double tol2 = 1.e-8;
5821 gp_Vec Vbs1( Pb1 - Ps1 ),Vbs2( Pb2 - Ps2 );
5822 if ( Vbs1.SquareMagnitude() > tol2 || Vbs2.SquareMagnitude() > tol2 ) {
5823 // Need node movement.
5825 // find X and Z axes to create trsf
5826 gp_Vec Zb( Pb1 - Pb2 ), Zs( Ps1 - Ps2 );
5828 if ( X.SquareMagnitude() <= gp::Resolution() * gp::Resolution() )
5830 X = gp_Ax2( gp::Origin(), Zb ).XDirection();
5833 gp_Ax3 toBordAx( Pb1, Zb, X );
5834 gp_Ax3 fromSideAx( Ps1, Zs, X );
5835 gp_Ax3 toGlobalAx( gp::Origin(), gp::DZ(), gp::DX() );
5837 gp_Trsf toBordSys, fromSide2Sys;
5838 toBordSys.SetTransformation( toBordAx );
5839 fromSide2Sys.SetTransformation( fromSideAx, toGlobalAx );
5840 fromSide2Sys.SetScaleFactor( Zs.Magnitude() / Zb.Magnitude() );
5843 for ( nBordIt = bordNodes.begin(); nBordIt != bordNodes.end(); nBordIt++ ) {
5844 const SMDS_MeshNode* n = *nBordIt;
5845 gp_XYZ xyz( n->X(),n->Y(),n->Z() );
5846 toBordSys.Transforms( xyz );
5847 fromSide2Sys.Transforms( xyz );
5848 nBordXYZ.insert( TNodeXYZMap::value_type( n, xyz ));
5852 // just insert nodes XYZ in the nBordXYZ map
5853 for ( nBordIt = bordNodes.begin(); nBordIt != bordNodes.end(); nBordIt++ ) {
5854 const SMDS_MeshNode* n = *nBordIt;
5855 nBordXYZ.insert( TNodeXYZMap::value_type( n, gp_XYZ( n->X(),n->Y(),n->Z() )));
5859 // 2. On the side 2, find the links most co-directed with the correspondent
5860 // links of the free border
5862 list< const SMDS_MeshElement* >& sideElems = eSide[ 1 ];
5863 list< const SMDS_MeshNode* >& sideNodes = nSide[ 1 ];
5864 sideNodes.push_back( theSideFirstNode );
5866 bool hasVolumes = false;
5867 LinkID_Gen aLinkID_Gen( GetMeshDS() );
5868 set<long> foundSideLinkIDs, checkedLinkIDs;
5869 SMDS_VolumeTool volume;
5870 //const SMDS_MeshNode* faceNodes[ 4 ];
5872 const SMDS_MeshNode* sideNode;
5873 const SMDS_MeshElement* sideElem;
5874 const SMDS_MeshNode* prevSideNode = theSideFirstNode;
5875 const SMDS_MeshNode* prevBordNode = theBordFirstNode;
5876 nBordIt = bordNodes.begin();
5878 // border node position and border link direction to compare with
5879 gp_XYZ bordPos = nBordXYZ[ *nBordIt ];
5880 gp_XYZ bordDir = bordPos - nBordXYZ[ prevBordNode ];
5881 // choose next side node by link direction or by closeness to
5882 // the current border node:
5883 bool searchByDir = ( *nBordIt != theBordLastNode );
5885 // find the next node on the Side 2
5887 double maxDot = -DBL_MAX, minDist = DBL_MAX;
5889 checkedLinkIDs.clear();
5890 gp_XYZ prevXYZ( prevSideNode->X(), prevSideNode->Y(), prevSideNode->Z() );
5892 // loop on inverse elements of current node (prevSideNode) on the Side 2
5893 SMDS_ElemIteratorPtr invElemIt = prevSideNode->GetInverseElementIterator();
5894 while ( invElemIt->more() )
5896 const SMDS_MeshElement* elem = invElemIt->next();
5897 // prepare data for a loop on links coming to prevSideNode, of a face or a volume
5898 int iPrevNode, iNode = 0, nbNodes = elem->NbNodes();
5899 vector< const SMDS_MeshNode* > faceNodes( nbNodes, (const SMDS_MeshNode*)0 );
5900 bool isVolume = volume.Set( elem );
5901 const SMDS_MeshNode** nodes = isVolume ? volume.GetNodes() : & faceNodes[0];
5902 if ( isVolume ) // --volume
5904 else if ( elem->GetType()==SMDSAbs_Face ) { // --face
5905 // retrieve all face nodes and find iPrevNode - an index of the prevSideNode
5906 if(elem->IsQuadratic()) {
5907 const SMDS_QuadraticFaceOfNodes* F =
5908 static_cast<const SMDS_QuadraticFaceOfNodes*>(elem);
5909 // use special nodes iterator
5910 SMDS_NodeIteratorPtr anIter = F->interlacedNodesIterator();
5911 while( anIter->more() ) {
5912 nodes[ iNode ] = anIter->next();
5913 if ( nodes[ iNode++ ] == prevSideNode )
5914 iPrevNode = iNode - 1;
5918 SMDS_ElemIteratorPtr nIt = elem->nodesIterator();
5919 while ( nIt->more() ) {
5920 nodes[ iNode ] = cast2Node( nIt->next() );
5921 if ( nodes[ iNode++ ] == prevSideNode )
5922 iPrevNode = iNode - 1;
5925 // there are 2 links to check
5930 // loop on links, to be precise, on the second node of links
5931 for ( iNode = 0; iNode < nbNodes; iNode++ ) {
5932 const SMDS_MeshNode* n = nodes[ iNode ];
5934 if ( !volume.IsLinked( n, prevSideNode ))
5938 if ( iNode ) // a node before prevSideNode
5939 n = nodes[ iPrevNode == 0 ? elem->NbNodes() - 1 : iPrevNode - 1 ];
5940 else // a node after prevSideNode
5941 n = nodes[ iPrevNode + 1 == elem->NbNodes() ? 0 : iPrevNode + 1 ];
5943 // check if this link was already used
5944 long iLink = aLinkID_Gen.GetLinkID( prevSideNode, n );
5945 bool isJustChecked = !checkedLinkIDs.insert( iLink ).second;
5946 if (!isJustChecked &&
5947 foundSideLinkIDs.find( iLink ) == foundSideLinkIDs.end() )
5949 // test a link geometrically
5950 gp_XYZ nextXYZ ( n->X(), n->Y(), n->Z() );
5951 bool linkIsBetter = false;
5952 double dot = 0.0, dist = 0.0;
5953 if ( searchByDir ) { // choose most co-directed link
5954 dot = bordDir * ( nextXYZ - prevXYZ ).Normalized();
5955 linkIsBetter = ( dot > maxDot );
5957 else { // choose link with the node closest to bordPos
5958 dist = ( nextXYZ - bordPos ).SquareModulus();
5959 linkIsBetter = ( dist < minDist );
5961 if ( linkIsBetter ) {
5970 } // loop on inverse elements of prevSideNode
5973 MESSAGE(" Cant find path by links of the Side 2 ");
5974 return SEW_BAD_SIDE_NODES;
5976 sideNodes.push_back( sideNode );
5977 sideElems.push_back( sideElem );
5978 foundSideLinkIDs.insert ( linkID );
5979 prevSideNode = sideNode;
5981 if ( *nBordIt == theBordLastNode )
5982 searchByDir = false;
5984 // find the next border link to compare with
5985 gp_XYZ sidePos( sideNode->X(), sideNode->Y(), sideNode->Z() );
5986 searchByDir = ( bordDir * ( sidePos - bordPos ) <= 0 );
5987 // move to next border node if sideNode is before forward border node (bordPos)
5988 while ( *nBordIt != theBordLastNode && !searchByDir ) {
5989 prevBordNode = *nBordIt;
5991 bordPos = nBordXYZ[ *nBordIt ];
5992 bordDir = bordPos - nBordXYZ[ prevBordNode ];
5993 searchByDir = ( bordDir * ( sidePos - bordPos ) <= 0 );
5997 while ( sideNode != theSideSecondNode );
5999 if ( hasVolumes && sideNodes.size () != bordNodes.size() && !toCreatePolyedrs) {
6000 MESSAGE("VOLUME SPLITTING IS FORBIDDEN");
6001 return SEW_VOLUMES_TO_SPLIT; // volume splitting is forbidden
6003 } // end nodes search on the side 2
6005 // ============================
6006 // sew the border to the side 2
6007 // ============================
6009 int nbNodes[] = { nSide[0].size(), nSide[1].size() };
6010 int maxNbNodes = Max( nbNodes[0], nbNodes[1] );
6012 TListOfListOfNodes nodeGroupsToMerge;
6013 if ( nbNodes[0] == nbNodes[1] ||
6014 ( theSideIsFreeBorder && !theSideThirdNode)) {
6016 // all nodes are to be merged
6018 for (nIt[0] = nSide[0].begin(), nIt[1] = nSide[1].begin();
6019 nIt[0] != nSide[0].end() && nIt[1] != nSide[1].end();
6020 nIt[0]++, nIt[1]++ )
6022 nodeGroupsToMerge.push_back( list<const SMDS_MeshNode*>() );
6023 nodeGroupsToMerge.back().push_back( *nIt[1] ); // to keep
6024 nodeGroupsToMerge.back().push_back( *nIt[0] ); // to remove
6029 // insert new nodes into the border and the side to get equal nb of segments
6031 // get normalized parameters of nodes on the borders
6032 //double param[ 2 ][ maxNbNodes ];
6034 param[0] = new double [ maxNbNodes ];
6035 param[1] = new double [ maxNbNodes ];
6037 for ( iBord = 0; iBord < 2; iBord++ ) { // loop on 2 borders
6038 list< const SMDS_MeshNode* >& nodes = nSide[ iBord ];
6039 list< const SMDS_MeshNode* >::iterator nIt = nodes.begin();
6040 const SMDS_MeshNode* nPrev = *nIt;
6041 double bordLength = 0;
6042 for ( iNode = 0; nIt != nodes.end(); nIt++, iNode++ ) { // loop on border nodes
6043 const SMDS_MeshNode* nCur = *nIt;
6044 gp_XYZ segment (nCur->X() - nPrev->X(),
6045 nCur->Y() - nPrev->Y(),
6046 nCur->Z() - nPrev->Z());
6047 double segmentLen = segment.Modulus();
6048 bordLength += segmentLen;
6049 param[ iBord ][ iNode ] = bordLength;
6052 // normalize within [0,1]
6053 for ( iNode = 0; iNode < nbNodes[ iBord ]; iNode++ ) {
6054 param[ iBord ][ iNode ] /= bordLength;
6058 // loop on border segments
6059 const SMDS_MeshNode *nPrev[ 2 ] = { 0, 0 };
6060 int i[ 2 ] = { 0, 0 };
6061 nIt[0] = nSide[0].begin(); eIt[0] = eSide[0].begin();
6062 nIt[1] = nSide[1].begin(); eIt[1] = eSide[1].begin();
6064 TElemOfNodeListMap insertMap;
6065 TElemOfNodeListMap::iterator insertMapIt;
6067 // key: elem to insert nodes into
6068 // value: 2 nodes to insert between + nodes to be inserted
6070 bool next[ 2 ] = { false, false };
6072 // find min adjacent segment length after sewing
6073 double nextParam = 10., prevParam = 0;
6074 for ( iBord = 0; iBord < 2; iBord++ ) { // loop on 2 borders
6075 if ( i[ iBord ] + 1 < nbNodes[ iBord ])
6076 nextParam = Min( nextParam, param[iBord][ i[iBord] + 1 ]);
6077 if ( i[ iBord ] > 0 )
6078 prevParam = Max( prevParam, param[iBord][ i[iBord] - 1 ]);
6080 double minParam = Min( param[ 0 ][ i[0] ], param[ 1 ][ i[1] ]);
6081 double maxParam = Max( param[ 0 ][ i[0] ], param[ 1 ][ i[1] ]);
6082 double minSegLen = Min( nextParam - minParam, maxParam - prevParam );
6084 // choose to insert or to merge nodes
6085 double du = param[ 1 ][ i[1] ] - param[ 0 ][ i[0] ];
6086 if ( Abs( du ) <= minSegLen * 0.2 ) {
6089 nodeGroupsToMerge.push_back( list<const SMDS_MeshNode*>() );
6090 const SMDS_MeshNode* n0 = *nIt[0];
6091 const SMDS_MeshNode* n1 = *nIt[1];
6092 nodeGroupsToMerge.back().push_back( n1 );
6093 nodeGroupsToMerge.back().push_back( n0 );
6094 // position of node of the border changes due to merge
6095 param[ 0 ][ i[0] ] += du;
6096 // move n1 for the sake of elem shape evaluation during insertion.
6097 // n1 will be removed by MergeNodes() anyway
6098 const_cast<SMDS_MeshNode*>( n0 )->setXYZ( n1->X(), n1->Y(), n1->Z() );
6099 next[0] = next[1] = true;
6104 int intoBord = ( du < 0 ) ? 0 : 1;
6105 const SMDS_MeshElement* elem = *eIt[ intoBord ];
6106 const SMDS_MeshNode* n1 = nPrev[ intoBord ];
6107 const SMDS_MeshNode* n2 = *nIt[ intoBord ];
6108 const SMDS_MeshNode* nIns = *nIt[ 1 - intoBord ];
6109 if ( intoBord == 1 ) {
6110 // move node of the border to be on a link of elem of the side
6111 gp_XYZ p1 (n1->X(), n1->Y(), n1->Z());
6112 gp_XYZ p2 (n2->X(), n2->Y(), n2->Z());
6113 double ratio = du / ( param[ 1 ][ i[1] ] - param[ 1 ][ i[1]-1 ]);
6114 gp_XYZ p = p2 * ( 1 - ratio ) + p1 * ratio;
6115 GetMeshDS()->MoveNode( nIns, p.X(), p.Y(), p.Z() );
6117 insertMapIt = insertMap.find( elem );
6118 bool notFound = ( insertMapIt == insertMap.end() );
6119 bool otherLink = ( !notFound && (*insertMapIt).second.front() != n1 );
6121 // insert into another link of the same element:
6122 // 1. perform insertion into the other link of the elem
6123 list<const SMDS_MeshNode*> & nodeList = (*insertMapIt).second;
6124 const SMDS_MeshNode* n12 = nodeList.front(); nodeList.pop_front();
6125 const SMDS_MeshNode* n22 = nodeList.front(); nodeList.pop_front();
6126 InsertNodesIntoLink( elem, n12, n22, nodeList, toCreatePolygons );
6127 // 2. perform insertion into the link of adjacent faces
6129 const SMDS_MeshElement* adjElem = findAdjacentFace( n12, n22, elem );
6131 InsertNodesIntoLink( adjElem, n12, n22, nodeList, toCreatePolygons );
6135 if (toCreatePolyedrs) {
6136 // perform insertion into the links of adjacent volumes
6137 UpdateVolumes(n12, n22, nodeList);
6139 // 3. find an element appeared on n1 and n2 after the insertion
6140 insertMap.erase( elem );
6141 elem = findAdjacentFace( n1, n2, 0 );
6143 if ( notFound || otherLink ) {
6144 // add element and nodes of the side into the insertMap
6145 insertMapIt = insertMap.insert
6146 ( TElemOfNodeListMap::value_type( elem, list<const SMDS_MeshNode*>() )).first;
6147 (*insertMapIt).second.push_back( n1 );
6148 (*insertMapIt).second.push_back( n2 );
6150 // add node to be inserted into elem
6151 (*insertMapIt).second.push_back( nIns );
6152 next[ 1 - intoBord ] = true;
6155 // go to the next segment
6156 for ( iBord = 0; iBord < 2; iBord++ ) { // loop on 2 borders
6157 if ( next[ iBord ] ) {
6158 if ( i[ iBord ] != 0 && eIt[ iBord ] != eSide[ iBord ].end())
6160 nPrev[ iBord ] = *nIt[ iBord ];
6161 nIt[ iBord ]++; i[ iBord ]++;
6165 while ( nIt[0] != nSide[0].end() && nIt[1] != nSide[1].end());
6167 // perform insertion of nodes into elements
6169 for (insertMapIt = insertMap.begin();
6170 insertMapIt != insertMap.end();
6173 const SMDS_MeshElement* elem = (*insertMapIt).first;
6174 list<const SMDS_MeshNode*> & nodeList = (*insertMapIt).second;
6175 const SMDS_MeshNode* n1 = nodeList.front(); nodeList.pop_front();
6176 const SMDS_MeshNode* n2 = nodeList.front(); nodeList.pop_front();
6178 InsertNodesIntoLink( elem, n1, n2, nodeList, toCreatePolygons );
6180 if ( !theSideIsFreeBorder ) {
6181 // look for and insert nodes into the faces adjacent to elem
6183 const SMDS_MeshElement* adjElem = findAdjacentFace( n1, n2, elem );
6185 InsertNodesIntoLink( adjElem, n1, n2, nodeList, toCreatePolygons );
6190 if (toCreatePolyedrs) {
6191 // perform insertion into the links of adjacent volumes
6192 UpdateVolumes(n1, n2, nodeList);
6198 } // end: insert new nodes
6200 MergeNodes ( nodeGroupsToMerge );
6205 //=======================================================================
6206 //function : InsertNodesIntoLink
6207 //purpose : insert theNodesToInsert into theFace between theBetweenNode1
6208 // and theBetweenNode2 and split theElement
6209 //=======================================================================
6211 void SMESH_MeshEditor::InsertNodesIntoLink(const SMDS_MeshElement* theFace,
6212 const SMDS_MeshNode* theBetweenNode1,
6213 const SMDS_MeshNode* theBetweenNode2,
6214 list<const SMDS_MeshNode*>& theNodesToInsert,
6215 const bool toCreatePoly)
6217 if ( theFace->GetType() != SMDSAbs_Face ) return;
6219 // find indices of 2 link nodes and of the rest nodes
6220 int iNode = 0, il1, il2, i3, i4;
6221 il1 = il2 = i3 = i4 = -1;
6222 //const SMDS_MeshNode* nodes[ theFace->NbNodes() ];
6223 vector<const SMDS_MeshNode*> nodes( theFace->NbNodes() );
6225 if(theFace->IsQuadratic()) {
6226 const SMDS_QuadraticFaceOfNodes* F =
6227 static_cast<const SMDS_QuadraticFaceOfNodes*>(theFace);
6228 // use special nodes iterator
6229 SMDS_NodeIteratorPtr anIter = F->interlacedNodesIterator();
6230 while( anIter->more() ) {
6231 const SMDS_MeshNode* n = anIter->next();
6232 if ( n == theBetweenNode1 )
6234 else if ( n == theBetweenNode2 )
6240 nodes[ iNode++ ] = n;
6244 SMDS_ElemIteratorPtr nodeIt = theFace->nodesIterator();
6245 while ( nodeIt->more() ) {
6246 const SMDS_MeshNode* n = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
6247 if ( n == theBetweenNode1 )
6249 else if ( n == theBetweenNode2 )
6255 nodes[ iNode++ ] = n;
6258 if ( il1 < 0 || il2 < 0 || i3 < 0 )
6261 // arrange link nodes to go one after another regarding the face orientation
6262 bool reverse = ( Abs( il2 - il1 ) == 1 ? il2 < il1 : il1 < il2 );
6263 list<const SMDS_MeshNode *> aNodesToInsert = theNodesToInsert;
6268 aNodesToInsert.reverse();
6270 // check that not link nodes of a quadrangles are in good order
6271 int nbFaceNodes = theFace->NbNodes();
6272 if ( nbFaceNodes == 4 && i4 - i3 != 1 ) {
6278 if (toCreatePoly || theFace->IsPoly()) {
6281 vector<const SMDS_MeshNode *> poly_nodes (nbFaceNodes + aNodesToInsert.size());
6283 // add nodes of face up to first node of link
6286 if(theFace->IsQuadratic()) {
6287 const SMDS_QuadraticFaceOfNodes* F =
6288 static_cast<const SMDS_QuadraticFaceOfNodes*>(theFace);
6289 // use special nodes iterator
6290 SMDS_NodeIteratorPtr anIter = F->interlacedNodesIterator();
6291 while( anIter->more() && !isFLN ) {
6292 const SMDS_MeshNode* n = anIter->next();
6293 poly_nodes[iNode++] = n;
6294 if (n == nodes[il1]) {
6298 // add nodes to insert
6299 list<const SMDS_MeshNode*>::iterator nIt = aNodesToInsert.begin();
6300 for (; nIt != aNodesToInsert.end(); nIt++) {
6301 poly_nodes[iNode++] = *nIt;
6303 // add nodes of face starting from last node of link
6304 while ( anIter->more() ) {
6305 poly_nodes[iNode++] = anIter->next();
6309 SMDS_ElemIteratorPtr nodeIt = theFace->nodesIterator();
6310 while ( nodeIt->more() && !isFLN ) {
6311 const SMDS_MeshNode* n = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
6312 poly_nodes[iNode++] = n;
6313 if (n == nodes[il1]) {
6317 // add nodes to insert
6318 list<const SMDS_MeshNode*>::iterator nIt = aNodesToInsert.begin();
6319 for (; nIt != aNodesToInsert.end(); nIt++) {
6320 poly_nodes[iNode++] = *nIt;
6322 // add nodes of face starting from last node of link
6323 while ( nodeIt->more() ) {
6324 const SMDS_MeshNode* n = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
6325 poly_nodes[iNode++] = n;
6329 // edit or replace the face
6330 SMESHDS_Mesh *aMesh = GetMeshDS();
6332 if (theFace->IsPoly()) {
6333 aMesh->ChangePolygonNodes(theFace, poly_nodes);
6336 int aShapeId = FindShape( theFace );
6338 SMDS_MeshElement* newElem = aMesh->AddPolygonalFace(poly_nodes);
6339 myLastCreatedElems.Append(newElem);
6340 if ( aShapeId && newElem )
6341 aMesh->SetMeshElementOnShape( newElem, aShapeId );
6343 aMesh->RemoveElement(theFace);
6348 if( !theFace->IsQuadratic() ) {
6350 // put aNodesToInsert between theBetweenNode1 and theBetweenNode2
6351 int nbLinkNodes = 2 + aNodesToInsert.size();
6352 //const SMDS_MeshNode* linkNodes[ nbLinkNodes ];
6353 vector<const SMDS_MeshNode*> linkNodes( nbLinkNodes );
6354 linkNodes[ 0 ] = nodes[ il1 ];
6355 linkNodes[ nbLinkNodes - 1 ] = nodes[ il2 ];
6356 list<const SMDS_MeshNode*>::iterator nIt = aNodesToInsert.begin();
6357 for ( iNode = 1; nIt != aNodesToInsert.end(); nIt++ ) {
6358 linkNodes[ iNode++ ] = *nIt;
6360 // decide how to split a quadrangle: compare possible variants
6361 // and choose which of splits to be a quadrangle
6362 int i1, i2, iSplit, nbSplits = nbLinkNodes - 1, iBestQuad;
6363 if ( nbFaceNodes == 3 ) {
6364 iBestQuad = nbSplits;
6367 else if ( nbFaceNodes == 4 ) {
6368 SMESH::Controls::NumericalFunctorPtr aCrit( new SMESH::Controls::AspectRatio);
6369 double aBestRate = DBL_MAX;
6370 for ( int iQuad = 0; iQuad < nbSplits; iQuad++ ) {
6372 double aBadRate = 0;
6373 // evaluate elements quality
6374 for ( iSplit = 0; iSplit < nbSplits; iSplit++ ) {
6375 if ( iSplit == iQuad ) {
6376 SMDS_FaceOfNodes quad (linkNodes[ i1++ ],
6380 aBadRate += getBadRate( &quad, aCrit );
6383 SMDS_FaceOfNodes tria (linkNodes[ i1++ ],
6385 nodes[ iSplit < iQuad ? i4 : i3 ]);
6386 aBadRate += getBadRate( &tria, aCrit );
6390 if ( aBadRate < aBestRate ) {
6392 aBestRate = aBadRate;
6397 // create new elements
6398 SMESHDS_Mesh *aMesh = GetMeshDS();
6399 int aShapeId = FindShape( theFace );
6402 for ( iSplit = 0; iSplit < nbSplits - 1; iSplit++ ) {
6403 SMDS_MeshElement* newElem = 0;
6404 if ( iSplit == iBestQuad )
6405 newElem = aMesh->AddFace (linkNodes[ i1++ ],
6410 newElem = aMesh->AddFace (linkNodes[ i1++ ],
6412 nodes[ iSplit < iBestQuad ? i4 : i3 ]);
6413 myLastCreatedElems.Append(newElem);
6414 if ( aShapeId && newElem )
6415 aMesh->SetMeshElementOnShape( newElem, aShapeId );
6418 // change nodes of theFace
6419 const SMDS_MeshNode* newNodes[ 4 ];
6420 newNodes[ 0 ] = linkNodes[ i1 ];
6421 newNodes[ 1 ] = linkNodes[ i2 ];
6422 newNodes[ 2 ] = nodes[ iSplit >= iBestQuad ? i3 : i4 ];
6423 newNodes[ 3 ] = nodes[ i4 ];
6424 aMesh->ChangeElementNodes( theFace, newNodes, iSplit == iBestQuad ? 4 : 3 );
6425 } // end if(!theFace->IsQuadratic())
6426 else { // theFace is quadratic
6427 // we have to split theFace on simple triangles and one simple quadrangle
6429 int nbshift = tmp*2;
6430 // shift nodes in nodes[] by nbshift
6432 for(i=0; i<nbshift; i++) {
6433 const SMDS_MeshNode* n = nodes[0];
6434 for(j=0; j<nbFaceNodes-1; j++) {
6435 nodes[j] = nodes[j+1];
6437 nodes[nbFaceNodes-1] = n;
6439 il1 = il1 - nbshift;
6440 // now have to insert nodes between n0 and n1 or n1 and n2 (see below)
6441 // n0 n1 n2 n0 n1 n2
6442 // +-----+-----+ +-----+-----+
6451 // create new elements
6452 SMESHDS_Mesh *aMesh = GetMeshDS();
6453 int aShapeId = FindShape( theFace );
6456 if(nbFaceNodes==6) { // quadratic triangle
6457 SMDS_MeshElement* newElem =
6458 aMesh->AddFace(nodes[3],nodes[4],nodes[5]);
6459 myLastCreatedElems.Append(newElem);
6460 if ( aShapeId && newElem )
6461 aMesh->SetMeshElementOnShape( newElem, aShapeId );
6462 if(theFace->IsMediumNode(nodes[il1])) {
6463 // create quadrangle
6464 newElem = aMesh->AddFace(nodes[0],nodes[1],nodes[3],nodes[5]);
6465 myLastCreatedElems.Append(newElem);
6466 if ( aShapeId && newElem )
6467 aMesh->SetMeshElementOnShape( newElem, aShapeId );
6473 // create quadrangle
6474 newElem = aMesh->AddFace(nodes[1],nodes[2],nodes[3],nodes[5]);
6475 myLastCreatedElems.Append(newElem);
6476 if ( aShapeId && newElem )
6477 aMesh->SetMeshElementOnShape( newElem, aShapeId );
6483 else { // nbFaceNodes==8 - quadratic quadrangle
6484 SMDS_MeshElement* newElem =
6485 aMesh->AddFace(nodes[3],nodes[4],nodes[5]);
6486 myLastCreatedElems.Append(newElem);
6487 if ( aShapeId && newElem )
6488 aMesh->SetMeshElementOnShape( newElem, aShapeId );
6489 newElem = aMesh->AddFace(nodes[5],nodes[6],nodes[7]);
6490 myLastCreatedElems.Append(newElem);
6491 if ( aShapeId && newElem )
6492 aMesh->SetMeshElementOnShape( newElem, aShapeId );
6493 newElem = aMesh->AddFace(nodes[5],nodes[7],nodes[3]);
6494 myLastCreatedElems.Append(newElem);
6495 if ( aShapeId && newElem )
6496 aMesh->SetMeshElementOnShape( newElem, aShapeId );
6497 if(theFace->IsMediumNode(nodes[il1])) {
6498 // create quadrangle
6499 newElem = aMesh->AddFace(nodes[0],nodes[1],nodes[3],nodes[7]);
6500 myLastCreatedElems.Append(newElem);
6501 if ( aShapeId && newElem )
6502 aMesh->SetMeshElementOnShape( newElem, aShapeId );
6508 // create quadrangle
6509 newElem = aMesh->AddFace(nodes[1],nodes[2],nodes[3],nodes[7]);
6510 myLastCreatedElems.Append(newElem);
6511 if ( aShapeId && newElem )
6512 aMesh->SetMeshElementOnShape( newElem, aShapeId );
6518 // create needed triangles using n1,n2,n3 and inserted nodes
6519 int nbn = 2 + aNodesToInsert.size();
6520 //const SMDS_MeshNode* aNodes[nbn];
6521 vector<const SMDS_MeshNode*> aNodes(nbn);
6522 aNodes[0] = nodes[n1];
6523 aNodes[nbn-1] = nodes[n2];
6524 list<const SMDS_MeshNode*>::iterator nIt = aNodesToInsert.begin();
6525 for ( iNode = 1; nIt != aNodesToInsert.end(); nIt++ ) {
6526 aNodes[iNode++] = *nIt;
6528 for(i=1; i<nbn; i++) {
6529 SMDS_MeshElement* newElem =
6530 aMesh->AddFace(aNodes[i-1],aNodes[i],nodes[n3]);
6531 myLastCreatedElems.Append(newElem);
6532 if ( aShapeId && newElem )
6533 aMesh->SetMeshElementOnShape( newElem, aShapeId );
6535 // remove old quadratic face
6536 aMesh->RemoveElement(theFace);
6540 //=======================================================================
6541 //function : UpdateVolumes
6543 //=======================================================================
6544 void SMESH_MeshEditor::UpdateVolumes (const SMDS_MeshNode* theBetweenNode1,
6545 const SMDS_MeshNode* theBetweenNode2,
6546 list<const SMDS_MeshNode*>& theNodesToInsert)
6548 myLastCreatedElems.Clear();
6549 myLastCreatedNodes.Clear();
6551 SMDS_ElemIteratorPtr invElemIt = theBetweenNode1->GetInverseElementIterator(SMDSAbs_Volume);
6552 while (invElemIt->more()) { // loop on inverse elements of theBetweenNode1
6553 const SMDS_MeshElement* elem = invElemIt->next();
6555 // check, if current volume has link theBetweenNode1 - theBetweenNode2
6556 SMDS_VolumeTool aVolume (elem);
6557 if (!aVolume.IsLinked(theBetweenNode1, theBetweenNode2))
6560 // insert new nodes in all faces of the volume, sharing link theBetweenNode1 - theBetweenNode2
6561 int iface, nbFaces = aVolume.NbFaces();
6562 vector<const SMDS_MeshNode *> poly_nodes;
6563 vector<int> quantities (nbFaces);
6565 for (iface = 0; iface < nbFaces; iface++) {
6566 int nbFaceNodes = aVolume.NbFaceNodes(iface), nbInserted = 0;
6567 // faceNodes will contain (nbFaceNodes + 1) nodes, last = first
6568 const SMDS_MeshNode** faceNodes = aVolume.GetFaceNodes(iface);
6570 for (int inode = 0; inode < nbFaceNodes; inode++) {
6571 poly_nodes.push_back(faceNodes[inode]);
6573 if (nbInserted == 0) {
6574 if (faceNodes[inode] == theBetweenNode1) {
6575 if (faceNodes[inode + 1] == theBetweenNode2) {
6576 nbInserted = theNodesToInsert.size();
6578 // add nodes to insert
6579 list<const SMDS_MeshNode*>::iterator nIt = theNodesToInsert.begin();
6580 for (; nIt != theNodesToInsert.end(); nIt++) {
6581 poly_nodes.push_back(*nIt);
6585 else if (faceNodes[inode] == theBetweenNode2) {
6586 if (faceNodes[inode + 1] == theBetweenNode1) {
6587 nbInserted = theNodesToInsert.size();
6589 // add nodes to insert in reversed order
6590 list<const SMDS_MeshNode*>::iterator nIt = theNodesToInsert.end();
6592 for (; nIt != theNodesToInsert.begin(); nIt--) {
6593 poly_nodes.push_back(*nIt);
6595 poly_nodes.push_back(*nIt);
6602 quantities[iface] = nbFaceNodes + nbInserted;
6605 // Replace or update the volume
6606 SMESHDS_Mesh *aMesh = GetMeshDS();
6608 if (elem->IsPoly()) {
6609 aMesh->ChangePolyhedronNodes(elem, poly_nodes, quantities);
6613 int aShapeId = FindShape( elem );
6615 SMDS_MeshElement* newElem =
6616 aMesh->AddPolyhedralVolume(poly_nodes, quantities);
6617 myLastCreatedElems.Append(newElem);
6618 if (aShapeId && newElem)
6619 aMesh->SetMeshElementOnShape(newElem, aShapeId);
6621 aMesh->RemoveElement(elem);
6626 //=======================================================================
6628 * \brief Convert elements contained in a submesh to quadratic
6629 * \retval int - nb of checked elements
6631 //=======================================================================
6633 int SMESH_MeshEditor::convertElemToQuadratic(SMESHDS_SubMesh * theSm,
6634 SMESH_MesherHelper& theHelper,
6635 const bool theForce3d)
6638 if( !theSm ) return nbElem;
6640 const bool notFromGroups = false;
6641 SMDS_ElemIteratorPtr ElemItr = theSm->GetElements();
6642 while(ElemItr->more())
6645 const SMDS_MeshElement* elem = ElemItr->next();
6646 if( !elem || elem->IsQuadratic() ) continue;
6648 int id = elem->GetID();
6649 int nbNodes = elem->NbNodes();
6650 vector<const SMDS_MeshNode *> aNds (nbNodes);
6652 for(int i = 0; i < nbNodes; i++)
6654 aNds[i] = elem->GetNode(i);
6656 SMDSAbs_ElementType aType = elem->GetType();
6658 GetMeshDS()->RemoveFreeElement(elem, theSm, notFromGroups);
6660 const SMDS_MeshElement* NewElem = 0;
6666 NewElem = theHelper.AddEdge(aNds[0], aNds[1], id, theForce3d);
6674 NewElem = theHelper.AddFace(aNds[0], aNds[1], aNds[2], id, theForce3d);
6677 NewElem = theHelper.AddFace(aNds[0], aNds[1], aNds[2], aNds[3], id, theForce3d);
6684 case SMDSAbs_Volume :
6689 NewElem = theHelper.AddVolume(aNds[0], aNds[1], aNds[2], aNds[3], id, theForce3d);
6692 NewElem = theHelper.AddVolume(aNds[0], aNds[1], aNds[2], aNds[3], aNds[4], aNds[5], id, theForce3d);
6695 NewElem = theHelper.AddVolume(aNds[0], aNds[1], aNds[2], aNds[3],
6696 aNds[4], aNds[5], aNds[6], aNds[7], id, theForce3d);
6706 ReplaceElemInGroups( elem, NewElem, GetMeshDS());
6708 theSm->AddElement( NewElem );
6713 //=======================================================================
6714 //function : ConvertToQuadratic
6716 //=======================================================================
6717 void SMESH_MeshEditor::ConvertToQuadratic(const bool theForce3d)
6719 SMESHDS_Mesh* meshDS = GetMeshDS();
6721 SMESH_MesherHelper aHelper(*myMesh);
6722 aHelper.SetIsQuadratic( true );
6723 const bool notFromGroups = false;
6725 int nbCheckedElems = 0;
6726 if ( myMesh->HasShapeToMesh() )
6728 if ( SMESH_subMesh *aSubMesh = myMesh->GetSubMeshContaining(myMesh->GetShapeToMesh()))
6730 SMESH_subMeshIteratorPtr smIt = aSubMesh->getDependsOnIterator(true,false);
6731 while ( smIt->more() ) {
6732 SMESH_subMesh* sm = smIt->next();
6733 if ( SMESHDS_SubMesh *smDS = sm->GetSubMeshDS() ) {
6734 aHelper.SetSubShape( sm->GetSubShape() );
6735 nbCheckedElems += convertElemToQuadratic(smDS, aHelper, theForce3d);
6740 int totalNbElems = meshDS->NbEdges() + meshDS->NbFaces() + meshDS->NbVolumes();
6741 if ( nbCheckedElems < totalNbElems ) // not all elements are in submeshes
6743 SMESHDS_SubMesh *smDS = 0;
6744 SMDS_EdgeIteratorPtr aEdgeItr = meshDS->edgesIterator();
6745 while(aEdgeItr->more())
6747 const SMDS_MeshEdge* edge = aEdgeItr->next();
6748 if(edge && !edge->IsQuadratic())
6750 int id = edge->GetID();
6751 const SMDS_MeshNode* n1 = edge->GetNode(0);
6752 const SMDS_MeshNode* n2 = edge->GetNode(1);
6754 meshDS->RemoveFreeElement(edge, smDS, notFromGroups);
6756 const SMDS_MeshEdge* NewEdge = aHelper.AddEdge(n1, n2, id, theForce3d);
6757 ReplaceElemInGroups( edge, NewEdge, GetMeshDS());
6760 SMDS_FaceIteratorPtr aFaceItr = meshDS->facesIterator();
6761 while(aFaceItr->more())
6763 const SMDS_MeshFace* face = aFaceItr->next();
6764 if(!face || face->IsQuadratic() ) continue;
6766 int id = face->GetID();
6767 int nbNodes = face->NbNodes();
6768 vector<const SMDS_MeshNode *> aNds (nbNodes);
6770 for(int i = 0; i < nbNodes; i++)
6772 aNds[i] = face->GetNode(i);
6775 meshDS->RemoveFreeElement(face, smDS, notFromGroups);
6777 SMDS_MeshFace * NewFace = 0;
6781 NewFace = aHelper.AddFace(aNds[0], aNds[1], aNds[2], id, theForce3d);
6784 NewFace = aHelper.AddFace(aNds[0], aNds[1], aNds[2], aNds[3], id, theForce3d);
6789 ReplaceElemInGroups( face, NewFace, GetMeshDS());
6791 SMDS_VolumeIteratorPtr aVolumeItr = meshDS->volumesIterator();
6792 while(aVolumeItr->more())
6794 const SMDS_MeshVolume* volume = aVolumeItr->next();
6795 if(!volume || volume->IsQuadratic() ) continue;
6797 int id = volume->GetID();
6798 int nbNodes = volume->NbNodes();
6799 vector<const SMDS_MeshNode *> aNds (nbNodes);
6801 for(int i = 0; i < nbNodes; i++)
6803 aNds[i] = volume->GetNode(i);
6806 meshDS->RemoveFreeElement(volume, smDS, notFromGroups);
6808 SMDS_MeshVolume * NewVolume = 0;
6812 NewVolume = aHelper.AddVolume(aNds[0], aNds[1], aNds[2],
6813 aNds[3], id, theForce3d );
6816 NewVolume = aHelper.AddVolume(aNds[0], aNds[1], aNds[2],
6817 aNds[3], aNds[4], aNds[5], id, theForce3d);
6820 NewVolume = aHelper.AddVolume(aNds[0], aNds[1], aNds[2], aNds[3],
6821 aNds[4], aNds[5], aNds[6], aNds[7], id, theForce3d);
6826 ReplaceElemInGroups(volume, NewVolume, meshDS);
6831 //=======================================================================
6833 * \brief Convert quadratic elements to linear ones and remove quadratic nodes
6834 * \retval int - nb of checked elements
6836 //=======================================================================
6838 int SMESH_MeshEditor::removeQuadElem(SMESHDS_SubMesh * theSm,
6839 SMDS_ElemIteratorPtr theItr,
6840 const int theShapeID)
6843 SMESHDS_Mesh* meshDS = GetMeshDS();
6844 const bool notFromGroups = false;
6846 while( theItr->more() )
6848 const SMDS_MeshElement* elem = theItr->next();
6850 if( elem && elem->IsQuadratic())
6852 int id = elem->GetID();
6853 int nbNodes = elem->NbNodes();
6854 vector<const SMDS_MeshNode *> aNds, mediumNodes;
6855 aNds.reserve( nbNodes );
6856 mediumNodes.reserve( nbNodes );
6858 for(int i = 0; i < nbNodes; i++)
6860 const SMDS_MeshNode* n = elem->GetNode(i);
6862 if( elem->IsMediumNode( n ) )
6863 mediumNodes.push_back( n );
6865 aNds.push_back( n );
6867 if( aNds.empty() ) continue;
6868 SMDSAbs_ElementType aType = elem->GetType();
6870 //remove old quadratic element
6871 meshDS->RemoveFreeElement( elem, theSm, notFromGroups );
6873 SMDS_MeshElement * NewElem = AddElement( aNds, aType, false, id );
6874 ReplaceElemInGroups(elem, NewElem, meshDS);
6875 if( theSm && NewElem )
6876 theSm->AddElement( NewElem );
6878 // remove medium nodes
6879 vector<const SMDS_MeshNode*>::iterator nIt = mediumNodes.begin();
6880 for ( ; nIt != mediumNodes.end(); ++nIt ) {
6881 const SMDS_MeshNode* n = *nIt;
6882 if ( n->NbInverseElements() == 0 ) {
6883 if ( n->GetPosition()->GetShapeId() != theShapeID )
6884 meshDS->RemoveFreeNode( n, meshDS->MeshElements
6885 ( n->GetPosition()->GetShapeId() ));
6887 meshDS->RemoveFreeNode( n, theSm );
6895 //=======================================================================
6896 //function : ConvertFromQuadratic
6898 //=======================================================================
6899 bool SMESH_MeshEditor::ConvertFromQuadratic()
6901 int nbCheckedElems = 0;
6902 if ( myMesh->HasShapeToMesh() )
6904 if ( SMESH_subMesh *aSubMesh = myMesh->GetSubMeshContaining(myMesh->GetShapeToMesh()))
6906 SMESH_subMeshIteratorPtr smIt = aSubMesh->getDependsOnIterator(true,false);
6907 while ( smIt->more() ) {
6908 SMESH_subMesh* sm = smIt->next();
6909 if ( SMESHDS_SubMesh *smDS = sm->GetSubMeshDS() )
6910 nbCheckedElems += removeQuadElem( smDS, smDS->GetElements(), sm->GetId() );
6916 GetMeshDS()->NbEdges() + GetMeshDS()->NbFaces() + GetMeshDS()->NbVolumes();
6917 if ( nbCheckedElems < totalNbElems ) // not all elements are in submeshes
6919 SMESHDS_SubMesh *aSM = 0;
6920 removeQuadElem( aSM, GetMeshDS()->elementsIterator(), 0 );
6926 //=======================================================================
6927 //function : SewSideElements
6929 //=======================================================================
6931 SMESH_MeshEditor::Sew_Error
6932 SMESH_MeshEditor::SewSideElements (TIDSortedElemSet& theSide1,
6933 TIDSortedElemSet& theSide2,
6934 const SMDS_MeshNode* theFirstNode1,
6935 const SMDS_MeshNode* theFirstNode2,
6936 const SMDS_MeshNode* theSecondNode1,
6937 const SMDS_MeshNode* theSecondNode2)
6939 myLastCreatedElems.Clear();
6940 myLastCreatedNodes.Clear();
6942 MESSAGE ("::::SewSideElements()");
6943 if ( theSide1.size() != theSide2.size() )
6944 return SEW_DIFF_NB_OF_ELEMENTS;
6946 Sew_Error aResult = SEW_OK;
6948 // 1. Build set of faces representing each side
6949 // 2. Find which nodes of the side 1 to merge with ones on the side 2
6950 // 3. Replace nodes in elements of the side 1 and remove replaced nodes
6952 // =======================================================================
6953 // 1. Build set of faces representing each side:
6954 // =======================================================================
6955 // a. build set of nodes belonging to faces
6956 // b. complete set of faces: find missing fices whose nodes are in set of nodes
6957 // c. create temporary faces representing side of volumes if correspondent
6958 // face does not exist
6960 SMESHDS_Mesh* aMesh = GetMeshDS();
6961 SMDS_Mesh aTmpFacesMesh;
6962 set<const SMDS_MeshElement*> faceSet1, faceSet2;
6963 set<const SMDS_MeshElement*> volSet1, volSet2;
6964 set<const SMDS_MeshNode*> nodeSet1, nodeSet2;
6965 set<const SMDS_MeshElement*> * faceSetPtr[] = { &faceSet1, &faceSet2 };
6966 set<const SMDS_MeshElement*> * volSetPtr[] = { &volSet1, &volSet2 };
6967 set<const SMDS_MeshNode*> * nodeSetPtr[] = { &nodeSet1, &nodeSet2 };
6968 TIDSortedElemSet * elemSetPtr[] = { &theSide1, &theSide2 };
6969 int iSide, iFace, iNode;
6971 for ( iSide = 0; iSide < 2; iSide++ ) {
6972 set<const SMDS_MeshNode*> * nodeSet = nodeSetPtr[ iSide ];
6973 TIDSortedElemSet * elemSet = elemSetPtr[ iSide ];
6974 set<const SMDS_MeshElement*> * faceSet = faceSetPtr[ iSide ];
6975 set<const SMDS_MeshElement*> * volSet = volSetPtr [ iSide ];
6976 set<const SMDS_MeshElement*>::iterator vIt;
6977 TIDSortedElemSet::iterator eIt;
6978 set<const SMDS_MeshNode*>::iterator nIt;
6980 // check that given nodes belong to given elements
6981 const SMDS_MeshNode* n1 = ( iSide == 0 ) ? theFirstNode1 : theFirstNode2;
6982 const SMDS_MeshNode* n2 = ( iSide == 0 ) ? theSecondNode1 : theSecondNode2;
6983 int firstIndex = -1, secondIndex = -1;
6984 for (eIt = elemSet->begin(); eIt != elemSet->end(); eIt++ ) {
6985 const SMDS_MeshElement* elem = *eIt;
6986 if ( firstIndex < 0 ) firstIndex = elem->GetNodeIndex( n1 );
6987 if ( secondIndex < 0 ) secondIndex = elem->GetNodeIndex( n2 );
6988 if ( firstIndex > -1 && secondIndex > -1 ) break;
6990 if ( firstIndex < 0 || secondIndex < 0 ) {
6991 // we can simply return until temporary faces created
6992 return (iSide == 0 ) ? SEW_BAD_SIDE1_NODES : SEW_BAD_SIDE2_NODES;
6995 // -----------------------------------------------------------
6996 // 1a. Collect nodes of existing faces
6997 // and build set of face nodes in order to detect missing
6998 // faces corresponing to sides of volumes
6999 // -----------------------------------------------------------
7001 set< set <const SMDS_MeshNode*> > setOfFaceNodeSet;
7003 // loop on the given element of a side
7004 for (eIt = elemSet->begin(); eIt != elemSet->end(); eIt++ ) {
7005 //const SMDS_MeshElement* elem = *eIt;
7006 const SMDS_MeshElement* elem = *eIt;
7007 if ( elem->GetType() == SMDSAbs_Face ) {
7008 faceSet->insert( elem );
7009 set <const SMDS_MeshNode*> faceNodeSet;
7010 SMDS_ElemIteratorPtr nodeIt = elem->nodesIterator();
7011 while ( nodeIt->more() ) {
7012 const SMDS_MeshNode* n = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
7013 nodeSet->insert( n );
7014 faceNodeSet.insert( n );
7016 setOfFaceNodeSet.insert( faceNodeSet );
7018 else if ( elem->GetType() == SMDSAbs_Volume )
7019 volSet->insert( elem );
7021 // ------------------------------------------------------------------------------
7022 // 1b. Complete set of faces: find missing fices whose nodes are in set of nodes
7023 // ------------------------------------------------------------------------------
7025 for ( nIt = nodeSet->begin(); nIt != nodeSet->end(); nIt++ ) { // loop on nodes of iSide
7026 SMDS_ElemIteratorPtr fIt = (*nIt)->GetInverseElementIterator(SMDSAbs_Face);
7027 while ( fIt->more() ) { // loop on faces sharing a node
7028 const SMDS_MeshElement* f = fIt->next();
7029 if ( faceSet->find( f ) == faceSet->end() ) {
7030 // check if all nodes are in nodeSet and
7031 // complete setOfFaceNodeSet if they are
7032 set <const SMDS_MeshNode*> faceNodeSet;
7033 SMDS_ElemIteratorPtr nodeIt = f->nodesIterator();
7034 bool allInSet = true;
7035 while ( nodeIt->more() && allInSet ) { // loop on nodes of a face
7036 const SMDS_MeshNode* n = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
7037 if ( nodeSet->find( n ) == nodeSet->end() )
7040 faceNodeSet.insert( n );
7043 faceSet->insert( f );
7044 setOfFaceNodeSet.insert( faceNodeSet );
7050 // -------------------------------------------------------------------------
7051 // 1c. Create temporary faces representing sides of volumes if correspondent
7052 // face does not exist
7053 // -------------------------------------------------------------------------
7055 if ( !volSet->empty() ) {
7056 //int nodeSetSize = nodeSet->size();
7058 // loop on given volumes
7059 for ( vIt = volSet->begin(); vIt != volSet->end(); vIt++ ) {
7060 SMDS_VolumeTool vol (*vIt);
7061 // loop on volume faces: find free faces
7062 // --------------------------------------
7063 list<const SMDS_MeshElement* > freeFaceList;
7064 for ( iFace = 0; iFace < vol.NbFaces(); iFace++ ) {
7065 if ( !vol.IsFreeFace( iFace ))
7067 // check if there is already a face with same nodes in a face set
7068 const SMDS_MeshElement* aFreeFace = 0;
7069 const SMDS_MeshNode** fNodes = vol.GetFaceNodes( iFace );
7070 int nbNodes = vol.NbFaceNodes( iFace );
7071 set <const SMDS_MeshNode*> faceNodeSet;
7072 vol.GetFaceNodes( iFace, faceNodeSet );
7073 bool isNewFace = setOfFaceNodeSet.insert( faceNodeSet ).second;
7075 // no such a face is given but it still can exist, check it
7076 if ( nbNodes == 3 ) {
7077 aFreeFace = aMesh->FindFace( fNodes[0],fNodes[1],fNodes[2] );
7079 else if ( nbNodes == 4 ) {
7080 aFreeFace = aMesh->FindFace( fNodes[0],fNodes[1],fNodes[2],fNodes[3] );
7083 vector<const SMDS_MeshNode *> poly_nodes ( fNodes, & fNodes[nbNodes]);
7084 aFreeFace = aMesh->FindFace(poly_nodes);
7088 // create a temporary face
7089 if ( nbNodes == 3 ) {
7090 aFreeFace = aTmpFacesMesh.AddFace( fNodes[0],fNodes[1],fNodes[2] );
7092 else if ( nbNodes == 4 ) {
7093 aFreeFace = aTmpFacesMesh.AddFace( fNodes[0],fNodes[1],fNodes[2],fNodes[3] );
7096 vector<const SMDS_MeshNode *> poly_nodes ( fNodes, & fNodes[nbNodes]);
7097 aFreeFace = aTmpFacesMesh.AddPolygonalFace(poly_nodes);
7101 freeFaceList.push_back( aFreeFace );
7103 } // loop on faces of a volume
7105 // choose one of several free faces
7106 // --------------------------------------
7107 if ( freeFaceList.size() > 1 ) {
7108 // choose a face having max nb of nodes shared by other elems of a side
7109 int maxNbNodes = -1/*, nbExcludedFaces = 0*/;
7110 list<const SMDS_MeshElement* >::iterator fIt = freeFaceList.begin();
7111 while ( fIt != freeFaceList.end() ) { // loop on free faces
7112 int nbSharedNodes = 0;
7113 SMDS_ElemIteratorPtr nodeIt = (*fIt)->nodesIterator();
7114 while ( nodeIt->more() ) { // loop on free face nodes
7115 const SMDS_MeshNode* n =
7116 static_cast<const SMDS_MeshNode*>( nodeIt->next() );
7117 SMDS_ElemIteratorPtr invElemIt = n->GetInverseElementIterator();
7118 while ( invElemIt->more() ) {
7119 const SMDS_MeshElement* e = invElemIt->next();
7120 if ( faceSet->find( e ) != faceSet->end() )
7122 if ( elemSet->find( e ) != elemSet->end() )
7126 if ( nbSharedNodes >= maxNbNodes ) {
7127 maxNbNodes = nbSharedNodes;
7131 freeFaceList.erase( fIt++ ); // here fIt++ occures before erase
7133 if ( freeFaceList.size() > 1 )
7135 // could not choose one face, use another way
7136 // choose a face most close to the bary center of the opposite side
7137 gp_XYZ aBC( 0., 0., 0. );
7138 set <const SMDS_MeshNode*> addedNodes;
7139 TIDSortedElemSet * elemSet2 = elemSetPtr[ 1 - iSide ];
7140 eIt = elemSet2->begin();
7141 for ( eIt = elemSet2->begin(); eIt != elemSet2->end(); eIt++ ) {
7142 SMDS_ElemIteratorPtr nodeIt = (*eIt)->nodesIterator();
7143 while ( nodeIt->more() ) { // loop on free face nodes
7144 const SMDS_MeshNode* n =
7145 static_cast<const SMDS_MeshNode*>( nodeIt->next() );
7146 if ( addedNodes.insert( n ).second )
7147 aBC += gp_XYZ( n->X(),n->Y(),n->Z() );
7150 aBC /= addedNodes.size();
7151 double minDist = DBL_MAX;
7152 fIt = freeFaceList.begin();
7153 while ( fIt != freeFaceList.end() ) { // loop on free faces
7155 SMDS_ElemIteratorPtr nodeIt = (*fIt)->nodesIterator();
7156 while ( nodeIt->more() ) { // loop on free face nodes
7157 const SMDS_MeshNode* n =
7158 static_cast<const SMDS_MeshNode*>( nodeIt->next() );
7159 gp_XYZ p( n->X(),n->Y(),n->Z() );
7160 dist += ( aBC - p ).SquareModulus();
7162 if ( dist < minDist ) {
7164 freeFaceList.erase( freeFaceList.begin(), fIt++ );
7167 fIt = freeFaceList.erase( fIt++ );
7170 } // choose one of several free faces of a volume
7172 if ( freeFaceList.size() == 1 ) {
7173 const SMDS_MeshElement* aFreeFace = freeFaceList.front();
7174 faceSet->insert( aFreeFace );
7175 // complete a node set with nodes of a found free face
7176 // for ( iNode = 0; iNode < ; iNode++ )
7177 // nodeSet->insert( fNodes[ iNode ] );
7180 } // loop on volumes of a side
7182 // // complete a set of faces if new nodes in a nodeSet appeared
7183 // // ----------------------------------------------------------
7184 // if ( nodeSetSize != nodeSet->size() ) {
7185 // for ( ; nIt != nodeSet->end(); nIt++ ) { // loop on nodes of iSide
7186 // SMDS_ElemIteratorPtr fIt = (*nIt)->GetInverseElementIterator(SMDSAbs_Face);
7187 // while ( fIt->more() ) { // loop on faces sharing a node
7188 // const SMDS_MeshElement* f = fIt->next();
7189 // if ( faceSet->find( f ) == faceSet->end() ) {
7190 // // check if all nodes are in nodeSet and
7191 // // complete setOfFaceNodeSet if they are
7192 // set <const SMDS_MeshNode*> faceNodeSet;
7193 // SMDS_ElemIteratorPtr nodeIt = f->nodesIterator();
7194 // bool allInSet = true;
7195 // while ( nodeIt->more() && allInSet ) { // loop on nodes of a face
7196 // const SMDS_MeshNode* n = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
7197 // if ( nodeSet->find( n ) == nodeSet->end() )
7198 // allInSet = false;
7200 // faceNodeSet.insert( n );
7202 // if ( allInSet ) {
7203 // faceSet->insert( f );
7204 // setOfFaceNodeSet.insert( faceNodeSet );
7210 } // Create temporary faces, if there are volumes given
7213 if ( faceSet1.size() != faceSet2.size() ) {
7214 // delete temporary faces: they are in reverseElements of actual nodes
7215 SMDS_FaceIteratorPtr tmpFaceIt = aTmpFacesMesh.facesIterator();
7216 while ( tmpFaceIt->more() )
7217 aTmpFacesMesh.RemoveElement( tmpFaceIt->next() );
7218 MESSAGE("Diff nb of faces");
7219 return SEW_TOPO_DIFF_SETS_OF_ELEMENTS;
7222 // ============================================================
7223 // 2. Find nodes to merge:
7224 // bind a node to remove to a node to put instead
7225 // ============================================================
7227 TNodeNodeMap nReplaceMap; // bind a node to remove to a node to put instead
7228 if ( theFirstNode1 != theFirstNode2 )
7229 nReplaceMap.insert( TNodeNodeMap::value_type( theFirstNode1, theFirstNode2 ));
7230 if ( theSecondNode1 != theSecondNode2 )
7231 nReplaceMap.insert( TNodeNodeMap::value_type( theSecondNode1, theSecondNode2 ));
7233 LinkID_Gen aLinkID_Gen( GetMeshDS() );
7234 set< long > linkIdSet; // links to process
7235 linkIdSet.insert( aLinkID_Gen.GetLinkID( theFirstNode1, theSecondNode1 ));
7237 typedef pair< const SMDS_MeshNode*, const SMDS_MeshNode* > NLink;
7238 list< NLink > linkList[2];
7239 linkList[0].push_back( NLink( theFirstNode1, theSecondNode1 ));
7240 linkList[1].push_back( NLink( theFirstNode2, theSecondNode2 ));
7241 // loop on links in linkList; find faces by links and append links
7242 // of the found faces to linkList
7243 list< NLink >::iterator linkIt[] = { linkList[0].begin(), linkList[1].begin() } ;
7244 for ( ; linkIt[0] != linkList[0].end(); linkIt[0]++, linkIt[1]++ ) {
7245 NLink link[] = { *linkIt[0], *linkIt[1] };
7246 long linkID = aLinkID_Gen.GetLinkID( link[0].first, link[0].second );
7247 if ( linkIdSet.find( linkID ) == linkIdSet.end() )
7250 // by links, find faces in the face sets,
7251 // and find indices of link nodes in the found faces;
7252 // in a face set, there is only one or no face sharing a link
7253 // ---------------------------------------------------------------
7255 const SMDS_MeshElement* face[] = { 0, 0 };
7256 //const SMDS_MeshNode* faceNodes[ 2 ][ 5 ];
7257 vector<const SMDS_MeshNode*> fnodes1(9);
7258 vector<const SMDS_MeshNode*> fnodes2(9);
7259 //const SMDS_MeshNode* notLinkNodes[ 2 ][ 2 ] = {{ 0, 0 },{ 0, 0 }} ;
7260 vector<const SMDS_MeshNode*> notLinkNodes1(6);
7261 vector<const SMDS_MeshNode*> notLinkNodes2(6);
7262 int iLinkNode[2][2];
7263 for ( iSide = 0; iSide < 2; iSide++ ) { // loop on 2 sides
7264 const SMDS_MeshNode* n1 = link[iSide].first;
7265 const SMDS_MeshNode* n2 = link[iSide].second;
7266 set<const SMDS_MeshElement*> * faceSet = faceSetPtr[ iSide ];
7267 set< const SMDS_MeshElement* > fMap;
7268 for ( int i = 0; i < 2; i++ ) { // loop on 2 nodes of a link
7269 const SMDS_MeshNode* n = i ? n1 : n2; // a node of a link
7270 SMDS_ElemIteratorPtr fIt = n->GetInverseElementIterator(SMDSAbs_Face);
7271 while ( fIt->more() ) { // loop on faces sharing a node
7272 const SMDS_MeshElement* f = fIt->next();
7273 if (faceSet->find( f ) != faceSet->end() && // f is in face set
7274 ! fMap.insert( f ).second ) // f encounters twice
7276 if ( face[ iSide ] ) {
7277 MESSAGE( "2 faces per link " );
7278 aResult = iSide ? SEW_BAD_SIDE2_NODES : SEW_BAD_SIDE1_NODES;
7282 faceSet->erase( f );
7283 // get face nodes and find ones of a link
7288 fnodes1.resize(f->NbNodes()+1);
7289 notLinkNodes1.resize(f->NbNodes()-2);
7292 fnodes2.resize(f->NbNodes()+1);
7293 notLinkNodes2.resize(f->NbNodes()-2);
7296 if(!f->IsQuadratic()) {
7297 SMDS_ElemIteratorPtr nIt = f->nodesIterator();
7298 while ( nIt->more() ) {
7299 const SMDS_MeshNode* n =
7300 static_cast<const SMDS_MeshNode*>( nIt->next() );
7302 iLinkNode[ iSide ][ 0 ] = iNode;
7304 else if ( n == n2 ) {
7305 iLinkNode[ iSide ][ 1 ] = iNode;
7307 //else if ( notLinkNodes[ iSide ][ 0 ] )
7308 // notLinkNodes[ iSide ][ 1 ] = n;
7310 // notLinkNodes[ iSide ][ 0 ] = n;
7314 notLinkNodes1[nbl] = n;
7315 //notLinkNodes1.push_back(n);
7317 notLinkNodes2[nbl] = n;
7318 //notLinkNodes2.push_back(n);
7320 //faceNodes[ iSide ][ iNode++ ] = n;
7322 fnodes1[iNode++] = n;
7325 fnodes2[iNode++] = n;
7329 else { // f->IsQuadratic()
7330 const SMDS_QuadraticFaceOfNodes* F =
7331 static_cast<const SMDS_QuadraticFaceOfNodes*>(f);
7332 // use special nodes iterator
7333 SMDS_NodeIteratorPtr anIter = F->interlacedNodesIterator();
7334 while ( anIter->more() ) {
7335 const SMDS_MeshNode* n =
7336 static_cast<const SMDS_MeshNode*>( anIter->next() );
7338 iLinkNode[ iSide ][ 0 ] = iNode;
7340 else if ( n == n2 ) {
7341 iLinkNode[ iSide ][ 1 ] = iNode;
7346 notLinkNodes1[nbl] = n;
7349 notLinkNodes2[nbl] = n;
7353 fnodes1[iNode++] = n;
7356 fnodes2[iNode++] = n;
7360 //faceNodes[ iSide ][ iNode ] = faceNodes[ iSide ][ 0 ];
7362 fnodes1[iNode] = fnodes1[0];
7365 fnodes2[iNode] = fnodes1[0];
7372 // check similarity of elements of the sides
7373 if (aResult == SEW_OK && ( face[0] && !face[1] ) || ( !face[0] && face[1] )) {
7374 MESSAGE("Correspondent face not found on side " << ( face[0] ? 1 : 0 ));
7375 if ( nReplaceMap.size() == 2 ) { // faces on input nodes not found
7376 aResult = ( face[0] ? SEW_BAD_SIDE2_NODES : SEW_BAD_SIDE1_NODES );
7379 aResult = SEW_TOPO_DIFF_SETS_OF_ELEMENTS;
7381 break; // do not return because it s necessary to remove tmp faces
7384 // set nodes to merge
7385 // -------------------
7387 if ( face[0] && face[1] ) {
7388 int nbNodes = face[0]->NbNodes();
7389 if ( nbNodes != face[1]->NbNodes() ) {
7390 MESSAGE("Diff nb of face nodes");
7391 aResult = SEW_TOPO_DIFF_SETS_OF_ELEMENTS;
7392 break; // do not return because it s necessary to remove tmp faces
7394 bool reverse[] = { false, false }; // order of notLinkNodes of quadrangle
7395 if ( nbNodes == 3 ) {
7396 //nReplaceMap.insert( TNodeNodeMap::value_type
7397 // ( notLinkNodes[0][0], notLinkNodes[1][0] ));
7398 nReplaceMap.insert( TNodeNodeMap::value_type
7399 ( notLinkNodes1[0], notLinkNodes2[0] ));
7402 for ( iSide = 0; iSide < 2; iSide++ ) { // loop on 2 sides
7403 // analyse link orientation in faces
7404 int i1 = iLinkNode[ iSide ][ 0 ];
7405 int i2 = iLinkNode[ iSide ][ 1 ];
7406 reverse[ iSide ] = Abs( i1 - i2 ) == 1 ? i1 > i2 : i2 > i1;
7407 // if notLinkNodes are the first and the last ones, then
7408 // their order does not correspond to the link orientation
7409 if (( i1 == 1 && i2 == 2 ) ||
7410 ( i1 == 2 && i2 == 1 ))
7411 reverse[ iSide ] = !reverse[ iSide ];
7413 if ( reverse[0] == reverse[1] ) {
7414 //nReplaceMap.insert( TNodeNodeMap::value_type
7415 // ( notLinkNodes[0][0], notLinkNodes[1][0] ));
7416 //nReplaceMap.insert( TNodeNodeMap::value_type
7417 // ( notLinkNodes[0][1], notLinkNodes[1][1] ));
7418 for(int nn=0; nn<nbNodes-2; nn++) {
7419 nReplaceMap.insert( TNodeNodeMap::value_type
7420 ( notLinkNodes1[nn], notLinkNodes2[nn] ));
7424 //nReplaceMap.insert( TNodeNodeMap::value_type
7425 // ( notLinkNodes[0][0], notLinkNodes[1][1] ));
7426 //nReplaceMap.insert( TNodeNodeMap::value_type
7427 // ( notLinkNodes[0][1], notLinkNodes[1][0] ));
7428 for(int nn=0; nn<nbNodes-2; nn++) {
7429 nReplaceMap.insert( TNodeNodeMap::value_type
7430 ( notLinkNodes1[nn], notLinkNodes2[nbNodes-3-nn] ));
7435 // add other links of the faces to linkList
7436 // -----------------------------------------
7438 //const SMDS_MeshNode** nodes = faceNodes[ 0 ];
7439 for ( iNode = 0; iNode < nbNodes; iNode++ ) {
7440 //linkID = aLinkID_Gen.GetLinkID( nodes[iNode], nodes[iNode+1] );
7441 linkID = aLinkID_Gen.GetLinkID( fnodes1[iNode], fnodes1[iNode+1] );
7442 pair< set<long>::iterator, bool > iter_isnew = linkIdSet.insert( linkID );
7443 if ( !iter_isnew.second ) { // already in a set: no need to process
7444 linkIdSet.erase( iter_isnew.first );
7446 else // new in set == encountered for the first time: add
7448 //const SMDS_MeshNode* n1 = nodes[ iNode ];
7449 //const SMDS_MeshNode* n2 = nodes[ iNode + 1];
7450 const SMDS_MeshNode* n1 = fnodes1[ iNode ];
7451 const SMDS_MeshNode* n2 = fnodes1[ iNode + 1];
7452 linkList[0].push_back ( NLink( n1, n2 ));
7453 linkList[1].push_back ( NLink( nReplaceMap[n1], nReplaceMap[n2] ));
7457 } // loop on link lists
7459 if ( aResult == SEW_OK &&
7460 ( linkIt[0] != linkList[0].end() ||
7461 !faceSetPtr[0]->empty() || !faceSetPtr[1]->empty() )) {
7462 MESSAGE( (linkIt[0] != linkList[0].end()) <<" "<< (faceSetPtr[0]->empty()) <<
7463 " " << (faceSetPtr[1]->empty()));
7464 aResult = SEW_TOPO_DIFF_SETS_OF_ELEMENTS;
7467 // ====================================================================
7468 // 3. Replace nodes in elements of the side 1 and remove replaced nodes
7469 // ====================================================================
7471 // delete temporary faces: they are in reverseElements of actual nodes
7472 SMDS_FaceIteratorPtr tmpFaceIt = aTmpFacesMesh.facesIterator();
7473 while ( tmpFaceIt->more() )
7474 aTmpFacesMesh.RemoveElement( tmpFaceIt->next() );
7476 if ( aResult != SEW_OK)
7479 list< int > nodeIDsToRemove/*, elemIDsToRemove*/;
7480 // loop on nodes replacement map
7481 TNodeNodeMap::iterator nReplaceMapIt = nReplaceMap.begin(), nnIt;
7482 for ( ; nReplaceMapIt != nReplaceMap.end(); nReplaceMapIt++ )
7483 if ( (*nReplaceMapIt).first != (*nReplaceMapIt).second ) {
7484 const SMDS_MeshNode* nToRemove = (*nReplaceMapIt).first;
7485 nodeIDsToRemove.push_back( nToRemove->GetID() );
7486 // loop on elements sharing nToRemove
7487 SMDS_ElemIteratorPtr invElemIt = nToRemove->GetInverseElementIterator();
7488 while ( invElemIt->more() ) {
7489 const SMDS_MeshElement* e = invElemIt->next();
7490 // get a new suite of nodes: make replacement
7491 int nbReplaced = 0, i = 0, nbNodes = e->NbNodes();
7492 vector< const SMDS_MeshNode*> nodes( nbNodes );
7493 SMDS_ElemIteratorPtr nIt = e->nodesIterator();
7494 while ( nIt->more() ) {
7495 const SMDS_MeshNode* n =
7496 static_cast<const SMDS_MeshNode*>( nIt->next() );
7497 nnIt = nReplaceMap.find( n );
7498 if ( nnIt != nReplaceMap.end() ) {
7504 // if ( nbReplaced == nbNodes && e->GetType() == SMDSAbs_Face )
7505 // elemIDsToRemove.push_back( e->GetID() );
7508 aMesh->ChangeElementNodes( e, & nodes[0], nbNodes );
7512 Remove( nodeIDsToRemove, true );
7517 //================================================================================
7519 * \brief Find corresponding nodes in two sets of faces
7520 * \param theSide1 - first face set
7521 * \param theSide2 - second first face
7522 * \param theFirstNode1 - a boundary node of set 1
7523 * \param theFirstNode2 - a node of set 2 corresponding to theFirstNode1
7524 * \param theSecondNode1 - a boundary node of set 1 linked with theFirstNode1
7525 * \param theSecondNode2 - a node of set 2 corresponding to theSecondNode1
7526 * \param nReplaceMap - output map of corresponding nodes
7527 * \retval bool - is a success or not
7529 //================================================================================
7532 //#define DEBUG_MATCHING_NODES
7535 SMESH_MeshEditor::Sew_Error
7536 SMESH_MeshEditor::FindMatchingNodes(set<const SMDS_MeshElement*>& theSide1,
7537 set<const SMDS_MeshElement*>& theSide2,
7538 const SMDS_MeshNode* theFirstNode1,
7539 const SMDS_MeshNode* theFirstNode2,
7540 const SMDS_MeshNode* theSecondNode1,
7541 const SMDS_MeshNode* theSecondNode2,
7542 TNodeNodeMap & nReplaceMap)
7544 set<const SMDS_MeshElement*> * faceSetPtr[] = { &theSide1, &theSide2 };
7546 nReplaceMap.clear();
7547 if ( theFirstNode1 != theFirstNode2 )
7548 nReplaceMap.insert( make_pair( theFirstNode1, theFirstNode2 ));
7549 if ( theSecondNode1 != theSecondNode2 )
7550 nReplaceMap.insert( make_pair( theSecondNode1, theSecondNode2 ));
7552 set< SMESH_TLink > linkSet; // set of nodes where order of nodes is ignored
7553 linkSet.insert( SMESH_TLink( theFirstNode1, theSecondNode1 ));
7555 list< NLink > linkList[2];
7556 linkList[0].push_back( NLink( theFirstNode1, theSecondNode1 ));
7557 linkList[1].push_back( NLink( theFirstNode2, theSecondNode2 ));
7559 // loop on links in linkList; find faces by links and append links
7560 // of the found faces to linkList
7561 list< NLink >::iterator linkIt[] = { linkList[0].begin(), linkList[1].begin() } ;
7562 for ( ; linkIt[0] != linkList[0].end(); linkIt[0]++, linkIt[1]++ ) {
7563 NLink link[] = { *linkIt[0], *linkIt[1] };
7564 if ( linkSet.find( link[0] ) == linkSet.end() )
7567 // by links, find faces in the face sets,
7568 // and find indices of link nodes in the found faces;
7569 // in a face set, there is only one or no face sharing a link
7570 // ---------------------------------------------------------------
7572 const SMDS_MeshElement* face[] = { 0, 0 };
7573 list<const SMDS_MeshNode*> notLinkNodes[2];
7574 //bool reverse[] = { false, false }; // order of notLinkNodes
7576 for ( int iSide = 0; iSide < 2; iSide++ ) // loop on 2 sides
7578 const SMDS_MeshNode* n1 = link[iSide].first;
7579 const SMDS_MeshNode* n2 = link[iSide].second;
7580 set<const SMDS_MeshElement*> * faceSet = faceSetPtr[ iSide ];
7581 set< const SMDS_MeshElement* > facesOfNode1;
7582 for ( int iNode = 0; iNode < 2; iNode++ ) // loop on 2 nodes of a link
7584 // during a loop of the first node, we find all faces around n1,
7585 // during a loop of the second node, we find one face sharing both n1 and n2
7586 const SMDS_MeshNode* n = iNode ? n1 : n2; // a node of a link
7587 SMDS_ElemIteratorPtr fIt = n->GetInverseElementIterator(SMDSAbs_Face);
7588 while ( fIt->more() ) { // loop on faces sharing a node
7589 const SMDS_MeshElement* f = fIt->next();
7590 if (faceSet->find( f ) != faceSet->end() && // f is in face set
7591 ! facesOfNode1.insert( f ).second ) // f encounters twice
7593 if ( face[ iSide ] ) {
7594 MESSAGE( "2 faces per link " );
7595 return ( iSide ? SEW_BAD_SIDE2_NODES : SEW_BAD_SIDE1_NODES );
7598 faceSet->erase( f );
7600 // get not link nodes
7601 int nbN = f->NbNodes();
7602 if ( f->IsQuadratic() )
7604 nbNodes[ iSide ] = nbN;
7605 list< const SMDS_MeshNode* > & nodes = notLinkNodes[ iSide ];
7606 int i1 = f->GetNodeIndex( n1 );
7607 int i2 = f->GetNodeIndex( n2 );
7608 int iEnd = nbN, iBeg = -1, iDelta = 1;
7609 bool reverse = ( Abs( i1 - i2 ) == 1 ? i1 > i2 : i2 > i1 );
7611 std::swap( iEnd, iBeg ); iDelta = -1;
7616 if ( i == iEnd ) i = iBeg + iDelta;
7617 if ( i == i1 ) break;
7618 nodes.push_back ( f->GetNode( i ) );
7624 // check similarity of elements of the sides
7625 if (( face[0] && !face[1] ) || ( !face[0] && face[1] )) {
7626 MESSAGE("Correspondent face not found on side " << ( face[0] ? 1 : 0 ));
7627 if ( nReplaceMap.size() == 2 ) { // faces on input nodes not found
7628 return ( face[0] ? SEW_BAD_SIDE2_NODES : SEW_BAD_SIDE1_NODES );
7631 return SEW_TOPO_DIFF_SETS_OF_ELEMENTS;
7635 // set nodes to merge
7636 // -------------------
7638 if ( face[0] && face[1] ) {
7639 if ( nbNodes[0] != nbNodes[1] ) {
7640 MESSAGE("Diff nb of face nodes");
7641 return SEW_TOPO_DIFF_SETS_OF_ELEMENTS;
7643 #ifdef DEBUG_MATCHING_NODES
7644 MESSAGE ( " Link 1: " << link[0].first->GetID() <<" "<< link[0].second->GetID()
7645 << " F 1: " << face[0] << "| Link 2: " << link[1].first->GetID() <<" "
7646 << link[1].second->GetID() << " F 2: " << face[1] << " | Bind: " ) ;
7648 int nbN = nbNodes[0];
7650 list<const SMDS_MeshNode*>::iterator n1 = notLinkNodes[0].begin();
7651 list<const SMDS_MeshNode*>::iterator n2 = notLinkNodes[1].begin();
7652 for ( int i = 0 ; i < nbN - 2; ++i ) {
7653 #ifdef DEBUG_MATCHING_NODES
7654 MESSAGE ( (*n1)->GetID() << " to " << (*n2)->GetID() );
7656 nReplaceMap.insert( make_pair( *(n1++), *(n2++) ));
7660 // add other links of the face 1 to linkList
7661 // -----------------------------------------
7663 const SMDS_MeshElement* f0 = face[0];
7664 const SMDS_MeshNode* n1 = f0->GetNode( nbN - 1 );
7665 for ( int i = 0; i < nbN; i++ )
7667 const SMDS_MeshNode* n2 = f0->GetNode( i );
7668 pair< set< SMESH_TLink >::iterator, bool > iter_isnew =
7669 linkSet.insert( SMESH_TLink( n1, n2 ));
7670 if ( !iter_isnew.second ) { // already in a set: no need to process
7671 linkSet.erase( iter_isnew.first );
7673 else // new in set == encountered for the first time: add
7675 #ifdef DEBUG_MATCHING_NODES
7676 MESSAGE ( "Add link 1: " << n1->GetID() << " " << n2->GetID() << " "
7677 << " | link 2: " << nReplaceMap[n1]->GetID() << " " << nReplaceMap[n2]->GetID() << " " );
7679 linkList[0].push_back ( NLink( n1, n2 ));
7680 linkList[1].push_back ( NLink( nReplaceMap[n1], nReplaceMap[n2] ));
7685 } // loop on link lists
7691 \brief Creates a hole in a mesh by doubling the nodes of some particular elements
7692 \param theNodes - identifiers of nodes to be doubled
7693 \param theModifiedElems - identifiers of elements to be updated by the new (doubled)
7694 nodes. If list of element identifiers is empty then nodes are doubled but
7695 they not assigned to elements
7696 \return TRUE if operation has been completed successfully, FALSE otherwise
7698 bool SMESH_MeshEditor::DoubleNodes( const std::list< int >& theListOfNodes,
7699 const std::list< int >& theListOfModifiedElems )
7701 myLastCreatedElems.Clear();
7702 myLastCreatedNodes.Clear();
7704 if ( theListOfNodes.size() == 0 )
7707 SMESHDS_Mesh* aMeshDS = GetMeshDS();
7711 // iterate through nodes and duplicate them
7713 std::map< const SMDS_MeshNode*, const SMDS_MeshNode* > anOldNodeToNewNode;
7715 std::list< int >::const_iterator aNodeIter;
7716 for ( aNodeIter = theListOfNodes.begin(); aNodeIter != theListOfNodes.end(); ++aNodeIter )
7718 int aCurr = *aNodeIter;
7719 SMDS_MeshNode* aNode = (SMDS_MeshNode*)aMeshDS->FindNode( aCurr );
7725 const SMDS_MeshNode* aNewNode = aMeshDS->AddNode( aNode->X(), aNode->Y(), aNode->Z() );
7728 anOldNodeToNewNode[ aNode ] = aNewNode;
7729 myLastCreatedNodes.Append( aNewNode );
7733 // Create map of new nodes for modified elements
7735 std::map< SMDS_MeshElement*, vector<const SMDS_MeshNode*> > anElemToNodes;
7737 std::list< int >::const_iterator anElemIter;
7738 for ( anElemIter = theListOfModifiedElems.begin();
7739 anElemIter != theListOfModifiedElems.end(); ++anElemIter )
7741 int aCurr = *anElemIter;
7742 SMDS_MeshElement* anElem = (SMDS_MeshElement*)aMeshDS->FindElement( aCurr );
7746 vector<const SMDS_MeshNode*> aNodeArr( anElem->NbNodes() );
7748 SMDS_ElemIteratorPtr anIter = anElem->nodesIterator();
7750 while ( anIter->more() )
7752 SMDS_MeshNode* aCurrNode = (SMDS_MeshNode*)anIter->next();
7753 if ( aCurr && anOldNodeToNewNode.find( aCurrNode ) != anOldNodeToNewNode.end() )
7755 const SMDS_MeshNode* aNewNode = anOldNodeToNewNode[ aCurrNode ];
7756 aNodeArr[ ind++ ] = aNewNode;
7759 aNodeArr[ ind++ ] = aCurrNode;
7761 anElemToNodes[ anElem ] = aNodeArr;
7764 // Change nodes of elements
7766 std::map< SMDS_MeshElement*, vector<const SMDS_MeshNode*> >::iterator
7767 anElemToNodesIter = anElemToNodes.begin();
7768 for ( ; anElemToNodesIter != anElemToNodes.end(); ++anElemToNodesIter )
7770 const SMDS_MeshElement* anElem = anElemToNodesIter->first;
7771 vector<const SMDS_MeshNode*> aNodeArr = anElemToNodesIter->second;
7773 aMeshDS->ChangeElementNodes( anElem, &aNodeArr[ 0 ], anElem->NbNodes() );