Salome HOME
Update copyrights
[modules/hexablock.git] / src / HEXABLOCK / HexElements.cxx
old mode 100755 (executable)
new mode 100644 (file)
index 730771f..087d093
@@ -1,12 +1,12 @@
 
 // C++ : Grilles
 
-// Copyright (C) 2009-2012  CEA/DEN, EDF R&D
+// Copyright (C) 2009-2019  CEA/DEN, EDF R&D
 //
 // 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 <cmath>
 #include <map>
 
-static bool db=false;
-
 BEGIN_NAMESPACE_HEXA
 
+static bool db=on_debug();
+
 // ====================================================== Constructeur
-Elements::Elements (Document* doc) : EltBase (doc)
+Elements::Elements (Document* doc) : EltBase (doc, EL_GRID)
 {
    glob  = Globale::getInstance ();
 
@@ -54,9 +54,13 @@ Elements::Elements (Document* doc) : EltBase (doc)
    cyl_dispo  = CYL_NOFILL;
    revo_lution = false;
    prism_vec   = false;
+   val_absolues = false;       // Version Hexa6
+   cyl_phimax   = M_PI/2;
+   under_v6     = true;
 }
 // ====================================================== Constructeur
-Elements::Elements (Document* doc, int nx, int ny, int nz) : EltBase (doc)
+Elements::Elements (Document* doc, int nx, int ny, int nz)
+        : EltBase (doc, EL_GRID)
 {
    glob  = Globale::getInstance ();
 
@@ -73,6 +77,8 @@ Elements::Elements (Document* doc, int nx, int ny, int nz) : EltBase (doc)
 
    resize (GR_CYLINDRIC, nx, ny, nz);
    cyl_closed = true;
+   grid_geom  = NULL;
+   cyl_phimax   = M_PI/2;
 }
 // ====================================================== Constructeur (clonage)
 Elements::Elements (Elements* orig) : EltBase (orig->el_root)
@@ -92,6 +98,7 @@ Elements::Elements (Elements* orig) : EltBase (orig->el_root)
 
    resize (orig->grid_type, orig->size_hx, orig->size_hy, orig->size_hz);
    cyl_closed = orig->cyl_closed;
+   cyl_phimax = orig->cyl_phimax;
 }
 // ====================================================== resize
 void Elements::resize (EnumGrid type, int nx, int ny, int nz, int nplus)
@@ -132,7 +139,7 @@ void Elements::resize (EnumGrid type, int nx, int ny, int nz, int nplus)
            return;
 
       case GR_JOINT :
-           nbr_orig = std::max (nx, 1); 
+           nbr_orig = std::max (nx, 1);
            gr_hauteur  = std::max (ny, 1);
            size_hx  = nbr_orig;
            size_hy  = 1;
@@ -147,7 +154,7 @@ void Elements::resize (EnumGrid type, int nx, int ny, int nz, int nplus)
 
       case GR_REPLACE :
            nbr_orig    = std::max (nx, 1);    // nb quads du pattern
-           gr_hauteur  = ny + 1;              // Hauteur des hexas 
+           gr_hauteur  = ny + 1;              // Hauteur des hexas
            pat_nbvertex  = std::max (nz, 1);    // nb vertex du pattern
            pat_nbedges   = std::max (nplus, 1); // nb edges du pattern
            size_hx  = nbr_orig;
@@ -177,180 +184,24 @@ void Elements::resize (EnumGrid type, int nx, int ny, int nz, int nplus)
            break;
       }
 
-   tab_hexa  .resize (nbr_hexas);
-   tab_quad  .resize (nbr_quads);
-   tab_edge  .resize (nbr_edges);
-   tab_vertex.resize (nbr_vertex);
+   tab_hexa  .resize (nbr_hexas,  NULL);
+   tab_quad  .resize (nbr_quads,  NULL);
+   tab_edge  .resize (nbr_edges,  NULL);
+   tab_vertex.resize (nbr_vertex, NULL);
 
    ker_vertex = nbr_vertex;
-
-   for (int nc=0 ; nc< nbr_hexas  ; nc++) tab_hexa   [nc] = NULL;
-   for (int nc=0 ; nc< nbr_quads  ; nc++) tab_quad   [nc] = NULL;
-   for (int nc=0 ; nc< nbr_edges  ; nc++) tab_edge   [nc] = NULL;
-   for (int nc=0 ; nc< nbr_vertex ; nc++) tab_vertex [nc] = NULL;
 }
-
-// ====================================================== makeCartesianGrid
-int Elements::makeCartesianGrid (Vertex* orig, Vector* v1, Vector* v2, 
-                   Vector* v3, int px, int py, int pz, int mx, int my, int mz)
+// ====================================================== saveVtk
+int Elements::saveVtk (cpchar nomfic)
 {
-   if (BadElement (orig) || BadElement(v1) || BadElement(v2) || BadElement(v3)
-                         || px<=0 || py<=0 || pz <= 0
-       || v1->getNorm () <= Epsil || v2->getNorm () <= Epsil
-       || v3->getNorm () <= Epsil)
-      {
-      setError ();
-      return HERR;
-      }
-
-   resize (GR_CARTESIAN, px+mx, py+my, pz+mz);
+   DumpStart ("saveVtk", nomfic);
 
-   makeCartesianNodes (orig, v1, v2, v3, px, py, pz, mx, my, mz);
-
-   fillGrid ();
-   return HOK;
-}
-// ====================================================== makeCylindricalGrid
-int Elements::makeCylindricalGrid (Vertex* c, Vector* b, Vector* h, 
-         double dr, double da, double dl, int nr, int na, int nl, bool fill)
-{
-   if (BadElement (c) || BadElement(b) || BadElement(h)
-       || nr<=0 || na<=0 || nl <= 0 || dr < Epsil || da < Epsil || dl < Epsil 
-       || b->getNorm () <= Epsil || h->getNorm () <= Epsil)
-      {
-      setError ();
-      return HERR;
-      }
-   resize (GR_CYLINDRIC, nr, na, nl);
-   cyl_closed = da >= 360.0;
-   makeCylindricalNodes (c, b, h, dr, da, dl, nr, na, nl, fill);
-   fillGrid ();
-   assoCylinder (c, h, da);
-   return HOK;
-}
-// ====================================================== makeSphericalGrid
-int Elements::makeSphericalGrid (Vertex* c, Vector* dv, int nb, double  k)
-{
-   if (BadElement (c) || BadElement(dv) 
-       || nb<=0 || k < Epsil 
-       || dv->getDx()<=Epsil || dv->getDy()<=Epsil || dv->getDz()<=Epsil)
+   pfile    vtk = fopen (nomfic, "w");
+   if (vtk==NULL)
       {
-      setError ();
+      DumpReturn (HERR);
       return HERR;
       }
-
-   resize (GR_SPHERIC, nb);
-
-   Vertex* i_node [HV_MAXI];    // Les noeuds de l'hexa englobant
-   Edge*   i_edge [HE_MAXI];    // Les noeuds de l'hexa englobant
-   Quad*   i_quad [HQ_MAXI];    // Les noeuds de l'hexa englobant
-
-   for (int nro=0 ; nro<HV_MAXI; nro++)
-       {
-       double dx = glob->CoordVertex (nro, dir_x) * dv->getDx();
-       double dy = glob->CoordVertex (nro, dir_y) * dv->getDy();
-       double dz = glob->CoordVertex (nro, dir_z) * dv->getDz();
-
-       i_node [nro] = el_root->addVertex (c->getX ()+dx, c->getY ()+dy, 
-                                          c->getZ ()+dz);
-       }
-
-   for (int nro=0 ; nro<HE_MAXI; nro++)
-       {
-       int v1 = glob->EdgeVertex (nro, V_AMONT);
-       int v2 = glob->EdgeVertex (nro, V_AVAL);
-       i_edge[nro] = newEdge (i_node[v1], i_node[v2]);
-
-       if (db)
-          {
-          char nm0[8], nm1 [8], nm2 [8];
-          printf (" %2d : %s = %s = [%s, %s] = [%d,%d] = [%s,%s]\n", nro, 
-                 glob->namofHexaEdge(nro), i_edge[nro]->getName(nm0), 
-                 glob->namofHexaVertex(v1), glob->namofHexaVertex(v2), v1, v2,
-                 i_node[v1]->getName(nm1), i_node[v2]->getName(nm2));
-          }
-       }
-        
-   for (int nro=0 ; nro<HQ_MAXI; nro++)
-       i_quad[nro] = newQuad (i_edge[glob->QuadEdge (nro, E_A)], 
-                              i_edge[glob->QuadEdge (nro, E_B)], 
-                              i_edge[glob->QuadEdge (nro, E_C)], 
-                              i_edge[glob->QuadEdge (nro, E_D)]);
-
-   tab_hexa.push_back (newHexa (i_quad[Q_A], i_quad[Q_B], i_quad[Q_C], 
-                                i_quad[Q_D], i_quad[Q_E], i_quad[Q_F]));
-   double lambda  = 1;
-   double dcell   = 1;
-   for (int niv=0; niv<gr_rayon ; niv++)
-       {
-       double lambda0 = lambda;
-       dcell  *= k;
-       lambda += dcell;
-       addStrate (i_quad, i_edge, i_node, c,  lambda/lambda0);
-       }
-       
-   return HOK;
-}
-// ====================================================== addStrate
-int Elements::addStrate (Quad* i_quad[], Edge* i_edge[], Vertex* i_node[], 
-                        Vertex* center,  double lambda)
-{
-   Vertex* e_node [HV_MAXI];    // Les noeuds de l'hexa englobant
-   Edge*   e_edge [HE_MAXI];    // Les noeuds de l'hexa englobant
-   Quad*   e_quad [HQ_MAXI];    // Les noeuds de l'hexa englobant
-
-   Edge*   d_edge [HV_MAXI];    // Les aretes diagonales (1 par sommet)
-   Quad*   d_quad [HE_MAXI];    // Les faces  diagonales (1 par arete)
-
-                                          // Les sommets
-                                          //  + les aretes diagonales
-   for (int nv=0 ; nv<HV_MAXI ; nv++)
-       {
-       double px0 = center->getX ();
-       double py0 = center->getY ();
-       double pz0 = center->getZ ();
-       e_node[nv] = el_root->addVertex (px0+lambda*(i_node[nv]->getX()-px0),
-                                        py0+lambda*(i_node[nv]->getY()-py0),
-                                        pz0+lambda*(i_node[nv]->getZ()-pz0));
-
-       d_edge[nv] = newEdge (i_node[nv], e_node[nv]);
-       }
-                                          // Les aretes exterieures
-                                          //  + les faces diagonales
-   for (int nro=0 ; nro<HE_MAXI ; nro++)
-       {
-       int nv0  = glob->EdgeVertex (nro, V_AMONT);
-       int nv1  = glob->EdgeVertex (nro, V_AVAL );
-       e_edge[nro] = newEdge (e_node [nv0], e_node [nv1]);
-       d_quad[nro] = newQuad (i_edge [nro], d_edge [nv0], 
-                              e_edge [nro], d_edge [nv1]);
-       }
-                                          // Les faces exterieures
-                                          //  + les hexas
-   Hexa* strate = NULL;
-   for (int nro=0 ; nro<HQ_MAXI ; nro++)
-       {
-       int ne0 = glob->QuadEdge (nro, E_A);
-       int ne1 = glob->QuadEdge (nro, E_B);
-       int ne2 = glob->QuadEdge (nro, E_C);
-       int ne3 = glob->QuadEdge (nro, E_D);
-
-       e_quad[nro] = newQuad (e_edge[ne0], e_edge[ne1], 
-                              e_edge[ne2], e_edge[ne3]);
-       strate  = newHexa (i_quad[nro], e_quad[nro], d_quad[ne0], 
-                          d_quad[ne2], d_quad[ne1], d_quad[ne3]);
-       tab_hexa.push_back (strate);
-       }
-
-   for (int nv=0 ; nv<HV_MAXI ; nv++) i_node  [nv] = e_node [nv];
-   for (int ns=0 ; ns<HE_MAXI ; ns++) i_edge  [ns] = e_edge [ns];
-   for (int nq=0 ; nq<HQ_MAXI ; nq++) i_quad  [nq] = e_quad [nq];
-
-   return HOK;
-}
-// ====================================================== saveVtk
-int Elements::saveVtk (cpchar nomfic)
-{
                                            // -- 1) Raz numerotation precedente
    for (int nro=0 ; nro<nbr_hexas ; nro++)
        {
@@ -358,6 +209,13 @@ int Elements::saveVtk (cpchar nomfic)
        if (cell!=NULL && cell->isHere())
            cell->razNodes ();
        }
+   for (int nro=0 ; nro<size_hplus ; nro++)
+       {
+       Hexa* cell = ker_hexa[nro];
+       if (cell!=NULL && cell->isHere())
+           cell->razNodes ();
+       }
+
 
    int nbnodes = 0;
    int nbcells = 0;
@@ -372,7 +230,16 @@ int Elements::saveVtk (cpchar nomfic)
           }
        }
 
-   pfile    vtk = fopen (nomfic, "w");
+   for (int nro=0 ; nro<size_hplus ; nro++)
+       {
+       Hexa* cell = ker_hexa[nro];
+       if (cell!=NULL && cell->isHere())
+          {
+          nbcells ++;
+          nbnodes += cell->countNodes ();
+          }
+       }
+
    fprintf (vtk, "# vtk DataFile Version 3.1\n");
    fprintf (vtk, "%s \n", nomfic);
    fprintf (vtk, "ASCII\n");
@@ -384,6 +251,12 @@ int Elements::saveVtk (cpchar nomfic)
    for (int nro=0 ; nro<nbr_hexas ; nro++)
        {
        Hexa* cell = tab_hexa[nro];
+       if (cell!=NULL && cell->isHere())
+          cell->printNodes (vtk, nronode);
+       }
+   for (int nro=0 ; nro<size_hplus ; nro++)
+       {
+       Hexa* cell = ker_hexa[nro];
        if (cell!=NULL && cell->isHere())
           cell->printNodes (vtk, nronode);
        }
@@ -397,19 +270,26 @@ int Elements::saveVtk (cpchar nomfic)
        if (cell!=NULL && cell->isHere())
           cell->printHexa (vtk);
        }
+   for (int nro=0 ; nro<size_hplus ; nro++)
+       {
+       Hexa* cell = ker_hexa[nro];
+       if (cell!=NULL && cell->isHere())
+          cell->printHexa (vtk);
+       }
 
    fprintf (vtk, "CELL_TYPES %d\n", nbcells);
    for (int nro=0 ; nro<nbcells ; nro++)
        fprintf (vtk, "%d\n", HE_MAXI);
 
    fclose (vtk);
+   DumpReturn (HOK);
    return HOK;
 }
 // ============ = = = = = = = = = Vertices ..........
 // ====================================================== newEdge
 Edge* Elements::newEdge (Vertex* v1, Vertex* v2)
 {
-   if (v1==NULL || v2==NULL) 
+   if (v1==NULL || v2==NULL)
       return NULL;
 
    Edge* elt = new Edge (v1, v2);
@@ -425,7 +305,7 @@ Quad* Elements::newQuad (Edge* e1, Edge* e2, Edge* e3, Edge* e4)
    return elt;
 }
 // ====================================================== newHexa
-Hexa* Elements::newHexa (Quad* qa, Quad* qb, Quad* qc, Quad* qd, 
+Hexa* Elements::newHexa (Quad* qa, Quad* qb, Quad* qc, Quad* qd,
                                              Quad* qe, Quad* qf)
 {
    if (qa==NULL || qb==NULL || qc==NULL|| qd==NULL || qe==NULL|| qf==NULL)
@@ -434,87 +314,6 @@ Hexa* Elements::newHexa (Quad* qa, Quad* qb, Quad* qc, Quad* qd,
    Hexa*  elt = new Hexa (qa, qb, qc, qd, qe, qf);
    return elt;
 }
-// ====================================================== joinQuads
-int Elements::joinQuads (Quads& orig, int nb, Vertex* v1, Vertex* v2, 
-                  Vertex* v3, Vertex* v4, Quad* cible)
-{
-   resize (GR_JOINT, orig.size(), nb);
-
-   el_root->markAll (IS_NONE);
-   db = on_debug();
-   // db = el_root->debug ();
-   
-   gr_hauteur  = nb;
-   nbr_orig = orig.size();
-
-   for (int nro=0 ; nro<nbr_orig ; nro++)
-       {
-       Quad* face = orig[nro];
-       if (face==NULL)
-          {
-          printf ("\n");
-          printf (" *** joinQuads : donnees incorrectes\n");
-          printf (" *** le %deme quadrangle de depart est NULL\n", nro);
-          setError ();
-          return HERR;
-          }
-       else if (face->getNbrParents()>1)
-          {
-          printf ("\n");
-          printf (" *** joinQuads : donnees incorrectes\n");
-          printf (" *** le %deme quadrangle de depart n'est pas une " 
-                  "face externe\n", nro);
-          face->dump ();
-          setError ();
-          return HERR;
-          }
-       orig [nro]->setMark (nro);
-       tab_orig.push_back (orig[nro]);
-       }
-
-   Edge* e_orig = tab_orig[0] -> findEdge (v1, v3);
-   Edge* e_dest = cible       -> findEdge (v2, v4);
-
-   if (e_orig==NULL)
-      {
-      printf ("\n");
-      printf (" *** joinQuads : donnees incorrectes\n");
-      printf (" *** Les vertex v1 et v3 passes en argument ne definissent\n");
-      printf (" *** pas une arete du 1er quadrangle de depart\n");
-      printf ("\n");
-      HexDump (v1);
-      HexDump (v3);
-      HexDump (orig[0]);
-      }
-
-   if (e_dest==NULL)
-      {
-      printf ("\n");
-      printf (" *** joinQuads : donnees incorrectes\n");
-      printf (" *** Les vertex v2 et v4 passes en argument ne definissent\n");
-      printf (" *** pas une arete du quadrangle cible\n");
-      printf ("\n");
-      HexDump (v2);
-      HexDump (v4);
-      HexDump (cible);
-      }
-
-   if (e_orig==NULL || e_dest==NULL)
-      {
-      setError ();
-      return HERR;
-      }
-
-   StrOrient orient (v1, v3, v2, v4);
-   int ier =  this   ->coupler (0, cible, &orient);
-   if (ier!=HOK)
-      {
-      setError ();
-      return HERR;
-      }
-   ier = orig[0]->coupler (cible, &orient, this);
-   return ier;
-}
 // ======================================================== coupler
 int Elements::coupler (int nquad, Quad* dest, StrOrient* orient)
 {
@@ -527,20 +326,20 @@ int Elements::coupler (int nquad, Quad* dest, StrOrient* orient)
    int n12 = orig->indexVertex (orient->v12);
    int n21 = dest->indexVertex (orient->v21);
    int n22 = dest->indexVertex (orient->v22);
-   
+
                // ---------------- Les 4 sommets initiaux
-   Vertex* vorig[QUAD4] = { orient->v11, orient->v12, 
+   Vertex* vorig[QUAD4] = { orient->v11, orient->v12,
                             orig->getVertex((n11+2) MODULO QUAD4),
                             orig->getVertex((n12+2) MODULO QUAD4) };
 
-   Vertex* vdest[QUAD4] = { orient->v21, orient->v22, 
+   Vertex* vdest[QUAD4] = { orient->v21, orient->v22,
                             dest->getVertex((n21+2) MODULO QUAD4),
                             dest->getVertex((n22+2) MODULO QUAD4) };
    if (db)
       {
       printf ("Quad nro %d : ", nquad);
-      orig->printName (" est couple avec "); 
-      dest->printName ("\n"); 
+      orig->printName (" est couple avec ");
+      dest->printName ("\n");
       printf ("Orientation : (");
       for (int ii=0 ; ii<QUAD4 ; ii++) printf("%s ", vorig[ii]->getName());
       printf (")\n");
@@ -563,28 +362,28 @@ int Elements::coupler (int nquad, Quad* dest, StrOrient* orient)
           double py0 = nd1->getY();
           double pz0 = nd1->getZ();
 
-          double dx = (nd2->getX() - nd1->getX()) / gr_hauteur;
-          double dy = (nd2->getY() - nd1->getY()) / gr_hauteur;
-          double dz = (nd2->getZ() - nd1->getZ()) / gr_hauteur;
+          double dxt = nd2->getX() - px0;
+          double dyt = nd2->getY() - py0;
+          double dzt = nd2->getZ() - pz0;
 
           nd1->setMark (indVertex (nquad, ns, 0));
 
           Vertex* nd = nd1;
           for (int nh=0 ; nh<gr_hauteur ; nh++)
               {
-              int nh1 = nh + 1;
+              double coeff = under_v6 ? ((double)(nh+1))/gr_hauteur : cum_values [nh];
               Vertex* ndp = nd;
               if (nh == gr_hauteur-1)
                  nd = nd2 ;
-              else 
-                 nd = el_root->addVertex (px0 + nh1*dx, py0 + nh1*dy,  
-                                          pz0 + nh1*dz);
+              else
+                 nd = el_root->addVertex (px0 + coeff*dxt, py0 + coeff*dyt,
+                                          pz0 + coeff*dzt);
               int nv  = indVertex (nquad, ns, nh);
               tab_vertex [nv] = nd;
               tab_edge   [nv] = el_root->addEdge (ndp, nd);
               if (db)
-                 printf (" Edge vertical nro %d = %s = (%s, %s)\n", nv, 
-                         tab_edge[nv]->getName(), 
+                 printf (" Edge vertical nro %d = %s = (%s, %s)\n", nv,
+                         tab_edge[nv]->getName(),
                          ndp->getName(), nd->getName());
               }
           }
@@ -625,7 +424,7 @@ int Elements::coupler (int nquad, Quad* dest, StrOrient* orient)
               eb  = tab_edge [nvb];
               if (nh==gr_hauteur-1)
                  ea = dest->findEdge (vdest[ns], vdest[next]);
-              else 
+              else
                  ea = el_root->addEdge (tab_vertex [nva], tab_vertex [nvb]);
 
               propagateAssociation (ec, ea, eb);
@@ -660,7 +459,7 @@ int Elements::coupler (int nquad, Quad* dest, StrOrient* orient)
        if (nh == gr_hauteur-1)
           fb = dest;
        else
-          fb = newQuad (tab_edge [nroEdgeH (nquad, E_A, nh)], 
+          fb = newQuad (tab_edge [nroEdgeH (nquad, E_A, nh)],
                         tab_edge [nroEdgeH (nquad, E_B, nh)],
                         tab_edge [nroEdgeH (nquad, E_C, nh)],
                         tab_edge [nroEdgeH (nquad, E_D, nh)]);
@@ -676,99 +475,6 @@ int Elements::coupler (int nquad, Quad* dest, StrOrient* orient)
        }
    return HOK;
 }
-// ====================================================== makeCartesianNodes
-int Elements::makeCartesianNodes (Vertex* orig, Vector* v1, Vector* v2, 
-                   Vector* v3, int px, int py, int pz, int mx, int my, int mz)
-{
-   double dx  = v1->getDx() + v2->getDx() + v3->getDx();
-   double dy  = v1->getDy() + v2->getDy() + v3->getDy();
-   double dz  = v1->getDz() + v2->getDz() + v3->getDz();
-
-   double px0 = orig->getX () - mx * dx;
-   double py0 = orig->getY () - my * dy;
-   double pz0 = orig->getZ () - mz * dz;
-
-   int nbre= 0;
-
-   for (int nz=0 ; nz<size_vz ; nz++)
-       for (int ny=0 ; ny<size_vy ; ny++)
-           for (int nx=0 ; nx<size_vx ; nx++)
-               {
-               Vertex* node = orig;
-               if (nx!=mx || ny!=my || nz!=mz)
-                  node = el_root->addVertex (px0 + nx*dx,  py0 + ny*dy,  
-                                                          pz0 + nz*dz);
-               setVertex (node, nx, ny, nz);
-               nbre++;
-               }
-   return HOK;
-}
-// ====================================================== makeCylindricalNodes
-int Elements::makeCylindricalNodes (Vertex* orig, Vector* base, Vector* haut, 
-            double dr, double da, double dl, int nr, int na, int nl, bool fill)
-{
-   int ier = makeBasicCylinder (dr, da, dl, nr, na, nl, fill);
-   if (ier!=HOK) 
-       return ier;
-
-   transfoVertices  (orig,  base, haut);
-   return HOK;
-}
-// ====================================================== transfoVertices
-void Elements::transfoVertices (Vertex* orig, Vector* base, Vector* haut)
-{
-   Vector* iprim = new Vector (base);
-   Vector* jprim = new Vector (base);
-   Vector* kprim = new Vector (haut);
-
-   int ier = kprim->renormer ();
-   if (ier!=HOK) 
-       return;
-
-   jprim->vectoriel (kprim, base);
-   ier = jprim->renormer ();
-   if (ier!=HOK) 
-       return;
-
-   iprim->vectoriel (jprim, kprim);
-   transfoVertices  (orig,  iprim, jprim, kprim);
-}
-// ====================================================== transfoVertices
-void Elements::transfoVertices (Vertex* orig, Vector* iprim, Vector* jprim, 
-                                Vector* kprim)
-{
-   double matrice[DIM3][DIM3]={{iprim->getDx(),jprim->getDx(),kprim->getDx()},
-                               {iprim->getDy(),jprim->getDy(),kprim->getDy()},
-                               {iprim->getDz(),jprim->getDz(),kprim->getDz()}};
-
-   double matkx = orig->getX();
-   double matky = orig->getY();
-   double matkz = orig->getZ();
-
-   int nbre = tab_vertex.size ();
-   for (int nro=0 ; nro<nbre ; nro++)
-       {
-       if (tab_vertex[nro] != NULL) 
-           tab_vertex[nro]->setMark (NO_USED);
-       }
-
-   for (int nro=0 ; nro<nbre ; nro++)
-       {
-       Vertex* node =  tab_vertex[nro];
-       if (node != NULL  && node->getMark() == NO_USED)
-          {
-          double point [DIM3] = {node->getX(), node->getY(), node->getZ()};
-          double result[DIM3] = {matkx, matky, matkz};
-
-          for (int ni=0 ; ni<DIM3; ni++)
-              for (int nj=0 ; nj<DIM3; nj++)
-                  result [ni] += matrice[ni][nj] * point[nj];
-               
-          node->setCoord (result[dir_x], result[dir_y], result[dir_z]);
-          node->setMark (IS_USED);
-          }
-       }
-}
 // ====================================================== transform
 int Elements::transform (Matrix* matrice)
 {
@@ -821,7 +527,7 @@ int Elements::cutHexas  (const Edges& t_edges, int nbcuts)
        Edge* arete   = t_edges [nro];
        v_amont [nro] = arete->getAmont ();
        v_aval  [nro] = arete->getAval  ();
-       if (db) 
+       if (db)
           {
           printf (" %3d : Edge = (", nro);
           v_amont[nro]->printName (", ");
@@ -842,14 +548,14 @@ int Elements::cutHexas  (const Edges& t_edges, int nbcuts)
               for (int nh=0 ; nh<nbcubes ; nh++)
                   {
                   Hexa* hexa = quad->getParent (nh);
-                  if (hexa->getMark () != IS_USED)  
+                  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 (" %3d : Quad = ", nbfaces);
@@ -869,37 +575,37 @@ int Elements::cutHexas  (const Edges& t_edges, int nbcuts)
    nbr_vertex    = nbnodes*(nbcuts+2);
    int nbpiliers = nbnodes*(nbcuts+1);            // aretes verticales
    int nbpoutres = nbcells*(nbcuts+2)*QUAD4;      // aretes horizontales
-   nbr_edges     = nbpoutres; 
+   nbr_edges     = nbpoutres;
    nbr_quads     = nbcells*(nbcuts+1)*QUAD4;      // faces Verticales
    nbr_hexas     = nbcells*(nbcuts+1);
 
                    // ------------------- Les noeuds et les aretes verticales
-   tab_quad.resize   (nbr_quads);
-   tab_edge.resize   (nbr_edges);
-   tab_hexa.resize   (nbr_hexas);
-   tab_vertex.resize (nbr_vertex);
+   tab_quad.resize   (nbr_quads,  NULL);
+   tab_edge.resize   (nbr_edges,  NULL);
+   tab_hexa.resize   (nbr_hexas,  NULL);
+   tab_vertex.resize (nbr_vertex, NULL);
    vector <Edge*>    tab_pilier (nbpiliers);
 
    int nbinter = nbcuts + 1;
    for (int ned=0; ned<nbnodes ; ned++)
        {
-       Shapes ass_shapes = t_edges[ned] -> getAssociations ();
        Edges  ass_edges;
        t_edges [ned]->remove ();
        Vertex* ndamont = v_amont [ned];
        Vertex* ndaval  = v_aval  [ned];
 
-       double dx = (ndaval->getX() - ndamont->getX()) / nbinter;
-       double dy = (ndaval->getY() - ndamont->getY()) / nbinter;
-       double dz = (ndaval->getZ() - ndamont->getZ()) / nbinter;
+       double lgx = ndaval->getX() - ndamont->getX();
+       double lgy = ndaval->getY() - ndamont->getY();
+       double lgz = ndaval->getZ() - ndamont->getZ();
 
        Vertex* nd0 = tab_vertex [ned] = ndamont;
        for (int nc=0; nc<nbcuts ; nc++)
            {
            int nc1 = nc+1;
-           Vertex* nd1 = el_root->addVertex (ndamont->getX() + nc1*dx, 
-                                             ndamont->getY() + nc1*dy, 
-                                             ndamont->getZ() + nc1*dz);
+           double coeff = under_v6 ? ((double)(nc+1))/nbinter : cum_values[nc]; 
+           Vertex* nd1 = el_root->addVertex (ndamont->getX() + coeff*lgx,
+                                             ndamont->getY() + coeff*lgy,
+                                             ndamont->getZ() + coeff*lgz);
            tab_vertex [nc1*nbnodes + ned] = nd1;
            tab_pilier [nc *nbnodes + ned] = newEdge (nd0, nd1);
            ass_edges.push_back (tab_pilier [nc *nbnodes + ned]);
@@ -909,7 +615,8 @@ int Elements::cutHexas  (const Edges& t_edges, int nbcuts)
        tab_pilier [nbcuts *nbnodes + ned] = newEdge (nd0, ndaval);
        ass_edges.push_back (tab_pilier[nbcuts *nbnodes + ned]);
        ndamont->setMark (ned);
-       cutAssociation (ass_shapes, ass_edges);
+       if (t_edges[ned]->isAssociated())
+           cutAssociation (t_edges[ned], ass_edges);
        }
                    // ------------------- Les aretes horizontales
                    // ------------------- Les faces verticales
@@ -930,9 +637,9 @@ int Elements::cutHexas  (const Edges& t_edges, int nbcuts)
                   {
                   tab_edge [nc*sizelig + nmur] = tab_edge [nc*sizelig + nmur0];
                   tab_quad [nc*sizelig + nmur] = tab_quad [nc *sizelig + nmur0];
-                  if (db) 
+                  if (db)
                      {
-                     printf (" %2d : %d quad_vertical [%02d] =", nro, ns, 
+                     printf (" %2d : %d quad_vertical [%02d] =", nro, ns,
                                                         nc*sizelig + nmur);
                      printf (" quad_vertical [%02d]\n",  nc*sizelig + nmur0);
                      }
@@ -963,7 +670,7 @@ int Elements::cutHexas  (const Edges& t_edges, int nbcuts)
                      {
                      v1 = vis_a_vis [vs1];
                      v2 = vis_a_vis [vs2];
-                     ed2 = toit->findEdge (v1, v2); 
+                     ed2 = toit->findEdge (v1, v2);
                      }
 
                   tab_edge [nc1*sizelig + nmur] = ed2;
@@ -971,9 +678,9 @@ int Elements::cutHexas  (const Edges& t_edges, int nbcuts)
                             tab_pilier [nc*nbnodes + nd1], ed2,
                             tab_pilier [nc*nbnodes + nd2]);
                   ed0 = ed2;
-                  if (db) 
+                  if (db)
                      {
-                     printf (" %2d : %d quad_vertical [%02d] = ", nro, ns, 
+                     printf (" %2d : %d quad_vertical [%02d] = ", nro, ns,
                                                          nc*sizelig + nmur);
                      PrintName (tab_quad [nc *sizelig + nmur]);
                      printf ("\n");
@@ -984,7 +691,7 @@ int Elements::cutHexas  (const Edges& t_edges, int nbcuts)
        }
                    // ------------------- Les faces horizontales
                    // ------------------- Les hexas
-   // Rappel : sizelig = nbcells*QUAD4 
+   // Rappel : sizelig = nbcells*QUAD4
    for (int nro=0; nro<nbcells ; nro++)
        {
        Quad* qa = q_amont [nro];
@@ -995,7 +702,7 @@ int Elements::cutHexas  (const Edges& t_edges, int nbcuts)
            if (nc<nbcuts)
               {
               int  edh   = (nc+1)*nbcells*QUAD4 + nro*QUAD4;
-              qb = newQuad (tab_edge[edh],   tab_edge[edh+1], 
+              qb = newQuad (tab_edge[edh],   tab_edge[edh+1],
                             tab_edge[edh+2], tab_edge[edh+3]);
               if (BadElement(qb))
                   return HERR;
@@ -1018,46 +725,46 @@ void clear_associations (std::vector<Quad*>& table)
 {
    int nbelts = table.size();
    for (int nro=0 ; nro<nbelts ; nro++)
-       clear_association (table[nro]);
+       {
+       Quad* elt = table [nro];
+       if (elt != NULL && elt->isValid())
+           elt->clearAssociation ();
+       }
 }
 // ====================================================== clear_associations
 void clear_associations (std::vector<Edge*>& table)
 {
    int nbelts = table.size();
    for (int nro=0 ; nro<nbelts ; nro++)
-       clear_association (table[nro]);
+       {
+       Edge* elt = table [nro];
+       if (elt != NULL && elt->isValid())
+           elt->clearAssociation ();
+       }
 }
 // ====================================================== clear_associations
 void clear_associations (std::vector<Vertex*>& table)
 {
    int nbelts = table.size();
    for (int nro=0 ; nro<nbelts ; nro++)
-       clear_association (table[nro]);
+       {
+       Vertex* elt = table [nro];
+       if (elt != NULL && elt->isValid())
+           elt->clearAssociation ();
+       }
 }
 // ====================================================== clearAssociation
 void Elements::clearAssociation  ()
 {
-// clear_associations (tab_hexa);
    clear_associations (tab_quad);
    clear_associations (tab_edge);
    clear_associations (tab_vertex);
 
-   //  clear_associations (tab_orig);
 
    clear_associations (ker_hquad);
    clear_associations (ker_vquad);
    clear_associations (ker_hedge);
    clear_associations (ker_vedge);
-
-/* ***********************************************
-   nbelts = tab_hexa.size();
-   for (int nro=0 ; nro<nbelts ; nro++)
-       {
-       Hexa* elt = tab_hexa[nro];
-       if (elt != NULL && elt->isValid())
-           elt->clearAssociation ();
-       }
-   *********************************************** */
 }
 // ============================================================ findVertex
 int Elements::findVertex (double vx, double vy, double vz)
@@ -1080,5 +787,18 @@ int Elements::findVertex (double vx, double vy, double vz)
        }
    return NOTHING;
 }
+// ============================================================ findQuad
+Quad* Elements::findQuad (Edge* e1, Edge* e2)
+{
+   int nbre = tab_quad.size();
+   for (int nro=0 ; nro<nbre ; nro++)
+       {
+       Quad* quad = tab_quad [nro];
+       if (quad != NULL && quad->isHere ()
+                        && quad->definedBy (e1, e2))
+              return quad;
+       }
+   return NULL;
+}
 
 END_NAMESPACE_HEXA