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