//=======================================================================
//function : CorrectWires
//purpose  :
//=======================================================================
  Standard_Boolean GEOMAlgo_Tools::CorrectWires(const TopoDS_Shape& aShape)
{
  Standard_Boolean bRet;
  TopoDS_Iterator aItF;
  TopExp_Explorer aExp;
  TopTools_MapOfShape aMF;
  GeomAdaptor_Surface aGAS;
  GeomAbs_SurfaceType aTS;
  TopLoc_Location aLoc;
  //
  bRet=Standard_False;
  //
  aExp.Init(aShape, TopAbs_FACE);
  for (; aExp.More(); aExp.Next()) {
    const TopoDS_Face& aF=*((TopoDS_Face*)&aExp.Current());
    if (aMF.Add(aF)) {
      const Handle(Geom_Surface)& aS=BRep_Tool::Surface(aF, aLoc);
      aGAS.Load(aS);
      aTS=aGAS.GetType();
      if (aTS==GeomAbs_Cylinder || aTS==GeomAbs_Plane) {
	aItF.Initialize(aF);
	for (; aItF.More(); aItF.Next()) {
	  const TopoDS_Wire& aW=*((TopoDS_Wire*)&aItF.Value());
	  if (CorrectWire(aW, aF)) {
	   bRet=Standard_True;
	  }
	}
      }
    }
  }
  return bRet;
}

//=======================================================================
//function : FillContainers
//purpose  :
//=======================================================================
void GEOMAlgo_Gluer2::FillContainers(const TopAbs_ShapeEnum aType)
{
  Standard_Boolean bHasImage, bToReverse;
  Standard_Integer i, aNbW;
  TopoDS_Shape aWnew, aEnew;
  TopoDS_Iterator aItS;
  BRep_Builder aBB;
  TopTools_IndexedMapOfShape aMW;
  TopTools_MapOfShape aMFence;
  //
  myErrorStatus=0;
  myWarningStatus=0;
  //
  TopExp::MapShapes(myArgument, aType, aMW);
  //
  aNbW=aMW.Extent();
  for (i=1; i<=aNbW; ++i) {
    const TopoDS_Shape& aW=aMW(i);
    //
    if (!aMFence.Add(aW)) {
      continue;
    }
    //
    bHasImage=HasImage(aW);
    if (!bHasImage) {
      continue;
    }
    //
    GEOMAlgo_Tools3D::MakeContainer(aType, aWnew);
    aWnew.Orientation(aW.Orientation());
    //
    aItS.Initialize(aW);
    for (; aItS.More(); aItS.Next()) {
      const TopoDS_Shape& aE=aItS.Value();
      if (myOrigins.IsBound(aE)) {
        aEnew=myOrigins.Find(aE);
        //
        bToReverse=GEOMAlgo_Tools3D::IsSplitToReverse(aEnew, aE, myContext);
        if (bToReverse) {
          aEnew.Reverse();
        }
        //
        aBB.Add(aWnew, aEnew);
      }
      else {
        aBB.Add(aWnew, aE);
      }
    }
    //
    //myImages / myOrigins
    TopTools_ListOfShape aLSD;
    //
    aLSD.Append(aW);
    myImages.Bind(aWnew, aLSD);
    myOrigins.Bind(aW, aWnew);
    //
  }//for (i=1; i<=aNbE; ++i) {
}

//=======================================================================
//function : FillBRepShapes
//purpose  :
//=======================================================================
void GEOMAlgo_Gluer2::FillBRepShapes(const TopAbs_ShapeEnum theType)
{
  Standard_Boolean bHasImage, bIsToWork;
  Standard_Integer i, aNbE;
  TopoDS_Iterator aItS;
  TopoDS_Shape aEnew;
  TopTools_IndexedMapOfShape aME;
  TopTools_MapOfShape aMFence;
  TopTools_ListIteratorOfListOfShape aItLS;
  //
  myErrorStatus=0;
  myWarningStatus=0;
  //
  TopExp::MapShapes(myArgument, theType, aME);
  //
  aNbE=aME.Extent();
  for (i=1; i<=aNbE; ++i) {
    const TopoDS_Shape& aE=aME(i);
    //
    if (!aMFence.Add(aE)) {
      continue;
    }
    //
    bIsToWork=myOriginsToWork.IsBound(aE);
    bHasImage=HasImage(aE);
    if (!bHasImage && !bIsToWork) {
      continue;
    }
    //
    MakeBRepShapes(aE, aEnew);
    //
    //myImages / myOrigins
    if (bIsToWork) {
      const TopoDS_Shape& aSkey=myOriginsToWork.Find(aE);
      const TopTools_ListOfShape& aLSD=myImagesToWork.Find(aSkey);
      //
      myImages.Bind(aEnew, aLSD);
      //
      aItLS.Initialize(aLSD);
      for (; aItLS.More(); aItLS.Next()) {
        const TopoDS_Shape& aEx=aItLS.Value();
        myOrigins.Bind(aEx, aEnew);
        //
        aMFence.Add(aEx);
      }
    }
    else {
      TopTools_ListOfShape aLSD;
      //
      aLSD.Append(aE);
      myImages.Bind(aEnew, aLSD);
      myOrigins.Bind(aE, aEnew);
    }
  }//for (i=1; i<=aNbF; ++i) {
}

int SMESH_MesherHelper::NbAncestors(const TopoDS_Shape& shape,
                                    const SMESH_Mesh&   mesh,
                                    TopAbs_ShapeEnum    ancestorType/*=TopAbs_SHAPE*/)
{
  TopTools_MapOfShape ancestors;
  TopTools_ListIteratorOfListOfShape ansIt( mesh.GetAncestors(shape) );
  for ( ; ansIt.More(); ansIt.Next() ) {
    if ( ancestorType == TopAbs_SHAPE || ansIt.Value().ShapeType() == ancestorType )
      ancestors.Add( ansIt.Value() );
  }
  return ancestors.Extent();
}

//=======================================================================
//function :  SupressFaces
//purpose  :
//=======================================================================
void SuppressFacesRec (const TopTools_SequenceOfShape& theShapesFaces,
                       const TopoDS_Shape&             theOriginalShape,
                       TopoDS_Shape&                   theOutShape)
{
  if ((theOriginalShape.ShapeType() != TopAbs_COMPOUND &&
       theOriginalShape.ShapeType() != TopAbs_COMPSOLID))
  {
    ShHealOper_RemoveFace aHealer (theOriginalShape);
    Standard_Boolean aResult = aHealer.Perform(theShapesFaces);

    if (aResult)
      theOutShape = aHealer.GetResultShape();
    else
      raiseNotDoneExeption(aHealer.GetErrorStatus());
  }
  else
  {
    BRep_Builder BB;
    TopoDS_Compound CC;
    BB.MakeCompound(CC);

    TopTools_MapOfShape mapShape;
    TopoDS_Iterator It (theOriginalShape, Standard_True, Standard_True);

    for (; It.More(); It.Next()) {
      TopoDS_Shape aShape_i = It.Value();
      if (mapShape.Add(aShape_i)) {
        // check, if current shape contains at least one of faces to be removed
        bool isFound = false;
        TopTools_IndexedMapOfShape aShapes_i;
        TopExp::MapShapes(aShape_i, aShapes_i);
        for (int i = 1; i <= theShapesFaces.Length() && !isFound; i++) {
          const TopoDS_Shape& aFace_i = theShapesFaces.Value(i);
          if (aShapes_i.Contains(aFace_i)) isFound = true;
        }
        if (isFound) {
          TopoDS_Shape anOutSh_i;
          SuppressFacesRec(theShapesFaces, aShape_i, anOutSh_i);
          if ( !anOutSh_i.IsNull() )
            BB.Add(CC, anOutSh_i);
        }
        else {
          // nothing to do
          BB.Add(CC, aShape_i);
        }
      }
    }
    theOutShape = CC;
  }
}

//=======================================================================
//function : MakeInternalWires
//purpose  : 
//=======================================================================
void MakeInternalWires(const TopTools_MapOfShape& theME,
		       TopTools_ListOfShape& theWires)
{
  TopTools_MapIteratorOfMapOfShape aItM;
  TopTools_MapOfShape aAddedMap;
  TopTools_ListIteratorOfListOfShape aItE;
  TopTools_IndexedDataMapOfShapeListOfShape aMVE;
  BRep_Builder aBB;
  //
  aItM.Initialize(theME);
  for (; aItM.More(); aItM.Next()) {
    const TopoDS_Shape& aE=aItM.Key();
    TopExp::MapShapesAndAncestors(aE, TopAbs_VERTEX, TopAbs_EDGE, aMVE);
  }
  //
  aItM.Initialize(theME);
  for (; aItM.More(); aItM.Next()) {
    TopoDS_Shape aEE=aItM.Key();
    if (!aAddedMap.Add(aEE)) {
      continue;
    }
    //
    // make a new shell
    TopoDS_Wire aW;
    aBB.MakeWire(aW);    
    aEE.Orientation(TopAbs_INTERNAL);
    aBB.Add(aW, aEE);
    //
    TopoDS_Iterator aItAdded (aW);
    for (; aItAdded.More(); aItAdded.Next()) {
      const TopoDS_Shape& aE =aItAdded.Value();
      //
      TopExp_Explorer aExp(aE, TopAbs_VERTEX);
      for (; aExp.More(); aExp.Next()) {
        const TopoDS_Shape& aV =aExp.Current();
	const TopTools_ListOfShape& aLE=aMVE.FindFromKey(aV);
	aItE.Initialize(aLE);
	for (; aItE.More(); aItE.Next()) { 
	  TopoDS_Shape aEL=aItE.Value();
	  if (aAddedMap.Add(aEL)){
	    aEL.Orientation(TopAbs_INTERNAL);
	    aBB.Add(aW, aEL);
	  }
	}
      }
    }
    theWires.Append(aW);
  }
}

void BRepOffsetAPI_MakeOffsetFix::MakeWire(TopoDS_Shape& wire)
{
    // get the edges of the wire and check which of the stored edges
    // serve as input of the wire
    TopTools_MapOfShape aMap;
    TopExp_Explorer xp(wire, TopAbs_EDGE);
    while (xp.More()) {
        aMap.Add(xp.Current());
        xp.Next();
    }

    std::list<TopoDS_Edge> edgeList;
    for (auto itLoc : myLocations) {
        const TopTools_ListOfShape& newShapes = mkOffset.Generated(itLoc.first);
        for (TopTools_ListIteratorOfListOfShape it(newShapes); it.More(); it.Next()) {
            TopoDS_Shape newShape = it.Value();

            if (aMap.Contains(newShape)) {
                newShape.Move(itLoc.second);
                edgeList.push_back(TopoDS::Edge(newShape));
            }
        }
    }

    if (!edgeList.empty()) {
        BRepBuilderAPI_MakeWire mkWire;
        mkWire.Add(edgeList.front());
        edgeList.pop_front();
        wire = mkWire.Wire();

        bool found = false;
        do {
            found = false;
            for (std::list<TopoDS_Edge>::iterator pE = edgeList.begin(); pE != edgeList.end(); ++pE) {
                mkWire.Add(*pE);
                if (mkWire.Error() != BRepBuilderAPI_DisconnectedWire) {
                    // edge added ==> remove it from list
                    found = true;
                    edgeList.erase(pE);
                    wire = mkWire.Wire();
                    break;
                }
            }
        }
        while (found);
    }
}

Standard_Boolean ShHealOper_Sewing::getWires(const TopoDS_Shape& theSewShape) const
{
  if(theSewShape.ShapeType() != TopAbs_COMPOUND)
    return Standard_False;
  
  Handle(TopTools_HSequenceOfShape) aSeqEdges = new TopTools_HSequenceOfShape;
  TopExp_Explorer aexpEdges(theSewShape,TopAbs_EDGE,TopAbs_WIRE);
  for ( ; aexpEdges.More(); aexpEdges.Next()) {
    aSeqEdges->Append(aexpEdges.Current());
  }
  if(aSeqEdges->Length() <2)
    return Standard_False;
  //get manifold wires from sewed edges.
  Standard_Real aTol = 0.;
  Standard_Boolean aShared = Standard_True;
  Handle(TopTools_HSequenceOfShape) aTmpWires = new TopTools_HSequenceOfShape;
  ShapeAnalysis_FreeBounds::ConnectEdgesToWires(aSeqEdges, aTol, aShared, aTmpWires);

  TopTools_MapOfShape aMapEdges;
  Standard_Integer i =1;
  for( ; i <= aSeqEdges->Length(); i++)
    aMapEdges.Add(aSeqEdges->Value(i));
  
  //remove free edges from result shape.
  TopoDS_Compound aNewComp;
  deleteFreeEdges(theSewShape,aMapEdges,aNewComp);
  
  //add new wires in the result shape.
  BRep_Builder aB;
  for( i =1; i <= aTmpWires->Length(); i++) {
    TopoDS_Iterator aite(aTmpWires->Value(i));
    Standard_Integer nbe =0;
    TopoDS_Shape aE;
    for( ; aite.More() && nbe < 3; aite.Next(),nbe++)
      aE = aite.Value();
    if(!nbe)
      continue;
    else if(nbe ==1)
      aB.Add(aNewComp,aE);
    else
      aB.Add(aNewComp,aTmpWires->Value(i));
  }
  
  myContext->Replace(theSewShape,aNewComp);
  return Standard_True;
}

//=======================================================================
//function : MapShapes
//purpose  :
//=======================================================================
void MapShapes(const TopoDS_Shape& theS,
               TopTools_MapOfShape& theM)
{
    theM.Add(theS);
    TopoDS_Iterator anIt;
    anIt.Initialize(theS);
    for (; anIt.More(); anIt.Next()) {
        const TopoDS_Shape& aSx=anIt.Value();
        MapShapes(aSx, theM);
    }
}

//=======================================================================
//function : FillImagesCompound
//purpose  :
//=======================================================================
void FillImagesCompound(const TopoDS_Shape& theS,
                        BRepAlgo_Image& theImages,
                        TopTools_MapOfShape& theMFP)
{
    Standard_Boolean bInterferred;
    TopAbs_ShapeEnum aTypeX;
    TopAbs_Orientation aOrX;
    TopoDS_Iterator aIt;
    BRep_Builder aBB;
    TopTools_ListIteratorOfListOfShape aItIm;
    //
    if (!theMFP.Add(theS)) {
        return;
    }
    //
    bInterferred=Standard_False;
    aIt.Initialize(theS);
    for (; aIt.More(); aIt.Next()) {
        const TopoDS_Shape& aSX=aIt.Value();
        aTypeX=aSX.ShapeType();
        if (aTypeX==TopAbs_COMPOUND) {
            FillImagesCompound(aSX, theImages, theMFP);
        }
        if (theImages.HasImage(aSX)) {
            bInterferred=Standard_True;
        }
    }
    if (!bInterferred) {
        return;
    }
    //
    TopoDS_Shape aCIm;
    GEOMAlgo_Tools3D::MakeContainer(TopAbs_COMPOUND, aCIm);
    //
    aIt.Initialize(theS);
    for (; aIt.More(); aIt.Next()) {
        const TopoDS_Shape& aSX=aIt.Value();
        aOrX=aSX.Orientation();
        if (theImages.HasImage(aSX)) {
            const TopTools_ListOfShape& aLFIm=theImages.Image(aSX);
            aItIm.Initialize(aLFIm);
            for (; aItIm.More(); aItIm.Next()) {
                TopoDS_Shape aSXIm=aItIm.Value();
                aSXIm.Orientation(aOrX);
                aBB.Add(aCIm, aSXIm);
            }
        }
        else {
            aBB.Add(aCIm, aSX);
        }
    }
    theImages.Bind(theS, aCIm);
}

//=======================================================================
//function : CompsolidToCompound
//purpose  :
//=======================================================================
TopoDS_Shape GEOMUtils::CompsolidToCompound (const TopoDS_Shape& theCompsolid)
{
  if (theCompsolid.ShapeType() != TopAbs_COMPSOLID) {
    return theCompsolid;
  }

  TopoDS_Compound aCompound;
  BRep_Builder B;
  B.MakeCompound(aCompound);

  TopTools_MapOfShape mapShape;
  TopoDS_Iterator It (theCompsolid, Standard_True, Standard_True);

  for (; It.More(); It.Next()) {
    TopoDS_Shape aShape_i = It.Value();
    if (mapShape.Add(aShape_i)) {
      B.Add(aCompound, aShape_i);
    }
  }

  return aCompound;
}

//modified by NIZNHY-PKV Thu Feb 16 12:25:16 2012f
//=======================================================================
//function : IsClosed
//purpose  :
//=======================================================================
Standard_Boolean IsClosed(const TopoDS_Edge& aE,
			  const TopoDS_Face& aF)
{
  Standard_Boolean bRet;
  //
  bRet=BRep_Tool::IsClosed(aE, aF);
  if (bRet) {
    TopTools_MapOfShape aM;
    TopExp_Explorer aExp(aF, TopAbs_EDGE);
    for (; aExp.More(); aExp.Next()) {
      const TopoDS_Shape& aEx=aExp.Current();
      //
      if (aM.Add(aEx)) {
	bRet=aEx.IsSame(aE);
	if (bRet) {
	  break;
	}
      }
    }
  }
  return bRet;
}

//=======================================================================
//function : AddSimpleShapes
//purpose  :
//=======================================================================
void GEOMUtils::AddSimpleShapes (const TopoDS_Shape& theShape, TopTools_ListOfShape& theList)
{
  if (theShape.ShapeType() != TopAbs_COMPOUND &&
      theShape.ShapeType() != TopAbs_COMPSOLID) {
    theList.Append(theShape);
    return;
  }

  TopTools_MapOfShape mapShape;
  TopoDS_Iterator It (theShape, Standard_True, Standard_True);

  for (; It.More(); It.Next()) {
    TopoDS_Shape aShape_i = It.Value();
    if (mapShape.Add(aShape_i)) {
      if (aShape_i.ShapeType() == TopAbs_COMPOUND ||
          aShape_i.ShapeType() == TopAbs_COMPSOLID) {
        AddSimpleShapes(aShape_i, theList);
      } else {
        theList.Append(aShape_i);
      }
    }
  }
}

//=======================================================================
//function : isGoodForChamfer
//purpose  :
//=======================================================================
static Standard_Boolean isGoodForChamfer (const TopoDS_Shape& theShape)
{
  if (theShape.ShapeType() == TopAbs_SHELL ||
      theShape.ShapeType() == TopAbs_SOLID ||
      theShape.ShapeType() == TopAbs_COMPSOLID) {
    return Standard_True;
  }

  if (theShape.ShapeType() == TopAbs_COMPOUND) {
    TopTools_MapOfShape mapShape;
    TopoDS_Iterator It (theShape, Standard_False, Standard_False);
    for (; It.More(); It.Next()) {
      if (mapShape.Add(It.Value())) {
        if (!isGoodForChamfer(It.Value())) {
          return Standard_False;
        }
      }
    }
    return Standard_True;
  }

  return Standard_False;
}

//=======================================================================
//function : Execute
//purpose  :
//=======================================================================
Standard_Integer GEOMImpl_Fillet1dDriver::Execute(TFunction_Logbook& log) const
{
  if (Label().IsNull()) return 0;
  Handle(GEOM_Function) aFunction = GEOM_Function::GetFunction(Label());

  GEOMImpl_IFillet1d aCI (aFunction);

  Handle(GEOM_Function) aRefShape = aCI.GetShape();
  TopoDS_Shape aShape = aRefShape->GetValue();
  if (aShape.IsNull())
    return 0;
  if (aShape.ShapeType() != TopAbs_WIRE)
    Standard_ConstructionError::Raise("Wrong arguments: polyline as wire must be given");

  TopoDS_Wire aWire = TopoDS::Wire(aShape);

  double rad = aCI.GetR();

  if ( rad < Precision::Confusion())
    return 0;

  // collect vertices for make fillet
  TopTools_ListOfShape aVertexList;
  TopTools_MapOfShape mapShape;
  int aLen = aCI.GetLength();
  if ( aLen > 0 ) {
    for (int ind = 1; ind <= aLen; ind++) {
      TopoDS_Shape aShapeVertex;
      if (GEOMImpl_ILocalOperations::GetSubShape
          (aWire, aCI.GetVertex(ind), aShapeVertex))
        if (mapShape.Add(aShapeVertex))
          aVertexList.Append( aShapeVertex );
    }
  } else { // get all vertices from wire
    TopExp_Explorer anExp( aWire, TopAbs_VERTEX );
    for ( ; anExp.More(); anExp.Next() ) {
      if (mapShape.Add(anExp.Current()))
        aVertexList.Append( anExp.Current() );
    }
  }
  if (aVertexList.IsEmpty())
    Standard_ConstructionError::Raise("Invalid input no vertices to make fillet");

  //INFO: this algorithm implemented in assumption that user can select both
  //  vertices of some edges to make fillet. In this case we should remember
  //  already modified initial edges to take care in next fillet step
  TopTools_DataMapOfShapeShape anEdgeToEdgeMap;

  //iterates on vertices, and make fillet on each couple of edges
  //collect result fillet edges in list
  TopTools_ListOfShape aListOfNewEdge;
  // remember relation between initial and modified map
  TopTools_IndexedDataMapOfShapeListOfShape aMapVToEdges;
  TopExp::MapShapesAndAncestors( aWire, TopAbs_VERTEX, TopAbs_EDGE, aMapVToEdges );
  TopTools_ListIteratorOfListOfShape anIt( aVertexList );
  for ( ; anIt.More(); anIt.Next() ) {
    TopoDS_Vertex aV = TopoDS::Vertex( anIt.Value() );
    if ( aV.IsNull() || !aMapVToEdges.Contains( aV ) )
      continue;
    const TopTools_ListOfShape& aVertexEdges = aMapVToEdges.FindFromKey( aV );
    if ( aVertexEdges.Extent() != 2 )
      continue; // no input data to make fillet
    TopoDS_Edge anEdge1 = TopoDS::Edge( aVertexEdges.First() );
    TopoDS_Edge anEdge2 = TopoDS::Edge( aVertexEdges.Last() );
    // check if initial edges already modified in previous fillet operation
    if ( anEdgeToEdgeMap.IsBound( anEdge1 ) ) anEdge1 = TopoDS::Edge(anEdgeToEdgeMap.Find( anEdge1 ));
    if ( anEdgeToEdgeMap.IsBound( anEdge2 ) ) anEdge2 = TopoDS::Edge(anEdgeToEdgeMap.Find( anEdge2 ));
    if ( anEdge1.IsNull() || anEdge2.IsNull() || anEdge1.IsSame( anEdge2 ) )
      continue; //no input data to make fillet

    // create plane on 2 edges
    gp_Pln aPlane;
    if ( !takePlane(anEdge1, anEdge2, aV, aPlane) )
      continue; // seems edges does not belong to same plane or parallel (fillet can not be build)

    GEOMImpl_Fillet1d aFilletAlgo(anEdge1, anEdge2, aPlane);
    if ( !aFilletAlgo.Perform(rad) )
      continue; // can not create fillet with given radius

    // take fillet result in given vertex
    TopoDS_Edge aModifE1, aModifE2;
    TopoDS_Edge aNewE = aFilletAlgo.Result(BRep_Tool::Pnt(aV), aModifE1, aModifE2);
    if (aNewE.IsNull())
      continue; // no result found

    // add  new created edges and take modified edges
    aListOfNewEdge.Append( aNewE );

    // check if face edges modified,
    // if yes, than map to original edges (from vertex-edges list), because edges can be modified before
    if (aModifE1.IsNull() || !anEdge1.IsSame( aModifE1 ))
      addEdgeRelation( anEdgeToEdgeMap, TopoDS::Edge(aVertexEdges.First()), aModifE1 );
    if (aModifE2.IsNull() || !anEdge2.IsSame( aModifE2 ))
      addEdgeRelation( anEdgeToEdgeMap, TopoDS::Edge(aVertexEdges.Last()), aModifE2 );
  }

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

bool NETGENPlugin_Mesher::fillNgMesh(netgen::OCCGeometry&           occgeom,
                                     netgen::Mesh&                  ngMesh,
                                     vector<SMDS_MeshNode*>&        nodeVec,
                                     const list< SMESH_subMesh* > & meshedSM)
{
  TNode2IdMap nodeNgIdMap;

  TopTools_MapOfShape visitedShapes;

  SMESH_MesherHelper helper (*_mesh);

  int faceID = occgeom.fmap.Extent();
  _faceDescriptors.clear();

  list< SMESH_subMesh* >::const_iterator smIt, smEnd = meshedSM.end();
  for ( smIt = meshedSM.begin(); smIt != smEnd; ++smIt )
  {
    SMESH_subMesh* sm = *smIt;
    if ( !visitedShapes.Add( sm->GetSubShape() ))
      continue;

    SMESHDS_SubMesh * smDS = sm->GetSubMeshDS();

    switch ( sm->GetSubShape().ShapeType() )
    {
    case TopAbs_EDGE: { // EDGE
      // ----------------------
      const TopoDS_Edge& geomEdge  = TopoDS::Edge( sm->GetSubShape() );

      // Add ng segments for each not meshed face the edge bounds
      TopTools_MapOfShape visitedAncestors;
      const TopTools_ListOfShape& ancestors = _mesh->GetAncestors( geomEdge );
      TopTools_ListIteratorOfListOfShape ancestorIt ( ancestors );
      for ( ; ancestorIt.More(); ancestorIt.Next() )
      {
        const TopoDS_Shape & ans = ancestorIt.Value();
        if ( ans.ShapeType() != TopAbs_FACE || !visitedAncestors.Add( ans ))
          continue;
        const TopoDS_Face& face = TopoDS::Face( ans );

        int faceID = occgeom.fmap.FindIndex( face );
        if ( faceID < 1 )
          continue; // meshed face

        // find out orientation of geomEdge within face
        bool isForwad = false;
        for ( TopExp_Explorer exp( face, TopAbs_EDGE ); exp.More(); exp.Next() ) {
          if ( geomEdge.IsSame( exp.Current() )) {
            isForwad = ( exp.Current().Orientation() == geomEdge.Orientation() );
            break;
          }
        }
        bool isQuad = smDS->GetElements()->next()->IsQuadratic();

        // get all nodes from geomEdge
        StdMeshers_FaceSide fSide( face, geomEdge, _mesh, isForwad, isQuad );
        const vector<UVPtStruct>& points = fSide.GetUVPtStruct();
        int i, nbSeg = fSide.NbSegments();

        double otherSeamParam = 0;
        helper.SetSubShape( face );
        bool isSeam = helper.IsRealSeam( geomEdge );
        if ( isSeam )
          otherSeamParam =
            helper.GetOtherParam( helper.GetPeriodicIndex() == 1 ? points[0].u : points[0].v );

        // add segments

        int prevNgId = ngNodeId( points[0].node, ngMesh, nodeNgIdMap );

        for ( i = 0; i < nbSeg; ++i )
        {
          const UVPtStruct& p1 = points[ i ];
          const UVPtStruct& p2 = points[ i+1 ];

          netgen::Segment seg;
          // ng node ids
          seg.p1 = prevNgId;
          seg.p2 = prevNgId = ngNodeId( p2.node, ngMesh, nodeNgIdMap );
          // node param on curve
          seg.epgeominfo[ 0 ].dist = p1.param;
          seg.epgeominfo[ 1 ].dist = p2.param;
          // uv on face
          seg.epgeominfo[ 0 ].u = p1.u;
          seg.epgeominfo[ 0 ].v = p1.v;
          seg.epgeominfo[ 1 ].u = p2.u;
          seg.epgeominfo[ 1 ].v = p2.v;

          //seg.epgeominfo[ iEnd ].edgenr = edgeID; //  = geom.emap.FindIndex(edge);
          seg.si = faceID;                   // = geom.fmap.FindIndex (face);
          seg.edgenr = ngMesh.GetNSeg() + 1; // segment id
          ngMesh.AddSegment (seg);

          if ( isSeam )
          {
            if ( helper.GetPeriodicIndex() == 1 ) {
              seg.epgeominfo[ 0 ].u = otherSeamParam;
              seg.epgeominfo[ 1 ].u = otherSeamParam;
              swap (seg.epgeominfo[0].v, seg.epgeominfo[1].v);
            } else {
//.........这里部分代码省略.........

//=============================================================================
bool NETGENPlugin_Mesher::Compute()
{
#ifdef WNT
  netgen::MeshingParameters& mparams = netgen::GlobalMeshingParameters();
#else
  netgen::MeshingParameters& mparams = netgen::mparam;
#endif  
  MESSAGE("Compute with:\n"
          " max size = " << mparams.maxh << "\n"
          " segments per edge = " << mparams.segmentsperedge);
  MESSAGE("\n"
          " growth rate = " << mparams.grading << "\n"
          " elements per radius = " << mparams.curvaturesafety << "\n"
          " second order = " << mparams.secondorder << "\n"
          " quad allowed = " << mparams.quad);

  SMESH_ComputeErrorPtr error = SMESH_ComputeError::New();
  nglib::Ng_Init();

  // -------------------------
  // Prepare OCC geometry
  // -------------------------

  netgen::OCCGeometry occgeo;
  list< SMESH_subMesh* > meshedSM;
  PrepareOCCgeometry( occgeo, _shape, *_mesh, &meshedSM );

  // -------------------------
  // Generate the mesh
  // -------------------------

  netgen::Mesh *ngMesh = NULL;

  SMESH_Comment comment;
  int err = 0;
  int nbInitNod = 0;
  int nbInitSeg = 0;
  int nbInitFac = 0;
  // vector of nodes in which node index == netgen ID
  vector< SMDS_MeshNode* > nodeVec;
  try
  {
    // ----------------
    // compute 1D mesh
    // ----------------
    // pass 1D simple parameters to NETGEN
    if ( _simpleHyp ) {
      if ( int nbSeg = _simpleHyp->GetNumberOfSegments() ) {
        // nb of segments
        mparams.segmentsperedge = nbSeg + 0.1;
        mparams.maxh = occgeo.boundingbox.Diam();
        mparams.grading = 0.01;
      }
      else {
        // segment length
        mparams.segmentsperedge = 1;
        mparams.maxh = _simpleHyp->GetLocalLength();
      }
    }
    // let netgen create ngMesh and calculate element size on not meshed shapes
    char *optstr = 0;
    int startWith = netgen::MESHCONST_ANALYSE;
    int endWith   = netgen::MESHCONST_ANALYSE;
    err = netgen::OCCGenerateMesh(occgeo, ngMesh, startWith, endWith, optstr);
    if (err) comment << "Error in netgen::OCCGenerateMesh() at MESHCONST_ANALYSE step";

    // fill ngMesh with nodes and elements of computed submeshes
    err = ! fillNgMesh(occgeo, *ngMesh, nodeVec, meshedSM);
    nbInitNod = ngMesh->GetNP();
    nbInitSeg = ngMesh->GetNSeg();
    nbInitFac = ngMesh->GetNSE();

    // compute mesh
    if (!err)
    {
      startWith = endWith = netgen::MESHCONST_MESHEDGES;
      err = netgen::OCCGenerateMesh(occgeo, ngMesh, startWith, endWith, optstr);
      if (err) comment << "Error in netgen::OCCGenerateMesh() at 1D mesh generation";
    }
    // ---------------------
    // compute surface mesh
    // ---------------------
    if (!err)
    {
      // pass 2D simple parameters to NETGEN
      if ( _simpleHyp ) {
        if ( double area = _simpleHyp->GetMaxElementArea() ) {
          // face area
          mparams.maxh = sqrt(2. * area/sqrt(3.0));
          mparams.grading = 0.4; // moderate size growth
        }
        else {
          // length from edges
          double length = 0;
          TopTools_MapOfShape tmpMap;
          for ( TopExp_Explorer exp( _shape, TopAbs_EDGE ); exp.More(); exp.Next() )
            if( tmpMap.Add(exp.Current()) )
              length += SMESH_Algo::EdgeLength( TopoDS::Edge( exp.Current() ));

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

//=======================================================================
// function: FillInternalVertices
// purpose:
//=======================================================================
void GEOMAlgo_Builder::FillInternalVertices()
{
  const NMTDS_ShapesDataStructure& aDS=*myPaveFiller->DS();
  NMTTools_PaveFiller* pPF=myPaveFiller;
  NMTDS_InterfPool* pIP=pPF->IP();
  const Handle(IntTools_Context)& aCtx= pPF->Context();
  //
  BOPTools_CArray1OfVSInterference& aVFs=pIP->VSInterferences();
  BOPTools_CArray1OfESInterference& aEFs=pIP->ESInterferences();
  const NMTTools_IndexedDataMapOfIndexedMapOfInteger& aMAV=pPF->AloneVertices();
  //
  Standard_Boolean bHasImage;
  Standard_Integer i, j, nF, aNbS, nV, nVSD, n1, n2, iFlag;
  Standard_Integer aNbVFs, aNbAVF, aNbEFs, aNbVC, aNbE, aNbV;
  Standard_Real aU1, aU2, aTol;
  NMTTools_IndexedDataMapOfIndexedMapOfInteger aMFMV;
  TopTools_MapOfShape aMFence;
  TopTools_ListIteratorOfListOfShape aIt, aItV;
  BRep_Builder aBB;
  //
  // 1. Collect face-vertex candidates [aMFMV]
  //
  // 1.1. VFs
  aNbVFs=aVFs.Extent();
  for (i=1; i<=aNbVFs; ++i) {
    const BOPTools_VSInterference& aVS=aVFs(i);
    aVS.Indices(n1, n2);
    nF=n2;
    nV=n1;
    if (aDS.Shape(n1).ShapeType()==TopAbs_FACE) {
      nF=n1;
      nV=n2;
    }
    nVSD=pPF->FindSDVertex(nV);
    if (nVSD) {
      nV=nVSD;
    }
    //
    UpdateCandidates(nF, nV, aMFMV);
  }
  //
  // 1.2 EFs
  aNbEFs=aEFs.Extent();
  for (i=1; i<=aNbEFs; ++i) {
    const BOPTools_ESInterference& aEF=aEFs(i);
    aEF.Indices(n1, n2);
    nV=aEF.NewShape();
    if (!nV) {
      continue;
    }
    const TopoDS_Shape& aSnew=aDS.Shape(nV);
    if (aSnew.ShapeType()!=TopAbs_VERTEX) {
      continue;
    }
    //
    nF=(aDS.Shape(n1).ShapeType()==TopAbs_FACE) ? n1 : n2;
    nVSD=pPF->FindSDVertex(nV);
    if (nVSD) {
      nV=nVSD;
    }
    UpdateCandidates(nF, nV, aMFMV);
  }
  //
  aNbS=aDS.NumberOfShapesOfTheObject();
  for (nF=1; nF<=aNbS; ++nF) {
    const TopoDS_Shape& aF=aDS.Shape(nF);
    //
    if (aF.ShapeType()!=TopAbs_FACE) {
      continue;
    }
    if (!aMFence.Add(aF)) {
      continue;
    }
    //
    const TopoDS_Face& aFF=TopoDS::Face(aF);
    aTol=BRep_Tool::Tolerance(aFF);
    //
    // 1.3 FFs
    if (aMAV.Contains(nF)) {
      const TColStd_IndexedMapOfInteger& aMAVF=aMAV.FindFromKey(nF);
      aNbAVF=aMAVF.Extent();
      for (j=1; j<=aNbAVF; ++j) {
        nV=aMAVF(j);
        nVSD=pPF->FindSDVertex(nV);
        if (nVSD) {
          nV=nVSD;
        }
        //
        UpdateCandidates(nF, nV, aMFMV);
      }
    }
    //
    // 1.4 Internal vertices of the face nF
    BooleanOperations_OnceExplorer aExp(aDS);
    aExp.Init(nF, TopAbs_VERTEX);
    for (; aExp.More(); aExp.Next()) {
//.........这里部分代码省略.........

bool StdMeshers_RadialPrism_3D::Evaluate(SMESH_Mesh& aMesh,
                                         const TopoDS_Shape& aShape,
                                         MapShapeNbElems& aResMap)
{
  // 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 ) {
    std::vector<int> aResVec(SMDSEntity_Last);
    for(int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aResVec[i] = 0;
    SMESH_subMesh * sm = aMesh.GetSubMesh(aShape);
    aResMap.insert(std::make_pair(sm,aResVec));
    SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
    smError.reset( new SMESH_ComputeError(COMPERR_ALGO_FAILED,"Submesh can not be evaluated",this));
    return false;
  }

  // Associate sub-shapes of the shells
  ProjectionUtils::TShapeShapeMap shape2ShapeMap;
  if ( !ProjectionUtils::FindSubShapeAssociation( outerShell, &aMesh,
                                                  innerShell, &aMesh,
                                                  shape2ShapeMap) ) {
    std::vector<int> aResVec(SMDSEntity_Last);
    for(int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aResVec[i] = 0;
    SMESH_subMesh * sm = aMesh.GetSubMesh(aShape);
    aResMap.insert(std::make_pair(sm,aResVec));
    SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
    smError.reset( new SMESH_ComputeError(COMPERR_ALGO_FAILED,"Submesh can not be evaluated",this));
    return false;
  }

  // get info for outer shell
  int nb0d_Out=0, nb2d_3_Out=0, nb2d_4_Out=0;
  //TopTools_SequenceOfShape FacesOut;
  for (TopExp_Explorer exp(outerShell, TopAbs_FACE); exp.More(); exp.Next()) {
    //FacesOut.Append(exp.Current());
    SMESH_subMesh *aSubMesh = aMesh.GetSubMesh(exp.Current());
    MapShapeNbElemsItr anIt = aResMap.find(aSubMesh);
    std::vector<int> aVec = (*anIt).second;
    nb0d_Out += aVec[SMDSEntity_Node];
    nb2d_3_Out += Max(aVec[SMDSEntity_Triangle],aVec[SMDSEntity_Quad_Triangle]);
    nb2d_4_Out += Max(aVec[SMDSEntity_Quadrangle],aVec[SMDSEntity_Quad_Quadrangle]);
  }
  int nb1d_Out = 0;
  TopTools_MapOfShape tmpMap;
  for (TopExp_Explorer exp(outerShell, TopAbs_EDGE); exp.More(); exp.Next()) {
    if( tmpMap.Contains( exp.Current() ) )
      continue;
    tmpMap.Add( exp.Current() );
    SMESH_subMesh *aSubMesh = aMesh.GetSubMesh(exp.Current());
    MapShapeNbElemsItr anIt = aResMap.find(aSubMesh);
    std::vector<int> aVec = (*anIt).second;
    nb0d_Out += aVec[SMDSEntity_Node];
    nb1d_Out += Max(aVec[SMDSEntity_Edge],aVec[SMDSEntity_Quad_Edge]);
  }
  tmpMap.Clear();
  for (TopExp_Explorer exp(outerShell, TopAbs_VERTEX); exp.More(); exp.Next()) {
    if( tmpMap.Contains( exp.Current() ) )
      continue;
    tmpMap.Add( exp.Current() );
    nb0d_Out++;
  }

  // get info for inner shell
  int nb0d_In=0, nb2d_3_In=0, nb2d_4_In=0;
  //TopTools_SequenceOfShape FacesIn;
  for (TopExp_Explorer exp(innerShell, TopAbs_FACE); exp.More(); exp.Next()) {
    //FacesIn.Append(exp.Current());
    SMESH_subMesh *aSubMesh = aMesh.GetSubMesh(exp.Current());
    MapShapeNbElemsItr anIt = aResMap.find(aSubMesh);
    std::vector<int> aVec = (*anIt).second;
    nb0d_In += aVec[SMDSEntity_Node];
    nb2d_3_In += Max(aVec[SMDSEntity_Triangle],aVec[SMDSEntity_Quad_Triangle]);
    nb2d_4_In += Max(aVec[SMDSEntity_Quadrangle],aVec[SMDSEntity_Quad_Quadrangle]);
  }
  int nb1d_In = 0;
  tmpMap.Clear();
  bool IsQuadratic = false;
  bool IsFirst = true;
  for (TopExp_Explorer exp(innerShell, TopAbs_EDGE); exp.More(); exp.Next()) {
    if( tmpMap.Contains( exp.Current() ) )
      continue;
    tmpMap.Add( exp.Current() );
    SMESH_subMesh *aSubMesh = aMesh.GetSubMesh(exp.Current());
    MapShapeNbElemsItr anIt = aResMap.find(aSubMesh);
    std::vector<int> aVec = (*anIt).second;
    nb0d_In += aVec[SMDSEntity_Node];
    nb1d_In += Max(aVec[SMDSEntity_Edge],aVec[SMDSEntity_Quad_Edge]);
    if(IsFirst) {
      IsQuadratic = (aVec[SMDSEntity_Quad_Edge] > aVec[SMDSEntity_Edge]);
      IsFirst = false;
    }
  }
  tmpMap.Clear();
  for (TopExp_Explorer exp(innerShell, TopAbs_VERTEX); exp.More(); exp.Next()) {
//.........这里部分代码省略.........

//=======================================================================
// function: FillSameDomainFaces
// purpose:
//=======================================================================
void GEOMAlgo_Builder::FillSameDomainFaces()
{
  Standard_Boolean bIsSDF, bHasImage1, bHasImage2, bForward;
  Standard_Integer i, j, aNbFF, nF1, nF2, aNbPBInOn, aNbC, aNbSE;
  Standard_Integer aNbF1, aNbF2, i2s, aNbSD;
  TopTools_MapOfShape aMFence;
  TopTools_ListOfShape aLX1, aLX2;
  TopTools_ListIteratorOfListOfShape aItF1, aItF2;
  NMTTools_ListOfCoupleOfShape aLCS;
  //
  const NMTDS_ShapesDataStructure& aDS=*myPaveFiller->DS();
  NMTTools_PaveFiller* pPF=myPaveFiller;
  NMTDS_InterfPool* pIP=pPF->IP();
  BOPTools_CArray1OfSSInterference& aFFs=pIP->SSInterferences();
  const Handle(IntTools_Context)& aCtx= pPF->Context();
  //
  //
  //mySameDomainShapes.Clear();
  //
  // 1. For each FF find among images of faces
  //    all pairs of same domain faces (SDF) [=> aLCS]
  aNbFF=aFFs.Extent();
  for (i=1; i<=aNbFF; ++i) {
    BOPTools_SSInterference&