Salome HOME
bos #18858 Use sphinx_rtd_theme as theme for SALOME documentation built with Sphinx
[modules/hexablock.git] / src / HEXABLOCK / HexElements_ter.cxx
old mode 100755 (executable)
new mode 100644 (file)
index cf1a35e..4d5a462
@@ -1,12 +1,11 @@
-
 // C++ : Table d'hexaedres (Evol Versions 3)
 
-// 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 "HexElements.hxx"
 
+#include "HexDocument.hxx"
 #include "HexVector.hxx"
 #include "HexVertex.hxx"
 #include "HexEdge.hxx"
 #include "HexHexa.hxx"
 #include "HexMatrix.hxx"
-#include "HexShape.hxx"
 #include "HexGlobale.hxx"
 
+#include "HexNewShape.hxx"
+#include "HexEdgeShape.hxx"
+#include "HexAssoEdge.hxx"
+
 #include <cmath>
 
 BEGIN_NAMESPACE_HEXA
 
-void geom_create_circle (double* milieu, double rayon, double* normale, 
-                         double* base, string& brep);
-void geom_create_sphere (double* milieu, double radius, string& brep);
-
-void translate_brep  (string& brep, double vdir[], string& trep);
-void transfo_brep  (string& brep, Matrix* matrice, string& trep);
-void geom_asso_point (Vertex* node);
-
 static bool db=false;
 
-// ====================================================== makeRind
-int Elements::makeRind (EnumGrid type, Vertex* center, Vector* vx, Vector* vz, 
-                        double radext, double radint, double radhole, 
-                       Vertex* plorig, double angle, int nr, int na, int nl)
-{
-   double phi1;
-   int ier = controlRind (type, center, vx, vz, radext, radint, radhole,
-                          plorig, angle, nr, na, nl, cyl_phi0, phi1);
-   if (ier!=HOK)
-      return ier;
-
-   resize (type, nr, na, nl);
-
-   cyl_radext  = radext;
-   cyl_radint  = radint;
-   cyl_radhole = radhole;
-   cyl_closed  = type==GR_HEMISPHERIC || type==GR_RIND;
-
-   double theta  = cyl_closed ? 2*M_PI : M_PI*angle/180;
-   cyl_dphi      = (phi1-cyl_phi0)/ size_hz;
-   cyl_dtheta    = theta / size_hy;
-
-   int nb_secteurs = cyl_closed ? size_vy-1 : size_vy;
-
-   for (int ny=0 ; ny<nb_secteurs ; ny++)
-       for (int nx=0 ; nx<size_vx ; nx++)
-           for (int nz=0 ; nz<size_vz ; nz++)
-               {
-               double px, py, pz;
-               getCylPoint (nx, ny, nz, px, py, pz);
-               Vertex* node = el_root->addVertex (px, py, pz);
-               setVertex (node, nx, ny, nz);
-               }
-   if (cyl_closed) 
-      for (int nx=0 ; nx<size_vx ; nx++)
-          for (int nz=0 ; nz<size_vz ; nz++)
-              {
-              Vertex* node = getVertexIJK (nx, 0, nz);
-              setVertex (node, nx, size_vy-1, nz);
-              }
-
-   transfoVertices (center, vx, vz);
-   fillGrid ();
-   assoCylinder (center, vz, angle);
-   return HOK;
-}
-// ====================================================== getCylPoint 
-int Elements::getCylPoint (int nr, int na, int nh, double& px, double& py,  
+// ====================================================== getCylPoint
+int Elements::getCylPoint (int nr, int na, int nh, double& px, double& py,
                            double& pz)
 {
    if (grid_type == GR_CYLINDRIC)
@@ -125,40 +74,40 @@ int Elements::getCylPoint (int nr, int na, int nh, double& px, double& py,
 
    return HOK;
 }
-// ====================================================== controlRind 
-int Elements::controlRind (EnumGrid type, Vertex* cx, Vector* vx, Vector* vz, 
+// ====================================================== controlRind
+int Elements::controlRind (EnumGrid type, Vertex* cx, Vector* vx, Vector* vz,
                            double rext, double rint, double rhole,
-                           Vertex* px, double angle, 
-                           int nrad, int nang, int nhaut, 
+                           Vertex* px, double angle,
+                           int nrad, int nang, int nhaut,
                            double &phi0, double &phi1)
 {
    const double Epsil1 = 1e-6;
    phi0  = phi1 = 0;
 
-   if (cx == NULL || vx == NULL || vz == NULL) 
+   if (cx == NULL || vx == NULL || vz == NULL)
       return HERR;
 
    if (nrad<=0 || nang<=0 || nhaut<=0)
       return HERR;
 
-   if (rext <= 0.0) 
+   if (rext <= 0.0)
       return HERR;
 
-   if (rint >= rext) 
+   if (rint >= rext)
       return HERR;
 
-   if (rhole > rint) 
+   if (rhole > rint)
       return HERR;
 
    double nvx = vx->getNorm();
    double nvz = vz->getNorm();
 
-   if (nvx < Epsil1 || nvz <  Epsil1) 
+   if (nvx < Epsil1 || nvz <  Epsil1)
       return HERR;
 
    double alpha = asin (rhole/rext);
    double beta  = -M_PI*DEMI;
-   if (type==GR_HEMISPHERIC || type==GR_PART_SPHERIC) 
+   if (type==GR_HEMISPHERIC || type==GR_PART_SPHERIC)
        alpha = 2*alpha;
 
    if (px!=NULL)
@@ -167,17 +116,17 @@ int Elements::controlRind (EnumGrid type, Vertex* cx, Vector* vx, Vector* vz,
       double oh = ((px->getX() - cx->getX()) * vz->getDx()
                 +  (px->getY() - cx->getY()) * vz->getDy()
                 +  (px->getZ() - cx->getZ()) * vz->getDz()) / nvz;
-      if (oh > rext) 
+      if (oh > rext)
          return HERR;
-      else if (oh > -rext) 
+      else if (oh > -rext)
          beta  = asin (oh/rext);
       }
 
    phi0 = std::max (alpha - M_PI*DEMI, beta);
    phi1 = M_PI*DEMI - alpha;
-   return HOK; 
+   return HOK;
 }
-// ====================================================== getHexas 
+// ====================================================== getHexas
 int Elements::getHexas (Hexas& liste)
 {
    liste.clear ();
@@ -189,139 +138,155 @@ int Elements::getHexas (Hexas& liste)
        }
    return HOK;
 }
-// ====================================================== assoCylinder 
+// ====================================================== assoCylinder
 void Elements::assoCylinder (Vertex* ori, Vector* normal, double angle)
 {
+   if (normal==NULL || ori==NULL)
+       return;
+
    RealVector tangles;
-   assoCylinders (ori, normal, angle, tangles);
+   Real3      vz, center;
+   normal->getCoord (vz);
+   ori   ->getPoint (center);
+   assoCylinders (center, vz, angle, tangles);
+}
+// ====================================================== assoCylinders
+void Elements::assoCylinders (Vertex* ori, Vector* normal, double angle,
+                              RealVector& t_angles)
+{
+   if (normal==NULL || ori==NULL)
+       return;
+
+   Real3 vz, center;
+   normal->getCoord (vz);
+   ori   ->getPoint (center);
+   assoCylinders (center, vz, angle, t_angles);
 }
-// ====================================================== assoCylinders 
-void Elements::assoCylinders (Vertex* ori, Vector* normal, double angle, 
+// ====================================================== assoCylinders
+void Elements::assoCylinders (double* ori, double* vk, double angle,
                               RealVector& t_angles)
 {
-   int    na      = t_angles.size();
-   bool   regul   = na == 0;
-   double alpha   = angle/size_hy;
+   char     name [12];
+   sprintf (name, "grid_%02d", el_id);
+   NewShape* geom = el_root->addShape (name, SH_CYLINDER);
+   geom -> openShape();
 
-   string brep;
-   Real3 vk = { normal->getDx(), normal->getDy(), normal->getDz() };
-   normer_vecteur (vk);
+   std::string brep;
+   // Real3 vk = { normal->getDx(), normal->getDy(), normal->getDz() };
+   // normer_vecteur (vk);
 
    for (int nz=0 ; nz<size_vz ; nz++)
        {
        for (int nx=0 ; nx<size_vx ; nx++)
            {
-           Vertex* pm = getVertexIJK (nx, 0, nz); 
-           Real3   om = { pm->getX() - ori->getX()
-                          pm->getY() - ori->getY()
-                          pm->getZ() - ori->getZ() };
+           Vertex* pm = getVertexIJK (nx, 0, nz);
+           Real3   om = { pm->getX() - ori[dir_x]
+                          pm->getY() - ori[dir_y]
+                          pm->getZ() - ori[dir_z] };
 
            double oh     = prod_scalaire (om, vk);
            double rayon  = 0;
            Real3  ph, hm;
            for (int dd=dir_x; dd<=dir_z ; dd++)
                {
-               ph [dd] = ori->getCoord(dd) + oh*vk[dd]; 
+               ph [dd] = ori[dd] + oh*vk[dd];
                hm [dd] = pm ->getCoord(dd) - ph[dd];
                rayon  += hm[dd] * hm[dd];
                }
 
            rayon = sqrt (rayon);
-           geom_create_circle (ph, rayon, vk, hm, brep);
+           int subid = geom->addCircle (ph, rayon, vk, hm);
 
-           double  pmax = 0;
            for (int ny=0 ; ny<size_hy ; ny++)
                {
-               double beta  = regul ? alpha : t_angles [ny]; 
-               double pmin = pmax;
-               pmax  += beta/360;
-               Edge* edge   =  getEdgeJ  (nx, ny, nz);
-               Shape* shape = new Shape (brep);
-               shape-> setBounds (pmin, pmax);
-               edge->addAssociation (shape);
+               double pmin = t_angles [ny]/360;
+               double pmax = t_angles [ny+1]/360;
+               Edge*  edge = getEdgeJ (nx, ny, nz);
+               geom->addAssociation (edge, subid, pmin, pmax);
+               if (db) std::cout << " assoCylinders : ny= " << ny 
+                                << ", nz= " << nz << " param = (" 
+                                << pmin << ", " << pmax  << ")\n";
                }
            }
        }
-   
    // Association automatique des vertex non associes -> bph
    // Traitement des faces spheriques
 
    Real3 vi = { -vk[dir_x],  -vk[dir_y],  -vk[dir_z] };
-   Real3 po = { ori->getX(), ori->getY(), ori->getZ() };
 
    switch (grid_type)
       {
       case GR_HEMISPHERIC  :    // Pour l'exterieur
       case GR_PART_SPHERIC :
-           assoRind (po, vi, size_vx-1);
+           assoRind (ori, vi, size_vx-1, geom);
            break;
       case GR_PART_RIND    :    // Exterieur + interieur
       case GR_RIND         :
-           assoRind (po, vi, 0);
-           assoRind (po, vi, size_vx-1);
+           assoRind (ori, vi, 0, geom);
+           assoRind (ori, vi, size_vx-1, geom);
            break;
       default :
            break;
       }
-   assoResiduelle ();  // Association des sommets residuels
+   geom->closeShape ();
 }
-// ====================================================== assoRind 
+// ====================================================== calcul_param
+double calcul_param (double* ori, double* vi, Vertex* vert)
+{
+   Real3  pnt, vj;
+   vert->getPoint (pnt);
+   calc_vecteur (ori, pnt, vj);
+   normer_vecteur  (vj);
+   double kos = prod_scalaire (vi, vj);
+   double alpha = acos (kos) / (2*M_PI);
+   return alpha;
+}
+// ====================================================== assoRind
 // Association des meridiennes
 // Creation sphere geometrique + association faces
-void Elements::assoRind (double* ori, double* vi, int nx)
+void Elements::assoRind (double* ori, double* vi, int nx, NewShape* geom)
 {
    Real3  vk;
-   Edges  contour;
-   string c_brep, s_brep;
-   Shape  shape (c_brep);
-   Shapes t_shape;
-   t_shape.push_back (&shape);
-
-   double radius = nx==0 ? cyl_radint : cyl_radext;
-   double paramin = (cyl_phi0 + M_PI/2) / (2*M_PI);
-   double paramax = paramin + size_hz*cyl_dphi/(2*M_PI);
-
-   paramin = std::max (paramin, 0.0);
-   paramax = std::min (paramax, 1.0);
-   geom_create_sphere (ori, radius, s_brep);
+   double radius  = nx==0 ? cyl_radint : cyl_radext;
+   int    sphid   = geom->addSphere (ori, radius);
 
    int nz1 = size_vz/2;
    int nb_secteurs = cyl_closed ? size_vy-1 : size_vy;
 
    for (int ny=0 ; ny<nb_secteurs ; ny++)
        {
-       Vertex* pm = getVertexIJK (nx, ny, nz1); 
+       Vertex* pm = getVertexIJK (nx, ny, nz1);
        Real3   vj = { pm->getX(), pm->getY(), pm->getZ() };
        prod_vectoriel (vi, vj, vk);
        double rayon = cyl_radint + nx*(cyl_radext-cyl_radint)/size_hx;
-       geom_create_circle (ori, rayon, vk, vi, c_brep);
-       shape.setBrep   (c_brep);
-       shape.setBounds (paramin, paramax);
-       contour.clear ();
+       int    subid = geom->addCircle (ori, rayon, vk, vi);
 
        for (int nz=0 ; nz<size_hz ; nz++)
            {
-           contour.push_back (getEdgeK  (nx, ny, nz));
-           Quad* quad = getQuadJK (nx, ny, nz);
-           if (quad != NULL)
-              {
-              Shape* sshape = new Shape (s_brep);
-              quad->addAssociation (sshape);
-              }
+           Quad* quad  = getQuadJK (nx, ny, nz);
+           Edge* edge  = getEdgeK  (nx, ny, nz);
+           Vertex* nd1 = edge->getVertex (V_AMONT);
+           Vertex* nd2 = edge->getVertex (V_AVAL);
+           double pmin = calcul_param (ori, vi, nd1);
+           double pmax = calcul_param (ori, vi, nd2);
+           std::cout << " assoRind : ny= " << ny << ", nz= " << nz 
+                    << " param = (" << pmin << ", " << pmax  << ")\n";
+
+           geom->addAssociation (edge, subid, pmin, pmax);
+           geom->addAssociation (quad, sphid);
            }
-       cutAssociation (t_shape, contour, false);
        }
 }
-// ====================================================== assoCircle 
+// ====================================================== assoCircle
 // ==== utilise pour les spheres carrees
-void Elements::assoCircle (double* center, Edge* ed1, Edge* ed2)
+void Elements::assoCircle (double* center, Edge* ed1, Edge* ed2, NewShape* geom)
 {
    Real3 oa,  ob, normal;
    Real3 pta, ptb, ptc, ptd;
-   string brep;
+   std::string brep;
 
 //  Les 2 edges dont les petits cotes d'un rectangle de rapport L/l=sqrt(2)
-//  Soit le cercle circonscrit a ce rectangle. 
+//  Soit le cercle circonscrit a ce rectangle.
 //    * la largeur est balayee par l'angle alpha
 //    * la longueur par l'angle beta = pi -alpha
 
@@ -342,112 +307,48 @@ void Elements::assoCircle (double* center, Edge* ed1, Edge* ed2)
    calc_vecteur   (center, pta, oa);
    calc_vecteur   (center, ptb, ob);
    prod_vectoriel (oa, ob, normal);
-   double rayon = calc_norme (oa);
-
-   geom_create_circle (center, rayon, normal, oa, brep);
 
-   Shape* asso1 = new Shape (brep);
-   Shape* asso2 = new Shape (brep);
+   double rayon = calc_norme (oa);
+   int    subid = geom->addCircle (center, rayon, normal, oa);
 
-   const double alpha = atan (sqrt(2)/2)/M_PI; //  angle proche de 70.528 degres
-   asso1->setBounds (0,   alpha);
-   asso2->setBounds (0.5, alpha + 0.5);
+   const double alpha = atan (sqrt(2.)/2)/M_PI; //  angle proche de 70.528 degres
+   // asso1->setBounds (0,   alpha);
+   // asso2->setBounds (0.5, alpha + 0.5);
 
-   ed1->addAssociation (asso1);
-   ed2->addAssociation (asso2);
+   geom->addAssociation (ed1, subid, 0.0, alpha);
+   geom->addAssociation (ed2, subid, 0.5, alpha+0.5);
 }
-// ====================================================== assoSphere 
-void Elements::assoSphere (Vertex* ori, Edge* t_edge[], Quad* t_quad[])
+// ====================================================== assoSphere
+void Elements::assoSphere (double* center, Edge* t_edge[], Quad* t_quad[])
 {
-   Real3 center, sommet;
-   ori->getPoint(center);
+   char     name [12];
+   sprintf (name, "grid_%02d", el_id);
+   NewShape* geom = el_root->addShape (name, SH_SPHERE);
+   geom -> openShape ();
 
-   assoCircle (center, t_edge [E_AC], t_edge [E_BD]);
-   assoCircle (center, t_edge [E_AD], t_edge [E_BC]);
-   assoCircle (center, t_edge [E_AE], t_edge [E_BF]);
-   assoCircle (center, t_edge [E_AF], t_edge [E_BE]);
-   assoCircle (center, t_edge [E_CE], t_edge [E_DF]);
-   assoCircle (center, t_edge [E_CF], t_edge [E_DE]);
+   Real3 sommet;
+
+   assoCircle (center, t_edge [E_AC], t_edge [E_BD], geom);
+   assoCircle (center, t_edge [E_AD], t_edge [E_BC], geom);
+   assoCircle (center, t_edge [E_AE], t_edge [E_BF], geom);
+   assoCircle (center, t_edge [E_AF], t_edge [E_BE], geom);
+   assoCircle (center, t_edge [E_CE], t_edge [E_DF], geom);
+   assoCircle (center, t_edge [E_CF], t_edge [E_DE], geom);
 
    t_edge[E_AC]->getVertex(V_AMONT)->getPoint (sommet);
-   double radius = calc_distance (center, sommet);;
+   double radius = calc_distance (center, sommet);
 
-   string brep;
-   geom_create_sphere (center, radius, brep);
+   int sphid = geom -> addSphere (center, radius);
+   geom->closeShape();
 
    for (int nf=0 ; nf < HQ_MAXI ; nf++)
-       {
-       Shape* shape = new Shape (brep);
-       t_quad [nf]->addAssociation (shape);
-       }
-
-   assoResiduelle ();  // Association des sommets residuels
-}
-// ====================================================== makeSphericalGrid
-int Elements::makeSphericalGrid (Vertex* c, double rayon, int nb, double  k)
-{
-   if (nb<=0 || rayon <=Epsil || k <= Epsil || BadElement (c))
-      {
-      setError ();
-      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) * rayon;
-       double dy = glob->CoordVertex (nro, dir_y) * rayon;
-       double dz = glob->CoordVertex (nro, dir_z) * rayon;
-
-       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);
-       }
-       
-   assoSphere (c, i_edge, i_quad);
-   return HOK;
+       t_quad[nf]->addAssociation (geom, sphid);
 }
 // ==================================================== propagateAssociation
 int Elements::propagateAssociation (Edge* orig, Edge* dest, Edge* dir)
 {
+    return HERR;
+#if 0
    if (revo_lution || orig==NULL || dest==NULL || dir==NULL)
       return HERR;
 
@@ -465,7 +366,7 @@ int Elements::propagateAssociation (Edge* orig, Edge* dest, Edge* dir)
    if (vo1==NULL || vd1==NULL)
       return HERR;
 
-   string  trep;
+   std::string  trep;
    Real3   pa, pb, vdir1, vdir2;
    calc_vecteur (vo1->getPoint (pa), vd1->getPoint (pb), vdir1);
    calc_vecteur (vo2->getPoint (pa), vd2->getPoint (pb), vdir2);
@@ -480,11 +381,11 @@ int Elements::propagateAssociation (Edge* orig, Edge* dest, Edge* dir)
           Shape* shape  = tab_shapes[nro];
           if (shape!=NULL)
              {
-             string brep   = shape->getBrep();
+             std::string brep   = shape->getBrep();
              translate_brep (brep, vdir1, trep);
-             Shape* tshape = new Shape (trep);
-             tshape->setBounds (shape->debut, shape->fin);
-             dest->addAssociation (tshape);
+             // Shape* tshape = new Shape (trep);
+             // tshape->setBounds (shape->getStart(), shape->getEnd());
+             // dest->addAssociation (tshape);
              }
           }
       }
@@ -495,85 +396,60 @@ int Elements::propagateAssociation (Edge* orig, Edge* dest, Edge* dir)
        Shape* shape = vo1->getAssociation ();
        if (shape!=NULL && vd1->getAssociation ()==NULL)
           {
-          string brep   = shape->getBrep();
+          std::string brep   = shape->getBrep();
           translate_brep (brep, vdir, trep);
-          Shape* tshape = new Shape (trep);
-          vd1->setAssociation (tshape);
+          // Shape* tshape = new Shape (trep);
+          // vd1->setAssociation (tshape);
           }
        vo1  = vo2;
        vd1  = vd2;
        vdir = vdir2;
        }
-
    return HOK;
+#endif
 }
 // ==================================================== prismAssociation
-int Elements::prismAssociation (Edge* orig, Edge* dest, int nh, Edge* dir)
+int Elements::prismAssociation (Edge* lorig, Edge* ldest, int nh)
 {
-   if (orig==NULL || dest==NULL || dir==NULL)
+   if (lorig==NULL || ldest==NULL)
       return HERR;
 
    updateMatrix (nh);
-   const Shapes& tab_shapes = orig->getAssociations ();
-   const Shapes& tab_dest   = dest->getAssociations ();
-   int   nbdest             = tab_dest.size();
-   int   nbshapes           = tab_shapes.size();
-   bool  on_edge            = nbshapes>0 && nbdest==0;
-
-   Vertex* vo1 = orig->commonVertex (dir);
-   Vertex* vd1 = dest->commonVertex (dir);
-
-   if (vo1==NULL || vd1==NULL)
-      return HERR;
-
-   string  trep;
-   Real3   porig, pdest, vdir;
-   vo1->getPoint (porig);
-   vd1->getPoint (pdest);
-   calc_vecteur  (porig, pdest, vdir);
-    
-   if (on_edge)
-      {
-      for (int nro=0 ; nro<nbshapes ; nro++)
-          {
-          Shape* shape  = tab_shapes[nro];
-          if (shape!=NULL)
-             {
-             string brep   = shape->getBrep();
-             //   translate_brep (brep, vdir, trep);
-             transfo_brep (brep, &gen_matrix, trep);
-             Shape* tshape = new Shape (trep);
-             tshape->setBounds (shape->debut, shape->fin);
-             dest->addAssociation (tshape);
-             }
-          }
-      }
+   int nbshapes = lorig->countAssociation ();
+   if (nbshapes==0)
+       return HOK;
 
-   for (int nro=V_AMONT ; nro<=V_AVAL ; nro++)
+   NewShape* new_shape = getShape ();
+   for (int nro=0 ; nro<nbshapes ; nro++)
        {
-       Shape* shape = vo1->getAssociation ();
-       if (shape!=NULL && vd1->getAssociation ()==NULL)
-          {
-          string brep   = shape->getBrep();
-          //  translate_brep (brep, vdir, trep);
-          transfo_brep (brep, &gen_matrix, trep);
-          Shape* tshape = new Shape (trep);
-          vd1->setAssociation (tshape);
-          }
-       vo1 = orig->opposedVertex (vo1);
-       vd1 = dest->opposedVertex (vd1);
+       AssoEdge*  asso  = lorig->getAssociation (nro);
+       EdgeShape* shape = asso ->getEdgeShape();
+       double     p1    = asso ->getStart();
+       double     p2    = asso ->getEnd  ();
+       int        subid = 0;
+       char       name [24];
+
+       sprintf (name, "0x%lx#%d", (unsigned long) shape, nh);
+       std::map<std::string,int>::iterator iter = map_shape.find (name);
+       if (iter != map_shape.end())
+          subid = iter->second;
+       else
+          subid = map_shape[name] = new_shape->transfoShape (cum_matrix, shape);
+
+       new_shape->addAssociation (ldest, subid, p1, p2);
+       printf (" prismAsso : %s -> %s(%d)  = %s [ %g , %g]\n",
+                 lorig->getName(), ldest->getName(), nh, name, p1, p2);
        }
    return HOK;
 }
 // ====================================================== assoResiduelle
 void Elements::assoResiduelle ()
 {
-   int nbre = tab_vertex.size();
-   for (int nv=0 ; nv<nbre ; nv++)
-       {
-       geom_asso_point (tab_vertex [nv]);
-       }
+   // int nbre = tab_vertex.size();
+   // for (int nv=0 ; nv<nbre ; nv++)
+       // {
+       // geom_asso_point (tab_vertex [nv]);
+       // }
 }
 // ====================================================== moveDisco
 void Elements::moveDisco (Hexa* hexa)
@@ -589,4 +465,206 @@ void Elements::moveDisco (Hexa* hexa)
        matrix.perform (tab_vertex [nv]);
        }
 }
+// =================================================== getStrate
+Hexa* Elements::getStrate (int couche, EnumHQuad nroface)
+{
+   Hexa* cell = NULL;
+   int   nro  = couche <= 0 ? 0 : (couche-1)*HQ_MAXI + nroface + 1;
+
+   if (nbr_hexas==0 || nro >= nbr_hexas)
+      cell = NULL;
+   else
+      cell = tab_hexa [nro];
+
+   return cell;
+}
+// ============================================================  setHexa
+void Elements::setHexa (Hexa* elt, int nro)
+{
+   if (nro >=0 && nro < nbr_hexas)
+       tab_hexa [nro] = elt;
+}
+// ============================================================  setQuad
+void Elements::setQuad (Quad* elt, int nro)
+{
+   if (nro >=0 && nro < nbr_quads)
+       tab_quad [nro] = elt;
+}
+// ============================================================  setEdge
+void Elements::setEdge (Edge* elt, int nro)
+{
+   if (nro >=0 && nro < nbr_edges)
+       tab_edge [nro] = elt;
+}
+// ============================================================  setVertex
+void Elements::setVertex (Vertex* elt, int nro)
+{
+   if (nro >=0 && nro < nbr_vertex)
+       tab_vertex [nro] = elt;
+}
+// -----------------------------------------------------------------------
+// -----------------------------------------------------------------------
+// ============================================================  getHexa
+Hexa* Elements::getHexa (int nro)
+{
+   DumpStart ("getHexa", nro);
+
+   Hexa* elt = NULL;
+   int  nombre=tab_hexa.size();
+   // if (nro >=0 && nro < nbr_hexas && el_status == HOK Abu 2010/05/06
+   if (nro >=0 && nro < nombre && el_status == HOK
+               && tab_hexa [nro] != NULL && tab_hexa [nro]->isValid())
+      elt = tab_hexa [nro];
+
+   DumpReturn (elt);
+   return elt;
+}
+// ============================================================  getQuad
+Quad* Elements::getQuad (int nro)
+{
+   DumpStart ("getQuad", nro);
+
+   Quad* elt = NULL;
+   if (nro >=0 && nro < nbr_quads && el_status == HOK
+               && tab_quad [nro] != NULL && tab_quad [nro]->isValid())
+      elt = tab_quad [nro];
+
+   DumpReturn (elt);
+   return elt;
+}
+// ============================================================  getEdge
+Edge* Elements::getEdge (int nro)
+{
+   DumpStart ("getEdge", nro);
+
+   Edge* elt = NULL;
+   if (nro >=0 && nro < nbr_edges && el_status == HOK
+               && tab_edge [nro] != NULL && tab_edge [nro]->isValid())
+      elt = tab_edge [nro];
+
+   DumpReturn (elt);
+   return elt;
+}
+// ============================================================  getVertex
+Vertex* Elements::getVertex (int nro)
+{
+   DumpStart ("getVertex", nro);
+
+   Vertex* elt = NULL;
+   if (nro >=0 && nro <  nbr_vertex && el_status == HOK
+               && tab_vertex [nro] != NULL && tab_vertex [nro]->isValid())
+      elt = tab_vertex [nro];
+
+   DumpReturn (elt);
+   return elt;
+}
+// ============================================================  indVertex
+int Elements::indVertex (int nquad, int nsommet, int nh)
+{
+   int nro = nsommet  + QUAD4*nquad + nbr_orig*QUAD4*(nh+1);
+   return nro;
+}
+// ============================================================  nroVertex
+int Elements::nroVertex (int nquad, int nsommet, int nh)
+{
+   int nro = nsommet  + QUAD4*nquad + nbr_orig*QUAD4*nh;
+   return nro;
+}
+// ============================================================  indVertex
+int Elements::indVertex (int ref_edge, int nh)
+{
+   int    nro = ref_edge + nbr_orig*QUAD4*nh;
+   return nro;
+}
+// ============================================================  nroEdgeH
+int Elements::nroEdgeH (int nvertex)
+{
+   return QUAD4*nbr_orig*gr_hauteur + nvertex;
+}
+// ============================================================  nroEdgeH
+int Elements::nroEdgeH (int nquad, int nsommet, int nh)
+{
+   return QUAD4*nbr_orig*gr_hauteur + indVertex (nquad, nsommet, nh);
+}
+// ============================================================  nroHexa
+int Elements::nroHexa (int nquad, int nh)
+{
+   int nro = gr_hauteur*nquad + nh;
+   return nro;
+}
+
+// ============================================================  addHexa
+void Elements::addHexa (Hexa* element)
+{
+   tab_hexa.push_back (element);
+   nbr_hexas ++;
+}
+// ============================================================  addQuad
+void Elements::addQuad (Quad* element)
+{
+   tab_quad.push_back (element);
+   nbr_quads ++;
+}
+// ============================================================  addEdge
+void Elements::addEdge (Edge* element)
+{
+   tab_edge.push_back (element);
+   nbr_edges ++;
+}
+// ============================================================  addVertex
+void Elements::addVertex (Vertex* element)
+{
+   tab_vertex.push_back (element);
+   nbr_vertex ++;
+}
+// ============================================================  findHexa
+int Elements::findHexa (Hexa* element)
+{
+   int nbre = tab_hexa.size();
+   for (int nro=0 ; nro<nbre ; nro++)
+       if (tab_hexa [nro] == element)
+          return nro;
+   return NOTHING;
+}
+// ============================================================  findQuad
+int Elements::findQuad (Quad* element)
+{
+   int nbre = tab_quad.size();
+   for (int nro=0 ; nro<nbre ; nro++)
+       if (tab_quad [nro] == element)
+          return nro;
+   return NOTHING;
+}
+// ============================================================  findEdge
+int Elements::findEdge (Edge* element)
+{
+   int nbre = tab_edge.size();
+   for (int nro=0 ; nro<nbre ; nro++)
+       if (tab_edge [nro] == element)
+          return nro;
+   return NOTHING;
+}
+// ============================================================  findVertex
+int Elements::findVertex (Vertex* element)
+{
+   int nbre = tab_vertex.size();
+   for (int nro=0 ; nro<nbre ; nro++)
+       if (tab_vertex [nro] == element)
+          return nro;
+   return NOTHING;
+}
+// ========================================================= saveVtk (avec nro)
+int Elements::saveVtk  (cpchar radical, int &nro)
+{
+   char num[8];
+   sprintf (num, "%d", nro);
+   nro ++;
+
+   std::string filename = radical;
+   filename += num;
+   filename += ".vtk";
+   int ier = saveVtk (filename.c_str());
+   return ier;
+}
+
 END_NAMESPACE_HEXA