PyObject* FemMeshPy::addEdge(PyObject *args)
{
    int n1,n2;
    if (!PyArg_ParseTuple(args, "ii",&n1,&n2))
        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);
        if (!node1 || !node2)
            throw std::runtime_error("Failed to get node of the given indices");
        SMDS_MeshEdge* edge = meshDS->AddEdge(node1, node2);
        if (!edge)
            throw std::runtime_error("Failed to add edge");
        return Py::new_reference_to(Py::Int(edge->GetID()));
    }
    catch (const std::exception& e) {
        PyErr_SetString(Base::BaseExceptionFreeCADError, e.what());
        return 0;
    }
}

PyObject* FemMeshPy::addFace(PyObject *args)
{
    int n1,n2,n3;
    if (!PyArg_ParseTuple(args, "iii",&n1,&n2,&n3))
        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);
        if (!node1 || !node2 || !node3)
            throw std::runtime_error("Failed to get node of the given indices");
        SMDS_MeshFace* face = meshDS->AddFace(node1, node2, node3);
        if (!face)
            throw std::runtime_error("Failed to add face");
        return Py::new_reference_to(Py::Int(face->GetID()));
    }
    catch (const std::exception& e) {
        PyErr_SetString(PyExc_Exception, e.what());
        return 0;
    }
}

TSideVector StdMeshers_FaceSide::GetFaceWires(const TopoDS_Face& theFace,
                                              SMESH_Mesh &       theMesh,
                                              const bool         theIgnoreMediumNodes,
                                              TError &           theError)
{
  TopoDS_Vertex V1;
  list< TopoDS_Edge > edges;
  list< int > nbEdgesInWires;
  int nbWires = SMESH_Block::GetOrderedEdges (theFace, V1, edges, nbEdgesInWires);

  // split list of all edges into separate wires
  TSideVector wires( nbWires );
  list< int >::iterator nbE = nbEdgesInWires.begin();
  list< TopoDS_Edge >::iterator from, to;
  from = to = edges.begin();
  for ( int iW = 0; iW < nbWires; ++iW )
  {
    std::advance( to, *nbE++ );
    list< TopoDS_Edge > wireEdges( from, to );
    // assure that there is a node on the first vertex
    // as StdMeshers_FaceSide::GetUVPtStruct() requires
    while ( !SMESH_Algo::VertexNode( TopExp::FirstVertex( wireEdges.front(), true),
                                     theMesh.GetMeshDS()))
    {
      wireEdges.splice(wireEdges.end(), wireEdges,
                       wireEdges.begin(), ++wireEdges.begin());
      if ( from->IsSame( wireEdges.front() )) {
        theError = TError
          ( new SMESH_ComputeError(COMPERR_BAD_INPUT_MESH,"No nodes on vertices"));
        return TSideVector(0);
      }
    }
    // find out side orientation, which is important if there are several wires (PAL19080) 
    bool isForward = true;
    if ( nbWires > 1 ) {
      TopExp_Explorer e( theFace, TopAbs_EDGE );
      while ( ! e.Current().IsSame( wireEdges.back() ))
        e.Next();
      isForward = ( e.Current().Orientation() == wireEdges.back().Orientation() );
    }

    StdMeshers_FaceSide* wire = new StdMeshers_FaceSide( theFace, wireEdges, &theMesh,
                                                         isForward, theIgnoreMediumNodes);
    wires[ iW ] = StdMeshers_FaceSidePtr( wire );
    from = to;
  }
  return wires;
}

bool EvaluatePentahedralMesh(SMESH_Mesh & aMesh,
                             const TopoDS_Shape & aShape,
                             MapShapeNbElems& aResMap)
{
  StdMeshers_Penta_3D anAlgo;
  bool bOK = anAlgo.Evaluate(aMesh, aShape, aResMap);

  //err = anAlgo.GetComputeError();
  //if ( !bOK && anAlgo.ErrorStatus() == 5 )
  if( !bOK ) {
    static StdMeshers_Prism_3D * aPrism3D = 0;
    if ( !aPrism3D ) {
      SMESH_Gen* gen = aMesh.GetGen();
      aPrism3D = new StdMeshers_Prism_3D( gen->GetANewId(), 0, gen );
    }
    SMESH_Hypothesis::Hypothesis_Status aStatus;
    if ( aPrism3D->CheckHypothesis( aMesh, aShape, aStatus ) ) {
      return aPrism3D->Evaluate(aMesh, aShape, aResMap);
    }
  }

  return bOK;
}

bool NETGENPlugin_NETGEN_3D::Evaluate(SMESH_Mesh& aMesh,
                                      const TopoDS_Shape& aShape,
                                      MapShapeNbElems& aResMap)
{
    int nbtri = 0, nbqua = 0;
    double fullArea = 0.0;
    for (TopExp_Explorer exp(aShape, TopAbs_FACE); exp.More(); exp.Next()) {
        TopoDS_Face F = TopoDS::Face( exp.Current() );
        SMESH_subMesh *sm = aMesh.GetSubMesh(F);
        MapShapeNbElemsItr anIt = aResMap.find(sm);
        if( anIt==aResMap.end() ) {
            SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
            smError.reset( new SMESH_ComputeError(COMPERR_ALGO_FAILED,"Submesh can not be evaluated",this));
            return false;
        }
        std::vector<int> aVec = (*anIt).second;
        nbtri += Max(aVec[SMDSEntity_Triangle],aVec[SMDSEntity_Quad_Triangle]);
        nbqua += Max(aVec[SMDSEntity_Quadrangle],aVec[SMDSEntity_Quad_Quadrangle]);
        GProp_GProps G;
        BRepGProp::SurfaceProperties(F,G);
        double anArea = G.Mass();
        fullArea += anArea;
    }

    // collect info from edges
    int nb0d_e = 0, nb1d_e = 0;
    bool IsQuadratic = false;
    bool IsFirst = true;
    TopTools_MapOfShape tmpMap;
    for (TopExp_Explorer exp(aShape, TopAbs_EDGE); exp.More(); exp.Next()) {
        TopoDS_Edge E = TopoDS::Edge(exp.Current());
        if( tmpMap.Contains(E) )
            continue;
        tmpMap.Add(E);
        SMESH_subMesh *aSubMesh = aMesh.GetSubMesh(exp.Current());
        MapShapeNbElemsItr anIt = aResMap.find(aSubMesh);
        if( anIt==aResMap.end() ) {
            SMESH_ComputeErrorPtr& smError = aSubMesh->GetComputeError();
            smError.reset( new SMESH_ComputeError(COMPERR_ALGO_FAILED,
                                                  "Submesh can not be evaluated",this));
            return false;
        }
        std::vector<int> aVec = (*anIt).second;
        nb0d_e += aVec[SMDSEntity_Node];
        nb1d_e += Max(aVec[SMDSEntity_Edge],aVec[SMDSEntity_Quad_Edge]);
        if(IsFirst) {
            IsQuadratic = (aVec[SMDSEntity_Quad_Edge] > aVec[SMDSEntity_Edge]);
            IsFirst = false;
        }
    }
    tmpMap.Clear();

    double ELen_face = sqrt(2.* ( fullArea/(nbtri+nbqua*2) ) / sqrt(3.0) );
    double ELen_vol = pow( 72, 1/6. ) * pow( _maxElementVolume, 1/3. );
    double ELen = Min(ELen_vol,ELen_face*2);

    GProp_GProps G;
    BRepGProp::VolumeProperties(aShape,G);
    double aVolume = G.Mass();
    double tetrVol = 0.1179*ELen*ELen*ELen;
    double CoeffQuality = 0.9;
    int nbVols = (int)aVolume/tetrVol/CoeffQuality;
    int nb1d_f = (nbtri*3 + nbqua*4 - nb1d_e) / 2;
    int nb1d_in = (int) ( nbVols*6 - nb1d_e - nb1d_f ) / 5;
    std::vector<int> aVec(SMDSEntity_Last);
    for(int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aVec[i]=0;
    if( IsQuadratic ) {
        aVec[SMDSEntity_Node] = nb1d_in/6 + 1 + nb1d_in;
        aVec[SMDSEntity_Quad_Tetra] = nbVols - nbqua*2;
        aVec[SMDSEntity_Quad_Pyramid] = nbqua;
    }
    else {
        aVec[SMDSEntity_Node] = nb1d_in/6 + 1;
        aVec[SMDSEntity_Tetra] = nbVols - nbqua*2;
        aVec[SMDSEntity_Pyramid] = nbqua;
    }
    SMESH_subMesh *sm = aMesh.GetSubMesh(aShape);
    aResMap.insert(std::make_pair(sm,aVec));

    return true;
}

bool NETGENPlugin_NETGEN_3D::compute(SMESH_Mesh&                     aMesh,
                                     SMESH_MesherHelper&             helper,
                                     vector< const SMDS_MeshNode* >& nodeVec,
                                     Ng_Mesh *                       Netgen_mesh)
{
  netgen::multithread.terminate = 0;

  netgen::Mesh* ngMesh = (netgen::Mesh*)Netgen_mesh;
  int Netgen_NbOfNodes = Ng_GetNP(Netgen_mesh);

#if NETGEN_VERSION < 5
  char *optstr = 0;
#endif
  int startWith = netgen::MESHCONST_MESHVOLUME;
  int endWith   = netgen::MESHCONST_OPTVOLUME;
  int err = 1;

  NETGENPlugin_Mesher aMesher( &aMesh, helper.GetSubShape(), /*isVolume=*/true );
  netgen::OCCGeometry occgeo;
  
  if ( _hypParameters )
  {
    aMesher.SetParameters( _hypParameters );
    if ( !_hypParameters->GetOptimize() )
      endWith = netgen::MESHCONST_MESHVOLUME;
  }
  else if ( _hypMaxElementVolume )
  {
    netgen::mparam.maxh = pow( 72, 1/6. ) * pow( _maxElementVolume, 1/3. );
    // limitVolumeSize( ngMesh, netgen::mparam.maxh ); // result is unpredictable
  }
  else if ( aMesh.HasShapeToMesh() )
  {
    aMesher.PrepareOCCgeometry( occgeo, helper.GetSubShape(), aMesh );
    netgen::mparam.maxh = occgeo.GetBoundingBox().Diam()/2;
  }
  else
  {
    netgen::Point3d pmin, pmax;
    ngMesh->GetBox (pmin, pmax);
    netgen::mparam.maxh = Dist(pmin, pmax)/2;
  }

  if ( !_hypParameters && aMesh.HasShapeToMesh() )
  {
    netgen::mparam.minh = aMesher.GetDefaultMinSize( helper.GetSubShape(), netgen::mparam.maxh );
  }

  try
  {
    OCC_CATCH_SIGNALS;

#if NETGEN_VERSION > 4
    ngMesh->CalcLocalH(netgen::mparam.grading);
    err = netgen::OCCGenerateMesh(occgeo, ngMesh, netgen::mparam, startWith, endWith);
#else
    ngMesh->CalcLocalH();
    err = netgen::OCCGenerateMesh(occgeo, ngMesh, startWith, endWith, optstr);
#endif
    if(netgen::multithread.terminate)
      return false;
    if ( err )
      error(SMESH_Comment("Error in netgen::OCCGenerateMesh() at ") << netgen::multithread.task);
  }
  catch (Standard_Failure& ex)
  {
    SMESH_Comment str("Exception in  netgen::OCCGenerateMesh()");
    str << " at " << netgen::multithread.task
        << ": " << ex.DynamicType()->Name();
    if ( ex.GetMessageString() && strlen( ex.GetMessageString() ))
      str << ": " << ex.GetMessageString();
    error(str);
  }
  catch (netgen::NgException exc)
  {
    SMESH_Comment str("NgException");
    if ( strlen( netgen::multithread.task ) > 0 )
      str << " at " << netgen::multithread.task;
    str << ": " << exc.What();
    error(str);
  }
  catch (...)
  {
    SMESH_Comment str("Exception in  netgen::OCCGenerateMesh()");
    if ( strlen( netgen::multithread.task ) > 0 )
      str << " at " << netgen::multithread.task;
    error(str);
  }

  int Netgen_NbOfNodesNew = Ng_GetNP(Netgen_mesh);
  int Netgen_NbOfTetra    = Ng_GetNE(Netgen_mesh);

  MESSAGE("End of Volume Mesh Generation. err=" << err <<
          ", nb new nodes: " << Netgen_NbOfNodesNew - Netgen_NbOfNodes <<
          ", nb tetra: " << Netgen_NbOfTetra);

  // -------------------------------------------------------------------
  // Feed back the SMESHDS with the generated Nodes and Volume Elements
  // -------------------------------------------------------------------

//.........这里部分代码省略.........

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 += ": ";
//.........这里部分代码省略.........

PyObject* FemMeshPy::addFace(PyObject *args)
{
    SMESH_Mesh* mesh = getFemMeshPtr()->getSMesh();
    SMESHDS_Mesh* meshDS = mesh->GetMeshDS();

    int n1,n2,n3;
    if (PyArg_ParseTuple(args, "iii",&n1,&n2,&n3))
    {
        // old form, deprecated
        try {
            const SMDS_MeshNode* node1 = meshDS->FindNode(n1);
            const SMDS_MeshNode* node2 = meshDS->FindNode(n2);
            const SMDS_MeshNode* node3 = meshDS->FindNode(n3);
            if (!node1 || !node2 || !node3)
                throw std::runtime_error("Failed to get node of the given indices");
            SMDS_MeshFace* face = meshDS->AddFace(node1, node2, node3);
            if (!face)
                throw std::runtime_error("Failed to add face");
            return Py::new_reference_to(Py::Long(face->GetID()));
        }
        catch (const std::exception& e) {
            PyErr_SetString(Base::BaseExceptionFreeCADError, e.what());
            return 0;
        }
    }
    PyErr_Clear();

    PyObject *obj;
    int ElementId=-1;
    if (PyArg_ParseTuple(args, "O!|i", &PyList_Type, &obj, &ElementId))
    {
        Py::List list(obj);
        std::vector<const SMDS_MeshNode*> Nodes;
        for (Py::List::iterator it = list.begin(); it != list.end(); ++it) {
#if PY_MAJOR_VERSION >= 3
            Py::Long NoNr(*it);
#else
            Py::Int NoNr(*it);
#endif
            const SMDS_MeshNode* node = meshDS->FindNode(NoNr);
            if (!node)
                throw std::runtime_error("Failed to get node of the given indices");
            Nodes.push_back(node);
        }

        SMDS_MeshFace* face=0;
        if(ElementId != -1) {
            switch(Nodes.size()){
                case 3:
                    face = meshDS->AddFaceWithID(Nodes[0],Nodes[1],Nodes[2],ElementId);
                    if (!face)
                        throw std::runtime_error("Failed to add triangular face with given ElementId");
                    break;
                case 4:
                    face = meshDS->AddFaceWithID(Nodes[0],Nodes[1],Nodes[2],Nodes[3],ElementId);
                    if (!face)
                        throw std::runtime_error("Failed to add face with given ElementId");
                    break;
                case 6:
                    face = meshDS->AddFaceWithID(Nodes[0],Nodes[1],Nodes[2],Nodes[3],Nodes[4],Nodes[5],ElementId);
                    if (!face)
                        throw std::runtime_error("Failed to add face with given ElementId");
                    break;
                case 8:
                    face = meshDS->AddFaceWithID(Nodes[0],Nodes[1],Nodes[2],Nodes[3],Nodes[4],Nodes[5],Nodes[6],Nodes[7],ElementId);
                    if (!face)
                        throw std::runtime_error("Failed to add face with given ElementId");
                    break;
                default:
                    throw std::runtime_error("Unknown node count, [3|4|6|8] are allowed"); //unknown face type
            }
        }else{
            switch(Nodes.size()){
                case 3:
                    face = meshDS->AddFace(Nodes[0],Nodes[1],Nodes[2]);
                    if (!face)
                        throw std::runtime_error("Failed to add triangular face");
                    break;
                case 4:
                    face = meshDS->AddFace(Nodes[0],Nodes[1],Nodes[2],Nodes[3]);
                    if (!face)
                        throw std::runtime_error("Failed to add face");
                    break;
                case 6:
                    face = meshDS->AddFace(Nodes[0],Nodes[1],Nodes[2],Nodes[3],Nodes[4],Nodes[5]);
                    if (!face)
                        throw std::runtime_error("Failed to add face");
                    break;
                case 8:
                    face = meshDS->AddFace(Nodes[0],Nodes[1],Nodes[2],Nodes[3],Nodes[4],Nodes[5],Nodes[6],Nodes[7]);
                    if (!face)
                        throw std::runtime_error("Failed to add face");
                    break;
                default:
                    throw std::runtime_error("Unknown node count, [4|5|6|8] are allowed"); //unknown face type
            }
        }

        return Py::new_reference_to(Py::Long(face->GetID()));

//.........这里部分代码省略.........

bool SMESH_Gen::IsGlobalHypothesis(const SMESH_Hypothesis* theHyp, SMESH_Mesh& aMesh)
{
  SMESH_HypoFilter filter( SMESH_HypoFilter::Is( theHyp ));
  return aMesh.GetHypothesis( aMesh.GetMeshDS()->ShapeToMesh(), filter, false );
}

bool StdMeshers_RadialPrism_3D::Compute(SMESH_Mesh& aMesh, const TopoDS_Shape& aShape)
{
  TopExp_Explorer exp;
  SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();

  myHelper = new SMESH_MesherHelper( aMesh );
  myHelper->IsQuadraticSubMesh( aShape );
  // to delete helper at exit from Compute()
  std::auto_ptr<SMESH_MesherHelper> helperDeleter( myHelper );

  // get 2 shells
  TopoDS_Solid solid = TopoDS::Solid( aShape );
  TopoDS_Shell outerShell = BRepClass3d::OuterShell( solid );
  TopoDS_Shape innerShell;
  int nbShells = 0;
  for ( TopoDS_Iterator It (solid); It.More(); It.Next(), ++nbShells )
    if ( !outerShell.IsSame( It.Value() ))
      innerShell = It.Value();
  if ( nbShells != 2 )
    return error(COMPERR_BAD_SHAPE, SMESH_Comment("Must be 2 shells but not ")<<nbShells);

  // ----------------------------------
  // Associate subshapes of the shells
  // ----------------------------------

  TAssocTool::TShapeShapeMap shape2ShapeMap;
  if ( !TAssocTool::FindSubShapeAssociation( outerShell, &aMesh,
                                             innerShell, &aMesh,
                                             shape2ShapeMap) )
    return error(COMPERR_BAD_SHAPE,"Topology of inner and outer shells seems different" );

  // ------------------
  // Make mesh
  // ------------------

  TNode2ColumnMap node2columnMap;
  myLayerPositions.clear();

  for ( exp.Init( outerShell, TopAbs_FACE ); exp.More(); exp.Next() )
  {
    // Corresponding subshapes
    TopoDS_Face outFace = TopoDS::Face( exp.Current() );
    TopoDS_Face inFace;
    if ( !shape2ShapeMap.IsBound( outFace )) {
      return error(SMESH_Comment("Corresponding inner face not found for face #" )
                   << meshDS->ShapeToIndex( outFace ));
    } else {
      inFace = TopoDS::Face( shape2ShapeMap( outFace ));
    }

    // Find matching nodes of in and out faces
    TNodeNodeMap nodeIn2OutMap;
    if ( ! TAssocTool::FindMatchingNodesOnFaces( inFace, &aMesh, outFace, &aMesh,
                                                 shape2ShapeMap, nodeIn2OutMap ))
      return error(COMPERR_BAD_INPUT_MESH,SMESH_Comment("Mesh on faces #")
                   << meshDS->ShapeToIndex( outFace ) << " and "
                   << meshDS->ShapeToIndex( inFace ) << " seems different" );

    // Create volumes

    SMDS_ElemIteratorPtr faceIt = meshDS->MeshElements( inFace )->GetElements();
    while ( faceIt->more() ) // loop on faces on inFace
    {
      const SMDS_MeshElement* face = faceIt->next();
      if ( !face || face->GetType() != SMDSAbs_Face )
        continue;
      int nbNodes = face->NbNodes();
      if ( face->IsQuadratic() )
        nbNodes /= 2;

      // find node columns for each node
      vector< const TNodeColumn* > columns( nbNodes );
      for ( int i = 0; i < nbNodes; ++i )
      {
        const SMDS_MeshNode* nIn = face->GetNode( i );
        TNode2ColumnMap::iterator n_col = node2columnMap.find( nIn );
        if ( n_col != node2columnMap.end() ) {
          columns[ i ] = & n_col->second;
        }
        else {
          TNodeNodeMap::iterator nInOut = nodeIn2OutMap.find( nIn );
          if ( nInOut == nodeIn2OutMap.end() )
            RETURN_BAD_RESULT("No matching node for "<< nIn->GetID() <<
                              " in face "<< face->GetID());
          columns[ i ] = makeNodeColumn( node2columnMap, nIn, nInOut->second );
        }
      }

      StdMeshers_Prism_3D::AddPrisms( columns, myHelper );
    }
  } // loop on faces of out shell

  return true;
}

bool StdMeshers_Projection_3D::Compute(SMESH_Mesh& aMesh, const TopoDS_Shape& aShape)
{
  if ( !_sourceHypo )
    return false;

  SMESH_Mesh * srcMesh = _sourceHypo->GetSourceMesh();
  SMESH_Mesh * tgtMesh = & aMesh;
  if ( !srcMesh )
    srcMesh = tgtMesh;

  SMESHDS_Mesh * srcMeshDS = srcMesh->GetMeshDS();
  SMESHDS_Mesh * tgtMeshDS = tgtMesh->GetMeshDS();

  // get shell from shape3D
  TopoDS_Shell srcShell, tgtShell;
  TopExp_Explorer exp( _sourceHypo->GetSource3DShape(), TopAbs_SHELL );
  int nbShell;
  for ( nbShell = 0; exp.More(); exp.Next(), ++nbShell )
    srcShell = TopoDS::Shell( exp.Current() );
  if ( nbShell != 1 )
    return error(COMPERR_BAD_SHAPE,
                 SMESH_Comment("Source shape must have 1 shell but not ") << nbShell);

  exp.Init( aShape, TopAbs_SHELL );
  for ( nbShell = 0; exp.More(); exp.Next(), ++nbShell )
    tgtShell = TopoDS::Shell( exp.Current() );
  if ( nbShell != 1 )
    return error(COMPERR_BAD_SHAPE,
                 SMESH_Comment("Target shape must have 1 shell but not ") << nbShell);

  // Check that shapes are blocks
  if ( TAssocTool::Count( tgtShell, TopAbs_FACE , 1 ) != 6 ||
       TAssocTool::Count( tgtShell, TopAbs_EDGE , 1 ) != 12 ||
       TAssocTool::Count( tgtShell, TopAbs_WIRE , 1 ) != 6 )
    return error(COMPERR_BAD_SHAPE, "Target shape is not a block");
  if ( TAssocTool::Count( srcShell, TopAbs_FACE , 1 ) != 6 ||
       TAssocTool::Count( srcShell, TopAbs_EDGE , 1 ) != 12 ||
       TAssocTool::Count( srcShell, TopAbs_WIRE , 1 ) != 6 )
    return error(COMPERR_BAD_SHAPE, "Source shape is not a block");

  // Assure that mesh on a source shape is computed

  SMESH_subMesh* srcSubMesh = srcMesh->GetSubMesh( _sourceHypo->GetSource3DShape() );
  //SMESH_subMesh* tgtSubMesh = tgtMesh->GetSubMesh( aShape );

  if ( tgtMesh == srcMesh && !aShape.IsSame( _sourceHypo->GetSource3DShape() )) {
    if ( !TAssocTool::MakeComputed( srcSubMesh ))
      return error(COMPERR_BAD_INPUT_MESH,"Source mesh not computed");
  }
  else {
    if ( !srcSubMesh->IsMeshComputed() )
      return error(COMPERR_BAD_INPUT_MESH,"Source mesh not computed");
  }

  // Find 2 pairs of corresponding vertices

  TopoDS_Vertex tgtV000, tgtV100, srcV000, srcV100;
  TAssocTool::TShapeShapeMap shape2ShapeMap;

  if ( _sourceHypo->HasVertexAssociation() )
  {
    tgtV000 = _sourceHypo->GetTargetVertex(1);
    tgtV100 = _sourceHypo->GetTargetVertex(2);
    srcV000 = _sourceHypo->GetSourceVertex(1);
    srcV100 = _sourceHypo->GetSourceVertex(2);
  }
  else
  {
    if ( !TAssocTool::FindSubShapeAssociation( tgtShell, tgtMesh, srcShell, srcMesh,
                                               shape2ShapeMap) )
      return error(COMPERR_BAD_SHAPE,"Topology of source and target shapes seems different" );

    exp.Init( tgtShell, TopAbs_EDGE );
    TopExp::Vertices( TopoDS::Edge( exp.Current() ), tgtV000, tgtV100 );

    if ( !shape2ShapeMap.IsBound( tgtV000 ) || !shape2ShapeMap.IsBound( tgtV100 ))
      return error("Association of subshapes failed" );
    srcV000 = TopoDS::Vertex( shape2ShapeMap( tgtV000 ));
    srcV100 = TopoDS::Vertex( shape2ShapeMap( tgtV100 ));
    if ( !TAssocTool::IsSubShape( srcV000, srcShell ) ||
         !TAssocTool::IsSubShape( srcV100, srcShell ))
      return error("Incorrect association of subshapes" );
  }

  // Load 2 SMESH_Block's with src and tgt shells

  SMESH_Block srcBlock, tgtBlock;
  TopTools_IndexedMapOfOrientedShape scrShapes, tgtShapes;
  if ( !tgtBlock.LoadBlockShapes( tgtShell, tgtV000, tgtV100, tgtShapes ))
    return error(COMPERR_BAD_SHAPE, "Can't detect block subshapes. Not a block?");

  if ( !srcBlock.LoadBlockShapes( srcShell, srcV000, srcV100, scrShapes ))
    return error(COMPERR_BAD_SHAPE, "Can't detect block subshapes. Not a block?");

  // Find matching nodes of src and tgt shells

  TNodeNodeMap src2tgtNodeMap;
  for ( int fId = SMESH_Block::ID_FirstF; fId < SMESH_Block::ID_Shell; ++fId )
  {
    // Corresponding subshapes
//.........这里部分代码省略.........

bool StdMeshers_RadialQuadrangle_1D2D::Evaluate(SMESH_Mesh& aMesh,
                                                const TopoDS_Shape& aShape,
                                                MapShapeNbElems& aResMap)
{
  if( aShape.ShapeType() != TopAbs_FACE ) {
    return false;
  }
  SMESH_subMesh * smf = aMesh.GetSubMesh(aShape);
  MapShapeNbElemsItr anIt = aResMap.find(smf);
  if( anIt != aResMap.end() ) {
    return false;
  }

  myLayerPositions.clear();
  gp_Pnt P0(0,0,0);
  gp_Pnt P1(100,0,0);
  computeLayerPositions(P0,P1);

  TopoDS_Edge E1,E2,E3;
  Handle(Geom_Curve) C1,C2,C3;
  double f1,l1,f2,l2,f3,l3;
  int nbe = 0;
  TopExp_Explorer exp;
  for ( exp.Init( aShape, TopAbs_EDGE ); exp.More(); exp.Next() ) {
    nbe++;
    TopoDS_Edge E = TopoDS::Edge( exp.Current() );
    if(nbe==1) {
      E1 = E;
      C1 = BRep_Tool::Curve(E,f1,l1);
    }
    else if(nbe==2) {
      E2 = E;
      C2 = BRep_Tool::Curve(E,f2,l2);
    }
    else if(nbe==3) {
      E3 = E;
      C3 = BRep_Tool::Curve(E,f3,l3);
    }
  }

  TopoDS_Edge CircEdge, LinEdge1, LinEdge2;
  int nb0d=0, nb2d_tria=0, nb2d_quad=0;
  bool isQuadratic = false;
  if(nbe==1) {
    // C1 must be a circle
    Handle(Geom_Circle) aCirc = Handle(Geom_Circle)::DownCast(C1);
    if( !aCirc.IsNull() ) {
      bool ok = _gen->Evaluate( aMesh, CircEdge, aResMap );
      if(ok) {
        SMESH_subMesh * sm = aMesh.GetSubMesh(CircEdge);
        MapShapeNbElemsItr anIt = aResMap.find(sm);
        vector<int> aVec = (*anIt).second;
        isQuadratic = aVec[SMDSEntity_Quad_Edge]>aVec[SMDSEntity_Edge];
        if(isQuadratic) {
          // main nodes
          nb0d = (aVec[SMDSEntity_Node]+1) * myLayerPositions.size();
          // radial medium nodes
          nb0d += (aVec[SMDSEntity_Node]+1) * (myLayerPositions.size()+1);
          // other medium nodes
          nb0d += (aVec[SMDSEntity_Node]+1) * myLayerPositions.size();
        }
        else {
          nb0d = (aVec[SMDSEntity_Node]+1) * myLayerPositions.size();
        }
        nb2d_tria = aVec[SMDSEntity_Node] + 1;
        nb2d_quad = nb0d;
      }
    }
  }
  else if(nbe==2) {
    // one curve must be a half of circle and other curve must be
    // a segment of line
    Handle(Geom_TrimmedCurve) tc = Handle(Geom_TrimmedCurve)::DownCast(C1);
    while( !tc.IsNull() ) {
      C1 = tc->BasisCurve();
      tc = Handle(Geom_TrimmedCurve)::DownCast(C1);
    }
    tc = Handle(Geom_TrimmedCurve)::DownCast(C2);
    while( !tc.IsNull() ) {
      C2 = tc->BasisCurve();
      tc = Handle(Geom_TrimmedCurve)::DownCast(C2);
    }
    Handle(Geom_Circle) aCirc = Handle(Geom_Circle)::DownCast(C1);
    Handle(Geom_Line) aLine = Handle(Geom_Line)::DownCast(C2);
    CircEdge = E1;
    LinEdge1 = E2;
    double fp = f1;
    double lp = l1;
    if( aCirc.IsNull() ) {
      aCirc = Handle(Geom_Circle)::DownCast(C2);
      CircEdge = E2;
      LinEdge1 = E1;
      fp = f2;
      lp = l2;
      aLine = Handle(Geom_Line)::DownCast(C3);
    }
    bool ok = !aCirc.IsNull() && !aLine.IsNull();
    if( fabs(fabs(lp-fp)-M_PI) > Precision::Confusion() ) {
      // not half of circle
      ok = false;
//.........这里部分代码省略.........

bool StdMeshers_RadialQuadrangle_1D2D::Compute(SMESH_Mesh&         aMesh,
                                               const TopoDS_Shape& aShape)
{
  TopExp_Explorer exp;
  SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();

  myHelper = new SMESH_MesherHelper( aMesh );
  myHelper->IsQuadraticSubMesh( aShape );
  // to delete helper at exit from Compute()
  auto_ptr<SMESH_MesherHelper> helperDeleter( myHelper );

  myLayerPositions.clear();

  TopoDS_Edge CircEdge, LinEdge1, LinEdge2;
  int nbe = analyseFace( aShape, CircEdge, LinEdge1, LinEdge2 );
  if( nbe>3 || nbe < 1 || CircEdge.IsNull() )
    return error(COMPERR_BAD_SHAPE);
  
  gp_Pnt P0,P1;
  // points for rotation
  TColgp_SequenceOfPnt Points;
  // angles for rotation
  TColStd_SequenceOfReal Angles;
  // Nodes1 and Nodes2 - nodes along radiuses
  // CNodes - nodes on circle edge
  vector< const SMDS_MeshNode* > Nodes1, Nodes2, CNodes;
  SMDS_MeshNode * NC;
  // parameters edge nodes on face
  TColgp_SequenceOfPnt2d Pnts2d1;
  gp_Pnt2d PC;

  int faceID = meshDS->ShapeToIndex(aShape);
  TopoDS_Face F = TopoDS::Face(aShape);
  Handle(Geom_Surface) S = BRep_Tool::Surface(F);

  if(nbe==1)
  {
    Handle(Geom_Circle) aCirc = Handle(Geom_Circle)::DownCast( getCurve( CircEdge ));

    bool ok = _gen->Compute( aMesh, CircEdge );
    if( !ok ) return false;
    map< double, const SMDS_MeshNode* > theNodes;
    ok = GetSortedNodesOnEdge(aMesh.GetMeshDS(),CircEdge,true,theNodes);
    if( !ok ) return false;

    CNodes.clear();
    map< double, const SMDS_MeshNode* >::iterator itn = theNodes.begin();
    const SMDS_MeshNode* NF = (*itn).second;
    CNodes.push_back( (*itn).second );
    double fang = (*itn).first;
    if ( itn != theNodes.end() ) {
      itn++;
      for(; itn != theNodes.end(); itn++ ) {
        CNodes.push_back( (*itn).second );
        double ang = (*itn).first - fang;
        if( ang>M_PI ) ang = ang - 2*M_PI;
        if( ang<-M_PI ) ang = ang + 2*M_PI;
        Angles.Append( ang ); 
      }
    }
    P1 = gp_Pnt( NF->X(), NF->Y(), NF->Z() );
    P0 = aCirc->Location();

    myLayerPositions.clear();
    computeLayerPositions(P0,P1);

    exp.Init( CircEdge, TopAbs_VERTEX );
    TopoDS_Vertex V1 = TopoDS::Vertex( exp.Current() );
    gp_Pnt2d p2dV = BRep_Tool::Parameters( V1, TopoDS::Face(aShape) );

    NC = meshDS->AddNode(P0.X(), P0.Y(), P0.Z());
    GeomAPI_ProjectPointOnSurf PPS(P0,S);
    double U0,V0;
    PPS.Parameters(1,U0,V0);
    meshDS->SetNodeOnFace(NC, faceID, U0, V0);
    PC = gp_Pnt2d(U0,V0);

    gp_Vec aVec(P0,P1);
    gp_Vec2d aVec2d(PC,p2dV);
    Nodes1.resize( myLayerPositions.size()+1 );
    Nodes2.resize( myLayerPositions.size()+1 );
    int i = 0;
    for(; i<myLayerPositions.size(); i++) {
      gp_Pnt P( P0.X() + aVec.X()*myLayerPositions[i],
                P0.Y() + aVec.Y()*myLayerPositions[i],
                P0.Z() + aVec.Z()*myLayerPositions[i] );
      Points.Append(P);
      SMDS_MeshNode * node = meshDS->AddNode(P.X(), P.Y(), P.Z());
      Nodes1[i] = node;
      Nodes2[i] = node;
      double U = PC.X() + aVec2d.X()*myLayerPositions[i];
      double V = PC.Y() + aVec2d.Y()*myLayerPositions[i];
      meshDS->SetNodeOnFace( node, faceID, U, V );
      Pnts2d1.Append(gp_Pnt2d(U,V));
    }
    Nodes1[Nodes1.size()-1] = NF;
    Nodes2[Nodes1.size()-1] = NF;
  }
  else if(nbe==2 && LinEdge1.Orientation() != TopAbs_INTERNAL )
  {
//.........这里部分代码省略.........

//.........这里部分代码省略.........
        MeshCore::MeshKernel kernel;
        kernel.Adopt(verts, faces, true);

        Mesh::MeshObject* meshdata = new Mesh::MeshObject();
        meshdata->swap(kernel);
        if (createSegm) {
            int index = 0;
            for (auto it : colorMap) {
                Mesh::Segment segm(meshdata, false);
                for (auto jt : it.second) {
                    segm.addIndices(meshSegments[jt]);
                }
                segm.save(true);
                std::stringstream str;
                str << "patch" << index++;
                segm.setName(str.str());
                meshdata->addSegment(segm);
            }
        }
        else {
            for (auto it : meshSegments) {
                meshdata->addSegment(it);
            }
        }
        return meshdata;
    }

#ifndef HAVE_SMESH
    throw Base::Exception("SMESH is not available on this platform");
#else
    std::list<SMESH_Hypothesis*> hypoth;

    SMESH_Gen* meshgen = SMESH_Gen::get();
    SMESH_Mesh* mesh = meshgen->CreateMesh(0, true);
    int hyp=0;

    switch (method) {
#if defined (HAVE_NETGEN)
    case Netgen: {
        NETGENPlugin_Hypothesis_2D* hyp2d = new NETGENPlugin_Hypothesis_2D(hyp++,0,meshgen);

        if (fineness >=0 && fineness < 5) {
            hyp2d->SetFineness(NETGENPlugin_Hypothesis_2D::Fineness(fineness));
        }
        // user defined values
        else {
            if (growthRate > 0)
                hyp2d->SetGrowthRate(growthRate);
            if (nbSegPerEdge > 0)
                hyp2d->SetNbSegPerEdge(nbSegPerEdge);
            if (nbSegPerRadius > 0)
                hyp2d->SetNbSegPerRadius(nbSegPerRadius);
        }

        hyp2d->SetQuadAllowed(allowquad);
        hyp2d->SetOptimize(optimize);
        hyp2d->SetSecondOrder(secondOrder); // apply bisecting to create four triangles out of one
        hypoth.push_back(hyp2d);

        NETGENPlugin_NETGEN_2D* alg2d = new NETGENPlugin_NETGEN_2D(hyp++,0,meshgen);
        hypoth.push_back(alg2d);
    } break;
#endif
#if defined (HAVE_MEFISTO)
    case Mefisto: {
        if (maxLength > 0) {

bool NETGENPlugin_NETGEN_2D_ONLY::Evaluate(SMESH_Mesh& aMesh,
        const TopoDS_Shape& aShape,
        MapShapeNbElems& aResMap)
{
    TopoDS_Face F = TopoDS::Face(aShape);
    if(F.IsNull())
        return false;

    // collect info from edges
    int nb0d = 0, nb1d = 0;
    bool IsQuadratic = false;
    bool IsFirst = true;
    double fullLen = 0.0;
    TopTools_MapOfShape tmpMap;
    for (TopExp_Explorer exp(F, TopAbs_EDGE); exp.More(); exp.Next()) {
        TopoDS_Edge E = TopoDS::Edge(exp.Current());
        if( tmpMap.Contains(E) )
            continue;
        tmpMap.Add(E);
        SMESH_subMesh *aSubMesh = aMesh.GetSubMesh(exp.Current());
        MapShapeNbElemsItr anIt = aResMap.find(aSubMesh);
        if( anIt==aResMap.end() ) {
            SMESH_subMesh *sm = aMesh.GetSubMesh(F);
            SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
            smError.reset( new SMESH_ComputeError(COMPERR_ALGO_FAILED,"Submesh can not be evaluated",this));
            return false;
        }
        std::vector<int> aVec = (*anIt).second;
        nb0d += aVec[SMDSEntity_Node];
        nb1d += Max(aVec[SMDSEntity_Edge],aVec[SMDSEntity_Quad_Edge]);
        double aLen = SMESH_Algo::EdgeLength(E);
        fullLen += aLen;
        if(IsFirst) {
            IsQuadratic = (aVec[SMDSEntity_Quad_Edge] > aVec[SMDSEntity_Edge]);
            IsFirst = false;
        }
    }
    tmpMap.Clear();

    // compute edge length
    double ELen = 0;
    if (_hypLengthFromEdges || (!_hypLengthFromEdges && !_hypMaxElementArea)) {
        if ( nb1d > 0 )
            ELen = fullLen / nb1d;
    }
    if ( _hypMaxElementArea ) {
        double maxArea = _hypMaxElementArea->GetMaxArea();
        ELen = sqrt(2. * maxArea/sqrt(3.0));
    }
    GProp_GProps G;
    BRepGProp::SurfaceProperties(F,G);
    double anArea = G.Mass();

    const int hugeNb = numeric_limits<int>::max()/10;
    if ( anArea / hugeNb > ELen*ELen )
    {
        SMESH_subMesh *sm = aMesh.GetSubMesh(F);
        SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
        smError.reset( new SMESH_ComputeError(COMPERR_ALGO_FAILED,"Submesh can not be evaluated.\nToo small element length",this));
        return false;
    }
    int nbFaces = (int) ( anArea / ( ELen*ELen*sqrt(3.) / 4 ) );
    int nbNodes = (int) ( ( nbFaces*3 - (nb1d-1)*2 ) / 6 + 1 );
    std::vector<int> aVec(SMDSEntity_Last);
    for(int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aVec[i]=0;
    if( IsQuadratic ) {
        aVec[SMDSEntity_Node] = nbNodes;
        aVec[SMDSEntity_Quad_Triangle] = nbFaces;
    }
    else {
        aVec[SMDSEntity_Node] = nbNodes;
        aVec[SMDSEntity_Triangle] = nbFaces;
    }
    SMESH_subMesh *sm = aMesh.GetSubMesh(F);
    aResMap.insert(std::make_pair(sm,aVec));

    return true;
}

static bool checkConformIgnoredAlgos(SMESH_Mesh&               aMesh,
                                     SMESH_subMesh*            aSubMesh,
                                     const SMESH_Algo*         aGlobIgnoAlgo,
                                     const SMESH_Algo*         aLocIgnoAlgo,
                                     bool &                    checkConform,
                                     set<SMESH_subMesh*>&      aCheckedMap,
                                     list< SMESH_Gen::TAlgoStateError > & theErrors)
{
  ASSERT( aSubMesh );
  if ( aSubMesh->GetSubShape().ShapeType() == TopAbs_VERTEX)
    return true;


  bool ret = true;

  const list<const SMESHDS_Hypothesis*>& listHyp =
    aMesh.GetMeshDS()->GetHypothesis( aSubMesh->GetSubShape() );
  list<const SMESHDS_Hypothesis*>::const_iterator it=listHyp.begin();
  for ( ; it != listHyp.end(); it++)
  {
    const SMESHDS_Hypothesis * aHyp = *it;
    if (aHyp->GetType() == SMESHDS_Hypothesis::PARAM_ALGO)
      continue;

    const SMESH_Algo* algo = dynamic_cast<const SMESH_Algo*> (aHyp);
    ASSERT ( algo );

    if ( aLocIgnoAlgo ) // algo is hidden by a local algo of upper dim
    {
      theErrors.push_back( SMESH_Gen::TAlgoStateError() );
      theErrors.back().Set( SMESH_Hypothesis::HYP_HIDDEN_ALGO, algo, false );
      INFOS( "Local <" << algo->GetName() << "> is hidden by local <"
            << aLocIgnoAlgo->GetName() << ">");
    }
    else
    {
      bool       isGlobal = (aMesh.IsMainShape( aSubMesh->GetSubShape() ));
      int             dim = algo->GetDim();
      int aMaxGlobIgnoDim = ( aGlobIgnoAlgo ? aGlobIgnoAlgo->GetDim() : -1 );
      bool    isNeededDim = ( aGlobIgnoAlgo ? aGlobIgnoAlgo->NeedLowerHyps( dim ) : false );

      if (( dim < aMaxGlobIgnoDim && !isNeededDim ) &&
          ( isGlobal || !aGlobIgnoAlgo->SupportSubmeshes() ))
      {
        // algo is hidden by a global algo
        theErrors.push_back( SMESH_Gen::TAlgoStateError() );
        theErrors.back().Set( SMESH_Hypothesis::HYP_HIDDEN_ALGO, algo, true );
        INFOS( ( isGlobal ? "Global" : "Local" )
              << " <" <<