TRUST 1.9.8
HPC thermohydraulic platform
Loading...
Searching...
No Matches
Octree.cpp
1/****************************************************************************
2* Copyright (c) 2026, CEA
3* All rights reserved.
4*
5* Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
6* 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
7* 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.
8* 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission.
9*
10* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
11* IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
12* OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
13*
14*****************************************************************************/
15
16#ifndef NbCasesParNoeuds
17#define NbCasesParNoeuds 512
18#endif
19
20#include <TRUSTVects.h>
21#include <Domaine.h>
22#include <Octree.h>
23
24Implemente_instanciable_sans_constructeur_32_64(OctreeRoot_32_64,"OctreeRoot",Objet_U);
25
26template <typename _SIZE_>
28{
29 int nb_octrees=Octree_32_64::nombre_d_octrees();
30 int i;
31 is << finl << "--------------------------------" << finl;
32 is << "niveau : " << niveau() << " taille : " << taille() << finl;
33 for(i=0; i<nb_octrees; i++)
34 {
35 if(les_octrees[i]!=0)
36 les_octrees[i]->printOn(is);
37 else
38 is << "vide" << " ";
39 }
40 return is << finl;
41}
42
43template <typename _SIZE_>
45{
47}
48
49template <typename _SIZE_>
51{
53}
54
55
56/*! @brief
57 *
58 * @param (OctreeLoc& loc)
59 * @param (int i)
60 */
61void OctreeLoc::construire(const OctreeLoc& loc, int i)
62{
63 double xmil=(loc.xmin+loc.xmax)/2;
64 double ymil=(loc.ymin+loc.ymax)/2;
65 double zmil=(loc.zmin+loc.zmax)/2;
66 construire(loc, i, xmil, ymil, zmil);
67}
68
69
70/*! @brief
71 *
72 * @param (OctreeLoc& loc)
73 * @param (int i)
74 * @param (double xmil)
75 * @param (double ymil)
76 * @param (double zmil)
77 */
78void OctreeLoc::construire(const OctreeLoc& loc, int i, double xmil, double ymil, double zmil)
79{
80 switch(i)
81 {
82 case GaucheArriereBas:
83 xmin=loc.xmin;
84 xmax=xmil;
85 ymin=loc.ymin;
86 ymax=ymil;
87 zmin=loc.zmin;
88 zmax=zmil;
89 break;
90 case GaucheArriereHaut:
91 xmin=loc.xmin;
92 xmax=xmil;
93 ymin=loc.ymin;
94 ymax=ymil;
95 zmin=zmil;
96 zmax=loc.zmax;
97 break;
98 case GaucheAvantBas:
99 xmin=loc.xmin;
100 xmax=xmil;
101 ymin=ymil;
102 ymax=loc.ymax;
103 zmin=loc.zmin;
104 zmax=zmil;
105 break;
106 case GaucheAvantHaut:
107 xmin=loc.xmin;
108 xmax=xmil;
109 ymin=ymil;
110 ymax=loc.ymax;
111 zmin=zmil;
112 zmax=loc.zmax;
113 break;
114 case DroitArriereBas:
115 xmin=xmil;
116 xmax=loc.xmax;
117 ymin=loc.ymin;
118 ymax=ymil;
119 zmin=loc.zmin;
120 zmax=zmil;
121 break;
122 case DroitArriereHaut:
123 xmin=xmil;
124 xmax=loc.xmax;
125 ymin=loc.ymin;
126 ymax=ymil;
127 zmin=zmil;
128 zmax=loc.zmax;
129 break;
130 case DroitAvantBas:
131 xmin=xmil;
132 xmax=loc.xmax;
133 ymin=ymil;
134 ymax=loc.ymax;
135 zmin=loc.zmin;
136 zmax=zmil;
137 break;
138 case DroitAvantHaut:
139 xmin=xmil;
140 xmax=loc.xmax;
141 ymin=ymil;
142 ymax=loc.ymax;
143 zmin=zmil;
144 zmax=loc.zmax;
145 break;
146 }
147}
148
149
150/*! @brief
151 *
152 * @param (double x)
153 * @param (double y)
154 * @param (double z)
155 * @return (int)
156 */
157int OctreeLoc::direction(double x, double y, double z) const
158{
159 double xmil=(xmin+xmax)/2;
160 double ymil=(ymin+ymax)/2;
161 double zmil=(zmin+zmax)/2;
162
163 if( (x<xmil) && (y<ymil) && (z<zmil))
164 return GaucheArriereBas;
165 if( (x<xmil) && (y<ymil) && !(z<zmil))
166 return GaucheArriereHaut;
167 if( (x<xmil) && !(y<ymil) && (z<zmil))
168 return GaucheAvantBas;
169 if( (x<xmil) && !(y<ymil) && !(z<zmil))
170 return GaucheAvantHaut;
171 if( !(x<xmil) && (y<ymil) && (z<zmil))
172 return DroitArriereBas;
173 if( !(x<xmil) && (y<ymil) && !(z<zmil))
174 return DroitArriereHaut;
175 if( !(x<xmil) && !(y<ymil) && (z<zmil))
176 return DroitAvantBas;
177 if( !(x<xmil) && !(y<ymil) && !(z<zmil))
178 return DroitAvantHaut;
179
180 Cerr << "Error in OctreeLoc::direction" << finl;
182 return -1;
183}
184
185
186/*! @brief
187 *
188 * @param (OctreeLoc& loc)
189 * @param (double x)
190 * @param (double y)
191 * @param (double z)
192 * @return (int)
193 */
194template <typename _SIZE_>
195typename Octree_32_64<_SIZE_>::int_t Octree_32_64<_SIZE_>::rang_elem_loc(const OctreeLoc& loc, double x, double y, double z) const
196{
197 int_t element=-1;
198 int j=1;
199 while(j<nombre_d_octrees())
200 {
201 if(les_octrees[j]!=0)
202 break;
203 else
204 j++;
205 }
206 if(j<nombre_d_octrees())
207 {
208 int i=loc.direction(x, y, z);
209 OctreeLoc loci;
210 loci.construire(loc, i);
211 if(les_octrees[i]==0)
212 return -1;
213 element = les_octrees[i]->rang_elem_loc(loci, x, y, z);
214 }
215 else
216 element = les_octrees[0]->rang_elem_loc(loc, x, y, z);
217 return element;
218}
219
220
221/*! @brief
222 *
223 * @param (OctreeLoc& loc)
224 * @param (int prems)
225 * @param (double x)
226 * @param (double y)
227 * @param (double z)
228 * @return (int)
229 */
230template <typename _SIZE_>
231typename Octree_32_64<_SIZE_>::int_t Octree_32_64<_SIZE_>::rang_elem_depuis_loc(const OctreeLoc& loc, int_t prems, double x, double y, double z) const
232{
233 int_t element=-1;
234 int j=1;
235 while(j<nombre_d_octrees())
236 {
237 if(les_octrees[j]!=0)
238 break;
239 else
240 j++;
241 }
242 if(j<nombre_d_octrees())
243 {
244 int i=loc.direction(x, y, z);
245 OctreeLoc loci;
246 loci.construire(loc, i);
247 if(les_octrees[i]==0)
248 return -1;
249 element = les_octrees[i]->rang_elem_depuis_loc(loci, prems, x, y, z);
250 }
251 else
252 element = les_octrees[0]->rang_elem_depuis_loc(loc, prems, x, y, z);
253 return element;
254}
255
256
257template <typename _SIZE_>
258void Octree_32_64<_SIZE_>::ranger_elem_1D(ArrOfInt& ok, int_t elem, int_t i, int nb_som_elem, const DoubleTab_t& coord, const IntTab_t& les_elems,
259 SmallArrOfTID_t& compteur, Vect_IntTab_t& SousTab, double xmil)
260{
261 for(int som=0; som<nb_som_elem; som++)
262 {
263 int_t sommet=les_elems(i,som);
264 assert(sommet>=0);
265 if(coord(sommet,0)<xmil)
266 {
267 if(!ok[GaucheArriereBas])
268 {
269 SousTab[GaucheArriereBas][compteur[GaucheArriereBas]]=i;
270 compteur[GaucheArriereBas]++;
271 ok[GaucheArriereBas]=1;
272 }
273 }
274 else
275 {
276 if(!ok[DroitArriereBas])
277 {
278 SousTab[DroitArriereBas][compteur[DroitArriereBas]]=i;
279 compteur[DroitArriereBas]++;
280 ok[DroitArriereBas]=1;
281 }
282 }
283 }
284}
285
286namespace // Anonymous namespace
287{
288
289template <typename _SIZE_>
290inline void range2D(double x, double y, double xmil, double ymil,
291 ArrOfInt& ok, TRUST_Vector<IntTab_T<_SIZE_>>& SousTab, SmallArrOfTID_T<_SIZE_>& compteur, _SIZE_ i)
292{
293 if(x<xmil)
294 {
295 if(y<ymil)
296 {
297 if(!ok[GaucheArriereBas])
298 {
299 SousTab[GaucheArriereBas][compteur[GaucheArriereBas]]=i;
300 compteur[GaucheArriereBas]++;
301 ok[GaucheArriereBas]=1;
302 }
303 }
304 else
305 {
306 if(!ok[GaucheAvantBas])
307 {
308 SousTab[GaucheAvantBas][compteur[GaucheAvantBas]]=i;
309 compteur[GaucheAvantBas]++;
310 ok[GaucheAvantBas]=1;
311 }
312 }
313 }
314 else
315 {
316 if(y<ymil)
317 {
318 if(!ok[DroitArriereBas])
319 {
320 SousTab[DroitArriereBas][compteur[DroitArriereBas]]=i;
321 compteur[DroitArriereBas]++;
322 ok[DroitArriereBas]=1;
323 }
324 }
325 else
326 {
327 if(!ok[DroitAvantBas])
328 {
329 SousTab[DroitAvantBas][compteur[DroitAvantBas]]=i;
330 compteur[DroitAvantBas]++;
331 ok[DroitAvantBas]=1;
332 }
333 }
334 }
335}
336
337template <typename _SIZE_>
338inline void range3D(double x, double y, double z,
339 double xmil, double ymil, double zmil,
340 ArrOfInt& ok, TRUST_Vector<IntTab_T<_SIZE_>>& SousTab, SmallArrOfTID_T<_SIZE_>& compteur, _SIZE_ i)
341{
342 if(x<xmil)
343 {
344 if(y<ymil)
345 {
346 if(z<zmil)
347 {
348 if(!ok[GaucheArriereBas])
349 {
350 SousTab[GaucheArriereBas][compteur[GaucheArriereBas]]=i;
351 compteur[GaucheArriereBas]++;
352 ok[GaucheArriereBas]=1;
353 }
354 }
355 else
356 {
357 if(!ok[GaucheArriereHaut])
358 {
359 SousTab[GaucheArriereHaut][compteur[GaucheArriereHaut]]=i;
360 compteur[GaucheArriereHaut]++;
361 ok[GaucheArriereHaut]=1;
362 }
363 }
364 }
365 else
366 {
367 if(z<zmil)
368 {
369 if(!ok[GaucheAvantBas])
370 {
371 SousTab[GaucheAvantBas][compteur[GaucheAvantBas]]=i;
372 compteur[GaucheAvantBas]++;
373 ok[GaucheAvantBas]=1;
374 }
375 }
376 else
377 {
378 if(!ok[GaucheAvantHaut])
379 {
380 SousTab[GaucheAvantHaut][compteur[GaucheAvantHaut]]=i;
381 compteur[GaucheAvantHaut]++;
382 ok[GaucheAvantHaut]=1;
383 }
384 }
385 }
386 }
387 else
388 {
389 if(y<ymil)
390 {
391 if(z<zmil)
392 {
393 if(!ok[DroitArriereBas])
394 {
395 SousTab[DroitArriereBas][compteur[DroitArriereBas]]=i;
396 compteur[DroitArriereBas]++;
397 ok[DroitArriereBas]=1;
398 }
399 }
400 else
401 {
402 if(!ok[DroitArriereHaut])
403 {
404 SousTab[DroitArriereHaut][compteur[DroitArriereHaut]]=i;
405 compteur[DroitArriereHaut]++;
406 ok[DroitArriereHaut]=1;
407 }
408 }
409 }
410 else
411 {
412 if(z<zmil)
413 {
414 if(!ok[DroitAvantBas])
415 {
416 SousTab[DroitAvantBas][compteur[DroitAvantBas]]=i;
417 compteur[DroitAvantBas]++;
418 ok[DroitAvantBas]=1;
419 }
420 }
421 else
422 {
423 if(!ok[DroitAvantHaut])
424 {
425 SousTab[DroitAvantHaut][compteur[DroitAvantHaut]]=i;
426 compteur[DroitAvantHaut]++;
427 ok[DroitAvantHaut]=1;
428 }
429 }
430 }
431 }
432}
433
434}// end anonymous namespace
435
436
437template <typename _SIZE_>
438void Octree_32_64<_SIZE_>::ranger_elem_2D(ArrOfInt& ok, int_t elem, int_t idx, int nb_som_elem, const DoubleTab_t& coord, const IntTab_t& elems,
439 SmallArrOfTID_t& compteur, Vect_IntTab_t& SousTab, double xmil, double ymil)
440{
441 double xmin=DMAXFLOAT, xmax=-DMAXFLOAT;
442 double ymin=xmin, ymax=xmax;
443 for(int som=0; som<nb_som_elem; som++)
444 {
445 int_t sommet=elems(idx,som);
446// assert(sommet>=0);
447// GF dans le cas dde polyedre le tableau elem est surdimensionne et peut contenir des -1
448 if (sommet>=0)
449 {
450 xmin=std::min(xmin, coord(sommet,0));
451 xmax=std::max(xmax, coord(sommet,0));
452 ymin=std::min(ymin, coord(sommet,1));
453 ymax=std::max(ymax, coord(sommet,1));
454 }
455 }
456 range2D(xmin, ymin, xmil, ymil, ok, SousTab, compteur, idx);
457 range2D(xmin, ymax, xmil, ymil, ok, SousTab, compteur, idx);
458 range2D(xmax, ymin, xmil, ymil, ok, SousTab, compteur, idx);
459 range2D(xmax, ymax, xmil, ymil, ok, SousTab, compteur, idx);
460}
461
462template <typename _SIZE_>
463void Octree_32_64<_SIZE_>::ranger_elem_3D(ArrOfInt& ok, int_t elem, int_t idx, int nb_som_elem, const DoubleTab_t& coord, const IntTab_t& elems,
464 SmallArrOfTID_t& compteur, Vect_IntTab_t& SousTab, double xmil, double ymil, double zmil)
465{
466 double epsilon=get_epsilon();
467 double xmin=DMAXFLOAT, xmax=-DMAXFLOAT;
468 double ymin=xmin, ymax=xmax;
469 double zmin=xmin, zmax=xmax;
470 for(int som=0; som<nb_som_elem; som++)
471 {
472 int_t sommet=elems(idx,som);
473// assert(sommet>=0);
474// GF dans le cas dde polyedre le tableau elem est surdimensionne et peut contenir des -1
475 if (sommet>=0)
476 {
477 xmin=std::min(xmin, coord(sommet,0))-epsilon;
478 xmax=std::max(xmax, coord(sommet,0))+epsilon;
479 ymin=std::min(ymin, coord(sommet,1))-epsilon;
480 ymax=std::max(ymax, coord(sommet,1))+epsilon;
481 zmin=std::min(zmin, coord(sommet,2))-epsilon;
482 zmax=std::max(zmax, coord(sommet,2))+epsilon;
483 }
484 }
485 range3D(xmin, ymin, zmin, xmil, ymil, zmil, ok, SousTab, compteur, idx);
486 range3D(xmin, ymax, zmin, xmil, ymil, zmil, ok, SousTab, compteur, idx);
487 range3D(xmax, ymin, zmin, xmil, ymil, zmil, ok, SousTab, compteur, idx);
488 range3D(xmax, ymax, zmin, xmil, ymil, zmil, ok, SousTab, compteur, idx);
489 range3D(xmin, ymin, zmax, xmil, ymil, zmil, ok, SousTab, compteur, idx);
490 range3D(xmin, ymax, zmax, xmil, ymil, zmil, ok, SousTab, compteur, idx);
491 range3D(xmax, ymin, zmax, xmil, ymil, zmil, ok, SousTab, compteur, idx);
492 range3D(xmax, ymax, zmax, xmil, ymil, zmil, ok, SousTab, compteur, idx);
493}
494
495/*! @brief
496 *
497 * @param (int nb_octrees)
498 * @param (ArrOfInt& Tab)
499 * @param (OctreeLoc& loc)
500 * @param (Octree* pe)
501 */
502template <typename _SIZE_>
503void Octree_32_64<_SIZE_>::construire(int nb_octrees, const ArrOfInt_t& Tab, const OctreeLoc& loc, Octree_t* pe)
504{
505 static int level=0;
506 int_t nb_elem=Tab.size_array();
507 pere=pe;
508 les_octrees=new Octree_32_64*[nb_octrees];
509
510 if(niveau()>level)
511 {
512 level=niveau();
513 Process::Journal() << "Octree: level=" << level
514 << " nb_elem=" << nb_elem << finl;
515 if(level > 100) Process::exit();
516 }
517
518 if(nb_elem<=NbCasesParNoeuds)
519 {
520 for(int i=0; i<nb_octrees; i++)
521 les_octrees[i]=0;
522 les_octrees[0]=new OctreeFloor_32_64<_SIZE_>(this, Tab, loc);
523 }
524 else
525 {
526 //else
527 const Domaine_t& dom=domaine();
528 Vect_IntTab_t SousTab(nb_octrees);
529 SmallArrOfTID_t compteur(nb_octrees);
530 int nb_som_elem=dom.nb_som_elem();
531 const IntTab_t& les_elems=dom.les_elems();
532 const DoubleTab_t& coord_sommets=dom.coord_sommets();
533 for(int i=0; i<nb_octrees; i++)
534 SousTab[i].resize(nb_elem);
535 double xmil=(loc.xmin+loc.xmax)/2;
536 double ymil=(loc.ymin+loc.ymax)/2;
537 double zmil=(loc.zmin+loc.zmax)/2;
538 ArrOfInt ok(nb_octrees);
539 for(int_t elem=0; elem<nb_elem; elem++)
540 {
541 int_t idx=Tab[elem];
542 ok=0;
543 switch(Objet_U::dimension)
544 {
545 case 1:
546 // ArrOfInt oks, int_t elem, int_t idx, int nb_som_elem, DoubleTab_t coords_som, IntTab_t les_elems,
547 // SmallArrOfTID_t compteur, ArrOfInt_t sous_tab, double xmil
548 ranger_elem_1D(ok, elem, idx , nb_som_elem, coord_sommets,
549 les_elems, compteur, SousTab, xmil);
550 break;
551 case 2:
552 ranger_elem_2D(ok, elem, idx , nb_som_elem, coord_sommets,
553 les_elems, compteur, SousTab, xmil, ymil);
554 break;
555 case 3:
556 ranger_elem_3D(ok, elem, idx ,nb_som_elem, coord_sommets,
557 les_elems, compteur, SousTab, xmil, ymil, zmil);
558 break;
559 default:
560 Cerr << "Space dimension " << Objet_U::dimension << " incorrect" << finl;
562 }
563 }
564
565 for(int i=0; i<nb_octrees; i++)
566 {
567 les_octrees[i]=0;
568 OctreeLoc loci;
569 loci.construire(loc, i, xmil, ymil, zmil);
570 SousTab[i].resize(compteur[i]);
571 if(compteur[i]==nb_elem)
572 les_octrees[i]=new OctreeFloor_32_64<_SIZE_>(this, SousTab[i], loci);
573 else if(compteur[i]>NbCasesParNoeuds)
574 les_octrees[i]=new Octree_32_64<_SIZE_>(nb_octrees, this, SousTab[i], loci);
575 else if(compteur[i]!=0)
576 les_octrees[i]=new OctreeFloor_32_64<_SIZE_>(this, SousTab[i], loci);
577 }
578 }
579}
580
581
582/*! @brief Detruit l'octree.
583 *
584 * Methode appelee par le destructeur
585 *
586 */
587template <typename _SIZE_>
589{
590 if(les_octrees==0)
591 return;
593 for(int i=0; i<nb_octrees; i++)
594 if(les_octrees[i])
595 delete les_octrees[i];
596 delete[] les_octrees;
597 les_octrees=0;
598}
599
600
601template <typename _SIZE_>
603{
604 int i=0;
605 if (pere==0) return i;
606 return (pere->niveau()+1);
607}
608
609
610/*! @brief Renvoie la taille de l'octree.
611 *
612 * @return (unsigned) la taille de l'octree
613 */
614template <typename _SIZE_>
616{
618 int_t la_taille = (int_t) sizeof(Octree);
619 for(int i=0; i<nb_octrees; i++)
620 if(les_octrees[i]!=0)
621 la_taille+= les_octrees[i]->taille();
622 return la_taille;
623}
624
625
626template <typename _SIZE_>
628{
629 this->detruire();
630 valid_=1;
631 reel_=reel_prec;
632 pere=0;
633 const Domaine_t& dom = domaine();
634
635 {
636 // Calcul du min et du max des coordonnees dans chaque direction
637 double min[3] = { 1e30, 1e30, 1e30 };
638 double max[3] = {-1e30,-1e30,-1e30 };
639 const DoubleTab_t& tab = dom.coord_sommets();
640 // tot sinon bouding box trop petite
641 // peut etre teste reel_ ...
642 const int_t n = tab.dimension_tot(0);
643 const int dim = tab.dimension_int(1);
644 assert(dim <= 3);
645 for (int_t i = 0; i < n; i++)
646 {
647 for (int j = 0; j < 3; j++)
648 {
649 if (j>=dim) break;
650 double x = tab(i,j);
651 if (x < min[j])
652 min[j] = x;
653 if (x > max[j])
654 max[j] = x;
655 }
656 }
657 loc.xmin = min[0];
658 loc.xmax = max[0];
659 if (axi)
660 {
661 loc.ymin = 0.;
662 loc.ymax = 2. * M_PI;
663 }
664 else
665 {
666 loc.ymin = min[1];
667 loc.ymax = max[1];
668 }
669 if (dim > 2)
670 {
671 loc.zmin = min[2];
672 loc.zmax = max[2];
673 }
674 else
675 {
676 loc.zmin = 0.;
677 loc.zmax = 1.;
678 }
679 }
680 int_t nb_elem=dom.nb_elem_tot();
681 if(reel_prec)
682 nb_elem=dom.nb_elem();
683
684 ArrOfInt_t SousTab(nb_elem);
685 for(int_t i=0; i<nb_elem; i++)
686 SousTab[i]=i;
688 Cerr << "Construction of the OctreeRoot OK " << finl;
689}
690
691
692/*! @brief
693 *
694 * @param (double x)
695 * @param (double y)
696 * @param (double z)
697 * @return (int)
698 * @throws Erreur dans OctreeRoot_32_64<_SIZE_>::rang_sommet
699 * @throws On a pas trouve de sommet a ces coordonnees
700 */
701template <typename _SIZE_>
703{
704 int_t elem=rang_elem(x, y, z);
705 if (elem != -1)
706 {
707 const Domaine_t& dom=domaine();
708 int nb_som_elem=dom.nb_som_elem();
709 double epsilon=this->get_epsilon();
710 for(int i=0; i<nb_som_elem; i++)
711 {
712 //for i
713 int_t sommet=dom.sommet_elem(elem,i);
714 if(sommet<0)
715 {
716 Cerr << "This is not a mesh !" << finl;
717 Cerr << dom << finl;
719 }
720 switch(Objet_U::dimension)
721 {
722 //switch
723 case 1:
724 if(std::fabs(dom.coord(sommet,0)-x)<epsilon)
725 return sommet;
726 break;
727 case 2:
728 if((std::fabs(dom.coord(sommet,0)-x)<epsilon) &&
729 (std::fabs(dom.coord(sommet,1)-y)<epsilon))
730 return sommet;
731 break;
732 case 3:
733 if((std::fabs(dom.coord(sommet,0)-x)<epsilon) &&
734 (std::fabs(dom.coord(sommet,1)-y)<epsilon) &&
735 (std::fabs(dom.coord(sommet,2)-z)<epsilon))
736 return sommet;
737 break;
738 default:
739 Cerr << "Error in OctreeRoot_32_64<_SIZE_>::rang_sommet" << finl;
741 return -1;
742 }
743 }
744 return -1;
745 }
746 else
747 return -1;
748}
749
750/*! @brief
751 *
752 * @param (double x)
753 * @param (double y)
754 * @param (double z)
755 * @return (int)
756 * @throws Erreur dans OctreeRoot_32_64<_SIZE_>::rang_arete
757 * @throws On a pas trouve d'arete a ces coordonnees
758 */
759template <typename _SIZE_>
761{
762 // Recherche l'element contenant le point x,y,z:
763 int_t elem=rang_elem(x, y, z);
764 const IntTab_t& elem_aretes=domaine().elem_aretes();
765 if (elem>=elem_aretes.dimension(0)) return -2;
766 if (elem != -1 && elem<elem_aretes.dimension(0))
767 {
768 // Boucle sur les aretes de l'element
769 double epsilon=this->get_epsilon();
770 const IntTab_t& aretes_som=domaine().aretes_som();
771 const DoubleTab_t& coord=domaine().coord_sommets();
772 int nb_aretes_elem = elem_aretes.dimension_int(1);
773 for(int i=0; i<nb_aretes_elem; i++)
774 {
775 // On boucle sur les aretes de l'element pour trouver celle
776 // dont le centre de gravite coincide avec le point x,y,z
777 int_t arete=elem_aretes(elem,i);
778 int_t s0=aretes_som(arete,0);
779 int_t s1=aretes_som(arete,1);
780 switch(Objet_U::dimension)
781 {
782 case 3:
783 if((std::fabs(0.5*(coord(s0,0)+coord(s1,0))-x)<epsilon) &&
784 (std::fabs(0.5*(coord(s0,1)+coord(s1,1))-y)<epsilon) &&
785 (std::fabs(0.5*(coord(s0,2)+coord(s1,2))-z)<epsilon))
786 return arete;
787 break;
788 default:
789 Cerr << "Error in OctreeRoot_32_64<_SIZE_>::rang_elem" << finl;
791 return -1;
792 }
793 }
794 return -1;
795 }
796 else
797 return -1;
798}
799
800/*! @brief
801 *
802 * @param (double x)
803 * @param (double y)
804 * @param (double z)
805 * @return (int)
806 */
807template <typename _SIZE_>
809{
810 double epsilon = this->get_epsilon();
811 if( (x<loc.xmin-epsilon) || (y<loc.ymin-epsilon) || (z<loc.zmin-epsilon)
812 || (x>loc.xmax+epsilon) || (y>loc.ymax+epsilon) || (z>loc.zmax+epsilon) )
813 return -1;
815}
816
817
818/*! @brief
819 *
820 * @param (int prems)
821 * @param (double x)
822 * @param (double y)
823 * @param (double z)
824 */
825template <typename _SIZE_>
827{
828 double epsilon = this->get_epsilon();
829 if( (x<loc.xmin-epsilon) || (y<loc.ymin-epsilon) || (z<loc.zmin-epsilon)
830 || (x>loc.xmax+epsilon) || (y>loc.ymax+epsilon) || (z>loc.zmax+epsilon) )
831 return -1;
832 return Octree_32_64<_SIZE_>::rang_elem_depuis_loc(loc, prems, x, y, z);
833}
834
835
836/*! @brief
837 *
838 * @param (DoubleTab& positions)
839 * @param (ArrOfInt& sommets)
840 * @return (int)
841 * @throws dimension d'espace non prevue
842 */
843template <typename _SIZE_>
845{
846 int sz=positions.dimension(0);
847 assert(positions.dimension(1)==Objet_U::dimension);
848 sommets.resize_array(sz);
849 for(int i=0; i<sz; i++)
850 {
851 double x=0, y=0, z=0;
852 switch(Objet_U::dimension)
853 {
854 case 1:
855 x=positions(i,0);
856 break;
857 case 2:
858 x=positions(i,0);
859 y=positions(i,1);
860 break;
861 case 3:
862 x=positions(i,0);
863 y=positions(i,1);
864 z=positions(i,2);
865 break;
866 default:
867 Cerr << "Error in OctreeRoot_32_64<_SIZE_>::rang_sommet" << finl;
869 }
870 sommets[i]=rang_sommet(x,y,z);
871 }
872 return 1;
873}
874
875/*! @brief
876 *
877 * @param (IntTab& elem_aretes) Definition des aretes par leur sommet
878 * @param (DoubleTab& positions)
879 * @param (ArrOfInt& aretes) Les aretes trouvees
880 * @return (int)
881 * @throws dimension d'espace non prevue
882 */
883template <typename _SIZE_>
885{
886 int sz=positions.dimension(0);
887 assert(positions.dimension(1)==Objet_U::dimension);
888 aretes.resize_array(sz);
889 for(int i=0; i<sz; i++)
890 {
891 double x=0, y=0, z=0;
892 switch(Objet_U::dimension)
893 {
894 case 1:
895 x=positions(i,0);
896 break;
897 case 2:
898 x=positions(i,0);
899 y=positions(i,1);
900 break;
901 case 3:
902 x=positions(i,0);
903 y=positions(i,1);
904 z=positions(i,2);
905 break;
906 default:
907 Cerr << "Error in OctreeRoot_32_64<_SIZE_>::rang_arete" << finl;
909 }
910 aretes[i]=rang_arete(x,y,z);
911 }
912 return 1;
913}
914
915/*! @brief
916 *
917 * @param (DoubleTab& positions)
918 * @param (ArrOfInt& elems)
919 * @return (int)
920 * @throws dimension d'espace non prevue
921 */
922template <typename _SIZE_>
924{
925 int sz=positions.dimension(0);
926 assert(positions.dimension(1)==Objet_U::dimension);
927 elems.resize_array(sz);
928 for(int i=0; i<sz; i++)
929 {
930 double x=0, y=0, z=0;
931 switch(Objet_U::dimension)
932 {
933 case 1:
934 x=positions(i,0);
935 break;
936 case 2:
937 x=positions(i,0);
938 y=positions(i,1);
939 break;
940 case 3:
941 x=positions(i,0);
942 y=positions(i,1);
943 z=positions(i,2);
944 break;
945 default:
946 Cerr << "Error in OctreeRoot_32_64<_SIZE_>::rang_elem" << finl;
948 }
949 elems[i]=rang_elem(x,y,z);
950 }
951 return 1;
952}
953
954
955/*! @brief
956 *
957 * @param (DoubleTab& positions)
958 * @param (ArrOfInt& prems)
959 * @param (ArrOfInt& elems)
960 * @throws dimension d'espace non prevue
961 */
962template <typename _SIZE_>
964{
965 int sz=positions.dimension(0);
966 assert(positions.dimension(1)==Objet_U::dimension);
967 elems.resize_array(sz);
968 for(int i=0; i<sz; i++)
969 {
970 double x=0, y=0, z=0;
971 switch(Objet_U::dimension)
972 {
973 case 1:
974 x=positions(i,0);
975 break;
976 case 2:
977 x=positions(i,0);
978 y=positions(i,1);
979 break;
980 case 3:
981 x=positions(i,0);
982 y=positions(i,1);
983 z=positions(i,2);
984 break;
985 default:
986 Cerr << "Error in OctreeRoot_32_64<_SIZE_>::rang_elem_depuis" << finl;
988 }
989 elems[i]=rang_elem_depuis(prems[i],x,y,z);
990 }
991 return 1;
992}
993
994
995/*! @brief
996 *
997 * @param (ArrOfInt& elements) we suppose we look for a small number of elements!
998 * @param (double x)
999 * @param (double y)
1000 * @param (double z)
1001 */
1002template <typename _SIZE_>
1003void OctreeRoot_32_64<_SIZE_>::rang_elems_sommet(SmallArrOfTID_t& elements, double x, double y, double z) const
1004{
1005 std::vector<int_t> els;
1006 bool fini=false;
1007 int i = 0;
1008 int_t el=0;
1009 while(!fini)
1010 {
1012 if(el==-1)
1013 fini=1;
1014 else
1015 {
1016 els.push_back(el);
1017 i++;
1018 el++;
1019 }
1020 }
1021 elements.resize_array(i);
1022 for(int j=0; j<i; j++)
1023 elements[j] = els[j];
1024}
1025
1026
1027/*! @brief Renvoie vrai si le domaine associe a l'octree est non nulle.
1028 *
1029 * @return (int) code de retour propage
1030 */
1031template <typename _SIZE_>
1033{
1034 if((!le_dom)||(valid_!=1))
1035 // L'Octree n'est pas construit ou le domaine est nulle
1036 return 0;
1037 else
1038 return 1;
1039}
1040
1041template <typename _SIZE_>
1043{
1044 if(valid_!=0)
1045 {
1046 Cerr << "Invalidation of the Octree" << finl;
1047 valid_ = 0;
1048 this->detruire();
1049 }
1050}
1051
1052
1053/*! @brief
1054 *
1055 * @param (OctreeLoc& loc)
1056 * @param (double x)
1057 * @param (double y)
1058 * @param (double z)
1059 * @return (int)
1060 */
1061template <typename _SIZE_>
1062typename OctreeFloor_32_64<_SIZE_>::int_t OctreeFloor_32_64<_SIZE_>::rang_elem_loc(const OctreeLoc& loc, double x, double y, double z) const
1063{
1064 int_t sz=num_elem.size_array();
1065 int_t element=-1;
1066 const OWN_PTR(Elem_geom_base_32_64<_SIZE_>)& elemgeom=this->domaine().type_elem();
1067 pos[0]=x;
1068 pos[1]=y;
1069 if(Objet_U::dimension>2) pos[2]=z;
1070 // ToDo Kokkos: exemple de virtual function facile a implementer:
1071 for (int_t ielem = 0; ielem < sz; ielem++)
1072 {
1073 if(elemgeom->contient(pos,num_elem[ielem]))
1074 {
1075 element=num_elem[ielem];
1076 break;
1077 }
1078 else
1079 {
1080 // Process::Journal() << "PAS TROUVE !!: " << finl;
1081 }
1082 }
1083 return element;
1084}
1085
1086/*! @brief
1087 *
1088 * @param (OctreeLoc& loc)
1089 * @param (int prems)
1090 * @param (double x)
1091 * @param (double y)
1092 * @param (double z)
1093 * @return (int)
1094 */
1095template <typename _SIZE_>
1096typename OctreeFloor_32_64<_SIZE_>::int_t OctreeFloor_32_64<_SIZE_>::rang_elem_depuis_loc(const OctreeLoc& loc, int_t prems, double x, double y, double z) const
1097{
1098 int_t sz=num_elem.size_array();
1099 int_t element=-1;
1100 const OWN_PTR(Elem_geom_base_32_64<_SIZE_>)& elemgeom = this->domaine().type_elem();
1101 pos[0]=x;
1102 pos[1]=y;
1103 if(Objet_U::dimension>2) pos[2]=z;
1104 bool trouve=false;
1105 int_t ielem=0;
1106 while((ielem < sz) && !trouve)
1107 {
1108 if(elemgeom->contient(pos,num_elem[ielem]))
1109 {
1110 element=num_elem[ielem];
1111 if(element >= prems)
1112 trouve=true;
1113 else
1114 ielem++;
1115 }
1116 else
1117 ielem++;
1118 }
1119 if(ielem>=sz)
1120 element=-1;
1121 return element;
1122}
1123
1124
1125/*! @brief
1126 *
1127 * @param (Octree* pe)
1128 * @param (ArrOfInt& Tab)
1129 * @param (OctreeLoc& loc)
1130 */
1131template <typename _SIZE_>
1133{
1134 assert(pe!=0);
1135 pere=pe;
1136 les_octrees=0;
1137 num_elem=Tab;
1138 // Ajout B.M: tri des elements dans l'ordre croissant de l'indice local.
1139 // Ainsi, lors des recherches de sommets ou faces de joint, on tombe
1140 // sur l'element reel, pas l'element virtuel
1141 num_elem.ordonne_array();
1142}
1143
1144
1145/*! @brief Renvoie la taille de l'octree.
1146 *
1147 * @return (unsigned) la taille de l'octree
1148 */
1149template <typename _SIZE_>
1151{
1152 return (int_t)sizeof(OctreeFloor_32_64)+num_elem.size_array()*(int_t)sizeof(int);
1153}
1154
1155
1156template class Octree_32_64<int>;
1157template class OctreeRoot_32_64<int>;
1158template class OctreeFloor_32_64<int>;
1159#if INT_is_64_ == 2
1160template class Octree_32_64<trustIdType>;
1161template class OctreeRoot_32_64<trustIdType>;
1162template class OctreeFloor_32_64<trustIdType>;
1163#endif
1164
int nb_som_elem() const
Renvoie le nombre de sommets des elements geometriques constituants le domaine.
Definition Domaine.h:474
int_t nb_elem_tot() const
Definition Domaine.h:132
IntTab_t & les_elems()
Definition Domaine.h:129
int_t nb_elem() const
Definition Domaine.h:131
const DoubleTab_t & coord_sommets() const
Definition Domaine.h:112
double coord(int_t i, int j) const
Definition Domaine.h:110
int_t sommet_elem(int_t i, int j) const
Renvoie le numero (global) du j-ieme sommet du i-ieme element.
Definition Domaine.h:136
Classe Elem_geom_base Cette classe est la classe de base pour la definition d'elements.
Class defining operators and methods for all reading operation in an input flow (file,...
Definition Entree.h:42
classe Objet_U Cette classe est la classe de base des Objets de TRUST
Definition Objet_U.h:73
static int dimension
Definition Objet_U.h:99
static int axi
Definition Objet_U.h:101
Classe OctreeFloor.
Definition Octree.h:164
Octree_32_64 ** les_octrees
Definition Octree.h:91
int_t rang_elem_depuis_loc(const OctreeLoc &, int_t prems, double x, double y=0, double z=0) const override
Definition Octree.cpp:1096
Octree_32_64< _SIZE_ > Octree_t
Definition Octree.h:167
ArrOfDouble pos
Definition Octree.h:189
ArrOfInt_T< _SIZE_ > ArrOfInt_t
Definition Octree.h:168
OctreeFloor_32_64(Octree_t *mon_pere, const ArrOfInt_t &val, const OctreeLoc &loc)
Definition Octree.h:170
void construire(Octree_t *, const ArrOfInt_t &, const OctreeLoc &)
Definition Octree.cpp:1132
int_t taille() const override
Renvoie la taille de l'octree.
Definition Octree.cpp:1150
Octree_32_64 * pere
Definition Octree.h:92
int_t rang_elem_loc(const OctreeLoc &, double x, double y=0, double z=0) const override
Definition Octree.cpp:1062
ArrOfInt_t num_elem
Definition Octree.h:188
Classe OctreeRoot.
Definition Octree.h:102
int_t rang_sommet(double x, double y=0, double z=0) const
Definition Octree.cpp:702
SmallArrOfTID_T< _SIZE_ > SmallArrOfTID_t
Definition Octree.h:109
int_t rang_elem(double x, double y=0, double z=0) const
Definition Octree.cpp:808
int_t rang_elem_depuis(int_t prems, double x, double y=0, double z=0) const
Definition Octree.cpp:826
IntTab_T< _SIZE_ > IntTab_t
Definition Octree.h:108
void construire(int reel=0)
Definition Octree.cpp:627
void rang_elems_sommet(SmallArrOfTID_t &, double x, double y=0, double z=0) const
Definition Octree.cpp:1003
int construit() const
Renvoie vrai si le domaine associe a l'octree est non nulle.
Definition Octree.cpp:1032
ArrOfInt_T< _SIZE_ > ArrOfInt_t
Definition Octree.h:107
OctreeLoc loc
Definition Octree.h:154
int_t rang_arete(double x, double y=0, double z=0) const
Definition Octree.cpp:760
Octree_32_64 * pere
Definition Octree.h:92
const Domaine_t & domaine() const override
Definition Octree.h:142
DoubleTab_T< _SIZE_ > DoubleTab_t
Definition Octree.h:110
Domaine_32_64< _SIZE_ > Domaine_t
Definition Octree.h:111
Classe Octree.
Definition Octree.h:54
Octree_32_64 ** les_octrees
Definition Octree.h:91
void detruire()
Detruit l'octree.
Definition Octree.cpp:588
virtual int_t rang_elem_depuis_loc(const OctreeLoc &, int_t prems, double x, double y=0, double z=0) const
Definition Octree.cpp:231
_SIZE_ int_t
Definition Octree.h:56
DoubleTab_T< _SIZE_ > DoubleTab_t
Definition Octree.h:60
virtual int_t taille() const
Renvoie la taille de l'octree.
Definition Octree.cpp:615
void ranger_elem_3D(ArrOfInt &oks, int_t elem, int_t idx, int nb_som_elem, const DoubleTab_t &coords_som, const IntTab_t &les_elems, SmallArrOfTID_t &compteur, Vect_IntTab_t &sous_tab, double xmil, double ymil, double zmil)
Definition Octree.cpp:463
IntTab_T< _SIZE_ > IntTab_t
Definition Octree.h:59
virtual int_t rang_elem_loc(const OctreeLoc &, double x, double y=0, double z=0) const
Definition Octree.cpp:195
TRUST_Vector< IntTab_t > Vect_IntTab_t
Definition Octree.h:61
static int nombre_d_octrees()
Definition Octree.h:76
int niveau() const
Definition Octree.cpp:602
ArrOfInt_T< _SIZE_ > ArrOfInt_t
Definition Octree.h:57
Octree_32_64()
Definition Octree.h:66
void ranger_elem_1D(ArrOfInt &oks, int_t elem, int_t idx, int nb_som_elem, const DoubleTab_t &coords_som, const IntTab_t &les_elems, SmallArrOfTID_t &compteur, Vect_IntTab_t &sous_tab, double xmil)
Definition Octree.cpp:258
SmallArrOfTID_T< _SIZE_ > SmallArrOfTID_t
Definition Octree.h:58
void construire(int, const ArrOfInt_t &, const OctreeLoc &, Octree_t *p=0)
Definition Octree.cpp:503
double get_epsilon() const
Definition Octree.h:89
virtual Sortie & printOn(Sortie &is) const
Definition Octree.cpp:27
Octree_32_64 * pere
Definition Octree.h:92
void ranger_elem_2D(ArrOfInt &oks, int_t elem, int_t idx, int nb_som_elem, const DoubleTab_t &coords_som, const IntTab_t &les_elems, SmallArrOfTID_t &compteur, Vect_IntTab_t &sous_tab, double xmil, double ymil)
Definition Octree.cpp:438
virtual Entree & readOn(Entree &is)
Definition Octree.h:79
virtual const Domaine_t & domaine() const
Definition Octree.h:72
Domaine_32_64< _SIZE_ > Domaine_t
Definition Octree.h:63
Octree_32_64< _SIZE_ > Octree_t
Definition Octree.h:64
static Sortie & Journal(int message_level=0)
Renvoie un objet statique de type Sortie qui sert de journal d'evenements.
Definition Process.cpp:588
static void exit(int exit_code=-1)
Routine de sortie de TRUST dans une region Kokkos.
Definition Process.cpp:455
Classe de base des flux de sortie.
Definition Sortie.h:52
_SIZE_ size_array() const
void resize_array(_SIZE_ new_size, RESIZE_OPTIONS opt=RESIZE_OPTIONS::COPY_INIT)
int dimension_int(int d) const
Definition TRUSTTab.tpp:152
_SIZE_ dimension_tot(int) const override
Definition TRUSTTab.tpp:160
_SIZE_ dimension(int d) const
Definition TRUSTTab.tpp:133
classe TRUST_Vector
void resize(int i)
Struct OctreeLoc.
Definition Octree.h:40
double xmax
Definition Octree.h:41
double zmin
Definition Octree.h:41
void construire(const OctreeLoc &loc, int i)
Definition Octree.cpp:61
double ymax
Definition Octree.h:41
double zmax
Definition Octree.h:41
double xmin
Definition Octree.h:41
double ymin
Definition Octree.h:41
int direction(double x, double y, double z) const
Definition Octree.cpp:157