TrioCFD 1.9.9_beta
TrioCFD documentation
Loading...
Searching...
No Matches
Pave.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#include <Pave.h>
17#include <math.h>
18#include <Domaine.h>
19
20Implemente_instanciable_32_64(Pave_32_64,"Pave",Domaine_32_64<_T_>);
21
22// XD mailler_base objet_lecture mailler_base INHERITS_BRACE Basic class to mesh.
23
24// XD pave mailler_base pave NO_BRACE Class to create a pave (block) with boundaries.
25// XD attr name chaine name REQ Name of the pave (block).
26// XD attr bloc bloc_pave bloc REQ Definition of the pave (block).
27// XD attr list_bord list_bord list_bord REQ Domain boundaries definition.
28
29// XD pave_64 pave pave_64 NO_BRACE Same as Pave for big (64b) domain
30
31// XD bloc_pave objet_lecture nul INHERITS_BRACE Class to create a pave.
32// XD attr Origine listf Origine OPT Keyword to define the pave (block) origin, that is to say one of the 8 block points
33// XD_CONT (or 4 in a 2D coordinate system).
34// XD attr longueurs listf longueurs OPT Keyword to define the block dimensions, that is to say knowing the origin,
35// XD_CONT length along the axes.
36// XD attr nombre_de_noeuds listentierf nombre_de_noeuds OPT Keyword to define the discretization (nodenumber) in each
37// XD_CONT direction.
38// XD attr facteurs listf facteurs OPT Keyword to define stretching factors for mesh discretization in each direction.
39// XD_CONT This is a real number which must be positive (by default 1.0). A stretching factor other than 1 allows
40// XD_CONT refinement on one edge in one direction.
41// XD attr symx rien symx OPT Keyword to define a block mesh that is symmetrical with respect to the YZ plane
42// XD_CONT (respectively Y-axis in 2D) passing through the block centre.
43// XD attr symy rien symy OPT Keyword to define a block mesh that is symmetrical with respect to the XZ plane
44// XD_CONT (respectively X-axis in 2D) passing through the block centre.
45// XD attr symz rien symz OPT Keyword defining a block mesh that is symmetrical with respect to the XY plane passing
46// XD_CONT through the block centre.
47// XD attr xtanh floattant xtanh OPT Keyword to generate mesh with tanh (hyperbolic tangent) variation in the
48// XD_CONT X-direction.
49// XD attr xtanh_dilatation entier(into=[-1,0,1]) xtanh_dilatation OPT Keyword to generate mesh with tanh (hyperbolic
50// XD_CONT tangent) variation in the X-direction. xtanh_dilatation: The value may be -1,0,1 (0 by default): 0: coarse
51// XD_CONT mesh at the middle of the channel and smaller near the walls -1: coarse mesh at the left side of the channel
52// XD_CONT and smaller at the right side 1: coarse mesh at the right side of the channel and smaller near the left side
53// XD_CONT of the channel.
54// XD attr xtanh_taille_premiere_maille floattant xtanh_taille_premiere_maille OPT Size of the first cell of the mesh
55// XD_CONT with tanh (hyperbolic tangent) variation in the X-direction.
56// XD attr ytanh floattant ytanh OPT Keyword to generate mesh with tanh (hyperbolic tangent) variation in the
57// XD_CONT Y-direction.
58// XD attr ytanh_dilatation entier(into=[-1,0,1]) ytanh_dilatation OPT Keyword to generate mesh with tanh (hyperbolic
59// XD_CONT tangent) variation in the Y-direction. ytanh_dilatation: The value may be -1,0,1 (0 by default): 0: coarse
60// XD_CONT mesh at the middle of the channel and smaller near the walls -1: coarse mesh at the bottom of the channel and
61// XD_CONT smaller near the top 1: coarse mesh at the top of the channel and smaller near the bottom.
62// XD attr ytanh_taille_premiere_maille floattant ytanh_taille_premiere_maille OPT Size of the first cell of the mesh
63// XD_CONT with tanh (hyperbolic tangent) variation in the Y-direction.
64// XD attr ztanh floattant ztanh OPT Keyword to generate mesh with tanh (hyperbolic tangent) variation in the
65// XD_CONT Z-direction.
66// XD attr ztanh_dilatation entier(into=[-1,0,1]) ztanh_dilatation OPT Keyword to generate mesh with tanh (hyperbolic
67// XD_CONT tangent) variation in the Z-direction. tanh_dilatation: The value may be -1,0,1 (0 by default): 0: coarse
68// XD_CONT mesh at the middle of the channel and smaller near the walls -1: coarse mesh at the back of the channel and
69// XD_CONT smaller near the front 1: coarse mesh at the front of the channel and smaller near the back.
70// XD attr ztanh_taille_premiere_maille floattant ztanh_taille_premiere_maille OPT Size of the first cell of the mesh
71// XD_CONT with tanh (hyperbolic tangent) variation in the Z-direction.
72
73// XD epsilon mailler_base epsilon NO_BRACE Two points will be confused if the distance between them is less than eps.
74// XD_CONT By default, eps is set to 1e-12. The keyword Epsilon allows an alternative value to be assigned to eps.
75// XD attr eps floattant eps REQ New value of precision.
76
77// XD domain mailler_base domain NO_BRACE Class to reuse a domain.
78// XD attr domain_name ref_domaine domain_name REQ Name of domain.
79
80// XD list_bloc_mailler listobj list_bloc_mailler BRACE mailler_base COMMA List of block mesh.
81
82template <typename _SIZE_>
83Sortie& Pave_32_64<_SIZE_>::printOn(Sortie& s ) const
84{
85 return Domaine_t::printOn(s) ;
86}
87
88/*! @brief Reads the specifications of a pave from an input stream.
89 *
90 * The format for reading a pave in the data set is as follows:
91 * Pave nom_pave
92 * {
93 * Origine OX OY (OZ)
94 * Longueurs LX LY (LZ)
95 * Nombre_de_noeuds NX NY (NZ)
96 * Facteurs Fx Fy (Fz)
97 * (Symx)
98 * (Symy)
99 * (Symz)
100 * }
101 * {
102 * (Bord) nom X = X0 Y0 <= Y <= Y1 Z0 <= Z <= Z1
103 * ...
104 * (Raccord) local homogene nom X = X0 Y0 <= Y <= Y1 Z0 <= Z <= Z1
105 * ...
106 * (Internes) nom X = X0 Y0 <= Y <= Y1 Z0 <= Z <= Z1
107 * ...
108 * (Joint) nom X = X0 Y0 <= Y <= Y1 Z0 <= Z <= Z1 PE_voisin
109 * ...
110 * (Groupe_Faces) nom X = X0 Y0 <= Y <= Y1 Z0 <= Z <= Z1
111 * ...
112 * }
113 *
114 * @param (Entree& is) an input stream
115 * @return (Entree&) the modified input stream
116 * @throws space dimension required for meshing
117 * @throws opening brace expected
118 * @throws Symy only makes sense for dimension >= 2
119 * @throws Symz only makes sense in dimension 3
120 * @throws progression factors must be positive
121 * @throws there must be at least two cells in x
122 * @throws opening brace expected before reading boundaries
123 * @throws unrecognised keyword
124 */
125template <typename _SIZE_>
126Entree& Pave_32_64<_SIZE_>::readOn(Entree& is)
127{
128 if(this->dimension<0)
129 {
130 Cerr << "You must give the dimension before to mesh" << finl;
131 Cerr << " The syntax is \"Dimension dim\" " << finl ;
132 Process::exit();
133 }
134 facteurs_.resize(this->dimension);
135 facteurs_ = 1.;
136 symetrique_.resize(this->dimension);
137 symetrique_=0;
138 origine_.resize(this->dimension);
139 longueurs_.resize(this->dimension);
140 nb_noeuds_.resize(this->dimension);
141 Motcle motlu;
142 int rang;
143
144 typer_();
145
146 is >> this->nom_;
147 Cerr << "Reading of the block " << this->nom_ << finl;
148 is >> motlu;
149 if (motlu!="{")
150 {
151 Cerr << "We expected a { after " << this->nom_ << finl;
152 Process::exit();
153 }
154 {
155 Motcles les_mots(18);
156 {
157 les_mots[0]="Origine";
158 les_mots[1]="Longueurs";
159 les_mots[2]="Nombre_de_noeuds";
160 les_mots[3]="Facteurs";
161 les_mots[4]="Symx";
162 les_mots[5]="Symy";
163 les_mots[6]="Symz";
164 les_mots[7]="xtanh";
165 les_mots[8]="xtanh_dilatation";
166 les_mots[9]="xtanh_taille_premiere_maille";
167 les_mots[10]="ytanh";
168 les_mots[11]="ytanh_dilatation";
169 les_mots[12]="ytanh_taille_premiere_maille";
170 les_mots[13]="ztanh";
171 les_mots[14]="ztanh_dilatation";
172 les_mots[15]="ztanh_taille_premiere_maille";
173 les_mots[16]="reprise_VEF";
174 les_mots[17]="}";
175 }
176 while(motlu != "}")
177 {
178 is >> motlu;
179 rang = les_mots.search(motlu);
180 switch(rang)
181 {
182 case 0:
183 for(int i=0; i< this->dimension; i++)
184 is >> origine_(i);
185 break;
186 case 1:
187 lire_longueurs(is);
188 break;
189 case 2:
190 lire_noeuds(is);
191 break;
192 case 3:
193 for(int i=0; i< this->dimension; i++)
194 is >> facteurs_(i);
195 break;
196 case 4: // Symx
197 symetrique_(0)=1;
198 break;
199 case 5: // Symy
200 if(this->dimension<2)
201 {
202 Cerr << "Symy has meaning only in dimension >= 2" << finl;
203 Process::exit();
204 }
205 symetrique_(1)=1;
206 break;
207 case 6: // Symz
208 if(this->dimension<3)
209 {
210 Cerr << "Symz has meaning only in dimension >= 3" << finl;
211 Process::exit();
212 }
213 symetrique_(2)=1;
214 break;
215 case 7:
216 {
217 is >> xa_tanh;
218 break;
219 }
220 case 8:
221 {
222 is >> xtanh_dilatation;
223 break;
224 }
225 case 9:
226 {
227 double x1, x_tmp;
228 is >> x1;
229 xa_tanh=.1;
230 double fac_sym;
231 assert(xtanh_dilatation!=-123);
232 if (xtanh_dilatation != 0 )
233 fac_sym=1.;
234 else
235 fac_sym=2.;
236 assert(!est_egal(xa_tanh,-123.));
237 for(int decimale=1; decimale<7; decimale++)
238 {
239 x_tmp=longueurs_(0)/fac_sym*(1.+tanh((-1.+fac_sym*1./((Mx-1)*1.))*atanh(xa_tanh))/xa_tanh);
240 while ( x_tmp > x1 )
241 {
242 xa_tanh+=pow(10.,-decimale);
243 xa_tanh = std::min(xa_tanh,1-Objet_U::precision_geom);
244 if ( std::fabs(xa_tanh) < Objet_U::precision_geom)
245 {
246 Cerr << "Error: The coefficient xa_tanh has a value of : " << xa_tanh << finl ;
247 Cerr << "So the mesh can't be generated with tanh (hyperbolic tangent) variation!" << finl ;
248 Cerr << "You must decrease either the xtanh_taille_premiere_maille size of the first cell of the mesh " << finl;
249 Cerr << "or the nombre_de_noeuds Nx node number in the X direction." << finl;
250 Process::exit();
251 }
252 x_tmp=longueurs_(0)/fac_sym*(1.+tanh((-1.+fac_sym*1./((Mx-1)*1.))*atanh(xa_tanh))/xa_tanh);
253 }
254 xa_tanh-=pow(10.,-decimale);
255 xa_tanh = std::min(xa_tanh,1-Objet_U::precision_geom);
256 if ( std::fabs(xa_tanh) < Objet_U::precision_geom)
257 {
258 Cerr << "Error: The coefficient xa_tanh has a value of : " << xa_tanh << finl ;
259 Cerr << "So the mesh can't be generated with tanh (hyperbolic tangent) variation!" << finl ;
260 Cerr << "You must decrease either the xtanh_taille_premiere_maille size of the first cell of the mesh " << finl;
261 Cerr << "or the nombre_de_noeuds Nx node number in the X direction." << finl;
262 Process::exit();
263 }
264 }
265 Cerr << "The coefficient xa_tanh has a value of : " << xa_tanh << finl ;
266 break;
267 }
268 case 10:
269 {
270 is >> a_tanh;
271 break;
272 }
273 case 11:
274 {
275 is >> tanh_dilatation;
276 break;
277 }
278 case 12:
279 {
280 double y1, y_tmp;
281 is >> y1;
282 a_tanh=.1;
283 double fac_sym;
284 assert(tanh_dilatation!=-123);
285 if (tanh_dilatation != 0 )
286 fac_sym=1.;
287 else
288 fac_sym=2.;
289 assert(!est_egal(a_tanh,-123.));
290 for(int decimale=1; decimale<7; decimale++)
291 {
292 y_tmp=longueurs_(1)/fac_sym*(1.+tanh((-1.+fac_sym*1./((My-1)*1.))*atanh(a_tanh))/a_tanh);
293 while ( y_tmp > y1 )
294 {
295 a_tanh+=pow(10.,-decimale);
296 a_tanh = std::min(a_tanh,1-Objet_U::precision_geom);
297 if ( std::fabs(a_tanh) < Objet_U::precision_geom)
298 {
299 Cerr << "Error: The coefficient ya_tanh has a value of : " << a_tanh << finl ;
300 Cerr << "So the mesh can't be generated with tanh (hyperbolic tangent) variation!" << finl ;
301 Cerr << "You must decrease either the ytanh_taille_premiere_maille size of the first cell of the mesh " << finl;
302 Cerr << "or the nombre_de_noeuds Ny node number in the Y direction." << finl;
303 Process::exit();
304 }
305 y_tmp=longueurs_(1)/fac_sym*(1.+tanh((-1.+fac_sym*1./((My-1)*1.))*atanh(a_tanh))/a_tanh);
306 }
307 a_tanh-=pow(10.,-decimale);
308 a_tanh = std::min(a_tanh,1-Objet_U::precision_geom);
309 if ( std::fabs(a_tanh) < Objet_U::precision_geom)
310 {
311 Cerr << "Error: The coefficient ya_tanh has a value of : " << a_tanh << finl ;
312 Cerr << "So the mesh can't be generated with tanh (hyperbolic tangent) variation!" << finl ;
313 Cerr << "You must decrease either the ytanh_taille_premiere_maille size of the first cell of the mesh " << finl;
314 Cerr << "or the nombre_de_noeuds Ny node number in the Y direction." << finl;
315 Process::exit();
316 }
317 }
318 Cerr << "The coefficient ya_tanh has a value of : " << a_tanh << finl ;
319 break;
320 }
321 case 13:
322 {
323 is >> za_tanh;
324 break;
325 }
326 case 14:
327 {
328 is >> ztanh_dilatation;
329 break;
330 }
331 case 15:
332 {
333 double z1, z_tmp;
334 is >> z1;
335 za_tanh=.1;
336 double fac_sym;
337 assert(ztanh_dilatation!=-123);
338 if (ztanh_dilatation != 0 )
339 fac_sym=1.;
340 else
341 fac_sym=2.;
342 assert(!est_egal(za_tanh,-123.));
343 for(int decimale=1; decimale<7; decimale++)
344 {
345 z_tmp=longueurs_(2)/fac_sym*(1.+tanh((-1.+fac_sym*1./((Mz-1)*1.))*atanh(za_tanh))/za_tanh);
346 while ( z_tmp > z1 )
347 {
348 za_tanh+=pow(10.,-decimale);
349 za_tanh = std::min(za_tanh,1-Objet_U::precision_geom);
350 if ( std::fabs(za_tanh) < Objet_U::precision_geom)
351 {
352 Cerr << "Error: The coefficient za_tanh has a value of : " << za_tanh << finl ;
353 Cerr << "So the mesh can't be generated with tanh (hyperbolic tangent) variation!" << finl ;
354 Cerr << "You must decrease either the ztanh_taille_premiere_maille size of the first cell of the mesh " << finl;
355 Cerr << "or the nombre_de_noeuds Nz node number in the Z direction." << finl;
356 Process::exit();
357 }
358 z_tmp=longueurs_(2)/fac_sym*(1.+tanh((-1.+fac_sym*1./((Mz-1)*1.))*atanh(za_tanh))/za_tanh);
359 }
360 za_tanh-=pow(10.,-decimale);
361 za_tanh = std::min(za_tanh,1-Objet_U::precision_geom);
362 if ( std::fabs(za_tanh) < Objet_U::precision_geom)
363 {
364 Cerr << "Error: The coefficient za_tanh has a value of : " << za_tanh << finl ;
365 Cerr << "So the mesh can't be generated with tanh (hyperbolic tangent) variation!" << finl ;
366 Cerr << "You must decrease either the ztanh_taille_premiere_maille size of the first cell of the mesh " << finl;
367 Cerr << "or the nombre_de_noeuds Nz node number in the Z direction." << finl;
368 Process::exit();
369 }
370 }
371 Cerr << "The coefficient za_tanh has a value of : " << za_tanh << finl ;
372 break;
373 }
374 case 16:
375 {
376 rep_VEF=true;
377 break;
378 }
379 case 17: // Fin
380 break;
381 default:
382 if (motlu == "tanh" || motlu == "tanh_dilatation" || motlu == "tanh_taille_premiere_maille")
383 Cerr << "Error: '" << motlu << "' keyword is obsolete since V1.7.9. We renamed it to 'Y"<< motlu <<"'" << finl;
384 else
385 {
386 Cerr << motlu << " is not understood " << finl;
387 Cerr << les_mots;
388 }
389 Process::exit();
390 }
391 }
392 pas_.resize(this->dimension);
393 if(min_array(facteurs_)<=0.)
394 {
395 Cerr << "The progression factors must be positives" << finl;
396 Process::exit();
397 }
398 switch(this->dimension)
399 {
400 case 2:
401 maille2D();
402 break;
403 case 3:
404 maille3D();
405 break;
406 }
407 }
408 if(Mx<2)
409 {
410 Cerr << "There must be at least two cells on x" << finl;
411 Process::exit();
412 }
413 if (this->axi)
414 {
415 // Coordinates are modulo 2 pi
416 double deux_pi=M_PI*2.0 ;
417 {
418 int_t nb_som=this->sommets_.dimension(0);
419 for(int_t i=0; i<nb_som; i++)
420 {
421 this->sommets_(i, 1) -= static_cast<double>((int_t) this->sommets_(i, 1)); // remove integral part
422 this->sommets_(i, 1) *= deux_pi;
423 }
424 }
425
426 longueurs_(1)*=deux_pi;
427 origine_(1)*=deux_pi;
428 }
429 {
430 is >> motlu;
431 if (motlu!="{")
432 {
433 Cerr << "We expected a { before reading the boundaries" << finl;
434 Process::exit();
435 }
436 Motcles les_mots(7);
437 {
438 les_mots[0]="Bord";
439 les_mots[1]="Raccord";
440 les_mots[2]="Internes";
441 les_mots[3]="Trou";
442 les_mots[4]="Joint";
443 les_mots[5]="Groupe_Faces";
444 les_mots[6]="}";
445 }
446 while(motlu != "}")
447 {
448 is >> motlu;
449 rang = les_mots.search(motlu);
450 switch(rang)
451 {
452 case 0: // Bord
453 {
454 Bord_32_64<_SIZE_>& newbord=this->mes_faces_bord_.add(Bord_32_64<_SIZE_>());
455 lire_front(is , newbord );
456 }
457 break;
458 case 1: // Raccord
459 {
460 OWN_PTR(Raccord_base_32_64<_SIZE_>)& newraccord=this->mes_faces_raccord_.add(OWN_PTR(Raccord_base_32_64<_SIZE_>)());
461 Nom type="Raccord_";
462 Nom local, homogene;
463 is >> local >> homogene;
464 type+=local;
465 type+="_";
466 type+=homogene;
467 if (type == "local")
468 {
469 Cerr << "Do not use sequential local connection!" << finl;
470 Process::exit();
471 }
472 newraccord.typer(type);
473 lire_front(is , newraccord.valeur() );
474 }
475 break;
476 case 2: // Plaques
477 {
478 Bord_Interne_32_64<_SIZE_>& faces=this->mes_bords_int_.add(Bord_Interne_32_64<_SIZE_>());
479 lire_front(is , faces );
480 }
481 break;
482 case 5: // Groupe de Faces
483 {
484 Groupe_Faces_32_64<_SIZE_>& faces=this->mes_groupes_faces_.add(Groupe_Faces_32_64<_SIZE_>());
485 lire_front(is , faces );
486 }
487 break;
488 case 3: // Trous
489 break;
490 case 4: // Joint
491 {
492 Joint_32_64<_SIZE_>& newjoint=this->mes_faces_joint_.add(Joint_32_64<_SIZE_>());
493 lire_front(is , newjoint );
494 }
495 break;
496 case 6: // Fin
497 break;
498 default:
499 Cerr << motlu << " is not understood " << finl;
500 Cerr << les_mots;
501 Process::exit();
502 }
503 }
504 }
505
506 return is;
507}
508
509///*! @brief Performs a 1D meshing of the pave using the parameter values read by ReadOn.
510// *
511// */
512//template <typename _SIZE_>
513//void Pave_32_64<_SIZE_>::maille1D()
514//{
515// Cerr << "Step of 1D mesh in progress... " << finl;
516// double epsilon_geom=Objet_U::precision_geom;
517// int i;
518// if(longueurs_(0)<0)
519// {
520// origine_(0)+=longueurs_(0);
521// longueurs_(0)=-longueurs_(0);
522// facteurs_(0)=1./facteurs_(0);
523// }
524// if( (!symetrique_(0)) || (std::fabs(facteurs_(0)-1.)<epsilon_geom) )
525// {
526// if(std::fabs(facteurs_(0)-1.)>epsilon_geom)
527// pas_(0)=longueurs_(0)*(facteurs_(0)-1)/(pow(facteurs_(0),Nx)-1);
528// else
529// pas_(0)=longueurs_(0)/Nx;
530// double dx=pas_(0);
531// double x=origine_(0);
532// for (i=0; i<Nx; i++)
533// {
534// coord_noeud(i)=x;
535// x+=dx;
536// dx*=facteurs_(0);
537// }
538// coord_noeud(Nx)=origine_(0)+longueurs_(0);
539// }
540// else
541// {
542// int Ix=Nx/2;
543// if( (Ix*2)==Nx )
544// pas_(0)=0.5*longueurs_(0)*(facteurs_(0)-1)/(pow(facteurs_(0),Ix)-1);
545// else
546// {
547// if(Ix)
548// {
549// Ix+=1;
550// pas_(0)=0.5*longueurs_(0)/( (pow(facteurs_(0),Ix-1)-1)/(facteurs_(0)-1)
551// +0.5*pow(facteurs_(0),Ix) );
552// }
553// else
554// {
555// Ix+=1;
556// pas_(0)=0;
557// }
558// }
559// double dx=pas_(0);
560// double x=0;
561// for (i=0; i<Ix; i++)
562// {
563// coord_noeud(i)=origine_(0)+x;
564// coord_noeud(Nx-i)=origine_(0)+longueurs_(0)-x;
565// x+=dx;
566// dx*=facteurs_(0);
567// }
568// if( (Ix*2)==Nx )
569// coord_noeud(Ix)=0.5*(coord_noeud(Ix+1)+coord_noeud(Ix-1));
570// }
571// Cerr << finl;
572// Cerr << "Step of mesh ended " << finl;
573//}
574
575
576/*! @brief Performs a 2D meshing of the box with the values of the parameters read by ReadOn.
577 *
578 */
579template <typename _SIZE_>
580void Pave_32_64<_SIZE_>::maille2D()
581{
582 Cerr << "Step of 2D mesh in progress... " << finl;
583 double epsilon_geom=Objet_U::precision_geom;
584 int i, j;
585 int pourcent=0;
586 int tmp;
587 if(longueurs_(0)<0)
588 {
589 origine_(0)+=longueurs_(0);
590 longueurs_(0)=-longueurs_(0);
591 facteurs_(0)=1./facteurs_(0);
592 }
593 if(longueurs_(1)<0)
594 {
595 origine_(1)+=longueurs_(1);
596 longueurs_(1)=-longueurs_(1);
597 facteurs_(1)=1./facteurs_(1);
598 }
599 if( (!symetrique_(0)) || (std::fabs(facteurs_(0)-1.)<epsilon_geom) )
600 {
601 assert(!est_egal(xa_tanh,-123.));
602 if ((xa_tanh>epsilon_geom)&&(xa_tanh<1.) )
603 {
604 /// Tanh mesh in x
605 Cerr << "In if( (xa_tanh>epsilon_geom)&&(xa_tanh<1.) )" << finl;
606 pourcent =0;
607 double fac_sym;
608 assert(xtanh_dilatation!=-123);
609 if (xtanh_dilatation != 0)
610 fac_sym=1.;
611 else
612 fac_sym=2.;
613
614 if (xtanh_dilatation != -1 )
615 {
616 for (i=0; i<Mx; i++)
617 {
618 Cerr << "i=" << i << finl;
619 tmp=(i*100)/Mx;
620 if(tmp > pourcent)
621 {
622 pourcent=tmp;
623 Cerr << "\rIt was created " << pourcent << "% of coordx" << flush;
624 }
625 assert(!est_egal(xa_tanh,-123.));
626 for (j=0; j<My; j++)
627 coord_noeud(i,j,0)=origine_(0)+longueurs_(0)/fac_sym*(1.+tanh((-1.+fac_sym*i/((Mx-1)*1.))*atanh(xa_tanh))/xa_tanh);
628 Cerr << "It was created " << pourcent << "% of coordx" << flush;
629 }
630 }
631 else
632 {
633 for (i=0; i<Mx; i++)
634 {
635 Cerr << "i=" << i << finl;
636 tmp=(i*100)/Mx;
637 if(tmp > pourcent)
638 {
639 pourcent=tmp;
640 Cerr << "\rIt was created " << pourcent << "% of coordx" << flush;
641 }
642 assert(!est_egal(xa_tanh,-123.));
643 for (j=0; j<My; j++)
644 coord_noeud(i,j,0)=origine_(0)+longueurs_(0)-longueurs_(0)/fac_sym*(1.+tanh((-1.+fac_sym*(Mx-1-i)/((Mx-1)*1.))*atanh(xa_tanh))/xa_tanh);
645 Cerr << "It was created " << pourcent << "% of coordx" << flush;
646 }
647 }
648 ///////////////////////////////
649 for (j=0; j<My; j++)
650 coord_noeud(0,j,0)=origine_(0);
651 }
652 else
653 {
654 double x;
655 for (i=0; i<Nx; i++)
656 {
657 tmp=(i*100)/Nx;
658 if(tmp > pourcent)
659 {
660 pourcent=tmp;
661 Cerr << "It was created " << pourcent << "% of coordx" << flush;
662 }
663 if (std::fabs(facteurs_(0)-1.)>epsilon_geom)
664 x=origine_(0)+longueurs_(0)*(pow(facteurs_(0),i)-1)/(pow(facteurs_(0),Nx)-1);
665 else
666 x=origine_(0)+longueurs_(0)*i/Nx;
667 for (j=0; j<My; j++)
668 coord_noeud(i,j,0)=x;
669 }
670 for (j=0; j<My; j++)
671 coord_noeud(Nx,j,0)=origine_(0)+longueurs_(0);
672 }
673 }
674 else
675 {
676 int Ix=Nx/2;
677 if( (Ix*2)==Nx )
678 pas_(0)=0.5*longueurs_(0)*(facteurs_(0)-1)/(pow(facteurs_(0),Ix)-1);
679 else
680 {
681 if(Ix)
682 {
683 Ix+=1;
684 pas_(0)=0.5*longueurs_(0)/( (pow(facteurs_(0),Ix-1)-1)/(facteurs_(0)-1)
685 +0.5*pow(facteurs_(0),Ix) );
686 }
687 else
688 {
689 Ix+=1;
690 pas_(0)=longueurs_(0);
691 }
692 }
693 double dx=pas_(0);
694 double x=0;
695 for (i=0; i<Ix; i++)
696 {
697 tmp=(i*100)/Ix;
698 if(tmp > pourcent)
699 {
700 pourcent=tmp;
701 Cerr << "It was created " << pourcent << "% of coordx" << flush;
702 }
703 for (j=0; j<My; j++)
704 {
705 coord_noeud(i,j,0)=origine_(0)+x;
706 coord_noeud(Nx-i,j,0)=origine_(0)+longueurs_(0)-x;
707 }
708 x+=dx;
709 dx*=facteurs_(0);
710 }
711 if( (Ix*2)==Nx )
712 for (j=0; j<My; j++)
713 coord_noeud(Ix,j,0)=0.5*(coord_noeud(Ix+1,j,0)
714 +coord_noeud(Ix-1,j,0));
715 }
716 pourcent=0;
717 //
718 if( (!symetrique_(1)) || (std::fabs(facteurs_(1)-1.)<epsilon_geom) )
719 {
720 assert(!est_egal(a_tanh,-123.));
721 if ((a_tanh>epsilon_geom)&&(a_tanh<1.) )
722 {
723 /// Tanh mesh in y
724 Cerr << "In if( (a_tanh>epsilon_geom)&&(a_tanh<1.) )" << finl;
725 pourcent =0;
726 double fac_sym;
727 assert(tanh_dilatation!=-123);
728 if (tanh_dilatation != 0)
729 fac_sym=1.;
730 else
731 fac_sym=2.;
732
733 if (tanh_dilatation != -1 )
734 {
735 for (j=0; j<My; j++)
736 {
737 Cerr << "j=" << j << finl;
738 tmp=(j*100)/My;
739 if(tmp > pourcent)
740 {
741 pourcent=tmp;
742 Cerr << "\rIt was created " << pourcent << "% of coordy" << flush;
743 }
744 assert(!est_egal(a_tanh,-123.));
745 for (i=0; i<Mx; i++)
746 coord_noeud(i,j,1)=origine_(1)+longueurs_(1)/fac_sym*(1.+tanh((-1.+fac_sym*j/((My-1)*1.))*atanh(a_tanh))/a_tanh);
747 Cerr << "coord_noeud(0,j,1)=" << coord_noeud(0,j,1) << finl;
748 Cerr << "It was created " << pourcent << "% of coordy" << flush;
749 }
750 }
751 else
752 {
753 for (j=0; j<My; j++)
754 {
755 Cerr << "j=" << j << finl;
756 tmp=(j*100)/My;
757 if(tmp > pourcent)
758 {
759 pourcent=tmp;
760 Cerr << "\rIt was created " << pourcent << "% of coordy" << flush;
761 }
762 assert(!est_egal(a_tanh,-123.));
763 for (i=0; i<Mx; i++)
764 coord_noeud(i,j,1)=origine_(1)+longueurs_(1)-longueurs_(1)/fac_sym*(1.+tanh((-1.+fac_sym*(My-1-j)/((My-1)*1.))*atanh(a_tanh))/a_tanh);
765 Cerr << "coord_noeud(0,j,1)=" << coord_noeud(0,j,1) << finl;
766 Cerr << "It was created " << pourcent << "% of coordy" << flush;
767 }
768 }
769 ///////////////////////////////
770 for (i=0; i<Mx; i++)
771 coord_noeud(i,0,1)=origine_(1);
772 }
773 else
774 {
775 double y;
776 for (j=0; j<Ny; j++)
777 {
778 tmp=(j*100)/Ny;
779 if(tmp > pourcent)
780 {
781 pourcent=tmp;
782 Cerr << "It was created " << pourcent << "% of coordy" << flush;
783 }
784 if (std::fabs(facteurs_(1)-1.)>epsilon_geom)
785 y=origine_(1)+longueurs_(1)*(pow(facteurs_(1),j)-1)/(pow(facteurs_(1),Ny)-1);
786 else
787 y=origine_(1)+longueurs_(1)*j/Ny;
788 for (i=0; i<Mx; i++)
789 coord_noeud(i,j,1)=y;
790 }
791 for (i=0; i<Mx; i++)
792 coord_noeud(i,j,1)=origine_(1)+longueurs_(1);
793 }
794 }
795 else
796 {
797 int Iy=Ny/2;
798 if( (Iy*2)==Ny )
799 pas_(1)=0.5*longueurs_(1)*(facteurs_(1)-1)/(pow(facteurs_(1),Iy)-1);
800 else
801 {
802 if(Iy)
803 {
804 Iy+=1;
805 pas_(1)=0.5*longueurs_(1)/( (pow(facteurs_(1),Iy-1)-1)/(facteurs_(1)-1)
806 +0.5*pow(facteurs_(1),Iy) );
807 }
808 else
809 {
810 Iy+=1;
811 pas_(1)=longueurs_(1);
812 }
813 }
814 double dy=pas_(1);
815 double y=0;
816 for (j=0; j<Iy; j++)
817 {
818 tmp=(j*100)/Iy;
819 if(tmp > pourcent)
820 {
821 pourcent=tmp;
822 Cerr << "It was created " << pourcent << "% of coordy" << flush;
823 }
824 for (i=0; i<Mx; i++)
825 {
826 coord_noeud(i,j,1)=origine_(1)+y;
827 coord_noeud(i,Ny-j,1)=origine_(1)+longueurs_(1)-y;
828 }
829 y+=dy;
830 dy*=facteurs_(1);
831 }
832 if( (Iy*2)==Ny )
833 for (i=0; i<Mx; i++)
834 coord_noeud(i,Iy,1)=0.5*(coord_noeud(i,Iy-1,1)
835 +coord_noeud(i,Iy+1,1));
836 }
837 Cerr << finl;
838 Cerr << "Step of mesh ended " << finl;
839}
840
841/*! @brief Performs a 3D meshing of the box with the values of the parameters read by ReadOn.
842 */
843template <typename _SIZE_>
844void Pave_32_64<_SIZE_>::maille3D()
845{
846 Cerr << "Step of 3D mesh in progress... " << finl;
847 double epsilon_geom=Objet_U::precision_geom;
848 int i, j, k;
849 int pourcent=0;
850 int tmp;
851 if(longueurs_(0)<0)
852 {
853 origine_(0)+=longueurs_(0);
854 longueurs_(0)=-longueurs_(0);
855 facteurs_(0)=1./facteurs_(0);
856 }
857 if(longueurs_(1)<0)
858 {
859 origine_(1)+=longueurs_(1);
860 longueurs_(1)=-longueurs_(1);
861 facteurs_(1)=1./facteurs_(1);
862 }
863 if(longueurs_(2)<0)
864 {
865 origine_(2)+=longueurs_(2);
866 longueurs_(2)=-longueurs_(2);
867 facteurs_(2)=1./facteurs_(2);
868 }
869 if( (!symetrique_(0)) || (std::fabs(facteurs_(0)-1.)<epsilon_geom) )
870 {
871 assert(!est_egal(xa_tanh,-123.));
872 if ((xa_tanh>epsilon_geom)&&(xa_tanh<1.) )
873 {
874 /// Tanh mesh in x
875 Cerr << "In if( (xa_tanh>epsilon)&&(xa_tanh<1.) )" << finl;
876 pourcent =0;
877 double fac_sym;
878 assert(xtanh_dilatation!=-123);
879 if (xtanh_dilatation != 0)
880 fac_sym=1.;
881 else
882 fac_sym=2.;
883
884 if (xtanh_dilatation != -1 )
885 {
886 for (i=0; i<Mx; i++)
887 {
888 Cerr << "i=" << i << finl;
889 tmp=(i*100)/Mx;
890 if(tmp > pourcent)
891 {
892 pourcent=tmp;
893 Cerr << "\rIt was created " << pourcent << "% of coordx" << flush;
894 }
895 assert(!est_egal(xa_tanh,-123.));
896 for (j=0; j<My; j++)
897 for (k=0; k<Mz; k++)
898 coord_noeud(i,j,k,0)=origine_(0)+longueurs_(0)/fac_sym*(1.+tanh((-1.+fac_sym*i/((Mx-1)*1.))*atanh(xa_tanh))/xa_tanh);
899 Cerr << "coord_noeud(i,0,0,0)=" << coord_noeud(i,0,0,0) << finl;
900 }
901 }
902 else
903 {
904 for (i=0; i<Mx; i++)
905 {
906 Cerr << "i=" << i << finl;
907 tmp=(i*100)/Mx;
908 if(tmp > pourcent)
909 {
910 pourcent=tmp;
911 Cerr << "It was created " << pourcent << "% of coordx" << flush;
912 }
913 assert(!est_egal(xa_tanh,-123.));
914 for (j=0; j<My; j++)
915 for (k=0; k<Mz; k++)
916 coord_noeud(i,j,k,0)=origine_(0)+longueurs_(0)-longueurs_(0)/fac_sym*(1.+tanh((-1.+fac_sym*(Mx-1-i)/((Mx-1)*1.))*atanh(xa_tanh))/xa_tanh);
917
918 Cerr << "coord_noeud(i,0,0,0)=" << coord_noeud(i,0,0,0) << finl;
919 }
920 }
921 for (j=0; j<My; j++)
922 for (k=0; k<Mz; k++)
923 coord_noeud(0,j,k,0)=origine_(0);
924 }
925 ///////////////////////////////
926 else
927 {
928
929 double x;
930 for (i=0; i<Nx; i++)
931 {
932 tmp=(i*100)/Nx;
933 if(tmp > pourcent)
934 {
935 pourcent=tmp;
936 Cerr << "It was created " << pourcent << "% of coordx" << flush;
937 }
938 if(std::fabs(facteurs_(0)-1.)>epsilon_geom)
939 x=origine_(0)+longueurs_(0)*(pow(facteurs_(0),i)-1)/(pow(facteurs_(0),Nx)-1);
940 else
941 x=origine_(0)+longueurs_(0)*i/Nx;
942 for (j=0; j<My; j++)
943 for (k=0; k<Mz; k++)
944 coord_noeud(i,j,k,0)=x;
945 }
946 for (j=0; j<My; j++)
947 for (k=0; k<Mz; k++)
948 coord_noeud(Nx,j,k,0)=origine_(0)+longueurs_(0);
949 }
950 }
951 else
952 {
953 int Ix=Nx/2;
954 if( (Ix*2)==Nx )
955 pas_(0)=0.5*longueurs_(0)*(facteurs_(0)-1)/(pow(facteurs_(0),Ix)-1);
956 else
957 {
958 if(Ix)
959 {
960 Ix+=1;
961 pas_(0)=0.5*longueurs_(0)/( (pow(facteurs_(0),Ix-1)-1)/(facteurs_(0)-1)
962 +0.5*pow(facteurs_(0),Ix) );
963 }
964 else
965 {
966 Ix+=1;
967 pas_(0)=longueurs_(0);
968 }
969 }
970 double dx=pas_(0);
971 double x=0;
972 for (i=0; i<Ix; i++)
973 {
974 tmp=(i*100)/Ix;
975 if(tmp > pourcent)
976 {
977 pourcent=tmp;
978 Cerr << "It was created " << pourcent << "% of coordx" << flush;
979 }
980 for (j=0; j<My; j++)
981 for (k=0; k<Mz; k++)
982 {
983 coord_noeud(i,j,k,0)=origine_(0)+x;
984 coord_noeud(Nx-i,j,k,0)=origine_(0)+longueurs_(0)-x;
985 }
986 x+=dx;
987 dx*=facteurs_(0);
988 }
989 if( (Ix*2)==Nx )
990 for (j=0; j<My; j++)
991 for (k=0; k<Mz; k++)
992 coord_noeud(Ix,j,k,0)=0.5*(coord_noeud(Ix+1,j,k,0)+
993 coord_noeud(Ix-1,j,k,0));
994 }
995 //
996 pourcent=0;
997 Cerr << finl;
998 if( (!symetrique_(1)) || (std::fabs(facteurs_(1)-1.)<epsilon_geom) )
999 {
1000 assert(!est_egal(a_tanh,-123.));
1001 if ((a_tanh>epsilon_geom)&&(a_tanh<1.) )
1002 {
1003 /// Tanh mesh in y
1004 Cerr << "In if( (a_tanh>epsilon)&&(a_tanh<1.) )" << finl;
1005 pourcent =0;
1006 double fac_sym;
1007 assert(tanh_dilatation!=-123);
1008 if (tanh_dilatation != 0)
1009 fac_sym=1.;
1010 else
1011 fac_sym=2.;
1012
1013 if (tanh_dilatation != -1 )
1014 {
1015 for (j=0; j<My; j++)
1016 {
1017 Cerr << "j=" << j << finl;
1018 tmp=(j*100)/My;
1019 if(tmp > pourcent)
1020 {
1021 pourcent=tmp;
1022 Cerr << "\rIt was created " << pourcent << "% of coordy" << flush;
1023 }
1024 assert(!est_egal(a_tanh,-123.));
1025 for (i=0; i<Mx; i++)
1026 for (k=0; k<Mz; k++)
1027 coord_noeud(i,j,k,1)=origine_(1)+longueurs_(1)/fac_sym*(1.+tanh((-1.+fac_sym*j/((My-1)*1.))*atanh(a_tanh))/a_tanh);
1028 Cerr << "coord_noeud(0,j,0,1)=" << coord_noeud(0,j,0,1) << finl;
1029 }
1030 }
1031 else
1032 {
1033 for (j=0; j<My; j++)
1034 {
1035 Cerr << "j=" << j << finl;
1036 tmp=(j*100)/My;
1037 if(tmp > pourcent)
1038 {
1039 pourcent=tmp;
1040 Cerr << "It was created " << pourcent << "% of coordy" << flush;
1041 }
1042 assert(!est_egal(a_tanh,-123.));
1043 for (i=0; i<Mx; i++)
1044 for (k=0; k<Mz; k++)
1045 coord_noeud(i,j,k,1)=origine_(1)+longueurs_(1)-longueurs_(1)/fac_sym*(1.+tanh((-1.+fac_sym*(My-1-j)/((My-1)*1.))*atanh(a_tanh))/a_tanh);
1046
1047 Cerr << "coord_noeud(0,j,0,1)=" << coord_noeud(0,j,0,1) << finl;
1048 }
1049 }
1050 for (i=0; i<Mx; i++)
1051 for (k=0; k<Mz; k++)
1052 coord_noeud(i,0,k,1)=origine_(1);
1053 }
1054 ///////////////////////////////
1055 else
1056 {
1057 double y;
1058 for (j=0; j<Ny; j++)
1059 {
1060 tmp=(j*100)/Ny;
1061 if(tmp > pourcent)
1062 {
1063 pourcent=tmp;
1064 Cerr << "It was created " << pourcent << "% of coordy" << flush;
1065 }
1066 if(std::fabs(facteurs_(1)-1.)>epsilon_geom)
1067 y=origine_(1)+longueurs_(1)*(pow(facteurs_(1),j)-1)/(pow(facteurs_(1),Ny)-1);
1068 else
1069 y=origine_(1)+longueurs_(1)*j/Ny;
1070
1071 for (i=0; i<Mx; i++)
1072 for (k=0; k<Mz; k++)
1073 coord_noeud(i,j,k,1)=y;
1074 }
1075 for (i=0; i<Mx; i++)
1076 for (k=0; k<Mz; k++)
1077 coord_noeud(i,j,k,1)=origine_(1)+longueurs_(1);
1078 }
1079 }
1080 else
1081 {
1082 int Iy=Ny/2;
1083 if( (Iy*2)==Ny )
1084 pas_(1)=0.5*longueurs_(1)*(facteurs_(1)-1)/(pow(facteurs_(1),Iy)-1);
1085 else
1086 {
1087 if(Iy)
1088 {
1089 Iy+=1;
1090 pas_(1)=0.5*longueurs_(1)/( (pow(facteurs_(1),Iy-1)-1)/(facteurs_(1)-1)
1091 +0.5*pow(facteurs_(1),Iy) );
1092 }
1093 else
1094 {
1095 Iy+=1;
1096 pas_(1)=longueurs_(1);
1097 }
1098 }
1099 double dy=pas_(1);
1100 double y=0;
1101 for (j=0; j<Iy; j++)
1102 {
1103 tmp=(j*100)/Iy;
1104 if(tmp > pourcent)
1105 {
1106 pourcent=tmp;
1107 Cerr << "It was created " << pourcent << "% of coordy" << flush;
1108 }
1109 for (i=0; i<Mx; i++)
1110 for (k=0; k<Mz; k++)
1111 {
1112 coord_noeud(i,j,k,1)=origine_(1)+y;
1113 coord_noeud(i,Ny-j,k,1)=origine_(1)+longueurs_(1)-y;
1114 }
1115 y+=dy;
1116 dy*=facteurs_(1);
1117 }
1118 if( (Iy*2)==Ny )
1119 for (i=0; i<Mx; i++)
1120 for (k=0; k<Mz; k++)
1121 coord_noeud(i,Iy,k,1)=0.5*(coord_noeud(i,Iy-1,k,1)
1122 +coord_noeud(i,Iy+1,k,1));
1123 }
1124 //
1125 pourcent=0;
1126 Cerr << finl;
1127 if( (!symetrique_(2)) ||( std::fabs(facteurs_(2)-1.)<epsilon_geom) )
1128 {
1129 assert(!est_egal(za_tanh,-123.));
1130 if ((za_tanh>epsilon_geom)&&(za_tanh<1.) )
1131 {
1132 /// Tanh mesh in z
1133 Cerr << "In if( (za_tanh>epsilon)&&(za_tanh<1.) )" << finl;
1134 pourcent =0;
1135 double fac_sym;
1136 assert(ztanh_dilatation!=-123);
1137 if (ztanh_dilatation != 0)
1138 fac_sym=1.;
1139 else
1140 fac_sym=2.;
1141
1142 if (ztanh_dilatation != -1 )
1143 {
1144 for (k=0; k<Mz; k++)
1145 {
1146 Cerr << "k=" << k << finl;
1147 tmp=(k*100)/Mz;
1148 if(tmp > pourcent)
1149 {
1150 pourcent=tmp;
1151 Cerr << "\rIt was created " << pourcent << "% of coordz" << flush;
1152 }
1153 assert(!est_egal(za_tanh,-123.));
1154 for (i=0; i<Mx; i++)
1155 for (j=0; j<My; j++)
1156 coord_noeud(i,j,k,2)=origine_(2)+longueurs_(2)/fac_sym*(1.+tanh((-1.+fac_sym*k/((Mz-1)*1.))*atanh(za_tanh))/za_tanh);
1157 Cerr << "coord_noeud(0,0,k,2)=" << coord_noeud(0,0,k,2) << finl;
1158 }
1159 }
1160 else
1161 {
1162 for (k=0; k<Mz; k++)
1163 {
1164 Cerr << "k=" << k << finl;
1165 tmp=(k*100)/Mz;
1166 if(tmp > pourcent)
1167 {
1168 pourcent=tmp;
1169 Cerr << "It was created " << pourcent << "% of coordz" << flush;
1170 }
1171 assert(!est_egal(za_tanh,-123.));
1172 for (i=0; i<Mx; i++)
1173 for (j=0; j<My; j++)
1174 coord_noeud(i,j,k,2)=origine_(2)+longueurs_(2)-longueurs_(2)/fac_sym*(1.+tanh((-1.+fac_sym*(Mz-1-k)/((Mz-1)*1.))*atanh(za_tanh))/za_tanh);
1175
1176 Cerr << "coord_noeud(0,0,k,2)=" << coord_noeud(0,0,k,2) << finl;
1177 }
1178 }
1179 for (i=0; i<Mx; i++)
1180 for (j=0; j<My; j++)
1181 coord_noeud(i,j,0,2)=origine_(2);
1182 }
1183 ///////////////////////////////
1184 else
1185 {
1186 double z;
1187 for (k=0; k<Nz; k++)
1188 {
1189 tmp=(k*100)/Nz;
1190 if(tmp > pourcent)
1191 {
1192 pourcent=tmp;
1193 Cerr << "It was created " << pourcent << "% of coordz" << flush;
1194 }
1195 if(std::fabs(facteurs_(2)-1.)>epsilon_geom)
1196 z=origine_(2)+longueurs_(2)*(pow(facteurs_(2),k)-1)/(pow(facteurs_(2),Nz)-1);
1197 else
1198 z=origine_(2)+longueurs_(2)*k/Nz;
1199 for (i=0; i<Mx; i++)
1200 for (j=0; j<My; j++)
1201 coord_noeud(i,j,k,2)=z;
1202 }
1203 for (i=0; i<Mx; i++)
1204 for (j=0; j<My; j++)
1205 coord_noeud(i,j,Nz,2)=origine_(2)+longueurs_(2);
1206 }
1207 }
1208 else
1209 {
1210 int Iz=Nz/2;
1211 if( (Iz*2)==Nz )
1212 pas_(2)=0.5*longueurs_(2)*(facteurs_(2)-1)/(pow(facteurs_(2),Iz)-1);
1213 else
1214 {
1215 if(Iz)
1216 {
1217 Iz+=1;
1218 pas_(2)=0.5*longueurs_(2)/( (pow(facteurs_(2),Iz-1)-1)/(facteurs_(2)-1)
1219 +0.5*pow(facteurs_(2),Iz) );
1220 }
1221 else
1222 {
1223 Iz+=1;
1224 pas_(2)=longueurs_(2);
1225 }
1226 }
1227 double dz=pas_(2);
1228 double z=0;
1229 for (k=0; k<Iz; k++)
1230 {
1231 tmp=(k*100)/Iz;
1232 if(tmp > pourcent)
1233 {
1234 pourcent=tmp;
1235 Cerr << "It was created " << pourcent << "% of coordz" << flush;
1236 }
1237 for (i=0; i<Mx; i++)
1238 for (j=0; j<My; j++)
1239 {
1240 coord_noeud(i,j,k,2)=origine_(2)+z;
1241 coord_noeud(i,j,Nz-k,2)=origine_(2)+longueurs_(2)-z;
1242 }
1243 z+=dz;
1244 dz*=facteurs_(2);
1245 }
1246 if( (Iz*2)==Nz )
1247 for (i=0; i<Mx; i++)
1248 for (j=0; j<My; j++)
1249 coord_noeud(i,j,Iz,2)=0.5*(coord_noeud(i,j,Iz+1,2)
1250 +coord_noeud(i,j,Iz-1,2));
1251 }
1252
1253
1254 // PaQ: to allow switching from a VDF mesh to a VEF mesh using tetraedriser homogene
1255
1256 if(rep_VEF)
1257 {
1258 for (j=1; j<My; j+=2)
1259 for (i=0; i<Mx; i++)
1260 for (k=0; k<Mz; k++)
1261 coord_noeud(i,j,k,1)=coord_noeud(i,j-1,k,1)+0.5*(coord_noeud(i,j+1,k,1)-coord_noeud(i,j-1,k,1));
1262 }
1263
1264
1265 Cerr << finl;
1266 Cerr << "Step of mesh ended " << finl;
1267}
1268
1269
1270/*! @brief Types the pave according to the space dimension and the requested coordinate system (cylindrical
1271 *
1272 * or Cartesian).
1273 * In dimension 1 : the pave is meshed with segments
1274 * In dimension 2 (axi): the pave is meshed with rectangles (axi)
1275 * In dimension 3 (axi): the pave is meshed with hexahedra (axi)
1276 */
1277template <typename _SIZE_>
1278void Pave_32_64<_SIZE_>::typer_()
1279{
1280 Nom typ;
1281 switch(this->dimension)
1282 {
1283 case 1:
1284 typ = "Segment";
1285 break;
1286 case 2:
1287 if (this->axi)
1288 typ = "Rectangle_Axi";
1289 else if (this->bidim_axi)
1290 typ = "Rectangle_2D_Axi";
1291 else
1292 typ = "Rectangle";
1293 break;
1294 case 3:
1295 if (this->axi)
1296 typ = "Hexaedre_Axi";
1297 else
1298 typ = "Hexaedre";
1299 break;
1300 default :
1301 throw;
1302 }
1303
1304 Domaine_t::typer(typ);
1305}
1306
1307/*! @brief Reads the lengths LX LY [LZ] from the data set input stream.
1308 *
1309 * Method called by readOn(Entree&)
1310 *
1311 * @param (Entree& is) an input stream
1312 * @throws The Length in axi is given as number of turns, must be between 0 and 1
1313 */
1314template <typename _SIZE_>
1315void Pave_32_64<_SIZE_>::lire_longueurs(Entree& is)
1316{
1317 int i;
1318 for(i=0; i< this->dimension; i++)
1319 is >> longueurs_(i);
1320 if(this->axi)
1321 {
1322 if(longueurs_(1) >1.+this->epsilon_)
1323 {
1324 Cerr << "The length \"Longueurs\" is in number of turns in axisymmetric" << finl;
1325 Cerr << "It must be between 0 and 1" << finl;
1326 Process::exit();
1327 }
1328 else if(longueurs_(1)==1.)
1329 tour_complet=true;
1330 }
1331}
1332
1333/*! @brief Reads the number of nodes from the data set input stream and builds the nodes.
1334 *
1335 * Method called by readOn(Entree&)
1336 *
1337 * @param (Entree& is) an input stream
1338 * @throws in axi coordinates, lengths must be read first
1339 * @throws unsupported space dimension
1340 */
1341template <typename _SIZE_>
1342void Pave_32_64<_SIZE_>::lire_noeuds(Entree& is)
1343{
1344 int i, j, k;
1345 if( (this->axi) && (longueurs_(0)==0.))
1346 {
1347 Cerr<< "You must first read the length in axisymmetric!" << finl;
1348 Process::exit();
1349 }
1350 for(i=0; i< this->dimension; i++)
1351 is >> nb_noeuds_(i);
1352
1353 // If the user is using 'Pave' instead of 'Pave_64', check no overflow:
1354 if (std::is_same<_SIZE_, int>::value)
1355 {
1356 trustIdType nb_nod = 1;
1357 for(int d=0; d < Objet_U::dimension; d++)
1358 nb_nod *= nb_noeuds_[d];
1359 if(8*nb_nod > std::numeric_limits<int>::max()) // 8* ... because later storage of elem connectivity will require 8 columns
1360 {
1361 Cerr << "ERROR: Pave: the total number of nodes you specified is bigger than the 32b limit!" << finl;
1362 Cerr << " If this is intended, you should switch to the 64 mode by using the 'Pave_64' (and Domaine_64 etc.) keyword." << finl;
1363 Process::exit();
1364 }
1365 }
1366
1367 if (min_array(nb_noeuds_)<2)
1368 {
1369 Cerr << "\nError: The number of nodes in directions Nx and Ny (and Nz) for 'Pave " << this->nom_ << "' must be greater than 1." << finl;
1370 Cerr << "If you want to define a unique cell in a given direction, set the number of nodes to 2 in that direction." << finl;
1371 Process::exit();
1372 }
1373 if(this->dimension==1)
1374 {
1375 Mx=nb_noeuds_(0);
1376 Nx=nb_noeuds_(0)-1;
1377 this->mes_elems_.resize(Nx,2);
1378 this->sommets_.resize(Mx);
1379 for (i=0; i<Nx; i++)
1380 {
1381 maille_sommet(i,0)=numero_sommet(i);
1382 maille_sommet(i,1)=numero_sommet(i+1);
1383 }
1384 }
1385 else if(this->dimension==2)
1386 {
1387 Mx=nb_noeuds_(0);
1388 My=nb_noeuds_(1);
1389 Nx=nb_noeuds_(0)-1;
1390 Ny=nb_noeuds_(1)-1;
1391 if(tour_complet) Ny++;
1392 this->mes_elems_.resize(Nx*Ny,4);
1393 this->sommets_.resize(Mx*My,2);
1394 for (i=0; i<Nx; i++)
1395 for (j=0; j<Ny; j++)
1396 {
1397 maille_sommet(i,j,0)=numero_sommet(i, j);
1398 maille_sommet(i,j,1)=numero_sommet(i+1, j);
1399 maille_sommet(i,j,2)=numero_sommet(i, j+1);
1400 maille_sommet(i,j,3)=numero_sommet(i+1, j+1);
1401 }
1402 }
1403 else if(this->dimension==3)
1404 {
1405 Mx=nb_noeuds_(0);
1406 My=nb_noeuds_(1);
1407 Mz=nb_noeuds_(2);
1408 Nx=nb_noeuds_(0)-1;
1409 Ny=nb_noeuds_(1)-1;
1410 Nz=nb_noeuds_(2)-1;
1411 if(tour_complet) Ny++;
1412 this->mes_elems_.resize(Nx*Ny*Nz,8);
1413 this->sommets_.resize(Mx*My*Mz,3);
1414 for (i=0; i<Nx; i++)
1415 for ( j=0; j<Ny; j++)
1416 for ( k=0; k<Nz; k++)
1417 {
1418 maille_sommet(i,j,k,0)=numero_sommet(i, j, k);
1419 maille_sommet(i,j,k,1)=numero_sommet(i+1, j, k);
1420 maille_sommet(i,j,k,2)=numero_sommet(i, j+1, k);
1421 maille_sommet(i,j,k,3)=numero_sommet(i+1, j+1, k);
1422 maille_sommet(i,j,k,4)=numero_sommet(i, j, k+1);
1423 maille_sommet(i,j,k,5)=numero_sommet(i+1, j, k+1);
1424 maille_sommet(i,j,k,6)=numero_sommet(i, j+1, k+1);
1425 maille_sommet(i,j,k,7)=numero_sommet(i+1, j+1, k+1);
1426 }
1427 }
1428 else
1429 {
1430 Cerr << "dimension = " << this->dimension << "not provided " << finl;
1431 Process::exit();
1432 }
1433}
1434
1435/*! @brief Reads the specifications of a boundary from the data set input stream and builds it.
1436 *
1437 * Format:
1438 * nom_front X = X0 Y0 <= Y <= Y1 Z0 <= Z <= Z1
1439 *
1440 * @param (Entree& is) an input stream
1441 * @param (Frontiere& front) the boundary being read
1442 * @throws keyword "X" expected
1443 * @throws keyword "=" expected
1444 * @throws invalid X endpoint
1445 * @throws keyword "X" or "Y" expected
1446 * @throws keyword "<=" expected
1447 * @throws keyword "Y" expected
1448 * @throws invalid Y endpoint
1449 * @throws there is no boundary in theta, you have meshed
1450 * a complete annulus
1451 * @throws keyword "X" or "Y" or "Z" expected
1452 * @throws keyword "Z" expected
1453 * @throws invalid Z endpoint
1454 */
1455template <typename _SIZE_>
1456void Pave_32_64<_SIZE_>::lire_front(Entree& is, Frontiere_t& front)
1457{
1458 int i, j, k;
1459 Nom nom_front;
1460 is >> nom_front;
1461 Cerr << "Reading of the boundary " << nom_front << finl;
1462 front.nommer(nom_front);
1463 front.typer_faces(this->elem_->type_face());
1464 bool internes=(sub_type(Bord_Interne_32_64<_SIZE_>, front) || sub_type(Joint_32_64<_SIZE_>, front));
1465 bool groupe_faces=sub_type(Groupe_Faces_32_64<_SIZE_>, front);
1466 if(this->dimension==1)
1467 {
1468 Nom X, egal;
1469 double coupe;
1470 is >> X >> egal >> coupe ;
1471 if( X!="X")
1472 {
1473 Cerr << "We expected a = X instead of" << X <<
1474 " before " << egal << " " << coupe << finl;
1475 Process::exit();
1476 }
1477 if( egal!="=")
1478 {
1479 Cerr << "We expected a = after " << X << finl;
1480 Process::exit();
1481 }
1482 front.dimensionner(0);
1483 IntTab_t som(1,1);
1484 if ( (coupe != origine_(0)) && (coupe != (origine_(0)+longueurs_(0))) )
1485 {
1486 Cerr << coupe << " is not an extremity" << finl;
1487 Cerr << "The extremities are : "
1488 << origine_(0) << " " << origine_(0)+longueurs_(0) << finl;
1489 Process::exit();
1490 }
1491 if(std::fabs(coupe - origine_(0))<this->epsilon_)
1492 som(0,0)=0;
1493 else
1494 som(0,0)=Mx;
1495 front.ajouter_faces(som);
1496 front.associer_domaine(*this);
1497 }
1498 else if(this->dimension==2)
1499 {
1500 Nom X, Y, egal, infegal1, infegal2;
1501 double coupe;
1502 double xmin, xmax;
1503 is >> X >> egal >> coupe ;
1504 if( (X!="X") && (X!="Y") )
1505 {
1506 Cerr << "We expected a = X or Y instead of" << X <<
1507 " before " << egal << " " << coupe << finl;
1508 Process::exit();
1509 }
1510 if( egal!="=")
1511 {
1512 Cerr << "We expected a = after " << X << finl;
1513 Process::exit();
1514 }
1515 is >> xmin >> infegal1 >> Y >> infegal2 >> xmax;
1516 if( infegal1!="<=")
1517 {
1518 Cerr << "We expected a <= after " << xmin << " and not " << infegal1 << finl;
1519 Process::exit();
1520 }
1521 if (X=="X")
1522 {
1523 if( Y!="Y" )
1524 {
1525 Cerr << "We expected a Y after " << infegal1 << finl;
1526 Process::exit();
1527 }
1528 if (this->axi)
1529 {
1530 double deux_pi=M_PI*2.0 ;
1531 xmin*=deux_pi;
1532 xmax*=deux_pi;
1533 }
1534 }
1535 else
1536 {
1537 if(tour_complet)
1538 {
1539 Cerr << "There is no boundary in teta! " << finl;
1540 Cerr << "You have meshed a complete crown! " << finl;
1541 Process::exit();
1542 }
1543 if( Y!="X" )
1544 {
1545 Cerr << "We expected a X after " << infegal1 << finl;
1546 Process::exit();
1547 }
1548 }
1549 if( infegal2!="<=")
1550 {
1551 Cerr << "We expected a <= after " << Y << " and not " << infegal2 <<finl;
1552 Process::exit();
1553 }
1554 front.dimensionner(0);
1555 IntTab_t som;
1556 if (X=="X")
1557 {
1558 if ( (std::fabs(origine_(0) - coupe)>this->epsilon_) &&
1559 (std::fabs(coupe - origine_(0)-longueurs_(0))>this->epsilon_) && (!internes) && (!groupe_faces))
1560 {
1561 Cerr << "X = " << coupe << " is not a boundary" << finl;
1562 Process::exit();
1563 }
1564 if ( ((std::fabs(origine_(0) - coupe)<this->epsilon_) ||
1565 (std::fabs(coupe - origine_(0)-longueurs_(0))<this->epsilon_)) && (internes))
1566 {
1567 Cerr << "X = " << coupe << " is a boundary" << finl;
1568 Process::exit();
1569 }
1570 int jmin, jmax;
1571 if (std::fabs(origine_(0) - coupe)<this->epsilon_)
1572 i=0;
1573 else
1574 i=Nx;
1575 if(internes || groupe_faces)
1576 for(i=0; coord_noeud(i,0,0)+this->epsilon_<coupe; i++) {};
1577 jmin=0;
1578 for(; coord_noeud(0,jmin,1)+this->epsilon_<xmin; jmin++) {};
1579 if(std::fabs(xmax-(origine_(1)+longueurs_(1)))<this->epsilon_)
1580 {
1581 jmax=My-1+tour_complet;
1582 }
1583 else
1584 {
1585 jmax = jmin;
1586 for(; coord_noeud(0,jmax,1)+this->epsilon_<xmax; jmax++) {};
1587 // if(jmax-jmin==0)
1588 // {
1589 // Cerr << "boundary named " << nom_front << " is empty!" << finl;
1590 // Process::exit();
1591 // }
1592 }
1593 som.resize(jmax-jmin,2);
1594 for(j=jmin; j<jmax; j++)
1595 {
1596 som(j-jmin,0)=numero_sommet(i, j);
1597 som(j-jmin,1)=numero_sommet(i, j+1);
1598 }
1599 }
1600 else
1601 {
1602 if (this->axi)
1603 {
1604 double deux_pi=M_PI*2.0 ;
1605 coupe*=deux_pi;
1606 }
1607 if ( (std::fabs(origine_(1) - coupe)>this->epsilon_) &&
1608 (std::fabs(coupe - origine_(1)-longueurs_(1))>this->epsilon_) && (!internes) && (!groupe_faces))
1609 {
1610 Cerr << "Y = " << coupe << " is not a boundary" << finl;
1611 Process::exit();
1612 }
1613 if ( ((std::fabs(origine_(1) - coupe)<this->epsilon_) ||
1614 (std::fabs(coupe - origine_(1)-longueurs_(1))<this->epsilon_)) && (internes))
1615 {
1616 Cerr << "Y = " << coupe <<" is a boundary" << finl;
1617 Process::exit();
1618 }
1619 int imin, imax;
1620 if (std::fabs(coupe-origine_(1))<this->epsilon_)
1621 j=0;
1622 else
1623 j=Ny;
1624 if(internes || groupe_faces)
1625 for(j=0; coord_noeud(0,j,1)+this->epsilon_<coupe; j++) {};
1626 imin=0;
1627 for(; coord_noeud(imin,0,0)+this->epsilon_<xmin; imin++) {};
1628 if(std::fabs(xmax-(origine_(0)+longueurs_(0)))<this->epsilon_)
1629 imax=Mx-1;
1630 else
1631 {
1632 imax = imin;
1633 for(; coord_noeud(imax,0,0)+this->epsilon_<xmax; imax++) {};
1634 // if(imax-imin==0)
1635 // {
1636 // Cerr << "boundary named " << nom_front << " is empty!" << finl;
1637 // Process::exit();
1638 // }
1639 }
1640 som.resize(imax-imin,2);
1641 for(i=imin; i<imax; i++)
1642 {
1643 som(i-imin,0)=numero_sommet(i, j);
1644 som(i-imin,1)=numero_sommet(i+1, j);
1645 }
1646 }
1647 front.ajouter_faces(som);
1648 front.associer_domaine(*this);
1649 }
1650 else if(this->dimension==3)
1651 {
1652 Nom X, Y, Z, egal, infegal1, infegal2, infegal3, infegal4;
1653 double coupe;
1654 double xmin, xmax, ymin, ymax;
1655 is >> X >> egal >> coupe ;
1656 if( (X!="X") && (X!="Y") && (X!="Z"))
1657 {
1658 Cerr << "We expected a = X or Y or Z instead of" << X <<
1659 " before " << egal << " " << coupe << finl;
1660 Process::exit();
1661 }
1662 if( egal!="=")
1663 {
1664 Cerr << "We expected a = after " << X << finl;
1665 Process::exit();
1666 }
1667 is >> xmin >> infegal1 >> Y >> infegal2 >> xmax;
1668 if( infegal1!="<=")
1669 {
1670 Cerr << "We expected a <= after " << xmin << " and not " << infegal1 <<finl;
1671 Process::exit();
1672 }
1673 if( infegal2!="<=")
1674 {
1675 Cerr << "We expected a <= after " << Y << " and not " << infegal2 << finl;
1676 Process::exit();
1677 }
1678 is >> ymin >> infegal3 >> Z >> infegal4 >> ymax;
1679 if( infegal3!="<=")
1680 {
1681 Cerr << "We expected a <= after " << ymin << " and not " << infegal3 << finl;
1682 Process::exit();
1683 }
1684 if( infegal4!="<=")
1685 {
1686 Cerr << "We expected a <= after " << Z << " and not " << infegal4 << finl;
1687 Process::exit();
1688 }
1689 if (X=="X")
1690 {
1691 if( Y!="Y" )
1692 {
1693 Cerr << "We expected a Y and not " << Y << finl;
1694 Process::exit();
1695 }
1696 if (this->axi)
1697 {
1698 double deux_pi=M_PI*2.0 ;
1699 xmin*=deux_pi;
1700 xmax*=deux_pi;
1701 }
1702 if( Z!="Z")
1703 {
1704 Cerr << "We expected a Z and not " << Z << finl;
1705 Process::exit();
1706 }
1707 }
1708 else if(X=="Y")
1709 {
1710 if(tour_complet)
1711 {
1712 Cerr << "There is no boundary in teta! " << finl;
1713 Cerr << "You have meshed a complete crown! " << finl;
1714 Process::exit();
1715 }
1716 if( Y!="X" )
1717 {
1718 Cerr << "We expected a X and not " << Y << finl;
1719 Process::exit();
1720 }
1721 if( Z!="Z")
1722 {
1723 Cerr << "We expected a Z and not " << Z << finl;
1724 Process::exit();
1725 }
1726 }
1727 else
1728 {
1729 if( Y!="X" )
1730 {
1731 Cerr << "We expected a X and not " << Y << finl;
1732 Process::exit();
1733 }
1734 if( Z!="Y" )
1735 {
1736 Cerr << "We expected a Z and not " << Z << finl;
1737 Process::exit();
1738 }
1739 if (this->axi)
1740 {
1741 double deux_pi=M_PI*2.0 ;
1742 ymin*=deux_pi;
1743 ymax*=deux_pi;
1744 }
1745 }
1746 front.dimensionner(0);
1747 IntTab_t som;
1748 if (X=="X")
1749 {
1750 if ( (std::fabs(origine_(0)-coupe)>this->epsilon_) &&
1751 (std::fabs(origine_(0)+longueurs_(0)-coupe)>this->epsilon_) && (!internes) && (!groupe_faces))
1752 {
1753 Cerr << "X = " << coupe << " is not a boundary " << finl;
1754 Process::exit();
1755 }
1756 if ( ((std::fabs(origine_(0)-coupe)<this->epsilon_) ||
1757 (std::fabs(origine_(0)+longueurs_(0)-coupe)<this->epsilon_)) && (internes) )
1758 {
1759 Cerr << "X = " << coupe << " is a boundary " << finl;
1760 Process::exit();
1761 }
1762 int jmin, jmax, kmin, kmax;
1763 if (std::fabs(origine_(0)-coupe)<this->epsilon_)
1764 i=0;
1765 else
1766 i=Nx;
1767 if(internes || groupe_faces)
1768 for(i=0; coord_noeud(i,0,0,0)+this->epsilon_<coupe; i++) {};
1769 jmin=kmin=0;
1770 for(; coord_noeud(0,jmin,0,1)+this->epsilon_<xmin; jmin++) {};
1771 for(; coord_noeud(0,0,kmin,2)+this->epsilon_<ymin; kmin++) {};
1772 jmax = jmin;
1773 kmax=kmin;
1774 if(std::fabs(xmax-(origine_(1)+longueurs_(1)))<this->epsilon_)
1775 {
1776 jmax=My-1+tour_complet;
1777 }
1778 else
1779 for(; coord_noeud(0,(jmax),0,1)+this->epsilon_<xmax; jmax++) {};
1780 for(; coord_noeud(0,0,(kmax),2)+this->epsilon_<ymax; kmax++) {};
1781 // if((jmax-jmin)*(kmax-kmin)==0)
1782 // {
1783 // Cerr << "boundary named " << nom_front << " is empty!" << finl;
1784 // Process::exit();
1785 // }
1786 som.resize((jmax-jmin)*(kmax-kmin),4);
1787 for(j=jmin; j<jmax; j++)
1788 for(k=kmin; k<kmax; k++)
1789 {
1790 som((k-kmin)*(jmax-jmin)+j-jmin,0)=numero_sommet(i, j, k);
1791 som((k-kmin)*(jmax-jmin)+j-jmin,1)=numero_sommet(i, j+1, k);
1792 som((k-kmin)*(jmax-jmin)+j-jmin,2)=numero_sommet(i, j, k+1);
1793 som((k-kmin)*(jmax-jmin)+j-jmin,3)=numero_sommet(i, j+1, k+1);
1794 }
1795 }
1796 else if(X=="Y")
1797 {
1798 if (this->axi)
1799 {
1800 double deux_pi=M_PI*2.0 ;
1801 coupe*=deux_pi;
1802 }
1803 if ( (std::fabs(origine_(1) - coupe)>this->epsilon_) &&
1804 (std::fabs(coupe - origine_(1)-longueurs_(1))>this->epsilon_) && (!internes) && (!groupe_faces) )
1805 {
1806 Cerr << "Y = " << coupe << " is not a boundary " << finl;
1807 Process::exit();
1808 }
1809 if ( ((std::fabs(origine_(1) - coupe)<this->epsilon_) ||
1810 (std::fabs(coupe - origine_(1)-longueurs_(1))<this->epsilon_)) && (internes))
1811 {
1812 Cerr << "Y = " << coupe << " is a boundary " << finl;
1813 Process::exit();
1814 }
1815 int imin, imax, kmin, kmax;
1816 if (std::fabs(coupe-origine_(1))<this->epsilon_)
1817 j=0;
1818 else
1819 j=Ny;
1820 if(internes || groupe_faces)
1821 for(j=0; coord_noeud(0,j,0,1)+this->epsilon_<coupe; j++) {};
1822 imin=kmin=0;
1823 for(; coord_noeud(imin,0,0,0)+this->epsilon_<xmin; imin++) {};
1824 for(; coord_noeud(0,0,kmin,2)+this->epsilon_<ymin; kmin++) {};
1825 imax = imin;
1826 kmax=kmin;
1827 for(; coord_noeud(imax,0,0,0)+this->epsilon_<xmax; imax++) {};
1828 for(; coord_noeud(0,0,(kmax),2)+this->epsilon_<ymax; kmax++) {};
1829 // if((imax-imin)*(kmax-kmin)==0)
1830 // {
1831 // Cerr << "boundary named " << nom_front << " is empty!" << finl;
1832 // Process::exit();
1833 // }
1834 som.resize((imax-imin)*(kmax-kmin),4);
1835 for(i=imin; i<imax; i++)
1836 for(k=kmin; k<kmax; k++)
1837 {
1838 som((k-kmin)*(imax-imin)+i-imin,0)=numero_sommet(i, j, k);
1839 som((k-kmin)*(imax-imin)+i-imin,1)=numero_sommet(i+1, j, k);
1840 som((k-kmin)*(imax-imin)+i-imin,2)=numero_sommet(i, j, k+1);
1841 som((k-kmin)*(imax-imin)+i-imin,3)=numero_sommet(i+1, j, k+1);
1842 }
1843 }
1844 else
1845 {
1846 if ( (std::fabs(origine_(2)-coupe)>this->epsilon_) &&
1847 (std::fabs(coupe-origine_(2)-longueurs_(2))>this->epsilon_) && (!internes) && (!groupe_faces))
1848 {
1849 Cerr << "Z = " << coupe << " is not a boundary " << finl;
1850 Process::exit();
1851 }
1852 if ( ((std::fabs(origine_(2)-coupe)<this->epsilon_) ||
1853 (std::fabs(coupe-origine_(2)-longueurs_(2))<this->epsilon_)) && (internes))
1854 {
1855 Cerr << "Z = " << coupe << " is a boundary " << finl;
1856 Process::exit();
1857 }
1858 int imin, imax, jmin, jmax;
1859 if (std::fabs(coupe-origine_(2))<this->epsilon_)
1860 k=0;
1861 else
1862 k=Nz;
1863 if(internes || groupe_faces)
1864 for(k=0; coord_noeud(0,0,k,3)+this->epsilon_<coupe; k++) {};
1865 imin=jmin=0;
1866 for(; coord_noeud(imin,0,0,0)+this->epsilon_<xmin; imin++) {};
1867 for(; coord_noeud(0,jmin,0,1)+this->epsilon_<ymin; jmin++) {};
1868 imax = imin;
1869 jmax=jmin;
1870 for(; coord_noeud(imax,0,0,0)+this->epsilon_<xmax; imax++) {};
1871 if(std::fabs(ymax-(origine_(1)+longueurs_(1)))<this->epsilon_)
1872 {
1873 jmax=My-1+tour_complet;
1874 }
1875 else
1876 for(; coord_noeud(0,jmax,0,1)+this->epsilon_<ymax; jmax++) {};
1877 // if((imax-imin)*(jmax-jmin)==0)
1878 // {
1879 // Cerr << "boundary named " << nom_front << " is empty!" << finl;
1880 // Process::exit();
1881 // }
1882 som.resize((imax-imin)*(jmax-jmin),4);
1883 for(i=imin; i<imax; i++)
1884 for(j=jmin; j<jmax; j++)
1885 {
1886 som((j-jmin)*(imax-imin)+i-imin,0)=numero_sommet(i, j, k);
1887 som((j-jmin)*(imax-imin)+i-imin,1)=numero_sommet(i+1, j, k);
1888 som((j-jmin)*(imax-imin)+i-imin,2)=numero_sommet(i, j+1, k);
1889 som((j-jmin)*(imax-imin)+i-imin,3)=numero_sommet(i+1, j+1, k);
1890 }
1891 }
1892 front.ajouter_faces(som);
1893 front.associer_domaine(*this);
1894 }
1895 else
1896 {
1897 Process::exit();
1898 }
1899 Cerr << "End of reading of the boundary : " << nom_front << finl;
1900}
1901
1902template class Pave_32_64<int>;
1903#if INT_is_64_ == 2
1904template class Pave_32_64<trustIdType>;
1905#endif
1906
Bord_32_64
Class Bord This class represents a boundary of a domain, it is a type of frontier.
Definition Bord.h:31
Bord_Interne_32_64
Class Bord_Interne The class serves to represent a set of faces that are internal.
Definition Bord_Interne.h:32
Domaine_32_64
class Domaine_32_64 A Domain is a mesh composed of a set of geometric elements of the same type.
Definition Domaine.h:62
Domaine_32_64::groupe_faces
Groupe_Faces_t & groupe_faces(int i)
Definition Domaine.h:221
Domaine_32_64::sommets_
DoubleTab_t sommets_
Definition Domaine.h:386
Domaine_32_64::typer
void typer(const Nom &)
Sets the element type of the domain using the name passed as parameter.
Definition Domaine.h:457
Domaine_32_64::mes_elems_
IntTab_t mes_elems_
Definition Domaine.h:391
Domaine_base::epsilon_
double epsilon_
Definition Domaine_base.h:115
Entree
Class defining operators and methods for all reading operation in an input flow (file,...
Definition Entree.h:42
Frontiere_32_64::ajouter_faces
void ajouter_faces(const IntTab_t &)
Adds one or more faces to the boundary; the face(s) are specified by an array.
Definition Frontiere.cpp:86
Frontiere_32_64::nommer
void nommer(const Nom &) override
Gives a name to the boundary.
Definition Frontiere.cpp:74
Frontiere_32_64::dimensionner
void dimensionner(int_t i)
Sizes the boundary, i.e. sets its number of faces.
Definition Frontiere.h:61
Frontiere_32_64::typer_faces
void typer_faces(const Motcle &)
Sets the type of the boundary faces.
Definition Frontiere.cpp:96
Frontiere_32_64::associer_domaine
void associer_domaine(const Domaine_t &)
Associates the boundary to the domain it belongs to.
Definition Frontiere.cpp:63
Frontiere_32_64::add
void add(const Frontiere_32_64 &)
Adds the vertices (and faces) of the given boundary to this object (Frontiere_32_64).
Definition Frontiere.cpp:173
Groupe_Faces_32_64
Groupe_Face class — represents a selection of faces read from a med file.
Definition Groupe_Faces.h:27
Joint_32_64
The Joint class is a Frontiere that contains the joint faces and vertices with the neighboring domain...
Definition Joint.h:34
Motcle
A character string (Nom) in uppercase.
Definition Motcle.h:26
Motcles
An array of Motcle objects.
Definition Motcle.h:63
Nom
class Nom: a character string for naming TRUST objects.
Definition Nom.h:31
Objet_U::Entree
friend class Entree
Definition Objet_U.h:71
Objet_U::dimension
static int dimension
Definition Objet_U.h:94
Objet_U::readOn
virtual Entree & readOn(Entree &)
Reads an Objet_U from an input stream. Virtual method to override.
Definition Objet_U.cpp:289
Objet_U::precision_geom
static double precision_geom
Definition Objet_U.h:81
Objet_U::bidim_axi
static int bidim_axi
Definition Objet_U.h:97
Objet_U::axi
static int axi
Definition Objet_U.h:96
Objet_U::printOn
virtual Sortie & printOn(Sortie &) const
Writes the object to an output stream. Virtual method to override.
Definition Objet_U.cpp:278
Pave_32_64
Class Pave A domain that is particularly easy to mesh!
Definition Pave.h:49
Pave_32_64::a_tanh
double a_tanh
Definition Pave.h:92
Pave_32_64::maille_sommet
int_t & maille_sommet(int i, int l)
Returns a reference to the index of the l-th vertex of the i-th cell (along X) of the pave.
Definition Pave.h:201
Pave_32_64::ztanh_dilatation
int ztanh_dilatation
Definition Pave.h:97
Pave_32_64::tanh_dilatation
int tanh_dilatation
Definition Pave.h:93
Pave_32_64::za_tanh
double za_tanh
Definition Pave.h:96
Pave_32_64::numero_sommet
int_t numero_sommet(int i)
Returns the index of the i-th vertex (along X).
Definition Pave.h:153
Pave_32_64::rep_VEF
bool rep_VEF
Definition Pave.h:98
Pave_32_64::Nx
int Nx
Definition Pave.h:91
Pave_32_64::facteurs_
DoubleVect facteurs_
Definition Pave.h:87
Pave_32_64::tour_complet
bool tour_complet
Definition Pave.h:99
Pave_32_64::Ny
int Ny
Definition Pave.h:91
Pave_32_64::maille3D
void maille3D()
Performs a 3D meshing of the box with the values of the parameters read by ReadOn.
Definition Pave.cpp:844
Pave_32_64::My
int My
Definition Pave.h:91
Pave_32_64::xa_tanh
double xa_tanh
Definition Pave.h:94
Pave_32_64::lire_noeuds
void lire_noeuds(Entree &is)
Reads the number of nodes from the data set input stream and builds the nodes.
Definition Pave.cpp:1342
Pave_32_64::longueurs_
DoubleVect longueurs_
Definition Pave.h:87
Pave_32_64::maille2D
void maille2D()
Performs a 1D meshing of the pave using the parameter values read by ReadOn.
Definition Pave.cpp:580
Pave_32_64::Frontiere_t
Frontiere_32_64< _SIZE_ > Frontiere_t
Definition Pave.h:58
Pave_32_64::Mz
int Mz
Definition Pave.h:91
Pave_32_64::coord_noeud
double & coord_noeud(int i)
Returns a reference to the coordinates of the i-th node.
Definition Pave.h:242
Pave_32_64::Mx
int Mx
Definition Pave.h:91
Pave_32_64::xtanh_dilatation
int xtanh_dilatation
Definition Pave.h:95
Pave_32_64::typer_
void typer_()
Types the pave according to the space dimension and the requested coordinate system (cylindrical.
Definition Pave.cpp:1278
Pave_32_64::IntTab_t
IntTab_T< _SIZE_ > IntTab_t
Definition Pave.h:54
Pave_32_64::origine_
DoubleVect origine_
Definition Pave.h:87
Pave_32_64::lire_front
void lire_front(Entree &, Frontiere_t &)
Reads the specifications of a boundary from the data set input stream and builds it.
Definition Pave.cpp:1456
Pave_32_64::nb_noeuds_
IntVect nb_noeuds_
Definition Pave.h:89
Pave_32_64::symetrique_
IntVect symetrique_
Definition Pave.h:88
Pave_32_64::lire_longueurs
void lire_longueurs(Entree &is)
Reads the lengths LX LY [LZ] from the data set input stream.
Definition Pave.cpp:1315
Pave_32_64::pas_
DoubleVect pas_
Definition Pave.h:87
Pave_32_64::Nz
int Nz
Definition Pave.h:91
Process::exit
static void exit(int exit_code=-1)
Exit routine for TRUST within a Kokkos region.
Definition Process.cpp:466
Raccord_base_32_64
Class Raccord_base This class is simply a boundary; it is the base class of the.
Definition Raccord_base.h:39
Sortie
Base class for output streams.
Definition Sortie.h:52
TRUSTTab::resize
void resize(_SIZE_ n, RESIZE_OPTIONS opt=RESIZE_OPTIONS::COPY_INIT)
Definition TRUSTTab.tpp:469