TrioCFD 1.9.9_beta
TrioCFD documentation
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// Functions of class 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 //This procedure corresponds to the Gaussian elimination method
56 test_flora = 0; //test to determine whether the matrix is invertible (0) or not (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<<"Error flora: non-invertible matrix at index "<<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<<"Error flora: non-invertible matrix at index"<<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 //This procedure corresponds to the Gaussian elimination method
104 test_flora = 1;//test to determine whether the matrix is invertible (1) or not (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<<"Error flora: non-invertible matrix at index "<<m<<finl;
133 test_flora = 1;
134 }
135 else //swap columns m and 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<<"Error flora: non-invertible matrix at index"<<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 //Builds the QR decomposition of A. The upper triangle R is stored in the upper triangle of A,
176 //except for the diagonal elements, which are stored in D. The orthogonal matrix Q is represented
177 //as a product of N-1 matrices Q(1), ..., Q(N-1) where Q(j) = Id-(uj*ujt)/cj. The i-th component of uj is 0
178 //for i=1, ...,j-1 and A(i,j) for i=j, ..., N. sing returns 0 if the decomposition is possible and 1 otherwise.
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 //Solves the system Rx=B, where R is the upper triangular matrix stored in A and D, from qrdcmp
226 //the result is stored in 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 //Solves the linear system 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<<"error"<<finl;
257}
258
259
260
261
262//biconjugate gradient method
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() is not parallel" << 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 computes -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 // For the convection treatment, standard polyhedra (which do not "see"
439 // any boundary conditions) are distinguished from non-standard polyhedra
440 // that have at least one face on the boundary.
441 // A standard polyhedron has n facets on which the convection scheme is applied.
442 // For a non-standard polyhedron with Dirichlet boundary conditions,
443 // part of the facets are carried by boundary faces.
444 // In short, the convection treatment for a polyhedron depends
445 // on the type (triangle, tetrahedron, ...) and the number of Dirichlet faces.
446
447 double psc;
448 int poly,face_adj,fa7,i,j,n_bord;
449 int num_face, rang ,itypcl;
450 int num10,num20;
451 int ncomp_ch_transporte, first;
452 if (transporte.nb_dim() == 1)
453 ncomp_ch_transporte=1;
454 else
455 ncomp_ch_transporte= transporte.dimension(1);
456
457 IntVect face(nfac);
458 DoubleVect vs(dimension);
459 DoubleVect vc(dimension);
460 DoubleTab vsom(nsom,dimension);
461 DoubleVect cc(dimension);
462 DoubleTab derive(1,1) ;
463 int N ,M , sing , cal_amont;
464 DoubleVect trans(ncomp_ch_transporte) ;
465 if (dimension == 2)
466 {
467 N = 5 ;
468 M = 8 ;
469 derive.resize(N,ncomp_ch_transporte) ;
470 }
471 else // (dimension == 3)
472 {
473 N = 9 ;
474 M = 15 ;
475 derive.resize(N,ncomp_ch_transporte) ;
476 }
477 // Reset the array used for stability time step computation
478 fluent_ = 0;
479
480 // Special treatment for periodic faces
481
482 int nb_faces_perio = 0;
483 // Loop to count the number of periodic faces
484 for (n_bord=0; n_bord<domaine_VEF.nb_front_Cl(); n_bord++)
485 {
486 const Cond_lim& la_cl = domaine_Cl_VEF.les_conditions_limites(n_bord);
487 if (sub_type(Periodique,la_cl.valeur()))
488 {
489 const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
490 nb_faces_perio += le_bord.nb_faces();
491 }
492 }
493
494 DoubleTab tab;
495 if (ncomp_ch_transporte == 1)
496 tab.resize(nb_faces_perio);
497 else
498 tab.resize(nb_faces_perio,ncomp_ch_transporte);
499
500 // Loop to fill tab
501 nb_faces_perio=0;
502 for (n_bord=0; n_bord<domaine_VEF.nb_front_Cl(); n_bord++)
503 {
504 const Cond_lim& la_cl = domaine_Cl_VEF.les_conditions_limites(n_bord);
505 if (sub_type(Periodique,la_cl.valeur()))
506 {
507 // const Periodique& la_cl_perio = ref_cast(Periodique, la_cl.valeur());
508 const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
509 int num1 = le_bord.num_premiere_face();
510 int num2 = num1 + le_bord.nb_faces();
511 for (num_face=num1; num_face<num2; num_face++)
512 {
513 if (ncomp_ch_transporte == 1)
514 tab(nb_faces_perio) = resu(num_face);
515 else
516 for (int comp=0; comp<ncomp_ch_transporte; comp++)
517 tab(nb_faces_perio,comp) = resu(num_face,comp);
518 nb_faces_perio++;
519 }
520 }
521 }
522
523
524 // Non-standard polyhedra are stored in 2 groups in the Domaine_VEF:
525 // - boundary and shared polyhedra
526 // - boundary and non-shared polyhedra
527 // Polyhedra are processed in the order they appear in the domain
528
529 // loop over polyhedra
530
531 int nlim = -1 ;
532 const IntTab& KEL=domaine_VEF.type_elem().KEL();
533 for (poly=0; poly<nb_elem_tot; poly++)
534 {
535
536 rang = rang_elem_non_std(poly);
537 if (rang==-1)
538 itypcl=0;
539 else
540 itypcl=domaine_Cl_VEF.type_elem_Cl(rang);
541
542 // compute face indices of the polyhedron
543 for (face_adj=0; face_adj<nfac; face_adj++)
544 face[face_adj]= elem_faces(poly,face_adj);
545
546 // compute velocity at the vertices of the polyhedra
547 for (j=0; j<dimension; j++)
548 {
549 vs[j] = la_vitesse.valeurs()(face[0],j);
550 for (i=1; i<nfac; i++)
551 vs[j]+= la_vitesse.valeurs()(face[i],j);
552 }
553 for (i=0; i<nsom; i++)
554 for (j=0; j<dimension; j++)
555 vsom(i,j) = vs[j] - dimension*la_vitesse.valeurs()(face[i],j);
556
557 // compute vc
558 domaine_VEF.type_elem().calcul_vc(face,vc,vs,vsom,vitesse(),
559 itypcl,porosite_face);
560
561
562 /* Cerr<<"premiere_face_int = "<<premiere_face_int<<", domaine_VEF.nb_faces_joint() = "<<domaine_VEF.nb_faces_joint()<<finl;
563 Cerr<<"premiere_face_std() = "<<domaine_VEF.premiere_face_std()<<finl;
564 Cerr<<"nb_elem() = "<<domaine_VEF.nb_elem()<<", nb_elem_tot() = "<<domaine_VEF.nb_elem_tot()<<finl;
565 Cerr<<"nb_faces() = "<<domaine_VEF.nb_faces()<<", nb_faces_tot() = "<<domaine_VEF.nb_faces_tot()<<finl;
566 Cerr<<"nb_faces_bord() = "<<domaine_VEF.nb_faces_bord()<<", premiere_face_bord() = "<<domaine_VEF.premiere_face_bord()<<finl;*/
567
568 cal_amont = 0 ;
569
570 int elem0,elem1,face_adj_glob;
571 for (face_adj=0; face_adj<nfac; face_adj++)
572 {
573 face_adj_glob = face[face_adj];
574 elem0 = face_voisins(face_adj_glob,0);
575 elem1 = face_voisins(face_adj_glob,1);
576 if ((elem0 == -1) || (elem1 == -1))
577 cal_amont++ ;
578 }
579
580 // compute reconstruction polynomial
581 if(dimension == 2 )
582 {
583 if ( cal_amont == 0) poly_DI_L2_2d(N,M,derive,poly,ncomp_ch_transporte,transporte,sing);
584 }
585 else
586 {
587 if ( cal_amont == 0) poly_DI_L2_3d(N,M,derive,poly,ncomp_ch_transporte,transporte,sing);
588 }
589
590 // Loop over facets of the polyhedron
591
592 for (fa7=0; fa7<nfa7; fa7++)
593 {
594 if (rang==-1)
595 for (i=0; i<dimension; i++)
596 cc[i] = facette_normales(poly,fa7,i);
597 else
598 for (i=0; i<dimension; i++)
599 cc[i] = normales_facettes_Cl(rang,fa7,i);
600
601 // Apply the convection scheme to each vertex of the facet
602
603 // reconstruction seulement wenn first = 0
604
605 first = -1 ;
606
607 // Process the vertex/vertices that are also vertices of the polyhedron
608 for (i=0; i<nb_som_facette-1; i++)
609 {
610 first++ ;
611
612 psc =0;
613 for (j=0; j<dimension; j++)
614 psc+= (vc(j)/double(nb_som_facette-1)+vsom(KEL(i+2,fa7),j))*cc[j];
615 psc /= nb_som_facette;
616
617 num10 = face[KEL(0,fa7)];
618 num20 = face[KEL(1,fa7)];
619
620 if ( cal_amont > 0)
621 {
622 convbis(psc,num10,num20,transporte,ncomp_ch_transporte,resu,fluent_);
623 }
624 else
625 {
626 if(dimension == 2 )
627 {
628 reconst_DI_L2_2d(derive,poly,psc,num10,num20,transporte,ncomp_ch_transporte,resu,fluent_,sing,
629 nlim );
630 }
631 else
632 {
633 reconst_DI_L2_3d(derive,poly,psc,num10,num20,transporte,ncomp_ch_transporte,resu,fluent_,sing,
634 first,trans,nlim);
635 }
636 }
637 }
638
639 }
640
641 } // end of loop
642
643 Cerr << " limitiert in " << nlim << " valeurs " << finl ;
644
645 int voisine;
646 nb_faces_perio = 0;
647 double diff1,diff2;
648
649 // Size the array of convective fluxes at the domain boundary
650 DoubleTab& flux_b = flux_bords_;
651 flux_b.resize(domaine_VEF.nb_faces_bord(),ncomp_ch_transporte);
652 flux_b = 0.;
653
654 // Loop over boundaries to process boundary conditions
655 // a convective term is accounted for only for
656 // Neumann_sortie_libre boundary conditions
657
658 for (n_bord=0; n_bord<domaine_VEF.nb_front_Cl(); n_bord++)
659 {
660
661 const Cond_lim& la_cl = domaine_Cl_VEF.les_conditions_limites(n_bord);
662
663 if (sub_type(Neumann_sortie_libre,la_cl.valeur()))
664 {
665 const Neumann_sortie_libre& la_sortie_libre = ref_cast(Neumann_sortie_libre, la_cl.valeur());
666 const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
667 int num1 = le_bord.num_premiere_face();
668 int num2 = num1 + le_bord.nb_faces();
669 for (num_face=num1; num_face<num2; num_face++)
670 {
671 psc =0;
672 for (i=0; i<dimension; i++)
673 psc += la_vitesse.valeurs()(num_face,i)*face_normales(num_face,i);
674 if (psc>0)
675 if (ncomp_ch_transporte == 1)
676 {
677 resu(num_face) -= psc*transporte(num_face);
678 flux_b(num_face,0) -= psc*transporte(num_face);
679 }
680 else
681 for (i=0; i<ncomp_ch_transporte; i++)
682 {
683 resu(num_face,i) -= psc*transporte(num_face,i);
684 flux_b(num_face,i) -= psc*transporte(num_face,i);
685 }
686 else
687 {
688 if (ncomp_ch_transporte == 1)
689 {
690 resu(num_face) -= psc*la_sortie_libre.val_ext(num_face-num1);
691 flux_b(num_face,0) -= psc*la_sortie_libre.val_ext(num_face-num1);
692 }
693 else
694 for (i=0; i<ncomp_ch_transporte; i++)
695 {
696 resu(num_face,i) -= psc*la_sortie_libre.val_ext(num_face-num1,i);
697 flux_b(num_face,i) -= psc*la_sortie_libre.val_ext(num_face-num1);
698 }
699 fluent_[num_face] -= psc;
700 }
701 }
702 }
703 else if (sub_type(Periodique,la_cl.valeur()))
704 {
705 const Periodique& la_cl_perio = ref_cast(Periodique, la_cl.valeur());
706 const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
707 int num1 = le_bord.num_premiere_face();
708 int num2 = num1 + le_bord.nb_faces();
709 IntVect fait(le_bord.nb_faces());
710 fait = 0;
711 for (num_face=num1; num_face<num2; num_face++)
712 {
713 if (fait[num_face-num1] == 0)
714 {
715 voisine = la_cl_perio.face_associee(num_face-num1) + num1;
716
717 if (ncomp_ch_transporte == 1)
718 {
719 diff1 = resu(num_face)-tab(nb_faces_perio);
720 diff2 = resu(voisine)-tab(nb_faces_perio+voisine-num_face);
721 resu(voisine) += diff1;
722 resu(num_face) += diff2;
723 flux_b(voisine,1) += diff1;
724 flux_b(num_face,0) += diff2;
725 }
726 else
727 for (int comp=0; comp<ncomp_ch_transporte; comp++)
728 {
729 diff1 = resu(num_face,comp)-tab(nb_faces_perio,comp);
730 diff2 = resu(voisine,comp)-tab(nb_faces_perio+voisine-num_face,comp);
731 resu(voisine,comp) += diff1;
732 resu(num_face,comp) += diff2;
733 flux_b(voisine,comp) += diff1;
734 flux_b(num_face,comp) += diff2;
735 }
736
737 fait[num_face-num1]= 1;
738 fait[voisine-num1] = 1;
739 }
740 nb_faces_perio++;
741 }
742 }
743 }
744 modifier_flux(*this);
745 return resu;
746
747}
748
749
750void Op_Conv_DI_L2_VEF_Face::reconst_DI_L2_2d(DoubleTab& derive ,int poly,
751 double psc,int num1,int num2,
752 const DoubleTab& transporte,
753 int ncomp,
754 DoubleTab& resu,
755 DoubleVect& tab_fluent, int sing, int& nlim) const
756{
757
758 //Cerr << " in reconst_DI_L2_2d " << sing << finl ;
759
760 const Domaine_VEF& domaine_VEF = le_dom_vef.valeur();
761
762 const IntTab& elem_faces = domaine_VEF.elem_faces();
763 // const IntTab& face_voisins = domaine_VEF.face_voisins();
764
765 const Domaine& domaine = domaine_VEF.domaine();
766
767 int nfac = domaine.nb_faces_elem();
768
769 const DoubleTab& xv = domaine_VEF.xv();
770 const DoubleTab& xp = domaine_VEF.xp();
771
772 double dt=equation().schema_temps().pas_de_temps();
773 const Champ_Inc_base& la_vitesse=vitesse_.valeur();
774
775 int i, j, face_adj , face_glob, numg=-1 ;
776 DoubleVect dist(dimension) ;
777 DoubleVect trans_c_g(ncomp) ;
778 DoubleVect vs(ncomp) ;
779 DoubleVect trans(ncomp) ;
780 DoubleTab coor_trans(dimension,dimension) ;
781
782 trans_c_g = 0. ;
783
784 for (j=0; j<ncomp; j++)
785 {
786 for(face_adj=0; face_adj<nfac; face_adj ++)
787 {
788 face_glob = elem_faces(poly, face_adj);
789
790 if (ncomp == 1) trans_c_g(0) += transporte(face_glob)/double(nfac) ;
791 else trans_c_g(j) += transporte(face_glob,j)/double(nfac) ;
792
793 if ((face_glob != num1) && (face_glob != num2)) numg = face_glob ;
794
795 }
796 }
797 int num3 = elem_faces(poly, 0 ) ;
798
799 //coor_trans(0,0) = cos(3.14159/4.) ;
800 //coor_trans(0,1) = cos(3.14159/4.) ;
801 //coor_trans(1,0) = cos(3.14159/2.-3.14159/4.) ;
802 //coor_trans(1,1) = cos(3.14159/2.+3.14159/4.) ;
803 coor_trans(0,0) = 1. ;
804 coor_trans(0,1) = 0. ;
805 coor_trans(1,0) = 0. ;
806 coor_trans(1,1) = 1. ;
807
808 for (i=0; i<dimension; i++) vs(i)= (la_vitesse.valeurs()(num1,i)+la_vitesse.valeurs()(num2,i)) / 2.;
809
810 dist = 0. ;
811
812 for (i=0; i<dimension; i++)
813 for (j=0; j<dimension; j++)
814 dist(i) += (2.*xp(poly,j) - xv(numg,j) - xv(num3,j) - vs(j)*dt/2. ) * coor_trans(i,j) ;
815
816 double dtrans_max, dtrans_min, dtrans_cen ;
817
818 int amont, aval ;
819
820 if (psc >= 0)
821 {
822 amont = num1;
823 aval = num2;
824 }
825 else
826 {
827 amont = num2;
828 aval = num1;
829 }
830
831 for (j=0; j<ncomp; j++)
832 {
833
834 if (sing == 0 )
835 {
836 if(ncomp == 1)
837 {
838 trans(0) = transporte(num3) + derive(4)*dist(0) * coor_trans(0,0)
839 + derive(3)*dist(1) * coor_trans(1,1)
840 + 1./2.*( derive(2)*dist(0)*dist(0) + derive(1)*dist(1)*dist(1) )
841 + derive(0)*dist(0)*dist(1) ;
842
843 dtrans_cen = transporte(amont) + transporte(aval) - trans_c_g(0) ;
844 dtrans_max = std::max( transporte(amont), dtrans_cen ) ;
845 dtrans_min = std::min( transporte(amont), dtrans_cen ) ;
846 }
847 else
848 {
849 trans(j) = transporte(num3,j) + derive(4,j)*dist(0) * coor_trans(0,0)
850 + derive(3,j)*dist(1) * coor_trans(1,1)
851 + 1./2.*( derive(2,j)*dist(0)*dist(0) + derive(1,j)*dist(1)*dist(1) )
852 + derive(0,j)*dist(0)*dist(1) ;
853
854
855 //dtrans_cen = transporte(num1,j) + transporte(num2,j) - trans_c_g(j) ;
856 //dtrans_cen = 1./2.*(transporte(amont,j) + transporte(aval,j)) ;
857 dtrans_cen = transporte(num1,j) + transporte(num2,j) - trans_c_g(j) ;
858
859 //dtrans_max = std::max( transporte(amont,j), transporte(aval,j)) ;
860 dtrans_max = std::max( transporte(amont,j) , dtrans_cen) ;
861 dtrans_min = std::min( transporte(amont,j) , dtrans_cen) ;
862 //dtrans_min = std::min( dtrans_min, dtrans_cen) ;
863
864 }
865
866 if (trans(j) > dtrans_max )
867 {
868 //Cerr << "lim_max "<< finl ;
869 nlim++ ;
870 trans(j) = dtrans_max ;
871 }
872 if (trans(j) < dtrans_min )
873 {
874 //Cerr << "lim_min "<< finl ;
875 nlim++ ;
876 trans(j) = dtrans_min ;
877 }
878 }
879 else
880 {
881 if(ncomp == 1) trans(0) = ( transporte(num1) + transporte(num2) ) - trans_c_g(0) ;
882 else trans(j) = ( transporte(num1,j) + transporte(num2,j) ) - trans_c_g(j) ;
883 Cerr << " singx != 0 " << finl ;
884 }
885 }
886
887 // double kwave = 1. ;
888 // double lo_x = xv(num1,0)+xv(num2,0)-xp(poly,0) ;
889 // double lo_y = xv(num1,1)+xv(num2,1)-xp(poly,1) ;
890
891 // double exact0= -sin(kwave*lo_x)*cos(kwave*lo_y);
892 // double exact1= cos(kwave*lo_x)*sin(kwave*lo_y);
893 // double exact0 = kwave*lo_x*lo_x ;
894 // double exact1 = kwave*lo_y*lo_y ;
895
896 // double moyenne0 = ( transporte(num1,0) + transporte(num2,0) ) - trans_c_g(0) ;
897 // double moyenne1 = ( transporte(num1,1) + transporte(num2,1) ) - trans_c_g(1) ;
898
899 // double xv1 = (xv(num1,0)+xv(num2,0)-xp(poly,0)) - (xv(num3,0) + dist(0) - coor_trans(0) ) ;
900 // double xv2 = (xv(num1,1)+xv(num2,1)-xp(poly,1)) - (xv(num3,1) + dist(1) - coor_trans(1) ) ;
901
902 // double moyenne0 = transporte(numg,0) ;
903 // double moyenne1 = transporte(numg,1) ;
904
905 // if (sing == 0 ) Cerr << exact0 << " " << moyenne0 << " " << trans(0) << " "
906 // << exact1 << " " << moyenne1 << " " << trans(1) << finl ;
907
908 // if (sing == 0 ) Cerr << moyenne0 << " " << trans(0) << " "
909 // << moyenne1 << " " << trans(1) << finl ;
910 // Cerr << numg << " " << num3 << finl ;
911
912 double flux;
913
914 if (psc >= 0)
915 {
916 tab_fluent[num2] += psc;
917 }
918 else
919 {
920 tab_fluent[num1] -= psc;
921 }
922
923 if (ncomp == 1)
924 {
925 flux = trans(0)*psc;
926 resu(num1) -= flux;
927 resu(num2) += flux;
928 }
929 else
930 {
931 for (i=0; i<ncomp; i++)
932 {
933 flux = trans(i)*psc;
934 resu(num1,i) -= flux;
935 resu(num2,i) += flux;
936 }
937 }
938
939}
940
941void Op_Conv_DI_L2_VEF_Face::reconst_DI_L2_3d(DoubleTab& derive, int poly,
942 double psc,int num1,int num2,
943 const DoubleTab& transporte,
944 int ncomp,
945 DoubleTab& resu,
946 DoubleVect& tab_fluent ,
947 int sing, int first, DoubleVect& trans,
948 int& nlim) const
949{
950
951 const Domaine_VEF& domaine_VEF = le_dom_vef.valeur();
952
953 const IntTab& elem_faces = domaine_VEF.elem_faces();
954 // const IntTab& face_voisins = domaine_VEF.face_voisins();
955
956 // const int nb_faces = domaine_VEF.nb_faces();
957 // const int nfa7 = domaine_VEF.type_elem().nb_facette();
958 // const int nb_elem = domaine_VEF.nb_elem();
959 // const int nb_elem_tot = domaine_VEF.nb_elem_tot();
960 // int premiere_face_std = domaine_VEF.premiere_face_std();
961
962 const Domaine& domaine = domaine_VEF.domaine();
963
964 int nfac = domaine.nb_faces_elem();
965 // int nsom = domaine.nb_som_elem();
966 // int nb_som_facette = domaine.type_elem()->nb_som_face();
967
968 const DoubleTab& xv = domaine_VEF.xv();
969 const DoubleTab& xp = domaine_VEF.xp();
970
971
972 double dt=equation().schema_temps().pas_de_temps();
973 const Champ_Inc_base& la_vitesse=vitesse_.valeur();
974
975 double dtrans_max, dtrans_min, dtrans_cen;//, dtrans_aval ;
976
977 int i, j, face_adj , face_glob ;
978
979 DoubleVect dist(dimension) ;
980 DoubleTab coor_trans(dimension,dimension) ;
981
982 DoubleVect trans_c_g(ncomp) ;
983 double flux;
984 IntVect face(nfac) ;
985
986 DoubleVect vs(dimension) ;
987
988 int amont, aval ;
989
990 if (psc >= 0)
991 {
992 amont = num1;
993 aval = num2;
994 }
995 else
996 {
997 amont = num2;
998 aval = num1;
999 }
1000
1001 if(first == 0)
1002 {
1003 trans_c_g = 0. ;
1004
1005 for(face_adj=0; face_adj<nfac; face_adj ++)
1006 {
1007 face_glob = elem_faces(poly, face_adj );
1008 face(face_adj) = face_glob ;
1009 for (i=0; i<ncomp; i++)
1010 trans_c_g(i) += transporte(face_glob,i)/double(nfac) ;
1011
1012 }
1013 int num3 = elem_faces(poly, 0 ) ;
1014
1015 coor_trans(0,0) = cos(3.14159/4.) ;
1016 coor_trans(0,1) = cos(3.14159/4.) ;
1017 coor_trans(0,2) = cos(3.14159/2.) ;
1018
1019 coor_trans(1,0) = cos(3.14159/2.-3.14159/4.) ;
1020 coor_trans(1,1) = cos(3.14159/2.+3.14159/4.) ;
1021 coor_trans(1,2) = cos(3.14159/2.) ;
1022
1023 coor_trans(2,0) = cos(3.14159/2.) ;
1024 coor_trans(2,1) = cos(3.14159/2.) ;
1025 coor_trans(2,2) = cos(0.) ;
1026
1027 //coor_trans = 0.;
1028 //coor_trans(0,0) = 1. ;
1029 //coor_trans(1,1) = 1. ;
1030 //coor_trans(2,2) = 1. ;
1031
1032 vs =0. ;
1033 for(face_adj=0; face_adj<nfac; face_adj ++)
1034 {
1035 face_glob = elem_faces(poly, face_adj );
1036 face(face_adj) = face_glob ;
1037 for (i=0; i<dimension; i++)
1038 vs(i) -= la_vitesse.valeurs()(face_glob,i)/double(nfac) ;
1039
1040 }
1041 for (i=0; i<dimension; i++) vs(i) += la_vitesse.valeurs()(num1,i) + la_vitesse.valeurs()(num2,i) ;
1042
1043 dist = 0. ;
1044
1045 for(face_adj=0; face_adj<nfac; face_adj ++)
1046 {
1047 face_glob = elem_faces(poly, face_adj);
1048 if ((face_glob == num1) || (face_glob == num2 ))
1049 for (i=0; i<dimension; i++)
1050 for (j=0; j<dimension; j++)
1051 dist(i) += xv(face_glob,j) * coor_trans(i,j) ;
1052 }
1053
1054 for(i=0; i<dimension; i++)
1055 for (j=0; j<dimension; j++)
1056 dist(i) -= (xp(poly,j) + xv(num3,j) + vs(j)*dt/2.) * coor_trans(i,j) ;
1057
1058 for (j=0; j<ncomp; j++)
1059 {
1060 if (sing == 0 )
1061 {
1062 if(ncomp == 1)
1063 {
1064 trans(j) = transporte(num3,j) + derive(8)*dist(0) * coor_trans(0,0)
1065 + derive(7)*dist(1) * coor_trans(1,1)
1066 + derive(6)*dist(2) * coor_trans(2,2)
1067 + 1./2.*( derive(5)*dist(0)*dist(0) + derive(4)*dist(1)*dist(1)
1068 + derive(3)*dist(2)*dist(2) )
1069 + derive(2)*dist(0)*dist(1) + derive(1)*dist(0)*dist(2)
1070 + derive(0)*dist(1)*dist(2) ;
1071
1072 dtrans_cen = transporte(amont) + transporte(aval) - trans_c_g(0) ;
1073 //dtrans_max = std::max( transporte(amont), transporte(aval) ) ;
1074 dtrans_max = std::max( transporte(amont), dtrans_cen ) ;
1075 dtrans_min = std::min( transporte(amont), dtrans_cen ) ;
1076 //dtrans_min = std::min( dtrans_min, dtrans_cen ) ;
1077 }
1078 else
1079 {
1080 trans(j) = transporte(num3,j) + derive(8,j)*dist(0) * coor_trans(0,0)
1081 + derive(7,j)*dist(1) * coor_trans(1,1)
1082 + derive(6,j)*dist(2) * coor_trans(2,2)
1083 + 1./2.*( derive(5,j)*dist(0)*dist(0) + derive(4,j)*dist(1)*dist(1)
1084 + derive(3,j)*dist(2)*dist(2) )
1085 + derive(2,j)*dist(0)*dist(1) + derive(1,j)*dist(0)*dist(2)
1086 + derive(0,j)*dist(1)*dist(2) ;
1087
1088 dtrans_cen = transporte(amont,j) + transporte(aval,j) - trans_c_g(j) ;
1089 dtrans_max = std::max( transporte(amont,j), dtrans_cen ) ;
1090
1091 dtrans_min = std::min( transporte(amont,j), dtrans_cen ) ;
1092
1093 }
1094
1095 if (trans(j) > dtrans_max )
1096 {
1097 //Cerr << " lim max in element " << poly << finl ;
1098 trans(j) = dtrans_max ;
1099 nlim++ ;
1100 }
1101 if (trans(j) < dtrans_min )
1102 {
1103 //Cerr << " lim min in element " << poly << finl ;
1104 trans(j) = dtrans_min ;
1105 nlim++ ;
1106 }
1107
1108 }
1109 else
1110 {
1111 if(ncomp == 1)
1112 {
1113 trans(0) = transporte(num1) + transporte(num2) - trans_c_g(0) ;
1114
1115 }
1116 else
1117 {
1118 trans(j) = transporte(num1,j) + transporte(num2,j) - trans_c_g(j) ;
1119 //Cerr << " singx != 0 " << finl ;
1120 }
1121 }
1122 }
1123 //double lo_x = xv(num1,0)+xv(num2,0)-xp(poly,0) ;
1124 //double lo_y = xv(num1,1)+xv(num2,1)-xp(poly,1) ;
1125 //double lo_z = xv(num1,2)+xv(num2,2)-xp(poly,2) ;
1126
1127 //double exact0= -sin(kwave*lo_x)*cos(kwave*lo_y);
1128 //double exact1= cos(kwave*lo_x)*sin(kwave*lo_y);
1129
1130 //double exact0 = lo_x*lo_x ;
1131 //double exact1 = lo_y*lo_y ;
1132 //double exact2 = lo_z*lo_z ;
1133
1134 //double moyenne0 = ( transporte(num1,0) + transporte(num2,0) ) - trans_c_g(0) ;
1135 //double moyenne1 = ( transporte(num1,1) + transporte(num2,1) ) - trans_c_g(1) ;
1136 //double moyenne2 = ( transporte(num1,2) + transporte(num2,2) ) - trans_c_g(2) ;
1137 //double moyenne0 = transporte(num1,0) ;
1138 //double moyenne1 = transporte(num1,1) ;
1139 //double moyenne2 = transporte(num1,2) ; ;
1140
1141 //double x0 = ( xv(num1,0) + xv(num2,0) ) - xp(poly,0) ;
1142 //double x1 = ( xv(num1,1) + xv(num2,1) ) - xp(poly,1) ;
1143 //double x2 = ( xv(num1,2) + xv(num2,2) ) - xp(poly,2) ;
1144 //double x0 = xv(num1,0) ;
1145 //double x1 = xv(num1,1) ;
1146 //double x2 = xv(num1,2) ;
1147
1148 //double y0 = xv(num3,0) + dist(0) ;
1149 //double y1 = xv(num3,1) + dist(1) ;
1150 //double y2 = xv(num3,2) + dist(2) ;
1151
1152 //if (sing == 0 ) Cerr << moyenne0 << " " <<trans(0) << " " << moyenne1 <<
1153 // " " << trans(1) << " " << moyenne2 << " " << trans(2) <<finl ;
1154 //Cerr << " exact0" << " " << exact0 << " " << "exact1" << " " << exact1
1155 // << " exact2" << " " << exact2 << finl ;
1156
1157 // Cerr << x0<< " " << y0 << " "<<x1 << " "<< y1<< " " << x2 << " "<< y2 << finl ;
1158 } // ende first = 0 und normal weiter
1159
1160 if (psc >= 0)
1161 {
1162 tab_fluent[num2] += psc;
1163 }
1164 else
1165 {
1166 tab_fluent[num1] -= psc;
1167 }
1168
1169 if (ncomp == 1)
1170 {
1171 flux = trans(0)*psc;
1172 resu(num1) -= flux;
1173 resu(num2) += flux;
1174 }
1175 else
1176 {
1177 for (i=0; i<ncomp; i++)
1178 {
1179 flux = trans(i)*psc;
1180 resu(num1,i) -= flux;
1181 resu(num2,i) += flux;
1182 }
1183 }
1184
1185}
1186void Op_Conv_DI_L2_VEF_Face::poly_DI_L2_2d(int N, int M, DoubleTab& derive ,int poly, int ncomp,
1187 const DoubleTab& transporte, int& sing) const
1188{
1189 //Cerr << " in poly_DI_L2_2d " << finl ;
1190 const Domaine_VEF& domaine_VEF = le_dom_vef.valeur();
1191
1192 const IntTab& elem_faces = domaine_VEF.elem_faces();
1193 const IntTab& face_voisins = domaine_VEF.face_voisins();
1194
1195 const Domaine& domaine = domaine_VEF.domaine();
1196
1197 int nfac = domaine.nb_faces_elem();
1198
1199 const DoubleTab& xv = domaine_VEF.xv();
1200 const DoubleTab& xp = domaine_VEF.xp();
1201
1202 int i, j, face_adj , face_glob , poly1 ;
1203 DoubleVect dist(dimension) ;
1204 DoubleTab coor_trans(dimension,dimension) ;
1205 IntVect face(nfac) ;
1206
1207 DoubleTab L(M,N) ;
1208 DoubleTab B(M,ncomp) ;
1209 DoubleVect dTransp_x(N) ;
1210 DoubleVect C(N) ;
1211 DoubleVect D(N) ;
1212 DoubleVect SM(N) ;
1213 double poid ;
1214
1215 for(face_adj=0; face_adj<nfac; face_adj ++)
1216 {
1217 face_glob = elem_faces(poly, face_adj );
1218 face(face_adj) = face_glob ;
1219 }
1220
1221 int num3 = elem_faces(poly, 0);
1222
1223 //coor_trans(0) = cos( (xv(num3,1) - xp(poly,1)) / (xv(num3,0) - xp(poly,0))) ;
1224 //coor_trans(1) = cos(1.- (xv(num3,1) - xp(poly,1)) / (xv(num3,0) - xp(poly,0))) ;
1225
1226 //coor_trans(0,0) = cos(3.14159/4.) ;
1227 //coor_trans(0,1) = cos(3.14159/4.) ;
1228 //coor_trans(1,0) = cos(3.14159/2.-3.14159/4.) ;
1229 //coor_trans(1,1) = cos(3.14159/2.+3.14159/4.) ;
1230 coor_trans(0,0) = 1. ;
1231 coor_trans(0,1) = 0. ;
1232 coor_trans(1,0) = 0. ;
1233 coor_trans(1,1) = 1. ;
1234 int row = -1 ;
1235
1236 for(int face_adj_poly =0; face_adj_poly < nfac; face_adj_poly ++)
1237 {
1238 poly1 = face_voisins(face[face_adj_poly],0);
1239 if (poly1 == poly ) poly1 = face_voisins(face[face_adj_poly], 1);
1240
1241 for(face_adj=0; face_adj<nfac; face_adj ++)
1242 {
1243 face_glob = elem_faces(poly1, face_adj);
1244 if (face_glob != num3 )
1245 {
1246 row++ ;
1247 dist = 0.;
1248 for (i=0; i<dimension; i++)
1249 for (j=0; j<dimension; j++)
1250 dist(i) += (xv(face_glob,j) - xv(num3,j)) * coor_trans(i,j) ;
1251
1252 poid = 0. ;
1253 for (i=0; i<dimension; i++)
1254 poid += ( (xv(face_glob,i)-xp(poly,i)) * (xv(face_glob,i)-xp(poly,i)) );
1255 poid = 1./sqrt(poid) ;
1256
1257 L(row,4) = poid * dist(0) * coor_trans(0,0) ;
1258 L(row,3) = poid * dist(1) * coor_trans(1,1) ;
1259 L(row,2) = poid * 1./2.*dist(0)*dist(0) ;
1260 L(row,1) = poid * 1./2.*dist(1)*dist(1) ;
1261 L(row,0) = poid * dist(0)*dist(1) ;
1262
1263 if (ncomp == 1)
1264 B(row,0) = poid * (transporte(face_glob) - transporte(num3) ) ;
1265 else for (j=0; j<ncomp; j++)
1266 B(row,j) = poid * (transporte(face_glob,j) - transporte(num3,j)) ;
1267
1268 }
1269 }
1270
1271 }
1272
1273 DoubleTab Lij(N,N) ;
1274 DoubleTab Bij(N,ncomp) ;
1275
1276 for (j=0; j<N; j++)
1277 for (i=0; i<ncomp; i++)
1278 for (int k=0; k<M; k++) Bij(j,i) += L(k,j) * B(k,i) ;
1279
1280 for (i=0; i<N; i++)
1281 for (j=0; j<N; j++)
1282 for (int k=0; k<M; k++) Lij(i,j) += L(k,i) * L(k,j) ;
1283
1284
1285 for (j=0; j<ncomp; j++)
1286 {
1287 for (i=0; i<N; i++) SM(i)=Bij(i,j) ;
1288
1289 qrsolv(Lij, N, SM, dTransp_x, sing, j, C, D);
1290
1291 if (sing == 0 )
1292 {
1293 if(ncomp == 1)
1294 {
1295 for (int val_N= 0; val_N < N ; val_N++ )
1296 derive(val_N) = dTransp_x(val_N) ;
1297 }
1298 else
1299 {
1300 for (int val_N= 0; val_N < N ; val_N++ )
1301 derive(val_N,j) = dTransp_x(val_N) ;
1302 }
1303 }
1304 }
1305}
1306
1307void Op_Conv_DI_L2_VEF_Face::poly_DI_L2_3d(int N, int M, DoubleTab& derive ,int poly, int ncomp,
1308 const DoubleTab& transporte ,int& sing) const
1309{
1310 // Cerr<<"Op_Conv_DI_L2_VEF_Face::poly_DI_L2_3d 0 "<<finl;
1311 const Domaine_VEF& domaine_VEF = le_dom_vef.valeur();
1312
1313 const IntTab& elem_faces = domaine_VEF.elem_faces();
1314 const IntTab& face_voisins = domaine_VEF.face_voisins();
1315
1316 const Domaine& domaine = domaine_VEF.domaine();
1317
1318 int nfac = domaine.nb_faces_elem();
1319
1320 const DoubleTab& xv = domaine_VEF.xv();
1321 const DoubleTab& xp = domaine_VEF.xp();
1322
1323 int i, j, face_adj , face_glob , poly1 ;
1324
1325 DoubleVect dist(dimension) ;
1326 DoubleTab coor_trans(dimension,dimension) ;
1327 IntVect face(nfac) ;
1328
1329 DoubleTab L(M,N) ;
1330 DoubleTab B(M,ncomp) ;
1331 DoubleVect dTransp_x(N) ;
1332 DoubleVect C(N) ;
1333 DoubleVect D(N) ;
1334 DoubleVect SM(N) ;
1335 double poid ;
1336 // Cerr<<"Op_Conv_DI_L2_VEF_Face::poly_DI_L2_3d 1 "<<finl;
1337
1338 for(face_adj=0; face_adj<nfac; face_adj ++)
1339 {
1340 face_glob = elem_faces(poly, face_adj );
1341 face(face_adj) = face_glob ;
1342 }
1343 // Cerr<<"Op_Conv_DI_L2_VEF_Face::poly_DI_L2_3d 2 "<<finl;
1344
1345 int num3 = elem_faces(poly, 0);
1346
1347 coor_trans(0,0) = cos(3.14159/4.) ;
1348 coor_trans(0,1) = cos(3.14159/4.) ;
1349 coor_trans(0,2) = cos(3.14159/2.) ;
1350
1351 coor_trans(1,0) = cos(3.14159/2.-3.14159/4.) ;
1352 coor_trans(1,1) = cos(3.14159/2.+3.14159/4.) ;
1353 coor_trans(1,2) = cos(3.14159/2.) ;
1354
1355 coor_trans(2,0) = cos(3.14159/2.) ;
1356 coor_trans(2,1) = cos(3.14159/2.) ;
1357 coor_trans(2,2) = cos(0.) ;
1358
1359 //coor_trans = 0.;
1360 //coor_trans(0,0) = 1. ;
1361 //coor_trans(1,1) = 1. ;
1362 //coor_trans(2,2) = 1. ;
1363
1364 int row = -1 ;
1365 // Cerr<<"Op_Conv_DI_L2_VEF_Face::poly_DI_L2_3d 3 "<<finl;
1366
1367 for(int face_adj_poly =0; face_adj_poly < nfac; face_adj_poly ++)
1368 {
1369 poly1 = face_voisins(face[face_adj_poly],0);
1370 // Cerr<<"avant poly = "<<poly<<", poly1 = "<<poly1<<", proc = "<<me()<<finl;
1371 // Cerr<<"face[face_adj_poly] = "<<face[face_adj_poly]<<finl;
1372 if (poly1 == poly ) poly1 = face_voisins(face[face_adj_poly], 1);
1373 // Cerr<<"apres poly = "<<poly<<", poly1 = "<<poly1<<", proc = "<<me()<<finl;
1374 for(face_adj=0; face_adj<nfac; face_adj ++)
1375 {
1376 face_glob = elem_faces(poly1, face_adj);
1377 if (face_glob != num3 )
1378 {
1379 row++ ;
1380 dist=0 ;
1381 for (i=0; i<dimension; i++)
1382 for (j=0; j<dimension; j++)
1383 dist(i) += (xv(face_glob,j) - xv(num3,j)) * coor_trans(i,j) ;
1384
1385 poid = 0. ;
1386 for (i=0; i<dimension; i++)
1387 poid += ( (xv(face_glob,i)-xp(poly,i)) * (xv(face_glob,i)-xp(poly,i)) );
1388 poid = 1./(poid) ;
1389
1390
1391 L(row,8) = poid*dist(0) * coor_trans(0,0) ;
1392 L(row,7) = poid*dist(1) * coor_trans(1,1) ;
1393 L(row,6) = poid*dist(2) * coor_trans(2,2) ;
1394 L(row,5) = poid*1./2.*dist(0)*dist(0) ;
1395 L(row,4) = poid*1./2.*dist(1)*dist(1) ;
1396 L(row,3) = poid*1./2.*dist(2)*dist(2) ;
1397 L(row,2) = poid*dist(0)*dist(1) ;
1398 L(row,1) = poid*dist(0)*dist(2) ;
1399 L(row,0) = poid*dist(1)*dist(2) ;
1400
1401 if (ncomp == 1)
1402 B(row,0) = poid * (transporte(face_glob) - transporte(num3) ) ;
1403 else for (j=0; j<ncomp; j++)
1404 B(row,j) = poid * (transporte(face_glob,j) - transporte(num3,j)) ;
1405 }
1406 }
1407 }
1408
1409 // Cerr<<"Op_Conv_DI_L2_VEF_Face::poly_DI_L2_3d 4 "<<finl;
1410
1411 DoubleTab Lij(N,N) ;
1412 DoubleTab Bij(N,ncomp) ;
1413
1414 for (j=0; j<N; j++)
1415 for (i=0; i<ncomp; i++)
1416 for (int k=0; k<M; k++) Bij(j,i) += L(k,j) * B(k,i) ;
1417
1418 // Cerr<<"Op_Conv_DI_L2_VEF_Face::poly_DI_L2_3d 5 "<<finl;
1419
1420 for (i=0; i<N; i++)
1421 for (j=0; j<N; j++)
1422 for (int k=0; k<M; k++) Lij(i,j) += L(k,i) * L(k,j) ;
1423 // Cerr<<"Op_Conv_DI_L2_VEF_Face::poly_DI_L2_3d 6 "<<finl;
1424
1425
1426 for (j=0; j<ncomp; j++)
1427 {
1428 for (i=0; i<N; i++) SM(i)=Bij(i,j) ;
1429
1430 qrsolv(Lij, N, SM, dTransp_x, sing, j, C, D);
1431
1432 if (sing == 0 )
1433 {
1434 if(ncomp == 1)
1435 {
1436 for (int val_N= 0; val_N < N ; val_N++ )
1437 derive(val_N) = dTransp_x(val_N) ;
1438 }
1439 else
1440 {
1441 for (int val_N= 0; val_N < N ; val_N++ )
1442 derive(val_N,j) = dTransp_x(val_N) ;
1443 }
1444 }
1445 }
1446
1447 // Cerr<<"Op_Conv_DI_L2_VEF_Face::poly_DI_L2_3d 7 "<<finl;
1448
1449}
1450
Champ_Inc_base
Class Champ_Inc_base.
Definition Champ_Inc_base.h:53
Champ_Inc_base::valeurs
DoubleTab & valeurs() override
Returns the array of field values at the current time.
Definition Champ_Inc_base.h:77
Champ_P1NC
Definition Champ_P1NC.h:35
Champ_base
class Champ_base This class is the base of the fields hierarchy.
Definition Champ_base.h:43
Cond_lim
class Cond_lim Generic class used to represent any class
Definition Cond_lim.h:31
Domaine_32_64::nb_faces_elem
int nb_faces_elem(int=0) const
Returns the number of faces of type i of the geometric elements that make up the domain.
Definition Domaine.h:484
Domaine_Cl_VEF
Definition Domaine_Cl_VEF.h:38
Domaine_Cl_VEF::type_elem_Cl
int type_elem_Cl(int i) const
Definition Domaine_Cl_VEF.h:60
Domaine_Cl_VEF::normales_facettes_Cl
DoubleTab & normales_facettes_Cl()
Definition Domaine_Cl_VEF.h:54
Domaine_Cl_dis_base::les_conditions_limites
const Cond_lim & les_conditions_limites(int) const
Returns the i-th boundary condition.
Definition Domaine_Cl_dis_base.cpp:386
Domaine_VEF
class Domaine_VEF
Definition Domaine_VEF.h:53
Domaine_VEF::rang_elem_non_std
IntVect & rang_elem_non_std()
Definition Domaine_VEF.h:85
Domaine_VEF::type_elem
const Elem_VEF_base & type_elem() const
Definition Domaine_VEF.h:74
Domaine_VEF::facette_normales
auto & facette_normales()
Definition Domaine_VEF.h:83
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
Returns the index of the i-th face of element num_elem; the face numbering convention is.
Definition Domaine_VF.h:542
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
Returns the neighbouring element of num_face in direction i.
Definition Domaine_VF.h:418
Domaine_VF::nb_faces_bord
int nb_faces_bord() const
Returns the number of faces on which boundary conditions are applied:
Definition Domaine_VF.h:512
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()
Returns the time scheme associated with the 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
Returns the reference to the equation pointed to by MorEqn::mon_equation.
Definition MorEqn.h:62
Neumann_sortie_libre
Neumann_sortie_libre This class represents an open boundary without imposed velocity.
Definition Neumann_sortie_libre.h:34
Neumann_sortie_libre::val_ext
double val_ext(int i) const override
Returns the value of the i-th component of the field imposed on the exterior of the boundary.
Definition Neumann_sortie_libre.cpp:120
Objet_U::dimension
static int dimension
Definition Objet_U.h:94
Objet_U::que_suis_je
const Nom & que_suis_je() const
Returns the string identifying the class.
Definition Objet_U.cpp:104
Objet_U::readOn
virtual Entree & readOn(Entree &)
Reads an Objet_U from an input stream. Virtual method to override.
Definition Objet_U.cpp:289
Objet_U::printOn
virtual Sortie & printOn(Sortie &) const
Writes the object to an output stream. Virtual method to override.
Definition Objet_U.cpp:278
Op_Conv_DI_L2_VEF_Face
class Op_Conv_DI_L2_VEF_Face
Definition Op_Conv_DI_L2_VEF_Face.h:33
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:1307
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:1186
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:941
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:750
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
class Periodique This class represents a periodic boundary condition.
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:108
Process::exit
static void exit(int exit_code=-1)
Exit routine for TRUST within a Kokkos region.
Definition Process.cpp:466
Process::is_sequential
static bool is_sequential()
Definition Process.cpp:113
Schema_Temps_base::pas_de_temps
double pas_de_temps() const
Returns the current time step (delta_t).
Definition Schema_Temps_base.h:394
Sortie
Base class for output streams.
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