Salome HOME
Updated copyright comment
[modules/hexablock.git] / src / HEXABLOCK / HexElements_bis.cxx
old mode 100755 (executable)
new mode 100644 (file)
index fd30415..9b26745
@@ -1,6 +1,25 @@
 
 // C++ : Table d'hexaedres
 
+// 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, 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
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+// Lesser General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License along with this library; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
+//
+// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+//
+
 #include "HexElements.hxx"
 
 #include "HexDocument.hxx"
 #include "HexEdge.hxx"
 
 #include "HexGlobale.hxx"
-#include "HexCylinder.hxx"
+#include "HexNewShape.hxx"
 
 #include <map>
 
-static bool db=false;
-
 BEGIN_NAMESPACE_HEXA
 
 // ====================================================== 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_type != GR_CARTESIAN && grid_type != GR_CYLINDRIC)
+   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_type != GR_CARTESIAN && grid_type != GR_CYLINDRIC)
+   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_type != GR_CARTESIAN && grid_type != GR_CYLINDRIC)
+   else if (grid_nocart)
       return NULL;
 
    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_type != GR_CARTESIAN && grid_type != GR_CYLINDRIC)
+   else if (grid_nocart)
       return NULL;
 
    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_type != GR_CARTESIAN && grid_type != GR_CYLINDRIC)
+   else if (grid_nocart)
       return NULL;
 
    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_type != GR_CARTESIAN && grid_type != GR_CYLINDRIC)
+   else if (grid_nocart)
       return NULL;
 
    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_type != GR_CARTESIAN && grid_type != GR_CYLINDRIC)
+   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_type != GR_CARTESIAN && grid_type != GR_CYLINDRIC)
+   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);
@@ -141,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)
@@ -151,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;
@@ -171,52 +222,66 @@ void Elements::remove ()
        if (tab_hexa[nh] != NULL)
            tab_hexa[nh]->remove();
 }
-// ====================================================== prismQuads
-int Elements::prismQuads (Quads& tstart, Vector* dir, int hauteur)
-{
-   el_root->markAll (NO_USED);
-   int nbcells   = tstart.size ();
-   nbr_vertex    = 0;
-   nbr_edges     = 0;
-
-   nbr_hexas = nbcells*hauteur;
-
-   tab_hexa.resize (nbr_hexas);
-   tab_quad.clear ();          // verticaux
-   tab_edge.clear ();
-   tab_pilier.clear ();
-   tab_vertex.clear ();
-
-   for (int nro=0 ; nro<nbcells ; nro++)
-       {
-       pushHexas (nro, tstart[nro], dir->getDx(), dir->getDy(), dir->getDz(), 
-                  hauteur);
-       }
-   return HOK;
-}
-// ====================================================== pushHexas
-int  Elements::pushHexas (int nro, Quad* qbase, double dx, double dy, 
-                          double dz, int hauteur, Quad* toit, int tcorr[])
+// ====================================================== prismHexas
+int  Elements::prismHexas (int nro, Quad* qbase, int hauteur)
 {
-   int ind_node [QUAD4];
+   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;
           Vertex* nd1 = NULL;
+          double beta = 0;
+          if (revo_lution)
+             {
+             Real3 centre, point, om;
+             revo_center->getPoint (centre);
+             vbase      ->getPoint (point);
+
+             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*revo_axe[dd];
+                 hm [dd] = point  [dd] - ph[dd];
+                 rayon  += hm[dd] * hm[dd];
+                 }
+             subid = grid_geom->addCircle (ph, sqrt(rayon), revo_axe, hm);
+             }
+
           for (int nh=0 ; nh<hauteur ; nh++)
               {
-              nd1 = el_root->addVertex (nd0->getX() + dx, nd0->getY() + dy,
-                                        nd0->getZ() + dz);
+              nd1 = el_root->addVertex (nd0->getX(), nd0->getY(), nd0->getZ());
+              updateMatrix (nh);
+              gen_matrix.perform   (nd1);
               tab_vertex.push_back (nd1);
-              tab_pilier.push_back (newEdge (nd0, nd1));
+              Edge* pilier = newEdge (nd0, nd1);
+              tab_pilier.push_back (pilier);
+              if (revo_lution)
+                 {
+                 double alpha = beta;
+                 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;
               }
           }
@@ -227,7 +292,7 @@ int  Elements::pushHexas (int nro, Quad* qbase, double dx, double dy,
    int ind_poutre [QUAD4];
    for (int ns=0 ; ns<QUAD4 ; ns++)
        {
-       Edge* ebase = qbase ->getEdge (ns);
+       Edge* ebase = qbase->getEdge (ns);
        int   indx  = ebase->getMark ();
        if (indx<0)
           {
@@ -240,7 +305,7 @@ int  Elements::pushHexas (int nro, Quad* qbase, double dx, double dy,
           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];
@@ -248,6 +313,8 @@ int  Elements::pushHexas (int nro, Quad* qbase, double dx, double dy,
               Quad* mur = newQuad (ed0, ed1, ed2, ed3);
               tab_edge.push_back (ed0);
               tab_quad.push_back (mur);
+              if (ebase->isAssociated ())
+                  prismAssociation (ebase, ed0, nh);
               }
           }
        ind_poutre [ns] = indx;
@@ -262,131 +329,94 @@ int  Elements::pushHexas (int nro, Quad* qbase, double dx, double dy,
    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;
        }
    return HOK;
 }
-// ====================================================== cutHexas
-int Elements::cutHexas0 (const Edges& t_edges, int nbcuts)
+// ====================================================== updateMatrix
+void Elements::updateMatrix (int nh)
 {
-                                       // 1) marquage des hexas
-   el_root->markAll (NO_USED);
-                                       // 2) Memo noeuds
-   vector <Quad*> q_amont;
-   vector <Quad*> q_aval;
-
-   int nbnodes = t_edges.size();
-   vector <Vertex*> v_amont (nbnodes);
-   vector <Vertex*> v_aval  (nbnodes);
-
-   std::map <Vertex*, Vertex*> vertex_aval;
-
-   for (int nro=0; nro<nbnodes ; nro++)
-       {
-       Edge*   arete = t_edges [nro];
-       Vertex* amont = arete->getAmont ();
-       v_amont [nro] = amont;
-       vertex_aval [amont] = arete->getAval  ();
-       int nbcells   = arete->getNbrParents ();
-
-       for (int nq=0 ; nq<nbcells ; nq++)
-           {
-           Quad* quad = arete->getParent (nq);
-           if (quad->getMark () != IS_USED)
-              {
-              quad->setMark (IS_USED);
-              int nbcubes = quad->getNbrParents ();
-              for (int nh=0 ; nh<nbcubes ; nh++)
-                  {
-                  Hexa* hexa = quad->getParent (nh);
-                  if (hexa->getMark () != IS_USED)  
-                     {
-                     hexa->setMark (IS_USED);
-                     int namont = hexa->getBase (v_amont[nro], arete);
-                     int naval  = glob->getOpposedQuad (namont);
-                     q_amont.push_back (hexa->getQuad  (namont));
-                     q_aval .push_back (hexa->getQuad  (naval ));
-                     if (db)
-                        {
-                        printf (" Quad = ");
-                        hexa->printName (", ");
-                        printf (" Faces = (");
-                        hexa->getQuad (namont)->printName (", ");
-                        hexa->getQuad (naval )->printName (")\n");
-                        }
-                     }
-                  }
-              }
-           }
-       }
-                   // ------------------- Dimensionnement
-   int hauteur = nbcuts + 1;
-   int nbcells = q_amont.size ();
-   nbr_hexas   = nbcells * hauteur;
-
-   tab_hexa.resize   (nbr_hexas);
-   tab_quad.clear ();          // verticaux
-   tab_edge.clear ();
-   tab_pilier.clear ();
-   tab_vertex.clear ();
+   if (revo_lution)
+      {
+      gen_matrix.defRotation (revo_center, revo_axe, gen_values[nh]);
+      cum_matrix.defRotation (revo_center, revo_axe, cum_values[nh]);
+      }
+   else
+      {
+      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;
+          trans [nc] = prism_dir [nc]*hauteur;
+          }
 
-   for (int ned=0; ned<nbnodes ; ned++)
-       t_edges [ned]->remove ();
+      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]);
 
-   for (int nro=0; nro<nbcells ; nro++)
-       {
-       int tcorr [QUAD4];
-       Quad* sol  = q_amont[nro];
-       Quad* toit = q_aval [nro];
+   nbelts = tab_pilier.size();
+   for (int nro=0 ; nro<nbelts ; nro++)
+       tab_edge.push_back (tab_pilier[nro]);
 
-       for (int nv=0; nv<QUAD4 ; nv++)
-           {
-           Vertex* dest = vertex_aval [sol->getVertex(nv)];
-           tcorr [nv] = toit->indexVertex (dest);
-           }
 
-       pushHexas (nro, sol, 0,0,0, hauteur, toit, tcorr);
-       }
-   return HOK;
+   nbr_hexas  = tab_hexa.size ();
+   nbr_edges  = tab_edge.size ();
+   nbr_quads  = tab_quad.size ();
+   nbr_vertex = tab_vertex.size ();
 }
-// ====================================================== makeCylinder
-int Elements::makeCylinder (Cylinder* cyl, Vector* vx, int nr, int na, int nl)
+// ====================================================== getShape
+NewShape* Elements::getShape()
 {
-   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, true);
-   fillGrid ();
-   return HOK;
+   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;
 }
-// ====================================================== makePipe
-int Elements::makePipe (Cylinder* cyl, Vector* vx, int nr, int na, int nl)
+// ====================================================== closeShape
+void Elements::closeShape()
 {
-   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 ();
-   return HOK;
-}
+   if (grid_geom==NULL)
+       return;
 
+   grid_geom -> closeShape();
+   grid_geom = NULL;
+}
 END_NAMESPACE_HEXA