src/StdMeshers/StdMeshers_Helper.cxx

   1 // File:      StdMeshers_Helper.cxx
   2 // Created:   15.02.06 15:22:41
   3 // Author:    Sergey KUUL
   4 // Copyright: Open CASCADE 2006
   5
   6
   7 #include "StdMeshers_Helper.hxx"
   8
   9 #include "TopTools_MapOfShape.hxx"
  10 #include "BRepTools.hxx"
  11 #include "BRepTools_WireExplorer.hxx"
  12 #include "SMDS_FacePosition.hxx"
  13 #include "SMDS_EdgePosition.hxx"
  14 #include "SMESH_subMesh.hxx"
  15 #include "Geom2d_Curve.hxx"
  16 #include "Geom_Curve.hxx"
  17 #include "Geom_Surface.hxx"
  18 #include "BRep_Tool.hxx"
  19 #include "gp_Pnt2d.hxx"
  20
  21
  22 //=======================================================================
  23 //function : CheckShape
  24 //purpose  :
  25 //=======================================================================
  26
  27 bool StdMeshers_Helper::IsQuadraticSubMesh(const TopoDS_Shape& aSh,
  28                                            const bool QuadMode)
  29 {
  30   SMESHDS_Mesh* meshDS = GetMesh()->GetMeshDS();
  31   myShapeID = meshDS->ShapeToIndex(aSh);
  32   myCreateQuadratic = false;
  33   if(QuadMode) {
  34     // we can create quadratic elements only if each elements
  35     // created on given shape is quadratic
  36     // also we have to fill myNLinkNodeMap
  37     myCreateQuadratic = true;
  38     if(aSh.ShapeType()!=TopAbs_FACE) {
  39       for (TopExp_Explorer exp(aSh, TopAbs_FACE); exp.More() && myCreateQuadratic; exp.Next()) {
  40         const TopoDS_Face& F = TopoDS::Face(exp.Current());
  41         SMDS_ElemIteratorPtr itf = GetMesh()->GetSubMesh(F)->GetSubMeshDS()->GetElements();
  42         while(itf->more()) {
  43           const SMDS_MeshElement* f = itf->next();
  44           if( f->GetType()==SMDSAbs_Face && !f->IsQuadratic() ) {
  45             myCreateQuadratic = false;
  46             break;
  47           }
  48           SMDS_ElemIteratorPtr itn = f->nodesIterator();
  49           if(f->NbNodes()==6) {
  50             const SMDS_MeshNode* Ns[6];
  51             int i = 0;
  52             while(itn->more()) {
  53               Ns[i++] = static_cast<const SMDS_MeshNode*>( itn->next() );
  54             }
  55             AddNLinkNode(Ns[0],Ns[1],Ns[3]);
  56             AddNLinkNode(Ns[1],Ns[2],Ns[4]);
  57             AddNLinkNode(Ns[2],Ns[0],Ns[5]);
  58           }
  59           else if(f->NbNodes()==8) {
  60             const SMDS_MeshNode* Ns[8];
  61             int i = 0;
  62             while(itn->more()) {
  63               Ns[i++] = static_cast<const SMDS_MeshNode*>( itn->next() );
  64             }
  65             AddNLinkNode(Ns[0],Ns[1],Ns[4]);
  66             AddNLinkNode(Ns[1],Ns[2],Ns[5]);
  67             AddNLinkNode(Ns[2],Ns[3],Ns[6]);
  68             AddNLinkNode(Ns[3],Ns[0],Ns[7]);
  69           }
  70         }
  71       }
  72     }
  73     else {
  74       TopTools_MapOfShape aMap;
  75       // check edges
  76       const TopoDS_Face& F = TopoDS::Face(aSh);
  77       const TopoDS_Wire& W = BRepTools::OuterWire(F);
  78       BRepTools_WireExplorer wexp (W, F);
  79       for (wexp.Init(W, F); wexp.More() && myCreateQuadratic; wexp.Next()) {
  80         const TopoDS_Edge& E = wexp.Current();
  81         if(aMap.Contains(E))
  82           continue;
  83         aMap.Add(E);
  84         SMDS_ElemIteratorPtr it = GetMesh()->GetSubMesh(E)->GetSubMeshDS()->GetElements();
  85         while(it->more()) {
  86           const SMDS_MeshElement* e = it->next();
  87           if( e->GetType()==SMDSAbs_Edge && !e->IsQuadratic() ) {
  88             myCreateQuadratic = false;
  89             break;
  90           }
  91           // fill NLinkNodeMap
  92           SMDS_ElemIteratorPtr nodeIt = e->nodesIterator();
  93           const SMDS_MeshNode* n1 = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
  94           const SMDS_MeshNode* n2 = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
  95           const SMDS_MeshNode* n3 = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
  96           NLink link(( n1 < n2 ? n1 : n2 ), ( n1 < n2 ? n2 : n1 ));
  97           myNLinkNodeMap.insert(NLinkNodeMap::value_type(link,n3));
  98           myNLinkNodeMap[link] = n3;
  99         }
 100       }
 101     }
 102   }
 103
 104   if(!myCreateQuadratic) {
 105     myNLinkNodeMap.clear();
 106   }
 107
 108   return myCreateQuadratic;
 109 }
 110
 111
 112 //=======================================================================
 113 //function : IsMedium
 114 //purpose  :
 115 //=======================================================================
 116
 117 bool StdMeshers_Helper::IsMedium(const SMDS_MeshNode* n)
 118 {
 119   SMDS_ElemIteratorPtr it = n->GetInverseElementIterator();
 120   while (it->more()) {
 121     const SMDS_MeshElement* elem = it->next();
 122     return elem->IsMediumNode(n);
 123   }
 124   return false;
 125 }
 126
 127
 128 //=======================================================================
 129 //function : AddNLinkNode
 130 //purpose  :
 131 //=======================================================================
 132 /*!
 133  * Auxilary function for filling myNLinkNodeMap
 134  */
 135 void StdMeshers_Helper::AddNLinkNode(const SMDS_MeshNode* n1,
 136                                      const SMDS_MeshNode* n2,
 137                                      const SMDS_MeshNode* n12)
 138 {
 139   NLink link(( n1 < n2 ? n1 : n2 ), ( n1 < n2 ? n2 : n1 ));
 140   ItNLinkNode itLN = myNLinkNodeMap.find( link );
 141   if ( itLN == myNLinkNodeMap.end() ) {
 142     // add new record to map
 143     myNLinkNodeMap.insert(NLinkNodeMap::value_type(link,n12));
 144   }
 145 }
 146
 147
 148 //=======================================================================
 149 //function : GetNodeUV
 150 //purpose  :
 151 //=======================================================================
 152 /*!
 153  * Auxilary function for GetMediumNode()
 154  */
 155 gp_XY StdMeshers_Helper::GetNodeUV(const TopoDS_Face& F,
 156                                    const SMDS_MeshNode* n)
 157 {
 158   gp_Pnt2d p2d;
 159   const SMDS_PositionPtr Pos = n->GetPosition();
 160   if(Pos->GetTypeOfPosition()==SMDS_TOP_FACE) {
 161 //cout<<"face"<<endl;
 162     // node has position on face
 163     const SMDS_FacePosition* fpos =
 164       static_cast<const SMDS_FacePosition*>(n->GetPosition().get());
 165     p2d = gp_Pnt2d(fpos->GetUParameter(),fpos->GetVParameter());
 166 //cout<<"fpos->GetUParameter()="<<fpos->GetUParameter()<<" fpos->GetVParameter()="<<fpos->GetVParameter()<<endl;
 167   }
 168   else if(Pos->GetTypeOfPosition()==SMDS_TOP_EDGE) {
 169     // node has position on edge => it is needed to find
 170     // corresponding edge from face, get pcurve for this
 171     // edge and recieve value from this pcurve
 172     const SMDS_EdgePosition* epos =
 173       static_cast<const SMDS_EdgePosition*>(n->GetPosition().get());
 174     double param = epos->GetUParameter();
 175     SMESHDS_Mesh* meshDS = GetMesh()->GetMeshDS();
 176     int edgeID = Pos->GetShapeId();
 177     const TopoDS_Edge& E = TopoDS::Edge(meshDS->IndexToShape(edgeID));
 178     double f, l;
 179     Handle(Geom2d_Curve) C2d = BRep_Tool::CurveOnSurface(E, F, f, l);
 180     p2d = C2d->Value(param);
 181   }
 182   else { // vertex position
 183     SMESHDS_Mesh * meshDS = GetMesh()->GetMeshDS();
 184     const TopoDS_Wire& W = BRepTools::OuterWire(F);
 185     BRepTools_WireExplorer wexp (W, F);
 186     TopoDS_Edge E;
 187     for (wexp.Init(W, F); wexp.More(); wexp.Next()) {
 188       E = wexp.Current();
 189       TopoDS_Vertex V1 = TopExp::FirstVertex(E);
 190       TopoDS_Vertex V2 = TopExp::LastVertex(E);
 191       if( meshDS->ShapeToIndex(V1) != Pos->GetShapeId() &&
 192           meshDS->ShapeToIndex(V2) != Pos->GetShapeId() ) continue;
 193       double f2, l2;
 194       Handle(Geom2d_Curve) C2d = BRep_Tool::CurveOnSurface(E, F, f2, l2);
 195       double f, l;
 196       TopLoc_Location L;
 197       Handle(Geom_Curve) C = BRep_Tool::Curve(E, L, f, l);
 198       gp_Pnt P = C->Value(f);
 199       if(!L.IsIdentity()) P = P.Transformed(L.Transformation());
 200       double tol = BRep_Tool::Tolerance(E);
 201       double dx = n->X() - P.X();
 202       double dy = n->Y() - P.Y();
 203       double dz = n->Z() - P.Z();
 204       double dist = sqrt(dx*dx + dy*dy + dz*dz);
 205       if(dist<tol) {
 206         p2d = C2d->Value(f2);
 207       }
 208       else {
 209         p2d = C2d->Value(l2);
 210       }
 211     }
 212   }
 213   return p2d.XY();
 214 }
 215
 216
 217 //=======================================================================
 218 //function : GetMediumNode
 219 //purpose  :
 220 //=======================================================================
 221 /*!
 222  * Special function for search or creation medium node
 223  */
 224 const SMDS_MeshNode* StdMeshers_Helper::GetMediumNode(const SMDS_MeshNode* n1,
 225                                                       const SMDS_MeshNode* n2,
 226                                                       const bool force3d)
 227 {
 228 //cout<<"n1: "<<n1;
 229 //cout<<"n2: "<<n2;
 230   NLink link(( n1 < n2 ? n1 : n2 ), ( n1 < n2 ? n2 : n1 ));
 231   ItNLinkNode itLN = myNLinkNodeMap.find( link );
 232   if ( itLN != myNLinkNodeMap.end() ) {
 233     return (*itLN).second;
 234   }
 235   else {
 236     // create medium node
 237     SMDS_MeshNode* n12;
 238     SMESHDS_Mesh* meshDS = GetMesh()->GetMeshDS();
 239     if(!force3d) {
 240       // we try to create medium node using UV parameters of
 241       // nodes, else - medium between corresponding 3d points
 242       const SMDS_PositionPtr Pos1 = n1->GetPosition();
 243       const SMDS_PositionPtr Pos2 = n2->GetPosition();
 244       int faceID = -1;
 245       if( Pos1->GetTypeOfPosition()==SMDS_TOP_FACE ) {
 246         faceID = Pos1->GetShapeId();
 247       }
 248       else if( Pos2->GetTypeOfPosition()==SMDS_TOP_FACE ) {
 249         faceID = Pos2->GetShapeId();
 250       }
 251       if(faceID>-1) {
 252         TopoDS_Face F = TopoDS::Face(meshDS->IndexToShape(faceID));
 253         gp_XY p1 = GetNodeUV(F,n1);
 254         gp_XY p2 = GetNodeUV(F,n2);
 255         double u = (p1.X()+p2.X())/2.;
 256         double v = (p1.Y()+p2.Y())/2.;
 257         Handle(Geom_Surface) S = BRep_Tool::Surface(F);
 258         gp_Pnt P = S->Value(u, v);
 259         n12 = meshDS->AddNode(P.X(), P.Y(), P.Z());
 260         meshDS->SetNodeOnFace(n12, faceID, u, v);
 261         myNLinkNodeMap.insert(NLinkNodeMap::value_type(link,n12));
 262         return n12;
 263       }
 264     }
 265     // 3d variant
 266     double x = ( n1->X() + n2->X() )/2.;
 267     double y = ( n1->Y() + n2->Y() )/2.;
 268     double z = ( n1->Z() + n2->Z() )/2.;
 269     n12 = meshDS->AddNode(x,y,z);
 270     meshDS->SetNodeInVolume(n12, myShapeID);
 271     myNLinkNodeMap.insert(NLinkNodeMap::value_type(link,n12));
 272     return n12;
 273   }
 274 }
 275
 276
 277 //=======================================================================
 278 //function : AddFace
 279 //purpose  :
 280 //=======================================================================
 281 /*!
 282  * Special function for creation quadratic triangle
 283  */
 284 SMDS_MeshFace* StdMeshers_Helper::AddFace(const SMDS_MeshNode* n1,
 285                                           const SMDS_MeshNode* n2,
 286                                           const SMDS_MeshNode* n3)
 287 {
 288   SMESHDS_Mesh * meshDS = GetMesh()->GetMeshDS();
 289   if(!myCreateQuadratic) {
 290     return  meshDS->AddFace(n1, n2, n3);
 291   }
 292
 293   const SMDS_MeshNode* n12 = GetMediumNode(n1,n2,false);
 294   const SMDS_MeshNode* n23 = GetMediumNode(n2,n3,false);
 295   const SMDS_MeshNode* n31 = GetMediumNode(n3,n1,false);
 296
 297   return  meshDS->AddFace(n1, n2, n3, n12, n23, n31);
 298 }
 299
 300
 301 //=======================================================================
 302 //function : AddFace
 303 //purpose  :
 304 //=======================================================================
 305 /*!
 306  * Special function for creation quadratic quadrangle
 307  */
 308 SMDS_MeshFace* StdMeshers_Helper::AddFace(const SMDS_MeshNode* n1,
 309                                           const SMDS_MeshNode* n2,
 310                                           const SMDS_MeshNode* n3,
 311                                           const SMDS_MeshNode* n4)
 312 {
 313   SMESHDS_Mesh * meshDS = GetMesh()->GetMeshDS();
 314   if(!myCreateQuadratic) {
 315     return  meshDS->AddFace(n1, n2, n3, n4);
 316   }
 317
 318   const SMDS_MeshNode* n12 = GetMediumNode(n1,n2,false);
 319   const SMDS_MeshNode* n23 = GetMediumNode(n2,n3,false);
 320   const SMDS_MeshNode* n34 = GetMediumNode(n3,n4,false);
 321   const SMDS_MeshNode* n41 = GetMediumNode(n4,n1,false);
 322
 323   return  meshDS->AddFace(n1, n2, n3, n4, n12, n23, n34, n41);
 324 }
 325
 326
 327 //=======================================================================
 328 //function : AddVolume
 329 //purpose  :
 330 //=======================================================================
 331 /*!
 332  * Special function for creation quadratic volume
 333  */
 334 SMDS_MeshVolume* StdMeshers_Helper::AddVolume(const SMDS_MeshNode* n1,
 335                                               const SMDS_MeshNode* n2,
 336                                               const SMDS_MeshNode* n3,
 337                                               const SMDS_MeshNode* n4,
 338                                               const SMDS_MeshNode* n5,
 339                                               const SMDS_MeshNode* n6)
 340 {
 341   SMESHDS_Mesh * meshDS = GetMesh()->GetMeshDS();
 342   if(!myCreateQuadratic) {
 343     return meshDS->AddVolume(n1, n2, n3, n4, n5, n6);
 344   }
 345
 346   const SMDS_MeshNode* n12 = GetMediumNode(n1,n2,true);
 347   const SMDS_MeshNode* n23 = GetMediumNode(n2,n3,true);
 348   const SMDS_MeshNode* n31 = GetMediumNode(n3,n1,true);
 349
 350   const SMDS_MeshNode* n45 = GetMediumNode(n4,n5,true);
 351   const SMDS_MeshNode* n56 = GetMediumNode(n5,n6,true);
 352   const SMDS_MeshNode* n64 = GetMediumNode(n6,n4,true);
 353
 354   const SMDS_MeshNode* n14 = GetMediumNode(n1,n4,true);
 355   const SMDS_MeshNode* n25 = GetMediumNode(n2,n5,true);
 356   const SMDS_MeshNode* n36 = GetMediumNode(n3,n6,true);
 357
 358   return meshDS->AddVolume(n1, n2, n3, n4, n5, n6,
 359                            n12, n23, n31, n45, n56, n64, n14, n25, n36);
 360 }
 361
 362
 363 //=======================================================================
 364 //function : AddVolume
 365 //purpose  :
 366 //=======================================================================
 367 /*!
 368  * Special function for creation quadratic volume
 369  */
 370 SMDS_MeshVolume* StdMeshers_Helper::AddVolume(const SMDS_MeshNode* n1,
 371                                               const SMDS_MeshNode* n2,
 372                                               const SMDS_MeshNode* n3,
 373                                               const SMDS_MeshNode* n4)
 374 {
 375   SMESHDS_Mesh * meshDS = GetMesh()->GetMeshDS();
 376   if(!myCreateQuadratic) {
 377     return meshDS->AddVolume(n1, n2, n3, n4);
 378   }
 379
 380   const SMDS_MeshNode* n12 = GetMediumNode(n1,n2,true);
 381   const SMDS_MeshNode* n23 = GetMediumNode(n2,n3,true);
 382   const SMDS_MeshNode* n31 = GetMediumNode(n3,n1,true);
 383
 384   const SMDS_MeshNode* n14 = GetMediumNode(n1,n4,true);
 385   const SMDS_MeshNode* n24 = GetMediumNode(n2,n4,true);
 386   const SMDS_MeshNode* n34 = GetMediumNode(n3,n4,true);
 387
 388   return meshDS->AddVolume(n1, n2, n3, n4, n12, n23, n31, n14, n24, n34);
 389 }
 390
 391
 392 //=======================================================================
 393 //function : AddVolume
 394 //purpose  :
 395 //=======================================================================
 396 /*!
 397  * Special function for creation quadratic volume
 398  */
 399 SMDS_MeshVolume* StdMeshers_Helper::AddVolume(const SMDS_MeshNode* n1,
 400                                               const SMDS_MeshNode* n2,
 401                                               const SMDS_MeshNode* n3,
 402                                               const SMDS_MeshNode* n4,
 403                                               const SMDS_MeshNode* n5,
 404                                               const SMDS_MeshNode* n6,
 405                                               const SMDS_MeshNode* n7,
 406                                               const SMDS_MeshNode* n8)
 407 {
 408   SMESHDS_Mesh * meshDS = GetMesh()->GetMeshDS();
 409   if(!myCreateQuadratic) {
 410     return meshDS->AddVolume(n1, n2, n3, n4, n5, n6, n7, n8);
 411   }
 412
 413   const SMDS_MeshNode* n12 = GetMediumNode(n1,n2,true);
 414   const SMDS_MeshNode* n23 = GetMediumNode(n2,n3,true);
 415   const SMDS_MeshNode* n34 = GetMediumNode(n3,n4,true);
 416   const SMDS_MeshNode* n41 = GetMediumNode(n4,n1,true);
 417
 418   const SMDS_MeshNode* n56 = GetMediumNode(n5,n6,true);
 419   const SMDS_MeshNode* n67 = GetMediumNode(n6,n7,true);
 420   const SMDS_MeshNode* n78 = GetMediumNode(n7,n8,true);
 421   const SMDS_MeshNode* n85 = GetMediumNode(n8,n5,true);
 422
 423   const SMDS_MeshNode* n15 = GetMediumNode(n1,n5,true);
 424   const SMDS_MeshNode* n26 = GetMediumNode(n2,n6,true);
 425   const SMDS_MeshNode* n37 = GetMediumNode(n3,n7,true);
 426   const SMDS_MeshNode* n48 = GetMediumNode(n4,n8,true);
 427
 428   return meshDS->AddVolume(n1, n2, n3, n4, n5, n6, n7, n8,
 429                            n12, n23, n34, n41, n56, n67,
 430                            n78, n85, n15, n26, n37, n48);
 431 }
 432
 433