smesh/src/StdMeshers/StdMeshers_Helper.cxx

// File:      StdMeshers_Helper.cxx
// Created:   15.02.06 15:22:41
// Author:    Sergey KUUL
// Copyright: Open CASCADE 2006


#include "StdMeshers_Helper.hxx"

#include "TopTools_MapOfShape.hxx"
#include "BRepTools.hxx"
#include "BRepTools_WireExplorer.hxx"
#include "SMDS_FacePosition.hxx"
#include "SMDS_EdgePosition.hxx"
#include "SMESH_subMesh.hxx"
#include "Geom2d_Curve.hxx"
#include "Geom_Curve.hxx"
#include "Geom_Surface.hxx"
#include "BRep_Tool.hxx"
#include "gp_Pnt2d.hxx"


//=======================================================================
//function : CheckShape
//purpose  : 
//=======================================================================

bool StdMeshers_Helper::IsQuadraticSubMesh(const TopoDS_Shape& aSh,
                                           const bool QuadMode)
{
  SMESHDS_Mesh* meshDS = GetMesh()->GetMeshDS();
  myShapeID = meshDS->ShapeToIndex(aSh);
  myCreateQuadratic = false;
  if(QuadMode) {
    // we can create quadratic elements only if each elements
    // created on given shape is quadratic
    // also we have to fill myNLinkNodeMap
    myCreateQuadratic = true;
    if(aSh.ShapeType()!=TopAbs_FACE) {
      for (TopExp_Explorer exp(aSh, TopAbs_FACE); exp.More() && myCreateQuadratic; exp.Next()) {
        const TopoDS_Face& F = TopoDS::Face(exp.Current());
        SMDS_ElemIteratorPtr itf = GetMesh()->GetSubMesh(F)->GetSubMeshDS()->GetElements();
        while(itf->more()) {
          const SMDS_MeshElement* f = itf->next();
          if( f->GetType()==SMDSAbs_Face && !f->IsQuadratic() ) {
            myCreateQuadratic = false;
            break;
          }
          SMDS_ElemIteratorPtr itn = f->nodesIterator();
          if(f->NbNodes()==6) {
            const SMDS_MeshNode* Ns[6];
            int i = 0;
            while(itn->more()) {
              Ns[i++] = static_cast<const SMDS_MeshNode*>( itn->next() );
            }
            AddNLinkNode(Ns[0],Ns[1],Ns[3]);
            AddNLinkNode(Ns[1],Ns[2],Ns[4]);
            AddNLinkNode(Ns[2],Ns[0],Ns[5]);
          }
          else if(f->NbNodes()==8) {
            const SMDS_MeshNode* Ns[8];
            int i = 0;
            while(itn->more()) {
              Ns[i++] = static_cast<const SMDS_MeshNode*>( itn->next() );
            }
            AddNLinkNode(Ns[0],Ns[1],Ns[4]);
            AddNLinkNode(Ns[1],Ns[2],Ns[5]);
            AddNLinkNode(Ns[2],Ns[3],Ns[6]);
            AddNLinkNode(Ns[3],Ns[0],Ns[7]);
          }
        }
      }
    }
    else {
      TopTools_MapOfShape aMap;
      // check edges
      const TopoDS_Face& F = TopoDS::Face(aSh);
      const TopoDS_Wire& W = BRepTools::OuterWire(F);
      BRepTools_WireExplorer wexp (W, F);
      for (wexp.Init(W, F); wexp.More() && myCreateQuadratic; wexp.Next()) {
        const TopoDS_Edge& E = wexp.Current();
        if(aMap.Contains(E))
          continue;
        aMap.Add(E);
        SMDS_ElemIteratorPtr it = GetMesh()->GetSubMesh(E)->GetSubMeshDS()->GetElements();
        while(it->more()) {
          const SMDS_MeshElement* e = it->next();
          if( e->GetType()==SMDSAbs_Edge && !e->IsQuadratic() ) {
            myCreateQuadratic = false;
            break;
          }
          // fill NLinkNodeMap
          SMDS_ElemIteratorPtr nodeIt = e->nodesIterator();
          const SMDS_MeshNode* n1 = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
          const SMDS_MeshNode* n2 = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
          const SMDS_MeshNode* n3 = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
          NLink link(( n1 < n2 ? n1 : n2 ), ( n1 < n2 ? n2 : n1 ));
          myNLinkNodeMap.insert(NLinkNodeMap::value_type(link,n3));
          myNLinkNodeMap[link] = n3;
        }
      }
    }
  }

  if(!myCreateQuadratic) {
    myNLinkNodeMap.clear();
  }

  return myCreateQuadratic;
}


//=======================================================================
//function : IsMedium
//purpose  : 
//=======================================================================

bool StdMeshers_Helper::IsMedium(const SMDS_MeshNode* n)
{
  SMDS_ElemIteratorPtr it = n->GetInverseElementIterator();
  while (it->more()) {
    const SMDS_MeshElement* elem = it->next();
    return elem->IsMediumNode(n);
  }
  return false;
}


//=======================================================================
//function : AddNLinkNode
//purpose  : 
//=======================================================================
/*!
 * Auxilary function for filling myNLinkNodeMap
 */
void StdMeshers_Helper::AddNLinkNode(const SMDS_MeshNode* n1,
                                     const SMDS_MeshNode* n2,
                                     const SMDS_MeshNode* n12)
{
  NLink link(( n1 < n2 ? n1 : n2 ), ( n1 < n2 ? n2 : n1 ));
  ItNLinkNode itLN = myNLinkNodeMap.find( link );
  if ( itLN == myNLinkNodeMap.end() ) {
    // add new record to map
    myNLinkNodeMap.insert(NLinkNodeMap::value_type(link,n12));
  }
}


//=======================================================================
//function : GetNodeUV
//purpose  : 
//=======================================================================
/*!
 * Auxilary function for GetMediumNode()
 */
gp_XY StdMeshers_Helper::GetNodeUV(const TopoDS_Face& F,
                                   const SMDS_MeshNode* n)
{
  gp_Pnt2d p2d;
  const SMDS_PositionPtr Pos = n->GetPosition();
  if(Pos->GetTypeOfPosition()==SMDS_TOP_FACE) {
//cout<<"face"<<endl;
    // node has position on face
    const SMDS_FacePosition* fpos =
      static_cast<const SMDS_FacePosition*>(n->GetPosition().get());
    p2d = gp_Pnt2d(fpos->GetUParameter(),fpos->GetVParameter());
//cout<<"fpos->GetUParameter()="<<fpos->GetUParameter()<<" fpos->GetVParameter()="<<fpos->GetVParameter()<<endl;
  }
  else if(Pos->GetTypeOfPosition()==SMDS_TOP_EDGE) {
    // node has position on edge => it is needed to find
    // corresponding edge from face, get pcurve for this
    // edge and recieve value from this pcurve
    const SMDS_EdgePosition* epos =
      static_cast<const SMDS_EdgePosition*>(n->GetPosition().get());
    double param = epos->GetUParameter();
    SMESHDS_Mesh* meshDS = GetMesh()->GetMeshDS();
    int edgeID = Pos->GetShapeId();
    const TopoDS_Edge& E = TopoDS::Edge(meshDS->IndexToShape(edgeID));
    double f, l;
    Handle(Geom2d_Curve) C2d = BRep_Tool::CurveOnSurface(E, F, f, l);
    p2d = C2d->Value(param);
  }
  else { // vertex position
    SMESHDS_Mesh * meshDS = GetMesh()->GetMeshDS();
    const TopoDS_Wire& W = BRepTools::OuterWire(F);
    BRepTools_WireExplorer wexp (W, F);
    TopoDS_Edge E;
    for (wexp.Init(W, F); wexp.More(); wexp.Next()) {
      E = wexp.Current();
      TopoDS_Vertex V1 = TopExp::FirstVertex(E);
      TopoDS_Vertex V2 = TopExp::LastVertex(E);
      if( meshDS->ShapeToIndex(V1) != Pos->GetShapeId() &&
          meshDS->ShapeToIndex(V2) != Pos->GetShapeId() ) continue;
      double f2, l2;
      Handle(Geom2d_Curve) C2d = BRep_Tool::CurveOnSurface(E, F, f2, l2);
      double f, l;
      TopLoc_Location L;
      Handle(Geom_Curve) C = BRep_Tool::Curve(E, L, f, l);
      gp_Pnt P = C->Value(f);                                                                          
      if(!L.IsIdentity()) P = P.Transformed(L.Transformation());
      double tol = BRep_Tool::Tolerance(E);
      double dx = n->X() - P.X();
      double dy = n->Y() - P.Y();
      double dz = n->Z() - P.Z();
      double dist = sqrt(dx*dx + dy*dy + dz*dz);
      if(dist<tol) {
        p2d = C2d->Value(f2);
      }
      else {
        p2d = C2d->Value(l2);
      }
    }
  }
  return p2d.XY();
}


//=======================================================================
//function : GetMediumNode
//purpose  : 
//=======================================================================
/*!
 * Special function for search or creation medium node
 */
const SMDS_MeshNode* StdMeshers_Helper::GetMediumNode(const SMDS_MeshNode* n1,
                                                      const SMDS_MeshNode* n2,
                                                      const bool force3d)
{
//cout<<"n1: "<<n1;
//cout<<"n2: "<<n2;
  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 = GetMesh()->GetMeshDS();
    if(!force3d) {
      // we try to create medium node using UV parameters of
      // nodes, else - medium between corresponding 3d points
      const SMDS_PositionPtr Pos1 = n1->GetPosition();
      const SMDS_PositionPtr Pos2 = n2->GetPosition();
      int faceID = -1;
      if( Pos1->GetTypeOfPosition()==SMDS_TOP_FACE ) {
        faceID = Pos1->GetShapeId();
      }
      else if( Pos2->GetTypeOfPosition()==SMDS_TOP_FACE ) {
        faceID = Pos2->GetShapeId();
      }
      if(faceID>-1) {
        TopoDS_Face F = TopoDS::Face(meshDS->IndexToShape(faceID));
        gp_XY p1 = GetNodeUV(F,n1);
        gp_XY p2 = GetNodeUV(F,n2);
        double u = (p1.X()+p2.X())/2.;
        double v = (p1.Y()+p2.Y())/2.;
        Handle(Geom_Surface) S = BRep_Tool::Surface(F);
        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;
      }
    }
    // 3d variant
    double x = ( n1->X() + n2->X() )/2.;
    double y = ( n1->Y() + n2->Y() )/2.;
    double z = ( n1->Z() + n2->Z() )/2.;
    n12 = meshDS->AddNode(x,y,z);
    meshDS->SetNodeInVolume(n12, myShapeID);
    myNLinkNodeMap.insert(NLinkNodeMap::value_type(link,n12));
    return n12;
  }
}


//=======================================================================
//function : AddFace
//purpose  : 
//=======================================================================
/*!
 * Special function for creation quadratic triangle
 */
SMDS_MeshFace* StdMeshers_Helper::AddFace(const SMDS_MeshNode* n1,
                                          const SMDS_MeshNode* n2,
                                          const SMDS_MeshNode* n3)
{
  SMESHDS_Mesh * meshDS = GetMesh()->GetMeshDS();
  if(!myCreateQuadratic) {
    return  meshDS->AddFace(n1, n2, n3);
  }

  const SMDS_MeshNode* n12 = GetMediumNode(n1,n2,false);
  const SMDS_MeshNode* n23 = GetMediumNode(n2,n3,false);
  const SMDS_MeshNode* n31 = GetMediumNode(n3,n1,false);

  return  meshDS->AddFace(n1, n2, n3, n12, n23, n31);
}


//=======================================================================
//function : AddFace
//purpose  : 
//=======================================================================
/*!
 * Special function for creation quadratic quadrangle
 */
SMDS_MeshFace* StdMeshers_Helper::AddFace(const SMDS_MeshNode* n1,
                                          const SMDS_MeshNode* n2,
                                          const SMDS_MeshNode* n3,
                                          const SMDS_MeshNode* n4)
{
  SMESHDS_Mesh * meshDS = GetMesh()->GetMeshDS();
  if(!myCreateQuadratic) {
    return  meshDS->AddFace(n1, n2, n3, n4);
  }

  const SMDS_MeshNode* n12 = GetMediumNode(n1,n2,false);
  const SMDS_MeshNode* n23 = GetMediumNode(n2,n3,false);
  const SMDS_MeshNode* n34 = GetMediumNode(n3,n4,false);
  const SMDS_MeshNode* n41 = GetMediumNode(n4,n1,false);

  return  meshDS->AddFace(n1, n2, n3, n4, n12, n23, n34, n41);
}


//=======================================================================
//function : AddVolume
//purpose  : 
//=======================================================================
/*!
 * Special function for creation quadratic volume
 */
SMDS_MeshVolume* StdMeshers_Helper::AddVolume(const SMDS_MeshNode* n1,
                                              const SMDS_MeshNode* n2,
                                              const SMDS_MeshNode* n3,
                                              const SMDS_MeshNode* n4,
                                              const SMDS_MeshNode* n5,
                                              const SMDS_MeshNode* n6)
{
  SMESHDS_Mesh * meshDS = GetMesh()->GetMeshDS();
  if(!myCreateQuadratic) {
    return meshDS->AddVolume(n1, n2, n3, n4, n5, n6);
  }

  const SMDS_MeshNode* n12 = GetMediumNode(n1,n2,true);
  const SMDS_MeshNode* n23 = GetMediumNode(n2,n3,true);
  const SMDS_MeshNode* n31 = GetMediumNode(n3,n1,true);

  const SMDS_MeshNode* n45 = GetMediumNode(n4,n5,true);
  const SMDS_MeshNode* n56 = GetMediumNode(n5,n6,true);
  const SMDS_MeshNode* n64 = GetMediumNode(n6,n4,true);

  const SMDS_MeshNode* n14 = GetMediumNode(n1,n4,true);
  const SMDS_MeshNode* n25 = GetMediumNode(n2,n5,true);
  const SMDS_MeshNode* n36 = GetMediumNode(n3,n6,true);

  return meshDS->AddVolume(n1, n2, n3, n4, n5, n6,
                           n12, n23, n31, n45, n56, n64, n14, n25, n36);
}


//=======================================================================
//function : AddVolume
//purpose  : 
//=======================================================================
/*!
 * Special function for creation quadratic volume
 */
SMDS_MeshVolume* StdMeshers_Helper::AddVolume(const SMDS_MeshNode* n1,
                                              const SMDS_MeshNode* n2,
                                              const SMDS_MeshNode* n3,
                                              const SMDS_MeshNode* n4)
{
  SMESHDS_Mesh * meshDS = GetMesh()->GetMeshDS();
  if(!myCreateQuadratic) {
    return meshDS->AddVolume(n1, n2, n3, n4);
  }

  const SMDS_MeshNode* n12 = GetMediumNode(n1,n2,true);
  const SMDS_MeshNode* n23 = GetMediumNode(n2,n3,true);
  const SMDS_MeshNode* n31 = GetMediumNode(n3,n1,true);

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

  return meshDS->AddVolume(n1, n2, n3, n4, n12, n23, n31, n14, n24, n34);
}


//=======================================================================
//function : AddVolume
//purpose  : 
//=======================================================================
/*!
 * Special function for creation quadratic volume
 */
SMDS_MeshVolume* StdMeshers_Helper::AddVolume(const SMDS_MeshNode* n1,
                                              const SMDS_MeshNode* n2,
                                              const SMDS_MeshNode* n3,
                                              const SMDS_MeshNode* n4,
                                              const SMDS_MeshNode* n5,
                                              const SMDS_MeshNode* n6,
                                              const SMDS_MeshNode* n7,
                                              const SMDS_MeshNode* n8)
{
  SMESHDS_Mesh * meshDS = GetMesh()->GetMeshDS();
  if(!myCreateQuadratic) {
    return meshDS->AddVolume(n1, n2, n3, n4, n5, n6, n7, n8);
  }

  const SMDS_MeshNode* n12 = GetMediumNode(n1,n2,true);
  const SMDS_MeshNode* n23 = GetMediumNode(n2,n3,true);
  const SMDS_MeshNode* n34 = GetMediumNode(n3,n4,true);
  const SMDS_MeshNode* n41 = GetMediumNode(n4,n1,true);

  const SMDS_MeshNode* n56 = GetMediumNode(n5,n6,true);
  const SMDS_MeshNode* n67 = GetMediumNode(n6,n7,true);
  const SMDS_MeshNode* n78 = GetMediumNode(n7,n8,true);
  const SMDS_MeshNode* n85 = GetMediumNode(n8,n5,true);

  const SMDS_MeshNode* n15 = GetMediumNode(n1,n5,true);
  const SMDS_MeshNode* n26 = GetMediumNode(n2,n6,true);
  const SMDS_MeshNode* n37 = GetMediumNode(n3,n7,true);
  const SMDS_MeshNode* n48 = GetMediumNode(n4,n8,true);

  return meshDS->AddVolume(n1, n2, n3, n4, n5, n6, n7, n8,
                           n12, n23, n34, n41, n56, n67,
                           n78, n85, n15, n26, n37, n48);
}