src/SMESH/SMESH_MesherHelper.cxx

   1 // Copyright (C) 2005  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
   2 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
   3 //
   4 // This library is free software; you can redistribute it and/or
   5 // modify it under the terms of the GNU Lesser General Public
   6 // License as published by the Free Software Foundation; either
   7 // version 2.1 of the License.
   8 //
   9 // This library is distributed in the hope that it will be useful
  10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  12 // Lesser General Public License for more details.
  13 //
  14 // You should have received a copy of the GNU Lesser General Public
  15 // License along with this library; if not, write to the Free Software
  16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  17 //
  18 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  19 //
  20 // File:      SMESH_MesherHelper.cxx
  21 // Created:   15.02.06 15:22:41
  22 // Author:    Sergey KUUL
  23 // Copyright: Open CASCADE 2006
  24
  25
  26 #include "SMESH_MesherHelper.hxx"
  27
  28 #include "SMDS_FacePosition.hxx"
  29 #include "SMDS_EdgePosition.hxx"
  30 #include "SMESH_MeshEditor.hxx"
  31
  32 #include <BRepAdaptor_Surface.hxx>
  33 #include <BRepTools.hxx>
  34 #include <BRep_Tool.hxx>
  35 #include <Geom2d_Curve.hxx>
  36 #include <Geom_Curve.hxx>
  37 #include <Geom_Surface.hxx>
  38 #include <TopExp_Explorer.hxx>
  39 #include <TopTools_MapOfShape.hxx>
  40 #include <gp_Pnt2d.hxx>
  41 #include <ShapeAnalysis.hxx>
  42
  43 #include <utilities.h>
  44
  45 #define RETURN_BAD_RESULT(msg) { MESSAGE(msg); return false; }
  46
  47 //================================================================================
  48 /*!
  49  * \brief Constructor
  50  */
  51 //================================================================================
  52
  53 SMESH_MesherHelper::SMESH_MesherHelper(SMESH_Mesh& theMesh)
  54   : myMesh(&theMesh), myShapeID(-1), myCreateQuadratic(false)
  55 {
  56   mySetElemOnShape = ( ! myMesh->HasShapeToMesh() );
  57 }
  58
  59 //=======================================================================
  60 //function : CheckShape
  61 //purpose  :
  62 //=======================================================================
  63
  64 bool SMESH_MesherHelper::IsQuadraticSubMesh(const TopoDS_Shape& aSh)
  65 {
  66   SMESHDS_Mesh* meshDS = GetMeshDS();
  67   // we can create quadratic elements only if all elements
  68   // created on subshapes of given shape are quadratic
  69   // also we have to fill myNLinkNodeMap
  70   myCreateQuadratic = true;
  71   mySeamShapeIds.clear();
  72   TopAbs_ShapeEnum subType( aSh.ShapeType()==TopAbs_FACE ? TopAbs_EDGE : TopAbs_FACE );
  73   SMDSAbs_ElementType elemType( subType==TopAbs_FACE ? SMDSAbs_Face : SMDSAbs_Edge );
  74
  75   TopExp_Explorer exp( aSh, subType );
  76   for (; exp.More() && myCreateQuadratic; exp.Next()) {
  77     if ( SMESHDS_SubMesh * subMesh = meshDS->MeshElements( exp.Current() )) {
  78       if ( SMDS_ElemIteratorPtr it = subMesh->GetElements() ) {
  79         while(it->more()) {
  80           const SMDS_MeshElement* e = it->next();
  81           if ( e->GetType() != elemType || !e->IsQuadratic() ) {
  82             myCreateQuadratic = false;
  83             break;
  84           }
  85           else {
  86             // fill NLinkNodeMap
  87             switch ( e->NbNodes() ) {
  88             case 3:
  89               AddNLinkNode(e->GetNode(0),e->GetNode(1),e->GetNode(2)); break;
  90             case 6:
  91               AddNLinkNode(e->GetNode(0),e->GetNode(1),e->GetNode(3));
  92               AddNLinkNode(e->GetNode(1),e->GetNode(2),e->GetNode(4));
  93               AddNLinkNode(e->GetNode(2),e->GetNode(0),e->GetNode(5)); break;
  94             case 8:
  95               AddNLinkNode(e->GetNode(0),e->GetNode(1),e->GetNode(4));
  96               AddNLinkNode(e->GetNode(1),e->GetNode(2),e->GetNode(5));
  97               AddNLinkNode(e->GetNode(2),e->GetNode(3),e->GetNode(6));
  98               AddNLinkNode(e->GetNode(3),e->GetNode(0),e->GetNode(7));
  99               break;
 100             default:
 101               myCreateQuadratic = false;
 102               break;
 103             }
 104           }
 105         }
 106       }
 107     }
 108   }
 109
 110   if(!myCreateQuadratic) {
 111     myNLinkNodeMap.clear();
 112   }
 113   SetSubShape( aSh );
 114
 115   return myCreateQuadratic;
 116 }
 117
 118 //================================================================================
 119 /*!
 120  * \brief Set geomerty to make elements on
 121   * \param aSh - geomertic shape
 122  */
 123 //================================================================================
 124
 125 void SMESH_MesherHelper::SetSubShape(const int aShID)
 126 {
 127   if ( aShID == myShapeID )
 128     return;
 129   if ( aShID > 1 )
 130     SetSubShape( GetMeshDS()->IndexToShape( aShID ));
 131   else
 132     SetSubShape( TopoDS_Shape() );
 133 }
 134
 135 //================================================================================
 136 /*!
 137  * \brief Set geomerty to make elements on
 138   * \param aSh - geomertic shape
 139  */
 140 //================================================================================
 141
 142 void SMESH_MesherHelper::SetSubShape(const TopoDS_Shape& aSh)
 143 {
 144   if ( myShape.IsSame( aSh ))
 145     return;
 146
 147   myShape = aSh;
 148   mySeamShapeIds.clear();
 149
 150   if ( myShape.IsNull() ) {
 151     myShapeID  = -1;
 152     return;
 153   }
 154   SMESHDS_Mesh* meshDS = GetMeshDS();
 155   myShapeID = meshDS->ShapeToIndex(aSh);
 156
 157   // treatment of periodic faces
 158   for ( TopExp_Explorer eF( aSh, TopAbs_FACE ); eF.More(); eF.Next() )
 159   {
 160     const TopoDS_Face& face = TopoDS::Face( eF.Current() );
 161     BRepAdaptor_Surface surface( face );
 162     if ( surface.IsUPeriodic() || surface.IsVPeriodic() )
 163     {
 164       // look for a seam edge
 165       for ( TopExp_Explorer exp( face, TopAbs_EDGE ); exp.More(); exp.Next()) {
 166         const TopoDS_Edge& edge = TopoDS::Edge( exp.Current() );
 167         if ( BRep_Tool::IsClosed( edge, face )) {
 168           // initialize myPar1, myPar2 and myParIndex
 169           if ( mySeamShapeIds.empty() ) {
 170             gp_Pnt2d uv1, uv2;
 171             BRep_Tool::UVPoints( edge, face, uv1, uv2 );
 172             if ( Abs( uv1.Coord(1) - uv2.Coord(1) ) < Abs( uv1.Coord(2) - uv2.Coord(2) ))
 173             {
 174               myParIndex = 1; // U periodic
 175               myPar1 = surface.FirstUParameter();
 176               myPar2 = surface.LastUParameter();
 177             }
 178             else {
 179               myParIndex = 2;  // V periodic
 180               myPar1 = surface.FirstVParameter();
 181               myPar2 = surface.LastVParameter();
 182             }
 183           }
 184           // store shapes indices
 185           mySeamShapeIds.insert( meshDS->ShapeToIndex( exp.Current() ));
 186           for ( TopExp_Explorer v( exp.Current(), TopAbs_VERTEX ); v.More(); v.Next() )
 187             mySeamShapeIds.insert( meshDS->ShapeToIndex( v.Current() ));
 188         }
 189       }
 190     }
 191   }
 192 }
 193
 194 //================================================================================
 195   /*!
 196    * \brief Check if inFaceNode argument is necessary for call GetNodeUV(F,..)
 197     * \param F - the face
 198     * \retval bool - return true if the face is periodic
 199    */
 200 //================================================================================
 201
 202 bool SMESH_MesherHelper::GetNodeUVneedInFaceNode(const TopoDS_Face& F) const
 203 {
 204   if ( F.IsNull() ) return !mySeamShapeIds.empty();
 205
 206   if ( !F.IsNull() && !myShape.IsNull() && myShape.IsSame( F ))
 207     return !mySeamShapeIds.empty();
 208
 209   Handle(Geom_Surface) aSurface = BRep_Tool::Surface( F );
 210   if ( !aSurface.IsNull() )
 211     return ( aSurface->IsUPeriodic() || aSurface->IsVPeriodic() );
 212
 213   return false;
 214 }
 215
 216 //=======================================================================
 217 //function : IsMedium
 218 //purpose  :
 219 //=======================================================================
 220
 221 bool SMESH_MesherHelper::IsMedium(const SMDS_MeshNode*      node,
 222                                  const SMDSAbs_ElementType typeToCheck)
 223 {
 224   return SMESH_MeshEditor::IsMedium( node, typeToCheck );
 225 }
 226
 227 //=======================================================================
 228 //function : AddNLinkNode
 229 //purpose  :
 230 //=======================================================================
 231 /*!
 232  * Auxilary function for filling myNLinkNodeMap
 233  */
 234 void SMESH_MesherHelper::AddNLinkNode(const SMDS_MeshNode* n1,
 235                                      const SMDS_MeshNode* n2,
 236                                      const SMDS_MeshNode* n12)
 237 {
 238   NLink link( n1, n2 );
 239   if ( n1 > n2 ) link = NLink( n2, n1 );
 240   // add new record to map
 241   myNLinkNodeMap.insert( make_pair(link,n12));
 242 }
 243
 244 //=======================================================================
 245 /*!
 246  * \brief Select UV on either of 2 pcurves of a seam edge, closest to the given UV
 247  * \param uv1 - UV on the seam
 248  * \param uv2 - UV within a face
 249  * \retval gp_Pnt2d - selected UV
 250  */
 251 //=======================================================================
 252
 253 gp_Pnt2d SMESH_MesherHelper::GetUVOnSeam( const gp_Pnt2d& uv1, const gp_Pnt2d& uv2 ) const
 254 {
 255   double p1 = uv1.Coord( myParIndex );
 256   double p2 = uv2.Coord( myParIndex );
 257   double p3 = ( Abs( p1 - myPar1 ) < Abs( p1 - myPar2 )) ? myPar2 : myPar1;
 258   if ( Abs( p2 - p1 ) > Abs( p2 - p3 ))
 259     p1 = p3;
 260   gp_Pnt2d result = uv1;
 261   result.SetCoord( myParIndex, p1 );
 262   return result;
 263 }
 264
 265 //=======================================================================
 266 /*!
 267  * \brief Return node UV on face
 268  * \param F - the face
 269  * \param n - the node
 270  * \param n2 - a node of element being created located inside a face
 271  * \retval gp_XY - resulting UV
 272  *
 273  * Auxilary function called form GetMediumNode()
 274  */
 275 //=======================================================================
 276
 277 gp_XY SMESH_MesherHelper::GetNodeUV(const TopoDS_Face&   F,
 278                                     const SMDS_MeshNode* n,
 279                                     const SMDS_MeshNode* n2) const
 280 {
 281   gp_Pnt2d uv( 1e100, 1e100 );
 282   const SMDS_PositionPtr Pos = n->GetPosition();
 283   if(Pos->GetTypeOfPosition()==SMDS_TOP_FACE)
 284   {
 285     // node has position on face
 286     const SMDS_FacePosition* fpos =
 287       static_cast<const SMDS_FacePosition*>(n->GetPosition().get());
 288     uv = gp_Pnt2d(fpos->GetUParameter(),fpos->GetVParameter());
 289   }
 290   else if(Pos->GetTypeOfPosition()==SMDS_TOP_EDGE)
 291   {
 292     // node has position on edge => it is needed to find
 293     // corresponding edge from face, get pcurve for this
 294     // edge and recieve value from this pcurve
 295     const SMDS_EdgePosition* epos =
 296       static_cast<const SMDS_EdgePosition*>(n->GetPosition().get());
 297     SMESHDS_Mesh* meshDS = GetMeshDS();
 298     int edgeID = Pos->GetShapeId();
 299     TopoDS_Edge E = TopoDS::Edge(meshDS->IndexToShape(edgeID));
 300     double f, l;
 301     TopLoc_Location loc;
 302     Handle(Geom2d_Curve) C2d = BRep_Tool::CurveOnSurface(E, F, f, l);
 303     uv = C2d->Value( epos->GetUParameter() );
 304     // for a node on a seam edge select one of UVs on 2 pcurves
 305     if ( n2 && mySeamShapeIds.find( edgeID ) != mySeamShapeIds.end() )
 306       uv = GetUVOnSeam( uv, GetNodeUV( F, n2, 0 ));
 307   }
 308   else if(Pos->GetTypeOfPosition()==SMDS_TOP_VERTEX)
 309   {
 310     int vertexID = n->GetPosition()->GetShapeId();
 311     const TopoDS_Vertex& V = TopoDS::Vertex(GetMeshDS()->IndexToShape(vertexID));
 312     uv = BRep_Tool::Parameters( V, F );
 313     if ( n2 && mySeamShapeIds.find( vertexID ) != mySeamShapeIds.end() )
 314       uv = GetUVOnSeam( uv, GetNodeUV( F, n2, 0 ));
 315   }
 316   return uv.XY();
 317 }
 318
 319 //=======================================================================
 320 /*!
 321  * \brief Return node U on edge
 322  * \param E - the Edge
 323  * \param n - the node
 324  * \retval double - resulting U
 325  *
 326  * Auxilary function called form GetMediumNode()
 327  */
 328 //=======================================================================
 329
 330 double SMESH_MesherHelper::GetNodeU(const TopoDS_Edge&   E,
 331                                     const SMDS_MeshNode* n)
 332 {
 333   double param = 0;
 334   const SMDS_PositionPtr Pos = n->GetPosition();
 335   if(Pos->GetTypeOfPosition()==SMDS_TOP_EDGE) {
 336     const SMDS_EdgePosition* epos =
 337       static_cast<const SMDS_EdgePosition*>(n->GetPosition().get());
 338     param =  epos->GetUParameter();
 339   }
 340   else if(Pos->GetTypeOfPosition()==SMDS_TOP_VERTEX) {
 341     SMESHDS_Mesh * meshDS = GetMeshDS();
 342     int vertexID = n->GetPosition()->GetShapeId();
 343     const TopoDS_Vertex& V = TopoDS::Vertex(meshDS->IndexToShape(vertexID));
 344     param =  BRep_Tool::Parameter( V, E );
 345   }
 346   return param;
 347 }
 348
 349 //=======================================================================
 350 //function : GetMediumNode
 351 //purpose  :
 352 //=======================================================================
 353 /*!
 354  * Special function for search or creation medium node
 355  */
 356 const SMDS_MeshNode* SMESH_MesherHelper::GetMediumNode(const SMDS_MeshNode* n1,
 357                                                        const SMDS_MeshNode* n2,
 358                                                        bool force3d)
 359 {
 360   TopAbs_ShapeEnum shapeType = myShape.IsNull() ? TopAbs_SHAPE : myShape.ShapeType();
 361
 362   NLink link(( n1 < n2 ? n1 : n2 ), ( n1 < n2 ? n2 : n1 ));
 363   ItNLinkNode itLN = myNLinkNodeMap.find( link );
 364   if ( itLN != myNLinkNodeMap.end() ) {
 365     return (*itLN).second;
 366   }
 367   else {
 368     // create medium node
 369     SMDS_MeshNode* n12;
 370     SMESHDS_Mesh* meshDS = GetMeshDS();
 371     int faceID = -1, edgeID = -1;
 372     const SMDS_PositionPtr Pos1 = n1->GetPosition();
 373     const SMDS_PositionPtr Pos2 = n2->GetPosition();
 374
 375     if( myShape.IsNull() )
 376     {
 377       if( Pos1->GetTypeOfPosition()==SMDS_TOP_FACE ) {
 378         faceID = Pos1->GetShapeId();
 379       }
 380       else if( Pos2->GetTypeOfPosition()==SMDS_TOP_FACE ) {
 381         faceID = Pos2->GetShapeId();
 382       }
 383
 384       if( Pos1->GetTypeOfPosition()==SMDS_TOP_EDGE ) {
 385         edgeID = Pos1->GetShapeId();
 386       }
 387       if( Pos2->GetTypeOfPosition()==SMDS_TOP_EDGE ) {
 388         edgeID = Pos2->GetShapeId();
 389       }
 390     }
 391
 392     if(!force3d) {
 393       // we try to create medium node using UV parameters of
 394       // nodes, else - medium between corresponding 3d points
 395       if(faceID>-1 || shapeType == TopAbs_FACE) {
 396         // obtaining a face and 2d points for nodes
 397         TopoDS_Face F;
 398         if( myShape.IsNull() )
 399           F = TopoDS::Face(meshDS->IndexToShape(faceID));
 400         else {
 401           F = TopoDS::Face(myShape);
 402           faceID = myShapeID;
 403         }
 404
 405         gp_XY p1 = GetNodeUV(F,n1,n2);
 406         gp_XY p2 = GetNodeUV(F,n2,n1);
 407
 408         //checking if surface is periodic
 409         Handle(Geom_Surface) S = BRep_Tool::Surface(F);
 410         Standard_Real UF,UL,VF,VL;
 411         S->Bounds(UF,UL,VF,VL);
 412
 413         Standard_Real u,v;
 414         Standard_Boolean isUPeriodic = S->IsUPeriodic();
 415         if(isUPeriodic) {
 416           Standard_Real UPeriod = S->UPeriod();
 417           Standard_Real p2x = p2.X()+ShapeAnalysis::AdjustByPeriod(p2.X(),p1.X(),UPeriod);
 418           Standard_Real pmid = (p1.X()+p2x)/2.;
 419           u = pmid+ShapeAnalysis::AdjustToPeriod(pmid,UF,UL);
 420         }
 421         else
 422           u= (p1.X()+p2.X())/2.;
 423
 424         Standard_Boolean isVPeriodic = S->IsVPeriodic();
 425         if(isVPeriodic) {
 426           Standard_Real VPeriod = S->VPeriod();
 427           Standard_Real p2y = p2.Y()+ShapeAnalysis::AdjustByPeriod(p2.Y(),p1.Y(),VPeriod);
 428           Standard_Real pmid = (p1.Y()+p2y)/2.;
 429           v = pmid+ShapeAnalysis::AdjustToPeriod(pmid,VF,VL);
 430         }
 431         else
 432           v = (p1.Y()+p2.Y())/2.;
 433
 434         gp_Pnt P = S->Value(u, v);
 435         n12 = meshDS->AddNode(P.X(), P.Y(), P.Z());
 436         meshDS->SetNodeOnFace(n12, faceID, u, v);
 437         myNLinkNodeMap.insert(NLinkNodeMap::value_type(link,n12));
 438         return n12;
 439       }
 440       if (edgeID>-1 || shapeType == TopAbs_EDGE) {
 441
 442         TopoDS_Edge E;
 443         if( myShape.IsNull() )
 444           E = TopoDS::Edge(meshDS->IndexToShape(edgeID));
 445         else {
 446           E = TopoDS::Edge(myShape);
 447           edgeID = myShapeID;
 448         }
 449
 450         double p1 = GetNodeU(E,n1);
 451         double p2 = GetNodeU(E,n2);
 452
 453         double f,l;
 454         Handle(Geom_Curve) C = BRep_Tool::Curve(E, f, l);
 455         if(!C.IsNull()) {
 456
 457           Standard_Boolean isPeriodic = C->IsPeriodic();
 458           double u;
 459           if(isPeriodic) {
 460             Standard_Real Period = C->Period();
 461             Standard_Real p = p2+ShapeAnalysis::AdjustByPeriod(p2,p1,Period);
 462             Standard_Real pmid = (p1+p)/2.;
 463             u = pmid+ShapeAnalysis::AdjustToPeriod(pmid,C->FirstParameter(),C->LastParameter());
 464           }
 465           else
 466             u = (p1+p2)/2.;
 467
 468           gp_Pnt P = C->Value( u );
 469           n12 = meshDS->AddNode(P.X(), P.Y(), P.Z());
 470           meshDS->SetNodeOnEdge(n12, edgeID, u);
 471           myNLinkNodeMap.insert(NLinkNodeMap::value_type(link,n12));
 472           return n12;
 473         }
 474       }
 475     }
 476     // 3d variant
 477     double x = ( n1->X() + n2->X() )/2.;
 478     double y = ( n1->Y() + n2->Y() )/2.;
 479     double z = ( n1->Z() + n2->Z() )/2.;
 480     n12 = meshDS->AddNode(x,y,z);
 481     if(edgeID>-1)
 482         meshDS->SetNodeOnEdge(n12, edgeID);
 483     else if(faceID>-1)
 484         meshDS->SetNodeOnFace(n12, faceID);
 485     else
 486       meshDS->SetNodeInVolume(n12, myShapeID);
 487     myNLinkNodeMap.insert(NLinkNodeMap::value_type(link,n12));
 488     return n12;
 489   }
 490 }
 491
 492 //=======================================================================
 493 /*!
 494  * Creates a node
 495  */
 496 //=======================================================================
 497
 498 SMDS_MeshNode* SMESH_MesherHelper::AddNode(double x, double y, double z, int ID)
 499 {
 500   SMESHDS_Mesh * meshDS = GetMeshDS();
 501   SMDS_MeshNode* node = 0;
 502   if ( ID )
 503     node = meshDS->AddNodeWithID( x, y, z, ID );
 504   else
 505     node = meshDS->AddNode( x, y, z );
 506   if ( mySetElemOnShape && myShapeID > 0 ) {
 507     switch ( myShape.ShapeType() ) {
 508     case TopAbs_SOLID:  meshDS->SetNodeInVolume( node, myShapeID); break;
 509     case TopAbs_SHELL:  meshDS->SetNodeInVolume( node, myShapeID); break;
 510     case TopAbs_FACE:   meshDS->SetNodeOnFace(   node, myShapeID); break;
 511     case TopAbs_EDGE:   meshDS->SetNodeOnEdge(   node, myShapeID); break;
 512     case TopAbs_VERTEX: meshDS->SetNodeOnVertex( node, myShapeID); break;
 513     default: ;
 514     }
 515   }
 516   return node;
 517 }
 518
 519 //=======================================================================
 520 /*!
 521  * Creates quadratic or linear edge
 522  */
 523 //=======================================================================
 524
 525 SMDS_MeshEdge* SMESH_MesherHelper::AddEdge(const SMDS_MeshNode* n1,
 526                                                 const SMDS_MeshNode* n2,
 527                                                 const int id,
 528                                                 const bool force3d)
 529 {
 530   SMESHDS_Mesh * meshDS = GetMeshDS();
 531
 532   SMDS_MeshEdge* edge = 0;
 533   if (myCreateQuadratic) {
 534     const SMDS_MeshNode* n12 = GetMediumNode(n1,n2,force3d);
 535     if(id)
 536       edge = meshDS->AddEdgeWithID(n1, n2, n12, id);
 537     else
 538       edge = meshDS->AddEdge(n1, n2, n12);
 539   }
 540   else {
 541     if(id)
 542       edge = meshDS->AddEdgeWithID(n1, n2, id);
 543     else
 544       edge = meshDS->AddEdge(n1, n2);
 545   }
 546
 547   if ( mySetElemOnShape && myShapeID > 0 )
 548     meshDS->SetMeshElementOnShape( edge, myShapeID );
 549
 550   return edge;
 551 }
 552
 553 //=======================================================================
 554 /*!
 555  * Creates quadratic or linear triangle
 556  */
 557 //=======================================================================
 558
 559 SMDS_MeshFace* SMESH_MesherHelper::AddFace(const SMDS_MeshNode* n1,
 560                                            const SMDS_MeshNode* n2,
 561                                            const SMDS_MeshNode* n3,
 562                                            const int id,
 563                                            const bool force3d)
 564 {
 565   SMESHDS_Mesh * meshDS = GetMeshDS();
 566   SMDS_MeshFace* elem = 0;
 567   if(!myCreateQuadratic) {
 568     if(id)
 569       elem = meshDS->AddFaceWithID(n1, n2, n3, id);
 570     else
 571       elem = meshDS->AddFace(n1, n2, n3);
 572   }
 573   else {
 574     const SMDS_MeshNode* n12 = GetMediumNode(n1,n2,force3d);
 575     const SMDS_MeshNode* n23 = GetMediumNode(n2,n3,force3d);
 576     const SMDS_MeshNode* n31 = GetMediumNode(n3,n1,force3d);
 577
 578     if(id)
 579       elem = meshDS->AddFaceWithID(n1, n2, n3, n12, n23, n31, id);
 580     else
 581       elem = meshDS->AddFace(n1, n2, n3, n12, n23, n31);
 582   }
 583   if ( mySetElemOnShape && myShapeID > 0 )
 584     meshDS->SetMeshElementOnShape( elem, myShapeID );
 585
 586   return elem;
 587 }
 588
 589 //=======================================================================
 590 /*!
 591  * Creates quadratic or linear quadrangle
 592  */
 593 //=======================================================================
 594
 595 SMDS_MeshFace* SMESH_MesherHelper::AddFace(const SMDS_MeshNode* n1,
 596                                            const SMDS_MeshNode* n2,
 597                                            const SMDS_MeshNode* n3,
 598                                            const SMDS_MeshNode* n4,
 599                                            const int id,
 600                                            const bool force3d)
 601 {
 602   SMESHDS_Mesh * meshDS = GetMeshDS();
 603   SMDS_MeshFace* elem = 0;
 604   if(!myCreateQuadratic) {
 605     if(id)
 606       elem = meshDS->AddFaceWithID(n1, n2, n3, n4, id);
 607     else
 608       elem = meshDS->AddFace(n1, n2, n3, n4);
 609   }
 610   else {
 611     const SMDS_MeshNode* n12 = GetMediumNode(n1,n2,force3d);
 612     const SMDS_MeshNode* n23 = GetMediumNode(n2,n3,force3d);
 613     const SMDS_MeshNode* n34 = GetMediumNode(n3,n4,force3d);
 614     const SMDS_MeshNode* n41 = GetMediumNode(n4,n1,force3d);
 615
 616     if(id)
 617       elem = meshDS->AddFaceWithID(n1, n2, n3, n4, n12, n23, n34, n41, id);
 618     else
 619       elem = meshDS->AddFace(n1, n2, n3, n4, n12, n23, n34, n41);
 620   }
 621   if ( mySetElemOnShape && myShapeID > 0 )
 622     meshDS->SetMeshElementOnShape( elem, myShapeID );
 623
 624   return elem;
 625 }
 626
 627 //=======================================================================
 628 /*!
 629  * Creates quadratic or linear volume
 630  */
 631 //=======================================================================
 632
 633 SMDS_MeshVolume* SMESH_MesherHelper::AddVolume(const SMDS_MeshNode* n1,
 634                                                const SMDS_MeshNode* n2,
 635                                                const SMDS_MeshNode* n3,
 636                                                const SMDS_MeshNode* n4,
 637                                                const SMDS_MeshNode* n5,
 638                                                const SMDS_MeshNode* n6,
 639                                                const int id,
 640                                                const bool force3d)
 641 {
 642   SMESHDS_Mesh * meshDS = GetMeshDS();
 643   SMDS_MeshVolume* elem = 0;
 644   if(!myCreateQuadratic) {
 645     if(id)
 646       elem = meshDS->AddVolumeWithID(n1, n2, n3, n4, n5, n6, id);
 647     else
 648       elem = meshDS->AddVolume(n1, n2, n3, n4, n5, n6);
 649   }
 650   else {
 651     const SMDS_MeshNode* n12 = GetMediumNode(n1,n2,force3d);
 652     const SMDS_MeshNode* n23 = GetMediumNode(n2,n3,force3d);
 653     const SMDS_MeshNode* n31 = GetMediumNode(n3,n1,force3d);
 654
 655     const SMDS_MeshNode* n45 = GetMediumNode(n4,n5,force3d);
 656     const SMDS_MeshNode* n56 = GetMediumNode(n5,n6,force3d);
 657     const SMDS_MeshNode* n64 = GetMediumNode(n6,n4,force3d);
 658
 659     const SMDS_MeshNode* n14 = GetMediumNode(n1,n4,force3d);
 660     const SMDS_MeshNode* n25 = GetMediumNode(n2,n5,force3d);
 661     const SMDS_MeshNode* n36 = GetMediumNode(n3,n6,force3d);
 662
 663     if(id)
 664       elem = meshDS->AddVolumeWithID(n1, n2, n3, n4, n5, n6,
 665                                      n12, n23, n31, n45, n56, n64, n14, n25, n36, id);
 666     else
 667       elem = meshDS->AddVolume(n1, n2, n3, n4, n5, n6,
 668                                n12, n23, n31, n45, n56, n64, n14, n25, n36);
 669   }
 670   if ( mySetElemOnShape && myShapeID > 0 )
 671     meshDS->SetMeshElementOnShape( elem, myShapeID );
 672
 673   return elem;
 674 }
 675
 676 //=======================================================================
 677 /*!
 678  * Creates quadratic or linear volume
 679  */
 680 //=======================================================================
 681
 682 SMDS_MeshVolume* SMESH_MesherHelper::AddVolume(const SMDS_MeshNode* n1,
 683                                                const SMDS_MeshNode* n2,
 684                                                const SMDS_MeshNode* n3,
 685                                                const SMDS_MeshNode* n4,
 686                                                const int id,
 687                                                const bool force3d)
 688 {
 689   SMESHDS_Mesh * meshDS = GetMeshDS();
 690   SMDS_MeshVolume* elem = 0;
 691   if(!myCreateQuadratic) {
 692     if(id)
 693       elem = meshDS->AddVolumeWithID(n1, n2, n3, n4, id);
 694     else
 695       elem = meshDS->AddVolume(n1, n2, n3, n4);
 696   }
 697   else {
 698     const SMDS_MeshNode* n12 = GetMediumNode(n1,n2,force3d);
 699     const SMDS_MeshNode* n23 = GetMediumNode(n2,n3,force3d);
 700     const SMDS_MeshNode* n31 = GetMediumNode(n3,n1,force3d);
 701
 702     const SMDS_MeshNode* n14 = GetMediumNode(n1,n4,force3d);
 703     const SMDS_MeshNode* n24 = GetMediumNode(n2,n4,force3d);
 704     const SMDS_MeshNode* n34 = GetMediumNode(n3,n4,force3d);
 705
 706     if(id)
 707       elem = meshDS->AddVolumeWithID(n1, n2, n3, n4, n12, n23, n31, n14, n24, n34, id);
 708     else
 709       elem = meshDS->AddVolume(n1, n2, n3, n4, n12, n23, n31, n14, n24, n34);
 710   }
 711   if ( mySetElemOnShape && myShapeID > 0 )
 712     meshDS->SetMeshElementOnShape( elem, myShapeID );
 713
 714   return elem;
 715 }
 716
 717 //=======================================================================
 718 /*!
 719  * Creates quadratic or linear pyramid
 720  */
 721 //=======================================================================
 722
 723 SMDS_MeshVolume* SMESH_MesherHelper::AddVolume(const SMDS_MeshNode* n1,
 724                                                const SMDS_MeshNode* n2,
 725                                                const SMDS_MeshNode* n3,
 726                                                const SMDS_MeshNode* n4,
 727                                                const SMDS_MeshNode* n5,
 728                                                const int id,
 729                                                const bool force3d)
 730 {
 731   SMDS_MeshVolume* elem = 0;
 732   if(!myCreateQuadratic) {
 733     if(id)
 734       elem = GetMeshDS()->AddVolumeWithID(n1, n2, n3, n4, n5, id);
 735     else
 736       elem = GetMeshDS()->AddVolume(n1, n2, n3, n4, n5);
 737   }
 738   else {
 739     const SMDS_MeshNode* n12 = GetMediumNode(n1,n2,force3d);
 740     const SMDS_MeshNode* n23 = GetMediumNode(n2,n3,force3d);
 741     const SMDS_MeshNode* n34 = GetMediumNode(n3,n4,force3d);
 742     const SMDS_MeshNode* n41 = GetMediumNode(n4,n1,force3d);
 743
 744     const SMDS_MeshNode* n15 = GetMediumNode(n1,n5,force3d);
 745     const SMDS_MeshNode* n25 = GetMediumNode(n2,n5,force3d);
 746     const SMDS_MeshNode* n35 = GetMediumNode(n3,n5,force3d);
 747     const SMDS_MeshNode* n45 = GetMediumNode(n4,n5,force3d);
 748
 749     if(id)
 750       elem = GetMeshDS()->AddVolumeWithID ( n1,  n2,  n3,  n4,  n5,
 751                                             n12, n23, n34, n41,
 752                                             n15, n25, n35, n45,
 753                                             id);
 754     else
 755       elem = GetMeshDS()->AddVolume( n1,  n2,  n3,  n4,  n5,
 756                                      n12, n23, n34, n41,
 757                                      n15, n25, n35, n45);
 758   }
 759   if ( mySetElemOnShape && myShapeID > 0 )
 760     GetMeshDS()->SetMeshElementOnShape( elem, myShapeID );
 761
 762   return elem;
 763 }
 764
 765 //=======================================================================
 766 /*!
 767  * Creates quadratic or linear hexahedron
 768  */
 769 //=======================================================================
 770
 771 SMDS_MeshVolume* SMESH_MesherHelper::AddVolume(const SMDS_MeshNode* n1,
 772                                                const SMDS_MeshNode* n2,
 773                                                const SMDS_MeshNode* n3,
 774                                                const SMDS_MeshNode* n4,
 775                                                const SMDS_MeshNode* n5,
 776                                                const SMDS_MeshNode* n6,
 777                                                const SMDS_MeshNode* n7,
 778                                                const SMDS_MeshNode* n8,
 779                                                const int id,
 780                                                const bool force3d)
 781 {
 782   SMESHDS_Mesh * meshDS = GetMeshDS();
 783   SMDS_MeshVolume* elem = 0;
 784   if(!myCreateQuadratic) {
 785     if(id)
 786       elem = meshDS->AddVolumeWithID(n1, n2, n3, n4, n5, n6, n7, n8, id);
 787     else
 788       elem = meshDS->AddVolume(n1, n2, n3, n4, n5, n6, n7, n8);
 789   }
 790   else {
 791     const SMDS_MeshNode* n12 = GetMediumNode(n1,n2,force3d);
 792     const SMDS_MeshNode* n23 = GetMediumNode(n2,n3,force3d);
 793     const SMDS_MeshNode* n34 = GetMediumNode(n3,n4,force3d);
 794     const SMDS_MeshNode* n41 = GetMediumNode(n4,n1,force3d);
 795
 796     const SMDS_MeshNode* n56 = GetMediumNode(n5,n6,force3d);
 797     const SMDS_MeshNode* n67 = GetMediumNode(n6,n7,force3d);
 798     const SMDS_MeshNode* n78 = GetMediumNode(n7,n8,force3d);
 799     const SMDS_MeshNode* n85 = GetMediumNode(n8,n5,force3d);
 800
 801     const SMDS_MeshNode* n15 = GetMediumNode(n1,n5,force3d);
 802     const SMDS_MeshNode* n26 = GetMediumNode(n2,n6,force3d);
 803     const SMDS_MeshNode* n37 = GetMediumNode(n3,n7,force3d);
 804     const SMDS_MeshNode* n48 = GetMediumNode(n4,n8,force3d);
 805
 806     if(id)
 807       elem = meshDS->AddVolumeWithID(n1, n2, n3, n4, n5, n6, n7, n8,
 808                                      n12, n23, n34, n41, n56, n67,
 809                                      n78, n85, n15, n26, n37, n48, id);
 810     else
 811       elem = meshDS->AddVolume(n1, n2, n3, n4, n5, n6, n7, n8,
 812                                n12, n23, n34, n41, n56, n67,
 813                                n78, n85, n15, n26, n37, n48);
 814   }
 815   if ( mySetElemOnShape && myShapeID > 0 )
 816     meshDS->SetMeshElementOnShape( elem, myShapeID );
 817
 818   return elem;
 819 }
 820
 821 //=======================================================================
 822   /*!
 823    * \brief Load nodes bound to face into a map of node columns
 824     * \param theParam2ColumnMap - map of node columns to fill
 825     * \param theFace - the face on which nodes are searched for
 826     * \param theBaseEdge - the edge nodes of which are columns' bases
 827     * \param theMesh - the mesh containing nodes
 828     * \retval bool - false if something is wrong
 829    *
 830    * The key of the map is a normalized parameter of each
 831    * base node on theBaseEdge.
 832    * This method works in supposition that nodes on the face
 833    * forms a rectangular grid and elements can be quardrangles or triangles
 834    */
 835 //=======================================================================
 836
 837 bool SMESH_MesherHelper::LoadNodeColumns(TParam2ColumnMap & theParam2ColumnMap,
 838                                          const TopoDS_Face& theFace,
 839                                          const TopoDS_Edge& theBaseEdge,
 840                                          SMESHDS_Mesh*      theMesh)
 841 {
 842   // get vertices of theBaseEdge
 843   TopoDS_Vertex vfb, vlb, vft; // first and last, bottom and top vertices
 844   TopoDS_Edge eFrw = TopoDS::Edge( theBaseEdge.Oriented( TopAbs_FORWARD ));
 845   TopExp::Vertices( eFrw, vfb, vlb );
 846
 847   // find the other edges of theFace and orientation of e1
 848   TopoDS_Edge e1, e2, eTop;
 849   bool rev1, CumOri = false;
 850   TopExp_Explorer exp( theFace, TopAbs_EDGE );
 851   int nbEdges = 0;
 852   for ( ; exp.More(); exp.Next() ) {
 853     if ( ++nbEdges > 4 ) {
 854       return false; // more than 4 edges in theFace
 855     }
 856     TopoDS_Edge e = TopoDS::Edge( exp.Current() );
 857     if ( theBaseEdge.IsSame( e ))
 858       continue;
 859     TopoDS_Vertex vCommon;
 860     if ( !TopExp::CommonVertex( theBaseEdge, e, vCommon ))
 861       eTop = e;
 862     else if ( vCommon.IsSame( vfb )) {
 863       e1 = e;
 864       vft = TopExp::LastVertex( e1, CumOri );
 865       rev1 = vfb.IsSame( vft );
 866       if ( rev1 )
 867         vft = TopExp::FirstVertex( e1, CumOri );
 868     }
 869     else
 870       e2 = e;
 871   }
 872   if ( nbEdges < 4 ) {
 873     return false; // less than 4 edges in theFace
 874   }
 875   if ( e2.IsNull() && vfb.IsSame( vlb ))
 876     e2 = e1;
 877
 878   // submeshes corresponding to shapes
 879   SMESHDS_SubMesh* smFace = theMesh->MeshElements( theFace );
 880   SMESHDS_SubMesh* smb = theMesh->MeshElements( theBaseEdge );
 881   SMESHDS_SubMesh* smt = theMesh->MeshElements( eTop );
 882   SMESHDS_SubMesh* sm1 = theMesh->MeshElements( e1 );
 883   SMESHDS_SubMesh* sm2 = theMesh->MeshElements( e2 );
 884   SMESHDS_SubMesh* smVfb = theMesh->MeshElements( vfb );
 885   SMESHDS_SubMesh* smVlb = theMesh->MeshElements( vlb );
 886   SMESHDS_SubMesh* smVft = theMesh->MeshElements( vft );
 887   if (!smFace || !smb || !smt || !sm1 || !sm2 || !smVfb || !smVlb || !smVft ) {
 888     RETURN_BAD_RESULT( "NULL submesh " <<smFace<<" "<<smb<<" "<<smt<<" "<<
 889                        sm1<<" "<<sm2<<" "<<smVfb<<" "<<smVlb<<" "<<smVft);
 890   }
 891   if ( smb->NbNodes() != smt->NbNodes() || sm1->NbNodes() != sm2->NbNodes() ) {
 892     RETURN_BAD_RESULT(" Diff nb of nodes on opposite edges" );
 893   }
 894   if (smVfb->NbNodes() != 1 || smVlb->NbNodes() != 1 || smVft->NbNodes() != 1) {
 895     RETURN_BAD_RESULT("Empty submesh of vertex");
 896   }
 897   // define whether mesh is quadratic
 898   bool isQuadraticMesh = false;
 899   SMDS_ElemIteratorPtr eIt = smFace->GetElements();
 900   if ( !eIt->more() ) {
 901     RETURN_BAD_RESULT("No elements on the face");
 902   }
 903   const SMDS_MeshElement* e = eIt->next();
 904   isQuadraticMesh = e->IsQuadratic();
 905
 906   if ( sm1->NbNodes() * smb->NbNodes() != smFace->NbNodes() ) {
 907     // check quadratic case
 908     if ( isQuadraticMesh ) {
 909       // what if there are quadrangles and triangles mixed?
 910 //       int n1 = sm1->NbNodes()/2;
 911 //       int n2 = smb->NbNodes()/2;
 912 //       int n3 = sm1->NbNodes() - n1;
 913 //       int n4 = smb->NbNodes() - n2;
 914 //       int nf = sm1->NbNodes()*smb->NbNodes() - n3*n4;
 915 //       if( nf != smFace->NbNodes() ) {
 916 //         MESSAGE( "Wrong nb face nodes: " <<
 917 //                 sm1->NbNodes()<<" "<<smb->NbNodes()<<" "<<smFace->NbNodes());
 918 //         return false;
 919 //       }
 920     }
 921     else {
 922       RETURN_BAD_RESULT( "Wrong nb face nodes: " <<
 923                          sm1->NbNodes()<<" "<<smb->NbNodes()<<" "<<smFace->NbNodes());
 924     }
 925   }
 926   // IJ size
 927   int vsize = sm1->NbNodes() + 2;
 928   int hsize = smb->NbNodes() + 2;
 929   if(isQuadraticMesh) {
 930     vsize = vsize - sm1->NbNodes()/2 -1;
 931     hsize = hsize - smb->NbNodes()/2 -1;
 932   }
 933
 934   // load nodes from theBaseEdge
 935
 936   set<const SMDS_MeshNode*> loadedNodes;
 937   const SMDS_MeshNode* nullNode = 0;
 938
 939   vector<const SMDS_MeshNode*> & nVecf = theParam2ColumnMap[ 0.];
 940   nVecf.resize( vsize, nullNode );
 941   loadedNodes.insert( nVecf[ 0 ] = smVfb->GetNodes()->next() );
 942
 943   vector<const SMDS_MeshNode*> & nVecl = theParam2ColumnMap[ 1.];
 944   nVecl.resize( vsize, nullNode );
 945   loadedNodes.insert( nVecl[ 0 ] = smVlb->GetNodes()->next() );
 946
 947   double f, l;
 948   BRep_Tool::Range( eFrw, f, l );
 949   double range = l - f;
 950   SMDS_NodeIteratorPtr nIt = smb->GetNodes();
 951   const SMDS_MeshNode* node;
 952   while ( nIt->more() ) {
 953     node = nIt->next();
 954     if(IsMedium(node))
 955       continue;
 956     const SMDS_EdgePosition* pos =
 957       dynamic_cast<const SMDS_EdgePosition*>( node->GetPosition().get() );
 958     if ( !pos ) {
 959       return false;
 960     }
 961     double u = ( pos->GetUParameter() - f ) / range;
 962     vector<const SMDS_MeshNode*> & nVec = theParam2ColumnMap[ u ];
 963     nVec.resize( vsize, nullNode );
 964     loadedNodes.insert( nVec[ 0 ] = node );
 965   }
 966   if ( theParam2ColumnMap.size() != hsize ) {
 967     RETURN_BAD_RESULT( "Wrong node positions on theBaseEdge" );
 968   }
 969
 970   // load nodes from e1
 971
 972   map< double, const SMDS_MeshNode*> sortedNodes; // sort by param on edge
 973   nIt = sm1->GetNodes();
 974   while ( nIt->more() ) {
 975     node = nIt->next();
 976     if(IsMedium(node))
 977       continue;
 978     const SMDS_EdgePosition* pos =
 979       dynamic_cast<const SMDS_EdgePosition*>( node->GetPosition().get() );
 980     if ( !pos ) {
 981       return false;
 982     }
 983     sortedNodes.insert( make_pair( pos->GetUParameter(), node ));
 984   }
 985   loadedNodes.insert( nVecf[ vsize - 1 ] = smVft->GetNodes()->next() );
 986   map< double, const SMDS_MeshNode*>::iterator u_n = sortedNodes.begin();
 987   int row  = rev1 ? vsize - 1 : 0;
 988   int dRow = rev1 ? -1 : +1;
 989   for ( ; u_n != sortedNodes.end(); u_n++ ) {
 990     row += dRow;
 991     loadedNodes.insert( nVecf[ row ] = u_n->second );
 992   }
 993
 994   // try to load the rest nodes
 995
 996   // get all faces from theFace
 997   TIDSortedElemSet allFaces, foundFaces;
 998   eIt = smFace->GetElements();
 999   while ( eIt->more() ) {
1000     const SMDS_MeshElement* e = eIt->next();
1001     if ( e->GetType() == SMDSAbs_Face )
1002       allFaces.insert( e );
1003   }
1004   // Starting from 2 neighbour nodes on theBaseEdge, look for a face
1005   // the nodes belong to, and between the nodes of the found face,
1006   // look for a not loaded node considering this node to be the next
1007   // in a column of the starting second node. Repeat, starting
1008   // from nodes next to the previous starting nodes in their columns,
1009   // and so on while a face can be found. Then go the the next pair
1010   // of nodes on theBaseEdge.
1011   TParam2ColumnMap::iterator par_nVec_1 = theParam2ColumnMap.begin();
1012   TParam2ColumnMap::iterator par_nVec_2 = par_nVec_1;
1013   // loop on columns
1014   int col = 0;
1015   for ( par_nVec_2++; par_nVec_2 != theParam2ColumnMap.end(); par_nVec_1++, par_nVec_2++ ) {
1016     col++;
1017     row = 0;
1018     const SMDS_MeshNode* n1 = par_nVec_1->second[ row ];
1019     const SMDS_MeshNode* n2 = par_nVec_2->second[ row ];
1020     const SMDS_MeshElement* face = 0;
1021     bool lastColOnClosedFace = ( nVecf[ row ] == n2 );
1022     do {
1023       // look for a face by 2 nodes
1024       face = SMESH_MeshEditor::FindFaceInSet( n1, n2, allFaces, foundFaces );
1025       if ( face ) {
1026         int nbFaceNodes = face->NbNodes();
1027         if ( face->IsQuadratic() )
1028           nbFaceNodes /= 2;
1029         if ( nbFaceNodes>4 ) {
1030           RETURN_BAD_RESULT(" Too many nodes in a face: " << nbFaceNodes );
1031         }
1032         // look for a not loaded node of the <face>
1033         bool found = false;
1034         const SMDS_MeshNode* n3 = 0; // a node defferent from n1 and n2
1035         for ( int i = 0; i < nbFaceNodes && !found; ++i ) {
1036           node = face->GetNode( i );
1037           found = loadedNodes.insert( node ).second;
1038           if ( !found && node != n1 && node != n2 )
1039             n3 = node;
1040         }
1041         if ( lastColOnClosedFace && row + 1 < vsize ) {
1042           node = nVecf[ row + 1 ];
1043           found = ( face->GetNodeIndex( node ) >= 0 );
1044         }
1045         if ( found ) {
1046           if ( ++row > vsize - 1 ) {
1047             RETURN_BAD_RESULT( "Too many nodes in column "<< col <<": "<< row+1);
1048           }
1049           par_nVec_2->second[ row ] = node;
1050           foundFaces.insert( face );
1051           n2 = node;
1052           if ( nbFaceNodes==4 ) {
1053             n1 = par_nVec_1->second[ row ];
1054           }
1055         }
1056         else if ( nbFaceNodes==3 && n3 == par_nVec_1->second[ row + 1 ] ) {
1057           n1 = n3;
1058         }
1059         else  {
1060           RETURN_BAD_RESULT( "Not quad mesh, column "<< col );
1061         }
1062       }
1063     }
1064     while ( face && n1 && n2 );
1065
1066     if ( row < vsize - 1 ) {
1067       MESSAGE( "Too few nodes in column "<< col <<": "<< row+1);
1068       MESSAGE( "Base node 1: "<< par_nVec_1->second[0]);
1069       MESSAGE( "Base node 2: "<< par_nVec_2->second[0]);
1070       if ( n1 ) { MESSAGE( "Current node 1: "<< n1); }
1071       else      { MESSAGE( "Current node 1: NULL");  }
1072       if ( n2 ) { MESSAGE( "Current node 2: "<< n2); }
1073       else      { MESSAGE( "Current node 2: NULL");  }
1074       MESSAGE( "first base node: "<< theParam2ColumnMap.begin()->second[0]);
1075       MESSAGE( "last base node: "<< theParam2ColumnMap.rbegin()->second[0]);
1076       return false;
1077     }
1078   } // loop on columns
1079
1080   return true;
1081 }