X-Git-Url: http://git.salome-platform.org/gitweb/?a=blobdiff_plain;ds=sidebyside;f=src%2FHEXABLOCK%2FHexElements.cxx;h=087d0939827492799d2fb25563eaa410da9fc747;hb=0ebe741f3ee8a69f147a5c02b41f33d3f3f02e0f;hp=dafc9245cb894e20bef6d46527c9d0f9d7de84b0;hpb=3362101497b1bc0ded71b74c0806ac06c64d49d3;p=modules%2Fhexablock.git diff --git a/src/HEXABLOCK/HexElements.cxx b/src/HEXABLOCK/HexElements.cxx old mode 100755 new mode 100644 index dafc924..087d093 --- a/src/HEXABLOCK/HexElements.cxx +++ b/src/HEXABLOCK/HexElements.cxx @@ -1,24 +1,42 @@ -// C++ : Table des noeuds +// C++ : Grilles + +// 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, 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 "HexVector.hxx" #include "HexVertex.hxx" #include "HexHexa.hxx" #include "HexEdge.hxx" - #include "HexGlobale.hxx" #include - -static bool db=true; +#include 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 (); @@ -26,16 +44,23 @@ Elements::Elements (Document* doc) : EltBase (doc) size_qx = size_ex = size_vx = size_hx = 0; size_qy = size_ey = size_vy = size_hy = 0; size_qz = size_ez = size_vz = size_hz = 0; - size_qplus = size_eplus = size_vplus = size_hplus = 0; + size_qvplus = size_qhplus = size_ehplus = size_evplus = size_hplus = 0; nbr_hexas = nbr_quads = nbr_edges = nbr_vertex = 0; ker_vertex = 0; cyl_closed = false; cyl_fill = false; + grid_nocart = true; 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 (); @@ -43,7 +68,7 @@ Elements::Elements (Document* doc, int nx, int ny, int nz) : EltBase (doc) size_qx = size_ex = size_vx = size_hx = 0; size_qy = size_ey = size_vy = size_hy = 0; size_qz = size_ez = size_vz = size_hz = 0; - size_qplus = size_eplus = size_vplus = size_hplus = 0; + size_qvplus = size_qhplus = size_ehplus = size_evplus = size_hplus = 0; nbr_hexas = nbr_quads = nbr_edges = nbr_vertex = 0; cyl_closed = true; @@ -52,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) @@ -66,22 +93,25 @@ Elements::Elements (Elements* orig) : EltBase (orig->el_root) size_qx = size_ex = size_vx = size_hx = 0; size_qy = size_ey = size_vy = size_hy = 0; size_qz = size_ez = size_vz = size_hz = 0; - size_qplus = size_eplus = size_vplus = size_hplus = 0; + size_qvplus = size_qhplus = size_ehplus = size_evplus = size_hplus = 0; nbr_hexas = nbr_quads = nbr_edges = nbr_vertex = 0; 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) +void Elements::resize (EnumGrid type, int nx, int ny, int nz, int nplus) { - grid_type = type; + grid_type = type; + grid_nocart = true; switch (grid_type) { case GR_CARTESIAN : case GR_CYLINDRIC : default : + grid_nocart = false; size_hx = std::max (nx, 1); size_hy = std::max (ny, 1); size_hz = std::max (nz, 1); @@ -109,7 +139,7 @@ void Elements::resize (EnumGrid type, int nx, int ny, int nz) 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; @@ -122,6 +152,21 @@ void Elements::resize (EnumGrid type, int nx, int ny, int nz) nbr_edges = 2*nbr_vertex; break; + case GR_REPLACE : + nbr_orig = std::max (nx, 1); // nb quads du pattern + 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; + size_hy = 1; + size_hz = gr_hauteur; + + nbr_hexas = nbr_orig * gr_hauteur; + nbr_vertex = pat_nbvertex * (gr_hauteur+1); + nbr_edges = pat_nbedges * (gr_hauteur+1) + pat_nbvertex*gr_hauteur; + nbr_quads = nbr_orig * (gr_hauteur+1) + pat_nbedges *gr_hauteur; + break; + case GR_BICYL : cyl_closed = true; size_hx = 1; @@ -139,148 +184,24 @@ void Elements::resize (EnumGrid type, int nx, int ny, int nz) 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) -{ - resize (GR_CARTESIAN, px+mx, py+my, pz+mz); - - 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) -{ - resize (GR_CYLINDRIC, nr, na, nl); - cyl_closed = da >= 360.0; - makeCylindricalNodes (c, b, h, dr, da, dl, nr, na, nl, fill); - fillGrid (); - return HOK; -} -// ====================================================== makeSphericalGrid -int Elements::makeSphericalGrid (Vertex* c, Vector* dv, int nb, double k) -{ - resize (GR_SPHERIC, nb); - - if (nb<0) - return HERR; - else if (dv->getDx()<=ZEROR || dv->getDy()<=ZEROR || dv->getDz()<=ZEROR) - return HERR; - - 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 ; nroCoordVertex (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 ; nroEdgeVertex (nro, V_AMONT)], - i_node[glob->EdgeVertex (nro, V_AVAL) ]); - - for (int nro=0 ; nroQuadEdge (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; nivgetX (); - 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 ; nroEdgeVertex (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 ; nroQuadEdge (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 ; nvisHere()) cell->razNodes (); } + for (int nro=0 ; nroisHere()) + cell->razNodes (); + } + int nbnodes = 0; int nbcells = 0; @@ -302,7 +230,16 @@ int Elements::saveVtk (cpchar nomfic) } } - pfile vtk = fopen (nomfic, "w"); + for (int nro=0 ; nroisHere()) + { + nbcells ++; + nbnodes += cell->countNodes (); + } + } + fprintf (vtk, "# vtk DataFile Version 3.1\n"); fprintf (vtk, "%s \n", nomfic); fprintf (vtk, "ASCII\n"); @@ -314,6 +251,12 @@ int Elements::saveVtk (cpchar nomfic) for (int nro=0 ; nroisHere()) + cell->printNodes (vtk, nronode); + } + for (int nro=0 ; nroisHere()) cell->printNodes (vtk, nronode); } @@ -327,19 +270,26 @@ int Elements::saveVtk (cpchar nomfic) if (cell!=NULL && cell->isHere()) cell->printHexa (vtk); } + for (int nro=0 ; nroisHere()) + cell->printHexa (vtk); + } fprintf (vtk, "CELL_TYPES %d\n", nbcells); for (int nro=0 ; nromarkAll (IS_NONE); - - gr_hauteur = nb; - nbr_orig = orig.size(); - - for (int nro=0 ; nrogetNbrParents()>1) - { - printf ("\n"); - printf (" *** joinQuads : donnees incorrectes\n"); - printf (" *** le %deme quadrangle de depart n'est pas une " - "face externe\n", nro); - face->dump (); - 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); - HexDump (e_orig); - HexDump (e_dest); - - 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) - return HERR; - - StrOrient orient (v1, v3, v2, v4); - this ->coupler (0, cible, &orient); - orig[0]->coupler (cible, &orient, this); - return HOK; -} // ======================================================== coupler -void Elements::coupler (int nquad, Quad* dest, StrOrient* orient) +int Elements::coupler (int nquad, Quad* dest, StrOrient* orient) { Quad* orig = tab_orig [nquad]; @@ -447,26 +326,29 @@ void 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); - - printf ("Quad nro %d : ", nquad); - orig->printName (" est couple avec "); - dest->printName (", "); - printf ("Orientation=(%d,%d) (%d,%d)\n", n11, n12, n21, n22); + // ---------------- 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"); + printf ("Orientation : ("); + for (int ii=0 ; iigetName()); + printf (")\n"); + printf (" -> ("); + for (int ii=0 ; iigetName()); + printf (")\n"); + } // ---------------- Les sommets + les aretes verticales -#define Put(e) if (e) { printf ( " ... " #e " = ") ; e->printName("\n"); } \ - else printf ( " ... " #e " = 0x0\n") -#define Putn(e,n) { printf ( " ... " #e "[%d] = ",n) ; if(e[n]) e[n]->printName("\n"); \ - else printf ( "0x0\n") - for (int ns=0 ; ns< QUAD4 ; ns++) { Vertex* nd1 = vorig[ns]; @@ -480,25 +362,29 @@ void 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 ; nhaddVertex (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(), + ndp->getName(), nd->getName()); } } else @@ -538,10 +424,13 @@ void 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); tab_quad [nva] = newQuad (ea, eb, ec, ed); + if (BadElement (tab_quad [nva])) + return HERR; tab_edge [nha] = ea; } } @@ -570,7 +459,7 @@ void 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)]); @@ -581,96 +470,11 @@ void Elements::coupler (int nquad, Quad* dest, StrOrient* orient) ff = tab_quad [indVertex (nquad, E_D, nh)]; tab_hexa [nroHexa(nquad,nh)] = newHexa (fa,fb,fc,fd,fe,ff); + if (BadElement (tab_hexa [nroHexa(nquad,nh)])) + return HERR; } -} -// ====================================================== 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 ; nzaddVertex (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; - - Vector* iprim = new Vector (base); - Vector* jprim = new Vector (base); - Vector* kprim = new Vector (haut); - - ier = kprim->renormer (); - if (ier!=HOK) - return ier; - - jprim->vectoriel (kprim, base); - ier = jprim->renormer (); - if (ier!=HOK) - return ier; - - iprim->vectoriel (jprim, kprim); - transfoVertices (orig, iprim, jprim, kprim); - return HOK; -} -// ====================================================== 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 ; nrosetMark (NO_USED); - } - - for (int nro=0 ; nrogetMark() == NO_USED) - { - double point [DIM3] = {node->getX(), node->getY(), node->getZ()}; - double result[DIM3] = {matkx, matky, matkz}; - - for (int ni=0 ; nisetCoord (result[dir_x], result[dir_y], result[dir_z]); - node->setMark (IS_USED); - } - } -} // ====================================================== transform int Elements::transform (Matrix* matrice) { @@ -689,6 +493,18 @@ int Elements::transform (Matrix* matrice) cell->moveNodes (matrice); } + if (nbr_hexas!=0 ) + return HOK; + // -- Cas pathologique : il n'y avait pas d'hexas + // -- On se rabat sur les sommets + + for (int nro=0 ; nroisHere()) + matrice->perform (node); + } + return HOK; } // ====================================================== cutHexas @@ -699,17 +515,28 @@ int Elements::cutHexas (const Edges& t_edges, int nbcuts) // 2) Memo noeuds vector q_amont; vector q_aval; + map vis_a_vis; int nbnodes = t_edges.size(); vector v_amont (nbnodes); vector v_aval (nbnodes); + int nbfaces = 0; for (int nro=0; nrogetAmont (); v_aval [nro] = arete->getAval (); - int nbcells = arete->getNbrParents (); + if (db) + { + printf (" %3d : Edge = (", nro); + v_amont[nro]->printName (", "); + v_aval [nro]->printName (")\n"); + } + + vis_a_vis [v_amont[nro]] = v_aval[nro]; + vis_a_vis [v_aval[nro]] = v_amont[nro]; + int nbcells = arete->getNbrParents (); for (int nq=0 ; nqgetParent (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 (" Quad = "); + printf (" %3d : Quad = ", nbfaces); hexa->printName (", "); printf (" Faces = ("); hexa->getQuad (namont)->printName (", "); hexa->getQuad (naval )->printName (")\n"); + nbfaces ++; } } } @@ -747,54 +575,59 @@ 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 tab_pilier (nbpiliers); int nbinter = nbcuts + 1; for (int ned=0; nedremove (); 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; ncaddVertex (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]); nd0 = nd1; } tab_vertex [nbinter*nbnodes + ned] = ndaval; tab_pilier [nbcuts *nbnodes + ned] = newEdge (nd0, ndaval); + ass_edges.push_back (tab_pilier[nbcuts *nbnodes + ned]); ndamont->setMark (ned); + if (t_edges[ned]->isAssociated()) + cutAssociation (t_edges[ned], ass_edges); } // ------------------- Les aretes horizontales // ------------------- Les faces verticales + HexDisplay (nbcells); int sizelig = nbcells*QUAD4; for (int nro=0; nrogetEdge (ns); int nmur = nro*QUAD4 + ns; int nmur0 = plinthe->getMark(); @@ -804,46 +637,61 @@ 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]; - printf (" quad_vertical [%02d] = ", nc*sizelig + nmur); - printf (" quad_vertical [%02d]\n", nc*sizelig + nmur0); - printf ("\n"); + if (db) + { + printf (" %2d : %d quad_vertical [%02d] =", nro, ns, + nc*sizelig + nmur); + printf (" quad_vertical [%02d]\n", nc*sizelig + nmur0); + } } tab_edge [nbinter*sizelig+nmur] = tab_edge[nbinter*sizelig+nmur0]; } else { plinthe->setMark (nmur); - Vertex* v1 = sol->getVertex (ns); - Vertex* v2 = sol->getVertex ((ns+1) MODULO QUAD4); - int nd1 = v1->getMark (); - int nd2 = v2->getMark (); + Vertex* vs1 = sol->getVertex (ns); + Vertex* vs2 = sol->getVertex ((ns+1) MODULO QUAD4); + int nd1 = vs1->getMark (); + int nd2 = vs2->getMark (); Edge* ed0 = tab_edge [nmur] = plinthe; - Edge* edtoit = toit->getEdge (ns); + Edge* ed2 = NULL; + Vertex* v1 = NULL; + Vertex* v2 = NULL; for (int nc=0 ; ncfindEdge (v1, v2); + } + tab_edge [nc1*sizelig + nmur] = ed2; tab_quad [nc *sizelig + nmur] = newQuad (ed0, tab_pilier [nc*nbnodes + nd1], ed2, tab_pilier [nc*nbnodes + nd2]); ed0 = ed2; - printf (" quad_vertical [%02d] = ", nc*sizelig + nmur); - PrintName (tab_quad [nc *sizelig + nmur]); - printf ("\n"); + if (db) + { + printf (" %2d : %d quad_vertical [%02d] = ", nro, ns, + nc*sizelig + nmur); + PrintName (tab_quad [nc *sizelig + nmur]); + printf ("\n"); + } } } } } // ------------------- Les faces horizontales // ------------------- Les hexas - // Rappel : sizelig = nbcells*QUAD4 + // Rappel : sizelig = nbcells*QUAD4 for (int nro=0; nro& table) +{ + int nbelts = table.size(); + for (int nro=0 ; nroisValid()) + elt->clearAssociation (); + } +} +// ====================================================== clear_associations +void clear_associations (std::vector& table) +{ + int nbelts = table.size(); + for (int nro=0 ; nroisValid()) + elt->clearAssociation (); + } +} +// ====================================================== clear_associations +void clear_associations (std::vector& table) +{ + int nbelts = table.size(); + for (int nro=0 ; nroisValid()) + elt->clearAssociation (); + } +} +// ====================================================== clearAssociation +void Elements::clearAssociation () +{ + clear_associations (tab_quad); + clear_associations (tab_edge); + clear_associations (tab_vertex); + + + clear_associations (ker_hquad); + clear_associations (ker_vquad); + clear_associations (ker_hedge); + clear_associations (ker_vedge); +} +// ============================================================ findVertex +int Elements::findVertex (double vx, double vy, double vz) +{ + double tol = el_root->getTolerance (); + double xmin = vx - tol; + double xmax = vx + tol; + double ymin = vy - tol; + double ymax = vy + tol; + double zmin = vz - tol; + double zmax = vz + tol; + + int nbre = tab_vertex.size(); + for (int nro=0 ; nroisHere () + && node->isin (xmin, xmax, ymin, ymax, zmin, zmax)) + return nro; + } + return NOTHING; +} +// ============================================================ findQuad +Quad* Elements::findQuad (Edge* e1, Edge* e2) +{ + int nbre = tab_quad.size(); + for (int nro=0 ; nroisHere () + && quad->definedBy (e1, e2)) + return quad; + } + return NULL; +} + END_NAMESPACE_HEXA