2 // C++ : Gestion des cylindres croises : construction
4 // Copyright (C) 2009-2013 CEA/DEN, EDF R&D
6 // This library is free software; you can redistribute it and/or
7 // modify it under the terms of the GNU Lesser General Public
8 // License as published by the Free Software Foundation; either
9 // version 2.1 of the License.
11 // This library is distributed in the hope that it will be useful,
12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14 // Lesser General Public License for more details.
16 // You should have received a copy of the GNU Lesser General Public
17 // License along with this library; if not, write to the Free Software
18 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
20 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
23 #include "HexCrossElements.hxx"
25 #include "HexDocument.hxx"
26 #include "HexVector.hxx"
27 #include "HexVertex.hxx"
28 #include "HexHexa.hxx"
29 #include "HexEdge.hxx"
31 #include "HexGlobale.hxx"
32 #include "HexCylinder.hxx"
33 #include "HexNewShape.hxx"
34 #include "HexBiCylinderShape.hxx"
35 #include "HexAssoEdge.hxx"
36 #include "HexEdgeShape.hxx"
38 static bool db = false;
42 static const int MaxLevel = 7;
43 static const double epaiss2 = 0.5;
44 static const double UnSur2pi = DEMI/M_PI;
46 double calcul_centre (Vertex* orig, Vertex* inter);
48 // ====================================================== createBigCyl
49 void CrossElements::createBigCyl ()
52 const int kv_mil = cyl_left != NO_PIPE ? 2 : 3 ;
54 enum { k0, k1, k2, k3, k4 };
57 z0 = big_hauteur[0] = - calcul_centre (cross_cyl2->getBase(), cross_center);
58 z4 = big_hauteur[4] = z0 + cross_cyl2->getHeight ();
60 big_hauteur[1] = getVertexIJK (CylSmall, NxExt, S_SE, kv_mil)->getZ();
61 z2 = big_hauteur[2] = getVertexIJK (CylSmall, NxExt, S_E, kv_mil)->getZ();
62 big_hauteur[3] = getVertexIJK (CylSmall, NxExt, S_NE, kv_mil)->getZ();
64 setVertex (CylBig, iv0, 0, k0, 0, 0, z0);
65 addSlice (CylBig, NxInt, k0, z0);
66 addSlice (CylBig, NxExt, k0, z0);
68 setVertex (CylBig, iv0, 0, k4, 0, 0, z4);
69 addSlice (CylBig, NxInt, k4, z4);
70 addSlice (CylBig, NxExt, k4, z4);
72 //------------------------------- Points intermediaires :
74 double xe1 = getVertexIJK (CylBig, NxExt, S_N, k0)->getX();
75 double ye1 = getVertexIJK (CylBig, NxExt, S_N, k0)->getY();
77 double xe2 = getVertexIJK (CylBig, NxExt, S_S, k0)->getX();
78 double ye2 = getVertexIJK (CylBig, NxExt, S_S, k0)->getY();
80 double xi1 = getVertexIJK (CylBig, NxInt, S_N, k0)->getX();
81 double yi1 = getVertexIJK (CylBig, NxInt, S_N, k0)->getY();
83 double xi2 = getVertexIJK (CylBig, NxInt, S_S, k0)->getX();
84 double yi2 = getVertexIJK (CylBig, NxInt, S_S, k0)->getY();
86 //------------------------------- Reprise des vertex du cylindre 1
88 //------------------------------- Centre
89 if (grid_type != GR_BIPIPE)
90 copyVertex (NxExt, S_S, 3, iv0, 0, k1);
92 //------------------------------- Creation Vertex Nord-Sud
93 for (int nk=k1 ; nk<k4 ; nk++)
95 setVertex (CylBig, NxExt, S_N, nk, xe1, ye1, big_hauteur[nk]);
96 setVertex (CylBig, NxExt, S_S, nk, xe2, ye2, big_hauteur[nk]);
97 setVertex (CylBig, NxInt, S_N, nk, xi1, yi1, big_hauteur[nk]);
98 setVertex (CylBig, NxInt, S_S, nk, xi2, yi2, big_hauteur[nk]);
101 //------------------------------- Face k1 externe :
103 copyVertex (NxExt, S_S, 5, NxExt, S_E, k1);
104 copyVertex (NxExt, S_SE, 5, NxExt, S_NE, k1);
105 copyVertex (NxExt, S_SW, 5, NxExt, S_SE, k1);
107 copyVertex (NxExt, S_S, 1, NxExt, S_W, k1);
108 copyVertex (NxExt, S_SE, 1, NxExt, S_NW, k1);
109 copyVertex (NxExt, S_SW, 1, NxExt, S_SW, k1);
111 //------------------------------- Face k1 interne :
113 copyVertex (NxExt, S_S, 2, NxInt, S_W, k1);
114 copyVertex (NxExt, S_SE, 2, NxInt, S_NW, k1);
115 copyVertex (NxExt, S_SW, 2, NxInt, S_SW, k1);
117 copyVertex (NxExt, S_SE, 3, NxInt, S_N, k1); // 16/02/2012
118 copyVertex (NxExt, S_SW, 3, NxInt, S_S, k1); // 16/02/2012
120 copyVertex (NxExt, S_S, 4, NxInt, S_E, k1);
121 copyVertex (NxExt, S_SE, 4, NxInt, S_NE, k1);
122 copyVertex (NxExt, S_SW, 4, NxInt, S_SE, k1);
124 // ------------------------------------------------------------ K2
125 //------------------------------- Centre
126 // copyVertex (iv0, 0, 3, iv0, 0, k2);
127 // ------------------------------ Face k2 externe :
130 copyVertex (NxExt, S_E, 5, NxExt, S_NE, k2);
131 copyVertex (NxExt, S_W, 5, NxExt, S_SE, k2);
133 copyVertex (NxExt, S_E, 1, NxExt, S_NW, k2);
134 copyVertex (NxExt, S_W, 1, NxExt, S_SW, k2);
136 //------------------------------- Face k2 interne :
138 copyVertex (NxExt, S_E, 2, NxInt, S_NW, k2);
139 copyVertex (NxExt, S_W, 2, NxInt, S_SW, k2);
141 copyVertex (NxExt, S_E, 4, NxInt, S_NE, k2);
142 copyVertex (NxExt, S_W, 4, NxInt, S_SE, k2);
144 if (cyl_left == NO_PIPE)
146 addVertex (CylBig, NxExt, S_W, k2, z2, cross_rayon [CylBig][NxExt]);
147 addVertex (CylBig, NxInt, S_W, k2, z2, cross_rayon [CylBig][NxInt]);
150 if (cyl_right == NO_PIPE)
152 addVertex (CylBig, NxExt, S_E, k2, z2, cross_rayon [CylBig][NxExt]);
153 addVertex (CylBig, NxInt, S_E, k2, z2, cross_rayon [CylBig][NxInt]);
156 // ------------------------------------------------------------
158 //------------------------------- Centre
159 if (grid_type != GR_BIPIPE)
160 copyVertex (NxExt, S_N, 3, iv0, 0, k3);
162 //------------------------------- Face k3 externe :
163 copyVertex (NxExt, S_N, 5, NxExt, S_E, k3);
164 copyVertex (NxExt, S_NE, 5, NxExt, S_NE, k3);
165 copyVertex (NxExt, S_NW, 5, NxExt, S_SE, k3);
167 copyVertex (NxExt, S_N, 1, NxExt, S_W, k3);
168 copyVertex (NxExt, S_NE, 1, NxExt, S_NW, k3);
169 copyVertex (NxExt, S_NW, 1, NxExt, S_SW, k3);
171 //------------------------------- Face k3 interne :
173 copyVertex (NxExt, S_N, 2, NxInt, S_W, k3);
174 copyVertex (NxExt, S_NE, 2, NxInt, S_NW, k3);
175 copyVertex (NxExt, S_NW, 2, NxInt, S_SW, k3);
177 copyVertex (NxExt, S_NE, 3, NxInt, S_N, k3); // 16/02/2012
178 copyVertex (NxExt, S_NW, 3, NxInt, S_S, k3); // 16/02/2012
180 copyVertex (NxExt, S_N, 4, NxInt, S_E, k3);
181 copyVertex (NxExt, S_NE, 4, NxInt, S_NE, k3);
182 copyVertex (NxExt, S_NW, 4, NxInt, S_SE, k3);
184 //------------------------------- Remplissage
186 if (grid_type == GR_BIPIPE)
188 for (int nj=0; nj<S_MAXI ; nj++)
190 Vertex* node = getVertexIJK (CylBig, NxInt, nj, 0);
191 double px = node->getX();
192 double py = node->getY();
193 for (int nk=1; nk<4 ; nk++)
195 node = getVertexIJK (CylBig, NxInt, nj, nk);
206 // ====================================================== createLittleCyl
207 void CrossElements::createLittleCyl ()
209 double c1 = calcul_centre (cross_cyl1->getBase(), cross_center);
210 // double cosalpha = cos (angle_inter[CylBig]);
211 double prayext = cross_rayon [CylBig][NxExt]; // * cosalpha;
212 double prayint = cross_rayon [CylBig][NxInt]; // * cosalpha;
214 double t_haut [MaxLevel] = { -c1, -prayext, -prayint,
216 cross_cyl1->getHeight () -c1 };
218 double rm = ( cross_rayon [CylSmall][NxExt]
219 + cross_rayon [CylSmall][NxInt]) / 2;
220 double rc = cross_rayon [CylBig] [NxInt];
222 double t_rayext [MaxLevel] = { -1, -1, rm, rc, rm, -1, -1 };
224 /* *******************************************************************
225 int nkdeb = at_left ? 0 : size_v1z/2;
226 int nkfin = at_right ? size_v1z : size_v1z/2;
227 ******************************************************************* */
232 case IS_HERE : nkdeb = 0;
234 case NO_PIPE : nkdeb = size_v1z/2;
236 case NO_CYL : nkdeb = 1;
241 case IS_HERE : nkfin = size_v1z;
243 case NO_PIPE : nkfin = size_v1z/2;
245 case NO_CYL : nkfin = size_v1z-1;
249 for (int nk = nkdeb ; nk<nkfin ; nk++)
251 double px = t_haut [nk];
252 if (grid_type != GR_BIPIPE)
253 addVertex (CylSmall, 0, 0, nk, px, 0);
255 addSlice (CylSmall, NxInt, nk, px);
256 addSlice (CylSmall, NxExt, nk, px, t_rayext[nk]);
262 // ====================================================== crossCylinders
263 int CrossElements::crossCylinders (Cylinder* lun, Cylinder* lautre, bool fill)
268 if (lun->getRadius() < lautre->getRadius())
279 double pr = cross_cyl1->getRadius ();
280 double gr = cross_cyl2->getRadius ();
283 cout << " **** crossCylinders : les deux rayons sont trop proches"
289 cross_center = cross_cyl2->interCylinder(cross_cyl1, at_left, at_right);
290 if (cross_center==NULL)
296 if (at_left) cyl_left = IS_HERE;
297 else if (grid_type == GR_BIPIPE) cyl_left = NO_PIPE;
298 else cyl_left = NO_CYL;
300 if (at_right) cyl_right = IS_HERE;
301 else if (grid_type == GR_BIPIPE) cyl_right = NO_PIPE;
302 else cyl_right = NO_CYL;
304 double cross_gray1 = cross_cyl1->getRadius ();
305 double cross_gray2 = cross_cyl2->getRadius ();
306 double cross_igray1 = cross_gray1 * epaiss2;
307 double cross_igray2 = cross_gray2 * epaiss2;
311 HexDisplay (cross_gray1);
312 HexDisplay (cross_gray2);
313 HexDisplay (cross_igray2);
314 HexDisplay (cross_igray1);
317 angle_inter [CylSmall] = M_PI/6;
318 angle_inter [CylBig] = asin (cross_gray1/cross_gray2);
319 angle_intermed = angle_inter [CylBig];
321 cross_rayon [CylSmall] [0] = cross_rayon [CylBig] [0] = 0;
323 cross_rayon [CylSmall] [NxInt] = cross_igray1;
324 cross_rayon [CylSmall] [NxExt] = cross_gray1;
325 cross_rayon [CylBig] [NxInt] = cross_igray2;
326 cross_rayon [CylBig] [NxExt] = cross_gray2;
331 double h1=0, h3=0, dx=1 , dy=1;
332 if (cyl_left!=NO_PIPE)
334 adjustLittleSlice (NxExt, 1, NxExt);
335 adjustLittleSlice (NxInt, 1, NxExt);
336 if (grid_type == GR_BIPIPE)
337 adjustLittleSlice (NxInt, 2, NxInt);
338 dx = getVertexIJK (CylBig, NxExt, S_NW, 1)->getX();
339 dy = getVertexIJK (CylBig, NxExt, S_NW, 1)->getY();
340 h1 = getVertexIJK (CylBig, NxExt, S_NW, 1)->getZ();
341 h3 = getVertexIJK (CylBig, NxExt, S_NW, 3)->getZ();
344 if (cyl_right!=NO_PIPE)
346 adjustLittleSlice (NxExt, 5, NxExt);
347 adjustLittleSlice (NxInt, 5, NxExt);
348 if (grid_type == GR_BIPIPE)
349 adjustLittleSlice (NxInt, 4, NxInt);
350 dx = -getVertexIJK (CylBig, NxExt, S_NE, 1)->getX();
351 dy = getVertexIJK (CylBig, NxExt, S_NE, 1)->getY();
352 h1 = getVertexIJK (CylBig, NxExt, S_NE, 1)->getZ();
353 h3 = getVertexIJK (CylBig, NxExt, S_NE, 3)->getZ();
356 for (int ny=0 ; ny<S_MAXI ; ny++)
358 if (ny != S_E && ny != S_W)
360 getVertexIJK (CylBig, NxExt, ny, 1)->setZ(h1);
361 getVertexIJK (CylBig, NxExt, ny, 3)->setZ(h3);
365 angle_intermed = atan (dy/dx);
366 Vector* iprim = new Vector (cross_cyl1->getDirection());
367 Vector* kprim = new Vector (cross_cyl2->getDirection());
368 Vector* jprim = new Vector (kprim);
372 jprim->vectoriel (kprim, iprim);
374 transfoVertices (cross_center, iprim, jprim, kprim);
377 iprim->getCoord (snorm);
378 kprim->getCoord (bnorm);
381 sprintf (name, "grid_%02d", el_id);
382 grid_geom = el_root->addShape (name, SH_INTER);
383 grid_geom -> openShape();
384 assoCylinder (CylSmall, snorm);
385 assoCylinder (CylBig, bnorm);
386 grid_geom -> closeShape();
388 if (cyl_left == IS_HERE)
390 assoIntersection (NxExt, 1, snorm, bnorm);
391 if (grid_type == GR_BIPIPE)
393 assoIntersection (NxInt, 2, snorm, bnorm);
397 if (cyl_right == IS_HERE)
399 assoIntersection (NxExt, NbrSlices1-1, snorm, bnorm);
400 if (grid_type == GR_BIPIPE)
402 assoIntersection (NxInt, NbrSlices1-2, snorm, bnorm);
406 if (cyl_left == IS_HERE)
408 adjustAsso (NxExt, S_N, 1, V_AVAL);
409 adjustAsso (NxExt, S_SW, 1, V_AMONT);
411 adjustAsso (NxExt, S_N, 3, V_AVAL);
412 adjustAsso (NxExt, S_SW, 3, V_AMONT);
415 if (cyl_right == IS_HERE)
417 adjustAsso (NxExt, S_NE, 1, V_AMONT);
418 adjustAsso (NxExt, S_S, 1, V_AVAL);
420 adjustAsso (NxExt, S_NE, 3, V_AMONT);
421 adjustAsso (NxExt, S_S, 3, V_AVAL);
424 // assoResiduelle ();
427 // ====================================================== copyVertex
428 void CrossElements::copyVertex (int i1, int j1, int k1, int i2, int j2, int k2)
430 Vertex* node = getVertexIJK (CylSmall, i1, j1, k1);
433 setVertex (node, CylBig, i2, j2, k2);
437 double rayon = cross_rayon [CylBig][i2];
438 double hauteur = big_hauteur [k2];
440 addVertex (CylBig, i2, j2, k2, hauteur, rayon);
442 // ===================================================== assoCylinder
443 void CrossElements::assoCylinder (int cyl, double* normal)
445 Real3 base, vec1, center, east, west, nordest;
447 if (cyl==CylSmall && cyl_left != IS_HERE)
450 Vertex* v_e = getVertexIJK (cyl, NxExt, S_E , nk);
451 Vertex* v_ne = getVertexIJK (cyl, NxExt, S_NE, nk);
452 Vertex* v_w = getVertexIJK (cyl, NxExt, S_W , nk);
454 v_e->getPoint (east);
455 v_w->getPoint (west);
456 calc_vecteur (west, east, base);
458 for (int nro=0 ; nro<DIM3 ; nro++)
459 center[nro] = (east[nro] + west[nro])/2;
461 v_ne->getPoint (nordest);
462 calc_vecteur (center, nordest, vec1);
464 double ps = prod_scalaire (base, vec1);
465 double pnorm = calc_norme(base) * calc_norme(vec1);
467 angle_inter [cyl] = acos (ps/pnorm);
469 if (cyl==CylBig || at_left)
471 assoSlice (cyl, base, normal, NxExt, 0);
472 if (grid_type == GR_BIPIPE)
474 assoSlice (cyl, base, normal, NxInt, 0);
475 assoSlice (cyl, base, normal, NxInt, 1);
479 if (cyl==CylBig || cyl_right == IS_HERE)
481 assoSlice (cyl, base, normal, NxExt, size_hz[cyl]);
482 if (grid_type == GR_BIPIPE)
484 assoSlice (cyl, base, normal, NxInt, size_hz[cyl]);
485 assoSlice (cyl, base, normal, NxInt, size_hz[cyl]-1);
490 for (int nz=1 ; nz<=3 ; nz++)
491 assoBigMiddle (base, normal, nz);
493 // ===================================================== assoBigMiddle
494 void CrossElements::assoBigMiddle (double* base, double* normal, int nzs)
496 Real3 center, pnt1, pnt2;
499 Vertex* v_n = getVertexIJK (CylBig, nx, S_N , nzs);
500 Vertex* v_s = getVertexIJK (CylBig, nx, S_S , nzs);
502 v_s->getPoint (pnt1);
503 v_n->getPoint (pnt2);
505 double rayon = calc_distance (pnt1, pnt2)/2;
506 // double alpha = angle_inter [CylBig];
509 // double h1 = cross_rayon[CylSmall][NxExt] * cos (angle_inter[CylSmall]);
510 // alpha = asin (h1/cross_rayon[CylBig][NxExt]);
513 for (int nro=0 ; nro<DIM3 ; nro++)
514 center[nro] = (pnt1[nro] + pnt2[nro])/2;
516 int subid = grid_geom->addCircle (center, rayon, normal, base);
518 if (cyl_right != IS_HERE)
520 assoArc (CylBig, NxExt, S_SE, nzs, subid);
521 assoArc (CylBig, NxExt, S_E, nzs, subid);
524 assoArc (CylBig, NxExt, S_NE, nzs, subid);
525 assoArc (CylBig, NxExt, S_N , nzs, subid);
527 if (cyl_left != IS_HERE)
529 assoArc (CylBig, NxExt, S_NW, nzs, subid);
530 assoArc (CylBig, NxExt, S_W , nzs, subid);
533 assoArc (CylBig, NxExt, S_SW, nzs, subid);
534 assoArc (CylBig, NxExt, S_S , nzs, subid);
536 // ====================================================== adjustLittleSlice
537 void CrossElements::adjustLittleSlice (int ni, int nk, int nibig)
539 Vertex* node = getVertexIJK (CylSmall, ni, S_N, nk);
543 double grayon2 = cross_rayon[CylBig][nibig] * cross_rayon[CylBig][nibig];
544 double prayon = cross_rayon[CylSmall][ni];
546 for (int nj=0; nj<S_MAXI ; nj++)
548 node = getVertexIJK (CylSmall, ni, nj, nk);
549 double angle = getAngle (CylSmall, nj);
550 double py = prayon * cos (angle);
551 double pz = prayon * sin (angle);
552 double px = sqrt (grayon2-py*py);
554 double qx = node->getX();
555 if (qx>=0) node->setX ( px);
556 else node->setX (-px);
561 // ===================================================== assoIntersection
562 int CrossElements::assoIntersection (int nxs, int nzs, double* snorm,
566 Real3 pse, psw, sorig, sbase;
567 Real3 pbe, pbw, borig, bbase;
569 cross_center->getPoint (center);
570 int nz = nzs < MiddleSlice1 ? 0 : NbrSlices1;
572 getVertexIJK (CylSmall, nxs, S_E , nz)->getPoint (pse);
573 getVertexIJK (CylSmall, nxs, S_W , nz)->getPoint (psw);
574 getVertexIJK (CylBig, nxs, S_E , 0) ->getPoint (pbe);
575 getVertexIJK (CylBig, nxs, S_W , 0) ->getPoint (pbw);
577 double srayon = calc_distance (pse, psw)/2;
578 double brayon = calc_distance (pbe, pbw)/2;
580 calc_milieu (psw, pse, sorig);
581 calc_milieu (pbw, pbe, borig);
582 calc_vecteur (psw, pse, sbase);
583 calc_vecteur (pbw, pbe, bbase);
585 double shaut = calc_distance (center, sorig);
586 double bhaut = calc_distance (center, borig)*2;
587 double* orig = nzs < MiddleSlice1 ? sorig : center; // Pb orientation
589 BiCylinderShape bicyl_shape (el_root);
590 int ier = bicyl_shape.defineCyls (borig, bnorm, bbase, brayon, bhaut,
591 orig, snorm, sbase, srayon, shaut);
595 for (int ny=S_E ; ny<=S_SE ; ny++)
597 Vertex* node = getVertexIJK (CylSmall, nxs, ny, nzs);
599 node->clearAssociation ();
602 for (int ny=S_E ; ny<=S_SE ; ny++)
604 Edge* edge = getEdgeJ (CylSmall, nxs, ny, nzs);
605 bicyl_shape.associate (edge);
610 // ===================================================== adjustAsso
611 void CrossElements::adjustAsso (int nx, int ny, int nz, int sens)
613 Edge* edge = getEdgeJ (CylBig, nx, ny, nz);
618 Vertex* vertex = edge->getVertex (sens);
619 AssoEdge* asso = edge->getAssociation (0);
623 vertex->getPoint (point);
625 EdgeShape* line = asso->getEdgeShape ();
626 double param = line->getParam (point);
631 asso->setStart (param);
633 asso->setEnd (param);