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_Downward.hxx"
30 #include "SMDS_EdgePosition.hxx"
31 #include "SMDS_FaceOfNodes.hxx"
32 #include "SMDS_FacePosition.hxx"
33 #include "SMDS_LinearEdge.hxx"
34 #include "SMDS_MeshGroup.hxx"
35 #include "SMDS_SetIterator.hxx"
36 #include "SMDS_SpacePosition.hxx"
37 #include "SMDS_VolumeTool.hxx"
38 #include "SMESHDS_Group.hxx"
39 #include "SMESHDS_Mesh.hxx"
40 #include "SMESH_Algo.hxx"
41 #include "SMESH_ControlsDef.hxx"
42 #include "SMESH_Group.hxx"
43 #include "SMESH_Mesh.hxx"
44 #include "SMESH_MeshAlgos.hxx"
45 #include "SMESH_MesherHelper.hxx"
46 #include "SMESH_OctreeNode.hxx"
47 #include "SMESH_subMesh.hxx"
49 #include <Basics_OCCTVersion.hxx>
51 #include "utilities.h"
54 #include <BRepAdaptor_Surface.hxx>
55 #include <BRepBuilderAPI_MakeEdge.hxx>
56 #include <BRepClass3d_SolidClassifier.hxx>
57 #include <BRep_Tool.hxx>
59 #include <Extrema_GenExtPS.hxx>
60 #include <Extrema_POnCurv.hxx>
61 #include <Extrema_POnSurf.hxx>
62 #include <Geom2d_Curve.hxx>
63 #include <GeomAdaptor_Surface.hxx>
64 #include <Geom_Curve.hxx>
65 #include <Geom_Surface.hxx>
66 #include <Precision.hxx>
67 #include <TColStd_ListOfInteger.hxx>
68 #include <TopAbs_State.hxx>
70 #include <TopExp_Explorer.hxx>
71 #include <TopTools_ListIteratorOfListOfShape.hxx>
72 #include <TopTools_ListOfShape.hxx>
73 #include <TopTools_SequenceOfShape.hxx>
75 #include <TopoDS_Edge.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 Return mesh DS
132 //================================================================================
134 SMESHDS_Mesh * SMESH_MeshEditor::GetMeshDS()
136 return myMesh->GetMeshDS();
140 //================================================================================
142 * \brief Clears myLastCreatedNodes and myLastCreatedElems
144 //================================================================================
146 void SMESH_MeshEditor::ClearLastCreated()
148 myLastCreatedNodes.Clear();
149 myLastCreatedElems.Clear();
152 //================================================================================
154 * \brief Initializes members by an existing element
155 * \param [in] elem - the source element
156 * \param [in] basicOnly - if true, does not set additional data of Ball and Polyhedron
158 //================================================================================
160 SMESH_MeshEditor::ElemFeatures&
161 SMESH_MeshEditor::ElemFeatures::Init( const SMDS_MeshElement* elem, bool basicOnly )
165 myType = elem->GetType();
166 if ( myType == SMDSAbs_Face || myType == SMDSAbs_Volume )
168 myIsPoly = elem->IsPoly();
171 myIsQuad = elem->IsQuadratic();
172 if ( myType == SMDSAbs_Volume && !basicOnly )
174 vector<int > quant = static_cast<const SMDS_VtkVolume* >( elem )->GetQuantities();
175 myPolyhedQuantities.swap( quant );
179 else if ( myType == SMDSAbs_Ball && !basicOnly )
181 myBallDiameter = static_cast<const SMDS_BallElement*>(elem)->GetDiameter();
187 //=======================================================================
191 //=======================================================================
194 SMESH_MeshEditor::AddElement(const vector<const SMDS_MeshNode*> & node,
195 const ElemFeatures& features)
197 SMDS_MeshElement* e = 0;
198 int nbnode = node.size();
199 SMESHDS_Mesh* mesh = GetMeshDS();
200 const int ID = features.myID;
202 switch ( features.myType ) {
204 if ( !features.myIsPoly ) {
206 if ( ID >= 1 ) e = mesh->AddFaceWithID(node[0], node[1], node[2], ID);
207 else e = mesh->AddFace (node[0], node[1], node[2] );
209 else if (nbnode == 4) {
210 if ( ID >= 1 ) e = mesh->AddFaceWithID(node[0], node[1], node[2], node[3], ID);
211 else e = mesh->AddFace (node[0], node[1], node[2], node[3] );
213 else if (nbnode == 6) {
214 if ( ID >= 1 ) e = mesh->AddFaceWithID(node[0], node[1], node[2], node[3],
215 node[4], node[5], ID);
216 else e = mesh->AddFace (node[0], node[1], node[2], node[3],
219 else if (nbnode == 7) {
220 if ( ID >= 1 ) e = mesh->AddFaceWithID(node[0], node[1], node[2], node[3],
221 node[4], node[5], node[6], ID);
222 else e = mesh->AddFace (node[0], node[1], node[2], node[3],
223 node[4], node[5], node[6] );
225 else if (nbnode == 8) {
226 if ( ID >= 1 ) e = mesh->AddFaceWithID(node[0], node[1], node[2], node[3],
227 node[4], node[5], node[6], node[7], ID);
228 else e = mesh->AddFace (node[0], node[1], node[2], node[3],
229 node[4], node[5], node[6], node[7] );
231 else if (nbnode == 9) {
232 if ( ID >= 1 ) e = mesh->AddFaceWithID(node[0], node[1], node[2], node[3],
233 node[4], node[5], node[6], node[7], node[8], ID);
234 else e = mesh->AddFace (node[0], node[1], node[2], node[3],
235 node[4], node[5], node[6], node[7], node[8] );
238 else if ( !features.myIsQuad )
240 if ( ID >= 1 ) e = mesh->AddPolygonalFaceWithID(node, ID);
241 else e = mesh->AddPolygonalFace (node );
243 else if ( nbnode % 2 == 0 ) // just a protection
245 if ( ID >= 1 ) e = mesh->AddQuadPolygonalFaceWithID(node, ID);
246 else e = mesh->AddQuadPolygonalFace (node );
251 if ( !features.myIsPoly ) {
253 if ( ID >= 1 ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3], ID);
254 else e = mesh->AddVolume (node[0], node[1], node[2], node[3] );
256 else if (nbnode == 5) {
257 if ( ID >= 1 ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3],
259 else e = mesh->AddVolume (node[0], node[1], node[2], node[3],
262 else if (nbnode == 6) {
263 if ( ID >= 1 ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3],
264 node[4], node[5], ID);
265 else e = mesh->AddVolume (node[0], node[1], node[2], node[3],
268 else if (nbnode == 8) {
269 if ( ID >= 1 ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3],
270 node[4], node[5], node[6], node[7], ID);
271 else e = mesh->AddVolume (node[0], node[1], node[2], node[3],
272 node[4], node[5], node[6], node[7] );
274 else if (nbnode == 10) {
275 if ( ID >= 1 ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3],
276 node[4], node[5], node[6], node[7],
277 node[8], node[9], ID);
278 else e = mesh->AddVolume (node[0], node[1], node[2], node[3],
279 node[4], node[5], node[6], node[7],
282 else if (nbnode == 12) {
283 if ( ID >= 1 ) e = mesh->AddVolumeWithID(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], ID);
286 else e = mesh->AddVolume (node[0], node[1], node[2], node[3],
287 node[4], node[5], node[6], node[7],
288 node[8], node[9], node[10], node[11] );
290 else if (nbnode == 13) {
291 if ( ID >= 1 ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3],
292 node[4], node[5], node[6], node[7],
293 node[8], node[9], node[10],node[11],
295 else e = mesh->AddVolume (node[0], node[1], node[2], node[3],
296 node[4], node[5], node[6], node[7],
297 node[8], node[9], node[10],node[11],
300 else if (nbnode == 15) {
301 if ( ID >= 1 ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3],
302 node[4], node[5], node[6], node[7],
303 node[8], node[9], node[10],node[11],
304 node[12],node[13],node[14],ID);
305 else e = mesh->AddVolume (node[0], node[1], node[2], node[3],
306 node[4], node[5], node[6], node[7],
307 node[8], node[9], node[10],node[11],
308 node[12],node[13],node[14] );
310 else if (nbnode == 20) {
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],ID);
316 else e = mesh->AddVolume (node[0], node[1], node[2], node[3],
317 node[4], node[5], node[6], node[7],
318 node[8], node[9], node[10],node[11],
319 node[12],node[13],node[14],node[15],
320 node[16],node[17],node[18],node[19] );
322 else if (nbnode == 27) {
323 if ( ID >= 1 ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3],
324 node[4], node[5], node[6], node[7],
325 node[8], node[9], node[10],node[11],
326 node[12],node[13],node[14],node[15],
327 node[16],node[17],node[18],node[19],
328 node[20],node[21],node[22],node[23],
329 node[24],node[25],node[26], ID);
330 else e = mesh->AddVolume (node[0], node[1], node[2], node[3],
331 node[4], node[5], node[6], node[7],
332 node[8], node[9], node[10],node[11],
333 node[12],node[13],node[14],node[15],
334 node[16],node[17],node[18],node[19],
335 node[20],node[21],node[22],node[23],
336 node[24],node[25],node[26] );
339 else if ( !features.myIsQuad )
341 if ( ID >= 1 ) e = mesh->AddPolyhedralVolumeWithID(node, features.myPolyhedQuantities, ID);
342 else e = mesh->AddPolyhedralVolume (node, features.myPolyhedQuantities );
346 // if ( ID >= 1 ) e = mesh->AddQuadPolyhedralVolumeWithID(node, features.myPolyhedQuantities,ID);
347 // else e = mesh->AddQuadPolyhedralVolume (node, features.myPolyhedQuantities );
353 if ( ID >= 1 ) e = mesh->AddEdgeWithID(node[0], node[1], ID);
354 else e = mesh->AddEdge (node[0], node[1] );
356 else if ( nbnode == 3 ) {
357 if ( ID >= 1 ) e = mesh->AddEdgeWithID(node[0], node[1], node[2], ID);
358 else e = mesh->AddEdge (node[0], node[1], node[2] );
362 case SMDSAbs_0DElement:
364 if ( ID >= 1 ) e = mesh->Add0DElementWithID(node[0], ID);
365 else e = mesh->Add0DElement (node[0] );
370 if ( ID >= 1 ) e = mesh->AddNodeWithID(node[0]->X(), node[0]->Y(), node[0]->Z(), ID);
371 else e = mesh->AddNode (node[0]->X(), node[0]->Y(), node[0]->Z() );
375 if ( ID >= 1 ) e = mesh->AddBallWithID(node[0], features.myBallDiameter, ID);
376 else e = mesh->AddBall (node[0], features.myBallDiameter );
381 if ( e ) myLastCreatedElems.Append( e );
385 //=======================================================================
389 //=======================================================================
391 SMDS_MeshElement* SMESH_MeshEditor::AddElement(const vector<int> & nodeIDs,
392 const ElemFeatures& features)
394 vector<const SMDS_MeshNode*> nodes;
395 nodes.reserve( nodeIDs.size() );
396 vector<int>::const_iterator id = nodeIDs.begin();
397 while ( id != nodeIDs.end() ) {
398 if ( const SMDS_MeshNode* node = GetMeshDS()->FindNode( *id++ ))
399 nodes.push_back( node );
403 return AddElement( nodes, features );
406 //=======================================================================
408 //purpose : Remove a node or an element.
409 // Modify a compute state of sub-meshes which become empty
410 //=======================================================================
412 int SMESH_MeshEditor::Remove (const list< int >& theIDs,
415 myLastCreatedElems.Clear();
416 myLastCreatedNodes.Clear();
418 SMESHDS_Mesh* aMesh = GetMeshDS();
419 set< SMESH_subMesh *> smmap;
422 list<int>::const_iterator it = theIDs.begin();
423 for ( ; it != theIDs.end(); it++ ) {
424 const SMDS_MeshElement * elem;
426 elem = aMesh->FindNode( *it );
428 elem = aMesh->FindElement( *it );
432 // Notify VERTEX sub-meshes about modification
434 const SMDS_MeshNode* node = cast2Node( elem );
435 if ( node->GetPosition()->GetTypeOfPosition() == SMDS_TOP_VERTEX )
436 if ( int aShapeID = node->getshapeId() )
437 if ( SMESH_subMesh * sm = GetMesh()->GetSubMeshContaining( aShapeID ) )
440 // Find sub-meshes to notify about modification
441 // SMDS_ElemIteratorPtr nodeIt = elem->nodesIterator();
442 // while ( nodeIt->more() ) {
443 // const SMDS_MeshNode* node = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
444 // const SMDS_PositionPtr& aPosition = node->GetPosition();
445 // if ( aPosition.get() ) {
446 // if ( int aShapeID = aPosition->GetShapeId() ) {
447 // if ( SMESH_subMesh * sm = GetMesh()->GetSubMeshContaining( aShapeID ) )
448 // smmap.insert( sm );
455 aMesh->RemoveNode( static_cast< const SMDS_MeshNode* >( elem ));
457 aMesh->RemoveElement( elem );
461 // Notify sub-meshes about modification
462 if ( !smmap.empty() ) {
463 set< SMESH_subMesh *>::iterator smIt;
464 for ( smIt = smmap.begin(); smIt != smmap.end(); smIt++ )
465 (*smIt)->ComputeStateEngine( SMESH_subMesh::MESH_ENTITY_REMOVED );
468 // // Check if the whole mesh becomes empty
469 // if ( SMESH_subMesh * sm = GetMesh()->GetSubMeshContaining( 1 ) )
470 // sm->ComputeStateEngine( SMESH_subMesh::CHECK_COMPUTE_STATE );
475 //================================================================================
477 * \brief Create 0D elements on all nodes of the given object except those
478 * nodes on which a 0D element already exists.
479 * \param elements - Elements on whose nodes to create 0D elements; if empty,
480 * the all mesh is treated
481 * \param all0DElems - returns all 0D elements found or created on nodes of \a elements
483 //================================================================================
485 void SMESH_MeshEditor::Create0DElementsOnAllNodes( const TIDSortedElemSet& elements,
486 TIDSortedElemSet& all0DElems )
488 SMDS_ElemIteratorPtr elemIt;
489 vector< const SMDS_MeshElement* > allNodes;
490 if ( elements.empty() )
492 allNodes.reserve( GetMeshDS()->NbNodes() );
493 elemIt = GetMeshDS()->elementsIterator( SMDSAbs_Node );
494 while ( elemIt->more() )
495 allNodes.push_back( elemIt->next() );
497 elemIt = elemSetIterator( allNodes );
501 elemIt = elemSetIterator( elements );
504 while ( elemIt->more() )
506 const SMDS_MeshElement* e = elemIt->next();
507 SMDS_ElemIteratorPtr nodeIt = e->nodesIterator();
508 while ( nodeIt->more() )
510 const SMDS_MeshNode* n = cast2Node( nodeIt->next() );
511 SMDS_ElemIteratorPtr it0D = n->GetInverseElementIterator( SMDSAbs_0DElement );
513 all0DElems.insert( it0D->next() );
515 myLastCreatedElems.Append( GetMeshDS()->Add0DElement( n ));
516 all0DElems.insert( myLastCreatedElems.Last() );
522 //=======================================================================
523 //function : FindShape
524 //purpose : Return an index of the shape theElem is on
525 // or zero if a shape not found
526 //=======================================================================
528 int SMESH_MeshEditor::FindShape (const SMDS_MeshElement * theElem)
530 myLastCreatedElems.Clear();
531 myLastCreatedNodes.Clear();
533 SMESHDS_Mesh * aMesh = GetMeshDS();
534 if ( aMesh->ShapeToMesh().IsNull() )
537 int aShapeID = theElem->getshapeId();
541 if ( SMESHDS_SubMesh * sm = aMesh->MeshElements( aShapeID ))
542 if ( sm->Contains( theElem ))
545 if ( theElem->GetType() == SMDSAbs_Node ) {
546 MESSAGE( ":( Error: invalid myShapeId of node " << theElem->GetID() );
549 MESSAGE( ":( Error: invalid myShapeId of element " << theElem->GetID() );
552 TopoDS_Shape aShape; // the shape a node of theElem is on
553 if ( theElem->GetType() != SMDSAbs_Node )
555 SMDS_ElemIteratorPtr nodeIt = theElem->nodesIterator();
556 while ( nodeIt->more() ) {
557 const SMDS_MeshNode* node = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
558 if ((aShapeID = node->getshapeId()) > 0) {
559 if ( SMESHDS_SubMesh * sm = aMesh->MeshElements( aShapeID ) ) {
560 if ( sm->Contains( theElem ))
562 if ( aShape.IsNull() )
563 aShape = aMesh->IndexToShape( aShapeID );
569 // None of nodes is on a proper shape,
570 // find the shape among ancestors of aShape on which a node is
571 if ( !aShape.IsNull() ) {
572 TopTools_ListIteratorOfListOfShape ancIt( GetMesh()->GetAncestors( aShape ));
573 for ( ; ancIt.More(); ancIt.Next() ) {
574 SMESHDS_SubMesh * sm = aMesh->MeshElements( ancIt.Value() );
575 if ( sm && sm->Contains( theElem ))
576 return aMesh->ShapeToIndex( ancIt.Value() );
581 SMESHDS_SubMeshIteratorPtr smIt = GetMeshDS()->SubMeshes();
582 while ( const SMESHDS_SubMesh* sm = smIt->next() )
583 if ( sm->Contains( theElem ))
590 //=======================================================================
591 //function : IsMedium
593 //=======================================================================
595 bool SMESH_MeshEditor::IsMedium(const SMDS_MeshNode* node,
596 const SMDSAbs_ElementType typeToCheck)
598 bool isMedium = false;
599 SMDS_ElemIteratorPtr it = node->GetInverseElementIterator(typeToCheck);
600 while (it->more() && !isMedium ) {
601 const SMDS_MeshElement* elem = it->next();
602 isMedium = elem->IsMediumNode(node);
607 //=======================================================================
608 //function : shiftNodesQuadTria
609 //purpose : Shift nodes in the array corresponded to quadratic triangle
610 // example: (0,1,2,3,4,5) -> (1,2,0,4,5,3)
611 //=======================================================================
613 static void shiftNodesQuadTria(vector< const SMDS_MeshNode* >& aNodes)
615 const SMDS_MeshNode* nd1 = aNodes[0];
616 aNodes[0] = aNodes[1];
617 aNodes[1] = aNodes[2];
619 const SMDS_MeshNode* nd2 = aNodes[3];
620 aNodes[3] = aNodes[4];
621 aNodes[4] = aNodes[5];
625 //=======================================================================
626 //function : nbEdgeConnectivity
627 //purpose : return number of the edges connected with the theNode.
628 // if theEdges has connections with the other type of the
629 // elements, return -1
630 //=======================================================================
632 static int nbEdgeConnectivity(const SMDS_MeshNode* theNode)
634 // SMDS_ElemIteratorPtr elemIt = theNode->GetInverseElementIterator();
636 // while(elemIt->more()) {
641 return theNode->NbInverseElements();
644 //=======================================================================
645 //function : getNodesFromTwoTria
647 //=======================================================================
649 static bool getNodesFromTwoTria(const SMDS_MeshElement * theTria1,
650 const SMDS_MeshElement * theTria2,
651 vector< const SMDS_MeshNode*>& N1,
652 vector< const SMDS_MeshNode*>& N2)
654 N1.assign( theTria1->begin_nodes(), theTria1->end_nodes() );
655 if ( N1.size() < 6 ) return false;
656 N2.assign( theTria2->begin_nodes(), theTria2->end_nodes() );
657 if ( N2.size() < 6 ) return false;
659 int sames[3] = {-1,-1,-1};
671 if(nbsames!=2) return false;
673 shiftNodesQuadTria(N1);
675 shiftNodesQuadTria(N1);
678 i = sames[0] + sames[1] + sames[2];
680 shiftNodesQuadTria(N2);
682 // now we receive following N1 and N2 (using numeration as in the image below)
683 // tria1 : (1 2 4 5 9 7) and tria2 : (3 4 2 8 9 6)
684 // i.e. first nodes from both arrays form a new diagonal
688 //=======================================================================
689 //function : InverseDiag
690 //purpose : Replace two neighbour triangles with ones built on the same 4 nodes
691 // but having other common link.
692 // Return False if args are improper
693 //=======================================================================
695 bool SMESH_MeshEditor::InverseDiag (const SMDS_MeshElement * theTria1,
696 const SMDS_MeshElement * theTria2 )
698 myLastCreatedElems.Clear();
699 myLastCreatedNodes.Clear();
701 if (!theTria1 || !theTria2)
704 const SMDS_VtkFace* F1 = dynamic_cast<const SMDS_VtkFace*>( theTria1 );
705 if (!F1) return false;
706 const SMDS_VtkFace* F2 = dynamic_cast<const SMDS_VtkFace*>( theTria2 );
707 if (!F2) return false;
708 if ((theTria1->GetEntityType() == SMDSEntity_Triangle) &&
709 (theTria2->GetEntityType() == SMDSEntity_Triangle)) {
711 // 1 +--+ A theTria1: ( 1 A B ) A->2 ( 1 2 B ) 1 +--+ A
712 // | /| theTria2: ( B A 2 ) B->1 ( 1 A 2 ) |\ |
716 // put nodes in array and find out indices of the same ones
717 const SMDS_MeshNode* aNodes [6];
718 int sameInd [] = { -1, -1, -1, -1, -1, -1 };
720 SMDS_ElemIteratorPtr it = theTria1->nodesIterator();
721 while ( it->more() ) {
722 aNodes[ i ] = static_cast<const SMDS_MeshNode*>( it->next() );
724 if ( i > 2 ) // theTria2
725 // find same node of theTria1
726 for ( int j = 0; j < 3; j++ )
727 if ( aNodes[ i ] == aNodes[ j ]) {
736 return false; // theTria1 is not a triangle
737 it = theTria2->nodesIterator();
739 if ( i == 6 && it->more() )
740 return false; // theTria2 is not a triangle
743 // find indices of 1,2 and of A,B in theTria1
744 int iA = -1, iB = 0, i1 = 0, i2 = 0;
745 for ( i = 0; i < 6; i++ ) {
746 if ( sameInd [ i ] == -1 ) {
751 if ( iA >= 0) iB = i;
755 // nodes 1 and 2 should not be the same
756 if ( aNodes[ i1 ] == aNodes[ i2 ] )
760 aNodes[ iA ] = aNodes[ i2 ];
762 aNodes[ sameInd[ iB ]] = aNodes[ i1 ];
764 GetMeshDS()->ChangeElementNodes( theTria1, aNodes, 3 );
765 GetMeshDS()->ChangeElementNodes( theTria2, &aNodes[ 3 ], 3 );
769 } // end if(F1 && F2)
771 // check case of quadratic faces
772 if (theTria1->GetEntityType() != SMDSEntity_Quad_Triangle &&
773 theTria1->GetEntityType() != SMDSEntity_BiQuad_Triangle)
775 if (theTria2->GetEntityType() != SMDSEntity_Quad_Triangle&&
776 theTria2->GetEntityType() != SMDSEntity_BiQuad_Triangle)
780 // 1 +--+--+ 2 theTria1: (1 2 4 5 9 7) or (2 4 1 9 7 5) or (4 1 2 7 5 9)
781 // | /| theTria2: (2 3 4 6 8 9) or (3 4 2 8 9 6) or (4 2 3 9 6 8)
789 vector< const SMDS_MeshNode* > N1;
790 vector< const SMDS_MeshNode* > N2;
791 if(!getNodesFromTwoTria(theTria1,theTria2,N1,N2))
793 // now we receive following N1 and N2 (using numeration as above image)
794 // tria1 : (1 2 4 5 9 7) and tria2 : (3 4 2 8 9 6)
795 // i.e. first nodes from both arrays determ new diagonal
797 vector< const SMDS_MeshNode*> N1new( N1.size() );
798 vector< const SMDS_MeshNode*> N2new( N2.size() );
799 N1new.back() = N1.back(); // central node of biquadratic
800 N2new.back() = N2.back();
801 N1new[0] = N1[0]; N2new[0] = N1[0];
802 N1new[1] = N2[0]; N2new[1] = N1[1];
803 N1new[2] = N2[1]; N2new[2] = N2[0];
804 N1new[3] = N1[4]; N2new[3] = N1[3];
805 N1new[4] = N2[3]; N2new[4] = N2[5];
806 N1new[5] = N1[5]; N2new[5] = N1[4];
807 // change nodes in faces
808 GetMeshDS()->ChangeElementNodes( theTria1, &N1new[0], N1new.size() );
809 GetMeshDS()->ChangeElementNodes( theTria2, &N2new[0], N2new.size() );
811 // move the central node of biquadratic triangle
812 SMESH_MesherHelper helper( *GetMesh() );
813 for ( int is2nd = 0; is2nd < 2; ++is2nd )
815 const SMDS_MeshElement* tria = is2nd ? theTria2 : theTria1;
816 vector< const SMDS_MeshNode*>& nodes = is2nd ? N2new : N1new;
817 if ( nodes.size() < 7 )
819 helper.SetSubShape( tria->getshapeId() );
820 const TopoDS_Face& F = TopoDS::Face( helper.GetSubShape() );
824 xyz = ( SMESH_TNodeXYZ( nodes[3] ) +
825 SMESH_TNodeXYZ( nodes[4] ) +
826 SMESH_TNodeXYZ( nodes[5] )) / 3.;
831 gp_XY uv = ( helper.GetNodeUV( F, nodes[3], nodes[2], &checkUV ) +
832 helper.GetNodeUV( F, nodes[4], nodes[0], &checkUV ) +
833 helper.GetNodeUV( F, nodes[5], nodes[1], &checkUV )) / 3.;
835 Handle(Geom_Surface) S = BRep_Tool::Surface(F,loc);
836 xyz = S->Value( uv.X(), uv.Y() );
837 xyz.Transform( loc );
838 if ( nodes[6]->GetPosition()->GetTypeOfPosition() == SMDS_TOP_FACE && // set UV
839 nodes[6]->getshapeId() > 0 )
840 GetMeshDS()->SetNodeOnFace( nodes[6], nodes[6]->getshapeId(), uv.X(), uv.Y() );
842 GetMeshDS()->MoveNode( nodes[6], xyz.X(), xyz.Y(), xyz.Z() );
847 //=======================================================================
848 //function : findTriangles
849 //purpose : find triangles sharing theNode1-theNode2 link
850 //=======================================================================
852 static bool findTriangles(const SMDS_MeshNode * theNode1,
853 const SMDS_MeshNode * theNode2,
854 const SMDS_MeshElement*& theTria1,
855 const SMDS_MeshElement*& theTria2)
857 if ( !theNode1 || !theNode2 ) return false;
859 theTria1 = theTria2 = 0;
861 set< const SMDS_MeshElement* > emap;
862 SMDS_ElemIteratorPtr it = theNode1->GetInverseElementIterator(SMDSAbs_Face);
864 const SMDS_MeshElement* elem = it->next();
865 if ( elem->NbCornerNodes() == 3 )
868 it = theNode2->GetInverseElementIterator(SMDSAbs_Face);
870 const SMDS_MeshElement* elem = it->next();
871 if ( emap.count( elem )) {
879 // theTria1 must be element with minimum ID
880 if ( theTria2->GetID() < theTria1->GetID() )
881 std::swap( theTria2, theTria1 );
889 //=======================================================================
890 //function : InverseDiag
891 //purpose : Replace two neighbour triangles sharing theNode1-theNode2 link
892 // with ones built on the same 4 nodes but having other common link.
893 // Return false if proper faces not found
894 //=======================================================================
896 bool SMESH_MeshEditor::InverseDiag (const SMDS_MeshNode * theNode1,
897 const SMDS_MeshNode * theNode2)
899 myLastCreatedElems.Clear();
900 myLastCreatedNodes.Clear();
902 const SMDS_MeshElement *tr1, *tr2;
903 if ( !findTriangles( theNode1, theNode2, tr1, tr2 ))
906 const SMDS_VtkFace* F1 = dynamic_cast<const SMDS_VtkFace*>( tr1 );
907 if (!F1) return false;
908 const SMDS_VtkFace* F2 = dynamic_cast<const SMDS_VtkFace*>( tr2 );
909 if (!F2) return false;
910 if ((tr1->GetEntityType() == SMDSEntity_Triangle) &&
911 (tr2->GetEntityType() == SMDSEntity_Triangle)) {
913 // 1 +--+ A tr1: ( 1 A B ) A->2 ( 1 2 B ) 1 +--+ A
914 // | /| tr2: ( B A 2 ) B->1 ( 1 A 2 ) |\ |
918 // put nodes in array
919 // and find indices of 1,2 and of A in tr1 and of B in tr2
920 int i, iA1 = 0, i1 = 0;
921 const SMDS_MeshNode* aNodes1 [3];
922 SMDS_ElemIteratorPtr it;
923 for (i = 0, it = tr1->nodesIterator(); it->more(); i++ ) {
924 aNodes1[ i ] = static_cast<const SMDS_MeshNode*>( it->next() );
925 if ( aNodes1[ i ] == theNode1 )
926 iA1 = i; // node A in tr1
927 else if ( aNodes1[ i ] != theNode2 )
931 const SMDS_MeshNode* aNodes2 [3];
932 for (i = 0, it = tr2->nodesIterator(); it->more(); i++ ) {
933 aNodes2[ i ] = static_cast<const SMDS_MeshNode*>( it->next() );
934 if ( aNodes2[ i ] == theNode2 )
935 iB2 = i; // node B in tr2
936 else if ( aNodes2[ i ] != theNode1 )
940 // nodes 1 and 2 should not be the same
941 if ( aNodes1[ i1 ] == aNodes2[ i2 ] )
945 aNodes1[ iA1 ] = aNodes2[ i2 ];
947 aNodes2[ iB2 ] = aNodes1[ i1 ];
949 GetMeshDS()->ChangeElementNodes( tr1, aNodes1, 3 );
950 GetMeshDS()->ChangeElementNodes( tr2, aNodes2, 3 );
955 // check case of quadratic faces
956 return InverseDiag(tr1,tr2);
959 //=======================================================================
960 //function : getQuadrangleNodes
961 //purpose : fill theQuadNodes - nodes of a quadrangle resulting from
962 // fusion of triangles tr1 and tr2 having shared link on
963 // theNode1 and theNode2
964 //=======================================================================
966 bool getQuadrangleNodes(const SMDS_MeshNode * theQuadNodes [],
967 const SMDS_MeshNode * theNode1,
968 const SMDS_MeshNode * theNode2,
969 const SMDS_MeshElement * tr1,
970 const SMDS_MeshElement * tr2 )
972 if( tr1->NbNodes() != tr2->NbNodes() )
974 // find the 4-th node to insert into tr1
975 const SMDS_MeshNode* n4 = 0;
976 SMDS_ElemIteratorPtr it = tr2->nodesIterator();
978 while ( !n4 && i<3 ) {
979 const SMDS_MeshNode * n = cast2Node( it->next() );
981 bool isDiag = ( n == theNode1 || n == theNode2 );
985 // Make an array of nodes to be in a quadrangle
986 int iNode = 0, iFirstDiag = -1;
987 it = tr1->nodesIterator();
990 const SMDS_MeshNode * n = cast2Node( it->next() );
992 bool isDiag = ( n == theNode1 || n == theNode2 );
994 if ( iFirstDiag < 0 )
996 else if ( iNode - iFirstDiag == 1 )
997 theQuadNodes[ iNode++ ] = n4; // insert the 4-th node between diagonal nodes
999 else if ( n == n4 ) {
1000 return false; // tr1 and tr2 should not have all the same nodes
1002 theQuadNodes[ iNode++ ] = n;
1004 if ( iNode == 3 ) // diagonal nodes have 0 and 2 indices
1005 theQuadNodes[ iNode ] = n4;
1010 //=======================================================================
1011 //function : DeleteDiag
1012 //purpose : Replace two neighbour triangles sharing theNode1-theNode2 link
1013 // with a quadrangle built on the same 4 nodes.
1014 // Return false if proper faces not found
1015 //=======================================================================
1017 bool SMESH_MeshEditor::DeleteDiag (const SMDS_MeshNode * theNode1,
1018 const SMDS_MeshNode * theNode2)
1020 myLastCreatedElems.Clear();
1021 myLastCreatedNodes.Clear();
1023 const SMDS_MeshElement *tr1, *tr2;
1024 if ( !findTriangles( theNode1, theNode2, tr1, tr2 ))
1027 const SMDS_VtkFace* F1 = dynamic_cast<const SMDS_VtkFace*>( tr1 );
1028 if (!F1) return false;
1029 const SMDS_VtkFace* F2 = dynamic_cast<const SMDS_VtkFace*>( tr2 );
1030 if (!F2) return false;
1031 SMESHDS_Mesh * aMesh = GetMeshDS();
1033 if ((tr1->GetEntityType() == SMDSEntity_Triangle) &&
1034 (tr2->GetEntityType() == SMDSEntity_Triangle)) {
1036 const SMDS_MeshNode* aNodes [ 4 ];
1037 if ( ! getQuadrangleNodes( aNodes, theNode1, theNode2, tr1, tr2 ))
1040 const SMDS_MeshElement* newElem = 0;
1041 newElem = aMesh->AddFace( aNodes[0], aNodes[1], aNodes[2], aNodes[3] );
1042 myLastCreatedElems.Append(newElem);
1043 AddToSameGroups( newElem, tr1, aMesh );
1044 int aShapeId = tr1->getshapeId();
1047 aMesh->SetMeshElementOnShape( newElem, aShapeId );
1049 aMesh->RemoveElement( tr1 );
1050 aMesh->RemoveElement( tr2 );
1055 // check case of quadratic faces
1056 if (tr1->GetEntityType() != SMDSEntity_Quad_Triangle)
1058 if (tr2->GetEntityType() != SMDSEntity_Quad_Triangle)
1062 // 1 +--+--+ 2 tr1: (1 2 4 5 9 7) or (2 4 1 9 7 5) or (4 1 2 7 5 9)
1063 // | /| tr2: (2 3 4 6 8 9) or (3 4 2 8 9 6) or (4 2 3 9 6 8)
1071 vector< const SMDS_MeshNode* > N1;
1072 vector< const SMDS_MeshNode* > N2;
1073 if(!getNodesFromTwoTria(tr1,tr2,N1,N2))
1075 // now we receive following N1 and N2 (using numeration as above image)
1076 // tria1 : (1 2 4 5 9 7) and tria2 : (3 4 2 8 9 6)
1077 // i.e. first nodes from both arrays determ new diagonal
1079 const SMDS_MeshNode* aNodes[8];
1089 const SMDS_MeshElement* newElem = 0;
1090 newElem = aMesh->AddFace( aNodes[0], aNodes[1], aNodes[2], aNodes[3],
1091 aNodes[4], aNodes[5], aNodes[6], aNodes[7]);
1092 myLastCreatedElems.Append(newElem);
1093 AddToSameGroups( newElem, tr1, aMesh );
1094 int aShapeId = tr1->getshapeId();
1097 aMesh->SetMeshElementOnShape( newElem, aShapeId );
1099 aMesh->RemoveElement( tr1 );
1100 aMesh->RemoveElement( tr2 );
1102 // remove middle node (9)
1103 GetMeshDS()->RemoveNode( N1[4] );
1108 //=======================================================================
1109 //function : Reorient
1110 //purpose : Reverse theElement orientation
1111 //=======================================================================
1113 bool SMESH_MeshEditor::Reorient (const SMDS_MeshElement * theElem)
1115 myLastCreatedElems.Clear();
1116 myLastCreatedNodes.Clear();
1120 SMDS_ElemIteratorPtr it = theElem->nodesIterator();
1121 if ( !it || !it->more() )
1124 const SMDSAbs_ElementType type = theElem->GetType();
1125 if ( type < SMDSAbs_Edge || type > SMDSAbs_Volume )
1128 const SMDSAbs_EntityType geomType = theElem->GetEntityType();
1129 if ( geomType == SMDSEntity_Polyhedra ) // polyhedron
1131 const SMDS_VtkVolume* aPolyedre =
1132 dynamic_cast<const SMDS_VtkVolume*>( theElem );
1134 MESSAGE("Warning: bad volumic element");
1137 const int nbFaces = aPolyedre->NbFaces();
1138 vector<const SMDS_MeshNode *> poly_nodes;
1139 vector<int> quantities (nbFaces);
1141 // reverse each face of the polyedre
1142 for (int iface = 1; iface <= nbFaces; iface++) {
1143 int inode, nbFaceNodes = aPolyedre->NbFaceNodes(iface);
1144 quantities[iface - 1] = nbFaceNodes;
1146 for (inode = nbFaceNodes; inode >= 1; inode--) {
1147 const SMDS_MeshNode* curNode = aPolyedre->GetFaceNode(iface, inode);
1148 poly_nodes.push_back(curNode);
1151 return GetMeshDS()->ChangePolyhedronNodes( theElem, poly_nodes, quantities );
1153 else // other elements
1155 vector<const SMDS_MeshNode*> nodes( theElem->begin_nodes(), theElem->end_nodes() );
1156 const std::vector<int>& interlace = SMDS_MeshCell::reverseSmdsOrder( geomType, nodes.size() );
1157 if ( interlace.empty() )
1159 std::reverse( nodes.begin(), nodes.end() ); // obsolete, just in case
1163 SMDS_MeshCell::applyInterlace( interlace, nodes );
1165 return GetMeshDS()->ChangeElementNodes( theElem, &nodes[0], nodes.size() );
1170 //================================================================================
1172 * \brief Reorient faces.
1173 * \param theFaces - the faces to reorient. If empty the whole mesh is meant
1174 * \param theDirection - desired direction of normal of \a theFace
1175 * \param theFace - one of \a theFaces that sould be oriented according to
1176 * \a theDirection and whose orientation defines orientation of other faces
1177 * \return number of reoriented faces.
1179 //================================================================================
1181 int SMESH_MeshEditor::Reorient2D (TIDSortedElemSet & theFaces,
1182 const gp_Dir& theDirection,
1183 const SMDS_MeshElement * theFace)
1186 if ( !theFace || theFace->GetType() != SMDSAbs_Face ) return nbReori;
1188 if ( theFaces.empty() )
1190 SMDS_FaceIteratorPtr fIt = GetMeshDS()->facesIterator(/*idInceasingOrder=*/true);
1191 while ( fIt->more() )
1192 theFaces.insert( theFaces.end(), fIt->next() );
1195 // orient theFace according to theDirection
1197 SMESH_MeshAlgos::FaceNormal( theFace, normal, /*normalized=*/false );
1198 if ( normal * theDirection.XYZ() < 0 )
1199 nbReori += Reorient( theFace );
1201 // Orient other faces
1203 set< const SMDS_MeshElement* > startFaces, visitedFaces;
1204 TIDSortedElemSet avoidSet;
1205 set< SMESH_TLink > checkedLinks;
1206 pair< set< SMESH_TLink >::iterator, bool > linkIt_isNew;
1208 if ( theFaces.size() > 1 )// leave 1 face to prevent finding not selected faces
1209 theFaces.erase( theFace );
1210 startFaces.insert( theFace );
1212 int nodeInd1, nodeInd2;
1213 const SMDS_MeshElement* otherFace;
1214 vector< const SMDS_MeshElement* > facesNearLink;
1215 vector< std::pair< int, int > > nodeIndsOfFace;
1217 set< const SMDS_MeshElement* >::iterator startFace = startFaces.begin();
1218 while ( !startFaces.empty() )
1220 startFace = startFaces.begin();
1221 theFace = *startFace;
1222 startFaces.erase( startFace );
1223 if ( !visitedFaces.insert( theFace ).second )
1227 avoidSet.insert(theFace);
1229 NLink link( theFace->GetNode( 0 ), (SMDS_MeshNode *) 0 );
1231 const int nbNodes = theFace->NbCornerNodes();
1232 for ( int i = 0; i < nbNodes; ++i ) // loop on links of theFace
1234 link.second = theFace->GetNode(( i+1 ) % nbNodes );
1235 linkIt_isNew = checkedLinks.insert( link );
1236 if ( !linkIt_isNew.second )
1238 // link has already been checked and won't be encountered more
1239 // if the group (theFaces) is manifold
1240 //checkedLinks.erase( linkIt_isNew.first );
1244 facesNearLink.clear();
1245 nodeIndsOfFace.clear();
1246 while (( otherFace = SMESH_MeshAlgos::FindFaceInSet( link.first, link.second,
1248 &nodeInd1, &nodeInd2 )))
1249 if ( otherFace != theFace)
1251 facesNearLink.push_back( otherFace );
1252 nodeIndsOfFace.push_back( make_pair( nodeInd1, nodeInd2 ));
1253 avoidSet.insert( otherFace );
1255 if ( facesNearLink.size() > 1 )
1257 // NON-MANIFOLD mesh shell !
1258 // select a face most co-directed with theFace,
1259 // other faces won't be visited this time
1261 SMESH_MeshAlgos::FaceNormal( theFace, NF, /*normalized=*/false );
1262 double proj, maxProj = -1;
1263 for ( size_t i = 0; i < facesNearLink.size(); ++i ) {
1264 SMESH_MeshAlgos::FaceNormal( facesNearLink[i], NOF, /*normalized=*/false );
1265 if (( proj = Abs( NF * NOF )) > maxProj ) {
1267 otherFace = facesNearLink[i];
1268 nodeInd1 = nodeIndsOfFace[i].first;
1269 nodeInd2 = nodeIndsOfFace[i].second;
1272 // not to visit rejected faces
1273 for ( size_t i = 0; i < facesNearLink.size(); ++i )
1274 if ( facesNearLink[i] != otherFace && theFaces.size() > 1 )
1275 visitedFaces.insert( facesNearLink[i] );
1277 else if ( facesNearLink.size() == 1 )
1279 otherFace = facesNearLink[0];
1280 nodeInd1 = nodeIndsOfFace.back().first;
1281 nodeInd2 = nodeIndsOfFace.back().second;
1283 if ( otherFace && otherFace != theFace)
1285 // link must be reverse in otherFace if orientation ot otherFace
1286 // is same as that of theFace
1287 if ( abs(nodeInd2-nodeInd1) == 1 ? nodeInd2 > nodeInd1 : nodeInd1 > nodeInd2 )
1289 nbReori += Reorient( otherFace );
1291 startFaces.insert( otherFace );
1294 std::swap( link.first, link.second ); // reverse the link
1300 //================================================================================
1302 * \brief Reorient faces basing on orientation of adjacent volumes.
1303 * \param theFaces - faces to reorient. If empty, all mesh faces are treated.
1304 * \param theVolumes - reference volumes.
1305 * \param theOutsideNormal - to orient faces to have their normal
1306 * pointing either \a outside or \a inside the adjacent volumes.
1307 * \return number of reoriented faces.
1309 //================================================================================
1311 int SMESH_MeshEditor::Reorient2DBy3D (TIDSortedElemSet & theFaces,
1312 TIDSortedElemSet & theVolumes,
1313 const bool theOutsideNormal)
1317 SMDS_ElemIteratorPtr faceIt;
1318 if ( theFaces.empty() )
1319 faceIt = GetMeshDS()->elementsIterator( SMDSAbs_Face );
1321 faceIt = elemSetIterator( theFaces );
1323 vector< const SMDS_MeshNode* > faceNodes;
1324 TIDSortedElemSet checkedVolumes;
1325 set< const SMDS_MeshNode* > faceNodesSet;
1326 SMDS_VolumeTool volumeTool;
1328 while ( faceIt->more() ) // loop on given faces
1330 const SMDS_MeshElement* face = faceIt->next();
1331 if ( face->GetType() != SMDSAbs_Face )
1334 const size_t nbCornersNodes = face->NbCornerNodes();
1335 faceNodes.assign( face->begin_nodes(), face->end_nodes() );
1337 checkedVolumes.clear();
1338 SMDS_ElemIteratorPtr vIt = faceNodes[ 0 ]->GetInverseElementIterator( SMDSAbs_Volume );
1339 while ( vIt->more() )
1341 const SMDS_MeshElement* volume = vIt->next();
1343 if ( !checkedVolumes.insert( volume ).second )
1345 if ( !theVolumes.empty() && !theVolumes.count( volume ))
1348 // is volume adjacent?
1349 bool allNodesCommon = true;
1350 for ( size_t iN = 1; iN < nbCornersNodes && allNodesCommon; ++iN )
1351 allNodesCommon = ( volume->GetNodeIndex( faceNodes[ iN ]) > -1 );
1352 if ( !allNodesCommon )
1355 // get nodes of a corresponding volume facet
1356 faceNodesSet.clear();
1357 faceNodesSet.insert( faceNodes.begin(), faceNodes.end() );
1358 volumeTool.Set( volume );
1359 int facetID = volumeTool.GetFaceIndex( faceNodesSet );
1360 if ( facetID < 0 ) continue;
1361 volumeTool.SetExternalNormal();
1362 const SMDS_MeshNode** facetNodes = volumeTool.GetFaceNodes( facetID );
1364 // compare order of faceNodes and facetNodes
1365 const int iQ = 1 + ( nbCornersNodes < faceNodes.size() );
1367 for ( int i = 0; i < 2; ++i )
1369 const SMDS_MeshNode* n = facetNodes[ i*iQ ];
1370 for ( size_t iN = 0; iN < nbCornersNodes; ++iN )
1371 if ( faceNodes[ iN ] == n )
1377 bool isOutside = Abs( iNN[0]-iNN[1] ) == 1 ? iNN[0] < iNN[1] : iNN[0] > iNN[1];
1378 if ( isOutside != theOutsideNormal )
1379 nbReori += Reorient( face );
1381 } // loop on given faces
1386 //=======================================================================
1387 //function : getBadRate
1389 //=======================================================================
1391 static double getBadRate (const SMDS_MeshElement* theElem,
1392 SMESH::Controls::NumericalFunctorPtr& theCrit)
1394 SMESH::Controls::TSequenceOfXYZ P;
1395 if ( !theElem || !theCrit->GetPoints( theElem, P ))
1397 return theCrit->GetBadRate( theCrit->GetValue( P ), theElem->NbNodes() );
1398 //return theCrit->GetBadRate( theCrit->GetValue( theElem->GetID() ), theElem->NbNodes() );
1401 //=======================================================================
1402 //function : QuadToTri
1403 //purpose : Cut quadrangles into triangles.
1404 // theCrit is used to select a diagonal to cut
1405 //=======================================================================
1407 bool SMESH_MeshEditor::QuadToTri (TIDSortedElemSet & theElems,
1408 SMESH::Controls::NumericalFunctorPtr theCrit)
1410 myLastCreatedElems.Clear();
1411 myLastCreatedNodes.Clear();
1413 if ( !theCrit.get() )
1416 SMESHDS_Mesh * aMesh = GetMeshDS();
1418 Handle(Geom_Surface) surface;
1419 SMESH_MesherHelper helper( *GetMesh() );
1421 TIDSortedElemSet::iterator itElem;
1422 for ( itElem = theElems.begin(); itElem != theElems.end(); itElem++ )
1424 const SMDS_MeshElement* elem = *itElem;
1425 if ( !elem || elem->GetType() != SMDSAbs_Face )
1427 if ( elem->NbCornerNodes() != 4 )
1430 // retrieve element nodes
1431 vector< const SMDS_MeshNode* > aNodes( elem->begin_nodes(), elem->end_nodes() );
1433 // compare two sets of possible triangles
1434 double aBadRate1, aBadRate2; // to what extent a set is bad
1435 SMDS_FaceOfNodes tr1 ( aNodes[0], aNodes[1], aNodes[2] );
1436 SMDS_FaceOfNodes tr2 ( aNodes[2], aNodes[3], aNodes[0] );
1437 aBadRate1 = getBadRate( &tr1, theCrit ) + getBadRate( &tr2, theCrit );
1439 SMDS_FaceOfNodes tr3 ( aNodes[1], aNodes[2], aNodes[3] );
1440 SMDS_FaceOfNodes tr4 ( aNodes[3], aNodes[0], aNodes[1] );
1441 aBadRate2 = getBadRate( &tr3, theCrit ) + getBadRate( &tr4, theCrit );
1443 const int aShapeId = FindShape( elem );
1444 const SMDS_MeshElement* newElem1 = 0;
1445 const SMDS_MeshElement* newElem2 = 0;
1447 if ( !elem->IsQuadratic() ) // split liner quadrangle
1449 // for MaxElementLength2D functor we return minimum diagonal for splitting,
1450 // because aBadRate1=2*len(diagonal 1-3); aBadRate2=2*len(diagonal 2-4)
1451 if ( aBadRate1 <= aBadRate2 ) {
1452 // tr1 + tr2 is better
1453 newElem1 = aMesh->AddFace( aNodes[2], aNodes[3], aNodes[0] );
1454 newElem2 = aMesh->AddFace( aNodes[2], aNodes[0], aNodes[1] );
1457 // tr3 + tr4 is better
1458 newElem1 = aMesh->AddFace( aNodes[3], aNodes[0], aNodes[1] );
1459 newElem2 = aMesh->AddFace( aNodes[3], aNodes[1], aNodes[2] );
1462 else // split quadratic quadrangle
1464 helper.SetIsQuadratic( true );
1465 helper.SetIsBiQuadratic( aNodes.size() == 9 );
1467 helper.AddTLinks( static_cast< const SMDS_MeshFace* >( elem ));
1468 if ( aNodes.size() == 9 )
1470 helper.SetIsBiQuadratic( true );
1471 if ( aBadRate1 <= aBadRate2 )
1472 helper.AddTLinkNode( aNodes[0], aNodes[2], aNodes[8] );
1474 helper.AddTLinkNode( aNodes[1], aNodes[3], aNodes[8] );
1476 // create a new element
1477 if ( aBadRate1 <= aBadRate2 ) {
1478 newElem1 = helper.AddFace( aNodes[2], aNodes[3], aNodes[0] );
1479 newElem2 = helper.AddFace( aNodes[2], aNodes[0], aNodes[1] );
1482 newElem1 = helper.AddFace( aNodes[3], aNodes[0], aNodes[1] );
1483 newElem2 = helper.AddFace( aNodes[3], aNodes[1], aNodes[2] );
1487 // care of a new element
1489 myLastCreatedElems.Append(newElem1);
1490 myLastCreatedElems.Append(newElem2);
1491 AddToSameGroups( newElem1, elem, aMesh );
1492 AddToSameGroups( newElem2, elem, aMesh );
1494 // put a new triangle on the same shape
1496 aMesh->SetMeshElementOnShape( newElem1, aShapeId );
1497 aMesh->SetMeshElementOnShape( newElem2, aShapeId );
1499 aMesh->RemoveElement( elem );
1504 //=======================================================================
1506 * \brief Split each of given quadrangles into 4 triangles.
1507 * \param theElems - The faces to be splitted. If empty all faces are split.
1509 //=======================================================================
1511 void SMESH_MeshEditor::QuadTo4Tri (TIDSortedElemSet & theElems)
1513 myLastCreatedElems.Clear();
1514 myLastCreatedNodes.Clear();
1516 SMESH_MesherHelper helper( *GetMesh() );
1517 helper.SetElementsOnShape( true );
1519 SMDS_ElemIteratorPtr faceIt;
1520 if ( theElems.empty() ) faceIt = GetMeshDS()->elementsIterator(SMDSAbs_Face);
1521 else faceIt = elemSetIterator( theElems );
1524 gp_XY uv [9]; uv[8] = gp_XY(0,0);
1526 vector< const SMDS_MeshNode* > nodes;
1527 SMESHDS_SubMesh* subMeshDS = 0;
1529 Handle(Geom_Surface) surface;
1530 TopLoc_Location loc;
1532 while ( faceIt->more() )
1534 const SMDS_MeshElement* quad = faceIt->next();
1535 if ( !quad || quad->NbCornerNodes() != 4 )
1538 // get a surface the quad is on
1540 if ( quad->getshapeId() < 1 )
1543 helper.SetSubShape( 0 );
1546 else if ( quad->getshapeId() != helper.GetSubShapeID() )
1548 helper.SetSubShape( quad->getshapeId() );
1549 if ( !helper.GetSubShape().IsNull() &&
1550 helper.GetSubShape().ShapeType() == TopAbs_FACE )
1552 F = TopoDS::Face( helper.GetSubShape() );
1553 surface = BRep_Tool::Surface( F, loc );
1554 subMeshDS = GetMeshDS()->MeshElements( quad->getshapeId() );
1558 helper.SetSubShape( 0 );
1563 // create a central node
1565 const SMDS_MeshNode* nCentral;
1566 nodes.assign( quad->begin_nodes(), quad->end_nodes() );
1568 if ( nodes.size() == 9 )
1570 nCentral = nodes.back();
1577 for ( ; iN < nodes.size(); ++iN )
1578 xyz[ iN ] = SMESH_TNodeXYZ( nodes[ iN ] );
1580 for ( ; iN < 8; ++iN ) // mid-side points of a linear qudrangle
1581 xyz[ iN ] = 0.5 * ( xyz[ iN - 4 ] + xyz[( iN - 3 )%4 ] );
1583 xyz[ 8 ] = helper.calcTFI( 0.5, 0.5,
1584 xyz[0], xyz[1], xyz[2], xyz[3],
1585 xyz[4], xyz[5], xyz[6], xyz[7] );
1589 for ( ; iN < nodes.size(); ++iN )
1590 uv[ iN ] = helper.GetNodeUV( F, nodes[iN], nodes[(iN+2)%4], &checkUV );
1592 for ( ; iN < 8; ++iN ) // UV of mid-side points of a linear qudrangle
1593 uv[ iN ] = helper.GetMiddleUV( surface, uv[ iN - 4 ], uv[( iN - 3 )%4 ] );
1595 uv[ 8 ] = helper.calcTFI( 0.5, 0.5,
1596 uv[0], uv[1], uv[2], uv[3],
1597 uv[4], uv[5], uv[6], uv[7] );
1599 gp_Pnt p = surface->Value( uv[8].X(), uv[8].Y() ).Transformed( loc );
1603 nCentral = helper.AddNode( xyz[8].X(), xyz[8].Y(), xyz[8].Z(), /*id=*/0,
1604 uv[8].X(), uv[8].Y() );
1605 myLastCreatedNodes.Append( nCentral );
1608 // create 4 triangles
1610 helper.SetIsQuadratic ( nodes.size() > 4 );
1611 helper.SetIsBiQuadratic( nodes.size() == 9 );
1612 if ( helper.GetIsQuadratic() )
1613 helper.AddTLinks( static_cast< const SMDS_MeshFace*>( quad ));
1615 GetMeshDS()->RemoveFreeElement( quad, subMeshDS, /*fromGroups=*/false );
1617 for ( int i = 0; i < 4; ++i )
1619 SMDS_MeshElement* tria = helper.AddFace( nodes[ i ],
1622 ReplaceElemInGroups( tria, quad, GetMeshDS() );
1623 myLastCreatedElems.Append( tria );
1628 //=======================================================================
1629 //function : BestSplit
1630 //purpose : Find better diagonal for cutting.
1631 //=======================================================================
1633 int SMESH_MeshEditor::BestSplit (const SMDS_MeshElement* theQuad,
1634 SMESH::Controls::NumericalFunctorPtr theCrit)
1636 myLastCreatedElems.Clear();
1637 myLastCreatedNodes.Clear();
1642 if (!theQuad || theQuad->GetType() != SMDSAbs_Face )
1645 if( theQuad->NbNodes()==4 ||
1646 (theQuad->NbNodes()==8 && theQuad->IsQuadratic()) ) {
1648 // retrieve element nodes
1649 const SMDS_MeshNode* aNodes [4];
1650 SMDS_ElemIteratorPtr itN = theQuad->nodesIterator();
1652 //while (itN->more())
1654 aNodes[ i++ ] = static_cast<const SMDS_MeshNode*>( itN->next() );
1656 // compare two sets of possible triangles
1657 double aBadRate1, aBadRate2; // to what extent a set is bad
1658 SMDS_FaceOfNodes tr1 ( aNodes[0], aNodes[1], aNodes[2] );
1659 SMDS_FaceOfNodes tr2 ( aNodes[2], aNodes[3], aNodes[0] );
1660 aBadRate1 = getBadRate( &tr1, theCrit ) + getBadRate( &tr2, theCrit );
1662 SMDS_FaceOfNodes tr3 ( aNodes[1], aNodes[2], aNodes[3] );
1663 SMDS_FaceOfNodes tr4 ( aNodes[3], aNodes[0], aNodes[1] );
1664 aBadRate2 = getBadRate( &tr3, theCrit ) + getBadRate( &tr4, theCrit );
1665 // for MaxElementLength2D functor we return minimum diagonal for splitting,
1666 // because aBadRate1=2*len(diagonal 1-3); aBadRate2=2*len(diagonal 2-4)
1667 if (aBadRate1 <= aBadRate2) // tr1 + tr2 is better
1668 return 1; // diagonal 1-3
1670 return 2; // diagonal 2-4
1677 // Methods of splitting volumes into tetra
1679 const int theHexTo5_1[5*4+1] =
1681 0, 1, 2, 5, 0, 4, 5, 7, 0, 2, 3, 7, 2, 5, 6, 7, 0, 5, 2, 7, -1
1683 const int theHexTo5_2[5*4+1] =
1685 1, 2, 3, 6, 1, 4, 5, 6, 0, 1, 3, 4, 3, 4, 6, 7, 1, 3, 4, 6, -1
1687 const int* theHexTo5[2] = { theHexTo5_1, theHexTo5_2 };
1689 const int theHexTo6_1[6*4+1] =
1691 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
1693 const int theHexTo6_2[6*4+1] =
1695 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
1697 const int theHexTo6_3[6*4+1] =
1699 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
1701 const int theHexTo6_4[6*4+1] =
1703 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
1705 const int* theHexTo6[4] = { theHexTo6_1, theHexTo6_2, theHexTo6_3, theHexTo6_4 };
1707 const int thePyraTo2_1[2*4+1] =
1709 0, 1, 2, 4, 0, 2, 3, 4, -1
1711 const int thePyraTo2_2[2*4+1] =
1713 1, 2, 3, 4, 1, 3, 0, 4, -1
1715 const int* thePyraTo2[2] = { thePyraTo2_1, thePyraTo2_2 };
1717 const int thePentaTo3_1[3*4+1] =
1719 0, 1, 2, 3, 1, 3, 4, 2, 2, 3, 4, 5, -1
1721 const int thePentaTo3_2[3*4+1] =
1723 1, 2, 0, 4, 2, 4, 5, 0, 0, 4, 5, 3, -1
1725 const int thePentaTo3_3[3*4+1] =
1727 2, 0, 1, 5, 0, 5, 3, 1, 1, 5, 3, 4, -1
1729 const int thePentaTo3_4[3*4+1] =
1731 0, 1, 2, 3, 1, 3, 4, 5, 2, 3, 1, 5, -1
1733 const int thePentaTo3_5[3*4+1] =
1735 1, 2, 0, 4, 2, 4, 5, 3, 0, 4, 2, 3, -1
1737 const int thePentaTo3_6[3*4+1] =
1739 2, 0, 1, 5, 0, 5, 3, 4, 1, 5, 0, 4, -1
1741 const int* thePentaTo3[6] = { thePentaTo3_1, thePentaTo3_2, thePentaTo3_3,
1742 thePentaTo3_4, thePentaTo3_5, thePentaTo3_6 };
1744 // Methods of splitting hexahedron into prisms
1746 const int theHexTo4Prisms_BT[6*4+1] = // bottom-top
1748 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
1750 const int theHexTo4Prisms_LR[6*4+1] = // left-right
1752 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
1754 const int theHexTo4Prisms_FB[6*4+1] = // front-back
1756 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
1759 const int theHexTo2Prisms_BT_1[6*2+1] =
1761 0, 1, 3, 4, 5, 7, 1, 2, 3, 5, 6, 7, -1
1763 const int theHexTo2Prisms_BT_2[6*2+1] =
1765 0, 1, 2, 4, 5, 6, 0, 2, 3, 4, 6, 7, -1
1767 const int* theHexTo2Prisms_BT[2] = { theHexTo2Prisms_BT_1, theHexTo2Prisms_BT_2 };
1769 const int theHexTo2Prisms_LR_1[6*2+1] =
1771 1, 0, 4, 2, 3, 7, 1, 4, 5, 2, 7, 6, -1
1773 const int theHexTo2Prisms_LR_2[6*2+1] =
1775 1, 0, 4, 2, 3, 7, 1, 4, 5, 2, 7, 6, -1
1777 const int* theHexTo2Prisms_LR[2] = { theHexTo2Prisms_LR_1, theHexTo2Prisms_LR_2 };
1779 const int theHexTo2Prisms_FB_1[6*2+1] =
1781 0, 3, 4, 1, 2, 5, 3, 7, 4, 2, 6, 5, -1
1783 const int theHexTo2Prisms_FB_2[6*2+1] =
1785 0, 3, 7, 1, 2, 7, 0, 7, 4, 1, 6, 5, -1
1787 const int* theHexTo2Prisms_FB[2] = { theHexTo2Prisms_FB_1, theHexTo2Prisms_FB_2 };
1790 struct TTriangleFacet //!< stores indices of three nodes of tetra facet
1793 TTriangleFacet(int n1, int n2, int n3): _n1(n1), _n2(n2), _n3(n3) {}
1794 bool contains(int n) const { return ( n == _n1 || n == _n2 || n == _n3 ); }
1795 bool hasAdjacentVol( const SMDS_MeshElement* elem,
1796 const SMDSAbs_GeometryType geom = SMDSGeom_TETRA) const;
1802 const int* _connectivity; //!< foursomes of tetra connectivy finished by -1
1803 bool _baryNode; //!< additional node is to be created at cell barycenter
1804 bool _ownConn; //!< to delete _connectivity in destructor
1805 map<int, const SMDS_MeshNode*> _faceBaryNode; //!< map face index to node at BC of face
1807 TSplitMethod( int nbTet=0, const int* conn=0, bool addNode=false)
1808 : _nbSplits(nbTet), _nbCorners(4), _connectivity(conn), _baryNode(addNode), _ownConn(false) {}
1809 ~TSplitMethod() { if ( _ownConn ) delete [] _connectivity; _connectivity = 0; }
1810 bool hasFacet( const TTriangleFacet& facet ) const
1812 if ( _nbCorners == 4 )
1814 const int* tetConn = _connectivity;
1815 for ( ; tetConn[0] >= 0; tetConn += 4 )
1816 if (( facet.contains( tetConn[0] ) +
1817 facet.contains( tetConn[1] ) +
1818 facet.contains( tetConn[2] ) +
1819 facet.contains( tetConn[3] )) == 3 )
1822 else // prism, _nbCorners == 6
1824 const int* prismConn = _connectivity;
1825 for ( ; prismConn[0] >= 0; prismConn += 6 )
1827 if (( facet.contains( prismConn[0] ) &&
1828 facet.contains( prismConn[1] ) &&
1829 facet.contains( prismConn[2] ))
1831 ( facet.contains( prismConn[3] ) &&
1832 facet.contains( prismConn[4] ) &&
1833 facet.contains( prismConn[5] )))
1841 //=======================================================================
1843 * \brief return TSplitMethod for the given element to split into tetrahedra
1845 //=======================================================================
1847 TSplitMethod getTetraSplitMethod( SMDS_VolumeTool& vol, const int theMethodFlags)
1849 const int iQ = vol.Element()->IsQuadratic() ? 2 : 1;
1851 // at HEXA_TO_24 method, each face of volume is split into triangles each based on
1852 // an edge and a face barycenter; tertaherdons are based on triangles and
1853 // a volume barycenter
1854 const bool is24TetMode = ( theMethodFlags == SMESH_MeshEditor::HEXA_TO_24 );
1856 // Find out how adjacent volumes are split
1858 vector < list< TTriangleFacet > > triaSplitsByFace( vol.NbFaces() ); // splits of each side
1859 int hasAdjacentSplits = 0, maxTetConnSize = 0;
1860 for ( int iF = 0; iF < vol.NbFaces(); ++iF )
1862 int nbNodes = vol.NbFaceNodes( iF ) / iQ;
1863 maxTetConnSize += 4 * ( nbNodes - (is24TetMode ? 0 : 2));
1864 if ( nbNodes < 4 ) continue;
1866 list< TTriangleFacet >& triaSplits = triaSplitsByFace[ iF ];
1867 const int* nInd = vol.GetFaceNodesIndices( iF );
1870 TTriangleFacet t012( nInd[0*iQ], nInd[1*iQ], nInd[2*iQ] );
1871 TTriangleFacet t123( nInd[1*iQ], nInd[2*iQ], nInd[3*iQ] );
1872 if ( t012.hasAdjacentVol( vol.Element() )) triaSplits.push_back( t012 );
1873 else if ( t123.hasAdjacentVol( vol.Element() )) triaSplits.push_back( t123 );
1877 int iCom = 0; // common node of triangle faces to split into
1878 for ( int iVar = 0; iVar < nbNodes; ++iVar, ++iCom )
1880 TTriangleFacet t012( nInd[ iQ * ( iCom )],
1881 nInd[ iQ * ( (iCom+1)%nbNodes )],
1882 nInd[ iQ * ( (iCom+2)%nbNodes )]);
1883 TTriangleFacet t023( nInd[ iQ * ( iCom )],
1884 nInd[ iQ * ( (iCom+2)%nbNodes )],
1885 nInd[ iQ * ( (iCom+3)%nbNodes )]);
1886 if ( t012.hasAdjacentVol( vol.Element() ) && t023.hasAdjacentVol( vol.Element() ))
1888 triaSplits.push_back( t012 );
1889 triaSplits.push_back( t023 );
1894 if ( !triaSplits.empty() )
1895 hasAdjacentSplits = true;
1898 // Among variants of split method select one compliant with adjacent volumes
1900 TSplitMethod method;
1901 if ( !vol.Element()->IsPoly() && !is24TetMode )
1903 int nbVariants = 2, nbTet = 0;
1904 const int** connVariants = 0;
1905 switch ( vol.Element()->GetEntityType() )
1907 case SMDSEntity_Hexa:
1908 case SMDSEntity_Quad_Hexa:
1909 case SMDSEntity_TriQuad_Hexa:
1910 if ( theMethodFlags == SMESH_MeshEditor::HEXA_TO_5 )
1911 connVariants = theHexTo5, nbTet = 5;
1913 connVariants = theHexTo6, nbTet = 6, nbVariants = 4;
1915 case SMDSEntity_Pyramid:
1916 case SMDSEntity_Quad_Pyramid:
1917 connVariants = thePyraTo2; nbTet = 2;
1919 case SMDSEntity_Penta:
1920 case SMDSEntity_Quad_Penta:
1921 connVariants = thePentaTo3; nbTet = 3; nbVariants = 6;
1926 for ( int variant = 0; variant < nbVariants && method._nbSplits == 0; ++variant )
1928 // check method compliancy with adjacent tetras,
1929 // all found splits must be among facets of tetras described by this method
1930 method = TSplitMethod( nbTet, connVariants[variant] );
1931 if ( hasAdjacentSplits && method._nbSplits > 0 )
1933 bool facetCreated = true;
1934 for ( size_t iF = 0; facetCreated && iF < triaSplitsByFace.size(); ++iF )
1936 list< TTriangleFacet >::const_iterator facet = triaSplitsByFace[iF].begin();
1937 for ( ; facetCreated && facet != triaSplitsByFace[iF].end(); ++facet )
1938 facetCreated = method.hasFacet( *facet );
1940 if ( !facetCreated )
1941 method = TSplitMethod(0); // incompatible method
1945 if ( method._nbSplits < 1 )
1947 // No standard method is applicable, use a generic solution:
1948 // each facet of a volume is split into triangles and
1949 // each of triangles and a volume barycenter form a tetrahedron.
1951 const bool isHex27 = ( vol.Element()->GetEntityType() == SMDSEntity_TriQuad_Hexa );
1953 int* connectivity = new int[ maxTetConnSize + 1 ];
1954 method._connectivity = connectivity;
1955 method._ownConn = true;
1956 method._baryNode = !isHex27; // to create central node or not
1959 int baryCenInd = vol.NbNodes() - int( isHex27 );
1960 for ( int iF = 0; iF < vol.NbFaces(); ++iF )
1962 const int nbNodes = vol.NbFaceNodes( iF ) / iQ;
1963 const int* nInd = vol.GetFaceNodesIndices( iF );
1964 // find common node of triangle facets of tetra to create
1965 int iCommon = 0; // index in linear numeration
1966 const list< TTriangleFacet >& triaSplits = triaSplitsByFace[ iF ];
1967 if ( !triaSplits.empty() )
1970 const TTriangleFacet* facet = &triaSplits.front();
1971 for ( ; iCommon < nbNodes-1 ; ++iCommon )
1972 if ( facet->contains( nInd[ iQ * iCommon ]) &&
1973 facet->contains( nInd[ iQ * ((iCommon+2)%nbNodes) ]))
1976 else if ( nbNodes > 3 && !is24TetMode )
1978 // find the best method of splitting into triangles by aspect ratio
1979 SMESH::Controls::NumericalFunctorPtr aspectRatio( new SMESH::Controls::AspectRatio);
1980 map< double, int > badness2iCommon;
1981 const SMDS_MeshNode** nodes = vol.GetFaceNodes( iF );
1982 int nbVariants = ( nbNodes == 4 ? 2 : nbNodes );
1983 for ( int iVar = 0; iVar < nbVariants; ++iVar, ++iCommon )
1986 for ( int iLast = iCommon+2; iLast < iCommon+nbNodes; ++iLast )
1988 SMDS_FaceOfNodes tria ( nodes[ iQ*( iCommon )],
1989 nodes[ iQ*((iLast-1)%nbNodes)],
1990 nodes[ iQ*((iLast )%nbNodes)]);
1991 badness += getBadRate( &tria, aspectRatio );
1993 badness2iCommon.insert( make_pair( badness, iCommon ));
1995 // use iCommon with lowest badness
1996 iCommon = badness2iCommon.begin()->second;
1998 if ( iCommon >= nbNodes )
1999 iCommon = 0; // something wrong
2001 // fill connectivity of tetrahedra based on a current face
2002 int nbTet = nbNodes - 2;
2003 if ( is24TetMode && nbNodes > 3 && triaSplits.empty())
2008 faceBaryCenInd = vol.GetCenterNodeIndex( iF );
2009 method._faceBaryNode[ iF ] = vol.GetNodes()[ faceBaryCenInd ];
2013 method._faceBaryNode[ iF ] = 0;
2014 faceBaryCenInd = baryCenInd + method._faceBaryNode.size();
2017 for ( int i = 0; i < nbTet; ++i )
2019 int i1 = i, i2 = (i+1) % nbNodes;
2020 if ( !vol.IsFaceExternal( iF )) swap( i1, i2 );
2021 connectivity[ connSize++ ] = nInd[ iQ * i1 ];
2022 connectivity[ connSize++ ] = nInd[ iQ * i2 ];
2023 connectivity[ connSize++ ] = faceBaryCenInd;
2024 connectivity[ connSize++ ] = baryCenInd;
2029 for ( int i = 0; i < nbTet; ++i )
2031 int i1 = (iCommon+1+i) % nbNodes, i2 = (iCommon+2+i) % nbNodes;
2032 if ( !vol.IsFaceExternal( iF )) swap( i1, i2 );
2033 connectivity[ connSize++ ] = nInd[ iQ * iCommon ];
2034 connectivity[ connSize++ ] = nInd[ iQ * i1 ];
2035 connectivity[ connSize++ ] = nInd[ iQ * i2 ];
2036 connectivity[ connSize++ ] = baryCenInd;
2039 method._nbSplits += nbTet;
2041 } // loop on volume faces
2043 connectivity[ connSize++ ] = -1;
2045 } // end of generic solution
2049 //=======================================================================
2051 * \brief return TSplitMethod to split haxhedron into prisms
2053 //=======================================================================
2055 TSplitMethod getPrismSplitMethod( SMDS_VolumeTool& vol,
2056 const int methodFlags,
2057 const int facetToSplit)
2059 // order of facets in HEX according to SMDS_VolumeTool::Hexa_F :
2061 const int iF = ( facetToSplit < 2 ) ? 0 : 1 + ( facetToSplit-2 ) % 2; // [0,1,2]
2063 if ( methodFlags == SMESH_MeshEditor::HEXA_TO_4_PRISMS )
2065 static TSplitMethod to4methods[4]; // order BT, LR, FB
2066 if ( to4methods[iF]._nbSplits == 0 )
2070 to4methods[iF]._connectivity = theHexTo4Prisms_BT;
2071 to4methods[iF]._faceBaryNode[ 0 ] = 0;
2072 to4methods[iF]._faceBaryNode[ 1 ] = 0;
2075 to4methods[iF]._connectivity = theHexTo4Prisms_LR;
2076 to4methods[iF]._faceBaryNode[ 2 ] = 0;
2077 to4methods[iF]._faceBaryNode[ 4 ] = 0;
2080 to4methods[iF]._connectivity = theHexTo4Prisms_FB;
2081 to4methods[iF]._faceBaryNode[ 3 ] = 0;
2082 to4methods[iF]._faceBaryNode[ 5 ] = 0;
2084 default: return to4methods[3];
2086 to4methods[iF]._nbSplits = 4;
2087 to4methods[iF]._nbCorners = 6;
2089 return to4methods[iF];
2091 // else if ( methodFlags == HEXA_TO_2_PRISMS )
2093 TSplitMethod method;
2095 const int iQ = vol.Element()->IsQuadratic() ? 2 : 1;
2097 const int nbVariants = 2, nbSplits = 2;
2098 const int** connVariants = 0;
2100 case 0: connVariants = theHexTo2Prisms_BT; break;
2101 case 1: connVariants = theHexTo2Prisms_LR; break;
2102 case 2: connVariants = theHexTo2Prisms_FB; break;
2103 default: return method;
2106 // look for prisms adjacent via facetToSplit and an opposite one
2107 for ( int is2nd = 0; is2nd < 2; ++is2nd )
2109 int iFacet = is2nd ? vol.GetOppFaceIndexOfHex( facetToSplit ) : facetToSplit;
2110 int nbNodes = vol.NbFaceNodes( iFacet ) / iQ;
2111 if ( nbNodes != 4 ) return method;
2113 const int* nInd = vol.GetFaceNodesIndices( iFacet );
2114 TTriangleFacet t012( nInd[0*iQ], nInd[1*iQ], nInd[2*iQ] );
2115 TTriangleFacet t123( nInd[1*iQ], nInd[2*iQ], nInd[3*iQ] );
2117 if ( t012.hasAdjacentVol( vol.Element(), SMDSGeom_PENTA ))
2119 else if ( t123.hasAdjacentVol( vol.Element(), SMDSGeom_PENTA ))
2124 // there are adjacent prism
2125 for ( int variant = 0; variant < nbVariants; ++variant )
2127 // check method compliancy with adjacent prisms,
2128 // the found prism facets must be among facets of prisms described by current method
2129 method._nbSplits = nbSplits;
2130 method._nbCorners = 6;
2131 method._connectivity = connVariants[ variant ];
2132 if ( method.hasFacet( *t ))
2137 // No adjacent prisms. Select a variant with a best aspect ratio.
2139 double badness[2] = { 0., 0. };
2140 static SMESH::Controls::NumericalFunctorPtr aspectRatio( new SMESH::Controls::AspectRatio);
2141 const SMDS_MeshNode** nodes = vol.GetNodes();
2142 for ( int variant = 0; variant < nbVariants; ++variant )
2143 for ( int is2nd = 0; is2nd < 2; ++is2nd )
2145 int iFacet = is2nd ? vol.GetOppFaceIndexOfHex( facetToSplit ) : facetToSplit;
2146 const int* nInd = vol.GetFaceNodesIndices( iFacet );
2148 method._connectivity = connVariants[ variant ];
2149 TTriangleFacet t012( nInd[0*iQ], nInd[1*iQ], nInd[2*iQ] );
2150 TTriangleFacet t123( nInd[1*iQ], nInd[2*iQ], nInd[3*iQ] );
2151 TTriangleFacet* t = ( method.hasFacet( t012 )) ? & t012 : & t123;
2153 SMDS_FaceOfNodes tria ( nodes[ t->_n1 ],
2156 badness[ variant ] += getBadRate( &tria, aspectRatio );
2158 const int iBetter = ( badness[1] < badness[0] && badness[0]-badness[1] > 0.1 * badness[0] );
2160 method._nbSplits = nbSplits;
2161 method._nbCorners = 6;
2162 method._connectivity = connVariants[ iBetter ];
2167 //================================================================================
2169 * \brief Check if there is a tetraherdon adjacent to the given element via this facet
2171 //================================================================================
2173 bool TTriangleFacet::hasAdjacentVol( const SMDS_MeshElement* elem,
2174 const SMDSAbs_GeometryType geom ) const
2176 // find the tetrahedron including the three nodes of facet
2177 const SMDS_MeshNode* n1 = elem->GetNode(_n1);
2178 const SMDS_MeshNode* n2 = elem->GetNode(_n2);
2179 const SMDS_MeshNode* n3 = elem->GetNode(_n3);
2180 SMDS_ElemIteratorPtr volIt1 = n1->GetInverseElementIterator(SMDSAbs_Volume);
2181 while ( volIt1->more() )
2183 const SMDS_MeshElement* v = volIt1->next();
2184 if ( v->GetGeomType() != geom )
2186 const int lastCornerInd = v->NbCornerNodes() - 1;
2187 if ( v->IsQuadratic() && v->GetNodeIndex( n1 ) > lastCornerInd )
2188 continue; // medium node not allowed
2189 const int ind2 = v->GetNodeIndex( n2 );
2190 if ( ind2 < 0 || lastCornerInd < ind2 )
2192 const int ind3 = v->GetNodeIndex( n3 );
2193 if ( ind3 < 0 || lastCornerInd < ind3 )
2200 //=======================================================================
2202 * \brief A key of a face of volume
2204 //=======================================================================
2206 struct TVolumeFaceKey: pair< pair< int, int>, pair< int, int> >
2208 TVolumeFaceKey( SMDS_VolumeTool& vol, int iF )
2210 TIDSortedNodeSet sortedNodes;
2211 const int iQ = vol.Element()->IsQuadratic() ? 2 : 1;
2212 int nbNodes = vol.NbFaceNodes( iF );
2213 const SMDS_MeshNode** fNodes = vol.GetFaceNodes( iF );
2214 for ( int i = 0; i < nbNodes; i += iQ )
2215 sortedNodes.insert( fNodes[i] );
2216 TIDSortedNodeSet::iterator n = sortedNodes.begin();
2217 first.first = (*(n++))->GetID();
2218 first.second = (*(n++))->GetID();
2219 second.first = (*(n++))->GetID();
2220 second.second = ( sortedNodes.size() > 3 ) ? (*(n++))->GetID() : 0;
2225 //=======================================================================
2226 //function : SplitVolumes
2227 //purpose : Split volume elements into tetrahedra or prisms.
2228 // If facet ID < 0, element is split into tetrahedra,
2229 // else a hexahedron is split into prisms so that the given facet is
2230 // split into triangles
2231 //=======================================================================
2233 void SMESH_MeshEditor::SplitVolumes (const TFacetOfElem & theElems,
2234 const int theMethodFlags)
2236 SMDS_VolumeTool volTool;
2237 SMESH_MesherHelper helper( *GetMesh()), fHelper(*GetMesh());
2238 fHelper.ToFixNodeParameters( true );
2240 SMESHDS_SubMesh* subMesh = 0;//GetMeshDS()->MeshElements(1);
2241 SMESHDS_SubMesh* fSubMesh = 0;//subMesh;
2243 SMESH_SequenceOfElemPtr newNodes, newElems;
2245 // map face of volume to it's baricenrtic node
2246 map< TVolumeFaceKey, const SMDS_MeshNode* > volFace2BaryNode;
2248 vector<const SMDS_MeshElement* > splitVols;
2250 TFacetOfElem::const_iterator elem2facet = theElems.begin();
2251 for ( ; elem2facet != theElems.end(); ++elem2facet )
2253 const SMDS_MeshElement* elem = elem2facet->first;
2254 const int facetToSplit = elem2facet->second;
2255 if ( elem->GetType() != SMDSAbs_Volume )
2257 const SMDSAbs_EntityType geomType = elem->GetEntityType();
2258 if ( geomType == SMDSEntity_Tetra || geomType == SMDSEntity_Quad_Tetra )
2261 if ( !volTool.Set( elem, /*ignoreCentralNodes=*/false )) continue; // strange...
2263 TSplitMethod splitMethod = ( facetToSplit < 0 ?
2264 getTetraSplitMethod( volTool, theMethodFlags ) :
2265 getPrismSplitMethod( volTool, theMethodFlags, facetToSplit ));
2266 if ( splitMethod._nbSplits < 1 ) continue;
2268 // find submesh to add new tetras to
2269 if ( !subMesh || !subMesh->Contains( elem ))
2271 int shapeID = FindShape( elem );
2272 helper.SetSubShape( shapeID ); // helper will add tetras to the found submesh
2273 subMesh = GetMeshDS()->MeshElements( shapeID );
2276 if ( elem->IsQuadratic() )
2279 // add quadratic links to the helper
2280 for ( int iF = 0; iF < volTool.NbFaces(); ++iF )
2282 const SMDS_MeshNode** fNodes = volTool.GetFaceNodes( iF );
2283 int nbN = volTool.NbFaceNodes( iF ) - bool( volTool.GetCenterNodeIndex(iF) > 0 );
2284 for ( int iN = 0; iN < nbN; iN += iQ )
2285 helper.AddTLinkNode( fNodes[iN], fNodes[iN+2], fNodes[iN+1] );
2287 helper.SetIsQuadratic( true );
2292 helper.SetIsQuadratic( false );
2294 vector<const SMDS_MeshNode*> nodes( volTool.GetNodes(),
2295 volTool.GetNodes() + elem->NbNodes() );
2296 helper.SetElementsOnShape( true );
2297 if ( splitMethod._baryNode )
2299 // make a node at barycenter
2300 volTool.GetBaryCenter( bc[0], bc[1], bc[2] );
2301 SMDS_MeshNode* gcNode = helper.AddNode( bc[0], bc[1], bc[2] );
2302 nodes.push_back( gcNode );
2303 newNodes.Append( gcNode );
2305 if ( !splitMethod._faceBaryNode.empty() )
2307 // make or find baricentric nodes of faces
2308 map<int, const SMDS_MeshNode*>::iterator iF_n = splitMethod._faceBaryNode.begin();
2309 for ( ; iF_n != splitMethod._faceBaryNode.end(); ++iF_n )
2311 map< TVolumeFaceKey, const SMDS_MeshNode* >::iterator f_n =
2312 volFace2BaryNode.insert
2313 ( make_pair( TVolumeFaceKey( volTool,iF_n->first ), iF_n->second )).first;
2316 volTool.GetFaceBaryCenter( iF_n->first, bc[0], bc[1], bc[2] );
2317 newNodes.Append( f_n->second = helper.AddNode( bc[0], bc[1], bc[2] ));
2319 nodes.push_back( iF_n->second = f_n->second );
2324 splitVols.resize( splitMethod._nbSplits ); // splits of a volume
2325 const int* volConn = splitMethod._connectivity;
2326 if ( splitMethod._nbCorners == 4 ) // tetra
2327 for ( int i = 0; i < splitMethod._nbSplits; ++i, volConn += splitMethod._nbCorners )
2328 newElems.Append( splitVols[ i ] = helper.AddVolume( nodes[ volConn[0] ],
2329 nodes[ volConn[1] ],
2330 nodes[ volConn[2] ],
2331 nodes[ volConn[3] ]));
2333 for ( int i = 0; i < splitMethod._nbSplits; ++i, volConn += splitMethod._nbCorners )
2334 newElems.Append( splitVols[ i ] = helper.AddVolume( nodes[ volConn[0] ],
2335 nodes[ volConn[1] ],
2336 nodes[ volConn[2] ],
2337 nodes[ volConn[3] ],
2338 nodes[ volConn[4] ],
2339 nodes[ volConn[5] ]));
2341 ReplaceElemInGroups( elem, splitVols, GetMeshDS() );
2343 // Split faces on sides of the split volume
2345 const SMDS_MeshNode** volNodes = volTool.GetNodes();
2346 for ( int iF = 0; iF < volTool.NbFaces(); ++iF )
2348 const int nbNodes = volTool.NbFaceNodes( iF ) / iQ;
2349 if ( nbNodes < 4 ) continue;
2351 // find an existing face
2352 vector<const SMDS_MeshNode*> fNodes( volTool.GetFaceNodes( iF ),
2353 volTool.GetFaceNodes( iF ) + volTool.NbFaceNodes( iF ));
2354 while ( const SMDS_MeshElement* face = GetMeshDS()->FindElement( fNodes, SMDSAbs_Face,
2355 /*noMedium=*/false))
2358 helper.SetElementsOnShape( false );
2359 vector< const SMDS_MeshElement* > triangles;
2361 // find submesh to add new triangles in
2362 if ( !fSubMesh || !fSubMesh->Contains( face ))
2364 int shapeID = FindShape( face );
2365 fSubMesh = GetMeshDS()->MeshElements( shapeID );
2367 map<int, const SMDS_MeshNode*>::iterator iF_n = splitMethod._faceBaryNode.find(iF);
2368 if ( iF_n != splitMethod._faceBaryNode.end() )
2370 const SMDS_MeshNode *baryNode = iF_n->second;
2371 for ( int iN = 0; iN < nbNodes*iQ; iN += iQ )
2373 const SMDS_MeshNode* n1 = fNodes[iN];
2374 const SMDS_MeshNode *n2 = fNodes[(iN+iQ)%(nbNodes*iQ)];
2375 const SMDS_MeshNode *n3 = baryNode;
2376 if ( !volTool.IsFaceExternal( iF ))
2378 triangles.push_back( helper.AddFace( n1,n2,n3 ));
2380 if ( fSubMesh ) // update position of the bary node on geometry
2383 subMesh->RemoveNode( baryNode, false );
2384 GetMeshDS()->SetNodeOnFace( baryNode, fSubMesh->GetID() );
2385 const TopoDS_Shape& s = GetMeshDS()->IndexToShape( fSubMesh->GetID() );
2386 if ( !s.IsNull() && s.ShapeType() == TopAbs_FACE )
2388 fHelper.SetSubShape( s );
2389 gp_XY uv( 1e100, 1e100 );
2391 if ( !fHelper.CheckNodeUV( TopoDS::Face( s ), baryNode,
2392 uv, /*tol=*/1e-7, /*force=*/true, distXYZ ) &&
2395 // node is too far from the surface
2396 GetMeshDS()->MoveNode( baryNode, distXYZ[1], distXYZ[2], distXYZ[3] );
2397 const_cast<SMDS_MeshNode*>( baryNode )->SetPosition
2398 ( SMDS_PositionPtr( new SMDS_FacePosition( uv.X(), uv.Y() )));
2405 // among possible triangles create ones discribed by split method
2406 const int* nInd = volTool.GetFaceNodesIndices( iF );
2407 int nbVariants = ( nbNodes == 4 ? 2 : nbNodes );
2408 int iCom = 0; // common node of triangle faces to split into
2409 list< TTriangleFacet > facets;
2410 for ( int iVar = 0; iVar < nbVariants; ++iVar, ++iCom )
2412 TTriangleFacet t012( nInd[ iQ * ( iCom )],
2413 nInd[ iQ * ( (iCom+1)%nbNodes )],
2414 nInd[ iQ * ( (iCom+2)%nbNodes )]);
2415 TTriangleFacet t023( nInd[ iQ * ( iCom )],
2416 nInd[ iQ * ( (iCom+2)%nbNodes )],
2417 nInd[ iQ * ( (iCom+3)%nbNodes )]);
2418 if ( splitMethod.hasFacet( t012 ) && splitMethod.hasFacet( t023 ))
2420 facets.push_back( t012 );
2421 facets.push_back( t023 );
2422 for ( int iLast = iCom+4; iLast < iCom+nbNodes; ++iLast )
2423 facets.push_back( TTriangleFacet( nInd[ iQ * ( iCom )],
2424 nInd[ iQ * ((iLast-1)%nbNodes )],
2425 nInd[ iQ * ((iLast )%nbNodes )]));
2429 list< TTriangleFacet >::iterator facet = facets.begin();
2430 if ( facet == facets.end() )
2432 for ( ; facet != facets.end(); ++facet )
2434 if ( !volTool.IsFaceExternal( iF ))
2435 swap( facet->_n2, facet->_n3 );
2436 triangles.push_back( helper.AddFace( volNodes[ facet->_n1 ],
2437 volNodes[ facet->_n2 ],
2438 volNodes[ facet->_n3 ]));
2441 for ( size_t i = 0; i < triangles.size(); ++i )
2443 if ( !triangles[ i ]) continue;
2445 fSubMesh->AddElement( triangles[ i ]);
2446 newElems.Append( triangles[ i ]);
2448 ReplaceElemInGroups( face, triangles, GetMeshDS() );
2449 GetMeshDS()->RemoveFreeElement( face, fSubMesh, /*fromGroups=*/false );
2451 } // while a face based on facet nodes exists
2452 } // loop on volume faces to split them into triangles
2454 GetMeshDS()->RemoveFreeElement( elem, subMesh, /*fromGroups=*/false );
2456 if ( geomType == SMDSEntity_TriQuad_Hexa )
2458 // remove medium nodes that could become free
2459 for ( int i = 20; i < volTool.NbNodes(); ++i )
2460 if ( volNodes[i]->NbInverseElements() == 0 )
2461 GetMeshDS()->RemoveNode( volNodes[i] );
2463 } // loop on volumes to split
2465 myLastCreatedNodes = newNodes;
2466 myLastCreatedElems = newElems;
2469 //=======================================================================
2470 //function : GetHexaFacetsToSplit
2471 //purpose : For hexahedra that will be split into prisms, finds facets to
2472 // split into triangles. Only hexahedra adjacent to the one closest
2473 // to theFacetNormal.Location() are returned.
2474 //param [in,out] theHexas - the hexahedra
2475 //param [in] theFacetNormal - facet normal
2476 //param [out] theFacets - the hexahedra and found facet IDs
2477 //=======================================================================
2479 void SMESH_MeshEditor::GetHexaFacetsToSplit( TIDSortedElemSet& theHexas,
2480 const gp_Ax1& theFacetNormal,
2481 TFacetOfElem & theFacets)
2483 #define THIS_METHOD "SMESH_MeshEditor::GetHexaFacetsToSplit(): "
2485 // Find a hexa closest to the location of theFacetNormal
2487 const SMDS_MeshElement* startHex;
2489 // get SMDS_ElemIteratorPtr on theHexas
2490 typedef const SMDS_MeshElement* TValue;
2491 typedef TIDSortedElemSet::iterator TSetIterator;
2492 typedef SMDS::SimpleAccessor<TValue,TSetIterator> TAccesor;
2493 typedef SMDS_MeshElement::GeomFilter TFilter;
2494 typedef SMDS_SetIterator < TValue, TSetIterator, TAccesor, TFilter > TElemSetIter;
2495 SMDS_ElemIteratorPtr elemIt = SMDS_ElemIteratorPtr
2496 ( new TElemSetIter( theHexas.begin(),
2498 SMDS_MeshElement::GeomFilter( SMDSGeom_HEXA )));
2500 SMESH_ElementSearcher* searcher =
2501 SMESH_MeshAlgos::GetElementSearcher( *myMesh->GetMeshDS(), elemIt );
2503 startHex = searcher->FindClosestTo( theFacetNormal.Location(), SMDSAbs_Volume );
2508 throw SALOME_Exception( THIS_METHOD "startHex not found");
2511 // Select a facet of startHex by theFacetNormal
2513 SMDS_VolumeTool vTool( startHex );
2514 double norm[3], dot, maxDot = 0;
2516 for ( int iF = 0; iF < vTool.NbFaces(); ++iF )
2517 if ( vTool.GetFaceNormal( iF, norm[0], norm[1], norm[2] ))
2519 dot = Abs( theFacetNormal.Direction().Dot( gp_Dir( norm[0], norm[1], norm[2] )));
2527 throw SALOME_Exception( THIS_METHOD "facet of startHex not found");
2529 // Fill theFacets starting from facetID of startHex
2531 // facets used for seach of volumes adjacent to already treated ones
2532 typedef pair< TFacetOfElem::iterator, int > TElemFacets;
2533 typedef map< TVolumeFaceKey, TElemFacets > TFacetMap;
2534 TFacetMap facetsToCheck;
2536 set<const SMDS_MeshNode*> facetNodes;
2537 const SMDS_MeshElement* curHex;
2539 const bool allHex = ((int) theHexas.size() == myMesh->NbHexas() );
2543 // move in two directions from startHex via facetID
2544 for ( int is2nd = 0; is2nd < 2; ++is2nd )
2547 int curFacet = facetID;
2548 if ( is2nd ) // do not treat startHex twice
2550 vTool.Set( curHex );
2551 if ( vTool.IsFreeFace( curFacet, &curHex ))
2557 vTool.GetFaceNodes( curFacet, facetNodes );
2558 vTool.Set( curHex );
2559 curFacet = vTool.GetFaceIndex( facetNodes );
2564 // store a facet to split
2565 if ( curHex->GetGeomType() != SMDSGeom_HEXA )
2567 theFacets.insert( make_pair( curHex, -1 ));
2570 if ( !allHex && !theHexas.count( curHex ))
2573 pair< TFacetOfElem::iterator, bool > facetIt2isNew =
2574 theFacets.insert( make_pair( curHex, curFacet ));
2575 if ( !facetIt2isNew.second )
2578 // remember not-to-split facets in facetsToCheck
2579 int oppFacet = vTool.GetOppFaceIndexOfHex( curFacet );
2580 for ( int iF = 0; iF < vTool.NbFaces(); ++iF )
2582 if ( iF == curFacet && iF == oppFacet )
2584 TVolumeFaceKey facetKey ( vTool, iF );
2585 TElemFacets elemFacet( facetIt2isNew.first, iF );
2586 pair< TFacetMap::iterator, bool > it2isnew =
2587 facetsToCheck.insert( make_pair( facetKey, elemFacet ));
2588 if ( !it2isnew.second )
2589 facetsToCheck.erase( it2isnew.first ); // adjacent hex already checked
2591 // pass to a volume adjacent via oppFacet
2592 if ( vTool.IsFreeFace( oppFacet, &curHex ))
2598 // get a new curFacet
2599 vTool.GetFaceNodes( oppFacet, facetNodes );
2600 vTool.Set( curHex );
2601 curFacet = vTool.GetFaceIndex( facetNodes, /*hint=*/curFacet );
2604 } // move in two directions from startHex via facetID
2606 // Find a new startHex by facetsToCheck
2610 TFacetMap::iterator fIt = facetsToCheck.begin();
2611 while ( !startHex && fIt != facetsToCheck.end() )
2613 const TElemFacets& elemFacets = fIt->second;
2614 const SMDS_MeshElement* hex = elemFacets.first->first;
2615 int splitFacet = elemFacets.first->second;
2616 int lateralFacet = elemFacets.second;
2617 facetsToCheck.erase( fIt );
2618 fIt = facetsToCheck.begin();
2621 if ( vTool.IsFreeFace( lateralFacet, &curHex ) ||
2622 curHex->GetGeomType() != SMDSGeom_HEXA )
2624 if ( !allHex && !theHexas.count( curHex ))
2629 // find a facet of startHex to split
2631 set<const SMDS_MeshNode*> lateralNodes;
2632 vTool.GetFaceNodes( lateralFacet, lateralNodes );
2633 vTool.GetFaceNodes( splitFacet, facetNodes );
2634 int oppLateralFacet = vTool.GetOppFaceIndexOfHex( lateralFacet );
2635 vTool.Set( startHex );
2636 lateralFacet = vTool.GetFaceIndex( lateralNodes, oppLateralFacet );
2638 // look for a facet of startHex having common nodes with facetNodes
2639 // but not lateralFacet
2640 for ( int iF = 0; iF < vTool.NbFaces(); ++iF )
2642 if ( iF == lateralFacet )
2644 int nbCommonNodes = 0;
2645 const SMDS_MeshNode** nn = vTool.GetFaceNodes( iF );
2646 for ( int iN = 0, nbN = vTool.NbFaceNodes( iF ); iN < nbN; ++iN )
2647 nbCommonNodes += facetNodes.count( nn[ iN ]);
2649 if ( nbCommonNodes >= 2 )
2656 throw SALOME_Exception( THIS_METHOD "facet of a new startHex not found");
2658 } // while ( startHex )
2665 //================================================================================
2667 * \brief Selects nodes of several elements according to a given interlace
2668 * \param [in] srcNodes - nodes to select from
2669 * \param [out] tgtNodesVec - array of nodes of several elements to fill in
2670 * \param [in] interlace - indices of nodes for all elements
2671 * \param [in] nbElems - nb of elements
2672 * \param [in] nbNodes - nb of nodes in each element
2673 * \param [in] mesh - the mesh
2674 * \param [out] elemQueue - a list to push elements found by the selected nodes
2675 * \param [in] type - type of elements to look for
2677 //================================================================================
2679 void selectNodes( const vector< const SMDS_MeshNode* >& srcNodes,
2680 vector< const SMDS_MeshNode* >* tgtNodesVec,
2681 const int* interlace,
2684 SMESHDS_Mesh* mesh = 0,
2685 list< const SMDS_MeshElement* >* elemQueue=0,
2686 SMDSAbs_ElementType type=SMDSAbs_All)
2688 for ( int iE = 0; iE < nbElems; ++iE )
2690 vector< const SMDS_MeshNode* >& elemNodes = tgtNodesVec[iE];
2691 const int* select = & interlace[iE*nbNodes];
2692 elemNodes.resize( nbNodes );
2693 for ( int iN = 0; iN < nbNodes; ++iN )
2694 elemNodes[iN] = srcNodes[ select[ iN ]];
2696 const SMDS_MeshElement* e;
2698 for ( int iE = 0; iE < nbElems; ++iE )
2699 if (( e = mesh->FindElement( tgtNodesVec[iE], type, /*noMedium=*/false)))
2700 elemQueue->push_back( e );
2704 //=======================================================================
2706 * Split bi-quadratic elements into linear ones without creation of additional nodes
2707 * - bi-quadratic triangle will be split into 3 linear quadrangles;
2708 * - bi-quadratic quadrangle will be split into 4 linear quadrangles;
2709 * - tri-quadratic hexahedron will be split into 8 linear hexahedra;
2710 * Quadratic elements of lower dimension adjacent to the split bi-quadratic element
2711 * will be split in order to keep the mesh conformal.
2712 * \param elems - elements to split
2714 //=======================================================================
2716 void SMESH_MeshEditor::SplitBiQuadraticIntoLinear(TIDSortedElemSet& theElems)
2718 vector< const SMDS_MeshNode* > elemNodes(27), subNodes[12], splitNodes[8];
2719 vector<const SMDS_MeshElement* > splitElems;
2720 list< const SMDS_MeshElement* > elemQueue;
2721 list< const SMDS_MeshElement* >::iterator elemIt;
2723 SMESHDS_Mesh * mesh = GetMeshDS();
2724 ElemFeatures *elemType, hexaType(SMDSAbs_Volume), quadType(SMDSAbs_Face), segType(SMDSAbs_Edge);
2725 int nbElems, nbNodes;
2727 TIDSortedElemSet::iterator elemSetIt = theElems.begin();
2728 for ( ; elemSetIt != theElems.end(); ++elemSetIt )
2731 elemQueue.push_back( *elemSetIt );
2732 for ( elemIt = elemQueue.begin(); elemIt != elemQueue.end(); ++elemIt )
2734 const SMDS_MeshElement* elem = *elemIt;
2735 switch( elem->GetEntityType() )
2737 case SMDSEntity_TriQuad_Hexa: // HEX27
2739 elemNodes.assign( elem->begin_nodes(), elem->end_nodes() );
2740 nbElems = nbNodes = 8;
2741 elemType = & hexaType;
2743 // get nodes for new elements
2744 static int vInd[8][8] = {{ 0,8,20,11, 16,21,26,24 },
2745 { 1,9,20,8, 17,22,26,21 },
2746 { 2,10,20,9, 18,23,26,22 },
2747 { 3,11,20,10, 19,24,26,23 },
2748 { 16,21,26,24, 4,12,25,15 },
2749 { 17,22,26,21, 5,13,25,12 },
2750 { 18,23,26,22, 6,14,25,13 },
2751 { 19,24,26,23, 7,15,25,14 }};
2752 selectNodes( elemNodes, & splitNodes[0], &vInd[0][0], nbElems, nbNodes );
2754 // add boundary faces to elemQueue
2755 static int fInd[6][9] = {{ 0,1,2,3, 8,9,10,11, 20 },
2756 { 4,5,6,7, 12,13,14,15, 25 },
2757 { 0,1,5,4, 8,17,12,16, 21 },
2758 { 1,2,6,5, 9,18,13,17, 22 },
2759 { 2,3,7,6, 10,19,14,18, 23 },
2760 { 3,0,4,7, 11,16,15,19, 24 }};
2761 selectNodes( elemNodes, & subNodes[0], &fInd[0][0], 6,9, mesh, &elemQueue, SMDSAbs_Face );
2763 // add boundary segments to elemQueue
2764 static int eInd[12][3] = {{ 0,1,8 }, { 1,2,9 }, { 2,3,10 }, { 3,0,11 },
2765 { 4,5,12}, { 5,6,13}, { 6,7,14 }, { 7,4,15 },
2766 { 0,4,16}, { 1,5,17}, { 2,6,18 }, { 3,7,19 }};
2767 selectNodes( elemNodes, & subNodes[0], &eInd[0][0], 12,3, mesh, &elemQueue, SMDSAbs_Edge );
2770 case SMDSEntity_BiQuad_Triangle: // TRIA7
2772 elemNodes.assign( elem->begin_nodes(), elem->end_nodes() );
2775 elemType = & quadType;
2777 // get nodes for new elements
2778 static int fInd[3][4] = {{ 0,3,6,5 }, { 1,4,6,3 }, { 2,5,6,4 }};
2779 selectNodes( elemNodes, & splitNodes[0], &fInd[0][0], nbElems, nbNodes );
2781 // add boundary segments to elemQueue
2782 static int eInd[3][3] = {{ 0,1,3 }, { 1,2,4 }, { 2,0,5 }};
2783 selectNodes( elemNodes, & subNodes[0], &eInd[0][0], 3,3, mesh, &elemQueue, SMDSAbs_Edge );
2786 case SMDSEntity_BiQuad_Quadrangle: // QUAD9
2788 elemNodes.assign( elem->begin_nodes(), elem->end_nodes() );
2791 elemType = & quadType;
2793 // get nodes for new elements
2794 static int fInd[4][4] = {{ 0,4,8,7 }, { 1,5,8,4 }, { 2,6,8,5 }, { 3,7,8,6 }};
2795 selectNodes( elemNodes, & splitNodes[0], &fInd[0][0], nbElems, nbNodes );
2797 // add boundary segments to elemQueue
2798 static int eInd[4][3] = {{ 0,1,4 }, { 1,2,5 }, { 2,3,6 }, { 3,0,7 }};
2799 selectNodes( elemNodes, & subNodes[0], &eInd[0][0], 4,3, mesh, &elemQueue, SMDSAbs_Edge );
2802 case SMDSEntity_Quad_Edge:
2804 if ( elemIt == elemQueue.begin() )
2805 continue; // an elem is in theElems
2806 elemNodes.assign( elem->begin_nodes(), elem->end_nodes() );
2809 elemType = & segType;
2811 // get nodes for new elements
2812 static int eInd[2][2] = {{ 0,2 }, { 2,1 }};
2813 selectNodes( elemNodes, & splitNodes[0], &eInd[0][0], nbElems, nbNodes );
2817 } // switch( elem->GetEntityType() )
2819 // Create new elements
2821 SMESHDS_SubMesh* subMesh = mesh->MeshElements( elem->getshapeId() );
2825 //elemType->SetID( elem->GetID() ); // create an elem with the same ID as a removed one
2826 mesh->RemoveFreeElement( elem, subMesh, /*fromGroups=*/false );
2827 //splitElems.push_back( AddElement( splitNodes[ 0 ], *elemType ));
2828 //elemType->SetID( -1 );
2830 for ( int iE = 0; iE < nbElems; ++iE )
2831 splitElems.push_back( AddElement( splitNodes[ iE ], *elemType ));
2834 ReplaceElemInGroups( elem, splitElems, mesh );
2837 for ( size_t i = 0; i < splitElems.size(); ++i )
2838 subMesh->AddElement( splitElems[i] );
2843 //=======================================================================
2844 //function : AddToSameGroups
2845 //purpose : add elemToAdd to the groups the elemInGroups belongs to
2846 //=======================================================================
2848 void SMESH_MeshEditor::AddToSameGroups (const SMDS_MeshElement* elemToAdd,
2849 const SMDS_MeshElement* elemInGroups,
2850 SMESHDS_Mesh * aMesh)
2852 const set<SMESHDS_GroupBase*>& groups = aMesh->GetGroups();
2853 if (!groups.empty()) {
2854 set<SMESHDS_GroupBase*>::const_iterator grIt = groups.begin();
2855 for ( ; grIt != groups.end(); grIt++ ) {
2856 SMESHDS_Group* group = dynamic_cast<SMESHDS_Group*>( *grIt );
2857 if ( group && group->Contains( elemInGroups ))
2858 group->SMDSGroup().Add( elemToAdd );
2864 //=======================================================================
2865 //function : RemoveElemFromGroups
2866 //purpose : Remove removeelem to the groups the elemInGroups belongs to
2867 //=======================================================================
2868 void SMESH_MeshEditor::RemoveElemFromGroups (const SMDS_MeshElement* removeelem,
2869 SMESHDS_Mesh * aMesh)
2871 const set<SMESHDS_GroupBase*>& groups = aMesh->GetGroups();
2872 if (!groups.empty())
2874 set<SMESHDS_GroupBase*>::const_iterator GrIt = groups.begin();
2875 for (; GrIt != groups.end(); GrIt++)
2877 SMESHDS_Group* grp = dynamic_cast<SMESHDS_Group*>(*GrIt);
2878 if (!grp || grp->IsEmpty()) continue;
2879 grp->SMDSGroup().Remove(removeelem);
2884 //================================================================================
2886 * \brief Replace elemToRm by elemToAdd in the all groups
2888 //================================================================================
2890 void SMESH_MeshEditor::ReplaceElemInGroups (const SMDS_MeshElement* elemToRm,
2891 const SMDS_MeshElement* elemToAdd,
2892 SMESHDS_Mesh * aMesh)
2894 const set<SMESHDS_GroupBase*>& groups = aMesh->GetGroups();
2895 if (!groups.empty()) {
2896 set<SMESHDS_GroupBase*>::const_iterator grIt = groups.begin();
2897 for ( ; grIt != groups.end(); grIt++ ) {
2898 SMESHDS_Group* group = dynamic_cast<SMESHDS_Group*>( *grIt );
2899 if ( group && group->SMDSGroup().Remove( elemToRm ) && elemToAdd )
2900 group->SMDSGroup().Add( elemToAdd );
2905 //================================================================================
2907 * \brief Replace elemToRm by elemToAdd in the all groups
2909 //================================================================================
2911 void SMESH_MeshEditor::ReplaceElemInGroups (const SMDS_MeshElement* elemToRm,
2912 const vector<const SMDS_MeshElement*>& elemToAdd,
2913 SMESHDS_Mesh * aMesh)
2915 const set<SMESHDS_GroupBase*>& groups = aMesh->GetGroups();
2916 if (!groups.empty())
2918 set<SMESHDS_GroupBase*>::const_iterator grIt = groups.begin();
2919 for ( ; grIt != groups.end(); grIt++ ) {
2920 SMESHDS_Group* group = dynamic_cast<SMESHDS_Group*>( *grIt );
2921 if ( group && group->SMDSGroup().Remove( elemToRm ) )
2922 for ( size_t i = 0; i < elemToAdd.size(); ++i )
2923 group->SMDSGroup().Add( elemToAdd[ i ] );
2928 //=======================================================================
2929 //function : QuadToTri
2930 //purpose : Cut quadrangles into triangles.
2931 // theCrit is used to select a diagonal to cut
2932 //=======================================================================
2934 bool SMESH_MeshEditor::QuadToTri (TIDSortedElemSet & theElems,
2935 const bool the13Diag)
2937 myLastCreatedElems.Clear();
2938 myLastCreatedNodes.Clear();
2940 SMESHDS_Mesh * aMesh = GetMeshDS();
2942 Handle(Geom_Surface) surface;
2943 SMESH_MesherHelper helper( *GetMesh() );
2945 TIDSortedElemSet::iterator itElem;
2946 for ( itElem = theElems.begin(); itElem != theElems.end(); itElem++ )
2948 const SMDS_MeshElement* elem = *itElem;
2949 if ( !elem || elem->GetGeomType() != SMDSGeom_QUADRANGLE )
2952 if ( elem->NbNodes() == 4 ) {
2953 // retrieve element nodes
2954 const SMDS_MeshNode* aNodes [4];
2955 SMDS_ElemIteratorPtr itN = elem->nodesIterator();
2957 while ( itN->more() )
2958 aNodes[ i++ ] = static_cast<const SMDS_MeshNode*>( itN->next() );
2960 int aShapeId = FindShape( elem );
2961 const SMDS_MeshElement* newElem1 = 0;
2962 const SMDS_MeshElement* newElem2 = 0;
2964 newElem1 = aMesh->AddFace( aNodes[2], aNodes[0], aNodes[1] );
2965 newElem2 = aMesh->AddFace( aNodes[2], aNodes[3], aNodes[0] );
2968 newElem1 = aMesh->AddFace( aNodes[3], aNodes[0], aNodes[1] );
2969 newElem2 = aMesh->AddFace( aNodes[3], aNodes[1], aNodes[2] );
2971 myLastCreatedElems.Append(newElem1);
2972 myLastCreatedElems.Append(newElem2);
2973 // put a new triangle on the same shape and add to the same groups
2976 aMesh->SetMeshElementOnShape( newElem1, aShapeId );
2977 aMesh->SetMeshElementOnShape( newElem2, aShapeId );
2979 AddToSameGroups( newElem1, elem, aMesh );
2980 AddToSameGroups( newElem2, elem, aMesh );
2981 aMesh->RemoveElement( elem );
2984 // Quadratic quadrangle
2986 else if ( elem->NbNodes() >= 8 )
2988 // get surface elem is on
2989 int aShapeId = FindShape( elem );
2990 if ( aShapeId != helper.GetSubShapeID() ) {