1 // Copyright (C) 2007-2016 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, or (at your option) any later version.
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
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_FacePosition.hxx"
33 #include "SMDS_SpacePosition.hxx"
34 #include "SMDS_MeshGroup.hxx"
35 #include "SMDS_LinearEdge.hxx"
36 #include "SMDS_Downward.hxx"
37 #include "SMDS_SetIterator.hxx"
39 #include "SMESHDS_Group.hxx"
40 #include "SMESHDS_Mesh.hxx"
42 #include "SMESH_Algo.hxx"
43 #include "SMESH_ControlsDef.hxx"
44 #include "SMESH_Group.hxx"
45 #include "SMESH_MeshAlgos.hxx"
46 #include "SMESH_MesherHelper.hxx"
47 #include "SMESH_OctreeNode.hxx"
48 #include "SMESH_subMesh.hxx"
50 #include <Basics_OCCTVersion.hxx>
52 #include "utilities.h"
55 #include <BRepAdaptor_Surface.hxx>
56 #include <BRepBuilderAPI_MakeEdge.hxx>
57 #include <BRepClass3d_SolidClassifier.hxx>
58 #include <BRep_Tool.hxx>
60 #include <Extrema_GenExtPS.hxx>
61 #include <Extrema_POnCurv.hxx>
62 #include <Extrema_POnSurf.hxx>
63 #include <Geom2d_Curve.hxx>
64 #include <GeomAdaptor_Surface.hxx>
65 #include <Geom_Curve.hxx>
66 #include <Geom_Surface.hxx>
67 #include <Precision.hxx>
68 #include <TColStd_ListOfInteger.hxx>
69 #include <TopAbs_State.hxx>
71 #include <TopExp_Explorer.hxx>
72 #include <TopTools_ListIteratorOfListOfShape.hxx>
73 #include <TopTools_ListOfShape.hxx>
74 #include <TopTools_SequenceOfShape.hxx>
76 #include <TopoDS_Face.hxx>
77 #include <TopoDS_Solid.hxx>
83 #include <gp_Trsf.hxx>
97 #include <boost/tuple/tuple.hpp>
99 #include <Standard_Failure.hxx>
100 #include <Standard_ErrorHandler.hxx>
102 #define cast2Node(elem) static_cast<const SMDS_MeshNode*>( elem )
105 using namespace SMESH::Controls;
109 template < class ELEM_SET >
110 SMDS_ElemIteratorPtr elemSetIterator( const ELEM_SET& elements )
112 typedef SMDS_SetIterator
113 < SMDS_pElement, typename ELEM_SET::const_iterator> TSetIterator;
114 return SMDS_ElemIteratorPtr( new TSetIterator( elements.begin(), elements.end() ));
118 //=======================================================================
119 //function : SMESH_MeshEditor
121 //=======================================================================
123 SMESH_MeshEditor::SMESH_MeshEditor( SMESH_Mesh* theMesh )
124 :myMesh( theMesh ) // theMesh may be NULL
128 //================================================================================
130 * \brief Clears myLastCreatedNodes and myLastCreatedElems
132 //================================================================================
134 void SMESH_MeshEditor::ClearLastCreated()
136 myLastCreatedNodes.Clear();
137 myLastCreatedElems.Clear();
140 //================================================================================
142 * \brief Initializes members by an existing element
143 * \param [in] elem - the source element
144 * \param [in] basicOnly - if true, does not set additional data of Ball and Polyhedron
146 //================================================================================
148 SMESH_MeshEditor::ElemFeatures&
149 SMESH_MeshEditor::ElemFeatures::Init( const SMDS_MeshElement* elem, bool basicOnly )
153 myType = elem->GetType();
154 if ( myType == SMDSAbs_Face || myType == SMDSAbs_Volume )
156 myIsPoly = elem->IsPoly();
159 myIsQuad = elem->IsQuadratic();
160 if ( myType == SMDSAbs_Volume && !basicOnly )
162 vector<int > quant = static_cast<const SMDS_VtkVolume* >( elem )->GetQuantities();
163 myPolyhedQuantities.swap( quant );
167 else if ( myType == SMDSAbs_Ball && !basicOnly )
169 myBallDiameter = static_cast<const SMDS_BallElement*>(elem)->GetDiameter();
175 //=======================================================================
179 //=======================================================================
182 SMESH_MeshEditor::AddElement(const vector<const SMDS_MeshNode*> & node,
183 const ElemFeatures& features)
185 SMDS_MeshElement* e = 0;
186 int nbnode = node.size();
187 SMESHDS_Mesh* mesh = GetMeshDS();
188 const int ID = features.myID;
190 switch ( features.myType ) {
192 if ( !features.myIsPoly ) {
194 if ( ID >= 1 ) e = mesh->AddFaceWithID(node[0], node[1], node[2], ID);
195 else e = mesh->AddFace (node[0], node[1], node[2] );
197 else if (nbnode == 4) {
198 if ( ID >= 1 ) e = mesh->AddFaceWithID(node[0], node[1], node[2], node[3], ID);
199 else e = mesh->AddFace (node[0], node[1], node[2], node[3] );
201 else if (nbnode == 6) {
202 if ( ID >= 1 ) e = mesh->AddFaceWithID(node[0], node[1], node[2], node[3],
203 node[4], node[5], ID);
204 else e = mesh->AddFace (node[0], node[1], node[2], node[3],
207 else if (nbnode == 7) {
208 if ( ID >= 1 ) e = mesh->AddFaceWithID(node[0], node[1], node[2], node[3],
209 node[4], node[5], node[6], ID);
210 else e = mesh->AddFace (node[0], node[1], node[2], node[3],
211 node[4], node[5], node[6] );
213 else if (nbnode == 8) {
214 if ( ID >= 1 ) e = mesh->AddFaceWithID(node[0], node[1], node[2], node[3],
215 node[4], node[5], node[6], node[7], ID);
216 else e = mesh->AddFace (node[0], node[1], node[2], node[3],
217 node[4], node[5], node[6], node[7] );
219 else if (nbnode == 9) {
220 if ( ID >= 1 ) e = mesh->AddFaceWithID(node[0], node[1], node[2], node[3],
221 node[4], node[5], node[6], node[7], node[8], ID);
222 else e = mesh->AddFace (node[0], node[1], node[2], node[3],
223 node[4], node[5], node[6], node[7], node[8] );
226 else if ( !features.myIsQuad )
228 if ( ID >= 1 ) e = mesh->AddPolygonalFaceWithID(node, ID);
229 else e = mesh->AddPolygonalFace (node );
231 else if ( nbnode % 2 == 0 ) // just a protection
233 if ( ID >= 1 ) e = mesh->AddQuadPolygonalFaceWithID(node, ID);
234 else e = mesh->AddQuadPolygonalFace (node );
239 if ( !features.myIsPoly ) {
241 if ( ID >= 1 ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3], ID);
242 else e = mesh->AddVolume (node[0], node[1], node[2], node[3] );
244 else if (nbnode == 5) {
245 if ( ID >= 1 ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3],
247 else e = mesh->AddVolume (node[0], node[1], node[2], node[3],
250 else if (nbnode == 6) {
251 if ( ID >= 1 ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3],
252 node[4], node[5], ID);
253 else e = mesh->AddVolume (node[0], node[1], node[2], node[3],
256 else if (nbnode == 8) {
257 if ( ID >= 1 ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3],
258 node[4], node[5], node[6], node[7], ID);
259 else e = mesh->AddVolume (node[0], node[1], node[2], node[3],
260 node[4], node[5], node[6], node[7] );
262 else if (nbnode == 10) {
263 if ( ID >= 1 ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3],
264 node[4], node[5], node[6], node[7],
265 node[8], node[9], ID);
266 else e = mesh->AddVolume (node[0], node[1], node[2], node[3],
267 node[4], node[5], node[6], node[7],
270 else if (nbnode == 12) {
271 if ( ID >= 1 ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3],
272 node[4], node[5], node[6], node[7],
273 node[8], node[9], node[10], node[11], ID);
274 else e = mesh->AddVolume (node[0], node[1], node[2], node[3],
275 node[4], node[5], node[6], node[7],
276 node[8], node[9], node[10], node[11] );
278 else if (nbnode == 13) {
279 if ( ID >= 1 ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3],
280 node[4], node[5], node[6], node[7],
281 node[8], node[9], node[10],node[11],
283 else e = mesh->AddVolume (node[0], node[1], node[2], node[3],
284 node[4], node[5], node[6], node[7],
285 node[8], node[9], node[10],node[11],
288 else if (nbnode == 15) {
289 if ( ID >= 1 ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3],
290 node[4], node[5], node[6], node[7],
291 node[8], node[9], node[10],node[11],
292 node[12],node[13],node[14],ID);
293 else e = mesh->AddVolume (node[0], node[1], node[2], node[3],
294 node[4], node[5], node[6], node[7],
295 node[8], node[9], node[10],node[11],
296 node[12],node[13],node[14] );
298 else if (nbnode == 20) {
299 if ( ID >= 1 ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3],
300 node[4], node[5], node[6], node[7],
301 node[8], node[9], node[10],node[11],
302 node[12],node[13],node[14],node[15],
303 node[16],node[17],node[18],node[19],ID);
304 else e = mesh->AddVolume (node[0], node[1], node[2], node[3],
305 node[4], node[5], node[6], node[7],
306 node[8], node[9], node[10],node[11],
307 node[12],node[13],node[14],node[15],
308 node[16],node[17],node[18],node[19] );
310 else if (nbnode == 27) {
311 if ( ID >= 1 ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3],
312 node[4], node[5], node[6], node[7],
313 node[8], node[9], node[10],node[11],
314 node[12],node[13],node[14],node[15],
315 node[16],node[17],node[18],node[19],
316 node[20],node[21],node[22],node[23],
317 node[24],node[25],node[26], ID);
318 else e = mesh->AddVolume (node[0], node[1], node[2], node[3],
319 node[4], node[5], node[6], node[7],
320 node[8], node[9], node[10],node[11],
321 node[12],node[13],node[14],node[15],
322 node[16],node[17],node[18],node[19],
323 node[20],node[21],node[22],node[23],
324 node[24],node[25],node[26] );
327 else if ( !features.myIsQuad )
329 if ( ID >= 1 ) e = mesh->AddPolyhedralVolumeWithID(node, features.myPolyhedQuantities, ID);
330 else e = mesh->AddPolyhedralVolume (node, features.myPolyhedQuantities );
334 // if ( ID >= 1 ) e = mesh->AddQuadPolyhedralVolumeWithID(node, features.myPolyhedQuantities,ID);
335 // else e = mesh->AddQuadPolyhedralVolume (node, features.myPolyhedQuantities );
341 if ( ID >= 1 ) e = mesh->AddEdgeWithID(node[0], node[1], ID);
342 else e = mesh->AddEdge (node[0], node[1] );
344 else if ( nbnode == 3 ) {
345 if ( ID >= 1 ) e = mesh->AddEdgeWithID(node[0], node[1], node[2], ID);
346 else e = mesh->AddEdge (node[0], node[1], node[2] );
350 case SMDSAbs_0DElement:
352 if ( ID >= 1 ) e = mesh->Add0DElementWithID(node[0], ID);
353 else e = mesh->Add0DElement (node[0] );
358 if ( ID >= 1 ) e = mesh->AddNodeWithID(node[0]->X(), node[0]->Y(), node[0]->Z(), ID);
359 else e = mesh->AddNode (node[0]->X(), node[0]->Y(), node[0]->Z() );
363 if ( ID >= 1 ) e = mesh->AddBallWithID(node[0], features.myBallDiameter, ID);
364 else e = mesh->AddBall (node[0], features.myBallDiameter );
369 if ( e ) myLastCreatedElems.Append( e );
373 //=======================================================================
377 //=======================================================================
379 SMDS_MeshElement* SMESH_MeshEditor::AddElement(const vector<int> & nodeIDs,
380 const ElemFeatures& features)
382 vector<const SMDS_MeshNode*> nodes;
383 nodes.reserve( nodeIDs.size() );
384 vector<int>::const_iterator id = nodeIDs.begin();
385 while ( id != nodeIDs.end() ) {
386 if ( const SMDS_MeshNode* node = GetMeshDS()->FindNode( *id++ ))
387 nodes.push_back( node );
391 return AddElement( nodes, features );
394 //=======================================================================
396 //purpose : Remove a node or an element.
397 // Modify a compute state of sub-meshes which become empty
398 //=======================================================================
400 int SMESH_MeshEditor::Remove (const list< int >& theIDs,
403 myLastCreatedElems.Clear();
404 myLastCreatedNodes.Clear();
406 SMESHDS_Mesh* aMesh = GetMeshDS();
407 set< SMESH_subMesh *> smmap;
410 list<int>::const_iterator it = theIDs.begin();
411 for ( ; it != theIDs.end(); it++ ) {
412 const SMDS_MeshElement * elem;
414 elem = aMesh->FindNode( *it );
416 elem = aMesh->FindElement( *it );
420 // Notify VERTEX sub-meshes about modification
422 const SMDS_MeshNode* node = cast2Node( elem );
423 if ( node->GetPosition()->GetTypeOfPosition() == SMDS_TOP_VERTEX )
424 if ( int aShapeID = node->getshapeId() )
425 if ( SMESH_subMesh * sm = GetMesh()->GetSubMeshContaining( aShapeID ) )
428 // Find sub-meshes to notify about modification
429 // SMDS_ElemIteratorPtr nodeIt = elem->nodesIterator();
430 // while ( nodeIt->more() ) {
431 // const SMDS_MeshNode* node = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
432 // const SMDS_PositionPtr& aPosition = node->GetPosition();
433 // if ( aPosition.get() ) {
434 // if ( int aShapeID = aPosition->GetShapeId() ) {
435 // if ( SMESH_subMesh * sm = GetMesh()->GetSubMeshContaining( aShapeID ) )
436 // smmap.insert( sm );
443 aMesh->RemoveNode( static_cast< const SMDS_MeshNode* >( elem ));
445 aMesh->RemoveElement( elem );
449 // Notify sub-meshes about modification
450 if ( !smmap.empty() ) {
451 set< SMESH_subMesh *>::iterator smIt;
452 for ( smIt = smmap.begin(); smIt != smmap.end(); smIt++ )
453 (*smIt)->ComputeStateEngine( SMESH_subMesh::MESH_ENTITY_REMOVED );
456 // // Check if the whole mesh becomes empty
457 // if ( SMESH_subMesh * sm = GetMesh()->GetSubMeshContaining( 1 ) )
458 // sm->ComputeStateEngine( SMESH_subMesh::CHECK_COMPUTE_STATE );
463 //================================================================================
465 * \brief Create 0D elements on all nodes of the given object except those
466 * nodes on which a 0D element already exists.
467 * \param elements - Elements on whose nodes to create 0D elements; if empty,
468 * the all mesh is treated
469 * \param all0DElems - returns all 0D elements found or created on nodes of \a elements
471 //================================================================================
473 void SMESH_MeshEditor::Create0DElementsOnAllNodes( const TIDSortedElemSet& elements,
474 TIDSortedElemSet& all0DElems )
476 SMDS_ElemIteratorPtr elemIt;
477 vector< const SMDS_MeshElement* > allNodes;
478 if ( elements.empty() )
480 allNodes.reserve( GetMeshDS()->NbNodes() );
481 elemIt = GetMeshDS()->elementsIterator( SMDSAbs_Node );
482 while ( elemIt->more() )
483 allNodes.push_back( elemIt->next() );
485 elemIt = elemSetIterator( allNodes );
489 elemIt = elemSetIterator( elements );
492 while ( elemIt->more() )
494 const SMDS_MeshElement* e = elemIt->next();
495 SMDS_ElemIteratorPtr nodeIt = e->nodesIterator();
496 while ( nodeIt->more() )
498 const SMDS_MeshNode* n = cast2Node( nodeIt->next() );
499 SMDS_ElemIteratorPtr it0D = n->GetInverseElementIterator( SMDSAbs_0DElement );
501 all0DElems.insert( it0D->next() );
503 myLastCreatedElems.Append( GetMeshDS()->Add0DElement( n ));
504 all0DElems.insert( myLastCreatedElems.Last() );
510 //=======================================================================
511 //function : FindShape
512 //purpose : Return an index of the shape theElem is on
513 // or zero if a shape not found
514 //=======================================================================
516 int SMESH_MeshEditor::FindShape (const SMDS_MeshElement * theElem)
518 myLastCreatedElems.Clear();
519 myLastCreatedNodes.Clear();
521 SMESHDS_Mesh * aMesh = GetMeshDS();
522 if ( aMesh->ShapeToMesh().IsNull() )
525 int aShapeID = theElem->getshapeId();
529 if ( SMESHDS_SubMesh * sm = aMesh->MeshElements( aShapeID ))
530 if ( sm->Contains( theElem ))
533 if ( theElem->GetType() == SMDSAbs_Node ) {
534 MESSAGE( ":( Error: invalid myShapeId of node " << theElem->GetID() );
537 MESSAGE( ":( Error: invalid myShapeId of element " << theElem->GetID() );
540 TopoDS_Shape aShape; // the shape a node of theElem is on
541 if ( theElem->GetType() != SMDSAbs_Node )
543 SMDS_ElemIteratorPtr nodeIt = theElem->nodesIterator();
544 while ( nodeIt->more() ) {
545 const SMDS_MeshNode* node = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
546 if ((aShapeID = node->getshapeId()) > 0) {
547 if ( SMESHDS_SubMesh * sm = aMesh->MeshElements( aShapeID ) ) {
548 if ( sm->Contains( theElem ))
550 if ( aShape.IsNull() )
551 aShape = aMesh->IndexToShape( aShapeID );
557 // None of nodes is on a proper shape,
558 // find the shape among ancestors of aShape on which a node is
559 if ( !aShape.IsNull() ) {
560 TopTools_ListIteratorOfListOfShape ancIt( GetMesh()->GetAncestors( aShape ));
561 for ( ; ancIt.More(); ancIt.Next() ) {
562 SMESHDS_SubMesh * sm = aMesh->MeshElements( ancIt.Value() );
563 if ( sm && sm->Contains( theElem ))
564 return aMesh->ShapeToIndex( ancIt.Value() );
569 SMESHDS_SubMeshIteratorPtr smIt = GetMeshDS()->SubMeshes();
570 while ( const SMESHDS_SubMesh* sm = smIt->next() )
571 if ( sm->Contains( theElem ))
578 //=======================================================================
579 //function : IsMedium
581 //=======================================================================
583 bool SMESH_MeshEditor::IsMedium(const SMDS_MeshNode* node,
584 const SMDSAbs_ElementType typeToCheck)
586 bool isMedium = false;
587 SMDS_ElemIteratorPtr it = node->GetInverseElementIterator(typeToCheck);
588 while (it->more() && !isMedium ) {
589 const SMDS_MeshElement* elem = it->next();
590 isMedium = elem->IsMediumNode(node);
595 //=======================================================================
596 //function : shiftNodesQuadTria
597 //purpose : Shift nodes in the array corresponded to quadratic triangle
598 // example: (0,1,2,3,4,5) -> (1,2,0,4,5,3)
599 //=======================================================================
601 static void shiftNodesQuadTria(vector< const SMDS_MeshNode* >& aNodes)
603 const SMDS_MeshNode* nd1 = aNodes[0];
604 aNodes[0] = aNodes[1];
605 aNodes[1] = aNodes[2];
607 const SMDS_MeshNode* nd2 = aNodes[3];
608 aNodes[3] = aNodes[4];
609 aNodes[4] = aNodes[5];
613 //=======================================================================
614 //function : nbEdgeConnectivity
615 //purpose : return number of the edges connected with the theNode.
616 // if theEdges has connections with the other type of the
617 // elements, return -1
618 //=======================================================================
620 static int nbEdgeConnectivity(const SMDS_MeshNode* theNode)
622 // SMDS_ElemIteratorPtr elemIt = theNode->GetInverseElementIterator();
624 // while(elemIt->more()) {
629 return theNode->NbInverseElements();
632 //=======================================================================
633 //function : getNodesFromTwoTria
635 //=======================================================================
637 static bool getNodesFromTwoTria(const SMDS_MeshElement * theTria1,
638 const SMDS_MeshElement * theTria2,
639 vector< const SMDS_MeshNode*>& N1,
640 vector< const SMDS_MeshNode*>& N2)
642 N1.assign( theTria1->begin_nodes(), theTria1->end_nodes() );
643 if ( N1.size() < 6 ) return false;
644 N2.assign( theTria2->begin_nodes(), theTria2->end_nodes() );
645 if ( N2.size() < 6 ) return false;
647 int sames[3] = {-1,-1,-1};
659 if(nbsames!=2) return false;
661 shiftNodesQuadTria(N1);
663 shiftNodesQuadTria(N1);
666 i = sames[0] + sames[1] + sames[2];
668 shiftNodesQuadTria(N2);
670 // now we receive following N1 and N2 (using numeration as in the image below)
671 // tria1 : (1 2 4 5 9 7) and tria2 : (3 4 2 8 9 6)
672 // i.e. first nodes from both arrays form a new diagonal
676 //=======================================================================
677 //function : InverseDiag
678 //purpose : Replace two neighbour triangles with ones built on the same 4 nodes
679 // but having other common link.
680 // Return False if args are improper
681 //=======================================================================
683 bool SMESH_MeshEditor::InverseDiag (const SMDS_MeshElement * theTria1,
684 const SMDS_MeshElement * theTria2 )
686 myLastCreatedElems.Clear();
687 myLastCreatedNodes.Clear();
689 if (!theTria1 || !theTria2)
692 const SMDS_VtkFace* F1 = dynamic_cast<const SMDS_VtkFace*>( theTria1 );
693 if (!F1) return false;
694 const SMDS_VtkFace* F2 = dynamic_cast<const SMDS_VtkFace*>( theTria2 );
695 if (!F2) return false;
696 if ((theTria1->GetEntityType() == SMDSEntity_Triangle) &&
697 (theTria2->GetEntityType() == SMDSEntity_Triangle)) {
699 // 1 +--+ A theTria1: ( 1 A B ) A->2 ( 1 2 B ) 1 +--+ A
700 // | /| theTria2: ( B A 2 ) B->1 ( 1 A 2 ) |\ |
704 // put nodes in array and find out indices of the same ones
705 const SMDS_MeshNode* aNodes [6];
706 int sameInd [] = { -1, -1, -1, -1, -1, -1 };
708 SMDS_ElemIteratorPtr it = theTria1->nodesIterator();
709 while ( it->more() ) {
710 aNodes[ i ] = static_cast<const SMDS_MeshNode*>( it->next() );
712 if ( i > 2 ) // theTria2
713 // find same node of theTria1
714 for ( int j = 0; j < 3; j++ )
715 if ( aNodes[ i ] == aNodes[ j ]) {
724 return false; // theTria1 is not a triangle
725 it = theTria2->nodesIterator();
727 if ( i == 6 && it->more() )
728 return false; // theTria2 is not a triangle
731 // find indices of 1,2 and of A,B in theTria1
732 int iA = -1, iB = 0, i1 = 0, i2 = 0;
733 for ( i = 0; i < 6; i++ ) {
734 if ( sameInd [ i ] == -1 ) {
739 if ( iA >= 0) iB = i;
743 // nodes 1 and 2 should not be the same
744 if ( aNodes[ i1 ] == aNodes[ i2 ] )
748 aNodes[ iA ] = aNodes[ i2 ];
750 aNodes[ sameInd[ iB ]] = aNodes[ i1 ];
752 GetMeshDS()->ChangeElementNodes( theTria1, aNodes, 3 );
753 GetMeshDS()->ChangeElementNodes( theTria2, &aNodes[ 3 ], 3 );
757 } // end if(F1 && F2)
759 // check case of quadratic faces
760 if (theTria1->GetEntityType() != SMDSEntity_Quad_Triangle &&
761 theTria1->GetEntityType() != SMDSEntity_BiQuad_Triangle)
763 if (theTria2->GetEntityType() != SMDSEntity_Quad_Triangle&&
764 theTria2->GetEntityType() != SMDSEntity_BiQuad_Triangle)
768 // 1 +--+--+ 2 theTria1: (1 2 4 5 9 7) or (2 4 1 9 7 5) or (4 1 2 7 5 9)
769 // | /| theTria2: (2 3 4 6 8 9) or (3 4 2 8 9 6) or (4 2 3 9 6 8)
777 vector< const SMDS_MeshNode* > N1;
778 vector< const SMDS_MeshNode* > N2;
779 if(!getNodesFromTwoTria(theTria1,theTria2,N1,N2))
781 // now we receive following N1 and N2 (using numeration as above image)
782 // tria1 : (1 2 4 5 9 7) and tria2 : (3 4 2 8 9 6)
783 // i.e. first nodes from both arrays determ new diagonal
785 vector< const SMDS_MeshNode*> N1new( N1.size() );
786 vector< const SMDS_MeshNode*> N2new( N2.size() );
787 N1new.back() = N1.back(); // central node of biquadratic
788 N2new.back() = N2.back();
789 N1new[0] = N1[0]; N2new[0] = N1[0];
790 N1new[1] = N2[0]; N2new[1] = N1[1];
791 N1new[2] = N2[1]; N2new[2] = N2[0];
792 N1new[3] = N1[4]; N2new[3] = N1[3];
793 N1new[4] = N2[3]; N2new[4] = N2[5];
794 N1new[5] = N1[5]; N2new[5] = N1[4];
795 // change nodes in faces
796 GetMeshDS()->ChangeElementNodes( theTria1, &N1new[0], N1new.size() );
797 GetMeshDS()->ChangeElementNodes( theTria2, &N2new[0], N2new.size() );
799 // move the central node of biquadratic triangle
800 SMESH_MesherHelper helper( *GetMesh() );
801 for ( int is2nd = 0; is2nd < 2; ++is2nd )
803 const SMDS_MeshElement* tria = is2nd ? theTria2 : theTria1;
804 vector< const SMDS_MeshNode*>& nodes = is2nd ? N2new : N1new;
805 if ( nodes.size() < 7 )
807 helper.SetSubShape( tria->getshapeId() );
808 const TopoDS_Face& F = TopoDS::Face( helper.GetSubShape() );
812 xyz = ( SMESH_TNodeXYZ( nodes[3] ) +
813 SMESH_TNodeXYZ( nodes[4] ) +
814 SMESH_TNodeXYZ( nodes[5] )) / 3.;
819 gp_XY uv = ( helper.GetNodeUV( F, nodes[3], nodes[2], &checkUV ) +
820 helper.GetNodeUV( F, nodes[4], nodes[0], &checkUV ) +
821 helper.GetNodeUV( F, nodes[5], nodes[1], &checkUV )) / 3.;
823 Handle(Geom_Surface) S = BRep_Tool::Surface(F,loc);
824 xyz = S->Value( uv.X(), uv.Y() );
825 xyz.Transform( loc );
826 if ( nodes[6]->GetPosition()->GetTypeOfPosition() == SMDS_TOP_FACE && // set UV
827 nodes[6]->getshapeId() > 0 )
828 GetMeshDS()->SetNodeOnFace( nodes[6], nodes[6]->getshapeId(), uv.X(), uv.Y() );
830 GetMeshDS()->MoveNode( nodes[6], xyz.X(), xyz.Y(), xyz.Z() );
835 //=======================================================================
836 //function : findTriangles
837 //purpose : find triangles sharing theNode1-theNode2 link
838 //=======================================================================
840 static bool findTriangles(const SMDS_MeshNode * theNode1,
841 const SMDS_MeshNode * theNode2,
842 const SMDS_MeshElement*& theTria1,
843 const SMDS_MeshElement*& theTria2)
845 if ( !theNode1 || !theNode2 ) return false;
847 theTria1 = theTria2 = 0;
849 set< const SMDS_MeshElement* > emap;
850 SMDS_ElemIteratorPtr it = theNode1->GetInverseElementIterator(SMDSAbs_Face);
852 const SMDS_MeshElement* elem = it->next();
853 if ( elem->NbCornerNodes() == 3 )
856 it = theNode2->GetInverseElementIterator(SMDSAbs_Face);
858 const SMDS_MeshElement* elem = it->next();
859 if ( emap.count( elem )) {
867 // theTria1 must be element with minimum ID
868 if ( theTria2->GetID() < theTria1->GetID() )
869 std::swap( theTria2, theTria1 );
877 //=======================================================================
878 //function : InverseDiag
879 //purpose : Replace two neighbour triangles sharing theNode1-theNode2 link
880 // with ones built on the same 4 nodes but having other common link.
881 // Return false if proper faces not found
882 //=======================================================================
884 bool SMESH_MeshEditor::InverseDiag (const SMDS_MeshNode * theNode1,
885 const SMDS_MeshNode * theNode2)
887 myLastCreatedElems.Clear();
888 myLastCreatedNodes.Clear();
890 const SMDS_MeshElement *tr1, *tr2;
891 if ( !findTriangles( theNode1, theNode2, tr1, tr2 ))
894 const SMDS_VtkFace* F1 = dynamic_cast<const SMDS_VtkFace*>( tr1 );
895 if (!F1) return false;
896 const SMDS_VtkFace* F2 = dynamic_cast<const SMDS_VtkFace*>( tr2 );
897 if (!F2) return false;
898 if ((tr1->GetEntityType() == SMDSEntity_Triangle) &&
899 (tr2->GetEntityType() == SMDSEntity_Triangle)) {
901 // 1 +--+ A tr1: ( 1 A B ) A->2 ( 1 2 B ) 1 +--+ A
902 // | /| tr2: ( B A 2 ) B->1 ( 1 A 2 ) |\ |
906 // put nodes in array
907 // and find indices of 1,2 and of A in tr1 and of B in tr2
908 int i, iA1 = 0, i1 = 0;
909 const SMDS_MeshNode* aNodes1 [3];
910 SMDS_ElemIteratorPtr it;
911 for (i = 0, it = tr1->nodesIterator(); it->more(); i++ ) {
912 aNodes1[ i ] = static_cast<const SMDS_MeshNode*>( it->next() );
913 if ( aNodes1[ i ] == theNode1 )
914 iA1 = i; // node A in tr1
915 else if ( aNodes1[ i ] != theNode2 )
919 const SMDS_MeshNode* aNodes2 [3];
920 for (i = 0, it = tr2->nodesIterator(); it->more(); i++ ) {
921 aNodes2[ i ] = static_cast<const SMDS_MeshNode*>( it->next() );
922 if ( aNodes2[ i ] == theNode2 )
923 iB2 = i; // node B in tr2
924 else if ( aNodes2[ i ] != theNode1 )
928 // nodes 1 and 2 should not be the same
929 if ( aNodes1[ i1 ] == aNodes2[ i2 ] )
933 aNodes1[ iA1 ] = aNodes2[ i2 ];
935 aNodes2[ iB2 ] = aNodes1[ i1 ];
937 GetMeshDS()->ChangeElementNodes( tr1, aNodes1, 3 );
938 GetMeshDS()->ChangeElementNodes( tr2, aNodes2, 3 );
943 // check case of quadratic faces
944 return InverseDiag(tr1,tr2);
947 //=======================================================================
948 //function : getQuadrangleNodes
949 //purpose : fill theQuadNodes - nodes of a quadrangle resulting from
950 // fusion of triangles tr1 and tr2 having shared link on
951 // theNode1 and theNode2
952 //=======================================================================
954 bool getQuadrangleNodes(const SMDS_MeshNode * theQuadNodes [],
955 const SMDS_MeshNode * theNode1,
956 const SMDS_MeshNode * theNode2,
957 const SMDS_MeshElement * tr1,
958 const SMDS_MeshElement * tr2 )
960 if( tr1->NbNodes() != tr2->NbNodes() )
962 // find the 4-th node to insert into tr1
963 const SMDS_MeshNode* n4 = 0;
964 SMDS_ElemIteratorPtr it = tr2->nodesIterator();
966 while ( !n4 && i<3 ) {
967 const SMDS_MeshNode * n = cast2Node( it->next() );
969 bool isDiag = ( n == theNode1 || n == theNode2 );
973 // Make an array of nodes to be in a quadrangle
974 int iNode = 0, iFirstDiag = -1;
975 it = tr1->nodesIterator();
978 const SMDS_MeshNode * n = cast2Node( it->next() );
980 bool isDiag = ( n == theNode1 || n == theNode2 );
982 if ( iFirstDiag < 0 )
984 else if ( iNode - iFirstDiag == 1 )
985 theQuadNodes[ iNode++ ] = n4; // insert the 4-th node between diagonal nodes
987 else if ( n == n4 ) {
988 return false; // tr1 and tr2 should not have all the same nodes
990 theQuadNodes[ iNode++ ] = n;
992 if ( iNode == 3 ) // diagonal nodes have 0 and 2 indices
993 theQuadNodes[ iNode ] = n4;
998 //=======================================================================
999 //function : DeleteDiag
1000 //purpose : Replace two neighbour triangles sharing theNode1-theNode2 link
1001 // with a quadrangle built on the same 4 nodes.
1002 // Return false if proper faces not found
1003 //=======================================================================
1005 bool SMESH_MeshEditor::DeleteDiag (const SMDS_MeshNode * theNode1,
1006 const SMDS_MeshNode * theNode2)
1008 myLastCreatedElems.Clear();
1009 myLastCreatedNodes.Clear();
1011 const SMDS_MeshElement *tr1, *tr2;
1012 if ( !findTriangles( theNode1, theNode2, tr1, tr2 ))
1015 const SMDS_VtkFace* F1 = dynamic_cast<const SMDS_VtkFace*>( tr1 );
1016 if (!F1) return false;
1017 const SMDS_VtkFace* F2 = dynamic_cast<const SMDS_VtkFace*>( tr2 );
1018 if (!F2) return false;
1019 SMESHDS_Mesh * aMesh = GetMeshDS();
1021 if ((tr1->GetEntityType() == SMDSEntity_Triangle) &&
1022 (tr2->GetEntityType() == SMDSEntity_Triangle)) {
1024 const SMDS_MeshNode* aNodes [ 4 ];
1025 if ( ! getQuadrangleNodes( aNodes, theNode1, theNode2, tr1, tr2 ))
1028 const SMDS_MeshElement* newElem = 0;
1029 newElem = aMesh->AddFace( aNodes[0], aNodes[1], aNodes[2], aNodes[3] );
1030 myLastCreatedElems.Append(newElem);
1031 AddToSameGroups( newElem, tr1, aMesh );
1032 int aShapeId = tr1->getshapeId();
1035 aMesh->SetMeshElementOnShape( newElem, aShapeId );
1037 aMesh->RemoveElement( tr1 );
1038 aMesh->RemoveElement( tr2 );
1043 // check case of quadratic faces
1044 if (tr1->GetEntityType() != SMDSEntity_Quad_Triangle)
1046 if (tr2->GetEntityType() != SMDSEntity_Quad_Triangle)
1050 // 1 +--+--+ 2 tr1: (1 2 4 5 9 7) or (2 4 1 9 7 5) or (4 1 2 7 5 9)
1051 // | /| tr2: (2 3 4 6 8 9) or (3 4 2 8 9 6) or (4 2 3 9 6 8)
1059 vector< const SMDS_MeshNode* > N1;
1060 vector< const SMDS_MeshNode* > N2;
1061 if(!getNodesFromTwoTria(tr1,tr2,N1,N2))
1063 // now we receive following N1 and N2 (using numeration as above image)
1064 // tria1 : (1 2 4 5 9 7) and tria2 : (3 4 2 8 9 6)
1065 // i.e. first nodes from both arrays determ new diagonal
1067 const SMDS_MeshNode* aNodes[8];
1077 const SMDS_MeshElement* newElem = 0;
1078 newElem = aMesh->AddFace( aNodes[0], aNodes[1], aNodes[2], aNodes[3],
1079 aNodes[4], aNodes[5], aNodes[6], aNodes[7]);
1080 myLastCreatedElems.Append(newElem);
1081 AddToSameGroups( newElem, tr1, aMesh );
1082 int aShapeId = tr1->getshapeId();
1085 aMesh->SetMeshElementOnShape( newElem, aShapeId );
1087 aMesh->RemoveElement( tr1 );
1088 aMesh->RemoveElement( tr2 );
1090 // remove middle node (9)
1091 GetMeshDS()->RemoveNode( N1[4] );
1096 //=======================================================================
1097 //function : Reorient
1098 //purpose : Reverse theElement orientation
1099 //=======================================================================
1101 bool SMESH_MeshEditor::Reorient (const SMDS_MeshElement * theElem)
1103 myLastCreatedElems.Clear();
1104 myLastCreatedNodes.Clear();
1108 SMDS_ElemIteratorPtr it = theElem->nodesIterator();
1109 if ( !it || !it->more() )
1112 const SMDSAbs_ElementType type = theElem->GetType();
1113 if ( type < SMDSAbs_Edge || type > SMDSAbs_Volume )
1116 const SMDSAbs_EntityType geomType = theElem->GetEntityType();
1117 if ( geomType == SMDSEntity_Polyhedra ) // polyhedron
1119 const SMDS_VtkVolume* aPolyedre =
1120 dynamic_cast<const SMDS_VtkVolume*>( theElem );
1122 MESSAGE("Warning: bad volumic element");
1125 const int nbFaces = aPolyedre->NbFaces();
1126 vector<const SMDS_MeshNode *> poly_nodes;
1127 vector<int> quantities (nbFaces);
1129 // reverse each face of the polyedre
1130 for (int iface = 1; iface <= nbFaces; iface++) {
1131 int inode, nbFaceNodes = aPolyedre->NbFaceNodes(iface);
1132 quantities[iface - 1] = nbFaceNodes;
1134 for (inode = nbFaceNodes; inode >= 1; inode--) {
1135 const SMDS_MeshNode* curNode = aPolyedre->GetFaceNode(iface, inode);
1136 poly_nodes.push_back(curNode);
1139 return GetMeshDS()->ChangePolyhedronNodes( theElem, poly_nodes, quantities );
1141 else // other elements
1143 vector<const SMDS_MeshNode*> nodes( theElem->begin_nodes(), theElem->end_nodes() );
1144 const std::vector<int>& interlace = SMDS_MeshCell::reverseSmdsOrder( geomType, nodes.size() );
1145 if ( interlace.empty() )
1147 std::reverse( nodes.begin(), nodes.end() ); // obsolete, just in case
1151 SMDS_MeshCell::applyInterlace( interlace, nodes );
1153 return GetMeshDS()->ChangeElementNodes( theElem, &nodes[0], nodes.size() );
1158 //================================================================================
1160 * \brief Reorient faces.
1161 * \param theFaces - the faces to reorient. If empty the whole mesh is meant
1162 * \param theDirection - desired direction of normal of \a theFace
1163 * \param theFace - one of \a theFaces that sould be oriented according to
1164 * \a theDirection and whose orientation defines orientation of other faces
1165 * \return number of reoriented faces.
1167 //================================================================================
1169 int SMESH_MeshEditor::Reorient2D (TIDSortedElemSet & theFaces,
1170 const gp_Dir& theDirection,
1171 const SMDS_MeshElement * theFace)
1174 if ( !theFace || theFace->GetType() != SMDSAbs_Face ) return nbReori;
1176 if ( theFaces.empty() )
1178 SMDS_FaceIteratorPtr fIt = GetMeshDS()->facesIterator(/*idInceasingOrder=*/true);
1179 while ( fIt->more() )
1180 theFaces.insert( theFaces.end(), fIt->next() );
1183 // orient theFace according to theDirection
1185 SMESH_MeshAlgos::FaceNormal( theFace, normal, /*normalized=*/false );
1186 if ( normal * theDirection.XYZ() < 0 )
1187 nbReori += Reorient( theFace );
1189 // Orient other faces
1191 set< const SMDS_MeshElement* > startFaces, visitedFaces;
1192 TIDSortedElemSet avoidSet;
1193 set< SMESH_TLink > checkedLinks;
1194 pair< set< SMESH_TLink >::iterator, bool > linkIt_isNew;
1196 if ( theFaces.size() > 1 )// leave 1 face to prevent finding not selected faces
1197 theFaces.erase( theFace );
1198 startFaces.insert( theFace );
1200 int nodeInd1, nodeInd2;
1201 const SMDS_MeshElement* otherFace;
1202 vector< const SMDS_MeshElement* > facesNearLink;
1203 vector< std::pair< int, int > > nodeIndsOfFace;
1205 set< const SMDS_MeshElement* >::iterator startFace = startFaces.begin();
1206 while ( !startFaces.empty() )
1208 startFace = startFaces.begin();
1209 theFace = *startFace;
1210 startFaces.erase( startFace );
1211 if ( !visitedFaces.insert( theFace ).second )
1215 avoidSet.insert(theFace);
1217 NLink link( theFace->GetNode( 0 ), (SMDS_MeshNode *) 0 );
1219 const int nbNodes = theFace->NbCornerNodes();
1220 for ( int i = 0; i < nbNodes; ++i ) // loop on links of theFace
1222 link.second = theFace->GetNode(( i+1 ) % nbNodes );
1223 linkIt_isNew = checkedLinks.insert( link );
1224 if ( !linkIt_isNew.second )
1226 // link has already been checked and won't be encountered more
1227 // if the group (theFaces) is manifold
1228 //checkedLinks.erase( linkIt_isNew.first );
1232 facesNearLink.clear();
1233 nodeIndsOfFace.clear();
1234 while (( otherFace = SMESH_MeshAlgos::FindFaceInSet( link.first, link.second,
1236 &nodeInd1, &nodeInd2 )))
1237 if ( otherFace != theFace)
1239 facesNearLink.push_back( otherFace );
1240 nodeIndsOfFace.push_back( make_pair( nodeInd1, nodeInd2 ));
1241 avoidSet.insert( otherFace );
1243 if ( facesNearLink.size() > 1 )
1245 // NON-MANIFOLD mesh shell !
1246 // select a face most co-directed with theFace,
1247 // other faces won't be visited this time
1249 SMESH_MeshAlgos::FaceNormal( theFace, NF, /*normalized=*/false );
1250 double proj, maxProj = -1;
1251 for ( size_t i = 0; i < facesNearLink.size(); ++i ) {
1252 SMESH_MeshAlgos::FaceNormal( facesNearLink[i], NOF, /*normalized=*/false );
1253 if (( proj = Abs( NF * NOF )) > maxProj ) {
1255 otherFace = facesNearLink[i];
1256 nodeInd1 = nodeIndsOfFace[i].first;
1257 nodeInd2 = nodeIndsOfFace[i].second;
1260 // not to visit rejected faces
1261 for ( size_t i = 0; i < facesNearLink.size(); ++i )
1262 if ( facesNearLink[i] != otherFace && theFaces.size() > 1 )
1263 visitedFaces.insert( facesNearLink[i] );
1265 else if ( facesNearLink.size() == 1 )
1267 otherFace = facesNearLink[0];
1268 nodeInd1 = nodeIndsOfFace.back().first;
1269 nodeInd2 = nodeIndsOfFace.back().second;
1271 if ( otherFace && otherFace != theFace)
1273 // link must be reverse in otherFace if orientation ot otherFace
1274 // is same as that of theFace
1275 if ( abs(nodeInd2-nodeInd1) == 1 ? nodeInd2 > nodeInd1 : nodeInd1 > nodeInd2 )
1277 nbReori += Reorient( otherFace );
1279 startFaces.insert( otherFace );
1282 std::swap( link.first, link.second ); // reverse the link
1288 //================================================================================
1290 * \brief Reorient faces basing on orientation of adjacent volumes.
1291 * \param theFaces - faces to reorient. If empty, all mesh faces are treated.
1292 * \param theVolumes - reference volumes.
1293 * \param theOutsideNormal - to orient faces to have their normal
1294 * pointing either \a outside or \a inside the adjacent volumes.
1295 * \return number of reoriented faces.
1297 //================================================================================
1299 int SMESH_MeshEditor::Reorient2DBy3D (TIDSortedElemSet & theFaces,
1300 TIDSortedElemSet & theVolumes,
1301 const bool theOutsideNormal)
1305 SMDS_ElemIteratorPtr faceIt;
1306 if ( theFaces.empty() )
1307 faceIt = GetMeshDS()->elementsIterator( SMDSAbs_Face );
1309 faceIt = elemSetIterator( theFaces );
1311 vector< const SMDS_MeshNode* > faceNodes;
1312 TIDSortedElemSet checkedVolumes;
1313 set< const SMDS_MeshNode* > faceNodesSet;
1314 SMDS_VolumeTool volumeTool;
1316 while ( faceIt->more() ) // loop on given faces
1318 const SMDS_MeshElement* face = faceIt->next();
1319 if ( face->GetType() != SMDSAbs_Face )
1322 const size_t nbCornersNodes = face->NbCornerNodes();
1323 faceNodes.assign( face->begin_nodes(), face->end_nodes() );
1325 checkedVolumes.clear();
1326 SMDS_ElemIteratorPtr vIt = faceNodes[ 0 ]->GetInverseElementIterator( SMDSAbs_Volume );
1327 while ( vIt->more() )
1329 const SMDS_MeshElement* volume = vIt->next();
1331 if ( !checkedVolumes.insert( volume ).second )
1333 if ( !theVolumes.empty() && !theVolumes.count( volume ))
1336 // is volume adjacent?
1337 bool allNodesCommon = true;
1338 for ( size_t iN = 1; iN < nbCornersNodes && allNodesCommon; ++iN )
1339 allNodesCommon = ( volume->GetNodeIndex( faceNodes[ iN ]) > -1 );
1340 if ( !allNodesCommon )
1343 // get nodes of a corresponding volume facet
1344 faceNodesSet.clear();
1345 faceNodesSet.insert( faceNodes.begin(), faceNodes.end() );
1346 volumeTool.Set( volume );
1347 int facetID = volumeTool.GetFaceIndex( faceNodesSet );
1348 if ( facetID < 0 ) continue;
1349 volumeTool.SetExternalNormal();
1350 const SMDS_MeshNode** facetNodes = volumeTool.GetFaceNodes( facetID );
1352 // compare order of faceNodes and facetNodes
1353 const int iQ = 1 + ( nbCornersNodes < faceNodes.size() );
1355 for ( int i = 0; i < 2; ++i )
1357 const SMDS_MeshNode* n = facetNodes[ i*iQ ];
1358 for ( size_t iN = 0; iN < nbCornersNodes; ++iN )
1359 if ( faceNodes[ iN ] == n )
1365 bool isOutside = Abs( iNN[0]-iNN[1] ) == 1 ? iNN[0] < iNN[1] : iNN[0] > iNN[1];
1366 if ( isOutside != theOutsideNormal )
1367 nbReori += Reorient( face );
1369 } // loop on given faces
1374 //=======================================================================
1375 //function : getBadRate
1377 //=======================================================================
1379 static double getBadRate (const SMDS_MeshElement* theElem,
1380 SMESH::Controls::NumericalFunctorPtr& theCrit)
1382 SMESH::Controls::TSequenceOfXYZ P;
1383 if ( !theElem || !theCrit->GetPoints( theElem, P ))
1385 return theCrit->GetBadRate( theCrit->GetValue( P ), theElem->NbNodes() );
1386 //return theCrit->GetBadRate( theCrit->GetValue( theElem->GetID() ), theElem->NbNodes() );
1389 //=======================================================================
1390 //function : QuadToTri
1391 //purpose : Cut quadrangles into triangles.
1392 // theCrit is used to select a diagonal to cut
1393 //=======================================================================
1395 bool SMESH_MeshEditor::QuadToTri (TIDSortedElemSet & theElems,
1396 SMESH::Controls::NumericalFunctorPtr theCrit)
1398 myLastCreatedElems.Clear();
1399 myLastCreatedNodes.Clear();
1401 if ( !theCrit.get() )
1404 SMESHDS_Mesh * aMesh = GetMeshDS();
1406 Handle(Geom_Surface) surface;
1407 SMESH_MesherHelper helper( *GetMesh() );
1409 TIDSortedElemSet::iterator itElem;
1410 for ( itElem = theElems.begin(); itElem != theElems.end(); itElem++ )
1412 const SMDS_MeshElement* elem = *itElem;
1413 if ( !elem || elem->GetType() != SMDSAbs_Face )
1415 if ( elem->NbCornerNodes() != 4 )
1418 // retrieve element nodes
1419 vector< const SMDS_MeshNode* > aNodes( elem->begin_nodes(), elem->end_nodes() );
1421 // compare two sets of possible triangles
1422 double aBadRate1, aBadRate2; // to what extent a set is bad
1423 SMDS_FaceOfNodes tr1 ( aNodes[0], aNodes[1], aNodes[2] );
1424 SMDS_FaceOfNodes tr2 ( aNodes[2], aNodes[3], aNodes[0] );
1425 aBadRate1 = getBadRate( &tr1, theCrit ) + getBadRate( &tr2, theCrit );
1427 SMDS_FaceOfNodes tr3 ( aNodes[1], aNodes[2], aNodes[3] );
1428 SMDS_FaceOfNodes tr4 ( aNodes[3], aNodes[0], aNodes[1] );
1429 aBadRate2 = getBadRate( &tr3, theCrit ) + getBadRate( &tr4, theCrit );
1431 const int aShapeId = FindShape( elem );
1432 const SMDS_MeshElement* newElem1 = 0;
1433 const SMDS_MeshElement* newElem2 = 0;
1435 if ( !elem->IsQuadratic() ) // split liner quadrangle
1437 // for MaxElementLength2D functor we return minimum diagonal for splitting,
1438 // because aBadRate1=2*len(diagonal 1-3); aBadRate2=2*len(diagonal 2-4)
1439 if ( aBadRate1 <= aBadRate2 ) {
1440 // tr1 + tr2 is better
1441 newElem1 = aMesh->AddFace( aNodes[2], aNodes[3], aNodes[0] );
1442 newElem2 = aMesh->AddFace( aNodes[2], aNodes[0], aNodes[1] );
1445 // tr3 + tr4 is better
1446 newElem1 = aMesh->AddFace( aNodes[3], aNodes[0], aNodes[1] );
1447 newElem2 = aMesh->AddFace( aNodes[3], aNodes[1], aNodes[2] );
1450 else // split quadratic quadrangle
1452 helper.SetIsQuadratic( true );
1453 helper.SetIsBiQuadratic( aNodes.size() == 9 );
1455 helper.AddTLinks( static_cast< const SMDS_MeshFace* >( elem ));
1456 if ( aNodes.size() == 9 )
1458 helper.SetIsBiQuadratic( true );
1459 if ( aBadRate1 <= aBadRate2 )
1460 helper.AddTLinkNode( aNodes[0], aNodes[2], aNodes[8] );
1462 helper.AddTLinkNode( aNodes[1], aNodes[3], aNodes[8] );
1464 // create a new element
1465 if ( aBadRate1 <= aBadRate2 ) {
1466 newElem1 = helper.AddFace( aNodes[2], aNodes[3], aNodes[0] );
1467 newElem2 = helper.AddFace( aNodes[2], aNodes[0], aNodes[1] );
1470 newElem1 = helper.AddFace( aNodes[3], aNodes[0], aNodes[1] );
1471 newElem2 = helper.AddFace( aNodes[3], aNodes[1], aNodes[2] );
1475 // care of a new element
1477 myLastCreatedElems.Append(newElem1);
1478 myLastCreatedElems.Append(newElem2);
1479 AddToSameGroups( newElem1, elem, aMesh );
1480 AddToSameGroups( newElem2, elem, aMesh );
1482 // put a new triangle on the same shape
1484 aMesh->SetMeshElementOnShape( newElem1, aShapeId );
1485 aMesh->SetMeshElementOnShape( newElem2, aShapeId );
1487 aMesh->RemoveElement( elem );
1492 //=======================================================================
1494 * \brief Split each of given quadrangles into 4 triangles.
1495 * \param theElems - The faces to be splitted. If empty all faces are split.
1497 //=======================================================================
1499 void SMESH_MeshEditor::QuadTo4Tri (TIDSortedElemSet & theElems)
1501 myLastCreatedElems.Clear();
1502 myLastCreatedNodes.Clear();
1504 SMESH_MesherHelper helper( *GetMesh() );
1505 helper.SetElementsOnShape( true );
1507 SMDS_ElemIteratorPtr faceIt;
1508 if ( theElems.empty() ) faceIt = GetMeshDS()->elementsIterator(SMDSAbs_Face);
1509 else faceIt = elemSetIterator( theElems );
1512 gp_XY uv [9]; uv[8] = gp_XY(0,0);
1514 vector< const SMDS_MeshNode* > nodes;
1515 SMESHDS_SubMesh* subMeshDS;
1517 Handle(Geom_Surface) surface;
1518 TopLoc_Location loc;
1520 while ( faceIt->more() )
1522 const SMDS_MeshElement* quad = faceIt->next();
1523 if ( !quad || quad->NbCornerNodes() != 4 )
1526 // get a surface the quad is on
1528 if ( quad->getshapeId() < 1 )
1531 helper.SetSubShape( 0 );
1534 else if ( quad->getshapeId() != helper.GetSubShapeID() )
1536 helper.SetSubShape( quad->getshapeId() );
1537 if ( !helper.GetSubShape().IsNull() &&
1538 helper.GetSubShape().ShapeType() == TopAbs_FACE )
1540 F = TopoDS::Face( helper.GetSubShape() );
1541 surface = BRep_Tool::Surface( F, loc );
1542 subMeshDS = GetMeshDS()->MeshElements( quad->getshapeId() );
1546 helper.SetSubShape( 0 );
1551 // create a central node
1553 const SMDS_MeshNode* nCentral;
1554 nodes.assign( quad->begin_nodes(), quad->end_nodes() );
1556 if ( nodes.size() == 9 )
1558 nCentral = nodes.back();
1565 for ( ; iN < nodes.size(); ++iN )
1566 xyz[ iN ] = SMESH_TNodeXYZ( nodes[ iN ] );
1568 for ( ; iN < 8; ++iN ) // mid-side points of a linear qudrangle
1569 xyz[ iN ] = 0.5 * ( xyz[ iN - 4 ] + xyz[( iN - 3 )%4 ] );
1571 xyz[ 8 ] = helper.calcTFI( 0.5, 0.5,
1572 xyz[0], xyz[1], xyz[2], xyz[3],
1573 xyz[4], xyz[5], xyz[6], xyz[7] );
1577 for ( ; iN < nodes.size(); ++iN )
1578 uv[ iN ] = helper.GetNodeUV( F, nodes[iN], nodes[(iN+2)%4], &checkUV );
1580 for ( ; iN < 8; ++iN ) // UV of mid-side points of a linear qudrangle
1581 uv[ iN ] = helper.GetMiddleUV( surface, uv[ iN - 4 ], uv[( iN - 3 )%4 ] );
1583 uv[ 8 ] = helper.calcTFI( 0.5, 0.5,
1584 uv[0], uv[1], uv[2], uv[3],
1585 uv[4], uv[5], uv[6], uv[7] );
1587 gp_Pnt p = surface->Value( uv[8].X(), uv[8].Y() ).Transformed( loc );
1591 nCentral = helper.AddNode( xyz[8].X(), xyz[8].Y(), xyz[8].Z(), /*id=*/0,
1592 uv[8].X(), uv[8].Y() );
1593 myLastCreatedNodes.Append( nCentral );
1596 // create 4 triangles
1598 helper.SetIsQuadratic ( nodes.size() > 4 );
1599 helper.SetIsBiQuadratic( nodes.size() == 9 );
1600 if ( helper.GetIsQuadratic() )
1601 helper.AddTLinks( static_cast< const SMDS_MeshFace*>( quad ));
1603 GetMeshDS()->RemoveFreeElement( quad, subMeshDS, /*fromGroups=*/false );
1605 for ( int i = 0; i < 4; ++i )
1607 SMDS_MeshElement* tria = helper.AddFace( nodes[ i ],
1610 ReplaceElemInGroups( tria, quad, GetMeshDS() );
1611 myLastCreatedElems.Append( tria );
1616 //=======================================================================
1617 //function : BestSplit
1618 //purpose : Find better diagonal for cutting.
1619 //=======================================================================
1621 int SMESH_MeshEditor::BestSplit (const SMDS_MeshElement* theQuad,
1622 SMESH::Controls::NumericalFunctorPtr theCrit)
1624 myLastCreatedElems.Clear();
1625 myLastCreatedNodes.Clear();
1630 if (!theQuad || theQuad->GetType() != SMDSAbs_Face )
1633 if( theQuad->NbNodes()==4 ||
1634 (theQuad->NbNodes()==8 && theQuad->IsQuadratic()) ) {
1636 // retrieve element nodes
1637 const SMDS_MeshNode* aNodes [4];
1638 SMDS_ElemIteratorPtr itN = theQuad->nodesIterator();
1640 //while (itN->more())
1642 aNodes[ i++ ] = static_cast<const SMDS_MeshNode*>( itN->next() );
1644 // compare two sets of possible triangles
1645 double aBadRate1, aBadRate2; // to what extent a set is bad
1646 SMDS_FaceOfNodes tr1 ( aNodes[0], aNodes[1], aNodes[2] );
1647 SMDS_FaceOfNodes tr2 ( aNodes[2], aNodes[3], aNodes[0] );
1648 aBadRate1 = getBadRate( &tr1, theCrit ) + getBadRate( &tr2, theCrit );
1650 SMDS_FaceOfNodes tr3 ( aNodes[1], aNodes[2], aNodes[3] );
1651 SMDS_FaceOfNodes tr4 ( aNodes[3], aNodes[0], aNodes[1] );
1652 aBadRate2 = getBadRate( &tr3, theCrit ) + getBadRate( &tr4, theCrit );
1653 // for MaxElementLength2D functor we return minimum diagonal for splitting,
1654 // because aBadRate1=2*len(diagonal 1-3); aBadRate2=2*len(diagonal 2-4)
1655 if (aBadRate1 <= aBadRate2) // tr1 + tr2 is better
1656 return 1; // diagonal 1-3
1658 return 2; // diagonal 2-4
1665 // Methods of splitting volumes into tetra
1667 const int theHexTo5_1[5*4+1] =
1669 0, 1, 2, 5, 0, 4, 5, 7, 0, 2, 3, 7, 2, 5, 6, 7, 0, 5, 2, 7, -1
1671 const int theHexTo5_2[5*4+1] =
1673 1, 2, 3, 6, 1, 4, 5, 6, 0, 1, 3, 4, 3, 4, 6, 7, 1, 3, 4, 6, -1
1675 const int* theHexTo5[2] = { theHexTo5_1, theHexTo5_2 };
1677 const int theHexTo6_1[6*4+1] =
1679 1, 5, 6, 0, 0, 1, 2, 6, 0, 4, 5, 6, 0, 4, 6, 7, 0, 2, 3, 6, 0, 3, 7, 6, -1
1681 const int theHexTo6_2[6*4+1] =
1683 2, 6, 7, 1, 1, 2, 3, 7, 1, 5, 6, 7, 1, 5, 7, 4, 1, 3, 0, 7, 1, 0, 4, 7, -1
1685 const int theHexTo6_3[6*4+1] =
1687 3, 7, 4, 2, 2, 3, 0, 4, 2, 6, 7, 4, 2, 6, 4, 5, 2, 0, 1, 4, 2, 1, 5, 4, -1
1689 const int theHexTo6_4[6*4+1] =
1691 0, 4, 5, 3, 3, 0, 1, 5, 3, 7, 4, 5, 3, 7, 5, 6, 3, 1, 2, 5, 3, 2, 6, 5, -1
1693 const int* theHexTo6[4] = { theHexTo6_1, theHexTo6_2, theHexTo6_3, theHexTo6_4 };
1695 const int thePyraTo2_1[2*4+1] =
1697 0, 1, 2, 4, 0, 2, 3, 4, -1
1699 const int thePyraTo2_2[2*4+1] =
1701 1, 2, 3, 4, 1, 3, 0, 4, -1
1703 const int* thePyraTo2[2] = { thePyraTo2_1, thePyraTo2_2 };
1705 const int thePentaTo3_1[3*4+1] =
1707 0, 1, 2, 3, 1, 3, 4, 2, 2, 3, 4, 5, -1
1709 const int thePentaTo3_2[3*4+1] =
1711 1, 2, 0, 4, 2, 4, 5, 0, 0, 4, 5, 3, -1
1713 const int thePentaTo3_3[3*4+1] =
1715 2, 0, 1, 5, 0, 5, 3, 1, 1, 5, 3, 4, -1
1717 const int thePentaTo3_4[3*4+1] =
1719 0, 1, 2, 3, 1, 3, 4, 5, 2, 3, 1, 5, -1
1721 const int thePentaTo3_5[3*4+1] =
1723 1, 2, 0, 4, 2, 4, 5, 3, 0, 4, 2, 3, -1
1725 const int thePentaTo3_6[3*4+1] =
1727 2, 0, 1, 5, 0, 5, 3, 4, 1, 5, 0, 4, -1
1729 const int* thePentaTo3[6] = { thePentaTo3_1, thePentaTo3_2, thePentaTo3_3,
1730 thePentaTo3_4, thePentaTo3_5, thePentaTo3_6 };
1732 // Methods of splitting hexahedron into prisms
1734 const int theHexTo4Prisms_BT[6*4+1] = // bottom-top
1736 0, 1, 8, 4, 5, 9, 1, 2, 8, 5, 6, 9, 2, 3, 8, 6, 7, 9, 3, 0, 8, 7, 4, 9, -1
1738 const int theHexTo4Prisms_LR[6*4+1] = // left-right
1740 1, 0, 8, 2, 3, 9, 0, 4, 8, 3, 7, 9, 4, 5, 8, 7, 6, 9, 5, 1, 8, 6, 2, 9, -1
1742 const int theHexTo4Prisms_FB[6*4+1] = // front-back
1744 0, 3, 9, 1, 2, 8, 3, 7, 9, 2, 6, 8, 7, 4, 9, 6, 5, 8, 4, 0, 9, 5, 1, 8, -1
1747 const int theHexTo2Prisms_BT_1[6*2+1] =
1749 0, 1, 3, 4, 5, 7, 1, 2, 3, 5, 6, 7, -1
1751 const int theHexTo2Prisms_BT_2[6*2+1] =
1753 0, 1, 2, 4, 5, 6, 0, 2, 3, 4, 6, 7, -1
1755 const int* theHexTo2Prisms_BT[2] = { theHexTo2Prisms_BT_1, theHexTo2Prisms_BT_2 };
1757 const int theHexTo2Prisms_LR_1[6*2+1] =
1759 1, 0, 4, 2, 3, 7, 1, 4, 5, 2, 7, 6, -1
1761 const int theHexTo2Prisms_LR_2[6*2+1] =
1763 1, 0, 4, 2, 3, 7, 1, 4, 5, 2, 7, 6, -1
1765 const int* theHexTo2Prisms_LR[2] = { theHexTo2Prisms_LR_1, theHexTo2Prisms_LR_2 };
1767 const int theHexTo2Prisms_FB_1[6*2+1] =
1769 0, 3, 4, 1, 2, 5, 3, 7, 4, 2, 6, 5, -1
1771 const int theHexTo2Prisms_FB_2[6*2+1] =
1773 0, 3, 7, 1, 2, 7, 0, 7, 4, 1, 6, 5, -1
1775 const int* theHexTo2Prisms_FB[2] = { theHexTo2Prisms_FB_1, theHexTo2Prisms_FB_2 };
1778 struct TTriangleFacet //!< stores indices of three nodes of tetra facet
1781 TTriangleFacet(int n1, int n2, int n3): _n1(n1), _n2(n2), _n3(n3) {}
1782 bool contains(int n) const { return ( n == _n1 || n == _n2 || n == _n3 ); }
1783 bool hasAdjacentVol( const SMDS_MeshElement* elem,
1784 const SMDSAbs_GeometryType geom = SMDSGeom_TETRA) const;
1790 const int* _connectivity; //!< foursomes of tetra connectivy finished by -1
1791 bool _baryNode; //!< additional node is to be created at cell barycenter
1792 bool _ownConn; //!< to delete _connectivity in destructor
1793 map<int, const SMDS_MeshNode*> _faceBaryNode; //!< map face index to node at BC of face
1795 TSplitMethod( int nbTet=0, const int* conn=0, bool addNode=false)
1796 : _nbSplits(nbTet), _nbCorners(4), _connectivity(conn), _baryNode(addNode), _ownConn(false) {}
1797 ~TSplitMethod() { if ( _ownConn ) delete [] _connectivity; _connectivity = 0; }
1798 bool hasFacet( const TTriangleFacet& facet ) const
1800 if ( _nbCorners == 4 )
1802 const int* tetConn = _connectivity;
1803 for ( ; tetConn[0] >= 0; tetConn += 4 )
1804 if (( facet.contains( tetConn[0] ) +
1805 facet.contains( tetConn[1] ) +
1806 facet.contains( tetConn[2] ) +
1807 facet.contains( tetConn[3] )) == 3 )
1810 else // prism, _nbCorners == 6
1812 const int* prismConn = _connectivity;
1813 for ( ; prismConn[0] >= 0; prismConn += 6 )
1815 if (( facet.contains( prismConn[0] ) &&
1816 facet.contains( prismConn[1] ) &&
1817 facet.contains( prismConn[2] ))
1819 ( facet.contains( prismConn[3] ) &&
1820 facet.contains( prismConn[4] ) &&
1821 facet.contains( prismConn[5] )))
1829 //=======================================================================
1831 * \brief return TSplitMethod for the given element to split into tetrahedra
1833 //=======================================================================
1835 TSplitMethod getTetraSplitMethod( SMDS_VolumeTool& vol, const int theMethodFlags)
1837 const int iQ = vol.Element()->IsQuadratic() ? 2 : 1;
1839 // at HEXA_TO_24 method, each face of volume is split into triangles each based on
1840 // an edge and a face barycenter; tertaherdons are based on triangles and
1841 // a volume barycenter
1842 const bool is24TetMode = ( theMethodFlags == SMESH_MeshEditor::HEXA_TO_24 );
1844 // Find out how adjacent volumes are split
1846 vector < list< TTriangleFacet > > triaSplitsByFace( vol.NbFaces() ); // splits of each side
1847 int hasAdjacentSplits = 0, maxTetConnSize = 0;
1848 for ( int iF = 0; iF < vol.NbFaces(); ++iF )
1850 int nbNodes = vol.NbFaceNodes( iF ) / iQ;
1851 maxTetConnSize += 4 * ( nbNodes - (is24TetMode ? 0 : 2));
1852 if ( nbNodes < 4 ) continue;
1854 list< TTriangleFacet >& triaSplits = triaSplitsByFace[ iF ];
1855 const int* nInd = vol.GetFaceNodesIndices( iF );
1858 TTriangleFacet t012( nInd[0*iQ], nInd[1*iQ], nInd[2*iQ] );
1859 TTriangleFacet t123( nInd[1*iQ], nInd[2*iQ], nInd[3*iQ] );
1860 if ( t012.hasAdjacentVol( vol.Element() )) triaSplits.push_back( t012 );
1861 else if ( t123.hasAdjacentVol( vol.Element() )) triaSplits.push_back( t123 );
1865 int iCom = 0; // common node of triangle faces to split into
1866 for ( int iVar = 0; iVar < nbNodes; ++iVar, ++iCom )
1868 TTriangleFacet t012( nInd[ iQ * ( iCom )],
1869 nInd[ iQ * ( (iCom+1)%nbNodes )],
1870 nInd[ iQ * ( (iCom+2)%nbNodes )]);
1871 TTriangleFacet t023( nInd[ iQ * ( iCom )],
1872 nInd[ iQ * ( (iCom+2)%nbNodes )],
1873 nInd[ iQ * ( (iCom+3)%nbNodes )]);
1874 if ( t012.hasAdjacentVol( vol.Element() ) && t023.hasAdjacentVol( vol.Element() ))
1876 triaSplits.push_back( t012 );
1877 triaSplits.push_back( t023 );
1882 if ( !triaSplits.empty() )
1883 hasAdjacentSplits = true;
1886 // Among variants of split method select one compliant with adjacent volumes
1888 TSplitMethod method;
1889 if ( !vol.Element()->IsPoly() && !is24TetMode )
1891 int nbVariants = 2, nbTet = 0;
1892 const int** connVariants = 0;
1893 switch ( vol.Element()->GetEntityType() )
1895 case SMDSEntity_Hexa:
1896 case SMDSEntity_Quad_Hexa:
1897 case SMDSEntity_TriQuad_Hexa:
1898 if ( theMethodFlags == SMESH_MeshEditor::HEXA_TO_5 )
1899 connVariants = theHexTo5, nbTet = 5;
1901 connVariants = theHexTo6, nbTet = 6, nbVariants = 4;
1903 case SMDSEntity_Pyramid:
1904 case SMDSEntity_Quad_Pyramid:
1905 connVariants = thePyraTo2; nbTet = 2;
1907 case SMDSEntity_Penta:
1908 case SMDSEntity_Quad_Penta:
1909 connVariants = thePentaTo3; nbTet = 3; nbVariants = 6;
1914 for ( int variant = 0; variant < nbVariants && method._nbSplits == 0; ++variant )
1916 // check method compliancy with adjacent tetras,
1917 // all found splits must be among facets of tetras described by this method
1918 method = TSplitMethod( nbTet, connVariants[variant] );
1919 if ( hasAdjacentSplits && method._nbSplits > 0 )
1921 bool facetCreated = true;
1922 for ( size_t iF = 0; facetCreated && iF < triaSplitsByFace.size(); ++iF )
1924 list< TTriangleFacet >::const_iterator facet = triaSplitsByFace[iF].begin();
1925 for ( ; facetCreated && facet != triaSplitsByFace[iF].end(); ++facet )
1926 facetCreated = method.hasFacet( *facet );
1928 if ( !facetCreated )
1929 method = TSplitMethod(0); // incompatible method
1933 if ( method._nbSplits < 1 )
1935 // No standard method is applicable, use a generic solution:
1936 // each facet of a volume is split into triangles and
1937 // each of triangles and a volume barycenter form a tetrahedron.
1939 const bool isHex27 = ( vol.Element()->GetEntityType() == SMDSEntity_TriQuad_Hexa );
1941 int* connectivity = new int[ maxTetConnSize + 1 ];
1942 method._connectivity = connectivity;
1943 method._ownConn = true;
1944 method._baryNode = !isHex27; // to create central node or not
1947 int baryCenInd = vol.NbNodes() - int( isHex27 );
1948 for ( int iF = 0; iF < vol.NbFaces(); ++iF )
1950 const int nbNodes = vol.NbFaceNodes( iF ) / iQ;
1951 const int* nInd = vol.GetFaceNodesIndices( iF );
1952 // find common node of triangle facets of tetra to create
1953 int iCommon = 0; // index in linear numeration
1954 const list< TTriangleFacet >& triaSplits = triaSplitsByFace[ iF ];
1955 if ( !triaSplits.empty() )
1958 const TTriangleFacet* facet = &triaSplits.front();
1959 for ( ; iCommon < nbNodes-1 ; ++iCommon )
1960 if ( facet->contains( nInd[ iQ * iCommon ]) &&
1961 facet->contains( nInd[ iQ * ((iCommon+2)%nbNodes) ]))
1964 else if ( nbNodes > 3 && !is24TetMode )
1966 // find the best method of splitting into triangles by aspect ratio
1967 SMESH::Controls::NumericalFunctorPtr aspectRatio( new SMESH::Controls::AspectRatio);
1968 map< double, int > badness2iCommon;
1969 const SMDS_MeshNode** nodes = vol.GetFaceNodes( iF );
1970 int nbVariants = ( nbNodes == 4 ? 2 : nbNodes );
1971 for ( int iVar = 0; iVar < nbVariants; ++iVar, ++iCommon )
1974 for ( int iLast = iCommon+2; iLast < iCommon+nbNodes; ++iLast )
1976 SMDS_FaceOfNodes tria ( nodes[ iQ*( iCommon )],
1977 nodes[ iQ*((iLast-1)%nbNodes)],
1978 nodes[ iQ*((iLast )%nbNodes)]);
1979 badness += getBadRate( &tria, aspectRatio );
1981 badness2iCommon.insert( make_pair( badness, iCommon ));
1983 // use iCommon with lowest badness
1984 iCommon = badness2iCommon.begin()->second;
1986 if ( iCommon >= nbNodes )
1987 iCommon = 0; // something wrong
1989 // fill connectivity of tetrahedra based on a current face
1990 int nbTet = nbNodes - 2;
1991 if ( is24TetMode && nbNodes > 3 && triaSplits.empty())
1996 faceBaryCenInd = vol.GetCenterNodeIndex( iF );
1997 method._faceBaryNode[ iF ] = vol.GetNodes()[ faceBaryCenInd ];
2001 method._faceBaryNode[ iF ] = 0;
2002 faceBaryCenInd = baryCenInd + method._faceBaryNode.size();
2005 for ( int i = 0; i < nbTet; ++i )
2007 int i1 = i, i2 = (i+1) % nbNodes;
2008 if ( !vol.IsFaceExternal( iF )) swap( i1, i2 );
2009 connectivity[ connSize++ ] = nInd[ iQ * i1 ];
2010 connectivity[ connSize++ ] = nInd[ iQ * i2 ];
2011 connectivity[ connSize++ ] = faceBaryCenInd;
2012 connectivity[ connSize++ ] = baryCenInd;
2017 for ( int i = 0; i < nbTet; ++i )
2019 int i1 = (iCommon+1+i) % nbNodes, i2 = (iCommon+2+i) % nbNodes;
2020 if ( !vol.IsFaceExternal( iF )) swap( i1, i2 );
2021 connectivity[ connSize++ ] = nInd[ iQ * iCommon ];
2022 connectivity[ connSize++ ] = nInd[ iQ * i1 ];
2023 connectivity[ connSize++ ] = nInd[ iQ * i2 ];
2024 connectivity[ connSize++ ] = baryCenInd;
2027 method._nbSplits += nbTet;
2029 } // loop on volume faces
2031 connectivity[ connSize++ ] = -1;
2033 } // end of generic solution
2037 //=======================================================================
2039 * \brief return TSplitMethod to split haxhedron into prisms
2041 //=======================================================================
2043 TSplitMethod getPrismSplitMethod( SMDS_VolumeTool& vol,
2044 const int methodFlags,
2045 const int facetToSplit)
2047 // order of facets in HEX according to SMDS_VolumeTool::Hexa_F :
2049 const int iF = ( facetToSplit < 2 ) ? 0 : 1 + ( facetToSplit-2 ) % 2; // [0,1,2]
2051 if ( methodFlags == SMESH_MeshEditor::HEXA_TO_4_PRISMS )
2053 static TSplitMethod to4methods[4]; // order BT, LR, FB
2054 if ( to4methods[iF]._nbSplits == 0 )
2058 to4methods[iF]._connectivity = theHexTo4Prisms_BT;
2059 to4methods[iF]._faceBaryNode[ 0 ] = 0;
2060 to4methods[iF]._faceBaryNode[ 1 ] = 0;
2063 to4methods[iF]._connectivity = theHexTo4Prisms_LR;
2064 to4methods[iF]._faceBaryNode[ 2 ] = 0;
2065 to4methods[iF]._faceBaryNode[ 4 ] = 0;
2068 to4methods[iF]._connectivity = theHexTo4Prisms_FB;
2069 to4methods[iF]._faceBaryNode[ 3 ] = 0;
2070 to4methods[iF]._faceBaryNode[ 5 ] = 0;
2072 default: return to4methods[3];
2074 to4methods[iF]._nbSplits = 4;
2075 to4methods[iF]._nbCorners = 6;
2077 return to4methods[iF];
2079 // else if ( methodFlags == HEXA_TO_2_PRISMS )
2081 TSplitMethod method;
2083 const int iQ = vol.Element()->IsQuadratic() ? 2 : 1;
2085 const int nbVariants = 2, nbSplits = 2;
2086 const int** connVariants = 0;
2088 case 0: connVariants = theHexTo2Prisms_BT; break;
2089 case 1: connVariants = theHexTo2Prisms_LR; break;
2090 case 2: connVariants = theHexTo2Prisms_FB; break;
2091 default: return method;
2094 // look for prisms adjacent via facetToSplit and an opposite one
2095 for ( int is2nd = 0; is2nd < 2; ++is2nd )
2097 int iFacet = is2nd ? vol.GetOppFaceIndexOfHex( facetToSplit ) : facetToSplit;
2098 int nbNodes = vol.NbFaceNodes( iFacet ) / iQ;
2099 if ( nbNodes != 4 ) return method;
2101 const int* nInd = vol.GetFaceNodesIndices( iFacet );
2102 TTriangleFacet t012( nInd[0*iQ], nInd[1*iQ], nInd[2*iQ] );
2103 TTriangleFacet t123( nInd[1*iQ], nInd[2*iQ], nInd[3*iQ] );
2105 if ( t012.hasAdjacentVol( vol.Element(), SMDSGeom_PENTA ))
2107 else if ( t123.hasAdjacentVol( vol.Element(), SMDSGeom_PENTA ))
2112 // there are adjacent prism
2113 for ( int variant = 0; variant < nbVariants; ++variant )
2115 // check method compliancy with adjacent prisms,
2116 // the found prism facets must be among facets of prisms described by current method
2117 method._nbSplits = nbSplits;
2118 method._nbCorners = 6;
2119 method._connectivity = connVariants[ variant ];
2120 if ( method.hasFacet( *t ))
2125 // No adjacent prisms. Select a variant with a best aspect ratio.
2127 double badness[2] = { 0, 0 };
2128 static SMESH::Controls::NumericalFunctorPtr aspectRatio( new SMESH::Controls::AspectRatio);
2129 const SMDS_MeshNode** nodes = vol.GetNodes();
2130 for ( int variant = 0; variant < nbVariants; ++variant )
2131 for ( int is2nd = 0; is2nd < 2; ++is2nd )
2133 int iFacet = is2nd ? vol.GetOppFaceIndexOfHex( facetToSplit ) : facetToSplit;
2134 const int* nInd = vol.GetFaceNodesIndices( iFacet );
2136 method._connectivity = connVariants[ variant ];
2137 TTriangleFacet t012( nInd[0*iQ], nInd[1*iQ], nInd[2*iQ] );
2138 TTriangleFacet t123( nInd[1*iQ], nInd[2*iQ], nInd[3*iQ] );
2139 TTriangleFacet* t = ( method.hasFacet( t012 )) ? & t012 : & t123;
2141 SMDS_FaceOfNodes tria ( nodes[ t->_n1 ],
2144 badness[ variant ] += getBadRate( &tria, aspectRatio );
2146 const int iBetter = ( badness[1] < badness[0] && badness[0]-badness[1] > 0.1 * badness[0] );
2148 method._nbSplits = nbSplits;
2149 method._nbCorners = 6;
2150 method._connectivity = connVariants[ iBetter ];
2155 //================================================================================
2157 * \brief Check if there is a tetraherdon adjacent to the given element via this facet
2159 //================================================================================
2161 bool TTriangleFacet::hasAdjacentVol( const SMDS_MeshElement* elem,
2162 const SMDSAbs_GeometryType geom ) const
2164 // find the tetrahedron including the three nodes of facet
2165 const SMDS_MeshNode* n1 = elem->GetNode(_n1);
2166 const SMDS_MeshNode* n2 = elem->GetNode(_n2);
2167 const SMDS_MeshNode* n3 = elem->GetNode(_n3);
2168 SMDS_ElemIteratorPtr volIt1 = n1->GetInverseElementIterator(SMDSAbs_Volume);
2169 while ( volIt1->more() )
2171 const SMDS_MeshElement* v = volIt1->next();
2172 if ( v->GetGeomType() != geom )
2174 const int lastCornerInd = v->NbCornerNodes() - 1;
2175 if ( v->IsQuadratic() && v->GetNodeIndex( n1 ) > lastCornerInd )
2176 continue; // medium node not allowed
2177 const int ind2 = v->GetNodeIndex( n2 );
2178 if ( ind2 < 0 || lastCornerInd < ind2 )
2180 const int ind3 = v->GetNodeIndex( n3 );
2181 if ( ind3 < 0 || lastCornerInd < ind3 )
2188 //=======================================================================
2190 * \brief A key of a face of volume
2192 //=======================================================================
2194 struct TVolumeFaceKey: pair< pair< int, int>, pair< int, int> >
2196 TVolumeFaceKey( SMDS_VolumeTool& vol, int iF )
2198 TIDSortedNodeSet sortedNodes;
2199 const int iQ = vol.Element()->IsQuadratic() ? 2 : 1;
2200 int nbNodes = vol.NbFaceNodes( iF );
2201 const SMDS_MeshNode** fNodes = vol.GetFaceNodes( iF );
2202 for ( int i = 0; i < nbNodes; i += iQ )
2203 sortedNodes.insert( fNodes[i] );
2204 TIDSortedNodeSet::iterator n = sortedNodes.begin();
2205 first.first = (*(n++))->GetID();
2206 first.second = (*(n++))->GetID();
2207 second.first = (*(n++))->GetID();
2208 second.second = ( sortedNodes.size() > 3 ) ? (*(n++))->GetID() : 0;
2213 //=======================================================================
2214 //function : SplitVolumes
2215 //purpose : Split volume elements into tetrahedra or prisms.
2216 // If facet ID < 0, element is split into tetrahedra,
2217 // else a hexahedron is split into prisms so that the given facet is
2218 // split into triangles
2219 //=======================================================================
2221 void SMESH_MeshEditor::SplitVolumes (const TFacetOfElem & theElems,
2222 const int theMethodFlags)
2224 SMDS_VolumeTool volTool;
2225 SMESH_MesherHelper helper( *GetMesh()), fHelper(*GetMesh());
2226 fHelper.ToFixNodeParameters( true );
2228 SMESHDS_SubMesh* subMesh = 0;//GetMeshDS()->MeshElements(1);
2229 SMESHDS_SubMesh* fSubMesh = 0;//subMesh;
2231 SMESH_SequenceOfElemPtr newNodes, newElems;
2233 // map face of volume to it's baricenrtic node
2234 map< TVolumeFaceKey, const SMDS_MeshNode* > volFace2BaryNode;
2236 vector<const SMDS_MeshElement* > splitVols;
2238 TFacetOfElem::const_iterator elem2facet = theElems.begin();
2239 for ( ; elem2facet != theElems.end(); ++elem2facet )
2241 const SMDS_MeshElement* elem = elem2facet->first;
2242 const int facetToSplit = elem2facet->second;
2243 if ( elem->GetType() != SMDSAbs_Volume )
2245 const SMDSAbs_EntityType geomType = elem->GetEntityType();
2246 if ( geomType == SMDSEntity_Tetra || geomType == SMDSEntity_Quad_Tetra )
2249 if ( !volTool.Set( elem, /*ignoreCentralNodes=*/false )) continue; // strange...
2251 TSplitMethod splitMethod = ( facetToSplit < 0 ?
2252 getTetraSplitMethod( volTool, theMethodFlags ) :
2253 getPrismSplitMethod( volTool, theMethodFlags, facetToSplit ));
2254 if ( splitMethod._nbSplits < 1 ) continue;
2256 // find submesh to add new tetras to
2257 if ( !subMesh || !subMesh->Contains( elem ))
2259 int shapeID = FindShape( elem );
2260 helper.SetSubShape( shapeID ); // helper will add tetras to the found submesh
2261 subMesh = GetMeshDS()->MeshElements( shapeID );
2264 if ( elem->IsQuadratic() )
2267 // add quadratic links to the helper
2268 for ( int iF = 0; iF < volTool.NbFaces(); ++iF )
2270 const SMDS_MeshNode** fNodes = volTool.GetFaceNodes( iF );
2271 int nbN = volTool.NbFaceNodes( iF ) - bool( volTool.GetCenterNodeIndex(iF) > 0 );
2272 for ( int iN = 0; iN < nbN; iN += iQ )
2273 helper.AddTLinkNode( fNodes[iN], fNodes[iN+2], fNodes[iN+1] );
2275 helper.SetIsQuadratic( true );
2280 helper.SetIsQuadratic( false );
2282 vector<const SMDS_MeshNode*> nodes( volTool.GetNodes(),
2283 volTool.GetNodes() + elem->NbNodes() );
2284 helper.SetElementsOnShape( true );
2285 if ( splitMethod._baryNode )
2287 // make a node at barycenter
2288 volTool.GetBaryCenter( bc[0], bc[1], bc[2] );
2289 SMDS_MeshNode* gcNode = helper.AddNode( bc[0], bc[1], bc[2] );
2290 nodes.push_back( gcNode );
2291 newNodes.Append( gcNode );
2293 if ( !splitMethod._faceBaryNode.empty() )
2295 // make or find baricentric nodes of faces
2296 map<int, const SMDS_MeshNode*>::iterator iF_n = splitMethod._faceBaryNode.begin();
2297 for ( ; iF_n != splitMethod._faceBaryNode.end(); ++iF_n )
2299 map< TVolumeFaceKey, const SMDS_MeshNode* >::iterator f_n =
2300 volFace2BaryNode.insert
2301 ( make_pair( TVolumeFaceKey( volTool,iF_n->first ), iF_n->second )).first;
2304 volTool.GetFaceBaryCenter( iF_n->first, bc[0], bc[1], bc[2] );
2305 newNodes.Append( f_n->second = helper.AddNode( bc[0], bc[1], bc[2] ));
2307 nodes.push_back( iF_n->second = f_n->second );
2312 splitVols.resize( splitMethod._nbSplits ); // splits of a volume
2313 const int* volConn = splitMethod._connectivity;
2314 if ( splitMethod._nbCorners == 4 ) // tetra
2315 for ( int i = 0; i < splitMethod._nbSplits; ++i, volConn += splitMethod._nbCorners )
2316 newElems.Append( splitVols[ i ] = helper.AddVolume( nodes[ volConn[0] ],
2317 nodes[ volConn[1] ],
2318 nodes[ volConn[2] ],
2319 nodes[ volConn[3] ]));
2321 for ( int i = 0; i < splitMethod._nbSplits; ++i, volConn += splitMethod._nbCorners )
2322 newElems.Append( splitVols[ i ] = helper.AddVolume( nodes[ volConn[0] ],
2323 nodes[ volConn[1] ],
2324 nodes[ volConn[2] ],
2325 nodes[ volConn[3] ],
2326 nodes[ volConn[4] ],
2327 nodes[ volConn[5] ]));
2329 ReplaceElemInGroups( elem, splitVols, GetMeshDS() );
2331 // Split faces on sides of the split volume
2333 const SMDS_MeshNode** volNodes = volTool.GetNodes();
2334 for ( int iF = 0; iF < volTool.NbFaces(); ++iF )
2336 const int nbNodes = volTool.NbFaceNodes( iF ) / iQ;
2337 if ( nbNodes < 4 ) continue;
2339 // find an existing face
2340 vector<const SMDS_MeshNode*> fNodes( volTool.GetFaceNodes( iF ),
2341 volTool.GetFaceNodes( iF ) + volTool.NbFaceNodes( iF ));
2342 while ( const SMDS_MeshElement* face = GetMeshDS()->FindElement( fNodes, SMDSAbs_Face,
2343 /*noMedium=*/false))
2346 helper.SetElementsOnShape( false );
2347 vector< const SMDS_MeshElement* > triangles;
2349 // find submesh to add new triangles in
2350 if ( !fSubMesh || !fSubMesh->Contains( face ))
2352 int shapeID = FindShape( face );
2353 fSubMesh = GetMeshDS()->MeshElements( shapeID );
2355 map<int, const SMDS_MeshNode*>::iterator iF_n = splitMethod._faceBaryNode.find(iF);
2356 if ( iF_n != splitMethod._faceBaryNode.end() )
2358 const SMDS_MeshNode *baryNode = iF_n->second;
2359 for ( int iN = 0; iN < nbNodes*iQ; iN += iQ )
2361 const SMDS_MeshNode* n1 = fNodes[iN];
2362 const SMDS_MeshNode *n2 = fNodes[(iN+iQ)%(nbNodes*iQ)];
2363 const SMDS_MeshNode *n3 = baryNode;
2364 if ( !volTool.IsFaceExternal( iF ))
2366 triangles.push_back( helper.AddFace( n1,n2,n3 ));
2368 if ( fSubMesh ) // update position of the bary node on geometry
2371 subMesh->RemoveNode( baryNode, false );
2372 GetMeshDS()->SetNodeOnFace( baryNode, fSubMesh->GetID() );
2373 const TopoDS_Shape& s = GetMeshDS()->IndexToShape( fSubMesh->GetID() );
2374 if ( !s.IsNull() && s.ShapeType() == TopAbs_FACE )
2376 fHelper.SetSubShape( s );
2377 gp_XY uv( 1e100, 1e100 );
2379 if ( !fHelper.CheckNodeUV( TopoDS::Face( s ), baryNode,
2380 uv, /*tol=*/1e-7, /*force=*/true, distXYZ ) &&
2383 // node is too far from the surface
2384 GetMeshDS()->MoveNode( baryNode, distXYZ[1], distXYZ[2], distXYZ[3] );
2385 const_cast<SMDS_MeshNode*>( baryNode )->SetPosition
2386 ( SMDS_PositionPtr( new SMDS_FacePosition( uv.X(), uv.Y() )));
2393 // among possible triangles create ones discribed by split method
2394 const int* nInd = volTool.GetFaceNodesIndices( iF );
2395 int nbVariants = ( nbNodes == 4 ? 2 : nbNodes );
2396 int iCom = 0; // common node of triangle faces to split into
2397 list< TTriangleFacet > facets;
2398 for ( int iVar = 0; iVar < nbVariants; ++iVar, ++iCom )
2400 TTriangleFacet t012( nInd[ iQ * ( iCom )],
2401 nInd[ iQ * ( (iCom+1)%nbNodes )],
2402 nInd[ iQ * ( (iCom+2)%nbNodes )]);
2403 TTriangleFacet t023( nInd[ iQ * ( iCom )],
2404 nInd[ iQ * ( (iCom+2)%nbNodes )],
2405 nInd[ iQ * ( (iCom+3)%nbNodes )]);
2406 if ( splitMethod.hasFacet( t012 ) && splitMethod.hasFacet( t023 ))
2408 facets.push_back( t012 );
2409 facets.push_back( t023 );
2410 for ( int iLast = iCom+4; iLast < iCom+nbNodes; ++iLast )
2411 facets.push_back( TTriangleFacet( nInd[ iQ * ( iCom )],
2412 nInd[ iQ * ((iLast-1)%nbNodes )],
2413 nInd[ iQ * ((iLast )%nbNodes )]));
2417 list< TTriangleFacet >::iterator facet = facets.begin();
2418 if ( facet == facets.end() )
2420 for ( ; facet != facets.end(); ++facet )
2422 if ( !volTool.IsFaceExternal( iF ))
2423 swap( facet->_n2, facet->_n3 );
2424 triangles.push_back( helper.AddFace( volNodes[ facet->_n1 ],
2425 volNodes[ facet->_n2 ],
2426 volNodes[ facet->_n3 ]));
2429 for ( size_t i = 0; i < triangles.size(); ++i )
2431 if ( !triangles[ i ]) continue;
2433 fSubMesh->AddElement( triangles[ i ]);
2434 newElems.Append( triangles[ i ]);
2436 ReplaceElemInGroups( face, triangles, GetMeshDS() );
2437 GetMeshDS()->RemoveFreeElement( face, fSubMesh, /*fromGroups=*/false );
2439 } // while a face based on facet nodes exists
2440 } // loop on volume faces to split them into triangles
2442 GetMeshDS()->RemoveFreeElement( elem, subMesh, /*fromGroups=*/false );
2444 if ( geomType == SMDSEntity_TriQuad_Hexa )
2446 // remove medium nodes that could become free
2447 for ( int i = 20; i < volTool.NbNodes(); ++i )
2448 if ( volNodes[i]->NbInverseElements() == 0 )
2449 GetMeshDS()->RemoveNode( volNodes[i] );
2451 } // loop on volumes to split
2453 myLastCreatedNodes = newNodes;
2454 myLastCreatedElems = newElems;
2457 //=======================================================================
2458 //function : GetHexaFacetsToSplit
2459 //purpose : For hexahedra that will be split into prisms, finds facets to
2460 // split into triangles. Only hexahedra adjacent to the one closest
2461 // to theFacetNormal.Location() are returned.
2462 //param [in,out] theHexas - the hexahedra
2463 //param [in] theFacetNormal - facet normal
2464 //param [out] theFacets - the hexahedra and found facet IDs
2465 //=======================================================================
2467 void SMESH_MeshEditor::GetHexaFacetsToSplit( TIDSortedElemSet& theHexas,
2468 const gp_Ax1& theFacetNormal,
2469 TFacetOfElem & theFacets)
2471 #define THIS_METHOD "SMESH_MeshEditor::GetHexaFacetsToSplit(): "
2473 // Find a hexa closest to the location of theFacetNormal
2475 const SMDS_MeshElement* startHex;
2477 // get SMDS_ElemIteratorPtr on theHexas
2478 typedef const SMDS_MeshElement* TValue;
2479 typedef TIDSortedElemSet::iterator TSetIterator;
2480 typedef SMDS::SimpleAccessor<TValue,TSetIterator> TAccesor;
2481 typedef SMDS_MeshElement::GeomFilter TFilter;
2482 typedef SMDS_SetIterator < TValue, TSetIterator, TAccesor, TFilter > TElemSetIter;
2483 SMDS_ElemIteratorPtr elemIt = SMDS_ElemIteratorPtr
2484 ( new TElemSetIter( theHexas.begin(),
2486 SMDS_MeshElement::GeomFilter( SMDSGeom_HEXA )));
2488 SMESH_ElementSearcher* searcher =
2489 SMESH_MeshAlgos::GetElementSearcher( *myMesh->GetMeshDS(), elemIt );
2491 startHex = searcher->FindClosestTo( theFacetNormal.Location(), SMDSAbs_Volume );
2496 throw SALOME_Exception( THIS_METHOD "startHex not found");
2499 // Select a facet of startHex by theFacetNormal
2501 SMDS_VolumeTool vTool( startHex );
2502 double norm[3], dot, maxDot = 0;
2504 for ( int iF = 0; iF < vTool.NbFaces(); ++iF )
2505 if ( vTool.GetFaceNormal( iF, norm[0], norm[1], norm[2] ))
2507 dot = Abs( theFacetNormal.Direction().Dot( gp_Dir( norm[0], norm[1], norm[2] )));
2515 throw SALOME_Exception( THIS_METHOD "facet of startHex not found");
2517 // Fill theFacets starting from facetID of startHex
2519 // facets used for seach of volumes adjacent to already treated ones
2520 typedef pair< TFacetOfElem::iterator, int > TElemFacets;
2521 typedef map< TVolumeFaceKey, TElemFacets > TFacetMap;
2522 TFacetMap facetsToCheck;
2524 set<const SMDS_MeshNode*> facetNodes;
2525 const SMDS_MeshElement* curHex;
2527 const bool allHex = ((int) theHexas.size() == myMesh->NbHexas() );
2531 // move in two directions from startHex via facetID
2532 for ( int is2nd = 0; is2nd < 2; ++is2nd )
2535 int curFacet = facetID;
2536 if ( is2nd ) // do not treat startHex twice
2538 vTool.Set( curHex );
2539 if ( vTool.IsFreeFace( curFacet, &curHex ))
2545 vTool.GetFaceNodes( curFacet, facetNodes );
2546 vTool.Set( curHex );
2547 curFacet = vTool.GetFaceIndex( facetNodes );
2552 // store a facet to split
2553 if ( curHex->GetGeomType() != SMDSGeom_HEXA )
2555 theFacets.insert( make_pair( curHex, -1 ));
2558 if ( !allHex && !theHexas.count( curHex ))
2561 pair< TFacetOfElem::iterator, bool > facetIt2isNew =
2562 theFacets.insert( make_pair( curHex, curFacet ));
2563 if ( !facetIt2isNew.second )
2566 // remember not-to-split facets in facetsToCheck
2567 int oppFacet = vTool.GetOppFaceIndexOfHex( curFacet );
2568 for ( int iF = 0; iF < vTool.NbFaces(); ++iF )
2570 if ( iF == curFacet && iF == oppFacet )
2572 TVolumeFaceKey facetKey ( vTool, iF );
2573 TElemFacets elemFacet( facetIt2isNew.first, iF );
2574 pair< TFacetMap::iterator, bool > it2isnew =
2575 facetsToCheck.insert( make_pair( facetKey, elemFacet ));
2576 if ( !it2isnew.second )
2577 facetsToCheck.erase( it2isnew.first ); // adjacent hex already checked
2579 // pass to a volume adjacent via oppFacet
2580 if ( vTool.IsFreeFace( oppFacet, &curHex ))
2586 // get a new curFacet
2587 vTool.GetFaceNodes( oppFacet, facetNodes );
2588 vTool.Set( curHex );
2589 curFacet = vTool.GetFaceIndex( facetNodes, /*hint=*/curFacet );
2592 } // move in two directions from startHex via facetID
2594 // Find a new startHex by facetsToCheck
2598 TFacetMap::iterator fIt = facetsToCheck.begin();
2599 while ( !startHex && fIt != facetsToCheck.end() )
2601 const TElemFacets& elemFacets = fIt->second;
2602 const SMDS_MeshElement* hex = elemFacets.first->first;
2603 int splitFacet = elemFacets.first->second;
2604 int lateralFacet = elemFacets.second;
2605 facetsToCheck.erase( fIt );
2606 fIt = facetsToCheck.begin();
2609 if ( vTool.IsFreeFace( lateralFacet, &curHex ) ||
2610 curHex->GetGeomType() != SMDSGeom_HEXA )
2612 if ( !allHex && !theHexas.count( curHex ))
2617 // find a facet of startHex to split
2619 set<const SMDS_MeshNode*> lateralNodes;
2620 vTool.GetFaceNodes( lateralFacet, lateralNodes );
2621 vTool.GetFaceNodes( splitFacet, facetNodes );
2622 int oppLateralFacet = vTool.GetOppFaceIndexOfHex( lateralFacet );
2623 vTool.Set( startHex );
2624 lateralFacet = vTool.GetFaceIndex( lateralNodes, oppLateralFacet );
2626 // look for a facet of startHex having common nodes with facetNodes
2627 // but not lateralFacet
2628 for ( int iF = 0; iF < vTool.NbFaces(); ++iF )
2630 if ( iF == lateralFacet )
2632 int nbCommonNodes = 0;
2633 const SMDS_MeshNode** nn = vTool.GetFaceNodes( iF );
2634 for ( int iN = 0, nbN = vTool.NbFaceNodes( iF ); iN < nbN; ++iN )
2635 nbCommonNodes += facetNodes.count( nn[ iN ]);
2637 if ( nbCommonNodes >= 2 )
2644 throw SALOME_Exception( THIS_METHOD "facet of a new startHex not found");
2646 } // while ( startHex )
2653 //================================================================================
2655 * \brief Selects nodes of several elements according to a given interlace
2656 * \param [in] srcNodes - nodes to select from
2657 * \param [out] tgtNodesVec - array of nodes of several elements to fill in
2658 * \param [in] interlace - indices of nodes for all elements
2659 * \param [in] nbElems - nb of elements
2660 * \param [in] nbNodes - nb of nodes in each element
2661 * \param [in] mesh - the mesh
2662 * \param [out] elemQueue - a list to push elements found by the selected nodes
2663 * \param [in] type - type of elements to look for
2665 //================================================================================
2667 void selectNodes( const vector< const SMDS_MeshNode* >& srcNodes,
2668 vector< const SMDS_MeshNode* >* tgtNodesVec,
2669 const int* interlace,
2672 SMESHDS_Mesh* mesh = 0,
2673 list< const SMDS_MeshElement* >* elemQueue=0,
2674 SMDSAbs_ElementType type=SMDSAbs_All)
2676 for ( int iE = 0; iE < nbElems; ++iE )
2678 vector< const SMDS_MeshNode* >& elemNodes = tgtNodesVec[iE];
2679 const int* select = & interlace[iE*nbNodes];
2680 elemNodes.resize( nbNodes );
2681 for ( int iN = 0; iN < nbNodes; ++iN )
2682 elemNodes[iN] = srcNodes[ select[ iN ]];
2684 const SMDS_MeshElement* e;
2686 for ( int iE = 0; iE < nbElems; ++iE )
2687 if (( e = mesh->FindElement( tgtNodesVec[iE], type, /*noMedium=*/false)))
2688 elemQueue->push_back( e );
2692 //=======================================================================
2694 * Split bi-quadratic elements into linear ones without creation of additional nodes
2695 * - bi-quadratic triangle will be split into 3 linear quadrangles;
2696 * - bi-quadratic quadrangle will be split into 4 linear quadrangles;
2697 * - tri-quadratic hexahedron will be split into 8 linear hexahedra;
2698 * Quadratic elements of lower dimension adjacent to the split bi-quadratic element
2699 * will be split in order to keep the mesh conformal.
2700 * \param elems - elements to split
2702 //=======================================================================
2704 void SMESH_MeshEditor::SplitBiQuadraticIntoLinear(TIDSortedElemSet& theElems)
2706 vector< const SMDS_MeshNode* > elemNodes(27), subNodes[12], splitNodes[8];
2707 vector<const SMDS_MeshElement* > splitElems;
2708 list< const SMDS_MeshElement* > elemQueue;
2709 list< const SMDS_MeshElement* >::iterator elemIt;
2711 SMESHDS_Mesh * mesh = GetMeshDS();
2712 ElemFeatures *elemType, hexaType(SMDSAbs_Volume), quadType(SMDSAbs_Face), segType(SMDSAbs_Edge);
2713 int nbElems, nbNodes;
2715 TIDSortedElemSet::iterator elemSetIt = theElems.begin();
2716 for ( ; elemSetIt != theElems.end(); ++elemSetIt )
2719 elemQueue.push_back( *elemSetIt );
2720 for ( elemIt = elemQueue.begin(); elemIt != elemQueue.end(); ++elemIt )
2722 const SMDS_MeshElement* elem = *elemIt;
2723 switch( elem->GetEntityType() )
2725 case SMDSEntity_TriQuad_Hexa: // HEX27
2727 elemNodes.assign( elem->begin_nodes(), elem->end_nodes() );
2728 nbElems = nbNodes = 8;
2729 elemType = & hexaType;
2731 // get nodes for new elements
2732 static int vInd[8][8] = {{ 0,8,20,11, 16,21,26,24 },
2733 { 1,9,20,8, 17,22,26,21 },
2734 { 2,10,20,9, 18,23,26,22 },
2735 { 3,11,20,10, 19,24,26,23 },
2736 { 16,21,26,24, 4,12,25,15 },
2737 { 17,22,26,21, 5,13,25,12 },
2738 { 18,23,26,22, 6,14,25,13 },
2739 { 19,24,26,23, 7,15,25,14 }};
2740 selectNodes( elemNodes, & splitNodes[0], &vInd[0][0], nbElems, nbNodes );
2742 // add boundary faces to elemQueue
2743 static int fInd[6][9] = {{ 0,1,2,3, 8,9,10,11, 20 },
2744 { 4,5,6,7, 12,13,14,15, 25 },
2745 { 0,1,5,4, 8,17,12,16, 21 },
2746 { 1,2,6,5, 9,18,13,17, 22 },
2747 { 2,3,7,6, 10,19,14,18, 23 },
2748 { 3,0,4,7, 11,16,15,19, 24 }};
2749 selectNodes( elemNodes, & subNodes[0], &fInd[0][0], 6,9, mesh, &elemQueue, SMDSAbs_Face );
2751 // add boundary segments to elemQueue
2752 static int eInd[12][3] = {{ 0,1,8 }, { 1,2,9 }, { 2,3,10 }, { 3,0,11 },
2753 { 4,5,12}, { 5,6,13}, { 6,7,14 }, { 7,4,15 },
2754 { 0,4,16}, { 1,5,17}, { 2,6,18 }, { 3,7,19 }};
2755 selectNodes( elemNodes, & subNodes[0], &eInd[0][0], 12,3, mesh, &elemQueue, SMDSAbs_Edge );
2758 case SMDSEntity_BiQuad_Triangle: // TRIA7
2760 elemNodes.assign( elem->begin_nodes(), elem->end_nodes() );
2763 elemType = & quadType;
2765 // get nodes for new elements
2766 static int fInd[3][4] = {{ 0,3,6,5 }, { 1,4,6,3 }, { 2,5,6,4 }};
2767 selectNodes( elemNodes, & splitNodes[0], &fInd[0][0], nbElems, nbNodes );
2769 // add boundary segments to elemQueue
2770 static int eInd[3][3] = {{ 0,1,3 }, { 1,2,4 }, { 2,0,5 }};
2771 selectNodes( elemNodes, & subNodes[0], &eInd[0][0], 3,3, mesh, &elemQueue, SMDSAbs_Edge );
2774 case SMDSEntity_BiQuad_Quadrangle: // QUAD9
2776 elemNodes.assign( elem->begin_nodes(), elem->end_nodes() );
2779 elemType = & quadType;
2781 // get nodes for new elements
2782 static int fInd[4][4] = {{ 0,4,8,7 }, { 1,5,8,4 }, { 2,6,8,5 }, { 3,7,8,6 }};
2783 selectNodes( elemNodes, & splitNodes[0], &fInd[0][0], nbElems, nbNodes );
2785 // add boundary segments to elemQueue
2786 static int eInd[4][3] = {{ 0,1,4 }, { 1,2,5 }, { 2,3,6 }, { 3,0,7 }};
2787 selectNodes( elemNodes, & subNodes[0], &eInd[0][0], 4,3, mesh, &elemQueue, SMDSAbs_Edge );
2790 case SMDSEntity_Quad_Edge:
2792 if ( elemIt == elemQueue.begin() )
2793 continue; // an elem is in theElems
2794 elemNodes.assign( elem->begin_nodes(), elem->end_nodes() );
2797 elemType = & segType;
2799 // get nodes for new elements
2800 static int eInd[2][2] = {{ 0,2 }, { 2,1 }};
2801 selectNodes( elemNodes, & splitNodes[0], &eInd[0][0], nbElems, nbNodes );
2805 } // switch( elem->GetEntityType() )
2807 // Create new elements
2809 SMESHDS_SubMesh* subMesh = mesh->MeshElements( elem->getshapeId() );
2813 //elemType->SetID( elem->GetID() ); // create an elem with the same ID as a removed one
2814 mesh->RemoveFreeElement( elem, subMesh, /*fromGroups=*/false );
2815 //splitElems.push_back( AddElement( splitNodes[ 0 ], *elemType ));
2816 //elemType->SetID( -1 );
2818 for ( int iE = 0; iE < nbElems; ++iE )
2819 splitElems.push_back( AddElement( splitNodes[ iE ], *elemType ));
2822 ReplaceElemInGroups( elem, splitElems, mesh );
2825 for ( size_t i = 0; i < splitElems.size(); ++i )
2826 subMesh->AddElement( splitElems[i] );
2831 //=======================================================================
2832 //function : AddToSameGroups
2833 //purpose : add elemToAdd to the groups the elemInGroups belongs to
2834 //=======================================================================
2836 void SMESH_MeshEditor::AddToSameGroups (const SMDS_MeshElement* elemToAdd,
2837 const SMDS_MeshElement* elemInGroups,
2838 SMESHDS_Mesh * aMesh)
2840 const set<SMESHDS_GroupBase*>& groups = aMesh->GetGroups();
2841 if (!groups.empty()) {
2842 set<SMESHDS_GroupBase*>::const_iterator grIt = groups.begin();
2843 for ( ; grIt != groups.end(); grIt++ ) {
2844 SMESHDS_Group* group = dynamic_cast<SMESHDS_Group*>( *grIt );
2845 if ( group && group->Contains( elemInGroups ))
2846 group->SMDSGroup().Add( elemToAdd );
2852 //=======================================================================
2853 //function : RemoveElemFromGroups
2854 //purpose : Remove removeelem to the groups the elemInGroups belongs to
2855 //=======================================================================
2856 void SMESH_MeshEditor::RemoveElemFromGroups (const SMDS_MeshElement* removeelem,
2857 SMESHDS_Mesh * aMesh)
2859 const set<SMESHDS_GroupBase*>& groups = aMesh->GetGroups();
2860 if (!groups.empty())
2862 set<SMESHDS_GroupBase*>::const_iterator GrIt = groups.begin();
2863 for (; GrIt != groups.end(); GrIt++)
2865 SMESHDS_Group* grp = dynamic_cast<SMESHDS_Group*>(*GrIt);
2866 if (!grp || grp->IsEmpty()) continue;
2867 grp->SMDSGroup().Remove(removeelem);
2872 //================================================================================
2874 * \brief Replace elemToRm by elemToAdd in the all groups
2876 //================================================================================
2878 void SMESH_MeshEditor::ReplaceElemInGroups (const SMDS_MeshElement* elemToRm,
2879 const SMDS_MeshElement* elemToAdd,
2880 SMESHDS_Mesh * aMesh)
2882 const set<SMESHDS_GroupBase*>& groups = aMesh->GetGroups();
2883 if (!groups.empty()) {
2884 set<SMESHDS_GroupBase*>::const_iterator grIt = groups.begin();
2885 for ( ; grIt != groups.end(); grIt++ ) {
2886 SMESHDS_Group* group = dynamic_cast<SMESHDS_Group*>( *grIt );
2887 if ( group && group->SMDSGroup().Remove( elemToRm ) && elemToAdd )
2888 group->SMDSGroup().Add( elemToAdd );
2893 //================================================================================
2895 * \brief Replace elemToRm by elemToAdd in the all groups
2897 //================================================================================
2899 void SMESH_MeshEditor::ReplaceElemInGroups (const SMDS_MeshElement* elemToRm,
2900 const vector<const SMDS_MeshElement*>& elemToAdd,
2901 SMESHDS_Mesh * aMesh)
2903 const set<SMESHDS_GroupBase*>& groups = aMesh->GetGroups();
2904 if (!groups.empty())
2906 set<SMESHDS_GroupBase*>::const_iterator grIt = groups.begin();
2907 for ( ; grIt != groups.end(); grIt++ ) {
2908 SMESHDS_Group* group = dynamic_cast<SMESHDS_Group*>( *grIt );
2909 if ( group && group->SMDSGroup().Remove( elemToRm ) )
2910 for ( size_t i = 0; i < elemToAdd.size(); ++i )
2911 group->SMDSGroup().Add( elemToAdd[ i ] );
2916 //=======================================================================
2917 //function : QuadToTri
2918 //purpose : Cut quadrangles into triangles.
2919 // theCrit is used to select a diagonal to cut
2920 //=======================================================================
2922 bool SMESH_MeshEditor::QuadToTri (TIDSortedElemSet & theElems,
2923 const bool the13Diag)
2925 myLastCreatedElems.Clear();
2926 myLastCreatedNodes.Clear();
2928 SMESHDS_Mesh * aMesh = GetMeshDS();
2930 Handle(Geom_Surface) surface;
2931 SMESH_MesherHelper helper( *GetMesh() );
2933 TIDSortedElemSet::iterator itElem;
2934 for ( itElem = theElems.begin(); itElem != theElems.end(); itElem++ )
2936 const SMDS_MeshElement* elem = *itElem;
2937 if ( !elem || elem->GetGeomType() != SMDSGeom_QUADRANGLE )
2940 if ( elem->NbNodes() == 4 ) {
2941 // retrieve element nodes
2942 const SMDS_MeshNode* aNodes [4];
2943 SMDS_ElemIteratorPtr itN = elem->nodesIterator();
2945 while ( itN->more() )
2946 aNodes[ i++ ] = static_cast<const SMDS_MeshNode*>( itN->next() );
2948 int aShapeId = FindShape( elem );
2949 const SMDS_MeshElement* newElem1 = 0;
2950 const SMDS_MeshElement* newElem2 = 0;
2952 newElem1 = aMesh->AddFace( aNodes[2], aNodes[0], aNodes[1] );
2953 newElem2 = aMesh->AddFace( aNodes[2], aNodes[3], aNodes[0] );
2956 newElem1 = aMesh->AddFace( aNodes[3], aNodes[0], aNodes[1] );
2957 newElem2 = aMesh->AddFace( aNodes[3], aNodes[1], aNodes[2] );
2959 myLastCreatedElems.Append(newElem1);
2960 myLastCreatedElems.Append(newElem2);
2961 // put a new triangle on the same shape and add to the same groups
2964 aMesh->SetMeshElementOnShape( newElem1, aShapeId );
2965 aMesh->SetMeshElementOnShape( newElem2, aShapeId );
2967 AddToSameGroups( newElem1, elem, aMesh );
2968 AddToSameGroups( newElem2, elem, aMesh );
2969 aMesh->RemoveElement( elem );
2972 // Quadratic quadrangle
2974 else if ( elem->NbNodes() >= 8 )
2976 // get surface elem is on
2977 int aShapeId = FindShape( elem );
2978 if ( aShapeId != helper.GetSubShapeID() ) {
2982 shape = aMesh->IndexToShape( aShapeId );
2983 if ( !shape.IsNull() && shape.ShapeType() == TopAbs_FACE ) {
2984 TopoDS_Face face = TopoDS::Face( shape );
2985 surface = BRep_Tool::Surface( face );
2986 if ( !surface.IsNull() )