src/SMESH/SMESH_Hexa_3D.cxx

   1 //  SMESH SMESH : implementaion of SMESH idl descriptions
   2 //
   3 //  Copyright (C) 2003  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
   4 //  CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
   5 //
   6 //  This library is free software; you can redistribute it and/or
   7 //  modify it under the terms of the GNU Lesser General Public
   8 //  License as published by the Free Software Foundation; either
   9 //  version 2.1 of the License.
  10 //
  11 //  This library is distributed in the hope that it will be useful,
  12 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
  13 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14 //  Lesser General Public License for more details.
  15 //
  16 //  You should have received a copy of the GNU Lesser General Public
  17 //  License along with this library; if not, write to the Free Software
  18 //  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  19 //
  20 //  See http://www.opencascade.org/SALOME/ or email : webmaster.salome@opencascade.org
  21 //
  22 //
  23 //
  24 //  File   : SMESH_Hexa_3D.cxx
  25 //  Author : Paul RASCLE, EDF
  26 //  Module : SMESH
  27 //  $Header$
  28
  29 using namespace std;
  30 using namespace std;
  31 #include "SMESH_Hexa_3D.hxx"
  32 #include "SMESH_Quadrangle_2D.hxx"
  33 #include "SMESH_Gen.hxx"
  34 #include "SMESH_Mesh.hxx"
  35
  36 #include "SMDS_MeshElement.hxx"
  37 #include "SMDS_MeshNode.hxx"
  38 #include "SMDS_FacePosition.hxx"
  39
  40 #include <TopExp.hxx>
  41 #include <TopTools_IndexedDataMapOfShapeListOfShape.hxx>
  42 #include <TopTools_ListOfShape.hxx>
  43 #include <TopTools_ListIteratorOfListOfShape.hxx>
  44 #include <TColStd_ListIteratorOfListOfInteger.hxx>
  45
  46 #include <BRep_Tool.hxx>
  47 #include <Geom_Surface.hxx>
  48 #include <Geom_Curve.hxx>
  49 #include <Geom2d_Curve.hxx>
  50 #include <Handle_Geom2d_Curve.hxx>
  51 #include <Handle_Geom_Curve.hxx>
  52
  53 #include "utilities.h"
  54
  55
  56 //=============================================================================
  57 /*!
  58  *
  59  */
  60 //=============================================================================
  61
  62 SMESH_Hexa_3D::SMESH_Hexa_3D(int hypId, int studyId,
  63                              SMESH_Gen* gen)
  64   : SMESH_3D_Algo(hypId, studyId, gen)
  65 {
  66   MESSAGE("SMESH_Hexa_3D::SMESH_Hexa_3D");
  67   _name = "Hexa_3D";
  68 //   _shapeType = TopAbs_SOLID;
  69   _shapeType = (1 << TopAbs_SHELL) | (1 << TopAbs_SOLID);// 1 bit /shape type
  70 //   MESSAGE("_shapeType octal " << oct << _shapeType);
  71   for (int i=0; i<6; i++) _quads[i] = 0;
  72 }
  73
  74 //=============================================================================
  75 /*!
  76  *
  77  */
  78 //=============================================================================
  79
  80 SMESH_Hexa_3D::~SMESH_Hexa_3D()
  81 {
  82   MESSAGE("SMESH_Hexa_3D::~SMESH_Hexa_3D");
  83 }
  84
  85 //=============================================================================
  86 /*!
  87  *
  88  */
  89 //=============================================================================
  90
  91 bool SMESH_Hexa_3D::CheckHypothesis(SMESH_Mesh& aMesh,
  92                                     const TopoDS_Shape& aShape)
  93 {
  94   MESSAGE("SMESH_Hexa_3D::CheckHypothesis");
  95
  96   bool isOk = true;
  97
  98   // nothing to check
  99
 100   return isOk;
 101 }
 102
 103 //=============================================================================
 104 /*!
 105  * Hexahedron mesh on hexaedron like form
 106  * -0.  - shape and face mesh verification
 107  * -1.  - identify faces and vertices of the "cube"
 108  * -2.  - Algorithm from:
 109  * "Application de l'interpolation transfinie à la création de maillages
 110  *  C0 ou G1 continus sur des triangles, quadrangles, tetraedres, pentaedres
 111  *  et hexaedres déformés."
 112  * Alain PERONNET - 8 janvier 1999
 113  */
 114 //=============================================================================
 115
 116 bool SMESH_Hexa_3D::Compute(SMESH_Mesh& aMesh,
 117                             const TopoDS_Shape& aShape)
 118   throw (SALOME_Exception)
 119 {
 120   MESSAGE("SMESH_Hexa_3D::Compute");
 121
 122   bool isOk = false;
 123   const Handle(SMESHDS_Mesh)& meshDS = aMesh.GetMeshDS();
 124   SMESH_subMesh* theSubMesh = aMesh.GetSubMesh(aShape);
 125   //const Handle(SMESHDS_SubMesh)& subMeshDS = theSubMesh->GetSubMeshDS();
 126
 127   // 0.  - shape and face mesh verification
 128   // 0.1 - shape must be a solid (or a shell) with 6 faces
 129   MESSAGE("---");
 130
 131   vector<SMESH_subMesh*> meshFaces;
 132   for (TopExp_Explorer exp(aShape,TopAbs_FACE);exp.More();exp.Next())
 133     {
 134       SMESH_subMesh* aSubMesh = aMesh.GetSubMeshContaining(exp.Current());
 135       ASSERT (aSubMesh);
 136       meshFaces.push_back(aSubMesh);
 137     }
 138   if (meshFaces.size() != 6)
 139     {
 140       SCRUTE(meshFaces.size());
 141       ASSERT(0);
 142       return false;
 143     }
 144
 145   // 0.2 - is each face meshed with Quadrangle_2D? (so, with a wire of 4 edges)
 146   MESSAGE("---");
 147
 148   for (int i=0; i<6; i++)
 149     {
 150       TopoDS_Shape aShape = meshFaces[i]->GetSubShape();
 151       SMESH_Algo* algo = _gen->GetAlgo(aMesh, aShape);
 152       string algoName = algo->GetName();
 153       if (algoName != "Quadrangle_2D")
 154         {
 155           // *** delete _quads
 156           SCRUTE(algoName);
 157           ASSERT(0);
 158           return false;
 159         }
 160       SMESH_Quadrangle_2D* quadAlgo =dynamic_cast<SMESH_Quadrangle_2D*> (algo);
 161       ASSERT(quadAlgo);
 162       try
 163         {
 164           _quads[i] = quadAlgo->CheckAnd2Dcompute(aMesh, aShape);
 165           // *** to delete after usage
 166         }
 167       catch (SALOME_Exception& S_ex)
 168         {
 169           // *** delete _quads
 170           // *** throw exception
 171           ASSERT(0);
 172         }
 173     }
 174
 175   // 1.  - identify faces and vertices of the "cube"
 176   // 1.1 - ancestor maps vertex->edges in the cube
 177   MESSAGE("---");
 178
 179   TopTools_IndexedDataMapOfShapeListOfShape MS;
 180   TopExp::MapShapesAndAncestors(aShape, TopAbs_VERTEX, TopAbs_EDGE, MS);
 181
 182   // 1.2 - first face is choosen as face Y=0 of the unit cube
 183   MESSAGE("---");
 184
 185   const TopoDS_Shape& aFace = meshFaces[0]->GetSubShape();
 186   const TopoDS_Face& F = TopoDS::Face(aFace);
 187
 188   // 1.3 - identify the 4 vertices of the face Y=0: V000, V100, V101, V001
 189   MESSAGE("---");
 190
 191   int i = 0;
 192   TopoDS_Edge E = _quads[0]->edge[i]; //edge will be Y=0,Z=0 on unit cube
 193   double f,l;
 194   Handle(Geom2d_Curve) C2d = BRep_Tool::CurveOnSurface(E,F,f,l);
 195   TopoDS_Vertex VFirst, VLast;
 196   TopExp::Vertices(E, VFirst, VLast); // corresponds to f and l
 197   bool isForward = (((l-f)*(_quads[0]->last[i] - _quads[0]->first[i])) > 0);
 198
 199   if (isForward)
 200     {
 201       _cube.V000 = VFirst; // will be (0,0,0) on the unit cube
 202       _cube.V100 = VLast;  // will be (1,0,0) on the unit cube
 203     }
 204   else
 205     {
 206       _cube.V000 = VLast;
 207       _cube.V100 = VFirst;
 208     }
 209
 210   i = 1;
 211   E = _quads[0]->edge[i];
 212   C2d = BRep_Tool::CurveOnSurface(E,F,f,l);
 213   TopExp::Vertices(E, VFirst, VLast);
 214   isForward = (((l-f)*(_quads[0]->last[i] - _quads[0]->first[i])) > 0);
 215   if (isForward) _cube.V101 = VLast;  // will be (1,0,1) on the unit cube
 216   else _cube.V101 = VFirst;
 217
 218   i = 2;
 219   E = _quads[0]->edge[i];
 220   C2d = BRep_Tool::CurveOnSurface(E,F,f,l);
 221   TopExp::Vertices(E, VFirst, VLast);
 222   isForward = (((l-f)*(_quads[0]->last[i] - _quads[0]->first[i])) > 0);
 223   if (isForward) _cube.V001 = VLast;  // will be (0,0,1) on the unit cube
 224   else _cube.V001 = VFirst;
 225
 226   // 1.4 - find edge X=0, Z=0 (ancestor of V000 not in face Y=0)
 227   //     - find edge X=1, Z=0 (ancestor of V100 not in face Y=0)
 228   //     - find edge X=1, Z=1 (ancestor of V101 not in face Y=0)
 229   //     - find edge X=0, Z=1 (ancestor of V001 not in face Y=0)
 230   MESSAGE("---");
 231
 232   TopoDS_Edge E_0Y0 = EdgeNotInFace(aMesh, aShape, F, _cube.V000, MS);
 233   ASSERT(! E_0Y0.IsNull());
 234
 235   TopoDS_Edge E_1Y0 = EdgeNotInFace(aMesh, aShape, F, _cube.V100, MS);
 236   ASSERT(! E_1Y0.IsNull());
 237
 238   TopoDS_Edge E_1Y1 = EdgeNotInFace(aMesh, aShape, F, _cube.V101, MS);
 239   ASSERT(! E_1Y1.IsNull());
 240
 241   TopoDS_Edge E_0Y1 = EdgeNotInFace(aMesh, aShape, F, _cube.V001, MS);
 242   ASSERT(! E_0Y1.IsNull());
 243
 244   // 1.5 - identify the 4 vertices in face Y=1: V010, V110, V111, V011
 245   MESSAGE("---");
 246
 247   TopExp::Vertices(E_0Y0, VFirst, VLast);
 248   if (VFirst.IsSame(_cube.V000)) _cube.V010 = VLast;
 249   else _cube.V010 = VFirst;
 250
 251   TopExp::Vertices(E_1Y0, VFirst, VLast);
 252   if (VFirst.IsSame(_cube.V100)) _cube.V110 = VLast;
 253   else _cube.V110 = VFirst;
 254
 255   TopExp::Vertices(E_1Y1, VFirst, VLast);
 256   if (VFirst.IsSame(_cube.V101)) _cube.V111 = VLast;
 257   else _cube.V111 = VFirst;
 258
 259   TopExp::Vertices(E_0Y1, VFirst, VLast);
 260   if (VFirst.IsSame(_cube.V001)) _cube.V011 = VLast;
 261   else _cube.V011 = VFirst;
 262
 263   // 1.6 - find remaining faces given 4 vertices
 264   MESSAGE("---");
 265
 266   _indY0 = 0;
 267   _cube.quad_Y0 = _quads[_indY0];
 268
 269   _indY1 = GetFaceIndex(aMesh, aShape, meshFaces,
 270                         _cube.V010,_cube.V011,_cube.V110,_cube.V111);
 271   _cube.quad_Y1 = _quads[_indY1];
 272
 273   _indZ0 = GetFaceIndex(aMesh, aShape, meshFaces,
 274                         _cube.V000,_cube.V010,_cube.V100,_cube.V110);
 275   _cube.quad_Z0 = _quads[_indZ0];
 276
 277   _indZ1 = GetFaceIndex(aMesh, aShape, meshFaces,
 278                   _cube.V001,_cube.V011,_cube.V101,_cube.V111);
 279   _cube.quad_Z1 = _quads[_indZ1];
 280
 281   _indX0 = GetFaceIndex(aMesh, aShape, meshFaces,
 282                         _cube.V000,_cube.V001,_cube.V010,_cube.V011);
 283   _cube.quad_X0 = _quads[_indX0];
 284
 285   _indX1 = GetFaceIndex(aMesh, aShape, meshFaces,
 286                         _cube.V100,_cube.V101,_cube.V110,_cube.V111);
 287   _cube.quad_X1 = _quads[_indX1];
 288
 289   MESSAGE("---");
 290
 291   // 1.7 - get convertion coefs from face 2D normalized to 3D normalized
 292
 293   Conv2DStruct cx0; // for face X=0
 294   Conv2DStruct cx1; // for face X=1
 295   Conv2DStruct cy0;
 296   Conv2DStruct cy1;
 297   Conv2DStruct cz0;
 298   Conv2DStruct cz1;
 299
 300   GetConv2DCoefs(*_cube.quad_X0, meshFaces[_indX0]->GetSubShape(),
 301                  _cube.V000,_cube.V010,_cube.V011,_cube.V001,
 302                  cx0);
 303   GetConv2DCoefs(*_cube.quad_X1, meshFaces[_indX1]->GetSubShape(),
 304                  _cube.V100,_cube.V110,_cube.V111,_cube.V101,
 305                  cx1);
 306   GetConv2DCoefs(*_cube.quad_Y0, meshFaces[_indY0]->GetSubShape(),
 307                  _cube.V000,_cube.V100,_cube.V101,_cube.V001,
 308                  cy0);
 309   GetConv2DCoefs(*_cube.quad_Y1, meshFaces[_indY1]->GetSubShape(),
 310                  _cube.V010,_cube.V110,_cube.V111,_cube.V011,
 311                  cy1);
 312   GetConv2DCoefs(*_cube.quad_Z0, meshFaces[_indZ0]->GetSubShape(),
 313                  _cube.V000,_cube.V100,_cube.V110,_cube.V010,
 314                  cz0);
 315   GetConv2DCoefs(*_cube.quad_Z1, meshFaces[_indZ1]->GetSubShape(),
 316                  _cube.V001,_cube.V101,_cube.V111,_cube.V011,
 317                  cz1);
 318
 319   // 1.8 - create a 3D structure for normalized values
 320
 321   MESSAGE("---");
 322   int nbx = _cube.quad_Y0->nbPts[0];
 323   int nby = _cube.quad_Y0->nbPts[1];
 324   int nbz;
 325   if (cx0.a1 != 0) nbz = _cube.quad_X0->nbPts[1];
 326   else  nbz = _cube.quad_X0->nbPts[0];
 327   //SCRUTE(nbx);
 328   //SCRUTE(nby);
 329   //SCRUTE(nbz);
 330   int nbxyz= nbx*nby*nbz;
 331   Point3DStruct* np = new Point3DStruct[nbxyz];
 332
 333   // 1.9 - store node indexes of faces
 334
 335   {
 336     const TopoDS_Face& F = TopoDS::Face(meshFaces[_indX0]->GetSubShape());
 337     const TColStd_ListOfInteger& indElt
 338       = aMesh.GetSubMesh(F)->GetSubMeshDS()->GetIDNodes();
 339     TColStd_ListIteratorOfListOfInteger itf(indElt);
 340
 341     faceQuadStruct* quad = _cube.quad_X0;
 342     int i=0; // j = x/face , k = y/face
 343     int nbdown = quad->nbPts[0];
 344     int nbright = quad->nbPts[1];
 345
 346     for (; itf.More(); itf.Next())
 347       {
 348         int nodeId = itf.Value();
 349         Handle (SMDS_MeshElement) elt = meshDS->FindNode(nodeId);
 350         Handle (SMDS_MeshNode) node = meshDS->GetNode(1, elt);
 351         Handle (SMDS_FacePosition) fpos
 352           = Handle (SMDS_FacePosition)::DownCast(node->GetPosition());
 353         double ri = fpos->GetUParameter();
 354         double rj = fpos->GetVParameter();
 355         int i1 = int(ri);
 356         int j1 = int(rj);
 357         int ij1 = j1*nbdown +i1;
 358         quad->uv_grid[ij1].nodeId = nodeId;
 359       }
 360
 361     for (int i1=0; i1<nbdown; i1++)
 362       for (int j1=0; j1<nbright; j1++)
 363         {
 364           int ij1 = j1*nbdown +i1;
 365           int j = cx0.ia*i1 + cx0.ib*j1 + cx0.ic; // j = x/face
 366           int k = cx0.ja*i1 + cx0.jb*j1 + cx0.jc; // k = y/face
 367           int ijk = k*nbx*nby + j*nbx + i;
 368           //MESSAGE(" "<<ij1<<" "<<i<<" "<<j<<" "<<ijk);
 369           np[ijk].nodeId = quad->uv_grid[ij1].nodeId;
 370           //SCRUTE(np[ijk].nodeId);
 371         }
 372   }
 373
 374   {
 375     const TopoDS_Face& F = TopoDS::Face(meshFaces[_indX1]->GetSubShape());
 376     const TColStd_ListOfInteger& indElt
 377       = aMesh.GetSubMesh(F)->GetSubMeshDS()->GetIDNodes();
 378     TColStd_ListIteratorOfListOfInteger itf(indElt);
 379
 380     faceQuadStruct* quad = _cube.quad_X1;
 381     int i=nbx-1; // j = x/face , k = y/face
 382     int nbdown = quad->nbPts[0];
 383     int nbright = quad->nbPts[1];
 384
 385     for (; itf.More(); itf.Next())
 386       {
 387         int nodeId = itf.Value();
 388         Handle (SMDS_MeshElement) elt = meshDS->FindNode(nodeId);
 389         Handle (SMDS_MeshNode) node = meshDS->GetNode(1, elt);
 390         Handle (SMDS_FacePosition) fpos
 391           = Handle (SMDS_FacePosition)::DownCast(node->GetPosition());
 392         double ri = fpos->GetUParameter();
 393         double rj = fpos->GetVParameter();
 394         int i1 = int(ri);
 395         int j1 = int(rj);
 396         int ij1 = j1*nbdown +i1;
 397         quad->uv_grid[ij1].nodeId = nodeId;
 398       }
 399
 400     for (int i1=0; i1<nbdown; i1++)
 401       for (int j1=0; j1<nbright; j1++)
 402         {
 403           int ij1 = j1*nbdown +i1;
 404           int j = cx1.ia*i1 + cx1.ib*j1 + cx1.ic; // j = x/face
 405           int k = cx1.ja*i1 + cx1.jb*j1 + cx1.jc; // k = y/face
 406           int ijk = k*nbx*nby + j*nbx + i;
 407           //MESSAGE(" "<<ij1<<" "<<i<<" "<<j<<" "<<ijk);
 408           np[ijk].nodeId = quad->uv_grid[ij1].nodeId;
 409           //SCRUTE(np[ijk].nodeId);
 410         }
 411   }
 412
 413   {
 414     const TopoDS_Face& F = TopoDS::Face(meshFaces[_indY0]->GetSubShape());
 415     const TColStd_ListOfInteger& indElt
 416       = aMesh.GetSubMesh(F)->GetSubMeshDS()->GetIDNodes();
 417     TColStd_ListIteratorOfListOfInteger itf(indElt);
 418
 419     faceQuadStruct* quad = _cube.quad_Y0;
 420     int j=0; // i = x/face , k = y/face
 421     int nbdown = quad->nbPts[0];
 422     int nbright = quad->nbPts[1];
 423
 424     for (; itf.More(); itf.Next())
 425       {
 426         int nodeId = itf.Value();
 427         Handle (SMDS_MeshElement) elt = meshDS->FindNode(nodeId);
 428         Handle (SMDS_MeshNode) node = meshDS->GetNode(1, elt);
 429         Handle (SMDS_FacePosition) fpos
 430           = Handle (SMDS_FacePosition)::DownCast(node->GetPosition());
 431         double ri = fpos->GetUParameter();
 432         double rj = fpos->GetVParameter();
 433         int i1 = int(ri);
 434         int j1 = int(rj);
 435         int ij1 = j1*nbdown +i1;
 436         quad->uv_grid[ij1].nodeId = nodeId;
 437       }
 438
 439     for (int i1=0; i1<nbdown; i1++)
 440       for (int j1=0; j1<nbright; j1++)
 441         {
 442           int ij1 = j1*nbdown +i1;
 443           int i = cy0.ia*i1 + cy0.ib*j1 + cy0.ic; // i = x/face
 444           int k = cy0.ja*i1 + cy0.jb*j1 + cy0.jc; // k = y/face
 445           int ijk = k*nbx*nby + j*nbx + i;
 446           //MESSAGE(" "<<ij1<<" "<<i<<" "<<j<<" "<<ijk);
 447           np[ijk].nodeId = quad->uv_grid[ij1].nodeId;
 448           //SCRUTE(np[ijk].nodeId);
 449         }
 450   }
 451
 452   {
 453     const TopoDS_Face& F = TopoDS::Face(meshFaces[_indY1]->GetSubShape());
 454     const TColStd_ListOfInteger& indElt
 455       = aMesh.GetSubMesh(F)->GetSubMeshDS()->GetIDNodes();
 456     TColStd_ListIteratorOfListOfInteger itf(indElt);
 457
 458     faceQuadStruct* quad = _cube.quad_Y1;
 459     int j=nby-1; // i = x/face , k = y/face
 460     int nbdown = quad->nbPts[0];
 461     int nbright = quad->nbPts[1];
 462
 463     for (; itf.More(); itf.Next())
 464       {
 465         int nodeId = itf.Value();
 466         Handle (SMDS_MeshElement) elt = meshDS->FindNode(nodeId);
 467         Handle (SMDS_MeshNode) node = meshDS->GetNode(1, elt);
 468         Handle (SMDS_FacePosition) fpos
 469           = Handle (SMDS_FacePosition)::DownCast(node->GetPosition());
 470         double ri = fpos->GetUParameter();
 471         double rj = fpos->GetVParameter();
 472         int i1 = int(ri);
 473         int j1 = int(rj);
 474         int ij1 = j1*nbdown +i1;
 475         quad->uv_grid[ij1].nodeId = nodeId;
 476       }
 477
 478     for (int i1=0; i1<nbdown; i1++)
 479       for (int j1=0; j1<nbright; j1++)
 480         {
 481           int ij1 = j1*nbdown +i1;
 482           int i = cy1.ia*i1 + cy1.ib*j1 + cy1.ic; // i = x/face
 483           int k = cy1.ja*i1 + cy1.jb*j1 + cy1.jc; // k = y/face
 484           int ijk = k*nbx*nby + j*nbx + i;
 485           //MESSAGE(" "<<ij1<<" "<<i<<" "<<j<<" "<<ijk);
 486           np[ijk].nodeId = quad->uv_grid[ij1].nodeId;
 487           //SCRUTE(np[ijk].nodeId);
 488         }
 489   }
 490
 491   {
 492     const TopoDS_Face& F = TopoDS::Face(meshFaces[_indZ0]->GetSubShape());
 493     const TColStd_ListOfInteger& indElt
 494       = aMesh.GetSubMesh(F)->GetSubMeshDS()->GetIDNodes();
 495     TColStd_ListIteratorOfListOfInteger itf(indElt);
 496
 497     faceQuadStruct* quad = _cube.quad_Z0;
 498     int k=0; // i = x/face , j = y/face
 499     int nbdown = quad->nbPts[0];
 500     int nbright = quad->nbPts[1];
 501
 502     for (; itf.More(); itf.Next())
 503       {
 504         int nodeId = itf.Value();
 505         Handle (SMDS_MeshElement) elt = meshDS->FindNode(nodeId);
 506         Handle (SMDS_MeshNode) node = meshDS->GetNode(1, elt);
 507         Handle (SMDS_FacePosition) fpos
 508           = Handle (SMDS_FacePosition)::DownCast(node->GetPosition());
 509         double ri = fpos->GetUParameter();
 510         double rj = fpos->GetVParameter();
 511         int i1 = int(ri);
 512         int j1 = int(rj);
 513         int ij1 = j1*nbdown +i1;
 514         quad->uv_grid[ij1].nodeId = nodeId;
 515       }
 516
 517     for (int i1=0; i1<nbdown; i1++)
 518       for (int j1=0; j1<nbright; j1++)
 519         {
 520           int ij1 = j1*nbdown +i1;
 521           int i = cz0.ia*i1 + cz0.ib*j1 + cz0.ic; // i = x/face
 522           int j = cz0.ja*i1 + cz0.jb*j1 + cz0.jc; // j = y/face
 523           int ijk = k*nbx*nby + j*nbx + i;
 524           //MESSAGE(" "<<ij1<<" "<<i<<" "<<j<<" "<<ijk);
 525           np[ijk].nodeId = quad->uv_grid[ij1].nodeId;
 526           //SCRUTE(np[ijk].nodeId);
 527         }
 528   }
 529
 530   {
 531     const TopoDS_Face& F = TopoDS::Face(meshFaces[_indZ1]->GetSubShape());
 532     const TColStd_ListOfInteger& indElt
 533       = aMesh.GetSubMesh(F)->GetSubMeshDS()->GetIDNodes();
 534     TColStd_ListIteratorOfListOfInteger itf(indElt);
 535
 536     faceQuadStruct* quad = _cube.quad_Z1;
 537     int k=nbz-1; // i = x/face , j = y/face
 538     int nbdown = quad->nbPts[0];
 539     int nbright = quad->nbPts[1];
 540
 541     for (; itf.More(); itf.Next())
 542       {
 543         int nodeId = itf.Value();
 544         Handle (SMDS_MeshElement) elt = meshDS->FindNode(nodeId);
 545         Handle (SMDS_MeshNode) node = meshDS->GetNode(1, elt);
 546         Handle (SMDS_FacePosition) fpos
 547           = Handle (SMDS_FacePosition)::DownCast(node->GetPosition());
 548         double ri = fpos->GetUParameter();
 549         double rj = fpos->GetVParameter();
 550         int i1 = int(ri);
 551         int j1 = int(rj);
 552         int ij1 = j1*nbdown +i1;
 553         quad->uv_grid[ij1].nodeId = nodeId;
 554       }
 555
 556     for (int i1=0; i1<nbdown; i1++)
 557       for (int j1=0; j1<nbright; j1++)
 558         {
 559           int ij1 = j1*nbdown +i1;
 560           int i = cz1.ia*i1 + cz1.ib*j1 + cz1.ic; // i = x/face
 561           int j = cz1.ja*i1 + cz1.jb*j1 + cz1.jc; // j = y/face
 562           int ijk = k*nbx*nby + j*nbx + i;
 563           //MESSAGE(" "<<ij1<<" "<<i<<" "<<j<<" "<<ijk);
 564           np[ijk].nodeId = quad->uv_grid[ij1].nodeId;
 565           //SCRUTE(np[ijk].nodeId);
 566         }
 567   }
 568
 569   // 2.0 - for each node of the cube:
 570   //       - get the 8 points 3D = 8 vertices of the cube
 571   //       - get the 12 points 3D on the 12 edges of the cube
 572   //       - get the 6 points 3D on the 6 faces with their ID
 573   //       - compute the point 3D
 574   //       - store the point 3D in SMESHDS, store its ID in 3D structure
 575
 576   TopoDS_Shell aShell;
 577   TopExp_Explorer exp(aShape,TopAbs_SHELL);
 578   if (exp.More())
 579     {
 580       aShell = TopoDS::Shell(exp.Current());
 581     }
 582   else
 583     {
 584       MESSAGE("no shell...");
 585       ASSERT(0);
 586     }
 587
 588   MESSAGE("---");
 589   Pt3 p000, p001, p010, p011, p100, p101, p110, p111;
 590   Pt3 px00, px01, px10, px11;
 591   Pt3 p0y0, p0y1, p1y0, p1y1;
 592   Pt3 p00z, p01z, p10z, p11z;
 593   Pt3 pxy0, pxy1, px0z, px1z, p0yz, p1yz;
 594
 595   GetPoint(p000,     0,     0,     0, nbx, nby, nbz, np, meshDS);
 596   GetPoint(p001,     0,     0, nbz-1, nbx, nby, nbz, np, meshDS);
 597   GetPoint(p010,     0, nby-1,     0, nbx, nby, nbz, np, meshDS);
 598   GetPoint(p011,     0, nby-1, nbz-1, nbx, nby, nbz, np, meshDS);
 599   GetPoint(p100, nbx-1,     0,     0, nbx, nby, nbz, np, meshDS);
 600   GetPoint(p101, nbx-1,     0, nbz-1, nbx, nby, nbz, np, meshDS);
 601   GetPoint(p110, nbx-1, nby-1,     0, nbx, nby, nbz, np, meshDS);
 602   GetPoint(p111, nbx-1, nby-1, nbz-1, nbx, nby, nbz, np, meshDS);
 603
 604   for (int i=1; i<nbx-1; i++)
 605     {
 606       for (int j=1; j<nby-1; j++)
 607         {
 608           for (int k=1; k<nbz-1; k++)
 609             {
 610               // *** seulement maillage regulier
 611               // 12 points on edges
 612               GetPoint(px00,     i,     0,     0, nbx, nby, nbz, np, meshDS);
 613               GetPoint(px01,     i,     0, nbz-1, nbx, nby, nbz, np, meshDS);
 614               GetPoint(px10,     i, nby-1,     0, nbx, nby, nbz, np, meshDS);
 615               GetPoint(px11,     i, nby-1, nbz-1, nbx, nby, nbz, np, meshDS);
 616
 617               GetPoint(p0y0,     0,     j,     0, nbx, nby, nbz, np, meshDS);
 618               GetPoint(p0y1,     0,     j, nbz-1, nbx, nby, nbz, np, meshDS);
 619               GetPoint(p1y0, nbx-1,     j,     0, nbx, nby, nbz, np, meshDS);
 620               GetPoint(p1y1, nbx-1,     j, nbz-1, nbx, nby, nbz, np, meshDS);
 621
 622               GetPoint(p00z,     0,     0,     k, nbx, nby, nbz, np, meshDS);
 623               GetPoint(p01z,     0, nby-1,     k, nbx, nby, nbz, np, meshDS);
 624               GetPoint(p10z, nbx-1,     0,     k, nbx, nby, nbz, np, meshDS);
 625               GetPoint(p11z, nbx-1, nby-1,     k, nbx, nby, nbz, np, meshDS);
 626
 627               // 12 points on faces
 628               GetPoint(pxy0,     i,     j,     0, nbx, nby, nbz, np, meshDS);
 629               GetPoint(pxy1,     i,     j, nbz-1, nbx, nby, nbz, np, meshDS);
 630               GetPoint(px0z,     i,     0,     k, nbx, nby, nbz, np, meshDS);
 631               GetPoint(px1z,     i, nby-1,     k, nbx, nby, nbz, np, meshDS);
 632               GetPoint(p0yz,     0,     j,     k, nbx, nby, nbz, np, meshDS);
 633               GetPoint(p1yz, nbx-1,     j,     k, nbx, nby, nbz, np, meshDS);
 634
 635               int ijk = k*nbx*nby + j*nbx + i;
 636               double x = double(i)/double(nbx-1); // *** seulement
 637               double y = double(j)/double(nby-1); // *** maillage
 638               double z = double(k)/double(nbz-1); // *** regulier
 639
 640               Pt3 X;
 641               for (int i=0; i<3; i++)
 642                 {
 643                   X[i] =
 644                      (1-x)*p0yz[i] + x*p1yz[i]
 645                     +(1-y)*px0z[i] + y*px1z[i]
 646                     +(1-z)*pxy0[i] + z*pxy1[i]
 647                     -(1-x)*((1-y)*p00z[i] + y*p01z[i])
 648                     -   x *((1-y)*p10z[i] + y*p11z[i])
 649                     -(1-y)*((1-z)*px00[i] + z*px01[i])
 650                     -   y *((1-z)*px10[i] + z*px11[i])
 651                     -(1-z)*((1-x)*p0y0[i] + x*p1y0[i])
 652                     -   z *((1-x)*p0y1[i] + x*p1y1[i])
 653                     +(1-x)*( (1-y)*((1-z)*p000[i] + z*p001[i])
 654                             +   y *((1-z)*p010[i] + z*p011[i]))
 655                     +   x *( (1-y)*((1-z)*p100[i] + z*p101[i])
 656                             +   y *((1-z)*p110[i] + z*p111[i]));
 657                 }
 658
 659
 660               int myNodeId = meshDS->AddNode(X[0],X[1],X[2]);
 661               //MESSAGE("point "<<myNodeId<<" "<<X[0]<<" "<<X[1]<<" "<<X[2]);
 662               np[ijk].nodeId = myNodeId;
 663               Handle (SMDS_MeshElement) elt = meshDS->FindNode(myNodeId);
 664               Handle (SMDS_MeshNode) node = meshDS->GetNode(1, elt);
 665
 666               //meshDS->SetNodeInVolume(node, TopoDS::Solid(aShape));
 667               meshDS->SetNodeInVolume(node, aShell);
 668             }
 669         }
 670     }
 671
 672   //2.1 - for each node of the cube (less 3 *1 Faces):
 673   //      - store hexahedron in SMESHDS
 674
 675   for (int i=0; i<nbx-1; i++)
 676     for (int j=0; j<nby-1; j++)
 677       for (int k=0; k<nbz-1; k++)
 678         {
 679           int n1 =    k *nbx*nby +    j *nbx + i;
 680           int n2 =    k *nbx*nby +    j *nbx + i+1;
 681           int n3 =    k *nbx*nby + (j+1)*nbx + i+1;
 682           int n4 =    k *nbx*nby + (j+1)*nbx + i;
 683           int n5 = (k+1)*nbx*nby +    j *nbx + i;
 684           int n6 = (k+1)*nbx*nby +    j *nbx + i+1;
 685           int n7 = (k+1)*nbx*nby + (j+1)*nbx + i+1;
 686           int n8 = (k+1)*nbx*nby + (j+1)*nbx + i;
 687
 688 //        MESSAGE(" "<<n1<<" "<<n2<<" "<<n3<<" "<<n4<<" "<<n5<<" "<<n6<<" "<<n7<<" "<<n8);
 689           //MESSAGE(" "<<np[n1].nodeId<<" "<<np[n2].nodeId<<" "<<np[n3].nodeId<<" "<<np[n4].nodeId<<" "<<np[n5].nodeId<<" "<<np[n6].nodeId<<" "<<np[n7].nodeId<<" "<<np[n8].nodeId);
 690
 691           int hexa = meshDS->AddVolume(np[n1].nodeId,
 692                                        np[n2].nodeId,
 693                                        np[n3].nodeId,
 694                                        np[n4].nodeId,
 695                                        np[n5].nodeId,
 696                                        np[n6].nodeId,
 697                                        np[n7].nodeId,
 698                                        np[n8].nodeId);
 699           Handle (SMDS_MeshElement) elt = meshDS->FindElement(hexa);
 700           meshDS->SetMeshElementOnShape(elt, aShell);
 701
 702           // *** 5 tetrahedres ... verifier orientations,
 703           //     mettre en coherence &vec quadrangles-> triangles
 704           //     choisir afficher 1 parmi edges, face et volumes
 705 //        int tetra1 = meshDS->AddVolume(np[n1].nodeId,
 706 //                                       np[n2].nodeId,
 707 //                                       np[n4].nodeId,
 708 //                                       np[n5].nodeId);
 709 //        int tetra2 = meshDS->AddVolume(np[n2].nodeId,
 710 //                                       np[n3].nodeId,
 711 //                                       np[n4].nodeId,
 712 //                                       np[n7].nodeId);
 713 //        int tetra3 = meshDS->AddVolume(np[n5].nodeId,
 714 //                                       np[n6].nodeId,
 715 //                                       np[n7].nodeId,
 716 //                                       np[n2].nodeId);
 717 //        int tetra4 = meshDS->AddVolume(np[n5].nodeId,
 718 //                                       np[n7].nodeId,
 719 //                                       np[n8].nodeId,
 720 //                                       np[n4].nodeId);
 721 //        int tetra5 = meshDS->AddVolume(np[n5].nodeId,
 722 //                                       np[n7].nodeId,
 723 //                                       np[n2].nodeId,
 724 //                                       np[n4].nodeId);
 725
 726         }
 727
 728   return true;
 729 }
 730
 731 //=============================================================================
 732 /*!
 733  *
 734  */
 735 //=============================================================================
 736
 737 void SMESH_Hexa_3D::GetPoint(Pt3 p,
 738                              int i, int j, int k,
 739                              int nbx, int nby, int nbz,
 740                              Point3DStruct *np,
 741                              const Handle(SMESHDS_Mesh)& meshDS)
 742 {
 743   int ijk = k*nbx*nby + j*nbx + i;
 744   int nodeId = np[ijk].nodeId;
 745   Handle (SMDS_MeshElement) elt = meshDS->FindNode(nodeId);
 746   Handle (SMDS_MeshNode) node = meshDS->GetNode(1, elt);
 747   p[0] = node->X();
 748   p[1] = node->Y();
 749   p[2] = node->Z();
 750   //MESSAGE(" "<<i<<" "<<j<<" "<<k<<" "<<p[0]<<" "<<p[1]<<" "<<p[2]);
 751 }
 752
 753 //=============================================================================
 754 /*!
 755  *
 756  */
 757 //=============================================================================
 758
 759 int SMESH_Hexa_3D::GetFaceIndex(SMESH_Mesh& aMesh,
 760                                 const TopoDS_Shape& aShape,
 761                                 const vector<SMESH_subMesh*>& meshFaces,
 762                                 const TopoDS_Vertex& V0,
 763                                 const TopoDS_Vertex& V1,
 764                                 const TopoDS_Vertex& V2,
 765                                 const TopoDS_Vertex& V3)
 766 {
 767   MESSAGE("SMESH_Hexa_3D::GetFaceIndex");
 768   int faceIndex = -1;
 769   for (int i=1; i<6; i++)
 770     {
 771       const TopoDS_Shape& aFace = meshFaces[i]->GetSubShape();
 772       //const TopoDS_Face& F = TopoDS::Face(aFace);
 773       TopTools_IndexedMapOfShape M;
 774       TopExp::MapShapes(aFace, TopAbs_VERTEX, M);
 775       bool verticesInShape = false;
 776       if (M.Contains(V0))
 777         if (M.Contains(V1))
 778           if (M.Contains(V2))
 779             if (M.Contains(V3)) verticesInShape = true;
 780       if (verticesInShape)
 781         {
 782           faceIndex = i;
 783           break;
 784         }
 785     }
 786   ASSERT(faceIndex>0);
 787   SCRUTE(faceIndex);
 788   return faceIndex;
 789 }
 790
 791 //=============================================================================
 792 /*!
 793  *
 794  */
 795 //=============================================================================
 796
 797 TopoDS_Edge
 798 SMESH_Hexa_3D::EdgeNotInFace(SMESH_Mesh& aMesh,
 799                              const TopoDS_Shape& aShape,
 800                              const TopoDS_Face& aFace,
 801                              const TopoDS_Vertex& aVertex,
 802                              const TopTools_IndexedDataMapOfShapeListOfShape& MS)
 803 {
 804   MESSAGE("SMESH_Hexa_3D::EdgeNotInFace");
 805   TopTools_IndexedDataMapOfShapeListOfShape MF;
 806   TopExp::MapShapesAndAncestors(aFace, TopAbs_VERTEX, TopAbs_EDGE, MF);
 807   const TopTools_ListOfShape& ancestorsInSolid = MS.FindFromKey(aVertex);
 808   const TopTools_ListOfShape& ancestorsInFace  = MF.FindFromKey(aVertex);
 809   SCRUTE(ancestorsInSolid.Extent());
 810   SCRUTE(ancestorsInFace.Extent());
 811   ASSERT(ancestorsInSolid.Extent() == 6); // 6 (edges doublees)
 812   ASSERT(ancestorsInFace.Extent() == 2);
 813
 814   TopoDS_Edge E;
 815   E.Nullify();
 816   TopTools_ListIteratorOfListOfShape its(ancestorsInSolid);
 817   for ( ; its.More();its.Next())
 818   {
 819     TopoDS_Shape ancestor = its.Value();
 820     TopTools_ListIteratorOfListOfShape itf(ancestorsInFace);
 821     bool isInFace = false;
 822     for ( ; itf.More();itf.Next())
 823       {
 824         TopoDS_Shape ancestorInFace = itf.Value();
 825         if (ancestorInFace.IsSame(ancestor))
 826           {
 827             isInFace = true;
 828             break;
 829           }
 830       }
 831     if (! isInFace)
 832       {
 833         E = TopoDS::Edge(ancestor);
 834         break;
 835       }
 836   }
 837   return E;
 838 }
 839
 840 //=============================================================================
 841 /*!
 842  *
 843  */
 844 //=============================================================================
 845
 846   void SMESH_Hexa_3D::GetConv2DCoefs(const faceQuadStruct& quad,
 847                                      const TopoDS_Shape& aShape,
 848                                      const TopoDS_Vertex& V0,
 849                                      const TopoDS_Vertex& V1,
 850                                      const TopoDS_Vertex& V2,
 851                                      const TopoDS_Vertex& V3,
 852                                      Conv2DStruct& conv)
 853 {
 854   MESSAGE("SMESH_Hexa_3D::GetConv2DCoefs");
 855   const TopoDS_Face& F = TopoDS::Face(aShape);
 856   TopoDS_Edge E = quad.edge[0];
 857   double f,l;
 858   Handle(Geom2d_Curve) C2d = BRep_Tool::CurveOnSurface(E,F,f,l);
 859   TopoDS_Vertex VFirst, VLast;
 860   TopExp::Vertices(E, VFirst, VLast); // corresponds to f and l
 861   bool isForward = (((l-f)*(quad.last[0] - quad.first[0])) > 0);
 862   TopoDS_Vertex VA, VB;
 863   if (isForward)
 864     {
 865       VA = VFirst; VB = VLast;
 866     }
 867   else
 868     {
 869      VA = VLast; VB = VFirst;
 870     }
 871   int a1,b1,c1,a2,b2,c2;
 872   if (VA.IsSame(V0))
 873     if (VB.IsSame(V1))
 874       {
 875         a1= 1; b1= 0; c1= 0; // x
 876         a2= 0; b2= 1; c2= 0; // y
 877       }
 878     else
 879       {
 880         ASSERT(VB.IsSame(V3));
 881         a1= 0; b1= 1; c1= 0; // y
 882         a2= 1; b2= 0; c2= 0; // x
 883       }
 884   if (VA.IsSame(V1))
 885     if (VB.IsSame(V2))
 886       {
 887         a1= 0; b1=-1; c1= 1; // 1-y
 888         a2= 1; b2= 0; c2= 0; // x
 889       }
 890     else
 891       {
 892         ASSERT(VB.IsSame(V0));
 893         a1=-1; b1= 0; c1= 1; // 1-x
 894         a2= 0; b2= 1; c2= 0; // y
 895       }
 896   if (VA.IsSame(V2))
 897     if (VB.IsSame(V3))
 898       {
 899         a1=-1; b1= 0; c1= 1; // 1-x
 900         a2= 0; b2=-1; c2= 1; // 1-y
 901       }
 902     else
 903       {
 904         ASSERT(VB.IsSame(V1));
 905         a1= 0; b1=-1; c1= 1; // 1-y
 906         a2=-1; b2= 0; c2= 1; // 1-x
 907       }
 908   if (VA.IsSame(V3))
 909     if (VB.IsSame(V0))
 910       {
 911         a1= 0; b1= 1; c1= 0; // y
 912         a2=-1; b2= 0; c2= 1; // 1-x
 913       }
 914     else
 915       {
 916         ASSERT(VB.IsSame(V2));
 917         a1= 1; b1= 0; c1= 0;  // x
 918         a2= 0; b2=-1; c2= 1;  // 1-y
 919       }
 920   MESSAGE("X = "<<c1<<"+ "<<a1<<"*x + "<<b1<<"*y");
 921   MESSAGE("Y = "<<c2<<"+ "<<a2<<"*x + "<<b2<<"*y");
 922   conv.a1=a1; conv.b1=b1; conv.c1=c1;
 923   conv.a2=a2; conv.b2=b2; conv.c2=c2;
 924
 925   int nbdown = quad.nbPts[0];
 926   int nbright = quad.nbPts[1];
 927   conv.ia = int(a1);
 928   conv.ib = int(b1);
 929   conv.ic = int(c1*a1*a1)*(nbdown-1) + int(c1*b1*b1)*(nbright-1);
 930   conv.ja = int(a2);
 931   conv.jb = int(b2);
 932   conv.jc = int(c2*a2*a2)*(nbdown-1) + int(c2*b2*b2)*(nbright-1);
 933   MESSAGE("I "<<conv.ia<<" "<<conv.ib<<" "<<conv.ic);
 934   MESSAGE("J "<<conv.ja<<" "<<conv.jb<<" "<<conv.jc);
 935 }
 936
 937 //=============================================================================
 938 /*!
 939  *
 940  */
 941 //=============================================================================
 942
 943 ostream & SMESH_Hexa_3D::SaveTo(ostream & save)
 944 {
 945   return save << this;
 946 }
 947
 948 //=============================================================================
 949 /*!
 950  *
 951  */
 952 //=============================================================================
 953
 954 istream & SMESH_Hexa_3D::LoadFrom(istream & load)
 955 {
 956   return load >> (*this);
 957 }
 958
 959 //=============================================================================
 960 /*!
 961  *
 962  */
 963 //=============================================================================
 964
 965 ostream & operator << (ostream & save, SMESH_Hexa_3D & hyp)
 966 {
 967   return save;
 968 }
 969
 970 //=============================================================================
 971 /*!
 972  *
 973  */
 974 //=============================================================================
 975
 976 istream & operator >> (istream & load, SMESH_Hexa_3D & hyp)
 977 {
 978   return load;
 979 }