+ SMDS_MeshEdge* edge = 0;
+ if (myCreateQuadratic) {
+ const SMDS_MeshNode* n12 = GetMediumNode(n1,n2,force3d);
+ if(id)
+ edge = meshDS->AddEdgeWithID(n1, n2, n12, id);
+ else
+ edge = meshDS->AddEdge(n1, n2, n12);
+ }
+ else {
+ if(id)
+ edge = meshDS->AddEdgeWithID(n1, n2, id);
+ else
+ edge = meshDS->AddEdge(n1, n2);
+ }
+
+ if ( mySetElemOnShape && myShapeID > 0 )
+ meshDS->SetMeshElementOnShape( edge, myShapeID );
+
+ return edge;
+}
+
+//=======================================================================
+//function : AddFace
+//purpose : Creates quadratic or linear triangle
+//=======================================================================
+
+SMDS_MeshFace* SMESH_MesherHelper::AddFace(const SMDS_MeshNode* n1,
+ const SMDS_MeshNode* n2,
+ const SMDS_MeshNode* n3,
+ const int id,
+ const bool force3d)
+{
+ SMESHDS_Mesh * meshDS = GetMeshDS();
+ SMDS_MeshFace* elem = 0;
+
+ if( n1==n2 || n2==n3 || n3==n1 )
+ return elem;
+
+ if(!myCreateQuadratic) {
+ if(id)
+ elem = meshDS->AddFaceWithID(n1, n2, n3, id);
+ else
+ elem = meshDS->AddFace(n1, n2, n3);
+ }
+ else {
+ const SMDS_MeshNode* n12 = GetMediumNode(n1,n2,force3d);
+ const SMDS_MeshNode* n23 = GetMediumNode(n2,n3,force3d);
+ const SMDS_MeshNode* n31 = GetMediumNode(n3,n1,force3d);
+ if(myCreateBiQuadratic)
+ {
+ const SMDS_MeshNode* nCenter = GetCentralNode(n1, n2, n3, n12, n23, n31, force3d);
+ if(id)
+ elem = meshDS->AddFaceWithID(n1, n2, n3, n12, n23, n31, nCenter, id);
+ else
+ elem = meshDS->AddFace(n1, n2, n3, n12, n23, n31, nCenter);
+ }
+ else
+ {
+ if(id)
+ elem = meshDS->AddFaceWithID(n1, n2, n3, n12, n23, n31, id);
+ else
+ elem = meshDS->AddFace(n1, n2, n3, n12, n23, n31);
+ }
+ }
+ if ( mySetElemOnShape && myShapeID > 0 )
+ meshDS->SetMeshElementOnShape( elem, myShapeID );
+
+ return elem;
+}
+
+//=======================================================================
+//function : AddFace
+//purpose : Creates bi-quadratic, quadratic or linear quadrangle
+//=======================================================================
+
+SMDS_MeshFace* SMESH_MesherHelper::AddFace(const SMDS_MeshNode* n1,
+ const SMDS_MeshNode* n2,
+ const SMDS_MeshNode* n3,
+ const SMDS_MeshNode* n4,
+ const int id,
+ const bool force3d)
+{
+ SMESHDS_Mesh * meshDS = GetMeshDS();
+ SMDS_MeshFace* elem = 0;
+
+ if( n1==n2 ) {
+ return AddFace(n1,n3,n4,id,force3d);
+ }
+ if( n1==n3 ) {
+ return AddFace(n1,n2,n4,id,force3d);
+ }
+ if( n1==n4 ) {
+ return AddFace(n1,n2,n3,id,force3d);
+ }
+ if( n2==n3 ) {
+ return AddFace(n1,n2,n4,id,force3d);
+ }
+ if( n2==n4 ) {
+ return AddFace(n1,n2,n3,id,force3d);
+ }
+ if( n3==n4 ) {
+ return AddFace(n1,n2,n3,id,force3d);
+ }
+
+ if(!myCreateQuadratic) {
+ if(id)
+ elem = meshDS->AddFaceWithID(n1, n2, n3, n4, id);
+ else
+ elem = meshDS->AddFace(n1, n2, n3, n4);
+ }
+ else {
+ const SMDS_MeshNode* n12 = GetMediumNode(n1,n2,force3d);
+ const SMDS_MeshNode* n23 = GetMediumNode(n2,n3,force3d);
+ const SMDS_MeshNode* n34 = GetMediumNode(n3,n4,force3d);
+ const SMDS_MeshNode* n41 = GetMediumNode(n4,n1,force3d);
+ if(myCreateBiQuadratic)
+ {
+ const SMDS_MeshNode* nCenter = GetCentralNode(n1, n2, n3, n4, n12, n23, n34, n41, force3d);
+ if(id)
+ elem = meshDS->AddFaceWithID(n1, n2, n3, n4, n12, n23, n34, n41, nCenter, id);
+ else
+ elem = meshDS->AddFace(n1, n2, n3, n4, n12, n23, n34, n41, nCenter);
+ }
+ else
+ {
+ if(id)
+ elem = meshDS->AddFaceWithID(n1, n2, n3, n4, n12, n23, n34, n41, id);
+ else
+ elem = meshDS->AddFace(n1, n2, n3, n4, n12, n23, n34, n41);
+ }
+ }
+ if ( mySetElemOnShape && myShapeID > 0 )
+ meshDS->SetMeshElementOnShape( elem, myShapeID );
+
+ return elem;
+}
+
+//=======================================================================
+//function : AddPolygonalFace
+//purpose : Creates polygon, with additional nodes in quadratic mesh
+//=======================================================================
+
+SMDS_MeshFace* SMESH_MesherHelper::AddPolygonalFace (const vector<const SMDS_MeshNode*>& nodes,
+ const int id,
+ const bool force3d)
+{
+ SMESHDS_Mesh * meshDS = GetMeshDS();
+ SMDS_MeshFace* elem = 0;
+
+ if(!myCreateQuadratic) {
+ if(id)
+ elem = meshDS->AddPolygonalFaceWithID(nodes, id);
+ else
+ elem = meshDS->AddPolygonalFace(nodes);
+ }
+ else {
+ vector<const SMDS_MeshNode*> newNodes;
+ for ( int i = 0; i < nodes.size(); ++i )
+ {
+ const SMDS_MeshNode* n1 = nodes[i];
+ const SMDS_MeshNode* n2 = nodes[(i+1)%nodes.size()];
+ const SMDS_MeshNode* n12 = GetMediumNode(n1,n2,force3d);
+ newNodes.push_back( n1 );
+ newNodes.push_back( n12 );
+ }
+ if(id)
+ elem = meshDS->AddPolygonalFaceWithID(newNodes, id);
+ else
+ elem = meshDS->AddPolygonalFace(newNodes);
+ }
+ if ( mySetElemOnShape && myShapeID > 0 )
+ meshDS->SetMeshElementOnShape( elem, myShapeID );
+
+ return elem;
+}
+
+//=======================================================================
+//function : AddVolume
+//purpose : Creates quadratic or linear prism
+//=======================================================================
+
+SMDS_MeshVolume* SMESH_MesherHelper::AddVolume(const SMDS_MeshNode* n1,
+ const SMDS_MeshNode* n2,
+ const SMDS_MeshNode* n3,
+ const SMDS_MeshNode* n4,
+ const SMDS_MeshNode* n5,
+ const SMDS_MeshNode* n6,
+ const int id,
+ const bool force3d)
+{
+ SMESHDS_Mesh * meshDS = GetMeshDS();
+ SMDS_MeshVolume* elem = 0;
+ if(!myCreateQuadratic) {
+ if(id)
+ elem = meshDS->AddVolumeWithID(n1, n2, n3, n4, n5, n6, id);
+ else
+ elem = meshDS->AddVolume(n1, n2, n3, n4, n5, n6);
+ }
+ else {
+ const SMDS_MeshNode* n12 = GetMediumNode(n1,n2,force3d);
+ const SMDS_MeshNode* n23 = GetMediumNode(n2,n3,force3d);
+ const SMDS_MeshNode* n31 = GetMediumNode(n3,n1,force3d);
+
+ const SMDS_MeshNode* n45 = GetMediumNode(n4,n5,force3d);
+ const SMDS_MeshNode* n56 = GetMediumNode(n5,n6,force3d);
+ const SMDS_MeshNode* n64 = GetMediumNode(n6,n4,force3d);
+
+ const SMDS_MeshNode* n14 = GetMediumNode(n1,n4,force3d);
+ const SMDS_MeshNode* n25 = GetMediumNode(n2,n5,force3d);
+ const SMDS_MeshNode* n36 = GetMediumNode(n3,n6,force3d);
+
+ if(id)
+ elem = meshDS->AddVolumeWithID(n1, n2, n3, n4, n5, n6,
+ n12, n23, n31, n45, n56, n64, n14, n25, n36, id);
+ else
+ elem = meshDS->AddVolume(n1, n2, n3, n4, n5, n6,
+ n12, n23, n31, n45, n56, n64, n14, n25, n36);
+ }
+ if ( mySetElemOnShape && myShapeID > 0 )
+ meshDS->SetMeshElementOnShape( elem, myShapeID );
+
+ return elem;
+}
+
+//=======================================================================
+//function : AddVolume
+//purpose : Creates quadratic or linear tetrahedron
+//=======================================================================
+
+SMDS_MeshVolume* SMESH_MesherHelper::AddVolume(const SMDS_MeshNode* n1,
+ const SMDS_MeshNode* n2,
+ const SMDS_MeshNode* n3,
+ const SMDS_MeshNode* n4,
+ const int id,
+ const bool force3d)
+{
+ SMESHDS_Mesh * meshDS = GetMeshDS();
+ SMDS_MeshVolume* elem = 0;
+ if(!myCreateQuadratic) {
+ if(id)
+ elem = meshDS->AddVolumeWithID(n1, n2, n3, n4, id);
+ else
+ elem = meshDS->AddVolume(n1, n2, n3, n4);
+ }
+ else {
+ const SMDS_MeshNode* n12 = GetMediumNode(n1,n2,force3d);
+ const SMDS_MeshNode* n23 = GetMediumNode(n2,n3,force3d);
+ const SMDS_MeshNode* n31 = GetMediumNode(n3,n1,force3d);
+
+ const SMDS_MeshNode* n14 = GetMediumNode(n1,n4,force3d);
+ const SMDS_MeshNode* n24 = GetMediumNode(n2,n4,force3d);
+ const SMDS_MeshNode* n34 = GetMediumNode(n3,n4,force3d);
+
+ if(id)
+ elem = meshDS->AddVolumeWithID(n1, n2, n3, n4, n12, n23, n31, n14, n24, n34, id);
+ else
+ elem = meshDS->AddVolume(n1, n2, n3, n4, n12, n23, n31, n14, n24, n34);
+ }
+ if ( mySetElemOnShape && myShapeID > 0 )
+ meshDS->SetMeshElementOnShape( elem, myShapeID );
+
+ return elem;
+}
+
+//=======================================================================
+//function : AddVolume
+//purpose : Creates quadratic or linear pyramid
+//=======================================================================
+
+SMDS_MeshVolume* SMESH_MesherHelper::AddVolume(const SMDS_MeshNode* n1,
+ const SMDS_MeshNode* n2,
+ const SMDS_MeshNode* n3,
+ const SMDS_MeshNode* n4,
+ const SMDS_MeshNode* n5,
+ const int id,
+ const bool force3d)
+{
+ SMDS_MeshVolume* elem = 0;
+ if(!myCreateQuadratic) {
+ if(id)
+ elem = GetMeshDS()->AddVolumeWithID(n1, n2, n3, n4, n5, id);
+ else
+ elem = GetMeshDS()->AddVolume(n1, n2, n3, n4, n5);
+ }
+ else {
+ const SMDS_MeshNode* n12 = GetMediumNode(n1,n2,force3d);
+ const SMDS_MeshNode* n23 = GetMediumNode(n2,n3,force3d);
+ const SMDS_MeshNode* n34 = GetMediumNode(n3,n4,force3d);
+ const SMDS_MeshNode* n41 = GetMediumNode(n4,n1,force3d);
+
+ const SMDS_MeshNode* n15 = GetMediumNode(n1,n5,force3d);
+ const SMDS_MeshNode* n25 = GetMediumNode(n2,n5,force3d);
+ const SMDS_MeshNode* n35 = GetMediumNode(n3,n5,force3d);
+ const SMDS_MeshNode* n45 = GetMediumNode(n4,n5,force3d);
+
+ if(id)
+ elem = GetMeshDS()->AddVolumeWithID ( n1, n2, n3, n4, n5,
+ n12, n23, n34, n41,
+ n15, n25, n35, n45,
+ id);
+ else
+ elem = GetMeshDS()->AddVolume( n1, n2, n3, n4, n5,
+ n12, n23, n34, n41,
+ n15, n25, n35, n45);
+ }
+ if ( mySetElemOnShape && myShapeID > 0 )
+ GetMeshDS()->SetMeshElementOnShape( elem, myShapeID );
+
+ return elem;
+}
+
+//=======================================================================
+//function : AddVolume
+//purpose : Creates bi-quadratic, quadratic or linear hexahedron
+//=======================================================================
+
+SMDS_MeshVolume* SMESH_MesherHelper::AddVolume(const SMDS_MeshNode* n1,
+ const SMDS_MeshNode* n2,
+ const SMDS_MeshNode* n3,
+ const SMDS_MeshNode* n4,
+ const SMDS_MeshNode* n5,
+ const SMDS_MeshNode* n6,
+ const SMDS_MeshNode* n7,
+ const SMDS_MeshNode* n8,
+ const int id,
+ const bool force3d)
+{
+ SMESHDS_Mesh * meshDS = GetMeshDS();
+ SMDS_MeshVolume* elem = 0;
+ if(!myCreateQuadratic) {
+ if(id)
+ elem = meshDS->AddVolumeWithID(n1, n2, n3, n4, n5, n6, n7, n8, id);
+ else
+ elem = meshDS->AddVolume(n1, n2, n3, n4, n5, n6, n7, n8);
+ }
+ else {
+ const SMDS_MeshNode* n12 = GetMediumNode(n1,n2,force3d);
+ const SMDS_MeshNode* n23 = GetMediumNode(n2,n3,force3d);
+ const SMDS_MeshNode* n34 = GetMediumNode(n3,n4,force3d);
+ const SMDS_MeshNode* n41 = GetMediumNode(n4,n1,force3d);
+
+ const SMDS_MeshNode* n56 = GetMediumNode(n5,n6,force3d);
+ const SMDS_MeshNode* n67 = GetMediumNode(n6,n7,force3d);
+ const SMDS_MeshNode* n78 = GetMediumNode(n7,n8,force3d);
+ const SMDS_MeshNode* n85 = GetMediumNode(n8,n5,force3d);
+
+ const SMDS_MeshNode* n15 = GetMediumNode(n1,n5,force3d);
+ const SMDS_MeshNode* n26 = GetMediumNode(n2,n6,force3d);
+ const SMDS_MeshNode* n37 = GetMediumNode(n3,n7,force3d);
+ const SMDS_MeshNode* n48 = GetMediumNode(n4,n8,force3d);
+ if(myCreateBiQuadratic)
+ {
+ const SMDS_MeshNode* n1234 = GetCentralNode(n1,n2,n3,n4,n12,n23,n34,n41,force3d);
+ const SMDS_MeshNode* n1256 = GetCentralNode(n1,n2,n5,n6,n12,n26,n56,n15,force3d);
+ const SMDS_MeshNode* n2367 = GetCentralNode(n2,n3,n6,n7,n23,n37,n67,n26,force3d);
+ const SMDS_MeshNode* n3478 = GetCentralNode(n3,n4,n7,n8,n34,n48,n78,n37,force3d);
+ const SMDS_MeshNode* n1458 = GetCentralNode(n1,n4,n5,n8,n41,n48,n15,n85,force3d);
+ const SMDS_MeshNode* n5678 = GetCentralNode(n5,n6,n7,n8,n56,n67,n78,n85,force3d);
+
+ vector<gp_XYZ> pointsOnShapes( SMESH_Block::ID_Shell );
+
+ pointsOnShapes[ SMESH_Block::ID_V000 ] = SMESH_TNodeXYZ( n4 );
+ pointsOnShapes[ SMESH_Block::ID_V100 ] = SMESH_TNodeXYZ( n8 );
+ pointsOnShapes[ SMESH_Block::ID_V010 ] = SMESH_TNodeXYZ( n3 );
+ pointsOnShapes[ SMESH_Block::ID_V110 ] = SMESH_TNodeXYZ( n7 );
+ pointsOnShapes[ SMESH_Block::ID_V001 ] = SMESH_TNodeXYZ( n1 );
+ pointsOnShapes[ SMESH_Block::ID_V101 ] = SMESH_TNodeXYZ( n5 );
+ pointsOnShapes[ SMESH_Block::ID_V011 ] = SMESH_TNodeXYZ( n2 );
+ pointsOnShapes[ SMESH_Block::ID_V111 ] = SMESH_TNodeXYZ( n6 );
+
+ pointsOnShapes[ SMESH_Block::ID_Ex00 ] = SMESH_TNodeXYZ( n48 );
+ pointsOnShapes[ SMESH_Block::ID_Ex10 ] = SMESH_TNodeXYZ( n37 );
+ pointsOnShapes[ SMESH_Block::ID_E0y0 ] = SMESH_TNodeXYZ( n15 );
+ pointsOnShapes[ SMESH_Block::ID_E1y0 ] = SMESH_TNodeXYZ( n26 );
+ pointsOnShapes[ SMESH_Block::ID_Ex01 ] = SMESH_TNodeXYZ( n34 );
+ pointsOnShapes[ SMESH_Block::ID_Ex11 ] = SMESH_TNodeXYZ( n78 );
+ pointsOnShapes[ SMESH_Block::ID_E0y1 ] = SMESH_TNodeXYZ( n12 );
+ pointsOnShapes[ SMESH_Block::ID_E1y1 ] = SMESH_TNodeXYZ( n56 );
+ pointsOnShapes[ SMESH_Block::ID_E00z ] = SMESH_TNodeXYZ( n41 );
+ pointsOnShapes[ SMESH_Block::ID_E10z ] = SMESH_TNodeXYZ( n85 );
+ pointsOnShapes[ SMESH_Block::ID_E01z ] = SMESH_TNodeXYZ( n23 );
+ pointsOnShapes[ SMESH_Block::ID_E11z ] = SMESH_TNodeXYZ( n67 );
+
+ pointsOnShapes[ SMESH_Block::ID_Fxy0 ] = SMESH_TNodeXYZ( n3478 );
+ pointsOnShapes[ SMESH_Block::ID_Fxy1 ] = SMESH_TNodeXYZ( n1256 );
+ pointsOnShapes[ SMESH_Block::ID_Fx0z ] = SMESH_TNodeXYZ( n1458 );
+ pointsOnShapes[ SMESH_Block::ID_Fx1z ] = SMESH_TNodeXYZ( n2367 );
+ pointsOnShapes[ SMESH_Block::ID_F0yz ] = SMESH_TNodeXYZ( n1234 );
+ pointsOnShapes[ SMESH_Block::ID_F1yz ] = SMESH_TNodeXYZ( n5678 );
+
+ gp_XYZ centerCube(0.5, 0.5, 0.5);
+ gp_XYZ nCenterElem;
+ SMESH_Block::ShellPoint( centerCube, pointsOnShapes, nCenterElem );
+ const SMDS_MeshNode* nCenter =
+ meshDS->AddNode( nCenterElem.X(), nCenterElem.Y(), nCenterElem.Z() );
+ meshDS->SetNodeInVolume( nCenter, myShapeID );
+
+ if(id)
+ elem = meshDS->AddVolumeWithID(n1, n2, n3, n4, n5, n6, n7, n8,
+ n12, n23, n34, n41, n56, n67,
+ n78, n85, n15, n26, n37, n48,
+ n1234, n1256, n2367, n3478, n1458, n5678, nCenter, id);
+ else
+ elem = meshDS->AddVolume(n1, n2, n3, n4, n5, n6, n7, n8,
+ n12, n23, n34, n41, n56, n67,
+ n78, n85, n15, n26, n37, n48,
+ n1234, n1256, n2367, n3478, n1458, n5678, nCenter);
+ }
+ else
+ {
+ if(id)
+ elem = meshDS->AddVolumeWithID(n1, n2, n3, n4, n5, n6, n7, n8,
+ n12, n23, n34, n41, n56, n67,
+ n78, n85, n15, n26, n37, n48, id);
+ else
+ elem = meshDS->AddVolume(n1, n2, n3, n4, n5, n6, n7, n8,
+ n12, n23, n34, n41, n56, n67,
+ n78, n85, n15, n26, n37, n48);
+ }
+ }
+ if ( mySetElemOnShape && myShapeID > 0 )
+ meshDS->SetMeshElementOnShape( elem, myShapeID );
+
+ return elem;
+}
+
+//=======================================================================
+//function : AddVolume
+//purpose : Creates LINEAR!!!!!!!!! octahedron
+//=======================================================================
+
+SMDS_MeshVolume* SMESH_MesherHelper::AddVolume(const SMDS_MeshNode* n1,
+ const SMDS_MeshNode* n2,
+ const SMDS_MeshNode* n3,
+ const SMDS_MeshNode* n4,
+ const SMDS_MeshNode* n5,
+ const SMDS_MeshNode* n6,
+ const SMDS_MeshNode* n7,
+ const SMDS_MeshNode* n8,
+ const SMDS_MeshNode* n9,
+ const SMDS_MeshNode* n10,
+ const SMDS_MeshNode* n11,
+ const SMDS_MeshNode* n12,
+ const int id,
+ bool force3d)
+{
+ SMESHDS_Mesh * meshDS = GetMeshDS();
+ SMDS_MeshVolume* elem = 0;
+ if(id)
+ elem = meshDS->AddVolumeWithID(n1,n2,n3,n4,n5,n6,n7,n8,n9,n10,n11,n12,id);
+ else
+ elem = meshDS->AddVolume(n1,n2,n3,n4,n5,n6,n7,n8,n9,n10,n11,n12);
+ if ( mySetElemOnShape && myShapeID > 0 )
+ meshDS->SetMeshElementOnShape( elem, myShapeID );
+ return elem;
+}
+
+//=======================================================================
+//function : AddPolyhedralVolume
+//purpose : Creates polyhedron. In quadratic mesh, adds medium nodes
+//=======================================================================
+
+SMDS_MeshVolume*
+SMESH_MesherHelper::AddPolyhedralVolume (const std::vector<const SMDS_MeshNode*>& nodes,
+ const std::vector<int>& quantities,
+ const int id,
+ const bool force3d)
+{
+ SMESHDS_Mesh * meshDS = GetMeshDS();
+ SMDS_MeshVolume* elem = 0;
+ if(!myCreateQuadratic)
+ {
+ if(id)
+ elem = meshDS->AddPolyhedralVolumeWithID(nodes, quantities, id);
+ else
+ elem = meshDS->AddPolyhedralVolume(nodes, quantities);
+ }
+ else
+ {
+ vector<const SMDS_MeshNode*> newNodes;
+ vector<int> newQuantities;
+ for ( int iFace=0, iN=0; iFace < quantities.size(); ++iFace)
+ {
+ int nbNodesInFace = quantities[iFace];
+ newQuantities.push_back(0);
+ for ( int i = 0; i < nbNodesInFace; ++i )
+ {
+ const SMDS_MeshNode* n1 = nodes[ iN + i ];
+ newNodes.push_back( n1 );
+ newQuantities.back()++;
+
+ const SMDS_MeshNode* n2 = nodes[ iN + ( i+1==nbNodesInFace ? 0 : i+1 )];
+// if ( n1->GetPosition()->GetTypeOfPosition() != SMDS_TOP_3DSPACE &&
+// n2->GetPosition()->GetTypeOfPosition() != SMDS_TOP_3DSPACE )
+ {
+ const SMDS_MeshNode* n12 = GetMediumNode(n1,n2,force3d);
+ newNodes.push_back( n12 );
+ newQuantities.back()++;
+ }
+ }
+ iN += nbNodesInFace;
+ }
+ if(id)
+ elem = meshDS->AddPolyhedralVolumeWithID( newNodes, newQuantities, id );
+ else
+ elem = meshDS->AddPolyhedralVolume( newNodes, newQuantities );
+ }
+ if ( mySetElemOnShape && myShapeID > 0 )
+ meshDS->SetMeshElementOnShape( elem, myShapeID );
+
+ return elem;
+}
+
+namespace
+{
+ //================================================================================
+ /*!
+ * \brief Check if a node belongs to any face of sub-mesh
+ */
+ //================================================================================
+
+ bool isNodeInSubMesh( const SMDS_MeshNode* n, const SMESHDS_SubMesh* sm )
+ {
+ SMDS_ElemIteratorPtr fIt = n->GetInverseElementIterator( SMDSAbs_Face );
+ while ( fIt->more() )
+ if ( sm->Contains( fIt->next() ))
+ return true;
+ return false;
+ }
+}
+
+//=======================================================================
+//function : IsSameElemGeometry
+//purpose : Returns true if all elements of a sub-mesh are of same shape
+//=======================================================================
+
+bool SMESH_MesherHelper::IsSameElemGeometry(const SMESHDS_SubMesh* smDS,
+ SMDSAbs_GeometryType shape,
+ const bool nullSubMeshRes)
+{
+ if ( !smDS ) return nullSubMeshRes;
+
+ SMDS_ElemIteratorPtr elemIt = smDS->GetElements();
+ while ( elemIt->more() ) {
+ const SMDS_MeshElement* e = elemIt->next();
+ if ( e->GetGeomType() != shape )
+ return false;
+ }
+ return true;
+}
+
+//=======================================================================
+//function : LoadNodeColumns
+//purpose : Load nodes bound to face into a map of node columns
+//=======================================================================
+
+bool SMESH_MesherHelper::LoadNodeColumns(TParam2ColumnMap & theParam2ColumnMap,
+ const TopoDS_Face& theFace,
+ const TopoDS_Edge& theBaseEdge,
+ SMESHDS_Mesh* theMesh,
+ SMESH_ProxyMesh* theProxyMesh)
+{
+ return LoadNodeColumns(theParam2ColumnMap,
+ theFace,
+ std::list<TopoDS_Edge>(1,theBaseEdge),
+ theMesh,
+ theProxyMesh);
+}
+
+//=======================================================================
+//function : LoadNodeColumns
+//purpose : Load nodes bound to face into a map of node columns
+//=======================================================================
+
+bool SMESH_MesherHelper::LoadNodeColumns(TParam2ColumnMap & theParam2ColumnMap,
+ const TopoDS_Face& theFace,
+ const std::list<TopoDS_Edge>& theBaseSide,
+ SMESHDS_Mesh* theMesh,
+ SMESH_ProxyMesh* theProxyMesh)
+{
+ // get a right sub-mesh of theFace
+
+ const SMESHDS_SubMesh* faceSubMesh = 0;
+ if ( theProxyMesh )
+ {
+ faceSubMesh = theProxyMesh->GetSubMesh( theFace );
+ if ( !faceSubMesh ||
+ faceSubMesh->NbElements() == 0 ||
+ theProxyMesh->IsTemporary( faceSubMesh->GetElements()->next() ))
+ {
+ // can use a proxy sub-mesh with not temporary elements only
+ faceSubMesh = 0;
+ theProxyMesh = 0;
+ }
+ }
+ if ( !faceSubMesh )
+ faceSubMesh = theMesh->MeshElements( theFace );
+ if ( !faceSubMesh || faceSubMesh->NbElements() == 0 )
+ return false;
+
+ if ( theParam2ColumnMap.empty() )
+ {
+ // get data of edges for normalization of params
+ vector< double > length;
+ double fullLen = 0;
+ list<TopoDS_Edge>::const_iterator edge;
+ {
+ for ( edge = theBaseSide.begin(); edge != theBaseSide.end(); ++edge )
+ {
+ double len = std::max( 1e-10, SMESH_Algo::EdgeLength( *edge ));
+ fullLen += len;
+ length.push_back( len );
+ }
+ }
+
+ // get nodes on theBaseEdge sorted by param on edge and initialize theParam2ColumnMap with them
+ edge = theBaseSide.begin();
+ for ( int iE = 0; edge != theBaseSide.end(); ++edge, ++iE )
+ {
+ map< double, const SMDS_MeshNode*> sortedBaseNN;
+ SMESH_Algo::GetSortedNodesOnEdge( theMesh, *edge,/*noMedium=*/true, sortedBaseNN);
+ if ( sortedBaseNN.empty() ) continue;
+
+ map< double, const SMDS_MeshNode*>::iterator u_n = sortedBaseNN.begin();
+ if ( theProxyMesh ) // from sortedBaseNN remove nodes not shared by faces of faceSubMesh
+ {
+ const SMDS_MeshNode* n1 = (++sortedBaseNN.begin())->second;
+ const SMDS_MeshNode* n2 = (++sortedBaseNN.rbegin())->second;
+ bool allNodesAreProxy = ( n1 != theProxyMesh->GetProxyNode( n1 ) &&
+ n2 != theProxyMesh->GetProxyNode( n2 ));
+ if ( allNodesAreProxy )
+ for ( u_n = sortedBaseNN.begin(); u_n != sortedBaseNN.end(); u_n++ )
+ u_n->second = theProxyMesh->GetProxyNode( u_n->second );
+
+ if ( u_n = sortedBaseNN.begin(), !isNodeInSubMesh( u_n->second, faceSubMesh ))
+ {
+ while ( ++u_n != sortedBaseNN.end() && !isNodeInSubMesh( u_n->second, faceSubMesh ));
+ sortedBaseNN.erase( sortedBaseNN.begin(), u_n );
+ }
+ if ( u_n = --sortedBaseNN.end(), !isNodeInSubMesh( u_n->second, faceSubMesh ))
+ {
+ while ( u_n != sortedBaseNN.begin() && !isNodeInSubMesh( (--u_n)->second, faceSubMesh ));
+ sortedBaseNN.erase( ++u_n, sortedBaseNN.end() );
+ }
+ if ( sortedBaseNN.empty() ) continue;
+ }
+
+ double f, l;
+ BRep_Tool::Range( *edge, f, l );
+ if ( edge->Orientation() == TopAbs_REVERSED ) std::swap( f, l );
+ const double coeff = 1. / ( l - f ) * length[iE] / fullLen;
+ const double prevPar = theParam2ColumnMap.empty() ? 0 : theParam2ColumnMap.rbegin()->first;
+ for ( u_n = sortedBaseNN.begin(); u_n != sortedBaseNN.end(); u_n++ )
+ {
+ double par = prevPar + coeff * ( u_n->first - f );
+ TParam2ColumnMap::iterator u2nn =
+ theParam2ColumnMap.insert( theParam2ColumnMap.end(), make_pair( par, TNodeColumn()));
+ u2nn->second.push_back( u_n->second );
+ }
+ }
+ if ( theParam2ColumnMap.empty() )
+ return false;
+ }
+
+ // nb rows of nodes
+ int prevNbRows = theParam2ColumnMap.begin()->second.size(); // current, at least 1 here
+ int expectedNbRows = faceSubMesh->NbElements() / ( theParam2ColumnMap.size()-1 ); // to be added
+
+ // fill theParam2ColumnMap column by column by passing from nodes on
+ // theBaseEdge up via mesh faces on theFace
+
+ TParam2ColumnMap::iterator par_nVec_1, par_nVec_2;
+ par_nVec_2 = theParam2ColumnMap.begin();
+ par_nVec_1 = par_nVec_2++;
+ TIDSortedElemSet emptySet, avoidSet;
+ for ( ; par_nVec_2 != theParam2ColumnMap.end(); ++par_nVec_1, ++par_nVec_2 )
+ {
+ vector<const SMDS_MeshNode*>& nCol1 = par_nVec_1->second;
+ vector<const SMDS_MeshNode*>& nCol2 = par_nVec_2->second;
+ nCol1.resize( prevNbRows + expectedNbRows );
+ nCol2.resize( prevNbRows + expectedNbRows );
+
+ int i1, i2, foundNbRows = 0;
+ const SMDS_MeshNode *n1 = nCol1[ prevNbRows-1 ];
+ const SMDS_MeshNode *n2 = nCol2[ prevNbRows-1 ];
+ // find face sharing node n1 and n2 and belonging to faceSubMesh
+ while ( const SMDS_MeshElement* face =
+ SMESH_MeshAlgos::FindFaceInSet( n1, n2, emptySet, avoidSet, &i1, &i2))
+ {
+ if ( faceSubMesh->Contains( face ))
+ {
+ int nbNodes = face->NbCornerNodes();
+ if ( nbNodes != 4 )
+ return false;
+ if ( foundNbRows + 1 > expectedNbRows )
+ return false;
+ n1 = face->GetNode( (i2+2) % 4 ); // opposite corner of quadrangle face
+ n2 = face->GetNode( (i1+2) % 4 );
+ nCol1[ prevNbRows + foundNbRows] = n1;
+ nCol2[ prevNbRows + foundNbRows] = n2;
+ ++foundNbRows;
+ }
+ avoidSet.insert( face );
+ }
+ if ( foundNbRows != expectedNbRows )
+ return false;
+ avoidSet.clear();
+ }
+ return ( theParam2ColumnMap.size() > 1 &&
+ theParam2ColumnMap.begin()->second.size() == prevNbRows + expectedNbRows );
+}
+
+namespace
+{
+ //================================================================================
+ /*!
+ * \brief Return true if a node is at a corner of a 2D structured mesh of FACE
+ */
+ //================================================================================
+
+ bool isCornerOfStructure( const SMDS_MeshNode* n,
+ const SMESHDS_SubMesh* faceSM,
+ SMESH_MesherHelper& faceAnalyser )
+ {
+ int nbFacesInSM = 0;
+ if ( n ) {
+ SMDS_ElemIteratorPtr fIt = n->GetInverseElementIterator( SMDSAbs_Face );
+ while ( fIt->more() )
+ nbFacesInSM += faceSM->Contains( fIt->next() );
+ }
+ if ( nbFacesInSM == 1 )
+ return true;
+
+ if ( nbFacesInSM == 2 && n->GetPosition()->GetTypeOfPosition() == SMDS_TOP_VERTEX )
+ {
+ return faceAnalyser.IsRealSeam( n->getshapeId() );
+ }
+ return false;
+ }
+}
+
+//=======================================================================
+//function : IsStructured
+//purpose : Return true if 2D mesh on FACE is structured
+//=======================================================================
+
+bool SMESH_MesherHelper::IsStructured( SMESH_subMesh* faceSM )
+{
+ SMESHDS_SubMesh* fSM = faceSM->GetSubMeshDS();
+ if ( !fSM || fSM->NbElements() == 0 )
+ return false;
+
+ list< TopoDS_Edge > edges;
+ list< int > nbEdgesInWires;
+ int nbWires = SMESH_Block::GetOrderedEdges( TopoDS::Face( faceSM->GetSubShape() ),
+ edges, nbEdgesInWires );
+ if ( nbWires != 1 /*|| nbEdgesInWires.front() != 4*/ ) // allow composite sides
+ return false;
+
+ // algo: find corners of a structure and then analyze nb of faces and
+ // length of structure sides
+
+ SMESHDS_Mesh* meshDS = faceSM->GetFather()->GetMeshDS();
+ SMESH_MesherHelper faceAnalyser( *faceSM->GetFather() );
+ faceAnalyser.SetSubShape( faceSM->GetSubShape() );
+
+ // rotate edges to get the first node being at corner
+ // (in principle it's not necessary but so far none SALOME algo can make
+ // such a structured mesh that all corner nodes are not on VERTEXes)
+ bool isCorner = false;
+ int nbRemainEdges = nbEdgesInWires.front();
+ do {
+ TopoDS_Vertex V = IthVertex( 0, edges.front() );
+ isCorner = isCornerOfStructure( SMESH_Algo::VertexNode( V, meshDS ),
+ fSM, faceAnalyser);
+ if ( !isCorner ) {
+ edges.splice( edges.end(), edges, edges.begin() );
+ --nbRemainEdges;
+ }
+ }
+ while ( !isCorner && nbRemainEdges > 0 );
+
+ if ( !isCorner )
+ return false;
+
+ // get all nodes from EDGEs
+ list< const SMDS_MeshNode* > nodes;
+ list< TopoDS_Edge >::iterator edge = edges.begin();
+ for ( ; edge != edges.end(); ++edge )
+ {
+ map< double, const SMDS_MeshNode* > u2Nodes;
+ if ( !SMESH_Algo::GetSortedNodesOnEdge( meshDS, *edge,
+ /*skipMedium=*/true, u2Nodes ))
+ return false;
+
+ list< const SMDS_MeshNode* > edgeNodes;
+ map< double, const SMDS_MeshNode* >::iterator u2n = u2Nodes.begin();
+ for ( ; u2n != u2Nodes.end(); ++u2n )
+ edgeNodes.push_back( u2n->second );
+ if ( edge->Orientation() == TopAbs_REVERSED )
+ edgeNodes.reverse();
+
+ if ( !nodes.empty() && nodes.back() == edgeNodes.front() )
+ edgeNodes.pop_front();
+ nodes.splice( nodes.end(), edgeNodes, edgeNodes.begin(), edgeNodes.end() );
+ }
+
+ // get length of structured sides
+ vector<int> nbEdgesInSide;
+ int nbEdges = 0;
+ list< const SMDS_MeshNode* >::iterator n = ++nodes.begin();
+ for ( ; n != nodes.end(); ++n )
+ {
+ ++nbEdges;
+ if ( isCornerOfStructure( *n, fSM, faceAnalyser )) {
+ nbEdgesInSide.push_back( nbEdges );
+ nbEdges = 0;
+ }
+ }
+
+ // checks
+ if ( nbEdgesInSide.size() != 4 )
+ return false;
+ if ( nbEdgesInSide[0] != nbEdgesInSide[2] )
+ return false;
+ if ( nbEdgesInSide[1] != nbEdgesInSide[3] )
+ return false;
+ if ( nbEdgesInSide[0] * nbEdgesInSide[1] != fSM->NbElements() )
+ return false;
+
+ return true;
+}
+
+//================================================================================
+/*!
+ * \brief Find out elements orientation on a geometrical face
+ * \param theFace - The face correctly oriented in the shape being meshed
+ * \retval bool - true if the face normal and the normal of first element
+ * in the correspoding submesh point in different directions
+ */
+//================================================================================
+
+bool SMESH_MesherHelper::IsReversedSubMesh (const TopoDS_Face& theFace)
+{
+ if ( theFace.IsNull() )
+ return false;
+
+ // find out orientation of a meshed face
+ int faceID = GetMeshDS()->ShapeToIndex( theFace );
+ TopoDS_Shape aMeshedFace = GetMeshDS()->IndexToShape( faceID );
+ bool isReversed = ( theFace.Orientation() != aMeshedFace.Orientation() );
+
+ const SMESHDS_SubMesh * aSubMeshDSFace = GetMeshDS()->MeshElements( faceID );
+ if ( !aSubMeshDSFace )
+ return isReversed;
+
+ // find an element with a good normal
+ gp_Vec Ne;
+ bool normalOK = false;
+ gp_XY uv;
+ SMDS_ElemIteratorPtr iteratorElem = aSubMeshDSFace->GetElements();
+ while ( !normalOK && iteratorElem->more() ) // loop on elements on theFace
+ {
+ const SMDS_MeshElement* elem = iteratorElem->next();
+ if ( elem && elem->NbCornerNodes() > 2 )
+ {
+ SMESH_TNodeXYZ nPnt[3];
+ SMDS_ElemIteratorPtr nodesIt = elem->nodesIterator();
+ int iNodeOnFace = 0, iPosDim = SMDS_TOP_VERTEX;
+ for ( int iN = 0; nodesIt->more() && iN < 3; ++iN) // loop on nodes
+ {
+ nPnt[ iN ] = nodesIt->next();
+ if ( nPnt[ iN ]._node->GetPosition()->GetTypeOfPosition() > iPosDim )
+ {
+ iNodeOnFace = iN;
+ iPosDim = nPnt[ iN ]._node->GetPosition()->GetTypeOfPosition();
+ }
+ }
+ // compute normal
+ gp_Vec v01( nPnt[0], nPnt[1] ), v02( nPnt[0], nPnt[2] );
+ if ( v01.SquareMagnitude() > RealSmall() &&
+ v02.SquareMagnitude() > RealSmall() )
+ {
+ Ne = v01 ^ v02;
+ if (( normalOK = ( Ne.SquareMagnitude() > RealSmall() )))
+ uv = GetNodeUV( theFace, nPnt[iNodeOnFace]._node, 0, &normalOK );
+ }
+ }
+ }
+ if ( !normalOK )
+ return isReversed;
+
+ // face normal at node position
+ TopLoc_Location loc;
+ Handle(Geom_Surface) surf = BRep_Tool::Surface( theFace, loc );
+ // if ( surf.IsNull() || surf->Continuity() < GeomAbs_C1 )
+ // some surfaces not detected as GeomAbs_C1 are nevertheless correct for meshing
+ if ( surf.IsNull() || surf->Continuity() < GeomAbs_C0 )
+ {
+ if (!surf.IsNull())
+ MESSAGE("surf->Continuity() < GeomAbs_C1 " << (surf->Continuity() < GeomAbs_C1));
+ return isReversed;
+ }
+ gp_Vec d1u, d1v; gp_Pnt p;
+ surf->D1( uv.X(), uv.Y(), p, d1u, d1v );
+ gp_Vec Nf = (d1u ^ d1v).Transformed( loc );
+
+ if ( theFace.Orientation() == TopAbs_REVERSED )
+ Nf.Reverse();
+
+ return Ne * Nf < 0.;
+}
+
+//=======================================================================
+//function : Count
+//purpose : Count nb of sub-shapes
+//=======================================================================
+
+int SMESH_MesherHelper::Count(const TopoDS_Shape& shape,
+ const TopAbs_ShapeEnum type,
+ const bool ignoreSame)
+{
+ if ( ignoreSame ) {
+ TopTools_IndexedMapOfShape map;
+ TopExp::MapShapes( shape, type, map );
+ return map.Extent();
+ }
+ else {
+ int nb = 0;
+ for ( TopExp_Explorer exp( shape, type ); exp.More(); exp.Next() )
+ ++nb;
+ return nb;
+ }
+}
+
+//=======================================================================
+//function : NbAncestors
+//purpose : Return number of unique ancestors of the shape
+//=======================================================================
+
+int SMESH_MesherHelper::NbAncestors(const TopoDS_Shape& shape,
+ const SMESH_Mesh& mesh,
+ TopAbs_ShapeEnum ancestorType/*=TopAbs_SHAPE*/)
+{
+ TopTools_MapOfShape ancestors;
+ TopTools_ListIteratorOfListOfShape ansIt( mesh.GetAncestors(shape) );
+ for ( ; ansIt.More(); ansIt.Next() ) {
+ if ( ancestorType == TopAbs_SHAPE || ansIt.Value().ShapeType() == ancestorType )
+ ancestors.Add( ansIt.Value() );
+ }
+ return ancestors.Extent();
+}
+
+//=======================================================================
+//function : GetSubShapeOri
+//purpose : Return orientation of sub-shape in the main shape
+//=======================================================================
+
+TopAbs_Orientation SMESH_MesherHelper::GetSubShapeOri(const TopoDS_Shape& shape,
+ const TopoDS_Shape& subShape)
+{
+ TopAbs_Orientation ori = TopAbs_Orientation(-1);
+ if ( !shape.IsNull() && !subShape.IsNull() )
+ {
+ TopExp_Explorer e( shape, subShape.ShapeType() );
+ if ( shape.Orientation() >= TopAbs_INTERNAL ) // TopAbs_INTERNAL or TopAbs_EXTERNAL
+ e.Init( shape.Oriented(TopAbs_FORWARD), subShape.ShapeType() );
+ for ( ; e.More(); e.Next())
+ if ( subShape.IsSame( e.Current() ))
+ break;
+ if ( e.More() )
+ ori = e.Current().Orientation();
+ }
+ return ori;
+}
+
+//=======================================================================
+//function : IsSubShape
+//purpose :
+//=======================================================================
+
+bool SMESH_MesherHelper::IsSubShape( const TopoDS_Shape& shape,
+ const TopoDS_Shape& mainShape )
+{
+ if ( !shape.IsNull() && !mainShape.IsNull() )
+ {
+ for ( TopExp_Explorer exp( mainShape, shape.ShapeType());
+ exp.More();
+ exp.Next() )
+ if ( shape.IsSame( exp.Current() ))
+ return true;
+ }
+ SCRUTE((shape.IsNull()));
+ SCRUTE((mainShape.IsNull()));
+ return false;
+}
+
+//=======================================================================
+//function : IsSubShape
+//purpose :
+//=======================================================================
+
+bool SMESH_MesherHelper::IsSubShape( const TopoDS_Shape& shape, SMESH_Mesh* aMesh )
+{
+ if ( shape.IsNull() || !aMesh )
+ return false;
+ return
+ aMesh->GetMeshDS()->ShapeToIndex( shape ) ||
+ // PAL16202
+ (shape.ShapeType() == TopAbs_COMPOUND && aMesh->GetMeshDS()->IsGroupOfSubShapes( shape ));
+}
+
+//================================================================================
+/*!
+ * \brief Return maximal tolerance of shape
+ */
+//================================================================================
+
+double SMESH_MesherHelper::MaxTolerance( const TopoDS_Shape& shape )
+{
+ double tol = Precision::Confusion();
+ TopExp_Explorer exp;
+ for ( exp.Init( shape, TopAbs_FACE ); exp.More(); exp.Next() )
+ tol = Max( tol, BRep_Tool::Tolerance( TopoDS::Face( exp.Current())));
+ for ( exp.Init( shape, TopAbs_EDGE ); exp.More(); exp.Next() )
+ tol = Max( tol, BRep_Tool::Tolerance( TopoDS::Edge( exp.Current())));
+ for ( exp.Init( shape, TopAbs_VERTEX ); exp.More(); exp.Next() )
+ tol = Max( tol, BRep_Tool::Tolerance( TopoDS::Vertex( exp.Current())));
+
+ return tol;
+}
+
+//================================================================================
+/*!
+ * \brief Check if the first and last vertices of an edge are the same
+ * \param anEdge - the edge to check
+ * \retval bool - true if same
+ */
+//================================================================================
+
+bool SMESH_MesherHelper::IsClosedEdge( const TopoDS_Edge& anEdge )
+{
+ if ( anEdge.Orientation() >= TopAbs_INTERNAL )
+ return IsClosedEdge( TopoDS::Edge( anEdge.Oriented( TopAbs_FORWARD )));
+ return TopExp::FirstVertex( anEdge ).IsSame( TopExp::LastVertex( anEdge ));
+}
+
+//================================================================================
+/*!
+ * \brief Wrapper over TopExp::FirstVertex() and TopExp::LastVertex() fixing them
+ * in the case of INTERNAL edge
+ */
+//================================================================================
+
+TopoDS_Vertex SMESH_MesherHelper::IthVertex( const bool is2nd,
+ TopoDS_Edge anEdge,
+ const bool CumOri )
+{
+ if ( anEdge.Orientation() >= TopAbs_INTERNAL )
+ anEdge.Orientation( TopAbs_FORWARD );
+
+ const TopAbs_Orientation tgtOri = is2nd ? TopAbs_REVERSED : TopAbs_FORWARD;
+ TopoDS_Iterator vIt( anEdge, CumOri );
+ while ( vIt.More() && vIt.Value().Orientation() != tgtOri )
+ vIt.Next();
+
+ return ( vIt.More() ? TopoDS::Vertex(vIt.Value()) : TopoDS_Vertex() );
+}
+
+//================================================================================
+/*!
+ * \brief Return type of shape contained in a group
+ * \param group - a shape of type TopAbs_COMPOUND
+ * \param avoidCompound - not to return TopAbs_COMPOUND
+ */
+//================================================================================
+
+TopAbs_ShapeEnum SMESH_MesherHelper::GetGroupType(const TopoDS_Shape& group,
+ const bool avoidCompound)
+{
+ if ( !group.IsNull() )
+ {
+ if ( group.ShapeType() != TopAbs_COMPOUND )
+ return group.ShapeType();
+
+ // iterate on a compound
+ TopoDS_Iterator it( group );
+ if ( it.More() )
+ return avoidCompound ? GetGroupType( it.Value() ) : it.Value().ShapeType();
+ }
+ return TopAbs_SHAPE;
+}
+
+//=======================================================================
+//function : IsQuadraticMesh
+//purpose : Check mesh without geometry for: if all elements on this shape are quadratic,
+// quadratic elements will be created.
+// Used then generated 3D mesh without geometry.
+//=======================================================================
+
+SMESH_MesherHelper:: MType SMESH_MesherHelper::IsQuadraticMesh()
+{
+ int NbAllEdgsAndFaces=0;
+ int NbQuadFacesAndEdgs=0;
+ int NbFacesAndEdges=0;
+ //All faces and edges
+ NbAllEdgsAndFaces = myMesh->NbEdges() + myMesh->NbFaces();
+ if ( NbAllEdgsAndFaces == 0 )
+ return SMESH_MesherHelper::LINEAR;
+
+ //Quadratic faces and edges
+ NbQuadFacesAndEdgs = myMesh->NbEdges(ORDER_QUADRATIC) + myMesh->NbFaces(ORDER_QUADRATIC);
+
+ //Linear faces and edges
+ NbFacesAndEdges = myMesh->NbEdges(ORDER_LINEAR) + myMesh->NbFaces(ORDER_LINEAR);
+
+ if (NbAllEdgsAndFaces == NbQuadFacesAndEdgs) {
+ //Quadratic mesh
+ return SMESH_MesherHelper::QUADRATIC;
+ }
+ else if (NbAllEdgsAndFaces == NbFacesAndEdges) {
+ //Linear mesh
+ return SMESH_MesherHelper::LINEAR;
+ }
+ else
+ //Mesh with both type of elements
+ return SMESH_MesherHelper::COMP;
+}
+
+//=======================================================================
+//function : GetOtherParam
+//purpose : Return an alternative parameter for a node on seam
+//=======================================================================
+
+double SMESH_MesherHelper::GetOtherParam(const double param) const
+{
+ int i = myParIndex & U_periodic ? 0 : 1;
+ return fabs(param-myPar1[i]) < fabs(param-myPar2[i]) ? myPar2[i] : myPar1[i];
+}
+
+namespace {
+
+ //=======================================================================
+ /*!
+ * \brief Iterator on ancestors of the given type
+ */
+ //=======================================================================
+
+ struct TAncestorsIterator : public SMDS_Iterator<const TopoDS_Shape*>
+ {
+ TopTools_ListIteratorOfListOfShape _ancIter;
+ TopAbs_ShapeEnum _type;
+ TopTools_MapOfShape _encountered;
+ TAncestorsIterator( const TopTools_ListOfShape& ancestors, TopAbs_ShapeEnum type)
+ : _ancIter( ancestors ), _type( type )
+ {
+ if ( _ancIter.More() ) {
+ if ( _ancIter.Value().ShapeType() != _type ) next();
+ else _encountered.Add( _ancIter.Value() );
+ }
+ }
+ virtual bool more()
+ {
+ return _ancIter.More();
+ }
+ virtual const TopoDS_Shape* next()
+ {
+ const TopoDS_Shape* s = _ancIter.More() ? & _ancIter.Value() : 0;
+ if ( _ancIter.More() )
+ for ( _ancIter.Next(); _ancIter.More(); _ancIter.Next())
+ if ( _ancIter.Value().ShapeType() == _type && _encountered.Add( _ancIter.Value() ))
+ break;
+ return s;
+ }
+ };
+
+} // namespace
+
+//=======================================================================
+/*!
+ * \brief Return iterator on ancestors of the given type
+ */
+//=======================================================================
+
+PShapeIteratorPtr SMESH_MesherHelper::GetAncestors(const TopoDS_Shape& shape,
+ const SMESH_Mesh& mesh,
+ TopAbs_ShapeEnum ancestorType)
+{
+ return PShapeIteratorPtr( new TAncestorsIterator( mesh.GetAncestors(shape), ancestorType));
+}
+
+//=======================================================================
+//function : GetCommonAncestor
+//purpose : Find a common ancestors of two shapes of the given type
+//=======================================================================
+
+TopoDS_Shape SMESH_MesherHelper::GetCommonAncestor(const TopoDS_Shape& shape1,
+ const TopoDS_Shape& shape2,
+ const SMESH_Mesh& mesh,
+ TopAbs_ShapeEnum ancestorType)
+{
+ TopoDS_Shape commonAnc;
+ if ( !shape1.IsNull() && !shape2.IsNull() )
+ {
+ PShapeIteratorPtr ancIt = GetAncestors( shape1, mesh, ancestorType );
+ while ( const TopoDS_Shape* anc = ancIt->next() )
+ if ( IsSubShape( shape2, *anc ))
+ {
+ commonAnc = *anc;
+ break;
+ }
+ }
+ return commonAnc;
+}
+
+//#include <Perf_Meter.hxx>
+
+//=======================================================================
+namespace { // Structures used by FixQuadraticElements()
+//=======================================================================
+
+#define __DMP__(txt) \
+ //cout << txt
+#define MSG(txt) __DMP__(txt<<endl)
+#define MSGBEG(txt) __DMP__(txt)
+
+ //const double straightTol2 = 1e-33; // to detect straing links
+ bool isStraightLink(double linkLen2, double middleNodeMove2)
+ {
+ // straight if <node move> < 1/15 * <link length>
+ return middleNodeMove2 < 1/15./15. * linkLen2;
+ }
+
+ struct QFace;
+ // ---------------------------------------
+ /*!
+ * \brief Quadratic link knowing its faces
+ */
+ struct QLink: public SMESH_TLink
+ {
+ const SMDS_MeshNode* _mediumNode;
+ mutable vector<const QFace* > _faces;
+ mutable gp_Vec _nodeMove;
+ mutable int _nbMoves;
+
+ QLink(const SMDS_MeshNode* n1, const SMDS_MeshNode* n2, const SMDS_MeshNode* nm):
+ SMESH_TLink( n1,n2 ), _mediumNode(nm), _nodeMove(0,0,0), _nbMoves(0) {
+ _faces.reserve(4);
+ //if ( MediumPos() != SMDS_TOP_3DSPACE )
+ _nodeMove = MediumPnt() - MiddlePnt();
+ }
+ void SetContinuesFaces() const;
+ const QFace* GetContinuesFace( const QFace* face ) const;
+ bool OnBoundary() const;
+ gp_XYZ MiddlePnt() const { return ( XYZ( node1() ) + XYZ( node2() )) / 2.; }
+ gp_XYZ MediumPnt() const { return XYZ( _mediumNode ); }
+
+ SMDS_TypeOfPosition MediumPos() const
+ { return _mediumNode->GetPosition()->GetTypeOfPosition(); }
+ SMDS_TypeOfPosition EndPos(bool isSecond) const
+ { return (isSecond ? node2() : node1())->GetPosition()->GetTypeOfPosition(); }
+ const SMDS_MeshNode* EndPosNode(SMDS_TypeOfPosition pos) const
+ { return EndPos(0) == pos ? node1() : EndPos(1) == pos ? node2() : 0; }
+
+ void Move(const gp_Vec& move, bool sum=false) const
+ { _nodeMove += move; _nbMoves += sum ? (_nbMoves==0) : 1; }
+ gp_XYZ Move() const { return _nodeMove.XYZ() / _nbMoves; }
+ bool IsMoved() const { return (_nbMoves > 0 /*&& !IsStraight()*/); }
+ bool IsStraight() const
+ { return isStraightLink( (XYZ(node1())-XYZ(node2())).SquareModulus(),
+ _nodeMove.SquareMagnitude());
+ }
+ bool operator<(const QLink& other) const {
+ return (node1()->GetID() == other.node1()->GetID() ?
+ node2()->GetID() < other.node2()->GetID() :
+ node1()->GetID() < other.node1()->GetID());
+ }
+// struct PtrComparator {
+// bool operator() (const QLink* l1, const QLink* l2 ) const { return *l1 < *l2; }
+// };
+ };
+ // ---------------------------------------------------------
+ /*!
+ * \brief Link in the chain of links; it connects two faces
+ */
+ struct TChainLink
+ {
+ const QLink* _qlink;
+ mutable const QFace* _qfaces[2];
+
+ TChainLink(const QLink* qlink=0):_qlink(qlink) {
+ _qfaces[0] = _qfaces[1] = 0;
+ }
+ void SetFace(const QFace* face) const { int iF = _qfaces[0] ? 1 : 0; _qfaces[iF]=face; }
+
+ bool IsBoundary() const { return !_qfaces[1]; }
+
+ void RemoveFace( const QFace* face ) const
+ { _qfaces[(face == _qfaces[1])] = 0; if (!_qfaces[0]) std::swap(_qfaces[0],_qfaces[1]); }
+
+ const QFace* NextFace( const QFace* f ) const
+ { return _qfaces[0]==f ? _qfaces[1] : _qfaces[0]; }
+
+ const SMDS_MeshNode* NextNode( const SMDS_MeshNode* n ) const
+ { return n == _qlink->node1() ? _qlink->node2() : _qlink->node1(); }
+
+ bool operator<(const TChainLink& other) const { return *_qlink < *other._qlink; }
+
+ operator bool() const { return (_qlink); }
+
+ const QLink* operator->() const { return _qlink; }
+
+ gp_Vec Normal() const;
+
+ bool IsStraight() const;
+ };
+ // --------------------------------------------------------------------
+ typedef list< TChainLink > TChain;
+ typedef set < TChainLink > TLinkSet;
+ typedef TLinkSet::const_iterator TLinkInSet;
+
+ const int theFirstStep = 5;
+
+ enum { ERR_OK, ERR_TRI, ERR_PRISM, ERR_UNKNOWN }; // errors of QFace::GetLinkChain()
+ // --------------------------------------------------------------------
+ /*!
+ * \brief Quadratic face shared by two volumes and bound by QLinks
+ */
+ struct QFace: public TIDSortedNodeSet
+ {
+ mutable const SMDS_MeshElement* _volumes[2];
+ mutable vector< const QLink* > _sides;
+ mutable bool _sideIsAdded[4]; // added in chain of links
+ gp_Vec _normal;
+#ifdef _DEBUG_
+ mutable const SMDS_MeshElement* _face;
+#endif
+
+ QFace( const vector< const QLink*>& links, const SMDS_MeshElement* face=0 );
+
+ void SetVolume(const SMDS_MeshElement* v) const { _volumes[ _volumes[0] ? 1 : 0 ] = v; }
+
+ int NbVolumes() const { return !_volumes[0] ? 0 : !_volumes[1] ? 1 : 2; }
+
+ void AddSelfToLinks() const {
+ for ( int i = 0; i < _sides.size(); ++i )
+ _sides[i]->_faces.push_back( this );
+ }
+ int LinkIndex( const QLink* side ) const {
+ for (int i=0; i<_sides.size(); ++i ) if ( _sides[i] == side ) return i;
+ return -1;
+ }
+ bool GetLinkChain( int iSide, TChain& chain, SMDS_TypeOfPosition pos, int& err) const;
+
+ bool GetLinkChain( TChainLink& link, TChain& chain, SMDS_TypeOfPosition pos, int& err) const
+ {
+ int i = LinkIndex( link._qlink );
+ if ( i < 0 ) return true;
+ _sideIsAdded[i] = true;
+ link.SetFace( this );
+ // continue from opposite link
+ return GetLinkChain( (i+2)%_sides.size(), chain, pos, err );
+ }
+ bool IsBoundary() const { return !_volumes[1]; }
+
+ bool Contains( const SMDS_MeshNode* node ) const { return count(node); }
+
+ bool IsSpoiled(const QLink* bentLink ) const;
+
+ TLinkInSet GetBoundaryLink( const TLinkSet& links,
+ const TChainLink& avoidLink,
+ TLinkInSet * notBoundaryLink = 0,
+ const SMDS_MeshNode* nodeToContain = 0,
+ bool * isAdjacentUsed = 0,
+ int nbRecursionsLeft = -1) const;
+
+ TLinkInSet GetLinkByNode( const TLinkSet& links,
+ const TChainLink& avoidLink,
+ const SMDS_MeshNode* nodeToContain) const;
+
+ const SMDS_MeshNode* GetNodeInFace() const {
+ for ( int iL = 0; iL < _sides.size(); ++iL )
+ if ( _sides[iL]->MediumPos() == SMDS_TOP_FACE ) return _sides[iL]->_mediumNode;
+ return 0;
+ }
+
+ gp_Vec LinkNorm(const int i, SMESH_MesherHelper* theFaceHelper=0) const;
+
+ double MoveByBoundary( const TChainLink& theLink,
+ const gp_Vec& theRefVec,
+ const TLinkSet& theLinks,
+ SMESH_MesherHelper* theFaceHelper=0,
+ const double thePrevLen=0,
+ const int theStep=theFirstStep,
+ gp_Vec* theLinkNorm=0,
+ double theSign=1.0) const;
+ };
+
+ //================================================================================
+ /*!
+ * \brief Dump QLink and QFace
+ */
+ ostream& operator << (ostream& out, const QLink& l)
+ {
+ out <<"QLink nodes: "
+ << l.node1()->GetID() << " - "
+ << l._mediumNode->GetID() << " - "
+ << l.node2()->GetID() << endl;
+ return out;
+ }
+ ostream& operator << (ostream& out, const QFace& f)
+ {
+ out <<"QFace nodes: "/*<< &f << " "*/;
+ for ( TIDSortedNodeSet::const_iterator n = f.begin(); n != f.end(); ++n )
+ out << (*n)->GetID() << " ";
+ out << " \tvolumes: "
+ << (f._volumes[0] ? f._volumes[0]->GetID() : 0) << " "
+ << (f._volumes[1] ? f._volumes[1]->GetID() : 0);
+ out << " \tNormal: "<< f._normal.X() <<", "<<f._normal.Y() <<", "<<f._normal.Z() << endl;
+ return out;
+ }
+
+ //================================================================================
+ /*!
+ * \brief Construct QFace from QLinks
+ */
+ //================================================================================
+
+ QFace::QFace( const vector< const QLink*>& links, const SMDS_MeshElement* face )
+ {
+ _volumes[0] = _volumes[1] = 0;
+ _sides = links;
+ _sideIsAdded[0]=_sideIsAdded[1]=_sideIsAdded[2]=_sideIsAdded[3]=false;
+ _normal.SetCoord(0,0,0);
+ for ( int i = 1; i < _sides.size(); ++i ) {
+ const QLink *l1 = _sides[i-1], *l2 = _sides[i];
+ insert( l1->node1() ); insert( l1->node2() );
+ // compute normal
+ gp_Vec v1( XYZ( l1->node2()), XYZ( l1->node1()));
+ gp_Vec v2( XYZ( l2->node1()), XYZ( l2->node2()));
+ if ( l1->node1() != l2->node1() && l1->node2() != l2->node2() )
+ v1.Reverse();
+ _normal += v1 ^ v2;
+ }
+ double normSqSize = _normal.SquareMagnitude();
+ if ( normSqSize > numeric_limits<double>::min() )
+ _normal /= sqrt( normSqSize );
+ else
+ _normal.SetCoord(1e-33,0,0);
+
+#ifdef _DEBUG_
+ _face = face;
+#endif
+ }
+ //================================================================================
+ /*!
+ * \brief Make up a chain of links
+ * \param iSide - link to add first
+ * \param chain - chain to fill in
+ * \param pos - postion of medium nodes the links should have
+ * \param error - out, specifies what is wrong
+ * \retval bool - false if valid chain can't be built; "valid" means that links
+ * of the chain belongs to rectangles bounding hexahedrons
+ */
+ //================================================================================
+
+ bool QFace::GetLinkChain( int iSide, TChain& chain, SMDS_TypeOfPosition pos, int& error) const
+ {
+ if ( iSide >= _sides.size() ) // wrong argument iSide
+ return false;
+ if ( _sideIsAdded[ iSide ]) // already in chain
+ return true;
+
+ if ( _sides.size() != 4 ) { // triangle - visit all my continous faces
+ MSGBEG( *this );
+ TLinkSet links;
+ list< const QFace* > faces( 1, this );
+ while ( !faces.empty() ) {
+ const QFace* face = faces.front();
+ for ( int i = 0; i < face->_sides.size(); ++i ) {
+ if ( !face->_sideIsAdded[i] && face->_sides[i] ) {
+ face->_sideIsAdded[i] = true;
+ // find a face side in the chain
+ TLinkInSet chLink = links.insert( TChainLink(face->_sides[i])).first;
+// TChain::iterator chLink = chain.begin();
+// for ( ; chLink != chain.end(); ++chLink )
+// if ( chLink->_qlink == face->_sides[i] )
+// break;
+// if ( chLink == chain.end() )
+// chLink = chain.insert( chain.begin(), TChainLink(face->_sides[i]));
+ // add a face to a chained link and put a continues face in the queue
+ chLink->SetFace( face );
+ if ( face->_sides[i]->MediumPos() == pos )
+ if ( const QFace* contFace = face->_sides[i]->GetContinuesFace( face ))
+ if ( contFace->_sides.size() == 3 )
+ faces.push_back( contFace );
+ }
+ }
+ faces.pop_front();
+ }
+ if ( error < ERR_TRI )
+ error = ERR_TRI;
+ chain.insert( chain.end(), links.begin(),links.end() );
+ return false;
+ }
+ _sideIsAdded[iSide] = true; // not to add this link to chain again
+ const QLink* link = _sides[iSide];
+ if ( !link)
+ return true;
+
+ // add link into chain
+ TChain::iterator chLink = chain.insert( chain.begin(), TChainLink(link));
+ chLink->SetFace( this );
+ MSGBEG( *this );
+
+ // propagate from a quadrangle to neighbour faces
+ if ( link->MediumPos() >= pos ) {
+ int nbLinkFaces = link->_faces.size();
+ if ( nbLinkFaces == 4 || (/*nbLinkFaces < 4 && */link->OnBoundary())) {
+ // hexahedral mesh or boundary quadrangles - goto a continous face
+ if ( const QFace* f = link->GetContinuesFace( this ))
+ if ( f->_sides.size() == 4 )
+ return f->GetLinkChain( *chLink, chain, pos, error );
+ }
+ else {
+ TChainLink chLink(link); // side face of prismatic mesh - visit all faces of iSide
+ for ( int i = 0; i < nbLinkFaces; ++i )
+ if ( link->_faces[i] )
+ link->_faces[i]->GetLinkChain( chLink, chain, pos, error );
+ if ( error < ERR_PRISM )
+ error = ERR_PRISM;
+ return false;
+ }
+ }
+ return true;
+ }
+
+ //================================================================================
+ /*!
+ * \brief Return a boundary link of the triangle face
+ * \param links - set of all links
+ * \param avoidLink - link not to return
+ * \param notBoundaryLink - out, neither the returned link nor avoidLink
+ * \param nodeToContain - node the returned link must contain; if provided, search
+ * also performed on adjacent faces
+ * \param isAdjacentUsed - returns true if link is found in adjacent faces
+ * \param nbRecursionsLeft - to limit recursion
+ */
+ //================================================================================
+
+ TLinkInSet QFace::GetBoundaryLink( const TLinkSet& links,
+ const TChainLink& avoidLink,
+ TLinkInSet * notBoundaryLink,
+ const SMDS_MeshNode* nodeToContain,
+ bool * isAdjacentUsed,
+ int nbRecursionsLeft) const
+ {
+ TLinkInSet linksEnd = links.end(), boundaryLink = linksEnd;
+
+ typedef list< pair< const QFace*, TLinkInSet > > TFaceLinkList;
+ TFaceLinkList adjacentFaces;
+
+ for ( int iL = 0; iL < _sides.size(); ++iL )
+ {
+ if ( avoidLink._qlink == _sides[iL] )
+ continue;
+ TLinkInSet link = links.find( _sides[iL] );
+ if ( link == linksEnd ) continue;
+ if ( (*link)->MediumPos() > SMDS_TOP_FACE )
+ continue; // We work on faces here, don't go inside a solid
+
+ // check link
+ if ( link->IsBoundary() ) {
+ if ( !nodeToContain ||
+ (*link)->node1() == nodeToContain ||
+ (*link)->node2() == nodeToContain )
+ {
+ boundaryLink = link;
+ if ( !notBoundaryLink ) break;
+ }
+ }
+ else if ( notBoundaryLink ) {
+ *notBoundaryLink = link;
+ if ( boundaryLink != linksEnd ) break;
+ }
+
+ if ( boundaryLink == linksEnd && nodeToContain ) // collect adjacent faces
+ if ( const QFace* adj = link->NextFace( this ))
+ if ( adj->Contains( nodeToContain ))
+ adjacentFaces.push_back( make_pair( adj, link ));
+ }
+
+ if ( isAdjacentUsed ) *isAdjacentUsed = false;
+ if ( boundaryLink == linksEnd && nodeToContain && nbRecursionsLeft) // check adjacent faces
+ {
+ if ( nbRecursionsLeft < 0 )
+ nbRecursionsLeft = nodeToContain->NbInverseElements();
+ TFaceLinkList::iterator adj = adjacentFaces.begin();
+ for ( ; boundaryLink == linksEnd && adj != adjacentFaces.end(); ++adj )
+ boundaryLink = adj->first->GetBoundaryLink( links, *(adj->second), 0, nodeToContain,
+ isAdjacentUsed, nbRecursionsLeft-1);
+ if ( isAdjacentUsed ) *isAdjacentUsed = true;
+ }
+ return boundaryLink;
+ }
+ //================================================================================
+ /*!
+ * \brief Return a link ending at the given node but not avoidLink
+ */
+ //================================================================================
+
+ TLinkInSet QFace::GetLinkByNode( const TLinkSet& links,
+ const TChainLink& avoidLink,
+ const SMDS_MeshNode* nodeToContain) const
+ {
+ for ( int i = 0; i < _sides.size(); ++i )
+ if ( avoidLink._qlink != _sides[i] &&
+ (_sides[i]->node1() == nodeToContain || _sides[i]->node2() == nodeToContain ))
+ return links.find( _sides[ i ]);
+ return links.end();
+ }
+
+ //================================================================================
+ /*!
+ * \brief Return normal to the i-th side pointing outside the face
+ */
+ //================================================================================
+
+ gp_Vec QFace::LinkNorm(const int i, SMESH_MesherHelper* /*uvHelper*/) const
+ {
+ gp_Vec norm, vecOut;
+// if ( uvHelper ) {
+// TopoDS_Face face = TopoDS::Face( uvHelper->GetSubShape());
+// const SMDS_MeshNode* inFaceNode = uvHelper->GetNodeUVneedInFaceNode() ? GetNodeInFace() : 0;
+// gp_XY uv1 = uvHelper->GetNodeUV( face, _sides[i]->node1(), inFaceNode );
+// gp_XY uv2 = uvHelper->GetNodeUV( face, _sides[i]->node2(), inFaceNode );
+// norm.SetCoord( uv1.Y() - uv2.Y(), uv2.X() - uv1.X(), 0 );
+
+// const QLink* otherLink = _sides[(i + 1) % _sides.size()];
+// const SMDS_MeshNode* otherNode =
+// otherLink->node1() == _sides[i]->node1() ? otherLink->node2() : otherLink->node1();
+// gp_XY pIn = uvHelper->GetNodeUV( face, otherNode, inFaceNode );
+// vecOut.SetCoord( uv1.X() - pIn.X(), uv1.Y() - pIn.Y(), 0 );
+// }
+// else {
+ norm = _normal ^ gp_Vec( XYZ(_sides[i]->node1()), XYZ(_sides[i]->node2()));
+ gp_XYZ pIn = ( XYZ( _sides[0]->node1() ) +
+ XYZ( _sides[0]->node2() ) +
+ XYZ( _sides[1]->node1() )) / 3.;
+ vecOut.SetXYZ( _sides[i]->MiddlePnt() - pIn );
+ //}
+ if ( norm * vecOut < 0 )
+ norm.Reverse();
+ double mag2 = norm.SquareMagnitude();
+ if ( mag2 > numeric_limits<double>::min() )
+ norm /= sqrt( mag2 );
+ return norm;
+ }
+ //================================================================================
+ /*!
+ * \brief Move medium node of theLink according to its distance from boundary
+ * \param theLink - link to fix
+ * \param theRefVec - movement of boundary
+ * \param theLinks - all adjacent links of continous triangles
+ * \param theFaceHelper - helper is not used so far
+ * \param thePrevLen - distance from the boundary
+ * \param theStep - number of steps till movement propagation limit
+ * \param theLinkNorm - out normal to theLink
+ * \param theSign - 1 or -1 depending on movement of boundary
+ * \retval double - distance from boundary to propagation limit or other boundary
+ */
+ //================================================================================
+
+ double QFace::MoveByBoundary( const TChainLink& theLink,
+ const gp_Vec& theRefVec,
+ const TLinkSet& theLinks,
+ SMESH_MesherHelper* theFaceHelper,
+ const double thePrevLen,
+ const int theStep,
+ gp_Vec* theLinkNorm,
+ double theSign) const
+ {
+ if ( !theStep )
+ return thePrevLen; // propagation limit reached
+
+ int iL; // index of theLink
+ for ( iL = 0; iL < _sides.size(); ++iL )
+ if ( theLink._qlink == _sides[ iL ])
+ break;
+
+ MSG(string(theStep,'.')<<" Ref( "<<theRefVec.X()<<","<<theRefVec.Y()<<","<<theRefVec.Z()<<" )"
+ <<" thePrevLen " << thePrevLen);
+ MSG(string(theStep,'.')<<" "<<*theLink._qlink);
+
+ gp_Vec linkNorm = -LinkNorm( iL/*, theFaceHelper*/ ); // normal to theLink
+ double refProj = theRefVec * linkNorm; // project movement vector to normal of theLink
+ if ( theStep == theFirstStep )
+ theSign = refProj < 0. ? -1. : 1.;
+ else if ( theSign * refProj < 0.4 * theRefVec.Magnitude())
+ return thePrevLen; // to propagate movement forward only, not in side dir or backward
+
+ int iL1 = (iL + 1) % 3, iL2 = (iL + 2) % 3; // indices of the two other links of triangle
+ TLinkInSet link1 = theLinks.find( _sides[iL1] );
+ TLinkInSet link2 = theLinks.find( _sides[iL2] );
+ if ( link1 == theLinks.end() || link2 == theLinks.end() )
+ return thePrevLen;
+ const QFace* f1 = link1->NextFace( this ); // adjacent faces
+ const QFace* f2 = link2->NextFace( this );
+
+ // propagate to adjacent faces till limit step or boundary
+ double len1 = thePrevLen + (theLink->MiddlePnt() - _sides[iL1]->MiddlePnt()).Modulus();
+ double len2 = thePrevLen + (theLink->MiddlePnt() - _sides[iL2]->MiddlePnt()).Modulus();
+ gp_Vec linkDir1(0,0,0); // initialize to avoid valgrind error ("Conditional jump...")
+ gp_Vec linkDir2(0,0,0);
+ try {
+ OCC_CATCH_SIGNALS;
+ if ( f1 && theLink->MediumPos() <= (*link1)->MediumPos() )
+ len1 = f1->MoveByBoundary
+ ( *link1, theRefVec, theLinks, theFaceHelper, len1, theStep-1, &linkDir1, theSign);
+ else
+ linkDir1 = LinkNorm( iL1/*, theFaceHelper*/ );
+ } catch (...) {
+ MSG( " --------------- EXCEPTION");
+ return thePrevLen;
+ }
+ try {
+ OCC_CATCH_SIGNALS;
+ if ( f2 && theLink->MediumPos() <= (*link2)->MediumPos() )
+ len2 = f2->MoveByBoundary
+ ( *link2, theRefVec, theLinks, theFaceHelper, len2, theStep-1, &linkDir2, theSign);
+ else
+ linkDir2 = LinkNorm( iL2/*, theFaceHelper*/ );
+ } catch (...) {
+ MSG( " --------------- EXCEPTION");
+ return thePrevLen;
+ }
+
+ double fullLen = 0;
+ if ( theStep != theFirstStep )
+ {
+ // choose chain length by direction of propagation most codirected with theRefVec
+ bool choose1 = ( theRefVec * linkDir1 * theSign > theRefVec * linkDir2 * theSign );
+ fullLen = choose1 ? len1 : len2;
+ double r = thePrevLen / fullLen;
+
+ gp_Vec move = linkNorm * refProj * ( 1 - r );
+ theLink->Move( move, true );
+
+ MSG(string(theStep,'.')<<" Move "<< theLink->_mediumNode->GetID()<<
+ " by " << refProj * ( 1 - r ) << " following " <<
+ (choose1 ? *link1->_qlink : *link2->_qlink));
+
+ if ( theLinkNorm ) *theLinkNorm = linkNorm;
+ }
+ return fullLen;
+ }
+
+ //================================================================================
+ /*!
+ * \brief Checks if the face is distorted due to bentLink
+ */
+ //================================================================================
+
+ bool QFace::IsSpoiled(const QLink* bentLink ) const
+ {
+ // code is valid for convex faces only
+ gp_XYZ gc(0,0,0);
+ for ( TIDSortedNodeSet::const_iterator n = begin(); n!=end(); ++n)
+ gc += XYZ( *n ) / size();
+ for (unsigned i = 0; i < _sides.size(); ++i )
+ {
+ if ( _sides[i] == bentLink ) continue;
+ gp_Vec linkNorm = _normal ^ gp_Vec( XYZ(_sides[i]->node1()), XYZ(_sides[i]->node2()));
+ gp_Vec vecOut( gc, _sides[i]->MiddlePnt() );
+ if ( linkNorm * vecOut < 0 )
+ linkNorm.Reverse();
+ double mag2 = linkNorm.SquareMagnitude();
+ if ( mag2 > numeric_limits<double>::min() )
+ linkNorm /= sqrt( mag2 );
+ gp_Vec vecBent ( _sides[i]->MiddlePnt(), bentLink->MediumPnt());
+ gp_Vec vecStraight( _sides[i]->MiddlePnt(), bentLink->MiddlePnt());
+ if ( vecBent * linkNorm > -0.1*vecStraight.Magnitude() )
+ return true;
+ }
+ return false;
+
+ }
+
+ //================================================================================
+ /*!
+ * \brief Find pairs of continues faces
+ */
+ //================================================================================
+
+ void QLink::SetContinuesFaces() const
+ {
+ // x0 x - QLink, [-|] - QFace, v - volume
+ // v0 | v1
+ // | Between _faces of link x2 two vertical faces are continues
+ // x1----x2-----x3 and two horizontal faces are continues. We set vertical faces
+ // | to _faces[0] and _faces[1] and horizontal faces to
+ // v2 | v3 _faces[2] and _faces[3] (or vise versa).
+ // x4
+
+ if ( _faces.empty() )
+ return;
+ int iFaceCont = -1, nbBoundary = 0, iBoundary[2]={-1,-1};
+ if ( _faces[0]->IsBoundary() )
+ iBoundary[ nbBoundary++ ] = 0;
+ for ( int iF = 1; iFaceCont < 0 && iF < _faces.size(); ++iF )
+ {
+ // look for a face bounding none of volumes bound by _faces[0]
+ bool sameVol = false;
+ int nbVol = _faces[iF]->NbVolumes();
+ for ( int iV = 0; !sameVol && iV < nbVol; ++iV )
+ sameVol = ( _faces[iF]->_volumes[iV] == _faces[0]->_volumes[0] ||
+ _faces[iF]->_volumes[iV] == _faces[0]->_volumes[1]);
+ if ( !sameVol )
+ iFaceCont = iF;
+ if ( _faces[iF]->IsBoundary() )
+ iBoundary[ nbBoundary++ ] = iF;
+ }
+ // Set continues faces: arrange _faces to have
+ // _faces[0] continues to _faces[1]
+ // _faces[2] continues to _faces[3]
+ if ( nbBoundary == 2 ) // bnd faces are continues
+ {
+ if (( iBoundary[0] < 2 ) != ( iBoundary[1] < 2 ))
+ {
+ int iNear0 = iBoundary[0] < 2 ? 1-iBoundary[0] : 5-iBoundary[0];
+ std::swap( _faces[ iBoundary[1] ], _faces[iNear0] );
+ }
+ }
+ else if ( iFaceCont > 0 ) // continues faces found
+ {
+ if ( iFaceCont != 1 )
+ std::swap( _faces[1], _faces[iFaceCont] );
+ }
+ else if ( _faces.size() > 1 ) // not found, set NULL by the first face
+ {
+ _faces.insert( ++_faces.begin(), 0 );
+ }
+ }
+ //================================================================================
+ /*!
+ * \brief Return a face continues to the given one
+ */
+ //================================================================================
+
+ const QFace* QLink::GetContinuesFace( const QFace* face ) const
+ {
+ for ( int i = 0; i < _faces.size(); ++i ) {
+ if ( _faces[i] == face ) {
+ int iF = i < 2 ? 1-i : 5-i;
+ return iF < _faces.size() ? _faces[iF] : 0;
+ }
+ }
+ return 0;
+ }
+ //================================================================================
+ /*!
+ * \brief True if link is on mesh boundary
+ */
+ //================================================================================
+
+ bool QLink::OnBoundary() const
+ {
+ for ( int i = 0; i < _faces.size(); ++i )
+ if (_faces[i] && _faces[i]->IsBoundary()) return true;
+ return false;
+ }
+ //================================================================================
+ /*!
+ * \brief Return normal of link of the chain
+ */
+ //================================================================================
+
+ gp_Vec TChainLink::Normal() const {
+ gp_Vec norm;
+ if (_qfaces[0]) norm = _qfaces[0]->_normal;
+ if (_qfaces[1]) norm += _qfaces[1]->_normal;
+ return norm;
+ }
+ //================================================================================
+ /*!
+ * \brief Test link curvature taking into account size of faces
+ */
+ //================================================================================
+
+ bool TChainLink::IsStraight() const
+ {
+ bool isStraight = _qlink->IsStraight();
+ if ( isStraight && _qfaces[0] && !_qfaces[1] )
+ {
+ int i = _qfaces[0]->LinkIndex( _qlink );
+ int iOpp = ( i + 2 ) % _qfaces[0]->_sides.size();
+ gp_XYZ mid1 = _qlink->MiddlePnt();
+ gp_XYZ mid2 = _qfaces[0]->_sides[ iOpp ]->MiddlePnt();
+ double faceSize2 = (mid1-mid2).SquareModulus();
+ isStraight = _qlink->_nodeMove.SquareMagnitude() < 1/10./10. * faceSize2;
+ }
+ return isStraight;
+ }
+
+ //================================================================================
+ /*!
+ * \brief Move medium nodes of vertical links of pentahedrons adjacent by side faces
+ */
+ //================================================================================
+
+ void fixPrism( TChain& allLinks )
+ {
+ // separate boundary links from internal ones
+ typedef set<const QLink*/*, QLink::PtrComparator*/> QLinkSet;
+ QLinkSet interLinks, bndLinks1, bndLink2;
+
+ bool isCurved = false;
+ for ( TChain::iterator lnk = allLinks.begin(); lnk != allLinks.end(); ++lnk ) {
+ if ( (*lnk)->OnBoundary() )
+ bndLinks1.insert( lnk->_qlink );
+ else
+ interLinks.insert( lnk->_qlink );
+ isCurved = isCurved || !lnk->IsStraight();
+ }
+ if ( !isCurved )
+ return; // no need to move
+
+ QLinkSet *curBndLinks = &bndLinks1, *newBndLinks = &bndLink2;
+
+ while ( !interLinks.empty() && !curBndLinks->empty() )
+ {
+ // propagate movement from boundary links to connected internal links
+ QLinkSet::iterator bnd = curBndLinks->begin(), bndEnd = curBndLinks->end();
+ for ( ; bnd != bndEnd; ++bnd )
+ {
+ const QLink* bndLink = *bnd;
+ for ( int i = 0; i < bndLink->_faces.size(); ++i ) // loop on faces of bndLink
+ {
+ const QFace* face = bndLink->_faces[i]; // quadrange lateral face of a prism
+ if ( !face ) continue;
+ // find and move internal link opposite to bndLink within the face
+ int interInd = ( face->LinkIndex( bndLink ) + 2 ) % face->_sides.size();
+ const QLink* interLink = face->_sides[ interInd ];
+ QLinkSet::iterator pInterLink = interLinks.find( interLink );
+ if ( pInterLink == interLinks.end() ) continue; // not internal link
+ interLink->Move( bndLink->_nodeMove );
+ // treated internal links become new boundary ones
+ interLinks.erase( pInterLink );
+ newBndLinks->insert( interLink );
+ }
+ }
+ curBndLinks->clear();
+ std::swap( curBndLinks, newBndLinks );
+ }
+ }
+
+ //================================================================================
+ /*!
+ * \brief Fix links of continues triangles near curved boundary
+ */
+ //================================================================================
+
+ void fixTriaNearBoundary( TChain & allLinks, SMESH_MesherHelper& /*helper*/)
+ {
+ if ( allLinks.empty() ) return;
+
+ TLinkSet linkSet( allLinks.begin(), allLinks.end());
+ TLinkInSet linkIt = linkSet.begin(), linksEnd = linkSet.end();
+
+ for ( linkIt = linkSet.begin(); linkIt != linksEnd; ++linkIt)
+ {
+ if ( linkIt->IsBoundary() && !linkIt->IsStraight() && linkIt->_qfaces[0])
+ {
+ // move iff a boundary link is bent towards inside of a face (issue 0021084)
+ const QFace* face = linkIt->_qfaces[0];
+ gp_XYZ pIn = ( face->_sides[0]->MiddlePnt() +
+ face->_sides[1]->MiddlePnt() +
+ face->_sides[2]->MiddlePnt() ) / 3.;
+ gp_XYZ insideDir( pIn - (*linkIt)->MiddlePnt());
+ bool linkBentInside = ((*linkIt)->_nodeMove.Dot( insideDir ) > 0 );
+ //if ( face->IsSpoiled( linkIt->_qlink ))
+ if ( linkBentInside )
+ face->MoveByBoundary( *linkIt, (*linkIt)->_nodeMove, linkSet );
+ }