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

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;
}

//.........这里部分代码省略.........
      else
        edge = meshDS->AddEdge(nodeVec.at(pinds[0]), nodeVec.at(pinds[1]),
                                nodeVec.at(pinds[2]));
      if (!edge)
      {
        if ( !comment.size() ) comment << "Cannot create a mesh edge";
        MESSAGE("Cannot create a mesh edge");
        nbSeg = nbFac = nbVol = isOK = 0;
        break;
      }
      if (!aEdge.IsNull())
        meshDS->SetMeshElementOnShape(edge, aEdge);
    }

    // create mesh faces along geometric faces
    for (i = nbInitFac+1; i <= nbFac/* && isOK*/; ++i )
    {
      const netgen::Element2d& elem = ngMesh->SurfaceElement(i);
      int aGeomFaceInd = elem.GetIndex();
      TopoDS_Face aFace;
      if (aGeomFaceInd > 0 && aGeomFaceInd <= occgeo.fmap.Extent())
        aFace = TopoDS::Face(occgeo.fmap(aGeomFaceInd));
      vector<SMDS_MeshNode*> nodes;
      for (int j=1; j <= elem.GetNP(); ++j)
      {
        int pind = elem.PNum(j);
        SMDS_MeshNode* node = nodeVec.at(pind);
        nodes.push_back(node);
        if (pind <= nbInitNod || pindMap.Contains(pind))
          continue;
        if (!aFace.IsNull())
        {
          const netgen::PointGeomInfo& pgi = elem.GeomInfoPi(j);
          meshDS->SetNodeOnFace(node, aFace, pgi.u, pgi.v);
          pindMap.Add(pind);
        }
      }
      SMDS_MeshFace* face = NULL;
      switch (elem.GetType())
      {
      case netgen::TRIG:
        face = meshDS->AddFace(nodes[0],nodes[1],nodes[2]);
        break;
      case netgen::QUAD:
        face = meshDS->AddFace(nodes[0],nodes[1],nodes[2],nodes[3]);
        break;
      case netgen::TRIG6:
        face = meshDS->AddFace(nodes[0],nodes[1],nodes[2],nodes[5],nodes[3],nodes[4]);
        break;
      case netgen::QUAD8:
        face = meshDS->AddFace(nodes[0],nodes[1],nodes[2],nodes[3],
                               nodes[4],nodes[7],nodes[5],nodes[6]);
        break;
      default:
        MESSAGE("NETGEN created a face of unexpected type, ignoring");
        continue;
      }
      if (!face)
      {
        if ( !comment.size() ) comment << "Cannot create a mesh face";
        MESSAGE("Cannot create a mesh face");
        nbSeg = nbFac = nbVol = isOK = 0;
        break;
      }
      if (!aFace.IsNull())
        meshDS->SetMeshElementOnShape(face, aFace);

/*!
 * 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_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;
}