X-Git-Url: http://git.salome-platform.org/gitweb/?a=blobdiff_plain;f=src%2FHEXABLOCK%2FHexElements.cxx;h=51ffae5e1d5ba4a6081ec5169db6d7acc72dade8;hb=5f7c68a2d0182cb6adf02b2b17032febac539364;hp=34479879c3845fcf1669b85ef44964c4315387c8;hpb=6924a056f811baefa30f31083b93b10f7dae3a35;p=modules%2Fhexablock.git diff --git a/src/HEXABLOCK/HexElements.cxx b/src/HEXABLOCK/HexElements.cxx index 3447987..51ffae5 100755 --- a/src/HEXABLOCK/HexElements.cxx +++ b/src/HEXABLOCK/HexElements.cxx @@ -1,12 +1,12 @@ // C++ : Grilles -// Copyright (C) 2009-2013 CEA/DEN, EDF R&D +// Copyright (C) 2009-2016 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 @@ -31,10 +31,10 @@ #include #include -static bool db=false; - BEGIN_NAMESPACE_HEXA +static bool db=on_debug(); + // ====================================================== Constructeur Elements::Elements (Document* doc) : EltBase (doc, EL_GRID) { @@ -54,6 +54,9 @@ Elements::Elements (Document* doc) : EltBase (doc, EL_GRID) 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) @@ -74,6 +77,8 @@ Elements::Elements (Document* doc, int nx, int ny, int nz) 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) @@ -93,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) @@ -178,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; - } + DumpStart ("saveVtk", nomfic); - 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) -{ - 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 ; 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); - 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 ; 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; @@ -373,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"); @@ -385,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); } @@ -398,12 +270,19 @@ 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); - db = on_debug(); - // db = el_root->debug (); - - 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 (); - 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) { @@ -564,22 +362,22 @@ 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 ; nhaddVertex (px0 + nh1*dx, py0 + nh1*dy, - pz0 + nh1*dz); + 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); @@ -677,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 ; 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; - - 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 ; 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) { @@ -875,32 +580,32 @@ int Elements::cutHexas (const Edges& t_edges, int nbcuts) 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; 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; 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]); @@ -910,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 @@ -1081,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 ; nroisHere () + && quad->definedBy (e1, e2)) + return quad; + } + return NULL; +} END_NAMESPACE_HEXA