2 // C++ : Gestion des cylindres croises
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/
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 "HexNewShape.hxx"
33 #include "HexBiCylinderShape.hxx"
35 /// static const double UnSur2pi = DEMI/M_PI;
39 static bool db = on_debug ();
40 static const double cos45 = cos (M_PI/4);
41 #define NameOf(x) ((x)?x->getName():"Null")
43 // ====================================================== Constructeur
44 BiCylinder::BiCylinder (Document* doc)
48 grid_type = GR_BIPIPE;
49 at_right = at_left = true;
51 tab_hexa.push_back (NULL);
52 tab_quad.push_back (NULL);
53 tab_edge.push_back (NULL);
54 tab_vertex.push_back (NULL);
55 nbr_vertex = nbr_edges = nbr_quads = nbr_hexas = 1;
59 // ====================================================== getHexaIJK
60 Hexa* BiCylinder::getHexaIJK (int cyl, int nx, int ny, int nz)
62 int key = getKey (cyl, nx, ny, nz);
63 int nro = map_hexa [key];
64 return tab_hexa [nro];
66 // ====================================================== getQuadIJ
67 Quad* BiCylinder::getQuadIJ (int cyl, int nx, int ny, int nz)
69 int key = getKey (cyl, dir_z, nx, ny, nz);
70 int nro = map_quad [key];
71 return tab_quad [nro];
73 // ====================================================== getQuadJK
74 Quad* BiCylinder::getQuadJK (int cyl, int nx, int ny, int nz)
76 int key = getKey (cyl, dir_x, nx, ny, nz);
77 int nro = map_quad [key];
78 return tab_quad [nro];
80 // ====================================================== getQuadIK
81 Quad* BiCylinder::getQuadIK (int cyl, int nx, int ny, int nz)
83 int key = getKey (cyl, dir_y, nx, ny, nz);
84 int nro = map_quad [key];
85 return tab_quad [nro];
87 // ====================================================== getEdgeI
88 Edge* BiCylinder::getEdgeI (int cyl, int nx, int ny, int nz)
90 int key = getKey (cyl, dir_x, nx, ny, nz);
91 int nro = map_edge [key];
92 return tab_edge [nro];
94 // ====================================================== getEdgeJ
95 Edge* BiCylinder::getEdgeJ (int cyl, int nx, int ny, int nz)
97 int key = getKey (cyl, dir_y, nx, ny, nz);
98 int nro = map_edge [key];
99 return tab_edge [nro];
101 // ====================================================== getEdgeK
102 Edge* BiCylinder::getEdgeK (int cyl, int nx, int ny, int nz)
104 int key = getKey (cyl, dir_z, nx, ny, nz);
105 int nro = map_edge [key];
106 return tab_edge [nro];
108 // ====================================================== getVertexIJK
109 Vertex* BiCylinder::getVertexIJK (int cyl, int nx, int ny, int nz)
111 int key = getKey (cyl, nx, ny, nz);
112 int nro = map_vertex [key];
115 printf ("getVertexIJK (%d, %d,%d,%d) = V%d = ", cyl, nx, ny, nz, nro);
116 if (tab_vertex[nro] == NULL)
119 printf ("%s = (%g, %g, %g)\n", NameOf (tab_vertex[nro]),
120 tab_vertex[nro]->getX(), tab_vertex[nro]->getY(),
121 tab_vertex[nro]->getZ());
123 return tab_vertex [nro];
125 // ------------------------------------------------------------------------
126 // ====================================================== addVertex
127 Vertex* BiCylinder::addVertex (double px, double py, double pz, int cyl,
128 int nx, int ny, int nz)
130 Vertex* node = el_root->addVertex (px, py, pz);
131 int key = getKey (cyl, nx, ny, nz);
132 tab_vertex.push_back (node);
133 map_vertex [key] = nbr_vertex;
137 printf (" tab_vertex [%d, %d,%d,%d] = V%d = ", cyl, nx, ny, nz,
142 printf ("%s = (%g, %g, %g)\n", NameOf (node), node->getX(),
143 node->getY(), node->getZ());
148 // ====================================================== addEdge
149 Edge* BiCylinder::addEdge (Vertex* v1, Vertex* v2, int cyl, int dir, int nx,
152 int key = getKey (cyl, dir, nx, ny, nz);
153 int nro = map_edge [key];
158 Edge* edge = tab_edge [nro];
160 printf ("pres_edge [%d,%d, %d,%d,%d] = E%d = ", cyl, dir, nx, ny, nz,
162 printf ("%s ((%s, %s)) = (%s,%s)\n", NameOf(edge),
163 NameOf(edge->getVertex(0)), NameOf(edge->getVertex(1)),
164 NameOf(v1), NameOf(v2));
165 if (edge->definedBy (v1,v2))
166 return tab_edge [nro];
168 printf (" **** Incoherence !!\n");
172 if (v1==NULL || v2==NULL)
175 Edge* edge = findEdge (v1, v2);
179 edge = newEdge (v1, v2);
180 map_edge [key] = nbr_edges;
181 tab_edge.push_back (edge);
186 printf (" tab_edge [%d,%d, %d,%d,%d] = E%d = ", cyl, dir, nx, ny, nz,
191 printf ("%s = (%s, %s)\n", NameOf(edge),
192 NameOf(edge->getVertex(0)),
193 NameOf(edge->getVertex(1)));
197 // ====================================================== addQuad
198 Quad* BiCylinder::addQuad (Edge* e1, Edge* e2, Edge* e3, Edge* e4, int cyl,
199 int dir, int nx, int ny, int nz)
201 int key = getKey (cyl, dir, nx, ny, nz);
202 int nro = map_quad [key];
207 Quad* quad = tab_quad [nro];
209 printf ("pres_quad [%d,%d, %d,%d,%d] = Q%d = ", cyl, dir, nx, ny, nz,
211 printf ("%s (%s, %s, %s, %s)\n", NameOf(quad),
212 NameOf(e1), NameOf(e2), NameOf(e3), NameOf(e4));
213 if (quad->definedBy (e1,e3))
214 return tab_quad [nro];
217 printf (" **** Incoherence !!\n");
218 printf ("%s = (%s, %s, %s, %s)\n", NameOf(quad),
219 NameOf(quad->getEdge(0)), NameOf(quad->getEdge(1)),
220 NameOf(quad->getEdge(1)), NameOf(quad->getEdge(2)));
221 printf ("%s = (%s, %s, %s, %s)\n", NameOf(quad),
222 NameOf(quad->getVertex(0)), NameOf(quad->getVertex(1)),
223 NameOf(quad->getVertex(1)), NameOf(quad->getVertex(2)));
228 Quad* quad = Elements::findQuad (e1, e2);
232 quad = newQuad (e1, e2, e3, e4);
233 map_quad [key] = nbr_quads;
234 tab_quad.push_back (quad);
239 printf (" tab_quad [%d,%d, %d,%d,%d] = Q%d = ", cyl, dir, nx, ny, nz,
244 printf ("%s = (%s, %s, %s, %s)\n", NameOf(quad),
245 NameOf(quad->getEdge(0)), NameOf(quad->getEdge(1)),
246 NameOf(quad->getEdge(2)), NameOf(quad->getEdge(3)));
250 // ====================================================== addHexa
251 Hexa* BiCylinder::addHexa (Quad* q1, Quad* q2, Quad* q3, Quad* q4, Quad* q5,
252 Quad* q6, int cyl, int nx, int ny, int nz)
254 int key = getKey (cyl, nx, ny, nz);
255 int nro = map_hexa [key];
258 if (tab_hexa [nro]->definedBy (q1,q2))
260 printf (" tab_hexa [%d, %d,%d,%d] = H%d est deja la\n ",
261 cyl, nx, ny, nz, nbr_hexas);
262 return tab_hexa [nro];
264 printf (" Incoherence : tab_hexa [%d, %d,%d,%d] = H%d = %s\n ",
265 cyl, nx, ny, nz, nbr_hexas, tab_hexa [nro]->getName ());
268 Hexa* hexa = newHexa (q1, q2, q3, q4, q5, q6);
269 map_hexa [key] = nbr_hexas;
270 tab_hexa.push_back (hexa);
274 printf (" tab_hexa [%d, %d,%d,%d] = H%d = ", cyl, nx, ny, nz, nbr_hexas);
278 printf ("%s = (%s, %s, %s, %s, %s, %s)\n", NameOf(hexa),
279 NameOf(hexa->getQuad(0)), NameOf(hexa->getQuad(1)),
280 NameOf(hexa->getQuad(2)), NameOf(hexa->getQuad(3)),
281 NameOf(hexa->getQuad(4)), NameOf(hexa->getQuad(5)));
286 // ------------------------------------------------------------------------
287 // ====================================================== findVertex
288 Vertex* BiCylinder::findVertex (double px, double py, double pz,
289 int nx, int ny, int nz)
291 for (it_map=map_vertex.begin() ; it_map!=map_vertex.end() ; ++it_map)
293 int nro = it_map->second;
294 Vertex* elt = tab_vertex[nro];
295 if (elt != NULL && elt->definedBy (px, py, pz))
297 int key = getKey (CylBig, nx, ny, nz);
298 map_vertex [key] = nro;
300 printf ("findVertex [Big,%d,%d,%d] = V%d = '%d' = %s\n",
301 nx, ny, nz, nro, it_map->first,
302 NameOf(tab_vertex[nro]));
306 printf ("**** Recherche du vertex (%g, %g, %g)\n", px, py, pz);
307 fatal_error ("HexBiCylinder : Vertex non trouve");
310 // ====================================================== findEdge
311 Edge* BiCylinder::findEdge (Vertex* v1, Vertex* v2)
313 int nbedges = tab_edge.size();
314 for (int nro = 0 ; nro<nbedges ; nro++)
316 Edge* elt = tab_edge[nro];
317 if (elt != NULL && elt->definedBy (v1, v2))
322 // ====================================================== findEdge
323 Edge* BiCylinder::findEdge (Vertex* v1, Vertex* v2, int dir, int nx,
326 for (it_map=map_edge.begin() ; it_map!=map_edge.end() ; ++it_map)
328 int nro = it_map->second;
329 Edge* elt = tab_edge[nro];
330 if (elt != NULL && elt->definedBy (v1, v2))
332 int key = getKey (CylBig, dir, nx, ny, nz);
333 map_edge [key] = nro;
335 printf ("findEdge [%d, %d,%d,%d] = E%d = '%d' = %s\n",
336 dir, nx, ny, nz, nro, it_map->first,
341 fatal_error ("HexBiCylinder : Edge non trouve");
344 // ====================================================== findQuad
345 Quad* BiCylinder::findQuad (Vertex* v1, Vertex* v2,
346 int dir, int nx, int ny, int nz)
348 for (it_map=map_quad.begin() ; it_map!=map_quad.end() ; ++it_map)
350 int nro = it_map->second;
351 Quad* elt = tab_quad[nro];
352 if (elt != NULL && elt->definedBy (v1, v2))
354 int key = getKey (CylBig, dir, nx, ny, nz);
355 map_quad [key] = nro;
357 printf ("findQuad [%d, %d,%d,%d] = E%d = '%d' = %s\n",
358 dir, nx, ny, nz, nro, it_map->first,
363 fatal_error ("HexBiCylinder : Quad non trouve");
366 // ====================================================== findQuad
367 Quad* BiCylinder::findQuad (Edge* e1, Edge* e2, int dir, int nx, int ny, int nz)
369 for (it_map=map_quad.begin() ; it_map!=map_quad.end() ; ++it_map)
371 int nro = it_map->second;
372 Quad* elt = tab_quad[nro];
373 if (elt != NULL && elt->definedBy (e1, e2))
375 int key = getKey (CylBig, dir, nx, ny, nz);
376 map_quad [key] = nro;
378 printf ("findQuad [%d, %d,%d,%d] = E%d = '%d' = %s\n",
379 dir, nx, ny, nz, nro, it_map->first, NameOf(elt));
383 fatal_error ("HexBiCylinder : Quad non trouve");
386 /*******************************************
387 // ====================================================== findQuad
388 Quad* BiCylinder::findQuad (Edge* e1, Edge* e2)
390 for (it_map=map_quad.begin() ; it_map!=map_quad.end() ; ++it_map)
392 int nro = it_map->second;
393 Quad* elt = tab_quad[nro];
394 if (elt != NULL && elt->definedBy (e1, e2))
396 int key = getKey (CylBig, dir, nx, ny, nz);
397 map_quad [key] = nro;
399 printf ("findQuad [%d, %d,%d,%d] = E%d = '%d' = %s\n",
400 dir, nx, ny, nz, nro, it_map->first, NameOf(elt));
404 fatal_error ("HexBiCylinder : Quad non trouve");
407 ********************************************/
408 // ====================================================== findHexa
409 Hexa* BiCylinder::findHexa (Quad* q1, Quad* q2, int nx, int ny, int nz)
411 for (it_map=map_hexa.begin() ; it_map!=map_hexa.end() ; ++it_map)
413 int nro = it_map->second;
414 Hexa* elt = tab_hexa[nro];
415 if (elt != NULL && elt->definedBy (q1, q2))
417 int key = getKey (CylBig, nx, ny, nz);
418 map_hexa [key] = nro;
420 printf ("findHexa [%d,%d,%d] = H%d = '%d'\n",
421 nx, ny, nz, nro, it_map->first);
425 fatal_error ("HexBiCylinder : Hexa non trouve");
428 // ====================================================== createLittlePipe
429 void BiCylinder::createLittlePipe ()
431 Real lg = cross_hauteur[CylSmall];
432 Real h1 = calc_distance (cross_center, cross_orismall);
433 Real h2 = cross_rayext[CylBig]*cos45;
434 Real h3 = (cross_rayext[CylBig] - cross_rayint[CylBig])*cos45;
437 // double t_haut [NbrVslices] = { -h1, -h2, -h2+h3, h2-h3, h2, lg-h1 };
441 addSlice (CylSmall, 0, 0, -h1, cross_rayint [CylSmall]);
442 addSlice (CylSmall, 1, 0, -h1, cross_rayext [CylSmall]);
443 addSlice (CylSmall, 1, 1, -h2, cross_rayext [CylSmall]);
444 addSlice (CylSmall, 0, 2, h3-h2, cross_rayint [CylSmall]);
448 addSlice (CylSmall, 0, 1, -h2, cross_rayint [CylSmall]);
450 addSlice (CylSmall, 1, 1, -h2, cross_rayext [CylSmall]);
451 addSlice (CylSmall, 0, 2, h3-h2, cross_rayint [CylSmall]);
455 addSlice (CylSmall, 2, 1, -h2, cross_rayext [CylBig]);
456 addSlice (CylSmall, 1, 2, h3-h2, cross_rayint [CylBig]);
457 addSlice (CylSmall, 1, 3, h2-h3, cross_rayint [CylBig]);
458 addSlice (CylSmall, 2, 4, h2, cross_rayext [CylBig]);
462 addSlice (CylSmall, 0, 3, h2-h3, cross_rayint [CylSmall]);
463 addSlice (CylSmall, 1, 4, h2, cross_rayext [CylSmall]);
464 addSlice (CylSmall, 0, 5, h5, cross_rayint [CylSmall]);
465 addSlice (CylSmall, 1, 5, h5, cross_rayext [CylSmall]);
469 addSlice (CylSmall, 0, 3, h2-h3, cross_rayint [CylSmall]);
470 addSlice (CylSmall, 1, 4, h2, cross_rayext [CylSmall]);
472 addSlice (CylSmall, 0, 4, h2, cross_rayint [CylSmall]);
474 //--------------------------------- Remplissage
478 fillSlice (CylSmall, 0,0, 0,2, 1,1, 1,0);
479 fillSlice (CylSmall, 0,2, 1,2, 2,1, 1,1);
482 addCube (CylSmall, 0,0,2);
487 fillSlice (CylSmall, 0,1, 0,2, 1,2, 2,1);
488 addCube (CylSmall, 1,0);
491 addCube ( CylSmall, 1, 2);
493 fillSlice (CylSmall, 1,2, 1,3, 2,4, 2,1, true);
496 addCube (CylSmall, 2,0, 3);
497 fillSlice (CylSmall, 0,2, 0,3, 1,3, 1,2);
502 fillSlice (CylSmall, 1,3, 0,3, 1,4, 2,4);
503 fillSlice (CylSmall, 0,3, 0,5, 1,5, 1,4);
505 addCube (CylSmall, 3, 0, 5);
509 fillSlice (CylSmall, 0,3, 0,4, 2,4, 1,3);
510 addCube (CylSmall, 3, 0);
513 addCube ( CylSmall, 3, 1);
515 // ====================================================== createBigPipe
516 void BiCylinder::createBigPipe ()
518 const int cyl = CylBig;
519 Real lg = cross_hauteur[cyl];
520 Real rext = cross_rayext [cyl];
521 Real rint = cross_rayint [cyl];
523 Real h1 = calc_distance (cross_center, cross_oribig);
524 Real h2 = rext * cos45;
526 Real dh = (rext - rint)*cos45;
528 addSlice (CylBig, 0, 0, -h1, rint);
529 addSlice (CylBig, 1, 0, -h1, rext);
530 addSlice (CylBig, 0, 1, -h2+dh, rint, true);
531 addSlice (CylBig, 1, 1, -h2, rext, true);
532 addSlice (CylBig, 0, 2, h2-dh, rint, true);
533 addSlice (CylBig, 1, 2, h2, rext, true);
534 addSlice (CylBig, 0, 3, h3, rint);
535 addSlice (CylBig, 1, 3, h3, rext);
540 fillSlice (CylBig, 0,0, 0,1, 1,1, 1,0);
543 fillSlice (CylBig, 0,2, 0,3, 1,3, 1,2);
545 // ====================================================== adjustLittleSlice
546 void BiCylinder::adjustLittleSlice (int ni, int nk)
548 Vertex* node = getVertexIJK (CylSmall, ni, 0, nk);
552 double grayon2 = cross_rayext[CylBig] * cross_rayext[CylBig];
553 double prayon = cross_rayext[CylSmall];
557 grayon2 = cross_rayint[CylBig] * cross_rayint[CylBig];
558 prayon = cross_rayint[CylSmall];
561 for (int nj=0; nj<NbrCotes ; nj++)
563 node = getVertexIJK (CylSmall, ni, nj, nk);
564 double angle = getAngle (nj);
565 double py = prayon * cos (angle);
566 double pz = prayon * sin (angle);
567 double px = sqrt (grayon2-py*py);
569 double qx = node->getX();
570 if (qx>=0) node->setX ( px);
571 else node->setX (-px);
576 // ====================================================== addSlice
577 void BiCylinder::addSlice (int cyl, int ni, int nk, double hauteur,
578 double rayon, bool find)
580 for (int nj=0 ; nj<NbrCotes ; nj++)
582 double theta = getAngle (nj);
583 double px = rayon*cos(theta);
584 double py = rayon*sin(theta);
585 // Find = forcement le gros cylindre
587 findVertex (px, py, hauteur, ni, nj, nk);
588 else if (cyl==CylBig)
589 addVertex (px, py, hauteur, CylBig, ni, nj, nk);
591 addVertex (hauteur, px, py, CylSmall, ni, nj, nk);
594 // ====================================================== addCarre
595 void BiCylinder::addCarre (int cyl, int nk, double hauteur, double rayon,
598 for (int nj=0 ; nj<NbrCotes ; nj++)
600 double theta = getAngle (nj);
601 double px = rayon*cos(theta);
602 double py = rayon*sin(theta);
604 findVertex (px, py, hauteur, 0, nj, nk);
605 else if (cyl==CylBig)
606 addVertex (px, py, hauteur, CylBig, 0, nj, nk);
608 addVertex (hauteur, px, py, cyl, 0, nj, nk);
611 // ====================================================== adjustLittleSquare
612 void BiCylinder::adjustLittleSquare (int nk)
614 double grayon2 = cross_rayext[CylBig] * cross_rayext[CylBig];
616 for (int nj=0 ; nj<NbrCotes ; nj++)
618 Vertex* node = getVertexIJK (CylSmall, 0, nj, nk);
619 double py = node->getY();
620 double px = sqrt (grayon2 - py*py);
622 if (node->getX() > 0) node->setX ( px);
623 else node->setX (-px);
626 // ====================================================== addCube
627 void BiCylinder::addCube (int cyl, int nk0, int nvi0, int k1)
629 /* *****************************************************************
630 H=bed +----bd-----+ bdf=G......0
634 E=bce +----bc-----+...|...bcf=F
638 D=ade...|...+----ad-|---+ adf=C.....3 | I
642 A=ace +----ac-----+ acf=B.....2 +-----> K
644 ***************************************************************** */
645 enum { nj0, nj1, nj2, nj3, ni0=0 };
646 const int nk1 = k1>0 ? k1 : nk0+1;
647 const int nvi1 = nvi0==1 ? 2
651 Vertex* vacf = getVertexIJK (cyl, nvi1, nj2, nk1);
652 Vertex* vadf = getVertexIJK (cyl, nvi1, nj3, nk1);
653 Vertex* vbdf = getVertexIJK (cyl, nvi1, nj0, nk1);
654 Vertex* vbcf = getVertexIJK (cyl, nvi1, nj1, nk1);
656 Edge* eaf = addEdge (vacf, vadf, cyl, dir_y, ni0, nj2, nk1);
657 Edge* edf = addEdge (vadf, vbdf, cyl, dir_y, ni0, nj3, nk1);
658 Edge* ebf = addEdge (vbcf, vbdf, cyl, dir_y, ni0, nj0, nk1);
659 Edge* ecf = addEdge (vacf, vbcf, cyl, dir_y, ni0, nj1, nk1);
661 Quad* qf = addQuad (eaf, edf, ebf, ecf, cyl, dir_z, ni0 , nj0, nk1);
663 Vertex* vace = getVertexIJK (cyl, nvi0, nj2, nk0);
664 Vertex* vade = getVertexIJK (cyl, nvi0, nj3, nk0);
665 Vertex* vbde = getVertexIJK (cyl, nvi0, nj0, nk0);
666 Vertex* vbce = getVertexIJK (cyl, nvi0, nj1, nk0);
668 Edge* eac = addEdge (vace, vacf, cyl, dir_z, ni0, nj2, nk0);
669 Edge* ead = addEdge (vade, vadf, cyl, dir_z, ni0, nj3, nk0);
670 Edge* ebd = addEdge (vbde, vbdf, cyl, dir_z, ni0, nj0, nk0);
671 Edge* ebc = addEdge (vbce, vbcf, cyl, dir_z, ni0, nj1, nk0);
673 Edge* eae = addEdge (vace, vade, cyl, dir_y, ni0, nj2, nk0);
674 Edge* ede = addEdge (vade, vbde, cyl, dir_y, ni0, nj3, nk0);
675 Edge* ebe = addEdge (vbce, vbde, cyl, dir_y, ni0, nj0, nk0);
676 Edge* ece = addEdge (vace, vbce, cyl, dir_y, ni0, nj1, nk0);
678 // Quad* qe = addQuad (eae, ede, ebe, ece, cyl, dir_z, ni0 , nj0, nk0);
679 Quad* qe = addQuad (eae, ece, ebe, ede, cyl, dir_z, ni0 , nj0, nk0);
681 // Quad* qa = addQuad (eac, eaf, ead, eae, cyl, dir_y, ni0 , nj2, nk0);
682 Quad* qa = addQuad (eac, eae, ead, eaf, cyl, dir_y, ni0 , nj2, nk0);
683 // Quad* qd = addQuad (ead, edf, ebd, ede, cyl, dir_x, ni0 , nj3, nk0);
684 Quad* qd = addQuad (ead, ede, ebd, edf, cyl, dir_x, ni0 , nj3, nk0);
685 Quad* qb = addQuad (ebc, ebf, ebd, ebe, cyl, dir_y, ni0 , nj0, nk0);
686 Quad* qc = addQuad (eac, ecf, ebc, ece, cyl, dir_x, ni0 , nj1, nk0);
688 addHexa (qa, qb, qc, qd, qe, qf, cyl, ni0, nj0, nk0);
690 // ====================================================== fillSlice
691 /* *****************************************************************
692 H=bde +----bd-----+ bdf=G
696 E=bce +----bc-----+...|...bcf=F
700 D=ade...|...+----ad-|---+ adf=C | I
704 A=ace +----ac-----+ acf=B +-----> K
706 ***************************************************************** */
707 void BiCylinder::fillSlice (int cyl, int nia, int nka, int nib, int nkb,
708 int nic, int nkc, int nid, int nkd,
711 for (int nj0=0; nj0<NbrCotes ; nj0++)
715 if (nj1>=NbrCotes) nj1=0;
716 Vertex* vace = getVertexIJK (cyl, nia, nj0, nka);
717 Vertex* vacf = getVertexIJK (cyl, nib, nj0, nkb);
718 Vertex* vadf = getVertexIJK (cyl, nic, nj0, nkc);
719 Vertex* vade = getVertexIJK (cyl, nid, nj0, nkd);
721 Vertex* vbce = getVertexIJK (cyl, nia, nj1, nka);
722 Vertex* vbcf = getVertexIJK (cyl, nib, nj1, nkb);
723 Vertex* vbdf = getVertexIJK (cyl, nic, nj1, nkc);
724 Vertex* vbde = getVertexIJK (cyl, nid, nj1, nkd);
726 Edge* eac = addEdge (vace, vacf, cyl, dir_z, nia, nj0, nka);
727 Edge* ead = addEdge (vadf, vade, cyl, dir_z, nid, nj0, nkd);
728 Edge* eae = addEdge (vade, vace, cyl, dir_x, nia, nj0, nka);
729 Edge* eaf = addEdge (vacf, vadf, cyl, dir_x, nib, nj0, nkb);
731 Edge* ebc = addEdge (vbce, vbcf, cyl, dir_z, nia, nj1, nka);
732 Edge* ebd = addEdge (vbdf, vbde, cyl, dir_z, nid, nj1, nkd);
733 Edge* ebe = addEdge (vbde, vbce, cyl, dir_x, nia, nj1, nka);
734 Edge* ebf = addEdge (vbcf, vbdf, cyl, dir_x, nib, nj1, nkb);
736 Edge* ece = addEdge (vace, vbce, cyl, dir_y, nia, nj0, nka);
737 Edge* ecf = addEdge (vacf, vbcf, cyl, dir_y, nib, nj0, nkb);
738 Edge* edf = addEdge (vadf, vbdf, cyl, dir_y, nic, nj0, nkc);
739 Edge* ede = addEdge (vade, vbde, cyl, dir_y, nid, nj0, nkd);
741 // Quad* qa = addQuad (eac, eaf, ead, eae, cyl, dir_y, nia+1, nj0, nka);
742 Quad* qa = addQuad (eac, eae, ead, eaf, cyl, dir_y, nia+1, nj0, nka);
743 Quad* qb = addQuad (ebc, ebf, ebd, ebe, cyl, dir_y, nia+1, nj1, nka);
745 Quad* qc = addQuad (eac, ecf, ebc, ece, cyl, dir_x, nia+1, nj0, nka);
746 // Quad* qd = addQuad (ead, edf, ebd, ede, cyl, dir_x, nid+1, nj0, nkd);
747 Quad* qd = addQuad (ead, ede, ebd, edf, cyl, dir_x, nid+1, nj0, nkd);
749 // Quad* qe = addQuad (eae, ede, ebe, ece, cyl, dir_z, nia+1, nj0, nka);
750 Quad* qe = addQuad (eae, ece, ebe, ede, cyl, dir_z, nia+1, nj0, nka);
751 Quad* qf = addQuad (eaf, edf, ebf, ecf, cyl, dir_z, nib+1, nj0, nkb);
753 addHexa (qa, qb, qc, qd, qe, qf, cyl, nia+1, nj0, nka);
756 // ====================================================== assoCylinders
757 void BiCylinder::assoCylinders ()
760 sprintf (name, "grid_%02d", el_id);
761 grid_geom = el_root->addShape (name, SH_INTER);
762 grid_geom -> openShape();
765 int imin = cyl_fill ? 1 : 0;
768 for (int ni=imin ; ni<2 ; ni++)
770 assoSlice (CylSmall, ni, 0, cross_dirsmall);
771 assoSlice (CylSmall, ni, s_kmax, cross_dirsmall);
773 for (int nk=0 ; nk<g_ksize ; nk++)
774 assoSlice (CylBig, ni, nk, cross_dirbig);
780 assoIntersection (1, 1);
782 assoIntersection (0, 2);
787 assoIntersection (1, 4);
789 assoIntersection (0, 3);
792 grid_geom -> closeShape();
794 // ====================================================== assoSlice
795 void BiCylinder::assoSlice (int cyl, int nx, int nzs, double* normal)
797 Real3 center, pnt0, pnt1, pnt2, vbase;
800 int ny2 = NbrCotes/2;
802 Vertex* v0 = getVertexIJK (cyl, nx, ny0 , nzs);
803 Vertex* v1 = getVertexIJK (cyl, nx, ny1 , nzs);
804 Vertex* v2 = getVertexIJK (cyl, nx, ny2 , nzs);
806 if (v0==NULL || v1==NULL || v2==NULL)
814 for (int nro=0 ; nro<DIM3 ; nro++)
816 center[nro] = (pnt0[nro] + pnt2[nro])/2;
817 vbase [nro] = (pnt0[nro] - pnt1[nro]);
818 rayon += carre (center [nro] - pnt0 [nro]);
821 int subid = grid_geom->addCircle (center, sqrt(rayon), normal, vbase);
822 for (int ny=0 ; ny<NbrCotes ; ny++)
824 int ny1 = (ny+1) MODULO NbrCotes;
825 Vertex* va = getVertexIJK (cyl, nx, ny , nzs);
826 Vertex* vb = getVertexIJK (cyl, nx, ny1, nzs);
828 assoArc (cyl, ny, va, vb, subid);
831 // ===================================================== assoArc
832 void BiCylinder::assoArc (int cyl, int ny, Vertex* v1, Vertex* v2, int subid)
834 const double Decal = 1.0 / NbrCotes;
835 const double Start = Decal / 2;
838 double posit = Start + ny*Decal;
839 Edge* edge = findEdge (v1, v2);
843 grid_geom->addAssociation (v1, subid, posit);
847 grid_geom->addAssociation (edge, subid, posit, posit+Decal);
851 grid_geom->addAssociation (edge, subid, posit, 1.0);
852 grid_geom->addAssociation (edge, subid, 0, Start);
855 // ===================================================== assoArc
856 void BiCylinder::assoArc (int cyl, int nx, int ny, int nz, int subid)
858 const double Decal = 1.0 / NbrCotes;
859 const double Start = Decal / 2;
862 double posit = Start + ny*Decal;
863 Edge* edge = getEdgeJ (cyl, nx, ny, nz);
868 Vertex* node = getVertexIJK (cyl, nx, ny, nz);
869 grid_geom->addAssociation (node, subid, posit);
873 grid_geom->addAssociation (edge, subid, posit, posit+Decal);
874 // node = getVertexIJK (cyl, nx, ny2, nz);
875 // grid_geom->addAssociation (node, subid, angle2*UnSur2pi);
879 grid_geom->addAssociation (edge, subid, posit, 1.0);
880 grid_geom->addAssociation (edge, subid, 0, Start);
883 // ===================================================== assoIntersection
884 int BiCylinder::assoIntersection (int nxs, int nzs)
887 Real3 pse, psw, sorig, sbase;
888 Real3 pbe, pbw, borig, bbase;
891 int nyd = NbrCotes/2;
892 int MiddleSlice1 = 3;
894 int nz = nzs < MiddleSlice1 ? 0 : 5;
896 Vertex* vse = getVertexIJK (CylSmall, nxs, ny0 , nz);
897 Vertex* vsw = getVertexIJK (CylSmall, nxs, nyd , nz);
898 Vertex* vbe = getVertexIJK (CylBig, nxs, ny0 , 0);
899 Vertex* vbw = getVertexIJK (CylBig, nxs, nyd , 0);
901 if (vse==NULL || vsw==NULL || vbe==NULL || vbw==NULL)
909 double srayon = calc_distance (pse, psw)/2;
910 double brayon = calc_distance (pbe, pbw)/2;
912 calc_milieu (psw, pse, sorig);
913 calc_milieu (pbw, pbe, borig);
914 calc_vecteur (psw, pse, sbase);
915 calc_vecteur (pbw, pbe, bbase);
917 double shaut = calc_distance (cross_center, sorig);
918 double bhaut = calc_distance (cross_center, borig)*2;
919 double* orig = nzs < MiddleSlice1 ? sorig : cross_center; // Pb orientation
921 BiCylinderShape bicyl_shape (el_root);
922 ier = bicyl_shape.defineCyls (borig, cross_dirbig, bbase, brayon, bhaut,
923 orig, cross_dirsmall, sbase, srayon, shaut);
927 for (int ny=0 ; ny<NbrCotes ; ny++)
929 Vertex* node = getVertexIJK (CylSmall, nxs, ny, nzs);
931 node->clearAssociation ();
932 // Edge* edge = getEdgeJ (CylSmall, nxs, ny, nzs);
933 // if (edge!=NULL) edge->clearAssociation ();
936 for (int ny=0 ; ny<NbrCotes ; ny++)
938 int ny1 = (ny+1) MODULO QUAD4;
939 Vertex* node0 = getVertexIJK (CylSmall, nxs, ny, nzs);
940 Vertex* node1 = getVertexIJK (CylSmall, nxs, ny1, nzs);
941 // Edge* edge = getEdgeJ (CylSmall, nxs, ny, nzs);
942 Edge* edge = findEdge (node0, node1);
943 bicyl_shape.associate (edge);
948 // ====================================================== makeCylinders
949 int BiCylinder::makeCylinders (Vertex* ori1, double* vz1, double r1, double h1,
950 Vertex* ori2, double* vz2, double r2, double h2)
952 grid_type = GR_BICYL;
953 int ier = makePipes (ori1,vz1, r1/2, r1,h1, ori2,vz2, r2/2, r2,h2, true);
956 // ====================================================== makePipes
957 int BiCylinder::makePipes (Vertex* ori1, double* vz1, double rint1,
958 double rext1, double h1, Vertex* ori2, double* vz2,
959 double rint2, double rext2, double h2, bool fill)
962 cross_hauteur [CylBig] = h2;
963 cross_rayext [CylBig ] = rext2;
964 cross_rayint [CylBig ] = rint2;
966 cross_hauteur [CylSmall] = h1;
967 cross_rayext [CylSmall] = rext1;
968 cross_rayint [CylSmall] = rint1;
970 ori1->getPoint (cross_orismall);
971 ori2->getPoint (cross_oribig );
972 copy_vecteur (vz2, cross_dirbig);
973 copy_vecteur (vz1, cross_dirsmall);
975 // Intersection des 2 cylindres :
976 int ier = checkInter (cross_oribig, vz2, rext2, h2,
977 cross_orismall, vz1, rext1, h1,
978 cross_center, at_left, at_right);
987 adjustLittleSlice (1, 1);
988 adjustLittleSlice (0, 2);
992 adjustLittleSlice (0, 1);
997 adjustLittleSlice (0, 3);
998 adjustLittleSlice (1, 4);
1002 adjustLittleSlice (0, 4);
1006 transfoVertices (cross_center, cross_dirsmall, cross_dirbig);
1009 // assoResiduelle ();
1010 // el_root->reorderQuads ();