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

SMDS_MeshVolume* SMESH_MesherHelper::AddVolume(const SMDS_MeshNode* n1,
                                               const SMDS_MeshNode* n2,
                                               const SMDS_MeshNode* n3,
                                               const SMDS_MeshNode* n4,
                                               const int id, 
					       const bool force3d)
{
  SMESHDS_Mesh * meshDS = GetMeshDS();
  SMDS_MeshVolume* elem = 0;
  if(!myCreateQuadratic) {
    if(id)
      elem = meshDS->AddVolumeWithID(n1, n2, n3, n4, id);
    else
      elem = meshDS->AddVolume(n1, n2, n3, n4);
  }
  else {
    const SMDS_MeshNode* n12 = GetMediumNode(n1,n2,force3d);
    const SMDS_MeshNode* n23 = GetMediumNode(n2,n3,force3d);
    const SMDS_MeshNode* n31 = GetMediumNode(n3,n1,force3d);

    const SMDS_MeshNode* n14 = GetMediumNode(n1,n4,force3d);
    const SMDS_MeshNode* n24 = GetMediumNode(n2,n4,force3d);
    const SMDS_MeshNode* n34 = GetMediumNode(n3,n4,force3d);

    if(id)
      elem = meshDS->AddVolumeWithID(n1, n2, n3, n4, n12, n23, n31, n14, n24, n34, id);
    else
      elem = meshDS->AddVolume(n1, n2, n3, n4, n12, n23, n31, n14, n24, n34);
  }
  if ( mySetElemOnShape && myShapeID > 0 )
    meshDS->SetMeshElementOnShape( elem, myShapeID );

  return elem;
}

SMDS_MeshEdge* SMESH_MesherHelper::AddEdge(const SMDS_MeshNode* n1,
                                                const SMDS_MeshNode* n2,
                                                const int id,
                                                const bool force3d)
{
  SMESHDS_Mesh * meshDS = GetMeshDS();
  
  SMDS_MeshEdge* edge = 0;
  if (myCreateQuadratic) {
    const SMDS_MeshNode* n12 = GetMediumNode(n1,n2,force3d);
    if(id)
      edge = meshDS->AddEdgeWithID(n1, n2, n12, id);
    else
      edge = meshDS->AddEdge(n1, n2, n12);
  }
  else {
    if(id)
      edge = meshDS->AddEdgeWithID(n1, n2, id);
    else
      edge = meshDS->AddEdge(n1, n2);
  }

  if ( mySetElemOnShape && myShapeID > 0 )
    meshDS->SetMeshElementOnShape( edge, myShapeID );

  return edge;
}

//.........这里部分代码省略.........
          if (aEdge.IsNull())
          {
            int aGeomEdgeInd = seg.epgeominfo[j].edgenr;
            if (aGeomEdgeInd > 0 && aGeomEdgeInd <= occgeo.emap.Extent())
              aEdge = TopoDS::Edge(occgeo.emap(aGeomEdgeInd));
          }
          param = seg.epgeominfo[j].dist;
          param2 += param;
        }
        else
          param = param2 * 0.5;
        if (pind <= nbInitNod || pindMap.Contains(pind))
          continue;
        if (!aEdge.IsNull())
        {
          meshDS->SetNodeOnEdge(nodeVec.at(pind), aEdge, param);
          pindMap.Add(pind);
        }
      }
      SMDS_MeshEdge* edge;
      if (nbp < 3) // second order ?
        edge = meshDS->AddEdge(nodeVec.at(pinds[0]), nodeVec.at(pinds[1]));
      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:

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

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

//.........这里部分代码省略.........
        double param;
        if(ori)
          param = fp + dp2*(1-myLayerPositions[i]);
        else
          param = cp + dp2*myLayerPositions[i];
        meshDS->SetNodeOnEdge(node, edgeID, param);
        P = gp_Pnt( P0.X() - aVec.X()*myLayerPositions[i],
                    P0.Y() - aVec.Y()*myLayerPositions[i],
                    P0.Z() - aVec.Z()*myLayerPositions[i] );
        node = meshDS->AddNode(P.X(), P.Y(), P.Z());
        Nodes2[i] = node;
        if(!ori)
          param = fp + dp2*(1-myLayerPositions[i]);
        else
          param = cp + dp2*myLayerPositions[i];
        meshDS->SetNodeOnEdge(node, edgeID, param);
        // parameters on face
        gp_Pnt2d P2d( PC.X() + V2d.X()*myLayerPositions[i],
                      PC.Y() + V2d.Y()*myLayerPositions[i] );
        Pnts2d1.Append(P2d);
      }
      Nodes1[ myLayerPositions.size() ] = NF;
      Nodes2[ myLayerPositions.size() ] = NL;
      // create 1D elements on edge
      vector< const SMDS_MeshNode* > tmpNodes;
      tmpNodes.resize(2*Nodes1.size()+1);
      for(i=0; i<Nodes2.size(); i++)
        tmpNodes[Nodes2.size()-i-1] = Nodes2[i];
      tmpNodes[Nodes2.size()] = NC;
      for(i=0; i<Nodes1.size(); i++)
        tmpNodes[Nodes2.size()+1+i] = Nodes1[i];
      for(i=1; i<tmpNodes.size(); i++) {
        SMDS_MeshEdge* ME = myHelper->AddEdge( tmpNodes[i-1], tmpNodes[i] );
        if(ME) meshDS->SetMeshElementOnShape(ME, edgeID);
      }
      markLinEdgeAsComputedByMe( LinEdge1, aMesh.GetSubMesh( F ));
    }
  }
  else // nbe==3 or ( nbe==2 && linEdge is INTERNAL )
  {
    if (nbe==2 && LinEdge1.Orientation() == TopAbs_INTERNAL )
      LinEdge2 = LinEdge1;

    // one curve must be a part of circle and other curves must be
    // segments of line
    double fp, lp;
    Handle(Geom_Circle) aCirc = Handle(Geom_Circle)::DownCast( getCurve( CircEdge ));
    Handle(Geom_Line) aLine1 = Handle(Geom_Line)::DownCast( getCurve( LinEdge1 ));
    Handle(Geom_Line) aLine2 = Handle(Geom_Line)::DownCast( getCurve( LinEdge2 ));
    if( aLine1.IsNull() || aLine2.IsNull() ) {
      // other curve not line
      return error(COMPERR_BAD_SHAPE);
    }

    bool linEdge1Computed = false;
    if ( SMESH_subMesh* sm1 = aMesh.GetSubMesh(LinEdge1))
      if( !sm1->IsEmpty() ) {
        if( isEdgeCompitaballyMeshed( LinEdge1, aMesh.GetSubMesh(F) ))
          linEdge1Computed = true;
        else
          return error("Invalid set of hypotheses");
      }
    bool linEdge2Computed = false;
    if ( SMESH_subMesh* sm2 = aMesh.GetSubMesh(LinEdge2))
      if( !sm2->IsEmpty() ) {
        if( isEdgeCompitaballyMeshed( LinEdge2, aMesh.GetSubMesh(F)  ))

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