C++ SMESHDS_Mesh类代码示例

PyObject* FemMeshPy::addQuad(PyObject *args)
{
    int n1,n2,n3,n4;
    if (!PyArg_ParseTuple(args, "iiii",&n1,&n2,&n3,&n4))
        return 0;

    try {
        SMESH_Mesh* mesh = getFemMeshPtr()->getSMesh();
        SMESHDS_Mesh* meshDS = mesh->GetMeshDS();
        const SMDS_MeshNode* node1 = meshDS->FindNode(n1);
        const SMDS_MeshNode* node2 = meshDS->FindNode(n2);
        const SMDS_MeshNode* node3 = meshDS->FindNode(n3);
        const SMDS_MeshNode* node4 = meshDS->FindNode(n4);
        if (!node1 || !node2 || !node3 || !node4)
            throw std::runtime_error("Failed to get node of the given indices");
        SMDS_MeshFace* face = meshDS->AddFace(node1, node2, node3, node4);
        if (!face)
            throw std::runtime_error("Failed to add quad");
        return Py::new_reference_to(Py::Long(face->GetID()));
    }
    catch (const std::exception& e) {
        PyErr_SetString(Base::BaseExceptionFreeCADError, e.what());
        return 0;
    }
}

TNodeColumn* StdMeshers_RadialPrism_3D::makeNodeColumn( TNode2ColumnMap&     n2ColMap,
                                                        const SMDS_MeshNode* outNode,
                                                        const SMDS_MeshNode* inNode)
{
  SMESHDS_Mesh * meshDS = myHelper->GetMeshDS();
  int shapeID = myHelper->GetSubShapeID();

  if ( myLayerPositions.empty() ) {
    gp_Pnt pIn = gpXYZ( inNode ), pOut = gpXYZ( outNode );
    computeLayerPositions( pIn, pOut );
  }
  int nbSegments = myLayerPositions.size() + 1;

  TNode2ColumnMap::iterator n_col =
    n2ColMap.insert( make_pair( outNode, TNodeColumn() )).first;
  TNodeColumn & column = n_col->second;
  column.resize( nbSegments + 1 );
  column.front() = outNode;
  column.back() =  inNode;

  gp_XYZ p1 = gpXYZ( outNode );
  gp_XYZ p2 = gpXYZ( inNode );
  for ( int z = 1; z < nbSegments; ++z )
  {
    double r = myLayerPositions[ z - 1 ];
    gp_XYZ p = ( 1 - r ) * p1 + r * p2;
    SMDS_MeshNode* n = meshDS->AddNode( p.X(), p.Y(), p.Z() );
    meshDS->SetNodeInVolume( n, shapeID );
    column[ z ] = n;
  }

  return & column;
}

PyObject* FemMeshPy::addVolume(PyObject *args)
{
    int n1,n2,n3,n4;
    if (!PyArg_ParseTuple(args, "iiii",&n1,&n2,&n3,&n4))
        return 0;

    try {
        SMESH_Mesh* mesh = getFemMeshPtr()->getSMesh();
        SMESHDS_Mesh* meshDS = mesh->GetMeshDS();
        const SMDS_MeshNode* node1 = meshDS->FindNode(n1);
        const SMDS_MeshNode* node2 = meshDS->FindNode(n2);
        const SMDS_MeshNode* node3 = meshDS->FindNode(n3);
        const SMDS_MeshNode* node4 = meshDS->FindNode(n4);
        if (!node1 || !node2 || !node3 || !node4)
            throw std::runtime_error("Failed to get node of the given indices");
        SMDS_MeshVolume* vol = meshDS->AddVolume(node1, node2, node3, node4);
        if (!vol)
            throw std::runtime_error("Failed to add volume");
        return Py::new_reference_to(Py::Int(vol->GetID()));
    }
    catch (const std::exception& e) {
        PyErr_SetString(PyExc_Exception, e.what());
        return 0;
    }
}

bool SMESH_MesherHelper::IsQuadraticSubMesh(const TopoDS_Shape& aSh)
{
  SMESHDS_Mesh* meshDS = GetMeshDS();
  // we can create quadratic elements only if all elements
  // created on subshapes of given shape are quadratic
  // also we have to fill myNLinkNodeMap
  myCreateQuadratic = true;
  mySeamShapeIds.clear();
  myDegenShapeIds.clear();
  TopAbs_ShapeEnum subType( aSh.ShapeType()==TopAbs_FACE ? TopAbs_EDGE : TopAbs_FACE );
  SMDSAbs_ElementType elemType( subType==TopAbs_FACE ? SMDSAbs_Face : SMDSAbs_Edge );

  int nbOldLinks = myNLinkNodeMap.size();

  TopExp_Explorer exp( aSh, subType );
  for (; exp.More() && myCreateQuadratic; exp.Next()) {
    if ( SMESHDS_SubMesh * subMesh = meshDS->MeshElements( exp.Current() )) {
      if ( SMDS_ElemIteratorPtr it = subMesh->GetElements() ) {
        while(it->more()) {
          const SMDS_MeshElement* e = it->next();
          if ( e->GetType() != elemType || !e->IsQuadratic() ) {
            myCreateQuadratic = false;
            break;
          }
          else {
            // fill NLinkNodeMap
            switch ( e->NbNodes() ) {
            case 3:
              AddNLinkNode(e->GetNode(0),e->GetNode(1),e->GetNode(2)); break;
            case 6:
              AddNLinkNode(e->GetNode(0),e->GetNode(1),e->GetNode(3));
              AddNLinkNode(e->GetNode(1),e->GetNode(2),e->GetNode(4));
              AddNLinkNode(e->GetNode(2),e->GetNode(0),e->GetNode(5)); break;
            case 8:
              AddNLinkNode(e->GetNode(0),e->GetNode(1),e->GetNode(4));
              AddNLinkNode(e->GetNode(1),e->GetNode(2),e->GetNode(5));
              AddNLinkNode(e->GetNode(2),e->GetNode(3),e->GetNode(6));
              AddNLinkNode(e->GetNode(3),e->GetNode(0),e->GetNode(7));
              break;
            default:
              myCreateQuadratic = false;
              break;
            }
          }
        }
      }
    }
  }

  if ( nbOldLinks == myNLinkNodeMap.size() )
    myCreateQuadratic = false;

  if(!myCreateQuadratic) {
    myNLinkNodeMap.clear();
  }
  SetSubShape( aSh );

  return myCreateQuadratic;
}

App::DocumentObjectExecReturn *FemMeshShapeNetgenObject::execute(void)
{
#ifdef FCWithNetgen

    Fem::FemMesh newMesh;

    Part::Feature *feat = Shape.getValue<Part::Feature*>();
    TopoDS_Shape shape = feat->Shape.getValue();

    NETGENPlugin_Mesher myNetGenMesher(newMesh.getSMesh(),shape,true);
    NETGENPlugin_Hypothesis* tet= new NETGENPlugin_Hypothesis(0,1,newMesh.getGenerator());
    tet->SetMaxSize(MaxSize.getValue());
    tet->SetSecondOrder(SecondOrder.getValue());
    tet->SetOptimize(Optimize.getValue());
    int iFineness = Fineness.getValue();
    tet->SetFineness((NETGENPlugin_Hypothesis::Fineness)iFineness);
    if(iFineness == 5){
        tet->SetGrowthRate(GrowthRate.getValue());
        tet->SetNbSegPerEdge(NbSegsPerEdge.getValue());
        tet->SetNbSegPerRadius(NbSegsPerRadius.getValue());
    }
    myNetGenMesher.SetParameters( tet);
    newMesh.getSMesh()->ShapeToMesh(shape);

    myNetGenMesher.Compute();

    // throw Base::RuntimeError("Compute Done\n");

    SMESHDS_Mesh* data = const_cast<SMESH_Mesh*>(newMesh.getSMesh())->GetMeshDS();
    const SMDS_MeshInfo& info = data->GetMeshInfo();
    int numFaces = data->NbFaces();
    int numNode = info.NbNodes();
  //int numTria = info.NbTriangles();
  //int numQuad = info.NbQuadrangles();
  //int numPoly = info.NbPolygons();
    int numVolu = info.NbVolumes();
  //int numTetr = info.NbTetras();
  //int numHexa = info.NbHexas();
  //int numPyrd = info.NbPyramids();
  //int numPris = info.NbPrisms();
  //int numHedr = info.NbPolyhedrons();

    Base::Console().Log("NetgenMesh: %i Nodes, %i Volumes, %i Faces\n",numNode,numVolu,numFaces);

      FemMesh.setValue(newMesh);
    return App::DocumentObject::StdReturn;
#else
    return new App::DocumentObjectExecReturn("The FEM module is built without NETGEN support. Meshing will not work!!!", this);
#endif
}

double SMESH_MesherHelper::GetNodeU(const TopoDS_Edge&   E,
                                    const SMDS_MeshNode* n)
{
  double param = 0;
  const SMDS_PositionPtr Pos = n->GetPosition();
  if(Pos->GetTypeOfPosition()==SMDS_TOP_EDGE) {
    const SMDS_EdgePosition* epos =
      static_cast<const SMDS_EdgePosition*>(n->GetPosition().get());
    param =  epos->GetUParameter();
  }
  else if(Pos->GetTypeOfPosition()==SMDS_TOP_VERTEX) {
    SMESHDS_Mesh * meshDS = GetMeshDS();
    int vertexID = n->GetPosition()->GetShapeId();
    const TopoDS_Vertex& V = TopoDS::Vertex(meshDS->IndexToShape(vertexID));
    param =  BRep_Tool::Parameter( V, E );
  }
  return param;
}

PyObject* FemMeshPy::addNode(PyObject *args)
{
    double x,y,z;
    int i = -1;
    if (PyArg_ParseTuple(args, "ddd",&x,&y,&z)){
        try {
            SMESH_Mesh* mesh = getFemMeshPtr()->getSMesh();
            SMESHDS_Mesh* meshDS = mesh->GetMeshDS();
            SMDS_MeshNode* node = meshDS->AddNode(x,y,z);
            if (!node)
                throw std::runtime_error("Failed to add node");
            return Py::new_reference_to(Py::Int(node->GetID()));
        }
        catch (const std::exception& e) {
            PyErr_SetString(Base::BaseExceptionFreeCADError, e.what());
            return 0;
        }
    }
    PyErr_Clear();

    if (PyArg_ParseTuple(args, "dddi",&x,&y,&z,&i)){
        try {
            SMESH_Mesh* mesh = getFemMeshPtr()->getSMesh();
            SMESHDS_Mesh* meshDS = mesh->GetMeshDS();
            SMDS_MeshNode* node = meshDS->AddNodeWithID(x,y,z,i);
            if (!node)
                throw std::runtime_error("Failed to add node");
            return Py::new_reference_to(Py::Int(node->GetID()));
        }
        catch (const std::exception& e) {
            PyErr_SetString(Base::BaseExceptionFreeCADError, e.what());
            return 0;
        }
    }
    PyErr_SetString(PyExc_TypeError, "addNode() accepts:\n"
        "-- addNode(x,y,z)\n"
        "-- addNode(x,y,z,ElemId)\n");
    return 0;

}

void SMESHDS_GroupOnGeom::SetShape( const TopoDS_Shape& theShape)
{
  SMESHDS_Mesh* aMesh = const_cast<SMESHDS_Mesh*>( GetMesh() );
  mySubMesh = aMesh->MeshElements( aMesh->AddCompoundSubmesh( theShape ));
  myShape   = theShape;
}

bool NETGENPlugin_NETGEN_3D::Compute(SMESH_Mesh&         aMesh,
                                     const TopoDS_Shape& aShape)
{
    MESSAGE("NETGENPlugin_NETGEN_3D::Compute with maxElmentsize = " << _maxElementVolume);

    SMESHDS_Mesh* meshDS = aMesh.GetMeshDS();

    const int invalid_ID = -1;

    SMESH::Controls::Area areaControl;
    SMESH::Controls::TSequenceOfXYZ nodesCoords;

    // -------------------------------------------------------------------
    // get triangles on aShell and make a map of nodes to Netgen node IDs
    // -------------------------------------------------------------------

    SMESH_MesherHelper helper(aMesh);
    SMESH_MesherHelper* myTool = &helper;
    bool _quadraticMesh = myTool->IsQuadraticSubMesh(aShape);

    typedef map< const SMDS_MeshNode*, int, TIDCompare > TNodeToIDMap;
    TNodeToIDMap nodeToNetgenID;
    list< const SMDS_MeshElement* > triangles;
    list< bool >                    isReversed; // orientation of triangles

    TopAbs_ShapeEnum mainType = aMesh.GetShapeToMesh().ShapeType();
    bool checkReverse = ( mainType == TopAbs_COMPOUND || mainType == TopAbs_COMPSOLID );

    // for the degeneraged edge: ignore all but one node on it;
    // map storing ids of degen edges and vertices and their netgen id:
    map< int, int* > degenShapeIdToPtrNgId;
    map< int, int* >::iterator shId_ngId;
    list< int > degenNgIds;

    StdMeshers_QuadToTriaAdaptor Adaptor;
    Adaptor.Compute(aMesh,aShape);

    for (TopExp_Explorer exp(aShape,TopAbs_FACE); exp.More(); exp.Next())
    {
        const TopoDS_Shape& aShapeFace = exp.Current();
        const SMESHDS_SubMesh * aSubMeshDSFace = meshDS->MeshElements( aShapeFace );
        if ( aSubMeshDSFace )
        {
            bool isRev = false;
            if ( checkReverse && helper.NbAncestors(aShapeFace, aMesh, aShape.ShapeType()) > 1 )
                // IsReversedSubMesh() can work wrong on strongly curved faces,
                // so we use it as less as possible
                isRev = SMESH_Algo::IsReversedSubMesh( TopoDS::Face(aShapeFace), meshDS );

            SMDS_ElemIteratorPtr iteratorElem = aSubMeshDSFace->GetElements();
            while ( iteratorElem->more() ) // loop on elements on a face
            {
                // check element
                const SMDS_MeshElement* elem = iteratorElem->next();
                if ( !elem )
                    return error( COMPERR_BAD_INPUT_MESH, "Null element encounters");
                bool isTraingle = ( elem->NbNodes()==3 || (_quadraticMesh && elem->NbNodes()==6 ));
                if ( !isTraingle ) {
                    //return error( COMPERR_BAD_INPUT_MESH,
                    //              SMESH_Comment("Not triangle element ")<<elem->GetID());
                    // using adaptor
                    const list<const SMDS_FaceOfNodes*>* faces = Adaptor.GetTriangles(elem);
                    if(faces==0) {
                        return error( COMPERR_BAD_INPUT_MESH,
                                      SMESH_Comment("Not triangles in adaptor for element ")<<elem->GetID());
                    }
                    list<const SMDS_FaceOfNodes*>::const_iterator itf = faces->begin();
                    for(; itf!=faces->end(); itf++ ) {
                        triangles.push_back( (*itf) );
                        isReversed.push_back( isRev );
                        // put triange's nodes to nodeToNetgenID map
                        SMDS_ElemIteratorPtr triangleNodesIt = (*itf)->nodesIterator();
                        while ( triangleNodesIt->more() ) {
                            const SMDS_MeshNode * node =
                                static_cast<const SMDS_MeshNode *>(triangleNodesIt->next());
                            if(myTool->IsMedium(node))
                                continue;
                            nodeToNetgenID.insert( make_pair( node, invalid_ID ));
                        }
                    }
                }
                else {
                    // keep a triangle
                    triangles.push_back( elem );
                    isReversed.push_back( isRev );
                    // put elem nodes to nodeToNetgenID map
                    SMDS_ElemIteratorPtr triangleNodesIt = elem->nodesIterator();
                    while ( triangleNodesIt->more() ) {
                        const SMDS_MeshNode * node =
                            static_cast<const SMDS_MeshNode *>(triangleNodesIt->next());
                        if(myTool->IsMedium(node))
                            continue;
                        nodeToNetgenID.insert( make_pair( node, invalid_ID ));
                    }
                }
#ifdef _DEBUG_
                // check if a trainge is degenerated
                areaControl.GetPoints( elem, nodesCoords );
                double area = areaControl.GetValue( nodesCoords );
                if ( area <= DBL_MIN ) {
//.........这里部分代码省略.........

bool NETGENPlugin_NETGEN_3D::Compute(SMESH_Mesh&         aMesh,
                                     const TopoDS_Shape& aShape)
{
  netgen::multithread.terminate = 0;
  netgen::multithread.task = "Volume meshing";
  _progressByTic = -1.;

  SMESHDS_Mesh* meshDS = aMesh.GetMeshDS();

  SMESH_MesherHelper helper(aMesh);
  bool _quadraticMesh = helper.IsQuadraticSubMesh(aShape);
  helper.SetElementsOnShape( true );

  int Netgen_NbOfNodes = 0;
  double Netgen_point[3];
  int Netgen_triangle[3];

  NETGENPlugin_NetgenLibWrapper ngLib;
  Ng_Mesh * Netgen_mesh = ngLib._ngMesh;

  // vector of nodes in which node index == netgen ID
  vector< const SMDS_MeshNode* > nodeVec;
  {
    const int invalid_ID = -1;

    SMESH::Controls::Area areaControl;
    SMESH::Controls::TSequenceOfXYZ nodesCoords;

    // maps nodes to ng ID
    typedef map< const SMDS_MeshNode*, int, TIDCompare > TNodeToIDMap;
    typedef TNodeToIDMap::value_type                     TN2ID;
    TNodeToIDMap nodeToNetgenID;

    // find internal shapes
    NETGENPlugin_Internals internals( aMesh, aShape, /*is3D=*/true );

    // ---------------------------------
    // Feed the Netgen with surface mesh
    // ---------------------------------

    TopAbs_ShapeEnum mainType = aMesh.GetShapeToMesh().ShapeType();
    bool checkReverse = ( mainType == TopAbs_COMPOUND || mainType == TopAbs_COMPSOLID );

    SMESH_ProxyMesh::Ptr proxyMesh( new SMESH_ProxyMesh( aMesh ));
    if ( _viscousLayersHyp )
    {
      netgen::multithread.percent = 3;
      proxyMesh = _viscousLayersHyp->Compute( aMesh, aShape );
      if ( !proxyMesh )
        return false;
    }
    if ( aMesh.NbQuadrangles() > 0 )
    {
      netgen::multithread.percent = 6;
      StdMeshers_QuadToTriaAdaptor* Adaptor = new StdMeshers_QuadToTriaAdaptor;
      Adaptor->Compute(aMesh,aShape,proxyMesh.get());
      proxyMesh.reset( Adaptor );
    }

    for ( TopExp_Explorer exFa( aShape, TopAbs_FACE ); exFa.More(); exFa.Next())
    {
      const TopoDS_Shape& aShapeFace = exFa.Current();
      int faceID = meshDS->ShapeToIndex( aShapeFace );
      bool isInternalFace = internals.isInternalShape( faceID );
      bool isRev = false;
      if ( checkReverse && !isInternalFace &&
           helper.NbAncestors(aShapeFace, aMesh, aShape.ShapeType()) > 1 )
        // IsReversedSubMesh() can work wrong on strongly curved faces,
        // so we use it as less as possible
        isRev = helper.IsReversedSubMesh( TopoDS::Face( aShapeFace ));

      const SMESHDS_SubMesh * aSubMeshDSFace = proxyMesh->GetSubMesh( aShapeFace );
      if ( !aSubMeshDSFace ) continue;
      SMDS_ElemIteratorPtr iteratorElem = aSubMeshDSFace->GetElements();
      while ( iteratorElem->more() ) // loop on elements on a geom face
      {
        // check mesh face
        const SMDS_MeshElement* elem = iteratorElem->next();
        if ( !elem )
          return error( COMPERR_BAD_INPUT_MESH, "Null element encounters");
        if ( elem->NbCornerNodes() != 3 )
          return error( COMPERR_BAD_INPUT_MESH, "Not triangle element encounters");

        // Add nodes of triangles and triangles them-selves to netgen mesh

        // add three nodes of triangle
        bool hasDegen = false;
        for ( int iN = 0; iN < 3; ++iN )
        {
          const SMDS_MeshNode* node = elem->GetNode( iN );
          const int shapeID = node->getshapeId();
          if ( node->GetPosition()->GetTypeOfPosition() == SMDS_TOP_EDGE &&
               helper.IsDegenShape( shapeID ))
          {
            // ignore all nodes on degeneraged edge and use node on its vertex instead
            TopoDS_Shape vertex = TopoDS_Iterator( meshDS->IndexToShape( shapeID )).Value();
            node = SMESH_Algo::VertexNode( TopoDS::Vertex( vertex ), meshDS );
            hasDegen = true;
          }
          int& ngID = nodeToNetgenID.insert(TN2ID( node, invalid_ID )).first->second;
//.........这里部分代码省略.........

bool NETGENPlugin_NETGEN_2D_ONLY::Compute(SMESH_Mesh&         aMesh,
                                          const TopoDS_Shape& aShape)
{
  MESSAGE("NETGENPlugin_NETGEN_2D_ONLY::Compute()");

  SMESHDS_Mesh* meshDS = aMesh.GetMeshDS();
  int faceID = meshDS->ShapeToIndex( aShape );

  SMESH_MesherHelper helper(aMesh);
  _quadraticMesh = helper.IsQuadraticSubMesh(aShape);
  helper.SetElementsOnShape( true );
  const bool ignoreMediumNodes = _quadraticMesh;
  
  // ------------------------
  // get all edges of a face
  // ------------------------
  const TopoDS_Face F = TopoDS::Face( aShape.Oriented( TopAbs_FORWARD ));
  TError problem;
  TSideVector wires = StdMeshers_FaceSide::GetFaceWires( F, aMesh, ignoreMediumNodes, problem );
  if ( problem && !problem->IsOK() )
    return error( problem );
  int nbWires = wires.size();
  if ( nbWires == 0 )
    return error( "Problem in StdMeshers_FaceSide::GetFaceWires()");
  if ( wires[0]->NbSegments() < 3 ) // ex: a circle with 2 segments
    return error(COMPERR_BAD_INPUT_MESH,
                 SMESH_Comment("Too few segments: ")<<wires[0]->NbSegments());

  // -------------------------
  // Make input netgen mesh
  // -------------------------

  Ng_Init();
  netgen::Mesh * ngMesh = new netgen::Mesh ();

  netgen::OCCGeometry occgeo;
  NETGENPlugin_Mesher::PrepareOCCgeometry( occgeo, F, aMesh );
  occgeo.fmap.Clear(); // face can be reversed, which is wrong in this case (issue 19978)
  occgeo.fmap.Add( F );

  vector< const SMDS_MeshNode* > nodeVec;
  problem = AddSegmentsToMesh( *ngMesh, occgeo, wires, helper, nodeVec );
  if ( problem && !problem->IsOK() ) {
    delete ngMesh; Ng_Exit();
    return error( problem );
  }

  // --------------------
  // compute edge length
  // --------------------

  double edgeLength = 0;
  if (_hypLengthFromEdges || (!_hypLengthFromEdges && !_hypMaxElementArea))
  {
    int nbSegments = 0;
    for ( int iW = 0; iW < nbWires; ++iW )
    {
      edgeLength += wires[ iW ]->Length();
      nbSegments += wires[ iW ]->NbSegments();
    }
    if ( nbSegments )
      edgeLength /= nbSegments;
  }
  if ( _hypMaxElementArea )
  {
    double maxArea = _hypMaxElementArea->GetMaxArea();
    edgeLength = sqrt(2. * maxArea/sqrt(3.0));
  }
  if ( edgeLength < DBL_MIN )
    edgeLength = occgeo.GetBoundingBox().Diam();

  //cout << " edgeLength = " << edgeLength << endl;

  netgen::mparam.maxh = edgeLength;
  netgen::mparam.quad = _hypQuadranglePreference ? 1 : 0;
  //ngMesh->SetGlobalH ( edgeLength );

  // -------------------------
  // Generate surface mesh
  // -------------------------

  char *optstr = 0;
  int startWith = MESHCONST_MESHSURFACE;
  int endWith   = MESHCONST_OPTSURFACE;
  int err = 1;

  try {
#if (OCC_VERSION_MAJOR << 16 | OCC_VERSION_MINOR << 8 | OCC_VERSION_MAINTENANCE) > 0x060100
    OCC_CATCH_SIGNALS;
#endif
#ifdef NETGEN_V5
    err = netgen::OCCGenerateMesh(occgeo, ngMesh,netgen::mparam, startWith, endWith);
#else
    err = netgen::OCCGenerateMesh(occgeo, ngMesh, startWith, endWith, optstr);
#endif
  }
  catch (Standard_Failure& ex) {
    string comment = ex.DynamicType()->Name();
    if ( ex.GetMessageString() && strlen( ex.GetMessageString() )) {
      comment += ": ";
//.........这里部分代码省略.........

//.........这里部分代码省略.........
      // ------------------------------------------------
      // sort list of sub-meshes according to mesh order
      // ------------------------------------------------
      smVec.assign( smWithAlgoSupportingSubmeshes[ aShapeDim ].begin(),
                    smWithAlgoSupportingSubmeshes[ aShapeDim ].end() );
      aMesh.SortByMeshOrder( smVec );

      // ------------------------------------------------------------
      // compute sub-meshes with local uni-dimensional algos under
      // sub-meshes with all-dimensional algos
      // ------------------------------------------------------------
      // start from lower shapes
      for ( size_t i = 0; i < smVec.size(); ++i )
      {
        sm = smVec[i];

        // get a shape the algo is assigned to
        if ( !GetAlgo( sm, & algoShape ))
          continue; // strange...

        // look for more local algos
        smIt = sm->getDependsOnIterator(!includeSelf, !complexShapeFirst);
        while ( smIt->more() )
        {
          SMESH_subMesh* smToCompute = smIt->next();

          const TopoDS_Shape& aSubShape = smToCompute->GetSubShape();
          const int aShapeDim = GetShapeDim( aSubShape );
          //if ( aSubShape.ShapeType() == TopAbs_VERTEX ) continue;
          if ( aShapeDim < 1 ) continue;

          // check for preview dimension limitations
          if ( aShapesId && GetShapeDim( aSubShape.ShapeType() ) > (int)aDim )
            continue;

          SMESH_HypoFilter filter( SMESH_HypoFilter::IsAlgo() );
          filter
            .And( SMESH_HypoFilter::IsApplicableTo( aSubShape ))
            .And( SMESH_HypoFilter::IsMoreLocalThan( algoShape, aMesh ));

          if ( SMESH_Algo* subAlgo = (SMESH_Algo*) aMesh.GetHypothesis( smToCompute, filter, true))
          {
            if ( ! subAlgo->NeedDiscreteBoundary() ) continue;
            SMESH_Hypothesis::Hypothesis_Status status;
            if ( subAlgo->CheckHypothesis( aMesh, aSubShape, status ))
              // mesh a lower smToCompute starting from vertices
              Compute( aMesh, aSubShape, aShapeOnly, /*anUpward=*/true, aDim, aShapesId );
          }
        }
      }
      // --------------------------------
      // apply the all-dimensional algos
      // --------------------------------
      for ( size_t i = 0; i < smVec.size(); ++i )
      {
        sm = smVec[i];
        if ( sm->GetComputeState() == SMESH_subMesh::READY_TO_COMPUTE)
        {
          const TopAbs_ShapeEnum shapeType = sm->GetSubShape().ShapeType();
          // check for preview dimension limitations
          if ( aShapesId && GetShapeDim( shapeType ) > (int)aDim )
            continue;

          if (_compute_canceled)
            return false;
          setCurrentSubMesh( sm );
          sm->ComputeStateEngine( computeEvent );
          setCurrentSubMesh( NULL );
          if ( aShapesId )
            aShapesId->insert( sm->GetId() );
        }
      }
    } // loop on shape dimensions

    // -----------------------------------------------
    // mesh the rest sub-shapes starting from vertices
    // -----------------------------------------------
    ret = Compute( aMesh, aShape, aShapeOnly, /*anUpward=*/true, aDim, aShapesId );
  }

  MESSAGE( "VSR - SMESH_Gen::Compute() finished, OK = " << ret);
  MEMOSTAT;

  SMESHDS_Mesh *myMesh = aMesh.GetMeshDS();
  MESSAGE("*** compactMesh after compute");
  myMesh->compactMesh();

  // fix quadratic mesh by bending iternal links near concave boundary
  if ( aShape.IsSame( aMesh.GetShapeToMesh() ) &&
       !aShapesId && // not preview
       ret ) // everything is OK
  {
    SMESH_MesherHelper aHelper( aMesh );
    if ( aHelper.IsQuadraticMesh() != SMESH_MesherHelper::LINEAR )
    {
      aHelper.FixQuadraticElements( sm->GetComputeError() );
    }
  }
  return ret;
}

/*!
 * Special function for search or creation medium node
 */
const SMDS_MeshNode* SMESH_MesherHelper::GetMediumNode(const SMDS_MeshNode* n1,
                                                       const SMDS_MeshNode* n2,
                                                       bool force3d)
{
  TopAbs_ShapeEnum shapeType = myShape.IsNull() ? TopAbs_SHAPE : myShape.ShapeType();

  NLink link(( n1 < n2 ? n1 : n2 ), ( n1 < n2 ? n2 : n1 ));
  ItNLinkNode itLN = myNLinkNodeMap.find( link );
  if ( itLN != myNLinkNodeMap.end() ) {
    return (*itLN).second;
  }
  else {
    // create medium node
    SMDS_MeshNode* n12;
    SMESHDS_Mesh* meshDS = GetMeshDS();
    int faceID = -1, edgeID = -1;
    const SMDS_PositionPtr Pos1 = n1->GetPosition();
    const SMDS_PositionPtr Pos2 = n2->GetPosition();
  
    if( myShape.IsNull() )
    {
      if( Pos1->GetTypeOfPosition()==SMDS_TOP_FACE ) {
        faceID = Pos1->GetShapeId();
      }
      else if( Pos2->GetTypeOfPosition()==SMDS_TOP_FACE ) {
        faceID = Pos2->GetShapeId();
      }

      if( Pos1->GetTypeOfPosition()==SMDS_TOP_EDGE ) {
        edgeID = Pos1->GetShapeId();
      }
      if( Pos2->GetTypeOfPosition()==SMDS_TOP_EDGE ) {
        edgeID = Pos2->GetShapeId();
      }
    }

    if(!force3d) {
      // we try to create medium node using UV parameters of
      // nodes, else - medium between corresponding 3d points
      if(faceID>-1 || shapeType == TopAbs_FACE) {
	// obtaining a face and 2d points for nodes
	TopoDS_Face F;
	if( myShape.IsNull() )
          F = TopoDS::Face(meshDS->IndexToShape(faceID));
	else {
          F = TopoDS::Face(myShape);
          faceID = myShapeID;
        }

	gp_XY p1 = GetNodeUV(F,n1,n2);
        gp_XY p2 = GetNodeUV(F,n2,n1);

	if ( IsDegenShape( Pos1->GetShapeId() ))
	  p1.SetCoord( myParIndex, p2.Coord( myParIndex ));
	else if ( IsDegenShape( Pos2->GetShapeId() ))
	  p2.SetCoord( myParIndex, p1.Coord( myParIndex ));

	//checking if surface is periodic
	Handle(Geom_Surface) S = BRep_Tool::Surface(F);
	Standard_Real UF,UL,VF,VL;
	S->Bounds(UF,UL,VF,VL);

	Standard_Real u,v;
	Standard_Boolean isUPeriodic = S->IsUPeriodic();
	if(isUPeriodic) {
	  Standard_Real UPeriod = S->UPeriod();
	  Standard_Real p2x = p2.X()+ShapeAnalysis::AdjustByPeriod(p2.X(),p1.X(),UPeriod);
	  Standard_Real pmid = (p1.X()+p2x)/2.;
	  u = pmid+ShapeAnalysis::AdjustToPeriod(pmid,UF,UL);
	}
	else 
	  u= (p1.X()+p2.X())/2.;

	Standard_Boolean isVPeriodic = S->IsVPeriodic();
	if(isVPeriodic) {
	  Standard_Real VPeriod = S->VPeriod();
	  Standard_Real p2y = p2.Y()+ShapeAnalysis::AdjustByPeriod(p2.Y(),p1.Y(),VPeriod);
	  Standard_Real pmid = (p1.Y()+p2y)/2.;
	  v = pmid+ShapeAnalysis::AdjustToPeriod(pmid,VF,VL);
	}
	else
	  v = (p1.Y()+p2.Y())/2.;

        gp_Pnt P = S->Value(u, v);
        n12 = meshDS->AddNode(P.X(), P.Y(), P.Z());
        meshDS->SetNodeOnFace(n12, faceID, u, v);
        myNLinkNodeMap.insert(NLinkNodeMap::value_type(link,n12));
        return n12;
      }
      if (edgeID>-1 || shapeType == TopAbs_EDGE) {

	TopoDS_Edge E;
	if( myShape.IsNull() )
          E = TopoDS::Edge(meshDS->IndexToShape(edgeID));
	else {
          E = TopoDS::Edge(myShape);
          edgeID = myShapeID;
//.........这里部分代码省略.........

bool StdMeshers_CompositeSegment_1D::Compute(SMESH_Mesh &         aMesh,
                                             const TopoDS_Shape & aShape)
{
  TopoDS_Edge edge = TopoDS::Edge( aShape );
  SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();

  // Get edges to be discretized as a whole
  TopoDS_Face nullFace;
  auto_ptr< StdMeshers_FaceSide > side( GetFaceSide(aMesh, edge, nullFace, true ));
  //side->dump("IN COMPOSITE SEG");

  if ( side->NbEdges() < 2 )
    return StdMeshers_Regular_1D::Compute( aMesh, aShape );

  // update segment lenght computed by StdMeshers_AutomaticLength
  const list <const SMESHDS_Hypothesis * > & hyps = GetUsedHypothesis(aMesh, aShape);
  if ( !hyps.empty() ) {
    StdMeshers_AutomaticLength * autoLenHyp = const_cast<StdMeshers_AutomaticLength *>
      (dynamic_cast <const StdMeshers_AutomaticLength * >(hyps.front()));
    if ( autoLenHyp )
      _value[ BEG_LENGTH_IND ]= autoLenHyp->GetLength( &aMesh, side->Length() );
  }

  // Compute node parameters
  auto_ptr< BRepAdaptor_CompCurve > C3d ( side->GetCurve3d() );
  double f = C3d->FirstParameter(), l = C3d->LastParameter();
  list< double > params;
  if ( !computeInternalParameters ( aMesh, *C3d, side->Length(), f, l, params, false ))
    return false;

  // Redistribute parameters near ends
  TopoDS_Vertex VFirst = side->FirstVertex();
  TopoDS_Vertex VLast  = side->LastVertex();
  redistributeNearVertices( aMesh, *C3d, side->Length(), params, VFirst, VLast );

  params.push_front(f);
  params.push_back(l);
  int nbNodes = params.size();

  // Create mesh

  const SMDS_MeshNode * nFirst = SMESH_Algo::VertexNode( VFirst, meshDS );
  const SMDS_MeshNode * nLast  = SMESH_Algo::VertexNode( VLast, meshDS );
  if (!nFirst)
    return error(COMPERR_BAD_INPUT_MESH, TComm("No node on vertex ")
                 <<meshDS->ShapeToIndex(VFirst));
  if (!nLast)
    return error(COMPERR_BAD_INPUT_MESH, TComm("No node on vertex ")
                 <<meshDS->ShapeToIndex(VLast));

  vector<const SMDS_MeshNode*> nodes( nbNodes, (const SMDS_MeshNode*)0 );
  nodes.front() = nFirst;
  nodes.back()  = nLast;

  // create internal nodes
  list< double >::iterator parIt = params.begin();
  double prevPar = *parIt;
  Standard_Real u;
  for ( int iN = 0; parIt != params.end(); ++iN, ++parIt)
  {
    if ( !nodes[ iN ] ) {
      gp_Pnt p = C3d->Value( *parIt );
      SMDS_MeshNode* n = meshDS->AddNode( p.X(), p.Y(), p.Z());
      C3d->Edge( *parIt, edge, u );
      meshDS->SetNodeOnEdge( n, edge, u );
//       cout << "new NODE: par="<<*parIt<<" ePar="<<u<<" e="<<edge.TShape().operator->()
//            << " " << n << endl;
      nodes[ iN ] = n;
    }
    // create edges
    if ( iN ) {
      double mPar = ( prevPar + *parIt )/2;
      if ( _quadraticMesh ) {
        // create medium node
        double segLen = GCPnts_AbscissaPoint::Length(*C3d, prevPar, *parIt);
        GCPnts_AbscissaPoint ruler( *C3d, segLen/2., prevPar );
        if ( ruler.IsDone() )
          mPar = ruler.Parameter();
        gp_Pnt p = C3d->Value( mPar );
        SMDS_MeshNode* n = meshDS->AddNode( p.X(), p.Y(), p.Z());
        //cout << "new NODE "<< n << endl;
        meshDS->SetNodeOnEdge( n, edge, u );
        SMDS_MeshEdge * seg = meshDS->AddEdge(nodes[ iN-1 ], nodes[ iN ], n);
        meshDS->SetMeshElementOnShape(seg, edge);
      }
      else {
        C3d->Edge( mPar, edge, u );
        SMDS_MeshEdge * seg = meshDS->AddEdge(nodes[ iN-1 ], nodes[ iN ]);
        meshDS->SetMeshElementOnShape(seg, edge);
      }
    }
    prevPar = *parIt;
  }

  // remove nodes on internal vertices
  for ( int iE = 1; iE < side->NbEdges(); ++iE )
  {
    TopoDS_Vertex V = side->FirstVertex( iE );
    while ( const SMDS_MeshNode * n = SMESH_Algo::VertexNode( V, meshDS ))
      meshDS->RemoveNode( n );
//.........这里部分代码省略.........

gp_XY SMESH_MesherHelper::GetNodeUV(const TopoDS_Face&   F,
                                    const SMDS_MeshNode* n,
                                    const SMDS_MeshNode* n2) const
{
  gp_Pnt2d uv( 1e100, 1e100 );
  const SMDS_PositionPtr Pos = n->GetPosition();
  if(Pos->GetTypeOfPosition()==SMDS_TOP_FACE)
  {
    // node has position on face
    const SMDS_FacePosition* fpos =
      static_cast<const SMDS_FacePosition*>(n->GetPosition().get());
    uv = gp_Pnt2d(fpos->GetUParameter(),fpos->GetVParameter());
  }
  else if(Pos->GetTypeOfPosition()==SMDS_TOP_EDGE)
  {
    // node has position on edge => it is needed to find
    // corresponding edge from face, get pcurve for this
    // edge and recieve value from this pcurve
    const SMDS_EdgePosition* epos =
      static_cast<const SMDS_EdgePosition*>(n->GetPosition().get());
    SMESHDS_Mesh* meshDS = GetMeshDS();
    int edgeID = Pos->GetShapeId();
    TopoDS_Edge E = TopoDS::Edge(meshDS->IndexToShape(edgeID));
    double f, l;
    Handle(Geom2d_Curve) C2d = BRep_Tool::CurveOnSurface(E, F, f, l);
    uv = C2d->Value( epos->GetUParameter() );
    // for a node on a seam edge select one of UVs on 2 pcurves
    if ( n2 && IsSeamShape( edgeID ) )
      uv = GetUVOnSeam( uv, GetNodeUV( F, n2, 0 ));
  }
  else if(Pos->GetTypeOfPosition()==SMDS_TOP_VERTEX)
  {
    if ( int vertexID = n->GetPosition()->GetShapeId() ) {
      bool ok = true;
      const TopoDS_Vertex& V = TopoDS::Vertex(GetMeshDS()->IndexToShape(vertexID));
      try {
        uv = BRep_Tool::Parameters( V, F );
      }
      catch (Standard_Failure& exc) {
        ok = false;
      }
      if ( !ok ) {
        for ( TopExp_Explorer vert(F,TopAbs_VERTEX); !ok && vert.More(); vert.Next() )
          ok = ( V == vert.Current() );
        if ( !ok ) {
#ifdef _DEBUG_
          MESSAGE ( "SMESH_MesherHelper::GetNodeUV(); Vertex " << vertexID
               << " not in face " << GetMeshDS()->ShapeToIndex( F ) );
#endif
          // get UV of a vertex closest to the node
          double dist = 1e100;
          gp_Pnt pn ( n->X(),n->Y(),n->Z() );
          for ( TopExp_Explorer vert(F,TopAbs_VERTEX); !ok && vert.More(); vert.Next() ) {
            TopoDS_Vertex curV = TopoDS::Vertex( vert.Current() );
            gp_Pnt p = BRep_Tool::Pnt( curV );
            double curDist = p.SquareDistance( pn );
            if ( curDist < dist ) {
              dist = curDist;
              uv = BRep_Tool::Parameters( curV, F );
              if ( dist < DBL_MIN ) break;
            }
          }
        }
        else {
          TopTools_ListIteratorOfListOfShape it( myMesh->GetAncestors( V ));
          for ( ; it.More(); it.Next() ) {
            if ( it.Value().ShapeType() == TopAbs_EDGE ) {
              const TopoDS_Edge & edge = TopoDS::Edge( it.Value() );
              double f,l;
              Handle(Geom2d_Curve) C2d = BRep_Tool::CurveOnSurface(edge, F, f, l);
              if ( !C2d.IsNull() ) {
                double u = ( V == TopExp::FirstVertex( edge ) ) ?  f : l;
                uv = C2d->Value( u );
                break;
              }
            }
          }
        }
      }
      if ( n2 && IsSeamShape( vertexID ) )
        uv = GetUVOnSeam( uv, GetNodeUV( F, n2, 0 ));
    }
  }
  return uv.XY();
}