X-Git-Url: http://git.salome-platform.org/gitweb/?a=blobdiff_plain;f=src%2FHEXABLOCK%2FHexElements_ter.cxx;h=1c129cb3bc49c32f958d7ca7517e8cb102848f4f;hb=6100168b28f7fe487a360409ef4f8e7c38df4450;hp=cf1a35ef194d17407095de80c82148e32861b189;hpb=477b7b7c83bbdc72658b0561639a253afb8ab62d;p=modules%2Fhexablock.git diff --git a/src/HEXABLOCK/HexElements_ter.cxx b/src/HEXABLOCK/HexElements_ter.cxx index cf1a35e..1c129cb 100755 --- a/src/HEXABLOCK/HexElements_ter.cxx +++ b/src/HEXABLOCK/HexElements_ter.cxx @@ -1,7 +1,6 @@ - // C++ : Table d'hexaedres (Evol Versions 3) -// Copyright (C) 2009-2012 CEA/DEN, EDF R&D +// Copyright (C) 2009-2013 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 @@ -22,76 +21,26 @@ #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 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 ; nyaddVertex (px, py, pz); - setVertex (node, nx, ny, nz); - } - if (cyl_closed) - for (int nx=0 ; nx= 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, +// ====================================================== assoCylinders +void Elements::assoCylinders (Vertex* ori, Vector* normal, double angle, RealVector& t_angles) { - int na = t_angles.size(); - bool regul = na == 0; - double alpha = angle/size_hy; + 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 (double* ori, double* vk, double angle, + RealVector& t_angles) +{ + 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); + // Real3 vk = { normal->getDx(), normal->getDy(), normal->getDz() }; + // normer_vecteur (vk); for (int nz=0 ; nzgetX() - 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 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) 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 ; nygetX(), 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 ; nzaddAssociation (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); + 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; // 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 (); + + Real3 sommet; - 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]); + 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 ; nroCoordVertex (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 ; 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; nivaddAssociation (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; @@ -482,9 +383,9 @@ int Elements::propagateAssociation (Edge* orig, Edge* dest, Edge* dir) { 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); } } } @@ -497,83 +398,58 @@ int Elements::propagateAssociation (Edge* orig, Edge* dest, Edge* dir) { 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 ; nrogetBrep(); - // 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 ; nrogetAssociation (); - 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); + map::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= 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