2 // C++ : Gestion des cylindres croises
4 // Copyright (C) 2009-2011 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/
21 // or email : webmaster.salome@opencascade.com
23 #include "HexBiCylinder.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"
34 static bool db = false;
38 static const double cos45 = cos (M_PI/4);
41 void geom_define_line (string& brep);
42 void geom_asso_point (double angle, Vertex* node);
44 void geom_create_circle (double* milieu, double rayon, double* normale,
45 double* base, string& brep);
46 int geom_create_cylcyl (double* borig, double* bnorm, double* bbase,
47 double bray, double bhaut,
48 double* sorig, double* snorm, double* sbase,
49 double sray, double shaut);
50 int geom_asso_cylcyl (Edge* edge);
52 // ====================================================== Constructeur
53 BiCylinder::BiCylinder (Document* doc)
61 at_right = at_left = true;
63 tab_hexa.push_back (NULL);
64 tab_quad.push_back (NULL);
65 tab_edge.push_back (NULL);
66 tab_vertex.push_back (NULL);
67 nbr_vertex = nbr_edges = nbr_quads = nbr_hexas = 1;
69 // ====================================================== getHexaIJK
70 Hexa* BiCylinder::getHexaIJK (int cyl, int nx, int ny, int nz)
72 int key = getKey (cyl, nx, ny, nz);
73 int nro = map_hexa [key];
74 return tab_hexa [nro];
76 // ====================================================== getQuadIJ
77 Quad* BiCylinder::getQuadIJ (int cyl, int nx, int ny, int nz)
79 int key = getKey (cyl, dir_z, nx, ny, nz);
80 int nro = map_quad [key];
81 return tab_quad [nro];
83 // ====================================================== getQuadJK
84 Quad* BiCylinder::getQuadJK (int cyl, int nx, int ny, int nz)
86 int key = getKey (cyl, dir_x, nx, ny, nz);
87 int nro = map_quad [key];
88 return tab_quad [nro];
90 // ====================================================== getQuadIK
91 Quad* BiCylinder::getQuadIK (int cyl, int nx, int ny, int nz)
93 int key = getKey (cyl, dir_y, nx, ny, nz);
94 int nro = map_quad [key];
95 return tab_quad [nro];
97 // ====================================================== getEdgeI
98 Edge* BiCylinder::getEdgeI (int cyl, int nx, int ny, int nz)
100 int key = getKey (cyl, dir_x, nx, ny, nz);
101 int nro = map_edge [key];
102 return tab_edge [nro];
104 // ====================================================== getEdgeJ
105 Edge* BiCylinder::getEdgeJ (int cyl, int nx, int ny, int nz)
107 int key = getKey (cyl, dir_y, nx, ny, nz);
108 int nro = map_edge [key];
109 return tab_edge [nro];
111 // ====================================================== getEdgeK
112 Edge* BiCylinder::getEdgeK (int cyl, int nx, int ny, int nz)
114 int key = getKey (cyl, dir_z, nx, ny, nz);
115 int nro = map_edge [key];
116 return tab_edge [nro];
118 // ====================================================== getVertexIJK
119 Vertex* BiCylinder::getVertexIJK (int cyl, int nx, int ny, int nz)
121 int key = getKey (cyl, nx, ny, nz);
122 int nro = map_vertex [key];
125 printf ("getVertexIJK (%d, %d,%d,%d) = V%d = ", cyl, nx, ny, nz, nro);
126 if (tab_vertex[nro] == NULL)
129 printf ("%s = (%g, %g, %g)\n", tab_vertex[nro]->getName(),
130 tab_vertex[nro]->getX(), tab_vertex[nro]->getY(),
131 tab_vertex[nro]->getZ());
133 return tab_vertex [nro];
135 // ------------------------------------------------------------------------
136 // ====================================================== addVertex
137 Vertex* BiCylinder::addVertex (double px, double py, double pz, int cyl,
138 int nx, int ny, int nz)
140 Vertex* node = el_root->addVertex (px, py, pz);
141 int key = getKey (cyl, nx, ny, nz);
142 tab_vertex.push_back (node);
143 map_vertex [key] = nbr_vertex;
147 printf (" tab_vertex [%d, %d,%d,%d] = V%d = ", cyl, nx, ny, nz,
152 printf ("%s = (%g, %g, %g)\n", node->getName(), node->getX(),
153 node->getY(), node->getZ());
158 // ====================================================== addEdge
159 Edge* BiCylinder::addEdge (Vertex* v1, Vertex* v2, int cyl, int dir, int nx,
162 int key = getKey (cyl, dir, nx, ny, nz);
163 int nro = map_edge [key];
168 Edge* edge = tab_edge [nro];
170 printf ("pres_edge [%d,%d, %d,%d,%d] = E%d = ", cyl, dir, nx, ny, nz,
172 printf ("%s ((%s, %s)) = (%s,%s)\n", edge->getName(),
173 edge->getVertex(0)->getName(), edge->getVertex(1)->getName(),
174 v1->getName(), v2->getName());
175 if (NOT edge->definedBy (v1,v2))
176 printf (" **** Incoherence !!\n");
178 return tab_edge [nro];
181 if (v1==NULL || v2==NULL)
184 Edge* edge = findEdge (v1, v2);
186 edge = newEdge (v1, v2);
188 map_edge [key] = nbr_edges;
189 tab_edge.push_back (edge);
193 printf (" tab_edge [%d,%d, %d,%d,%d] = E%d = ", cyl, dir, nx, ny, nz,
198 printf ("%s = (%s, %s)\n", edge->getName(),
199 edge->getVertex(0)->getName(),
200 edge->getVertex(1)->getName());
205 // ====================================================== addQuad
206 Quad* BiCylinder::addQuad (Edge* e1, Edge* e2, Edge* e3, Edge* e4, int cyl,
207 int dir, int nx, int ny, int nz)
209 int key = getKey (cyl, dir, nx, ny, nz);
210 int nro = map_quad [key];
212 return tab_quad [nro];
214 Quad* quad = newQuad (e1, e2, e3, e4);
215 map_quad [key] = nbr_quads;
216 tab_quad.push_back (quad);
220 printf (" tab_quad [%d,%d, %d,%d,%d] = Q%d = ", cyl, dir, nx, ny, nz,
225 printf ("%s = (%s, %s, %s, %s)\n", quad->getName(),
226 quad->getEdge(0)->getName(), quad->getEdge(1)->getName(),
227 quad->getEdge(2)->getName(), quad->getEdge(3)->getName());
232 // ====================================================== addHexa
233 Hexa* BiCylinder::addHexa (Quad* q1, Quad* q2, Quad* q3, Quad* q4, Quad* q5,
234 Quad* q6, int cyl, int nx, int ny, int nz)
236 int key = getKey (cyl, nx, ny, nz);
237 int nro = map_hexa [key];
240 printf (" tab_hexa [%d, %d,%d,%d] = H%d est deja la\n ",
241 cyl, nx, ny, nz, nbr_hexas);
242 return tab_hexa [nro];
245 Hexa* hexa = newHexa (q1, q2, q3, q4, q5, q6);
246 map_hexa [key] = nbr_hexas;
247 tab_hexa.push_back (hexa);
251 printf (" tab_hexa [%d, %d,%d,%d] = H%d = ", cyl, nx, ny, nz, nbr_hexas);
255 printf ("%s = (%s, %s, %s, %s, %s, %s)\n", hexa->getName(),
256 hexa->getQuad(0)->getName(), hexa->getQuad(1)->getName(),
257 hexa->getQuad(2)->getName(), hexa->getQuad(3)->getName(),
258 hexa->getQuad(4)->getName(), hexa->getQuad(5)->getName());
263 // ------------------------------------------------------------------------
264 // ====================================================== findVertex
265 Vertex* BiCylinder::findVertex (double px, double py, double pz,
266 int nx, int ny, int nz)
268 for (it_map=map_vertex.begin() ; it_map!=map_vertex.end() ; ++it_map)
270 int nro = it_map->second;
271 Vertex* elt = tab_vertex[nro];
272 if (elt != NULL && elt->definedBy (px, py, pz))
274 int key = getKey (CylBig, nx, ny, nz);
275 map_vertex [key] = nro;
277 printf ("findVertex [Big,%d,%d,%d] = V%d = '%d' = %s\n",
278 nx, ny, nz, nro, it_map->first,
279 tab_vertex[nro]->getName());
283 printf ("**** Recherche du vertex (%g, %g, %g)\n", px, py, pz);
284 fatal_error ("HexBiCylinder : Vertex non trouve");
287 // ====================================================== findEdge
288 Edge* BiCylinder::findEdge (Vertex* v1, Vertex* v2)
290 int nbedges = tab_edge.size();
291 for (int nro = 0 ; nro<nbedges ; nro++)
293 Edge* elt = tab_edge[nro];
294 if (elt != NULL && elt->definedBy (v1, v2))
299 // ====================================================== findEdge
300 Edge* BiCylinder::findEdge (Vertex* v1, Vertex* v2, int dir, int nx,
303 for (it_map=map_edge.begin() ; it_map!=map_edge.end() ; ++it_map)
305 int nro = it_map->second;
306 Edge* elt = tab_edge[nro];
307 if (elt != NULL && elt->definedBy (v1, v2))
309 int key = getKey (CylBig, dir, nx, ny, nz);
310 map_edge [key] = nro;
312 printf ("findEdge [%d, %d,%d,%d] = E%d = '%d' = %s\n",
313 dir, nx, ny, nz, nro, it_map->first,
318 fatal_error ("HexBiCylinder : Edge non trouve");
321 // ====================================================== findQuad
322 Quad* BiCylinder::findQuad (Vertex* v1, Vertex* v2,
323 int dir, int nx, int ny, int nz)
325 for (it_map=map_quad.begin() ; it_map!=map_quad.end() ; ++it_map)
327 int nro = it_map->second;
328 Quad* elt = tab_quad[nro];
329 if (elt != NULL && elt->definedBy (v1, v2))
331 int key = getKey (CylBig, dir, nx, ny, nz);
332 map_quad [key] = nro;
334 printf ("findQuad [%d, %d,%d,%d] = E%d = '%d' = %s\n",
335 dir, nx, ny, nz, nro, it_map->first,
340 fatal_error ("HexBiCylinder : Quad non trouve");
343 // ====================================================== findQuad
344 Quad* BiCylinder::findQuad (Edge* e1, Edge* e2, int dir, int nx, int ny, int nz)
346 for (it_map=map_quad.begin() ; it_map!=map_quad.end() ; ++it_map)
348 int nro = it_map->second;
349 Quad* elt = tab_quad[nro];
350 if (elt != NULL && elt->definedBy (e1, e2))
352 int key = getKey (CylBig, dir, nx, ny, nz);
353 map_quad [key] = nro;
355 printf ("findQuad [%d, %d,%d,%d] = E%d = '%d' = %s\n",
356 dir, nx, ny, nz, nro, it_map->first, elt->getName());
360 fatal_error ("HexBiCylinder : Quad non trouve");
363 // ====================================================== findHexa
364 Hexa* BiCylinder::findHexa (Quad* q1, Quad* q2, int nx, int ny, int nz)
366 for (it_map=map_hexa.begin() ; it_map!=map_hexa.end() ; ++it_map)
368 int nro = it_map->second;
369 Hexa* elt = tab_hexa[nro];
370 if (elt != NULL && elt->definedBy (q1, q2))
372 int key = getKey (CylBig, nx, ny, nz);
373 map_hexa [key] = nro;
375 printf ("findHexa [%d,%d,%d] = H%d = '%d'\n",
376 nx, ny, nz, nro, it_map->first);
380 fatal_error ("HexBiCylinder : Hexa non trouve");
383 // ====================================================== crossCylinders
384 int BiCylinder::crossCylinders (Cylinder* lun, Cylinder* lautre)
386 if (lun->getRadius() < lautre->getRadius())
397 int ier = cross_cyl2->interCylinder (cross_cyl1, at_left, at_right,
403 cross_rayext [CylSmall] = cross_cyl1->getRadius ();
404 cross_rayext [CylBig] = cross_cyl2->getRadius ();
405 cross_rayint [CylSmall] = cross_rayext [CylSmall] / 2;
406 cross_rayint [CylBig ] = cross_rayext [CylBig ] / 2;
407 cross_hauteur [CylSmall] = cross_cyl1->getHeight ();
408 cross_hauteur [CylBig] = cross_cyl2->getHeight ();
412 adjustLittleSlice (1, 1);
413 adjustLittleSlice (0, 2);
414 adjustLittleSlice (0, 3);
415 adjustLittleSlice (1, 4);
417 Vector* iprim = new Vector (cross_cyl1->getDirection());
418 Vector* kprim = new Vector (cross_cyl2->getDirection());
419 Vector* jprim = new Vector (kprim);
423 jprim->vectoriel (kprim, iprim);
425 // transfoVertices (cross_center, iprim, jprim, kprim);
427 // Real3 snorm, bnorm;
428 // iprim->getCoord (snorm);
429 // kprim->getCoord (bnorm);
431 Real3 bnorm = {0, 0, 1};
432 Real3 snorm = {1, 0, 0};
433 assoCylinders (snorm, bnorm);
435 /*********************************************
438 assoIntersection (NxExt, 1, snorm, bnorm);
439 if (grid_type == GR_BIPIPE)
441 assoIntersection (NxInt, 2, snorm, bnorm);
447 assoIntersection (NxExt, NbrSlices1-1, snorm, bnorm);
448 if (grid_type == GR_BIPIPE)
450 assoIntersection (NxInt, NbrSlices1-2, snorm, bnorm);
454 ******************************************* */
458 // ====================================================== createLittleCyl
459 void BiCylinder::createLittleCyl ()
462 Vertex* vbase = cross_cyl1->getBase();
464 Real lg = cross_hauteur[CylSmall];
465 Real h1 = calc_distance (cross_center, vbase->getPoint (base));
466 Real h2 = cross_rayext[CylBig]*cos45;
467 Real h3 = (cross_rayext[CylBig] - cross_rayint[CylBig])*cos45;
469 double t_haut [NbrVslices] = { -h1, -h2, -h2+h3, h2-h3, h2, lg-h1 };
471 for (int nk=0; nk<NbrVslices ; nk++)
476 addSlice (CylSmall, 0, nk, t_haut[nk], cross_rayint [CylSmall]);
477 addSlice (CylSmall, 1, nk, t_haut[nk], cross_rayext [CylSmall]);
480 addSlice (CylSmall, 1, nk, t_haut[nk], cross_rayext [CylSmall]);
481 addSlice (CylSmall, 2, nk, t_haut[nk], cross_rayext [CylBig]);
484 addSlice (CylSmall, 0, nk, t_haut[nk], cross_rayint [CylSmall]);
485 addSlice (CylSmall, 1, nk, t_haut[nk], cross_rayint [CylBig]);
490 fillSlice (CylSmall, 0,0, 0,2, 1,1, 1,0);
491 fillSlice (CylSmall, 0,2, 1,2, 2,1, 1,1); // OK
492 // fillSlice (CylSmall, 1,1, 0,2, 1,2, 2,1); // Test
493 fillSlice (CylSmall, 1,2, 1,3, 2,4, 2,1, true);
494 fillSlice (CylSmall, 1,3, 0,3, 1,4, 2,4);
495 fillSlice (CylSmall, 0,3, 0,5, 1,5, 1,4);
498 // ====================================================== createBigCyl
499 void BiCylinder::createBigCyl ()
501 const int cyl = CylBig;
503 Vertex* vbase = cross_cyl2->getBase();
504 Real lg = cross_hauteur[cyl];
505 Real rext = cross_rayext [cyl];
506 Real rint = cross_rayint [cyl];
507 Real h1 = calc_distance (cross_center, vbase->getPoint (base));
508 Real h2 = rext * cos45;
510 Real dh = (rext - rint)*cos45;
512 addSlice (CylBig, 0, 0, -h1, rint);
513 addSlice (CylBig, 1, 0, -h1, rext);
514 addSlice (CylBig, 0, 1, -h2+dh, rint, true);
515 addSlice (CylBig, 1, 1, -h2, rext, true);
516 addSlice (CylBig, 0, 2, h2-dh, rint, true);
517 addSlice (CylBig, 1, 2, h2, rext, true);
518 addSlice (CylBig, 0, 3, h3, rint);
519 addSlice (CylBig, 1, 3, h3, rext);
522 fillSlice (CylBig, 0,0, 0,1, 1,1, 1,0);
523 fillSlice (CylBig, 0,2, 0,3, 1,3, 1,2);
525 // ====================================================== adjustLittleSlice
526 void BiCylinder::adjustLittleSlice (int ni, int nk)
528 Vertex* node = getVertexIJK (CylSmall, ni, 0, nk);
532 double grayon2 = cross_rayext[CylBig] * cross_rayext[CylBig];
533 double prayon = cross_rayext[CylSmall];
536 grayon2 = cross_rayint[CylBig] * cross_rayint[CylBig];
537 prayon = cross_rayint[CylSmall];
540 for (int nj=0; nj<NbrCotes ; nj++)
542 node = getVertexIJK (CylSmall, ni, nj, nk);
543 double angle = getAngle (nj);
544 double py = prayon * cos (angle);
545 double pz = prayon * sin (angle);
546 double px = sqrt (grayon2-py*py);
548 double qx = node->getX();
549 if (qx>=0) node->setX ( px);
550 else node->setX (-px);
555 // ====================================================== addSlice
556 void BiCylinder::addSlice (int cyl, int ni, int nk, double hauteur,
557 double rayon, bool find)
559 for (int nj=0 ; nj<NbrCotes ; nj++)
561 double theta = getAngle (nj);
562 double px = rayon*cos(theta);
563 double py = rayon*sin(theta);
564 // Find = forcement le gros cylindre
566 findVertex (px, py, hauteur, ni, nj, nk);
567 else if (cyl==CylBig)
568 addVertex (px, py, hauteur, CylBig, ni, nj, nk);
570 addVertex (hauteur, px, py, CylSmall, ni, nj, nk);
573 // ====================================================== fillSlice
574 /* *****************************************************************
575 H=bed +----bd-----+ bdf=G
579 E=bce +----bc-----+...|...bcf=F
583 D=ade...|...+----ad-|---+ adf=C | I
587 A=ace +----ac-----+ acf=B +-----> K
589 ***************************************************************** */
590 void BiCylinder::fillSlice (int cyl, int nia, int nka, int nib, int nkb,
591 int nic, int nkc, int nid, int nkd,
594 for (int nj0=0; nj0<NbrCotes ; nj0++)
598 if (nj1>=NbrCotes) nj1=0;
599 Vertex* vace = getVertexIJK (cyl, nia, nj0, nka);
600 Vertex* vacf = getVertexIJK (cyl, nib, nj0, nkb);
601 Vertex* vadf = getVertexIJK (cyl, nic, nj0, nkc);
602 Vertex* vade = getVertexIJK (cyl, nid, nj0, nkd);
604 Vertex* vbce = getVertexIJK (cyl, nia, nj1, nka);
605 Vertex* vbcf = getVertexIJK (cyl, nib, nj1, nkb);
606 Vertex* vbdf = getVertexIJK (cyl, nic, nj1, nkc);
607 Vertex* vbde = getVertexIJK (cyl, nid, nj1, nkd);
609 /* *******************
618 ******************* */
620 Edge* eac = addEdge (vace, vacf, cyl, dir_z, nia, nj0, nka);
621 Edge* ead = addEdge (vade, vadf, cyl, dir_z, nid, nj0, nkd);
622 Edge* eae = addEdge (vace, vade, cyl, dir_x, nia, nj0, nka);
623 Edge* eaf = addEdge (vacf, vadf, cyl, dir_x, nib, nj0, nkb);
625 Edge* ebc = addEdge (vbce, vbcf, cyl, dir_z, nia, nj1, nka);
626 Edge* ebd = addEdge (vbde, vbdf, cyl, dir_z, nid, nj1, nkd);
627 Edge* ebe = addEdge (vbce, vbde, cyl, dir_x, nia, nj1, nka);
628 Edge* ebf = addEdge (vbcf, vbdf, cyl, dir_x, nib, nj1, nkb);
630 Edge* ece = addEdge (vace, vbce, cyl, dir_y, nia, nj0, nka);
631 Edge* ecf = addEdge (vacf, vbcf, cyl, dir_y, nib, nj0, nkb);
632 Edge* edf = addEdge (vadf, vbdf, cyl, dir_y, nic, nj0, nkc);
633 Edge* ede = addEdge (vade, vbde, cyl, dir_y, nid, nj0, nkd);
635 Quad* qa = addQuad (eac, eaf, ead, eae, cyl, dir_y, nia , nj0, nka);
636 Quad* qb = addQuad (ebc, ebf, ebd, ebe, cyl, dir_y, nia , nj1, nka);
637 Quad* qc = addQuad (eac, ecf, ebc, ece, cyl, dir_x, nia , nj0, nka);
638 Quad* qd = addQuad (ead, edf, ebd, ede, cyl, dir_x, nid , nj0, nkd);
639 Quad* qe = addQuad (eae, ede, ebe, ece, cyl, dir_z, nia , nj0, nka);
640 Quad* qf = addQuad (eaf, edf, ebf, ecf, cyl, dir_z, nib , nj0, nkb);
642 addHexa (qa, qb, qc, qd, qe, qf, cyl, nia, nj0, nka);
645 // ====================================================== assoCylinders
646 void BiCylinder::assoCylinders (double* snormal, double* gnormal)
648 assoSlice (CylSmall, 0, 0, snormal);
649 assoSlice (CylSmall, 1, 0, snormal);
650 assoSlice (CylSmall, 0, 5, snormal);
651 assoSlice (CylSmall, 1, 5, snormal);
653 for (int nk=0 ; nk<4 ; nk++)
654 for (int ni=0 ; ni<2 ; ni++)
655 assoSlice (CylBig, ni, nk, gnormal);
657 assoIntersection (1, 1, snormal, gnormal);
658 assoIntersection (0, 2, snormal, gnormal);
659 assoIntersection (0, 3, snormal, gnormal);
660 assoIntersection (1, 4, snormal, gnormal);
663 // ====================================================== assoSlice
664 void BiCylinder::assoSlice (int cyl, int nx, int nzs, double* normal)
666 Real3 center, pnt1, pnt2, vbase;
669 int nyd = NbrCotes/2;
671 Vertex* v0 = getVertexIJK (cyl, nx, ny0 , nzs);
672 Vertex* vd = getVertexIJK (cyl, nx, nyd , nzs);
674 if (vd==NULL || v0==NULL)
681 for (int nro=0 ; nro<DIM3 ; nro++)
683 center[nro] = (pnt1[nro] + pnt2[nro])/2;
684 vbase [nro] = (pnt1[nro] - center[nro]);
685 rayon += vbase [nro]*vbase[nro];
687 rayon = sqrt (rayon);
695 geom_create_circle (center, rayon, normal, vbase, brep);
696 geom_define_line (brep); // pour geom_asso_point
698 for (int ny=0 ; ny<NbrCotes ; ny++)
700 assoArc (cyl, nx, ny, nzs, brep, rayon);
703 // ===================================================== assoArc
704 void BiCylinder::assoArc (int cyl, int nx, int ny, int nz, string& brep,
707 static const double Alpha = 1.0 / NbrCotes;
708 static const double Theta = 2*M_PI/ NbrCotes;
709 int ny1 = (ny+1) MODULO NbrCotes;
711 Edge* edge = getEdgeJ (cyl, nx, ny, nz);
712 Vertex* node0 = getVertexIJK (cyl, nx, ny, nz);
713 Vertex* node1 = getVertexIJK (cyl, nx, ny1, nz);
716 printf ("AssoArc : Edge = %s = (%s,%s) -> (%s,%s)\n", edge->getName(),
717 edge->getVertex(0)->getName(), edge->getVertex(1)->getName(),
718 node0->getName(), node1->getName());
720 // Shape* shape = new Shape (brep);
722 // shape->setBounds (ny*Alpha, (ny+1)*Alpha);
723 // edge ->addAssociation (shape);
725 geom_asso_point ( ny *Theta*rayon, node0);
726 geom_asso_point ((ny+1)*Theta*rayon, node1);
728 // ===================================================== assoIntersection
729 int BiCylinder::assoIntersection (int nxs, int nzs, double* snorm,
732 Real3 X_center = {0, 0, 0}; // provisoire, a la place de cross_center
733 Real3 pse, psw, sorig, sbase;
734 Real3 pbe, pbw, borig, bbase;
737 int nyd = NbrCotes/2;
738 int MiddleSlice1 = 3;
740 int nz = nzs < MiddleSlice1 ? 0 : 5;
742 getVertexIJK (CylSmall, nxs, ny0 , nz)->getPoint (pse);
743 getVertexIJK (CylSmall, nxs, nyd , nz)->getPoint (psw);
744 getVertexIJK (CylBig, nxs, ny0 , 0) ->getPoint (pbe);
745 getVertexIJK (CylBig, nxs, nyd , 0) ->getPoint (pbw);
747 double srayon = calc_distance (pse, psw)/2;
748 double brayon = calc_distance (pbe, pbw)/2;
750 calc_milieu (psw, pse, sorig);
751 calc_milieu (pbw, pbe, borig);
752 calc_vecteur (psw, pse, sbase);
753 calc_vecteur (pbw, pbe, bbase);
755 double shaut = calc_distance (X_center, sorig);
756 double bhaut = calc_distance (X_center, borig)*2;
757 double* orig = nzs < MiddleSlice1 ? sorig : X_center; // Pb orientation
769 int ier = geom_create_cylcyl (borig, bnorm, bbase, brayon, bhaut,
770 orig, snorm, sbase, srayon, shaut);
774 for (int ny=0 ; ny<NbrCotes ; ny++)
776 Vertex* node = getVertexIJK (CylSmall, nxs, ny, nzs);
778 node->clearAssociation ();
781 for (int ny=0 ; ny<NbrCotes ; ny++)
783 Edge* edge = getEdgeJ (CylSmall, nxs, ny, nzs);
784 geom_asso_cylcyl (edge);
789 // ======================================================== test_bicylinder
790 BiCylinder* test_bicylinder (Document* docu, int option)
792 Vertex* ori1 = docu->addVertex ( 0,0,0);
793 Vertex* ori2 = docu->addVertex (-5,0,5);
794 Vector* vz = docu->addVector ( 0,0,1);
795 Vector* vx = docu->addVector ( 1,0,0);
797 double r1 = 2*sqrt (2.0);
798 double r2 = 3*sqrt (2.0);
802 Cylinder* cyl1 = docu->addCylinder (ori1, vz, r1, l1);
803 Cylinder* cyl2 = docu->addCylinder (ori2, vx, r2, l2);
805 BiCylinder* grid = new BiCylinder (docu);
806 grid->crossCylinders (cyl1, cyl2);