TRUST 1.9.8
HPC thermohydraulic platform
Loading...
Searching...
No Matches
Op_Conv_DI_L2_VEF_Face.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
17#include <Op_Conv_DI_L2_VEF_Face.h>
18#include <Champ_P1NC.h>
19#include <Schema_Temps_base.h>
20#include <Periodique.h>
21#include <Neumann_sortie_libre.h>
22
23Implemente_instanciable(Op_Conv_DI_L2_VEF_Face,"Op_Conv_DI_L2_VEF_P1NC",Op_Conv_VEF_base);
24// XD convection_di_l2 convection_deriv di_l2 NO_BRACE Only for VEF discretization.
25
26//// printOn
27//
28
29Sortie& Op_Conv_DI_L2_VEF_Face::printOn(Sortie& s ) const
30{
31 return s << que_suis_je() ;
32}
33
34//// readOn
35//
36
37Entree& Op_Conv_DI_L2_VEF_Face::readOn(Entree& s )
38{
39 return s ;
40}
41
42//
43// Fonctions de la classe Op_Conv_DI_L2_VEF_Face
44//
45
46void Op_Conv_DI_L2_VEF_Face::associer_vitesse(const Champ_base& vit)
47{
48 const Champ_P1NC& inco = ref_cast(Champ_P1NC,vit);
49 vitesse_= inco;
50}
51
52
53void flora(DoubleTab A, int& N , DoubleVect B, DoubleVect& U, int& test_flora)
54{
55 //cette procedure correspond a la methode d'elimination de Gauss
56 test_flora = 0; //test pour savoir si la matrice est inversible(0) ou non(1)
57 int N1 = N-1;
58 int m, m1, i, i1, k, l;
59 double quo,SU;
60 for(m=0; m<N1; m++)
61 {
62 m1 = m+1;
63 if(std::fabs(A(m,m)) > 1e-20)
64 {
65 for(i=m1; i<N; i++)
66 {
67 quo = A(i,m)/A(m,m);
68 for(k=m; k<N; k++)
69 A(i,k) =A(i,k)-quo*A(m,k);
70 B(i) = B(i)-quo*B(m);
71 }
72 }
73 else
74 {
75 Cerr<<"Erreur flora: matrice non inversible a l'indice "<<m<<finl;
76 test_flora = 1;
77 }
78
79 }
80 if(std::fabs(A(N-1,N-1)) >= 1e-20)
81 {
82 U(N-1) = B(N-1)/A(N-1,N-1);
83 for(l=0; l<N1; l++)
84 {
85 i = N1-l-1;
86 SU = 0;
87 i1 = i+1;
88 for(k=i1; k<N; k++)
89 SU = SU+A(i,k)*U(k);
90 U(i) = (B(i)-SU)/A(i,i);
91 }
92 }
93 else
94 {
95 Cerr<<"Erreur flora: matrice non inversible a l'indice"<<N-1<<finl;
96 test_flora = 1;
97 }
98}
99
100
101void flora_p(DoubleTab& A, int& N, DoubleVect& B, DoubleVect& U, int& test_flora)
102{
103 //cette procedure correspond a la methode d'elimination de Gauss
104 test_flora = 1;//test pour savoir si la matrice est inversible(1) ou non(0)
105 int N1 = N-1;
106 int m, m1, i, j, i1, k, l;
107 int test=0;
108 double quo, SU, x, y;
109 for(m=0; m<N1; m++)
110 {
111 m1 = m+1;
112 if(std::fabs(A(m,m)) >= 1.e-10)
113 {
114 for(i=m1; i<N; i++)
115 {
116 quo = A(i,m)/A(m,m);
117 for(k=m; k<N; k++)
118 A(i,k) =A(i,k)-quo*A(m,k);
119 B(i) = B(i)-quo*B(m);
120 }
121 }
122 else
123 {
124 j=m+1;
125 while(j<N && test == 0)
126 {
127 if(std::fabs(A(m,j)) >= 1.e-10) test =1;
128 j++;
129 }
130 if(j == N)
131 {
132 Cerr<<"Erreur flora: matrice non inversible a l'indice "<<m<<finl;
133 test_flora = 1;
134 }
135 else //echange des colonnes m et j
136 {
137 for(i=0; i<N; i++)
138 {
139 x=A(i,m);
140 A(i,m)=A(i,j);
141 A(i,j)=x;
142 y=B(m);
143 B(m)=B(j);
144 B(j)=y;
145 }
146 }
147 test_flora = 0;
148
149 }
150
151 }
152 if(std::fabs(A(N-1,N-1)) >= 1.e-10)
153 {
154 U(N-1) = B(N-1)/A(N-1,N-1);
155 for(l=0; l<N1; l++)
156 {
157 i = N1-l-1;
158 SU = 0;
159 i1 = i+1;
160 for(k=i1; k<N; k++)
161 SU = SU+A(i,k)*U(k);
162 U(i) = (B(i)-SU)/A(i,i);
163 }
164 }
165 else
166 {
167 //Cerr<<"Erreur flora: matrice non inversible a l'indice"<<N-1<<finl;
168 test_flora = 0;
169 }
170}
171
172
173void qrdcmp(DoubleTab& A, int& N, DoubleVect& C, DoubleVect& D, int& sing)
174{
175 //construit la decomposition QR de A. Le triangle superieur R est stockee dans le triangle superieur de A,
176 //exepte les elements de la diagonale, qui sont stockes dans D. La matrice orthogonale Q est representee
177 //comme un produit de N-1 matrice Q(1), ..., Q(N-1) ou Q(j) = Id-(uj*ujt)/cj. La ieme composante de uj est 0
178 //pour i=1, ...,j-1 et A(i,j) pour i=j, ..., N. sing retourne 0 si la decomposition est possible et 1 sinon.
179
180 int i, j, k;
181 double scale, sigma, sum, tau;
182 sing =0;
183 for(k=0; k<N-1; k++)
184 {
185 scale = 0.;
186 for(i=k; i<N; i++)
187 if(scale < std::fabs(A(i,k))) scale = std::fabs(A(i,k));
188 if(std::fabs(scale)<1.e-15)
189 {
190 //cas singulier
191 Cerr << " huhu " << finl ;
192 sing = 1;
193 C(k) = 0;
194 D(k) = 0;
195 return ;
196 }
197 else
198 {
199 for(i=k; i<N; i++) A(i,k) /= scale;
200 for(sum=0.0,i=k; i<N; i++) sum += (A(i,k) * A(i,k));
201 if(A(k,k)>0) sigma = sqrt(sum);
202 else sigma = -1*sqrt(sum);
203 A(k,k) += sigma;
204 C(k) = sigma*A(k,k);
205 D(k) = -1*scale*sigma;
206 for(j=k+1; j<N; j++)
207 {
208 for(sum=0.0,i=k; i<N; i++) sum += A(i,k)*A(i,j);
209 tau = sum/C(k);
210 for(i=k; i<N; i++) A(i,j) -= tau*A(i,k);
211 }
212 }
213 }
214 D(N-1) = A(N-1,N-1);
215 if(std::fabs(D(N-1)) <1.e-12)
216 {
217 Cerr << " hoho " << finl ;
218 sing =1;
219 }
220}
221
222
223void rsolv(DoubleTab& A, int& N, DoubleVect& D, DoubleVect& B)
224{
225 //resout le systeme Rx=B, ou R est triangulaire superieure stockee dans A et D, provenant de qrdcmp
226 //le resultat est stocke dans B.
227 int i, j;
228 double sum;
229 B(N-1) /= D(N-1);
230 for(i=N-2; i>=0; i--)
231 {
232 for(sum=0.0,j=i+1; j<N; j++) sum += A(i,j)*B(j);
233 B(i) = (B(i)-sum)/D(i);
234 }
235}
236
237void qrsolv( DoubleTab& A, int& N, DoubleVect& B, DoubleVect& X, int& sing,
238 int& ncomp, DoubleVect& C, DoubleVect& D)
239{
240 //resout le systeme lineaire Ax=B
241 // DoubleVect C(N), D(N);
242 int i, j;
243 double sum, tau;
244 if(ncomp == 0 ) qrdcmp(A, N, C, D, sing);
245 if(sing == 0)
246 {
247 for(j=0; j<N-1; j++)
248 {
249 for(sum=0.0,i=j; i<N; i++) sum += A(i,j)*B(i);
250 tau = sum/C(j);
251 for(i=j; i<N; i++) B(i) -= tau*A(i,j);
252 }
253 rsolv(A, N, D, B);//resout Rx=QtB
254 for(i=0; i<N; i++) X(i) = B(i);
255 }
256 //else Cerr<<"erreur"<<finl;
257}
258
259
260
261
262//methode du gradient biconjugue
263void gradient_biconjugue(DoubleTab A, int n, DoubleVect b, DoubleVect& x, int& sing, int& niter)
264{
265 if (Process::is_parallel())
266 {
267 Cerr << "OpVEF_DI_L2.cpp: gradient_biconjugue() n'est pas parallele" << finl;
268 assert(0);
269 Process::exit();
270 }
271
272 double seuil ;
273 double dnew = 1. ;
274
275 double dold, alfa, beta ;
276 niter = 0 ;
277 sing = 0 ;
278
279 DoubleVect r(n) ;
280 DoubleVect r_tilda(n);
281 DoubleVect p(n) ;
282 DoubleVect p_tilda(n);
283 DoubleVect q(n) ;
284 //DoubleVect u(n);
285 double r_norme, b_norme = norme_array(b);
286
287 int i,j;
288
289 if(b_norme > 1.e-7 )
290 seuil = 1.e-5/b_norme;
291 else seuil = 1.e-10;
292
293 r = 0. ;
294 p = 0. ;
295 q = 0. ;
296
297 for(i=0; i<n; i++)
298 p(i) = b(i);
299
300 int nmax = 50 ;
301
302 for(i=0; i<n; i++)
303 for(j=0; j<n; j++)
304 {
305 r(i) = p(i)-A(i,j)*x(j) ;
306 r_tilda(i) = r(i) ;
307 }
308
309 p = 0 ;
310 dold = dotproduct_array(r_tilda, r);
311 r_norme = norme_array(r);
312
313 if(sqrt(dold) > seuil)
314 {
315 while ( ( r_norme > seuil ) && (niter++ < nmax) )
316 {
317 assert(Process::is_sequential()); // B.M. code visiblement faux en parallele
318 dnew = dotproduct_array(r_tilda, r);
319
320 if(dold == 0.) niter = nmax ;
321 else
322 {
323 beta = dnew/dold ;
324 for(i=0; i<n; i++)
325 {
326 p(i) = r(i)+beta*p(i) ;
327 p_tilda(i) = r_tilda(i)+beta*p_tilda(i) ;
328 }
329
330 q = 0. ;
331 for(i=0; i<n; i++)
332 for(j=0; j<n; j++)
333 q(i) += A(i,j)*p(j) ;
334
335 beta = dotproduct_array(p_tilda, q) ;
336 //Cerr<<"beta :"<<beta<<finl ;
337 if(beta == 0.) niter = nmax ;
338 else
339 {
340 alfa = dnew / beta ;
341
342 x.ajoute_sans_ech_esp_virt(alfa, p);
343 r.ajoute_sans_ech_esp_virt(-alfa, q);
344
345 q = 0. ;
346 for(i=0; i<n; i++)
347 for(j=0; j<n; j++)
348 q(i) += A(j,i)*p_tilda(j) ;
349
350 r_tilda.ajoute_sans_ech_esp_virt(-alfa, q);
351
352 dold = dnew;
353 //assert(dnew >= 0.);
354 r_norme = norme_array(r);
355 //Cerr << " niter " << niter << " dnew " << dnew <<finl ;
356 //Cerr << " r " << r.norme() << finl;
357 //Cerr << " x " << x << finl;
358 }
359 }
360 }
361 }
362 if ( niter >= nmax) sing = 1 ;
363 //Cerr<<"niter :"<<niter<<finl;
364}
365
366
367// convbis correspond au calcul de -1*terme_convection
368
369void convbis(double psc,int num1,int num2,
370 const DoubleTab& transporte, int ncomp,
371 DoubleTab& resu, DoubleVect& fluent)
372{
373 int comp,amont;
374 double flux;
375
376 if (psc >= 0)
377 {
378 amont = num1;
379 fluent[num2] += psc;
380 }
381 else
382 {
383 amont = num2;
384 fluent[num1] -= psc;
385 }
386
387 if (ncomp == 1)
388 {
389 flux = transporte(amont)*psc;
390 resu(num1) -= flux;
391 resu(num2) += flux;
392 }
393 else
394 for (comp=0; comp<ncomp; comp++)
395 {
396 flux = transporte(amont,comp)*psc;
397 resu(num1,comp) -= flux;
398 resu(num2,comp) += flux;
399 }
400}
401
402
403DoubleTab& Op_Conv_DI_L2_VEF_Face::ajouter(const DoubleTab& transporte,
404 DoubleTab& resu) const
405{
406 const Domaine_Cl_VEF& domaine_Cl_VEF = la_zcl_vef.valeur();
407 const Domaine_VEF& domaine_VEF = le_dom_vef.valeur();
408 const Champ_Inc_base& la_vitesse=vitesse_.valeur();
409
410 const IntTab& elem_faces = domaine_VEF.elem_faces();
411 const DoubleTab& face_normales = domaine_VEF.face_normales();
412 const auto& facette_normales = domaine_VEF.facette_normales();
413 // const DoubleVect& volumes_entrelaces = domaine_VEF.volumes_entrelaces();
414 const Domaine& domaine = domaine_VEF.domaine();
415 // const int nb_faces = domaine_VEF.nb_faces();
416 const int nfa7 = domaine_VEF.type_elem().nb_facette();
417 // const int nb_elem = domaine_VEF.nb_elem();
418 const int nb_elem_tot = domaine_VEF.nb_elem_tot();
419 const IntVect& rang_elem_non_std = domaine_VEF.rang_elem_non_std();
420 const IntTab& face_voisins = domaine_VEF.face_voisins();
421 const DoubleVect& porosite_face = equation().milieu().porosite_face();
422
423 /* const IntTab& face_voisins = domaine_VEF.face_voisins();
424 int jjj;
425 for(jjj = 0;jjj<domaine_VEF.nb_faces_tot();jjj++)
426 Cerr<<"face_voisins : face = "<<jjj<<" : "<<face_voisins(jjj,0)<<", "<<face_voisins(jjj,1)<<finl;*/
427
428
429 const DoubleTab& normales_facettes_Cl = domaine_Cl_VEF.normales_facettes_Cl();
430 // const DoubleVect& volumes_entrelaces_Cl = domaine_Cl_VEF.volumes_entrelaces_Cl();
431
432 int nfac = domaine.nb_faces_elem();
433 int nsom = domaine.nb_som_elem();
434 int nb_som_facette = domaine.type_elem()->nb_som_face();
435
436 // int premiere_face_int = domaine_VEF.premiere_face_int();
437
438 // Pour le traitement de la convection on distingue les polyedres
439 // standard qui ne "voient" pas les conditions aux limites et les
440 // polyedres non standard qui ont au moins une face sur le bord.
441 // Un polyedre standard a n facettes sur lesquelles on applique le
442 // schema de convection.
443 // Pour un polyedre non standard qui porte des conditions aux limites
444 // de Dirichlet, une partie des facettes sont portees par les faces.
445 // En bref pour un polyedre le traitement de la convection depend
446 // du type (triangle, tetraedre ...) et du nombre de faces de Dirichlet.
447
448 double psc;
449 int poly,face_adj,fa7,i,j,n_bord;
450 int num_face, rang ,itypcl;
451 int num10,num20;
452 int ncomp_ch_transporte, first;
453 if (transporte.nb_dim() == 1)
454 ncomp_ch_transporte=1;
455 else
456 ncomp_ch_transporte= transporte.dimension(1);
457
458 IntVect face(nfac);
459 DoubleVect vs(dimension);
460 DoubleVect vc(dimension);
461 DoubleTab vsom(nsom,dimension);
462 DoubleVect cc(dimension);
463 DoubleTab derive(1,1) ;
464 int N ,M , sing , cal_amont;
465 DoubleVect trans(ncomp_ch_transporte) ;
466 if (dimension == 2)
467 {
468 N = 5 ;
469 M = 8 ;
470 derive.resize(N,ncomp_ch_transporte) ;
471 }
472 else // (dimension == 3)
473 {
474 N = 9 ;
475 M = 15 ;
476 derive.resize(N,ncomp_ch_transporte) ;
477 }
478 // On remet a zero le tableau qui sert pour
479 // le calcul du pas de temps de stabilite
480 fluent_ = 0;
481
482 // Traitement particulier pour les faces de periodicite
483
484 int nb_faces_perio = 0;
485 // Boucle pour compter le nombre de faces de periodicite
486 for (n_bord=0; n_bord<domaine_VEF.nb_front_Cl(); n_bord++)
487 {
488 const Cond_lim& la_cl = domaine_Cl_VEF.les_conditions_limites(n_bord);
489 if (sub_type(Periodique,la_cl.valeur()))
490 {
491 const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
492 nb_faces_perio += le_bord.nb_faces();
493 }
494 }
495
496 DoubleTab tab;
497 if (ncomp_ch_transporte == 1)
498 tab.resize(nb_faces_perio);
499 else
500 tab.resize(nb_faces_perio,ncomp_ch_transporte);
501
502 // Boucle pour remplir tab
503 nb_faces_perio=0;
504 for (n_bord=0; n_bord<domaine_VEF.nb_front_Cl(); n_bord++)
505 {
506 const Cond_lim& la_cl = domaine_Cl_VEF.les_conditions_limites(n_bord);
507 if (sub_type(Periodique,la_cl.valeur()))
508 {
509 // const Periodique& la_cl_perio = ref_cast(Periodique, la_cl.valeur());
510 const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
511 int num1 = le_bord.num_premiere_face();
512 int num2 = num1 + le_bord.nb_faces();
513 for (num_face=num1; num_face<num2; num_face++)
514 {
515 if (ncomp_ch_transporte == 1)
516 tab(nb_faces_perio) = resu(num_face);
517 else
518 for (int comp=0; comp<ncomp_ch_transporte; comp++)
519 tab(nb_faces_perio,comp) = resu(num_face,comp);
520 nb_faces_perio++;
521 }
522 }
523 }
524
525
526 // Les polyedres non standard sont ranges en 2 groupes dans le Domaine_VEF:
527 // - polyedres bords et joints
528 // - polyedres bords et non joints
529 // On traite les polyedres en suivant l'ordre dans lequel ils figurent
530 // dans le domaine
531
532 // boucle sur les polys
533
534 int nlim = -1 ;
535 const IntTab& KEL=domaine_VEF.type_elem().KEL();
536 for (poly=0; poly<nb_elem_tot; poly++)
537 {
538
539 rang = rang_elem_non_std(poly);
540 if (rang==-1)
541 itypcl=0;
542 else
543 itypcl=domaine_Cl_VEF.type_elem_Cl(rang);
544
545 // calcul des numeros des faces du polyedre
546 for (face_adj=0; face_adj<nfac; face_adj++)
547 face[face_adj]= elem_faces(poly,face_adj);
548
549 // calcul de la vitesse aux sommets des polyedres
550 for (j=0; j<dimension; j++)
551 {
552 vs[j] = la_vitesse.valeurs()(face[0],j);
553 for (i=1; i<nfac; i++)
554 vs[j]+= la_vitesse.valeurs()(face[i],j);
555 }
556 for (i=0; i<nsom; i++)
557 for (j=0; j<dimension; j++)
558 vsom(i,j) = vs[j] - dimension*la_vitesse.valeurs()(face[i],j);
559
560 // calcul de vc
561 domaine_VEF.type_elem().calcul_vc(face,vc,vs,vsom,vitesse(),
562 itypcl,porosite_face);
563
564
565 /* Cerr<<"premiere_face_int = "<<premiere_face_int<<", domaine_VEF.nb_faces_joint() = "<<domaine_VEF.nb_faces_joint()<<finl;
566 Cerr<<"premiere_face_std() = "<<domaine_VEF.premiere_face_std()<<finl;
567 Cerr<<"nb_elem() = "<<domaine_VEF.nb_elem()<<", nb_elem_tot() = "<<domaine_VEF.nb_elem_tot()<<finl;
568 Cerr<<"nb_faces() = "<<domaine_VEF.nb_faces()<<", nb_faces_tot() = "<<domaine_VEF.nb_faces_tot()<<finl;
569 Cerr<<"nb_faces_bord() = "<<domaine_VEF.nb_faces_bord()<<", premiere_face_bord() = "<<domaine_VEF.premiere_face_bord()<<finl;*/
570
571 cal_amont = 0 ;
572
573 int elem0,elem1,face_adj_glob;
574 for (face_adj=0; face_adj<nfac; face_adj++)
575 {
576 face_adj_glob = face[face_adj];
577 elem0 = face_voisins(face_adj_glob,0);
578 elem1 = face_voisins(face_adj_glob,1);
579 if ((elem0 == -1) || (elem1 == -1))
580 cal_amont++ ;
581 }
582
583 // calcul polynom de reconstruction
584 if(dimension == 2 )
585 {
586 if ( cal_amont == 0) poly_DI_L2_2d(N,M,derive,poly,ncomp_ch_transporte,transporte,sing);
587 }
588 else
589 {
590 if ( cal_amont == 0) poly_DI_L2_3d(N,M,derive,poly,ncomp_ch_transporte,transporte,sing);
591 }
592
593 // Boucle sur les facettes du polyedre
594
595 for (fa7=0; fa7<nfa7; fa7++)
596 {
597 if (rang==-1)
598 for (i=0; i<dimension; i++)
599 cc[i] = facette_normales(poly,fa7,i);
600 else
601 for (i=0; i<dimension; i++)
602 cc[i] = normales_facettes_Cl(rang,fa7,i);
603
604 // On applique le schema de convection a chaque sommet de la facette
605
606 // reconstruction seulement wenn first = 0
607
608 first = -1 ;
609
610 // On traite le ou les sommets qui sont aussi des sommets du polyedre
611 for (i=0; i<nb_som_facette-1; i++)
612 {
613 first++ ;
614
615 psc =0;
616 for (j=0; j<dimension; j++)
617 psc+= (vc(j)/double(nb_som_facette-1)+vsom(KEL(i+2,fa7),j))*cc[j];
618 psc /= nb_som_facette;
619
620 num10 = face[KEL(0,fa7)];
621 num20 = face[KEL(1,fa7)];
622
623 if ( cal_amont > 0)
624 {
625 convbis(psc,num10,num20,transporte,ncomp_ch_transporte,resu,fluent_);
626 }
627 else
628 {
629 if(dimension == 2 )
630 {
631 reconst_DI_L2_2d(derive,poly,psc,num10,num20,transporte,ncomp_ch_transporte,resu,fluent_,sing,
632 nlim );
633 }
634 else
635 {
636 reconst_DI_L2_3d(derive,poly,psc,num10,num20,transporte,ncomp_ch_transporte,resu,fluent_,sing,
637 first,trans,nlim);
638 }
639 }
640 }
641
642 }
643
644 } // fin de la boucle
645
646 Cerr << " limitiert in " << nlim << " valeurs " << finl ;
647
648 int voisine;
649 nb_faces_perio = 0;
650 double diff1,diff2;
651
652 // Dimensionnement du tableau des flux convectifs au bord du domaine
653 // de calcul
654 DoubleTab& flux_b = flux_bords_;
655 flux_b.resize(domaine_VEF.nb_faces_bord(),ncomp_ch_transporte);
656 flux_b = 0.;
657
658 // Boucle sur les bords pour traiter les conditions aux limites
659 // il y a prise en compte d'un terme de convection pour les
660 // conditions aux limites de Neumann_sortie_libre seulement
661
662 for (n_bord=0; n_bord<domaine_VEF.nb_front_Cl(); n_bord++)
663 {
664
665 const Cond_lim& la_cl = domaine_Cl_VEF.les_conditions_limites(n_bord);
666
667 if (sub_type(Neumann_sortie_libre,la_cl.valeur()))
668 {
669 const Neumann_sortie_libre& la_sortie_libre = ref_cast(Neumann_sortie_libre, la_cl.valeur());
670 const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
671 int num1 = le_bord.num_premiere_face();
672 int num2 = num1 + le_bord.nb_faces();
673 for (num_face=num1; num_face<num2; num_face++)
674 {
675 psc =0;
676 for (i=0; i<dimension; i++)
677 psc += la_vitesse.valeurs()(num_face,i)*face_normales(num_face,i);
678 if (psc>0)
679 if (ncomp_ch_transporte == 1)
680 {
681 resu(num_face) -= psc*transporte(num_face);
682 flux_b(num_face,0) -= psc*transporte(num_face);
683 }
684 else
685 for (i=0; i<ncomp_ch_transporte; i++)
686 {
687 resu(num_face,i) -= psc*transporte(num_face,i);
688 flux_b(num_face,i) -= psc*transporte(num_face,i);
689 }
690 else
691 {
692 if (ncomp_ch_transporte == 1)
693 {
694 resu(num_face) -= psc*la_sortie_libre.val_ext(num_face-num1);
695 flux_b(num_face,0) -= psc*la_sortie_libre.val_ext(num_face-num1);
696 }
697 else
698 for (i=0; i<ncomp_ch_transporte; i++)
699 {
700 resu(num_face,i) -= psc*la_sortie_libre.val_ext(num_face-num1,i);
701 flux_b(num_face,i) -= psc*la_sortie_libre.val_ext(num_face-num1);
702 }
703 fluent_[num_face] -= psc;
704 }
705 }
706 }
707 else if (sub_type(Periodique,la_cl.valeur()))
708 {
709 const Periodique& la_cl_perio = ref_cast(Periodique, la_cl.valeur());
710 const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
711 int num1 = le_bord.num_premiere_face();
712 int num2 = num1 + le_bord.nb_faces();
713 IntVect fait(le_bord.nb_faces());
714 fait = 0;
715 for (num_face=num1; num_face<num2; num_face++)
716 {
717 if (fait[num_face-num1] == 0)
718 {
719 voisine = la_cl_perio.face_associee(num_face-num1) + num1;
720
721 if (ncomp_ch_transporte == 1)
722 {
723 diff1 = resu(num_face)-tab(nb_faces_perio);
724 diff2 = resu(voisine)-tab(nb_faces_perio+voisine-num_face);
725 resu(voisine) += diff1;
726 resu(num_face) += diff2;
727 flux_b(voisine,1) += diff1;
728 flux_b(num_face,0) += diff2;
729 }
730 else
731 for (int comp=0; comp<ncomp_ch_transporte; comp++)
732 {
733 diff1 = resu(num_face,comp)-tab(nb_faces_perio,comp);
734 diff2 = resu(voisine,comp)-tab(nb_faces_perio+voisine-num_face,comp);
735 resu(voisine,comp) += diff1;
736 resu(num_face,comp) += diff2;
737 flux_b(voisine,comp) += diff1;
738 flux_b(num_face,comp) += diff2;
739 }
740
741 fait[num_face-num1]= 1;
742 fait[voisine-num1] = 1;
743 }
744 nb_faces_perio++;
745 }
746 }
747 }
748 modifier_flux(*this);
749 return resu;
750
751}
752
753
754void Op_Conv_DI_L2_VEF_Face::reconst_DI_L2_2d(DoubleTab& derive ,int poly,
755 double psc,int num1,int num2,
756 const DoubleTab& transporte,
757 int ncomp,
758 DoubleTab& resu,
759 DoubleVect& tab_fluent, int sing, int& nlim) const
760{
761
762 //Cerr << " in reconst_DI_L2_2d " << sing << finl ;
763
764 const Domaine_VEF& domaine_VEF = le_dom_vef.valeur();
765
766 const IntTab& elem_faces = domaine_VEF.elem_faces();
767 // const IntTab& face_voisins = domaine_VEF.face_voisins();
768
769 const Domaine& domaine = domaine_VEF.domaine();
770
771 int nfac = domaine.nb_faces_elem();
772
773 const DoubleTab& xv = domaine_VEF.xv();
774 const DoubleTab& xp = domaine_VEF.xp();
775
776 double dt=equation().schema_temps().pas_de_temps();
777 const Champ_Inc_base& la_vitesse=vitesse_.valeur();
778
779 int i, j, face_adj , face_glob, numg=-1 ;
780 DoubleVect dist(dimension) ;
781 DoubleVect trans_c_g(ncomp) ;
782 DoubleVect vs(ncomp) ;
783 DoubleVect trans(ncomp) ;
784 DoubleTab coor_trans(dimension,dimension) ;
785
786 trans_c_g = 0. ;
787
788 for (j=0; j<ncomp; j++)
789 {
790 for(face_adj=0; face_adj<nfac; face_adj ++)
791 {
792 face_glob = elem_faces(poly, face_adj);
793
794 if (ncomp == 1) trans_c_g(0) += transporte(face_glob)/double(nfac) ;
795 else trans_c_g(j) += transporte(face_glob,j)/double(nfac) ;
796
797 if ((face_glob != num1) && (face_glob != num2)) numg = face_glob ;
798
799 }
800 }
801 int num3 = elem_faces(poly, 0 ) ;
802
803 //coor_trans(0,0) = cos(3.14159/4.) ;
804 //coor_trans(0,1) = cos(3.14159/4.) ;
805 //coor_trans(1,0) = cos(3.14159/2.-3.14159/4.) ;
806 //coor_trans(1,1) = cos(3.14159/2.+3.14159/4.) ;
807 coor_trans(0,0) = 1. ;
808 coor_trans(0,1) = 0. ;
809 coor_trans(1,0) = 0. ;
810 coor_trans(1,1) = 1. ;
811
812 for (i=0; i<dimension; i++) vs(i)= (la_vitesse.valeurs()(num1,i)+la_vitesse.valeurs()(num2,i)) / 2.;
813
814 dist = 0. ;
815
816 for (i=0; i<dimension; i++)
817 for (j=0; j<dimension; j++)
818 dist(i) += (2.*xp(poly,j) - xv(numg,j) - xv(num3,j) - vs(j)*dt/2. ) * coor_trans(i,j) ;
819
820 double dtrans_max, dtrans_min, dtrans_cen ;
821
822 int amont, aval ;
823
824 if (psc >= 0)
825 {
826 amont = num1;
827 aval = num2;
828 }
829 else
830 {
831 amont = num2;
832 aval = num1;
833 }
834
835 for (j=0; j<ncomp; j++)
836 {
837
838 if (sing == 0 )
839 {
840 if(ncomp == 1)
841 {
842 trans(0) = transporte(num3) + derive(4)*dist(0) * coor_trans(0,0)
843 + derive(3)*dist(1) * coor_trans(1,1)
844 + 1./2.*( derive(2)*dist(0)*dist(0) + derive(1)*dist(1)*dist(1) )
845 + derive(0)*dist(0)*dist(1) ;
846
847 dtrans_cen = transporte(amont) + transporte(aval) - trans_c_g(0) ;
848 dtrans_max = std::max( transporte(amont), dtrans_cen ) ;
849 dtrans_min = std::min( transporte(amont), dtrans_cen ) ;
850 }
851 else
852 {
853 trans(j) = transporte(num3,j) + derive(4,j)*dist(0) * coor_trans(0,0)
854 + derive(3,j)*dist(1) * coor_trans(1,1)
855 + 1./2.*( derive(2,j)*dist(0)*dist(0) + derive(1,j)*dist(1)*dist(1) )
856 + derive(0,j)*dist(0)*dist(1) ;
857
858
859 //dtrans_cen = transporte(num1,j) + transporte(num2,j) - trans_c_g(j) ;
860 //dtrans_cen = 1./2.*(transporte(amont,j) + transporte(aval,j)) ;
861 dtrans_cen = transporte(num1,j) + transporte(num2,j) - trans_c_g(j) ;
862
863 //dtrans_max = std::max( transporte(amont,j), transporte(aval,j)) ;
864 dtrans_max = std::max( transporte(amont,j) , dtrans_cen) ;
865 dtrans_min = std::min( transporte(amont,j) , dtrans_cen) ;
866 //dtrans_min = std::min( dtrans_min, dtrans_cen) ;
867
868 }
869
870 if (trans(j) > dtrans_max )
871 {
872 //Cerr << "lim_max "<< finl ;
873 nlim++ ;
874 trans(j) = dtrans_max ;
875 }
876 if (trans(j) < dtrans_min )
877 {
878 //Cerr << "lim_min "<< finl ;
879 nlim++ ;
880 trans(j) = dtrans_min ;
881 }
882 }
883 else
884 {
885 if(ncomp == 1) trans(0) = ( transporte(num1) + transporte(num2) ) - trans_c_g(0) ;
886 else trans(j) = ( transporte(num1,j) + transporte(num2,j) ) - trans_c_g(j) ;
887 Cerr << " singx != 0 " << finl ;
888 }
889 }
890
891 // double kwave = 1. ;
892 // double lo_x = xv(num1,0)+xv(num2,0)-xp(poly,0) ;
893 // double lo_y = xv(num1,1)+xv(num2,1)-xp(poly,1) ;
894
895 // double exact0= -sin(kwave*lo_x)*cos(kwave*lo_y);
896 // double exact1= cos(kwave*lo_x)*sin(kwave*lo_y);
897 // double exact0 = kwave*lo_x*lo_x ;
898 // double exact1 = kwave*lo_y*lo_y ;
899
900 // double moyenne0 = ( transporte(num1,0) + transporte(num2,0) ) - trans_c_g(0) ;
901 // double moyenne1 = ( transporte(num1,1) + transporte(num2,1) ) - trans_c_g(1) ;
902
903 // double xv1 = (xv(num1,0)+xv(num2,0)-xp(poly,0)) - (xv(num3,0) + dist(0) - coor_trans(0) ) ;
904 // double xv2 = (xv(num1,1)+xv(num2,1)-xp(poly,1)) - (xv(num3,1) + dist(1) - coor_trans(1) ) ;
905
906 // double moyenne0 = transporte(numg,0) ;
907 // double moyenne1 = transporte(numg,1) ;
908
909 // if (sing == 0 ) Cerr << exact0 << " " << moyenne0 << " " << trans(0) << " "
910 // << exact1 << " " << moyenne1 << " " << trans(1) << finl ;
911
912 // if (sing == 0 ) Cerr << moyenne0 << " " << trans(0) << " "
913 // << moyenne1 << " " << trans(1) << finl ;
914 // Cerr << numg << " " << num3 << finl ;
915
916 double flux;
917
918 if (psc >= 0)
919 {
920 tab_fluent[num2] += psc;
921 }
922 else
923 {
924 tab_fluent[num1] -= psc;
925 }
926
927 if (ncomp == 1)
928 {
929 flux = trans(0)*psc;
930 resu(num1) -= flux;
931 resu(num2) += flux;
932 }
933 else
934 {
935 for (i=0; i<ncomp; i++)
936 {
937 flux = trans(i)*psc;
938 resu(num1,i) -= flux;
939 resu(num2,i) += flux;
940 }
941 }
942
943}
944
945void Op_Conv_DI_L2_VEF_Face::reconst_DI_L2_3d(DoubleTab& derive, int poly,
946 double psc,int num1,int num2,
947 const DoubleTab& transporte,
948 int ncomp,
949 DoubleTab& resu,
950 DoubleVect& tab_fluent ,
951 int sing, int first, DoubleVect& trans,
952 int& nlim) const
953{
954
955 const Domaine_VEF& domaine_VEF = le_dom_vef.valeur();
956
957 const IntTab& elem_faces = domaine_VEF.elem_faces();
958 // const IntTab& face_voisins = domaine_VEF.face_voisins();
959
960 // const int nb_faces = domaine_VEF.nb_faces();
961 // const int nfa7 = domaine_VEF.type_elem().nb_facette();
962 // const int nb_elem = domaine_VEF.nb_elem();
963 // const int nb_elem_tot = domaine_VEF.nb_elem_tot();
964 // int premiere_face_std = domaine_VEF.premiere_face_std();
965
966 const Domaine& domaine = domaine_VEF.domaine();
967
968 int nfac = domaine.nb_faces_elem();
969 // int nsom = domaine.nb_som_elem();
970 // int nb_som_facette = domaine.type_elem()->nb_som_face();
971
972 const DoubleTab& xv = domaine_VEF.xv();
973 const DoubleTab& xp = domaine_VEF.xp();
974
975
976 double dt=equation().schema_temps().pas_de_temps();
977 const Champ_Inc_base& la_vitesse=vitesse_.valeur();
978
979 double dtrans_max, dtrans_min, dtrans_cen;//, dtrans_aval ;
980
981 int i, j, face_adj , face_glob ;
982
983 DoubleVect dist(dimension) ;
984 DoubleTab coor_trans(dimension,dimension) ;
985
986 DoubleVect trans_c_g(ncomp) ;
987 double flux;
988 IntVect face(nfac) ;
989
990 DoubleVect vs(dimension) ;
991
992 int amont, aval ;
993
994 if (psc >= 0)
995 {
996 amont = num1;
997 aval = num2;
998 }
999 else
1000 {
1001 amont = num2;
1002 aval = num1;
1003 }
1004
1005 if(first == 0)
1006 {
1007 trans_c_g = 0. ;
1008
1009 for(face_adj=0; face_adj<nfac; face_adj ++)
1010 {
1011 face_glob = elem_faces(poly, face_adj );
1012 face(face_adj) = face_glob ;
1013 for (i=0; i<ncomp; i++)
1014 trans_c_g(i) += transporte(face_glob,i)/double(nfac) ;
1015
1016 }
1017 int num3 = elem_faces(poly, 0 ) ;
1018
1019 coor_trans(0,0) = cos(3.14159/4.) ;
1020 coor_trans(0,1) = cos(3.14159/4.) ;
1021 coor_trans(0,2) = cos(3.14159/2.) ;
1022
1023 coor_trans(1,0) = cos(3.14159/2.-3.14159/4.) ;
1024 coor_trans(1,1) = cos(3.14159/2.+3.14159/4.) ;
1025 coor_trans(1,2) = cos(3.14159/2.) ;
1026
1027 coor_trans(2,0) = cos(3.14159/2.) ;
1028 coor_trans(2,1) = cos(3.14159/2.) ;
1029 coor_trans(2,2) = cos(0.) ;
1030
1031 //coor_trans = 0.;
1032 //coor_trans(0,0) = 1. ;
1033 //coor_trans(1,1) = 1. ;
1034 //coor_trans(2,2) = 1. ;
1035
1036 vs =0. ;
1037 for(face_adj=0; face_adj<nfac; face_adj ++)
1038 {
1039 face_glob = elem_faces(poly, face_adj );
1040 face(face_adj) = face_glob ;
1041 for (i=0; i<dimension; i++)
1042 vs(i) -= la_vitesse.valeurs()(face_glob,i)/double(nfac) ;
1043
1044 }
1045 for (i=0; i<dimension; i++) vs(i) += la_vitesse.valeurs()(num1,i) + la_vitesse.valeurs()(num2,i) ;
1046
1047 dist = 0. ;
1048
1049 for(face_adj=0; face_adj<nfac; face_adj ++)
1050 {
1051 face_glob = elem_faces(poly, face_adj);
1052 if ((face_glob == num1) || (face_glob == num2 ))
1053 for (i=0; i<dimension; i++)
1054 for (j=0; j<dimension; j++)
1055 dist(i) += xv(face_glob,j) * coor_trans(i,j) ;
1056 }
1057
1058 for(i=0; i<dimension; i++)
1059 for (j=0; j<dimension; j++)
1060 dist(i) -= (xp(poly,j) + xv(num3,j) + vs(j)*dt/2.) * coor_trans(i,j) ;
1061
1062 for (j=0; j<ncomp; j++)
1063 {
1064 if (sing == 0 )
1065 {
1066 if(ncomp == 1)
1067 {
1068 trans(j) = transporte(num3,j) + derive(8)*dist(0) * coor_trans(0,0)
1069 + derive(7)*dist(1) * coor_trans(1,1)
1070 + derive(6)*dist(2) * coor_trans(2,2)
1071 + 1./2.*( derive(5)*dist(0)*dist(0) + derive(4)*dist(1)*dist(1)
1072 + derive(3)*dist(2)*dist(2) )
1073 + derive(2)*dist(0)*dist(1) + derive(1)*dist(0)*dist(2)
1074 + derive(0)*dist(1)*dist(2) ;
1075
1076 dtrans_cen = transporte(amont) + transporte(aval) - trans_c_g(0) ;
1077 //dtrans_max = std::max( transporte(amont), transporte(aval) ) ;
1078 dtrans_max = std::max( transporte(amont), dtrans_cen ) ;
1079 dtrans_min = std::min( transporte(amont), dtrans_cen ) ;
1080 //dtrans_min = std::min( dtrans_min, dtrans_cen ) ;
1081 }
1082 else
1083 {
1084 trans(j) = transporte(num3,j) + derive(8,j)*dist(0) * coor_trans(0,0)
1085 + derive(7,j)*dist(1) * coor_trans(1,1)
1086 + derive(6,j)*dist(2) * coor_trans(2,2)
1087 + 1./2.*( derive(5,j)*dist(0)*dist(0) + derive(4,j)*dist(1)*dist(1)
1088 + derive(3,j)*dist(2)*dist(2) )
1089 + derive(2,j)*dist(0)*dist(1) + derive(1,j)*dist(0)*dist(2)
1090 + derive(0,j)*dist(1)*dist(2) ;
1091
1092 dtrans_cen = transporte(amont,j) + transporte(aval,j) - trans_c_g(j) ;
1093 dtrans_max = std::max( transporte(amont,j), dtrans_cen ) ;
1094
1095 dtrans_min = std::min( transporte(amont,j), dtrans_cen ) ;
1096
1097 }
1098
1099 if (trans(j) > dtrans_max )
1100 {
1101 //Cerr << " lim max in element " << poly << finl ;
1102 trans(j) = dtrans_max ;
1103 nlim++ ;
1104 }
1105 if (trans(j) < dtrans_min )
1106 {
1107 //Cerr << " lim min in element " << poly << finl ;
1108 trans(j) = dtrans_min ;
1109 nlim++ ;
1110 }
1111
1112 }
1113 else
1114 {
1115 if(ncomp == 1)
1116 {
1117 trans(0) = transporte(num1) + transporte(num2) - trans_c_g(0) ;
1118
1119 }
1120 else
1121 {
1122 trans(j) = transporte(num1,j) + transporte(num2,j) - trans_c_g(j) ;
1123 //Cerr << " singx != 0 " << finl ;
1124 }
1125 }
1126 }
1127 //double lo_x = xv(num1,0)+xv(num2,0)-xp(poly,0) ;
1128 //double lo_y = xv(num1,1)+xv(num2,1)-xp(poly,1) ;
1129 //double lo_z = xv(num1,2)+xv(num2,2)-xp(poly,2) ;
1130
1131 //double exact0= -sin(kwave*lo_x)*cos(kwave*lo_y);
1132 //double exact1= cos(kwave*lo_x)*sin(kwave*lo_y);
1133
1134 //double exact0 = lo_x*lo_x ;
1135 //double exact1 = lo_y*lo_y ;
1136 //double exact2 = lo_z*lo_z ;
1137
1138 //double moyenne0 = ( transporte(num1,0) + transporte(num2,0) ) - trans_c_g(0) ;
1139 //double moyenne1 = ( transporte(num1,1) + transporte(num2,1) ) - trans_c_g(1) ;
1140 //double moyenne2 = ( transporte(num1,2) + transporte(num2,2) ) - trans_c_g(2) ;
1141 //double moyenne0 = transporte(num1,0) ;
1142 //double moyenne1 = transporte(num1,1) ;
1143 //double moyenne2 = transporte(num1,2) ; ;
1144
1145 //double x0 = ( xv(num1,0) + xv(num2,0) ) - xp(poly,0) ;
1146 //double x1 = ( xv(num1,1) + xv(num2,1) ) - xp(poly,1) ;
1147 //double x2 = ( xv(num1,2) + xv(num2,2) ) - xp(poly,2) ;
1148 //double x0 = xv(num1,0) ;
1149 //double x1 = xv(num1,1) ;
1150 //double x2 = xv(num1,2) ;
1151
1152 //double y0 = xv(num3,0) + dist(0) ;
1153 //double y1 = xv(num3,1) + dist(1) ;
1154 //double y2 = xv(num3,2) + dist(2) ;
1155
1156 //if (sing == 0 ) Cerr << moyenne0 << " " <<trans(0) << " " << moyenne1 <<
1157 // " " << trans(1) << " " << moyenne2 << " " << trans(2) <<finl ;
1158 //Cerr << " exact0" << " " << exact0 << " " << "exact1" << " " << exact1
1159 // << " exact2" << " " << exact2 << finl ;
1160
1161 // Cerr << x0<< " " << y0 << " "<<x1 << " "<< y1<< " " << x2 << " "<< y2 << finl ;
1162 } // ende first = 0 und normal weiter
1163
1164 if (psc >= 0)
1165 {
1166 tab_fluent[num2] += psc;
1167 }
1168 else
1169 {
1170 tab_fluent[num1] -= psc;
1171 }
1172
1173 if (ncomp == 1)
1174 {
1175 flux = trans(0)*psc;
1176 resu(num1) -= flux;
1177 resu(num2) += flux;
1178 }
1179 else
1180 {
1181 for (i=0; i<ncomp; i++)
1182 {
1183 flux = trans(i)*psc;
1184 resu(num1,i) -= flux;
1185 resu(num2,i) += flux;
1186 }
1187 }
1188
1189}
1190void Op_Conv_DI_L2_VEF_Face::poly_DI_L2_2d(int N, int M, DoubleTab& derive ,int poly, int ncomp,
1191 const DoubleTab& transporte, int& sing) const
1192{
1193 //Cerr << " in poly_DI_L2_2d " << finl ;
1194 const Domaine_VEF& domaine_VEF = le_dom_vef.valeur();
1195
1196 const IntTab& elem_faces = domaine_VEF.elem_faces();
1197 const IntTab& face_voisins = domaine_VEF.face_voisins();
1198
1199 const Domaine& domaine = domaine_VEF.domaine();
1200
1201 int nfac = domaine.nb_faces_elem();
1202
1203 const DoubleTab& xv = domaine_VEF.xv();
1204 const DoubleTab& xp = domaine_VEF.xp();
1205
1206 int i, j, face_adj , face_glob , poly1 ;
1207 DoubleVect dist(dimension) ;
1208 DoubleTab coor_trans(dimension,dimension) ;
1209 IntVect face(nfac) ;
1210
1211 DoubleTab L(M,N) ;
1212 DoubleTab B(M,ncomp) ;
1213 DoubleVect dTransp_x(N) ;
1214 DoubleVect C(N) ;
1215 DoubleVect D(N) ;
1216 DoubleVect SM(N) ;
1217 double poid ;
1218
1219 for(face_adj=0; face_adj<nfac; face_adj ++)
1220 {
1221 face_glob = elem_faces(poly, face_adj );
1222 face(face_adj) = face_glob ;
1223 }
1224
1225 int num3 = elem_faces(poly, 0);
1226
1227 //coor_trans(0) = cos( (xv(num3,1) - xp(poly,1)) / (xv(num3,0) - xp(poly,0))) ;
1228 //coor_trans(1) = cos(1.- (xv(num3,1) - xp(poly,1)) / (xv(num3,0) - xp(poly,0))) ;
1229
1230 //coor_trans(0,0) = cos(3.14159/4.) ;
1231 //coor_trans(0,1) = cos(3.14159/4.) ;
1232 //coor_trans(1,0) = cos(3.14159/2.-3.14159/4.) ;
1233 //coor_trans(1,1) = cos(3.14159/2.+3.14159/4.) ;
1234 coor_trans(0,0) = 1. ;
1235 coor_trans(0,1) = 0. ;
1236 coor_trans(1,0) = 0. ;
1237 coor_trans(1,1) = 1. ;
1238 int row = -1 ;
1239
1240 for(int face_adj_poly =0; face_adj_poly < nfac; face_adj_poly ++)
1241 {
1242 poly1 = face_voisins(face[face_adj_poly],0);
1243 if (poly1 == poly ) poly1 = face_voisins(face[face_adj_poly], 1);
1244
1245 for(face_adj=0; face_adj<nfac; face_adj ++)
1246 {
1247 face_glob = elem_faces(poly1, face_adj);
1248 if (face_glob != num3 )
1249 {
1250 row++ ;
1251 dist = 0.;
1252 for (i=0; i<dimension; i++)
1253 for (j=0; j<dimension; j++)
1254 dist(i) += (xv(face_glob,j) - xv(num3,j)) * coor_trans(i,j) ;
1255
1256 poid = 0. ;
1257 for (i=0; i<dimension; i++)
1258 poid += ( (xv(face_glob,i)-xp(poly,i)) * (xv(face_glob,i)-xp(poly,i)) );
1259 poid = 1./sqrt(poid) ;
1260
1261 L(row,4) = poid * dist(0) * coor_trans(0,0) ;
1262 L(row,3) = poid * dist(1) * coor_trans(1,1) ;
1263 L(row,2) = poid * 1./2.*dist(0)*dist(0) ;
1264 L(row,1) = poid * 1./2.*dist(1)*dist(1) ;
1265 L(row,0) = poid * dist(0)*dist(1) ;
1266
1267 if (ncomp == 1)
1268 B(row,0) = poid * (transporte(face_glob) - transporte(num3) ) ;
1269 else for (j=0; j<ncomp; j++)
1270 B(row,j) = poid * (transporte(face_glob,j) - transporte(num3,j)) ;
1271
1272 }
1273 }
1274
1275 }
1276
1277 DoubleTab Lij(N,N) ;
1278 DoubleTab Bij(N,ncomp) ;
1279
1280 for (j=0; j<N; j++)
1281 for (i=0; i<ncomp; i++)
1282 for (int k=0; k<M; k++) Bij(j,i) += L(k,j) * B(k,i) ;
1283
1284 for (i=0; i<N; i++)
1285 for (j=0; j<N; j++)
1286 for (int k=0; k<M; k++) Lij(i,j) += L(k,i) * L(k,j) ;
1287
1288
1289 for (j=0; j<ncomp; j++)
1290 {
1291 for (i=0; i<N; i++) SM(i)=Bij(i,j) ;
1292
1293 qrsolv(Lij, N, SM, dTransp_x, sing, j, C, D);
1294
1295 if (sing == 0 )
1296 {
1297 if(ncomp == 1)
1298 {
1299 for (int val_N= 0; val_N < N ; val_N++ )
1300 derive(val_N) = dTransp_x(val_N) ;
1301 }
1302 else
1303 {
1304 for (int val_N= 0; val_N < N ; val_N++ )
1305 derive(val_N,j) = dTransp_x(val_N) ;
1306 }
1307 }
1308 }
1309}
1310
1311void Op_Conv_DI_L2_VEF_Face::poly_DI_L2_3d(int N, int M, DoubleTab& derive ,int poly, int ncomp,
1312 const DoubleTab& transporte ,int& sing) const
1313{
1314 // Cerr<<"Op_Conv_DI_L2_VEF_Face::poly_DI_L2_3d 0 "<<finl;
1315 const Domaine_VEF& domaine_VEF = le_dom_vef.valeur();
1316
1317 const IntTab& elem_faces = domaine_VEF.elem_faces();
1318 const IntTab& face_voisins = domaine_VEF.face_voisins();
1319
1320 const Domaine& domaine = domaine_VEF.domaine();
1321
1322 int nfac = domaine.nb_faces_elem();
1323
1324 const DoubleTab& xv = domaine_VEF.xv();
1325 const DoubleTab& xp = domaine_VEF.xp();
1326
1327 int i, j, face_adj , face_glob , poly1 ;
1328
1329 DoubleVect dist(dimension) ;
1330 DoubleTab coor_trans(dimension,dimension) ;
1331 IntVect face(nfac) ;
1332
1333 DoubleTab L(M,N) ;
1334 DoubleTab B(M,ncomp) ;
1335 DoubleVect dTransp_x(N) ;
1336 DoubleVect C(N) ;
1337 DoubleVect D(N) ;
1338 DoubleVect SM(N) ;
1339 double poid ;
1340 // Cerr<<"Op_Conv_DI_L2_VEF_Face::poly_DI_L2_3d 1 "<<finl;
1341
1342 for(face_adj=0; face_adj<nfac; face_adj ++)
1343 {
1344 face_glob = elem_faces(poly, face_adj );
1345 face(face_adj) = face_glob ;
1346 }
1347 // Cerr<<"Op_Conv_DI_L2_VEF_Face::poly_DI_L2_3d 2 "<<finl;
1348
1349 int num3 = elem_faces(poly, 0);
1350
1351 coor_trans(0,0) = cos(3.14159/4.) ;
1352 coor_trans(0,1) = cos(3.14159/4.) ;
1353 coor_trans(0,2) = cos(3.14159/2.) ;
1354
1355 coor_trans(1,0) = cos(3.14159/2.-3.14159/4.) ;
1356 coor_trans(1,1) = cos(3.14159/2.+3.14159/4.) ;
1357 coor_trans(1,2) = cos(3.14159/2.) ;
1358
1359 coor_trans(2,0) = cos(3.14159/2.) ;
1360 coor_trans(2,1) = cos(3.14159/2.) ;
1361 coor_trans(2,2) = cos(0.) ;
1362
1363 //coor_trans = 0.;
1364 //coor_trans(0,0) = 1. ;
1365 //coor_trans(1,1) = 1. ;
1366 //coor_trans(2,2) = 1. ;
1367
1368 int row = -1 ;
1369 // Cerr<<"Op_Conv_DI_L2_VEF_Face::poly_DI_L2_3d 3 "<<finl;
1370
1371 for(int face_adj_poly =0; face_adj_poly < nfac; face_adj_poly ++)
1372 {
1373 poly1 = face_voisins(face[face_adj_poly],0);
1374 // Cerr<<"avant poly = "<<poly<<", poly1 = "<<poly1<<", proc = "<<me()<<finl;
1375 // Cerr<<"face[face_adj_poly] = "<<face[face_adj_poly]<<finl;
1376 if (poly1 == poly ) poly1 = face_voisins(face[face_adj_poly], 1);
1377 // Cerr<<"apres poly = "<<poly<<", poly1 = "<<poly1<<", proc = "<<me()<<finl;
1378 for(face_adj=0; face_adj<nfac; face_adj ++)
1379 {
1380 face_glob = elem_faces(poly1, face_adj);
1381 if (face_glob != num3 )
1382 {
1383 row++ ;
1384 dist=0 ;
1385 for (i=0; i<dimension; i++)
1386 for (j=0; j<dimension; j++)
1387 dist(i) += (xv(face_glob,j) - xv(num3,j)) * coor_trans(i,j) ;
1388
1389 poid = 0. ;
1390 for (i=0; i<dimension; i++)
1391 poid += ( (xv(face_glob,i)-xp(poly,i)) * (xv(face_glob,i)-xp(poly,i)) );
1392 poid = 1./(poid) ;
1393
1394
1395 L(row,8) = poid*dist(0) * coor_trans(0,0) ;
1396 L(row,7) = poid*dist(1) * coor_trans(1,1) ;
1397 L(row,6) = poid*dist(2) * coor_trans(2,2) ;
1398 L(row,5) = poid*1./2.*dist(0)*dist(0) ;
1399 L(row,4) = poid*1./2.*dist(1)*dist(1) ;
1400 L(row,3) = poid*1./2.*dist(2)*dist(2) ;
1401 L(row,2) = poid*dist(0)*dist(1) ;
1402 L(row,1) = poid*dist(0)*dist(2) ;
1403 L(row,0) = poid*dist(1)*dist(2) ;
1404
1405 if (ncomp == 1)
1406 B(row,0) = poid * (transporte(face_glob) - transporte(num3) ) ;
1407 else for (j=0; j<ncomp; j++)
1408 B(row,j) = poid * (transporte(face_glob,j) - transporte(num3,j)) ;
1409 }
1410 }
1411 }
1412
1413 // Cerr<<"Op_Conv_DI_L2_VEF_Face::poly_DI_L2_3d 4 "<<finl;
1414
1415 DoubleTab Lij(N,N) ;
1416 DoubleTab Bij(N,ncomp) ;
1417
1418 for (j=0; j<N; j++)
1419 for (i=0; i<ncomp; i++)
1420 for (int k=0; k<M; k++) Bij(j,i) += L(k,j) * B(k,i) ;
1421
1422 // Cerr<<"Op_Conv_DI_L2_VEF_Face::poly_DI_L2_3d 5 "<<finl;
1423
1424 for (i=0; i<N; i++)
1425 for (j=0; j<N; j++)
1426 for (int k=0; k<M; k++) Lij(i,j) += L(k,i) * L(k,j) ;
1427 // Cerr<<"Op_Conv_DI_L2_VEF_Face::poly_DI_L2_3d 6 "<<finl;
1428
1429
1430 for (j=0; j<ncomp; j++)
1431 {
1432 for (i=0; i<N; i++) SM(i)=Bij(i,j) ;
1433
1434 qrsolv(Lij, N, SM, dTransp_x, sing, j, C, D);
1435
1436 if (sing == 0 )
1437 {
1438 if(ncomp == 1)
1439 {
1440 for (int val_N= 0; val_N < N ; val_N++ )
1441 derive(val_N) = dTransp_x(val_N) ;
1442 }
1443 else
1444 {
1445 for (int val_N= 0; val_N < N ; val_N++ )
1446 derive(val_N,j) = dTransp_x(val_N) ;
1447 }
1448 }
1449 }
1450
1451 // Cerr<<"Op_Conv_DI_L2_VEF_Face::poly_DI_L2_3d 7 "<<finl;
1452
1453}
1454
Champ_Inc_base
Classe Champ_Inc_base.
Definition Champ_Inc_base.h:53
Champ_Inc_base::valeurs
DoubleTab & valeurs() override
Renvoie le tableau des valeurs du champ au temps courant.
Definition Champ_Inc_base.h:77
Champ_P1NC
Definition Champ_P1NC.h:35
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
Domaine_32_64::nb_faces_elem
int nb_faces_elem(int=0) const
Renvoie le nombre de face de type i des elements geometriques constituants le domaine.
Definition Domaine.h:484
Domaine_Cl_VEF
Definition Domaine_Cl_VEF.h:40
Domaine_Cl_VEF::type_elem_Cl
int type_elem_Cl(int i) const
Definition Domaine_Cl_VEF.h:62
Domaine_Cl_VEF::normales_facettes_Cl
DoubleTab & normales_facettes_Cl()
Definition Domaine_Cl_VEF.h:56
Domaine_Cl_dis_base::les_conditions_limites
const Cond_lim & les_conditions_limites(int) const
Renvoie la i-ieme condition aux limites.
Definition Domaine_Cl_dis_base.cpp:386
Domaine_VEF
class Domaine_VEF
Definition Domaine_VEF.h:54
Domaine_VEF::rang_elem_non_std
IntVect & rang_elem_non_std()
Definition Domaine_VEF.h:86
Domaine_VEF::type_elem
const Elem_VEF_base & type_elem() const
Definition Domaine_VEF.h:75
Domaine_VEF::facette_normales
auto & facette_normales()
Definition Domaine_VEF.h:84
Domaine_VF::face_normales
virtual double face_normales(int face, int comp) const
Definition Domaine_VF.h:47
Domaine_VF::xv
double xv(int num_face, int k) const
Definition Domaine_VF.h:76
Domaine_VF::elem_faces
int elem_faces(int i, int j) const
renvoie le numero de le ieme face de la maille num_elem la facon dont ces faces sont numerotees est
Definition Domaine_VF.h:543
Domaine_VF::xp
double xp(int num_elem, int k) const
Definition Domaine_VF.h:77
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
Domaine_dis_base::nb_elem_tot
int nb_elem_tot() const
Definition Domaine_dis_base.h:50
Domaine_dis_base::nb_front_Cl
int nb_front_Cl() const
Definition Domaine_dis_base.h:53
Domaine_dis_base::domaine
const Domaine & domaine() const
Definition Domaine_dis_base.h:46
Elem_VEF_base::calcul_vc
virtual void calcul_vc(const ArrOfInt &, ArrOfDouble &, const ArrOfDouble &, const DoubleTab &, const Champ_Inc_base &, int, const DoubleVect &) const =0
Elem_VEF_base::KEL
const IntTab & KEL() const
Definition Elem_VEF_base.h:31
Elem_VEF_base::nb_facette
virtual int nb_facette() const =0
Entree
Class defining operators and methods for all reading operation in an input flow (file,...
Definition Entree.h:42
Equation_base::milieu
virtual const Milieu_base & milieu() const =0
Equation_base::schema_temps
Schema_Temps_base & schema_temps()
Renvoie le schema en temps associe a l'equation.
Definition Equation_base.cpp:865
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
Milieu_base::porosite_face
DoubleVect & porosite_face()
Definition Milieu_base.h:62
MorEqn::equation
const Equation_base & equation() const
Renvoie la reference sur l'equation pointe par MorEqn::mon_equation.
Definition MorEqn.h:62
Neumann_sortie_libre
classe Neumann_sortie_libre Cette classe represente une frontiere ouverte sans vitesse imposee
Definition Neumann_sortie_libre.h:34
Neumann_sortie_libre::val_ext
double val_ext(int i) const override
Renvoie la valeur de la i-eme composante du champ impose a l'exterieur de la frontiere.
Definition Neumann_sortie_libre.cpp:120
Objet_U::dimension
static int dimension
Definition Objet_U.h:99
Objet_U::que_suis_je
const Nom & que_suis_je() const
renvoie la chaine identifiant la classe.
Definition Objet_U.cpp:104
Objet_U::readOn
virtual Entree & readOn(Entree &)
Lecture d'un Objet_U sur un flot d'entree Methode a surcharger.
Definition Objet_U.cpp:293
Objet_U::printOn
virtual Sortie & printOn(Sortie &) const
Ecriture de l'objet sur un flot de sortie Methode a surcharger.
Definition Objet_U.cpp:282
Op_Conv_DI_L2_VEF_Face
class Op_Conv_DI_L2_VEF_Face
Definition Op_Conv_DI_L2_VEF_Face.h:34
Op_Conv_DI_L2_VEF_Face::poly_DI_L2_3d
void poly_DI_L2_3d(int, int, DoubleTab &, int, int, const DoubleTab &, int &) const
Definition Op_Conv_DI_L2_VEF_Face.cpp:1311
Op_Conv_DI_L2_VEF_Face::associer_vitesse
void associer_vitesse(const Champ_base &) override
Definition Op_Conv_DI_L2_VEF_Face.cpp:46
Op_Conv_DI_L2_VEF_Face::ajouter
DoubleTab & ajouter(const DoubleTab &, DoubleTab &) const override
Definition Op_Conv_DI_L2_VEF_Face.cpp:403
Op_Conv_DI_L2_VEF_Face::poly_DI_L2_2d
void poly_DI_L2_2d(int, int, DoubleTab &, int, int, const DoubleTab &, int &) const
Definition Op_Conv_DI_L2_VEF_Face.cpp:1190
Op_Conv_DI_L2_VEF_Face::reconst_DI_L2_3d
void reconst_DI_L2_3d(DoubleTab &, int, double, int, int, const DoubleTab &, int, DoubleTab &, DoubleVect &, int, int, DoubleVect &, int &) const
Definition Op_Conv_DI_L2_VEF_Face.cpp:945
Op_Conv_DI_L2_VEF_Face::reconst_DI_L2_2d
void reconst_DI_L2_2d(DoubleTab &, int, double, int, int, const DoubleTab &, int, DoubleTab &, DoubleVect &, int, int &) const
Definition Op_Conv_DI_L2_VEF_Face.cpp:754
Op_Conv_VEF_base
class Op_Conv_VEF_base
Definition Op_Conv_VEF_base.h:37
Op_Conv_VEF_base::fluent_
DoubleVect fluent_
Definition Op_Conv_VEF_base.h:61
Op_Conv_VEF_base::vitesse
const Champ_Inc_base & vitesse() const
Definition Op_Conv_VEF_base.h:41
Op_VEF_Face::modifier_flux
void modifier_flux(const Operateur_base &) const
Definition Op_VEF_Face.cpp:338
Operateur_base::flux_bords_
DoubleTab flux_bords_
Definition Operateur_base.h:140
Periodique
classe Periodique Cette classe represente une condition aux limites periodique.
Definition Periodique.h:31
Periodique::face_associee
int face_associee(int i) const
Definition Periodique.h:35
Process::is_parallel
static bool is_parallel()
Definition Process.cpp:110
Process::exit
static void exit(int exit_code=-1)
Routine de sortie de TRUST dans une region Kokkos.
Definition Process.cpp:455
Process::is_sequential
static bool is_sequential()
Definition Process.cpp:115
Schema_Temps_base::pas_de_temps
double pas_de_temps() const
Renvoie le pas de temps (delta_t) courant.
Definition Schema_Temps_base.h:393
Sortie
Classe de base des flux de sortie.
Definition Sortie.h:52
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
TRUSTVect::ajoute_sans_ech_esp_virt
void ajoute_sans_ech_esp_virt(_SCALAR_TYPE_ alpha, const TRUSTVect &y, Mp_vect_options opt=VECT_REAL_ITEMS)
Definition TRUSTVect.tpp:450