+ return true;
+}
+
+//=======================================================================
+//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<int, TopAbs_ShapeEnum>
+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;
+
+ if (( myShapeID == n1->getshapeId() || myShapeID == n2->getshapeId() ) && myShapeID > 0 )
+ {
+ shapeType = myShape.ShapeType();
+ shapeID = myShapeID;
+ }
+ else if ( n1->getshapeId() == n2->getshapeId() )
+ {
+ shapeID = n2->getshapeId();
+ shape = GetSubShapeByNode( n1, GetMeshDS() );
+ }
+ else
+ {
+ const SMDS_TypeOfPosition Pos1 = n1->GetPosition()->GetTypeOfPosition();
+ const SMDS_TypeOfPosition Pos2 = n2->GetPosition()->GetTypeOfPosition();
+
+ if ( Pos1 == SMDS_TOP_3DSPACE || Pos2 == SMDS_TOP_3DSPACE )
+ {
+ }
+ else if ( Pos1 == SMDS_TOP_FACE || Pos2 == SMDS_TOP_FACE )
+ {
+ if ( Pos1 != SMDS_TOP_FACE || Pos2 != SMDS_TOP_FACE )
+ {
+ if ( Pos1 != SMDS_TOP_FACE ) std::swap( n1,n2 );
+ TopoDS_Shape F = GetSubShapeByNode( n1, GetMeshDS() );
+ TopoDS_Shape S = GetSubShapeByNode( n2, GetMeshDS() );
+ if ( IsSubShape( S, F ))
+ {
+ shapeType = TopAbs_FACE;
+ shapeID = n1->getshapeId();
+ }
+ }
+ }
+ else if ( Pos1 == SMDS_TOP_EDGE && Pos2 == SMDS_TOP_EDGE )
+ {
+ TopoDS_Shape E1 = GetSubShapeByNode( n1, GetMeshDS() );
+ TopoDS_Shape E2 = GetSubShapeByNode( n2, GetMeshDS() );
+ shape = GetCommonAncestor( E1, E2, *myMesh, TopAbs_FACE );
+ }
+ else if ( Pos1 == SMDS_TOP_VERTEX && Pos2 == SMDS_TOP_VERTEX )
+ {
+ TopoDS_Shape V1 = GetSubShapeByNode( n1, GetMeshDS() );
+ TopoDS_Shape V2 = GetSubShapeByNode( n2, GetMeshDS() );
+ shape = GetCommonAncestor( V1, V2, *myMesh, TopAbs_EDGE );
+ if ( shape.IsNull() ) shape = GetCommonAncestor( V1, V2, *myMesh, TopAbs_FACE );
+ }
+ else // VERTEX and EDGE
+ {
+ if ( Pos1 != SMDS_TOP_VERTEX ) std::swap( n1,n2 );
+ TopoDS_Shape V = GetSubShapeByNode( n1, GetMeshDS() );
+ TopoDS_Shape E = GetSubShapeByNode( n2, GetMeshDS() );
+ if ( IsSubShape( V, E ))
+ shape = E;
+ else
+ shape = GetCommonAncestor( V, E, *myMesh, TopAbs_FACE );
+ }
+ }
+
+ if ( !shape.IsNull() )
+ {
+ if ( shapeID < 1 )
+ shapeID = GetMeshDS()->ShapeToIndex( shape );
+ shapeType = shape.ShapeType();
+ }
+ 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<TBiQuad, SMDS_MeshNode* >::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 a new medium node between given ones
+//=======================================================================
+
+const SMDS_MeshNode* SMESH_MesherHelper::GetMediumNode(const SMDS_MeshNode* n1,
+ const SMDS_MeshNode* n2,
+ bool force3d)
+{
+ // Find existing node
+
+ SMESH_TLink link(n1,n2);
+ ItTLinkNode itLN = myTLinkNodeMap.find( link );
+ if ( itLN != myTLinkNodeMap.end() ) {
+ return (*itLN).second;
+ }
+
+ // Create medium node
+
+ SMDS_MeshNode* n12;
+ SMESHDS_Mesh* meshDS = GetMeshDS();
+
+ if ( IsSeamShape( n1->getshapeId() ))
+ // to get a correct UV of a node on seam, the second node must have checked UV
+ std::swap( n1, n2 );
+
+ // get type of shape for the new medium node
+ int faceID = -1, edgeID = -1;
+ TopoDS_Edge E; double u [2];
+ TopoDS_Face F; gp_XY uv[2];
+ bool uvOK[2] = { false, false };
+
+ pair<int, TopAbs_ShapeEnum> pos = GetMediumPos( n1, n2, mySetElemOnShape );
+
+ // get positions of the given nodes on shapes
+ if ( pos.second == TopAbs_FACE )
+ {
+ F = TopoDS::Face(meshDS->IndexToShape( faceID = pos.first ));
+ uv[0] = GetNodeUV(F,n1,n2, force3d ? 0 : &uvOK[0]);
+ uv[1] = GetNodeUV(F,n2,n1, force3d ? 0 : &uvOK[1]);
+ }
+ else if ( pos.second == TopAbs_EDGE )
+ {
+ const SMDS_PositionPtr Pos1 = n1->GetPosition();
+ const SMDS_PositionPtr Pos2 = n2->GetPosition();
+ if ( Pos1->GetTypeOfPosition()==SMDS_TOP_EDGE &&
+ Pos2->GetTypeOfPosition()==SMDS_TOP_EDGE &&
+ n1->getshapeId() != n2->getshapeId() )
+ {
+ // issue 0021006
+ return getMediumNodeOnComposedWire(n1,n2,force3d);
+ }
+ E = TopoDS::Edge(meshDS->IndexToShape( edgeID = pos.first ));
+ u[0] = GetNodeU(E,n1,n2, force3d ? 0 : &uvOK[0]);
+ u[1] = GetNodeU(E,n2,n1, force3d ? 0 : &uvOK[1]);
+ }
+
+ if ( !force3d & uvOK[0] && uvOK[1] )
+ {
+ // we try to create medium node using UV parameters of
+ // nodes, else - medium between corresponding 3d points
+ if( ! F.IsNull() )
+ {
+ //if ( uvOK[0] && uvOK[1] )
+ {
+ if ( IsDegenShape( n1->getshapeId() )) {
+ if ( myParIndex & U_periodic ) uv[0].SetCoord( 1, uv[1].Coord( 1 ));
+ else uv[0].SetCoord( 2, uv[1].Coord( 2 ));
+ }
+ else if ( IsDegenShape( n2->getshapeId() )) {
+ if ( myParIndex & U_periodic ) uv[1].SetCoord( 1, uv[0].Coord( 1 ));
+ else uv[1].SetCoord( 2, uv[0].Coord( 2 ));
+ }
+
+ TopLoc_Location loc;
+ Handle(Geom_Surface) S = BRep_Tool::Surface(F,loc);
+ 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());
+ if ( mySetElemOnShape )
+ meshDS->SetNodeOnFace(n12, faceID, UV.X(), UV.Y());
+ myTLinkNodeMap.insert(make_pair(link,n12));
+ return n12;
+ }
+ }
+ else if ( !E.IsNull() )
+ {
+ double f,l;
+ Handle(Geom_Curve) C = BRep_Tool::Curve(E, f, l);
+ if(!C.IsNull())
+ {
+ Standard_Boolean isPeriodic = C->IsPeriodic();
+ double U;
+ if(isPeriodic) {
+ Standard_Real Period = C->Period();
+ Standard_Real p = u[1]+ShapeAnalysis::AdjustByPeriod(u[1],u[0],Period);
+ Standard_Real pmid = (u[0]+p)/2.;
+ U = pmid+ShapeAnalysis::AdjustToPeriod(pmid,C->FirstParameter(),C->LastParameter());
+ }
+ else
+ U = (u[0]+u[1])/2.;
+
+ gp_Pnt P = C->Value( U );
+ n12 = meshDS->AddNode(P.X(), P.Y(), P.Z());
+ if ( mySetElemOnShape )
+ meshDS->SetNodeOnEdge(n12, edgeID, U);
+ myTLinkNodeMap.insert(make_pair(link,n12));
+ return n12;
+ }
+ }
+ }
+
+ // 3d variant
+ double x = ( n1->X() + n2->X() )/2.;
+ double y = ( n1->Y() + n2->Y() )/2.;