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 iRepl[ nbRepl++ ] = iCur;
4851 curNodes[ iCur ] = n;
4852 bool isUnique = nodeSet.insert( n ).second;
4854 uniqueNodes[ iUnique++ ] = n;
4858 // Analyse element topology after replacement
4861 int nbUniqueNodes = nodeSet.size();
4862 if ( nbNodes != nbUniqueNodes ) { // some nodes stick
4863 // Polygons and Polyhedral volumes
4864 if (elem->IsPoly()) {
4866 if (elem->GetType() == SMDSAbs_Face) {
4868 vector<const SMDS_MeshNode *> face_nodes (nbNodes);
4870 for (; inode < nbNodes; inode++) {
4871 face_nodes[inode] = curNodes[inode];
4874 vector<const SMDS_MeshNode *> polygons_nodes;
4875 vector<int> quantities;
4876 int nbNew = SimplifyFace(face_nodes, polygons_nodes, quantities);
4880 for (int iface = 0; iface < nbNew - 1; iface++) {
4881 int nbNodes = quantities[iface];
4882 vector<const SMDS_MeshNode *> poly_nodes (nbNodes);
4883 for (int ii = 0; ii < nbNodes; ii++, inode++) {
4884 poly_nodes[ii] = polygons_nodes[inode];
4886 SMDS_MeshElement* newElem = aMesh->AddPolygonalFace(poly_nodes);
4887 myLastCreatedElems.Append(newElem);
4889 aMesh->SetMeshElementOnShape(newElem, aShapeId);
4891 aMesh->ChangeElementNodes(elem, &polygons_nodes[inode], quantities[nbNew - 1]);
4894 rmElemIds.push_back(elem->GetID());
4898 else if (elem->GetType() == SMDSAbs_Volume) {
4899 // Polyhedral volume
4900 if (nbUniqueNodes < 4) {
4901 rmElemIds.push_back(elem->GetID());
4904 // each face has to be analized in order to check volume validity
4905 const SMDS_PolyhedralVolumeOfNodes* aPolyedre =
4906 static_cast<const SMDS_PolyhedralVolumeOfNodes*>( elem );
4908 int nbFaces = aPolyedre->NbFaces();
4910 vector<const SMDS_MeshNode *> poly_nodes;
4911 vector<int> quantities;
4913 for (int iface = 1; iface <= nbFaces; iface++) {
4914 int nbFaceNodes = aPolyedre->NbFaceNodes(iface);
4915 vector<const SMDS_MeshNode *> faceNodes (nbFaceNodes);
4917 for (int inode = 1; inode <= nbFaceNodes; inode++) {
4918 const SMDS_MeshNode * faceNode = aPolyedre->GetFaceNode(iface, inode);
4919 TNodeNodeMap::iterator nnIt = nodeNodeMap.find(faceNode);
4920 if (nnIt != nodeNodeMap.end()) { // faceNode sticks
4921 faceNode = (*nnIt).second;
4923 faceNodes[inode - 1] = faceNode;
4926 SimplifyFace(faceNodes, poly_nodes, quantities);
4929 if (quantities.size() > 3) {
4930 // to be done: remove coincident faces
4933 if (quantities.size() > 3)
4934 aMesh->ChangePolyhedronNodes(elem, poly_nodes, quantities);
4936 rmElemIds.push_back(elem->GetID());
4940 rmElemIds.push_back(elem->GetID());
4951 switch ( nbNodes ) {
4952 case 2: ///////////////////////////////////// EDGE
4953 isOk = false; break;
4954 case 3: ///////////////////////////////////// TRIANGLE
4955 isOk = false; break;
4957 if ( elem->GetType() == SMDSAbs_Volume ) // TETRAHEDRON
4959 else { //////////////////////////////////// QUADRANGLE
4960 if ( nbUniqueNodes < 3 )
4962 else if ( nbRepl == 2 && iRepl[ 1 ] - iRepl[ 0 ] == 2 )
4963 isOk = false; // opposite nodes stick
4966 case 6: ///////////////////////////////////// PENTAHEDRON
4967 if ( nbUniqueNodes == 4 ) {
4968 // ---------------------------------> tetrahedron
4970 iRepl[ 0 ] > 2 && iRepl[ 1 ] > 2 && iRepl[ 2 ] > 2 ) {
4971 // all top nodes stick: reverse a bottom
4972 uniqueNodes[ 0 ] = curNodes [ 1 ];
4973 uniqueNodes[ 1 ] = curNodes [ 0 ];
4975 else if (nbRepl == 3 &&
4976 iRepl[ 0 ] < 3 && iRepl[ 1 ] < 3 && iRepl[ 2 ] < 3 ) {
4977 // all bottom nodes stick: set a top before
4978 uniqueNodes[ 3 ] = uniqueNodes [ 0 ];
4979 uniqueNodes[ 0 ] = curNodes [ 3 ];
4980 uniqueNodes[ 1 ] = curNodes [ 4 ];
4981 uniqueNodes[ 2 ] = curNodes [ 5 ];
4983 else if (nbRepl == 4 &&
4984 iRepl[ 2 ] - iRepl [ 0 ] == 3 && iRepl[ 3 ] - iRepl [ 1 ] == 3 ) {
4985 // a lateral face turns into a line: reverse a bottom
4986 uniqueNodes[ 0 ] = curNodes [ 1 ];
4987 uniqueNodes[ 1 ] = curNodes [ 0 ];
4992 else if ( nbUniqueNodes == 5 ) {
4993 // PENTAHEDRON --------------------> 2 tetrahedrons
4994 if ( nbRepl == 2 && iRepl[ 1 ] - iRepl [ 0 ] == 3 ) {
4995 // a bottom node sticks with a linked top one
4997 SMDS_MeshElement* newElem =
4998 aMesh->AddVolume(curNodes[ 3 ],
5001 curNodes[ iRepl[ 0 ] == 2 ? 1 : 2 ]);
5002 myLastCreatedElems.Append(newElem);
5004 aMesh->SetMeshElementOnShape( newElem, aShapeId );
5005 // 2. : reverse a bottom
5006 uniqueNodes[ 0 ] = curNodes [ 1 ];
5007 uniqueNodes[ 1 ] = curNodes [ 0 ];
5017 if(elem->IsQuadratic()) { // Quadratic quadrangle
5030 if( iRepl[0]==0 && iRepl[1]==1 && iRepl[2]==4 ) {
5031 uniqueNodes[0] = curNodes[0];
5032 uniqueNodes[1] = curNodes[2];
5033 uniqueNodes[2] = curNodes[3];
5034 uniqueNodes[3] = curNodes[5];
5035 uniqueNodes[4] = curNodes[6];
5036 uniqueNodes[5] = curNodes[7];
5039 if( iRepl[0]==0 && iRepl[1]==3 && iRepl[2]==7 ) {
5040 uniqueNodes[0] = curNodes[0];
5041 uniqueNodes[1] = curNodes[1];
5042 uniqueNodes[2] = curNodes[2];
5043 uniqueNodes[3] = curNodes[4];
5044 uniqueNodes[4] = curNodes[5];
5045 uniqueNodes[5] = curNodes[6];
5048 if( iRepl[0]==0 && iRepl[1]==4 && iRepl[2]==7 ) {
5049 uniqueNodes[0] = curNodes[1];
5050 uniqueNodes[1] = curNodes[2];
5051 uniqueNodes[2] = curNodes[3];
5052 uniqueNodes[3] = curNodes[5];
5053 uniqueNodes[4] = curNodes[6];
5054 uniqueNodes[5] = curNodes[0];
5057 if( iRepl[0]==1 && iRepl[1]==2 && iRepl[2]==5 ) {
5058 uniqueNodes[0] = curNodes[0];
5059 uniqueNodes[1] = curNodes[1];
5060 uniqueNodes[2] = curNodes[3];
5061 uniqueNodes[3] = curNodes[4];
5062 uniqueNodes[4] = curNodes[6];
5063 uniqueNodes[5] = curNodes[7];
5066 if( iRepl[0]==1 && iRepl[1]==4 && iRepl[2]==5 ) {
5067 uniqueNodes[0] = curNodes[0];
5068 uniqueNodes[1] = curNodes[2];
5069 uniqueNodes[2] = curNodes[3];
5070 uniqueNodes[3] = curNodes[1];
5071 uniqueNodes[4] = curNodes[6];
5072 uniqueNodes[5] = curNodes[7];
5075 if( iRepl[0]==2 && iRepl[1]==3 && iRepl[2]==6 ) {
5076 uniqueNodes[0] = curNodes[0];
5077 uniqueNodes[1] = curNodes[1];
5078 uniqueNodes[2] = curNodes[2];
5079 uniqueNodes[3] = curNodes[4];
5080 uniqueNodes[4] = curNodes[5];
5081 uniqueNodes[5] = curNodes[7];
5084 if( iRepl[0]==2 && iRepl[1]==5 && iRepl[2]==6 ) {
5085 uniqueNodes[0] = curNodes[0];
5086 uniqueNodes[1] = curNodes[1];
5087 uniqueNodes[2] = curNodes[3];
5088 uniqueNodes[3] = curNodes[4];
5089 uniqueNodes[4] = curNodes[2];
5090 uniqueNodes[5] = curNodes[7];
5093 if( iRepl[0]==3 && iRepl[1]==6 && iRepl[2]==7 ) {
5094 uniqueNodes[0] = curNodes[0];
5095 uniqueNodes[1] = curNodes[1];
5096 uniqueNodes[2] = curNodes[2];
5097 uniqueNodes[3] = curNodes[4];
5098 uniqueNodes[4] = curNodes[5];
5099 uniqueNodes[5] = curNodes[3];
5105 //////////////////////////////////// HEXAHEDRON
5107 SMDS_VolumeTool hexa (elem);
5108 hexa.SetExternalNormal();
5109 if ( nbUniqueNodes == 4 && nbRepl == 6 ) {
5110 //////////////////////// ---> tetrahedron
5111 for ( int iFace = 0; iFace < 6; iFace++ ) {
5112 const int *ind = hexa.GetFaceNodesIndices( iFace ); // indices of face nodes
5113 if (curNodes[ind[ 0 ]] == curNodes[ind[ 1 ]] &&
5114 curNodes[ind[ 0 ]] == curNodes[ind[ 2 ]] &&
5115 curNodes[ind[ 0 ]] == curNodes[ind[ 3 ]] ) {
5116 // one face turns into a point ...
5117 int iOppFace = hexa.GetOppFaceIndex( iFace );
5118 ind = hexa.GetFaceNodesIndices( iOppFace );
5120 iUnique = 2; // reverse a tetrahedron bottom
5121 for ( iCur = 0; iCur < 4 && nbStick < 2; iCur++ ) {
5122 if ( curNodes[ind[ iCur ]] == curNodes[ind[ iCur + 1 ]] )
5124 else if ( iUnique >= 0 )
5125 uniqueNodes[ iUnique-- ] = curNodes[ind[ iCur ]];
5127 if ( nbStick == 1 ) {
5128 // ... and the opposite one - into a triangle.
5130 ind = hexa.GetFaceNodesIndices( iFace );
5131 uniqueNodes[ 3 ] = curNodes[ind[ 0 ]];
5138 else if (nbUniqueNodes == 5 && nbRepl == 4 ) {
5139 //////////////////// HEXAHEDRON ---> 2 tetrahedrons
5140 for ( int iFace = 0; iFace < 6; iFace++ ) {
5141 const int *ind = hexa.GetFaceNodesIndices( iFace ); // indices of face nodes
5142 if (curNodes[ind[ 0 ]] == curNodes[ind[ 1 ]] &&
5143 curNodes[ind[ 0 ]] == curNodes[ind[ 2 ]] &&
5144 curNodes[ind[ 0 ]] == curNodes[ind[ 3 ]] ) {
5145 // one face turns into a point ...
5146 int iOppFace = hexa.GetOppFaceIndex( iFace );
5147 ind = hexa.GetFaceNodesIndices( iOppFace );
5149 iUnique = 2; // reverse a tetrahedron 1 bottom
5150 for ( iCur = 0; iCur < 4 && nbStick == 0; iCur++ ) {
5151 if ( curNodes[ind[ iCur ]] == curNodes[ind[ iCur + 1 ]] )
5153 else if ( iUnique >= 0 )
5154 uniqueNodes[ iUnique-- ] = curNodes[ind[ iCur ]];
5156 if ( nbStick == 0 ) {
5157 // ... and the opposite one is a quadrangle
5159 const int* indTop = hexa.GetFaceNodesIndices( iFace );
5160 uniqueNodes[ 3 ] = curNodes[indTop[ 0 ]];
5163 SMDS_MeshElement* newElem =
5164 aMesh->AddVolume(curNodes[ind[ 0 ]],
5167 curNodes[indTop[ 0 ]]);
5168 myLastCreatedElems.Append(newElem);
5170 aMesh->SetMeshElementOnShape( newElem, aShapeId );
5177 else if ( nbUniqueNodes == 6 && nbRepl == 4 ) {
5178 ////////////////// HEXAHEDRON ---> 2 tetrahedrons or 1 prism
5179 // find indices of quad and tri faces
5180 int iQuadFace[ 6 ], iTriFace[ 6 ], nbQuad = 0, nbTri = 0, iFace;
5181 for ( iFace = 0; iFace < 6; iFace++ ) {
5182 const int *ind = hexa.GetFaceNodesIndices( iFace ); // indices of face nodes
5184 for ( iCur = 0; iCur < 4; iCur++ )
5185 nodeSet.insert( curNodes[ind[ iCur ]] );
5186 nbUniqueNodes = nodeSet.size();
5187 if ( nbUniqueNodes == 3 )
5188 iTriFace[ nbTri++ ] = iFace;
5189 else if ( nbUniqueNodes == 4 )
5190 iQuadFace[ nbQuad++ ] = iFace;
5192 if (nbQuad == 2 && nbTri == 4 &&
5193 hexa.GetOppFaceIndex( iQuadFace[ 0 ] ) == iQuadFace[ 1 ]) {
5194 // 2 opposite quadrangles stuck with a diagonal;
5195 // sample groups of merged indices: (0-4)(2-6)
5196 // --------------------------------------------> 2 tetrahedrons
5197 const int *ind1 = hexa.GetFaceNodesIndices( iQuadFace[ 0 ]); // indices of quad1 nodes
5198 const int *ind2 = hexa.GetFaceNodesIndices( iQuadFace[ 1 ]);
5199 int i0, i1d, i2, i3d, i0t, i2t; // d-daigonal, t-top
5200 if (curNodes[ind1[ 0 ]] == curNodes[ind2[ 0 ]] &&
5201 curNodes[ind1[ 2 ]] == curNodes[ind2[ 2 ]]) {
5202 // stuck with 0-2 diagonal
5210 else if (curNodes[ind1[ 1 ]] == curNodes[ind2[ 3 ]] &&
5211 curNodes[ind1[ 3 ]] == curNodes[ind2[ 1 ]]) {
5212 // stuck with 1-3 diagonal
5224 uniqueNodes[ 0 ] = curNodes [ i0 ];
5225 uniqueNodes[ 1 ] = curNodes [ i1d ];
5226 uniqueNodes[ 2 ] = curNodes [ i3d ];
5227 uniqueNodes[ 3 ] = curNodes [ i0t ];
5230 SMDS_MeshElement* newElem = aMesh->AddVolume(curNodes[ i1d ],
5234 myLastCreatedElems.Append(newElem);
5236 aMesh->SetMeshElementOnShape( newElem, aShapeId );
5239 else if (( nbTri == 2 && nbQuad == 3 ) || // merged (0-4)(1-5)
5240 ( nbTri == 4 && nbQuad == 2 )) { // merged (7-4)(1-5)
5241 // --------------------------------------------> prism
5242 // find 2 opposite triangles
5244 for ( iFace = 0; iFace + 1 < nbTri; iFace++ ) {
5245 if ( hexa.GetOppFaceIndex( iTriFace[ iFace ] ) == iTriFace[ iFace + 1 ]) {
5246 // find indices of kept and replaced nodes
5247 // and fill unique nodes of 2 opposite triangles
5248 const int *ind1 = hexa.GetFaceNodesIndices( iTriFace[ iFace ]);
5249 const int *ind2 = hexa.GetFaceNodesIndices( iTriFace[ iFace + 1 ]);
5250 const SMDS_MeshNode** hexanodes = hexa.GetNodes();
5251 // fill unique nodes
5254 for ( iCur = 0; iCur < 4 && isOk; iCur++ ) {
5255 const SMDS_MeshNode* n = curNodes[ind1[ iCur ]];
5256 const SMDS_MeshNode* nInit = hexanodes[ind1[ iCur ]];
5258 // iCur of a linked node of the opposite face (make normals co-directed):
5259 int iCurOpp = ( iCur == 1 || iCur == 3 ) ? 4 - iCur : iCur;
5260 // check that correspondent corners of triangles are linked
5261 if ( !hexa.IsLinked( ind1[ iCur ], ind2[ iCurOpp ] ))
5264 uniqueNodes[ iUnique ] = n;
5265 uniqueNodes[ iUnique + 3 ] = curNodes[ind2[ iCurOpp ]];
5274 } // if ( nbUniqueNodes == 6 && nbRepl == 4 )
5280 } // switch ( nbNodes )
5282 } // if ( nbNodes != nbUniqueNodes ) // some nodes stick
5285 if (elem->IsPoly() && elem->GetType() == SMDSAbs_Volume) {
5286 // Change nodes of polyedre
5287 const SMDS_PolyhedralVolumeOfNodes* aPolyedre =
5288 static_cast<const SMDS_PolyhedralVolumeOfNodes*>( elem );
5290 int nbFaces = aPolyedre->NbFaces();
5292 vector<const SMDS_MeshNode *> poly_nodes;
5293 vector<int> quantities (nbFaces);
5295 for (int iface = 1; iface <= nbFaces; iface++) {
5296 int inode, nbFaceNodes = aPolyedre->NbFaceNodes(iface);
5297 quantities[iface - 1] = nbFaceNodes;
5299 for (inode = 1; inode <= nbFaceNodes; inode++) {
5300 const SMDS_MeshNode* curNode = aPolyedre->GetFaceNode(iface, inode);
5302 TNodeNodeMap::iterator nnIt = nodeNodeMap.find( curNode );
5303 if (nnIt != nodeNodeMap.end()) { // curNode sticks
5304 curNode = (*nnIt).second;
5306 poly_nodes.push_back(curNode);
5309 aMesh->ChangePolyhedronNodes( elem, poly_nodes, quantities );
5313 // Change regular element or polygon
5314 aMesh->ChangeElementNodes( elem, & uniqueNodes[0], nbUniqueNodes );
5318 // Remove invalid regular element or invalid polygon
5319 rmElemIds.push_back( elem->GetID() );
5322 } // loop on elements
5324 // Remove equal nodes and bad elements
5326 Remove( rmNodeIds, true );
5327 Remove( rmElemIds, false );
5332 // ========================================================
5333 // class : SortableElement
5334 // purpose : allow sorting elements basing on their nodes
5335 // ========================================================
5336 class SortableElement : public set <const SMDS_MeshElement*>
5340 SortableElement( const SMDS_MeshElement* theElem )
5343 SMDS_ElemIteratorPtr nodeIt = theElem->nodesIterator();
5344 while ( nodeIt->more() )
5345 this->insert( nodeIt->next() );
5348 const SMDS_MeshElement* Get() const
5351 void Set(const SMDS_MeshElement* e) const
5356 mutable const SMDS_MeshElement* myElem;
5359 //=======================================================================
5360 //function : FindEqualElements
5361 //purpose : Return list of group of elements built on the same nodes.
5362 // Search among theElements or in the whole mesh if theElements is empty
5363 //=======================================================================
5364 void SMESH_MeshEditor::FindEqualElements(set<const SMDS_MeshElement*> & theElements,
5365 TListOfListOfElementsID & theGroupsOfElementsID)
5367 myLastCreatedElems.Clear();
5368 myLastCreatedNodes.Clear();
5370 typedef set<const SMDS_MeshElement*> TElemsSet;
5371 typedef map< SortableElement, int > TMapOfNodeSet;
5372 typedef list<int> TGroupOfElems;
5375 if ( theElements.empty() )
5376 { // get all elements in the mesh
5377 SMDS_ElemIteratorPtr eIt = GetMeshDS()->elementsIterator();
5378 while ( eIt->more() )
5379 elems.insert( elems.end(), eIt->next());
5382 elems = theElements;
5384 vector< TGroupOfElems > arrayOfGroups;
5385 TGroupOfElems groupOfElems;
5386 TMapOfNodeSet mapOfNodeSet;
5388 TElemsSet::iterator elemIt = elems.begin();
5389 for ( int i = 0, j=0; elemIt != elems.end(); ++elemIt, ++j ) {
5390 const SMDS_MeshElement* curElem = *elemIt;
5391 SortableElement SE(curElem);
5394 pair< TMapOfNodeSet::iterator, bool> pp = mapOfNodeSet.insert(make_pair(SE, i));
5395 if( !(pp.second) ) {
5396 TMapOfNodeSet::iterator& itSE = pp.first;
5397 ind = (*itSE).second;
5398 arrayOfGroups[ind].push_back(curElem->GetID());
5401 groupOfElems.clear();
5402 groupOfElems.push_back(curElem->GetID());
5403 arrayOfGroups.push_back(groupOfElems);
5408 vector< TGroupOfElems >::iterator groupIt = arrayOfGroups.begin();
5409 for ( ; groupIt != arrayOfGroups.end(); ++groupIt ) {
5410 groupOfElems = *groupIt;
5411 if ( groupOfElems.size() > 1 ) {
5412 groupOfElems.sort();
5413 theGroupsOfElementsID.push_back(groupOfElems);
5418 //=======================================================================
5419 //function : MergeElements
5420 //purpose : In each given group, substitute all elements by the first one.
5421 //=======================================================================
5423 void SMESH_MeshEditor::MergeElements(TListOfListOfElementsID & theGroupsOfElementsID)
5425 myLastCreatedElems.Clear();
5426 myLastCreatedNodes.Clear();
5428 typedef list<int> TListOfIDs;
5429 TListOfIDs rmElemIds; // IDs of elems to remove
5431 SMESHDS_Mesh* aMesh = GetMeshDS();
5433 TListOfListOfElementsID::iterator groupsIt = theGroupsOfElementsID.begin();
5434 while ( groupsIt != theGroupsOfElementsID.end() ) {
5435 TListOfIDs& aGroupOfElemID = *groupsIt;
5436 aGroupOfElemID.sort();
5437 int elemIDToKeep = aGroupOfElemID.front();
5438 const SMDS_MeshElement* elemToKeep = aMesh->FindElement(elemIDToKeep);
5439 aGroupOfElemID.pop_front();
5440 TListOfIDs::iterator idIt = aGroupOfElemID.begin();
5441 while ( idIt != aGroupOfElemID.end() ) {
5442 int elemIDToRemove = *idIt;
5443 const SMDS_MeshElement* elemToRemove = aMesh->FindElement(elemIDToRemove);
5444 // add the kept element in groups of removed one (PAL15188)
5445 AddToSameGroups( elemToKeep, elemToRemove, aMesh );
5446 rmElemIds.push_back( elemIDToRemove );
5452 Remove( rmElemIds, false );
5455 //=======================================================================
5456 //function : MergeEqualElements
5457 //purpose : Remove all but one of elements built on the same nodes.
5458 //=======================================================================
5460 void SMESH_MeshEditor::MergeEqualElements()
5462 set<const SMDS_MeshElement*> aMeshElements; /* empty input -
5463 to merge equal elements in the whole mesh */
5464 TListOfListOfElementsID aGroupsOfElementsID;
5465 FindEqualElements(aMeshElements, aGroupsOfElementsID);
5466 MergeElements(aGroupsOfElementsID);
5469 //=======================================================================
5470 //function : FindFaceInSet
5471 //purpose : Return a face having linked nodes n1 and n2 and which is
5472 // - not in avoidSet,
5473 // - in elemSet provided that !elemSet.empty()
5474 //=======================================================================
5476 const SMDS_MeshElement*
5477 SMESH_MeshEditor::FindFaceInSet(const SMDS_MeshNode* n1,
5478 const SMDS_MeshNode* n2,
5479 const TIDSortedElemSet& elemSet,
5480 const TIDSortedElemSet& avoidSet)
5483 SMDS_ElemIteratorPtr invElemIt = n1->GetInverseElementIterator(SMDSAbs_Face);
5484 while ( invElemIt->more() ) { // loop on inverse elements of n1
5485 const SMDS_MeshElement* elem = invElemIt->next();
5486 if (avoidSet.find( elem ) != avoidSet.end() )
5488 if ( !elemSet.empty() && elemSet.find( elem ) == elemSet.end())
5490 // get face nodes and find index of n1
5491 int i1, nbN = elem->NbNodes(), iNode = 0;
5492 //const SMDS_MeshNode* faceNodes[ nbN ], *n;
5493 vector<const SMDS_MeshNode*> faceNodes( nbN );
5494 const SMDS_MeshNode* n;
5495 SMDS_ElemIteratorPtr nIt = elem->nodesIterator();
5496 while ( nIt->more() ) {
5497 faceNodes[ iNode ] = static_cast<const SMDS_MeshNode*>( nIt->next() );
5498 if ( faceNodes[ iNode++ ] == n1 )
5501 // find a n2 linked to n1
5502 if(!elem->IsQuadratic()) {
5503 for ( iNode = 0; iNode < 2; iNode++ ) {
5504 if ( iNode ) // node before n1
5505 n = faceNodes[ i1 == 0 ? nbN - 1 : i1 - 1 ];
5506 else // node after n1
5507 n = faceNodes[ i1 + 1 == nbN ? 0 : i1 + 1 ];
5512 else { // analysis for quadratic elements
5513 bool IsFind = false;
5514 // check using only corner nodes
5515 for ( iNode = 0; iNode < 2; iNode++ ) {
5516 if ( iNode ) // node before n1
5517 n = faceNodes[ i1 == 0 ? nbN/2 - 1 : i1 - 1 ];
5518 else // node after n1
5519 n = faceNodes[ i1 + 1 == nbN/2 ? 0 : i1 + 1 ];
5527 // check using all nodes
5528 const SMDS_QuadraticFaceOfNodes* F =
5529 static_cast<const SMDS_QuadraticFaceOfNodes*>(elem);
5530 // use special nodes iterator
5532 SMDS_NodeIteratorPtr anIter = F->interlacedNodesIterator();
5533 while ( anIter->more() ) {
5534 faceNodes[iNode] = static_cast<const SMDS_MeshNode*>(anIter->next());
5535 if ( faceNodes[ iNode++ ] == n1 )
5538 for ( iNode = 0; iNode < 2; iNode++ ) {
5539 if ( iNode ) // node before n1
5540 n = faceNodes[ i1 == 0 ? nbN - 1 : i1 - 1 ];
5541 else // node after n1
5542 n = faceNodes[ i1 + 1 == nbN ? 0 : i1 + 1 ];
5548 } // end analysis for quadratic elements
5553 //=======================================================================
5554 //function : findAdjacentFace
5556 //=======================================================================
5558 static const SMDS_MeshElement* findAdjacentFace(const SMDS_MeshNode* n1,
5559 const SMDS_MeshNode* n2,
5560 const SMDS_MeshElement* elem)
5562 TIDSortedElemSet elemSet, avoidSet;
5564 avoidSet.insert ( elem );
5565 return SMESH_MeshEditor::FindFaceInSet( n1, n2, elemSet, avoidSet );
5568 //=======================================================================
5569 //function : FindFreeBorder
5571 //=======================================================================
5573 #define ControlFreeBorder SMESH::Controls::FreeEdges::IsFreeEdge
5575 bool SMESH_MeshEditor::FindFreeBorder (const SMDS_MeshNode* theFirstNode,
5576 const SMDS_MeshNode* theSecondNode,
5577 const SMDS_MeshNode* theLastNode,
5578 list< const SMDS_MeshNode* > & theNodes,
5579 list< const SMDS_MeshElement* >& theFaces)
5581 if ( !theFirstNode || !theSecondNode )
5583 // find border face between theFirstNode and theSecondNode
5584 const SMDS_MeshElement* curElem = findAdjacentFace( theFirstNode, theSecondNode, 0 );
5588 theFaces.push_back( curElem );
5589 theNodes.push_back( theFirstNode );
5590 theNodes.push_back( theSecondNode );
5592 //vector<const SMDS_MeshNode*> nodes;
5593 const SMDS_MeshNode *nIgnore = theFirstNode, *nStart = theSecondNode;
5594 set < const SMDS_MeshElement* > foundElems;
5595 bool needTheLast = ( theLastNode != 0 );
5597 while ( nStart != theLastNode ) {
5598 if ( nStart == theFirstNode )
5599 return !needTheLast;
5601 // find all free border faces sharing form nStart
5603 list< const SMDS_MeshElement* > curElemList;
5604 list< const SMDS_MeshNode* > nStartList;
5605 SMDS_ElemIteratorPtr invElemIt = nStart->GetInverseElementIterator(SMDSAbs_Face);
5606 while ( invElemIt->more() ) {
5607 const SMDS_MeshElement* e = invElemIt->next();
5608 if ( e == curElem || foundElems.insert( e ).second ) {
5610 int iNode = 0, nbNodes = e->NbNodes();
5611 //const SMDS_MeshNode* nodes[nbNodes+1];
5612 vector<const SMDS_MeshNode*> nodes(nbNodes+1);
5614 if(e->IsQuadratic()) {
5615 const SMDS_QuadraticFaceOfNodes* F =
5616 static_cast<const SMDS_QuadraticFaceOfNodes*>(e);
5617 // use special nodes iterator
5618 SMDS_NodeIteratorPtr anIter = F->interlacedNodesIterator();
5619 while( anIter->more() ) {
5620 nodes[ iNode++ ] = anIter->next();
5624 SMDS_ElemIteratorPtr nIt = e->nodesIterator();
5625 while ( nIt->more() )
5626 nodes[ iNode++ ] = static_cast<const SMDS_MeshNode*>( nIt->next() );
5628 nodes[ iNode ] = nodes[ 0 ];
5630 for ( iNode = 0; iNode < nbNodes; iNode++ )
5631 if (((nodes[ iNode ] == nStart && nodes[ iNode + 1] != nIgnore ) ||
5632 (nodes[ iNode + 1] == nStart && nodes[ iNode ] != nIgnore )) &&
5633 ControlFreeBorder( &nodes[ iNode ], e->GetID() ))
5635 nStartList.push_back( nodes[ iNode + ( nodes[ iNode ] == nStart ? 1 : 0 )]);
5636 curElemList.push_back( e );
5640 // analyse the found
5642 int nbNewBorders = curElemList.size();
5643 if ( nbNewBorders == 0 ) {
5644 // no free border furthermore
5645 return !needTheLast;
5647 else if ( nbNewBorders == 1 ) {
5648 // one more element found
5650 nStart = nStartList.front();
5651 curElem = curElemList.front();
5652 theFaces.push_back( curElem );
5653 theNodes.push_back( nStart );
5656 // several continuations found
5657 list< const SMDS_MeshElement* >::iterator curElemIt;
5658 list< const SMDS_MeshNode* >::iterator nStartIt;
5659 // check if one of them reached the last node
5660 if ( needTheLast ) {
5661 for (curElemIt = curElemList.begin(), nStartIt = nStartList.begin();
5662 curElemIt!= curElemList.end();
5663 curElemIt++, nStartIt++ )
5664 if ( *nStartIt == theLastNode ) {
5665 theFaces.push_back( *curElemIt );
5666 theNodes.push_back( *nStartIt );
5670 // find the best free border by the continuations
5671 list<const SMDS_MeshNode*> contNodes[ 2 ], *cNL;
5672 list<const SMDS_MeshElement*> contFaces[ 2 ], *cFL;
5673 for (curElemIt = curElemList.begin(), nStartIt = nStartList.begin();
5674 curElemIt!= curElemList.end();
5675 curElemIt++, nStartIt++ )
5677 cNL = & contNodes[ contNodes[0].empty() ? 0 : 1 ];
5678 cFL = & contFaces[ contFaces[0].empty() ? 0 : 1 ];
5679 // find one more free border
5680 if ( ! FindFreeBorder( nStart, *nStartIt, theLastNode, *cNL, *cFL )) {
5684 else if ( !contNodes[0].empty() && !contNodes[1].empty() ) {
5685 // choice: clear a worse one
5686 int iLongest = ( contNodes[0].size() < contNodes[1].size() ? 1 : 0 );
5687 int iWorse = ( needTheLast ? 1 - iLongest : iLongest );
5688 contNodes[ iWorse ].clear();
5689 contFaces[ iWorse ].clear();
5692 if ( contNodes[0].empty() && contNodes[1].empty() )
5695 // append the best free border
5696 cNL = & contNodes[ contNodes[0].empty() ? 1 : 0 ];
5697 cFL = & contFaces[ contFaces[0].empty() ? 1 : 0 ];
5698 theNodes.pop_back(); // remove nIgnore
5699 theNodes.pop_back(); // remove nStart
5700 theFaces.pop_back(); // remove curElem
5701 list< const SMDS_MeshNode* >::iterator nIt = cNL->begin();
5702 list< const SMDS_MeshElement* >::iterator fIt = cFL->begin();
5703 for ( ; nIt != cNL->end(); nIt++ ) theNodes.push_back( *nIt );
5704 for ( ; fIt != cFL->end(); fIt++ ) theFaces.push_back( *fIt );
5707 } // several continuations found
5708 } // while ( nStart != theLastNode )
5713 //=======================================================================
5714 //function : CheckFreeBorderNodes
5715 //purpose : Return true if the tree nodes are on a free border
5716 //=======================================================================
5718 bool SMESH_MeshEditor::CheckFreeBorderNodes(const SMDS_MeshNode* theNode1,
5719 const SMDS_MeshNode* theNode2,
5720 const SMDS_MeshNode* theNode3)
5722 list< const SMDS_MeshNode* > nodes;
5723 list< const SMDS_MeshElement* > faces;
5724 return FindFreeBorder( theNode1, theNode2, theNode3, nodes, faces);
5727 //=======================================================================
5728 //function : SewFreeBorder
5730 //=======================================================================
5732 SMESH_MeshEditor::Sew_Error
5733 SMESH_MeshEditor::SewFreeBorder (const SMDS_MeshNode* theBordFirstNode,
5734 const SMDS_MeshNode* theBordSecondNode,
5735 const SMDS_MeshNode* theBordLastNode,
5736 const SMDS_MeshNode* theSideFirstNode,
5737 const SMDS_MeshNode* theSideSecondNode,
5738 const SMDS_MeshNode* theSideThirdNode,
5739 const bool theSideIsFreeBorder,
5740 const bool toCreatePolygons,
5741 const bool toCreatePolyedrs)
5743 myLastCreatedElems.Clear();
5744 myLastCreatedNodes.Clear();
5746 MESSAGE("::SewFreeBorder()");
5747 Sew_Error aResult = SEW_OK;
5749 // ====================================
5750 // find side nodes and elements
5751 // ====================================
5753 list< const SMDS_MeshNode* > nSide[ 2 ];
5754 list< const SMDS_MeshElement* > eSide[ 2 ];
5755 list< const SMDS_MeshNode* >::iterator nIt[ 2 ];
5756 list< const SMDS_MeshElement* >::iterator eIt[ 2 ];
5760 if (!FindFreeBorder(theBordFirstNode,theBordSecondNode,theBordLastNode,
5761 nSide[0], eSide[0])) {
5762 MESSAGE(" Free Border 1 not found " );
5763 aResult = SEW_BORDER1_NOT_FOUND;
5765 if (theSideIsFreeBorder) {
5768 if (!FindFreeBorder(theSideFirstNode, theSideSecondNode, theSideThirdNode,
5769 nSide[1], eSide[1])) {
5770 MESSAGE(" Free Border 2 not found " );
5771 aResult = ( aResult != SEW_OK ? SEW_BOTH_BORDERS_NOT_FOUND : SEW_BORDER2_NOT_FOUND );
5774 if ( aResult != SEW_OK )
5777 if (!theSideIsFreeBorder) {
5781 // -------------------------------------------------------------------------
5783 // 1. If nodes to merge are not coincident, move nodes of the free border
5784 // from the coord sys defined by the direction from the first to last
5785 // nodes of the border to the correspondent sys of the side 2
5786 // 2. On the side 2, find the links most co-directed with the correspondent
5787 // links of the free border
5788 // -------------------------------------------------------------------------
5790 // 1. Since sewing may brake if there are volumes to split on the side 2,
5791 // we wont move nodes but just compute new coordinates for them
5792 typedef map<const SMDS_MeshNode*, gp_XYZ> TNodeXYZMap;
5793 TNodeXYZMap nBordXYZ;
5794 list< const SMDS_MeshNode* >& bordNodes = nSide[ 0 ];
5795 list< const SMDS_MeshNode* >::iterator nBordIt;
5797 gp_XYZ Pb1( theBordFirstNode->X(), theBordFirstNode->Y(), theBordFirstNode->Z() );
5798 gp_XYZ Pb2( theBordLastNode->X(), theBordLastNode->Y(), theBordLastNode->Z() );
5799 gp_XYZ Ps1( theSideFirstNode->X(), theSideFirstNode->Y(), theSideFirstNode->Z() );
5800 gp_XYZ Ps2( theSideSecondNode->X(), theSideSecondNode->Y(), theSideSecondNode->Z() );
5801 double tol2 = 1.e-8;
5802 gp_Vec Vbs1( Pb1 - Ps1 ),Vbs2( Pb2 - Ps2 );
5803 if ( Vbs1.SquareMagnitude() > tol2 || Vbs2.SquareMagnitude() > tol2 ) {
5804 // Need node movement.
5806 // find X and Z axes to create trsf
5807 gp_Vec Zb( Pb1 - Pb2 ), Zs( Ps1 - Ps2 );
5809 if ( X.SquareMagnitude() <= gp::Resolution() * gp::Resolution() )
5811 X = gp_Ax2( gp::Origin(), Zb ).XDirection();
5814 gp_Ax3 toBordAx( Pb1, Zb, X );
5815 gp_Ax3 fromSideAx( Ps1, Zs, X );
5816 gp_Ax3 toGlobalAx( gp::Origin(), gp::DZ(), gp::DX() );
5818 gp_Trsf toBordSys, fromSide2Sys;
5819 toBordSys.SetTransformation( toBordAx );
5820 fromSide2Sys.SetTransformation( fromSideAx, toGlobalAx );
5821 fromSide2Sys.SetScaleFactor( Zs.Magnitude() / Zb.Magnitude() );
5824 for ( nBordIt = bordNodes.begin(); nBordIt != bordNodes.end(); nBordIt++ ) {
5825 const SMDS_MeshNode* n = *nBordIt;
5826 gp_XYZ xyz( n->X(),n->Y(),n->Z() );
5827 toBordSys.Transforms( xyz );
5828 fromSide2Sys.Transforms( xyz );
5829 nBordXYZ.insert( TNodeXYZMap::value_type( n, xyz ));
5833 // just insert nodes XYZ in the nBordXYZ map
5834 for ( nBordIt = bordNodes.begin(); nBordIt != bordNodes.end(); nBordIt++ ) {
5835 const SMDS_MeshNode* n = *nBordIt;
5836 nBordXYZ.insert( TNodeXYZMap::value_type( n, gp_XYZ( n->X(),n->Y(),n->Z() )));
5840 // 2. On the side 2, find the links most co-directed with the correspondent
5841 // links of the free border
5843 list< const SMDS_MeshElement* >& sideElems = eSide[ 1 ];
5844 list< const SMDS_MeshNode* >& sideNodes = nSide[ 1 ];
5845 sideNodes.push_back( theSideFirstNode );
5847 bool hasVolumes = false;
5848 LinkID_Gen aLinkID_Gen( GetMeshDS() );
5849 set<long> foundSideLinkIDs, checkedLinkIDs;
5850 SMDS_VolumeTool volume;
5851 //const SMDS_MeshNode* faceNodes[ 4 ];
5853 const SMDS_MeshNode* sideNode;
5854 const SMDS_MeshElement* sideElem;
5855 const SMDS_MeshNode* prevSideNode = theSideFirstNode;
5856 const SMDS_MeshNode* prevBordNode = theBordFirstNode;
5857 nBordIt = bordNodes.begin();
5859 // border node position and border link direction to compare with
5860 gp_XYZ bordPos = nBordXYZ[ *nBordIt ];
5861 gp_XYZ bordDir = bordPos - nBordXYZ[ prevBordNode ];
5862 // choose next side node by link direction or by closeness to
5863 // the current border node:
5864 bool searchByDir = ( *nBordIt != theBordLastNode );
5866 // find the next node on the Side 2
5868 double maxDot = -DBL_MAX, minDist = DBL_MAX;
5870 checkedLinkIDs.clear();
5871 gp_XYZ prevXYZ( prevSideNode->X(), prevSideNode->Y(), prevSideNode->Z() );
5873 // loop on inverse elements of current node (prevSideNode) on the Side 2
5874 SMDS_ElemIteratorPtr invElemIt = prevSideNode->GetInverseElementIterator();
5875 while ( invElemIt->more() )
5877 const SMDS_MeshElement* elem = invElemIt->next();
5878 // prepare data for a loop on links coming to prevSideNode, of a face or a volume
5879 int iPrevNode, iNode = 0, nbNodes = elem->NbNodes();
5880 vector< const SMDS_MeshNode* > faceNodes( nbNodes, (const SMDS_MeshNode*)0 );
5881 bool isVolume = volume.Set( elem );
5882 const SMDS_MeshNode** nodes = isVolume ? volume.GetNodes() : & faceNodes[0];
5883 if ( isVolume ) // --volume
5885 else if ( elem->GetType()==SMDSAbs_Face ) { // --face
5886 // retrieve all face nodes and find iPrevNode - an index of the prevSideNode
5887 if(elem->IsQuadratic()) {
5888 const SMDS_QuadraticFaceOfNodes* F =
5889 static_cast<const SMDS_QuadraticFaceOfNodes*>(elem);
5890 // use special nodes iterator
5891 SMDS_NodeIteratorPtr anIter = F->interlacedNodesIterator();
5892 while( anIter->more() ) {
5893 nodes[ iNode ] = anIter->next();
5894 if ( nodes[ iNode++ ] == prevSideNode )
5895 iPrevNode = iNode - 1;
5899 SMDS_ElemIteratorPtr nIt = elem->nodesIterator();
5900 while ( nIt->more() ) {
5901 nodes[ iNode ] = cast2Node( nIt->next() );
5902 if ( nodes[ iNode++ ] == prevSideNode )
5903 iPrevNode = iNode - 1;
5906 // there are 2 links to check
5911 // loop on links, to be precise, on the second node of links
5912 for ( iNode = 0; iNode < nbNodes; iNode++ ) {
5913 const SMDS_MeshNode* n = nodes[ iNode ];
5915 if ( !volume.IsLinked( n, prevSideNode ))
5919 if ( iNode ) // a node before prevSideNode
5920 n = nodes[ iPrevNode == 0 ? elem->NbNodes() - 1 : iPrevNode - 1 ];
5921 else // a node after prevSideNode
5922 n = nodes[ iPrevNode + 1 == elem->NbNodes() ? 0 : iPrevNode + 1 ];
5924 // check if this link was already used
5925 long iLink = aLinkID_Gen.GetLinkID( prevSideNode, n );
5926 bool isJustChecked = !checkedLinkIDs.insert( iLink ).second;
5927 if (!isJustChecked &&
5928 foundSideLinkIDs.find( iLink ) == foundSideLinkIDs.end() )
5930 // test a link geometrically
5931 gp_XYZ nextXYZ ( n->X(), n->Y(), n->Z() );
5932 bool linkIsBetter = false;
5933 double dot = 0.0, dist = 0.0;
5934 if ( searchByDir ) { // choose most co-directed link
5935 dot = bordDir * ( nextXYZ - prevXYZ ).Normalized();
5936 linkIsBetter = ( dot > maxDot );
5938 else { // choose link with the node closest to bordPos
5939 dist = ( nextXYZ - bordPos ).SquareModulus();
5940 linkIsBetter = ( dist < minDist );
5942 if ( linkIsBetter ) {
5951 } // loop on inverse elements of prevSideNode
5954 MESSAGE(" Cant find path by links of the Side 2 ");
5955 return SEW_BAD_SIDE_NODES;
5957 sideNodes.push_back( sideNode );
5958 sideElems.push_back( sideElem );
5959 foundSideLinkIDs.insert ( linkID );
5960 prevSideNode = sideNode;
5962 if ( *nBordIt == theBordLastNode )
5963 searchByDir = false;
5965 // find the next border link to compare with
5966 gp_XYZ sidePos( sideNode->X(), sideNode->Y(), sideNode->Z() );
5967 searchByDir = ( bordDir * ( sidePos - bordPos ) <= 0 );
5968 // move to next border node if sideNode is before forward border node (bordPos)
5969 while ( *nBordIt != theBordLastNode && !searchByDir ) {
5970 prevBordNode = *nBordIt;
5972 bordPos = nBordXYZ[ *nBordIt ];
5973 bordDir = bordPos - nBordXYZ[ prevBordNode ];
5974 searchByDir = ( bordDir * ( sidePos - bordPos ) <= 0 );
5978 while ( sideNode != theSideSecondNode );
5980 if ( hasVolumes && sideNodes.size () != bordNodes.size() && !toCreatePolyedrs) {
5981 MESSAGE("VOLUME SPLITTING IS FORBIDDEN");
5982 return SEW_VOLUMES_TO_SPLIT; // volume splitting is forbidden
5984 } // end nodes search on the side 2
5986 // ============================
5987 // sew the border to the side 2
5988 // ============================
5990 int nbNodes[] = { nSide[0].size(), nSide[1].size() };
5991 int maxNbNodes = Max( nbNodes[0], nbNodes[1] );
5993 TListOfListOfNodes nodeGroupsToMerge;
5994 if ( nbNodes[0] == nbNodes[1] ||
5995 ( theSideIsFreeBorder && !theSideThirdNode)) {
5997 // all nodes are to be merged
5999 for (nIt[0] = nSide[0].begin(), nIt[1] = nSide[1].begin();
6000 nIt[0] != nSide[0].end() && nIt[1] != nSide[1].end();
6001 nIt[0]++, nIt[1]++ )
6003 nodeGroupsToMerge.push_back( list<const SMDS_MeshNode*>() );
6004 nodeGroupsToMerge.back().push_back( *nIt[1] ); // to keep
6005 nodeGroupsToMerge.back().push_back( *nIt[0] ); // to remove
6010 // insert new nodes into the border and the side to get equal nb of segments
6012 // get normalized parameters of nodes on the borders
6013 //double param[ 2 ][ maxNbNodes ];
6015 param[0] = new double [ maxNbNodes ];
6016 param[1] = new double [ maxNbNodes ];
6018 for ( iBord = 0; iBord < 2; iBord++ ) { // loop on 2 borders
6019 list< const SMDS_MeshNode* >& nodes = nSide[ iBord ];
6020 list< const SMDS_MeshNode* >::iterator nIt = nodes.begin();
6021 const SMDS_MeshNode* nPrev = *nIt;
6022 double bordLength = 0;
6023 for ( iNode = 0; nIt != nodes.end(); nIt++, iNode++ ) { // loop on border nodes
6024 const SMDS_MeshNode* nCur = *nIt;
6025 gp_XYZ segment (nCur->X() - nPrev->X(),
6026 nCur->Y() - nPrev->Y(),
6027 nCur->Z() - nPrev->Z());
6028 double segmentLen = segment.Modulus();
6029 bordLength += segmentLen;
6030 param[ iBord ][ iNode ] = bordLength;
6033 // normalize within [0,1]
6034 for ( iNode = 0; iNode < nbNodes[ iBord ]; iNode++ ) {
6035 param[ iBord ][ iNode ] /= bordLength;
6039 // loop on border segments
6040 const SMDS_MeshNode *nPrev[ 2 ] = { 0, 0 };
6041 int i[ 2 ] = { 0, 0 };
6042 nIt[0] = nSide[0].begin(); eIt[0] = eSide[0].begin();
6043 nIt[1] = nSide[1].begin(); eIt[1] = eSide[1].begin();
6045 TElemOfNodeListMap insertMap;
6046 TElemOfNodeListMap::iterator insertMapIt;
6048 // key: elem to insert nodes into
6049 // value: 2 nodes to insert between + nodes to be inserted
6051 bool next[ 2 ] = { false, false };
6053 // find min adjacent segment length after sewing
6054 double nextParam = 10., prevParam = 0;
6055 for ( iBord = 0; iBord < 2; iBord++ ) { // loop on 2 borders
6056 if ( i[ iBord ] + 1 < nbNodes[ iBord ])
6057 nextParam = Min( nextParam, param[iBord][ i[iBord] + 1 ]);
6058 if ( i[ iBord ] > 0 )
6059 prevParam = Max( prevParam, param[iBord][ i[iBord] - 1 ]);
6061 double minParam = Min( param[ 0 ][ i[0] ], param[ 1 ][ i[1] ]);
6062 double maxParam = Max( param[ 0 ][ i[0] ], param[ 1 ][ i[1] ]);
6063 double minSegLen = Min( nextParam - minParam, maxParam - prevParam );
6065 // choose to insert or to merge nodes
6066 double du = param[ 1 ][ i[1] ] - param[ 0 ][ i[0] ];
6067 if ( Abs( du ) <= minSegLen * 0.2 ) {
6070 nodeGroupsToMerge.push_back( list<const SMDS_MeshNode*>() );
6071 const SMDS_MeshNode* n0 = *nIt[0];
6072 const SMDS_MeshNode* n1 = *nIt[1];
6073 nodeGroupsToMerge.back().push_back( n1 );
6074 nodeGroupsToMerge.back().push_back( n0 );
6075 // position of node of the border changes due to merge
6076 param[ 0 ][ i[0] ] += du;
6077 // move n1 for the sake of elem shape evaluation during insertion.
6078 // n1 will be removed by MergeNodes() anyway
6079 const_cast<SMDS_MeshNode*>( n0 )->setXYZ( n1->X(), n1->Y(), n1->Z() );
6080 next[0] = next[1] = true;
6085 int intoBord = ( du < 0 ) ? 0 : 1;
6086 const SMDS_MeshElement* elem = *eIt[ intoBord ];
6087 const SMDS_MeshNode* n1 = nPrev[ intoBord ];
6088 const SMDS_MeshNode* n2 = *nIt[ intoBord ];
6089 const SMDS_MeshNode* nIns = *nIt[ 1 - intoBord ];
6090 if ( intoBord == 1 ) {
6091 // move node of the border to be on a link of elem of the side
6092 gp_XYZ p1 (n1->X(), n1->Y(), n1->Z());
6093 gp_XYZ p2 (n2->X(), n2->Y(), n2->Z());
6094 double ratio = du / ( param[ 1 ][ i[1] ] - param[ 1 ][ i[1]-1 ]);
6095 gp_XYZ p = p2 * ( 1 - ratio ) + p1 * ratio;
6096 GetMeshDS()->MoveNode( nIns, p.X(), p.Y(), p.Z() );
6098 insertMapIt = insertMap.find( elem );
6099 bool notFound = ( insertMapIt == insertMap.end() );
6100 bool otherLink = ( !notFound && (*insertMapIt).second.front() != n1 );
6102 // insert into another link of the same element:
6103 // 1. perform insertion into the other link of the elem
6104 list<const SMDS_MeshNode*> & nodeList = (*insertMapIt).second;
6105 const SMDS_MeshNode* n12 = nodeList.front(); nodeList.pop_front();
6106 const SMDS_MeshNode* n22 = nodeList.front(); nodeList.pop_front();
6107 InsertNodesIntoLink( elem, n12, n22, nodeList, toCreatePolygons );
6108 // 2. perform insertion into the link of adjacent faces
6110 const SMDS_MeshElement* adjElem = findAdjacentFace( n12, n22, elem );
6112 InsertNodesIntoLink( adjElem, n12, n22, nodeList, toCreatePolygons );
6116 if (toCreatePolyedrs) {
6117 // perform insertion into the links of adjacent volumes
6118 UpdateVolumes(n12, n22, nodeList);
6120 // 3. find an element appeared on n1 and n2 after the insertion
6121 insertMap.erase( elem );
6122 elem = findAdjacentFace( n1, n2, 0 );
6124 if ( notFound || otherLink ) {
6125 // add element and nodes of the side into the insertMap
6126 insertMapIt = insertMap.insert
6127 ( TElemOfNodeListMap::value_type( elem, list<const SMDS_MeshNode*>() )).first;
6128 (*insertMapIt).second.push_back( n1 );
6129 (*insertMapIt).second.push_back( n2 );
6131 // add node to be inserted into elem
6132 (*insertMapIt).second.push_back( nIns );
6133 next[ 1 - intoBord ] = true;
6136 // go to the next segment
6137 for ( iBord = 0; iBord < 2; iBord++ ) { // loop on 2 borders
6138 if ( next[ iBord ] ) {
6139 if ( i[ iBord ] != 0 && eIt[ iBord ] != eSide[ iBord ].end())
6141 nPrev[ iBord ] = *nIt[ iBord ];
6142 nIt[ iBord ]++; i[ iBord ]++;
6146 while ( nIt[0] != nSide[0].end() && nIt[1] != nSide[1].end());
6148 // perform insertion of nodes into elements
6150 for (insertMapIt = insertMap.begin();
6151 insertMapIt != insertMap.end();
6154 const SMDS_MeshElement* elem = (*insertMapIt).first;
6155 list<const SMDS_MeshNode*> & nodeList = (*insertMapIt).second;
6156 const SMDS_MeshNode* n1 = nodeList.front(); nodeList.pop_front();
6157 const SMDS_MeshNode* n2 = nodeList.front(); nodeList.pop_front();
6159 InsertNodesIntoLink( elem, n1, n2, nodeList, toCreatePolygons );
6161 if ( !theSideIsFreeBorder ) {
6162 // look for and insert nodes into the faces adjacent to elem
6164 const SMDS_MeshElement* adjElem = findAdjacentFace( n1, n2, elem );
6166 InsertNodesIntoLink( adjElem, n1, n2, nodeList, toCreatePolygons );
6171 if (toCreatePolyedrs) {
6172 // perform insertion into the links of adjacent volumes
6173 UpdateVolumes(n1, n2, nodeList);
6179 } // end: insert new nodes
6181 MergeNodes ( nodeGroupsToMerge );
6186 //=======================================================================
6187 //function : InsertNodesIntoLink
6188 //purpose : insert theNodesToInsert into theFace between theBetweenNode1
6189 // and theBetweenNode2 and split theElement
6190 //=======================================================================
6192 void SMESH_MeshEditor::InsertNodesIntoLink(const SMDS_MeshElement* theFace,
6193 const SMDS_MeshNode* theBetweenNode1,
6194 const SMDS_MeshNode* theBetweenNode2,
6195 list<const SMDS_MeshNode*>& theNodesToInsert,
6196 const bool toCreatePoly)
6198 if ( theFace->GetType() != SMDSAbs_Face ) return;
6200 // find indices of 2 link nodes and of the rest nodes
6201 int iNode = 0, il1, il2, i3, i4;
6202 il1 = il2 = i3 = i4 = -1;
6203 //const SMDS_MeshNode* nodes[ theFace->NbNodes() ];
6204 vector<const SMDS_MeshNode*> nodes( theFace->NbNodes() );
6206 if(theFace->IsQuadratic()) {
6207 const SMDS_QuadraticFaceOfNodes* F =
6208 static_cast<const SMDS_QuadraticFaceOfNodes*>(theFace);
6209 // use special nodes iterator
6210 SMDS_NodeIteratorPtr anIter = F->interlacedNodesIterator();
6211 while( anIter->more() ) {
6212 const SMDS_MeshNode* n = anIter->next();
6213 if ( n == theBetweenNode1 )
6215 else if ( n == theBetweenNode2 )
6221 nodes[ iNode++ ] = n;
6225 SMDS_ElemIteratorPtr nodeIt = theFace->nodesIterator();
6226 while ( nodeIt->more() ) {
6227 const SMDS_MeshNode* n = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
6228 if ( n == theBetweenNode1 )
6230 else if ( n == theBetweenNode2 )
6236 nodes[ iNode++ ] = n;
6239 if ( il1 < 0 || il2 < 0 || i3 < 0 )
6242 // arrange link nodes to go one after another regarding the face orientation
6243 bool reverse = ( Abs( il2 - il1 ) == 1 ? il2 < il1 : il1 < il2 );
6244 list<const SMDS_MeshNode *> aNodesToInsert = theNodesToInsert;
6249 aNodesToInsert.reverse();
6251 // check that not link nodes of a quadrangles are in good order
6252 int nbFaceNodes = theFace->NbNodes();
6253 if ( nbFaceNodes == 4 && i4 - i3 != 1 ) {
6259 if (toCreatePoly || theFace->IsPoly()) {
6262 vector<const SMDS_MeshNode *> poly_nodes (nbFaceNodes + aNodesToInsert.size());
6264 // add nodes of face up to first node of link
6267 if(theFace->IsQuadratic()) {
6268 const SMDS_QuadraticFaceOfNodes* F =
6269 static_cast<const SMDS_QuadraticFaceOfNodes*>(theFace);
6270 // use special nodes iterator
6271 SMDS_NodeIteratorPtr anIter = F->interlacedNodesIterator();
6272 while( anIter->more() && !isFLN ) {
6273 const SMDS_MeshNode* n = anIter->next();
6274 poly_nodes[iNode++] = n;
6275 if (n == nodes[il1]) {
6279 // add nodes to insert
6280 list<const SMDS_MeshNode*>::iterator nIt = aNodesToInsert.begin();
6281 for (; nIt != aNodesToInsert.end(); nIt++) {
6282 poly_nodes[iNode++] = *nIt;
6284 // add nodes of face starting from last node of link
6285 while ( anIter->more() ) {
6286 poly_nodes[iNode++] = anIter->next();
6290 SMDS_ElemIteratorPtr nodeIt = theFace->nodesIterator();
6291 while ( nodeIt->more() && !isFLN ) {
6292 const SMDS_MeshNode* n = static_cast<const SMDS_MeshNode*>( nodeIt->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 ( nodeIt->more() ) {
6305 const SMDS_MeshNode* n = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
6306 poly_nodes[iNode++] = n;
6310 // edit or replace the face
6311 SMESHDS_Mesh *aMesh = GetMeshDS();
6313 if (theFace->IsPoly()) {
6314 aMesh->ChangePolygonNodes(theFace, poly_nodes);
6317 int aShapeId = FindShape( theFace );
6319 SMDS_MeshElement* newElem = aMesh->AddPolygonalFace(poly_nodes);
6320 myLastCreatedElems.Append(newElem);
6321 if ( aShapeId && newElem )
6322 aMesh->SetMeshElementOnShape( newElem, aShapeId );
6324 aMesh->RemoveElement(theFace);
6329 if( !theFace->IsQuadratic() ) {
6331 // put aNodesToInsert between theBetweenNode1 and theBetweenNode2
6332 int nbLinkNodes = 2 + aNodesToInsert.size();
6333 //const SMDS_MeshNode* linkNodes[ nbLinkNodes ];
6334 vector<const SMDS_MeshNode*> linkNodes( nbLinkNodes );
6335 linkNodes[ 0 ] = nodes[ il1 ];
6336 linkNodes[ nbLinkNodes - 1 ] = nodes[ il2 ];
6337 list<const SMDS_MeshNode*>::iterator nIt = aNodesToInsert.begin();
6338 for ( iNode = 1; nIt != aNodesToInsert.end(); nIt++ ) {
6339 linkNodes[ iNode++ ] = *nIt;
6341 // decide how to split a quadrangle: compare possible variants
6342 // and choose which of splits to be a quadrangle
6343 int i1, i2, iSplit, nbSplits = nbLinkNodes - 1, iBestQuad;
6344 if ( nbFaceNodes == 3 ) {
6345 iBestQuad = nbSplits;
6348 else if ( nbFaceNodes == 4 ) {
6349 SMESH::Controls::NumericalFunctorPtr aCrit( new SMESH::Controls::AspectRatio);
6350 double aBestRate = DBL_MAX;
6351 for ( int iQuad = 0; iQuad < nbSplits; iQuad++ ) {
6353 double aBadRate = 0;
6354 // evaluate elements quality
6355 for ( iSplit = 0; iSplit < nbSplits; iSplit++ ) {
6356 if ( iSplit == iQuad ) {
6357 SMDS_FaceOfNodes quad (linkNodes[ i1++ ],
6361 aBadRate += getBadRate( &quad, aCrit );
6364 SMDS_FaceOfNodes tria (linkNodes[ i1++ ],
6366 nodes[ iSplit < iQuad ? i4 : i3 ]);
6367 aBadRate += getBadRate( &tria, aCrit );
6371 if ( aBadRate < aBestRate ) {
6373 aBestRate = aBadRate;
6378 // create new elements
6379 SMESHDS_Mesh *aMesh = GetMeshDS();
6380 int aShapeId = FindShape( theFace );
6383 for ( iSplit = 0; iSplit < nbSplits - 1; iSplit++ ) {
6384 SMDS_MeshElement* newElem = 0;
6385 if ( iSplit == iBestQuad )
6386 newElem = aMesh->AddFace (linkNodes[ i1++ ],
6391 newElem = aMesh->AddFace (linkNodes[ i1++ ],
6393 nodes[ iSplit < iBestQuad ? i4 : i3 ]);
6394 myLastCreatedElems.Append(newElem);
6395 if ( aShapeId && newElem )
6396 aMesh->SetMeshElementOnShape( newElem, aShapeId );
6399 // change nodes of theFace
6400 const SMDS_MeshNode* newNodes[ 4 ];
6401 newNodes[ 0 ] = linkNodes[ i1 ];
6402 newNodes[ 1 ] = linkNodes[ i2 ];
6403 newNodes[ 2 ] = nodes[ iSplit >= iBestQuad ? i3 : i4 ];
6404 newNodes[ 3 ] = nodes[ i4 ];
6405 aMesh->ChangeElementNodes( theFace, newNodes, iSplit == iBestQuad ? 4 : 3 );
6406 } // end if(!theFace->IsQuadratic())
6407 else { // theFace is quadratic
6408 // we have to split theFace on simple triangles and one simple quadrangle
6410 int nbshift = tmp*2;
6411 // shift nodes in nodes[] by nbshift
6413 for(i=0; i<nbshift; i++) {
6414 const SMDS_MeshNode* n = nodes[0];
6415 for(j=0; j<nbFaceNodes-1; j++) {
6416 nodes[j] = nodes[j+1];
6418 nodes[nbFaceNodes-1] = n;
6420 il1 = il1 - nbshift;
6421 // now have to insert nodes between n0 and n1 or n1 and n2 (see below)
6422 // n0 n1 n2 n0 n1 n2
6423 // +-----+-----+ +-----+-----+
6432 // create new elements
6433 SMESHDS_Mesh *aMesh = GetMeshDS();
6434 int aShapeId = FindShape( theFace );
6437 if(nbFaceNodes==6) { // quadratic triangle
6438 SMDS_MeshElement* newElem =
6439 aMesh->AddFace(nodes[3],nodes[4],nodes[5]);
6440 myLastCreatedElems.Append(newElem);
6441 if ( aShapeId && newElem )
6442 aMesh->SetMeshElementOnShape( newElem, aShapeId );
6443 if(theFace->IsMediumNode(nodes[il1])) {
6444 // create quadrangle
6445 newElem = aMesh->AddFace(nodes[0],nodes[1],nodes[3],nodes[5]);
6446 myLastCreatedElems.Append(newElem);
6447 if ( aShapeId && newElem )
6448 aMesh->SetMeshElementOnShape( newElem, aShapeId );
6454 // create quadrangle
6455 newElem = aMesh->AddFace(nodes[1],nodes[2],nodes[3],nodes[5]);
6456 myLastCreatedElems.Append(newElem);
6457 if ( aShapeId && newElem )
6458 aMesh->SetMeshElementOnShape( newElem, aShapeId );
6464 else { // nbFaceNodes==8 - quadratic quadrangle
6465 SMDS_MeshElement* newElem =
6466 aMesh->AddFace(nodes[3],nodes[4],nodes[5]);
6467 myLastCreatedElems.Append(newElem);
6468 if ( aShapeId && newElem )
6469 aMesh->SetMeshElementOnShape( newElem, aShapeId );
6470 newElem = aMesh->AddFace(nodes[5],nodes[6],nodes[7]);
6471 myLastCreatedElems.Append(newElem);
6472 if ( aShapeId && newElem )
6473 aMesh->SetMeshElementOnShape( newElem, aShapeId );
6474 newElem = aMesh->AddFace(nodes[5],nodes[7],nodes[3]);
6475 myLastCreatedElems.Append(newElem);
6476 if ( aShapeId && newElem )
6477 aMesh->SetMeshElementOnShape( newElem, aShapeId );
6478 if(theFace->IsMediumNode(nodes[il1])) {
6479 // create quadrangle
6480 newElem = aMesh->AddFace(nodes[0],nodes[1],nodes[3],nodes[7]);
6481 myLastCreatedElems.Append(newElem);
6482 if ( aShapeId && newElem )
6483 aMesh->SetMeshElementOnShape( newElem, aShapeId );
6489 // create quadrangle
6490 newElem = aMesh->AddFace(nodes[1],nodes[2],nodes[3],nodes[7]);
6491 myLastCreatedElems.Append(newElem);
6492 if ( aShapeId && newElem )
6493 aMesh->SetMeshElementOnShape( newElem, aShapeId );
6499 // create needed triangles using n1,n2,n3 and inserted nodes
6500 int nbn = 2 + aNodesToInsert.size();
6501 //const SMDS_MeshNode* aNodes[nbn];
6502 vector<const SMDS_MeshNode*> aNodes(nbn);
6503 aNodes[0] = nodes[n1];
6504 aNodes[nbn-1] = nodes[n2];
6505 list<const SMDS_MeshNode*>::iterator nIt = aNodesToInsert.begin();
6506 for ( iNode = 1; nIt != aNodesToInsert.end(); nIt++ ) {
6507 aNodes[iNode++] = *nIt;
6509 for(i=1; i<nbn; i++) {
6510 SMDS_MeshElement* newElem =
6511 aMesh->AddFace(aNodes[i-1],aNodes[i],nodes[n3]);
6512 myLastCreatedElems.Append(newElem);
6513 if ( aShapeId && newElem )
6514 aMesh->SetMeshElementOnShape( newElem, aShapeId );
6516 // remove old quadratic face
6517 aMesh->RemoveElement(theFace);
6521 //=======================================================================
6522 //function : UpdateVolumes
6524 //=======================================================================
6525 void SMESH_MeshEditor::UpdateVolumes (const SMDS_MeshNode* theBetweenNode1,
6526 const SMDS_MeshNode* theBetweenNode2,
6527 list<const SMDS_MeshNode*>& theNodesToInsert)
6529 myLastCreatedElems.Clear();
6530 myLastCreatedNodes.Clear();
6532 SMDS_ElemIteratorPtr invElemIt = theBetweenNode1->GetInverseElementIterator(SMDSAbs_Volume);
6533 while (invElemIt->more()) { // loop on inverse elements of theBetweenNode1
6534 const SMDS_MeshElement* elem = invElemIt->next();
6536 // check, if current volume has link theBetweenNode1 - theBetweenNode2
6537 SMDS_VolumeTool aVolume (elem);
6538 if (!aVolume.IsLinked(theBetweenNode1, theBetweenNode2))
6541 // insert new nodes in all faces of the volume, sharing link theBetweenNode1 - theBetweenNode2
6542 int iface, nbFaces = aVolume.NbFaces();
6543 vector<const SMDS_MeshNode *> poly_nodes;
6544 vector<int> quantities (nbFaces);
6546 for (iface = 0; iface < nbFaces; iface++) {
6547 int nbFaceNodes = aVolume.NbFaceNodes(iface), nbInserted = 0;
6548 // faceNodes will contain (nbFaceNodes + 1) nodes, last = first
6549 const SMDS_MeshNode** faceNodes = aVolume.GetFaceNodes(iface);
6551 for (int inode = 0; inode < nbFaceNodes; inode++) {
6552 poly_nodes.push_back(faceNodes[inode]);
6554 if (nbInserted == 0) {
6555 if (faceNodes[inode] == theBetweenNode1) {
6556 if (faceNodes[inode + 1] == theBetweenNode2) {
6557 nbInserted = theNodesToInsert.size();
6559 // add nodes to insert
6560 list<const SMDS_MeshNode*>::iterator nIt = theNodesToInsert.begin();
6561 for (; nIt != theNodesToInsert.end(); nIt++) {
6562 poly_nodes.push_back(*nIt);
6566 else if (faceNodes[inode] == theBetweenNode2) {
6567 if (faceNodes[inode + 1] == theBetweenNode1) {
6568 nbInserted = theNodesToInsert.size();
6570 // add nodes to insert in reversed order
6571 list<const SMDS_MeshNode*>::iterator nIt = theNodesToInsert.end();
6573 for (; nIt != theNodesToInsert.begin(); nIt--) {
6574 poly_nodes.push_back(*nIt);
6576 poly_nodes.push_back(*nIt);
6583 quantities[iface] = nbFaceNodes + nbInserted;
6586 // Replace or update the volume
6587 SMESHDS_Mesh *aMesh = GetMeshDS();
6589 if (elem->IsPoly()) {
6590 aMesh->ChangePolyhedronNodes(elem, poly_nodes, quantities);
6594 int aShapeId = FindShape( elem );
6596 SMDS_MeshElement* newElem =
6597 aMesh->AddPolyhedralVolume(poly_nodes, quantities);
6598 myLastCreatedElems.Append(newElem);
6599 if (aShapeId && newElem)
6600 aMesh->SetMeshElementOnShape(newElem, aShapeId);
6602 aMesh->RemoveElement(elem);
6607 //=======================================================================
6609 * \brief Convert elements contained in a submesh to quadratic
6610 * \retval int - nb of checked elements
6612 //=======================================================================
6614 int SMESH_MeshEditor::convertElemToQuadratic(SMESHDS_SubMesh * theSm,
6615 SMESH_MesherHelper& theHelper,
6616 const bool theForce3d)
6619 if( !theSm ) return nbElem;
6621 const bool notFromGroups = false;
6622 SMDS_ElemIteratorPtr ElemItr = theSm->GetElements();
6623 while(ElemItr->more())
6626 const SMDS_MeshElement* elem = ElemItr->next();
6627 if( !elem || elem->IsQuadratic() ) continue;
6629 int id = elem->GetID();
6630 int nbNodes = elem->NbNodes();
6631 vector<const SMDS_MeshNode *> aNds (nbNodes);
6633 for(int i = 0; i < nbNodes; i++)
6635 aNds[i] = elem->GetNode(i);
6637 SMDSAbs_ElementType aType = elem->GetType();
6639 GetMeshDS()->RemoveFreeElement(elem, theSm, notFromGroups);
6641 const SMDS_MeshElement* NewElem = 0;
6647 NewElem = theHelper.AddEdge(aNds[0], aNds[1], id, theForce3d);
6655 NewElem = theHelper.AddFace(aNds[0], aNds[1], aNds[2], id, theForce3d);
6658 NewElem = theHelper.AddFace(aNds[0], aNds[1], aNds[2], aNds[3], id, theForce3d);
6665 case SMDSAbs_Volume :
6670 NewElem = theHelper.AddVolume(aNds[0], aNds[1], aNds[2], aNds[3], id, theForce3d);
6673 NewElem = theHelper.AddVolume(aNds[0], aNds[1], aNds[2], aNds[3], aNds[4], aNds[5], id, theForce3d);
6676 NewElem = theHelper.AddVolume(aNds[0], aNds[1], aNds[2], aNds[3],
6677 aNds[4], aNds[5], aNds[6], aNds[7], id, theForce3d);
6687 ReplaceElemInGroups( elem, NewElem, GetMeshDS());
6689 theSm->AddElement( NewElem );
6694 //=======================================================================
6695 //function : ConvertToQuadratic
6697 //=======================================================================
6698 void SMESH_MeshEditor::ConvertToQuadratic(const bool theForce3d)
6700 SMESHDS_Mesh* meshDS = GetMeshDS();
6702 SMESH_MesherHelper aHelper(*myMesh);
6703 aHelper.SetIsQuadratic( true );
6704 const bool notFromGroups = false;
6706 int nbCheckedElems = 0;
6707 if ( myMesh->HasShapeToMesh() )
6709 if ( SMESH_subMesh *aSubMesh = myMesh->GetSubMeshContaining(myMesh->GetShapeToMesh()))
6711 SMESH_subMeshIteratorPtr smIt = aSubMesh->getDependsOnIterator(true,false);
6712 while ( smIt->more() ) {
6713 SMESH_subMesh* sm = smIt->next();
6714 if ( SMESHDS_SubMesh *smDS = sm->GetSubMeshDS() ) {
6715 aHelper.SetSubShape( sm->GetSubShape() );
6716 nbCheckedElems += convertElemToQuadratic(smDS, aHelper, theForce3d);
6721 int totalNbElems = meshDS->NbEdges() + meshDS->NbFaces() + meshDS->NbVolumes();
6722 if ( nbCheckedElems < totalNbElems ) // not all elements are in submeshes
6724 SMESHDS_SubMesh *smDS = 0;
6725 SMDS_EdgeIteratorPtr aEdgeItr = meshDS->edgesIterator();
6726 while(aEdgeItr->more())
6728 const SMDS_MeshEdge* edge = aEdgeItr->next();
6729 if(edge && !edge->IsQuadratic())
6731 int id = edge->GetID();
6732 const SMDS_MeshNode* n1 = edge->GetNode(0);
6733 const SMDS_MeshNode* n2 = edge->GetNode(1);
6735 meshDS->RemoveFreeElement(edge, smDS, notFromGroups);
6737 const SMDS_MeshEdge* NewEdge = aHelper.AddEdge(n1, n2, id, theForce3d);
6738 ReplaceElemInGroups( edge, NewEdge, GetMeshDS());
6741 SMDS_FaceIteratorPtr aFaceItr = meshDS->facesIterator();
6742 while(aFaceItr->more())
6744 const SMDS_MeshFace* face = aFaceItr->next();
6745 if(!face || face->IsQuadratic() ) continue;
6747 int id = face->GetID();
6748 int nbNodes = face->NbNodes();
6749 vector<const SMDS_MeshNode *> aNds (nbNodes);
6751 for(int i = 0; i < nbNodes; i++)
6753 aNds[i] = face->GetNode(i);
6756 meshDS->RemoveFreeElement(face, smDS, notFromGroups);
6758 SMDS_MeshFace * NewFace = 0;
6762 NewFace = aHelper.AddFace(aNds[0], aNds[1], aNds[2], id, theForce3d);
6765 NewFace = aHelper.AddFace(aNds[0], aNds[1], aNds[2], aNds[3], id, theForce3d);
6770 ReplaceElemInGroups( face, NewFace, GetMeshDS());
6772 SMDS_VolumeIteratorPtr aVolumeItr = meshDS->volumesIterator();
6773 while(aVolumeItr->more())
6775 const SMDS_MeshVolume* volume = aVolumeItr->next();
6776 if(!volume || volume->IsQuadratic() ) continue;
6778 int id = volume->GetID();
6779 int nbNodes = volume->NbNodes();
6780 vector<const SMDS_MeshNode *> aNds (nbNodes);
6782 for(int i = 0; i < nbNodes; i++)
6784 aNds[i] = volume->GetNode(i);
6787 meshDS->RemoveFreeElement(volume, smDS, notFromGroups);
6789 SMDS_MeshVolume * NewVolume = 0;
6793 NewVolume = aHelper.AddVolume(aNds[0], aNds[1], aNds[2],
6794 aNds[3], id, theForce3d );
6797 NewVolume = aHelper.AddVolume(aNds[0], aNds[1], aNds[2],
6798 aNds[3], aNds[4], aNds[5], id, theForce3d);
6801 NewVolume = aHelper.AddVolume(aNds[0], aNds[1], aNds[2], aNds[3],
6802 aNds[4], aNds[5], aNds[6], aNds[7], id, theForce3d);
6807 ReplaceElemInGroups(volume, NewVolume, meshDS);
6812 //=======================================================================
6814 * \brief Convert quadratic elements to linear ones and remove quadratic nodes
6815 * \retval int - nb of checked elements
6817 //=======================================================================
6819 int SMESH_MeshEditor::removeQuadElem(SMESHDS_SubMesh * theSm,
6820 SMDS_ElemIteratorPtr theItr,
6821 const int theShapeID)
6824 SMESHDS_Mesh* meshDS = GetMeshDS();
6825 const bool notFromGroups = false;
6827 while( theItr->more() )
6829 const SMDS_MeshElement* elem = theItr->next();
6831 if( elem && elem->IsQuadratic())
6833 int id = elem->GetID();
6834 int nbNodes = elem->NbNodes();
6835 vector<const SMDS_MeshNode *> aNds, mediumNodes;
6836 aNds.reserve( nbNodes );
6837 mediumNodes.reserve( nbNodes );
6839 for(int i = 0; i < nbNodes; i++)
6841 const SMDS_MeshNode* n = elem->GetNode(i);
6843 if( elem->IsMediumNode( n ) )
6844 mediumNodes.push_back( n );
6846 aNds.push_back( n );
6848 if( aNds.empty() ) continue;
6849 SMDSAbs_ElementType aType = elem->GetType();
6851 //remove old quadratic element
6852 meshDS->RemoveFreeElement( elem, theSm, notFromGroups );
6854 SMDS_MeshElement * NewElem = AddElement( aNds, aType, false, id );
6855 ReplaceElemInGroups(elem, NewElem, meshDS);
6856 if( theSm && NewElem )
6857 theSm->AddElement( NewElem );
6859 // remove medium nodes
6860 vector<const SMDS_MeshNode*>::iterator nIt = mediumNodes.begin();
6861 for ( ; nIt != mediumNodes.end(); ++nIt ) {
6862 const SMDS_MeshNode* n = *nIt;
6863 if ( n->NbInverseElements() == 0 ) {
6864 if ( n->GetPosition()->GetShapeId() != theShapeID )
6865 meshDS->RemoveFreeNode( n, meshDS->MeshElements
6866 ( n->GetPosition()->GetShapeId() ));
6868 meshDS->RemoveFreeNode( n, theSm );
6876 //=======================================================================
6877 //function : ConvertFromQuadratic
6879 //=======================================================================
6880 bool SMESH_MeshEditor::ConvertFromQuadratic()
6882 int nbCheckedElems = 0;
6883 if ( myMesh->HasShapeToMesh() )
6885 if ( SMESH_subMesh *aSubMesh = myMesh->GetSubMeshContaining(myMesh->GetShapeToMesh()))
6887 SMESH_subMeshIteratorPtr smIt = aSubMesh->getDependsOnIterator(true,false);
6888 while ( smIt->more() ) {
6889 SMESH_subMesh* sm = smIt->next();
6890 if ( SMESHDS_SubMesh *smDS = sm->GetSubMeshDS() )
6891 nbCheckedElems += removeQuadElem( smDS, smDS->GetElements(), sm->GetId() );
6897 GetMeshDS()->NbEdges() + GetMeshDS()->NbFaces() + GetMeshDS()->NbVolumes();
6898 if ( nbCheckedElems < totalNbElems ) // not all elements are in submeshes
6900 SMESHDS_SubMesh *aSM = 0;
6901 removeQuadElem( aSM, GetMeshDS()->elementsIterator(), 0 );
6907 //=======================================================================
6908 //function : SewSideElements
6910 //=======================================================================
6912 SMESH_MeshEditor::Sew_Error
6913 SMESH_MeshEditor::SewSideElements (TIDSortedElemSet& theSide1,
6914 TIDSortedElemSet& theSide2,
6915 const SMDS_MeshNode* theFirstNode1,
6916 const SMDS_MeshNode* theFirstNode2,
6917 const SMDS_MeshNode* theSecondNode1,
6918 const SMDS_MeshNode* theSecondNode2)
6920 myLastCreatedElems.Clear();
6921 myLastCreatedNodes.Clear();
6923 MESSAGE ("::::SewSideElements()");
6924 if ( theSide1.size() != theSide2.size() )
6925 return SEW_DIFF_NB_OF_ELEMENTS;
6927 Sew_Error aResult = SEW_OK;
6929 // 1. Build set of faces representing each side
6930 // 2. Find which nodes of the side 1 to merge with ones on the side 2
6931 // 3. Replace nodes in elements of the side 1 and remove replaced nodes
6933 // =======================================================================
6934 // 1. Build set of faces representing each side:
6935 // =======================================================================
6936 // a. build set of nodes belonging to faces
6937 // b. complete set of faces: find missing fices whose nodes are in set of nodes
6938 // c. create temporary faces representing side of volumes if correspondent
6939 // face does not exist
6941 SMESHDS_Mesh* aMesh = GetMeshDS();
6942 SMDS_Mesh aTmpFacesMesh;
6943 set<const SMDS_MeshElement*> faceSet1, faceSet2;
6944 set<const SMDS_MeshElement*> volSet1, volSet2;
6945 set<const SMDS_MeshNode*> nodeSet1, nodeSet2;
6946 set<const SMDS_MeshElement*> * faceSetPtr[] = { &faceSet1, &faceSet2 };
6947 set<const SMDS_MeshElement*> * volSetPtr[] = { &volSet1, &volSet2 };
6948 set<const SMDS_MeshNode*> * nodeSetPtr[] = { &nodeSet1, &nodeSet2 };
6949 TIDSortedElemSet * elemSetPtr[] = { &theSide1, &theSide2 };
6950 int iSide, iFace, iNode;
6952 for ( iSide = 0; iSide < 2; iSide++ ) {
6953 set<const SMDS_MeshNode*> * nodeSet = nodeSetPtr[ iSide ];
6954 TIDSortedElemSet * elemSet = elemSetPtr[ iSide ];
6955 set<const SMDS_MeshElement*> * faceSet = faceSetPtr[ iSide ];
6956 set<const SMDS_MeshElement*> * volSet = volSetPtr [ iSide ];
6957 set<const SMDS_MeshElement*>::iterator vIt;
6958 TIDSortedElemSet::iterator eIt;
6959 set<const SMDS_MeshNode*>::iterator nIt;
6961 // check that given nodes belong to given elements
6962 const SMDS_MeshNode* n1 = ( iSide == 0 ) ? theFirstNode1 : theFirstNode2;
6963 const SMDS_MeshNode* n2 = ( iSide == 0 ) ? theSecondNode1 : theSecondNode2;
6964 int firstIndex = -1, secondIndex = -1;
6965 for (eIt = elemSet->begin(); eIt != elemSet->end(); eIt++ ) {
6966 const SMDS_MeshElement* elem = *eIt;
6967 if ( firstIndex < 0 ) firstIndex = elem->GetNodeIndex( n1 );
6968 if ( secondIndex < 0 ) secondIndex = elem->GetNodeIndex( n2 );
6969 if ( firstIndex > -1 && secondIndex > -1 ) break;
6971 if ( firstIndex < 0 || secondIndex < 0 ) {
6972 // we can simply return until temporary faces created
6973 return (iSide == 0 ) ? SEW_BAD_SIDE1_NODES : SEW_BAD_SIDE2_NODES;
6976 // -----------------------------------------------------------
6977 // 1a. Collect nodes of existing faces
6978 // and build set of face nodes in order to detect missing
6979 // faces corresponing to sides of volumes
6980 // -----------------------------------------------------------
6982 set< set <const SMDS_MeshNode*> > setOfFaceNodeSet;
6984 // loop on the given element of a side
6985 for (eIt = elemSet->begin(); eIt != elemSet->end(); eIt++ ) {
6986 //const SMDS_MeshElement* elem = *eIt;
6987 const SMDS_MeshElement* elem = *eIt;
6988 if ( elem->GetType() == SMDSAbs_Face ) {
6989 faceSet->insert( elem );
6990 set <const SMDS_MeshNode*> faceNodeSet;
6991 SMDS_ElemIteratorPtr nodeIt = elem->nodesIterator();
6992 while ( nodeIt->more() ) {
6993 const SMDS_MeshNode* n = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
6994 nodeSet->insert( n );
6995 faceNodeSet.insert( n );
6997 setOfFaceNodeSet.insert( faceNodeSet );
6999 else if ( elem->GetType() == SMDSAbs_Volume )
7000 volSet->insert( elem );
7002 // ------------------------------------------------------------------------------
7003 // 1b. Complete set of faces: find missing fices whose nodes are in set of nodes
7004 // ------------------------------------------------------------------------------
7006 for ( nIt = nodeSet->begin(); nIt != nodeSet->end(); nIt++ ) { // loop on nodes of iSide
7007 SMDS_ElemIteratorPtr fIt = (*nIt)->GetInverseElementIterator(SMDSAbs_Face);
7008 while ( fIt->more() ) { // loop on faces sharing a node
7009 const SMDS_MeshElement* f = fIt->next();
7010 if ( faceSet->find( f ) == faceSet->end() ) {
7011 // check if all nodes are in nodeSet and
7012 // complete setOfFaceNodeSet if they are
7013 set <const SMDS_MeshNode*> faceNodeSet;
7014 SMDS_ElemIteratorPtr nodeIt = f->nodesIterator();
7015 bool allInSet = true;
7016 while ( nodeIt->more() && allInSet ) { // loop on nodes of a face
7017 const SMDS_MeshNode* n = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
7018 if ( nodeSet->find( n ) == nodeSet->end() )
7021 faceNodeSet.insert( n );
7024 faceSet->insert( f );
7025 setOfFaceNodeSet.insert( faceNodeSet );
7031 // -------------------------------------------------------------------------
7032 // 1c. Create temporary faces representing sides of volumes if correspondent
7033 // face does not exist
7034 // -------------------------------------------------------------------------
7036 if ( !volSet->empty() ) {
7037 //int nodeSetSize = nodeSet->size();
7039 // loop on given volumes
7040 for ( vIt = volSet->begin(); vIt != volSet->end(); vIt++ ) {
7041 SMDS_VolumeTool vol (*vIt);
7042 // loop on volume faces: find free faces
7043 // --------------------------------------
7044 list<const SMDS_MeshElement* > freeFaceList;
7045 for ( iFace = 0; iFace < vol.NbFaces(); iFace++ ) {
7046 if ( !vol.IsFreeFace( iFace ))
7048 // check if there is already a face with same nodes in a face set
7049 const SMDS_MeshElement* aFreeFace = 0;
7050 const SMDS_MeshNode** fNodes = vol.GetFaceNodes( iFace );
7051 int nbNodes = vol.NbFaceNodes( iFace );
7052 set <const SMDS_MeshNode*> faceNodeSet;
7053 vol.GetFaceNodes( iFace, faceNodeSet );
7054 bool isNewFace = setOfFaceNodeSet.insert( faceNodeSet ).second;
7056 // no such a face is given but it still can exist, check it
7057 if ( nbNodes == 3 ) {
7058 aFreeFace = aMesh->FindFace( fNodes[0],fNodes[1],fNodes[2] );
7060 else if ( nbNodes == 4 ) {
7061 aFreeFace = aMesh->FindFace( fNodes[0],fNodes[1],fNodes[2],fNodes[3] );
7064 vector<const SMDS_MeshNode *> poly_nodes ( fNodes, & fNodes[nbNodes]);
7065 aFreeFace = aMesh->FindFace(poly_nodes);
7069 // create a temporary face
7070 if ( nbNodes == 3 ) {
7071 aFreeFace = aTmpFacesMesh.AddFace( fNodes[0],fNodes[1],fNodes[2] );
7073 else if ( nbNodes == 4 ) {
7074 aFreeFace = aTmpFacesMesh.AddFace( fNodes[0],fNodes[1],fNodes[2],fNodes[3] );
7077 vector<const SMDS_MeshNode *> poly_nodes ( fNodes, & fNodes[nbNodes]);
7078 aFreeFace = aTmpFacesMesh.AddPolygonalFace(poly_nodes);
7082 freeFaceList.push_back( aFreeFace );
7084 } // loop on faces of a volume
7086 // choose one of several free faces
7087 // --------------------------------------
7088 if ( freeFaceList.size() > 1 ) {
7089 // choose a face having max nb of nodes shared by other elems of a side
7090 int maxNbNodes = -1/*, nbExcludedFaces = 0*/;
7091 list<const SMDS_MeshElement* >::iterator fIt = freeFaceList.begin();
7092 while ( fIt != freeFaceList.end() ) { // loop on free faces
7093 int nbSharedNodes = 0;
7094 SMDS_ElemIteratorPtr nodeIt = (*fIt)->nodesIterator();
7095 while ( nodeIt->more() ) { // loop on free face nodes
7096 const SMDS_MeshNode* n =
7097 static_cast<const SMDS_MeshNode*>( nodeIt->next() );
7098 SMDS_ElemIteratorPtr invElemIt = n->GetInverseElementIterator();
7099 while ( invElemIt->more() ) {
7100 const SMDS_MeshElement* e = invElemIt->next();
7101 if ( faceSet->find( e ) != faceSet->end() )
7103 if ( elemSet->find( e ) != elemSet->end() )
7107 if ( nbSharedNodes >= maxNbNodes ) {
7108 maxNbNodes = nbSharedNodes;
7112 freeFaceList.erase( fIt++ ); // here fIt++ occures before erase
7114 if ( freeFaceList.size() > 1 )
7116 // could not choose one face, use another way
7117 // choose a face most close to the bary center of the opposite side
7118 gp_XYZ aBC( 0., 0., 0. );
7119 set <const SMDS_MeshNode*> addedNodes;
7120 TIDSortedElemSet * elemSet2 = elemSetPtr[ 1 - iSide ];
7121 eIt = elemSet2->begin();
7122 for ( eIt = elemSet2->begin(); eIt != elemSet2->end(); eIt++ ) {
7123 SMDS_ElemIteratorPtr nodeIt = (*eIt)->nodesIterator();
7124 while ( nodeIt->more() ) { // loop on free face nodes
7125 const SMDS_MeshNode* n =
7126 static_cast<const SMDS_MeshNode*>( nodeIt->next() );
7127 if ( addedNodes.insert( n ).second )
7128 aBC += gp_XYZ( n->X(),n->Y(),n->Z() );
7131 aBC /= addedNodes.size();
7132 double minDist = DBL_MAX;
7133 fIt = freeFaceList.begin();
7134 while ( fIt != freeFaceList.end() ) { // loop on free faces
7136 SMDS_ElemIteratorPtr nodeIt = (*fIt)->nodesIterator();
7137 while ( nodeIt->more() ) { // loop on free face nodes
7138 const SMDS_MeshNode* n =
7139 static_cast<const SMDS_MeshNode*>( nodeIt->next() );
7140 gp_XYZ p( n->X(),n->Y(),n->Z() );
7141 dist += ( aBC - p ).SquareModulus();
7143 if ( dist < minDist ) {
7145 freeFaceList.erase( freeFaceList.begin(), fIt++ );
7148 fIt = freeFaceList.erase( fIt++ );
7151 } // choose one of several free faces of a volume
7153 if ( freeFaceList.size() == 1 ) {
7154 const SMDS_MeshElement* aFreeFace = freeFaceList.front();
7155 faceSet->insert( aFreeFace );
7156 // complete a node set with nodes of a found free face
7157 // for ( iNode = 0; iNode < ; iNode++ )
7158 // nodeSet->insert( fNodes[ iNode ] );
7161 } // loop on volumes of a side
7163 // // complete a set of faces if new nodes in a nodeSet appeared
7164 // // ----------------------------------------------------------
7165 // if ( nodeSetSize != nodeSet->size() ) {
7166 // for ( ; nIt != nodeSet->end(); nIt++ ) { // loop on nodes of iSide
7167 // SMDS_ElemIteratorPtr fIt = (*nIt)->GetInverseElementIterator(SMDSAbs_Face);
7168 // while ( fIt->more() ) { // loop on faces sharing a node
7169 // const SMDS_MeshElement* f = fIt->next();
7170 // if ( faceSet->find( f ) == faceSet->end() ) {
7171 // // check if all nodes are in nodeSet and
7172 // // complete setOfFaceNodeSet if they are
7173 // set <const SMDS_MeshNode*> faceNodeSet;
7174 // SMDS_ElemIteratorPtr nodeIt = f->nodesIterator();
7175 // bool allInSet = true;
7176 // while ( nodeIt->more() && allInSet ) { // loop on nodes of a face
7177 // const SMDS_MeshNode* n = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
7178 // if ( nodeSet->find( n ) == nodeSet->end() )
7179 // allInSet = false;
7181 // faceNodeSet.insert( n );
7183 // if ( allInSet ) {
7184 // faceSet->insert( f );
7185 // setOfFaceNodeSet.insert( faceNodeSet );
7191 } // Create temporary faces, if there are volumes given
7194 if ( faceSet1.size() != faceSet2.size() ) {
7195 // delete temporary faces: they are in reverseElements of actual nodes
7196 SMDS_FaceIteratorPtr tmpFaceIt = aTmpFacesMesh.facesIterator();
7197 while ( tmpFaceIt->more() )
7198 aTmpFacesMesh.RemoveElement( tmpFaceIt->next() );
7199 MESSAGE("Diff nb of faces");
7200 return SEW_TOPO_DIFF_SETS_OF_ELEMENTS;
7203 // ============================================================
7204 // 2. Find nodes to merge:
7205 // bind a node to remove to a node to put instead
7206 // ============================================================
7208 TNodeNodeMap nReplaceMap; // bind a node to remove to a node to put instead
7209 if ( theFirstNode1 != theFirstNode2 )
7210 nReplaceMap.insert( TNodeNodeMap::value_type( theFirstNode1, theFirstNode2 ));
7211 if ( theSecondNode1 != theSecondNode2 )
7212 nReplaceMap.insert( TNodeNodeMap::value_type( theSecondNode1, theSecondNode2 ));
7214 LinkID_Gen aLinkID_Gen( GetMeshDS() );
7215 set< long > linkIdSet; // links to process
7216 linkIdSet.insert( aLinkID_Gen.GetLinkID( theFirstNode1, theSecondNode1 ));
7218 typedef pair< const SMDS_MeshNode*, const SMDS_MeshNode* > NLink;
7219 list< NLink > linkList[2];
7220 linkList[0].push_back( NLink( theFirstNode1, theSecondNode1 ));
7221 linkList[1].push_back( NLink( theFirstNode2, theSecondNode2 ));
7222 // loop on links in linkList; find faces by links and append links
7223 // of the found faces to linkList
7224 list< NLink >::iterator linkIt[] = { linkList[0].begin(), linkList[1].begin() } ;
7225 for ( ; linkIt[0] != linkList[0].end(); linkIt[0]++, linkIt[1]++ ) {
7226 NLink link[] = { *linkIt[0], *linkIt[1] };
7227 long linkID = aLinkID_Gen.GetLinkID( link[0].first, link[0].second );
7228 if ( linkIdSet.find( linkID ) == linkIdSet.end() )
7231 // by links, find faces in the face sets,
7232 // and find indices of link nodes in the found faces;
7233 // in a face set, there is only one or no face sharing a link
7234 // ---------------------------------------------------------------
7236 const SMDS_MeshElement* face[] = { 0, 0 };
7237 //const SMDS_MeshNode* faceNodes[ 2 ][ 5 ];
7238 vector<const SMDS_MeshNode*> fnodes1(9);
7239 vector<const SMDS_MeshNode*> fnodes2(9);
7240 //const SMDS_MeshNode* notLinkNodes[ 2 ][ 2 ] = {{ 0, 0 },{ 0, 0 }} ;
7241 vector<const SMDS_MeshNode*> notLinkNodes1(6);
7242 vector<const SMDS_MeshNode*> notLinkNodes2(6);
7243 int iLinkNode[2][2];
7244 for ( iSide = 0; iSide < 2; iSide++ ) { // loop on 2 sides
7245 const SMDS_MeshNode* n1 = link[iSide].first;
7246 const SMDS_MeshNode* n2 = link[iSide].second;
7247 set<const SMDS_MeshElement*> * faceSet = faceSetPtr[ iSide ];
7248 set< const SMDS_MeshElement* > fMap;
7249 for ( int i = 0; i < 2; i++ ) { // loop on 2 nodes of a link
7250 const SMDS_MeshNode* n = i ? n1 : n2; // a node of a link
7251 SMDS_ElemIteratorPtr fIt = n->GetInverseElementIterator(SMDSAbs_Face);
7252 while ( fIt->more() ) { // loop on faces sharing a node
7253 const SMDS_MeshElement* f = fIt->next();
7254 if (faceSet->find( f ) != faceSet->end() && // f is in face set
7255 ! fMap.insert( f ).second ) // f encounters twice
7257 if ( face[ iSide ] ) {
7258 MESSAGE( "2 faces per link " );
7259 aResult = iSide ? SEW_BAD_SIDE2_NODES : SEW_BAD_SIDE1_NODES;
7263 faceSet->erase( f );
7264 // get face nodes and find ones of a link
7269 fnodes1.resize(f->NbNodes()+1);
7270 notLinkNodes1.resize(f->NbNodes()-2);
7273 fnodes2.resize(f->NbNodes()+1);
7274 notLinkNodes2.resize(f->NbNodes()-2);
7277 if(!f->IsQuadratic()) {
7278 SMDS_ElemIteratorPtr nIt = f->nodesIterator();
7279 while ( nIt->more() ) {
7280 const SMDS_MeshNode* n =
7281 static_cast<const SMDS_MeshNode*>( nIt->next() );
7283 iLinkNode[ iSide ][ 0 ] = iNode;
7285 else if ( n == n2 ) {
7286 iLinkNode[ iSide ][ 1 ] = iNode;
7288 //else if ( notLinkNodes[ iSide ][ 0 ] )
7289 // notLinkNodes[ iSide ][ 1 ] = n;
7291 // notLinkNodes[ iSide ][ 0 ] = n;
7295 notLinkNodes1[nbl] = n;
7296 //notLinkNodes1.push_back(n);
7298 notLinkNodes2[nbl] = n;
7299 //notLinkNodes2.push_back(n);
7301 //faceNodes[ iSide ][ iNode++ ] = n;
7303 fnodes1[iNode++] = n;
7306 fnodes2[iNode++] = n;
7310 else { // f->IsQuadratic()
7311 const SMDS_QuadraticFaceOfNodes* F =
7312 static_cast<const SMDS_QuadraticFaceOfNodes*>(f);
7313 // use special nodes iterator
7314 SMDS_NodeIteratorPtr anIter = F->interlacedNodesIterator();
7315 while ( anIter->more() ) {
7316 const SMDS_MeshNode* n =
7317 static_cast<const SMDS_MeshNode*>( anIter->next() );
7319 iLinkNode[ iSide ][ 0 ] = iNode;
7321 else if ( n == n2 ) {
7322 iLinkNode[ iSide ][ 1 ] = iNode;
7327 notLinkNodes1[nbl] = n;
7330 notLinkNodes2[nbl] = n;
7334 fnodes1[iNode++] = n;
7337 fnodes2[iNode++] = n;
7341 //faceNodes[ iSide ][ iNode ] = faceNodes[ iSide ][ 0 ];
7343 fnodes1[iNode] = fnodes1[0];
7346 fnodes2[iNode] = fnodes1[0];
7353 // check similarity of elements of the sides
7354 if (aResult == SEW_OK && ( face[0] && !face[1] ) || ( !face[0] && face[1] )) {
7355 MESSAGE("Correspondent face not found on side " << ( face[0] ? 1 : 0 ));
7356 if ( nReplaceMap.size() == 2 ) { // faces on input nodes not found
7357 aResult = ( face[0] ? SEW_BAD_SIDE2_NODES : SEW_BAD_SIDE1_NODES );
7360 aResult = SEW_TOPO_DIFF_SETS_OF_ELEMENTS;
7362 break; // do not return because it s necessary to remove tmp faces
7365 // set nodes to merge
7366 // -------------------
7368 if ( face[0] && face[1] ) {
7369 int nbNodes = face[0]->NbNodes();
7370 if ( nbNodes != face[1]->NbNodes() ) {
7371 MESSAGE("Diff nb of face nodes");
7372 aResult = SEW_TOPO_DIFF_SETS_OF_ELEMENTS;
7373 break; // do not return because it s necessary to remove tmp faces
7375 bool reverse[] = { false, false }; // order of notLinkNodes of quadrangle
7376 if ( nbNodes == 3 ) {
7377 //nReplaceMap.insert( TNodeNodeMap::value_type
7378 // ( notLinkNodes[0][0], notLinkNodes[1][0] ));
7379 nReplaceMap.insert( TNodeNodeMap::value_type
7380 ( notLinkNodes1[0], notLinkNodes2[0] ));
7383 for ( iSide = 0; iSide < 2; iSide++ ) { // loop on 2 sides
7384 // analyse link orientation in faces
7385 int i1 = iLinkNode[ iSide ][ 0 ];
7386 int i2 = iLinkNode[ iSide ][ 1 ];
7387 reverse[ iSide ] = Abs( i1 - i2 ) == 1 ? i1 > i2 : i2 > i1;
7388 // if notLinkNodes are the first and the last ones, then
7389 // their order does not correspond to the link orientation
7390 if (( i1 == 1 && i2 == 2 ) ||
7391 ( i1 == 2 && i2 == 1 ))
7392 reverse[ iSide ] = !reverse[ iSide ];
7394 if ( reverse[0] == reverse[1] ) {
7395 //nReplaceMap.insert( TNodeNodeMap::value_type
7396 // ( notLinkNodes[0][0], notLinkNodes[1][0] ));
7397 //nReplaceMap.insert( TNodeNodeMap::value_type
7398 // ( notLinkNodes[0][1], notLinkNodes[1][1] ));
7399 for(int nn=0; nn<nbNodes-2; nn++) {
7400 nReplaceMap.insert( TNodeNodeMap::value_type
7401 ( notLinkNodes1[nn], notLinkNodes2[nn] ));
7405 //nReplaceMap.insert( TNodeNodeMap::value_type
7406 // ( notLinkNodes[0][0], notLinkNodes[1][1] ));
7407 //nReplaceMap.insert( TNodeNodeMap::value_type
7408 // ( notLinkNodes[0][1], notLinkNodes[1][0] ));
7409 for(int nn=0; nn<nbNodes-2; nn++) {
7410 nReplaceMap.insert( TNodeNodeMap::value_type
7411 ( notLinkNodes1[nn], notLinkNodes2[nbNodes-3-nn] ));
7416 // add other links of the faces to linkList
7417 // -----------------------------------------
7419 //const SMDS_MeshNode** nodes = faceNodes[ 0 ];
7420 for ( iNode = 0; iNode < nbNodes; iNode++ ) {
7421 //linkID = aLinkID_Gen.GetLinkID( nodes[iNode], nodes[iNode+1] );
7422 linkID = aLinkID_Gen.GetLinkID( fnodes1[iNode], fnodes1[iNode+1] );
7423 pair< set<long>::iterator, bool > iter_isnew = linkIdSet.insert( linkID );
7424 if ( !iter_isnew.second ) { // already in a set: no need to process
7425 linkIdSet.erase( iter_isnew.first );
7427 else // new in set == encountered for the first time: add
7429 //const SMDS_MeshNode* n1 = nodes[ iNode ];
7430 //const SMDS_MeshNode* n2 = nodes[ iNode + 1];
7431 const SMDS_MeshNode* n1 = fnodes1[ iNode ];
7432 const SMDS_MeshNode* n2 = fnodes1[ iNode + 1];
7433 linkList[0].push_back ( NLink( n1, n2 ));
7434 linkList[1].push_back ( NLink( nReplaceMap[n1], nReplaceMap[n2] ));
7438 } // loop on link lists
7440 if ( aResult == SEW_OK &&
7441 ( linkIt[0] != linkList[0].end() ||
7442 !faceSetPtr[0]->empty() || !faceSetPtr[1]->empty() )) {
7443 MESSAGE( (linkIt[0] != linkList[0].end()) <<" "<< (faceSetPtr[0]->empty()) <<
7444 " " << (faceSetPtr[1]->empty()));
7445 aResult = SEW_TOPO_DIFF_SETS_OF_ELEMENTS;
7448 // ====================================================================
7449 // 3. Replace nodes in elements of the side 1 and remove replaced nodes
7450 // ====================================================================
7452 // delete temporary faces: they are in reverseElements of actual nodes
7453 SMDS_FaceIteratorPtr tmpFaceIt = aTmpFacesMesh.facesIterator();
7454 while ( tmpFaceIt->more() )
7455 aTmpFacesMesh.RemoveElement( tmpFaceIt->next() );
7457 if ( aResult != SEW_OK)
7460 list< int > nodeIDsToRemove/*, elemIDsToRemove*/;
7461 // loop on nodes replacement map
7462 TNodeNodeMap::iterator nReplaceMapIt = nReplaceMap.begin(), nnIt;
7463 for ( ; nReplaceMapIt != nReplaceMap.end(); nReplaceMapIt++ )
7464 if ( (*nReplaceMapIt).first != (*nReplaceMapIt).second ) {
7465 const SMDS_MeshNode* nToRemove = (*nReplaceMapIt).first;
7466 nodeIDsToRemove.push_back( nToRemove->GetID() );
7467 // loop on elements sharing nToRemove
7468 SMDS_ElemIteratorPtr invElemIt = nToRemove->GetInverseElementIterator();
7469 while ( invElemIt->more() ) {
7470 const SMDS_MeshElement* e = invElemIt->next();
7471 // get a new suite of nodes: make replacement
7472 int nbReplaced = 0, i = 0, nbNodes = e->NbNodes();
7473 vector< const SMDS_MeshNode*> nodes( nbNodes );
7474 SMDS_ElemIteratorPtr nIt = e->nodesIterator();
7475 while ( nIt->more() ) {
7476 const SMDS_MeshNode* n =
7477 static_cast<const SMDS_MeshNode*>( nIt->next() );
7478 nnIt = nReplaceMap.find( n );
7479 if ( nnIt != nReplaceMap.end() ) {
7485 // if ( nbReplaced == nbNodes && e->GetType() == SMDSAbs_Face )
7486 // elemIDsToRemove.push_back( e->GetID() );
7489 aMesh->ChangeElementNodes( e, & nodes[0], nbNodes );
7493 Remove( nodeIDsToRemove, true );
7498 //================================================================================
7500 * \brief Find corresponding nodes in two sets of faces
7501 * \param theSide1 - first face set
7502 * \param theSide2 - second first face
7503 * \param theFirstNode1 - a boundary node of set 1
7504 * \param theFirstNode2 - a node of set 2 corresponding to theFirstNode1
7505 * \param theSecondNode1 - a boundary node of set 1 linked with theFirstNode1
7506 * \param theSecondNode2 - a node of set 2 corresponding to theSecondNode1
7507 * \param nReplaceMap - output map of corresponding nodes
7508 * \retval bool - is a success or not
7510 //================================================================================
7513 //#define DEBUG_MATCHING_NODES
7516 SMESH_MeshEditor::Sew_Error
7517 SMESH_MeshEditor::FindMatchingNodes(set<const SMDS_MeshElement*>& theSide1,
7518 set<const SMDS_MeshElement*>& theSide2,
7519 const SMDS_MeshNode* theFirstNode1,
7520 const SMDS_MeshNode* theFirstNode2,
7521 const SMDS_MeshNode* theSecondNode1,
7522 const SMDS_MeshNode* theSecondNode2,
7523 TNodeNodeMap & nReplaceMap)
7525 set<const SMDS_MeshElement*> * faceSetPtr[] = { &theSide1, &theSide2 };
7527 nReplaceMap.clear();
7528 if ( theFirstNode1 != theFirstNode2 )
7529 nReplaceMap.insert( make_pair( theFirstNode1, theFirstNode2 ));
7530 if ( theSecondNode1 != theSecondNode2 )
7531 nReplaceMap.insert( make_pair( theSecondNode1, theSecondNode2 ));
7533 set< SMESH_TLink > linkSet; // set of nodes where order of nodes is ignored
7534 linkSet.insert( SMESH_TLink( theFirstNode1, theSecondNode1 ));
7536 list< NLink > linkList[2];
7537 linkList[0].push_back( NLink( theFirstNode1, theSecondNode1 ));
7538 linkList[1].push_back( NLink( theFirstNode2, theSecondNode2 ));
7540 // loop on links in linkList; find faces by links and append links
7541 // of the found faces to linkList
7542 list< NLink >::iterator linkIt[] = { linkList[0].begin(), linkList[1].begin() } ;
7543 for ( ; linkIt[0] != linkList[0].end(); linkIt[0]++, linkIt[1]++ ) {
7544 NLink link[] = { *linkIt[0], *linkIt[1] };
7545 if ( linkSet.find( link[0] ) == linkSet.end() )
7548 // by links, find faces in the face sets,
7549 // and find indices of link nodes in the found faces;
7550 // in a face set, there is only one or no face sharing a link
7551 // ---------------------------------------------------------------
7553 const SMDS_MeshElement* face[] = { 0, 0 };
7554 list<const SMDS_MeshNode*> notLinkNodes[2];
7555 //bool reverse[] = { false, false }; // order of notLinkNodes
7557 for ( int iSide = 0; iSide < 2; iSide++ ) // loop on 2 sides
7559 const SMDS_MeshNode* n1 = link[iSide].first;
7560 const SMDS_MeshNode* n2 = link[iSide].second;
7561 set<const SMDS_MeshElement*> * faceSet = faceSetPtr[ iSide ];
7562 set< const SMDS_MeshElement* > facesOfNode1;
7563 for ( int iNode = 0; iNode < 2; iNode++ ) // loop on 2 nodes of a link
7565 // during a loop of the first node, we find all faces around n1,
7566 // during a loop of the second node, we find one face sharing both n1 and n2
7567 const SMDS_MeshNode* n = iNode ? n1 : n2; // a node of a link
7568 SMDS_ElemIteratorPtr fIt = n->GetInverseElementIterator(SMDSAbs_Face);
7569 while ( fIt->more() ) { // loop on faces sharing a node
7570 const SMDS_MeshElement* f = fIt->next();
7571 if (faceSet->find( f ) != faceSet->end() && // f is in face set
7572 ! facesOfNode1.insert( f ).second ) // f encounters twice
7574 if ( face[ iSide ] ) {
7575 MESSAGE( "2 faces per link " );
7576 return ( iSide ? SEW_BAD_SIDE2_NODES : SEW_BAD_SIDE1_NODES );
7579 faceSet->erase( f );
7581 // get not link nodes
7582 int nbN = f->NbNodes();
7583 if ( f->IsQuadratic() )
7585 nbNodes[ iSide ] = nbN;
7586 list< const SMDS_MeshNode* > & nodes = notLinkNodes[ iSide ];
7587 int i1 = f->GetNodeIndex( n1 );
7588 int i2 = f->GetNodeIndex( n2 );
7589 int iEnd = nbN, iBeg = -1, iDelta = 1;
7590 bool reverse = ( Abs( i1 - i2 ) == 1 ? i1 > i2 : i2 > i1 );
7592 std::swap( iEnd, iBeg ); iDelta = -1;
7597 if ( i == iEnd ) i = iBeg + iDelta;
7598 if ( i == i1 ) break;
7599 nodes.push_back ( f->GetNode( i ) );
7605 // check similarity of elements of the sides
7606 if (( face[0] && !face[1] ) || ( !face[0] && face[1] )) {
7607 MESSAGE("Correspondent face not found on side " << ( face[0] ? 1 : 0 ));
7608 if ( nReplaceMap.size() == 2 ) { // faces on input nodes not found
7609 return ( face[0] ? SEW_BAD_SIDE2_NODES : SEW_BAD_SIDE1_NODES );
7612 return SEW_TOPO_DIFF_SETS_OF_ELEMENTS;
7616 // set nodes to merge
7617 // -------------------
7619 if ( face[0] && face[1] ) {
7620 if ( nbNodes[0] != nbNodes[1] ) {
7621 MESSAGE("Diff nb of face nodes");
7622 return SEW_TOPO_DIFF_SETS_OF_ELEMENTS;
7624 #ifdef DEBUG_MATCHING_NODES
7625 MESSAGE ( " Link 1: " << link[0].first->GetID() <<" "<< link[0].second->GetID()
7626 << " F 1: " << face[0] << "| Link 2: " << link[1].first->GetID() <<" "
7627 << link[1].second->GetID() << " F 2: " << face[1] << " | Bind: " ) ;
7629 int nbN = nbNodes[0];
7631 list<const SMDS_MeshNode*>::iterator n1 = notLinkNodes[0].begin();
7632 list<const SMDS_MeshNode*>::iterator n2 = notLinkNodes[1].begin();
7633 for ( int i = 0 ; i < nbN - 2; ++i ) {
7634 #ifdef DEBUG_MATCHING_NODES
7635 MESSAGE ( (*n1)->GetID() << " to " << (*n2)->GetID() );
7637 nReplaceMap.insert( make_pair( *(n1++), *(n2++) ));
7641 // add other links of the face 1 to linkList
7642 // -----------------------------------------
7644 const SMDS_MeshElement* f0 = face[0];
7645 const SMDS_MeshNode* n1 = f0->GetNode( nbN - 1 );
7646 for ( int i = 0; i < nbN; i++ )
7648 const SMDS_MeshNode* n2 = f0->GetNode( i );
7649 pair< set< SMESH_TLink >::iterator, bool > iter_isnew =
7650 linkSet.insert( SMESH_TLink( n1, n2 ));
7651 if ( !iter_isnew.second ) { // already in a set: no need to process
7652 linkSet.erase( iter_isnew.first );
7654 else // new in set == encountered for the first time: add
7656 #ifdef DEBUG_MATCHING_NODES
7657 MESSAGE ( "Add link 1: " << n1->GetID() << " " << n2->GetID() << " "
7658 << " | link 2: " << nReplaceMap[n1]->GetID() << " " << nReplaceMap[n2]->GetID() << " " );
7660 linkList[0].push_back ( NLink( n1, n2 ));
7661 linkList[1].push_back ( NLink( nReplaceMap[n1], nReplaceMap[n2] ));
7666 } // loop on link lists
7672 \brief Creates a hole in a mesh by doubling the nodes of some particular elements
7673 \param theNodes - identifiers of nodes to be doubled
7674 \param theModifiedElems - identifiers of elements to be updated by the new (doubled)
7675 nodes. If list of element identifiers is empty then nodes are doubled but
7676 they not assigned to elements
7677 \return TRUE if operation has been completed successfully, FALSE otherwise
7679 bool SMESH_MeshEditor::DoubleNodes( const std::list< int >& theListOfNodes,
7680 const std::list< int >& theListOfModifiedElems )
7682 myLastCreatedElems.Clear();
7683 myLastCreatedNodes.Clear();
7685 if ( theListOfNodes.size() == 0 )
7688 SMESHDS_Mesh* aMeshDS = GetMeshDS();
7692 // iterate through nodes and duplicate them
7694 std::map< const SMDS_MeshNode*, const SMDS_MeshNode* > anOldNodeToNewNode;
7696 std::list< int >::const_iterator aNodeIter;
7697 for ( aNodeIter = theListOfNodes.begin(); aNodeIter != theListOfNodes.end(); ++aNodeIter )
7699 int aCurr = *aNodeIter;
7700 SMDS_MeshNode* aNode = (SMDS_MeshNode*)aMeshDS->FindNode( aCurr );
7706 const SMDS_MeshNode* aNewNode = aMeshDS->AddNode( aNode->X(), aNode->Y(), aNode->Z() );
7709 anOldNodeToNewNode[ aNode ] = aNewNode;
7710 myLastCreatedNodes.Append( aNewNode );
7714 // Create map of new nodes for modified elements
7716 std::map< SMDS_MeshElement*, vector<const SMDS_MeshNode*> > anElemToNodes;
7718 std::list< int >::const_iterator anElemIter;
7719 for ( anElemIter = theListOfModifiedElems.begin();
7720 anElemIter != theListOfModifiedElems.end(); ++anElemIter )
7722 int aCurr = *anElemIter;
7723 SMDS_MeshElement* anElem = (SMDS_MeshElement*)aMeshDS->FindElement( aCurr );
7727 vector<const SMDS_MeshNode*> aNodeArr( anElem->NbNodes() );
7729 SMDS_ElemIteratorPtr anIter = anElem->nodesIterator();
7731 while ( anIter->more() )
7733 SMDS_MeshNode* aCurrNode = (SMDS_MeshNode*)anIter->next();
7734 if ( aCurr && anOldNodeToNewNode.find( aCurrNode ) != anOldNodeToNewNode.end() )
7736 const SMDS_MeshNode* aNewNode = anOldNodeToNewNode[ aCurrNode ];
7737 aNodeArr[ ind++ ] = aNewNode;
7740 aNodeArr[ ind++ ] = aCurrNode;
7742 anElemToNodes[ anElem ] = aNodeArr;
7745 // Change nodes of elements
7747 std::map< SMDS_MeshElement*, vector<const SMDS_MeshNode*> >::iterator
7748 anElemToNodesIter = anElemToNodes.begin();
7749 for ( ; anElemToNodesIter != anElemToNodes.end(); ++anElemToNodesIter )
7751 const SMDS_MeshElement* anElem = anElemToNodesIter->first;
7752 vector<const SMDS_MeshNode*> aNodeArr = anElemToNodesIter->second;
7754 aMeshDS->ChangeElementNodes( anElem, &aNodeArr[ 0 ], anElem->NbNodes() );