C++ SMESHDS_Mesh::AddNode方法代码示例

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::addNode(PyObject *args)
{
    double x,y,z;
    if (!PyArg_ParseTuple(args, "ddd",&x,&y,&z))
        return 0;

    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(PyExc_Exception, e.what());
        return 0;
    }
}

SMDS_MeshNode* SMESH_MesherHelper::AddNode(double x, double y, double z, int ID)
{
  SMESHDS_Mesh * meshDS = GetMeshDS();
  SMDS_MeshNode* node = 0;
  if ( ID )
    node = meshDS->AddNodeWithID( x, y, z, ID );
  else
    node = meshDS->AddNode( x, y, z );
  if ( mySetElemOnShape && myShapeID > 0 ) {
    switch ( myShape.ShapeType() ) {
    case TopAbs_SOLID:  meshDS->SetNodeInVolume( node, myShapeID); break;
    case TopAbs_SHELL:  meshDS->SetNodeInVolume( node, myShapeID); break;
    case TopAbs_FACE:   meshDS->SetNodeOnFace(   node, myShapeID); break;
    case TopAbs_EDGE:   meshDS->SetNodeOnEdge(   node, myShapeID); break;
    case TopAbs_VERTEX: meshDS->SetNodeOnVertex( node, myShapeID); break;
    default: ;
    }
  }
  return node;
}

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;

}

//.........这里部分代码省略.........
  SMESHDS_Mesh* meshDS = _mesh->GetMeshDS();
  bool isOK = ( !err && (_isVolume ? (nbVol > 0) : (nbFac > 0)) );
  if ( true /*isOK*/ ) // get whatever built
  {
    // map of nodes assigned to submeshes
    NCollection_Map<int> pindMap;
    // create and insert nodes into nodeVec
    nodeVec.resize( nbNod + 1 );
    int i;
    for (i = nbInitNod+1; i <= nbNod /*&& isOK*/; ++i )
    {
      const netgen::MeshPoint& ngPoint = ngMesh->Point(i);
      SMDS_MeshNode* node = NULL;
      bool newNodeOnVertex = false;
      TopoDS_Vertex aVert;
      if (i-nbInitNod <= occgeo.vmap.Extent())
      {
        // point on vertex
        aVert = TopoDS::Vertex(occgeo.vmap(i-nbInitNod));
        SMESHDS_SubMesh * submesh = meshDS->MeshElements(aVert);
        if (submesh)
        {
          SMDS_NodeIteratorPtr it = submesh->GetNodes();
          if (it->more())
          {
            node = const_cast<SMDS_MeshNode*> (it->next());
            pindMap.Add(i);
          }
        }
        if (!node)
          newNodeOnVertex = true;
      }
      if (!node)
        node = meshDS->AddNode(ngPoint.X(), ngPoint.Y(), ngPoint.Z());
      if (!node)
      {
        MESSAGE("Cannot create a mesh node");
        if ( !comment.size() ) comment << "Cannot create a mesh node";
        nbSeg = nbFac = nbVol = isOK = 0;
        break;
      }
      nodeVec.at(i) = node;
      if (newNodeOnVertex)
      {
        // point on vertex
        meshDS->SetNodeOnVertex(node, aVert);
        pindMap.Add(i);
      }
    }

    // create mesh segments along geometric edges
    NCollection_Map<Link> linkMap;
    for (i = nbInitSeg+1; i <= nbSeg/* && isOK*/; ++i )
    {
      const netgen::Segment& seg = ngMesh->LineSegment(i);
      Link link(seg.p1, seg.p2);
      if (linkMap.Contains(link))
        continue;
      linkMap.Add(link);
      TopoDS_Edge aEdge;
      int pinds[3] = { seg.p1, seg.p2, seg.pmid };
      int nbp = 0;
      double param2 = 0;
      for (int j=0; j < 3; ++j)
      {
        int pind = pinds[j];

//.........这里部分代码省略.........
                 Netgen_triangle[2] != Netgen_triangle[1] ))
        {
            Ng_AddSurfaceElement(Netgen_mesh, NG_TRIG, Netgen_triangle);
        }
    }

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

    Ng_Meshing_Parameters Netgen_param;

    Netgen_param.secondorder = Netgen_param2ndOrder;
    Netgen_param.fineness = Netgen_paramFine;
    Netgen_param.maxh = Netgen_paramSize;

    Ng_Result status;

    try {
#if (OCC_VERSION_MAJOR << 16 | OCC_VERSION_MINOR << 8 | OCC_VERSION_MAINTENANCE) > 0x060100
        OCC_CATCH_SIGNALS;
#endif
        status = Ng_GenerateVolumeMesh(Netgen_mesh, &Netgen_param);
    }
    catch (Standard_Failure& exc) {
        error(COMPERR_OCC_EXCEPTION, exc.GetMessageString());
        status = NG_VOLUME_FAILURE;
    }
    catch (...) {
        error("Exception in Ng_GenerateVolumeMesh()");
        status = NG_VOLUME_FAILURE;
    }
    if ( GetComputeError()->IsOK() ) {
        switch ( status ) {
        case NG_SURFACE_INPUT_ERROR:
            error( status, "NG_SURFACE_INPUT_ERROR");
        case NG_VOLUME_FAILURE:
            error( status, "NG_VOLUME_FAILURE");
        case NG_STL_INPUT_ERROR:
            error( status, "NG_STL_INPUT_ERROR");
        case NG_SURFACE_FAILURE:
            error( status, "NG_SURFACE_FAILURE");
        case NG_FILE_NOT_FOUND:
            error( status, "NG_FILE_NOT_FOUND");
        };
    }

    int Netgen_NbOfNodesNew = Ng_GetNP(Netgen_mesh);

    int Netgen_NbOfTetra = Ng_GetNE(Netgen_mesh);

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

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

    bool isOK = ( /*status == NG_OK &&*/ Netgen_NbOfTetra > 0 );// get whatever built
    if ( isOK )
    {
        // vector of nodes in which node index == netgen ID
        vector< const SMDS_MeshNode* > nodeVec ( Netgen_NbOfNodesNew + 1 );
        // insert old nodes into nodeVec
        for ( n_id = nodeToNetgenID.begin(); n_id != nodeToNetgenID.end(); ++n_id ) {
            nodeVec.at( n_id->second ) = n_id->first;
        }
        // create and insert new nodes into nodeVec
        int nodeIndex = Netgen_NbOfNodes + 1;
        int shapeID = meshDS->ShapeToIndex( aShape );
        for ( ; nodeIndex <= Netgen_NbOfNodesNew; ++nodeIndex )
        {
            Ng_GetPoint( Netgen_mesh, nodeIndex, Netgen_point );
            SMDS_MeshNode * node = meshDS->AddNode(Netgen_point[0],
                                                   Netgen_point[1],
                                                   Netgen_point[2]);
            meshDS->SetNodeInVolume(node, shapeID);
            nodeVec.at(nodeIndex) = node;
        }

        // create tetrahedrons
        for ( int elemIndex = 1; elemIndex <= Netgen_NbOfTetra; ++elemIndex )
        {
            Ng_GetVolumeElement(Netgen_mesh, elemIndex, Netgen_tetrahedron);
            SMDS_MeshVolume * elt = myTool->AddVolume (nodeVec.at( Netgen_tetrahedron[0] ),
                                    nodeVec.at( Netgen_tetrahedron[1] ),
                                    nodeVec.at( Netgen_tetrahedron[2] ),
                                    nodeVec.at( Netgen_tetrahedron[3] ));
            meshDS->SetMeshElementOnShape(elt, shapeID );
        }
    }

    Ng_DeleteMesh(Netgen_mesh);
    Ng_Exit();

    NETGENPlugin_Mesher::RemoveTmpFiles();

    return (status == NG_OK);
}

/*!
 * 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;
//.........这里部分代码省略.........

//.........这里部分代码省略.........
  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 += ": ";
      comment += ex.GetMessageString();
    }
    error(COMPERR_OCC_EXCEPTION, comment);
  }
  catch (NgException exc) {
    error( SMESH_Comment("NgException: ") << exc.What() );
  }
  catch (...) {
    error(COMPERR_EXCEPTION,"Exception in netgen::OCCGenerateMesh()");
  }

  // ----------------------------------------------------
  // Fill the SMESHDS with the generated nodes and faces
  // ----------------------------------------------------

  int nbNodes = ngMesh->GetNP();
  int nbFaces = ngMesh->GetNSE();

  int nbInputNodes = nodeVec.size();
  nodeVec.resize( nbNodes, 0 );

  // add nodes
  for ( int i = nbInputNodes + 1; i <= nbNodes; ++i )
  {
    const MeshPoint& ngPoint = ngMesh->Point(i);
    SMDS_MeshNode * node = meshDS->AddNode(ngPoint(0), ngPoint(1), ngPoint(2));
    nodeVec[ i-1 ] = node;
  }

  // create faces
  bool reverse = ( aShape.Orientation() == TopAbs_REVERSED );
  for ( int i = 1; i <= nbFaces ; ++i )
  {
    const Element2d& elem = ngMesh->SurfaceElement(i);
    vector<const SMDS_MeshNode*> nodes( elem.GetNP() );
    for (int j=1; j <= elem.GetNP(); ++j)
    {
      int pind = elem.PNum(j);
      const SMDS_MeshNode* node = nodeVec.at(pind-1);
      if ( reverse )
        nodes[ nodes.size()-j ] = node;
      else
        nodes[ j-1 ] = node;
      if ( node->GetPosition()->GetTypeOfPosition() == SMDS_TOP_3DSPACE )
      {
        const PointGeomInfo& pgi = elem.GeomInfoPi(j);
        meshDS->SetNodeOnFace((SMDS_MeshNode*)node, faceID, pgi.u, pgi.v);
      }
    }
    SMDS_MeshFace* face = 0;
    if ( elem.GetType() == TRIG )
      face = helper.AddFace(nodes[0],nodes[1],nodes[2]);
    else
      face = helper.AddFace(nodes[0],nodes[1],nodes[2],nodes[3]);
  }

  Ng_DeleteMesh((nglib::Ng_Mesh*)ngMesh);
  Ng_Exit();

  NETGENPlugin_Mesher::RemoveTmpFiles();

  return !err;
}

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 );
//.........这里部分代码省略.........

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

  SMDS_VolumeTool volTool;
  SMESH_MesherHelper helper( *tgtMesh );
  helper.IsQuadraticSubMesh( aShape );

  SMESHDS_SubMesh* srcSMDS = srcSubMesh->GetSubMeshDS();
  SMDS_ElemIteratorPtr volIt = srcSMDS->GetElements();
  while ( volIt->more() ) // loop on source volumes
  {
    const SMDS_MeshElement* srcVol = volIt->next();
    if ( !srcVol || srcVol->GetType() != SMDSAbs_Volume )
        continue;
    int nbNodes = srcVol->NbNodes();
    SMDS_VolumeTool::VolumeType  volType = volTool.GetType( nbNodes );
    if ( srcVol->IsQuadratic() )
      nbNodes = volTool.NbCornerNodes( volType );

    // Find or create a new tgt node for each node of a src volume

    vector< const SMDS_MeshNode* > nodes( nbNodes );
    for ( int i = 0; i < nbNodes; ++i )
    {
      const SMDS_MeshNode* srcNode = srcVol->GetNode( i );
      const SMDS_MeshNode* tgtNode = 0;
      TNodeNodeMap::iterator sN_tN = src2tgtNodeMap.find( srcNode );
      if ( sN_tN != src2tgtNodeMap.end() ) // found
      {
        tgtNode = sN_tN->second;
      }
      else // Create a new tgt node
      {
        // compute normalized parameters of source node in srcBlock
        gp_Pnt srcCoord = gpXYZ( srcNode );
        gp_XYZ srcParam;
        if ( !srcBlock.ComputeParameters( srcCoord, srcParam ))
          return error(SMESH_Comment("Can't compute normalized parameters ")
                       << "for source node " << srcNode->GetID());
        // compute coordinates of target node by srcParam
        gp_XYZ tgtXYZ;
        if ( !tgtBlock.ShellPoint( srcParam, tgtXYZ ))
          return error("Can't compute coordinates by normalized parameters");
        // add node
        SMDS_MeshNode* newNode = tgtMeshDS->AddNode( tgtXYZ.X(), tgtXYZ.Y(), tgtXYZ.Z() );
        tgtMeshDS->SetNodeInVolume( newNode, helper.GetSubShapeID() );
        tgtNode = newNode;
        src2tgtNodeMap.insert( make_pair( srcNode, tgtNode ));
      }
      nodes[ i ] = tgtNode;
    }

    // Create a new volume

    SMDS_MeshVolume * tgtVol = 0;
    int id = 0, force3d = false;
    switch ( volType ) {
    case SMDS_VolumeTool::TETRA     :
    case SMDS_VolumeTool::QUAD_TETRA:
      tgtVol = helper.AddVolume( nodes[0],
                                 nodes[1],
                                 nodes[2],
                                 nodes[3], id, force3d); break;
    case SMDS_VolumeTool::PYRAM     :
    case SMDS_VolumeTool::QUAD_PYRAM:
      tgtVol = helper.AddVolume( nodes[0],
                                 nodes[1],
                                 nodes[2],
                                 nodes[3],
                                 nodes[4], id, force3d); break;
    case SMDS_VolumeTool::PENTA     :
    case SMDS_VolumeTool::QUAD_PENTA:
      tgtVol = helper.AddVolume( nodes[0],
                                 nodes[1],
                                 nodes[2],
                                 nodes[3],
                                 nodes[4],
                                 nodes[5], id, force3d); break;
    case SMDS_VolumeTool::HEXA      :
    case SMDS_VolumeTool::QUAD_HEXA :
      tgtVol = helper.AddVolume( nodes[0],
                                 nodes[1],
                                 nodes[2],
                                 nodes[3],
                                 nodes[4],
                                 nodes[5],
                                 nodes[6],
                                 nodes[7], id, force3d); break;
    default: // polyhedron
      const SMDS_PolyhedralVolumeOfNodes * poly =
        dynamic_cast<const SMDS_PolyhedralVolumeOfNodes*>( srcVol );
      if ( !poly )
        RETURN_BAD_RESULT("Unexpected volume type");
      tgtVol = tgtMeshDS->AddPolyhedralVolume( nodes, poly->GetQuanities() );
    }
    if ( tgtVol ) {
      tgtMeshDS->SetMeshElementOnShape( tgtVol, helper.GetSubShapeID() );
    }
  } // loop on volumes of src shell

  return true;
}

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 )
  {
//.........这里部分代码省略.........

//.........这里部分代码省略.........
            {
                err = 1;
                str << "Exception in  netgen::OCCGenerateMesh()"
                    << " at " << netgen::multithread.task
                    << ": " << ex.DynamicType()->Name();
                if ( ex.GetMessageString() && strlen( ex.GetMessageString() ))
                    str << ": " << ex.GetMessageString();
            }
            catch (...) {
                err = 1;
                str << "Exception in  netgen::OCCGenerateMesh()"
                    << " at " << netgen::multithread.task;
            }
            if ( err )
            {
                if ( aMesher.FixFaceMesh( occgeom, *ngMesh, 1 ))
                    break;
                if ( iLoop == LOC_SIZE )
                {
                    netgen::mparam.minh = netgen::mparam.maxh;
                    netgen::mparam.maxh = 0;
                    for ( int iW = 0; iW < wires.size(); ++iW )
                    {
                        StdMeshers_FaceSidePtr wire = wires[ iW ];
                        const vector<UVPtStruct>& uvPtVec = wire->GetUVPtStruct();
                        for ( size_t iP = 1; iP < uvPtVec.size(); ++iP )
                        {
                            SMESH_TNodeXYZ   p( uvPtVec[ iP ].node );
                            netgen::Point3d np( p.X(),p.Y(),p.Z());
                            double segLen = p.Distance( uvPtVec[ iP-1 ].node );
                            double   size = ngMesh->GetH( np );
                            netgen::mparam.minh = Min( netgen::mparam.minh, size );
                            netgen::mparam.maxh = Max( netgen::mparam.maxh, segLen );
                        }
                    }
                    //cerr << "min " << netgen::mparam.minh << " max " << netgen::mparam.maxh << endl;
                    netgen::mparam.minh *= 0.9;
                    netgen::mparam.maxh *= 1.1;
                    continue;
                }
                else
                {
                    faceErr.reset( new SMESH_ComputeError( COMPERR_ALGO_FAILED, str ));
                }
            }


            // ----------------------------------------------------
            // Fill the SMESHDS with the generated nodes and faces
            // ----------------------------------------------------

            int nbNodes = ngMesh->GetNP();
            int nbFaces = ngMesh->GetNSE();

            int nbInputNodes = nodeVec.size()-1;
            nodeVec.resize( nbNodes+1, 0 );

            // add nodes
            for ( int ngID = nbInputNodes + 1; ngID <= nbNodes; ++ngID )
            {
                const MeshPoint& ngPoint = ngMesh->Point( ngID );
                SMDS_MeshNode * node = meshDS->AddNode(ngPoint(0), ngPoint(1), ngPoint(2));
                nodeVec[ ngID ] = node;
            }

            // create faces
            int i,j;
            vector<const SMDS_MeshNode*> nodes;
            for ( i = 1; i <= nbFaces ; ++i )
            {
                const Element2d& elem = ngMesh->SurfaceElement(i);
                nodes.resize( elem.GetNP() );
                for (j=1; j <= elem.GetNP(); ++j)
                {
                    int pind = elem.PNum(j);
                    if ( pind < 1 )
                        break;
                    nodes[ j-1 ] = nodeVec[ pind ];
                    if ( nodes[ j-1 ]->GetPosition()->GetTypeOfPosition() == SMDS_TOP_3DSPACE )
                    {
                        const PointGeomInfo& pgi = elem.GeomInfoPi(j);
                        meshDS->SetNodeOnFace( nodes[ j-1 ], faceID, pgi.u, pgi.v);
                    }
                }
                if ( j > elem.GetNP() )
                {
                    SMDS_MeshFace* face = 0;
                    if ( elem.GetType() == TRIG )
                        face = helper.AddFace(nodes[0],nodes[1],nodes[2]);
                    else
                        face = helper.AddFace(nodes[0],nodes[1],nodes[2],nodes[3]);
                }
            }

            break;
        } // two attempts
    } // loop on FACEs

    return true;
}

//.........这里部分代码省略.........
        GetPoint(p00z, 0, 0, k, nbx, nby, nbz, np, meshDS);
        GetPoint(p01z, 0, nby - 1, k, nbx, nby, nbz, np, meshDS);
        GetPoint(p10z, nbx - 1, 0, k, nbx, nby, nbz, np, meshDS);
        GetPoint(p11z, nbx - 1, nby - 1, k, nbx, nby, nbz, np, meshDS);

        // 12 points on faces
        GetPoint(pxy0, i, j, 0, nbx, nby, nbz, np, meshDS);
        GetPoint(pxy1, i, j, nbz - 1, nbx, nby, nbz, np, meshDS);
        GetPoint(px0z, i, 0, k, nbx, nby, nbz, np, meshDS);
        GetPoint(px1z, i, nby - 1, k, nbx, nby, nbz, np, meshDS);
        GetPoint(p0yz, 0, j, k, nbx, nby, nbz, np, meshDS);
        GetPoint(p1yz, nbx - 1, j, k, nbx, nby, nbz, np, meshDS);

        int ijk = k * nbx * nby + j * nbx + i;
        double x = double (i) / double (nbx - 1);       // *** seulement
        double y = double (j) / double (nby - 1);       // *** maillage
        double z = double (k) / double (nbz - 1);       // *** regulier

        Pt3 X;
        for (int i = 0; i < 3; i++) {
          X[i] = (1 - x) * p0yz[i] + x * p1yz[i]
                 + (1 - y) * px0z[i] + y * px1z[i]
                 + (1 - z) * pxy0[i] + z * pxy1[i]
                 - (1 - x) * ((1 - y) * p00z[i] + y * p01z[i])
                 - x * ((1 - y) * p10z[i] + y * p11z[i])
                 - (1 - y) * ((1 - z) * px00[i] + z * px01[i])
                 - y * ((1 - z) * px10[i] + z * px11[i])
                 - (1 - z) * ((1 - x) * p0y0[i] + x * p1y0[i])
                 - z * ((1 - x) * p0y1[i] + x * p1y1[i])
                 + (1 - x) * ((1 - y) * ((1 - z) * p000[i] + z * p001[i])
                 + y * ((1 - z) * p010[i] + z * p011[i]))
                 + x * ((1 - y) * ((1 - z) * p100[i] + z * p101[i])
                 + y * ((1 - z) * p110[i] + z * p111[i]));
        }

        SMDS_MeshNode * node = meshDS->AddNode(X[0], X[1], X[2]);
        np[ijk].node = node;
        meshDS->SetNodeInVolume(node, shapeID);
      }
    }
  }

  // find orientation of furute volumes according to MED convention
  vector< bool > forward( nbx * nby );
  SMDS_VolumeTool vTool;
  for (int i = 0; i < nbx - 1; i++) {
    for (int j = 0; j < nby - 1; j++) {
      int n1 = j * nbx + i;
      int n2 = j * nbx + i + 1;
      int n3 = (j + 1) * nbx + i + 1;
      int n4 = (j + 1) * nbx + i;
      int n5 = nbx * nby + j * nbx + i;
      int n6 = nbx * nby + j * nbx + i + 1;
      int n7 = nbx * nby + (j + 1) * nbx + i + 1;
      int n8 = nbx * nby + (j + 1) * nbx + i;

      SMDS_VolumeOfNodes tmpVol (np[n1].node,np[n2].node,np[n3].node,np[n4].node,
                                 np[n5].node,np[n6].node,np[n7].node,np[n8].node);
      vTool.Set( &tmpVol );
      forward[ n1 ] = vTool.IsForward();
    }
  }

  //2.1 - for each node of the cube (less 3 *1 Faces):
  //      - store hexahedron in SMESHDS
  MESSAGE("Storing hexahedron into the DS");
  for (int i = 0; i < nbx - 1; i++) {
    for (int j = 0; j < nby - 1; j++) {
      bool isForw = forward.at( j * nbx + i );
      for (int k = 0; k < nbz - 1; k++) {
        int n1 = k * nbx * nby + j * nbx + i;
        int n2 = k * nbx * nby + j * nbx + i + 1;
        int n3 = k * nbx * nby + (j + 1) * nbx + i + 1;
        int n4 = k * nbx * nby + (j + 1) * nbx + i;
        int n5 = (k + 1) * nbx * nby + j * nbx + i;
        int n6 = (k + 1) * nbx * nby + j * nbx + i + 1;
        int n7 = (k + 1) * nbx * nby + (j + 1) * nbx + i + 1;
        int n8 = (k + 1) * nbx * nby + (j + 1) * nbx + i;

        SMDS_MeshVolume * elt;
        if ( isForw ) {
          elt = aTool.AddVolume(np[n1].node, np[n2].node,
                                np[n3].node, np[n4].node,
                                np[n5].node, np[n6].node,
                                np[n7].node, np[n8].node);
        }
        else {
          elt = aTool.AddVolume(np[n1].node, np[n4].node,
                                np[n3].node, np[n2].node,
                                np[n5].node, np[n8].node,
                                np[n7].node, np[n6].node);
        }
        
        meshDS->SetMeshElementOnShape(elt, shapeID);
      }
    }
  }
  if ( np ) delete [] np;
  return ClearAndReturn( aQuads, true );
}