X-Git-Url: http://git.salome-platform.org/gitweb/?p=modules%2Fsmesh.git;a=blobdiff_plain;f=src%2FSMESH%2FSMESH_MesherHelper.cxx;h=e1f9709707c252af600fac0d909ea4bc6e5de50e;hp=60d202870cda4ff3f7ee87df0616f084a24d19b1;hb=373c03904b8e3fc5490ff4e17716f0cdcb39c03c;hpb=b02b941bd87bd80d7f573289f6b961d1a1f00167 diff --git a/src/SMESH/SMESH_MesherHelper.cxx b/src/SMESH/SMESH_MesherHelper.cxx index 60d202870..e1f970970 100644 --- a/src/SMESH/SMESH_MesherHelper.cxx +++ b/src/SMESH/SMESH_MesherHelper.cxx @@ -31,14 +31,13 @@ #include "SMDS_FacePosition.hxx" #include "SMDS_IteratorOnIterators.hxx" #include "SMDS_VolumeTool.hxx" +#include "SMESH_Block.hxx" #include "SMESH_ProxyMesh.hxx" #include "SMESH_subMesh.hxx" #include #include -#include #include -#include #include #include #include @@ -82,7 +81,11 @@ namespace { //================================================================================ SMESH_MesherHelper::SMESH_MesherHelper(SMESH_Mesh& theMesh) - : myParIndex(0), myMesh(&theMesh), myShapeID(0), myCreateQuadratic(false), + : myParIndex(0), + myMesh(&theMesh), + myShapeID(0), + myCreateQuadratic(false), + myCreateBiQuadratic(false), myFixNodeParameters(false) { myPar1[0] = myPar2[0] = myPar1[1] = myPar2[1] = 0; @@ -199,7 +202,7 @@ bool SMESH_MesherHelper::IsQuadraticSubMesh(const TopoDS_Shape& aSh) //======================================================================= //function : SetSubShape -//purpose : Set geomerty to make elements on +//purpose : Set geometry to make elements on //======================================================================= void SMESH_MesherHelper::SetSubShape(const int aShID) @@ -214,7 +217,7 @@ void SMESH_MesherHelper::SetSubShape(const int aShID) //======================================================================= //function : SetSubShape -//purpose : Set geomerty to create elements on +//purpose : Set geometry to create elements on //======================================================================= void SMESH_MesherHelper::SetSubShape(const TopoDS_Shape& aSh) @@ -283,6 +286,18 @@ void SMESH_MesherHelper::SetSubShape(const TopoDS_Shape& aSh) myDegenShapeIds.insert( meshDS->ShapeToIndex( v.Current() )); } } + if ( !myDegenShapeIds.empty() && !myParIndex ) { + if ( surface->IsUPeriodic() || surface->IsUClosed() ) { + myParIndex |= U_periodic; + myPar1[0] = surf.FirstUParameter(); + myPar2[0] = surf.LastUParameter(); + } + else if ( surface->IsVPeriodic() || surface->IsVClosed() ) { + myParIndex |= V_periodic; + myPar1[1] = surf.FirstVParameter(); + myPar2[1] = surf.LastVParameter(); + } + } } } } @@ -923,11 +938,18 @@ bool SMESH_MesherHelper::CheckNodeU(const TopoDS_Edge& E, //function : GetMediumPos //purpose : Return index and type of the shape (EDGE or FACE only) to // set a medium node on +//param : useCurSubShape - if true, returns the shape set via SetSubShape() +// if any //======================================================================= -std::pair SMESH_MesherHelper::GetMediumPos(const SMDS_MeshNode* n1, - const SMDS_MeshNode* n2) +std::pair +SMESH_MesherHelper::GetMediumPos(const SMDS_MeshNode* n1, + const SMDS_MeshNode* n2, + const bool useCurSubShape) { + if ( useCurSubShape && !myShape.IsNull() ) + return std::make_pair( myShapeID, myShape.ShapeType() ); + TopAbs_ShapeEnum shapeType = TopAbs_SHAPE; int shapeID = -1; TopoDS_Shape shape; @@ -998,9 +1020,157 @@ std::pair SMESH_MesherHelper::GetMediumPos(const SMDS_Mes return make_pair( shapeID, shapeType ); } +//======================================================================= +//function : GetCentralNode +//purpose : Return existing or create a new central node for a quardilateral +// quadratic face given its 8 nodes. +//@param : force3d - true means node creation in between the given nodes, +// else node position is found on a geometrical face if any. +//======================================================================= + +const SMDS_MeshNode* SMESH_MesherHelper::GetCentralNode(const SMDS_MeshNode* n1, + const SMDS_MeshNode* n2, + const SMDS_MeshNode* n3, + const SMDS_MeshNode* n4, + const SMDS_MeshNode* n12, + const SMDS_MeshNode* n23, + const SMDS_MeshNode* n34, + const SMDS_MeshNode* n41, + bool force3d) +{ + SMDS_MeshNode *centralNode = 0; // central node to return + + // Find an existing central node + + TBiQuad keyOfMap(n1,n2,n3,n4); + std::map::iterator itMapCentralNode; + itMapCentralNode = myMapWithCentralNode.find( keyOfMap ); + if ( itMapCentralNode != myMapWithCentralNode.end() ) + { + return (*itMapCentralNode).second; + } + + // Get type of shape for the new central node + + TopAbs_ShapeEnum shapeType = TopAbs_SHAPE; + int shapeID = -1; + int faceID = -1; + TopoDS_Shape shape; + TopTools_ListIteratorOfListOfShape it; + + std::map< int, int > faceId2nbNodes; + std::map< int, int > ::iterator itMapWithIdFace; + + SMESHDS_Mesh* meshDS = GetMeshDS(); + + // check if a face lie on a FACE, i.e. its all corner nodes lie either on the FACE or + // on sub-shapes of the FACE + if ( GetMesh()->HasShapeToMesh() ) + { + const SMDS_MeshNode* nodes[] = { n1, n2, n3, n4 }; + for(int i = 0; i < 4; i++) + { + shape = GetSubShapeByNode( nodes[i], meshDS ); + if ( shape.IsNull() ) break; + if ( shape.ShapeType() == TopAbs_SOLID ) + { + shapeID = nodes[i]->getshapeId(); + shapeType = TopAbs_SOLID; + break; + } + if ( shape.ShapeType() == TopAbs_FACE ) + { + faceID = nodes[i]->getshapeId(); + itMapWithIdFace = faceId2nbNodes.insert( std::make_pair( faceID, 0 ) ).first; + itMapWithIdFace->second++; + } + else + { + PShapeIteratorPtr it = GetAncestors(shape, *GetMesh(), TopAbs_FACE ); + while ( const TopoDS_Shape* face = it->next() ) + { + faceID = meshDS->ShapeToIndex( *face ); + itMapWithIdFace = faceId2nbNodes.insert( std::make_pair( faceID, 0 ) ).first; + itMapWithIdFace->second++; + } + } + } + } + if ( shapeID < 1 && !faceId2nbNodes.empty() ) // SOLID not found + { + // find ID of the FACE the four corner nodes belong to + itMapWithIdFace = faceId2nbNodes.begin(); + for ( ; itMapWithIdFace != faceId2nbNodes.end(); ++itMapWithIdFace) + { + if ( itMapWithIdFace->second == 4 ) + { + shapeType = TopAbs_FACE; + faceID = (*itMapWithIdFace).first; + break; + } + } + } + + TopoDS_Face F; + if ( shapeType == TopAbs_FACE ) + { + F = TopoDS::Face( meshDS->IndexToShape( faceID )); + } + + // Create a node + + gp_XY uvAvg; + gp_Pnt P; + if ( !F.IsNull() ) + { + if ( !force3d ) + { + uvAvg = calcTFI (0.5, 0.5, + GetNodeUV(F,n1,n3), GetNodeUV(F,n2,n4), + GetNodeUV(F,n3,n1), GetNodeUV(F,n4,n2), + GetNodeUV(F,n12,n3), GetNodeUV(F,n23,n4), + GetNodeUV(F,n34,n2), GetNodeUV(F,n41,n2)); + TopLoc_Location loc; + Handle( Geom_Surface ) S = BRep_Tool::Surface( F, loc ); + P = S->Value( uvAvg.X(), uvAvg.Y() ).Transformed( loc ); + centralNode = meshDS->AddNode( P.X(), P.Y(), P.Z() ); + if ( mySetElemOnShape ) + meshDS->SetNodeOnFace( centralNode, faceID, uvAvg.X(), uvAvg.Y() ); + myMapWithCentralNode.insert( std::make_pair( keyOfMap, centralNode ) ); + return centralNode; + } + } + + P = ( SMESH_TNodeXYZ( n1 ) + + SMESH_TNodeXYZ( n2 ) + + SMESH_TNodeXYZ( n3 ) + + SMESH_TNodeXYZ( n4 ) ) / 4; + centralNode = meshDS->AddNode( P.X(), P.Y(), P.Z() ); + + if ( mySetElemOnShape ) + { + if ( !F.IsNull() ) + { + uvAvg = (GetNodeUV(F,n1,n3) + + GetNodeUV(F,n2,n4) + + GetNodeUV(F,n3,n1) + + GetNodeUV(F,n4,n2)) / 4; + CheckNodeUV( F, centralNode, uvAvg, 2*BRep_Tool::Tolerance( F ), /*force=*/true); + meshDS->SetNodeOnFace( centralNode, faceID, uvAvg.X(), uvAvg.Y() ); + } + else if ( shapeID > 0 ) + meshDS->SetNodeInVolume( centralNode, shapeID ); + else if ( myShapeID > 0 ) + meshDS->SetMeshElementOnShape( centralNode, myShapeID ); + } + + myMapWithCentralNode.insert( std::make_pair( keyOfMap, centralNode ) ); + return centralNode; +} + //======================================================================= //function : GetMediumNode -//purpose : Return existing or create new medium nodes between given ones +//purpose : Return existing or create a new medium node between given ones //======================================================================= const SMDS_MeshNode* SMESH_MesherHelper::GetMediumNode(const SMDS_MeshNode* n1, @@ -1030,7 +1200,7 @@ const SMDS_MeshNode* SMESH_MesherHelper::GetMediumNode(const SMDS_MeshNode* n1, TopoDS_Face F; gp_XY uv[2]; bool uvOK[2] = { false, false }; - pair pos = GetMediumPos( n1, n2 ); + pair pos = GetMediumPos( n1, n2, mySetElemOnShape ); // get positions of the given nodes on shapes if ( pos.second == TopAbs_FACE ) @@ -1077,7 +1247,8 @@ const SMDS_MeshNode* SMESH_MesherHelper::GetMediumNode(const SMDS_MeshNode* n1, gp_XY UV = GetMiddleUV( S, uv[0], uv[1] ); gp_Pnt P = S->Value( UV.X(), UV.Y() ).Transformed(loc); n12 = meshDS->AddNode(P.X(), P.Y(), P.Z()); - meshDS->SetNodeOnFace(n12, faceID, UV.X(), UV.Y()); + if ( mySetElemOnShape ) + meshDS->SetNodeOnFace(n12, faceID, UV.X(), UV.Y()); myTLinkNodeMap.insert(make_pair(link,n12)); return n12; } @@ -1101,7 +1272,8 @@ const SMDS_MeshNode* SMESH_MesherHelper::GetMediumNode(const SMDS_MeshNode* n1, gp_Pnt P = C->Value( U ); n12 = meshDS->AddNode(P.X(), P.Y(), P.Z()); - meshDS->SetNodeOnEdge(n12, edgeID, U); + if ( mySetElemOnShape ) + meshDS->SetNodeOnEdge(n12, edgeID, U); myTLinkNodeMap.insert(make_pair(link,n12)); return n12; } @@ -1114,21 +1286,24 @@ const SMDS_MeshNode* SMESH_MesherHelper::GetMediumNode(const SMDS_MeshNode* n1, double z = ( n1->Z() + n2->Z() )/2.; n12 = meshDS->AddNode(x,y,z); - if ( !F.IsNull() ) - { - gp_XY UV = ( uv[0] + uv[1] ) / 2.; - CheckNodeUV( F, n12, UV, 2*BRep_Tool::Tolerance( F ), /*force=*/true); - meshDS->SetNodeOnFace(n12, faceID, UV.X(), UV.Y() ); - } - else if ( !E.IsNull() ) + if ( mySetElemOnShape ) { - double U = ( u[0] + u[1] ) / 2.; - CheckNodeU( E, n12, U, 2*BRep_Tool::Tolerance( E ), /*force=*/true); - meshDS->SetNodeOnEdge(n12, edgeID, U); - } - else if ( myShapeID > 0 ) - { - meshDS->SetNodeInVolume(n12, myShapeID); + if ( !F.IsNull() ) + { + gp_XY UV = ( uv[0] + uv[1] ) / 2.; + CheckNodeUV( F, n12, UV, 2*BRep_Tool::Tolerance( F ), /*force=*/true); + meshDS->SetNodeOnFace(n12, faceID, UV.X(), UV.Y() ); + } + else if ( !E.IsNull() ) + { + double U = ( u[0] + u[1] ) / 2.; + CheckNodeU( E, n12, U, 2*BRep_Tool::Tolerance( E ), /*force=*/true); + meshDS->SetNodeOnEdge(n12, edgeID, U); + } + else if ( myShapeID > 0 ) + { + meshDS->SetMeshElementOnShape(n12, myShapeID); + } } myTLinkNodeMap.insert( make_pair( link, n12 )); @@ -1200,7 +1375,8 @@ const SMDS_MeshNode* SMESH_MesherHelper::getMediumNodeOnComposedWire(const SMDS_ GetMeshDS()->MoveNode( n12, p.X(), p.Y(), p.Z() ); } - GetMeshDS()->SetNodeOnEdge(n12, edges[iOkEdge], u); + if ( mySetElemOnShape ) + GetMeshDS()->SetNodeOnEdge(n12, edges[iOkEdge], u); myTLinkNodeMap.insert( make_pair( SMESH_TLink(n1,n2), n12 )); @@ -1308,7 +1484,7 @@ SMDS_MeshFace* SMESH_MesherHelper::AddFace(const SMDS_MeshNode* n1, //======================================================================= //function : AddFace -//purpose : Creates quadratic or linear quadrangle +//purpose : Creates bi-quadratic, quadratic or linear quadrangle //======================================================================= SMDS_MeshFace* SMESH_MesherHelper::AddFace(const SMDS_MeshNode* n1, @@ -1351,11 +1527,21 @@ SMDS_MeshFace* SMESH_MesherHelper::AddFace(const SMDS_MeshNode* n1, 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(id) - elem = meshDS->AddFaceWithID(n1, n2, n3, n4, n12, n23, n34, n41, id); + 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 - elem = meshDS->AddFace(n1, n2, n3, n4, n12, n23, n34, n41); + { + 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 ); @@ -1539,7 +1725,7 @@ SMDS_MeshVolume* SMESH_MesherHelper::AddVolume(const SMDS_MeshNode* n1, //======================================================================= //function : AddVolume -//purpose : Creates quadratic or linear hexahedron +//purpose : Creates bi-quadratic, quadratic or linear hexahedron //======================================================================= SMDS_MeshVolume* SMESH_MesherHelper::AddVolume(const SMDS_MeshNode* n1, @@ -1576,15 +1762,75 @@ SMDS_MeshVolume* SMESH_MesherHelper::AddVolume(const SMDS_MeshNode* n1, 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(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); + 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 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 - elem = meshDS->AddVolume(n1, n2, n3, n4, n5, n6, n7, n8, - n12, n23, n34, n41, n56, n67, - n78, n85, n15, n26, n37, n48); + { + 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 ); @@ -1697,6 +1943,26 @@ namespace } } +//======================================================================= +//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 @@ -1726,7 +1992,7 @@ bool SMESH_MesherHelper::LoadNodeColumns(TParam2ColumnMap & theParam2 SMESHDS_Mesh* theMesh, SMESH_ProxyMesh* theProxyMesh) { - // get a right submesh of theFace + // get a right sub-mesh of theFace const SMESHDS_SubMesh* faceSubMesh = 0; if ( theProxyMesh ) @@ -1746,70 +2012,74 @@ bool SMESH_MesherHelper::LoadNodeColumns(TParam2ColumnMap & theParam2 if ( !faceSubMesh || faceSubMesh->NbElements() == 0 ) return false; - // get data of edges for normalization of params - - vector< double > length; - double fullLen = 0; - list::const_iterator edge; + if ( theParam2ColumnMap.empty() ) { - for ( edge = theBaseSide.begin(); edge != theBaseSide.end(); ++edge ) + // get data of edges for normalization of params + + vector< double > length; + double fullLen = 0; + list::const_iterator edge; { - double len = std::max( 1e-10, SMESH_Algo::EdgeLength( *edge )); - fullLen += len; - length.push_back( len ); + 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*> sortedBaseNodes; - SMESH_Algo::GetSortedNodesOnEdge( theMesh, *edge,/*noMedium=*/true, sortedBaseNodes); - if ( sortedBaseNodes.empty() ) continue; - - map< double, const SMDS_MeshNode*>::iterator u_n = sortedBaseNodes.begin(); - if ( theProxyMesh ) // from sortedBaseNodes remove nodes not shared by faces of faceSubMesh + // 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 ) { - const SMDS_MeshNode* n1 = sortedBaseNodes.begin()->second; - const SMDS_MeshNode* n2 = sortedBaseNodes.rbegin()->second; - bool allNodesAreProxy = ( n1 != theProxyMesh->GetProxyNode( n1 ) && - n2 != theProxyMesh->GetProxyNode( n2 )); - if ( allNodesAreProxy ) - for ( u_n = sortedBaseNodes.begin(); u_n != sortedBaseNodes.end(); u_n++ ) - u_n->second = theProxyMesh->GetProxyNode( u_n->second ); + map< double, const SMDS_MeshNode*> sortedBaseNN; + SMESH_Algo::GetSortedNodesOnEdge( theMesh, *edge,/*noMedium=*/true, sortedBaseNN); + if ( sortedBaseNN.empty() ) continue; - if ( u_n = sortedBaseNodes.begin(), !isNodeInSubMesh( u_n->second, faceSubMesh )) + map< double, const SMDS_MeshNode*>::iterator u_n = sortedBaseNN.begin(); + if ( theProxyMesh ) // from sortedBaseNN remove nodes not shared by faces of faceSubMesh { - while ( ++u_n != sortedBaseNodes.end() && !isNodeInSubMesh( u_n->second, faceSubMesh )); - sortedBaseNodes.erase( sortedBaseNodes.begin(), u_n ); + 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 ); + } + else 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; } - else if ( u_n = --sortedBaseNodes.end(), !isNodeInSubMesh( u_n->second, faceSubMesh )) + + 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++ ) { - while ( u_n != sortedBaseNodes.begin() && !isNodeInSubMesh( (--u_n)->second, faceSubMesh )); - sortedBaseNodes.erase( ++u_n, sortedBaseNodes.end() ); + 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 ( sortedBaseNodes.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 = sortedBaseNodes.begin(); u_n != sortedBaseNodes.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; } - if ( theParam2ColumnMap.empty() ) - return false; - - int nbRows = 1 + faceSubMesh->NbElements() / ( theParam2ColumnMap.size()-1 ); + // 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 @@ -1822,35 +2092,253 @@ bool SMESH_MesherHelper::LoadNodeColumns(TParam2ColumnMap & theParam2 { vector& nCol1 = par_nVec_1->second; vector& nCol2 = par_nVec_2->second; - nCol1.resize( nbRows ); - nCol2.resize( nbRows ); + nCol1.resize( prevNbRows + expectedNbRows ); + nCol2.resize( prevNbRows + expectedNbRows ); - int i1, i2, iRow = 0; - const SMDS_MeshNode *n1 = nCol1[0], *n2 = nCol2[0]; + 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_MeshEditor::FindFaceInSet( n1, n2, emptySet, avoidSet, &i1, &i2)) { if ( faceSubMesh->Contains( face )) { - int nbNodes = face->IsQuadratic() ? face->NbNodes()/2 : face->NbNodes(); + 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 ); - if ( ++iRow >= nbRows ) - return false; - nCol1[ iRow ] = n1; - nCol2[ iRow ] = n2; - avoidSet.clear(); + nCol1[ prevNbRows + foundNbRows] = n1; + nCol2[ prevNbRows + foundNbRows] = n2; + ++foundNbRows; } avoidSet.insert( face ); } - // set a real height - nCol1.resize( iRow + 1 ); - nCol2.resize( iRow + 1 ); + 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 ) + 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 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(); + for ( int iN = 0; nodesIt->more() && iN < 3; ++iN) // loop on nodes + nPnt[ iN ] = nodesIt->next(); + + // 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[0]._node, nPnt[2]._node, &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; } - return theParam2ColumnMap.size() > 1 && theParam2ColumnMap.begin()->second.size() > 1; } //======================================================================= @@ -2244,7 +2732,7 @@ namespace { // Structures used by FixQuadraticElements() enum { ERR_OK, ERR_TRI, ERR_PRISM, ERR_UNKNOWN }; // errors of QFace::GetLinkChain() // -------------------------------------------------------------------- /*! - * \brief Face shared by two volumes and bound by QLinks + * \brief Quadratic face shared by two volumes and bound by QLinks */ struct QFace: public TIDSortedNodeSet { @@ -3428,7 +3916,7 @@ namespace { // Structures used by FixQuadraticElements() */ //======================================================================= -void SMESH_MesherHelper::FixQuadraticElements(SMESH_ComputeErrorPtr& error, +void SMESH_MesherHelper::FixQuadraticElements(SMESH_ComputeErrorPtr& compError, bool volumeOnly) { // setenv NO_FixQuadraticElements to know if FixQuadraticElements() is guilty of bad conversion @@ -3462,7 +3950,8 @@ void SMESH_MesherHelper::FixQuadraticElements(SMESH_ComputeErrorPtr& error, #endif SMESH_MesherHelper h(*myMesh); h.SetSubShape( s.Current() ); - h.FixQuadraticElements( error, false ); + h.ToFixNodeParameters(true); + h.FixQuadraticElements( compError, false ); } } // fix nodes on geom faces @@ -3473,12 +3962,12 @@ void SMESH_MesherHelper::FixQuadraticElements(SMESH_ComputeErrorPtr& error, MSG("FIX FACE " << nbfaces-- << " #" << GetMeshDS()->ShapeToIndex(fIt.Key())); SMESH_MesherHelper h(*myMesh); h.SetSubShape( fIt.Key() ); - h.FixQuadraticElements( error, true); h.ToFixNodeParameters(true); + h.FixQuadraticElements( compError, true); } //perf_print_all_meters(1); - if ( error && error->myName == EDITERR_NO_MEDIUM_ON_GEOM ) - error->myComment = "during conversion to quadratic, " + if ( compError && compError->myName == EDITERR_NO_MEDIUM_ON_GEOM ) + compError->myComment = "during conversion to quadratic, " "some medium nodes were not placed on geometry to avoid distorting elements"; return; } @@ -3521,7 +4010,7 @@ void SMESH_MesherHelper::FixQuadraticElements(SMESH_ComputeErrorPtr& error, // Issue 0020982 // Move medium nodes to the link middle for elements whose corner nodes // are out of geometrical boundary to fix distorted elements. - force3DOutOfBoundary( *this, error ); + force3DOutOfBoundary( *this, compError ); if ( elemType == SMDSAbs_Volume ) { @@ -3595,7 +4084,9 @@ void SMESH_MesherHelper::FixQuadraticElements(SMESH_ComputeErrorPtr& error, QLink link( face->GetNode(iN), face->GetNode((iN+1)%nbN), face->GetNode(iN+nbN) ); pLink = links.insert( link ).first; faceLinks[ iN ] = & *pLink; - if ( !isCurved ) + if ( !isCurved && + link.node1()->GetPosition()->GetTypeOfPosition() < 2 && + link.node2()->GetPosition()->GetTypeOfPosition() < 2 ) isCurved = !link.IsStraight(); } // store QFace @@ -3610,11 +4101,13 @@ void SMESH_MesherHelper::FixQuadraticElements(SMESH_ComputeErrorPtr& error, // 3. Compute displacement of medium nodes // --------------------------------------- - // two loops on QFaces: the first is to treat boundary links, the second is for internal ones + // two loops on QFaces: the first is to treat boundary links, the second is for internal ones. TopLoc_Location loc; - // not treat boundary of volumic submesh + bool checkUV; + // not to treat boundary of volumic sub-mesh. int isInside = ( elemType == SMDSAbs_Volume && volumeOnly ) ? 1 : 0; - for ( ; isInside < 2; ++isInside ) { + for ( ; isInside < 2; ++isInside ) + { MSG( "--------------- LOOP (inside=" << isInside << ") ------------------"); SMDS_TypeOfPosition pos = isInside ? SMDS_TOP_3DSPACE : SMDS_TOP_FACE; SMDS_TypeOfPosition bndPos = isInside ? SMDS_TOP_FACE : SMDS_TOP_EDGE; @@ -3689,7 +4182,6 @@ void SMESH_MesherHelper::FixQuadraticElements(SMESH_ComputeErrorPtr& error, gp_Vec move1 = chain.back ()->_nodeMove; TopoDS_Face face; - bool checkUV = true; if ( !isInside ) { // compute node displacement of end links of chain in parametric space of face @@ -3802,10 +4294,131 @@ void SMESH_MesherHelper::FixQuadraticElements(SMESH_ComputeErrorPtr& error, // 4. Move nodes // ------------- + TIDSortedElemSet biQuadQuas, triQuadHexa; + const SMDS_MeshElement *biQuadQua, *triQuadHex; + const bool toFixCentralNodes = ( myMesh->NbBiQuadQuadrangles() + + myMesh->NbTriQuadraticHexas() ); + for ( pLink = links.begin(); pLink != links.end(); ++pLink ) { - if ( pLink->IsMoved() ) { + if ( pLink->IsMoved() ) + { gp_Pnt p = pLink->MiddlePnt() + pLink->Move(); GetMeshDS()->MoveNode( pLink->_mediumNode, p.X(), p.Y(), p.Z()); + + // collect bi-quadratic elements + if ( toFixCentralNodes ) + { + biQuadQua = triQuadHex = 0; + SMDS_ElemIteratorPtr eIt = pLink->_mediumNode->GetInverseElementIterator(); + while ( eIt->more() ) + { + const SMDS_MeshElement* e = eIt->next(); + SMDSAbs_EntityType type = e->GetEntityType(); + if ( type == SMDSEntity_BiQuad_Quadrangle ) + biQuadQuas.insert( e ); + else if ( type == SMDSEntity_TriQuad_Hexa ) + triQuadHexa.insert( e ); + } + } + } + } + + // Fix positions of central nodes of bi-tri-quadratic elements + + // treat bi-quad quadrangles + { + vector< const SMDS_MeshNode* > nodes( 9 ); + gp_XY uv[ 9 ]; + //TIDSortedNodeSet checkedNodes; + TIDSortedElemSet::iterator quadIt = biQuadQuas.begin(); + for ( ; quadIt != biQuadQuas.end(); ++quadIt ) + { + const SMDS_MeshElement* quad = *quadIt; + // nodes + nodes.clear(); + nodes.assign( quad->begin_nodes(), quad->end_nodes() ); + // FACE + TopoDS_Shape S = GetSubShapeByNode( nodes.back(), GetMeshDS() ); + if ( S.IsNull() || S.ShapeType() != TopAbs_FACE ) continue; + const TopoDS_Face& F = TopoDS::Face( S ); + Handle( Geom_Surface ) surf = BRep_Tool::Surface( F, loc ); + const double tol = BRep_Tool::Tolerance( F ); + // UV + for ( int i = 0; i < 8; ++i ) + { + uv[ i ] = GetNodeUV( F, nodes[i], nodes[8], &checkUV ); + // as this method is used after mesh generation, UV of nodes is not + // updated according to bending links, so we update + if ( i > 3 && nodes[i]->GetPosition()->GetTypeOfPosition() == SMDS_TOP_FACE ) + CheckNodeUV( F, nodes[i], uv[ i ], 2*tol, /*force=*/true ); + } + // move central node + gp_XY uvCent = calcTFI (0.5, 0.5, uv[0],uv[1],uv[2],uv[3],uv[4],uv[5],uv[6],uv[7] ); + gp_Pnt p = surf->Value( uvCent.X(), uvCent.Y() ).Transformed( loc ); + GetMeshDS()->MoveNode( nodes[8], p.X(), p.Y(), p.Z()); + } + } + + // treat tri-quadratic hexahedra + { + SMDS_VolumeTool volExp; + TIDSortedElemSet::iterator hexIt = triQuadHexa.begin(); + for ( ; hexIt != triQuadHexa.end(); ++hexIt ) + { + volExp.Set( *hexIt, /*ignoreCentralNodes=*/false ); + + // fix nodes central in sides + for ( int iQuad = 0; iQuad < volExp.NbFaces(); ++iQuad ) + { + const SMDS_MeshNode** quadNodes = volExp.GetFaceNodes( iQuad ); + if ( quadNodes[8]->GetPosition()->GetTypeOfPosition() == SMDS_TOP_3DSPACE ) + { + gp_XYZ p = calcTFI( 0.5, 0.5, + SMESH_TNodeXYZ( quadNodes[0] ), SMESH_TNodeXYZ( quadNodes[2] ), + SMESH_TNodeXYZ( quadNodes[4] ), SMESH_TNodeXYZ( quadNodes[6] ), + SMESH_TNodeXYZ( quadNodes[1] ), SMESH_TNodeXYZ( quadNodes[3] ), + SMESH_TNodeXYZ( quadNodes[5] ), SMESH_TNodeXYZ( quadNodes[7] )); + GetMeshDS()->MoveNode( quadNodes[8], p.X(), p.Y(), p.Z()); + } + } + + // fix the volume central node + vector pointsOnShapes( SMESH_Block::ID_Shell ); + const SMDS_MeshNode** hexNodes = volExp.GetNodes(); + + pointsOnShapes[ SMESH_Block::ID_V000 ] = SMESH_TNodeXYZ( hexNodes[ 0 ] ); + pointsOnShapes[ SMESH_Block::ID_V100 ] = SMESH_TNodeXYZ( hexNodes[ 3 ] ); + pointsOnShapes[ SMESH_Block::ID_V010 ] = SMESH_TNodeXYZ( hexNodes[ 1 ] ); + pointsOnShapes[ SMESH_Block::ID_V110 ] = SMESH_TNodeXYZ( hexNodes[ 2 ] ); + pointsOnShapes[ SMESH_Block::ID_V001 ] = SMESH_TNodeXYZ( hexNodes[ 4 ] ); + pointsOnShapes[ SMESH_Block::ID_V101 ] = SMESH_TNodeXYZ( hexNodes[ 7 ] ); + pointsOnShapes[ SMESH_Block::ID_V011 ] = SMESH_TNodeXYZ( hexNodes[ 5 ] ); + pointsOnShapes[ SMESH_Block::ID_V111 ] = SMESH_TNodeXYZ( hexNodes[ 6 ] ); + + pointsOnShapes[ SMESH_Block::ID_Ex00 ] = SMESH_TNodeXYZ( hexNodes[ 11 ] ); + pointsOnShapes[ SMESH_Block::ID_Ex10 ] = SMESH_TNodeXYZ( hexNodes[ 9 ] ); + pointsOnShapes[ SMESH_Block::ID_E0y0 ] = SMESH_TNodeXYZ( hexNodes[ 8 ] ); + pointsOnShapes[ SMESH_Block::ID_E1y0 ] = SMESH_TNodeXYZ( hexNodes[ 10 ] ); + pointsOnShapes[ SMESH_Block::ID_Ex01 ] = SMESH_TNodeXYZ( hexNodes[ 15 ] ); + pointsOnShapes[ SMESH_Block::ID_Ex11 ] = SMESH_TNodeXYZ( hexNodes[ 13 ] ); + pointsOnShapes[ SMESH_Block::ID_E0y1 ] = SMESH_TNodeXYZ( hexNodes[ 12 ] ); + pointsOnShapes[ SMESH_Block::ID_E1y1 ] = SMESH_TNodeXYZ( hexNodes[ 14 ] ); + pointsOnShapes[ SMESH_Block::ID_E00z ] = SMESH_TNodeXYZ( hexNodes[ 16 ] ); + pointsOnShapes[ SMESH_Block::ID_E10z ] = SMESH_TNodeXYZ( hexNodes[ 19 ] ); + pointsOnShapes[ SMESH_Block::ID_E01z ] = SMESH_TNodeXYZ( hexNodes[ 17 ] ); + pointsOnShapes[ SMESH_Block::ID_E11z ] = SMESH_TNodeXYZ( hexNodes[ 18 ] ); + + pointsOnShapes[ SMESH_Block::ID_Fxy0 ] = SMESH_TNodeXYZ( hexNodes[ 20 ] ); + pointsOnShapes[ SMESH_Block::ID_Fxy1 ] = SMESH_TNodeXYZ( hexNodes[ 25 ] ); + pointsOnShapes[ SMESH_Block::ID_Fx0z ] = SMESH_TNodeXYZ( hexNodes[ 21 ] ); + pointsOnShapes[ SMESH_Block::ID_Fx1z ] = SMESH_TNodeXYZ( hexNodes[ 23 ] ); + pointsOnShapes[ SMESH_Block::ID_F0yz ] = SMESH_TNodeXYZ( hexNodes[ 24 ] ); + pointsOnShapes[ SMESH_Block::ID_F1yz ] = SMESH_TNodeXYZ( hexNodes[ 22 ] ); + + gp_XYZ nCenterParams(0.5, 0.5, 0.5), nCenterCoords; + SMESH_Block::ShellPoint( nCenterParams, pointsOnShapes, nCenterCoords ); + GetMeshDS()->MoveNode( hexNodes[26], + nCenterCoords.X(), nCenterCoords.Y(), nCenterCoords.Z()); } }