TRUST 1.9.8
HPC thermohydraulic platform
Loading...
Searching...
No Matches
Iterateur_PolyMAC_CDO_Elem.h
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 Iterateur_PolyMAC_CDO_Elem_included
17#define Iterateur_PolyMAC_CDO_Elem_included
18
19#include <Equation_base.h>
20#include <Milieu_base.h>
21#include <Champ_Uniforme.h>
22#include <Probleme_base.h>
23#include <Schema_Temps_base.h>
24#include <Operateur_base.h>
25#include <Operateur_Diff_base.h>
26#include <Op_Conv_PolyMAC_CDO_iterateur_base.h>
27#include <EcrFicPartage.h>
28#include <communications.h>
29#include <TRUSTTrav.h>
30
31template <class _TYPE_>
32class Iterateur_PolyMAC_CDO_Elem : public Iterateur_PolyMAC_CDO_base
33{
34 inline int duplique() const override
35 {
36 Iterateur_PolyMAC_CDO_Elem* xxx = new Iterateur_PolyMAC_CDO_Elem(*this);
37 if(!xxx) Process::exit("Not enough memory ");
38 return xxx->numero();
39 }
40
41 inline unsigned taille_memoire() const override { throw; }
42
43public:
44 inline Iterateur_PolyMAC_CDO_Elem() { }
45 inline Iterateur_PolyMAC_CDO_Elem(const Iterateur_PolyMAC_CDO_Elem<_TYPE_>& );
46 inline Evaluateur_PolyMAC_CDO& evaluateur() override;
47 inline const Evaluateur_PolyMAC_CDO& evaluateur() const override;
48 DoubleTab& calculer(const DoubleTab& , DoubleTab& ) const override;
49 DoubleTab& ajouter(const DoubleTab&, DoubleTab& ) const override;
50 void calculer_flux_bord(const DoubleTab&) const override;
51 void contribuer_au_second_membre(DoubleTab& ) const override;
52 void ajouter_contribution(const DoubleTab&, Matrice_Morse& ) const override;
53 void ajouter_contribution_vitesse(const DoubleTab&, Matrice_Morse& ) const override;
54 inline void completer_() override;
55 int impr(Sortie& os) const override;
56 void modifier_flux() const;
57
58protected:
59 _TYPE_ flux_evaluateur;
60 DoubleTab& ajouter_bords(const DoubleTab& , DoubleTab& ) const;
61 DoubleTab& ajouter_bords(const DoubleTab& , DoubleTab& , int ) const;
62 DoubleTab& ajouter_interne(const DoubleTab& , DoubleTab& ) const;
63 DoubleTab& ajouter_interne(const DoubleTab& , DoubleTab& , int ) const;
64 void contribuer_au_second_membre_bords(DoubleTab& ) const;
65 void contribuer_au_second_membre_bords(DoubleTab& , int ) const;
66 void contribuer_au_second_membre_interne(DoubleTab& ) const;
67 void contribuer_au_second_membre_interne(DoubleTab&, int ) const;
68 void ajouter_contribution_bords(const DoubleTab&, Matrice_Morse& ) const;
69 void ajouter_contribution_bords(const DoubleTab&, Matrice_Morse&, int ) const;
70 void ajouter_contribution_interne(const DoubleTab&, Matrice_Morse& ) const;
71 void ajouter_contribution_interne(const DoubleTab&, Matrice_Morse&, int ) const;
72 void ajouter_contribution_interne_vitesse(const DoubleTab&, Matrice_Morse& ) const;
73 void ajouter_contribution_bords_vitesse(const DoubleTab&, Matrice_Morse& ) const;
74 const Milieu_base& milieu() const;
75 IntTab elem;
76 mutable SFichier Flux; // Impression .out
77};
78
83
89template <class _TYPE_> inline const Evaluateur_PolyMAC_CDO& Iterateur_PolyMAC_CDO_Elem<_TYPE_>::evaluateur() const
90{
91 Evaluateur_PolyMAC_CDO& eval = (Evaluateur_PolyMAC_CDO&) flux_evaluateur;
92 return eval;
93}
94
95template <class _TYPE_> inline void Iterateur_PolyMAC_CDO_Elem<_TYPE_>::completer_()
96{
97 elem.ref(le_domaine->face_voisins());
98}
99
100template <class _TYPE_>
101const Milieu_base& Iterateur_PolyMAC_CDO_Elem<_TYPE_>::milieu() const
102{
103 return (la_zcl->equation()).milieu();
104}
105
106template <class _TYPE_>
107DoubleTab& Iterateur_PolyMAC_CDO_Elem<_TYPE_>::ajouter(const DoubleTab& donne,
108 DoubleTab& resu) const
109{
110 ((_TYPE_&) flux_evaluateur).mettre_a_jour();
111 assert(donne.nb_dim() < 3);
112 assert(la_zcl);
113 assert(le_domaine);
114 int ncomp=1;
115 if (donne.nb_dim() == 2)
116 ncomp=donne.dimension(1);
117 DoubleTab& flux_bords=op_base->flux_bords();
118 flux_bords.resize(le_domaine->nb_faces_bord(),ncomp);
119 flux_bords=0;
120 /* modif b.m.: on va faire += sur des items virtuels, initialiser les cases */
121 /* sinon risque que les cases soient invalides ou non initialisees */
122 // {
123 // int n = resu.size_array() - resu.size();
124 // double *data = resu.addr() + resu.size();
125 // for (; n; n--, data++)
126 // *data = 0.;
127 // }
128 if( ncomp == 1) /* cas scalaire */
129 {
130 ajouter_bords(donne, resu) ;
131 ajouter_interne(donne, resu) ;
132 }
133 else /* cas vectoriel */
134 {
135 ajouter_bords(donne, resu, ncomp) ;
136 ajouter_interne(donne, resu, ncomp) ;
137 }
138 modifier_flux() ;
139 return resu;
140}
141template <class _TYPE_> DoubleTab& Iterateur_PolyMAC_CDO_Elem<_TYPE_>::ajouter_bords(const DoubleTab& donnee,
142 DoubleTab& resu) const
143{
144 int elem1, elem2;
145 int ndeb, nfin;
146 int face;
147 int num_cl=0;
148 double flux;
149 int nb_front_Cl=le_domaine->nb_front_Cl();
150 DoubleTab& flux_bords=op_base->flux_bords();
151 for (; num_cl<nb_front_Cl; num_cl++)
152 {
153 /* pour chaque Condition Limite on regarde son type */
154 const Cond_lim& la_cl = la_zcl->les_conditions_limites(num_cl);
155 const Front_VF& frontiere_dis = ref_cast(Front_VF,la_cl->frontiere_dis());
156 ndeb = frontiere_dis.num_premiere_face();
157 nfin = ndeb + frontiere_dis.nb_faces();
158 /* Test en bidim axi */
159 if (bidim_axi && !sub_type(Symetrie,la_cl.valeur()))
160 {
161 if (nfin>ndeb && est_egal(le_domaine->face_surfaces()[ndeb],0))
162 {
163 Cerr << "Error in the definition of the boundary conditions." << finl;
164 Cerr << "The axis of revolution for this 2D calculation is along Y." << finl;
165 Cerr << "So you must specify symmetry for the boundary " << frontiere_dis.le_nom() << finl;
166 exit();
167 }
168 }
169 switch(type_cl(la_cl))
170 {
171 case symetrie :
172 if (flux_evaluateur.calculer_flux_faces_symetrie())
173 {
174 const Symetrie& cl =(const Symetrie&) (la_cl.valeur());
175 for (face=ndeb; face<nfin; face++)
176 {
177 flux = flux_evaluateur.flux_face(donnee, face, cl, ndeb);
178 if ( (elem1=elem(face,0)) > -1)
179 {
180 resu[elem1]+=flux;
181 flux_bords(face,0)+=flux;
182 }
183 if ( (elem2=elem(face,1)) > -1)
184 {
185 resu[elem2]-=flux;
186 flux_bords(face,0)-=flux;
187 }
188 }
189 }
190 break;
191 case sortie_libre :
192 if (flux_evaluateur.calculer_flux_faces_sortie_libre())
193 {
194 const Neumann_sortie_libre& cl =(const Neumann_sortie_libre&) (la_cl.valeur());
195 for (face=ndeb; face<nfin; face++)
196 {
197 flux = flux_evaluateur.flux_face(donnee, face, cl, ndeb);
198 if ( (elem1=elem(face,0)) > -1)
199 {
200 resu[elem1]+=flux;
201 flux_bords(face,0)+=flux;
202 }
203 if ( (elem2=elem(face,1)) > -1)
204 {
205 resu[elem2]-=flux;
206 flux_bords(face,0)-=flux;
207 }
208 }
209 }
210 break;
211 case entree_fluide :
212 if (flux_evaluateur.calculer_flux_faces_entree_fluide())
213 {
214 const Dirichlet_entree_fluide& cl =(const Dirichlet_entree_fluide&) (la_cl.valeur());
215 for (face=ndeb; face<nfin; face++)
216 {
217 flux = flux_evaluateur.flux_face(donnee, face, cl, ndeb);
218 if ( (elem1=elem(face,0)) > -1)
219 {
220 resu[elem1]+=flux;
221 flux_bords(face,0)+=flux;
222 }
223 if ( (elem2=elem(face,1)) > -1)
224 {
225 resu[elem2]-=flux;
226 flux_bords(face,0)-=flux;
227 }
228 }
229 }
230 break;
231 case paroi_fixe :
232 if (flux_evaluateur.calculer_flux_faces_paroi_fixe())
233 {
234 const Dirichlet_paroi_fixe& cl =(const Dirichlet_paroi_fixe&) (la_cl.valeur());
235 for (face=ndeb; face<nfin; face++)
236 {
237 flux = flux_evaluateur.flux_face(donnee, face, cl, ndeb);
238 if ( (elem1=elem(face,0)) > -1)
239 {
240 resu[elem1]+=flux;
241 flux_bords(face,0)+=flux;
242 }
243 if ( (elem2=elem(face,1)) > -1)
244 {
245 resu[elem2]-=flux;
246 flux_bords(face,0)-=flux;
247 }
248 }
249 }
250 break;
251 case paroi_defilante :
252 if (flux_evaluateur.calculer_flux_faces_paroi_defilante())
253 {
254 const Dirichlet_paroi_defilante& cl =(const Dirichlet_paroi_defilante&) (la_cl.valeur());
255 for (face=ndeb; face<nfin; face++)
256 {
257 flux = flux_evaluateur.flux_face(donnee, face, cl, ndeb);
258 if ( (elem1=elem(face,0)) > -1)
259 {
260 resu[elem1]+=flux;
261 flux_bords(face,0)+=flux;
262 }
263 if ( (elem2=elem(face,1)) > -1)
264 {
265 resu[elem2]-=flux;
266 flux_bords(face,0)-=flux;
267 }
268 }
269 }
270 break;
271 case paroi_adiabatique :
272 if (flux_evaluateur.calculer_flux_faces_paroi_adiabatique())
273 {
274 const Neumann_paroi_adiabatique& cl =(const Neumann_paroi_adiabatique&) (la_cl.valeur());
275 for (face=ndeb; face<nfin; face++)
276 {
277 /* on initialise elem1 elem2 et on fait planter */
278 elem1=-1;
279 elem2=-1;
280 assert(0);
281 exit();
282 flux = flux_evaluateur.flux_face(donnee, face, cl, ndeb);
283 resu[elem1]+=flux;
284 flux_bords(face,0)+=flux;
285 resu[elem2]-=flux;
286 flux_bords(face,0)-=flux;
287 }
288 }
289 break;
290 case paroi :
291 if (flux_evaluateur.calculer_flux_faces_paroi())
292 {
293 const Neumann_paroi& cl =(const Neumann_paroi&) (la_cl.valeur());
294 for (face=ndeb; face<nfin; face++)
295 {
296 flux = flux_evaluateur.flux_face(donnee, face, cl, ndeb);
297 if ( (elem1=elem(face,0)) > -1)
298 {
299 resu[elem1]+=flux;
300 flux_bords(face,0)+=flux;
301 }
302 if ( (elem2=elem(face,1)) > -1)
303 {
304 resu[elem2]-=flux;
305 flux_bords(face,0)-=flux;
306 }
307 }
308 }
309 break;
310 case echange_externe_impose :
311 if (flux_evaluateur.calculer_flux_faces_echange_externe_impose())
312 {
313 const Echange_externe_impose& cl =(const Echange_externe_impose&) (la_cl.valeur());
314
315 int boundary_index=-1;
316 if (le_domaine->front_VF(num_cl).le_nom() == frontiere_dis.le_nom())
317 boundary_index=num_cl;
318
319 for (face=ndeb; face<nfin; face++)
320 {
321 int local_face=le_domaine->front_VF(boundary_index).num_local_face(face);
322 flux = flux_evaluateur.flux_face(donnee, boundary_index,face,local_face, cl, ndeb);
323 if ( (elem1=elem(face,0)) > -1)
324 {
325 resu[elem1]+=flux;
326 flux_bords(face,0)+=flux;
327 }
328 if ( (elem2=elem(face,1)) > -1)
329 {
330 resu[elem2]-=flux;
331 flux_bords(face,0)-=flux;
332 }
333 }
334 }
335 break;
336 case echange_global_impose :
337 if (flux_evaluateur.calculer_flux_faces_echange_global_impose())
338 {
339 const Echange_global_impose& cl =(const Echange_global_impose&) (la_cl.valeur());
340 for (face=ndeb; face<nfin; face++)
341 {
342 flux = flux_evaluateur.flux_face(donnee, face, cl, ndeb);
343 if ( (elem1=elem(face,0)) > -1)
344 {
345 resu[elem1]+=flux;
346 flux_bords(face,0)+=flux;
347 }
348 if ( (elem2=elem(face,1)) > -1)
349 {
350 resu[elem2]-=flux;
351 flux_bords(face,0)-=flux;
352 }
353 }
354 }
355 break;
356 case periodique :
357 if (flux_evaluateur.calculer_flux_faces_periodique())
358 {
359 const Periodique& cl =(const Periodique&) (la_cl.valeur());
360 for (face=ndeb; face<nfin; face++)
361 {
362 flux = flux_evaluateur.flux_face(donnee, face, cl, ndeb);
363 if ( (elem1=elem(face,0)) > -1)
364 {
365 resu[elem1]+=0.5*flux;
366 if ( face < (ndeb+frontiere_dis.nb_faces()/2) )
367 flux_bords(face,0)+=flux;
368 }
369 if ( (elem2=elem(face,1)) > -1)
370 {
371 resu[elem2]-=0.5*flux;
372 if ( (ndeb+frontiere_dis.nb_faces()/2) <= face )
373 flux_bords(face,0)-=flux;
374 }
375 }
376 }
377 break;
378 case scalaire_impose_paroi :
379 break;
380 /*
381 case nouvelle_Cl_PolyMAC_CDO :
382 if (flux_evaluateur.calculer_flux_faces_echange_global_impose()){
383 const Nouvelle_Cl_PolyMAC_CDO& cl =(const Nouvelle_Cl_PolyMAC_CDO&) (la_cl.valeur());
384 for (face=ndeb; face<nfin; face++) {
385 if ( (elem1=elem(face,0)) > -1)
386 resu[elem1]+=flux_evaluateur.flux_face(donnee, face, cl, ndeb);
387 if ( (elem2=elem(face,1)) > -1)
388 resu[elem2]-=flux_evaluateur.flux_face(donnee, face, cl, ndeb);
389 }
390 }
391 break;
392 */
393 default :
394 Cerr << "On ne reconnait pas la condition limite : " << la_cl.valeur();
395 Cerr << "Dans Iterateur_PolyMAC_CDO_Elem<_TYPE_>::ajouter_bords"<<finl;
396 exit();
397 break;
398 }
399 }
400 return resu;
401}
402template <class _TYPE_> void Iterateur_PolyMAC_CDO_Elem<_TYPE_>::calculer_flux_bord(const DoubleTab& donnee) const
403{
404 ((_TYPE_&) flux_evaluateur).mettre_a_jour();
405 assert(donnee.nb_dim() < 3);
406 assert(la_zcl);
407 assert(le_domaine);
408 int ncomp=1;
409 if (donnee.nb_dim() == 2)
410 ncomp=donnee.dimension(1);
411 DoubleTab& flux_bords=op_base->flux_bords();
412 flux_bords.resize(le_domaine->nb_faces_bord(),ncomp);
413 flux_bords=0;
414
415 if( ncomp != 1) Process ::exit();/* cas scalaire */
416
417
418 int ndeb, nfin;
419 int face;
420 int num_cl=0;
421 double flux;
422 int nb_front_Cl=le_domaine->nb_front_Cl();
423 for (; num_cl<nb_front_Cl; num_cl++)
424 {
425 /* pour chaque Condition Limite on regarde son type */
426 const Cond_lim& la_cl = la_zcl->les_conditions_limites(num_cl);
427 const Front_VF& frontiere_dis = ref_cast(Front_VF,la_cl->frontiere_dis());
428 ndeb = frontiere_dis.num_premiere_face();
429 nfin = ndeb + frontiere_dis.nb_faces();
430 /* Test en bidim axi */
431 if (bidim_axi && !sub_type(Symetrie,la_cl.valeur()))
432 {
433 if (nfin>ndeb && est_egal(le_domaine->face_surfaces()[ndeb],0))
434 {
435 Cerr << "Error in the definition of the boundary conditions." << finl;
436 Cerr << "The axis of revolution for this 2D calculation is along Y." << finl;
437 Cerr << "So you must specify symmetry for the boundary " << frontiere_dis.le_nom() << finl;
438 exit();
439 }
440 }
441 switch(type_cl(la_cl))
442 {
443 case symetrie :
444 if (flux_evaluateur.calculer_flux_faces_symetrie())
445 {
446 const Symetrie& cl =(const Symetrie&) (la_cl.valeur());
447 for (face=ndeb; face<nfin; face++)
448 {
449 flux = flux_evaluateur.flux_face(donnee, face, cl, ndeb);
450 if ( (elem(face,0)) > -1)
451 {
452 flux_bords(face,0)+=flux;
453 }
454 if ( (elem(face,1)) > -1)
455 {
456 flux_bords(face,0)-=flux;
457 }
458 }
459 }
460 break;
461 case sortie_libre :
462 if (flux_evaluateur.calculer_flux_faces_sortie_libre())
463 {
464 const Neumann_sortie_libre& cl =(const Neumann_sortie_libre&) (la_cl.valeur());
465 for (face=ndeb; face<nfin; face++)
466 {
467 flux = flux_evaluateur.flux_face(donnee, face, cl, ndeb);
468 if ( (elem(face,0)) > -1)
469 {
470 flux_bords(face,0)+=flux;
471 }
472 if ( (elem(face,1)) > -1)
473 {
474 flux_bords(face,0)-=flux;
475 }
476 }
477 }
478 break;
479 case entree_fluide :
480 if (flux_evaluateur.calculer_flux_faces_entree_fluide())
481 {
482 const Dirichlet_entree_fluide& cl =(const Dirichlet_entree_fluide&) (la_cl.valeur());
483 for (face=ndeb; face<nfin; face++)
484 {
485 flux = flux_evaluateur.flux_face(donnee, face, cl, ndeb);
486 if ( (elem(face,0)) > -1)
487 {
488 flux_bords(face,0)+=flux;
489 }
490 if ( (elem(face,1)) > -1)
491 {
492 flux_bords(face,0)-=flux;
493 }
494 }
495 }
496 break;
497 case paroi_fixe :
498 if (flux_evaluateur.calculer_flux_faces_paroi_fixe())
499 {
500 const Dirichlet_paroi_fixe& cl =(const Dirichlet_paroi_fixe&) (la_cl.valeur());
501 for (face=ndeb; face<nfin; face++)
502 {
503 flux = flux_evaluateur.flux_face(donnee, face, cl, ndeb);
504 if ( (elem(face,0)) > -1)
505 {
506 flux_bords(face,0)+=flux;
507 }
508 if ( (elem(face,1)) > -1)
509 {
510 flux_bords(face,0)-=flux;
511 }
512 }
513 }
514 break;
515 case paroi_defilante :
516 if (flux_evaluateur.calculer_flux_faces_paroi_defilante())
517 {
518 const Dirichlet_paroi_defilante& cl =(const Dirichlet_paroi_defilante&) (la_cl.valeur());
519 for (face=ndeb; face<nfin; face++)
520 {
521 flux = flux_evaluateur.flux_face(donnee, face, cl, ndeb);
522 if ( (elem(face,0)) > -1)
523 {
524 flux_bords(face,0)+=flux;
525 }
526 if ( (elem(face,1)) > -1)
527 {
528 flux_bords(face,0)-=flux;
529 }
530 }
531 }
532 break;
533 case paroi_adiabatique :
534 if (flux_evaluateur.calculer_flux_faces_paroi_adiabatique())
535 {
536 const Neumann_paroi_adiabatique& cl =(const Neumann_paroi_adiabatique&) (la_cl.valeur());
537 for (face=ndeb; face<nfin; face++)
538 {
539 assert(0);
540 exit();
541 flux = flux_evaluateur.flux_face(donnee, face, cl, ndeb);
542 flux_bords(face,0)+=flux;
543 flux_bords(face,0)-=flux;
544 }
545 }
546 break;
547 case paroi :
548 if (flux_evaluateur.calculer_flux_faces_paroi())
549 {
550 const Neumann_paroi& cl =(const Neumann_paroi&) (la_cl.valeur());
551 for (face=ndeb; face<nfin; face++)
552 {
553 flux = flux_evaluateur.flux_face(donnee, face, cl, ndeb);
554 if ( (elem(face,0)) > -1)
555 {
556 flux_bords(face,0)+=flux;
557 }
558 if ( (elem(face,1)) > -1)
559 {
560 flux_bords(face,0)-=flux;
561 }
562 }
563 }
564 break;
565 case echange_externe_impose :
566 if (flux_evaluateur.calculer_flux_faces_echange_externe_impose())
567 {
568 const Echange_externe_impose& cl =(const Echange_externe_impose&) (la_cl.valeur());
569
570 int boundary_index=-1;
571 if (le_domaine->front_VF(num_cl).le_nom() == frontiere_dis.le_nom())
572 boundary_index=num_cl;
573
574 for (face=ndeb; face<nfin; face++)
575 {
576 int local_face=le_domaine->front_VF(boundary_index).num_local_face(face);
577 flux = flux_evaluateur.flux_face(donnee, boundary_index,face,local_face, cl, ndeb);
578 if ( (elem(face,0)) > -1)
579 {
580 flux_bords(face,0)+=flux;
581 }
582 if ( (elem(face,1)) > -1)
583 {
584 flux_bords(face,0)-=flux;
585 }
586 }
587 }
588 break;
589 case echange_global_impose :
590 if (flux_evaluateur.calculer_flux_faces_echange_global_impose())
591 {
592 const Echange_global_impose& cl =(const Echange_global_impose&) (la_cl.valeur());
593 for (face=ndeb; face<nfin; face++)
594 {
595 flux = flux_evaluateur.flux_face(donnee, face, cl, ndeb);
596 if ( (elem(face,0)) > -1)
597 {
598 flux_bords(face,0)+=flux;
599 }
600 if ( (elem(face,1)) > -1)
601 {
602 flux_bords(face,0)-=flux;
603 }
604 }
605 }
606 break;
607 case periodique :
608 if (flux_evaluateur.calculer_flux_faces_periodique())
609 {
610 const Periodique& cl =(const Periodique&) (la_cl.valeur());
611 for (face=ndeb; face<nfin; face++)
612 {
613 flux = flux_evaluateur.flux_face(donnee, face, cl, ndeb);
614 if ( (elem(face,0)) > -1)
615 {
616 if ( face < (ndeb+frontiere_dis.nb_faces()/2) )
617 flux_bords(face,0)+=flux;
618 }
619 if ( (elem(face,1)) > -1)
620 {
621 if ( (ndeb+frontiere_dis.nb_faces()/2) <= face )
622 flux_bords(face,0)-=flux;
623 }
624 }
625 }
626 break;
627 case scalaire_impose_paroi :
628 break;
629 /*
630 case nouvelle_Cl_PolyMAC_CDO :
631 if (flux_evaluateur.calculer_flux_faces_echange_global_impose()){
632 const Nouvelle_Cl_PolyMAC_CDO& cl =(const Nouvelle_Cl_PolyMAC_CDO&) (la_cl.valeur());
633 for (face=ndeb; face<nfin; face++) {
634 if ( (int elem1=elem(face,0)) > -1)
635 resu[elem1]+=flux_evaluateur.flux_face(donnee, face, cl, ndeb);
636 if ( (int elem2=elem(face,1)) > -1)
637 resu[elem2]-=flux_evaluateur.flux_face(donnee, face, cl, ndeb);
638 }
639 }
640 break;
641 */
642 default :
643 Cerr << "On ne reconnait pas la condition limite : " << la_cl.valeur();
644 Cerr << "Dans Iterateur_PolyMAC_CDO_Elem<_TYPE_>::ajouter_bords"<<finl;
645 exit();
646 break;
647 }
648 }
649 modifier_flux() ;
650}
651
652template <class _TYPE_> DoubleTab& Iterateur_PolyMAC_CDO_Elem<_TYPE_>::ajouter_bords(const DoubleTab& donnee,
653 DoubleTab& resu,int ncomp) const
654{
655 int elem1, elem2;
656 int ndeb, nfin;
657 int face,k;
658 DoubleVect flux(ncomp);
659 int num_cl=0;
660 int nb_front_Cl=le_domaine->nb_front_Cl();
661 DoubleTab& flux_bords=op_base->flux_bords();
662 for (; num_cl<nb_front_Cl; num_cl++)
663 {
664 const Cond_lim& la_cl = la_zcl->les_conditions_limites(num_cl);
665 const Front_VF& frontiere_dis = ref_cast(Front_VF,la_cl->frontiere_dis());
666 ndeb = frontiere_dis.num_premiere_face();
667 nfin = ndeb + frontiere_dis.nb_faces();
668 /* Test en bidim axi */
669 if (bidim_axi && !sub_type(Symetrie,la_cl.valeur()))
670 {
671 if (nfin>ndeb && est_egal(le_domaine->face_surfaces()[ndeb],0))
672 {
673 Cerr << "Error in the definition of the boundary conditions." << finl;
674 Cerr << "The axis of revolution for this 2D calculation is along Y." << finl;
675 Cerr << "So you must specify symmetry for the boundary " << frontiere_dis.le_nom() << finl;
676 exit();
677 }
678 }
679 switch(type_cl(la_cl))
680 {
681 case symetrie :
682 if (flux_evaluateur.calculer_flux_faces_symetrie())
683 {
684 const Symetrie& cl =(const Symetrie&) (la_cl.valeur());
685 for (face=ndeb; face<nfin; face++)
686 {
687 flux_evaluateur.flux_face(donnee, face, cl, ndeb, flux);
688 if ( (elem1=elem(face,0)) > -1)
689 for (k=0; k<ncomp; k++)
690 {
691 resu(elem1,k) +=flux(k);
692 flux_bords(face,k)+=flux(k);
693 }
694 if ( (elem2=elem(face,1)) > -1)
695 for (k=0; k<ncomp; k++)
696 {
697 resu(elem2,k) -=flux(k);
698 flux_bords(face,k)-=flux(k);
699 }
700 }
701 }
702 break;
703 case sortie_libre :
704 if (flux_evaluateur.calculer_flux_faces_sortie_libre())
705 {
706 const Neumann_sortie_libre& cl =(const Neumann_sortie_libre&) (la_cl.valeur());
707 for (face=ndeb; face<nfin; face++)
708 {
709 flux_evaluateur.flux_face(donnee, face, cl, ndeb, flux);
710 if ( (elem1=elem(face,0)) > -1)
711 for (k=0; k<ncomp; k++)
712 {
713 resu(elem1,k) +=flux(k);
714 flux_bords(face,k)+=flux(k);
715 }
716 if ( (elem2=elem(face,1)) > -1)
717 for (k=0; k<ncomp; k++)
718 {
719 resu(elem2,k) -=flux(k);
720 flux_bords(face,k)-=flux(k);
721 }
722 }
723 }
724 break;
725 case entree_fluide :
726 if (flux_evaluateur.calculer_flux_faces_entree_fluide())
727 {
728 const Dirichlet_entree_fluide& cl =(const Dirichlet_entree_fluide&) (la_cl.valeur());
729 for (face=ndeb; face<nfin; face++)
730 {
731 flux_evaluateur.flux_face(donnee, face, cl, ndeb, flux);
732 if ( (elem1=elem(face,0)) > -1)
733 for (k=0; k<ncomp; k++)
734 {
735 resu(elem1,k) +=flux(k);
736 flux_bords(face,k)+=flux(k);
737 }
738 if ( (elem2=elem(face,1)) > -1)
739 for (k=0; k<ncomp; k++)
740 {
741 resu(elem2,k) -=flux(k);
742 flux_bords(face,k)-=flux(k);
743 }
744 }
745 }
746 break;
747 case paroi_fixe :
748 if (flux_evaluateur.calculer_flux_faces_paroi_fixe())
749 {
750 const Dirichlet_paroi_fixe& cl =(const Dirichlet_paroi_fixe&) (la_cl.valeur());
751 for (face=ndeb; face<nfin; face++)
752 {
753 flux_evaluateur.flux_face(donnee, face, cl, ndeb, flux);
754 if ( (elem1=elem(face,0)) > -1)
755 for (k=0; k<ncomp; k++)
756 {
757 resu(elem1,k) +=flux(k);
758 flux_bords(face,k)+=flux(k);
759 }
760 if ( (elem2=elem(face,1)) > -1)
761 for (k=0; k<ncomp; k++)
762 {
763 resu(elem2,k) -=flux(k);
764 flux_bords(face,k)-=flux(k);
765 }
766 }
767 }
768 break;
769 case paroi_defilante :
770 if (flux_evaluateur.calculer_flux_faces_paroi_defilante())
771 {
772 const Dirichlet_paroi_defilante& cl =(const Dirichlet_paroi_defilante&) (la_cl.valeur());
773 for (face=ndeb; face<nfin; face++)
774 {
775 flux_evaluateur.flux_face(donnee, face, cl, ndeb, flux);
776 if ( (elem1=elem(face,0)) > -1)
777 for (k=0; k<ncomp; k++)
778 {
779 resu(elem1,k) +=flux(k);
780 flux_bords(face,k)+=flux(k);
781 }
782 if ( (elem2=elem(face,1)) > -1)
783 for (k=0; k<ncomp; k++)
784 {
785 resu(elem2,k) -=flux(k);
786 flux_bords(face,k)-=flux(k);
787 }
788 }
789 }
790 break;
791 case paroi_adiabatique :
792 if (flux_evaluateur.calculer_flux_faces_paroi_adiabatique())
793 {
794 const Neumann_paroi_adiabatique& cl =(const Neumann_paroi_adiabatique&) (la_cl.valeur());
795 for (face=ndeb; face<nfin; face++)
796 {
797 flux_evaluateur.flux_face(donnee, face, cl, ndeb, flux);
798 if ( (elem1=elem(face,0)) > -1)
799 for (k=0; k<ncomp; k++)
800 {
801 resu(elem1,k) +=flux(k);
802 flux_bords(face,k)+=flux(k);
803 }
804 if ( (elem2=elem(face,1)) > -1)
805 for (k=0; k<ncomp; k++)
806 {
807 resu(elem2,k) -=flux(k);
808 flux_bords(face,k)-=flux(k);
809 }
810 }
811 }
812 break;
813 case paroi :
814 if (flux_evaluateur.calculer_flux_faces_paroi())
815 {
816 const Neumann_paroi& cl =(const Neumann_paroi&) (la_cl.valeur());
817 for (face=ndeb; face<nfin; face++)
818 {
819 flux_evaluateur.flux_face(donnee, face, cl, ndeb, flux);
820 /* on initialise elem1 elem2 et on fait planter */
821 elem1=-1;
822 elem2=-1;
823 assert(0);
824 exit();
825 for (k=0; k<ncomp; k++)
826 {
827 resu(elem1,k) +=flux(k);
828 flux_bords(face,k)+=flux(k);
829 }
830 for (k=0; k<ncomp; k++)
831 {
832 resu(elem2,k) -=flux(k);
833 flux_bords(face,k)-=flux(k);
834 }
835 }
836 }
837 break;
838 case echange_externe_impose :
839 if (flux_evaluateur.calculer_flux_faces_echange_externe_impose())
840 {
841 const Echange_externe_impose& cl =(const Echange_externe_impose&) (la_cl.valeur());
842
843 int boundary_index=-1;
844 if (le_domaine->front_VF(num_cl).le_nom() == frontiere_dis.le_nom())
845 boundary_index=num_cl;
846
847 for (face=ndeb; face<nfin; face++)
848 {
849 int local_face=le_domaine->front_VF(boundary_index).num_local_face(face);
850 flux_evaluateur.flux_face(donnee, boundary_index, face, local_face, cl, ndeb, flux);
851 if ( (elem1=elem(face,0)) > -1)
852 for (k=0; k<ncomp; k++)
853 {
854 resu(elem1,k) +=flux(k);
855 flux_bords(face,k)+=flux(k);
856 }
857 if ( (elem2=elem(face,1)) > -1)
858 for (k=0; k<ncomp; k++)
859 {
860 resu(elem2,k) -=flux(k);
861 flux_bords(face,k)-=flux(k);
862 }
863 }
864 }
865 break;
866 case echange_global_impose :
867 if (flux_evaluateur.calculer_flux_faces_echange_global_impose())
868 {
869 const Echange_global_impose& cl =(const Echange_global_impose&) (la_cl.valeur());
870 for (face=ndeb; face<nfin; face++)
871 {
872 flux_evaluateur.flux_face(donnee, face, cl, ndeb, flux);
873 if ( (elem1=elem(face,0)) > -1)
874 for (k=0; k<ncomp; k++)
875 {
876 resu(elem1,k) +=flux(k);
877 flux_bords(face,k)+=flux(k);
878 }
879 if ( (elem2=elem(face,1)) > -1)
880 for (k=0; k<ncomp; k++)
881 {
882 resu(elem2,k) -=flux(k);
883 flux_bords(face,k)-=flux(k);
884 }
885 }
886 }
887 break;
888 case periodique :
889 if (flux_evaluateur.calculer_flux_faces_periodique())
890 {
891 const Periodique& cl =(const Periodique&) (la_cl.valeur());
892 for (face=ndeb; face<nfin; face++)
893 {
894 flux_evaluateur.flux_face(donnee, face, cl, ndeb, flux);
895 if ( (elem1=elem(face,0)) > -1)
896 for (k=0; k<ncomp; k++)
897 {
898 resu(elem1,k) +=0.5*flux(k);
899 if ( face < (ndeb+frontiere_dis.nb_faces()/2) )
900 flux_bords(face,k)+=flux(k);
901 }
902 if ( (elem2=elem(face,1)) > -1)
903 for (k=0; k<ncomp; k++)
904 {
905 resu(elem2,k) -=0.5*flux(k);
906 if ( (ndeb+frontiere_dis.nb_faces()/2) <= face )
907 flux_bords(face,k)-=flux(k);
908 }
909 }
910 }
911 break;
912 /*
913 case nouvelle_Cl_PolyMAC_CDO :
914 if (flux_evaluateur.calculer_flux_faces_echange_global_impose()){
915 const Nouvelle_Cl_PolyMAC_CDO& cl =(const Nouvelle_Cl_PolyMAC_CDO&) (la_cl.valeur());
916 for (face=ndeb; face<nfin; face++) {
917 flux_evaluateur.flux_face(donnee, face, cl, ndeb, flux);
918 if ( (elem1=elem(face,0)) > -1)
919 for (k=0; k<ncomp; k++)
920 resu(elem1,k) +=flux(k);
921 if ( (elem2=elem(face,1)) > -1)
922 for (k=0; k<ncomp; k++)
923 resu(elem2,k) -=flux(k);
924 }
925 }
926 break;
927 */
928 default :
929 Cerr << "On ne reconnait pas la condition limite : " << la_cl.valeur();
930 exit();
931 break;
932 }
933 }
934 return resu;
935}
936
937template <class _TYPE_> DoubleTab& Iterateur_PolyMAC_CDO_Elem<_TYPE_>::ajouter_interne(const DoubleTab& donnee,
938 DoubleTab& resu) const
939{
940 const Domaine_PolyMAC_CDO& domaine_PolyMAC_CDO = le_domaine.valeur();
941 double flux;
942 int face;
943 int ndeb = domaine_PolyMAC_CDO.premiere_face_int();
944 int nfin = domaine_PolyMAC_CDO.nb_faces();
945 for (face=ndeb; face<nfin; face++)
946 {
947 flux=flux_evaluateur.flux_faces_interne(donnee, face);
948 resu[elem(face,0)]+=flux;
949 resu[elem(face,1)]-=flux;
950 }
951 return resu;
952}
953template <class _TYPE_> DoubleTab& Iterateur_PolyMAC_CDO_Elem<_TYPE_>::ajouter_interne(const DoubleTab& donnee,
954 DoubleTab& resu,int ncomp) const
955{
956 const Domaine_PolyMAC_CDO& domaine_PolyMAC_CDO = le_domaine.valeur();
957 DoubleVect flux(ncomp);
958 int face,k;
959 int elem0,elem1;
960 int ndeb = domaine_PolyMAC_CDO.premiere_face_int();
961 int nfin = domaine_PolyMAC_CDO.nb_faces();
962 for (face=ndeb; face<nfin; face++)
963 {
964 flux_evaluateur.flux_faces_interne(donnee, face, flux);
965 elem0 = elem(face,0);
966 elem1 = elem(face,1);
967 for (k=0; k<ncomp; k++)
968 {
969 resu(elem0,k)+=flux(k);
970 resu(elem1,k)-=flux(k);
971 }
972 }
973 return resu;
974}
975template <class _TYPE_> DoubleTab& Iterateur_PolyMAC_CDO_Elem<_TYPE_>::calculer(const DoubleTab& inco, DoubleTab& resu) const
976{
977 operator_egal(resu, 0., VECT_REAL_ITEMS);
978 return ajouter(inco,resu);
979}
980template <class _TYPE_> void Iterateur_PolyMAC_CDO_Elem<_TYPE_>::modifier_flux() const
981{
982 if (op_base->equation().inconnue().le_nom()=="temperature"
983 && !( sub_type(Operateur_Diff_base,op_base.valeur()) && ref_cast(Operateur_Diff_base,op_base.valeur()).diffusivite().le_nom() == "conductivite" ) )
984 {
985 DoubleTab& flux_bords=op_base->flux_bords();
986 const Domaine_PolyMAC_CDO& le_domaine_vdf=ref_cast(Domaine_PolyMAC_CDO,op_base->equation().domaine_dis());
987 const Champ_base& rho = (op_base->equation()).milieu().masse_volumique();
988 const Champ_Don_base& Cp = (op_base->equation()).milieu().capacite_calorifique();
989 const IntTab& face_voisins=le_domaine_vdf.face_voisins();
990 int rho_uniforme=(sub_type(Champ_Uniforme,rho) ? 1:0);
991 int cp_uniforme=(sub_type(Champ_Uniforme,Cp) ? 1:0);
992 int is_rho_u=op_base->equation().probleme().is_dilatable();
993 if (is_rho_u)
994 {
995 const Operateur_base& op=op_base.valeur();
996 is_rho_u=0;
997 if (sub_type(Op_Conv_PolyMAC_CDO_iterateur_base,op))
998 if (ref_cast(Op_Conv_PolyMAC_CDO_iterateur_base,op).vitesse().le_nom()=="rho_u")
999 is_rho_u=1;
1000 }
1001 double Cp_=0,rho_=0;
1002 const int nb_faces_bords=le_domaine_vdf.nb_faces_bord();
1003 for (int face = 0; face < nb_faces_bords; face++)
1004 {
1005 int num_elem = face_voisins(face, 0);
1006 if (num_elem == -1)
1007 num_elem = face_voisins(face, 1);
1008 if (cp_uniforme)
1009 Cp_ = Cp.valeurs()(0, 0);
1010 else if (Cp.nb_comp() == 1)
1011 Cp_ = Cp.valeurs()(num_elem);
1012 else
1013 Cp_ = Cp.valeurs()(num_elem, 0);
1014 if (rho_uniforme)
1015 rho_ = rho.valeurs()(0, 0);
1016 else if (rho.nb_comp() == 1)
1017 rho_ = rho.valeurs()(num_elem);
1018 else
1019 rho_ = rho.valeurs()(num_elem, 0);
1020 /* si on est en QC temperature on a calcule div(rhou * T) */
1021 /* il ne faut pas remultiplier par rho */
1022 if (is_rho_u)
1023 rho_ = 1;
1024 flux_bords(face, 0) *= (rho_ * Cp_);
1025 }
1026 }
1027}
1028
1029template <class _TYPE_> int Iterateur_PolyMAC_CDO_Elem<_TYPE_>::impr(Sortie& os) const
1030{
1031 const Domaine& madomaine=le_domaine->domaine();
1032 const Domaine_PolyMAC_CDO& zpoly=ref_cast(Domaine_PolyMAC_CDO,op_base->equation().domaine_dis());
1033 const int impr_bord=(madomaine.bords_a_imprimer().est_vide() ? 0:1);
1034 const Schema_Temps_base& sch = la_zcl->equation().probleme().schema_temps();
1035 double temps = sch.temps_courant();
1036 const DoubleTab& inco = la_zcl->equation().inconnue().valeurs();
1037 DoubleTab& flux_bords=op_base->flux_bords();
1038 int ncomp=1;
1039 if (inco.nb_dim() == 2)
1040 ncomp=inco.dimension(1);
1041 flux_bords.resize(zpoly.premiere_face_int(),ncomp);
1042 DoubleVect bilan(flux_bords.dimension(1));
1043 int k,face;
1044 int nb_front_Cl=le_domaine->nb_front_Cl();
1045 DoubleTrav flux_bords2( 3, nb_front_Cl , flux_bords.dimension(1)) ;
1046 flux_bords2=0;
1047 /*flux_bord(k) -> flux_bords2(0,num_cl,k) */
1048 /*flux_bord_perio1(k) -> flux_bords2(1,num_cl,k) */
1049 /*flux_bord_perio2(k) -> flux_bords2(2,num_cl,k) */
1050 for (int num_cl=0; num_cl<nb_front_Cl; num_cl++)
1051 {
1052 const Cond_lim& la_cl = la_zcl->les_conditions_limites(num_cl);
1053 const Front_VF& frontiere_dis = ref_cast(Front_VF,la_cl->frontiere_dis());
1054 int ndeb = frontiere_dis.num_premiere_face();
1055 int nfin = ndeb + frontiere_dis.nb_faces();
1056 int periodicite = (type_cl(la_cl)==periodique?1:0);
1057 for (face=ndeb; face<nfin; face++)
1058 for(k=0; k<flux_bords.dimension(1); k++)
1059 {
1060 flux_bords2(0,num_cl,k)+=flux_bords(face, k);
1061 if(periodicite)
1062 {
1063 if( face < (ndeb+frontiere_dis.nb_faces()/2) )
1064 flux_bords2(1,num_cl,k)+=flux_bords(face, k);
1065 else
1066 flux_bords2(2,num_cl,k)+=flux_bords(face, k);
1067 }
1068 }
1069 } /* fin for num_cl */
1070 mp_sum_for_each_item(flux_bords2);
1071 if (je_suis_maitre())
1072 {
1073 op_base->ouvrir_fichier(Flux,"",1);
1074 Flux.add_col(temps);
1075 for (int num_cl=0; num_cl<nb_front_Cl; num_cl++)
1076 {
1077 const Cond_lim& la_cl = la_zcl->les_conditions_limites(num_cl);
1078 int periodicite = (type_cl(la_cl)==periodique?1:0);
1079 for(k=0; k<flux_bords.dimension(1); k++)
1080 {
1081 bilan(k)+=flux_bords2(0,num_cl,k);
1082 if(periodicite)
1083 {
1084 Flux.add_col(flux_bords2(1,num_cl,k));
1085 Flux.add_col(flux_bords2(2,num_cl,k));
1086 }
1087 else
1088 Flux.add_col(flux_bords2(0,num_cl,k));
1089 }
1090 }
1091 for(k=0; k<flux_bords.dimension(1); k++)
1092 Flux.add_col(bilan(k));
1093 Flux << finl;
1094 }
1095 const LIST(Nom)& Liste_Bords_a_imprimer = le_domaine->domaine().bords_a_imprimer();
1096 if (!Liste_Bords_a_imprimer.est_vide())
1097 {
1098 EcrFicPartage Flux_face;
1099 op_base->ouvrir_fichier_partage(Flux_face,"",impr_bord);
1100 for (int num_cl=0; num_cl<nb_front_Cl; num_cl++)
1101 {
1102 const Frontiere_dis_base& la_fr = la_zcl->les_conditions_limites(num_cl)->frontiere_dis();
1103 const Cond_lim& la_cl = la_zcl->les_conditions_limites(num_cl);
1104 const Front_VF& frontiere_dis = ref_cast(Front_VF,la_cl->frontiere_dis());
1105 int ndeb = frontiere_dis.num_premiere_face();
1106 int nfin = ndeb + frontiere_dis.nb_faces();
1107 if (madomaine.bords_a_imprimer().contient(la_fr.le_nom()))
1108 {
1109 Flux_face << "# Flux par face sur " << la_fr.le_nom() << " au temps " << temps << " : " << finl;
1110 for (face=ndeb; face<nfin; face++)
1111 {
1112 if (dimension == 2)
1113 Flux_face << "# Face a x= " << le_domaine->xv(face,0) << " y= " << le_domaine->xv(face,1) << " : ";
1114 else if (dimension == 3)
1115 Flux_face << "# Face a x= " << le_domaine->xv(face,0) << " y= " << le_domaine->xv(face,1) << " z= " << le_domaine->xv(face,2) << " : ";
1116 for(k=0; k<flux_bords.dimension(1); k++)
1117 Flux_face << flux_bords(face, k) << " ";
1118 Flux_face << finl;
1119 }
1120 Flux_face.syncfile();
1121 }
1122 }
1123 }
1124 return 1;
1125}
1126
1127template <class _TYPE_> void Iterateur_PolyMAC_CDO_Elem<_TYPE_>::contribuer_au_second_membre(DoubleTab& resu) const
1128{
1129 ((_TYPE_&) flux_evaluateur).mettre_a_jour();
1130 assert(resu.nb_dim() < 3);
1131 assert(la_zcl);
1132 assert(le_domaine);
1133 int ncomp=1;
1134 if (resu.nb_dim() == 2)
1135 ncomp=resu.dimension(1);
1136 assert(op_base->flux_bords().dimension(0)==le_domaine->nb_faces_bord()); /* resize deja fait */
1137 if( ncomp == 1) /* cas scalaire */
1138 {
1139 contribuer_au_second_membre_bords(resu) ;
1140 contribuer_au_second_membre_interne(resu) ;
1141 }
1142 else /* cas vectoriel */
1143 {
1144 contribuer_au_second_membre_bords(resu, ncomp) ;
1145 contribuer_au_second_membre_interne(resu, ncomp) ;
1146 }
1147}
1148template <class _TYPE_> void Iterateur_PolyMAC_CDO_Elem<_TYPE_>::contribuer_au_second_membre_bords(DoubleTab& resu) const
1149{
1150 int elem1, elem2;
1151 int ndeb, nfin;
1152 int face;
1153 int num_cl=0;
1154 double flux;
1155 int nb_front_Cl=le_domaine->nb_front_Cl();
1156 DoubleTab& flux_bords=op_base->flux_bords();
1157 for (; num_cl<nb_front_Cl; num_cl++)
1158 {
1159 /* pour chaque Condition Limite on regarde son type */
1160 const Cond_lim& la_cl = la_zcl->les_conditions_limites(num_cl);
1161 const Front_VF& frontiere_dis = ref_cast(Front_VF,la_cl->frontiere_dis());
1162 ndeb = frontiere_dis.num_premiere_face();
1163 nfin = ndeb + frontiere_dis.nb_faces();
1164 switch(type_cl(la_cl))
1165 {
1166 case symetrie :
1167 if (flux_evaluateur.calculer_flux_faces_symetrie())
1168 {
1169 const Symetrie& cl =(const Symetrie&) (la_cl.valeur());
1170 for (face=ndeb; face<nfin; face++)
1171 {
1172 flux = flux_evaluateur.secmem_face(face, cl, ndeb);
1173 if ( (elem1=elem(face,0)) > -1)
1174 {
1175 resu[elem1]+=flux;
1176 flux_bords(face,0)+=flux;
1177 }
1178 if ( (elem2=elem(face,1)) > -1)
1179 {
1180 resu[elem2]-=flux;
1181 flux_bords(face,0)-=flux;
1182 }
1183 }
1184 }
1185 break;
1186 case sortie_libre :
1187 if (flux_evaluateur.calculer_flux_faces_sortie_libre())
1188 {
1189 const Neumann_sortie_libre& cl =(const Neumann_sortie_libre&) (la_cl.valeur());
1190 for (face=ndeb; face<nfin; face++)
1191 {
1192 flux = flux_evaluateur.secmem_face(face, cl, ndeb);
1193 if ( (elem1=elem(face,0)) > -1)
1194 {
1195 resu[elem1]+=flux;
1196 flux_bords(face,0)+=flux;
1197 }
1198 if ( (elem2=elem(face,1)) > -1)
1199 {
1200 resu[elem2]-=flux;
1201 flux_bords(face,0)-=flux;
1202 }
1203 }
1204 }
1205 break;
1206 case entree_fluide :
1207 if (flux_evaluateur.calculer_flux_faces_entree_fluide())
1208 {
1209 const Dirichlet_entree_fluide& cl =(const Dirichlet_entree_fluide&) (la_cl.valeur());
1210 for (face=ndeb; face<nfin; face++)
1211 {
1212 flux = flux_evaluateur.secmem_face(face, cl, ndeb);
1213 if ( (elem1=elem(face,0)) > -1)
1214 {
1215 resu[elem1]+=flux;
1216 flux_bords(face,0)+=flux;
1217 }
1218 if ( (elem2=elem(face,1)) > -1)
1219 {
1220 resu[elem2]-=flux;
1221 flux_bords(face,0)-=flux;
1222 }
1223 }
1224 }
1225 break;
1226 case paroi_fixe :
1227 if (flux_evaluateur.calculer_flux_faces_paroi_fixe())
1228 {
1229 const Dirichlet_paroi_fixe& cl =(const Dirichlet_paroi_fixe&) (la_cl.valeur());
1230 for (face=ndeb; face<nfin; face++)
1231 {
1232 flux = flux_evaluateur.secmem_face(face, cl, ndeb);
1233 if ( (elem1=elem(face,0)) > -1)
1234 {
1235 resu[elem1]+=flux;
1236 flux_bords(face,0)+=flux;
1237 }
1238 if ( (elem2=elem(face,1)) > -1)
1239 {
1240 resu[elem2]-=flux;
1241 flux_bords(face,0)-=flux;
1242 }
1243 }
1244 }
1245 break;
1246 case paroi_defilante :
1247 if (flux_evaluateur.calculer_flux_faces_paroi_defilante())
1248 {
1249 const Dirichlet_paroi_defilante& cl =(const Dirichlet_paroi_defilante&) (la_cl.valeur());
1250 for (face=ndeb; face<nfin; face++)
1251 {
1252 flux = flux_evaluateur.secmem_face(face, cl, ndeb);
1253 if ( (elem1=elem(face,0)) > -1)
1254 {
1255 resu[elem1]+=flux;
1256 flux_bords(face,0)+=flux;
1257 }
1258 if ( (elem2=elem(face,1)) > -1)
1259 {
1260 resu[elem2]-=flux;
1261 flux_bords(face,0)-=flux;
1262 }
1263 }
1264 }
1265 break;
1266 case paroi_adiabatique :
1267 if (flux_evaluateur.calculer_flux_faces_paroi_adiabatique())
1268 {
1269 const Neumann_paroi_adiabatique& cl =(const Neumann_paroi_adiabatique&) (la_cl.valeur());
1270 for (face=ndeb; face<nfin; face++)
1271 {
1272 flux = flux_evaluateur.secmem_face(face, cl, ndeb);
1273 if ( (elem1=elem(face,0)) > -1)
1274 {
1275 resu[elem1]+=flux;
1276 flux_bords(face,0)+=flux;
1277 }
1278 if ( (elem2=elem(face,1)) > -1)
1279 {
1280 resu[elem2]-=flux;
1281 flux_bords(face,0)-=flux;
1282 }
1283 }
1284 }
1285 break;
1286 case paroi :
1287 if (flux_evaluateur.calculer_flux_faces_paroi())
1288 {
1289 const Neumann_paroi& cl =(const Neumann_paroi&) (la_cl.valeur());
1290 for (face=ndeb; face<nfin; face++)
1291 {
1292 flux = flux_evaluateur.secmem_face(face, cl, ndeb);
1293 if ( (elem1=elem(face,0)) > -1)
1294 {
1295 resu[elem1]+=flux;
1296 flux_bords(face,0)+=flux;
1297 }
1298 if ( (elem2=elem(face,1)) > -1)
1299 {
1300 resu[elem2]-=flux;
1301 flux_bords(face,0)-=flux;
1302 }
1303 }
1304 }
1305 break;
1306 case echange_externe_impose :
1307 if (flux_evaluateur.calculer_flux_faces_echange_externe_impose())
1308 {
1309 const Echange_externe_impose& cl =(const Echange_externe_impose&) (la_cl.valeur());
1310
1311 int boundary_index=-1;
1312 if (le_domaine->front_VF(num_cl).le_nom() == frontiere_dis.le_nom())
1313 boundary_index=num_cl;
1314
1315 for (face=ndeb; face<nfin; face++)
1316 {
1317 int local_face=le_domaine->front_VF(boundary_index).num_local_face(face);
1318 flux = flux_evaluateur.secmem_face(boundary_index,face,local_face, cl, ndeb);
1319 if ( (elem1=elem(face,0)) > -1)
1320 {
1321 resu[elem1]+=flux;
1322 flux_bords(face,0)+=flux;
1323 }
1324 if ( (elem2=elem(face,1)) > -1)
1325 {
1326 resu[elem2]-=flux;
1327 flux_bords(face,0)-=flux;
1328 }
1329 }
1330 }
1331 break;
1332 case echange_global_impose :
1333 if (flux_evaluateur.calculer_flux_faces_echange_global_impose())
1334 {
1335 const Echange_global_impose& cl =(const Echange_global_impose&) (la_cl.valeur());
1336 for (face=ndeb; face<nfin; face++)
1337 {
1338 flux = flux_evaluateur.secmem_face(face, cl, ndeb);
1339 if ( (elem1=elem(face,0)) > -1)
1340 {
1341 resu[elem1]+=flux;
1342 flux_bords(face,0)+=flux;
1343 }
1344 if ( (elem2=elem(face,1)) > -1)
1345 {
1346 resu[elem2]-=flux;
1347 flux_bords(face,0)-=flux;
1348 }
1349 }
1350 }
1351 break;
1352 case periodique :
1353 if (flux_evaluateur.calculer_flux_faces_periodique())
1354 {
1355 const Periodique& cl =(const Periodique&) (la_cl.valeur());
1356 for (face=ndeb; face<nfin; face++)
1357 {
1358 if ( (elem1=elem(face,0)) > -1)
1359 resu[elem1]+=0.5*flux_evaluateur.secmem_face(face, cl, ndeb);
1360 if ( (elem2=elem(face,1)) > -1)
1361 resu[elem2]-=0.5*flux_evaluateur.secmem_face(face, cl, ndeb);
1362 }
1363 }
1364 break;
1365 case scalaire_impose_paroi :
1366 break;
1367 /*
1368 case nouvelle_Cl_PolyMAC_CDO :
1369 if (flux_evaluateur.calculer_flux_faces_echange_global_impose()){
1370 const Nouvelle_Cl_PolyMAC_CDO& cl =(const Nouvelle_Cl_PolyMAC_CDO&) (la_cl.valeur());
1371 for (face=ndeb; face<nfin; face++) {
1372 if ( (elem1=elem(face,0)) > -1)
1373 resu[elem1]+=flux_evaluateur.secmem_face(face, cl, ndeb);
1374 if ( (elem2=elem(face,1)) > -1)
1375 resu[elem2]-=flux_evaluateur.secmem_face(face, cl, ndeb);
1376 }
1377 }
1378 break;
1379 */
1380 default :
1381 Cerr << "On ne reconnait pas la condition limite : " << la_cl.valeur();
1382 exit();
1383 break;
1384 }
1385 }
1386}
1387
1388template <class _TYPE_> void Iterateur_PolyMAC_CDO_Elem<_TYPE_>::contribuer_au_second_membre_bords(DoubleTab& resu,int ncomp) const
1389{
1390 int elem1, elem2;
1391 int ndeb, nfin;
1392 int face,k;
1393 DoubleVect flux(ncomp);
1394 int num_cl=0;
1395 int nb_front_Cl=le_domaine->nb_front_Cl();
1396 DoubleTab& flux_bords=op_base->flux_bords();
1397 for (; num_cl<nb_front_Cl; num_cl++)
1398 {
1399 const Cond_lim& la_cl = la_zcl->les_conditions_limites(num_cl);
1400 const Front_VF& frontiere_dis = ref_cast(Front_VF,la_cl->frontiere_dis());
1401 ndeb = frontiere_dis.num_premiere_face();
1402 nfin = ndeb + frontiere_dis.nb_faces();
1403 switch(type_cl(la_cl))
1404 {
1405 case symetrie :
1406 if (flux_evaluateur.calculer_flux_faces_symetrie())
1407 {
1408 const Symetrie& cl =(const Symetrie&) (la_cl.valeur());
1409 for (face=ndeb; face<nfin; face++)
1410 {
1411 flux_evaluateur.secmem_face(face, cl, ndeb, flux);
1412 if ( (elem1=elem(face,0)) > -1)
1413 for (k=0; k<ncomp; k++)
1414 {
1415 resu(elem1,k) +=flux(k);
1416 flux_bords(face,k)+=flux(k);
1417 }
1418 if ( (elem2=elem(face,1)) > -1)
1419 for (k=0; k<ncomp; k++)
1420 {
1421 resu(elem2,k) -=flux(k);
1422 flux_bords(face,k)-=flux(k);
1423 }
1424 }
1425 }
1426 break;
1427 case sortie_libre :
1428 if (flux_evaluateur.calculer_flux_faces_sortie_libre())
1429 {
1430 const Neumann_sortie_libre& cl =(const Neumann_sortie_libre&) (la_cl.valeur());
1431 for (face=ndeb; face<nfin; face++)
1432 {
1433 flux_evaluateur.secmem_face(face, cl, ndeb, flux);
1434 if ( (elem1=elem(face,0)) > -1)
1435 for (k=0; k<ncomp; k++)
1436 {
1437 resu(elem1,k) +=flux(k);
1438 flux_bords(face,k)+=flux(k);
1439 }
1440 if ( (elem2=elem(face,1)) > -1)
1441 for (k=0; k<ncomp; k++)
1442 {
1443 resu(elem2,k) -=flux(k);
1444 flux_bords(face,k)-=flux(k);
1445 }
1446 }
1447 }
1448 break;
1449 case entree_fluide :
1450 if (flux_evaluateur.calculer_flux_faces_entree_fluide())
1451 {
1452 const Dirichlet_entree_fluide& cl =(const Dirichlet_entree_fluide&) (la_cl.valeur());
1453 for (face=ndeb; face<nfin; face++)
1454 {
1455 flux_evaluateur.secmem_face(face, cl, ndeb, flux);
1456 if ( (elem1=elem(face,0)) > -1)
1457 for (k=0; k<ncomp; k++)
1458 {
1459 resu(elem1,k) +=flux(k);
1460 flux_bords(face,k)+=flux(k);
1461 }
1462 if ( (elem2=elem(face,1)) > -1)
1463 for (k=0; k<ncomp; k++)
1464 {
1465 resu(elem2,k) -=flux(k);
1466 flux_bords(face,k)-=flux(k);
1467 }
1468 }
1469 }
1470 break;
1471 case paroi_fixe :
1472 if (flux_evaluateur.calculer_flux_faces_paroi_fixe())
1473 {
1474 const Dirichlet_paroi_fixe& cl =(const Dirichlet_paroi_fixe&) (la_cl.valeur());
1475 for (face=ndeb; face<nfin; face++)
1476 {
1477 flux_evaluateur.secmem_face(face, cl, ndeb, flux);
1478 if ( (elem1=elem(face,0)) > -1)
1479 for (k=0; k<ncomp; k++)
1480 {
1481 resu(elem1,k) +=flux(k);
1482 flux_bords(face,k)+=flux(k);
1483 }
1484 if ( (elem2=elem(face,1)) > -1)
1485 for (k=0; k<ncomp; k++)
1486 {
1487 resu(elem2,k) -=flux(k);
1488 flux_bords(face,k)-=flux(k);
1489 }
1490 }
1491 }
1492 break;
1493 case paroi_defilante :
1494 if (flux_evaluateur.calculer_flux_faces_paroi_defilante())
1495 {
1496 const Dirichlet_paroi_defilante& cl =(const Dirichlet_paroi_defilante&) (la_cl.valeur());
1497 for (face=ndeb; face<nfin; face++)
1498 {
1499 flux_evaluateur.secmem_face(face, cl, ndeb, flux);
1500 if ( (elem1=elem(face,0)) > -1)
1501 for (k=0; k<ncomp; k++)
1502 {
1503 resu(elem1,k) +=flux(k);
1504 flux_bords(face,k)+=flux(k);
1505 }
1506 if ( (elem2=elem(face,1)) > -1)
1507 for (k=0; k<ncomp; k++)
1508 {
1509 resu(elem2,k) -=flux(k);
1510 flux_bords(face,k)-=flux(k);
1511 }
1512 }
1513 }
1514 break;
1515 case paroi_adiabatique :
1516 if (flux_evaluateur.calculer_flux_faces_paroi_adiabatique())
1517 {
1518 const Neumann_paroi_adiabatique& cl =(const Neumann_paroi_adiabatique&) (la_cl.valeur());
1519 for (face=ndeb; face<nfin; face++)
1520 {
1521 flux_evaluateur.secmem_face(face, cl, ndeb, flux);
1522 if ( (elem1=elem(face,0)) > -1)
1523 for (k=0; k<ncomp; k++)
1524 {
1525 resu(elem1,k) +=flux(k);
1526 flux_bords(face,k)+=flux(k);
1527 }
1528 if ( (elem2=elem(face,1)) > -1)
1529 for (k=0; k<ncomp; k++)
1530 {
1531 resu(elem2,k) -=flux(k);
1532 flux_bords(face,k)-=flux(k);
1533 }
1534 }
1535 }
1536 break;
1537 case paroi :
1538 if (flux_evaluateur.calculer_flux_faces_paroi())
1539 {
1540 const Neumann_paroi& cl =(const Neumann_paroi&) (la_cl.valeur());
1541 for (face=ndeb; face<nfin; face++)
1542 {
1543 flux_evaluateur.secmem_face(face, cl, ndeb, flux);
1544 if ( (elem1=elem(face,0)) > -1)
1545 for (k=0; k<ncomp; k++)
1546 {
1547 resu(elem1,k) +=flux(k);
1548 flux_bords(face,k)+=flux(k);
1549 }
1550 if ( (elem2=elem(face,1)) > -1)
1551 for (k=0; k<ncomp; k++)
1552 {
1553 resu(elem2,k) -=flux(k);
1554 flux_bords(face,k)-=flux(k);
1555 }
1556 }
1557 }
1558 break;
1559 case echange_externe_impose :
1560 if (flux_evaluateur.calculer_flux_faces_echange_externe_impose())
1561 {
1562 const Echange_externe_impose& cl =(const Echange_externe_impose&) (la_cl.valeur());
1563
1564 int boundary_index=-1;
1565 if (le_domaine->front_VF(num_cl).le_nom() == frontiere_dis.le_nom())
1566 boundary_index=num_cl;
1567
1568 for (face=ndeb; face<nfin; face++)
1569 {
1570 int local_face=le_domaine->front_VF(boundary_index).num_local_face(face);
1571 flux_evaluateur.secmem_face(boundary_index,face,local_face, cl, ndeb, flux);
1572 if ( (elem1=elem(face,0)) > -1)
1573 for (k=0; k<ncomp; k++)
1574 {
1575 resu(elem1,k) +=flux(k);
1576 flux_bords(face,k)+=flux(k);
1577 }
1578 if ( (elem2=elem(face,1)) > -1)
1579 for (k=0; k<ncomp; k++)
1580 {
1581 resu(elem2,k) -=flux(k);
1582 flux_bords(face,k)-=flux(k);
1583 }
1584 }
1585 }
1586 break;
1587 case echange_global_impose :
1588 if (flux_evaluateur.calculer_flux_faces_echange_global_impose())
1589 {
1590 const Echange_global_impose& cl =(const Echange_global_impose&) (la_cl.valeur());
1591 for (face=ndeb; face<nfin; face++)
1592 {
1593 flux_evaluateur.secmem_face(face, cl, ndeb, flux);
1594 if ( (elem1=elem(face,0)) > -1)
1595 for (k=0; k<ncomp; k++)
1596 {
1597 resu(elem1,k) +=flux(k);
1598 flux_bords(face,k)+=flux(k);
1599 }
1600 if ( (elem2=elem(face,1)) > -1)
1601 for (k=0; k<ncomp; k++)
1602 {
1603 resu(elem2,k) -=flux(k);
1604 flux_bords(face,k)-=flux(k);
1605 }
1606 }
1607 }
1608 break;
1609 case periodique :
1610 if (flux_evaluateur.calculer_flux_faces_periodique())
1611 {
1612 const Periodique& cl =(const Periodique&) (la_cl.valeur());
1613 for (face=ndeb; face<nfin; face++)
1614 {
1615 flux_evaluateur.secmem_face(face, cl, ndeb, flux);
1616 if ( (elem1=elem(face,0)) > -1)
1617 for (k=0; k<ncomp; k++)
1618 {
1619 resu(elem1,k) +=0.5*flux(k);
1620 flux_bords(face,k)+=0.5*flux(k);
1621 }
1622 if ( (elem2=elem(face,1)) > -1)
1623 for (k=0; k<ncomp; k++)
1624 {
1625 resu(elem2,k) -=0.5*flux(k);
1626 flux_bords(face,k)-=0.5*flux(k);
1627 }
1628 }
1629 }
1630 break;
1631 /*
1632 case nouvelle_Cl_PolyMAC_CDO :
1633 if (flux_evaluateur.calculer_flux_faces_echange_global_impose()){
1634 const Nouvelle_Cl_PolyMAC_CDO& cl =(const Nouvelle_Cl_PolyMAC_CDO&) (la_cl.valeur());
1635 for (face=ndeb; face<nfin; face++) {
1636 flux_evaluateur.secmem_face(face, cl, ndeb, flux);
1637 if ( (elem1=elem(face,0)) > -1)
1638 for (k=0; k<ncomp; k++)
1639 resu(elem1,k) +=flux(k);
1640 if ( (elem2=elem(face,1)) > -1)
1641 for (k=0; k<ncomp; k++)
1642 resu(elem2,k) -=flux(k);
1643 }
1644 }
1645 break;
1646 */
1647 default :
1648 Cerr << "On ne reconnait pas la condition limite : " << la_cl.valeur();
1649 exit();
1650 break;
1651 }
1652 }
1653}
1654
1655template <class _TYPE_> void Iterateur_PolyMAC_CDO_Elem<_TYPE_>::contribuer_au_second_membre_interne(DoubleTab& resu) const
1656{
1657 const Domaine_PolyMAC_CDO& domaine_PolyMAC_CDO = le_domaine.valeur();
1658 double flux;
1659 int face;
1660 int ndeb=domaine_PolyMAC_CDO.premiere_face_int();
1661 int nfin=domaine_PolyMAC_CDO.nb_faces();
1662 for (face=ndeb; face<nfin; face++)
1663 {
1664 flux=flux_evaluateur.secmem_faces_interne(face);
1665 resu[elem(face,0)]+=flux;
1666 resu[elem(face,1)]-=flux;
1667 }
1668}
1669template <class _TYPE_> void Iterateur_PolyMAC_CDO_Elem<_TYPE_>::contribuer_au_second_membre_interne( DoubleTab& resu,int ncomp) const
1670{
1671 const Domaine_PolyMAC_CDO& domaine_PolyMAC_CDO = le_domaine.valeur();
1672 DoubleVect flux(ncomp);
1673 int face,k;
1674 int elem0,elem1;
1675 int ndeb=domaine_PolyMAC_CDO.premiere_face_int();
1676 int nfin=domaine_PolyMAC_CDO.nb_faces();
1677 for (face=ndeb; face<nfin; face++)
1678 {
1679 flux_evaluateur.secmem_faces_interne(face, flux);
1680 elem0 = elem(face,0);
1681 elem1 = elem(face,1);
1682 for (k=0; k<ncomp; k++)
1683 {
1684 resu(elem0,k)+=flux(k);
1685 resu(elem1,k)-=flux(k);
1686 }
1687 }
1688}
1689
1690template <class _TYPE_> void Iterateur_PolyMAC_CDO_Elem<_TYPE_>::ajouter_contribution(const DoubleTab& inco, Matrice_Morse& matrice) const
1691{
1692 ((_TYPE_&) flux_evaluateur).mettre_a_jour();
1693 assert(inco.nb_dim() < 3);
1694 assert(la_zcl);
1695 assert(le_domaine);
1696 int ncomp=1;
1697 if (inco.nb_dim() == 2)
1698 ncomp=inco.dimension(1);
1699 DoubleTab& flux_bords=op_base->flux_bords();
1700 flux_bords.resize(le_domaine->nb_faces_bord(),ncomp);
1701 flux_bords=0;
1702 if( ncomp == 1) /* cas scalaire */
1703 {
1704 ajouter_contribution_bords(inco, matrice) ;
1705 ajouter_contribution_interne(inco, matrice) ;
1706 }
1707 else /* cas vectoriel */
1708 {
1709 ajouter_contribution_bords(inco, matrice, ncomp) ;
1710 ajouter_contribution_interne(inco, matrice, ncomp) ;
1711 }
1712}
1713template <class _TYPE_> void Iterateur_PolyMAC_CDO_Elem<_TYPE_>::ajouter_contribution_bords(const DoubleTab& inco, Matrice_Morse& matrice ) const
1714{
1715 int elem1, elem2;
1716 double aii=0, ajj=0;
1717 int ndeb, nfin;
1718 int face;
1719 int num_cl=0;
1720 int nb_front_Cl=le_domaine->nb_front_Cl();
1721 for (; num_cl<nb_front_Cl; num_cl++)
1722 {
1723 /* pour chaque Condition Limite on regarde son type */
1724 const Cond_lim& la_cl = la_zcl->les_conditions_limites(num_cl);
1725 const Front_VF& frontiere_dis = ref_cast(Front_VF,la_cl->frontiere_dis());
1726 ndeb = frontiere_dis.num_premiere_face();
1727 nfin = ndeb + frontiere_dis.nb_faces();
1728 switch(type_cl(la_cl))
1729 {
1730 case symetrie :
1731 if (flux_evaluateur.calculer_flux_faces_symetrie())
1732 {
1733 const Symetrie& cl =(const Symetrie&) (la_cl.valeur());
1734 for (face=ndeb; face<nfin; face++)
1735 {
1736 flux_evaluateur.coeffs_face(face,ndeb, cl, aii, ajj);
1737 if ( (elem1=elem(face,0)) > -1)
1738 {
1739 matrice(elem1,elem1)+=aii;
1740 }
1741 if ( (elem2=elem(face,1)) > -1)
1742 {
1743 matrice(elem2,elem2)+=ajj;
1744 }
1745 }
1746 }
1747 break;
1748 case sortie_libre :
1749 if (flux_evaluateur.calculer_flux_faces_sortie_libre())
1750 {
1751 const Neumann_sortie_libre& cl =(const Neumann_sortie_libre&) (la_cl.valeur());
1752 for (face=ndeb; face<nfin; face++)
1753 {
1754 flux_evaluateur.coeffs_face(face,ndeb, cl, aii, ajj);
1755 if ( (elem1=elem(face,0)) > -1)
1756 {
1757 matrice(elem1,elem1)+=aii;
1758 }
1759 if ( (elem2=elem(face,1)) > -1)
1760 {
1761 matrice(elem2,elem2)+=ajj;
1762 }
1763 }
1764 }
1765 break;
1766 case entree_fluide :
1767 if (flux_evaluateur.calculer_flux_faces_entree_fluide())
1768 {
1769 const Dirichlet_entree_fluide& cl =(const Dirichlet_entree_fluide&) (la_cl.valeur());
1770 for (face=ndeb; face<nfin; face++)
1771 {
1772 flux_evaluateur.coeffs_face(face,ndeb, cl, aii, ajj);
1773 if ( (elem1=elem(face,0)) > -1)
1774 {
1775 matrice(elem1,elem1)+=aii;
1776 }
1777 if ( (elem2=elem(face,1)) > -1)
1778 {
1779 matrice(elem2,elem2)+=ajj;
1780 }
1781 }
1782 }
1783 break;
1784 case paroi_fixe :
1785 if (flux_evaluateur.calculer_flux_faces_paroi_fixe())
1786 {
1787 const Dirichlet_paroi_fixe& cl =(const Dirichlet_paroi_fixe&) (la_cl.valeur());
1788 for (face=ndeb; face<nfin; face++)
1789 {
1790 flux_evaluateur.coeffs_face(face,ndeb, cl, aii, ajj);
1791 if ( (elem1=elem(face,0)) > -1)
1792 {
1793 matrice(elem1,elem1)+=aii;
1794 }
1795 if ( (elem2=elem(face,1)) > -1)
1796 {
1797 matrice(elem2,elem2)+=ajj;
1798 }
1799 }
1800 }
1801 break;
1802 case paroi_defilante :
1803 if (flux_evaluateur.calculer_flux_faces_paroi_defilante())
1804 {
1805 const Dirichlet_paroi_defilante& cl =(const Dirichlet_paroi_defilante&) (la_cl.valeur());
1806 for (face=ndeb; face<nfin; face++)
1807 {
1808 flux_evaluateur.coeffs_face(face,ndeb, cl, aii, ajj);
1809 if ( (elem1=elem(face,0)) > -1)
1810 {
1811 matrice(elem1,elem1)+=aii;
1812 }
1813 if ( (elem2=elem(face,1)) > -1)
1814 {
1815 matrice(elem2,elem2)+=ajj;
1816 }
1817 }
1818 }
1819 break;
1820 case paroi_adiabatique :
1821 if (flux_evaluateur.calculer_flux_faces_paroi_adiabatique())
1822 {
1823 const Neumann_paroi_adiabatique& cl =(const Neumann_paroi_adiabatique&) (la_cl.valeur());
1824 for (face=ndeb; face<nfin; face++)
1825 {
1826 flux_evaluateur.coeffs_face(face,ndeb, cl, aii, ajj);
1827 if ( (elem1=elem(face,0)) > -1)
1828 {
1829 matrice(elem1,elem1)+=aii;
1830 }
1831 if ( (elem2=elem(face,1)) > -1)
1832 {
1833 matrice(elem2,elem2)+=ajj;
1834 }
1835 }
1836 }
1837 break;
1838 case paroi :
1839 if (flux_evaluateur.calculer_flux_faces_paroi())
1840 {
1841 const Neumann_paroi& cl =(const Neumann_paroi&) (la_cl.valeur());
1842 for (face=ndeb; face<nfin; face++)
1843 {
1844 flux_evaluateur.coeffs_face(face,ndeb, cl, aii, ajj);
1845 if ( (elem1=elem(face,0)) > -1)
1846 {
1847 matrice(elem1,elem1)+=aii;
1848 }
1849 if ( (elem2=elem(face,1)) > -1)
1850 {
1851 matrice(elem2,elem2)+=ajj;
1852 }
1853 }
1854 }
1855 break;
1856 case echange_externe_impose :
1857 if (flux_evaluateur.calculer_flux_faces_echange_externe_impose())
1858 {
1859 const Echange_externe_impose& cl =(const Echange_externe_impose&) (la_cl.valeur());
1860
1861 int boundary_index=-1;
1862 if (le_domaine->front_VF(num_cl).le_nom() == frontiere_dis.le_nom())
1863 boundary_index=num_cl;
1864
1865 for (face=ndeb; face<nfin; face++)
1866 {
1867 int local_face=le_domaine->front_VF(boundary_index).num_local_face(face);
1868 flux_evaluateur.coeffs_face(boundary_index,face,local_face,ndeb, cl, aii, ajj);
1869 if ( (elem1=elem(face,0)) > -1)
1870 {
1871 matrice(elem1,elem1)+=aii;
1872 }
1873 if ( (elem2=elem(face,1)) > -1)
1874 {
1875 matrice(elem2,elem2)+=ajj;
1876 }
1877 }
1878 }
1879 break;
1880 case echange_global_impose :
1881 if (flux_evaluateur.calculer_flux_faces_echange_global_impose())
1882 {
1883 const Echange_global_impose& cl =(const Echange_global_impose&) (la_cl.valeur());
1884 for (face=ndeb; face<nfin; face++)
1885 {
1886 flux_evaluateur.coeffs_face(face,ndeb, cl, aii, ajj);
1887 if ( (elem1=elem(face,0)) > -1)
1888 {
1889 matrice(elem1,elem1)+=aii;
1890 }
1891 if ( (elem2=elem(face,1)) > -1)
1892 {
1893 matrice(elem2,elem2)+=ajj;
1894 }
1895 }
1896 }
1897 break;
1898 case periodique :
1899 if (flux_evaluateur.calculer_flux_faces_periodique())
1900 {
1901 const Periodique& cl =(const Periodique&) (la_cl.valeur());
1902 for (face=ndeb; face<nfin; face++)
1903 {
1904 /* GF Qui a mis l'appel au coeffs faces_internes ??? . Cest faux ...*/
1905 /* flux_evaluateur.coeffs_faces_interne(face, aii, ajj); */
1906 flux_evaluateur.coeffs_face(face,ndeb, cl, aii, ajj);
1907 elem1 = elem(face,0);
1908 elem2 = elem(face,1);
1909 matrice(elem1,elem1)+=0.5*aii;
1910 matrice(elem1,elem2)-=0.5*ajj;
1911 matrice(elem2,elem2)+=0.5*ajj;
1912 matrice(elem2,elem1)-=0.5*aii;
1913 }
1914 }
1915 break;
1916 case scalaire_impose_paroi :
1917 break;
1918 default :
1919 Cerr << "On ne reconnait pas la condition limite : " << la_cl.valeur();
1920 exit();
1921 break;
1922 }
1923 }
1924}
1925
1926template <class _TYPE_> void Iterateur_PolyMAC_CDO_Elem<_TYPE_>::ajouter_contribution_bords(const DoubleTab& inco, Matrice_Morse& matrice ,int ncomp) const
1927{
1928 int elem1, elem2;
1929 DoubleVect aii(ncomp), ajj(ncomp);
1930 int ndeb, nfin;
1931 int face,i;
1932 int num_cl=0;
1933 int nb_front_Cl=le_domaine->nb_front_Cl();
1934 for (; num_cl<nb_front_Cl; num_cl++)
1935 {
1936 const Cond_lim& la_cl = la_zcl->les_conditions_limites(num_cl);
1937 const Front_VF& frontiere_dis = ref_cast(Front_VF,la_cl->frontiere_dis());
1938 ndeb = frontiere_dis.num_premiere_face();
1939 nfin = ndeb + frontiere_dis.nb_faces();
1940 switch(type_cl(la_cl))
1941 {
1942 case symetrie :
1943 if (flux_evaluateur.calculer_flux_faces_symetrie())
1944 {
1945 const Symetrie& cl =(const Symetrie&) (la_cl.valeur());
1946 for (face=ndeb; face<nfin; face++)
1947 {
1948 elem1 = elem(face,0);
1949 elem2 = elem(face,1);
1950 flux_evaluateur.coeffs_face(face,ndeb, cl, aii, ajj);
1951 if ( (elem1=elem(face,0)) > -1)
1952 {
1953 for (i=0; i<ncomp; i++)
1954 {
1955 matrice(elem1*ncomp+i,elem1*ncomp+i)+=aii(i);
1956 }
1957 }
1958 if ( (elem2=elem(face,1)) > -1)
1959 {
1960 for (i=0; i<ncomp; i++)
1961 {
1962 matrice(elem2*ncomp+i,elem2*ncomp+i)+=ajj(i);
1963 }
1964 }
1965 }
1966 }
1967 break;
1968 case sortie_libre :
1969 if (flux_evaluateur.calculer_flux_faces_sortie_libre())
1970 {
1971 const Neumann_sortie_libre& cl =(const Neumann_sortie_libre&) (la_cl.valeur());
1972 for (face=ndeb; face<nfin; face++)
1973 {
1974 flux_evaluateur.coeffs_face(face,ndeb, cl, aii, ajj);
1975 elem1 = elem(face,0);
1976 elem2 = elem(face,1);
1977 flux_evaluateur.coeffs_face(face,ndeb, cl, aii, ajj);
1978 if ( (elem1=elem(face,0)) > -1)
1979 {
1980 for (i=0; i<ncomp; i++)
1981 {
1982 matrice(elem1*ncomp+i,elem1*ncomp+i)+=aii(i);
1983 }
1984 }
1985 if ( (elem2=elem(face,1)) > -1)
1986 {
1987 for (i=0; i<ncomp; i++)
1988 {
1989 matrice(elem2*ncomp+i,elem2*ncomp+i)+=ajj(i);
1990 }
1991 }
1992 }
1993 }
1994 break;
1995 case entree_fluide :
1996 if (flux_evaluateur.calculer_flux_faces_entree_fluide())
1997 {
1998 const Dirichlet_entree_fluide& cl =(const Dirichlet_entree_fluide&) (la_cl.valeur());
1999 for (face=ndeb; face<nfin; face++)
2000 {
2001 flux_evaluateur.coeffs_face(face,ndeb, cl, aii, ajj);
2002 elem1 = elem(face,0);
2003 elem2 = elem(face,1);
2004 flux_evaluateur.coeffs_face(face,ndeb, cl, aii, ajj);
2005 if ( (elem1=elem(face,0)) > -1)
2006 {
2007 for (i=0; i<ncomp; i++)
2008 {
2009 matrice(elem1*ncomp+i,elem1*ncomp+i)+=aii(i);
2010 }
2011 }
2012 if ( (elem2=elem(face,1)) > -1)
2013 {
2014 for (i=0; i<ncomp; i++)
2015 {
2016 matrice(elem2*ncomp+i,elem2*ncomp+i)+=ajj(i);
2017 }
2018 }
2019 }
2020 }
2021 break;
2022 case paroi_fixe :
2023 if (flux_evaluateur.calculer_flux_faces_paroi_fixe())
2024 {
2025 const Dirichlet_paroi_fixe& cl =(const Dirichlet_paroi_fixe&) (la_cl.valeur());
2026 for (face=ndeb; face<nfin; face++)
2027 {
2028 flux_evaluateur.coeffs_face(face,ndeb, cl, aii, ajj);
2029 elem1 = elem(face,0);
2030 elem2 = elem(face,1);
2031 flux_evaluateur.coeffs_face(face,ndeb, cl, aii, ajj);
2032 if ( (elem1=elem(face,0)) > -1)
2033 {
2034 for (i=0; i<ncomp; i++)
2035 {
2036 matrice(elem1*ncomp+i,elem1*ncomp+i)+=aii(i);
2037 }
2038 }
2039 if ( (elem2=elem(face,1)) > -1)
2040 {
2041 for (i=0; i<ncomp; i++)
2042 {
2043 matrice(elem2*ncomp+i,elem2*ncomp+i)+=ajj(i);
2044 }
2045 }
2046 }
2047 }
2048 break;
2049 case paroi_defilante :
2050 if (flux_evaluateur.calculer_flux_faces_paroi_defilante())
2051 {
2052 const Dirichlet_paroi_defilante& cl =(const Dirichlet_paroi_defilante&) (la_cl.valeur());
2053 for (face=ndeb; face<nfin; face++)
2054 {
2055 flux_evaluateur.coeffs_face(face,ndeb, cl, aii, ajj);
2056 elem1 = elem(face,0);
2057 elem2 = elem(face,1);
2058 flux_evaluateur.coeffs_face(face,ndeb, cl, aii, ajj);
2059 if ( (elem1=elem(face,0)) > -1)
2060 {
2061 for (i=0; i<ncomp; i++)
2062 {
2063 matrice(elem1*ncomp+i,elem1*ncomp+i)+=aii(i);
2064 }
2065 }
2066 if ( (elem2=elem(face,1)) > -1)
2067 {
2068 for (i=0; i<ncomp; i++)
2069 {
2070 matrice(elem2*ncomp+i,elem2*ncomp+i)+=ajj(i);
2071 }
2072 }
2073 }
2074 }
2075 break;
2076 case paroi_adiabatique :
2077 if (flux_evaluateur.calculer_flux_faces_paroi_adiabatique())
2078 {
2079 const Neumann_paroi_adiabatique& cl =(const Neumann_paroi_adiabatique&) (la_cl.valeur());
2080 for (face=ndeb; face<nfin; face++)
2081 {
2082 flux_evaluateur.coeffs_face(face,ndeb, cl, aii, ajj);
2083 elem1 = elem(face,0);
2084 elem2 = elem(face,1);
2085 flux_evaluateur.coeffs_face(face,ndeb, cl, aii, ajj);
2086 if ( (elem1=elem(face,0)) > -1)
2087 {
2088 for (i=0; i<ncomp; i++)
2089 {
2090 matrice(elem1*ncomp+i,elem1*ncomp+i)+=aii(i);
2091 }
2092 }
2093 if ( (elem2=elem(face,1)) > -1)
2094 {
2095 for (i=0; i<ncomp; i++)
2096 {
2097 matrice(elem2*ncomp+i,elem2*ncomp+i)+=ajj(i);
2098 }
2099 }
2100 }
2101 }
2102 break;
2103 case paroi :
2104 if (flux_evaluateur.calculer_flux_faces_paroi())
2105 {
2106 const Neumann_paroi& cl =(const Neumann_paroi&) (la_cl.valeur());
2107 for (face=ndeb; face<nfin; face++)
2108 {
2109 flux_evaluateur.coeffs_face(face,ndeb, cl, aii, ajj);
2110 elem1 = elem(face,0);
2111 elem2 = elem(face,1);
2112 flux_evaluateur.coeffs_face(face,ndeb, cl, aii, ajj);
2113 if ( (elem1=elem(face,0)) > -1)
2114 {
2115 for (i=0; i<ncomp; i++)
2116 {
2117 matrice(elem1*ncomp+i,elem1*ncomp+i)+=aii(i);
2118 }
2119 }
2120 if ( (elem2=elem(face,1)) > -1)
2121 {
2122 for (i=0; i<ncomp; i++)
2123 {
2124 matrice(elem2*ncomp+i,elem2*ncomp+i)+=ajj(i);
2125 }
2126 }
2127 }
2128 }
2129 break;
2130 case echange_externe_impose :
2131 if (flux_evaluateur.calculer_flux_faces_echange_externe_impose())
2132 {
2133 const Echange_externe_impose& cl =(const Echange_externe_impose&) (la_cl.valeur());
2134
2135 int boundary_index=-1;
2136 if (le_domaine->front_VF(num_cl).le_nom() == frontiere_dis.le_nom())
2137 boundary_index=num_cl;
2138
2139 for (face=ndeb; face<nfin; face++)
2140 {
2141 int local_face=le_domaine->front_VF(boundary_index).num_local_face(face);
2142 flux_evaluateur.coeffs_face(boundary_index,face,local_face,ndeb, cl, aii, ajj);
2143 elem1 = elem(face,0);
2144 elem2 = elem(face,1);
2145 if ( (elem1=elem(face,0)) > -1)
2146 {
2147 for (i=0; i<ncomp; i++)
2148 {
2149 matrice(elem1*ncomp+i,elem1*ncomp+i)+=aii(i);
2150 }
2151 }
2152 if ( (elem2=elem(face,1)) > -1)
2153 {
2154 for (i=0; i<ncomp; i++)
2155 {
2156 matrice(elem2*ncomp+i,elem2*ncomp+i)+=aii(i);
2157 }
2158 }
2159 }
2160 }
2161 break;
2162 case echange_global_impose :
2163 if (flux_evaluateur.calculer_flux_faces_echange_global_impose())
2164 {
2165 const Echange_global_impose& cl =(const Echange_global_impose&) (la_cl.valeur());
2166 for (face=ndeb; face<nfin; face++)
2167 {
2168 flux_evaluateur.coeffs_face(face,ndeb, cl, aii, ajj);
2169 elem1 = elem(face,0);
2170 elem2 = elem(face,1);
2171 flux_evaluateur.coeffs_face(face,ndeb, cl, aii, ajj);
2172 if ( (elem1=elem(face,0)) > -1)
2173 {
2174 for (i=0; i<ncomp; i++)
2175 {
2176 matrice(elem1*ncomp+i,elem1*ncomp+i)+=aii(i);
2177 }
2178 }
2179 if ( (elem2=elem(face,1)) > -1)
2180 {
2181 for (i=0; i<ncomp; i++)
2182 {
2183 matrice(elem2*ncomp+i,elem2*ncomp+i)+=ajj(i);
2184 }
2185 }
2186 }
2187 }
2188 break;
2189 case periodique :
2190 if (flux_evaluateur.calculer_flux_faces_periodique())
2191 {
2192 const Periodique& cl =(const Periodique&) (la_cl.valeur());
2193 for (face=ndeb; face<nfin; face++)
2194 {
2195 flux_evaluateur.coeffs_face(face,ndeb, cl, aii, ajj);
2196 elem1 = elem(face,0);
2197 elem2 = elem(face,1);
2198 flux_evaluateur.coeffs_face(face,ndeb, cl, aii, ajj);
2199 {
2200 for (i=0; i<ncomp; i++)
2201 {
2202 int n1=elem1*ncomp+i;
2203 int n2=elem2*ncomp+i;
2204 matrice(n1,n1)+=0.5*aii(i);
2205 matrice(n1,n2)-=0.5*ajj(i);
2206 matrice(n2,n2)+=0.5*ajj(i);
2207 matrice(n2,n1)-=0.5*aii(i);
2208 }
2209 }
2210 }
2211 }
2212 break;
2213 case scalaire_impose_paroi :
2214 // on n'a rien a convecter
2215 break;
2216 default :
2217 Cerr << "On ne reconnait pas la condition limite : " << la_cl.valeur();
2218 exit();
2219 break;
2220 }
2221 }
2222}
2223
2224template <class _TYPE_> void Iterateur_PolyMAC_CDO_Elem<_TYPE_>::ajouter_contribution_interne(const DoubleTab& inco, Matrice_Morse& matrice ) const
2225{
2226 const Domaine_PolyMAC_CDO& domaine_PolyMAC_CDO = le_domaine.valeur();
2227 int face;
2228 double aii=0, ajj=0;
2229 int elem1,elem2;
2230 int ndeb=domaine_PolyMAC_CDO.premiere_face_int();
2231 int nfin=domaine_PolyMAC_CDO.nb_faces();
2232 for (face=ndeb; face<nfin; face++)
2233 {
2234 elem1 = elem(face,0);
2235 elem2 = elem(face,1);
2236 flux_evaluateur.coeffs_faces_interne(face, aii, ajj);
2237 matrice(elem1,elem1)+=aii;
2238 matrice(elem1,elem2)-=ajj;
2239 matrice(elem2,elem2)+=ajj;
2240 matrice(elem2,elem1)-=aii;
2241 }
2242}
2243template <class _TYPE_> void Iterateur_PolyMAC_CDO_Elem<_TYPE_>::ajouter_contribution_interne(const DoubleTab& inco, Matrice_Morse& matrice ,int ncomp) const
2244{
2245 const Domaine_PolyMAC_CDO& domaine_PolyMAC_CDO = le_domaine.valeur();
2246 int face,i;
2247 DoubleVect aii(ncomp), ajj(ncomp);
2248 int elem1,elem2;
2249 int ndeb=domaine_PolyMAC_CDO.premiere_face_int();
2250 int nfin=domaine_PolyMAC_CDO.nb_faces();
2251 for (face=ndeb; face<nfin; face++)
2252 {
2253 elem1 = elem(face,0);
2254 elem2 = elem(face,1);
2255 flux_evaluateur.coeffs_faces_interne(face, aii, ajj);
2256 for (i=0; i<ncomp; i++)
2257 {
2258 int i0=elem1*ncomp+i,j0=elem2*ncomp+i;
2259 matrice(i0,i0)+=aii(i) ;
2260 matrice(i0,j0)-=ajj(i) ;
2261 matrice(j0,j0)+=ajj(i) ;
2262 matrice(j0,i0)-=aii(i) ;
2263 }
2264 }
2265}
2266
2267template <class _TYPE_> void Iterateur_PolyMAC_CDO_Elem<_TYPE_>::ajouter_contribution_vitesse(const DoubleTab& inco, Matrice_Morse& matrice) const
2268{
2269 ((_TYPE_&) flux_evaluateur).mettre_a_jour();
2270 assert(inco.nb_dim() < 3);
2271 assert(la_zcl);
2272 assert(le_domaine);
2273 int ncomp=1;
2274 if (inco.nb_dim() == 2)
2275 ncomp=inco.dimension(1);
2276 // DoubleTab& flux_bords=op_base->flux_bords();
2277 // flux_bords.resize(le_domaine->nb_faces_bord(),ncomp);
2278 // flux_bords=0;
2279 if( ncomp == 1) /* cas scalaire */
2280 {
2281 ajouter_contribution_bords_vitesse(inco, matrice);
2282 ajouter_contribution_interne_vitesse(inco, matrice);
2283 }
2284 else /* cas vectoriel */
2285 {
2286 abort();
2287 // ajouter_contribution_bords(inco, matrice, ncomp) ;
2288 // ajouter_contribution_interne(inco, matrice, ncomp) ;
2289 }
2290}
2291
2292template <class _TYPE_> void Iterateur_PolyMAC_CDO_Elem<_TYPE_>::ajouter_contribution_interne_vitesse(const DoubleTab& inco, Matrice_Morse& matrice) const
2293{
2294 const Domaine_PolyMAC_CDO& domaine_PolyMAC_CDO = le_domaine.valeur();
2295 double aef = 0;
2296 const int ndeb = domaine_PolyMAC_CDO.premiere_face_int();
2297 const int nfin = domaine_PolyMAC_CDO.nb_faces();
2298 for (int f = ndeb; f < nfin; f++)
2299 {
2300 const int e1 = elem(f, 0);
2301 const int e2 = elem(f, 1);
2302 aef = flux_evaluateur.coeffs_faces_interne_bloc_vitesse(inco, f);
2303 matrice(e1, f) += aef;
2304 matrice(e2, f) -= aef;
2305 }
2306}
2307
2308template <class _TYPE_> void Iterateur_PolyMAC_CDO_Elem<_TYPE_>::ajouter_contribution_bords_vitesse(const DoubleTab& inco, Matrice_Morse& matrice ) const
2309{
2310 int e1, e2;
2311 double aef=0;
2312 int ndeb, nfin;
2313 int nb_front_Cl=le_domaine->nb_front_Cl();
2314 for (int num_cl = 0; num_cl<nb_front_Cl; num_cl++)
2315 {
2316 /* pour chaque Condition Limite on regarde son type */
2317 const Cond_lim& la_cl = la_zcl->les_conditions_limites(num_cl);
2318 const Front_VF& frontiere_dis = ref_cast(Front_VF,la_cl->frontiere_dis());
2319 ndeb = frontiere_dis.num_premiere_face();
2320 nfin = ndeb + frontiere_dis.nb_faces();
2321 switch(type_cl(la_cl))
2322 {
2323 case symetrie :
2324 if (flux_evaluateur.calculer_flux_faces_symetrie())
2325 {
2326 const Symetrie& cl =(const Symetrie&) (la_cl.valeur());
2327 for (int f = ndeb; f < nfin; f++)
2328 {
2329 aef = flux_evaluateur.coeffs_face_bloc_vitesse(inco, f, cl, ndeb);
2330 if ( (e1 = elem(f, 0)) > -1) matrice(e1, f) += aef;
2331 if ( (e2 = elem(f, 1)) > -1) matrice(e2, f) -= aef;
2332 }
2333 }
2334 break;
2335 case sortie_libre :
2336 if (flux_evaluateur.calculer_flux_faces_sortie_libre())
2337 {
2338 const Neumann_sortie_libre& cl =(const Neumann_sortie_libre&) (la_cl.valeur());
2339 for (int f = ndeb; f < nfin; f++)
2340 {
2341 aef = flux_evaluateur.coeffs_face_bloc_vitesse(inco, f, cl, ndeb);
2342 if ( (e1 = elem(f, 0)) > -1) matrice(e1, f) += aef;
2343 if ( (e2 = elem(f, 1)) > -1) matrice(e2, f) -= aef;
2344 }
2345 }
2346 break;
2347 case entree_fluide :
2348 if (flux_evaluateur.calculer_flux_faces_entree_fluide())
2349 {
2350 const Dirichlet_entree_fluide& cl =(const Dirichlet_entree_fluide&) (la_cl.valeur());
2351 for (int f = ndeb; f < nfin; f++)
2352 {
2353 aef = flux_evaluateur.coeffs_face_bloc_vitesse(inco, f, cl, ndeb);
2354 if ( (e1 = elem(f, 0)) > -1) matrice(e1, f) += aef;
2355 if ( (e2 = elem(f, 1)) > -1) matrice(e2, f) -= aef;
2356 }
2357 }
2358 break;
2359 case paroi_fixe :
2360 if (flux_evaluateur.calculer_flux_faces_paroi_fixe())
2361 {
2362 const Dirichlet_paroi_fixe& cl =(const Dirichlet_paroi_fixe&) (la_cl.valeur());
2363 for (int f = ndeb; f < nfin; f++)
2364 {
2365 aef = flux_evaluateur.coeffs_face_bloc_vitesse(inco, f, cl, ndeb);
2366 if ( (e1 = elem(f, 0)) > -1) matrice(e1, f) += aef;
2367 if ( (e2 = elem(f, 1)) > -1) matrice(e2, f) -= aef;
2368 }
2369 }
2370 break;
2371 case paroi_defilante :
2372 if (flux_evaluateur.calculer_flux_faces_paroi_defilante())
2373 {
2374 const Dirichlet_paroi_defilante& cl =(const Dirichlet_paroi_defilante&) (la_cl.valeur());
2375 for (int f = ndeb; f < nfin; f++)
2376 {
2377 aef = flux_evaluateur.coeffs_face_bloc_vitesse(inco, f, cl, ndeb);
2378 if ( (e1 = elem(f, 0)) > -1) matrice(e1, f) += aef;
2379 if ( (e2 = elem(f, 1)) > -1) matrice(e2, f) -= aef;
2380 }
2381 }
2382 break;
2383 case paroi_adiabatique :
2384 if (flux_evaluateur.calculer_flux_faces_paroi_adiabatique())
2385 {
2386 const Neumann_paroi_adiabatique& cl =(const Neumann_paroi_adiabatique&) (la_cl.valeur());
2387 for (int f = ndeb; f < nfin; f++)
2388 {
2389 aef = flux_evaluateur.coeffs_face_bloc_vitesse(inco, f, cl, ndeb);
2390 if ( (e1 = elem(f, 0)) > -1) matrice(e1, f) += aef;
2391 if ( (e2 = elem(f, 1)) > -1) matrice(e2, f) -= aef;
2392 }
2393 }
2394 break;
2395 case paroi :
2396 if (flux_evaluateur.calculer_flux_faces_paroi())
2397 {
2398 const Neumann_paroi& cl =(const Neumann_paroi&) (la_cl.valeur());
2399 for (int f = ndeb; f < nfin; f++)
2400 {
2401 aef = flux_evaluateur.coeffs_face_bloc_vitesse(inco, f, cl, ndeb);
2402 if ( (e1 = elem(f, 0)) > -1) matrice(e1, f) += aef;
2403 if ( (e2 = elem(f, 1)) > -1) matrice(e2, f) -= aef;
2404 }
2405 }
2406 break;
2407 case echange_externe_impose :
2408 if (flux_evaluateur.calculer_flux_faces_echange_externe_impose())
2409 {
2410 const Echange_externe_impose& cl =(const Echange_externe_impose&) (la_cl.valeur());
2411 int boundary_index=-1;
2412 if (le_domaine->front_VF(num_cl).le_nom() == frontiere_dis.le_nom())
2413 boundary_index=num_cl;
2414
2415 for (int f = ndeb; f < nfin; f++)
2416 {
2417 int local_face=le_domaine->front_VF(boundary_index).num_local_face(f);
2418 aef = flux_evaluateur.coeffs_face_bloc_vitesse(inco, boundary_index, f, local_face, cl, ndeb);
2419 if ( (e1 = elem(f, 0)) > -1) matrice(e1, f) += aef;
2420 if ( (e2 = elem(f, 1)) > -1) matrice(e1, f) -= aef;
2421 }
2422 }
2423 break;
2424 case echange_global_impose :
2425 if (flux_evaluateur.calculer_flux_faces_echange_global_impose())
2426 {
2427 const Echange_global_impose& cl =(const Echange_global_impose&) (la_cl.valeur());
2428 for (int f = ndeb; f < nfin; f++)
2429 {
2430 aef = flux_evaluateur.coeffs_face_bloc_vitesse(inco, f, cl, ndeb);
2431 if ( (e1 = elem(f, 0)) > -1) matrice(e1, f) += aef;
2432 if ( (e2 = elem(f, 1)) > -1) matrice(e2, f) -= aef;
2433 }
2434 }
2435 break;
2436 case periodique :
2437 if (flux_evaluateur.calculer_flux_faces_periodique())
2438 {
2439 const Periodique& cl =(const Periodique&) (la_cl.valeur());
2440 for (int f = ndeb; f < nfin; f++)
2441 {
2442 aef = flux_evaluateur.coeffs_face_bloc_vitesse(inco, f, cl, ndeb);
2443 if ( (e1 = elem(f, 0)) > -1)
2444 {
2445 if ( f < (ndeb+frontiere_dis.nb_faces()/2) )
2446 matrice(e1, f) += aef;
2447 }
2448 if ( (e2 = elem(f, 1)) > -1)
2449 {
2450 if ( (ndeb+frontiere_dis.nb_faces()/2) <= f )
2451 matrice(e2, f) -= aef;
2452 }
2453 }
2454 }
2455 break;
2456 default :
2457 Cerr << "On ne reconnait pas la condition limite : " << la_cl.valeur();
2458 exit();
2459 break;
2460 }
2461
2462 }
2463}
2464
2465#endif /* Iterateur_PolyMAC_CDO_Elem_included */
Champ_Don_base
classe Champ_Don_base classe de base des Champs donnes (non calcules)
Definition Champ_Don_base.h:32
Champ_Don_base::valeurs
DoubleTab & valeurs() override
Surcharge Champ_base::valeurs() Renvoie le tableau des valeurs.
Definition Champ_Don_base.h:49
Champ_Proto::valeurs
virtual DoubleTab & valeurs()=0
Champ_Uniforme
classe Champ_Uniforme Represente un champ constant dans l'espace et dans le temps.
Definition Champ_Uniforme.h:26
Champ_base
classe Champ_base Cette classe est la base de la hierarchie des champs.
Definition Champ_base.h:43
Cond_lim
classe Cond_lim Classe generique servant a representer n'importe quelle classe
Definition Cond_lim.h:31
Dirichlet_entree_fluide
classe Dirichlet_entree_fluide Cette classe represente une condition aux limite imposant une grandeur
Definition Dirichlet_entree_fluide_leaves.h:33
Dirichlet_paroi_defilante
classe Dirichlet_paroi_defilante Impose la vitesse de paroi dnas une equation de type Navier_Stokes.
Definition Dirichlet_paroi_defilante.h:26
Dirichlet_paroi_fixe
classe Dirichlet_paroi_fixe Represente une paroi immobile dans une equation de type Navier_Stokes.
Definition Dirichlet_paroi_fixe.h:26
Domaine_PolyMAC_CDO
Definition Domaine_PolyMAC_CDO.h:31
Domaine_VF::nb_faces
int nb_faces() const
renvoie le nombre global de faces.
Definition Domaine_VF.h:471
Domaine_VF::premiere_face_int
int premiere_face_int() const
une face est interne ssi elle separe deux elements.
Definition Domaine_VF.h:463
Domaine_VF::face_voisins
int face_voisins(int num_face, int i) const
renvoie l'element voisin de numface dans la direction i.
Definition Domaine_VF.h:418
Domaine_VF::nb_faces_bord
int nb_faces_bord() const
renvoie le nombre de faces sur lesquelles sont appliquees les conditions limites :
Definition Domaine_VF.h:513
Echange_externe_impose
Classe Echange_externe_impose: Cette classe represente le cas particulier de la classe.
Definition Echange_externe_impose.h:46
Echange_global_impose
Classe Echange_global_impose Cette classe represente le cas particulier de la classe.
Definition Echange_global_impose.h:39
EcrFicPartage::syncfile
Sortie & syncfile() override
Provoque l'ecriture sur disque des donnees accumulees sur les differents processeurs depuis le dernie...
Definition EcrFicPartage.cpp:173
Evaluateur_PolyMAC_CDO
class Evaluateur_PolyMAC_CDO
Definition Evaluateur_PolyMAC_CDO.h:36
Field_base::nb_comp
virtual int nb_comp() const
Definition Field_base.h:56
Front_VF
class Front_VF
Definition Front_VF.h:36
Front_VF::nb_faces
int nb_faces() const
Definition Front_VF.h:53
Front_VF::num_premiere_face
int num_premiere_face() const
Definition Front_VF.h:63
Frontiere_dis_base
classe Frontiere_dis_base Classe representant une frontiere discretisee.
Definition Frontiere_dis_base.h:34
Frontiere_dis_base::le_nom
const Nom & le_nom() const override
Renvoie le nom de la frontiere geometrique.
Definition Frontiere_dis_base.cpp:81
Iterateur_PolyMAC_CDO_Elem
Definition Iterateur_PolyMAC_CDO_Elem.h:33
Iterateur_PolyMAC_CDO_Elem::ajouter_bords
DoubleTab & ajouter_bords(const DoubleTab &, DoubleTab &) const
Definition Iterateur_PolyMAC_CDO_Elem.h:141
Iterateur_PolyMAC_CDO_Elem::Flux
SFichier Flux
Definition Iterateur_PolyMAC_CDO_Elem.h:76
Iterateur_PolyMAC_CDO_Elem::impr
int impr(Sortie &os) const override
Definition Iterateur_PolyMAC_CDO_Elem.h:1029
Iterateur_PolyMAC_CDO_Elem::ajouter_interne
DoubleTab & ajouter_interne(const DoubleTab &, DoubleTab &) const
Definition Iterateur_PolyMAC_CDO_Elem.h:937
Iterateur_PolyMAC_CDO_Elem::Iterateur_PolyMAC_CDO_Elem
Iterateur_PolyMAC_CDO_Elem()
Definition Iterateur_PolyMAC_CDO_Elem.h:44
Iterateur_PolyMAC_CDO_Elem::contribuer_au_second_membre
void contribuer_au_second_membre(DoubleTab &) const override
Definition Iterateur_PolyMAC_CDO_Elem.h:1127
Iterateur_PolyMAC_CDO_Elem::ajouter
DoubleTab & ajouter(const DoubleTab &, DoubleTab &) const override
Definition Iterateur_PolyMAC_CDO_Elem.h:107
Iterateur_PolyMAC_CDO_Elem::evaluateur
Evaluateur_PolyMAC_CDO & evaluateur() override
Definition Iterateur_PolyMAC_CDO_Elem.h:84
Iterateur_PolyMAC_CDO_Elem::ajouter_contribution
void ajouter_contribution(const DoubleTab &, Matrice_Morse &) const override
Definition Iterateur_PolyMAC_CDO_Elem.h:1690
Iterateur_PolyMAC_CDO_Elem::milieu
const Milieu_base & milieu() const
Definition Iterateur_PolyMAC_CDO_Elem.h:101
Iterateur_PolyMAC_CDO_Elem::calculer_flux_bord
void calculer_flux_bord(const DoubleTab &) const override
Definition Iterateur_PolyMAC_CDO_Elem.h:402
Iterateur_PolyMAC_CDO_Elem::ajouter_contribution_interne
void ajouter_contribution_interne(const DoubleTab &, Matrice_Morse &) const
Definition Iterateur_PolyMAC_CDO_Elem.h:2224
Iterateur_PolyMAC_CDO_Elem::flux_evaluateur
_TYPE_ flux_evaluateur
Definition Iterateur_PolyMAC_CDO_Elem.h:59
Iterateur_PolyMAC_CDO_Elem::contribuer_au_second_membre_bords
void contribuer_au_second_membre_bords(DoubleTab &) const
Definition Iterateur_PolyMAC_CDO_Elem.h:1148
Iterateur_PolyMAC_CDO_Elem::calculer
DoubleTab & calculer(const DoubleTab &, DoubleTab &) const override
Definition Iterateur_PolyMAC_CDO_Elem.h:975
Iterateur_PolyMAC_CDO_Elem::contribuer_au_second_membre_interne
void contribuer_au_second_membre_interne(DoubleTab &) const
Definition Iterateur_PolyMAC_CDO_Elem.h:1655
Iterateur_PolyMAC_CDO_Elem::ajouter_contribution_bords_vitesse
void ajouter_contribution_bords_vitesse(const DoubleTab &, Matrice_Morse &) const
Definition Iterateur_PolyMAC_CDO_Elem.h:2308
Iterateur_PolyMAC_CDO_Elem::modifier_flux
void modifier_flux() const
Definition Iterateur_PolyMAC_CDO_Elem.h:980
Iterateur_PolyMAC_CDO_Elem::ajouter_contribution_interne_vitesse
void ajouter_contribution_interne_vitesse(const DoubleTab &, Matrice_Morse &) const
Definition Iterateur_PolyMAC_CDO_Elem.h:2292
Iterateur_PolyMAC_CDO_Elem::ajouter_contribution_vitesse
void ajouter_contribution_vitesse(const DoubleTab &, Matrice_Morse &) const override
Definition Iterateur_PolyMAC_CDO_Elem.h:2267
Iterateur_PolyMAC_CDO_Elem::elem
IntTab elem
Definition Iterateur_PolyMAC_CDO_Elem.h:75
Iterateur_PolyMAC_CDO_Elem::completer_
void completer_() override
Definition Iterateur_PolyMAC_CDO_Elem.h:95
Iterateur_PolyMAC_CDO_Elem::ajouter_contribution_bords
void ajouter_contribution_bords(const DoubleTab &, Matrice_Morse &) const
Definition Iterateur_PolyMAC_CDO_Elem.h:1713
Iterateur_PolyMAC_CDO_base
Definition Iterateur_PolyMAC_CDO_base.h:57
Iterateur_PolyMAC_CDO_base::type_cl
Type_Cl_PolyMAC_CDO type_cl(const Cond_lim &) const
Definition Iterateur_PolyMAC_CDO_base.h:85
Matrice_Morse
Classe Matrice_Morse Represente une matrice M (creuse), non necessairement carree.
Definition Matrice_Morse.h:50
Milieu_base
classe Milieu_base Cette classe est la base de la hierarchie des milieux (physiques)
Definition Milieu_base.h:50
Neumann_paroi_adiabatique
Classe Neumann_paroi_adiabatique Cette condition limite correspond a une paroi adiabatique dans une.
Definition Neumann_paroi_adiabatique.h:29
Neumann_paroi
Classe Neumann_paroi Cette condition limite correspond a un flux impose pour l'equation de.
Definition Neumann_paroi.h:29
Neumann_sortie_libre
classe Neumann_sortie_libre Cette classe represente une frontiere ouverte sans vitesse imposee
Definition Neumann_sortie_libre.h:34
Nom
class Nom Une chaine de caractere pour nommer les objets de TRUST
Definition Nom.h:31
Objet_U::duplique
virtual int duplique() const =0
Objet_U::taille_memoire
virtual unsigned taille_memoire() const =0
Objet_U::numero
int numero() const
Renvoie l'indice de l'objet dans Memoire::data.
Definition Objet_U.cpp:268
Objet_U::dimension
static int dimension
Definition Objet_U.h:99
Objet_U::Sortie
friend class Sortie
Definition Objet_U.h:75
Objet_U::le_nom
virtual const Nom & le_nom() const
Donne le nom de l'Objet_U Methode a surcharger : renvoie "neant" dans cette implementation.
Definition Objet_U.cpp:319
Objet_U::bidim_axi
static int bidim_axi
Definition Objet_U.h:102
Op_Conv_PolyMAC_CDO_iterateur_base
Definition Op_Conv_PolyMAC_CDO_iterateur_base.h:24
Operateur_Diff_base
classe Operateur_Diff_base Cette classe est la base de la hierarchie des operateurs representant
Definition Operateur_Diff_base.h:37
Operateur_base
classe Operateur_base Classe est la base de la hierarchie des objets representant un
Definition Operateur_base.h:48
Periodique
classe Periodique Cette classe represente une condition aux limites periodique.
Definition Periodique.h:31
Process::mp_sum_for_each_item
static void mp_sum_for_each_item(TRUSTArray< _TYPE_ > &x, int n=-1)
Definition Process.cpp:193
Process::abort
static void abort()
Routine de sortie de Trio-U sur une erreur abort().
Definition Process.cpp:570
Process::exit
static void exit(int exit_code=-1)
Routine de sortie de TRUST dans une region Kokkos.
Definition Process.cpp:455
Process::je_suis_maitre
static int je_suis_maitre()
renvoie 1 si on est sur le processeur maitre du groupe courant (c'est a dire me() == 0),...
Definition Process.cpp:86
SFichier
Cette classe est a la classe C++ ofstream ce que la classe Sortie est a la classe C++ ostream Elle re...
Definition SFichier.h:27
Schema_Temps_base
class Schema_Temps_base
Definition Schema_Temps_base.h:67
Schema_Temps_base::temps_courant
double temps_courant() const
Renvoie le temps courant.
Definition Schema_Temps_base.h:445
Sortie
Classe de base des flux de sortie.
Definition Sortie.h:52
Symetrie
classe Symetrie Sur les faces de symetrie on a les proprietes suivantes:
Definition Symetrie.h:37
TRUSTTab::nb_dim
int nb_dim() const
Definition TRUSTTab.h:199
TRUSTTab::resize
void resize(_SIZE_ n, RESIZE_OPTIONS opt=RESIZE_OPTIONS::COPY_INIT)
Definition TRUSTTab.tpp:469
TRUSTTab::dimension
_SIZE_ dimension(int d) const
Definition TRUSTTab.tpp:133