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->GetNodeWrap( iAfter ));
2054 linkedNodes.insert( elem->GetNodeWrap( 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->GetNodeWrap( 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);
2777 vector<bool> issimple(nbNodes);
2779 for ( iNode = 0; iNode < nbNodes; iNode++ ) {
2780 TNodeOfNodeListMapItr nnIt = newNodesItVec[ iNode ];
2781 const SMDS_MeshNode* node = nnIt->first;
2782 const list< const SMDS_MeshNode* > & listNewNodes = nnIt->second;
2783 if ( listNewNodes.empty() )
2786 issimple[iNode] = (listNewNodes.size()==nbSteps);
2788 itNN[ iNode ] = listNewNodes.begin();
2789 prevNod[ iNode ] = node;
2790 nextNod[ iNode ] = listNewNodes.front();
2791 //cout<<"iNode="<<iNode<<endl;
2792 //cout<<" prevNod[iNode]="<< prevNod[iNode]<<" nextNod[iNode]="<< nextNod[iNode]<<endl;
2793 if ( prevNod[ iNode ] != nextNod [ iNode ])
2794 iNotSameNode = iNode;
2798 sames[nbSame++] = iNode;
2801 //cout<<"1 nbSame="<<nbSame<<endl;
2802 if ( nbSame == nbNodes || nbSame > 2) {
2803 MESSAGE( " Too many same nodes of element " << elem->GetID() );
2807 // if( elem->IsQuadratic() && nbSame>0 ) {
2808 // MESSAGE( "Can not rotate quadratic element " << elem->GetID() );
2812 int iBeforeSame = 0, iAfterSame = 0, iOpposSame = 0;
2814 iBeforeSame = ( iSameNode == 0 ? nbNodes - 1 : iSameNode - 1 );
2815 iAfterSame = ( iSameNode + 1 == nbNodes ? 0 : iSameNode + 1 );
2816 iOpposSame = ( iSameNode - 2 < 0 ? iSameNode + 2 : iSameNode - 2 );
2820 //cout<<" prevNod[0]="<< prevNod[0]<<" prevNod[1]="<< prevNod[1]
2821 // <<" prevNod[2]="<< prevNod[2]<<" prevNod[3]="<< prevNod[4]
2822 // <<" prevNod[4]="<< prevNod[4]<<" prevNod[5]="<< prevNod[5]
2823 // <<" prevNod[6]="<< prevNod[6]<<" prevNod[7]="<< prevNod[7]<<endl;
2825 // check element orientation
2827 if ( nbNodes > 2 && !isReverse( prevNod, nextNod, nbNodes, iNotSameNode )) {
2828 //MESSAGE("Reversed elem " << elem );
2832 std::swap( iBeforeSame, iAfterSame );
2835 // make new elements
2836 for (int iStep = 0; iStep < nbSteps; iStep++ ) {
2838 for ( iNode = 0; iNode < nbNodes; iNode++ ) {
2839 if(issimple[iNode]) {
2840 nextNod[ iNode ] = *itNN[ iNode ];
2844 if( elem->GetType()==SMDSAbs_Node ) {
2845 // we have to use two nodes
2846 midlNod[ iNode ] = *itNN[ iNode ];
2848 nextNod[ iNode ] = *itNN[ iNode ];
2851 else if(!elem->IsQuadratic() || elem->IsMediumNode(prevNod[iNode]) ) {
2852 // we have to use each second node
2854 nextNod[ iNode ] = *itNN[ iNode ];
2858 // we have to use two nodes
2859 midlNod[ iNode ] = *itNN[ iNode ];
2861 nextNod[ iNode ] = *itNN[ iNode ];
2866 SMDS_MeshElement* aNewElem = 0;
2867 if(!elem->IsPoly()) {
2868 switch ( nbNodes ) {
2872 if ( nbSame == 0 ) {
2874 aNewElem = aMesh->AddEdge( prevNod[ 0 ], nextNod[ 0 ] );
2876 aNewElem = aMesh->AddEdge( prevNod[ 0 ], nextNod[ 0 ], midlNod[ 0 ] );
2882 aNewElem = aMesh->AddFace(prevNod[ 0 ], prevNod[ 1 ],
2883 nextNod[ 1 ], nextNod[ 0 ] );
2885 aNewElem = aMesh->AddFace(prevNod[ 0 ], prevNod[ 1 ],
2886 nextNod[ iNotSameNode ] );
2890 case 3: { // TRIANGLE or quadratic edge
2891 if(elem->GetType() == SMDSAbs_Face) { // TRIANGLE
2893 if ( nbSame == 0 ) // --- pentahedron
2894 aNewElem = aMesh->AddVolume (prevNod[ i0 ], prevNod[ 1 ], prevNod[ i2 ],
2895 nextNod[ i0 ], nextNod[ 1 ], nextNod[ i2 ] );
2897 else if ( nbSame == 1 ) // --- pyramid
2898 aNewElem = aMesh->AddVolume (prevNod[ iBeforeSame ], prevNod[ iAfterSame ],
2899 nextNod[ iAfterSame ], nextNod[ iBeforeSame ],
2900 nextNod[ iSameNode ]);
2902 else // 2 same nodes: --- tetrahedron
2903 aNewElem = aMesh->AddVolume (prevNod[ i0 ], prevNod[ 1 ], prevNod[ i2 ],
2904 nextNod[ iNotSameNode ]);
2906 else { // quadratic edge
2907 if(nbSame==0) { // quadratic quadrangle
2908 aNewElem = aMesh->AddFace(prevNod[0], nextNod[0], nextNod[1], prevNod[1],
2909 midlNod[0], nextNod[2], midlNod[1], prevNod[2]);
2911 else if(nbSame==1) { // quadratic triangle
2913 return; // medium node on axis
2914 else if(sames[0]==0) {
2915 aNewElem = aMesh->AddFace(prevNod[0], nextNod[1], prevNod[1],
2916 nextNod[2], midlNod[1], prevNod[2]);
2918 else { // sames[0]==1
2919 aNewElem = aMesh->AddFace(prevNod[0], nextNod[0], prevNod[1],
2920 midlNod[0], nextNod[2], prevNod[2]);
2928 case 4: { // QUADRANGLE
2930 if ( nbSame == 0 ) // --- hexahedron
2931 aNewElem = aMesh->AddVolume (prevNod[ i0 ], prevNod[ 1 ], prevNod[ i2 ], prevNod[ 3 ],
2932 nextNod[ i0 ], nextNod[ 1 ], nextNod[ i2 ], nextNod[ 3 ]);
2934 else if ( nbSame == 1 ) { // --- pyramid + pentahedron
2935 aNewElem = aMesh->AddVolume (prevNod[ iBeforeSame ], prevNod[ iAfterSame ],
2936 nextNod[ iAfterSame ], nextNod[ iBeforeSame ],
2937 nextNod[ iSameNode ]);
2938 newElems.push_back( aNewElem );
2939 aNewElem = aMesh->AddVolume (prevNod[ iAfterSame ], prevNod[ iOpposSame ],
2940 prevNod[ iBeforeSame ], nextNod[ iAfterSame ],
2941 nextNod[ iOpposSame ], nextNod[ iBeforeSame ] );
2943 else if ( nbSame == 2 ) { // pentahedron
2944 if ( prevNod[ iBeforeSame ] == nextNod[ iBeforeSame ] )
2945 // iBeforeSame is same too
2946 aNewElem = aMesh->AddVolume (prevNod[ iBeforeSame ], prevNod[ iOpposSame ],
2947 nextNod[ iOpposSame ], prevNod[ iSameNode ],
2948 prevNod[ iAfterSame ], nextNod[ iAfterSame ]);
2950 // iAfterSame is same too
2951 aNewElem = aMesh->AddVolume (prevNod[ iSameNode ], prevNod[ iBeforeSame ],
2952 nextNod[ iBeforeSame ], prevNod[ iAfterSame ],
2953 prevNod[ iOpposSame ], nextNod[ iOpposSame ]);
2957 case 6: { // quadratic triangle
2958 // create pentahedron with 15 nodes
2959 if(i0>0) { // reversed case
2960 aNewElem = aMesh->AddVolume (prevNod[0], prevNod[2], prevNod[1],
2961 nextNod[0], nextNod[2], nextNod[1],
2962 prevNod[5], prevNod[4], prevNod[3],
2963 nextNod[5], nextNod[4], nextNod[3],
2964 midlNod[0], midlNod[2], midlNod[1]);
2966 else { // not reversed case
2967 aNewElem = aMesh->AddVolume (prevNod[0], prevNod[1], prevNod[2],
2968 nextNod[0], nextNod[1], nextNod[2],
2969 prevNod[3], prevNod[4], prevNod[5],
2970 nextNod[3], nextNod[4], nextNod[5],
2971 midlNod[0], midlNod[1], midlNod[2]);
2975 case 8: { // quadratic quadrangle
2976 // create hexahedron with 20 nodes
2977 if(i0>0) { // reversed case
2978 aNewElem = aMesh->AddVolume (prevNod[0], prevNod[3], prevNod[2], prevNod[1],
2979 nextNod[0], nextNod[3], nextNod[2], nextNod[1],
2980 prevNod[7], prevNod[6], prevNod[5], prevNod[4],
2981 nextNod[7], nextNod[6], nextNod[5], nextNod[4],
2982 midlNod[0], midlNod[3], midlNod[2], midlNod[1]);
2984 else { // not reversed case
2985 aNewElem = aMesh->AddVolume (prevNod[0], prevNod[1], prevNod[2], prevNod[3],
2986 nextNod[0], nextNod[1], nextNod[2], nextNod[3],
2987 prevNod[4], prevNod[5], prevNod[6], prevNod[7],
2988 nextNod[4], nextNod[5], nextNod[6], nextNod[7],
2989 midlNod[0], midlNod[1], midlNod[2], midlNod[3]);
2994 // realized for extrusion only
2995 //vector<const SMDS_MeshNode*> polyedre_nodes (nbNodes*2 + 4*nbNodes);
2996 //vector<int> quantities (nbNodes + 2);
2998 //quantities[0] = nbNodes; // bottom of prism
2999 //for (int inode = 0; inode < nbNodes; inode++) {
3000 // polyedre_nodes[inode] = prevNod[inode];
3003 //quantities[1] = nbNodes; // top of prism
3004 //for (int inode = 0; inode < nbNodes; inode++) {
3005 // polyedre_nodes[nbNodes + inode] = nextNod[inode];
3008 //for (int iface = 0; iface < nbNodes; iface++) {
3009 // quantities[iface + 2] = 4;
3010 // int inextface = (iface == nbNodes - 1) ? 0 : iface + 1;
3011 // polyedre_nodes[2*nbNodes + 4*iface + 0] = prevNod[iface];
3012 // polyedre_nodes[2*nbNodes + 4*iface + 1] = prevNod[inextface];
3013 // polyedre_nodes[2*nbNodes + 4*iface + 2] = nextNod[inextface];
3014 // polyedre_nodes[2*nbNodes + 4*iface + 3] = nextNod[iface];
3016 //aNewElem = aMesh->AddPolyhedralVolume (polyedre_nodes, quantities);
3023 // realized for extrusion only
3024 vector<const SMDS_MeshNode*> polyedre_nodes (nbNodes*2 + 4*nbNodes);
3025 vector<int> quantities (nbNodes + 2);
3027 quantities[0] = nbNodes; // bottom of prism
3028 for (int inode = 0; inode < nbNodes; inode++) {
3029 polyedre_nodes[inode] = prevNod[inode];
3032 quantities[1] = nbNodes; // top of prism
3033 for (int inode = 0; inode < nbNodes; inode++) {
3034 polyedre_nodes[nbNodes + inode] = nextNod[inode];
3037 for (int iface = 0; iface < nbNodes; iface++) {
3038 quantities[iface + 2] = 4;
3039 int inextface = (iface == nbNodes - 1) ? 0 : iface + 1;
3040 polyedre_nodes[2*nbNodes + 4*iface + 0] = prevNod[iface];
3041 polyedre_nodes[2*nbNodes + 4*iface + 1] = prevNod[inextface];
3042 polyedre_nodes[2*nbNodes + 4*iface + 2] = nextNod[inextface];
3043 polyedre_nodes[2*nbNodes + 4*iface + 3] = nextNod[iface];
3045 aNewElem = aMesh->AddPolyhedralVolume (polyedre_nodes, quantities);
3049 newElems.push_back( aNewElem );
3050 myLastCreatedElems.Append(aNewElem);
3051 srcElements.Append( elem );
3054 // set new prev nodes
3055 for ( iNode = 0; iNode < nbNodes; iNode++ )
3056 prevNod[ iNode ] = nextNod[ iNode ];
3061 //=======================================================================
3063 * \brief Create 1D and 2D elements around swept elements
3064 * \param mapNewNodes - source nodes and ones generated from them
3065 * \param newElemsMap - source elements and ones generated from them
3066 * \param elemNewNodesMap - nodes generated from each node of each element
3067 * \param elemSet - all swept elements
3068 * \param nbSteps - number of sweeping steps
3069 * \param srcElements - to append elem for each generated element
3071 //=======================================================================
3073 void SMESH_MeshEditor::makeWalls (TNodeOfNodeListMap & mapNewNodes,
3074 TElemOfElemListMap & newElemsMap,
3075 TElemOfVecOfNnlmiMap & elemNewNodesMap,
3076 TIDSortedElemSet& elemSet,
3078 SMESH_SequenceOfElemPtr& srcElements)
3080 ASSERT( newElemsMap.size() == elemNewNodesMap.size() );
3081 SMESHDS_Mesh* aMesh = GetMeshDS();
3083 // Find nodes belonging to only one initial element - sweep them to get edges.
3085 TNodeOfNodeListMapItr nList = mapNewNodes.begin();
3086 for ( ; nList != mapNewNodes.end(); nList++ ) {
3087 const SMDS_MeshNode* node =
3088 static_cast<const SMDS_MeshNode*>( nList->first );
3089 SMDS_ElemIteratorPtr eIt = node->GetInverseElementIterator();
3090 int nbInitElems = 0;
3091 const SMDS_MeshElement* el = 0;
3092 SMDSAbs_ElementType highType = SMDSAbs_Edge; // count most complex elements only
3093 while ( eIt->more() && nbInitElems < 2 ) {
3095 SMDSAbs_ElementType type = el->GetType();
3096 if ( type == SMDSAbs_Volume || type < highType ) continue;
3097 if ( type > highType ) {
3101 if ( elemSet.find(el) != elemSet.end() )
3104 if ( nbInitElems < 2 ) {
3105 bool NotCreateEdge = el && el->IsQuadratic() && el->IsMediumNode(node);
3106 if(!NotCreateEdge) {
3107 vector<TNodeOfNodeListMapItr> newNodesItVec( 1, nList );
3108 list<const SMDS_MeshElement*> newEdges;
3109 sweepElement( node, newNodesItVec, newEdges, nbSteps, srcElements );
3114 // Make a ceiling for each element ie an equal element of last new nodes.
3115 // Find free links of faces - make edges and sweep them into faces.
3117 TElemOfElemListMap::iterator itElem = newElemsMap.begin();
3118 TElemOfVecOfNnlmiMap::iterator itElemNodes = elemNewNodesMap.begin();
3119 for ( ; itElem != newElemsMap.end(); itElem++, itElemNodes++ ) {
3120 const SMDS_MeshElement* elem = itElem->first;
3121 vector<TNodeOfNodeListMapItr>& vecNewNodes = itElemNodes->second;
3123 if ( elem->GetType() == SMDSAbs_Edge ) {
3124 // create a ceiling edge
3125 if (!elem->IsQuadratic()) {
3126 if ( !aMesh->FindEdge( vecNewNodes[ 0 ]->second.back(),
3127 vecNewNodes[ 1 ]->second.back())) {
3128 myLastCreatedElems.Append(aMesh->AddEdge(vecNewNodes[ 0 ]->second.back(),
3129 vecNewNodes[ 1 ]->second.back()));
3130 srcElements.Append( myLastCreatedElems.Last() );
3134 if ( !aMesh->FindEdge( vecNewNodes[ 0 ]->second.back(),
3135 vecNewNodes[ 1 ]->second.back(),
3136 vecNewNodes[ 2 ]->second.back())) {
3137 myLastCreatedElems.Append(aMesh->AddEdge(vecNewNodes[ 0 ]->second.back(),
3138 vecNewNodes[ 1 ]->second.back(),
3139 vecNewNodes[ 2 ]->second.back()));
3140 srcElements.Append( myLastCreatedElems.Last() );
3144 if ( elem->GetType() != SMDSAbs_Face )
3147 if(itElem->second.size()==0) continue;
3149 bool hasFreeLinks = false;
3151 TIDSortedElemSet avoidSet;
3152 avoidSet.insert( elem );
3154 set<const SMDS_MeshNode*> aFaceLastNodes;
3155 int iNode, nbNodes = vecNewNodes.size();
3156 if(!elem->IsQuadratic()) {
3157 // loop on the face nodes
3158 for ( iNode = 0; iNode < nbNodes; iNode++ ) {
3159 aFaceLastNodes.insert( vecNewNodes[ iNode ]->second.back() );
3160 // look for free links of the face
3161 int iNext = ( iNode + 1 == nbNodes ) ? 0 : iNode + 1;
3162 const SMDS_MeshNode* n1 = vecNewNodes[ iNode ]->first;
3163 const SMDS_MeshNode* n2 = vecNewNodes[ iNext ]->first;
3164 // check if a link is free
3165 if ( ! SMESH_MeshEditor::FindFaceInSet ( n1, n2, elemSet, avoidSet )) {
3166 hasFreeLinks = true;
3167 // make an edge and a ceiling for a new edge
3168 if ( !aMesh->FindEdge( n1, n2 )) {
3169 myLastCreatedElems.Append(aMesh->AddEdge( n1, n2 )); // free link edge
3170 srcElements.Append( myLastCreatedElems.Last() );
3172 n1 = vecNewNodes[ iNode ]->second.back();
3173 n2 = vecNewNodes[ iNext ]->second.back();
3174 if ( !aMesh->FindEdge( n1, n2 )) {
3175 myLastCreatedElems.Append(aMesh->AddEdge( n1, n2 )); // ceiling edge
3176 srcElements.Append( myLastCreatedElems.Last() );
3181 else { // elem is quadratic face
3182 int nbn = nbNodes/2;
3183 for ( iNode = 0; iNode < nbn; iNode++ ) {
3184 aFaceLastNodes.insert( vecNewNodes[ iNode ]->second.back() );
3185 int iNext = ( iNode + 1 == nbn ) ? 0 : iNode + 1;
3186 const SMDS_MeshNode* n1 = vecNewNodes[ iNode ]->first;
3187 const SMDS_MeshNode* n2 = vecNewNodes[ iNext ]->first;
3188 // check if a link is free
3189 if ( ! SMESH_MeshEditor::FindFaceInSet ( n1, n2, elemSet, avoidSet )) {
3190 hasFreeLinks = true;
3191 // make an edge and a ceiling for a new edge
3193 const SMDS_MeshNode* n3 = vecNewNodes[ iNode+nbn ]->first;
3194 if ( !aMesh->FindEdge( n1, n2, n3 )) {
3195 myLastCreatedElems.Append(aMesh->AddEdge( n1, n2, n3 )); // free link edge
3196 srcElements.Append( myLastCreatedElems.Last() );
3198 n1 = vecNewNodes[ iNode ]->second.back();
3199 n2 = vecNewNodes[ iNext ]->second.back();
3200 n3 = vecNewNodes[ iNode+nbn ]->second.back();
3201 if ( !aMesh->FindEdge( n1, n2, n3 )) {
3202 myLastCreatedElems.Append(aMesh->AddEdge( n1, n2, n3 )); // ceiling edge
3203 srcElements.Append( myLastCreatedElems.Last() );
3207 for ( iNode = nbn; iNode < 2*nbn; iNode++ ) {
3208 aFaceLastNodes.insert( vecNewNodes[ iNode ]->second.back() );
3212 // sweep free links into faces
3214 if ( hasFreeLinks ) {
3215 list<const SMDS_MeshElement*> & newVolumes = itElem->second;
3216 int iVol, volNb, nbVolumesByStep = newVolumes.size() / nbSteps;
3218 set<const SMDS_MeshNode*> initNodeSet, topNodeSet, faceNodeSet;
3219 for ( iNode = 0; iNode < nbNodes; iNode++ ) {
3220 initNodeSet.insert( vecNewNodes[ iNode ]->first );
3221 topNodeSet .insert( vecNewNodes[ iNode ]->second.back() );
3223 for ( volNb = 0; volNb < nbVolumesByStep; volNb++ ) {
3224 list<const SMDS_MeshElement*>::iterator v = newVolumes.begin();
3226 while ( iVol++ < volNb ) v++;
3227 // find indices of free faces of a volume and their source edges
3228 list< int > freeInd;
3229 list< const SMDS_MeshElement* > srcEdges; // source edges of free faces
3230 SMDS_VolumeTool vTool( *v );
3231 int iF, nbF = vTool.NbFaces();
3232 for ( iF = 0; iF < nbF; iF ++ ) {
3233 if (vTool.IsFreeFace( iF ) &&
3234 vTool.GetFaceNodes( iF, faceNodeSet ) &&
3235 initNodeSet != faceNodeSet) // except an initial face
3237 if ( nbSteps == 1 && faceNodeSet == topNodeSet )
3239 freeInd.push_back( iF );
3240 // find source edge of a free face iF
3241 vector<const SMDS_MeshNode*> commonNodes; // shared by the initial and free faces
3242 commonNodes.resize( initNodeSet.size(), NULL ); // avoid spoiling memory
3243 std::set_intersection( faceNodeSet.begin(), faceNodeSet.end(),
3244 initNodeSet.begin(), initNodeSet.end(),
3245 commonNodes.begin());
3246 if ( (*v)->IsQuadratic() )
3247 srcEdges.push_back(aMesh->FindEdge (commonNodes[0],commonNodes[1],commonNodes[2]));
3249 srcEdges.push_back(aMesh->FindEdge (commonNodes[0],commonNodes[1]));
3251 if ( !srcEdges.back() )
3253 cout << "SMESH_MeshEditor::makeWalls(), no source edge found for a free face #"
3254 << iF << " of volume #" << vTool.ID() << endl;
3259 if ( freeInd.empty() )
3262 // create faces for all steps;
3263 // if such a face has been already created by sweep of edge,
3264 // assure that its orientation is OK
3265 for ( int iStep = 0; iStep < nbSteps; iStep++ ) {
3267 vTool.SetExternalNormal();
3268 list< int >::iterator ind = freeInd.begin();
3269 list< const SMDS_MeshElement* >::iterator srcEdge = srcEdges.begin();
3270 for ( ; ind != freeInd.end(); ++ind, ++srcEdge ) // loop on free faces
3272 const SMDS_MeshNode** nodes = vTool.GetFaceNodes( *ind );
3273 int nbn = vTool.NbFaceNodes( *ind );
3275 case 3: { ///// triangle
3276 const SMDS_MeshFace * f = aMesh->FindFace( nodes[ 0 ], nodes[ 1 ], nodes[ 2 ]);
3278 myLastCreatedElems.Append(aMesh->AddFace( nodes[ 0 ], nodes[ 1 ], nodes[ 2 ] ));
3279 else if ( nodes[ 1 ] != f->GetNodeWrap( f->GetNodeIndex( nodes[ 0 ] ) + 1 ))
3280 aMesh->ChangeElementNodes( f, nodes, nbn );
3283 case 4: { ///// quadrangle
3284 const SMDS_MeshFace * f = aMesh->FindFace( nodes[ 0 ], nodes[ 1 ], nodes[ 2 ], nodes[ 3 ]);
3286 myLastCreatedElems.Append(aMesh->AddFace( nodes[ 0 ], nodes[ 1 ], nodes[ 2 ], nodes[ 3 ] ));
3287 else if ( nodes[ 1 ] != f->GetNodeWrap( f->GetNodeIndex( nodes[ 0 ] ) + 1 ))
3288 aMesh->ChangeElementNodes( f, nodes, nbn );
3292 if( (*v)->IsQuadratic() ) {
3293 if(nbn==6) { /////// quadratic triangle
3294 const SMDS_MeshFace * f = aMesh->FindFace( nodes[0], nodes[2], nodes[4],
3295 nodes[1], nodes[3], nodes[5] );
3297 myLastCreatedElems.Append(aMesh->AddFace(nodes[0], nodes[2], nodes[4],
3298 nodes[1], nodes[3], nodes[5]));
3299 else if ( nodes[ 2 ] != f->GetNodeWrap( f->GetNodeIndex( nodes[ 0 ] ) + 1 ))
3300 aMesh->ChangeElementNodes( f, nodes, nbn );
3302 else { /////// quadratic quadrangle
3303 const SMDS_MeshFace * f = aMesh->FindFace( nodes[0], nodes[2], nodes[4], nodes[6],
3304 nodes[1], nodes[3], nodes[5], nodes[7] );
3306 myLastCreatedElems.Append(aMesh->AddFace(nodes[0], nodes[2], nodes[4], nodes[6],
3307 nodes[1], nodes[3], nodes[5], nodes[7]));
3308 else if ( nodes[ 2 ] != f->GetNodeWrap( f->GetNodeIndex( nodes[ 0 ] ) + 1 ))
3309 aMesh->ChangeElementNodes( f, nodes, nbn );
3312 else { //////// polygon
3313 vector<const SMDS_MeshNode*> polygon_nodes ( nodes, &nodes[nbn] );
3314 const SMDS_MeshFace * f = aMesh->FindFace( polygon_nodes );
3316 myLastCreatedElems.Append(aMesh->AddPolygonalFace(polygon_nodes));
3317 else if ( nodes[ 1 ] != f->GetNodeWrap( f->GetNodeIndex( nodes[ 0 ] ) + 1 ))
3318 aMesh->ChangeElementNodes( f, nodes, nbn );
3321 while ( srcElements.Length() < myLastCreatedElems.Length() )
3322 srcElements.Append( *srcEdge );
3324 } // loop on free faces
3326 // go to the next volume
3328 while ( iVol++ < nbVolumesByStep ) v++;
3331 } // sweep free links into faces
3333 // Make a ceiling face with a normal external to a volume
3335 SMDS_VolumeTool lastVol( itElem->second.back() );
3337 int iF = lastVol.GetFaceIndex( aFaceLastNodes );
3339 lastVol.SetExternalNormal();
3340 const SMDS_MeshNode** nodes = lastVol.GetFaceNodes( iF );
3341 int nbn = lastVol.NbFaceNodes( iF );
3344 if (!hasFreeLinks ||
3345 !aMesh->FindFace( nodes[ 0 ], nodes[ 1 ], nodes[ 2 ]))
3346 myLastCreatedElems.Append(aMesh->AddFace( nodes[ 0 ], nodes[ 1 ], nodes[ 2 ] ));
3349 if (!hasFreeLinks ||
3350 !aMesh->FindFace( nodes[ 0 ], nodes[ 1 ], nodes[ 2 ], nodes[ 3 ]))
3351 myLastCreatedElems.Append(aMesh->AddFace( nodes[ 0 ], nodes[ 1 ], nodes[ 2 ], nodes[ 3 ] ));
3354 if(itElem->second.back()->IsQuadratic()) {
3356 if (!hasFreeLinks ||
3357 !aMesh->FindFace(nodes[0], nodes[2], nodes[4],
3358 nodes[1], nodes[3], nodes[5]) ) {
3359 myLastCreatedElems.Append(aMesh->AddFace(nodes[0], nodes[2], nodes[4],
3360 nodes[1], nodes[3], nodes[5]));
3364 if (!hasFreeLinks ||
3365 !aMesh->FindFace(nodes[0], nodes[2], nodes[4], nodes[6],
3366 nodes[1], nodes[3], nodes[5], nodes[7]) )
3367 myLastCreatedElems.Append(aMesh->AddFace(nodes[0], nodes[2], nodes[4], nodes[6],
3368 nodes[1], nodes[3], nodes[5], nodes[7]));
3372 vector<const SMDS_MeshNode*> polygon_nodes ( nodes, &nodes[nbn] );
3373 if (!hasFreeLinks || !aMesh->FindFace(polygon_nodes))
3374 myLastCreatedElems.Append(aMesh->AddPolygonalFace(polygon_nodes));
3378 while ( srcElements.Length() < myLastCreatedElems.Length() )
3379 srcElements.Append( myLastCreatedElems.Last() );
3381 } // loop on swept elements
3384 //=======================================================================
3385 //function : RotationSweep
3387 //=======================================================================
3389 SMESH_MeshEditor::PGroupIDs
3390 SMESH_MeshEditor::RotationSweep(TIDSortedElemSet & theElems,
3391 const gp_Ax1& theAxis,
3392 const double theAngle,
3393 const int theNbSteps,
3394 const double theTol,
3395 const bool theMakeGroups,
3396 const bool theMakeWalls)
3398 myLastCreatedElems.Clear();
3399 myLastCreatedNodes.Clear();
3401 // source elements for each generated one
3402 SMESH_SequenceOfElemPtr srcElems, srcNodes;
3404 MESSAGE( "RotationSweep()");
3406 aTrsf.SetRotation( theAxis, theAngle );
3408 aTrsf2.SetRotation( theAxis, theAngle/2. );
3410 gp_Lin aLine( theAxis );
3411 double aSqTol = theTol * theTol;
3413 SMESHDS_Mesh* aMesh = GetMeshDS();
3415 TNodeOfNodeListMap mapNewNodes;
3416 TElemOfVecOfNnlmiMap mapElemNewNodes;
3417 TElemOfElemListMap newElemsMap;
3420 TIDSortedElemSet::iterator itElem;
3421 for ( itElem = theElems.begin(); itElem != theElems.end(); itElem++ ) {
3422 const SMDS_MeshElement* elem = *itElem;
3423 if ( !elem || elem->GetType() == SMDSAbs_Volume )
3425 vector<TNodeOfNodeListMapItr> & newNodesItVec = mapElemNewNodes[ elem ];
3426 newNodesItVec.reserve( elem->NbNodes() );
3428 // loop on elem nodes
3429 SMDS_ElemIteratorPtr itN = elem->nodesIterator();
3430 while ( itN->more() )
3432 // check if a node has been already sweeped
3433 const SMDS_MeshNode* node = cast2Node( itN->next() );
3434 TNodeOfNodeListMapItr nIt = mapNewNodes.find( node );
3435 if ( nIt == mapNewNodes.end() ) {
3436 nIt = mapNewNodes.insert( make_pair( node, list<const SMDS_MeshNode*>() )).first;
3437 list<const SMDS_MeshNode*>& listNewNodes = nIt->second;
3440 gp_XYZ aXYZ( node->X(), node->Y(), node->Z() );
3442 aXYZ.Coord( coord[0], coord[1], coord[2] );
3443 bool isOnAxis = ( aLine.SquareDistance( aXYZ ) <= aSqTol );
3444 const SMDS_MeshNode * newNode = node;
3445 for ( int i = 0; i < theNbSteps; i++ ) {
3447 if( elem->IsQuadratic() && !elem->IsMediumNode(node) ) {
3449 aTrsf2.Transforms( coord[0], coord[1], coord[2] );
3450 //aTrsf.Transforms( coord[0], coord[1], coord[2] );
3451 newNode = aMesh->AddNode( coord[0], coord[1], coord[2] );
3452 myLastCreatedNodes.Append(newNode);
3453 srcNodes.Append( node );
3454 listNewNodes.push_back( newNode );
3455 aTrsf2.Transforms( coord[0], coord[1], coord[2] );
3456 //aTrsf.Transforms( coord[0], coord[1], coord[2] );
3459 aTrsf.Transforms( coord[0], coord[1], coord[2] );
3461 newNode = aMesh->AddNode( coord[0], coord[1], coord[2] );
3462 myLastCreatedNodes.Append(newNode);
3463 srcNodes.Append( node );
3465 listNewNodes.push_back( newNode );
3469 // if current elem is quadratic and current node is not medium
3470 // we have to check - may be it is needed to insert additional nodes
3471 if( elem->IsQuadratic() && !elem->IsMediumNode(node) ) {
3472 list< const SMDS_MeshNode* > & listNewNodes = nIt->second;
3473 if(listNewNodes.size()==theNbSteps) {
3474 listNewNodes.clear();
3476 gp_XYZ aXYZ( node->X(), node->Y(), node->Z() );
3478 aXYZ.Coord( coord[0], coord[1], coord[2] );
3479 const SMDS_MeshNode * newNode = node;
3480 for(int i = 0; i<theNbSteps; i++) {
3481 aTrsf2.Transforms( coord[0], coord[1], coord[2] );
3482 newNode = aMesh->AddNode( coord[0], coord[1], coord[2] );
3483 myLastCreatedNodes.Append(newNode);
3484 listNewNodes.push_back( newNode );
3485 srcNodes.Append( node );
3486 aTrsf2.Transforms( coord[0], coord[1], coord[2] );
3487 newNode = aMesh->AddNode( coord[0], coord[1], coord[2] );
3488 myLastCreatedNodes.Append(newNode);
3489 srcNodes.Append( node );
3490 listNewNodes.push_back( newNode );
3495 newNodesItVec.push_back( nIt );
3497 // make new elements
3498 sweepElement( elem, newNodesItVec, newElemsMap[elem], theNbSteps, srcElems );
3502 makeWalls( mapNewNodes, newElemsMap, mapElemNewNodes, theElems, theNbSteps, srcElems );
3504 PGroupIDs newGroupIDs;
3505 if ( theMakeGroups )
3506 newGroupIDs = generateGroups( srcNodes, srcElems, "rotated");
3512 //=======================================================================
3513 //function : CreateNode
3515 //=======================================================================
3516 const SMDS_MeshNode* SMESH_MeshEditor::CreateNode(const double x,
3519 const double tolnode,
3520 SMESH_SequenceOfNode& aNodes)
3522 myLastCreatedElems.Clear();
3523 myLastCreatedNodes.Clear();
3526 SMESHDS_Mesh * aMesh = myMesh->GetMeshDS();
3528 // try to search in sequence of existing nodes
3529 // if aNodes.Length()>0 we 'nave to use given sequence
3530 // else - use all nodes of mesh
3531 if(aNodes.Length()>0) {
3533 for(i=1; i<=aNodes.Length(); i++) {
3534 gp_Pnt P2(aNodes.Value(i)->X(),aNodes.Value(i)->Y(),aNodes.Value(i)->Z());
3535 if(P1.Distance(P2)<tolnode)
3536 return aNodes.Value(i);
3540 SMDS_NodeIteratorPtr itn = aMesh->nodesIterator();
3541 while(itn->more()) {
3542 const SMDS_MeshNode* aN = static_cast<const SMDS_MeshNode*> (itn->next());
3543 gp_Pnt P2(aN->X(),aN->Y(),aN->Z());
3544 if(P1.Distance(P2)<tolnode)
3549 // create new node and return it
3550 const SMDS_MeshNode* NewNode = aMesh->AddNode(x,y,z);
3551 myLastCreatedNodes.Append(NewNode);
3556 //=======================================================================
3557 //function : ExtrusionSweep
3559 //=======================================================================
3561 SMESH_MeshEditor::PGroupIDs
3562 SMESH_MeshEditor::ExtrusionSweep (TIDSortedElemSet & theElems,
3563 const gp_Vec& theStep,
3564 const int theNbSteps,
3565 TElemOfElemListMap& newElemsMap,
3566 const bool theMakeGroups,
3568 const double theTolerance)
3570 ExtrusParam aParams;
3571 aParams.myDir = gp_Dir(theStep);
3572 aParams.myNodes.Clear();
3573 aParams.mySteps = new TColStd_HSequenceOfReal;
3575 for(i=1; i<=theNbSteps; i++)
3576 aParams.mySteps->Append(theStep.Magnitude());
3579 ExtrusionSweep(theElems,aParams,newElemsMap,theMakeGroups,theFlags,theTolerance);
3583 //=======================================================================
3584 //function : ExtrusionSweep
3586 //=======================================================================
3588 SMESH_MeshEditor::PGroupIDs
3589 SMESH_MeshEditor::ExtrusionSweep (TIDSortedElemSet & theElems,
3590 ExtrusParam& theParams,
3591 TElemOfElemListMap& newElemsMap,
3592 const bool theMakeGroups,
3594 const double theTolerance)
3596 myLastCreatedElems.Clear();
3597 myLastCreatedNodes.Clear();
3599 // source elements for each generated one
3600 SMESH_SequenceOfElemPtr srcElems, srcNodes;
3602 SMESHDS_Mesh* aMesh = GetMeshDS();
3604 int nbsteps = theParams.mySteps->Length();
3606 TNodeOfNodeListMap mapNewNodes;
3607 //TNodeOfNodeVecMap mapNewNodes;
3608 TElemOfVecOfNnlmiMap mapElemNewNodes;
3609 //TElemOfVecOfMapNodesMap mapElemNewNodes;
3612 TIDSortedElemSet::iterator itElem;
3613 for ( itElem = theElems.begin(); itElem != theElems.end(); itElem++ ) {
3614 // check element type
3615 const SMDS_MeshElement* elem = *itElem;
3616 if ( !elem || elem->GetType() == SMDSAbs_Volume )
3619 vector<TNodeOfNodeListMapItr> & newNodesItVec = mapElemNewNodes[ elem ];
3620 //vector<TNodeOfNodeVecMapItr> & newNodesItVec = mapElemNewNodes[ elem ];
3621 newNodesItVec.reserve( elem->NbNodes() );
3623 // loop on elem nodes
3624 SMDS_ElemIteratorPtr itN = elem->nodesIterator();
3625 while ( itN->more() )
3627 // check if a node has been already sweeped
3628 const SMDS_MeshNode* node = cast2Node( itN->next() );
3629 TNodeOfNodeListMap::iterator nIt = mapNewNodes.find( node );
3630 //TNodeOfNodeVecMap::iterator nIt = mapNewNodes.find( node );
3631 if ( nIt == mapNewNodes.end() ) {
3632 nIt = mapNewNodes.insert( make_pair( node, list<const SMDS_MeshNode*>() )).first;
3633 //nIt = mapNewNodes.insert( make_pair( node, vector<const SMDS_MeshNode*>() )).first;
3634 list<const SMDS_MeshNode*>& listNewNodes = nIt->second;
3635 //vector<const SMDS_MeshNode*>& vecNewNodes = nIt->second;
3636 //vecNewNodes.reserve(nbsteps);
3639 double coord[] = { node->X(), node->Y(), node->Z() };
3640 //int nbsteps = theParams.mySteps->Length();
3641 for ( int i = 0; i < nbsteps; i++ ) {
3642 if( elem->IsQuadratic() && !elem->IsMediumNode(node) ) {
3643 // create additional node
3644 double x = coord[0] + theParams.myDir.X()*theParams.mySteps->Value(i+1)/2.;
3645 double y = coord[1] + theParams.myDir.Y()*theParams.mySteps->Value(i+1)/2.;
3646 double z = coord[2] + theParams.myDir.Z()*theParams.mySteps->Value(i+1)/2.;
3647 if( theFlags & EXTRUSION_FLAG_SEW ) {
3648 const SMDS_MeshNode * newNode = CreateNode(x, y, z,
3649 theTolerance, theParams.myNodes);
3650 listNewNodes.push_back( newNode );
3653 const SMDS_MeshNode * newNode = aMesh->AddNode(x, y, z);
3654 myLastCreatedNodes.Append(newNode);
3655 srcNodes.Append( node );
3656 listNewNodes.push_back( newNode );
3659 //aTrsf.Transforms( coord[0], coord[1], coord[2] );
3660 coord[0] = coord[0] + theParams.myDir.X()*theParams.mySteps->Value(i+1);
3661 coord[1] = coord[1] + theParams.myDir.Y()*theParams.mySteps->Value(i+1);
3662 coord[2] = coord[2] + theParams.myDir.Z()*theParams.mySteps->Value(i+1);
3663 if( theFlags & EXTRUSION_FLAG_SEW ) {
3664 const SMDS_MeshNode * newNode = CreateNode(coord[0], coord[1], coord[2],
3665 theTolerance, theParams.myNodes);
3666 listNewNodes.push_back( newNode );
3667 //vecNewNodes[i]=newNode;
3670 const SMDS_MeshNode * newNode = aMesh->AddNode( coord[0], coord[1], coord[2] );
3671 myLastCreatedNodes.Append(newNode);
3672 srcNodes.Append( node );
3673 listNewNodes.push_back( newNode );
3674 //vecNewNodes[i]=newNode;
3679 // if current elem is quadratic and current node is not medium
3680 // we have to check - may be it is needed to insert additional nodes
3681 if( elem->IsQuadratic() && !elem->IsMediumNode(node) ) {
3682 list< const SMDS_MeshNode* > & listNewNodes = nIt->second;
3683 if(listNewNodes.size()==nbsteps) {
3684 listNewNodes.clear();
3685 double coord[] = { node->X(), node->Y(), node->Z() };
3686 for ( int i = 0; i < nbsteps; i++ ) {
3687 double x = coord[0] + theParams.myDir.X()*theParams.mySteps->Value(i+1);
3688 double y = coord[1] + theParams.myDir.Y()*theParams.mySteps->Value(i+1);
3689 double z = coord[2] + theParams.myDir.Z()*theParams.mySteps->Value(i+1);
3690 if( theFlags & EXTRUSION_FLAG_SEW ) {
3691 const SMDS_MeshNode * newNode = CreateNode(x, y, z,
3692 theTolerance, theParams.myNodes);
3693 listNewNodes.push_back( newNode );
3696 const SMDS_MeshNode * newNode = aMesh->AddNode(x, y, z);
3697 myLastCreatedNodes.Append(newNode);
3698 srcNodes.Append( node );
3699 listNewNodes.push_back( newNode );
3701 coord[0] = coord[0] + theParams.myDir.X()*theParams.mySteps->Value(i+1);
3702 coord[1] = coord[1] + theParams.myDir.Y()*theParams.mySteps->Value(i+1);
3703 coord[2] = coord[2] + theParams.myDir.Z()*theParams.mySteps->Value(i+1);
3704 if( theFlags & EXTRUSION_FLAG_SEW ) {
3705 const SMDS_MeshNode * newNode = CreateNode(coord[0], coord[1], coord[2],
3706 theTolerance, theParams.myNodes);
3707 listNewNodes.push_back( newNode );
3710 const SMDS_MeshNode * newNode = aMesh->AddNode( coord[0], coord[1], coord[2] );
3711 myLastCreatedNodes.Append(newNode);
3712 srcNodes.Append( node );
3713 listNewNodes.push_back( newNode );
3719 newNodesItVec.push_back( nIt );
3721 // make new elements
3722 sweepElement( elem, newNodesItVec, newElemsMap[elem], nbsteps, srcElems );
3725 if( theFlags & EXTRUSION_FLAG_BOUNDARY ) {
3726 makeWalls( mapNewNodes, newElemsMap, mapElemNewNodes, theElems, nbsteps, srcElems );
3728 PGroupIDs newGroupIDs;
3729 if ( theMakeGroups )
3730 newGroupIDs = generateGroups( srcNodes, srcElems, "extruded");
3736 //=======================================================================
3737 //class : SMESH_MeshEditor_PathPoint
3738 //purpose : auxiliary class
3739 //=======================================================================
3740 class SMESH_MeshEditor_PathPoint {
3742 SMESH_MeshEditor_PathPoint() {
3743 myPnt.SetCoord(99., 99., 99.);
3744 myTgt.SetCoord(1.,0.,0.);
3748 void SetPnt(const gp_Pnt& aP3D){
3751 void SetTangent(const gp_Dir& aTgt){
3754 void SetAngle(const double& aBeta){
3757 void SetParameter(const double& aPrm){
3760 const gp_Pnt& Pnt()const{
3763 const gp_Dir& Tangent()const{
3766 double Angle()const{
3769 double Parameter()const{
3780 //=======================================================================
3781 //function : ExtrusionAlongTrack
3783 //=======================================================================
3784 SMESH_MeshEditor::Extrusion_Error
3785 SMESH_MeshEditor::ExtrusionAlongTrack (TIDSortedElemSet & theElements,
3786 SMESH_subMesh* theTrack,
3787 const SMDS_MeshNode* theN1,
3788 const bool theHasAngles,
3789 list<double>& theAngles,
3790 const bool theHasRefPoint,
3791 const gp_Pnt& theRefPoint,
3792 const bool theMakeGroups)
3794 myLastCreatedElems.Clear();
3795 myLastCreatedNodes.Clear();
3797 // source elements for each generated one
3798 SMESH_SequenceOfElemPtr srcElems, srcNodes;
3800 int j, aNbTP, aNbE, aNb;
3801 double aT1, aT2, aT, aAngle, aX, aY, aZ;
3802 std::list<double> aPrms;
3803 std::list<double>::iterator aItD;
3804 TIDSortedElemSet::iterator itElem;
3806 Standard_Real aTx1, aTx2, aL2, aTolVec, aTolVec2;
3810 Handle(Geom_Curve) aC3D;
3811 TopoDS_Edge aTrackEdge;
3812 TopoDS_Vertex aV1, aV2;
3814 SMDS_ElemIteratorPtr aItE;
3815 SMDS_NodeIteratorPtr aItN;
3816 SMDSAbs_ElementType aTypeE;
3818 TNodeOfNodeListMap mapNewNodes;
3819 TElemOfVecOfNnlmiMap mapElemNewNodes;
3820 TElemOfElemListMap newElemsMap;
3823 aTolVec2=aTolVec*aTolVec;
3826 aNbE = theElements.size();
3829 return EXTR_NO_ELEMENTS;
3831 // 1.1 Track Pattern
3834 SMESHDS_SubMesh* pSubMeshDS=theTrack->GetSubMeshDS();
3836 aItE = pSubMeshDS->GetElements();
3837 while ( aItE->more() ) {
3838 const SMDS_MeshElement* pE = aItE->next();
3839 aTypeE = pE->GetType();
3840 // Pattern must contain links only
3841 if ( aTypeE != SMDSAbs_Edge )
3842 return EXTR_PATH_NOT_EDGE;
3845 const TopoDS_Shape& aS = theTrack->GetSubShape();
3846 // Sub shape for the Pattern must be an Edge
3847 if ( aS.ShapeType() != TopAbs_EDGE )
3848 return EXTR_BAD_PATH_SHAPE;
3850 aTrackEdge = TopoDS::Edge( aS );
3851 // the Edge must not be degenerated
3852 if ( BRep_Tool::Degenerated( aTrackEdge ) )
3853 return EXTR_BAD_PATH_SHAPE;
3855 TopExp::Vertices( aTrackEdge, aV1, aV2 );
3856 aT1=BRep_Tool::Parameter( aV1, aTrackEdge );
3857 aT2=BRep_Tool::Parameter( aV2, aTrackEdge );
3859 aItN = theTrack->GetFather()->GetSubMesh( aV1 )->GetSubMeshDS()->GetNodes();
3860 const SMDS_MeshNode* aN1 = aItN->next();
3862 aItN = theTrack->GetFather()->GetSubMesh( aV2 )->GetSubMeshDS()->GetNodes();
3863 const SMDS_MeshNode* aN2 = aItN->next();
3865 // starting node must be aN1 or aN2
3866 if ( !( aN1 == theN1 || aN2 == theN1 ) )
3867 return EXTR_BAD_STARTING_NODE;
3869 aNbTP = pSubMeshDS->NbNodes() + 2;
3872 vector<double> aAngles( aNbTP );
3874 for ( j=0; j < aNbTP; ++j ) {
3878 if ( theHasAngles ) {
3879 aItD = theAngles.begin();
3880 for ( j=1; (aItD != theAngles.end()) && (j<aNbTP); ++aItD, ++j ) {
3882 aAngles[j] = aAngle;
3886 // 2. Collect parameters on the track edge
3887 aPrms.push_back( aT1 );
3888 aPrms.push_back( aT2 );
3890 aItN = pSubMeshDS->GetNodes();
3891 while ( aItN->more() ) {
3892 const SMDS_MeshNode* pNode = aItN->next();
3893 const SMDS_EdgePosition* pEPos =
3894 static_cast<const SMDS_EdgePosition*>( pNode->GetPosition().get() );
3895 aT = pEPos->GetUParameter();
3896 aPrms.push_back( aT );
3901 if ( aN1 == theN1 ) {
3913 SMESH_MeshEditor_PathPoint aPP;
3914 vector<SMESH_MeshEditor_PathPoint> aPPs( aNbTP );
3916 aC3D = BRep_Tool::Curve( aTrackEdge, aTx1, aTx2 );
3918 aItD = aPrms.begin();
3919 for ( j=0; aItD != aPrms.end(); ++aItD, ++j ) {
3921 aC3D->D1( aT, aP3D, aVec );
3922 aL2 = aVec.SquareMagnitude();
3923 if ( aL2 < aTolVec2 )
3924 return EXTR_CANT_GET_TANGENT;
3926 gp_Dir aTgt( aVec );
3927 aAngle = aAngles[j];
3930 aPP.SetTangent( aTgt );
3931 aPP.SetAngle( aAngle );
3932 aPP.SetParameter( aT );
3936 // 3. Center of rotation aV0
3938 if ( !theHasRefPoint ) {
3940 aGC.SetCoord( 0.,0.,0. );
3942 itElem = theElements.begin();
3943 for ( ; itElem != theElements.end(); itElem++ ) {
3944 const SMDS_MeshElement* elem = *itElem;
3946 SMDS_ElemIteratorPtr itN = elem->nodesIterator();
3947 while ( itN->more() ) {
3948 const SMDS_MeshNode* node = static_cast<const SMDS_MeshNode*>( itN->next() );
3953 if ( mapNewNodes.find( node ) == mapNewNodes.end() ) {
3954 list<const SMDS_MeshNode*> aLNx;
3955 mapNewNodes[node] = aLNx;
3957 gp_XYZ aXYZ( aX, aY, aZ );
3965 } // if (!theHasRefPoint) {
3966 mapNewNodes.clear();
3968 // 4. Processing the elements
3969 SMESHDS_Mesh* aMesh = GetMeshDS();
3971 for ( itElem = theElements.begin(); itElem != theElements.end(); itElem++ ) {
3972 // check element type
3973 const SMDS_MeshElement* elem = *itElem;
3974 aTypeE = elem->GetType();
3975 if ( !elem || ( aTypeE != SMDSAbs_Face && aTypeE != SMDSAbs_Edge ) )
3978 vector<TNodeOfNodeListMapItr> & newNodesItVec = mapElemNewNodes[ elem ];
3979 newNodesItVec.reserve( elem->NbNodes() );
3981 // loop on elem nodes
3983 SMDS_ElemIteratorPtr itN = elem->nodesIterator();
3984 while ( itN->more() )
3987 // check if a node has been already processed
3988 const SMDS_MeshNode* node =
3989 static_cast<const SMDS_MeshNode*>( itN->next() );
3990 TNodeOfNodeListMap::iterator nIt = mapNewNodes.find( node );
3991 if ( nIt == mapNewNodes.end() ) {
3992 nIt = mapNewNodes.insert( make_pair( node, list<const SMDS_MeshNode*>() )).first;
3993 list<const SMDS_MeshNode*>& listNewNodes = nIt->second;
3996 aX = node->X(); aY = node->Y(); aZ = node->Z();
3998 Standard_Real aAngle1x, aAngleT1T0, aTolAng;
3999 gp_Pnt aP0x, aP1x, aPN0, aPN1, aV0x, aV1x;
4000 gp_Ax1 anAx1, anAxT1T0;
4001 gp_Dir aDT1x, aDT0x, aDT1T0;
4006 aPN0.SetCoord(aX, aY, aZ);
4008 const SMESH_MeshEditor_PathPoint& aPP0 = aPPs[0];
4010 aDT0x= aPP0.Tangent();
4012 for ( j = 1; j < aNbTP; ++j ) {
4013 const SMESH_MeshEditor_PathPoint& aPP1 = aPPs[j];
4015 aDT1x = aPP1.Tangent();
4016 aAngle1x = aPP1.Angle();
4018 gp_Trsf aTrsf, aTrsfRot, aTrsfRotT1T0;
4020 gp_Vec aV01x( aP0x, aP1x );
4021 aTrsf.SetTranslation( aV01x );
4024 aV1x = aV0x.Transformed( aTrsf );
4025 aPN1 = aPN0.Transformed( aTrsf );
4027 // rotation 1 [ T1,T0 ]
4028 aAngleT1T0=-aDT1x.Angle( aDT0x );
4029 if (fabs(aAngleT1T0) > aTolAng) {
4031 anAxT1T0.SetLocation( aV1x );
4032 anAxT1T0.SetDirection( aDT1T0 );
4033 aTrsfRotT1T0.SetRotation( anAxT1T0, aAngleT1T0 );
4035 aPN1 = aPN1.Transformed( aTrsfRotT1T0 );
4039 if ( theHasAngles ) {
4040 anAx1.SetLocation( aV1x );
4041 anAx1.SetDirection( aDT1x );
4042 aTrsfRot.SetRotation( anAx1, aAngle1x );
4044 aPN1 = aPN1.Transformed( aTrsfRot );
4048 if( elem->IsQuadratic() && !elem->IsMediumNode(node) ) {
4049 // create additional node
4050 double x = ( aPN1.X() + aPN0.X() )/2.;
4051 double y = ( aPN1.Y() + aPN0.Y() )/2.;
4052 double z = ( aPN1.Z() + aPN0.Z() )/2.;
4053 const SMDS_MeshNode* newNode = aMesh->AddNode(x,y,z);
4054 myLastCreatedNodes.Append(newNode);
4055 srcNodes.Append( node );
4056 listNewNodes.push_back( newNode );
4061 const SMDS_MeshNode* newNode = aMesh->AddNode( aX, aY, aZ );
4062 myLastCreatedNodes.Append(newNode);
4063 srcNodes.Append( node );
4064 listNewNodes.push_back( newNode );
4074 // if current elem is quadratic and current node is not medium
4075 // we have to check - may be it is needed to insert additional nodes
4076 if( elem->IsQuadratic() && !elem->IsMediumNode(node) ) {
4077 list< const SMDS_MeshNode* > & listNewNodes = nIt->second;
4078 if(listNewNodes.size()==aNbTP-1) {
4079 vector<const SMDS_MeshNode*> aNodes(2*(aNbTP-1));
4080 gp_XYZ P(node->X(), node->Y(), node->Z());
4081 list< const SMDS_MeshNode* >::iterator it = listNewNodes.begin();
4083 for(i=0; i<aNbTP-1; i++) {
4084 const SMDS_MeshNode* N = *it;
4085 double x = ( N->X() + P.X() )/2.;
4086 double y = ( N->Y() + P.Y() )/2.;
4087 double z = ( N->Z() + P.Z() )/2.;
4088 const SMDS_MeshNode* newN = aMesh->AddNode(x,y,z);
4089 srcNodes.Append( node );
4090 myLastCreatedNodes.Append(newN);
4093 P = gp_XYZ(N->X(),N->Y(),N->Z());
4095 listNewNodes.clear();
4096 for(i=0; i<2*(aNbTP-1); i++) {
4097 listNewNodes.push_back(aNodes[i]);
4103 newNodesItVec.push_back( nIt );
4105 // make new elements
4106 //sweepElement( aMesh, elem, newNodesItVec, newElemsMap[elem],
4107 // newNodesItVec[0]->second.size(), myLastCreatedElems );
4108 sweepElement( elem, newNodesItVec, newElemsMap[elem], aNbTP-1, srcElems );
4111 makeWalls( mapNewNodes, newElemsMap, mapElemNewNodes, theElements, aNbTP-1, srcElems );
4113 if ( theMakeGroups )
4114 generateGroups( srcNodes, srcElems, "extruded");
4119 //=======================================================================
4120 //function : Transform
4122 //=======================================================================
4124 SMESH_MeshEditor::PGroupIDs
4125 SMESH_MeshEditor::Transform (TIDSortedElemSet & theElems,
4126 const gp_Trsf& theTrsf,
4128 const bool theMakeGroups,
4129 SMESH_Mesh* theTargetMesh)
4131 myLastCreatedElems.Clear();
4132 myLastCreatedNodes.Clear();
4134 bool needReverse = false;
4135 string groupPostfix;
4136 switch ( theTrsf.Form() ) {
4141 groupPostfix = "mirrored";
4144 groupPostfix = "rotated";
4146 case gp_Translation:
4147 groupPostfix = "translated";
4150 groupPostfix = "scaled";
4153 needReverse = false;
4154 groupPostfix = "transformed";
4157 SMESH_MeshEditor targetMeshEditor( theTargetMesh );
4158 SMESHDS_Mesh* aTgtMesh = theTargetMesh ? theTargetMesh->GetMeshDS() : 0;
4159 SMESHDS_Mesh* aMesh = GetMeshDS();
4162 // map old node to new one
4163 TNodeNodeMap nodeMap;
4165 // elements sharing moved nodes; those of them which have all
4166 // nodes mirrored but are not in theElems are to be reversed
4167 TIDSortedElemSet inverseElemSet;
4169 // source elements for each generated one
4170 SMESH_SequenceOfElemPtr srcElems, srcNodes;
4173 TIDSortedElemSet::iterator itElem;
4174 for ( itElem = theElems.begin(); itElem != theElems.end(); itElem++ ) {
4175 const SMDS_MeshElement* elem = *itElem;
4179 // loop on elem nodes
4180 SMDS_ElemIteratorPtr itN = elem->nodesIterator();
4181 while ( itN->more() ) {
4183 // check if a node has been already transformed
4184 const SMDS_MeshNode* node = cast2Node( itN->next() );
4185 pair<TNodeNodeMap::iterator,bool> n2n_isnew =
4186 nodeMap.insert( make_pair ( node, node ));
4187 if ( !n2n_isnew.second )
4191 coord[0] = node->X();
4192 coord[1] = node->Y();
4193 coord[2] = node->Z();
4194 theTrsf.Transforms( coord[0], coord[1], coord[2] );
4195 if ( theTargetMesh ) {
4196 const SMDS_MeshNode * newNode = aTgtMesh->AddNode( coord[0], coord[1], coord[2] );
4197 n2n_isnew.first->second = newNode;
4198 myLastCreatedNodes.Append(newNode);
4199 srcNodes.Append( node );
4201 else if ( theCopy ) {
4202 const SMDS_MeshNode * newNode = aMesh->AddNode( coord[0], coord[1], coord[2] );
4203 n2n_isnew.first->second = newNode;
4204 myLastCreatedNodes.Append(newNode);
4205 srcNodes.Append( node );
4208 aMesh->MoveNode( node, coord[0], coord[1], coord[2] );
4209 // node position on shape becomes invalid
4210 const_cast< SMDS_MeshNode* > ( node )->SetPosition
4211 ( SMDS_SpacePosition::originSpacePosition() );
4214 // keep inverse elements
4215 if ( !theCopy && !theTargetMesh && needReverse ) {
4216 SMDS_ElemIteratorPtr invElemIt = node->GetInverseElementIterator();
4217 while ( invElemIt->more() ) {
4218 const SMDS_MeshElement* iel = invElemIt->next();
4219 inverseElemSet.insert( iel );
4225 // either create new elements or reverse mirrored ones
4226 if ( !theCopy && !needReverse && !theTargetMesh )
4229 TIDSortedElemSet::iterator invElemIt = inverseElemSet.begin();
4230 for ( ; invElemIt != inverseElemSet.end(); invElemIt++ )
4231 theElems.insert( *invElemIt );
4233 // replicate or reverse elements
4236 REV_TETRA = 0, // = nbNodes - 4
4237 REV_PYRAMID = 1, // = nbNodes - 4
4238 REV_PENTA = 2, // = nbNodes - 4
4240 REV_HEXA = 4, // = nbNodes - 4
4244 { 2, 1, 0, 3, 4, 0, 0, 0 }, // REV_TETRA
4245 { 2, 1, 0, 3, 4, 0, 0, 0 }, // REV_PYRAMID
4246 { 2, 1, 0, 5, 4, 3, 0, 0 }, // REV_PENTA
4247 { 2, 1, 0, 3, 0, 0, 0, 0 }, // REV_FACE
4248 { 2, 1, 0, 3, 6, 5, 4, 7 }, // REV_HEXA
4249 { 0, 1, 2, 3, 4, 5, 6, 7 } // FORWARD
4252 for ( itElem = theElems.begin(); itElem != theElems.end(); itElem++ )
4254 const SMDS_MeshElement* elem = *itElem;
4255 if ( !elem || elem->GetType() == SMDSAbs_Node )
4258 int nbNodes = elem->NbNodes();
4259 int elemType = elem->GetType();
4261 if (elem->IsPoly()) {
4262 // Polygon or Polyhedral Volume
4263 switch ( elemType ) {
4266 vector<const SMDS_MeshNode*> poly_nodes (nbNodes);
4268 SMDS_ElemIteratorPtr itN = elem->nodesIterator();
4269 while (itN->more()) {
4270 const SMDS_MeshNode* node =
4271 static_cast<const SMDS_MeshNode*>(itN->next());
4272 TNodeNodeMap::iterator nodeMapIt = nodeMap.find(node);
4273 if (nodeMapIt == nodeMap.end())
4274 break; // not all nodes transformed
4276 // reverse mirrored faces and volumes
4277 poly_nodes[nbNodes - iNode - 1] = (*nodeMapIt).second;
4279 poly_nodes[iNode] = (*nodeMapIt).second;
4283 if ( iNode != nbNodes )
4284 continue; // not all nodes transformed
4286 if ( theTargetMesh ) {
4287 myLastCreatedElems.Append(aTgtMesh->AddPolygonalFace(poly_nodes));
4288 srcElems.Append( elem );
4290 else if ( theCopy ) {
4291 myLastCreatedElems.Append(aMesh->AddPolygonalFace(poly_nodes));
4292 srcElems.Append( elem );
4295 aMesh->ChangePolygonNodes(elem, poly_nodes);
4299 case SMDSAbs_Volume:
4301 // ATTENTION: Reversing is not yet done!!!
4302 const SMDS_PolyhedralVolumeOfNodes* aPolyedre =
4303 dynamic_cast<const SMDS_PolyhedralVolumeOfNodes*>( elem );
4305 MESSAGE("Warning: bad volumic element");
4309 vector<const SMDS_MeshNode*> poly_nodes;
4310 vector<int> quantities;
4312 bool allTransformed = true;
4313 int nbFaces = aPolyedre->NbFaces();
4314 for (int iface = 1; iface <= nbFaces && allTransformed; iface++) {
4315 int nbFaceNodes = aPolyedre->NbFaceNodes(iface);
4316 for (int inode = 1; inode <= nbFaceNodes && allTransformed; inode++) {
4317 const SMDS_MeshNode* node = aPolyedre->GetFaceNode(iface, inode);
4318 TNodeNodeMap::iterator nodeMapIt = nodeMap.find(node);
4319 if (nodeMapIt == nodeMap.end()) {
4320 allTransformed = false; // not all nodes transformed
4322 poly_nodes.push_back((*nodeMapIt).second);
4325 quantities.push_back(nbFaceNodes);
4327 if ( !allTransformed )
4328 continue; // not all nodes transformed
4330 if ( theTargetMesh ) {
4331 myLastCreatedElems.Append(aTgtMesh->AddPolyhedralVolume(poly_nodes, quantities));
4332 srcElems.Append( elem );
4334 else if ( theCopy ) {
4335 myLastCreatedElems.Append(aMesh->AddPolyhedralVolume(poly_nodes, quantities));
4336 srcElems.Append( elem );
4339 aMesh->ChangePolyhedronNodes(elem, poly_nodes, quantities);
4349 int* i = index[ FORWARD ];
4350 if ( needReverse && nbNodes > 2) // reverse mirrored faces and volumes
4351 if ( elemType == SMDSAbs_Face )
4352 i = index[ REV_FACE ];
4354 i = index[ nbNodes - 4 ];
4356 if(elem->IsQuadratic()) {
4357 static int anIds[] = {0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19};
4360 if(nbNodes==3) { // quadratic edge
4361 static int anIds[] = {1,0,2};
4364 else if(nbNodes==6) { // quadratic triangle
4365 static int anIds[] = {0,2,1,5,4,3};
4368 else if(nbNodes==8) { // quadratic quadrangle
4369 static int anIds[] = {0,3,2,1,7,6,5,4};
4372 else if(nbNodes==10) { // quadratic tetrahedron of 10 nodes
4373 static int anIds[] = {0,2,1,3,6,5,4,7,9,8};
4376 else if(nbNodes==13) { // quadratic pyramid of 13 nodes
4377 static int anIds[] = {0,3,2,1,4,8,7,6,5,9,12,11,10};
4380 else if(nbNodes==15) { // quadratic pentahedron with 15 nodes
4381 static int anIds[] = {0,2,1,3,5,4,8,7,6,11,10,9,12,14,13};
4384 else { // nbNodes==20 - quadratic hexahedron with 20 nodes
4385 static int anIds[] = {0,3,2,1,4,7,6,5,11,10,9,8,15,14,13,12,16,19,18,17};
4391 // find transformed nodes
4392 vector<const SMDS_MeshNode*> nodes(nbNodes);
4394 SMDS_ElemIteratorPtr itN = elem->nodesIterator();
4395 while ( itN->more() ) {
4396 const SMDS_MeshNode* node =
4397 static_cast<const SMDS_MeshNode*>( itN->next() );
4398 TNodeNodeMap::iterator nodeMapIt = nodeMap.find( node );
4399 if ( nodeMapIt == nodeMap.end() )
4400 break; // not all nodes transformed
4401 nodes[ i [ iNode++ ]] = (*nodeMapIt).second;
4403 if ( iNode != nbNodes )
4404 continue; // not all nodes transformed
4406 if ( theTargetMesh ) {
4407 if ( SMDS_MeshElement* copy =
4408 targetMeshEditor.AddElement( nodes, elem->GetType(), elem->IsPoly() )) {
4409 myLastCreatedElems.Append( copy );
4410 srcElems.Append( elem );
4413 else if ( theCopy ) {
4414 if ( SMDS_MeshElement* copy = AddElement( nodes, elem->GetType(), elem->IsPoly() )) {
4415 myLastCreatedElems.Append( copy );
4416 srcElems.Append( elem );
4420 // reverse element as it was reversed by transformation
4422 aMesh->ChangeElementNodes( elem, &nodes[0], nbNodes );
4426 PGroupIDs newGroupIDs;
4428 if ( theMakeGroups && theCopy ||
4429 theMakeGroups && theTargetMesh )
4430 newGroupIDs = generateGroups( srcNodes, srcElems, groupPostfix, theTargetMesh );
4435 //=======================================================================
4437 * \brief Create groups of elements made during transformation
4438 * \param nodeGens - nodes making corresponding myLastCreatedNodes
4439 * \param elemGens - elements making corresponding myLastCreatedElems
4440 * \param postfix - to append to names of new groups
4442 //=======================================================================
4444 SMESH_MeshEditor::PGroupIDs
4445 SMESH_MeshEditor::generateGroups(const SMESH_SequenceOfElemPtr& nodeGens,
4446 const SMESH_SequenceOfElemPtr& elemGens,
4447 const std::string& postfix,
4448 SMESH_Mesh* targetMesh)
4450 PGroupIDs newGroupIDs( new list<int> );
4451 SMESH_Mesh* mesh = targetMesh ? targetMesh : GetMesh();
4453 // Sort existing groups by types and collect their names
4455 // to store an old group and a generated new one
4456 typedef pair< SMESHDS_GroupBase*, SMDS_MeshGroup* > TOldNewGroup;
4457 vector< list< TOldNewGroup > > groupsByType( SMDSAbs_NbElementTypes );
4459 set< string > groupNames;
4461 SMDS_MeshGroup* nullNewGroup = (SMDS_MeshGroup*) 0;
4462 SMESH_Mesh::GroupIteratorPtr groupIt = GetMesh()->GetGroups();
4463 while ( groupIt->more() ) {
4464 SMESH_Group * group = groupIt->next();
4465 if ( !group ) continue;
4466 SMESHDS_GroupBase* groupDS = group->GetGroupDS();
4467 if ( !groupDS || groupDS->IsEmpty() ) continue;
4468 groupNames.insert( group->GetName() );
4469 groupDS->SetStoreName( group->GetName() );
4470 groupsByType[ groupDS->GetType() ].push_back( make_pair( groupDS, nullNewGroup ));
4475 // loop on nodes and elements
4476 for ( int isNodes = 0; isNodes < 2; ++isNodes )
4478 const SMESH_SequenceOfElemPtr& gens = isNodes ? nodeGens : elemGens;
4479 const SMESH_SequenceOfElemPtr& elems = isNodes ? myLastCreatedNodes : myLastCreatedElems;
4480 if ( gens.Length() != elems.Length() )
4481 throw SALOME_Exception(LOCALIZED("invalid args"));
4483 // loop on created elements
4484 for (int iElem = 1; iElem <= elems.Length(); ++iElem )
4486 const SMDS_MeshElement* sourceElem = gens( iElem );
4487 if ( !sourceElem ) {
4488 MESSAGE("generateGroups(): NULL source element");
4491 list< TOldNewGroup > & groupsOldNew = groupsByType[ sourceElem->GetType() ];
4492 if ( groupsOldNew.empty() ) {
4493 while ( iElem < gens.Length() && gens( iElem+1 ) == sourceElem )
4494 ++iElem; // skip all elements made by sourceElem
4497 // collect all elements made by sourceElem
4498 list< const SMDS_MeshElement* > resultElems;
4499 if ( const SMDS_MeshElement* resElem = elems( iElem ))
4500 if ( resElem != sourceElem )
4501 resultElems.push_back( resElem );
4502 while ( iElem < gens.Length() && gens( iElem+1 ) == sourceElem )
4503 if ( const SMDS_MeshElement* resElem = elems( ++iElem ))
4504 if ( resElem != sourceElem )
4505 resultElems.push_back( resElem );
4506 // do not generate element groups from node ones
4507 if ( sourceElem->GetType() == SMDSAbs_Node &&
4508 elems( iElem )->GetType() != SMDSAbs_Node )
4511 // add resultElems to groups made by ones the sourceElem belongs to
4512 list< TOldNewGroup >::iterator gOldNew, gLast = groupsOldNew.end();
4513 for ( gOldNew = groupsOldNew.begin(); gOldNew != gLast; ++gOldNew )
4515 SMESHDS_GroupBase* oldGroup = gOldNew->first;
4516 if ( oldGroup->Contains( sourceElem )) // sourceElem in oldGroup
4518 SMDS_MeshGroup* & newGroup = gOldNew->second;
4519 if ( !newGroup )// create a new group
4522 string name = oldGroup->GetStoreName();
4523 if ( !targetMesh ) {
4527 while ( !groupNames.insert( name ).second ) // name exists
4533 TCollection_AsciiString nbStr(nb+1);
4534 name.resize( name.rfind('_')+1 );
4535 name += nbStr.ToCString();
4542 SMESH_Group* group = mesh->AddGroup( resultElems.back()->GetType(),
4544 SMESHDS_Group* groupDS = static_cast<SMESHDS_Group*>(group->GetGroupDS());
4545 newGroup = & groupDS->SMDSGroup();
4546 newGroupIDs->push_back( id );
4549 // fill in a new group
4550 list< const SMDS_MeshElement* >::iterator resLast = resultElems.end(), resElemIt;
4551 for ( resElemIt = resultElems.begin(); resElemIt != resLast; ++resElemIt )
4552 newGroup->Add( *resElemIt );
4555 } // loop on created elements
4556 }// loop on nodes and elements
4561 //=======================================================================
4562 //function : FindCoincidentNodes
4563 //purpose : Return list of group of nodes close to each other within theTolerance
4564 // Search among theNodes or in the whole mesh if theNodes is empty using
4565 // an Octree algorithm
4566 //=======================================================================
4568 void SMESH_MeshEditor::FindCoincidentNodes (set<const SMDS_MeshNode*> & theNodes,
4569 const double theTolerance,
4570 TListOfListOfNodes & theGroupsOfNodes)
4572 myLastCreatedElems.Clear();
4573 myLastCreatedNodes.Clear();
4575 set<const SMDS_MeshNode*> nodes;
4576 if ( theNodes.empty() )
4577 { // get all nodes in the mesh
4578 SMDS_NodeIteratorPtr nIt = GetMeshDS()->nodesIterator();
4579 while ( nIt->more() )
4580 nodes.insert( nodes.end(),nIt->next());
4584 SMESH_OctreeNode::FindCoincidentNodes ( nodes, &theGroupsOfNodes, theTolerance);
4588 //=======================================================================
4590 * \brief Implementation of search for the node closest to point
4592 //=======================================================================
4594 struct SMESH_NodeSearcherImpl: public SMESH_NodeSearcher
4597 * \brief Constructor
4599 SMESH_NodeSearcherImpl( const SMESHDS_Mesh* theMesh )
4601 set<const SMDS_MeshNode*> nodes;
4603 SMDS_NodeIteratorPtr nIt = theMesh->nodesIterator();
4604 while ( nIt->more() )
4605 nodes.insert( nodes.end(), nIt->next() );
4607 myOctreeNode = new SMESH_OctreeNode(nodes) ;
4610 * \brief Do it's job
4612 const SMDS_MeshNode* FindClosestTo( const gp_Pnt& thePnt )
4614 SMDS_MeshNode tgtNode( thePnt.X(), thePnt.Y(), thePnt.Z() );
4615 list<const SMDS_MeshNode*> nodes;
4616 const double precision = 1e-6;
4617 //myOctreeNode->NodesAround( &tgtNode, &nodes, precision );
4619 double minSqDist = DBL_MAX;
4621 if ( nodes.empty() ) // get all nodes of OctreeNode's closest to thePnt
4623 // sort leafs by their distance from thePnt
4624 typedef map< double, SMESH_OctreeNode* > TDistTreeMap;
4625 TDistTreeMap treeMap;
4626 list< SMESH_OctreeNode* > treeList;
4627 list< SMESH_OctreeNode* >::iterator trIt;
4628 treeList.push_back( myOctreeNode );
4629 for ( trIt = treeList.begin(); trIt != treeList.end(); ++trIt)
4631 SMESH_OctreeNode* tree = *trIt;
4632 if ( !tree->isLeaf() ) { // put children to the queue
4633 SMESH_OctreeNodeIteratorPtr cIt = tree->GetChildrenIterator();
4634 while ( cIt->more() )
4635 treeList.push_back( cIt->next() );
4637 else if ( tree->NbNodes() ) { // put tree to treeMap
4638 tree->getBox( box );
4639 double sqDist = thePnt.SquareDistance( 0.5 * ( box.CornerMin() + box.CornerMax() ));
4640 pair<TDistTreeMap::iterator,bool> it_in = treeMap.insert( make_pair( sqDist, tree ));
4641 if ( !it_in.second ) // not unique distance to box center
4642 treeMap.insert( it_in.first, make_pair( sqDist - 1e-13*treeMap.size(), tree ));
4645 // find distance after which there is no sense to check tree's
4646 double sqLimit = DBL_MAX;
4647 TDistTreeMap::iterator sqDist_tree = treeMap.begin();
4648 if ( treeMap.size() > 5 ) {
4649 SMESH_OctreeNode* closestTree = sqDist_tree->second;
4650 closestTree->getBox( box );
4651 double limit = sqrt( sqDist_tree->first ) + sqrt ( box.SquareExtent() );
4652 sqLimit = limit * limit;
4654 // get all nodes from trees
4655 for ( ; sqDist_tree != treeMap.end(); ++sqDist_tree) {
4656 if ( sqDist_tree->first > sqLimit )
4658 SMESH_OctreeNode* tree = sqDist_tree->second;
4659 tree->NodesAround( tree->GetNodeIterator()->next(), &nodes );
4662 // find closest among nodes
4663 minSqDist = DBL_MAX;
4664 const SMDS_MeshNode* closestNode = 0;
4665 list<const SMDS_MeshNode*>::iterator nIt = nodes.begin();
4666 for ( ; nIt != nodes.end(); ++nIt ) {
4667 double sqDist = thePnt.SquareDistance( TNodeXYZ( *nIt ) );
4668 if ( minSqDist > sqDist ) {
4678 ~SMESH_NodeSearcherImpl() { delete myOctreeNode; }
4680 SMESH_OctreeNode* myOctreeNode;
4683 //=======================================================================
4685 * \brief Return SMESH_NodeSearcher
4687 //=======================================================================
4689 SMESH_NodeSearcher* SMESH_MeshEditor::GetNodeSearcher()
4691 return new SMESH_NodeSearcherImpl( GetMeshDS() );
4694 //=======================================================================
4695 //function : SimplifyFace
4697 //=======================================================================
4698 int SMESH_MeshEditor::SimplifyFace (const vector<const SMDS_MeshNode *> faceNodes,
4699 vector<const SMDS_MeshNode *>& poly_nodes,
4700 vector<int>& quantities) const
4702 int nbNodes = faceNodes.size();
4707 set<const SMDS_MeshNode*> nodeSet;
4709 // get simple seq of nodes
4710 //const SMDS_MeshNode* simpleNodes[ nbNodes ];
4711 vector<const SMDS_MeshNode*> simpleNodes( nbNodes );
4712 int iSimple = 0, nbUnique = 0;
4714 simpleNodes[iSimple++] = faceNodes[0];
4716 for (int iCur = 1; iCur < nbNodes; iCur++) {
4717 if (faceNodes[iCur] != simpleNodes[iSimple - 1]) {
4718 simpleNodes[iSimple++] = faceNodes[iCur];
4719 if (nodeSet.insert( faceNodes[iCur] ).second)
4723 int nbSimple = iSimple;
4724 if (simpleNodes[nbSimple - 1] == simpleNodes[0]) {
4734 bool foundLoop = (nbSimple > nbUnique);
4737 set<const SMDS_MeshNode*> loopSet;
4738 for (iSimple = 0; iSimple < nbSimple && !foundLoop; iSimple++) {
4739 const SMDS_MeshNode* n = simpleNodes[iSimple];
4740 if (!loopSet.insert( n ).second) {
4744 int iC = 0, curLast = iSimple;
4745 for (; iC < curLast; iC++) {
4746 if (simpleNodes[iC] == n) break;
4748 int loopLen = curLast - iC;
4750 // create sub-element
4752 quantities.push_back(loopLen);
4753 for (; iC < curLast; iC++) {
4754 poly_nodes.push_back(simpleNodes[iC]);
4757 // shift the rest nodes (place from the first loop position)
4758 for (iC = curLast + 1; iC < nbSimple; iC++) {
4759 simpleNodes[iC - loopLen] = simpleNodes[iC];
4761 nbSimple -= loopLen;
4764 } // for (iSimple = 0; iSimple < nbSimple; iSimple++)
4765 } // while (foundLoop)
4769 quantities.push_back(iSimple);
4770 for (int i = 0; i < iSimple; i++)
4771 poly_nodes.push_back(simpleNodes[i]);
4777 //=======================================================================
4778 //function : MergeNodes
4779 //purpose : In each group, the cdr of nodes are substituted by the first one
4781 //=======================================================================
4783 void SMESH_MeshEditor::MergeNodes (TListOfListOfNodes & theGroupsOfNodes)
4785 myLastCreatedElems.Clear();
4786 myLastCreatedNodes.Clear();
4788 SMESHDS_Mesh* aMesh = GetMeshDS();
4790 TNodeNodeMap nodeNodeMap; // node to replace - new node
4791 set<const SMDS_MeshElement*> elems; // all elements with changed nodes
4792 list< int > rmElemIds, rmNodeIds;
4794 // Fill nodeNodeMap and elems
4796 TListOfListOfNodes::iterator grIt = theGroupsOfNodes.begin();
4797 for ( ; grIt != theGroupsOfNodes.end(); grIt++ ) {
4798 list<const SMDS_MeshNode*>& nodes = *grIt;
4799 list<const SMDS_MeshNode*>::iterator nIt = nodes.begin();
4800 const SMDS_MeshNode* nToKeep = *nIt;
4801 for ( ++nIt; nIt != nodes.end(); nIt++ ) {
4802 const SMDS_MeshNode* nToRemove = *nIt;
4803 nodeNodeMap.insert( TNodeNodeMap::value_type( nToRemove, nToKeep ));
4804 if ( nToRemove != nToKeep ) {
4805 rmNodeIds.push_back( nToRemove->GetID() );
4806 AddToSameGroups( nToKeep, nToRemove, aMesh );
4809 SMDS_ElemIteratorPtr invElemIt = nToRemove->GetInverseElementIterator();
4810 while ( invElemIt->more() ) {
4811 const SMDS_MeshElement* elem = invElemIt->next();
4816 // Change element nodes or remove an element
4818 set<const SMDS_MeshElement*>::iterator eIt = elems.begin();
4819 for ( ; eIt != elems.end(); eIt++ ) {
4820 const SMDS_MeshElement* elem = *eIt;
4821 int nbNodes = elem->NbNodes();
4822 int aShapeId = FindShape( elem );
4824 set<const SMDS_MeshNode*> nodeSet;
4825 vector< const SMDS_MeshNode*> curNodes( nbNodes ), uniqueNodes( nbNodes );
4826 int iUnique = 0, iCur = 0, nbRepl = 0;
4827 vector<int> iRepl( nbNodes );
4829 // get new seq of nodes
4830 SMDS_ElemIteratorPtr itN = elem->nodesIterator();
4831 while ( itN->more() ) {
4832 const SMDS_MeshNode* n =
4833 static_cast<const SMDS_MeshNode*>( itN->next() );
4835 TNodeNodeMap::iterator nnIt = nodeNodeMap.find( n );
4836 if ( nnIt != nodeNodeMap.end() ) { // n sticks
4838 // BUG 0020185: begin
4840 bool stopRecur = false;
4841 set<const SMDS_MeshNode*> nodesRecur;
4842 nodesRecur.insert(n);
4843 while (!stopRecur) {
4844 TNodeNodeMap::iterator nnIt_i = nodeNodeMap.find( n );
4845 if ( nnIt_i != nodeNodeMap.end() ) { // n sticks
4846 n = (*nnIt_i).second;
4847 if (!nodesRecur.insert(n).second) {
4848 // error: recursive dependancy
4857 iRepl[ nbRepl++ ] = iCur;
4859 curNodes[ iCur ] = n;
4860 bool isUnique = nodeSet.insert( n ).second;
4862 uniqueNodes[ iUnique++ ] = n;
4866 // Analyse element topology after replacement
4869 int nbUniqueNodes = nodeSet.size();
4870 if ( nbNodes != nbUniqueNodes ) { // some nodes stick
4871 // Polygons and Polyhedral volumes
4872 if (elem->IsPoly()) {
4874 if (elem->GetType() == SMDSAbs_Face) {
4876 vector<const SMDS_MeshNode *> face_nodes (nbNodes);
4878 for (; inode < nbNodes; inode++) {
4879 face_nodes[inode] = curNodes[inode];
4882 vector<const SMDS_MeshNode *> polygons_nodes;
4883 vector<int> quantities;
4884 int nbNew = SimplifyFace(face_nodes, polygons_nodes, quantities);
4888 for (int iface = 0; iface < nbNew - 1; iface++) {
4889 int nbNodes = quantities[iface];
4890 vector<const SMDS_MeshNode *> poly_nodes (nbNodes);
4891 for (int ii = 0; ii < nbNodes; ii++, inode++) {
4892 poly_nodes[ii] = polygons_nodes[inode];
4894 SMDS_MeshElement* newElem = aMesh->AddPolygonalFace(poly_nodes);
4895 myLastCreatedElems.Append(newElem);
4897 aMesh->SetMeshElementOnShape(newElem, aShapeId);
4899 aMesh->ChangeElementNodes(elem, &polygons_nodes[inode], quantities[nbNew - 1]);
4902 rmElemIds.push_back(elem->GetID());
4906 else if (elem->GetType() == SMDSAbs_Volume) {
4907 // Polyhedral volume
4908 if (nbUniqueNodes < 4) {
4909 rmElemIds.push_back(elem->GetID());
4912 // each face has to be analized in order to check volume validity
4913 const SMDS_PolyhedralVolumeOfNodes* aPolyedre =
4914 static_cast<const SMDS_PolyhedralVolumeOfNodes*>( elem );
4916 int nbFaces = aPolyedre->NbFaces();
4918 vector<const SMDS_MeshNode *> poly_nodes;
4919 vector<int> quantities;
4921 for (int iface = 1; iface <= nbFaces; iface++) {
4922 int nbFaceNodes = aPolyedre->NbFaceNodes(iface);
4923 vector<const SMDS_MeshNode *> faceNodes (nbFaceNodes);
4925 for (int inode = 1; inode <= nbFaceNodes; inode++) {
4926 const SMDS_MeshNode * faceNode = aPolyedre->GetFaceNode(iface, inode);
4927 TNodeNodeMap::iterator nnIt = nodeNodeMap.find(faceNode);
4928 if (nnIt != nodeNodeMap.end()) { // faceNode sticks
4929 faceNode = (*nnIt).second;
4931 faceNodes[inode - 1] = faceNode;
4934 SimplifyFace(faceNodes, poly_nodes, quantities);
4937 if (quantities.size() > 3) {
4938 // to be done: remove coincident faces
4941 if (quantities.size() > 3)
4942 aMesh->ChangePolyhedronNodes(elem, poly_nodes, quantities);
4944 rmElemIds.push_back(elem->GetID());
4948 rmElemIds.push_back(elem->GetID());
4959 switch ( nbNodes ) {
4960 case 2: ///////////////////////////////////// EDGE
4961 isOk = false; break;
4962 case 3: ///////////////////////////////////// TRIANGLE
4963 isOk = false; break;
4965 if ( elem->GetType() == SMDSAbs_Volume ) // TETRAHEDRON
4967 else { //////////////////////////////////// QUADRANGLE
4968 if ( nbUniqueNodes < 3 )
4970 else if ( nbRepl == 2 && iRepl[ 1 ] - iRepl[ 0 ] == 2 )
4971 isOk = false; // opposite nodes stick
4974 case 6: ///////////////////////////////////// PENTAHEDRON
4975 if ( nbUniqueNodes == 4 ) {
4976 // ---------------------------------> tetrahedron
4978 iRepl[ 0 ] > 2 && iRepl[ 1 ] > 2 && iRepl[ 2 ] > 2 ) {
4979 // all top nodes stick: reverse a bottom
4980 uniqueNodes[ 0 ] = curNodes [ 1 ];
4981 uniqueNodes[ 1 ] = curNodes [ 0 ];
4983 else if (nbRepl == 3 &&
4984 iRepl[ 0 ] < 3 && iRepl[ 1 ] < 3 && iRepl[ 2 ] < 3 ) {
4985 // all bottom nodes stick: set a top before
4986 uniqueNodes[ 3 ] = uniqueNodes [ 0 ];
4987 uniqueNodes[ 0 ] = curNodes [ 3 ];
4988 uniqueNodes[ 1 ] = curNodes [ 4 ];
4989 uniqueNodes[ 2 ] = curNodes [ 5 ];
4991 else if (nbRepl == 4 &&
4992 iRepl[ 2 ] - iRepl [ 0 ] == 3 && iRepl[ 3 ] - iRepl [ 1 ] == 3 ) {
4993 // a lateral face turns into a line: reverse a bottom
4994 uniqueNodes[ 0 ] = curNodes [ 1 ];
4995 uniqueNodes[ 1 ] = curNodes [ 0 ];
5000 else if ( nbUniqueNodes == 5 ) {
5001 // PENTAHEDRON --------------------> 2 tetrahedrons
5002 if ( nbRepl == 2 && iRepl[ 1 ] - iRepl [ 0 ] == 3 ) {
5003 // a bottom node sticks with a linked top one
5005 SMDS_MeshElement* newElem =
5006 aMesh->AddVolume(curNodes[ 3 ],
5009 curNodes[ iRepl[ 0 ] == 2 ? 1 : 2 ]);
5010 myLastCreatedElems.Append(newElem);
5012 aMesh->SetMeshElementOnShape( newElem, aShapeId );
5013 // 2. : reverse a bottom
5014 uniqueNodes[ 0 ] = curNodes [ 1 ];
5015 uniqueNodes[ 1 ] = curNodes [ 0 ];
5025 if(elem->IsQuadratic()) { // Quadratic quadrangle
5038 if( iRepl[0]==0 && iRepl[1]==1 && iRepl[2]==4 ) {
5039 uniqueNodes[0] = curNodes[0];
5040 uniqueNodes[1] = curNodes[2];
5041 uniqueNodes[2] = curNodes[3];
5042 uniqueNodes[3] = curNodes[5];
5043 uniqueNodes[4] = curNodes[6];
5044 uniqueNodes[5] = curNodes[7];
5047 if( iRepl[0]==0 && iRepl[1]==3 && iRepl[2]==7 ) {
5048 uniqueNodes[0] = curNodes[0];
5049 uniqueNodes[1] = curNodes[1];
5050 uniqueNodes[2] = curNodes[2];
5051 uniqueNodes[3] = curNodes[4];
5052 uniqueNodes[4] = curNodes[5];
5053 uniqueNodes[5] = curNodes[6];
5056 if( iRepl[0]==0 && iRepl[1]==4 && iRepl[2]==7 ) {
5057 uniqueNodes[0] = curNodes[1];
5058 uniqueNodes[1] = curNodes[2];
5059 uniqueNodes[2] = curNodes[3];
5060 uniqueNodes[3] = curNodes[5];
5061 uniqueNodes[4] = curNodes[6];
5062 uniqueNodes[5] = curNodes[0];
5065 if( iRepl[0]==1 && iRepl[1]==2 && iRepl[2]==5 ) {
5066 uniqueNodes[0] = curNodes[0];
5067 uniqueNodes[1] = curNodes[1];
5068 uniqueNodes[2] = curNodes[3];
5069 uniqueNodes[3] = curNodes[4];
5070 uniqueNodes[4] = curNodes[6];
5071 uniqueNodes[5] = curNodes[7];
5074 if( iRepl[0]==1 && iRepl[1]==4 && iRepl[2]==5 ) {
5075 uniqueNodes[0] = curNodes[0];
5076 uniqueNodes[1] = curNodes[2];
5077 uniqueNodes[2] = curNodes[3];
5078 uniqueNodes[3] = curNodes[1];
5079 uniqueNodes[4] = curNodes[6];
5080 uniqueNodes[5] = curNodes[7];
5083 if( iRepl[0]==2 && iRepl[1]==3 && iRepl[2]==6 ) {
5084 uniqueNodes[0] = curNodes[0];
5085 uniqueNodes[1] = curNodes[1];
5086 uniqueNodes[2] = curNodes[2];
5087 uniqueNodes[3] = curNodes[4];
5088 uniqueNodes[4] = curNodes[5];
5089 uniqueNodes[5] = curNodes[7];
5092 if( iRepl[0]==2 && iRepl[1]==5 && iRepl[2]==6 ) {
5093 uniqueNodes[0] = curNodes[0];
5094 uniqueNodes[1] = curNodes[1];
5095 uniqueNodes[2] = curNodes[3];
5096 uniqueNodes[3] = curNodes[4];
5097 uniqueNodes[4] = curNodes[2];
5098 uniqueNodes[5] = curNodes[7];
5101 if( iRepl[0]==3 && iRepl[1]==6 && iRepl[2]==7 ) {
5102 uniqueNodes[0] = curNodes[0];
5103 uniqueNodes[1] = curNodes[1];
5104 uniqueNodes[2] = curNodes[2];
5105 uniqueNodes[3] = curNodes[4];
5106 uniqueNodes[4] = curNodes[5];
5107 uniqueNodes[5] = curNodes[3];
5113 //////////////////////////////////// HEXAHEDRON
5115 SMDS_VolumeTool hexa (elem);
5116 hexa.SetExternalNormal();
5117 if ( nbUniqueNodes == 4 && nbRepl == 6 ) {
5118 //////////////////////// ---> tetrahedron
5119 for ( int iFace = 0; iFace < 6; iFace++ ) {
5120 const int *ind = hexa.GetFaceNodesIndices( iFace ); // indices of face nodes
5121 if (curNodes[ind[ 0 ]] == curNodes[ind[ 1 ]] &&
5122 curNodes[ind[ 0 ]] == curNodes[ind[ 2 ]] &&
5123 curNodes[ind[ 0 ]] == curNodes[ind[ 3 ]] ) {
5124 // one face turns into a point ...
5125 int iOppFace = hexa.GetOppFaceIndex( iFace );
5126 ind = hexa.GetFaceNodesIndices( iOppFace );
5128 iUnique = 2; // reverse a tetrahedron bottom
5129 for ( iCur = 0; iCur < 4 && nbStick < 2; iCur++ ) {
5130 if ( curNodes[ind[ iCur ]] == curNodes[ind[ iCur + 1 ]] )
5132 else if ( iUnique >= 0 )
5133 uniqueNodes[ iUnique-- ] = curNodes[ind[ iCur ]];
5135 if ( nbStick == 1 ) {
5136 // ... and the opposite one - into a triangle.
5138 ind = hexa.GetFaceNodesIndices( iFace );
5139 uniqueNodes[ 3 ] = curNodes[ind[ 0 ]];
5146 else if (nbUniqueNodes == 5 && nbRepl == 4 ) {
5147 //////////////////// HEXAHEDRON ---> 2 tetrahedrons
5148 for ( int iFace = 0; iFace < 6; iFace++ ) {
5149 const int *ind = hexa.GetFaceNodesIndices( iFace ); // indices of face nodes
5150 if (curNodes[ind[ 0 ]] == curNodes[ind[ 1 ]] &&
5151 curNodes[ind[ 0 ]] == curNodes[ind[ 2 ]] &&
5152 curNodes[ind[ 0 ]] == curNodes[ind[ 3 ]] ) {
5153 // one face turns into a point ...
5154 int iOppFace = hexa.GetOppFaceIndex( iFace );
5155 ind = hexa.GetFaceNodesIndices( iOppFace );
5157 iUnique = 2; // reverse a tetrahedron 1 bottom
5158 for ( iCur = 0; iCur < 4 && nbStick == 0; iCur++ ) {
5159 if ( curNodes[ind[ iCur ]] == curNodes[ind[ iCur + 1 ]] )
5161 else if ( iUnique >= 0 )
5162 uniqueNodes[ iUnique-- ] = curNodes[ind[ iCur ]];
5164 if ( nbStick == 0 ) {
5165 // ... and the opposite one is a quadrangle
5167 const int* indTop = hexa.GetFaceNodesIndices( iFace );
5168 uniqueNodes[ 3 ] = curNodes[indTop[ 0 ]];
5171 SMDS_MeshElement* newElem =
5172 aMesh->AddVolume(curNodes[ind[ 0 ]],
5175 curNodes[indTop[ 0 ]]);
5176 myLastCreatedElems.Append(newElem);
5178 aMesh->SetMeshElementOnShape( newElem, aShapeId );
5185 else if ( nbUniqueNodes == 6 && nbRepl == 4 ) {
5186 ////////////////// HEXAHEDRON ---> 2 tetrahedrons or 1 prism
5187 // find indices of quad and tri faces
5188 int iQuadFace[ 6 ], iTriFace[ 6 ], nbQuad = 0, nbTri = 0, iFace;
5189 for ( iFace = 0; iFace < 6; iFace++ ) {
5190 const int *ind = hexa.GetFaceNodesIndices( iFace ); // indices of face nodes
5192 for ( iCur = 0; iCur < 4; iCur++ )
5193 nodeSet.insert( curNodes[ind[ iCur ]] );
5194 nbUniqueNodes = nodeSet.size();
5195 if ( nbUniqueNodes == 3 )
5196 iTriFace[ nbTri++ ] = iFace;
5197 else if ( nbUniqueNodes == 4 )
5198 iQuadFace[ nbQuad++ ] = iFace;
5200 if (nbQuad == 2 && nbTri == 4 &&
5201 hexa.GetOppFaceIndex( iQuadFace[ 0 ] ) == iQuadFace[ 1 ]) {
5202 // 2 opposite quadrangles stuck with a diagonal;
5203 // sample groups of merged indices: (0-4)(2-6)
5204 // --------------------------------------------> 2 tetrahedrons
5205 const int *ind1 = hexa.GetFaceNodesIndices( iQuadFace[ 0 ]); // indices of quad1 nodes
5206 const int *ind2 = hexa.GetFaceNodesIndices( iQuadFace[ 1 ]);
5207 int i0, i1d, i2, i3d, i0t, i2t; // d-daigonal, t-top
5208 if (curNodes[ind1[ 0 ]] == curNodes[ind2[ 0 ]] &&
5209 curNodes[ind1[ 2 ]] == curNodes[ind2[ 2 ]]) {
5210 // stuck with 0-2 diagonal
5218 else if (curNodes[ind1[ 1 ]] == curNodes[ind2[ 3 ]] &&
5219 curNodes[ind1[ 3 ]] == curNodes[ind2[ 1 ]]) {
5220 // stuck with 1-3 diagonal
5232 uniqueNodes[ 0 ] = curNodes [ i0 ];
5233 uniqueNodes[ 1 ] = curNodes [ i1d ];
5234 uniqueNodes[ 2 ] = curNodes [ i3d ];
5235 uniqueNodes[ 3 ] = curNodes [ i0t ];
5238 SMDS_MeshElement* newElem = aMesh->AddVolume(curNodes[ i1d ],
5242 myLastCreatedElems.Append(newElem);
5244 aMesh->SetMeshElementOnShape( newElem, aShapeId );
5247 else if (( nbTri == 2 && nbQuad == 3 ) || // merged (0-4)(1-5)
5248 ( nbTri == 4 && nbQuad == 2 )) { // merged (7-4)(1-5)
5249 // --------------------------------------------> prism
5250 // find 2 opposite triangles
5252 for ( iFace = 0; iFace + 1 < nbTri; iFace++ ) {
5253 if ( hexa.GetOppFaceIndex( iTriFace[ iFace ] ) == iTriFace[ iFace + 1 ]) {
5254 // find indices of kept and replaced nodes
5255 // and fill unique nodes of 2 opposite triangles
5256 const int *ind1 = hexa.GetFaceNodesIndices( iTriFace[ iFace ]);
5257 const int *ind2 = hexa.GetFaceNodesIndices( iTriFace[ iFace + 1 ]);
5258 const SMDS_MeshNode** hexanodes = hexa.GetNodes();
5259 // fill unique nodes
5262 for ( iCur = 0; iCur < 4 && isOk; iCur++ ) {
5263 const SMDS_MeshNode* n = curNodes[ind1[ iCur ]];
5264 const SMDS_MeshNode* nInit = hexanodes[ind1[ iCur ]];
5266 // iCur of a linked node of the opposite face (make normals co-directed):
5267 int iCurOpp = ( iCur == 1 || iCur == 3 ) ? 4 - iCur : iCur;
5268 // check that correspondent corners of triangles are linked
5269 if ( !hexa.IsLinked( ind1[ iCur ], ind2[ iCurOpp ] ))
5272 uniqueNodes[ iUnique ] = n;
5273 uniqueNodes[ iUnique + 3 ] = curNodes[ind2[ iCurOpp ]];
5282 } // if ( nbUniqueNodes == 6 && nbRepl == 4 )
5288 } // switch ( nbNodes )
5290 } // if ( nbNodes != nbUniqueNodes ) // some nodes stick
5293 if (elem->IsPoly() && elem->GetType() == SMDSAbs_Volume) {
5294 // Change nodes of polyedre
5295 const SMDS_PolyhedralVolumeOfNodes* aPolyedre =
5296 static_cast<const SMDS_PolyhedralVolumeOfNodes*>( elem );
5298 int nbFaces = aPolyedre->NbFaces();
5300 vector<const SMDS_MeshNode *> poly_nodes;
5301 vector<int> quantities (nbFaces);
5303 for (int iface = 1; iface <= nbFaces; iface++) {
5304 int inode, nbFaceNodes = aPolyedre->NbFaceNodes(iface);
5305 quantities[iface - 1] = nbFaceNodes;
5307 for (inode = 1; inode <= nbFaceNodes; inode++) {
5308 const SMDS_MeshNode* curNode = aPolyedre->GetFaceNode(iface, inode);
5310 TNodeNodeMap::iterator nnIt = nodeNodeMap.find( curNode );
5311 if (nnIt != nodeNodeMap.end()) { // curNode sticks
5312 curNode = (*nnIt).second;
5314 poly_nodes.push_back(curNode);
5317 aMesh->ChangePolyhedronNodes( elem, poly_nodes, quantities );
5321 // Change regular element or polygon
5322 aMesh->ChangeElementNodes( elem, & uniqueNodes[0], nbUniqueNodes );
5326 // Remove invalid regular element or invalid polygon
5327 rmElemIds.push_back( elem->GetID() );
5330 } // loop on elements
5332 // Remove equal nodes and bad elements
5334 Remove( rmNodeIds, true );
5335 Remove( rmElemIds, false );
5340 // ========================================================
5341 // class : SortableElement
5342 // purpose : allow sorting elements basing on their nodes
5343 // ========================================================
5344 class SortableElement : public set <const SMDS_MeshElement*>
5348 SortableElement( const SMDS_MeshElement* theElem )
5351 SMDS_ElemIteratorPtr nodeIt = theElem->nodesIterator();
5352 while ( nodeIt->more() )
5353 this->insert( nodeIt->next() );
5356 const SMDS_MeshElement* Get() const
5359 void Set(const SMDS_MeshElement* e) const
5364 mutable const SMDS_MeshElement* myElem;
5367 //=======================================================================
5368 //function : FindEqualElements
5369 //purpose : Return list of group of elements built on the same nodes.
5370 // Search among theElements or in the whole mesh if theElements is empty
5371 //=======================================================================
5372 void SMESH_MeshEditor::FindEqualElements(set<const SMDS_MeshElement*> & theElements,
5373 TListOfListOfElementsID & theGroupsOfElementsID)
5375 myLastCreatedElems.Clear();
5376 myLastCreatedNodes.Clear();
5378 typedef set<const SMDS_MeshElement*> TElemsSet;
5379 typedef map< SortableElement, int > TMapOfNodeSet;
5380 typedef list<int> TGroupOfElems;
5383 if ( theElements.empty() )
5384 { // get all elements in the mesh
5385 SMDS_ElemIteratorPtr eIt = GetMeshDS()->elementsIterator();
5386 while ( eIt->more() )
5387 elems.insert( elems.end(), eIt->next());
5390 elems = theElements;
5392 vector< TGroupOfElems > arrayOfGroups;
5393 TGroupOfElems groupOfElems;
5394 TMapOfNodeSet mapOfNodeSet;
5396 TElemsSet::iterator elemIt = elems.begin();
5397 for ( int i = 0, j=0; elemIt != elems.end(); ++elemIt, ++j ) {
5398 const SMDS_MeshElement* curElem = *elemIt;
5399 SortableElement SE(curElem);
5402 pair< TMapOfNodeSet::iterator, bool> pp = mapOfNodeSet.insert(make_pair(SE, i));
5403 if( !(pp.second) ) {
5404 TMapOfNodeSet::iterator& itSE = pp.first;
5405 ind = (*itSE).second;
5406 arrayOfGroups[ind].push_back(curElem->GetID());
5409 groupOfElems.clear();
5410 groupOfElems.push_back(curElem->GetID());
5411 arrayOfGroups.push_back(groupOfElems);
5416 vector< TGroupOfElems >::iterator groupIt = arrayOfGroups.begin();
5417 for ( ; groupIt != arrayOfGroups.end(); ++groupIt ) {
5418 groupOfElems = *groupIt;
5419 if ( groupOfElems.size() > 1 ) {
5420 groupOfElems.sort();
5421 theGroupsOfElementsID.push_back(groupOfElems);
5426 //=======================================================================
5427 //function : MergeElements
5428 //purpose : In each given group, substitute all elements by the first one.
5429 //=======================================================================
5431 void SMESH_MeshEditor::MergeElements(TListOfListOfElementsID & theGroupsOfElementsID)
5433 myLastCreatedElems.Clear();
5434 myLastCreatedNodes.Clear();
5436 typedef list<int> TListOfIDs;
5437 TListOfIDs rmElemIds; // IDs of elems to remove
5439 SMESHDS_Mesh* aMesh = GetMeshDS();
5441 TListOfListOfElementsID::iterator groupsIt = theGroupsOfElementsID.begin();
5442 while ( groupsIt != theGroupsOfElementsID.end() ) {
5443 TListOfIDs& aGroupOfElemID = *groupsIt;
5444 aGroupOfElemID.sort();
5445 int elemIDToKeep = aGroupOfElemID.front();
5446 const SMDS_MeshElement* elemToKeep = aMesh->FindElement(elemIDToKeep);
5447 aGroupOfElemID.pop_front();
5448 TListOfIDs::iterator idIt = aGroupOfElemID.begin();
5449 while ( idIt != aGroupOfElemID.end() ) {
5450 int elemIDToRemove = *idIt;
5451 const SMDS_MeshElement* elemToRemove = aMesh->FindElement(elemIDToRemove);
5452 // add the kept element in groups of removed one (PAL15188)
5453 AddToSameGroups( elemToKeep, elemToRemove, aMesh );
5454 rmElemIds.push_back( elemIDToRemove );
5460 Remove( rmElemIds, false );
5463 //=======================================================================
5464 //function : MergeEqualElements
5465 //purpose : Remove all but one of elements built on the same nodes.
5466 //=======================================================================
5468 void SMESH_MeshEditor::MergeEqualElements()
5470 set<const SMDS_MeshElement*> aMeshElements; /* empty input -
5471 to merge equal elements in the whole mesh */
5472 TListOfListOfElementsID aGroupsOfElementsID;
5473 FindEqualElements(aMeshElements, aGroupsOfElementsID);
5474 MergeElements(aGroupsOfElementsID);
5477 //=======================================================================
5478 //function : FindFaceInSet
5479 //purpose : Return a face having linked nodes n1 and n2 and which is
5480 // - not in avoidSet,
5481 // - in elemSet provided that !elemSet.empty()
5482 //=======================================================================
5484 const SMDS_MeshElement*
5485 SMESH_MeshEditor::FindFaceInSet(const SMDS_MeshNode* n1,
5486 const SMDS_MeshNode* n2,
5487 const TIDSortedElemSet& elemSet,
5488 const TIDSortedElemSet& avoidSet)
5491 SMDS_ElemIteratorPtr invElemIt = n1->GetInverseElementIterator(SMDSAbs_Face);
5492 while ( invElemIt->more() ) { // loop on inverse elements of n1
5493 const SMDS_MeshElement* elem = invElemIt->next();
5494 if (avoidSet.find( elem ) != avoidSet.end() )
5496 if ( !elemSet.empty() && elemSet.find( elem ) == elemSet.end())
5498 // get face nodes and find index of n1
5499 int i1, nbN = elem->NbNodes(), iNode = 0;
5500 //const SMDS_MeshNode* faceNodes[ nbN ], *n;
5501 vector<const SMDS_MeshNode*> faceNodes( nbN );
5502 const SMDS_MeshNode* n;
5503 SMDS_ElemIteratorPtr nIt = elem->nodesIterator();
5504 while ( nIt->more() ) {
5505 faceNodes[ iNode ] = static_cast<const SMDS_MeshNode*>( nIt->next() );
5506 if ( faceNodes[ iNode++ ] == n1 )
5509 // find a n2 linked to n1
5510 if(!elem->IsQuadratic()) {
5511 for ( iNode = 0; iNode < 2; iNode++ ) {
5512 if ( iNode ) // node before n1
5513 n = faceNodes[ i1 == 0 ? nbN - 1 : i1 - 1 ];
5514 else // node after n1
5515 n = faceNodes[ i1 + 1 == nbN ? 0 : i1 + 1 ];
5520 else { // analysis for quadratic elements
5521 bool IsFind = false;
5522 // check using only corner nodes
5523 for ( iNode = 0; iNode < 2; iNode++ ) {
5524 if ( iNode ) // node before n1
5525 n = faceNodes[ i1 == 0 ? nbN/2 - 1 : i1 - 1 ];
5526 else // node after n1
5527 n = faceNodes[ i1 + 1 == nbN/2 ? 0 : i1 + 1 ];
5535 // check using all nodes
5536 const SMDS_QuadraticFaceOfNodes* F =
5537 static_cast<const SMDS_QuadraticFaceOfNodes*>(elem);
5538 // use special nodes iterator
5540 SMDS_NodeIteratorPtr anIter = F->interlacedNodesIterator();
5541 while ( anIter->more() ) {
5542 faceNodes[iNode] = static_cast<const SMDS_MeshNode*>(anIter->next());
5543 if ( faceNodes[ iNode++ ] == n1 )
5546 for ( iNode = 0; iNode < 2; iNode++ ) {
5547 if ( iNode ) // node before n1
5548 n = faceNodes[ i1 == 0 ? nbN - 1 : i1 - 1 ];
5549 else // node after n1
5550 n = faceNodes[ i1 + 1 == nbN ? 0 : i1 + 1 ];
5556 } // end analysis for quadratic elements
5561 //=======================================================================
5562 //function : findAdjacentFace
5564 //=======================================================================
5566 static const SMDS_MeshElement* findAdjacentFace(const SMDS_MeshNode* n1,
5567 const SMDS_MeshNode* n2,
5568 const SMDS_MeshElement* elem)
5570 TIDSortedElemSet elemSet, avoidSet;
5572 avoidSet.insert ( elem );
5573 return SMESH_MeshEditor::FindFaceInSet( n1, n2, elemSet, avoidSet );
5576 //=======================================================================
5577 //function : FindFreeBorder
5579 //=======================================================================
5581 #define ControlFreeBorder SMESH::Controls::FreeEdges::IsFreeEdge
5583 bool SMESH_MeshEditor::FindFreeBorder (const SMDS_MeshNode* theFirstNode,
5584 const SMDS_MeshNode* theSecondNode,
5585 const SMDS_MeshNode* theLastNode,
5586 list< const SMDS_MeshNode* > & theNodes,
5587 list< const SMDS_MeshElement* >& theFaces)
5589 if ( !theFirstNode || !theSecondNode )
5591 // find border face between theFirstNode and theSecondNode
5592 const SMDS_MeshElement* curElem = findAdjacentFace( theFirstNode, theSecondNode, 0 );
5596 theFaces.push_back( curElem );
5597 theNodes.push_back( theFirstNode );
5598 theNodes.push_back( theSecondNode );
5600 //vector<const SMDS_MeshNode*> nodes;
5601 const SMDS_MeshNode *nIgnore = theFirstNode, *nStart = theSecondNode;
5602 set < const SMDS_MeshElement* > foundElems;
5603 bool needTheLast = ( theLastNode != 0 );
5605 while ( nStart != theLastNode ) {
5606 if ( nStart == theFirstNode )
5607 return !needTheLast;
5609 // find all free border faces sharing form nStart
5611 list< const SMDS_MeshElement* > curElemList;
5612 list< const SMDS_MeshNode* > nStartList;
5613 SMDS_ElemIteratorPtr invElemIt = nStart->GetInverseElementIterator(SMDSAbs_Face);
5614 while ( invElemIt->more() ) {
5615 const SMDS_MeshElement* e = invElemIt->next();
5616 if ( e == curElem || foundElems.insert( e ).second ) {
5618 int iNode = 0, nbNodes = e->NbNodes();
5619 //const SMDS_MeshNode* nodes[nbNodes+1];
5620 vector<const SMDS_MeshNode*> nodes(nbNodes+1);
5622 if(e->IsQuadratic()) {
5623 const SMDS_QuadraticFaceOfNodes* F =
5624 static_cast<const SMDS_QuadraticFaceOfNodes*>(e);
5625 // use special nodes iterator
5626 SMDS_NodeIteratorPtr anIter = F->interlacedNodesIterator();
5627 while( anIter->more() ) {
5628 nodes[ iNode++ ] = anIter->next();
5632 SMDS_ElemIteratorPtr nIt = e->nodesIterator();
5633 while ( nIt->more() )
5634 nodes[ iNode++ ] = static_cast<const SMDS_MeshNode*>( nIt->next() );
5636 nodes[ iNode ] = nodes[ 0 ];
5638 for ( iNode = 0; iNode < nbNodes; iNode++ )
5639 if (((nodes[ iNode ] == nStart && nodes[ iNode + 1] != nIgnore ) ||
5640 (nodes[ iNode + 1] == nStart && nodes[ iNode ] != nIgnore )) &&
5641 ControlFreeBorder( &nodes[ iNode ], e->GetID() ))
5643 nStartList.push_back( nodes[ iNode + ( nodes[ iNode ] == nStart ? 1 : 0 )]);
5644 curElemList.push_back( e );
5648 // analyse the found
5650 int nbNewBorders = curElemList.size();
5651 if ( nbNewBorders == 0 ) {
5652 // no free border furthermore
5653 return !needTheLast;
5655 else if ( nbNewBorders == 1 ) {
5656 // one more element found
5658 nStart = nStartList.front();
5659 curElem = curElemList.front();
5660 theFaces.push_back( curElem );
5661 theNodes.push_back( nStart );
5664 // several continuations found
5665 list< const SMDS_MeshElement* >::iterator curElemIt;
5666 list< const SMDS_MeshNode* >::iterator nStartIt;
5667 // check if one of them reached the last node
5668 if ( needTheLast ) {
5669 for (curElemIt = curElemList.begin(), nStartIt = nStartList.begin();
5670 curElemIt!= curElemList.end();
5671 curElemIt++, nStartIt++ )
5672 if ( *nStartIt == theLastNode ) {
5673 theFaces.push_back( *curElemIt );
5674 theNodes.push_back( *nStartIt );
5678 // find the best free border by the continuations
5679 list<const SMDS_MeshNode*> contNodes[ 2 ], *cNL;
5680 list<const SMDS_MeshElement*> contFaces[ 2 ], *cFL;
5681 for (curElemIt = curElemList.begin(), nStartIt = nStartList.begin();
5682 curElemIt!= curElemList.end();
5683 curElemIt++, nStartIt++ )
5685 cNL = & contNodes[ contNodes[0].empty() ? 0 : 1 ];
5686 cFL = & contFaces[ contFaces[0].empty() ? 0 : 1 ];
5687 // find one more free border
5688 if ( ! FindFreeBorder( nStart, *nStartIt, theLastNode, *cNL, *cFL )) {
5692 else if ( !contNodes[0].empty() && !contNodes[1].empty() ) {
5693 // choice: clear a worse one
5694 int iLongest = ( contNodes[0].size() < contNodes[1].size() ? 1 : 0 );
5695 int iWorse = ( needTheLast ? 1 - iLongest : iLongest );
5696 contNodes[ iWorse ].clear();
5697 contFaces[ iWorse ].clear();
5700 if ( contNodes[0].empty() && contNodes[1].empty() )
5703 // append the best free border
5704 cNL = & contNodes[ contNodes[0].empty() ? 1 : 0 ];
5705 cFL = & contFaces[ contFaces[0].empty() ? 1 : 0 ];
5706 theNodes.pop_back(); // remove nIgnore
5707 theNodes.pop_back(); // remove nStart
5708 theFaces.pop_back(); // remove curElem
5709 list< const SMDS_MeshNode* >::iterator nIt = cNL->begin();
5710 list< const SMDS_MeshElement* >::iterator fIt = cFL->begin();
5711 for ( ; nIt != cNL->end(); nIt++ ) theNodes.push_back( *nIt );
5712 for ( ; fIt != cFL->end(); fIt++ ) theFaces.push_back( *fIt );
5715 } // several continuations found
5716 } // while ( nStart != theLastNode )
5721 //=======================================================================
5722 //function : CheckFreeBorderNodes
5723 //purpose : Return true if the tree nodes are on a free border
5724 //=======================================================================
5726 bool SMESH_MeshEditor::CheckFreeBorderNodes(const SMDS_MeshNode* theNode1,
5727 const SMDS_MeshNode* theNode2,
5728 const SMDS_MeshNode* theNode3)
5730 list< const SMDS_MeshNode* > nodes;
5731 list< const SMDS_MeshElement* > faces;
5732 return FindFreeBorder( theNode1, theNode2, theNode3, nodes, faces);
5735 //=======================================================================
5736 //function : SewFreeBorder
5738 //=======================================================================
5740 SMESH_MeshEditor::Sew_Error
5741 SMESH_MeshEditor::SewFreeBorder (const SMDS_MeshNode* theBordFirstNode,
5742 const SMDS_MeshNode* theBordSecondNode,
5743 const SMDS_MeshNode* theBordLastNode,
5744 const SMDS_MeshNode* theSideFirstNode,
5745 const SMDS_MeshNode* theSideSecondNode,
5746 const SMDS_MeshNode* theSideThirdNode,
5747 const bool theSideIsFreeBorder,
5748 const bool toCreatePolygons,
5749 const bool toCreatePolyedrs)
5751 myLastCreatedElems.Clear();
5752 myLastCreatedNodes.Clear();
5754 MESSAGE("::SewFreeBorder()");
5755 Sew_Error aResult = SEW_OK;
5757 // ====================================
5758 // find side nodes and elements
5759 // ====================================
5761 list< const SMDS_MeshNode* > nSide[ 2 ];
5762 list< const SMDS_MeshElement* > eSide[ 2 ];
5763 list< const SMDS_MeshNode* >::iterator nIt[ 2 ];
5764 list< const SMDS_MeshElement* >::iterator eIt[ 2 ];
5768 if (!FindFreeBorder(theBordFirstNode,theBordSecondNode,theBordLastNode,
5769 nSide[0], eSide[0])) {
5770 MESSAGE(" Free Border 1 not found " );
5771 aResult = SEW_BORDER1_NOT_FOUND;
5773 if (theSideIsFreeBorder) {
5776 if (!FindFreeBorder(theSideFirstNode, theSideSecondNode, theSideThirdNode,
5777 nSide[1], eSide[1])) {
5778 MESSAGE(" Free Border 2 not found " );
5779 aResult = ( aResult != SEW_OK ? SEW_BOTH_BORDERS_NOT_FOUND : SEW_BORDER2_NOT_FOUND );
5782 if ( aResult != SEW_OK )
5785 if (!theSideIsFreeBorder) {
5789 // -------------------------------------------------------------------------
5791 // 1. If nodes to merge are not coincident, move nodes of the free border
5792 // from the coord sys defined by the direction from the first to last
5793 // nodes of the border to the correspondent sys of the side 2
5794 // 2. On the side 2, find the links most co-directed with the correspondent
5795 // links of the free border
5796 // -------------------------------------------------------------------------
5798 // 1. Since sewing may brake if there are volumes to split on the side 2,
5799 // we wont move nodes but just compute new coordinates for them
5800 typedef map<const SMDS_MeshNode*, gp_XYZ> TNodeXYZMap;
5801 TNodeXYZMap nBordXYZ;
5802 list< const SMDS_MeshNode* >& bordNodes = nSide[ 0 ];
5803 list< const SMDS_MeshNode* >::iterator nBordIt;
5805 gp_XYZ Pb1( theBordFirstNode->X(), theBordFirstNode->Y(), theBordFirstNode->Z() );
5806 gp_XYZ Pb2( theBordLastNode->X(), theBordLastNode->Y(), theBordLastNode->Z() );
5807 gp_XYZ Ps1( theSideFirstNode->X(), theSideFirstNode->Y(), theSideFirstNode->Z() );
5808 gp_XYZ Ps2( theSideSecondNode->X(), theSideSecondNode->Y(), theSideSecondNode->Z() );
5809 double tol2 = 1.e-8;
5810 gp_Vec Vbs1( Pb1 - Ps1 ),Vbs2( Pb2 - Ps2 );
5811 if ( Vbs1.SquareMagnitude() > tol2 || Vbs2.SquareMagnitude() > tol2 ) {
5812 // Need node movement.
5814 // find X and Z axes to create trsf
5815 gp_Vec Zb( Pb1 - Pb2 ), Zs( Ps1 - Ps2 );
5817 if ( X.SquareMagnitude() <= gp::Resolution() * gp::Resolution() )
5819 X = gp_Ax2( gp::Origin(), Zb ).XDirection();
5822 gp_Ax3 toBordAx( Pb1, Zb, X );
5823 gp_Ax3 fromSideAx( Ps1, Zs, X );
5824 gp_Ax3 toGlobalAx( gp::Origin(), gp::DZ(), gp::DX() );
5826 gp_Trsf toBordSys, fromSide2Sys;
5827 toBordSys.SetTransformation( toBordAx );
5828 fromSide2Sys.SetTransformation( fromSideAx, toGlobalAx );
5829 fromSide2Sys.SetScaleFactor( Zs.Magnitude() / Zb.Magnitude() );
5832 for ( nBordIt = bordNodes.begin(); nBordIt != bordNodes.end(); nBordIt++ ) {
5833 const SMDS_MeshNode* n = *nBordIt;
5834 gp_XYZ xyz( n->X(),n->Y(),n->Z() );
5835 toBordSys.Transforms( xyz );
5836 fromSide2Sys.Transforms( xyz );
5837 nBordXYZ.insert( TNodeXYZMap::value_type( n, xyz ));
5841 // just insert nodes XYZ in the nBordXYZ map
5842 for ( nBordIt = bordNodes.begin(); nBordIt != bordNodes.end(); nBordIt++ ) {
5843 const SMDS_MeshNode* n = *nBordIt;
5844 nBordXYZ.insert( TNodeXYZMap::value_type( n, gp_XYZ( n->X(),n->Y(),n->Z() )));
5848 // 2. On the side 2, find the links most co-directed with the correspondent
5849 // links of the free border
5851 list< const SMDS_MeshElement* >& sideElems = eSide[ 1 ];
5852 list< const SMDS_MeshNode* >& sideNodes = nSide[ 1 ];
5853 sideNodes.push_back( theSideFirstNode );
5855 bool hasVolumes = false;
5856 LinkID_Gen aLinkID_Gen( GetMeshDS() );
5857 set<long> foundSideLinkIDs, checkedLinkIDs;
5858 SMDS_VolumeTool volume;
5859 //const SMDS_MeshNode* faceNodes[ 4 ];
5861 const SMDS_MeshNode* sideNode;
5862 const SMDS_MeshElement* sideElem;
5863 const SMDS_MeshNode* prevSideNode = theSideFirstNode;
5864 const SMDS_MeshNode* prevBordNode = theBordFirstNode;
5865 nBordIt = bordNodes.begin();
5867 // border node position and border link direction to compare with
5868 gp_XYZ bordPos = nBordXYZ[ *nBordIt ];
5869 gp_XYZ bordDir = bordPos - nBordXYZ[ prevBordNode ];
5870 // choose next side node by link direction or by closeness to
5871 // the current border node:
5872 bool searchByDir = ( *nBordIt != theBordLastNode );
5874 // find the next node on the Side 2
5876 double maxDot = -DBL_MAX, minDist = DBL_MAX;
5878 checkedLinkIDs.clear();
5879 gp_XYZ prevXYZ( prevSideNode->X(), prevSideNode->Y(), prevSideNode->Z() );
5881 // loop on inverse elements of current node (prevSideNode) on the Side 2
5882 SMDS_ElemIteratorPtr invElemIt = prevSideNode->GetInverseElementIterator();
5883 while ( invElemIt->more() )
5885 const SMDS_MeshElement* elem = invElemIt->next();
5886 // prepare data for a loop on links coming to prevSideNode, of a face or a volume
5887 int iPrevNode, iNode = 0, nbNodes = elem->NbNodes();
5888 vector< const SMDS_MeshNode* > faceNodes( nbNodes, (const SMDS_MeshNode*)0 );
5889 bool isVolume = volume.Set( elem );
5890 const SMDS_MeshNode** nodes = isVolume ? volume.GetNodes() : & faceNodes[0];
5891 if ( isVolume ) // --volume
5893 else if ( elem->GetType()==SMDSAbs_Face ) { // --face
5894 // retrieve all face nodes and find iPrevNode - an index of the prevSideNode
5895 if(elem->IsQuadratic()) {
5896 const SMDS_QuadraticFaceOfNodes* F =
5897 static_cast<const SMDS_QuadraticFaceOfNodes*>(elem);
5898 // use special nodes iterator
5899 SMDS_NodeIteratorPtr anIter = F->interlacedNodesIterator();
5900 while( anIter->more() ) {
5901 nodes[ iNode ] = anIter->next();
5902 if ( nodes[ iNode++ ] == prevSideNode )
5903 iPrevNode = iNode - 1;
5907 SMDS_ElemIteratorPtr nIt = elem->nodesIterator();
5908 while ( nIt->more() ) {
5909 nodes[ iNode ] = cast2Node( nIt->next() );
5910 if ( nodes[ iNode++ ] == prevSideNode )
5911 iPrevNode = iNode - 1;
5914 // there are 2 links to check
5919 // loop on links, to be precise, on the second node of links
5920 for ( iNode = 0; iNode < nbNodes; iNode++ ) {
5921 const SMDS_MeshNode* n = nodes[ iNode ];
5923 if ( !volume.IsLinked( n, prevSideNode ))
5927 if ( iNode ) // a node before prevSideNode
5928 n = nodes[ iPrevNode == 0 ? elem->NbNodes() - 1 : iPrevNode - 1 ];
5929 else // a node after prevSideNode
5930 n = nodes[ iPrevNode + 1 == elem->NbNodes() ? 0 : iPrevNode + 1 ];
5932 // check if this link was already used
5933 long iLink = aLinkID_Gen.GetLinkID( prevSideNode, n );
5934 bool isJustChecked = !checkedLinkIDs.insert( iLink ).second;
5935 if (!isJustChecked &&
5936 foundSideLinkIDs.find( iLink ) == foundSideLinkIDs.end() )
5938 // test a link geometrically
5939 gp_XYZ nextXYZ ( n->X(), n->Y(), n->Z() );
5940 bool linkIsBetter = false;
5941 double dot = 0.0, dist = 0.0;
5942 if ( searchByDir ) { // choose most co-directed link
5943 dot = bordDir * ( nextXYZ - prevXYZ ).Normalized();
5944 linkIsBetter = ( dot > maxDot );
5946 else { // choose link with the node closest to bordPos
5947 dist = ( nextXYZ - bordPos ).SquareModulus();
5948 linkIsBetter = ( dist < minDist );
5950 if ( linkIsBetter ) {
5959 } // loop on inverse elements of prevSideNode
5962 MESSAGE(" Cant find path by links of the Side 2 ");
5963 return SEW_BAD_SIDE_NODES;
5965 sideNodes.push_back( sideNode );
5966 sideElems.push_back( sideElem );
5967 foundSideLinkIDs.insert ( linkID );
5968 prevSideNode = sideNode;
5970 if ( *nBordIt == theBordLastNode )
5971 searchByDir = false;
5973 // find the next border link to compare with
5974 gp_XYZ sidePos( sideNode->X(), sideNode->Y(), sideNode->Z() );
5975 searchByDir = ( bordDir * ( sidePos - bordPos ) <= 0 );
5976 // move to next border node if sideNode is before forward border node (bordPos)
5977 while ( *nBordIt != theBordLastNode && !searchByDir ) {
5978 prevBordNode = *nBordIt;
5980 bordPos = nBordXYZ[ *nBordIt ];
5981 bordDir = bordPos - nBordXYZ[ prevBordNode ];
5982 searchByDir = ( bordDir * ( sidePos - bordPos ) <= 0 );
5986 while ( sideNode != theSideSecondNode );
5988 if ( hasVolumes && sideNodes.size () != bordNodes.size() && !toCreatePolyedrs) {
5989 MESSAGE("VOLUME SPLITTING IS FORBIDDEN");
5990 return SEW_VOLUMES_TO_SPLIT; // volume splitting is forbidden
5992 } // end nodes search on the side 2
5994 // ============================
5995 // sew the border to the side 2
5996 // ============================
5998 int nbNodes[] = { nSide[0].size(), nSide[1].size() };
5999 int maxNbNodes = Max( nbNodes[0], nbNodes[1] );
6001 TListOfListOfNodes nodeGroupsToMerge;
6002 if ( nbNodes[0] == nbNodes[1] ||
6003 ( theSideIsFreeBorder && !theSideThirdNode)) {
6005 // all nodes are to be merged
6007 for (nIt[0] = nSide[0].begin(), nIt[1] = nSide[1].begin();
6008 nIt[0] != nSide[0].end() && nIt[1] != nSide[1].end();
6009 nIt[0]++, nIt[1]++ )
6011 nodeGroupsToMerge.push_back( list<const SMDS_MeshNode*>() );
6012 nodeGroupsToMerge.back().push_back( *nIt[1] ); // to keep
6013 nodeGroupsToMerge.back().push_back( *nIt[0] ); // to remove
6018 // insert new nodes into the border and the side to get equal nb of segments
6020 // get normalized parameters of nodes on the borders
6021 //double param[ 2 ][ maxNbNodes ];
6023 param[0] = new double [ maxNbNodes ];
6024 param[1] = new double [ maxNbNodes ];
6026 for ( iBord = 0; iBord < 2; iBord++ ) { // loop on 2 borders
6027 list< const SMDS_MeshNode* >& nodes = nSide[ iBord ];
6028 list< const SMDS_MeshNode* >::iterator nIt = nodes.begin();
6029 const SMDS_MeshNode* nPrev = *nIt;
6030 double bordLength = 0;
6031 for ( iNode = 0; nIt != nodes.end(); nIt++, iNode++ ) { // loop on border nodes
6032 const SMDS_MeshNode* nCur = *nIt;
6033 gp_XYZ segment (nCur->X() - nPrev->X(),
6034 nCur->Y() - nPrev->Y(),
6035 nCur->Z() - nPrev->Z());
6036 double segmentLen = segment.Modulus();
6037 bordLength += segmentLen;
6038 param[ iBord ][ iNode ] = bordLength;
6041 // normalize within [0,1]
6042 for ( iNode = 0; iNode < nbNodes[ iBord ]; iNode++ ) {
6043 param[ iBord ][ iNode ] /= bordLength;
6047 // loop on border segments
6048 const SMDS_MeshNode *nPrev[ 2 ] = { 0, 0 };
6049 int i[ 2 ] = { 0, 0 };
6050 nIt[0] = nSide[0].begin(); eIt[0] = eSide[0].begin();
6051 nIt[1] = nSide[1].begin(); eIt[1] = eSide[1].begin();
6053 TElemOfNodeListMap insertMap;
6054 TElemOfNodeListMap::iterator insertMapIt;
6056 // key: elem to insert nodes into
6057 // value: 2 nodes to insert between + nodes to be inserted
6059 bool next[ 2 ] = { false, false };
6061 // find min adjacent segment length after sewing
6062 double nextParam = 10., prevParam = 0;
6063 for ( iBord = 0; iBord < 2; iBord++ ) { // loop on 2 borders
6064 if ( i[ iBord ] + 1 < nbNodes[ iBord ])
6065 nextParam = Min( nextParam, param[iBord][ i[iBord] + 1 ]);
6066 if ( i[ iBord ] > 0 )
6067 prevParam = Max( prevParam, param[iBord][ i[iBord] - 1 ]);
6069 double minParam = Min( param[ 0 ][ i[0] ], param[ 1 ][ i[1] ]);
6070 double maxParam = Max( param[ 0 ][ i[0] ], param[ 1 ][ i[1] ]);
6071 double minSegLen = Min( nextParam - minParam, maxParam - prevParam );
6073 // choose to insert or to merge nodes
6074 double du = param[ 1 ][ i[1] ] - param[ 0 ][ i[0] ];
6075 if ( Abs( du ) <= minSegLen * 0.2 ) {
6078 nodeGroupsToMerge.push_back( list<const SMDS_MeshNode*>() );
6079 const SMDS_MeshNode* n0 = *nIt[0];
6080 const SMDS_MeshNode* n1 = *nIt[1];
6081 nodeGroupsToMerge.back().push_back( n1 );
6082 nodeGroupsToMerge.back().push_back( n0 );
6083 // position of node of the border changes due to merge
6084 param[ 0 ][ i[0] ] += du;
6085 // move n1 for the sake of elem shape evaluation during insertion.
6086 // n1 will be removed by MergeNodes() anyway
6087 const_cast<SMDS_MeshNode*>( n0 )->setXYZ( n1->X(), n1->Y(), n1->Z() );
6088 next[0] = next[1] = true;
6093 int intoBord = ( du < 0 ) ? 0 : 1;
6094 const SMDS_MeshElement* elem = *eIt[ intoBord ];
6095 const SMDS_MeshNode* n1 = nPrev[ intoBord ];
6096 const SMDS_MeshNode* n2 = *nIt[ intoBord ];
6097 const SMDS_MeshNode* nIns = *nIt[ 1 - intoBord ];
6098 if ( intoBord == 1 ) {
6099 // move node of the border to be on a link of elem of the side
6100 gp_XYZ p1 (n1->X(), n1->Y(), n1->Z());
6101 gp_XYZ p2 (n2->X(), n2->Y(), n2->Z());
6102 double ratio = du / ( param[ 1 ][ i[1] ] - param[ 1 ][ i[1]-1 ]);
6103 gp_XYZ p = p2 * ( 1 - ratio ) + p1 * ratio;
6104 GetMeshDS()->MoveNode( nIns, p.X(), p.Y(), p.Z() );
6106 insertMapIt = insertMap.find( elem );
6107 bool notFound = ( insertMapIt == insertMap.end() );
6108 bool otherLink = ( !notFound && (*insertMapIt).second.front() != n1 );
6110 // insert into another link of the same element:
6111 // 1. perform insertion into the other link of the elem
6112 list<const SMDS_MeshNode*> & nodeList = (*insertMapIt).second;
6113 const SMDS_MeshNode* n12 = nodeList.front(); nodeList.pop_front();
6114 const SMDS_MeshNode* n22 = nodeList.front(); nodeList.pop_front();
6115 InsertNodesIntoLink( elem, n12, n22, nodeList, toCreatePolygons );
6116 // 2. perform insertion into the link of adjacent faces
6118 const SMDS_MeshElement* adjElem = findAdjacentFace( n12, n22, elem );
6120 InsertNodesIntoLink( adjElem, n12, n22, nodeList, toCreatePolygons );
6124 if (toCreatePolyedrs) {
6125 // perform insertion into the links of adjacent volumes
6126 UpdateVolumes(n12, n22, nodeList);
6128 // 3. find an element appeared on n1 and n2 after the insertion
6129 insertMap.erase( elem );
6130 elem = findAdjacentFace( n1, n2, 0 );
6132 if ( notFound || otherLink ) {
6133 // add element and nodes of the side into the insertMap
6134 insertMapIt = insertMap.insert
6135 ( TElemOfNodeListMap::value_type( elem, list<const SMDS_MeshNode*>() )).first;
6136 (*insertMapIt).second.push_back( n1 );
6137 (*insertMapIt).second.push_back( n2 );
6139 // add node to be inserted into elem
6140 (*insertMapIt).second.push_back( nIns );
6141 next[ 1 - intoBord ] = true;
6144 // go to the next segment
6145 for ( iBord = 0; iBord < 2; iBord++ ) { // loop on 2 borders
6146 if ( next[ iBord ] ) {
6147 if ( i[ iBord ] != 0 && eIt[ iBord ] != eSide[ iBord ].end())
6149 nPrev[ iBord ] = *nIt[ iBord ];
6150 nIt[ iBord ]++; i[ iBord ]++;
6154 while ( nIt[0] != nSide[0].end() && nIt[1] != nSide[1].end());
6156 // perform insertion of nodes into elements
6158 for (insertMapIt = insertMap.begin();
6159 insertMapIt != insertMap.end();
6162 const SMDS_MeshElement* elem = (*insertMapIt).first;
6163 list<const SMDS_MeshNode*> & nodeList = (*insertMapIt).second;
6164 const SMDS_MeshNode* n1 = nodeList.front(); nodeList.pop_front();
6165 const SMDS_MeshNode* n2 = nodeList.front(); nodeList.pop_front();
6167 InsertNodesIntoLink( elem, n1, n2, nodeList, toCreatePolygons );
6169 if ( !theSideIsFreeBorder ) {
6170 // look for and insert nodes into the faces adjacent to elem
6172 const SMDS_MeshElement* adjElem = findAdjacentFace( n1, n2, elem );
6174 InsertNodesIntoLink( adjElem, n1, n2, nodeList, toCreatePolygons );
6179 if (toCreatePolyedrs) {
6180 // perform insertion into the links of adjacent volumes
6181 UpdateVolumes(n1, n2, nodeList);
6187 } // end: insert new nodes
6189 MergeNodes ( nodeGroupsToMerge );
6194 //=======================================================================
6195 //function : InsertNodesIntoLink
6196 //purpose : insert theNodesToInsert into theFace between theBetweenNode1
6197 // and theBetweenNode2 and split theElement
6198 //=======================================================================
6200 void SMESH_MeshEditor::InsertNodesIntoLink(const SMDS_MeshElement* theFace,
6201 const SMDS_MeshNode* theBetweenNode1,
6202 const SMDS_MeshNode* theBetweenNode2,
6203 list<const SMDS_MeshNode*>& theNodesToInsert,
6204 const bool toCreatePoly)
6206 if ( theFace->GetType() != SMDSAbs_Face ) return;
6208 // find indices of 2 link nodes and of the rest nodes
6209 int iNode = 0, il1, il2, i3, i4;
6210 il1 = il2 = i3 = i4 = -1;
6211 //const SMDS_MeshNode* nodes[ theFace->NbNodes() ];
6212 vector<const SMDS_MeshNode*> nodes( theFace->NbNodes() );
6214 if(theFace->IsQuadratic()) {
6215 const SMDS_QuadraticFaceOfNodes* F =
6216 static_cast<const SMDS_QuadraticFaceOfNodes*>(theFace);
6217 // use special nodes iterator
6218 SMDS_NodeIteratorPtr anIter = F->interlacedNodesIterator();
6219 while( anIter->more() ) {
6220 const SMDS_MeshNode* n = anIter->next();
6221 if ( n == theBetweenNode1 )
6223 else if ( n == theBetweenNode2 )
6229 nodes[ iNode++ ] = n;
6233 SMDS_ElemIteratorPtr nodeIt = theFace->nodesIterator();
6234 while ( nodeIt->more() ) {
6235 const SMDS_MeshNode* n = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
6236 if ( n == theBetweenNode1 )
6238 else if ( n == theBetweenNode2 )
6244 nodes[ iNode++ ] = n;
6247 if ( il1 < 0 || il2 < 0 || i3 < 0 )
6250 // arrange link nodes to go one after another regarding the face orientation
6251 bool reverse = ( Abs( il2 - il1 ) == 1 ? il2 < il1 : il1 < il2 );
6252 list<const SMDS_MeshNode *> aNodesToInsert = theNodesToInsert;
6257 aNodesToInsert.reverse();
6259 // check that not link nodes of a quadrangles are in good order
6260 int nbFaceNodes = theFace->NbNodes();
6261 if ( nbFaceNodes == 4 && i4 - i3 != 1 ) {
6267 if (toCreatePoly || theFace->IsPoly()) {
6270 vector<const SMDS_MeshNode *> poly_nodes (nbFaceNodes + aNodesToInsert.size());
6272 // add nodes of face up to first node of link
6275 if(theFace->IsQuadratic()) {
6276 const SMDS_QuadraticFaceOfNodes* F =
6277 static_cast<const SMDS_QuadraticFaceOfNodes*>(theFace);
6278 // use special nodes iterator
6279 SMDS_NodeIteratorPtr anIter = F->interlacedNodesIterator();
6280 while( anIter->more() && !isFLN ) {
6281 const SMDS_MeshNode* n = anIter->next();
6282 poly_nodes[iNode++] = n;
6283 if (n == nodes[il1]) {
6287 // add nodes to insert
6288 list<const SMDS_MeshNode*>::iterator nIt = aNodesToInsert.begin();
6289 for (; nIt != aNodesToInsert.end(); nIt++) {
6290 poly_nodes[iNode++] = *nIt;
6292 // add nodes of face starting from last node of link
6293 while ( anIter->more() ) {
6294 poly_nodes[iNode++] = anIter->next();
6298 SMDS_ElemIteratorPtr nodeIt = theFace->nodesIterator();
6299 while ( nodeIt->more() && !isFLN ) {
6300 const SMDS_MeshNode* n = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
6301 poly_nodes[iNode++] = n;
6302 if (n == nodes[il1]) {
6306 // add nodes to insert
6307 list<const SMDS_MeshNode*>::iterator nIt = aNodesToInsert.begin();
6308 for (; nIt != aNodesToInsert.end(); nIt++) {
6309 poly_nodes[iNode++] = *nIt;
6311 // add nodes of face starting from last node of link
6312 while ( nodeIt->more() ) {
6313 const SMDS_MeshNode* n = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
6314 poly_nodes[iNode++] = n;
6318 // edit or replace the face
6319 SMESHDS_Mesh *aMesh = GetMeshDS();
6321 if (theFace->IsPoly()) {
6322 aMesh->ChangePolygonNodes(theFace, poly_nodes);
6325 int aShapeId = FindShape( theFace );
6327 SMDS_MeshElement* newElem = aMesh->AddPolygonalFace(poly_nodes);
6328 myLastCreatedElems.Append(newElem);
6329 if ( aShapeId && newElem )
6330 aMesh->SetMeshElementOnShape( newElem, aShapeId );
6332 aMesh->RemoveElement(theFace);
6337 if( !theFace->IsQuadratic() ) {
6339 // put aNodesToInsert between theBetweenNode1 and theBetweenNode2
6340 int nbLinkNodes = 2 + aNodesToInsert.size();
6341 //const SMDS_MeshNode* linkNodes[ nbLinkNodes ];
6342 vector<const SMDS_MeshNode*> linkNodes( nbLinkNodes );
6343 linkNodes[ 0 ] = nodes[ il1 ];
6344 linkNodes[ nbLinkNodes - 1 ] = nodes[ il2 ];
6345 list<const SMDS_MeshNode*>::iterator nIt = aNodesToInsert.begin();
6346 for ( iNode = 1; nIt != aNodesToInsert.end(); nIt++ ) {
6347 linkNodes[ iNode++ ] = *nIt;
6349 // decide how to split a quadrangle: compare possible variants
6350 // and choose which of splits to be a quadrangle
6351 int i1, i2, iSplit, nbSplits = nbLinkNodes - 1, iBestQuad;
6352 if ( nbFaceNodes == 3 ) {
6353 iBestQuad = nbSplits;
6356 else if ( nbFaceNodes == 4 ) {
6357 SMESH::Controls::NumericalFunctorPtr aCrit( new SMESH::Controls::AspectRatio);
6358 double aBestRate = DBL_MAX;
6359 for ( int iQuad = 0; iQuad < nbSplits; iQuad++ ) {
6361 double aBadRate = 0;
6362 // evaluate elements quality
6363 for ( iSplit = 0; iSplit < nbSplits; iSplit++ ) {
6364 if ( iSplit == iQuad ) {
6365 SMDS_FaceOfNodes quad (linkNodes[ i1++ ],
6369 aBadRate += getBadRate( &quad, aCrit );
6372 SMDS_FaceOfNodes tria (linkNodes[ i1++ ],
6374 nodes[ iSplit < iQuad ? i4 : i3 ]);
6375 aBadRate += getBadRate( &tria, aCrit );
6379 if ( aBadRate < aBestRate ) {
6381 aBestRate = aBadRate;
6386 // create new elements
6387 SMESHDS_Mesh *aMesh = GetMeshDS();
6388 int aShapeId = FindShape( theFace );
6391 for ( iSplit = 0; iSplit < nbSplits - 1; iSplit++ ) {
6392 SMDS_MeshElement* newElem = 0;
6393 if ( iSplit == iBestQuad )
6394 newElem = aMesh->AddFace (linkNodes[ i1++ ],
6399 newElem = aMesh->AddFace (linkNodes[ i1++ ],
6401 nodes[ iSplit < iBestQuad ? i4 : i3 ]);
6402 myLastCreatedElems.Append(newElem);
6403 if ( aShapeId && newElem )
6404 aMesh->SetMeshElementOnShape( newElem, aShapeId );
6407 // change nodes of theFace
6408 const SMDS_MeshNode* newNodes[ 4 ];
6409 newNodes[ 0 ] = linkNodes[ i1 ];
6410 newNodes[ 1 ] = linkNodes[ i2 ];
6411 newNodes[ 2 ] = nodes[ iSplit >= iBestQuad ? i3 : i4 ];
6412 newNodes[ 3 ] = nodes[ i4 ];
6413 aMesh->ChangeElementNodes( theFace, newNodes, iSplit == iBestQuad ? 4 : 3 );
6414 } // end if(!theFace->IsQuadratic())
6415 else { // theFace is quadratic
6416 // we have to split theFace on simple triangles and one simple quadrangle
6418 int nbshift = tmp*2;
6419 // shift nodes in nodes[] by nbshift
6421 for(i=0; i<nbshift; i++) {
6422 const SMDS_MeshNode* n = nodes[0];
6423 for(j=0; j<nbFaceNodes-1; j++) {
6424 nodes[j] = nodes[j+1];
6426 nodes[nbFaceNodes-1] = n;
6428 il1 = il1 - nbshift;
6429 // now have to insert nodes between n0 and n1 or n1 and n2 (see below)
6430 // n0 n1 n2 n0 n1 n2
6431 // +-----+-----+ +-----+-----+
6440 // create new elements
6441 SMESHDS_Mesh *aMesh = GetMeshDS();
6442 int aShapeId = FindShape( theFace );
6445 if(nbFaceNodes==6) { // quadratic triangle
6446 SMDS_MeshElement* newElem =
6447 aMesh->AddFace(nodes[3],nodes[4],nodes[5]);
6448 myLastCreatedElems.Append(newElem);
6449 if ( aShapeId && newElem )
6450 aMesh->SetMeshElementOnShape( newElem, aShapeId );
6451 if(theFace->IsMediumNode(nodes[il1])) {
6452 // create quadrangle
6453 newElem = aMesh->AddFace(nodes[0],nodes[1],nodes[3],nodes[5]);
6454 myLastCreatedElems.Append(newElem);
6455 if ( aShapeId && newElem )
6456 aMesh->SetMeshElementOnShape( newElem, aShapeId );
6462 // create quadrangle
6463 newElem = aMesh->AddFace(nodes[1],nodes[2],nodes[3],nodes[5]);
6464 myLastCreatedElems.Append(newElem);
6465 if ( aShapeId && newElem )
6466 aMesh->SetMeshElementOnShape( newElem, aShapeId );
6472 else { // nbFaceNodes==8 - quadratic quadrangle
6473 SMDS_MeshElement* newElem =
6474 aMesh->AddFace(nodes[3],nodes[4],nodes[5]);
6475 myLastCreatedElems.Append(newElem);
6476 if ( aShapeId && newElem )
6477 aMesh->SetMeshElementOnShape( newElem, aShapeId );
6478 newElem = aMesh->AddFace(nodes[5],nodes[6],nodes[7]);
6479 myLastCreatedElems.Append(newElem);
6480 if ( aShapeId && newElem )
6481 aMesh->SetMeshElementOnShape( newElem, aShapeId );
6482 newElem = aMesh->AddFace(nodes[5],nodes[7],nodes[3]);
6483 myLastCreatedElems.Append(newElem);
6484 if ( aShapeId && newElem )
6485 aMesh->SetMeshElementOnShape( newElem, aShapeId );
6486 if(theFace->IsMediumNode(nodes[il1])) {
6487 // create quadrangle
6488 newElem = aMesh->AddFace(nodes[0],nodes[1],nodes[3],nodes[7]);
6489 myLastCreatedElems.Append(newElem);
6490 if ( aShapeId && newElem )
6491 aMesh->SetMeshElementOnShape( newElem, aShapeId );
6497 // create quadrangle
6498 newElem = aMesh->AddFace(nodes[1],nodes[2],nodes[3],nodes[7]);
6499 myLastCreatedElems.Append(newElem);
6500 if ( aShapeId && newElem )
6501 aMesh->SetMeshElementOnShape( newElem, aShapeId );
6507 // create needed triangles using n1,n2,n3 and inserted nodes
6508 int nbn = 2 + aNodesToInsert.size();
6509 //const SMDS_MeshNode* aNodes[nbn];
6510 vector<const SMDS_MeshNode*> aNodes(nbn);
6511 aNodes[0] = nodes[n1];
6512 aNodes[nbn-1] = nodes[n2];
6513 list<const SMDS_MeshNode*>::iterator nIt = aNodesToInsert.begin();
6514 for ( iNode = 1; nIt != aNodesToInsert.end(); nIt++ ) {
6515 aNodes[iNode++] = *nIt;
6517 for(i=1; i<nbn; i++) {
6518 SMDS_MeshElement* newElem =
6519 aMesh->AddFace(aNodes[i-1],aNodes[i],nodes[n3]);
6520 myLastCreatedElems.Append(newElem);
6521 if ( aShapeId && newElem )
6522 aMesh->SetMeshElementOnShape( newElem, aShapeId );
6524 // remove old quadratic face
6525 aMesh->RemoveElement(theFace);
6529 //=======================================================================
6530 //function : UpdateVolumes
6532 //=======================================================================
6533 void SMESH_MeshEditor::UpdateVolumes (const SMDS_MeshNode* theBetweenNode1,
6534 const SMDS_MeshNode* theBetweenNode2,
6535 list<const SMDS_MeshNode*>& theNodesToInsert)
6537 myLastCreatedElems.Clear();
6538 myLastCreatedNodes.Clear();
6540 SMDS_ElemIteratorPtr invElemIt = theBetweenNode1->GetInverseElementIterator(SMDSAbs_Volume);
6541 while (invElemIt->more()) { // loop on inverse elements of theBetweenNode1
6542 const SMDS_MeshElement* elem = invElemIt->next();
6544 // check, if current volume has link theBetweenNode1 - theBetweenNode2
6545 SMDS_VolumeTool aVolume (elem);
6546 if (!aVolume.IsLinked(theBetweenNode1, theBetweenNode2))
6549 // insert new nodes in all faces of the volume, sharing link theBetweenNode1 - theBetweenNode2
6550 int iface, nbFaces = aVolume.NbFaces();
6551 vector<const SMDS_MeshNode *> poly_nodes;
6552 vector<int> quantities (nbFaces);
6554 for (iface = 0; iface < nbFaces; iface++) {
6555 int nbFaceNodes = aVolume.NbFaceNodes(iface), nbInserted = 0;
6556 // faceNodes will contain (nbFaceNodes + 1) nodes, last = first
6557 const SMDS_MeshNode** faceNodes = aVolume.GetFaceNodes(iface);
6559 for (int inode = 0; inode < nbFaceNodes; inode++) {
6560 poly_nodes.push_back(faceNodes[inode]);
6562 if (nbInserted == 0) {
6563 if (faceNodes[inode] == theBetweenNode1) {
6564 if (faceNodes[inode + 1] == theBetweenNode2) {
6565 nbInserted = theNodesToInsert.size();
6567 // add nodes to insert
6568 list<const SMDS_MeshNode*>::iterator nIt = theNodesToInsert.begin();
6569 for (; nIt != theNodesToInsert.end(); nIt++) {
6570 poly_nodes.push_back(*nIt);
6574 else if (faceNodes[inode] == theBetweenNode2) {
6575 if (faceNodes[inode + 1] == theBetweenNode1) {
6576 nbInserted = theNodesToInsert.size();
6578 // add nodes to insert in reversed order
6579 list<const SMDS_MeshNode*>::iterator nIt = theNodesToInsert.end();
6581 for (; nIt != theNodesToInsert.begin(); nIt--) {
6582 poly_nodes.push_back(*nIt);
6584 poly_nodes.push_back(*nIt);
6591 quantities[iface] = nbFaceNodes + nbInserted;
6594 // Replace or update the volume
6595 SMESHDS_Mesh *aMesh = GetMeshDS();
6597 if (elem->IsPoly()) {
6598 aMesh->ChangePolyhedronNodes(elem, poly_nodes, quantities);
6602 int aShapeId = FindShape( elem );
6604 SMDS_MeshElement* newElem =
6605 aMesh->AddPolyhedralVolume(poly_nodes, quantities);
6606 myLastCreatedElems.Append(newElem);
6607 if (aShapeId && newElem)
6608 aMesh->SetMeshElementOnShape(newElem, aShapeId);
6610 aMesh->RemoveElement(elem);
6615 //=======================================================================
6617 * \brief Convert elements contained in a submesh to quadratic
6618 * \retval int - nb of checked elements
6620 //=======================================================================
6622 int SMESH_MeshEditor::convertElemToQuadratic(SMESHDS_SubMesh * theSm,
6623 SMESH_MesherHelper& theHelper,
6624 const bool theForce3d)
6627 if( !theSm ) return nbElem;
6629 const bool notFromGroups = false;
6630 SMDS_ElemIteratorPtr ElemItr = theSm->GetElements();
6631 while(ElemItr->more())
6634 const SMDS_MeshElement* elem = ElemItr->next();
6635 if( !elem || elem->IsQuadratic() ) continue;
6637 int id = elem->GetID();
6638 int nbNodes = elem->NbNodes();
6639 vector<const SMDS_MeshNode *> aNds (nbNodes);
6641 for(int i = 0; i < nbNodes; i++)
6643 aNds[i] = elem->GetNode(i);
6645 SMDSAbs_ElementType aType = elem->GetType();
6647 GetMeshDS()->RemoveFreeElement(elem, theSm, notFromGroups);
6649 const SMDS_MeshElement* NewElem = 0;
6655 NewElem = theHelper.AddEdge(aNds[0], aNds[1], id, theForce3d);
6663 NewElem = theHelper.AddFace(aNds[0], aNds[1], aNds[2], id, theForce3d);
6666 NewElem = theHelper.AddFace(aNds[0], aNds[1], aNds[2], aNds[3], id, theForce3d);
6673 case SMDSAbs_Volume :
6678 NewElem = theHelper.AddVolume(aNds[0], aNds[1], aNds[2], aNds[3], id, theForce3d);
6681 NewElem = theHelper.AddVolume(aNds[0], aNds[1], aNds[2], aNds[3], aNds[4], aNds[5], id, theForce3d);
6684 NewElem = theHelper.AddVolume(aNds[0], aNds[1], aNds[2], aNds[3],
6685 aNds[4], aNds[5], aNds[6], aNds[7], id, theForce3d);
6695 ReplaceElemInGroups( elem, NewElem, GetMeshDS());
6697 theSm->AddElement( NewElem );
6702 //=======================================================================
6703 //function : ConvertToQuadratic
6705 //=======================================================================
6706 void SMESH_MeshEditor::ConvertToQuadratic(const bool theForce3d)
6708 SMESHDS_Mesh* meshDS = GetMeshDS();
6710 SMESH_MesherHelper aHelper(*myMesh);
6711 aHelper.SetIsQuadratic( true );
6712 const bool notFromGroups = false;
6714 int nbCheckedElems = 0;
6715 if ( myMesh->HasShapeToMesh() )
6717 if ( SMESH_subMesh *aSubMesh = myMesh->GetSubMeshContaining(myMesh->GetShapeToMesh()))
6719 SMESH_subMeshIteratorPtr smIt = aSubMesh->getDependsOnIterator(true,false);
6720 while ( smIt->more() ) {
6721 SMESH_subMesh* sm = smIt->next();
6722 if ( SMESHDS_SubMesh *smDS = sm->GetSubMeshDS() ) {
6723 aHelper.SetSubShape( sm->GetSubShape() );
6724 if ( !theForce3d) aHelper.SetCheckNodePosition(true);
6725 nbCheckedElems += convertElemToQuadratic(smDS, aHelper, theForce3d);
6730 int totalNbElems = meshDS->NbEdges() + meshDS->NbFaces() + meshDS->NbVolumes();
6731 if ( nbCheckedElems < totalNbElems ) // not all elements are in submeshes
6733 SMESHDS_SubMesh *smDS = 0;
6734 SMDS_EdgeIteratorPtr aEdgeItr = meshDS->edgesIterator();
6735 while(aEdgeItr->more())
6737 const SMDS_MeshEdge* edge = aEdgeItr->next();
6738 if(edge && !edge->IsQuadratic())
6740 int id = edge->GetID();
6741 const SMDS_MeshNode* n1 = edge->GetNode(0);
6742 const SMDS_MeshNode* n2 = edge->GetNode(1);
6744 meshDS->RemoveFreeElement(edge, smDS, notFromGroups);
6746 const SMDS_MeshEdge* NewEdge = aHelper.AddEdge(n1, n2, id, theForce3d);
6747 ReplaceElemInGroups( edge, NewEdge, GetMeshDS());
6750 SMDS_FaceIteratorPtr aFaceItr = meshDS->facesIterator();
6751 while(aFaceItr->more())
6753 const SMDS_MeshFace* face = aFaceItr->next();
6754 if(!face || face->IsQuadratic() ) continue;
6756 int id = face->GetID();
6757 int nbNodes = face->NbNodes();
6758 vector<const SMDS_MeshNode *> aNds (nbNodes);
6760 for(int i = 0; i < nbNodes; i++)
6762 aNds[i] = face->GetNode(i);
6765 meshDS->RemoveFreeElement(face, smDS, notFromGroups);
6767 SMDS_MeshFace * NewFace = 0;
6771 NewFace = aHelper.AddFace(aNds[0], aNds[1], aNds[2], id, theForce3d);
6774 NewFace = aHelper.AddFace(aNds[0], aNds[1], aNds[2], aNds[3], id, theForce3d);
6779 ReplaceElemInGroups( face, NewFace, GetMeshDS());
6781 SMDS_VolumeIteratorPtr aVolumeItr = meshDS->volumesIterator();
6782 while(aVolumeItr->more())
6784 const SMDS_MeshVolume* volume = aVolumeItr->next();
6785 if(!volume || volume->IsQuadratic() ) continue;
6787 int id = volume->GetID();
6788 int nbNodes = volume->NbNodes();
6789 vector<const SMDS_MeshNode *> aNds (nbNodes);
6791 for(int i = 0; i < nbNodes; i++)
6793 aNds[i] = volume->GetNode(i);
6796 meshDS->RemoveFreeElement(volume, smDS, notFromGroups);
6798 SMDS_MeshVolume * NewVolume = 0;
6802 NewVolume = aHelper.AddVolume(aNds[0], aNds[1], aNds[2],
6803 aNds[3], id, theForce3d );
6806 NewVolume = aHelper.AddVolume(aNds[0], aNds[1], aNds[2],
6807 aNds[3], aNds[4], aNds[5], id, theForce3d);
6810 NewVolume = aHelper.AddVolume(aNds[0], aNds[1], aNds[2], aNds[3],
6811 aNds[4], aNds[5], aNds[6], aNds[7], id, theForce3d);
6816 ReplaceElemInGroups(volume, NewVolume, meshDS);
6819 if ( !theForce3d ) {
6820 aHelper.SetSubShape(0); // apply to the whole mesh
6821 aHelper.FixQuadraticElements();
6825 //=======================================================================
6827 * \brief Convert quadratic elements to linear ones and remove quadratic nodes
6828 * \retval int - nb of checked elements
6830 //=======================================================================
6832 int SMESH_MeshEditor::removeQuadElem(SMESHDS_SubMesh * theSm,
6833 SMDS_ElemIteratorPtr theItr,
6834 const int theShapeID)
6837 SMESHDS_Mesh* meshDS = GetMeshDS();
6838 const bool notFromGroups = false;
6840 while( theItr->more() )
6842 const SMDS_MeshElement* elem = theItr->next();
6844 if( elem && elem->IsQuadratic())
6846 int id = elem->GetID();
6847 int nbNodes = elem->NbNodes();
6848 vector<const SMDS_MeshNode *> aNds, mediumNodes;
6849 aNds.reserve( nbNodes );
6850 mediumNodes.reserve( nbNodes );
6852 for(int i = 0; i < nbNodes; i++)
6854 const SMDS_MeshNode* n = elem->GetNode(i);
6856 if( elem->IsMediumNode( n ) )
6857 mediumNodes.push_back( n );
6859 aNds.push_back( n );
6861 if( aNds.empty() ) continue;
6862 SMDSAbs_ElementType aType = elem->GetType();
6864 //remove old quadratic element
6865 meshDS->RemoveFreeElement( elem, theSm, notFromGroups );
6867 SMDS_MeshElement * NewElem = AddElement( aNds, aType, false, id );
6868 ReplaceElemInGroups(elem, NewElem, meshDS);
6869 if( theSm && NewElem )
6870 theSm->AddElement( NewElem );
6872 // remove medium nodes
6873 vector<const SMDS_MeshNode*>::iterator nIt = mediumNodes.begin();
6874 for ( ; nIt != mediumNodes.end(); ++nIt ) {
6875 const SMDS_MeshNode* n = *nIt;
6876 if ( n->NbInverseElements() == 0 ) {
6877 if ( n->GetPosition()->GetShapeId() != theShapeID )
6878 meshDS->RemoveFreeNode( n, meshDS->MeshElements
6879 ( n->GetPosition()->GetShapeId() ));
6881 meshDS->RemoveFreeNode( n, theSm );
6889 //=======================================================================
6890 //function : ConvertFromQuadratic
6892 //=======================================================================
6893 bool SMESH_MeshEditor::ConvertFromQuadratic()
6895 int nbCheckedElems = 0;
6896 if ( myMesh->HasShapeToMesh() )
6898 if ( SMESH_subMesh *aSubMesh = myMesh->GetSubMeshContaining(myMesh->GetShapeToMesh()))
6900 SMESH_subMeshIteratorPtr smIt = aSubMesh->getDependsOnIterator(true,false);
6901 while ( smIt->more() ) {
6902 SMESH_subMesh* sm = smIt->next();
6903 if ( SMESHDS_SubMesh *smDS = sm->GetSubMeshDS() )
6904 nbCheckedElems += removeQuadElem( smDS, smDS->GetElements(), sm->GetId() );
6910 GetMeshDS()->NbEdges() + GetMeshDS()->NbFaces() + GetMeshDS()->NbVolumes();
6911 if ( nbCheckedElems < totalNbElems ) // not all elements are in submeshes
6913 SMESHDS_SubMesh *aSM = 0;
6914 removeQuadElem( aSM, GetMeshDS()->elementsIterator(), 0 );
6920 //=======================================================================
6921 //function : SewSideElements
6923 //=======================================================================
6925 SMESH_MeshEditor::Sew_Error
6926 SMESH_MeshEditor::SewSideElements (TIDSortedElemSet& theSide1,
6927 TIDSortedElemSet& theSide2,
6928 const SMDS_MeshNode* theFirstNode1,
6929 const SMDS_MeshNode* theFirstNode2,
6930 const SMDS_MeshNode* theSecondNode1,
6931 const SMDS_MeshNode* theSecondNode2)
6933 myLastCreatedElems.Clear();
6934 myLastCreatedNodes.Clear();
6936 MESSAGE ("::::SewSideElements()");
6937 if ( theSide1.size() != theSide2.size() )
6938 return SEW_DIFF_NB_OF_ELEMENTS;
6940 Sew_Error aResult = SEW_OK;
6942 // 1. Build set of faces representing each side
6943 // 2. Find which nodes of the side 1 to merge with ones on the side 2
6944 // 3. Replace nodes in elements of the side 1 and remove replaced nodes
6946 // =======================================================================
6947 // 1. Build set of faces representing each side:
6948 // =======================================================================
6949 // a. build set of nodes belonging to faces
6950 // b. complete set of faces: find missing fices whose nodes are in set of nodes
6951 // c. create temporary faces representing side of volumes if correspondent
6952 // face does not exist
6954 SMESHDS_Mesh* aMesh = GetMeshDS();
6955 SMDS_Mesh aTmpFacesMesh;
6956 set<const SMDS_MeshElement*> faceSet1, faceSet2;
6957 set<const SMDS_MeshElement*> volSet1, volSet2;
6958 set<const SMDS_MeshNode*> nodeSet1, nodeSet2;
6959 set<const SMDS_MeshElement*> * faceSetPtr[] = { &faceSet1, &faceSet2 };
6960 set<const SMDS_MeshElement*> * volSetPtr[] = { &volSet1, &volSet2 };
6961 set<const SMDS_MeshNode*> * nodeSetPtr[] = { &nodeSet1, &nodeSet2 };
6962 TIDSortedElemSet * elemSetPtr[] = { &theSide1, &theSide2 };
6963 int iSide, iFace, iNode;
6965 for ( iSide = 0; iSide < 2; iSide++ ) {
6966 set<const SMDS_MeshNode*> * nodeSet = nodeSetPtr[ iSide ];
6967 TIDSortedElemSet * elemSet = elemSetPtr[ iSide ];
6968 set<const SMDS_MeshElement*> * faceSet = faceSetPtr[ iSide ];
6969 set<const SMDS_MeshElement*> * volSet = volSetPtr [ iSide ];
6970 set<const SMDS_MeshElement*>::iterator vIt;
6971 TIDSortedElemSet::iterator eIt;
6972 set<const SMDS_MeshNode*>::iterator nIt;
6974 // check that given nodes belong to given elements
6975 const SMDS_MeshNode* n1 = ( iSide == 0 ) ? theFirstNode1 : theFirstNode2;
6976 const SMDS_MeshNode* n2 = ( iSide == 0 ) ? theSecondNode1 : theSecondNode2;
6977 int firstIndex = -1, secondIndex = -1;
6978 for (eIt = elemSet->begin(); eIt != elemSet->end(); eIt++ ) {
6979 const SMDS_MeshElement* elem = *eIt;
6980 if ( firstIndex < 0 ) firstIndex = elem->GetNodeIndex( n1 );
6981 if ( secondIndex < 0 ) secondIndex = elem->GetNodeIndex( n2 );
6982 if ( firstIndex > -1 && secondIndex > -1 ) break;
6984 if ( firstIndex < 0 || secondIndex < 0 ) {
6985 // we can simply return until temporary faces created
6986 return (iSide == 0 ) ? SEW_BAD_SIDE1_NODES : SEW_BAD_SIDE2_NODES;
6989 // -----------------------------------------------------------
6990 // 1a. Collect nodes of existing faces
6991 // and build set of face nodes in order to detect missing
6992 // faces corresponing to sides of volumes
6993 // -----------------------------------------------------------
6995 set< set <const SMDS_MeshNode*> > setOfFaceNodeSet;
6997 // loop on the given element of a side
6998 for (eIt = elemSet->begin(); eIt != elemSet->end(); eIt++ ) {
6999 //const SMDS_MeshElement* elem = *eIt;
7000 const SMDS_MeshElement* elem = *eIt;
7001 if ( elem->GetType() == SMDSAbs_Face ) {
7002 faceSet->insert( elem );
7003 set <const SMDS_MeshNode*> faceNodeSet;
7004 SMDS_ElemIteratorPtr nodeIt = elem->nodesIterator();
7005 while ( nodeIt->more() ) {
7006 const SMDS_MeshNode* n = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
7007 nodeSet->insert( n );
7008 faceNodeSet.insert( n );
7010 setOfFaceNodeSet.insert( faceNodeSet );
7012 else if ( elem->GetType() == SMDSAbs_Volume )
7013 volSet->insert( elem );
7015 // ------------------------------------------------------------------------------
7016 // 1b. Complete set of faces: find missing fices whose nodes are in set of nodes
7017 // ------------------------------------------------------------------------------
7019 for ( nIt = nodeSet->begin(); nIt != nodeSet->end(); nIt++ ) { // loop on nodes of iSide
7020 SMDS_ElemIteratorPtr fIt = (*nIt)->GetInverseElementIterator(SMDSAbs_Face);
7021 while ( fIt->more() ) { // loop on faces sharing a node
7022 const SMDS_MeshElement* f = fIt->next();
7023 if ( faceSet->find( f ) == faceSet->end() ) {
7024 // check if all nodes are in nodeSet and
7025 // complete setOfFaceNodeSet if they are
7026 set <const SMDS_MeshNode*> faceNodeSet;
7027 SMDS_ElemIteratorPtr nodeIt = f->nodesIterator();
7028 bool allInSet = true;
7029 while ( nodeIt->more() && allInSet ) { // loop on nodes of a face
7030 const SMDS_MeshNode* n = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
7031 if ( nodeSet->find( n ) == nodeSet->end() )
7034 faceNodeSet.insert( n );
7037 faceSet->insert( f );
7038 setOfFaceNodeSet.insert( faceNodeSet );
7044 // -------------------------------------------------------------------------
7045 // 1c. Create temporary faces representing sides of volumes if correspondent
7046 // face does not exist
7047 // -------------------------------------------------------------------------
7049 if ( !volSet->empty() ) {
7050 //int nodeSetSize = nodeSet->size();
7052 // loop on given volumes
7053 for ( vIt = volSet->begin(); vIt != volSet->end(); vIt++ ) {
7054 SMDS_VolumeTool vol (*vIt);
7055 // loop on volume faces: find free faces
7056 // --------------------------------------
7057 list<const SMDS_MeshElement* > freeFaceList;
7058 for ( iFace = 0; iFace < vol.NbFaces(); iFace++ ) {
7059 if ( !vol.IsFreeFace( iFace ))
7061 // check if there is already a face with same nodes in a face set
7062 const SMDS_MeshElement* aFreeFace = 0;
7063 const SMDS_MeshNode** fNodes = vol.GetFaceNodes( iFace );
7064 int nbNodes = vol.NbFaceNodes( iFace );
7065 set <const SMDS_MeshNode*> faceNodeSet;
7066 vol.GetFaceNodes( iFace, faceNodeSet );
7067 bool isNewFace = setOfFaceNodeSet.insert( faceNodeSet ).second;
7069 // no such a face is given but it still can exist, check it
7070 if ( nbNodes == 3 ) {
7071 aFreeFace = aMesh->FindFace( fNodes[0],fNodes[1],fNodes[2] );
7073 else if ( nbNodes == 4 ) {
7074 aFreeFace = aMesh->FindFace( fNodes[0],fNodes[1],fNodes[2],fNodes[3] );
7077 vector<const SMDS_MeshNode *> poly_nodes ( fNodes, & fNodes[nbNodes]);
7078 aFreeFace = aMesh->FindFace(poly_nodes);
7082 // create a temporary face
7083 if ( nbNodes == 3 ) {
7084 aFreeFace = aTmpFacesMesh.AddFace( fNodes[0],fNodes[1],fNodes[2] );
7086 else if ( nbNodes == 4 ) {
7087 aFreeFace = aTmpFacesMesh.AddFace( fNodes[0],fNodes[1],fNodes[2],fNodes[3] );
7090 vector<const SMDS_MeshNode *> poly_nodes ( fNodes, & fNodes[nbNodes]);
7091 aFreeFace = aTmpFacesMesh.AddPolygonalFace(poly_nodes);
7095 freeFaceList.push_back( aFreeFace );
7097 } // loop on faces of a volume
7099 // choose one of several free faces
7100 // --------------------------------------
7101 if ( freeFaceList.size() > 1 ) {
7102 // choose a face having max nb of nodes shared by other elems of a side
7103 int maxNbNodes = -1/*, nbExcludedFaces = 0*/;
7104 list<const SMDS_MeshElement* >::iterator fIt = freeFaceList.begin();
7105 while ( fIt != freeFaceList.end() ) { // loop on free faces
7106 int nbSharedNodes = 0;
7107 SMDS_ElemIteratorPtr nodeIt = (*fIt)->nodesIterator();
7108 while ( nodeIt->more() ) { // loop on free face nodes
7109 const SMDS_MeshNode* n =
7110 static_cast<const SMDS_MeshNode*>( nodeIt->next() );
7111 SMDS_ElemIteratorPtr invElemIt = n->GetInverseElementIterator();
7112 while ( invElemIt->more() ) {
7113 const SMDS_MeshElement* e = invElemIt->next();
7114 if ( faceSet->find( e ) != faceSet->end() )
7116 if ( elemSet->find( e ) != elemSet->end() )
7120 if ( nbSharedNodes >= maxNbNodes ) {
7121 maxNbNodes = nbSharedNodes;
7125 freeFaceList.erase( fIt++ ); // here fIt++ occures before erase
7127 if ( freeFaceList.size() > 1 )
7129 // could not choose one face, use another way
7130 // choose a face most close to the bary center of the opposite side
7131 gp_XYZ aBC( 0., 0., 0. );
7132 set <const SMDS_MeshNode*> addedNodes;
7133 TIDSortedElemSet * elemSet2 = elemSetPtr[ 1 - iSide ];
7134 eIt = elemSet2->begin();
7135 for ( eIt = elemSet2->begin(); eIt != elemSet2->end(); eIt++ ) {
7136 SMDS_ElemIteratorPtr nodeIt = (*eIt)->nodesIterator();
7137 while ( nodeIt->more() ) { // loop on free face nodes
7138 const SMDS_MeshNode* n =
7139 static_cast<const SMDS_MeshNode*>( nodeIt->next() );
7140 if ( addedNodes.insert( n ).second )
7141 aBC += gp_XYZ( n->X(),n->Y(),n->Z() );
7144 aBC /= addedNodes.size();
7145 double minDist = DBL_MAX;
7146 fIt = freeFaceList.begin();
7147 while ( fIt != freeFaceList.end() ) { // loop on free faces
7149 SMDS_ElemIteratorPtr nodeIt = (*fIt)->nodesIterator();
7150 while ( nodeIt->more() ) { // loop on free face nodes
7151 const SMDS_MeshNode* n =
7152 static_cast<const SMDS_MeshNode*>( nodeIt->next() );
7153 gp_XYZ p( n->X(),n->Y(),n->Z() );
7154 dist += ( aBC - p ).SquareModulus();
7156 if ( dist < minDist ) {
7158 freeFaceList.erase( freeFaceList.begin(), fIt++ );
7161 fIt = freeFaceList.erase( fIt++ );
7164 } // choose one of several free faces of a volume
7166 if ( freeFaceList.size() == 1 ) {
7167 const SMDS_MeshElement* aFreeFace = freeFaceList.front();
7168 faceSet->insert( aFreeFace );
7169 // complete a node set with nodes of a found free face
7170 // for ( iNode = 0; iNode < ; iNode++ )
7171 // nodeSet->insert( fNodes[ iNode ] );
7174 } // loop on volumes of a side
7176 // // complete a set of faces if new nodes in a nodeSet appeared
7177 // // ----------------------------------------------------------
7178 // if ( nodeSetSize != nodeSet->size() ) {
7179 // for ( ; nIt != nodeSet->end(); nIt++ ) { // loop on nodes of iSide
7180 // SMDS_ElemIteratorPtr fIt = (*nIt)->GetInverseElementIterator(SMDSAbs_Face);
7181 // while ( fIt->more() ) { // loop on faces sharing a node
7182 // const SMDS_MeshElement* f = fIt->next();
7183 // if ( faceSet->find( f ) == faceSet->end() ) {
7184 // // check if all nodes are in nodeSet and
7185 // // complete setOfFaceNodeSet if they are
7186 // set <const SMDS_MeshNode*> faceNodeSet;
7187 // SMDS_ElemIteratorPtr nodeIt = f->nodesIterator();
7188 // bool allInSet = true;
7189 // while ( nodeIt->more() && allInSet ) { // loop on nodes of a face
7190 // const SMDS_MeshNode* n = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
7191 // if ( nodeSet->find( n ) == nodeSet->end() )
7192 // allInSet = false;
7194 // faceNodeSet.insert( n );
7196 // if ( allInSet ) {
7197 // faceSet->insert( f );
7198 // setOfFaceNodeSet.insert( faceNodeSet );
7204 } // Create temporary faces, if there are volumes given
7207 if ( faceSet1.size() != faceSet2.size() ) {
7208 // delete temporary faces: they are in reverseElements of actual nodes
7209 SMDS_FaceIteratorPtr tmpFaceIt = aTmpFacesMesh.facesIterator();
7210 while ( tmpFaceIt->more() )
7211 aTmpFacesMesh.RemoveElement( tmpFaceIt->next() );
7212 MESSAGE("Diff nb of faces");
7213 return SEW_TOPO_DIFF_SETS_OF_ELEMENTS;
7216 // ============================================================
7217 // 2. Find nodes to merge:
7218 // bind a node to remove to a node to put instead
7219 // ============================================================
7221 TNodeNodeMap nReplaceMap; // bind a node to remove to a node to put instead
7222 if ( theFirstNode1 != theFirstNode2 )
7223 nReplaceMap.insert( TNodeNodeMap::value_type( theFirstNode1, theFirstNode2 ));
7224 if ( theSecondNode1 != theSecondNode2 )
7225 nReplaceMap.insert( TNodeNodeMap::value_type( theSecondNode1, theSecondNode2 ));
7227 LinkID_Gen aLinkID_Gen( GetMeshDS() );
7228 set< long > linkIdSet; // links to process
7229 linkIdSet.insert( aLinkID_Gen.GetLinkID( theFirstNode1, theSecondNode1 ));
7231 typedef pair< const SMDS_MeshNode*, const SMDS_MeshNode* > NLink;
7232 list< NLink > linkList[2];
7233 linkList[0].push_back( NLink( theFirstNode1, theSecondNode1 ));
7234 linkList[1].push_back( NLink( theFirstNode2, theSecondNode2 ));
7235 // loop on links in linkList; find faces by links and append links
7236 // of the found faces to linkList
7237 list< NLink >::iterator linkIt[] = { linkList[0].begin(), linkList[1].begin() } ;
7238 for ( ; linkIt[0] != linkList[0].end(); linkIt[0]++, linkIt[1]++ ) {
7239 NLink link[] = { *linkIt[0], *linkIt[1] };
7240 long linkID = aLinkID_Gen.GetLinkID( link[0].first, link[0].second );
7241 if ( linkIdSet.find( linkID ) == linkIdSet.end() )
7244 // by links, find faces in the face sets,
7245 // and find indices of link nodes in the found faces;
7246 // in a face set, there is only one or no face sharing a link
7247 // ---------------------------------------------------------------
7249 const SMDS_MeshElement* face[] = { 0, 0 };
7250 //const SMDS_MeshNode* faceNodes[ 2 ][ 5 ];
7251 vector<const SMDS_MeshNode*> fnodes1(9);
7252 vector<const SMDS_MeshNode*> fnodes2(9);
7253 //const SMDS_MeshNode* notLinkNodes[ 2 ][ 2 ] = {{ 0, 0 },{ 0, 0 }} ;
7254 vector<const SMDS_MeshNode*> notLinkNodes1(6);
7255 vector<const SMDS_MeshNode*> notLinkNodes2(6);
7256 int iLinkNode[2][2];
7257 for ( iSide = 0; iSide < 2; iSide++ ) { // loop on 2 sides
7258 const SMDS_MeshNode* n1 = link[iSide].first;
7259 const SMDS_MeshNode* n2 = link[iSide].second;
7260 set<const SMDS_MeshElement*> * faceSet = faceSetPtr[ iSide ];
7261 set< const SMDS_MeshElement* > fMap;
7262 for ( int i = 0; i < 2; i++ ) { // loop on 2 nodes of a link
7263 const SMDS_MeshNode* n = i ? n1 : n2; // a node of a link
7264 SMDS_ElemIteratorPtr fIt = n->GetInverseElementIterator(SMDSAbs_Face);
7265 while ( fIt->more() ) { // loop on faces sharing a node
7266 const SMDS_MeshElement* f = fIt->next();
7267 if (faceSet->find( f ) != faceSet->end() && // f is in face set
7268 ! fMap.insert( f ).second ) // f encounters twice
7270 if ( face[ iSide ] ) {
7271 MESSAGE( "2 faces per link " );
7272 aResult = iSide ? SEW_BAD_SIDE2_NODES : SEW_BAD_SIDE1_NODES;
7276 faceSet->erase( f );
7277 // get face nodes and find ones of a link
7282 fnodes1.resize(f->NbNodes()+1);
7283 notLinkNodes1.resize(f->NbNodes()-2);
7286 fnodes2.resize(f->NbNodes()+1);
7287 notLinkNodes2.resize(f->NbNodes()-2);
7290 if(!f->IsQuadratic()) {
7291 SMDS_ElemIteratorPtr nIt = f->nodesIterator();
7292 while ( nIt->more() ) {
7293 const SMDS_MeshNode* n =
7294 static_cast<const SMDS_MeshNode*>( nIt->next() );
7296 iLinkNode[ iSide ][ 0 ] = iNode;
7298 else if ( n == n2 ) {
7299 iLinkNode[ iSide ][ 1 ] = iNode;
7301 //else if ( notLinkNodes[ iSide ][ 0 ] )
7302 // notLinkNodes[ iSide ][ 1 ] = n;
7304 // notLinkNodes[ iSide ][ 0 ] = n;
7308 notLinkNodes1[nbl] = n;
7309 //notLinkNodes1.push_back(n);
7311 notLinkNodes2[nbl] = n;
7312 //notLinkNodes2.push_back(n);
7314 //faceNodes[ iSide ][ iNode++ ] = n;
7316 fnodes1[iNode++] = n;
7319 fnodes2[iNode++] = n;
7323 else { // f->IsQuadratic()
7324 const SMDS_QuadraticFaceOfNodes* F =
7325 static_cast<const SMDS_QuadraticFaceOfNodes*>(f);
7326 // use special nodes iterator
7327 SMDS_NodeIteratorPtr anIter = F->interlacedNodesIterator();
7328 while ( anIter->more() ) {
7329 const SMDS_MeshNode* n =
7330 static_cast<const SMDS_MeshNode*>( anIter->next() );
7332 iLinkNode[ iSide ][ 0 ] = iNode;
7334 else if ( n == n2 ) {
7335 iLinkNode[ iSide ][ 1 ] = iNode;
7340 notLinkNodes1[nbl] = n;
7343 notLinkNodes2[nbl] = n;
7347 fnodes1[iNode++] = n;
7350 fnodes2[iNode++] = n;
7354 //faceNodes[ iSide ][ iNode ] = faceNodes[ iSide ][ 0 ];
7356 fnodes1[iNode] = fnodes1[0];
7359 fnodes2[iNode] = fnodes1[0];
7366 // check similarity of elements of the sides
7367 if (aResult == SEW_OK && ( face[0] && !face[1] ) || ( !face[0] && face[1] )) {
7368 MESSAGE("Correspondent face not found on side " << ( face[0] ? 1 : 0 ));
7369 if ( nReplaceMap.size() == 2 ) { // faces on input nodes not found
7370 aResult = ( face[0] ? SEW_BAD_SIDE2_NODES : SEW_BAD_SIDE1_NODES );
7373 aResult = SEW_TOPO_DIFF_SETS_OF_ELEMENTS;
7375 break; // do not return because it s necessary to remove tmp faces
7378 // set nodes to merge
7379 // -------------------
7381 if ( face[0] && face[1] ) {
7382 int nbNodes = face[0]->NbNodes();
7383 if ( nbNodes != face[1]->NbNodes() ) {
7384 MESSAGE("Diff nb of face nodes");
7385 aResult = SEW_TOPO_DIFF_SETS_OF_ELEMENTS;
7386 break; // do not return because it s necessary to remove tmp faces
7388 bool reverse[] = { false, false }; // order of notLinkNodes of quadrangle
7389 if ( nbNodes == 3 ) {
7390 //nReplaceMap.insert( TNodeNodeMap::value_type
7391 // ( notLinkNodes[0][0], notLinkNodes[1][0] ));
7392 nReplaceMap.insert( TNodeNodeMap::value_type
7393 ( notLinkNodes1[0], notLinkNodes2[0] ));
7396 for ( iSide = 0; iSide < 2; iSide++ ) { // loop on 2 sides
7397 // analyse link orientation in faces
7398 int i1 = iLinkNode[ iSide ][ 0 ];
7399 int i2 = iLinkNode[ iSide ][ 1 ];
7400 reverse[ iSide ] = Abs( i1 - i2 ) == 1 ? i1 > i2 : i2 > i1;
7401 // if notLinkNodes are the first and the last ones, then
7402 // their order does not correspond to the link orientation
7403 if (( i1 == 1 && i2 == 2 ) ||
7404 ( i1 == 2 && i2 == 1 ))
7405 reverse[ iSide ] = !reverse[ iSide ];
7407 if ( reverse[0] == reverse[1] ) {
7408 //nReplaceMap.insert( TNodeNodeMap::value_type
7409 // ( notLinkNodes[0][0], notLinkNodes[1][0] ));
7410 //nReplaceMap.insert( TNodeNodeMap::value_type
7411 // ( notLinkNodes[0][1], notLinkNodes[1][1] ));
7412 for(int nn=0; nn<nbNodes-2; nn++) {
7413 nReplaceMap.insert( TNodeNodeMap::value_type
7414 ( notLinkNodes1[nn], notLinkNodes2[nn] ));
7418 //nReplaceMap.insert( TNodeNodeMap::value_type
7419 // ( notLinkNodes[0][0], notLinkNodes[1][1] ));
7420 //nReplaceMap.insert( TNodeNodeMap::value_type
7421 // ( notLinkNodes[0][1], notLinkNodes[1][0] ));
7422 for(int nn=0; nn<nbNodes-2; nn++) {
7423 nReplaceMap.insert( TNodeNodeMap::value_type
7424 ( notLinkNodes1[nn], notLinkNodes2[nbNodes-3-nn] ));
7429 // add other links of the faces to linkList
7430 // -----------------------------------------
7432 //const SMDS_MeshNode** nodes = faceNodes[ 0 ];
7433 for ( iNode = 0; iNode < nbNodes; iNode++ ) {
7434 //linkID = aLinkID_Gen.GetLinkID( nodes[iNode], nodes[iNode+1] );
7435 linkID = aLinkID_Gen.GetLinkID( fnodes1[iNode], fnodes1[iNode+1] );
7436 pair< set<long>::iterator, bool > iter_isnew = linkIdSet.insert( linkID );
7437 if ( !iter_isnew.second ) { // already in a set: no need to process
7438 linkIdSet.erase( iter_isnew.first );
7440 else // new in set == encountered for the first time: add
7442 //const SMDS_MeshNode* n1 = nodes[ iNode ];
7443 //const SMDS_MeshNode* n2 = nodes[ iNode + 1];
7444 const SMDS_MeshNode* n1 = fnodes1[ iNode ];
7445 const SMDS_MeshNode* n2 = fnodes1[ iNode + 1];
7446 linkList[0].push_back ( NLink( n1, n2 ));
7447 linkList[1].push_back ( NLink( nReplaceMap[n1], nReplaceMap[n2] ));
7451 } // loop on link lists
7453 if ( aResult == SEW_OK &&
7454 ( linkIt[0] != linkList[0].end() ||
7455 !faceSetPtr[0]->empty() || !faceSetPtr[1]->empty() )) {
7456 MESSAGE( (linkIt[0] != linkList[0].end()) <<" "<< (faceSetPtr[0]->empty()) <<
7457 " " << (faceSetPtr[1]->empty()));
7458 aResult = SEW_TOPO_DIFF_SETS_OF_ELEMENTS;
7461 // ====================================================================
7462 // 3. Replace nodes in elements of the side 1 and remove replaced nodes
7463 // ====================================================================
7465 // delete temporary faces: they are in reverseElements of actual nodes
7466 SMDS_FaceIteratorPtr tmpFaceIt = aTmpFacesMesh.facesIterator();
7467 while ( tmpFaceIt->more() )
7468 aTmpFacesMesh.RemoveElement( tmpFaceIt->next() );
7470 if ( aResult != SEW_OK)
7473 list< int > nodeIDsToRemove/*, elemIDsToRemove*/;
7474 // loop on nodes replacement map
7475 TNodeNodeMap::iterator nReplaceMapIt = nReplaceMap.begin(), nnIt;
7476 for ( ; nReplaceMapIt != nReplaceMap.end(); nReplaceMapIt++ )
7477 if ( (*nReplaceMapIt).first != (*nReplaceMapIt).second ) {
7478 const SMDS_MeshNode* nToRemove = (*nReplaceMapIt).first;
7479 nodeIDsToRemove.push_back( nToRemove->GetID() );
7480 // loop on elements sharing nToRemove
7481 SMDS_ElemIteratorPtr invElemIt = nToRemove->GetInverseElementIterator();
7482 while ( invElemIt->more() ) {
7483 const SMDS_MeshElement* e = invElemIt->next();
7484 // get a new suite of nodes: make replacement
7485 int nbReplaced = 0, i = 0, nbNodes = e->NbNodes();
7486 vector< const SMDS_MeshNode*> nodes( nbNodes );
7487 SMDS_ElemIteratorPtr nIt = e->nodesIterator();
7488 while ( nIt->more() ) {
7489 const SMDS_MeshNode* n =
7490 static_cast<const SMDS_MeshNode*>( nIt->next() );
7491 nnIt = nReplaceMap.find( n );
7492 if ( nnIt != nReplaceMap.end() ) {
7498 // if ( nbReplaced == nbNodes && e->GetType() == SMDSAbs_Face )
7499 // elemIDsToRemove.push_back( e->GetID() );
7502 aMesh->ChangeElementNodes( e, & nodes[0], nbNodes );
7506 Remove( nodeIDsToRemove, true );
7511 //================================================================================
7513 * \brief Find corresponding nodes in two sets of faces
7514 * \param theSide1 - first face set
7515 * \param theSide2 - second first face
7516 * \param theFirstNode1 - a boundary node of set 1
7517 * \param theFirstNode2 - a node of set 2 corresponding to theFirstNode1
7518 * \param theSecondNode1 - a boundary node of set 1 linked with theFirstNode1
7519 * \param theSecondNode2 - a node of set 2 corresponding to theSecondNode1
7520 * \param nReplaceMap - output map of corresponding nodes
7521 * \retval bool - is a success or not
7523 //================================================================================
7526 //#define DEBUG_MATCHING_NODES
7529 SMESH_MeshEditor::Sew_Error
7530 SMESH_MeshEditor::FindMatchingNodes(set<const SMDS_MeshElement*>& theSide1,
7531 set<const SMDS_MeshElement*>& theSide2,
7532 const SMDS_MeshNode* theFirstNode1,
7533 const SMDS_MeshNode* theFirstNode2,
7534 const SMDS_MeshNode* theSecondNode1,
7535 const SMDS_MeshNode* theSecondNode2,
7536 TNodeNodeMap & nReplaceMap)
7538 set<const SMDS_MeshElement*> * faceSetPtr[] = { &theSide1, &theSide2 };
7540 nReplaceMap.clear();
7541 if ( theFirstNode1 != theFirstNode2 )
7542 nReplaceMap.insert( make_pair( theFirstNode1, theFirstNode2 ));
7543 if ( theSecondNode1 != theSecondNode2 )
7544 nReplaceMap.insert( make_pair( theSecondNode1, theSecondNode2 ));
7546 set< SMESH_TLink > linkSet; // set of nodes where order of nodes is ignored
7547 linkSet.insert( SMESH_TLink( theFirstNode1, theSecondNode1 ));
7549 list< NLink > linkList[2];
7550 linkList[0].push_back( NLink( theFirstNode1, theSecondNode1 ));
7551 linkList[1].push_back( NLink( theFirstNode2, theSecondNode2 ));
7553 // loop on links in linkList; find faces by links and append links
7554 // of the found faces to linkList
7555 list< NLink >::iterator linkIt[] = { linkList[0].begin(), linkList[1].begin() } ;
7556 for ( ; linkIt[0] != linkList[0].end(); linkIt[0]++, linkIt[1]++ ) {
7557 NLink link[] = { *linkIt[0], *linkIt[1] };
7558 if ( linkSet.find( link[0] ) == linkSet.end() )
7561 // by links, find faces in the face sets,
7562 // and find indices of link nodes in the found faces;
7563 // in a face set, there is only one or no face sharing a link
7564 // ---------------------------------------------------------------
7566 const SMDS_MeshElement* face[] = { 0, 0 };
7567 list<const SMDS_MeshNode*> notLinkNodes[2];
7568 //bool reverse[] = { false, false }; // order of notLinkNodes
7570 for ( int iSide = 0; iSide < 2; iSide++ ) // loop on 2 sides
7572 const SMDS_MeshNode* n1 = link[iSide].first;
7573 const SMDS_MeshNode* n2 = link[iSide].second;
7574 set<const SMDS_MeshElement*> * faceSet = faceSetPtr[ iSide ];
7575 set< const SMDS_MeshElement* > facesOfNode1;
7576 for ( int iNode = 0; iNode < 2; iNode++ ) // loop on 2 nodes of a link
7578 // during a loop of the first node, we find all faces around n1,
7579 // during a loop of the second node, we find one face sharing both n1 and n2
7580 const SMDS_MeshNode* n = iNode ? n1 : n2; // a node of a link
7581 SMDS_ElemIteratorPtr fIt = n->GetInverseElementIterator(SMDSAbs_Face);
7582 while ( fIt->more() ) { // loop on faces sharing a node
7583 const SMDS_MeshElement* f = fIt->next();
7584 if (faceSet->find( f ) != faceSet->end() && // f is in face set
7585 ! facesOfNode1.insert( f ).second ) // f encounters twice
7587 if ( face[ iSide ] ) {
7588 MESSAGE( "2 faces per link " );
7589 return ( iSide ? SEW_BAD_SIDE2_NODES : SEW_BAD_SIDE1_NODES );
7592 faceSet->erase( f );
7594 // get not link nodes
7595 int nbN = f->NbNodes();
7596 if ( f->IsQuadratic() )
7598 nbNodes[ iSide ] = nbN;
7599 list< const SMDS_MeshNode* > & nodes = notLinkNodes[ iSide ];
7600 int i1 = f->GetNodeIndex( n1 );
7601 int i2 = f->GetNodeIndex( n2 );
7602 int iEnd = nbN, iBeg = -1, iDelta = 1;
7603 bool reverse = ( Abs( i1 - i2 ) == 1 ? i1 > i2 : i2 > i1 );
7605 std::swap( iEnd, iBeg ); iDelta = -1;
7610 if ( i == iEnd ) i = iBeg + iDelta;
7611 if ( i == i1 ) break;
7612 nodes.push_back ( f->GetNode( i ) );
7618 // check similarity of elements of the sides
7619 if (( face[0] && !face[1] ) || ( !face[0] && face[1] )) {
7620 MESSAGE("Correspondent face not found on side " << ( face[0] ? 1 : 0 ));
7621 if ( nReplaceMap.size() == 2 ) { // faces on input nodes not found
7622 return ( face[0] ? SEW_BAD_SIDE2_NODES : SEW_BAD_SIDE1_NODES );
7625 return SEW_TOPO_DIFF_SETS_OF_ELEMENTS;
7629 // set nodes to merge
7630 // -------------------
7632 if ( face[0] && face[1] ) {
7633 if ( nbNodes[0] != nbNodes[1] ) {
7634 MESSAGE("Diff nb of face nodes");
7635 return SEW_TOPO_DIFF_SETS_OF_ELEMENTS;
7637 #ifdef DEBUG_MATCHING_NODES
7638 MESSAGE ( " Link 1: " << link[0].first->GetID() <<" "<< link[0].second->GetID()
7639 << " F 1: " << face[0] << "| Link 2: " << link[1].first->GetID() <<" "
7640 << link[1].second->GetID() << " F 2: " << face[1] << " | Bind: " ) ;
7642 int nbN = nbNodes[0];
7644 list<const SMDS_MeshNode*>::iterator n1 = notLinkNodes[0].begin();
7645 list<const SMDS_MeshNode*>::iterator n2 = notLinkNodes[1].begin();
7646 for ( int i = 0 ; i < nbN - 2; ++i ) {
7647 #ifdef DEBUG_MATCHING_NODES
7648 MESSAGE ( (*n1)->GetID() << " to " << (*n2)->GetID() );
7650 nReplaceMap.insert( make_pair( *(n1++), *(n2++) ));
7654 // add other links of the face 1 to linkList
7655 // -----------------------------------------
7657 const SMDS_MeshElement* f0 = face[0];
7658 const SMDS_MeshNode* n1 = f0->GetNode( nbN - 1 );
7659 for ( int i = 0; i < nbN; i++ )
7661 const SMDS_MeshNode* n2 = f0->GetNode( i );
7662 pair< set< SMESH_TLink >::iterator, bool > iter_isnew =
7663 linkSet.insert( SMESH_TLink( n1, n2 ));
7664 if ( !iter_isnew.second ) { // already in a set: no need to process
7665 linkSet.erase( iter_isnew.first );
7667 else // new in set == encountered for the first time: add
7669 #ifdef DEBUG_MATCHING_NODES
7670 MESSAGE ( "Add link 1: " << n1->GetID() << " " << n2->GetID() << " "
7671 << " | link 2: " << nReplaceMap[n1]->GetID() << " " << nReplaceMap[n2]->GetID() << " " );
7673 linkList[0].push_back ( NLink( n1, n2 ));
7674 linkList[1].push_back ( NLink( nReplaceMap[n1], nReplaceMap[n2] ));
7679 } // loop on link lists
7685 \brief Creates a hole in a mesh by doubling the nodes of some particular elements
7686 \param theNodes - identifiers of nodes to be doubled
7687 \param theModifiedElems - identifiers of elements to be updated by the new (doubled)
7688 nodes. If list of element identifiers is empty then nodes are doubled but
7689 they not assigned to elements
7690 \return TRUE if operation has been completed successfully, FALSE otherwise
7692 bool SMESH_MeshEditor::DoubleNodes( const std::list< int >& theListOfNodes,
7693 const std::list< int >& theListOfModifiedElems )
7695 myLastCreatedElems.Clear();
7696 myLastCreatedNodes.Clear();
7698 if ( theListOfNodes.size() == 0 )
7701 SMESHDS_Mesh* aMeshDS = GetMeshDS();
7705 // iterate through nodes and duplicate them
7707 std::map< const SMDS_MeshNode*, const SMDS_MeshNode* > anOldNodeToNewNode;
7709 std::list< int >::const_iterator aNodeIter;
7710 for ( aNodeIter = theListOfNodes.begin(); aNodeIter != theListOfNodes.end(); ++aNodeIter )
7712 int aCurr = *aNodeIter;
7713 SMDS_MeshNode* aNode = (SMDS_MeshNode*)aMeshDS->FindNode( aCurr );
7719 const SMDS_MeshNode* aNewNode = aMeshDS->AddNode( aNode->X(), aNode->Y(), aNode->Z() );
7722 anOldNodeToNewNode[ aNode ] = aNewNode;
7723 myLastCreatedNodes.Append( aNewNode );
7727 // Create map of new nodes for modified elements
7729 std::map< SMDS_MeshElement*, vector<const SMDS_MeshNode*> > anElemToNodes;
7731 std::list< int >::const_iterator anElemIter;
7732 for ( anElemIter = theListOfModifiedElems.begin();
7733 anElemIter != theListOfModifiedElems.end(); ++anElemIter )
7735 int aCurr = *anElemIter;
7736 SMDS_MeshElement* anElem = (SMDS_MeshElement*)aMeshDS->FindElement( aCurr );
7740 vector<const SMDS_MeshNode*> aNodeArr( anElem->NbNodes() );
7742 SMDS_ElemIteratorPtr anIter = anElem->nodesIterator();
7744 while ( anIter->more() )
7746 SMDS_MeshNode* aCurrNode = (SMDS_MeshNode*)anIter->next();
7747 if ( aCurr && anOldNodeToNewNode.find( aCurrNode ) != anOldNodeToNewNode.end() )
7749 const SMDS_MeshNode* aNewNode = anOldNodeToNewNode[ aCurrNode ];
7750 aNodeArr[ ind++ ] = aNewNode;
7753 aNodeArr[ ind++ ] = aCurrNode;
7755 anElemToNodes[ anElem ] = aNodeArr;
7758 // Change nodes of elements
7760 std::map< SMDS_MeshElement*, vector<const SMDS_MeshNode*> >::iterator
7761 anElemToNodesIter = anElemToNodes.begin();
7762 for ( ; anElemToNodesIter != anElemToNodes.end(); ++anElemToNodesIter )
7764 const SMDS_MeshElement* anElem = anElemToNodesIter->first;
7765 vector<const SMDS_MeshNode*> aNodeArr = anElemToNodesIter->second;
7767 aMeshDS->ChangeElementNodes( anElem, &aNodeArr[ 0 ], anElem->NbNodes() );