Salome HOME
Updated copyright comment
[modules/hexablock.git] / src / HEXABLOCK / HexElements_bis.cxx
old mode 100755 (executable)
new mode 100644 (file)
index 0bfc615..9b26745
@@ -1,12 +1,12 @@
 
 // C++ : Table d'hexaedres
 
-// Copyright (C) 2009-2012  CEA/DEN, EDF R&D
+// Copyright (C) 2009-2024  CEA, EDF
 //
 // This library is free software; you can redistribute it and/or
 // modify it under the terms of the GNU Lesser General Public
 // License as published by the Free Software Foundation; either
-// version 2.1 of the License.
+// version 2.1 of the License, or (at your option) any later version.
 //
 // This library is distributed in the hope that it will be useful,
 // but WITHOUT ANY WARRANTY; without even the implied warranty of
 #include "HexEdge.hxx"
 
 #include "HexGlobale.hxx"
-#include "HexCylinder.hxx"
-#include "HexShape.hxx"
+#include "HexNewShape.hxx"
 
 #include <map>
 
 BEGIN_NAMESPACE_HEXA
 
-void geom_dump_asso     (Edge* edge);
-void geom_create_circle (double* milieu, double rayon, double* normale, 
-                         double* base, string& brep);
-
 // ====================================================== getHexaIJK
 Hexa* Elements::getHexaIJK (int nx, int ny, int nz)
 {
+   if (isBad())
+      return NULL;
    if (nx<0 || nx>=size_hx ||  ny<0 || ny>=size_hy || nz<0 || nz>=size_hz)
       return NULL;
    else if (grid_nocart)
       return NULL;
 
-   int nro = nx + size_hx*ny + size_hx*size_hy*nz; 
+   int nro = nx + size_hx*ny + size_hx*size_hy*nz;
 
-   return tab_hexa [nro]; 
+   DumpStart  ("getHexaIJK", nx << ny << nz);
+   DumpReturn (tab_hexa [nro]);
+   return      tab_hexa [nro];
 }
 // ====================================================== getQuadIJ
 Quad* Elements::getQuadIJ (int nx, int ny, int nz)
 {
+   if (isBad())
+      return NULL;
    if (nx<0 || nx>=size_qx ||  ny<0 || ny>=size_qy || nz<0 || nz>=size_qz)
       return NULL;
    else if (grid_nocart)
       return NULL;
 
-   int nro = nx + size_qx*ny + size_qx*size_qy*nz 
+   int nro = nx + size_qx*ny + size_qx*size_qy*nz
                 + size_qx*size_qy*size_qz*dir_z;
-   return tab_quad [nro]; 
+
+   DumpStart  ("getQuadIJ", nx << ny << nz);
+   DumpReturn (tab_quad [nro]);
+   return tab_quad [nro];
 }
 // ====================================================== getQuadJK
 Quad* Elements::getQuadJK (int nx, int ny, int nz)
 {
+   if (isBad())
+      return NULL;
    if (nx<0 || nx>=size_qx ||  ny<0 || ny>=size_qy || nz<0 || nz>=size_qz)
       return NULL;
    else if (grid_nocart)
@@ -74,11 +80,15 @@ Quad* Elements::getQuadJK (int nx, int ny, int nz)
 
    int nro = nx + size_qx*ny + size_qx*size_qy*nz; // + dir_x*...
 
-   return tab_quad [nro]; 
+   DumpStart  ("getQuadJK", nx << ny << nz);
+   DumpReturn (tab_quad [nro]);
+   return tab_quad [nro];
 }
 // ====================================================== getQuadIK
 Quad* Elements::getQuadIK (int nx, int ny, int nz)
 {
+   if (isBad())
+      return NULL;
    if (nx<0 || nx>=size_qx ||  ny<0 || ny>=size_qy || nz<0 || nz>=size_qz)
       return NULL;
    else if (grid_nocart)
@@ -86,11 +96,15 @@ Quad* Elements::getQuadIK (int nx, int ny, int nz)
 
    int nro = nx + size_qx*ny + size_qx*size_qy*nz + size_qx*size_qy*size_qz;
 
-   return tab_quad [nro]; 
+   DumpStart  ("getQuadIK", nx << ny << nz);
+   DumpReturn (tab_quad [nro]);
+   return tab_quad [nro];
 }
 // ====================================================== getEdgeI
 Edge* Elements::getEdgeI (int nx, int ny, int nz)
 {
+   if (isBad())
+      return NULL;
    if (nx<0 || nx>=size_ex ||  ny<0 || ny>=size_ey || nz<0 || nz>=size_ez)
       return NULL;
    else if (grid_nocart)
@@ -98,11 +112,15 @@ Edge* Elements::getEdgeI (int nx, int ny, int nz)
 
    int nro = nx + size_ex*ny + size_ex*size_ey*nz;
 
-   return tab_edge [nro]; 
+   DumpStart  ("getEdgeI", nx << ny << nz);
+   DumpReturn (tab_edge [nro]);
+   return tab_edge [nro];
 }
 // ====================================================== getEdgeJ
 Edge* Elements::getEdgeJ (int nx, int ny, int nz)
 {
+   if (isBad())
+      return NULL;
    if (nx<0 || nx>=size_ex ||  ny<0 || ny>=size_ey || nz<0 || nz>=size_ez)
       return NULL;
    else if (grid_nocart)
@@ -110,47 +128,58 @@ Edge* Elements::getEdgeJ (int nx, int ny, int nz)
 
    int nro = nx + size_ex*ny + size_ex*size_ey*nz + size_ex*size_ey*size_ez;
 
-   return tab_edge [nro]; 
+   DumpStart  ("getEdgeJ", nx << ny << nz);
+   DumpReturn (tab_edge [nro]);
+   return tab_edge [nro];
 }
 // ====================================================== getEdgeK
 Edge* Elements::getEdgeK (int nx, int ny, int nz)
 {
+   if (isBad())
+      return NULL;
    if (nx<0 || nx>=size_ex ||  ny<0 || ny>=size_ey || nz<0 || nz>=size_ez)
       return NULL;
    else if (grid_nocart)
       return NULL;
 
-   int nro = nx + size_ex*ny + size_ex*size_ey*nz 
+   int nro = nx + size_ex*ny + size_ex*size_ey*nz
                 + size_ex*size_ey*size_ez*dir_z;
-   return tab_edge [nro]; 
+
+   DumpStart  ("getEdgeK", nx << ny << nz);
+   DumpReturn (tab_edge [nro]);
+   return tab_edge [nro];
 }
 // ====================================================== getVertexIJK
 Vertex* Elements::getVertexIJK (int nx, int ny, int nz)
 {
+   if (isBad())
+      return NULL;
    if (nx<0 || nx>=size_vx ||  ny<0 || ny>=size_vy || nz<0 || nz>=size_vz)
       return NULL;
    else if (grid_nocart)
       return NULL;
 
-   int nro = nx + size_vx*ny + size_vx*size_vy*nz; 
+   int nro = nx + size_vx*ny + size_vx*size_vy*nz;
 
-   return tab_vertex [nro]; 
+   DumpStart  ("getVertexIJK", nx << ny << nz);
+   DumpReturn (tab_vertex [nro]);
+   return tab_vertex [nro];
 }
 // ====================================================== setVertex
 void Elements::setVertex (Vertex* elt, int nx, int ny, int nz)
 {
-   if (   nx < 0 || nx >= size_vx || ny < 0 || ny >= size_vy 
-       || nz < 0 || nz >= size_vz) return; 
+   if (   nx < 0 || nx >= size_vx || ny < 0 || ny >= size_vy
+       || nz < 0 || nz >= size_vz) return;
 
    int nro = nx + size_vx*ny + size_vx*size_vy*nz;
    tab_vertex [nro] = elt;
 }
 // ====================================================== setVertex (2)
-void Elements::setVertex (int nx, int ny, int nz, double px, double py, 
+void Elements::setVertex (int nx, int ny, int nz, double px, double py,
                                                   double pz)
 {
-   if (   nx < 0 || nx >= size_vx || ny < 0 || ny >= size_vy 
-       || nz < 0 || nz >= size_vz) return; 
+   if (   nx < 0 || nx >= size_vx || ny < 0 || ny >= size_vy
+       || nz < 0 || nz >= size_vz) return;
 
    Vertex*    node = el_root->addVertex (px, py, pz);
    setVertex (node, nx, ny, nz);
@@ -163,7 +192,7 @@ void Elements::setEdge (Edge* elt, EnumCoord dir, int nx, int ny, int nz)
       return;
 
    int nro = nx + size_ex*ny + size_ex*size_ey*nz + size_ex*size_ey*size_ez*dir;
-   tab_edge [nro] = elt; 
+   tab_edge [nro] = elt;
 }
 // ====================================================== setQuad
 void Elements::setQuad (Quad* elt, EnumCoord dir, int nx, int ny, int nz)
@@ -173,13 +202,13 @@ void Elements::setQuad (Quad* elt, EnumCoord dir, int nx, int ny, int nz)
       return;
 
    int nro = nx + size_ex*ny + size_ex*size_ey*nz + size_ex*size_ey*size_ez*dir;
-   tab_quad [nro] = elt; 
+   tab_quad [nro] = elt;
 }
 // ====================================================== setHexa
 void Elements::setHexa (Hexa* elt, int nx, int ny, int nz)
 {
-   if (   nx < 0 || nx >= size_hx || ny < 0 || ny >= size_hy 
-       || nz < 0 || nz >= size_hz) return; 
+   if (   nx < 0 || nx >= size_hx || ny < 0 || ny >= size_hy
+       || nz < 0 || nz >= size_hz) return;
 
    int nro = nx + size_hx*ny + size_hx*size_hy*nz;
    tab_hexa [nro] = elt;
@@ -193,154 +222,23 @@ void Elements::remove ()
        if (tab_hexa[nh] != NULL)
            tab_hexa[nh]->remove();
 }
-// ====================================================== makeCylinder
-int Elements::makeCylinder (Cylinder* cyl, Vector* vx, int nr, int na, int nl)
-{
-   Vertex* orig = cyl->getBase ();
-   Vector* dir  = cyl->getDirection ();
-   double  ray  = cyl->getRadius ();
-   double  haut = cyl->getHeight ();
-   
-   resize (GR_CYLINDRIC, nr, na, nl);
-   cyl_closed = true;
-   makeCylindricalNodes (orig, vx, dir, ray/(nr+1), 360, haut/nl, 
-                         nr, na, nl, true);
-   fillGrid ();
-   assoCylinder (orig, dir, 360);
-   return HOK;
-}
-// ====================================================== makePipe
-int Elements::makePipe (Cylinder* cyl, Vector* vx, int nr, int na, int nl)
-{
-   Vertex* orig = cyl->getBase ();
-   Vector* dir  = cyl->getDirection ();
-   double  ray  = cyl->getRadius ();
-   double  haut = cyl->getHeight ();
-   
-   resize (GR_CYLINDRIC, nr, na, nl);
-   cyl_closed = true;
-   makeCylindricalNodes (orig, vx, dir, ray, 360, haut, nr, na, nl, false);
-   fillGrid ();
-   assoCylinder (orig, dir, 360);
-   return HOK;
-}
-// 
-// ---------------------------------------- prism Quads
-//
-// ====================================================== prismQuads
-int Elements::prismQuads (Quads& tstart, Vector* dir, int nbiter)
-{
-   el_root->markAll (NO_USED);
-   int nbcells   = tstart.size ();
-   nbr_vertex    = 0;
-   nbr_edges     = 0;
-
-   nbr_hexas = nbcells*nbiter;
-
-   tab_hexa.resize (nbr_hexas);
-   tab_quad.clear ();          // verticaux
-   ker_hquad.clear ();         // Horizontaux
-   tab_edge.clear ();
-   tab_pilier.clear ();
-   tab_vertex.clear ();
-
-   revo_lution = false;
-   prism_vec   = false;
-   gen_matrix.defTranslation (dir);
-
-   for (int nro=0 ; nro<nbcells ; nro++)
-       {
-       prismHexas (nro, tstart[nro], nbiter);
-       }
-
-   endPrism ();
-   return HOK;
-}
-// ====================================================== prismQuadsVec
-int Elements::prismQuadsVec (Quads& tstart, Vector* dir, RealVector& tlen, 
-                             int mode)
-{
-   int nbiter = tlen.size();
-   if (nbiter==0)
-      return HERR;
-
-   el_root->markAll (NO_USED);
-   int nbcells   = tstart.size ();
-   nbr_vertex    = 0;
-   nbr_edges     = 0;
-
-   nbr_hexas = nbcells*nbiter;
-
-   tab_hexa.resize (nbr_hexas);
-   tab_quad.clear ();          // verticaux
-   ker_hquad.clear ();         // Horizontaux
-   tab_edge.clear ();
-   tab_pilier.clear ();
-   tab_vertex.clear ();
-
-   revo_lution = false;
-   prism_vec   = true;
-   dir->getCoord  (prism_dir);
-   normer_vecteur (prism_dir);
-   gen_values = tlen;
-
-   for (int nro=0 ; nro<nbcells ; nro++)
-       {
-       prismHexas (nro, tstart[nro], nbiter);
-       }
-
-   endPrism ();
-   return HOK;
-}
-// ======================================================== revolutionQuads
-int Elements::revolutionQuads (Quads& start, Vertex* center, Vector* axis,
-                               RealVector &angles)
-{
-   int nbiter = angles.size();
-   if (center==NULL  || axis==NULL || nbiter==0)
-      return HERR;
-
-   el_root->markAll (NO_USED);
-   int nbcells   = start.size ();
-   nbr_vertex    = 0;
-   nbr_edges     = 0;
-
-   nbr_hexas   = nbcells*nbiter;
-
-   tab_hexa.resize (nbr_hexas);
-   tab_quad.clear ();          // verticaux
-   ker_hquad.clear ();         // Horizontaux
-   tab_edge.clear ();
-   tab_pilier.clear ();
-   tab_vertex.clear ();
-
-   revo_lution  = true;
-   prism_vec    = false;
-   revo_axis    = axis;
-   revo_center  = center;
-   gen_values = angles;
-
-   for (int nro=0 ; nro<nbcells ; nro++)
-       {
-       prismHexas (nro, start[nro], nbiter);
-       }
-
-   endPrism ();
-   return HOK;
-}
 // ====================================================== prismHexas
 int  Elements::prismHexas (int nro, Quad* qbase, int hauteur)
 {
-   int ind_node [QUAD4];
-   string c_rep;
+   int   ind_node [QUAD4];
+   int   subid = 0;
+   Real3 koord, transfo;
 
            // ----------------------------- Vertex + aretes verticales
    for (int ns=0 ; ns<QUAD4 ; ns++)
        {
-       Vertex* vbase = qbase ->getVertex (ns);
+       Vertex* vbase = qbase->getVertex (ns);
        int     indx  = vbase->getMark ();
        if (indx<0)
           {
+          bool asso = vbase->isAssociated();
+          vbase->getAssoCoord (koord);
+
           indx = nbr_vertex++;
           vbase->setMark (indx);
           Vertex* nd0 = vbase;
@@ -348,33 +246,21 @@ int  Elements::prismHexas (int nro, Quad* qbase, int hauteur)
           double beta = 0;
           if (revo_lution)
              {
-             Real3 centre, vk, point, om;
+             Real3 centre, point, om;
              revo_center->getPoint (centre);
              vbase      ->getPoint (point);
-             revo_axis  ->getCoord (vk);
-             normer_vecteur (vk);
 
-             calc_vecteur   (centre, point, om);
-             double oh     = prod_scalaire (om, vk);
-             double rayon  = 0;
+             calc_vecteur  (centre, point, om);
+             double oh    = prod_scalaire (om, revo_axe);
+             double rayon = 0;
              Real3  ph, hm;
              for (int dd=dir_x; dd<=dir_z ; dd++)
                  {
-                 ph [dd] = centre [dd] + oh*vk[dd]; 
+                 ph [dd] = centre [dd] + oh*revo_axe[dd];
                  hm [dd] = point  [dd] - ph[dd];
                  rayon  += hm[dd] * hm[dd];
                  }
-             rayon = sqrt (rayon);
-/********************************
-             PutCoord (centre);
-             PutCoord (point);
-             PutData  (oh);
-             PutCoord (ph);
-             PutData  (rayon);
-             PutCoord (vk);
-             PutCoord (hm);
-********************************/
-             geom_create_circle (ph, rayon, vk, hm, c_rep);
+             subid = grid_geom->addCircle (ph, sqrt(rayon), revo_axe, hm);
              }
 
           for (int nh=0 ; nh<hauteur ; nh++)
@@ -388,11 +274,13 @@ int  Elements::prismHexas (int nro, Quad* qbase, int hauteur)
               if (revo_lution)
                  {
                  double alpha = beta;
-                 beta = alpha + gen_values[nh]; 
-                 Shape* shape = new Shape (c_rep);
-                 shape->setBounds (alpha/360, beta/360);
-                 pilier->addAssociation (shape);
-                 //      geom_dump_asso (pilier);
+                 beta = alpha + gen_values[nh];
+                 grid_geom->addAssociation (pilier, subid, alpha/360, beta/360);
+                 }
+              if (asso)
+                 {
+                 cum_matrix.perform  (koord, transfo);
+                 nd1->setAssociation (transfo);
                  }
               nd0 = nd1;
               }
@@ -417,7 +305,7 @@ int  Elements::prismHexas (int nro, Quad* qbase, int hauteur)
           for (int nh=0 ; nh<hauteur ; nh++)
               {
               ed2 = ed0;
-              ed0 = newEdge (tab_vertex [nd1*hauteur + nh], 
+              ed0 = newEdge (tab_vertex [nd1*hauteur + nh],
                              tab_vertex [nd2*hauteur + nh]);
               ed1 = tab_pilier [nd1*hauteur + nh];
               ed3 = tab_pilier [nd2*hauteur + nh];
@@ -425,7 +313,8 @@ int  Elements::prismHexas (int nro, Quad* qbase, int hauteur)
               Quad* mur = newQuad (ed0, ed1, ed2, ed3);
               tab_edge.push_back (ed0);
               tab_quad.push_back (mur);
-              prismAssociation (ed2, ed0, nh, ed1);
+              if (ebase->isAssociated ())
+                  prismAssociation (ebase, ed0, nh);
               }
           }
        ind_poutre [ns] = indx;
@@ -440,14 +329,14 @@ int  Elements::prismHexas (int nro, Quad* qbase, int hauteur)
    int nv3 = hauteur*ind_poutre [3];
    for (int nh=0 ; nh<hauteur ; nh++)
        {
-       qb = newQuad (tab_edge [nh+nv0], tab_edge [nh+nv1], 
+       qb = newQuad (tab_edge [nh+nv0], tab_edge [nh+nv1],
                      tab_edge [nh+nv2], tab_edge [nh+nv3]);
        qc = tab_quad [nh + nv0];
        qd = tab_quad [nh + nv2];
        qe = tab_quad [nh + nv1];
        qf = tab_quad [nh + nv3];
 
-// *** tab_hexa [nh*hauteur + nro] = newHexa (qa, qb, qc, qd, qe, qf); Abu 
+// *** tab_hexa [nh*hauteur + nro] = newHexa (qa, qb, qc, qd, qe, qf); Abu
        tab_hexa [nro*hauteur + nh] = newHexa (qa, qb, qc, qd, qe, qf);
        ker_hquad.push_back (qb);
        qa = qb;
@@ -455,25 +344,46 @@ int  Elements::prismHexas (int nro, Quad* qbase, int hauteur)
    return HOK;
 }
 // ====================================================== updateMatrix
-void Elements::updateMatrix (int hauteur)
+void Elements::updateMatrix (int nh)
 {
    if (revo_lution)
       {
-      gen_matrix.defRotation (revo_center, revo_axis, gen_values[hauteur]);
+      gen_matrix.defRotation (revo_center, revo_axe, gen_values[nh]);
+      cum_matrix.defRotation (revo_center, revo_axe, cum_values[nh]);
       }
-   else if (prism_vec)
+   else
       {
-      double h0 = hauteur>0 ?  gen_values[hauteur-1] : 0;
-      double dh = gen_values[hauteur] - h0;
-      Real3 decal;
+      double hauteur = (nh+1)*prism_len;
+      double dh      = prism_len;
+      if (prism_vec)
+         {
+         if (under_v6)
+            {
+            hauteur = gen_values[nh];
+            dh      = nh>0 ? hauteur-gen_values[nh-1] : hauteur;
+            }
+         else
+            {
+            hauteur = cum_values[nh];
+            dh      = gen_values[nh];
+            }
+         }
+      Real3 trans, decal;
       for (int nc=dir_x ; nc<=dir_z ; nc++)
-          decal [nc] = prism_dir [nc]*dh; 
+          {
+          decal [nc] = prism_dir [nc]*dh;
+          trans [nc] = prism_dir [nc]*hauteur;
+          }
+
       gen_matrix.defTranslation (decal);
+      cum_matrix.defTranslation (trans);
       }
 }
 // ====================================================== endPrism
 void Elements::endPrism ()
 {
+   closeShape();
+
    int nbelts = ker_hquad.size();
    for (int nro=0 ; nro<nbelts ; nro++)
        tab_quad.push_back (ker_hquad[nro]);
@@ -488,4 +398,25 @@ void Elements::endPrism ()
    nbr_quads  = tab_quad.size ();
    nbr_vertex = tab_vertex.size ();
 }
+// ====================================================== getShape
+NewShape* Elements::getShape()
+{
+   if (grid_geom==NULL)
+      {
+      std::string name = "extrud_" + el_name;
+      grid_geom   = el_root->addShape (name.c_str(), SH_EXTRUD);
+      grid_geom -> openShape();
+      }
+
+   return grid_geom;
+}
+// ====================================================== closeShape
+void Elements::closeShape()
+{
+   if (grid_geom==NULL)
+       return;
+
+   grid_geom -> closeShape();
+   grid_geom = NULL;
+}
 END_NAMESPACE_HEXA