165 DoubleTab& resu)
const
168 const Domaine_VEF& domaine_VEF = le_dom_vef.valeur();
170 const IntTab& elem_faces = domaine_VEF.
elem_faces();
171 const DoubleTab& face_normales = domaine_VEF.
face_normales();
174 const Domaine& domaine = domaine_VEF.
domaine();
175 const int nb_faces = domaine_VEF.
nb_faces();
177 const int nb_elem = domaine_VEF.
nb_elem();
180 const IntTab& face_voisins = domaine_VEF.
face_voisins();
182 const DoubleVect& volumes = domaine_VEF.
volumes();
183 const DoubleTab& xv = domaine_VEF.
xv();
184 const DoubleTab& xg = domaine_VEF.
xp();
185 const DoubleTab& coord = domaine.coord_sommets();
187 const IntTab& les_Polys = domaine.les_elems();
193 int nfac = domaine.nb_faces_elem();
194 int nsom = domaine.nb_som_elem();
195 int nb_som_facette = domaine.type_elem()->nb_som_face();
210 if ((nom_elem==
"Tetra_VEF")||(nom_elem==
"Tri_VEF"))
224 int poly,poly1,poly2,face_adj,fa7,i,j,n_bord;
225 int num_face, rang ,itypcl;
226 int num10,num20,num3,num_som;
227 int ncomp_ch_transporte;
229 if (transporte.
nb_dim() == 1)
230 ncomp_ch_transporte=1;
232 ncomp_ch_transporte= transporte.
dimension(1);
234 int fac,elem1,elem2,comp0;
235 int nb_faces_ = domaine_VEF.
nb_faces();
238 DoubleVect flux(ncomp_ch_transporte);
239 DoubleVect fluxsom(ncomp_ch_transporte);
240 DoubleVect fluxg(ncomp_ch_transporte);
244 int nb_faces_perio = 0;
245 for (n_bord=0; n_bord<domaine_VEF.
nb_front_Cl(); n_bord++)
252 int num2 = num1 + le_bord.
nb_faces();
253 for (num_face=num1; num_face<num2; num_face++)
259 if (ncomp_ch_transporte == 1)
260 tab.
resize(nb_faces_perio);
262 tab.
resize(nb_faces_perio,ncomp_ch_transporte);
265 for (n_bord=0; n_bord<domaine_VEF.
nb_front_Cl(); n_bord++)
273 int num2 = num1 + le_bord.
nb_faces();
274 for (num_face=num1; num_face<num2; num_face++)
276 if (ncomp_ch_transporte == 1)
277 tab(nb_faces_perio) = resu(num_face);
279 for (
int comp=0; comp<ncomp_ch_transporte; comp++)
280 tab(nb_faces_perio,comp) = resu(num_face,comp);
293 DoubleTab gradient_elem(0, ncomp_ch_transporte,
dimension);
299 for (fac=0; fac< premiere_face_int; fac++)
301 elem1=face_voisins(fac,0);
302 if(ncomp_ch_transporte==1)
305 gradient_elem(elem1, 0, i) +=
306 face_normales(fac,i)*transporte(fac);
309 for (comp0=0; comp0<ncomp_ch_transporte; comp0++)
311 gradient_elem(elem1, comp0, i) +=
312 face_normales(fac,i)*transporte(fac,comp0);
316 for (; fac<nb_faces_; fac++)
318 elem1=face_voisins(fac,0);
319 elem2=face_voisins(fac,1);
320 if(ncomp_ch_transporte==1)
323 gradient_elem(elem1, 0, i) +=
324 face_normales(fac,i)*transporte(fac);
325 gradient_elem(elem2, 0, i) -=
326 face_normales(fac,i)*transporte(fac);
329 for (comp0=0; comp0<ncomp_ch_transporte; comp0++)
332 gradient_elem(elem1, comp0, i) +=
333 face_normales(fac,i)*transporte(fac,comp0);
334 gradient_elem(elem2, comp0, i) -=
335 face_normales(fac,i)*transporte(fac,comp0);
340 for (
int elem=0; elem<nb_elem; elem++)
341 for (comp0=0; comp0<ncomp_ch_transporte; comp0++)
343 gradient_elem(elem,comp0,i) /= volumes(elem);
371 const IntTab& KEL=type_elemvef.
KEL();
372 for (poly=0; poly<nb_elem_tot; poly++)
376 rang = rang_elem_non_std(poly);
384 for (face_adj=0; face_adj<nfac; face_adj++)
386 face(face_adj)= elem_faces(poly,face_adj);
397 vs(j) = la_vitesse.
valeurs()(face(0),j)*porosite_face(face(0));
398 for (i=1; i<nfac; i++)
399 vs(j)+= la_vitesse.
valeurs()(face(i),j)*porosite_face(face(i));
405 for (j=0; j<nsom; j++)
416 for (j=0; j<nsom; j++)
418 num_som = domaine.sommet_elem(poly,j);
433 for (fa7=0; fa7<nfa7; fa7++)
438 num10 = face(KEL(0,fa7));
439 num20 = face(KEL(1,fa7));
440 num3 = face(KEL(2,fa7));
444 poly1 = face_voisins(num10,0);
447 poly1 = face_voisins(num10,1);
450 poly2 = face_voisins(num20,0);
453 poly2 = face_voisins(num20,1);
456 scom = les_Polys(poly,KEL(2,fa7));
461 rx0(i) = xv(num20,i)-xv(num10,i);
467 cc[i] = facette_normales(poly, fa7, i);
470 cc[i] = normales_facettes_Cl(rang,fa7,i);
478 for (i=0; i<nb_som_facette-1; i++)
487 psc+=((vsom(KEL(i+2,fa7),j) + la_vitesse.
valeurs()(num3,j) * porosite_face(num3)))*cc[j];
490 coord_som(j)=coord(scom,j);
492 convkschemas_centre(
K,ncomp_ch_transporte,
dimension,poly,poly1,poly2,num10,num20,psc,transporte,
493 fluent_,flux,rx0,gradient_elem,coord_som);
508 if (ncomp_ch_transporte == 1)
512 xm = 0.5 *(coord(scom,j)+xv(num3,j));
517 fluxsom(0) += gradient_elem(poly,0,j)*(coord(scom,j)-xm);
521 fluxsom(0) += (teta*gradient_elem(poly,0,j) + (1.- teta)*gradient_elem(poly1,0,j))*(coord(scom,j)-xm);
529 fluxsom(0) += gradient_elem(poly,0,j)*(coord(scom,j)-xm);
533 fluxsom(0) += (teta*gradient_elem(poly,0,j) + (1.- teta)*gradient_elem(poly2,0,j))*(coord(scom,j)-xm);
537 fluxsom(0) += flux(0);
543 for (comp0=0; comp0<ncomp_ch_transporte; comp0++)
545 fluxsom(comp0) = flux(comp0);
548 xm = 0.5 *(coord(scom,j)+xv(num3,j));
553 fluxsom(comp0) += gradient_elem(poly,comp0,j)*(coord(scom,j)-xm);
557 fluxsom(comp0) += (teta*gradient_elem(poly,comp0,j) + (1.-teta)*gradient_elem(poly1,comp0,j))*(coord(scom,j)-xm);
565 fluxsom(comp0) += gradient_elem(poly,comp0,j)*(coord(scom,j)-xm);
569 fluxsom(comp0) += (teta*gradient_elem(poly,comp0,j) + (1.-teta)*gradient_elem(poly2,comp0,j))*(coord(scom,j)-xm);
575 fluxsom(comp0) *= psc;
585 if (ncomp_ch_transporte == 1)
589 xm = 0.5 *(coord(scom,j)+xv(num3,j));
594 fluxg(0) += gradient_elem(poly,0,j)*(xg(poly,j)-xm);
598 fluxg(0) += (teta*gradient_elem(poly,0,j) + (1.- teta)*gradient_elem(poly1,0,j))*(xg(poly,j)-xm);
606 fluxg(0) += gradient_elem(poly,0,j)*(xg(poly,j)-xm);
610 fluxg(0) += (teta*gradient_elem(poly,0,j) + (1.- teta)*gradient_elem(poly2,0,j))*(xg(poly,j)-xm);
621 for (comp0=0; comp0<ncomp_ch_transporte; comp0++)
623 fluxg(comp0) = flux(comp0);
626 xm = 0.5 *(coord(scom,j)+xv(num3,j));
631 fluxg(comp0) += gradient_elem(poly,comp0,j)*(xg(poly,j)-xm);
635 fluxg(comp0) += (teta*gradient_elem(poly,comp0,j) + (1.-teta)*gradient_elem(poly1,comp0,j))*(xg(poly,j)-xm);
643 fluxg(comp0) += gradient_elem(poly,comp0,j)*(xg(poly,j)-xm);
647 fluxg(comp0) += (teta*gradient_elem(poly,comp0,j) + (1.-teta)*gradient_elem(poly2,comp0,j))*(xg(poly,j)-xm);
662 if (ncomp_ch_transporte == 1)
664 resu(num10) -=flux(0)*psc;
665 resu(num20) += flux(0)*psc;
670 for (comp0=0; comp0<ncomp_ch_transporte; comp0++)
672 resu(num10,comp0) -= flux(comp0)*psc;
673 resu(num20,comp0) += flux(comp0)*psc;
681 for (poly=0; poly<nb_elem_tot; poly++)
684 for (face_adj=0; face_adj<nfac; face_adj++)
685 if(face_adj<nb_faces)
break;
688 rang = rang_elem_non_std(poly);
696 for (face_adj=0; face_adj<nfac; face_adj++)
698 face(face_adj)= elem_faces(poly,face_adj);
708 vs(j) = la_vitesse.
valeurs()(face(0),j)*porosite_face(face(0));
709 for (i=1; i<nfac; i++)
710 vs(j)+= la_vitesse.
valeurs()(face(i),j)*porosite_face(face(i));
716 for (j=0; j<nsom; j++)
727 for (j=0; j<nsom; j++)
729 num_som = domaine.sommet_elem(poly,j);
744 for (fa7=0; fa7<nfa7; fa7++)
749 num10 = face(KEL(0,fa7));
750 num20 = face(KEL(1,fa7));
751 num3 = face(KEL(2,fa7));
755 poly1 = face_voisins(num10,0);
758 poly1 = face_voisins(num10,1);
761 poly2 = face_voisins(num20,0);
764 poly2 = face_voisins(num20,1);
767 scom = les_Polys(poly,KEL(2,fa7));
772 rx0(i) = xv(num20,i)-xv(num10,i);
778 cc[i] = facette_normales(poly, fa7, i);
781 cc[i] = normales_facettes_Cl(rang,fa7,i);
789 for (i=0; i<nb_som_facette-1; i++)
798 psc+=((vsom(KEL(i+2,fa7),j) + la_vitesse.
valeurs()(num3,j) * porosite_face(num3)))*cc[j];
801 coord_som(j)=coord(scom,j);
802 convkschemas_centre(
K,ncomp_ch_transporte,
dimension,poly,poly1,poly2,num10,num20,psc,transporte,
803 fluent_,flux,rx0,gradient_elem,coord_som);
818 if (ncomp_ch_transporte == 1)
822 xm = 0.5 *(coord(scom,j)+xv(num3,j));
827 fluxsom(0) += gradient_elem(poly,0,j)*(coord(scom,j)-xm);
831 fluxsom(0) += (teta*gradient_elem(poly,0,j) + (1.- teta)*gradient_elem(poly1,0,j))*(coord(scom,j)-xm);
839 fluxsom(0) += gradient_elem(poly,0,j)*(coord(scom,j)-xm);
843 fluxsom(0) += (teta*gradient_elem(poly,0,j) + (1.- teta)*gradient_elem(poly2,0,j))*(coord(scom,j)-xm);
847 fluxsom(0) += flux(0);
853 for (comp0=0; comp0<ncomp_ch_transporte; comp0++)
855 fluxsom(comp0) = flux(comp0);
858 xm = 0.5 *(coord(scom,j)+xv(num3,j));
863 fluxsom(comp0) += gradient_elem(poly,comp0,j)*(coord(scom,j)-xm);
867 fluxsom(comp0) += (teta*gradient_elem(poly,comp0,j) + (1.-teta)*gradient_elem(poly1,comp0,j))*(coord(scom,j)-xm);
875 fluxsom(comp0) += gradient_elem(poly,comp0,j)*(coord(scom,j)-xm);
879 fluxsom(comp0) += (teta*gradient_elem(poly,comp0,j) + (1.-teta)*gradient_elem(poly2,comp0,j))*(coord(scom,j)-xm);
885 fluxsom(comp0) *= psc;
895 if (ncomp_ch_transporte == 1)
899 xm = 0.5 *(coord(scom,j)+xv(num3,j));
904 fluxg(0) += gradient_elem(poly,0,j)*(xg(poly,j)-xm);
908 fluxg(0) += (teta*gradient_elem(poly,0,j) + (1.- teta)*gradient_elem(poly1,0,j))*(xg(poly,j)-xm);
916 fluxg(0) += gradient_elem(poly,0,j)*(xg(poly,j)-xm);
920 fluxg(0) += (teta*gradient_elem(poly,0,j) + (1.- teta)*gradient_elem(poly2,0,j))*(xg(poly,j)-xm);
931 for (comp0=0; comp0<ncomp_ch_transporte; comp0++)
933 fluxg(comp0) = flux(comp0);
936 xm = 0.5 *(coord(scom,j)+xv(num3,j));
941 fluxg(comp0) += gradient_elem(poly,comp0,j)*(xg(poly,j)-xm);
945 fluxg(comp0) += (teta*gradient_elem(poly,comp0,j) + (1.-teta)*gradient_elem(poly1,comp0,j))*(xg(poly,j)-xm);
953 fluxg(comp0) += gradient_elem(poly,comp0,j)*(xg(poly,j)-xm);
957 fluxg(comp0) += (teta*gradient_elem(poly,comp0,j) + (1.-teta)*gradient_elem(poly2,comp0,j))*(xg(poly,j)-xm);
972 if (ncomp_ch_transporte == 1)
974 resu(num10) -= flux(0)*psc;
975 resu(num20) += flux(0)*psc;
980 for (comp0=0; comp0<ncomp_ch_transporte; comp0++)
982 resu(num10,comp0) -= flux(comp0)*psc;
983 resu(num20,comp0) +=flux(comp0)*psc;
1004 for (n_bord=0; n_bord<domaine_VEF.
nb_front_Cl(); n_bord++)
1013 int num2 = num1 + le_bord.
nb_faces();
1014 for (num_face=num1; num_face<num2; num_face++)
1018 psc += la_vitesse.
valeurs()(num_face,i)*face_normales(num_face,i)*porosite_face(num_face);
1020 if (ncomp_ch_transporte == 1)
1022 resu(num_face) -= psc*transporte(num_face);
1023 flux_b(num_face,0) -= psc*transporte(num_face);
1026 for (i=0; i<ncomp_ch_transporte; i++)
1028 resu(num_face,i) -= psc*transporte(num_face,i);
1029 flux_b(num_face,i) -= psc*transporte(num_face,i);
1033 if (ncomp_ch_transporte == 1)
1035 resu(num_face) -= psc*la_sortie_libre.
val_ext(num_face-num1);
1036 flux_b(num_face,0) -= psc*la_sortie_libre.
val_ext(num_face-num1);
1039 for (i=0; i<ncomp_ch_transporte; i++)
1041 resu(num_face,i) -= psc*la_sortie_libre.
val_ext(num_face-num1,i);
1042 flux_b(num_face,i) -= psc*la_sortie_libre.
val_ext(num_face-num1,i);
1048 else if (sub_type(
Periodique,la_cl.valeur()))
1053 int num2 = num1 + le_bord.
nb_faces();
1056 for (num_face=num1; num_face<num2; num_face++)
1058 if (fait[num_face-num1] == 0)
1062 if (ncomp_ch_transporte == 1)
1064 diff1 = resu(num_face)-tab(nb_faces_perio);
1065 diff2 = resu(voisine)-tab(nb_faces_perio+voisine-num_face);
1066 resu(voisine) += diff1;
1067 resu(num_face) += diff2;
1068 flux_b(voisine,0) += diff1;
1069 flux_b(num_face,0) += diff2;
1072 for (
int comp=0; comp<ncomp_ch_transporte; comp++)
1074 diff1 = resu(num_face,comp)-tab(nb_faces_perio,comp);
1075 diff2 = resu(voisine,comp)-tab(nb_faces_perio+voisine-num_face,comp);
1076 resu(voisine,comp) += diff1;
1077 resu(num_face,comp) += diff2;
1078 flux_b(num_face,comp) += diff1;
1079 flux_b(num_face,comp) += diff2;
1082 fait[num_face-num1]= 1;
1083 fait[voisine-num1] = 1;