src/SMESH/SMESH_Hexa_3D.cxx

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