16#include <Transport_Interfaces_FT_Disc.h>
17#include <Probleme_FT_Disc_gen.h>
20#include <EcritureLectureSpecial.h>
21#include <Discret_Thyd.h>
22#include <Schema_Temps_base.h>
24#include <Domaine_VF.h>
25#include <Domaine_VDF.h>
26#include <Sous_Domaine.h>
27#include <Champ_Face_VDF.h>
28#include <Champ_P1NC.h>
29#include <Champ_Fonc_P1NC.h>
30#include <Fluide_Diphasique.h>
31#include <Fluide_Incompressible.h>
32#include <Perf_counters.h>
33#include <Paroi_FT_disc.h>
34#include <Champ_Fonc_Face_VDF.h>
36#include <LecFicDiffuse.h>
37#include <Sauvegarde_Reprise_Maillage_FT.h>
40#include <Format_Post_Lata.h>
41#include <Connex_components.h>
42#include <Connex_components_FT.h>
43#include <Loi_horaire.h>
44#include <Transport_Marqueur_FT.h>
45#include <communications.h>
46#include <Convection_Diffusion_Temperature_FT_Disc.h>
48#include <Champ_Fonc_P0_base.h>
49#include <EFichierBin.h>
55#include <Dirichlet_paroi_fixe.h>
56#include <Dirichlet_paroi_defilante.h>
57#include <Echange_contact_VDF_FT_Disc.h>
58#include <TRUST_2_PDI.h>
71static void eval_vitesse(
double x,
double y,
double z,
double t,
73 double& vx,
double& vy,
double& vz)
75 int i0=0, i1=1, i2=2, i3=3;
94static void integrer_vitesse_imposee(
96 double temps,
double dt,
double& x,
double& y,
double& z)
104 parser_vx,parser_vy,parser_vz,vx,vy,vz);
106 eval_vitesse(x + vx * dt * 0.5,
110 parser_vx,parser_vy,parser_vz,vx1,vy1,vz1);
112 eval_vitesse(x + vx1 * dt * 0.5,
116 parser_vx,parser_vy,parser_vz,vx2,vy2,vz2);
118 eval_vitesse(x + vx2 * dt,
122 parser_vx,parser_vy,parser_vz,vx3,vy3,vz3);
124 x += (vx + 2.*vx1 + 2.*vx2 + vx3) / 6. * dt;
125 y += (vy + 2.*vy1 + 2.*vy2 + vy3) / 6. * dt;
126 z += (vz + 2.*vz1 + 2.*vz2 + vz3) / 6. * dt;
138static void calculer_normale_sommets_interface(
const Maillage_FT_Disc& maillage,
146 const IntTab& facettes = maillage.
facettes();
147 const int nsommets_faces = facettes.
line_size();
149 normale.
resize(nsom, dim);
151 double n[3]= {0,0,0};
153 for (
int i = 0; i < nfaces; i++)
160 const double surface = surface_facettes[i];
161 for (
int k = 0; k < dim; k++)
162 n[k] = normale_facettes(i,k) * surface;
164 for (
int j = 0; j < nsommets_faces; j++)
166 const int sommet = facettes(i,j);
167 for (
int k = 0; k < dim; k++)
168 normale(sommet, k) += n[k];
192 int index = fullname.
find(
".");
193 std::string fname = index ==-1 ? fullname.
getString() : fullname.
getString().substr(0, index);
194 Nom nom_fich_xyz(fname);
197 nom_fich_xyz +=
"_maillage_interface";
198 nom_fich_xyz +=
".xyz";
206 Cerr <<
"Transport_Interfaces_FT_Disc::init_save_file Backup file for the FT mesh has already been initialized!" << finl;
213 fic_front_sauv_->ouvrir(nom_fich_xyz);
220 if (!fic_front_sauv_)
222 Cerr <<
"Transport_Interfaces_FT_Disc::close_save_file Backup file for the FT mesh has not been initialized!" << finl;
228 fic_front_sauv_.valeur() <<
Nom(
"fin");
229 (fic_front_sauv_.valeur()).flush();
230 (fic_front_sauv_.valeur()).syncfile();
239 std::vector<YAML_data> data = indicatrice_cache->data_a_sauvegarder();
242 std::string name =
pb_name_.getString() +
"_" +
ident_.getString() +
"_indicatrice_cache_tag";
245 data.push_back(ictag);
255 bytes += indicatrice_cache->sauvegarder(os);
268 std::string name =
pb_name_.getString() +
"_" +
ident_.getString() +
"_indicatrice_cache_tag";
277 Cerr <<
"Transport_Interfaces_FT_Disc_interne::sauvegarder : WARNING ! PDI format is not supported for the FT mesh so we save it in xyz format " << finl;
285 xyz_os << mon_ident << finl;
294 xyz_os << mon_ident << finl;
311 Cerr <<
"Transport_Interfaces_FT_Disc_interne::reprendre" << finl;
317 Nom name =
pb_name_ +
"_" + indicatrice_cache->le_nom();
327 if (!indicatrice_cache)
330 indicatrice_cache.typer(type);
336 Cerr <<
"Transport_Interfaces_FT_Disc_interne::reprendre : WARNING ! PDI format is not supported for the FT mesh so we use xyz format" << finl;
342 xyz_is->ouvrir(nomfic);
345 Cerr <<
"Error during the opening of the restart file : " << nomfic << finl;
349 Cerr <<
"Reading interface in file " << nomfic << finl;
353 avancer_fichier(xyz_is, mon_ident);
358 avancer_fichier(xyz_is, mon_ident);
380 interpolation_repere_local_ = 0;
388Transport_Interfaces_FT_Disc::~Transport_Interfaces_FT_Disc()
390 delete variables_internes_;
427 param.
ajouter(
"VOFlike_correction_volume",&variables_internes_->VOFlike_correction_volume);
428 param.
ajouter(
"nb_lissage_correction_volume",&variables_internes_->nb_lissage_correction_volume);
429 param.
ajouter(
"nb_iterations_correction_volume",&variables_internes_->nb_iterations_correction_volume);
432 param.
ajouter(
"parcours_interface",&variables_internes_->parcours_interface_);
435 param.
ajouter(
"sous_domaine_volume_impose",&variables_internes_->nom_domaine_volume_impose_);
436 param.
ajouter_flag(
"interpolation_repere_local", &interpolation_repere_local_);
438 param.
ajouter(
"n_iterations_interpolation_ibc",&variables_internes_->n_iterations_interpolation_ibc);
440 param.
ajouter(
"nombre_facettes_retenues_par_cellule",&variables_internes_->nomb_fa7_accepted);
441 param.
ajouter(
"seuil_convergence_uzawa",&variables_internes_->seuil_uzawa) ;
442 param.
ajouter(
"nb_iteration_max_uzawa",&variables_internes_->nb_iter_uzawa) ;
445 param.
ajouter(
"temps_debut",&temps_debut_);
446 param.
ajouter(
"vitesse_imposee_regularisee", &variables_internes_->vimp_regul) ;
449 param.
ajouter(
"collision_model_fpi",&collision_model_);
456 if (un_mot==
"maillage")
462 else if (un_mot==
"methode_transport")
465 motcles2[0] =
"vitesse_imposee";
466 motcles2[1] =
"loi_horaire";
467 motcles2[2] =
"vitesse_interpolee";
470 Cerr << un_mot <<
" " << motlu << finl;
471 int rang = motcles2.
search(motlu);
478 <<
" analytical expressions vx(x,y,[z,]t) vy(...) [vz(...)]" << finl;
480 variables_internes_->methode_transport =
488 Nom& expression = variables_internes_->expression_vitesse_imposee[i];
491 Cerr <<
" Component " << i <<
" of the velocity : " << expression << finl;
497 Cerr <<
" Imposed motion by the time scheme : ";
501 Cerr << nom_loi << finl;
508 Cerr <<
" Transport velocity interpolated from the velocity \n"
509 <<
" of the following Navier_Stokes equation :" << finl;
511 variables_internes_->methode_transport =
516 Cerr << nom_eq << finl;
517 variables_internes_->refequation_vitesse_transport =
519 probleme_base_.valeur()).equation_hydraulique(nom_eq);
521 const Equation_base& eqn_hydraulique = variables_internes_->refequation_vitesse_transport.valeur();
534 Cerr <<
"Transport_Interfaces_FT_Disc::lire\n"
535 <<
"The options for methode_transport are :\n"
542 else if (un_mot==
"n_iterations_distance")
546 variables_internes_->n_iterations_distance = n;
548 Cerr <<
"Iterations number of the smoothing tool \"distance_interface\": " << n << finl;
551 else if (un_mot==
"iterations_correction_volume")
555 variables_internes_->iterations_correction_volume = n;
557 variables_internes_->VOFlike_correction_volume = 1;
558 variables_internes_->nb_lissage_correction_volume = n;
559 variables_internes_->nb_iterations_correction_volume = n;
562 Cerr <<
"Obsolete Keyword is read Iterations_correction_volume : " << n << finl;
563 Cerr <<
"For future compatibility, it is recommended to switch to the new syntax : " << finl;
565 Cerr <<
" VOFlike_correction_volume 1 # a flag (0 or 1) for activation #" << finl;
567 Cerr <<
" VOFlike_correction_volume 0 # a flag (0 or 1) for activation #" << finl;
568 Cerr <<
" nb_lissage_correction_volume " << n <<
" # Select 0 or N the number of smoothing to apply to avoid potential spikes due to volume correction. #" << finl;
569 Cerr <<
" nb_iterations_correction_volume " << n <<
" # to iterate on the volume correction until seuil is reached. #" << finl;
573 Cerr <<
"The value of nb_iterations_correction_volume is read from bloc remaillage at keyword: nb_iter_correction_volume" << finl;
578 else if (un_mot==
"volume_impose_phase_1")
582 variables_internes_->volume_impose_phase_1 = x;
584 Cerr <<
"volume_impose_phase_1 : " << x << finl;
587 else if (un_mot==
"methode_interpolation_v")
590 motcles2[0] =
"valeur_a_elem";
591 motcles2[1] =
"vdf_lineaire";
592 motcles2[2] =
"mean_volumic_velocity";
596 Cerr << un_mot <<
" " << motlu << finl;
597 int rang = motcles2.
search(motlu);
601 variables_internes_->methode_interpolation_v =
605 variables_internes_->methode_interpolation_v =
609 variables_internes_->methode_interpolation_v =
611 if (interpolation_repere_local_)
613 Cerr <<
"ERROR : interpolation_repere_local should not be employed"
614 " with methode_interpolation_v::MEAN_VOLUMIC_VELOCITY !!!!!!!!!!!" << finl;
620 Cerr <<
"Transport_Interfaces_FT_Disc::lire\n"
621 <<
"The options for " << un_mot <<
" are :\n"
627 else if (un_mot==
"parcours_interface")
629 is >> variables_internes_->parcours_interface_;
632 else if (un_mot==
"injecteur_interfaces")
638 Cerr <<
"Reading of the interfaces to inject in the file " << nom_fichier << finl;
640 ArrOfDouble& temps = variables_internes_->injection_interfaces_temps_;
641 ArrOfInt& phase = variables_internes_->injection_interfaces_phase_;
642 Noms& expr = variables_internes_->injection_interfaces_expressions_;
650 temps.append_array(t);
656 Cerr <<
"Time=" << t <<
" phase=" << ph <<
" expression=" << une_expr << finl;
659 envoyer_broadcast(variables_internes_->injection_interfaces_temps_, 0);
660 envoyer_broadcast(variables_internes_->injection_interfaces_phase_, 0);
661 envoyer_broadcast(variables_internes_->injection_interfaces_expressions_, 0);
665 else if (un_mot==
"interpolation_champ_face")
668 motcles2[0] =
"base";
669 motcles2[1] =
"lineaire";
672 int rang = motcles2.
search(motbis);
676 variables_internes_->interpolation_champ_face =
681 variables_internes_->interpolation_champ_face =
685 mots.add(
"vitesse_fluide_explicite");
686 mots.add(
"extrapolation_diphasique_solide");
687 mots.add(
"extrapolation_solide");
688 mots.add(
"distance_projete_face");
693 Motcle accouverte =
"{" , accfermee =
"}" ;
694 if (mot2 == accouverte)
698 while (mot2 != accfermee)
700 int rang2 = mots.
search(mot2);
705 variables_internes_->vf_explicite = 1 ;
706 Cerr <<
"Lecture du type de vitesse fluide pour le calcul de v_imp : vf_explicite " << finl;
711 variables_internes_->is_extra_diphasique_solide = 1 ;
712 variables_internes_->is_extra_solide = 0 ;
713 Cerr <<
"Lecture du type d'extrapolation considere : diphasique_solide " << finl;
718 variables_internes_->is_extra_solide = 1 ;
719 variables_internes_->is_extra_diphasique_solide = 1 ;
720 Cerr <<
"Lecture du type d'extrapolation considere : solide " << finl;
725 variables_internes_->is_distance_projete_face = 1 ;
726 Cerr <<
"Lecture du type de distance consideree dans l'interpolation : distance reelle entre projete et face " << finl;
735 Cerr <<
"Erreur a la lecture des parametres de l'interpolation lineaire : " << finl;
736 Cerr <<
"On attendait : " << accouverte <<
" , on a eu " << mot2 << finl;
742 Cerr <<
"Transport_Interfaces_FT_Disc::lire\n"
743 <<
" les options de interpolation_champ_face sont :\n"
748 else if (un_mot==
"type_vitesse_imposee")
751 motcles2[0] =
"uniforme";
752 motcles2[1] =
"analytique";
755 int rang = motcles2.
search(motbis);
765 Cerr <<
"Transport_Interfaces_FT_Disc::lire\n"
766 <<
" les options de type_vitesse_imposee sont :\n"
771 else if (un_mot==
"distance_IBC_faces")
774 motcles2[0] =
"initiale";
775 motcles2[1] =
"modifiee";
779 int rang = motcles2.
search(motbis);
791 Cerr <<
"Transport_Interfaces_FT_Disc::lire\n"
792 <<
" les options de distance_IBC_faces sont :\n"
797 else if (un_mot==
"distance_projete_faces")
800 motcles2[0] =
"initiale";
801 motcles2[1] =
"modifiee";
802 motcles2[2] =
"simplifiee";
806 int rang = motcles2.
search(motbis);
819 Cerr <<
"Transport_Interfaces_FT_Disc::lire\n"
820 <<
" les options de distance_projete_faces sont :\n"
825 else if (un_mot==
"type_indic_faces")
828 motcles2[0] =
"standard";
829 motcles2[1] =
"modifiee";
830 motcles2[2] =
"ai_based";
833 Cerr << un_mot <<
" " << motlu << finl;
834 int rang = motcles2.
search(motlu);
841 Cerr <<
" Standard interpolation of indicatrice to faces." << finl;
848 <<
" analytical expressions vx(x,y,[z,]t) vy(...) [vz(...)]" << finl;
857 variables_internes_->modified_indic_faces_position = 0. ;
858 variables_internes_->modified_indic_faces_thickness= 1. ;
860 mots.add(
"position");
861 mots.add(
"thickness");
865 Motcle accouverte =
"{" , accfermee =
"}" ;
866 if (mot == accouverte)
870 while (mot != accfermee)
872 int rang2 = mots.
search(mot);
877 is >> variables_internes_->modified_indic_faces_position ;
878 Cerr <<
"Lecture de la position de l interface modifiee " << finl;
883 is >> variables_internes_->modified_indic_faces_thickness ;
884 if ( variables_internes_->modified_indic_faces_thickness < 0.)
886 Cerr <<
"L epaisseur de l interface doit etre positive ou nulle!!" << finl;
889 Cerr <<
"Lecture de l epaisseur de l interface modifiee " << finl;
898 Cerr <<
"Erreur a la lecture des parametres de l'indicatrice modifiee " << finl;
899 Cerr <<
"On attendait : " << accouverte << finl;
904 Cerr <<
"L indicatrice face sera calculee a partir de la distance. Position : d=" << variables_internes_->modified_indic_faces_position <<
"h ; Epaisseur : "<< variables_internes_->modified_indic_faces_thickness <<
"h" <<finl;
912 Cerr <<
" The interpolation of indicatrice to faces is based on the interfacial area"
913 <<
" and on the normal to the interface." << finl;
917 Cerr <<
"Transport_Interfaces_FT_Disc::lire\n"
918 <<
" les options de type_indic_faces sont :\n"
935 const Conds_lim& les_cl = le_dom_Cl_dis->les_conditions_limites();
936 const int n = les_cl.size();
938 for (i = 0; i < n; i++)
945 Cerr <<
"Error in Transport_Interfaces_FT_Disc::verif_Cl():\n"
946 <<
" Boundary conditions for Transport_Interfaces_FT_Disc must be\n"
947 <<
" of type Paroi_FT_disc (example : \"Paroi_FT_disc symetrie\")" << finl;
954int facettes_trier_retirer_doublons(IntTab& fa7)
956 using line = std::array<int, _DIM_>;
957 line* chunked_fa7 =
reinterpret_cast<line*
>(fa7.
addr());
960 auto same_fa7 = [&](line A, line B)
963 std::sort(A.begin(), A.end());
964 std::sort(B.begin(), B.end());
966 for (
int i = 0; i < dim; i++)
976 std::sort(chunked_fa7, chunked_fa7 + nb_facettes, [&](
const line& fa1,
const line& fa2)
978 return (same_fa7(fa1,fa2) < 0);
983 for (
int i = 1; i < nb_facettes; i++)
985 line& fa1 = chunked_fa7[i];
986 line& fa2 = chunked_fa7[i-1];
987 if(same_fa7(fa1,fa2)!=0)
990 for(
int j=0; j<dim; j++)
991 fa7(count,j)=fa7(i,j);
1003 ArrOfInt phase_specifiee;
1005 int default_phase = -1;
1008 int reverse_normal = 0;
1014 Cerr <<
"Error in Transport_Interfaces_FT_Disc::lire_cond_init, expected { after fichier_geom" << finl;
1019 motcles[0] =
"fichier_geom";
1020 motcles[1] =
"point_phase";
1021 motcles[2] =
"default_phase";
1022 motcles[3] =
"lata_dump";
1023 motcles[4] =
"nom_domaine";
1024 motcles[5] =
"reverse_normal";
1030 const int rang = motcles.
search(motlu);
1035 Cerr <<
" Reading .geom file in following file: " << filename << finl;
1042 is >> phase_specifiee[i];
1043 Cerr <<
" Reading point in phase " << phase_specifiee[i];
1044 if (phase_specifiee[i] != 0 && phase_specifiee[i] != 1)
1046 Cerr <<
" Error reading point_phase : expected 0 or 1" << finl;
1052 Cerr <<
" " << points(i,j);
1058 is >> default_phase;
1059 if (default_phase != 0 && default_phase != 1)
1061 Cerr <<
" Error reading default_phase : expected 0 or 1" << finl;
1064 Cerr <<
" Default phase : " << default_phase << finl;
1068 Cerr <<
" Dumping connex components and indicator function in lata file : "
1069 << lata_file << finl;
1073 Cerr <<
" Reading interface data in domain: " << domain_name << finl;
1077 Cerr <<
" Reverse_normal : swap nodes to reverse the normal vector of the surface mesh" << finl;
1080 Cerr <<
" Unknown keyword " << motlu
1081 <<
"\n Known keywords are " << motcles << finl;
1087 if (filename !=
"??" && domain_name !=
"??")
1089 Cerr <<
"Error: you specified both a .geom file name AND a domain name" << finl;
1092 else if (filename ==
"??" && domain_name ==
"??")
1094 Cerr <<
"Error: you must specify a FICHIER_GEOM or a NOM_DOMAINE" << finl;
1097 if (filename !=
"??")
1099 Cerr <<
"Reading geometry file and building interface mesh" << finl;
1100 if (is_a_binary_file(filename))
1114 Cerr <<
"Reading interface in existing domain: " << domain_name << finl;
1117 Cerr <<
"Error : object " << domain_name <<
" is not a domain" << finl;
1124 Cerr <<
"Reversing mesh normal vectors" << finl;
1125 Domaine& domaine = ref_dom.valeur();
1126 IntTab& elems = domaine.les_elems();
1127 const int nb_elem = elems.dimension(0);
1128 const int nb_som_elem = elems.line_size();
1129 if (nb_som_elem != 2 && nb_som_elem != 3)
1131 Cerr <<
"Error: mesh has wrong dimension (must be segments in 2D, triangles in 3D)" << finl;
1134 const int j0 = nb_som_elem - 2;
1135 for (
int i = 0; i < nb_elem; i++)
1137 const int s0 = elems(i, j0);
1138 const int s1 = elems(i, j0+1);
1140 elems(i, j0+1) = s0;
1149 Domaine& domaine = ref_dom.valeur();
1150 IntTab& fa7 = domaine.les_elems();
1153 const int nb_facettes = fa7.
dimension(0);
1158 count = facettes_trier_retirer_doublons<2>(fa7);
1160 count = facettes_trier_retirer_doublons<3>(fa7);
1163 Cerr<<
"End correction SM 12/08 the removed faces number is: "<< nb_facettes <<
" - "<< count+1 <<
" = " << nb_facettes-count-1 << finl;
1168 Cerr <<
"Building interface mesh" << finl;
1172 Cerr <<
"Extracting connex components and assigning indicator function." << finl;
1181 for (
int i = 0; i < nb_elem; i++)
1184 if (index_elem[i] >= 0)
1190 const IntTab& elem_faces = domaine_vf.
elem_faces();
1192 const int nb_local_connex_components = search_connex_components_local(elem_faces, faces_elem, num_compo);
1193 const int nb_connex_components = compute_global_connex_components(num_compo, nb_local_connex_components);
1196 Cerr <<
" found " << nb_connex_components <<
" connex components" << finl;
1198 ArrOfInt phase_of_component(nb_connex_components);
1199 phase_of_component= default_phase;
1203 const int nb_pts = phase_specifiee.
size_array();
1206 ArrOfInt elem_points;
1208 for (
int i = 0; i < nb_pts; i++)
1210 int composante_connexe = -2;
1211 if (elem_points[i] >= 0)
1212 composante_connexe = num_compo[elem_points[i]];
1215 composante_connexe = (int)
mp_max(composante_connexe);
1216 if (composante_connexe == -2)
1218 Cerr <<
"Error : point " << i <<
" is not in the domain" << finl;
1222 else if (composante_connexe == -1)
1224 Cerr <<
"Error : point " << i <<
" is in connex_component -1"
1225 <<
" (this point is in a mesh cell containing the interface)."
1226 <<
" Use the dump_lata keyword to find a point away from the interface" << finl;
1232 const int phase = phase_specifiee[i];
1233 Cerr <<
"Assigning phase " << phase <<
" to component " << composante_connexe << finl;
1234 phase_of_component[composante_connexe] = phase;
1239 DoubleTab& indic = variables_internes_->indicatrice_cache->valeurs();
1240 for (
int i = 0; i < nb_elem; i++)
1242 const int compo = num_compo[i];
1245 const int phase = phase_of_component[compo];
1252 if (lata_file !=
"??")
1254 Cerr <<
"Writing lata file" << finl;
1256 const Domaine& un_dom = domaine_vf.
domaine();
1257 constexpr double TEMPS = 0.;
1258 constexpr int FIRST_POST = 1;
1263 DoubleTab data(nb_elem);
1264 for (
int i = 0; i < nb_elem; i++)
1265 data(i) = num_compo(i);
1272 Nom nom_champ(
"connex_component");
1273 lata.
ecrire_champ(un_dom, unites, noms_compo, 1, TEMPS, nom_champ, nom_dom,
"elem",
"scalar",
1275 nom_champ =
"indicatrice";
1276 lata.
ecrire_champ(un_dom, unites, noms_compo, 1, TEMPS, nom_champ, nom_dom,
"elem",
"scalar",
1278 nom_champ =
"distance";
1279 lata.
ecrire_champ(un_dom, unites, noms_compo, 1, TEMPS, nom_champ, nom_dom,
"elem",
"scalar",
1286 if (phase_of_component.
size_array() > 0 && min_array(phase_of_component) < 0)
1288 Cerr <<
"Error: some connex components of the indicator function were not initialized\n"
1289 <<
" you must either specify more points or specify a default_phase.\n"
1290 <<
" With the lata_dump keyword, you can watch connex components and find where\n"
1291 <<
" points should be added." << finl;
1309 Cerr <<
"Reading initial condition" << finl;
1311 motcles[0] =
"fonction";
1312 motcles[1] =
"fichier_geom";
1313 motcles[2] =
"fonction_ignorer_collision";
1314 motcles[3] =
"reprise";
1319 Cerr <<
"Erreur dans Transport_Interfaces_FT_Disc::lire_cond_init\n";
1320 Cerr <<
" On attendait {\n On a trouve " << motlu << finl;
1323 const Motcle virgule =
",";
1331 int rang = motcles.
search(motlu);
1347 if (nom_phase==
"ajout_phase1")
1351 else if (nom_phase==
"ajout_phase0")
1357 Cerr <<
"Transport_Interfaces_FT::lire_cond_init :\n";
1358 Cerr <<
" The keyword " << nom_phase <<
" is not understood.\n";
1359 Cerr <<
" Allowed keywords are : ajout_phase0 or ajout_phase1" << finl;
1366 if (
probleme().reprise_effectuee())
1368 Cerr <<
" Interface not build since a restarting is expected." << finl;
1372 Cerr <<
" Interface construction : " << expression << finl;
1377 int ignorer_collision = (rang==2);
1379 variables_internes_->indicatrice_cache->valeurs(),
1381 variables_internes_->distance_interface_sommets,
1385 if (ignorer_collision)
1386 Cerr <<
"Warning: the interface is colliding with an existing interface (ignorer_collision=1)." << finl;
1389 Cerr <<
"Error: the interface is colliding with an existing interface." << finl;
1409 Cerr <<
"A .xyz file is waited to restart the calculation for the interface equation." << finl;
1417 is2 >> tmp >> version;
1423 Cerr <<
"Transport_Interfaces_FT::lire_cond_init :\n";
1424 Cerr <<
" The keyword " << motlu <<
" is not understood here.\n";
1425 Cerr <<
" Allowed keywords are :\n" << motcles << finl;
1429 Cerr<<
"lectureF "<<motlu<<finl;
1431 while (motlu==virgule);
1435 Cerr <<
"Error for method Transport_Interfaces_FT_Disc::lire_cond_init\n";
1436 Cerr <<
" A } was expected\n It has been found " << motlu << finl;
1445 assert(!probleme_base_);
1449 Cerr <<
"Error for method Transport_Interfaces_FT_Disc::associer_pb_base\n";
1450 Cerr <<
" The Transport_Interfaces_FT_Disc equation must be associated to a problem\n";
1451 Cerr <<
" of type Probleme_FT_Disc_gen" << finl;
1452 Cerr <<
" or type Pb_Thermohydraulique_Especes_QC" << finl;
1453 Cerr <<
" or type Pb_Thermohydraulique_Especes_Turbulent_QC" << finl;
1456 probleme_base_ = un_probleme;
1460 variables_internes_->set_pb_name(pb_name);
1474 variables_internes_->set_ident(
le_nom());
1482 fieldname =
"INDICATRICE";
1483 fieldname += suffix;
1486 1 , nb_valeurs_temps,
1489 indicatrice_->associer_eqn(*
this);
1494 fieldname =
"INDICATRICE_CACHE";
1495 fieldname += suffix;
1500 variables_internes_->indicatrice_cache);
1501 variables_internes_->indicatrice_cache->associer_eqn(*
this);
1503 variables_internes_->indicatrice_cache->PDI_save_type(
true);
1507 fieldname =
"INDICATRICE_FACES";
1508 fieldname += suffix;
1511 1 , nb_valeurs_temps,
1513 indicatrice_faces_);
1514 indicatrice_faces_->associer_eqn(*
this);
1517 fieldname =
"VITESSE_FILTREE";
1518 fieldname += suffix;
1523 variables_internes_->vitesse_filtree);
1524 variables_internes_->vitesse_filtree->associer_eqn(*
this);
1528 fieldname =
"FLUX_TMP";
1529 fieldname += suffix;
1534 variables_internes_->tmp_flux);
1538 fieldname =
"INDEX_ELEM";
1543 variables_internes_->index_element);
1546 fieldname =
"NELEM_PAR_DIRECTION";
1551 variables_internes_->nelem_par_direction);
1552 champs_compris_.ajoute_champ(variables_internes_->nelem_par_direction);
1554 fieldname =
"DISTANCE_INTERFACE_ELEM";
1555 fieldname += suffix;
1560 variables_internes_->distance_interface);
1561 champs_compris_.ajoute_champ(variables_internes_->distance_interface);
1563 fieldname =
"DISTANCE_INTERFACE_FACE";
1564 fieldname += suffix;
1569 variables_internes_->distance_interface_faces);
1570 champs_compris_.ajoute_champ(variables_internes_->distance_interface_faces);
1572 fieldname =
"DISTANCE_INTERFACE_FACE_COR";
1573 fieldname += suffix;
1578 variables_internes_->distance_interface_faces_corrigee);
1579 champs_compris_.ajoute_champ(variables_internes_->distance_interface_faces_corrigee);
1581 fieldname =
"DISTANCE_INTERFACE_FACE_DIF";
1582 fieldname += suffix;
1587 variables_internes_->distance_interface_faces_difference);
1588 champs_compris_.ajoute_champ(variables_internes_->distance_interface_faces_difference);
1590 fieldname =
"VITESSE_IMP_INTERP";
1591 fieldname += suffix;
1596 vitesse_imp_interp_);
1603 DoubleTab& d = variables_internes_->distance_interface_sommets;
1608 variables_internes_->distance_interface_sommets = 0.;
1610 fieldname =
"NORMALE_INTERFACE";
1611 fieldname += suffix;
1616 variables_internes_->normale_interface);
1621 fieldname =
"SURFACE_INTERFACE";
1622 fieldname += suffix;
1627 variables_internes_->surface_interface);
1634 fieldname =
"TIME_DERIVATIVE_INTERFACE";
1635 fieldname += suffix;
1654 variables_internes_->connectivite_frontieres_.associer_domaine_vf(domaine_vf);
1668 variables_internes_->algorithmes_transport_.typer(
"Algorithmes_Transport_FT_Disc");
1672 Cerr <<
"Discretizing of the boundary conditions... ";
1675 Cerr <<
" OK " << finl;
1677 le_dom_Cl_dis->associer_eqn(*
this);
1678 le_dom_Cl_dis->associer_inconnue(
inconnue());
1690 Journal() <<
"Transport_Interfaces_FT_Disc::remailler_interface " <<
le_nom() << finl;
1694 Champ_base& indicatrice = variables_internes_->indicatrice_cache.valeur();
1703 if (not(preparer_calcul_anticipated_done_))
1709 Process::Journal()<<
"Transport_Interfaces_FT_Disc::preparer_calcul does nothing because it has been anticipated in the problem"<<finl;
1718 preparer_calcul_anticipated_done_ =
true;
1722 compute_nb_particles_tot();
1737 if (collision_model_)
1742 collision_model_.valeur().preparer_calcul(domain_vdf, nb_particles_tot_,
1743 ns, *
this, schema_comm_FT);
1749 le_dom_Cl_dis->initialiser(temps);
1751 if (
probleme().reprise_effectuee())
1753 Cerr <<
"Calculation restarting : no initial meshing" << finl;
1757 Cerr <<
"Initial condition remeshing." << finl;
1764 indicatrice_->changer_temps(temps);
1765 variables_internes_->indicatrice_cache->changer_temps(temps);
1766 indicatrice_faces_->changer_temps(temps);
1767 variables_internes_->vitesse_filtree->changer_temps(temps);
1768 variables_internes_->tmp_flux->changer_temps(temps);
1769 variables_internes_->index_element->changer_temps(temps);
1770 variables_internes_->nelem_par_direction->changer_temps(temps);
1771 variables_internes_->distance_interface->changer_temps(temps);
1772 variables_internes_->distance_interface_faces->changer_temps(temps);
1773 variables_internes_->distance_interface_faces_corrigee->changer_temps(temps);
1774 variables_internes_->distance_interface_faces_difference->changer_temps(temps);
1775 vitesse_imp_interp_->changer_temps(temps);
1776 variables_internes_->normale_interface->changer_temps(temps);
1777 variables_internes_->surface_interface->changer_temps(temps);
1780 variables_internes_->set_checkpoint_fname(checkpoint_fname);
1790 Cerr <<
"Error for the method Transport_Interfaces_FT_Disc::preparer_calcul\n"
1791 <<
"The transport method has not indicated for the equation\n"
1810 variables_internes_->init_save_file();
1815 variables_internes_->close_save_file();
1866 if (tag == variables_internes_->distance_normale_cache_tag)
1868 Cerr <<
"WARNING : Unneeded call to Transport_Interfaces_FT_Disc::update_normale_distance_interface" << finl;
1869 Cerr <<
" mesh tag: " << tag <<
", cache tag: "<< variables_internes_->distance_normale_cache_tag << finl;
1873 DoubleTab& distance = variables_internes_->distance_interface->valeurs();
1874 DoubleTab& normale = variables_internes_->normale_interface->valeurs();
1878 variables_internes_->n_iterations_distance);
1879 variables_internes_->distance_normale_cache_tag = tag;
1894 if (tag == variables_internes_->indicatrice_cache_tag)
1896 Cerr <<
"WARNING : Unneeded call to Transport_Interfaces_FT_Disc::update_indicatrice. tag = " << tag << finl;
1897 Cerr <<
" indicatrice is already up to date" << finl;
1902 if (tag != variables_internes_->distance_normale_cache_tag)
1904 Cerr <<
"ERROR : in Transport_Interfaces_FT_Disc::update_indicatrice : distance_normale not up to date" << finl;
1905 Cerr <<
" Call update_indicatrice_normale_distance before or update_normale_distance_interface to do both at the same time" << finl;
1906 Cerr <<
" mesh tag: " << tag <<
", cache tag: "<< variables_internes_->distance_normale_cache_tag << finl;
1913 Cerr <<
"ERROR : in Transport_Interfaces_FT_Disc::update_indicatrice : maillage pas parcouru" << finl;
1917 DoubleVect& valeurs_indicatrice = variables_internes_->indicatrice_cache->valeurs();
1921 variables_internes_->indicatrice_cache_tag = tag;
1936 if (tag != variables_internes_->indicatrice_cache_tag)
1938 Cerr <<
"Error: indicatrice not up to date" << finl;
1952 if (tag != variables_internes_->indicatrice_cache_tag)
1954 Cerr <<
"Error: indicatrice not up to date in get_indicatrice. " << finl;
1957 return variables_internes_->indicatrice_cache.valeur();
1963 if (tag != variables_internes_->distance_normale_cache_tag)
1965 Cerr <<
"Error: indicatrice not up to date in get_distance_interface. " << finl;
1968 return variables_internes_->distance_interface.valeur();
1974 if (tag != variables_internes_->distance_normale_cache_tag)
1976 Cerr <<
"Error: indicatrice not up to date in get_normale_interface. " << finl;
1979 return variables_internes_->normale_interface.valeur();
1993void interpoler_vitesse_point_vdf(
const Champ_base& champ_vitesse,
1994 const FTd_vecteur3& coord_som,
1996 FTd_vecteur3& vitesse,
1997 const int future_or_past)
1999 const DoubleTab& valeurs_v = (bool)(future_or_past) ? champ_vitesse.
futur() : champ_vitesse.
valeurs();
2001 const IntTab& elem_faces = domaine_vf.
elem_faces();
2002 const DoubleTab& xv = domaine_vf.
xv();
2003 const DoubleTab& xp = domaine_vf.
xp();
2004 const IntTab& face_voisins = domaine_vf.
face_voisins();
2033 for (compo = 0; compo < dim; compo++)
2042 int elem_voisins[3] = { -1, -1, -1 };
2044 int indice_local_face_voisine[3] = { -1, -1, -1 };
2048 int elem_diagonal = -1;
2052 for (direction = 0; direction < dim; direction++)
2055 const double x = coord_som[direction];
2057 if (direction == compo)
2062 const int face_inf_compo = elem_faces(element, compo);
2063 const int face_sup_compo = elem_faces(element, compo + dim);
2064 const double xmin = xv(face_inf_compo, direction);
2065 const double xmax = xv(face_sup_compo, direction);
2066 a = (xmax - x) / (xmax - xmin);
2074 const double centre_elem = xp(element, direction);
2077 const int i_face_voisine = direction + ((x < centre_elem) ? 0 : dim);
2079 const int face_voisine = elem_faces(element, i_face_voisine);
2081 const int elem_voisin =
2082 face_voisins(face_voisine, 0) + face_voisins(face_voisine, 1) - element;
2083 if (elem_voisin >= 0)
2086 const double centre_voisin = xp(elem_voisin, direction);
2087 a = (centre_voisin - x) / (centre_voisin - centre_elem);
2088 elem_voisins[direction] = elem_voisin;
2089 indice_local_face_voisine[direction] = i_face_voisine;
2102 coef[direction] = a;
2107 const int direction1 = (compo + 1) % 3;
2108 const int direction2 = (compo + 2) % 3;
2109 const int elem_voisin1 = elem_voisins[direction1];
2110 const int elem_voisin2 = elem_voisins[direction2];
2112 if (elem_voisin1 >= 0 && elem_voisin2 >= 0)
2115 int i_face_voisine, face_voisine, elem_diagonal1, elem_diagonal2;
2117 i_face_voisine = indice_local_face_voisine[direction2];
2118 face_voisine = elem_faces(elem_voisin1, i_face_voisine);
2120 elem_diagonal1 = face_voisins(face_voisine, 0)
2121 + face_voisins(face_voisine, 1) - elem_voisin1;
2122 if (elem_diagonal1 >= 0)
2126 i_face_voisine = indice_local_face_voisine[direction1];
2127 face_voisine = elem_faces(elem_voisin2, i_face_voisine);
2129 elem_diagonal2 = face_voisins(face_voisine, 0)
2130 + face_voisins(face_voisine, 1) - elem_voisin2;
2131 if (elem_diagonal1 == elem_diagonal2)
2134 elem_diagonal = elem_diagonal1;
2148 if (elem_diagonal < 0)
2151 elem_voisins[direction1] = -1;
2152 elem_voisins[direction2] = -1;
2164 const int direction1 = 1 - compo;
2165 int element_voisin = elem_voisins[direction1];
2166 if (element_voisin < 0)
2167 element_voisin = element;
2168 const int f00 = elem_faces(element, compo);
2169 const int f10 = elem_faces(element, compo + dim);
2170 const int f01 = elem_faces(element_voisin, compo);
2171 const int f11 = elem_faces(element_voisin, compo + dim);
2173 double ci = coef[compo];
2175 double cj = coef[direction1];
2177 ci * cj * valeurs_v(f00)
2178 + (1.-ci) * cj * valeurs_v(f10)
2179 + ci * (1.-cj) * valeurs_v(f01)
2180 + (1.-ci) * (1.-cj) * valeurs_v(f11);
2182 && (xv(f00,0) <DMINFLOAT)
2183 && ((fabs(valeurs_v(f00)-valeurs_v(f10))>DMINFLOAT) || (fabs(valeurs_v(f01)-valeurs_v(f11))>DMINFLOAT)))
2185 Cerr <<
"In bidim_axi, when interpolating u_r within the first cell, we use the value on the symetry axis u_r(r=0)=0." << finl;
2186 Cerr <<
"We take a simple mean on that and the value at the other face. But for a divergence-free field, neglecting dv/dy, " << finl;
2187 Cerr <<
"it would be better to assume a velocity as u(x) = x_1/x * u_1 (if x!=0). GB 2020/03/05." << finl;
2188 const double x = coord_som[0];
2189 Cerr <<
"Here, the difference is "<< (1.-ci) <<
" vs. " << xv(f10, 0)/x << finl;
2190 Cerr <<
"u1= " << valeurs_v(f10) <<
" direction1 : " << direction1 << finl;
2191 Cerr <<
"interpoler_vitesse_point_vdf of Transport_Interface..cpp not exiting but interpolation is adapted" << finl;
2192 Cerr <<
"Former: " << vitesse[compo];
2193 vitesse[compo] = xv(f10, 0)/x * (
2195 + (1.-cj) * valeurs_v(f11));
2196 Cerr <<
" New: " << vitesse[compo] << finl;
2203 const int direction1 = (compo + 1) % 3;
2204 const int direction2 = (compo + 2) % 3;
2205 int elem00 = element;
2206 int elem10 = elem_voisins[direction1];
2207 int elem01 = elem_voisins[direction2];
2208 int elem11 = elem_diagonal;
2209 if (elem10 < 0) elem10 = element;
2210 if (elem01 < 0) elem01 = element;
2211 if (elem11 < 0) elem11 = element;
2213 const int f000 = elem_faces(elem00, compo);
2214 const int f100 = elem_faces(elem00, compo + dim);
2215 const int f010 = elem_faces(elem10, compo);
2216 const int f110 = elem_faces(elem10, compo + dim);
2217 const int f001 = elem_faces(elem01, compo);
2218 const int f101 = elem_faces(elem01, compo + dim);
2219 const int f011 = elem_faces(elem11, compo);
2220 const int f111 = elem_faces(elem11, compo + dim);
2222 double ci = coef[compo];
2223 double cj = coef[direction1];
2224 double ck = coef[direction2];
2227 ( ci * cj * valeurs_v(f000)
2228 + (1.-ci) * cj * valeurs_v(f100)
2229 + ci * (1.-cj) * valeurs_v(f010)
2230 + (1.-ci) * (1.-cj) * valeurs_v(f110)) * ck
2231 +( ci * cj * valeurs_v(f001)
2232 + (1.-ci) * cj * valeurs_v(f101)
2233 + ci * (1.-cj) * valeurs_v(f011)
2234 + (1.-ci) * (1.-cj) * valeurs_v(f111)) * (1. - ck);
2257void interpoler_simple_vitesse_point_vdf(
const Champ_base& champ_vitesse,
2258 const FTd_vecteur3& coord_som,
2260 FTd_vecteur3& vitesse,
2261 const int future_or_past)
2263 const DoubleTab& valeurs_v = (bool)(future_or_past) ? champ_vitesse.
futur() : champ_vitesse.
valeurs();
2265 const IntTab& elem_faces = domaine_vf.
elem_faces();
2266 const DoubleTab& xv = domaine_vf.
xv();
2288 for (compo = 0; compo < dim; compo++)
2296 int direction = compo;
2298 const double x = coord_som[direction];
2302 const int face_inf_compo = elem_faces(element, compo);
2303 const int face_sup_compo = elem_faces(element, compo + dim);
2304 const double xmin = xv(face_inf_compo, direction);
2305 const double xmax = xv(face_sup_compo, direction);
2306 coef = (xmax - x) / (xmax - xmin);
2313 vitesse[compo] = coef * valeurs_v(face_inf_compo) + (1.-coef) * valeurs_v(face_sup_compo);
2320 const DoubleVect& volumes = domaine_vf.
volumes();
2321 int nd, nb_nd = indicatrice.
size();
2322 assert(nb_nd==domaine_vf.
nb_elem());
2324 double integrale = 0.;
2326 for (nd=0 ; nd<nb_nd ; nd++)
2328 integrale += indicatrice(nd) * volumes(nd);
2329 integrale_ph0 += (1 - indicatrice(nd))*volumes(nd);
2351 DoubleTab& vitesse_noeuds,
2354 const bool la_roue_de_vitesse_a_deja_tournee)
const
2356 switch(variables_internes_->methode_interpolation_v)
2365 Champ_base& champ_filtre = variables_internes_->vitesse_filtree.valeur();
2366 champ_vitesse_interp = & champ_filtre;
2374 const DoubleTab& val_champ_vitesse = (bool)(explicit_u_NS_) ? champ_vitesse.
futur() : champ_vitesse.
valeurs();
2375 champ_filtre.
valeurs() = val_champ_vitesse;
2381 if (explicit_u_NS_ && la_roue_de_vitesse_a_deja_tournee) champ_filtre.
futur() = champ_filtre.
valeurs();
2384 Process::exit(
"Code never verified. I guess something is missing like : champ_filtre.valeurs() = val_champ_vitesse");
2389 champ_vitesse_interp = &champ_vitesse;
2393 Cerr <<
"Error for the method Transport_Interfaces_FT_Disc::calculer_vitesse_transport_interpolee\n";
2394 Cerr <<
"The interpolation from a field " << champ_vitesse.
que_suis_je();
2395 Cerr <<
"is not developped." << finl;
2402 DoubleTab& les_positions = variables_internes_->doubletab_pos;
2403 IntVect& les_elements = variables_internes_->intvect_elements;
2404 DoubleTab& les_vitesses = variables_internes_->doubletab_vitesses;
2406 const DoubleTab& pos = maillage.
sommets();
2408 const int nb_pos_tot = pos.
dimension(0);
2410 les_positions.
resize(nb_pos_tot, dim);
2411 les_elements.
resize(nb_pos_tot);
2413 int nb_positions = 0;
2414 for (i = 0; i < nb_pos_tot; i++)
2416 const int num_elem = elem[i];
2419 for (j = 0; j < dim; j++)
2420 les_positions(nb_positions, j) = pos(i, j);
2421 les_elements(nb_positions) = num_elem;
2425 les_positions.
resize(nb_positions, dim);
2426 les_elements.
resize(nb_positions);
2427 les_vitesses.
resize(nb_positions, dim);
2429 if (explicit_u_NS_ && la_roue_de_vitesse_a_deja_tournee)
2433 champ_vitesse_interp->
valeur_aux_elems(les_positions, les_elements, les_vitesses);
2436 vitesse_noeuds.
resize(nb_pos_tot, dim);
2438 for (i = 0; i < nb_pos_tot; i++)
2442 for (j = 0; j < dim; j++)
2443 vitesse_noeuds(i, j) = les_vitesses(nb_positions, j);
2455 Cerr <<
"Error for the method Transport_Interfaces_FT_Disc::calculer_vitesse_transport_interpolee\n"
2456 <<
"the interpolation VDF_LINEAIRE is valid only for a VDF discretization with a Champ_face field\n"
2457 <<
" (type for the current field: " << champ_vitesse.
que_suis_je() << finl;
2460 const DoubleTab& pos = maillage.
sommets();
2462 const int nb_pos_tot = pos.
dimension(0);
2466 FTd_vecteur3 vitesse;
2467 vitesse_noeuds.
resize(nb_pos_tot, dim);
2468 for (i = 0; i < nb_pos_tot; i++)
2470 const int element = elem[i];
2474 for (j = 0; j < dim; j++)
2475 coord[j] = pos(i,j);
2479 interpoler_vitesse_point_vdf(champ_vitesse, coord, element, vitesse, explicit_u_NS_ && la_roue_de_vitesse_a_deja_tournee);
2484 interpoler_simple_vitesse_point_vdf(champ_vitesse, coord, element, vitesse, explicit_u_NS_ && la_roue_de_vitesse_a_deja_tournee);
2486 for (j = 0; j < dim; j++)
2487 vitesse_noeuds(i,j) = vitesse[j];
2495 const DoubleTab& vertices = maillage.
sommets();
2496 const int nb_vertices_tot = vertices.
dimension(0);
2499 ArrOfInt id_number_fa7(nb_fa7);
2501 int n = search_connex_components_local_FT(maillage, id_number_fa7);
2502 const int nb_particles_tot=compute_global_connex_components_FT(maillage, id_number_fa7, n);
2503 const DoubleTab& phase_indicator_function = indicatrice_->valeurs();
2504 const ArrOfInt& sommets_elem = maillage.
sommet_elem();
2505 const IntTab& facettes = maillage.
facettes();
2506 DoubleTab Particles_velocity(nb_particles_tot,
dimension);
2507 Particles_velocity=0;
2508 DoubleVect V_compo_elem(nb_particles_tot);
2510 ArrOfDouble Particles_surfaces(nb_particles_tot);
2511 Particles_surfaces=0;
2512 Equation_base& eqn_hydraulique = variables_internes_->refequation_vitesse_transport.valeur();
2514 const DoubleTab& tab_velocity=champ_vitesse.
valeurs();
2517 const DoubleVect& volumes_maille = domain_vdf.
volumes();
2518 const IntTab& elem_faces = domain_vdf.
elem_faces();
2537 for (
int elem=0; elem<domain_vdf.
nb_elem(); elem++)
2541 if (phase_indicator_function(elem)==id_solid_phase)
2543 const int compo= particles_eulerian_id_number(elem);
2544 V_compo_elem(compo)+=volumes_maille(elem)*(1-phase_indicator_function(elem));
2545 for (
int dim=0; dim<
dimension; dim++) Particles_velocity(compo,dim)+=0.5*
2546 (tab_velocity(elem_faces(elem,dim))
2547 +tab_velocity(elem_faces(elem,dim+
dimension))
2548 )*volumes_maille(elem);
2554 DoubleVect s_vcompo;
2556 tab_divide_any_shape(Particles_velocity, s_vcompo);
2560 DoubleTab Positions_compo(nb_particles_tot_,
dimension);
2562 for (
int fa7=0; fa7<nb_fa7; fa7++)
2566 const int compo = id_number_fa7(fa7);
2567 const double s_fa7 = surface_fa7(fa7);
2568 Particles_surfaces(compo)+=s_fa7;
2571 for (
int k = 0; k < vertices.
dimension(1); k++)
2573 int som = facettes(fa7, k);
2574 Positions_compo(compo, dim) += s_fa7 * vertices(som, dim)/
dimension;
2581 DoubleVect s_scompo;
2583 tab_divide_any_shape(Positions_compo, s_scompo);
2589 IntVect id_number_lagrangian_vertices;
2591 RESIZE_OPTIONS::NOCOPY_NOINIT);
2592 id_number_lagrangian_vertices = -1;
2594 const int dim = vitesse_noeuds.
dimension(1);
2595 for (
int iface = 0; iface < nb_fa7; iface++)
2597 const int id_number = id_number_fa7[iface];
2598 for (
int j = 0; j < dim; j++)
2599 id_number_lagrangian_vertices[facettes(iface, j)] = id_number;
2606 for (
int som = 0; som < nb_vertices_tot; som++)
2608 if (sommets_elem[som] >= 0)
2610 for (
int dim = 0; dim <
dimension; dim++)
2611 vitesse_noeuds(som, dim) =
2612 Particles_velocity(id_number_lagrangian_vertices(som), dim);
2621 Cerr <<
"Transport_Interfaces_FT_Disc::calculer_vitesse_transport_interpolee\n"
2622 <<
" interpolation method not developped" << finl;
2632 DoubleTab& val_scal_noeuds,
2636 switch(variables_internes_->methode_interpolation_v)
2642 champ_scal_interp=champ_scal;
2652 DoubleTab scal_filtre_val(champ_scal.
valeurs());
2655 ref_cast(
Champ_P1NC, champ_scal).filtrer_L2(scal_filtre_val);
2658 champ_scal_interp->valeurs() = scal_filtre_val;
2667 Cerr<<
"champ_scal type ="<<champ_scal.
que_suis_je()<<finl;
2668 Cerr <<
"Erreur dans Transport_Interfaces_FT_Disc::calculer_scalaire_interpole\n";
2669 Cerr <<
" L'interpolation a partir d'un champ " << champ_scal.
que_suis_je();
2670 Cerr <<
" n'est pas codee." << finl;
2672 Cerr<<
"champ_scal type ="<<champ_scal.
que_suis_je()<<finl;
2673 Cerr <<
"Error for the method Transport_Interfaces_FT_Disc::calculer_scalaire_interpole\n";
2674 Cerr <<
" The interpolation from a field " << champ_scal.
que_suis_je();
2675 Cerr <<
" is not developped." << finl;
2682 DoubleTab& les_positions = variables_internes_->doubletab_pos;
2683 IntVect& les_elements = variables_internes_->intvect_elements;
2686 const DoubleTab& pos = maillage.
sommets();
2688 const int nb_pos_tot = pos.
dimension(0);
2690 les_positions.
resize(nb_pos_tot, dim);
2691 les_elements.
resize(nb_pos_tot);
2693 int nb_positions = 0;
2694 for (i = 0; i < nb_pos_tot; i++)
2696 const int num_elem = elem[i];
2699 for (j = 0; j < dim; j++)
2700 les_positions(nb_positions, j) = pos(i, j);
2701 les_elements(nb_positions) = num_elem;
2705 les_positions.
resize(nb_positions, dim);
2706 les_elements.
resize(nb_positions);
2708 champ_scal_interp->valeur_aux_elems(les_positions, les_elements,val_scal_noeuds);
2721 Cerr <<
"Transport_Interfaces_FT_Disc::calculer_scalaire_interpole\n"
2722 <<
" the interpolation method is not coded" << finl;
2740 Cerr<<
" On ne doit pas resoudre Transport_Interfaces_FT_Disc en implicite "<<finl;
2744static void init_parser_v_impose(
const Noms& expression_vitesse,
Parser& parser_x,
Parser& parser_y,
Parser& parser_z,
double temps)
2748 std::string sx(expression_vitesse[0]);
2756 parser_x.
setVar(
"z", 0.);
2757 parser_x.
setVar(
"t", temps);
2759 std::string sy(expression_vitesse[1]);
2767 parser_y.
setVar(
"z", 0.);
2768 parser_y.
setVar(
"t", temps);
2770 Nom unused_expr(
"0");
2771 std::string sz(dimension3
2772 ? expression_vitesse[2]
2782 parser_z.
setVar(
"z", 0.);
2783 parser_z.
setVar(
"t", temps);
2805 DoubleTab& vpoint,
const DoubleTab& rho_faces,
2806 DoubleTab& source_val,
const double temps,
2807 const double dt,
const int is_explicite,
2810 if (!(temps>temps_debut_))
2816 DoubleTab vit_imposee;
2822 const IntTab& face_voisins = mon_dom_dis.
face_voisins();
2846 if(variables_internes_->vimp_regul && is_explicite)
2849 DoubleTab vitesse(vpoint0);
2851 vitesse += inco_val ;
2853 for (
int i=0 ; i < n; i++)
2855 if (indicatrice_faces(i) > 0. )
2857 double f = indicatrice_faces(i);
2858 for (
int j = 0; j < m; j++)
2861 vit_imposee(i,j) = f*vit_imposee(i,j) + (1.-f)*vitesse(i,j);
2863 vit_imposee(i,j) = f*vit_imposee(i,j) + (1.-f)*vitesse(i,j)/rho_faces(i);
2865 vitesse_imp_interp_->valeurs()(i,j)= vit_imposee( i,j ) ;
2872 calcul_source(inco_val,vpoint0,rho_faces,source_val,vit_imposee,indicatrice_faces,
2873 is_QC,dt,is_explicite,eta);
2876 const DoubleVect& volumes_entrelaces = ref_cast(
Domaine_VF,mon_dom_dis).volumes_entrelaces();
2880 DoubleTab termes_sources_face(vpoint);
2881 termes_sources_face=0.;
2882 modifie_source(termes_sources_face,source_val,rho_faces,n,m,is_QC,volumes_entrelaces,le_solveur_masse);
2885 for (i = 0; i < n; i++)
2886 for (j = 0; j < m; j++)
2887 vpoint(i, j) += termes_sources_face(i,j);
2895 Cerr <<
"Error for the method Transport_Interfaces_FT_Disc::modifier_vpoint_pour_imposer_vit\n"
2896 <<
" The transport equation is not of type \"methode_transport vitesse_imposee\"."
2897 <<
" or \"methode_transport loi_horaire\"."
2902 Debog::verifier(
"Transport_Interfaces_FT_Disc::calculer_vitesse_imposee vpoint",vpoint);
2907 const IntTab& face_voisins)
2909 DoubleTab& indicatrice_faces = indicatrice_faces_->valeurs();
2911 for (
int i = 0; i < nfaces; i++)
2913 const int elem0 = face_voisins(i, 0);
2914 const int elem1 = face_voisins(i, 1);
2915 indicatrice_faces(i)= 0.;
2918 indicatrice_faces(i) = indicatrice(elem0);
2920 indicatrice_faces(i) += indicatrice(elem1);
2921 if (elem0 >= 0 && elem1 >= 0)
2922 indicatrice_faces(i) *= 0.5;
2923 double bmax=1.-1e-9;
2924 if(indicatrice_faces(i) <= 1.0e-9) indicatrice_faces(i) = 0. ;
2925 else if(indicatrice_faces(i) >= bmax) indicatrice_faces(i) = 1. ;
2928 switch(variables_internes_->type_indic_faces_)
2941 const DoubleVect& face_surfaces = domaine_vdf.
face_surfaces();
2943 double& position = variables_internes_->modified_indic_faces_position;
2944 double& thickness = variables_internes_->modified_indic_faces_thickness;
2945 for (
int i = 0; i < nfaces; i++)
2947 double h=volumes_entrelaces(i)/face_surfaces(i);
2948 if (dist_face(i) > (position+thickness/2.)*h || indicatrice_faces(i)==1.)
2949 indicatrice_faces(i)=1.;
2950 else if (dist_face(i) >= (position-thickness/2.)*h && thickness !=0.)
2951 indicatrice_faces(i) = (dist_face(i) - position*h)/(thickness*h) + 0.5 ;
2953 indicatrice_faces(i)=0.;
2965 const Equation_base& eqn_hydraulique = variables_internes_->refequation_vitesse_transport.valeur();
2974 const int vef = (dim == 2);
2977 Cerr <<
"Code never applied or checked in VEF. You should read the algo first and assess it!" << finl;
2982 for (
int face = 0; face < nfaces; face++)
2984 double indic_face = 0.;
2986 for (v = 0; v < 2; v++)
2988 const int elem = face_voisins(face, v);
2992 const double indic = indicatrice[elem];
2994 if (indic <=5e-3 || indic >= 1.-5e-3)
3002 const double ai= interfacial_area(elem);
3003 if (fabs(ai)>DMINFLOAT)
3008 for (
int j = 0; j < dim; j++)
3011 const double nx = normale_elements(elem, j);
3018 Cerr <<
"Never tested. To be verified. It should depend on a scalar product with the vect (xp-xv)" << finl;
3026 const int dir = orientation[face];
3027 const double nx = normale_elements(elem, dir);
3044 tmp = indicatrice(elem);
3046 double bmax=1.-1e-9;
3047 if(tmp <= 1.0e-9) tmp = 0. ;
3048 else if(tmp >= bmax) tmp = 1. ;
3056 Cerr <<
" WTF, c'est impossible" << finl;
3064 const int elem_voisin = face_voisins(face, 1-v);
3065 const double indic = indicatrice[elem_voisin];
3081 indicatrice_faces_->valeurs()(face) = indic_face;
3086 Cerr <<
"Interpolation option AI_BASED in Transport_Interfaces_FT_Disc is not available "
3087 <<
"in Transport_Interfaces_FT_Disc::calcul_indicatrice_faces if we do not "
3088 <<
"have a Navier_Stokes_FT_Disc equation..." << finl;
3094 Cerr <<
"Transport_Interfaces_FT_Disc::calcul_indicatrice_faces\n"
3095 <<
" unknown case?" << finl;
3100 Debog::verifier(
"Transport_Interfaces_FT_Disc::calcul_indicatrice_faces indicatrice_faces",indicatrice_faces);
3105 return variables_internes_->vimp_regul;
3109 return indicatrice_faces_;
3116 const IntTab& face_voisins = mon_dom_dis.
face_voisins();
3118 return indicatrice_faces_;
3130 const DoubleTab& vpoint,
3131 const DoubleTab& rho_faces,
3132 DoubleTab& source_val,
const DoubleTab& vit_imposee,
const DoubleTab& indicatrice_faces,
3135 const int is_explicite,
3139 DoubleTab terme_explicite(vpoint);
3143 if (indicatrice_faces.
dimension(0) == nfaces )
3147 terme_explicite *= dt;
3151 terme_explicite *= 0.;
3154 for (
int i = 0; i < nfaces; i++)
3157 double indic = (indicatrice_faces(i) > 0. ? 1.0 : 0.0) ;
3160 if (variables_internes_->vimp_regul && !is_explicite ) indic = (indicatrice_faces(i) ==1. ? 1.0 : 0.0);
3162 double rho_face = rho_faces(i);
3165 for (
int j = 0; j < dim; j++)
3167 double increment_inco = inco_val(i,j) + terme_explicite(i,j);
3169 tsource = (vit_imposee(i,j) - increment_inco) / dt * indic * rho_face * c;
3171 tsource = (vit_imposee(i,j) * rho_face - increment_inco) / dt * indic * c;
3173 source_val(i,j) = tsource;
3180 Cerr <<
"Erreur dans Transport_Interfaces_FT_Disc::calcul_source\n" << finl;
3181 Cerr <<
"Dimension de indicatrice_face differente de dimension de vpoint.\n" << finl;
3195 file+=
"_Force_totale_sur_";
3196 else if( type==
"force_totale" )
3198 file+=
"_Friction_totale_sur_";
3199 else if( type==
"Friction" )
3201 file+=
"_Friction_conv_diff_sur_";
3202 else if( type==
"Pressure" )
3204 file+=
"_Friction_Pression_sur_";
3205 else if ( type ==
"moment")
3207 file+=
"_Moment_total_sur_";
3208 else if ( type==
"particles_trajectory" )
3210 file+=
"_particles_trajectory_";
3211 else if ( type==
"mean_rms_particles_velocity")
3213 file+=
"_mean_rms_particles_velocity_";
3214 else if (type==
"particles_data")
3216 file+=
"_particles_data_";
3217 else if (type==
"particle_hydrodynamic_forces")
3220 file+=
"_particle_hydrodynamic_forces_";
3222 else if( type==
"Stokes_theoretical_forces" )
3225 file+=
"_Stokes_theoretical_forces_";
3227 else if (type==
"particle_hydrodynamic_forces_bed")
3230 file+=
"_particle_hydrodynamic_forces_bed_";
3232 else if (type==
"fluid_to_particle_heat_transfer")
3235 file+=
"_fluid_to_particle_heat_transfer_";
3237 else if (type==
"fluid_to_particle_heat_transfer_bed")
3240 file+=
"_fluid_to_particle_heat_transfer_bed_";
3242 else if(type==
"Friction_diff")
3244 file+=
"_Friction_diff_sur_" ;
3247 Cerr <<
"The file " << type <<
" is not understood by Transport_Interfaces_FT_Disc::ouvrir_fichier. "
3251 const int rang=files.search(type);
3266 fic << (
Nom)
"# Printing " << (type==
"moment"?
"of the drag moment exerted":
"of the drag exerted");
3267 fic <<
" by the fluid on the interface " << equation.
le_nom();
3268 fic <<
" " << (type==
"moment"?
"[N.m]":
"[N]") << finl;
3285 if (nb_compo>1) ch+=
"Fx";
3286 if (nb_compo>=2) ch+=espace+
"Fy";
3287 if (nb_compo>=3) ch+=espace+
"Fz";
3294 fic <<
"#########################################" << finl;
3295 fic <<
"# Position - Velocity - Collision force #" << finl;
3296 fic <<
"#########################################" << finl;
3297 fic <<
"# Time [s]" << finl;
3298 fic <<
"# Position of the gravity center of the particle [m] (px py pz)" << finl;
3299 fic <<
"# Velocity of the gravity center of the particle [m/s] (vx vy vz)" << finl;
3300 fic <<
"# Collision force discretized on the particle volume [N] (fcx fcy fcz)" << finl;
3302 fic <<
"# Time" << espace <<
"px py pz" << espace <<
"vx vy vz" << espace <<
"fcx fcy fcz" << finl;
3308 fic <<
"#####################################################" << finl;
3309 fic <<
"# Average velocity - Average velocity squared - RMS #" << finl;
3310 fic <<
"#####################################################" << finl;
3311 fic <<
"# Time [s]" << finl;
3312 fic <<
"# Average velocity of purely solid cells. For each purely solid cell,"
3313 " the velocity at gravity center is computed as the average velocity of"
3314 " the opposing faces weighted by its volume. [m/s] (vx_av vy_av vz_av)" << finl;
3315 fic <<
"# Average velocity squared of purely solid cells. [m^2/s^2] (vx2_av vy2_av vz2_av)" << finl;
3316 fic <<
"# Once the average velocity and the average velocity squared is known,"
3317 " the RMS is computed as sqrt(abs(vi_av^2 - vi2_av)) with i in {x,y,z}."
3318 " [-] (rmsx rmsy rmsz)" << finl;
3320 fic <<
"# Time" << espace <<
"vx_av vy_av vz_av" << espace <<
"vx2_av vy2_av vz2_av"
3321 << espace <<
"rmsx rmsy rmsz" << finl;
3327 fic <<
"############################################################" << finl;
3328 fic <<
"# Position - Velocity - Collision force - Collision number #" << finl;
3329 fic <<
"############################################################" << finl;
3330 fic <<
"# Time [s] - Total number of collision " << finl;
3331 fic <<
"# Position of the gravity center of the particle [m] (px py pz)" << finl;
3332 fic <<
"# Velocity of the gravity center of the particle [m/s] (vx vy vz)" << finl;
3333 fic <<
"# Collision force discretized on the particle volume [N] (fcx fcy fcz)" << finl;
3335 fic <<
"# Time" <<espace <<
"total collision number" << finl;
3336 fic <<
"# particle_id px py pz vx vy vz fcx fcy fcz number_of_particles_in_collision"<< finl;
3341 fic <<
"###################################" << finl;
3342 fic <<
"# Hydrodynamic force computation #" << finl;
3343 fic <<
"###################################" << finl;
3345 fic <<
"# Time [s]"<< espace <<
"Particle surface [m^2]" << espace <<
3346 "Pressure force [N] (fpx fpy fpz)" << espace <<
"Friction force [N] (ffx ffy ffz)" << finl;
3351 fic <<
"#####################################################################" << finl;
3352 fic <<
"# Hydrodynamic force computation - Stokes theoretical configuration #" << finl;
3353 fic <<
"#####################################################################" << finl;
3354 fic <<
"# Time [s]"<< finl;
3355 fic <<
"# Stokes theoretical PRESSURE FORCE computed from the integration, on the lagrangian mesh,"
3356 " of the discretized analytical solution [N] (fpx_th fpy_th fpz_th)" << finl;
3357 fic <<
"# Stokes PRESSURE FORCE computed with the developed method on the theoretical"
3358 " pressure field discretized on the eulerian mesh [N] (fpx_th_interp fpy_th_interp"
3359 " fpz_th_interp)" << finl;
3360 fic <<
"# Stokes theoretical FRICTION FORCE computed from the integration,"
3361 " on the lagrangian mesh, of the discretized analytical solution [N]"
3362 " (ffx_th ffy_th ffz_th)" << finl;
3363 fic <<
"# Stokes FRICTION FORCE computed with the developed method on the theoretical"
3364 " velocity field discretized on the eulerian mesh [N] (ffx_th_interp"
3365 " ffy_th_interp ffz_th_interp)" << finl;
3367 fic <<
"# Time" << espace <<
"fpx_th fpy_th fpz_th" << espace <<
3368 "fpx_th_interp fpy_th_interp fpz_th_interp" << espace <<
"ffx_th ffy_th ffz_th" <<
3369 espace <<
"ffx_th_interp ffy_th_interp ffz_th_interp" << finl;
3374 fic <<
"#########################################################" << finl;
3375 fic <<
"# Hydrodynamic force computation in a particle assembly #" << finl;
3376 fic <<
"#########################################################" << finl;
3378 fic <<
"# Time [s]" << espace <<
"particle_id" << espace <<
"Pressure force [N] (fpx fpy fpz)"
3379 << espace <<
"Friction force [N] (ffx ffy ffz)" << espace <<
3380 "Percentage of facets for which forces were computable" << espace <<
3381 "Average fluid velocity in P2" <<
"Percentage of purely fluid cells in P2"<< finl;
3386 fic <<
"#########################" << finl;
3387 fic <<
"# Heat flux computation #" << finl;
3388 fic <<
"#########################" << finl;
3389 fic <<
"# Time [s]" << finl;
3390 fic <<
"# Computation of the heat flux received by the particle from the surrounding fluid. [W] (phi)" << finl;
3392 fic <<
"# Time" << espace <<
"phi" << finl;
3397 fic <<
"################################################" << finl;
3398 fic <<
"# Heat flux computation in a particle assembly #" << finl;
3399 fic <<
"################################################" << finl;
3400 fic <<
"# Time [s]" << finl;
3401 fic <<
"# Computation of the heat flux received by the particle from the surrounding fluid. [W] (phi_i), where i stands for the particle number." << finl;
3402 fic <<
"# Average temperature of purely fluid cells in P2. [K] (T_i)" << finl;
3404 fic <<
"Time" << espace <<
"phi_0 T_0 ... phi_N T_N" << finl;
3411 os.
ouvrir(file,ios::app);
3414 os.
setf(ios::scientific);
3417 const int n,
const int m,
const int is_QC,
3421 for (
int face=0; face<n; face++)
3422 for (
int dim=0; dim<m; dim++)
3423 termes_sources_face(face,dim)=vol_entrelaces(face)*source_val(face,dim);
3426 un_solv_masse.
appliquer(termes_sources_face);
3430 for (
int i = 0; i < n; i++)
3432 const double rho_face = rho_faces(i);
3434 for (
int j = 0; j < m; j++)
3435 termes_sources_face(i,j) = termes_sources_face(i,j) / rho_face;
3444 const int nbdim1 = source_val.
line_size() == 1;
3451 DoubleTab termes_sources_face(source_val);
3452 DoubleTab termes_pressure_face(pressure_part);
3453 DoubleTab termes_friction_face(friction_part);
3454 DoubleTab termes_diff_face(diff_part);
3458 const DoubleVect& vol_entrelaces = ref_cast(
Domaine_VF,mon_dom_dis).volumes_entrelaces();
3460 ArrOfInt sequential_items_flags;
3464 for (
int face=0; face<n; face++)
3466 double indic = (indicatrice_faces(face) > 0. ? 1.0 : 0.0);
3467 double coef = vol_entrelaces(face)*indic;
3468 for (
int dim=0; dim<m; dim++)
3470 termes_sources_face(face,dim)=source_val(face,dim)*coef;
3471 termes_pressure_face(face,dim)=pressure_part(face,dim)*coef;
3472 termes_friction_face(face,dim)=friction_part(face,dim)*coef;
3473 termes_diff_face(face,dim)=diff_part(face,dim)*coef;
3477 if (sequential_items_flags[face])
3483 values(0,j) -= termes_sources_face(face,0);
3484 values(1,j) -= termes_pressure_face(face,0);
3485 values(2,j) -= termes_friction_face(face,0);
3486 values(3,j) -= termes_diff_face(face,0);
3492 values(0,j) -= termes_sources_face(face,j);
3493 values(1,j) -= termes_pressure_face(face,j);
3494 values(2,j) -= termes_friction_face(face,j);
3495 values(3,j) -= termes_diff_face(face,j);
3510 ouvrir_fichier(Force,
"force_totale",1,*
this);
3516 Force << espace << values(0,k);
3518 const int impr_mom = 1 ;
3521 ouvrir_fichier(Pressure,
"Pressure",impr_mom,*
this);
3527 Pressure << espace << values(1,k);
3531 ouvrir_fichier(Friction,
"Friction",impr_mom,*
this);
3537 Friction << espace << values(2,k);
3541 ouvrir_fichier(Friction2,
"Friction_diff",impr_mom,*
this);
3546 Friction2 << espace << values(3,k);
3560 ouvrir_fichier(Force,
"force",1,*
this);
3564 Force << espace << force_[k];
3571 ouvrir_fichier(Moment,
"moment",impr_mom,*
this);
3573 for(
int k=0; k<moment_.size_array(); k++)
3574 Moment << espace << moment_[k];
3582 if (nb_particles_tot_<dim_max_impr)
3585 ouvrir_fichier(Particles_trajectory,
"particles_trajectory",1,*
this);
3587 for (
int particle=0; particle<nb_particles_tot_; particle++)
3589 Particles_trajectory << espace;
3592 Particles_trajectory << espace;
3596 Particles_trajectory << espace;
3598 if (collision_model_)
3600 const DoubleTab& lagrangian_contact_forces=
3601 collision_model_.valeur().get_lagrangian_contact_forces();
3603 Particles_trajectory << espace << lagrangian_contact_forces(particle,dim);
3605 Particles_trajectory << finl;
3611 SFichier Moy_Rms_Vitesse_Particule;
3612 ouvrir_fichier(Moy_Rms_Vitesse_Particule,
"mean_rms_particles_velocity",1,*
this);
3620 for (
int particle=0; particle<nb_particles_tot_; particle++)
3622 Moy_Rms_Vitesse_Particule << particle <<
" ";
3623 Moy_Rms_Vitesse_Particule << espace;
3624 Moy_Rms_Vitesse_Particule << particles_purely_solid_mesh_volume(particle)<< espace;
3627 Moy_Rms_Vitesse_Particule << espace << mean_velocity(particle,dim);
3628 Moy_Rms_Vitesse_Particule << espace;
3631 Moy_Rms_Vitesse_Particule << espace << mean_squared_velocity(particle,dim);
3632 Moy_Rms_Vitesse_Particule << espace;
3635 Moy_Rms_Vitesse_Particule << espace << rms_velocity(particle,dim);
3637 Moy_Rms_Vitesse_Particule << finl;
3642 ouvrir_fichier(Particles_data,
"particles_data",1,*
this);
3643 if (collision_model_)
3645 const int collision_number=collision_model_.valeur().get_collision_number();
3647 collision_number<< finl;
3649 for (
int particle = 0; particle < nb_particles_tot_; particle++)
3651 Particles_data << particle <<
" ";
3658 if (collision_model_)
3660 const DoubleTab& lagrangian_contact_forces=
3661 collision_model_.valeur().get_lagrangian_contact_forces();
3662 const ArrOfDouble particles_collision_number=collision_model_.valeur().get_particles_collision_number();
3665 Particles_data << lagrangian_contact_forces(particle,d) << espace ;
3667 Particles_data << particles_collision_number(particle) << espace ;
3669 Particles_data << finl;
3679 if (nb_particles_tot_<dim_max_impr)
3682 ouvrir_fichier(particle_hydro_forces,
"particle_hydrodynamic_forces",1,mon_eq);
3685 for (
int compo=0; compo<nb_particles_tot_; compo++)
3687 particle_hydro_forces << espace << total_surface_interf(compo) << espace;
3688 for (
int dim=0; dim<
dimension; dim++) particle_hydro_forces <<
3689 espace << total_pressure_force(compo,dim);
3690 particle_hydro_forces << espace;
3691 for (
int dim=0; dim<
dimension; dim++) particle_hydro_forces <<
3692 espace << total_friction_force(compo,dim);
3693 particle_hydro_forces << espace;
3696 particle_hydro_forces << finl;
3700 SFichier particle_hydro_forces_bed;
3701 ouvrir_fichier(particle_hydro_forces_bed,
"particle_hydrodynamic_forces_bed",1,mon_eq);
3702 particle_hydro_forces_bed <<
"TIME ";
3704 particle_hydro_forces_bed << finl;
3708 for (
int particle = 0; particle < nb_particles_tot_; particle++)
3711 particle_hydro_forces_bed << particle ;
3712 particle_hydro_forces_bed << espace << total_surface_interf(particle) << espace;
3713 for (dim=0; dim<
dimension; dim++) particle_hydro_forces_bed << espace <<
3714 total_pressure_force(particle,dim);
3715 particle_hydro_forces_bed << espace;
3716 for (dim=0; dim<
dimension; dim++) particle_hydro_forces_bed << espace <<
3717 total_friction_force(particle,dim);
3718 particle_hydro_forces_bed << espace;
3719 for (dim=0; dim<
dimension; dim++) particle_hydro_forces_bed << espace <<
3720 Nb_fa7_ok_prop(particle,dim);
3721 for (dim=0; dim<
dimension; dim++) particle_hydro_forces_bed << espace <<
3722 U_P2_moy(particle,dim);
3723 particle_hydro_forces_bed << espace << Prop_indic_fluide_P2(particle);
3724 particle_hydro_forces_bed << finl;
3735 get_pressure_force_Stokes_th();
3737 get_friction_force_Stokes_th();
3738 int nb_particles_tot = total_pressure_force.
dimension(0);
3740 if (nb_particles_tot<dim_max_impr)
3743 ouvrir_fichier(Force_Particule_th,
"Stokes_theoretical_forces",1,mon_eq);
3745 for (
int particle=0; particle<nb_particles_tot; particle++)
3747 Force_Particule_th << espace;
3748 for (
int dim=0; dim<
dimension; dim++) Force_Particule_th << espace <<
3749 total_pressure_force_Stokes(particle,dim);
3750 Force_Particule_th << espace;
3751 for (
int dim=0; dim<
dimension; dim++) Force_Particule_th << espace <<
3752 total_pressure_force(particle,dim);
3753 Force_Particule_th << espace;
3754 for (
int dim=0; dim<
dimension; dim++) Force_Particule_th << espace <<
3755 total_friction_force_Stokes(particle,dim);
3756 Force_Particule_th << espace;
3757 for (
int dim=0; dim<
dimension; dim++) Force_Particule_th << espace <<
3758 total_friction_force(particle,dim);
3759 Force_Particule_th << finl;
3767 const DoubleVect& total_heat_transfer =
3769 if (nb_particles_tot_ < dim_max_impr)
3772 ouvrir_fichier(Heat_Transfer,
3773 "fluid_to_particle_heat_transfer", 1, mon_eq);
3775 for (
int particle = 0; particle < nb_particles_tot_;
3778 Heat_Transfer << espace << total_heat_transfer(particle);
3780 Heat_Transfer << finl;
3785 ouvrir_fichier(Heat_Transfer_bed,
3786 "fluid_to_particle_heat_transfer_bed", 1, mon_eq);
3788 for (
int particle = 0; particle < nb_particles_tot_;
3791 Heat_Transfer_bed << espace << total_heat_transfer(particle);
3793 Heat_Transfer_bed << finl;
3810 DoubleTab tab_critere(present);
3812 tab_critere = present;
3813 tab_critere -=passe;
3819 const DoubleTab& rho_faces,DoubleTab& source_val,
3820 const int is_explicite,
const double eta)
3836 if ( !is_explicite )
3838 for (
int i = 0; i<n; i++)
3840 double indic = (indicatrice_faces(i) > 0. ? 1.0 : 0.0);
3841 for (
int j = 0; j < m; j++)
3844 source_val(i,j) -= vpoint(i,j) * indic * rho_faces(i) * c;
3846 source_val(i,j) -= vpoint(i,j) * indic * c;
3851 DoubleTab termes_sources_face(vpoint);
3852 const DoubleVect& vol_entrelaces = ref_cast(
Domaine_VF,mon_dom_dis).volumes_entrelaces();
3854 for (
int face=0; face<n; face++)
3855 for (
int dim=0; dim<m; dim++)
3857 double indic = (indicatrice_faces(face) > 0. ? 1.0 : 0.0);
3858 termes_sources_face(face,dim)=vol_entrelaces(face)*vol_entrelaces(face)*source_val(face,dim)*indic;
3862 le_solveur_masse.
appliquer(termes_sources_face);
3871 const ArrOfDouble& centre_gravite = domaine.cg_moments();
3880 ArrOfInt sequential_items_flags;
3884 for (
int i = 0; i < n; i++)
3888 if (sequential_items_flags[i])
3894 dforce[j] = -termes_sources_face(i);
3895 force_[j] += dforce[j];
3900 dforce[j] = -termes_sources_face(i,j);
3901 force_[j] += dforce[j];
3907 xgr[j] = centre_faces(i,j) - centre_gravite[j];
3910 moment_[0] += dforce[1]*xgr[0] - dforce[0]*xgr[1];
3913 moment_[0] += dforce[2]*xgr[1] - dforce[1]*xgr[2];
3914 moment_[1] += dforce[0]*xgr[2] - dforce[2]*xgr[0];
3915 moment_[2] += dforce[1]*xgr[0] - dforce[0]*xgr[1];
3932 const double temps=le_schema_en_temps->temps_courant();
3934 const int dimension3 = (dim==3);
3935 Parser parser_x, parser_y, parser_z;
3937 const IntTab& face_voisins = domaine_vf.
face_voisins();
3938 const int nfaces = face_voisins.
dimension(0);
3939 const DoubleTab& xv = domaine_vf.
xv();
3945 vit_ibc.
resize(nfaces,1);
3947 vit_ibc.
resize(nfaces,dim);
3949 init_parser_v_impose(variables_internes_->expression_vitesse_imposee,
3950 parser_x, parser_y, parser_z,temps);
3951 for (
int i = 0; i < nfaces; i++)
3953 for (
int j = 0; j < dim; j++)
3955 const double coord = xv(i,j);
3956 parser_x.
setVar(j, coord);
3957 parser_y.
setVar(j, coord);
3959 parser_z.
setVar(j, coord);
3967 vit_ibc(i,0) = parser_x.
eval();
3970 vit_ibc(i,0) = parser_y.
eval();
3973 vit_ibc(i,0) = parser_z.
eval();
3979 vit_ibc(i,0) = parser_x.
eval();
3980 vit_ibc(i,1) = parser_y.
eval();
3982 vit_ibc(i,2) = parser_z.
eval();
3992 const DoubleTab& vitesse,
3993 const DoubleTab& vpoint,
3998 const int dimension3 = (dim==3);
4000 const IntTab& face_voisins = mon_dom_dis.
face_voisins();
4001 const int nfaces = face_voisins.
dimension(0);
4002 const DoubleTab& xv = ref_cast(
Domaine_VF, mon_dom_dis).xv();
4008 vitesse_imp.
resize(nfaces,1);
4010 vitesse_imp.
resize(nfaces,dim);
4014 ArrOfDouble coord(dim);
4015 ArrOfDouble v_imp(dim);
4016 for (
int i = 0; i < nfaces; i++)
4019 for (
int j = 0; j < dim; j++)
4021 v_imp = variables_internes_->loi_horaire_->vitesse(temps+dt,coord);
4026 for (
int j = 0; j < dim; j++)
4027 vitesse_imp(i,j) = v_imp[j];
4032 switch(variables_internes_->interpolation_champ_face)
4036 Parser parser_x, parser_y, parser_z;
4037 init_parser_v_impose(variables_internes_->expression_vitesse_imposee,
4038 parser_x, parser_y, parser_z,temps);
4039 for (
int i = 0; i < nfaces; i++)
4042 for (
int j = 0; j < dim; j++)
4045 const double coord = xv(i,j);
4046 parser_x.
setVar(j, coord);
4047 parser_y.
setVar(j, coord);
4049 parser_z.
setVar(j, coord);
4057 vitesse_imp(i,0) = parser_x.
eval();
4060 vitesse_imp(i,0) = parser_y.
eval();
4063 vitesse_imp(i,0) = parser_z.
eval();
4071 vitesse_imp(i,0) = parser_x.
eval();
4072 vitesse_imp(i,1) = parser_y.
eval();
4074 vitesse_imp(i,2) = parser_z.
eval();
4082 if ( !variables_internes_->vf_explicite )
4089 dt_loc = (1./dt_loc)*(le_schema_en_temps->facteur_securite_pas());
4092 dt_loc = std::min(dt_loc,dt);
4093 Cerr <<
"Transport_Interfaces_FT_Disc::calculer_vitesse_imposee avec dt_loc : "<<dt_loc<<finl;
4098 const int phase = 0;
4099 vitesse_imp = vpoint;
4100 vitesse_imp*= dt_loc;
4101 vitesse_imp += vitesse;
4103 DoubleTab gradient(vitesse);
4104 interpoler_vitesse_face(dist_face,phase,variables_internes_->n_iterations_interpolation_ibc,vitesse_imp,gradient,temps,dt);
4105 Debog::verifier(
"Transport_Interfaces_FT_Disc::calcul_vitesse vitesse_imp apres interp. ",vitesse_imp);
4109 Cerr <<
"Transport_Interfaces_FT_Disc::calculer_vitesse_imposee\n"
4110 <<
" interpolation lineaire non codee en VEF" << finl;
4117 Cerr <<
"Transport_Interfaces_FT_Disc::calculer_vitesse_imposee\n"
4118 <<
" methode d'interpolation non codee" << finl;
4125 Cerr <<
"Error for the method Transport_Interfaces_FT_Disc::calcul_vitesse\n"
4126 <<
" The transport equation is not of type \"methode_transport vitesse_imposee\"."
4127 <<
" or \"methode_transport loi_horaire\"."
4146 const DoubleTab& dist_elem,
4147 const DoubleTab& normale_elem,
4148 DoubleTab& dist_face)
const
4150 statistics().create_custom_counter(
"Calculer_distance_interface",3,
"FrontTracking");
4151 statistics().begin_count(
"Calculer_distance_interface",statistics().get_last_opened_counter_level()+1);
4153 static const double distance_faces_invalides = -1.e30;
4156 const IntTab& elem_faces = domaine_vf.
elem_faces();
4157 const DoubleTab& xv = domaine_vf.
xv();
4158 const DoubleTab& xp = domaine_vf.
xp();
4160 ArrOfInt ncontrib(nb_faces);
4165 double centre[3] = {0., 0., 0.};
4166 double normale[3] = {0., 0., 0.};
4170 const int nb_faces_elem = elem_faces.
line_size();
4172 for (elem = 0; elem < nb_elem_tot; elem++)
4174 const double d1 = dist_elem(elem);
4176 if (d1 < distance_faces_invalides)
4179 for (i = 0; i < dim; i++)
4181 centre[i] = xp(elem, i);
4182 normale[i] = normale_elem(elem, i);
4185 for (i = 0; i < nb_faces_elem; i++)
4187 const int face = elem_faces(elem, i);
4190 if (face < nb_faces)
4194 for (j = 0; j < dim; j++)
4196 double position_centre_face = xv(face, j);
4197 d2 += (position_centre_face - centre[j]) * normale[j];
4199 dist_face(face) += d1 + d2;
4205 const double valeur_invalide = distance_faces_invalides * 1.1;
4206#ifdef __INTEL_COMPILER
4209 for (i = 0; i < nb_faces; i++)
4211 const int n = ncontrib[i];
4215 dist_face(i) = valeur_invalide;
4218 statistics().end_count(
"Calculer_distance_interface");
4226 if (tag == variables_internes_->distance_faces_cache_tag)
4227 return variables_internes_->distance_interface_faces.valeur();
4229 variables_internes_->distance_faces_cache_tag = tag;
4230 assert(tag == variables_internes_->distance_normale_cache_tag);
4233 DoubleTab& dist_face = variables_internes_->distance_interface_faces->valeurs();
4236 return variables_internes_->distance_interface_faces.valeur();
4240inline void check(
const DoubleTab& v_imp,
int& i_face,
double v,
int v_est_initialise)
4242 if (v_est_initialise && v_imp(i_face)!= v)
4244 Cerr <<
"=====================================================================================" << finl;
4245 Cerr <<
"You defined a non-uniform function for the keyword: methode_transport vitesse_imposee" << finl;
4246 Cerr <<
"Please add:" << finl;
4247 Cerr <<
"type_vitesse_imposee analytique" << finl;
4253 v_est_initialise = 1;
4274 const DoubleTab& distance_interface_faces,
4276 const int stencil_width,
4278 DoubleTab& gradient,
4279 const double t,
const double dt)
4288 const IntTab& elem_faces = domaine_vf.
elem_faces();
4290 const int nfaces = faces_elem.
dimension(0);
4291 const int nb_elem = domaine_vf.
nb_elem() ;
4292 const int nb_elem_tot = domaine_vf.
nb_elem_tot() ;
4293 const IntVect& orientation = domaine_vdf.
orientation();
4294 const DoubleTab& xv = domaine_vf.
xv();
4295 const DoubleTab& xp = domaine_vf.
xp();
4296 assert(champ.dimension(0) == nfaces);
4297 assert(gradient.dimension(0) == nfaces);
4298 assert(distance_interface_faces.
dimension(0) == nfaces);
4303 Parser parser_x, parser_y, parser_z;
4304 init_parser_v_impose(variables_internes_->expression_vitesse_imposee,parser_x, parser_y, parser_z, t);
4306 const double invalid_test = -1.e30;
4307 const double invalid_value = -2.e30;
4308 gradient = invalid_value;
4309 const double indic_phase = (phase == 0) ? 0. : 1.;
4315 DoubleTab dist_face_cor(distance_interface_faces) ;
4317 DoubleTab typeface(champ) ;
4328 IntTab trav(nb_elem,dim) ;
4329 if( nb_elem < nb_elem_tot )
4335 for(
int i_elem=0 ; i_elem<nb_elem ; i_elem++)
4337 for(
int i=0 ; i<dim ; i++ )
4338 xe(i)=xp(i_elem,i) ;
4340 if( indicatrice(i_elem) != 0. && indicatrice(i_elem) != 1. )
4342 for(
int dir = 0 ; dir < dim ; dir++ )
4344 int a = elem_faces(i_elem, dir) ;
4345 int b = elem_faces(i_elem, dir+dim) ;
4346 const double pas = std::fabs(domaine_vdf.
distance_face(a,b,dir)) ;
4349 trav(i_elem,dir) = traverse ;
4352 else if( indicatrice(i_elem) == 0. )
4355 for(
int k=0; k<dim; k++)
4356 trav(i_elem,k) = -1 ;
4358 else if( indicatrice(i_elem) == 1. )
4361 for(
int k=0; k<dim; k++)
4362 trav(i_elem,k) = -2 ;
4374 int cpt_halfref = 0 ;
4375 DoubleTab typefacejoint ;
4379 for(
int njoint=0; njoint<nbjoints; njoint++)
4381 const Joint& joint_temp = domaine_vf.
joint(njoint);
4383 const int nb_faces = indices_faces_joint.
dimension(0);
4384 for (
int j = 0; j < nb_faces; j++)
4386 int face_de_joint = indices_faces_joint(j, 1);
4392 for( i_face = 0 ; i_face < nfaces ; i_face++ )
4394 if( indicatrice_face(i_face) == 0. || indicatrice_face(i_face) == 1. )
4396 typeface(i_face) = indicatrice_face(i_face) ;
4400 typefacejoint = typeface ;
4402 for(
int i=0 ; i<FACEJOINT.
size() ; i++ )
4404 typeface(FACEJOINT[i]) = std::max(typefacejoint(FACEJOINT[i]),typeface(FACEJOINT[i])) ;
4405 if( typefacejoint(FACEJOINT[i]) < 0 && typeface(FACEJOINT[i]) > -1 )
4410 for( i_face = 0 ; i_face < nfaces ; i_face++ )
4413 if( indicatrice_face(i_face) != 0. && indicatrice_face(i_face) != 1. && indicatrice_face(i_face) != 0.5 )
4416 int reference_trouvee = 0 ;
4418 while( i_voisin < 2 && reference_trouvee == 0)
4420 const int elem_voisin = faces_elem(i_face,i_voisin) ;
4422 if( elem_voisin > -1 &&
4423 ( indicatrice(elem_voisin) == 0. || indicatrice(elem_voisin) == 1. ) )
4425 reference_trouvee = 1 ;
4426 typeface(i_face) = indicatrice(elem_voisin) ;
4433 typefacejoint = typeface ;
4435 for(
int i=0 ; i<FACEJOINT.
size() ; i++ )
4437 typeface(FACEJOINT[i]) = std::max(typefacejoint(FACEJOINT[i]),typeface(FACEJOINT[i])) ;
4438 if( typefacejoint(FACEJOINT[i]) < 0 && typeface(FACEJOINT[i]) > -1
4439 && indicatrice_face(FACEJOINT[i]) != 0. && indicatrice_face(FACEJOINT[i]) != 1.
4440 && indicatrice_face(FACEJOINT[i]) != 0.5 )
4445 for( i_face = 0 ; i_face < nfaces ; i_face++ )
4447 if( indicatrice_face(i_face) == 0.5 )
4450 if( faces_elem(i_face,0) > -1 )
4451 element = faces_elem(i_face,0) ;
4453 element = faces_elem(i_face,1) ;
4455 if( indicatrice(element) == 0. || indicatrice(element) == 1. )
4457 typeface(i_face) = 0.5 ;
4462 typefacejoint = typeface ;
4464 for(
int i=0 ; i<FACEJOINT.
size() ; i++ )
4466 typeface(FACEJOINT[i]) = std::max(typefacejoint(FACEJOINT[i]),typeface(FACEJOINT[i])) ;
4467 if( typefacejoint(FACEJOINT[i]) < 0 && typeface(FACEJOINT[i]) > -1 && indicatrice_face(FACEJOINT[i]) == 0.5 )
4486 int cpt_rest = nfaces - cpt_ref - cpt_vref - cpt_halfref ;
4490 double proc_stop = -1. ;
4497 if( proc_stop < 0. )
4500 for( i_face = 0 ; i_face<nfaces ; i_face++)
4502 if( typeface( i_face ) < 0 )
4504 const int ori = orientation[i_face] ;
4505 const int elem0 = faces_elem(i_face,0) ;
4506 const int elem1 = faces_elem(i_face,1) ;
4509 if( elem0 > -1 && elem0 < nb_elem_tot )
4511 else if( elem1 > -1 && elem1 < nb_elem_tot )
4518 int nb = trav(elem,ori) ;
4521 int j_face = elem_faces(elem, ori) + elem_faces(elem, ori+dim) - i_face ;
4523 if( typeface(j_face) > -1 )
4526 if( fmod( nb, 2. ) > 0 )
4531 typeface(i_face) = 1. - typeface(j_face) ;
4538 typeface(i_face) = typeface(j_face) ;
4544 elem = elem0 + elem1 - elem ;
4545 if( elem > -1 && elem < nb_elem_tot )
4547 nb = trav(elem,ori) ;
4548 j_face = elem_faces(elem, ori) + elem_faces(elem, ori+dim) - i_face ;
4550 if( typeface(j_face) > -1 )
4553 if( fmod( nb, 2. ) > 0 )
4558 typeface(i_face) = 1. - typeface(j_face) ;
4565 typeface(i_face) = typeface(j_face) ;
4574 typefacejoint = typeface ;
4576 for(
int i=0 ; i<FACEJOINT.
size() ; i++ )
4578 typeface(FACEJOINT[i]) = std::max(typefacejoint(FACEJOINT[i]),typeface(FACEJOINT[i])) ;
4579 if( typefacejoint(FACEJOINT[i]) < 0 && typeface(FACEJOINT[i]) > -1
4580 && indicatrice_face(FACEJOINT[i]) != 0. && indicatrice_face(FACEJOINT[i]) != 1.
4581 && indicatrice_face(FACEJOINT[i]) != 0.5 )
4586 if( cpt_diph == cpt_rest )
4590 stop =
mp_sum( proc_stop ) ;
4594 for( i_face = 0 ; i_face < nfaces ; i_face++ )
4596 if( indicatrice_face(i_face) != 0. && indicatrice_face(i_face) != 1. )
4598 int elem_voisin = 0 ;
4603 const int ori = orientation[i_face] ;
4604 if( faces_elem(i_face,0) > -1 )
4606 j_face = elem_faces(faces_elem(i_face,0), ori) + elem_faces(faces_elem(i_face,0), ori+dim) - i_face ;
4607 dxa = std::fabs(domaine_vdf.
distance_face(i_face,j_face,ori)) ;
4609 if( faces_elem(i_face,1) > -1 )
4611 j_face = elem_faces(faces_elem(i_face,1), ori) + elem_faces(faces_elem(i_face,1), ori+dim) - i_face ;
4612 dxb = std::fabs(domaine_vdf.
distance_face(i_face,j_face,ori)) ;
4616 elem_voisin = faces_elem(i_face,1) ;
4621 elem_voisin = faces_elem(i_face,0) ;
4625 for(
int k=0 ; k<dim ; k++ )
4629 int a = elem_faces(elem_voisin, k) ;
4630 int b = elem_faces(elem_voisin, k+dim) ;
4632 Lref += 0.25*pas*pas ;
4635 const double ratio = std::fabs(distance_interface_faces(i_face))/sqrt(Lref) ;
4636 const double tol_distface = 1e-3 ;
4637 const int test1 = ( ratio < tol_distface ) ;
4639 typeface(i_face) = 0.5 ;
4642 typefacejoint = typeface ;
4644 for(
int i=0 ; i<FACEJOINT.
size() ; i++ )
4646 typeface(FACEJOINT[i]) = typefacejoint(FACEJOINT[i]) ;
4659 for( i_face=0 ; i_face<nfaces ; i_face++ )
4661 if( distance_interface_faces(i_face) > invalid_test )
4662 dist_face_cor(i_face) = (2.*typeface(i_face) - 1.)*std::fabs(distance_interface_faces(i_face)) ;
4664 dist_face_cor(i_face) = distance_interface_faces(i_face) ;
4669 for( i_face=0 ; i_face<nfaces ; i_face++ )
4671 dist_face_cor(i_face) = distance_interface_faces(i_face) ;
4678 for( i_face=0 ; i_face<nfaces ; i_face++ )
4680 variables_internes_->distance_interface_faces_corrigee->valeurs()(i_face) = dist_face_cor(i_face) ;
4681 variables_internes_->distance_interface_faces_difference->valeurs()(i_face) = dist_face_cor(i_face) - distance_interface_faces(i_face) ;
4687 int cpt_face_fluide_signefaux = 0 ;
4688 int cpt_face_solide_signefaux = 0 ;
4689 int cpt_face_diphasique_signefaux = 0 ;
4690 int cpt_face_diphasique_signemodifiee = 0 ;
4691 int cpt_face_diphasique_signemodifiee_dnulle = 0 ;
4692 int nfaces_diphasique = 0 ;
4693 int nfaces_fluide = 0 ;
4694 int nfaces_solide = 0 ;
4696 for( i_face=0 ; i_face<nfaces ; i_face++ )
4698 if( indicatrice_face( i_face ) == 0. )
4700 else if( indicatrice_face( i_face ) == 1. )
4703 nfaces_diphasique++;
4705 if( indicatrice_face( i_face ) == 0. && dist_face_cor( i_face ) >= 0. )
4707 cpt_face_fluide_signefaux++ ;
4709 else if( indicatrice_face( i_face ) == 1. && dist_face_cor( i_face ) <= 0. && dist_face_cor( i_face ) > invalid_test )
4711 cpt_face_solide_signefaux++ ;
4713 else if( indicatrice_face( i_face ) != 1. && indicatrice_face( i_face ) != 0. )
4716 if( faces_elem(i_face,0) > -1 )
4717 element = faces_elem(i_face,0) ;
4719 element = faces_elem(i_face,1) ;
4721 if( indicatrice( element ) == 0. && dist_face_cor( i_face ) >= 0. )
4723 cpt_face_diphasique_signefaux++ ;
4725 else if( indicatrice( element ) == 1. && dist_face_cor( i_face ) <= 0. && dist_face_cor( i_face ) > invalid_test )
4727 cpt_face_diphasique_signefaux++ ;
4731 if( std::fabs(dist_face_cor(i_face) - distance_interface_faces(i_face)) != 0. )
4733 cpt_face_diphasique_signemodifiee++ ;
4735 if( dist_face_cor( i_face ) == 0. )
4737 cpt_face_diphasique_signemodifiee_dnulle++ ;
4742 Cerr <<
" -----------------------------------------------------"<< finl ;
4743 Cerr <<
" -----------------------------------------------------"<< finl ;
4744 Cerr <<
" -----------------------------------------------------"<< finl ;
4745 Cerr <<
" Nombre de faces : "<< nfaces <<finl ;
4746 Cerr <<
" Nombre de faces fluides : "<<nfaces_fluide<<finl;
4747 Cerr <<
" Nombre de faces fluides avec mauvais signe : "<<cpt_face_fluide_signefaux<<finl;
4748 Cerr <<
" Nombre de faces solides : "<<nfaces_solide<<finl;
4749 Cerr <<
" Nombre de faces solides avec mauvais signe : "<<cpt_face_solide_signefaux<<finl;
4750 Cerr <<
" Nombre de faces diphasiques : "<<nfaces_diphasique<<finl;
4751 Cerr <<
" Nombre de faces diphasiques avec mauvais signe : "<<cpt_face_diphasique_signefaux<<finl;
4752 Cerr <<
" Nombre de faces diphasiques avec signe modifie : "<<cpt_face_diphasique_signemodifiee<<finl;
4753 Cerr <<
" Nombre de faces diphasiques avec signe modifie et distance nulle : "<<cpt_face_diphasique_signemodifiee_dnulle<<finl;
4754 Cerr <<
" -----------------------------------------------------"<< finl ;
4755 Cerr <<
" -----------------------------------------------------"<< finl ;
4756 Cerr <<
" -----------------------------------------------------"<< finl ;
4761 DoubleTab v_imp(champ) ;
4763 Cerr <<
" INITIALISATION DE Vimp a Vsolide pour les faces solides "<<finl ;
4764 for ( i_face = 0; i_face < nfaces; i_face++)
4766 const double d = dist_face_cor(i_face) ;
4767 if ( d>0. || (indicatrice_face(i_face) == 1. && d < -1.e20) )
4769 for (
int j = 0; j < dim; j++)
4771 const double coord = xv(i_face,j);
4772 parser_x.
setVar(j, coord) ;
4773 parser_y.
setVar(j, coord) ;
4775 parser_z.
setVar(j, coord) ;
4782 v_imp(i_face) = parser_x.
eval() ;
4785 v_imp(i_face) = parser_y.
eval() ;
4788 v_imp(i_face) = parser_z.
eval() ;
4799 Cerr <<
"Transport_Interfaces_FT_Disc_interne::UNIFORME " << finl;
4800 double vx=0,vy=0,vz=0;
4802 for ( i_face = 0; i_face < nfaces; i_face++)
4807 double d = dist_face_cor(i_face) ;
4808 if (indicatrice_face(i_face) == indic_phase && d > invalid_test)
4811 for (
int j = 0; j < dim; j++)
4813 const double coord = xv(i_face,j);
4814 parser_x.
setVar(j, coord);
4815 parser_y.
setVar(j, coord);
4817 parser_z.
setVar(j, coord);
4827 v_imp(i_face) = parser_x.
eval() ;
4828 check(v_imp,i_face,vx,ix);
4831 v_imp(i_face) = parser_y.
eval() ;
4832 check(v_imp,i_face,vy,iy);
4835 v_imp(i_face) = parser_z.
eval() ;
4836 check(v_imp,i_face,vz,iz);
4839 gradient(i_face) = (champ(i_face) - v_imp(i_face)) / d ;
4849 Cerr <<
"Transport_Interfaces_FT_Disc_interne::ANALYTIQUE "<< finl ;
4853 double indic_pos=variables_internes_->modified_indic_faces_position;
4854 double indic_thi=variables_internes_->modified_indic_faces_thickness;
4855 if (indic_pos-0.5*indic_thi < -1.)
4857 Cerr <<
"-------------------------------------------------------------------------" << finl;
4858 Cerr <<
"-------------------------------------------------------------------------" << finl;
4859 Cerr <<
"The datafile ask to use type_vitesse_imposee=ANALYTIQUE and" << finl;
4860 Cerr <<
"type_indic_faces modifiee with a too large transition to reach zero" << finl;
4861 Cerr <<
"There is a big risk that fluid gradients cannot be computed" << finl;
4862 Cerr <<
"because fluid points (Indic_faces=0) are too far from the interface" << finl;
4863 Cerr <<
"To fix that, just do one of the following things:" << finl;
4864 Cerr <<
"If the imposed velocity field is uniform, use type_vitesse_imposee=UNIFORM" << finl;
4865 Cerr <<
"If not, use \"distance_projete_faces simplifiee\" " << finl;
4866 Cerr <<
"-------------------------------------------------------------------------" << finl;
4867 Cerr <<
"-------------------------------------------------------------------------" << finl;
4875 Cerr <<
"Transport_Interfaces_FT_Disc_interne::PROJETE_SIMPLIFIE, calcul a partir des distances et normales" << finl;
4884 for (i_face = 0; i_face < nfaces; i_face++)
4886 double d = dist_face_cor(i_face) ;
4887 const int diphasique_tmp = (indicatrice_face(i_face) != 0. &&
4888 (indicatrice_face(i_face) != 1. || (indicatrice_face(i_face) == 1. && d<=0. && d> invalid_test)));
4889 const int diphasique = ( variables_internes_->is_extra_diphasique_solide ? diphasique_tmp : (diphasique_tmp && d <= 0.) ) ;
4890 const int fluide = ( indicatrice_face(i_face) == 0. ) ;
4891 const int solide = ( variables_internes_->is_extra_solide && (d>0. || (indicatrice_face(i_face) == 1. && d<-1.e20)) ) ;
4892 const int monophasique = ( fluide || solide ) ;
4894 if (diphasique || (monophasique && d> invalid_test ))
4896 DoubleTab coord_projete(dim) ;
4897 DoubleTab normale_face(dim) ;
4898 DoubleTab coord_face(dim) ;
4899 coord_projete = 0. ;
4902 double cpt_normale = 0.;
4904 for (
int i=0 ; i<dim ; i++)
4905 coord_face(i) = xv(i_face, i) ;
4908 for(
int i=0 ; i<2; i++)
4910 const int elem_voisin = faces_elem(i_face,i) ;
4911 if (elem_voisin > -1)
4913 for (
int j = 0; j < dim; j++)
4914 normale_face(j) += normale_elem(elem_voisin, j);
4918 assert(cpt_normale!=0.);
4919 normale_face /=cpt_normale;
4920 double norme_normale = 0.;
4921 for (
int i=0 ; i<dim ; i++)
4922 norme_normale +=normale_face(i)*normale_face(i);
4923 assert(norme_normale!=0.);
4924 normale_face /= sqrt(norme_normale);
4926 for (
int i=0 ; i<dim ; i++)
4927 coord_projete(i) = coord_face(i) - d*normale_face(i);
4929 for (
int j = 0; j < dim; j++)
4931 parser_x.
setVar(j, coord_projete(j)) ;
4932 parser_y.
setVar(j, coord_projete(j)) ;
4934 parser_z.
setVar(j, coord_projete(j)) ;
4941 v_imp(i_face) = parser_x.
eval() ;
4944 v_imp(i_face) = parser_y.
eval() ;
4947 v_imp(i_face) = parser_z.
eval() ;
4969 const int exec_planfa7existant = 0 ;
4985 int nb_fa7_accepted = variables_internes_->nomb_fa7_accepted;
4986 int max_fa7 = nb_fa7_accepted ;
4991 IntList OutElem, OutFa7 ;
4994 calcul_OutElemFa7(maillage,indicatrice,nb_elem,nb_fa7_accepted,OutElem,OutFa7) ;
5000 calcul_maxfa7(maillage,indicatrice,nb_elem,max_fa7,exec_planfa7existant) ;
5012 DoubleTab Tab100( nb_elem, nb_facettes_in_elem ) ;
5013 DoubleTab Tab101( nb_elem, nb_facettes_in_elem ) ;
5014 DoubleTab Tab102( nb_elem, nb_facettes_in_elem ) ;
5015 DoubleTab Tab103( nb_elem, nb_facettes_in_elem ) ;
5016 DoubleTab Tab110( nb_elem, nb_facettes_in_elem ) ;
5017 DoubleTab Tab111( nb_elem, nb_facettes_in_elem ) ;
5018 DoubleTab Tab112( nb_elem, nb_facettes_in_elem ) ;
5019 IntTab Tab12( nb_elem, nb_facettes_in_elem ) ;
5021 IntTab CptFacette( nb_elem ) ;
5025 DoubleTab Vertex(nfaces,dim) ;
5026 DoubleTab PPP(nfaces,dim) ;
5029 if( nb_elem < nb_elem_tot )
5056 DoubleTab Barycentre(nb_facettes, 3) ;
5059 Cerr <<
"STOKAGE DES FACETTES " << finl ;
5062 StockageFa7(maillage,CptFacette,Tab100,Tab101,Tab102,Tab103,
5063 Tab110,Tab111,Tab112,Tab12,Barycentre,indicatrice,
5064 OutElem,fa7,exec_planfa7existant) ;
5070 if( OutFa7.
size() != OutElem.
size()*nb_fa7_accepted )
5072 Cerr <<
"ERREUR DE DIMENSIONNEMENT " << finl ;
5073 Cerr <<
" OutFa7.size() "<<OutFa7.
size()<<finl ;
5074 Cerr <<
" OutElem.size()*nb_fa7_accepted "<<OutElem.
size()*nb_fa7_accepted<<finl ;
5080 IntTab TabOutFa7(OutElem.
size(),nb_fa7_accepted) ;
5081 for(
int i=0 ; i< OutElem.
size() ; i++)
5084 CptFacette( OutElem[i] ) = nb_fa7_accepted ;
5085 for(
int j=0 ; j< nb_fa7_accepted ; j++)
5087 TabOutFa7(i,j) = OutFa7[j+i*nb_fa7_accepted] ;
5090 StockageFa7(maillage,Tab100,Tab101,Tab102,Tab103,Tab110,Tab111,Tab112,
5091 Tab12,Barycentre,OutElem,TabOutFa7,fa7) ;
5104 if( nb_elem < nb_elem_tot )
5106 Cerr <<
" RENUMEROTATION DES FACETTES VIRTUELLES " << finl ;
5107 RenumFa7(Barycentre,Tab110,Tab111,Tab112,Tab12,CptFacette,
5108 nb_facettes,nb_facettes_dim) ;
5127 Cerr <<
" IDENTIFICATICATION DU VERTEX LE PLUS PROCHE POUR LES FACES CONCERNEES PAR LA PROJECTION " << finl ;
5128 Cerr <<
" CAS DES FACES APPARTENANT A AU MOINS UN ELEMENT TRAVERSE : ETAPE 1 " <<finl ;
5131 Cerr <<
" CAS DES FACES APPARTENANT A AU MOINS UN ELEMENT TRAVERSE : ETAPE 2 " <<finl ;
5134 Cerr <<
" CAS DES FACES MONOPHASIQUE : ETAPE 3 " <<finl ;
5137 Cerr <<
" FIN IDENTIFICATION " << finl ;
5140 int nb_proj_modif = 0 ;
5141 int nb_proj_modif_tot = 0 ;
5142 Cerr <<
"CALCUL DU PROJETE POUR LES FACES QUI APPARTIENNENT " << finl ;
5143 Cerr <<
"A UN ELEMENT TRAVERSE PAR L INTERFACE" << finl;
5145 Tab101,Tab102,Tab103,Tab12,CptFacette,v_imp,PPP,parser_x,parser_y,parser_z) ;
5148 Cerr <<
" Nombre de projete modifie pour les faces diphasiques : "<<nb_proj_modif_tot<<
" au temps "<<
schema_temps().
temps_courant()<<finl ;
5152 Cerr <<
"CALCUL DU PROJETE POUR LES FACES FLUIDE " << finl ;
5154 t,dt,Tab100,Tab101,Tab102,Tab103,Tab12,CptFacette,v_imp,PPP,parser_x,parser_y,parser_z) ;
5158 Cerr <<
" Nombre de projete modifie pour les faces fluides : "<<nb_proj_modif_tot<<
" au temps "<<
schema_temps().
temps_courant()<<finl ;
5159 Debog::verifier(
"Transport_Interfaces_FT_Disc::interpoler_vitesse_face vimp",v_imp);
5163 for (i_face = 0; i_face < nfaces; i_face++)
5165 double d = dist_face_cor(i_face) ;
5166 if (indicatrice_face(i_face) == indic_phase && d > invalid_test)
5168 gradient(i_face) = (champ(i_face) - v_imp(i_face))/d ;
5171 gradient.echange_espace_virtuel();
5172 Debog::verifier(
"Transport_Interfaces_FT_Disc::interpoler_vitesse_face gradient avant etalement",gradient);
5178 DoubleTab gradient_old;
5179 Cerr <<
"Nombre d'iterations de la methode interpoler_vitesse_face = " << stencil_width << finl;
5180 for (
int iteration = 0; iteration < stencil_width; iteration++)
5184 gradient_old = gradient;
5188 for (i_face = 0; i_face < nfaces; i_face++)
5190 const double d = dist_face_cor(i_face) ;
5191 if (indicatrice_face(i_face) != indic_phase && d > invalid_test)
5197 const int ori=orientation(i_face);
5198 const int elem0 = faces_elem(i_face, 0);
5199 const int elem1 = faces_elem(i_face, 1);
5206 for (
int direction=0; direction <
dimension; direction++)
5213 const int voisin0 = elem_faces(elem0, ori) + elem_faces(elem0, ori+
dimension) - i_face;
5214 if( (elem_faces(elem0, ori) != i_face) && (elem_faces(elem0, ori+
dimension) != i_face) )
5216 Cerr <<
"Attention, il y a un probleme sur la numerotation des faces. On calcule en effet" << finl ;
5217 Cerr <<
"elem_faces(elem0, ori) = " <<elem_faces(elem0, ori)<< finl ;
5218 Cerr <<
"elem_faces(elem0, ori+dimension) = " <<elem_faces(elem0, ori+
dimension)<< finl ;
5219 Cerr <<
"i_face = " <<i_face<< finl;
5220 Cerr <<
"tandis qu'au moins un des deux calculs 'elem_faces(elem0, ori)' ou 'elem_faces(elem0, ori+dimension)'"<<finl ;
5221 Cerr <<
"doit correspondre au nuemro de la face i_face"<<finl ;
5222 Cerr <<
"Attention, verifier les conditions aux limites."<<finl ;
5223 Cerr <<
"En effet, le codage suivant ne permet pas de prendre en compte des conditions" << finl ;
5224 Cerr <<
"de periodicite avec une vitesse solide imposee analytique." << finl ;
5227 const double grad0 = gradient_old(voisin0);
5228 if (grad0 > invalid_test)
5236 const int voisin1 = elem_faces(elem1, ori) + elem_faces(elem1, ori+
dimension) - i_face;
5237 if( (elem_faces(elem1, ori) != i_face) && (elem_faces(elem1, ori+
dimension) != i_face) )
5239 Cerr <<
"Attention, il y a un probleme sur la numerotation des faces. On calcule en effet" << finl ;
5240 Cerr <<
"elem_faces(elem1, ori) = " <<elem_faces(elem1, ori)<< finl ;
5241 Cerr <<
"elem_faces(elem1, ori+dimension) = " <<elem_faces(elem1, ori+
dimension)<< finl ;
5242 Cerr <<
"i_face = " <<i_face<< finl;
5243 Cerr <<
"tandis qu'au moins un des deux calculs 'elem_faces(elem1, ori)' ou 'elem_faces(elem1, ori+dimension)'"<<finl ;
5244 Cerr <<
"doit correspondre au nuemro de la face i_face"<<finl ;
5245 Cerr <<
"Attention, verifier les conditions aux limites."<<finl ;
5246 Cerr <<
"En effet, le codage suivant ne permet pas de prendre en compte des conditions" << finl ;
5247 Cerr <<
"de periodicite avec une vitesse solide imposee analytique." << finl ;
5250 const double grad1 = gradient_old(voisin1);
5251 if (grad1 > invalid_test)
5265 const int face0=elem_faces(elem,direction);
5266 const int face1=elem_faces(elem,direction+
dimension);
5267 const int elem_voisin0=faces_elem(face0,0)+faces_elem(face0,1)-elem;
5268 const int elem_voisin1=faces_elem(face1,0)+faces_elem(face1,1)-elem;
5272 double grad_voisin0=0.;
5275 double grad_voisin1=0.;
5277 if (elem_voisin0 >= 0)
5279 for (
int compo = 0; compo<
dimension; compo++)
5281 double c_gravite=xv(i_face,compo);
5282 double c_face_voisin1=xv(elem_faces(elem_voisin0,ori),compo);
5283 double c_face_voisin2=xv(elem_faces(elem_voisin0,ori+
dimension),compo);
5284 a_carre+=(c_gravite-c_face_voisin1)*(c_gravite-c_face_voisin1);
5285 b_carre+=(c_gravite-c_face_voisin2)*(c_gravite-c_face_voisin2);
5287 if (a_carre<b_carre)
5288 grad_voisin0=gradient_old(elem_faces(elem_voisin0,ori));
5290 grad_voisin0=gradient_old(elem_faces(elem_voisin0,ori+
dimension));
5291 if (grad_voisin0 > invalid_test)
5293 somme += grad_voisin0;
5297 if (elem_voisin1 >= 0)
5299 for (
int compo = 0; compo<
dimension; compo++)
5301 double c_gravite=xv(i_face,compo);
5302 double c_face_voisin3=xv(elem_faces(elem_voisin1,ori),compo);
5303 double c_face_voisin4=xv(elem_faces(elem_voisin1,ori+
dimension),compo);
5304 c_carre+=(c_gravite-c_face_voisin3)*(c_gravite-c_face_voisin3);
5305 d_carre+=(c_gravite-c_face_voisin4)*(c_gravite-c_face_voisin4);
5307 if (c_carre<d_carre)
5308 grad_voisin1=gradient_old(elem_faces(elem_voisin1,ori));
5310 grad_voisin1=gradient_old(elem_faces(elem_voisin1,ori+
dimension));
5311 if (grad_voisin1 > invalid_test)
5313 somme += grad_voisin1;
5321 gradient(i_face) = somme / coeff;
5324 gradient.echange_espace_virtuel();
5327 Debog::verifier(
"Transport_Interfaces_FT_Disc::interpoler_vitesse_face gradient",gradient);
5330 for (i_face = 0; i_face < nfaces; i_face++)
5332 if (indicatrice_face(i_face) != indic_phase)
5334 double d = dist_face_cor(i_face) ;
5335 const double grad = gradient[i_face];
5336 const double indic = indicatrice_face(i_face) ;
5338 const int interpolation = ( d < 0. ) ;
5339 const int extra_diphasique = ( variables_internes_->is_extra_diphasique_solide && indic != 0. && indic != 1. && d > 0. ) ;
5340 const int extra_solide = ( variables_internes_->is_extra_solide && indic == 1. && (d>0. || d < -1.e20) ) ;
5341 const int extrapolation = ( extra_diphasique || extra_solide ) ;
5345 for (
int j = 0; j < dim; j++)
5347 const double coord = xv(i_face,j);
5348 parser_x.
setVar(j, coord);
5349 parser_y.
setVar(j, coord);
5351 parser_z.
setVar(j, coord);
5356 v_imp(i_face) = 0. ;
5360 v_imp(i_face) = parser_x.
eval();
5363 v_imp(i_face) = parser_y.
eval();
5366 v_imp(i_face) = parser_z.
eval();
5369 champ( i_face ) = v_imp(i_face) ;
5372 if (grad > invalid_test && d > invalid_test && ( interpolation || extrapolation ) )
5374 champ(i_face) += d*grad ;
5383 champ(i_face) = v_imp(i_face) ;
5384 if (grad > invalid_test && d > invalid_test && ( interpolation || extrapolation ) )
5386 champ(i_face) += d * grad ;
5389 vitesse_imp_interp_->valeurs()(i_face)= champ( i_face ) ;
5393 vitesse_imp_interp_->valeurs()(i_face)=-3e30;
5396 champ.echange_espace_virtuel() ;
5397 Debog::verifier(
"Transport_Interfaces_FT_Disc::interpoler_vitesse_face champ",champ );
5403 const int dim,
const int ori,
5409 const ArrOfInt& index_elem = intersection.
index_elem() ;
5411 const IntTab& facettes = maillage.
facettes() ;
5413 const DoubleTab& sommets = maillage.
sommets() ;
5416 int index = index_elem[elem] ;
5418 const double tol_theta = 1e-4 ;
5419 double I_sur_arete_fa7 = 0. ;
5420 double I_sur_arete_vertex = 0. ;
5421 int I_dans_fa7 = 0 ;
5431 if( std::fabs( normale_facettes(i_facette,ori) ) > 0. )
5433 DoubleTab A(dim), B(dim), C(dim), U(dim), I(dim);
5439 for(
int i = 0 ; i<dim ; i++)
5444 for(
int i=0 ; i<dim ; i++)
5446 U(i) = normale_facettes(i_facette,i) ;
5447 A(i) = sommets(facettes(i_facette, 0),i) ;
5448 B(i) = sommets(facettes(i_facette, 1),i) ;
5450 C(i) = sommets(facettes(i_facette, 2),i) ;
5453 for(
int i =0 ; i<dim ; i++)
5456 I(ori) += (A(i) - I(i))*U(i)/U(ori) ;
5461 double FIscalU = (I(ori)-xe(ori))*U(ori) ;
5462 double FI = std::fabs( I(ori)-xe(ori) ) ;
5463 double theta = acos(FIscalU/FI) ;
5464 double theta_rel = std::fabs( theta - (2.0*atan(1.)) ) / (2.0*atan(1.)) ;
5465 double ecart_elem = std::fabs(xe(ori)-I(ori)) ;
5467 if( theta_rel >= tol_theta && ecart_elem <= 0.5*dx )
5471 double AIscalAB = (I(0)-A(0))*(B(0)-A(0)) + (I(1)-A(1))*(B(1)-A(1)) ;
5472 double ABscalAB = (B(0)-A(0))*(B(0)-A(0)) + (B(1)-A(1))*(B(1)-A(1)) ;
5473 double xx = AIscalAB / ABscalAB ;
5475 const int test1 = ( std::fabs( xx ) <= precision ) ;
5476 const int test2 = ( std::fabs( xx - 1. ) <= precision ) ;
5477 const int test3 = ( xx > precision ) ;
5478 const int test4 = ( xx < 1.-precision ) ;
5479 if( test1 || test2 )
5482 I_sur_arete_vertex += 1./2. ;
5484 else if( test3 && test4 )
5492 DoubleTab ABvectAC(dim) ;
5493 ABvectAC(0) = (B(1)-A(1))*(C(2)-A(2))-(B(2)-A(2))*(C(1)-A(1)) ;
5494 ABvectAC(1) = (B(2)-A(2))*(C(0)-A(0))-(B(0)-A(0))*(C(2)-A(2)) ;
5495 ABvectAC(2) = (B(0)-A(0))*(C(1)-A(1))-(B(1)-A(1))*(C(0)-A(0)) ;
5496 double ABvectACscalU = ABvectAC(0)*U(0) + ABvectAC(1)*U(1) + ABvectAC(2)*U(2) ;
5497 double AireABC = ABvectACscalU / 2. ;
5499 DoubleTab IBvectIC(dim) ;
5500 IBvectIC(0) = (B(1)-I(1))*(C(2)-I(2))-(B(2)-I(2))*(C(1)-I(1)) ;
5501 IBvectIC(1) = (B(2)-I(2))*(C(0)-I(0))-(B(0)-I(0))*(C(2)-I(2)) ;
5502 IBvectIC(2) = (B(0)-I(0))*(C(1)-I(1))-(B(1)-I(1))*(C(0)-I(0)) ;
5503 double IBvectICscalU = IBvectIC(0)*U(0) + IBvectIC(1)*U(1) + IBvectIC(2)*U(2) ;
5504 double AireIBC = IBvectICscalU / 2. ;
5506 DoubleTab ICvectIA(dim) ;
5507 ICvectIA(0) = (C(1)-I(1))*(A(2)-I(2))-(C(2)-I(2))*(A(1)-I(1)) ;
5508 ICvectIA(1) = (C(2)-I(2))*(A(0)-I(0))-(C(0)-I(0))*(A(2)-I(2)) ;
5509 ICvectIA(2) = (C(0)-I(0))*(A(1)-I(1))-(C(1)-I(1))*(A(0)-I(0)) ;
5510 double ICvectIAscalU = ICvectIA(0)*U(0) + ICvectIA(1)*U(1) + ICvectIA(2)*U(2) ;
5511 double AireICA = ICvectIAscalU / 2. ;
5513 DoubleTab IAvectIB(dim) ;
5514 IAvectIB(0) = (A(1)-I(1))*(B(2)-I(2))-(A(2)-I(2))*(B(1)-I(1)) ;
5515 IAvectIB(1) = (A(2)-I(2))*(B(0)-I(0))-(A(0)-I(0))*(B(2)-I(2)) ;
5516 IAvectIB(2) = (A(0)-I(0))*(B(1)-I(1))-(A(1)-I(1))*(B(0)-I(0)) ;
5517 double IAvectIBscalU = IAvectIB(0)*U(0) + IAvectIB(1)*U(1) + IAvectIB(2)*U(2) ;
5518 double AireIAB = IAvectIBscalU / 2. ;
5520 double alpha = AireIBC/AireABC ;
5521 double beta = AireICA/AireABC ;
5522 double gamma = AireIAB/AireABC ;
5524 int test1 = ( std::fabs(alpha-1.)<=precision && std::fabs(beta)<=precision && std::fabs(gamma)<=precision ) ;
5525 int test2 = ( std::fabs(beta-1.)<=precision && std::fabs(alpha)<=precision && std::fabs(gamma)<=precision ) ;
5526 int test3 = ( std::fabs(gamma-1.)<=precision && std::fabs(alpha)<=precision && std::fabs(beta)<=precision ) ;
5527 int test4 = ( std::fabs(alpha)<=precision && precision<beta && beta<1.-precision && precision<gamma && gamma<1.-precision ) ;
5528 int test5 = ( std::fabs(beta)<=precision && precision<alpha && alpha<1.-precision && precision<gamma && gamma<1.-precision ) ;
5529 int test6 = ( std::fabs(gamma)<=precision && precision<beta && beta<1.-precision && precision<alpha && alpha<1.-precision ) ;
5530 int test7 = ( precision<alpha && alpha<1.-precision) ;
5531 int test8 = ( precision<beta && beta<1.-precision) ;
5532 int test9 = ( precision<gamma && gamma<1.-precision) ;
5534 if( test1 || test2 || test3 )
5538 I_sur_arete_vertex += 1e-5 ;
5540 else if( test4 || test5 || test6 )
5543 I_sur_arete_fa7 += 1./2. ;
5545 else if( test7 && test8 && test9 )
5553 Cerr <<
" attention seules les dimension 2 et 3 sont traitees "<<finl ;
5561 traverse = int(ceil(I_sur_arete_vertex)) + int(ceil(I_sur_arete_fa7)) + I_dans_fa7 ;
5567 const int nb_elem,
int& nb_fa7_accepted,
5568 IntList& OutElem, IntList& OutFa7 )
5571 const ArrOfInt& index_elem = intersection.
index_elem() ;
5576 for (
int i_elem = 0; i_elem < nb_elem ; i_elem++)
5579 if( indicatrice(i_elem) != 0. && indicatrice(i_elem) != 1. )
5581 int index = index_elem[i_elem] ;
5595 E.
add(surfaces[i_facette]) ;
5604 if( F.
size() > nb_fa7_accepted )
5613 while( i<nb_fa7_accepted )
5616 double max_data = -1e30 ;
5621 if( E[l]>max_data && !Fa7Elem.
contient(F[l]) )
5630 Fa7Elem.
add( F[ll] ) ;
5631 OutFa7.
add( F[ll] ) ;
5641 const DoubleTab& indicatrice_face, DoubleTab& Vertex )
5648 const int nb_elem = domaine_vf.
nb_elem() ;
5649 const int nfaces = faces_elem.
dimension(0) ;
5651 const ArrOfInt& index_elem = intersection.
index_elem() ;
5652 const IntTab& facettes = maillage.
facettes() ;
5653 const DoubleTabFT& sommets = maillage.
sommets() ;
5654 const DoubleTab& xv = domaine_vf.
xv();
5655 DoubleTab coord(dim) ;
5658 for(
int i_face=0 ; i_face<nfaces ; i_face++ )
5661 double dist_fa7 = 1e+30 ;
5662 int ppp_calcule = 0 ;
5665 if( indicatrice_face(i_face) == 0.5 )
5668 if( faces_elem(i_face,0) > -1 )
5669 voisin = faces_elem(i_face,0) ;
5671 voisin = faces_elem(i_face,1) ;
5675 if( indicatrice(voisin) == 0. || indicatrice(voisin) == 1. )
5677 for(
int j=0; j<dim; j++)
5679 Vertex(i_face,j) = xv(i_face,j) ;
5685 if (indicatrice_face(i_face) != 0. && indicatrice_face(i_face) != 1. && ppp_calcule == 0 )
5688 for(
int i=0 ; i<2; i++)
5690 const int elem_voisin = faces_elem(i_face,i) ;
5692 if( elem_voisin<nb_elem && elem_voisin>-1 && indicatrice(elem_voisin) != 0. && indicatrice(elem_voisin) != 1. )
5694 int index = index_elem[elem_voisin] ;
5701 double dist_min = 1e+30 ;
5703 for(
int som=0; som<facettes.
dimension(1); som++)
5706 double dist_som = 0. ;
5707 for(
int j=0; j<dim; j++)
5710 double sj = sommets(facettes(i_facette, som), j) ;
5712 double xj = xv(i_face, j) ;
5714 dist_som += (xj - sj)*(xj - sj) ;
5717 dist_som = sqrt(dist_som) ;
5720 if( dist_som < dist_min )
5722 dist_min = dist_som ;
5723 for(
int j=0; j<dim; j++)
5725 coord(j) = sommets(facettes(i_facette, som), j) ;
5730 if( dist_min < dist_fa7 )
5732 dist_fa7 = dist_min ;
5733 for(
int j=0; j<dim; j++)
5735 Vertex(i_face,j) = coord(j) ;
5749 DoubleTab& Vertex, DoubleTab& PPP )
5756 const IntTab& elem_faces = domaine_vf.
elem_faces();
5757 const int nfaces = faces_elem.
dimension(0) ;
5758 const DoubleTab& xv = domaine_vf.
xv();
5760 for(
int i_face=0 ; i_face<nfaces ; i_face++ )
5764 int ppp_calcule = 0 ;
5767 if( indicatrice_face(i_face) == 0.5 )
5770 if( faces_elem(i_face,0) > -1 )
5771 voisin = faces_elem(i_face,0) ;
5773 voisin = faces_elem(i_face,1) ;
5778 if( indicatrice(voisin) == 0. || indicatrice(voisin) == 1. )
5780 for(
int j=0; j<dim; j++)
5782 PPP(i_face,j) = xv(i_face,j) ;
5788 if (indicatrice_face(i_face) != 0. && indicatrice_face(i_face) != 1. && ppp_calcule == 0 )
5790 for(
int i=0 ; i<dim ; i++)
5792 d += (Vertex(i_face,i)-xv(i_face,i))*(Vertex(i_face,i)-xv(i_face,i)) ;
5796 for(
int k=0 ; k<dim ; k++)
5798 PPP(i_face,k) = Vertex(i_face,k) ;
5803 for(
int i=0; i<2; i++ )
5805 const int elem = faces_elem(i_face,i) ;
5810 for (
int face_loc=0 ; face_loc<2*dim ; face_loc++)
5812 const int autre_face = elem_faces(elem,face_loc) ;
5813 if (indicatrice_face(autre_face) != 0. && (indicatrice_face(autre_face) != 1.
5814 || (indicatrice_face(autre_face) == 1. && d<0.)))
5817 for(
int j=0 ; j<dim ; j++)
5819 dd += (Vertex(autre_face,j)-xv(i_face,j))*(Vertex(autre_face,j)-xv(i_face,j)) ;
5825 for(
int k=0 ; k<dim ; k++)
5827 PPP(i_face,k) = Vertex(autre_face,k) ;
5846 const IntTab& elem_faces = domaine_vf.
elem_faces();
5847 const int nfaces = faces_elem.
dimension(0) ;
5848 const DoubleTab& xv = domaine_vf.
xv();
5850 for(
int i_face=0 ; i_face<nfaces ; i_face++ )
5853 double dmin = 1.e30 ;
5855 if (indicatrice_face(i_face) == 0. || indicatrice_face(i_face) == 1.)
5860 for(
int ii=0; ii<2; ii++ )
5862 const int elem = faces_elem(i_face,ii) ;
5867 for (
int face_loc=0 ; face_loc<2*dim ; face_loc++)
5869 const int face = elem_faces(elem,face_loc) ;
5871 if (indicatrice_face(face) != 0. && indicatrice_face(face) != 1.)
5874 for(
int i=0 ; i<2; i++)
5876 const int elem_voisin = faces_elem(face,i) ;
5880 if( elem_voisin>-1 && indicatrice(elem_voisin) != 0. && indicatrice(elem_voisin) != 1. )
5883 for (
int autre_face_loc=0 ; autre_face_loc<2*dim ; autre_face_loc++)
5885 const int autre_face = elem_faces(elem_voisin,autre_face_loc) ;
5887 for(
int j=0 ; j<dim ; j++)
5889 dd += (PPP(autre_face,j)-xv(i_face,j))*(PPP(autre_face,j)-xv(i_face,j)) ;
5895 for(
int k=0 ; k<dim ; k++)
5897 PPP(i_face,k) = PPP(autre_face,k) ;
5913 const DoubleTab& indicatrice,
5914 const int nb_elem,
int& max_fa7,
5915 const int exec_planfa7existant)
5918 const ArrOfInt& index_elem = intersection.
index_elem() ;
5922 for (
int i_elem = 0; i_elem < nb_elem ; i_elem++)
5925 if( indicatrice(i_elem) != 0. && indicatrice(i_elem) != 1. )
5927 int index = index_elem[i_elem] ;
5930 DoubleList A, B, C, D ;
5938 int test_liste = 0 ;
5941 if( !(A.
est_vide()) && exec_planfa7existant == 1 )
5948 if( exec_planfa7existant == 1 )
5950 double aa=0., bb=0., cc=0., dd=0. ;
5973 DoubleTab& Tab111, DoubleTab& Tab112,
5974 IntTab& Tab12, IntTab& CptFacette,
5975 const int nb_facettes,
5976 const int nb_facettes_dim )
5983 for (
int i_elem = nn; i_elem < nn_tot ; i_elem++)
5986 for(
int cpt=0 ; cpt < CptFacette(i_elem) ; cpt++ )
5988 int bary_trouve = 0 ;
5990 while( ii < nb_facettes && bary_trouve==0 )
5992 bool test_0 = (Barycentre(ii,0)-Tab110(i_elem,cpt))*(Barycentre(ii,0)-Tab110(i_elem,cpt))<eps ;
5993 bool test_1 = (Barycentre(ii,1)-Tab111(i_elem,cpt))*(Barycentre(ii,1)-Tab111(i_elem,cpt))<eps ;
5994 bool test_2 = (Barycentre(ii,2)-Tab112(i_elem,cpt))*(Barycentre(ii,2)-Tab112(i_elem,cpt))<eps ;
5996 if( test_0 && test_1 && test_2 )
5999 Tab12( i_elem, cpt ) = ii ;
6003 if( bary_trouve==0 && ii==nb_facettes )
6008 if( Tab12( i_elem, cpt ) < nb_facettes )
6009 Tab12( i_elem, cpt )+=nb_facettes_dim ;
6016 IntTab& CptFacette, DoubleTab& Tab100,
6017 DoubleTab& Tab101, DoubleTab& Tab102,
6018 DoubleTab& Tab103, DoubleTab& Tab110,
6019 DoubleTab& Tab111, DoubleTab& Tab112,
6020 IntTab& Tab12, DoubleTab& Barycentre,
6021 const DoubleTab& indicatrice,
6022 IntList& OutElem, ArrOfBit& fa7,
6023 const int exec_planfa7existant )
6026 const ArrOfInt& index_elem = intersection.
index_elem() ;
6032 for (
int i_elem = 0; i_elem < nn ; i_elem++)
6034 if( indicatrice(i_elem) != 0. && indicatrice(i_elem) != 1. && !OutElem.
contient(i_elem) )
6037 int index = index_elem[i_elem] ;
6038 CptFacette(i_elem) = 0 ;
6041 DoubleList A, B, C, D ;
6048 int test_liste = 0 ;
6049 if( !(A.
est_vide()) && exec_planfa7existant == 1 )
6058 Tab102(i_elem,cpt),Tab103(i_elem,cpt)) ;
6059 if( exec_planfa7existant == 1 )
6061 A.
add(Tab100(i_elem,cpt)) ;
6062 B.
add(Tab101(i_elem,cpt)) ;
6063 C.
add(Tab102(i_elem,cpt)) ;
6064 D.
add(Tab103(i_elem,cpt)) ;
6066 Tab12(i_elem,cpt) = i_facette ;
6074 BaryFa7(maillage,i_facette,Barycentre) ;
6077 Tab110(i_elem,cpt) = Barycentre(i_facette,0) ;
6078 Tab111(i_elem,cpt) = Barycentre(i_facette,1) ;
6079 Tab112(i_elem,cpt) = Barycentre(i_facette,2) ;
6086 CptFacette(i_elem) = cpt ;
6089 Cerr<<
"Attention, seulement "<< cpt <<
" facettes ont ete stokees "<< finl ;
6090 Cerr<<
"dans l'element "<< i_elem <<
" qui en contient plus."<< finl ;
6091 Cerr<<
"Regarder la finesse du maillage Lagrangien par rapport au maillage Eulerien" <<finl ;
6099 DoubleTab& Tab101, DoubleTab& Tab102, DoubleTab& Tab103,
6100 DoubleTab& Tab110,DoubleTab& Tab111, DoubleTab& Tab112,
6101 IntTab& Tab12, DoubleTab& Barycentre, IntList& OutElem,
6102 IntTab& TabOutFa7, ArrOfBit& fa7 )
6107 while( cpt_elem < OutElem.
size() )
6109 for(
int k=0 ; k<TabOutFa7.
line_size() ; k++ )
6111 const int i_facette = TabOutFa7(cpt_elem,k) ;
6113 Tab101(OutElem[cpt_elem],k),Tab102(OutElem[cpt_elem],k),
6114 Tab103(OutElem[cpt_elem],k)) ;
6115 Tab12(OutElem[cpt_elem],k) = i_facette ;
6123 BaryFa7(maillage,i_facette,Barycentre) ;
6126 Tab110(OutElem[cpt_elem],k) = Barycentre(i_facette,0) ;
6127 Tab111(OutElem[cpt_elem],k) = Barycentre(i_facette,1) ;
6128 Tab112(OutElem[cpt_elem],k) = Barycentre(i_facette,2) ;
6138 const DoubleTabFT& sommets = maillage.
sommets() ;
6139 const IntTab& facettes = maillage.
facettes() ;
6142 DoubleTab CoordSom( facettes.
line_size(), dim ) ;
6145 for (
int i=0 ; i<facettes.
line_size() ; i ++ )
6147 Bary(i_facette,i) = 0. ;
6148 Som(i) = facettes(i_facette, i) ;
6149 for (
int j = 0 ; j < dim ; j++ )
6151 CoordSom(i,j) = sommets(Som(i),j) ;
6155 for (
int k = 0 ; k< dim ; k++ )
6157 for(
int l = 0 ; l< facettes.
line_size() ; l++ )
6159 Bary(i_facette,k) += CoordSom(l,k) / double(facettes.
line_size()) ;
6163 Bary(i_facette,2) = 0. ;
6168 DoubleList B,DoubleList C,DoubleList D,
6169 const int i_facette,
int& test_liste )
6174 const DoubleTabFT& sommets = maillage.
sommets() ;
6175 const IntTab& facettes = maillage.
facettes() ;
6176 double x1, y1, z1 = 0. ;
6177 double x2, y2, z2 = 0. ;
6178 double x3 = 0., y3 = 0., z3 = 0. ;
6180 const int s1 = facettes(i_facette, 0) ;
6181 x1 = sommets(s1,0) ;
6182 y1 = sommets(s1,1) ;
6183 const int s2 = facettes(i_facette, 1) ;
6184 x2 = sommets(s2,0) ;
6185 y2 = sommets(s2,1) ;
6186 double eps = 1e-20 ;
6192 z1 = sommets(s1,2) ;
6193 z2 = sommets(s2,2) ;
6194 const int s3 = facettes(i_facette, 2) ;
6195 x3 = sommets(s3,0) ;
6196 y3 = sommets(s3,1) ;
6197 z3 = sommets(s3,2) ;
6201 while( i< A.
size() && test_liste==0 )
6203 double plan1 = D[i] + A[i] * x1 + B[i] * y1 + C[i] * z1 ;
6204 double plan2 = D[i] + A[i] * x2 + B[i] * y2 + C[i] * z2 ;
6206 plan3 = D[i] + A[i] * x3 + B[i] * y3 + C[i] * z3 ;
6209 if( std::fabs(plan1)<eps && std::fabs(plan2)<eps && std::fabs(plan3)<eps )
6217 double& a,
double& b,
double& c,
double& d )
6220 const DoubleTabFT& sommets = maillage.
sommets() ;
6221 const IntTab& facettes = maillage.
facettes() ;
6224 double x, y, z = 0. ;
6225 double inverse_norme = 0. ;
6227 const int s1 = facettes(i_facette, 0) ;
6229 a = normale_facettes(i_facette, 0) ;
6230 b = normale_facettes(i_facette, 1) ;
6236 c = normale_facettes(i_facette, 2) ;
6239 inverse_norme = 1. / sqrt( a*a + b*b + c*c ) ;
6241 a *= inverse_norme ;
6242 b *= inverse_norme ;
6243 c *= inverse_norme ;
6244 d = - (a * x + b * y + c * z) ;
6249 double& b,
double& c,
double& d )
6253 const DoubleTabFT& sommets = maillage.
sommets() ;
6254 const IntTab& facettes = maillage.
facettes() ;
6257 double x, y, z = 0. ;
6258 double inverse_norme = 0. ;
6260 const int s1 = facettes(i_facette, 0) ;
6262 a = normale_facettes(i_facette, 0) ;
6263 b = normale_facettes(i_facette, 1) ;
6269 c = normale_facettes(i_facette, 2) ;
6272 inverse_norme = 1. / sqrt( a*a + b*b + c*c ) ;
6274 a *= inverse_norme ;
6275 b *= inverse_norme ;
6276 c *= inverse_norme ;
6277 d = - (a * x + b * y + c * z) ;
6281 const int voisin1,
const DoubleTab& indicatrice,
double& tol )
6287 const IntTab& elem_faces = domaine_vf.
elem_faces();
6292 if( voisin0 > -1 && indicatrice(voisin0) != 0. && indicatrice(voisin0) != 1. )
6295 const int face_voisine = elem_faces(voisin0,ori) + elem_faces(voisin0,ori+dim) - i_face ;
6296 L(ori) = std::fabs(domaine_vdf.
distance_face(i_face,face_voisine,ori)) ;
6298 if( voisin1 > -1 && indicatrice(voisin1) != 0. && indicatrice(voisin1) != 1. )
6301 const int face_voisine = elem_faces(voisin1,ori) + elem_faces(voisin1,ori+dim) - i_face ;
6302 double xx = std::fabs(domaine_vdf.
distance_face(i_face,face_voisine,ori)) ;
6303 L(ori) = std::max( L(ori), xx ) ;
6306 for(
int k = 0 ; k<dim ; k++ )
6310 const int face0 = elem_faces(voisin,k) ;
6311 const int face1 = elem_faces(voisin,k+dim) ;
6312 L(k) = std::fabs(domaine_vdf.
distance_face(face0,face1,k))/2. ;
6315 for(
int k = 0 ; k<dim ; k++ )
6323 const int voisin1,
const DoubleTab& indicatrice_face,
6324 const DoubleTab& indicatrice,
double& tol )
6330 const IntVect& orientation = domaine_vdf.
orientation() ;
6331 const IntTab& elem_faces = domaine_vf.
elem_faces() ;
6343 const int face_voisine = elem_faces(voisin0,ori) + elem_faces(voisin0,ori+dim) - i_face ;
6344 if( indicatrice_face(face_voisine) != 0. && indicatrice_face(face_voisine) != 1. )
6346 L(ori) = std::fabs(domaine_vdf.
distance_face(i_face,face_voisine,ori)) ;
6348 const int voisin_voisin = faces_elem(face_voisine,0) + faces_elem(face_voisine,1) - voisin0 ;
6349 if( voisin_voisin > -1 )
6351 const int face_voisine_voisine = elem_faces(voisin_voisin,ori) + elem_faces(voisin_voisin,ori+dim) - face_voisine ;
6352 L(ori) += std::fabs(domaine_vdf.
distance_face(face_voisine,face_voisine_voisine,ori)) ;
6354 for(
int k = 0 ; k<dim ; k++ )
6358 const int face0 = elem_faces(voisin0,k) ;
6359 const int face1 = elem_faces(voisin0,k+dim) ;
6360 L(k) = std::fabs(domaine_vdf.
distance_face(face0,face1,k))/2. ;
6363 for(
int k = 0 ; k<dim ; k++ )
6374 const int face_voisine = elem_faces(voisin1,ori) + elem_faces(voisin1,ori+dim) -i_face ;
6375 if( indicatrice_face(face_voisine) != 0. && indicatrice_face(face_voisine) != 1. )
6377 L(ori) = std::fabs(domaine_vdf.
distance_face(i_face,face_voisine,ori)) ;
6379 const int voisin_voisin = faces_elem(face_voisine,0) + faces_elem(face_voisine,1) - voisin1 ;
6380 if( voisin_voisin > -1 )
6382 const int face_voisine_voisine = elem_faces(voisin_voisin,ori) + elem_faces(voisin_voisin,ori+dim) - face_voisine ;
6383 L(ori) += std::fabs(domaine_vdf.
distance_face(face_voisine,face_voisine_voisine,ori)) ;
6385 for(
int k = 0 ; k<dim ; k++ )
6389 const int face0 = elem_faces(voisin1,k) ;
6390 const int face1 = elem_faces(voisin1,k+dim) ;
6391 L(k) = std::fabs(domaine_vdf.
distance_face(face0,face1,k))/2. ;
6394 for(
int k = 0 ; k<dim ; k++ )
6401 double tol_long = std::max(tol0,tol1) ;
6413 const int face_voisine = elem_faces(voisin0,ori) + elem_faces(voisin0,ori+dim) - i_face ;
6414 L(ori) = std::fabs(domaine_vdf.
distance_face(i_face,face_voisine,ori)) ;
6419 const int face_voisine = elem_faces(voisin1,ori) + elem_faces(voisin1,ori+dim) - i_face ;
6420 double xx = std::fabs(domaine_vdf.
distance_face(i_face,face_voisine,ori)) ;
6421 L(ori) = std::max( L(ori), xx ) ;
6424 DoubleTab xx_max(dim) ;
6426 for (
int i=0 ; i<2*dim ; i++)
6428 const int autre_face = elem_faces(voisin,i) ;
6429 const int face_voisine = elem_faces(voisin,ori) + elem_faces(voisin,ori+dim) - i_face ;
6431 if( autre_face != i_face && autre_face != face_voisine
6432 && indicatrice_face(autre_face) != 0. && indicatrice_face(autre_face) != 1. )
6434 const int voisin_voisin = faces_elem(autre_face,0) + faces_elem(autre_face,1) - voisin ;
6435 const int ori_face = orientation[autre_face] ;
6436 const int autre_face_voisine = elem_faces(voisin_voisin,ori_face) + elem_faces(voisin_voisin,ori_face+dim) - autre_face ;
6437 double xx = std::fabs(domaine_vdf.
distance_face(autre_face,autre_face_voisine,ori_face)) ;
6438 if( xx_max(ori_face) < xx )
6439 xx_max(ori_face) = xx ;
6443 DoubleTab TOL_TMP(dim-1) ;
6446 for(
int k=0 ; k<dim ; k++)
6450 for(
int l=0 ; l<dim ; l++)
6452 if( l != k && l != ori )
6454 const int face0 = elem_faces(voisin,k) ;
6455 const int face1 = elem_faces(voisin,k+dim) ;
6456 double yy = 0.5*std::fabs(domaine_vdf.
distance_face(face0,face1,k)) ;
6457 double a = 0.5*yy + xx_max(k) ;
6458 const int face3 = elem_faces(voisin,l) ;
6459 const int face4 = elem_faces(voisin,l+dim) ;
6460 double b = 0.5*std::fabs(domaine_vdf.
distance_face(face3,face4,l)) ;
6461 TOL_TMP(cpt) = L(ori)*L(ori) + a*a + b*b ;
6467 double tol_trans = TOL_TMP(0) ;
6469 tol_trans = std::max(tol_trans,TOL_TMP(1)) ;
6470 tol_trans = sqrt( tol_trans ) ;
6474 tol = std::max(tol_long,tol_trans) ;
6479 const DoubleTab& V, DoubleTab& coord_projete,
int& cpt )
6481 const int dim = coord_projete.
dimension(0) ;
6482 const double dist_IBC = std::fabs(d) ;
6484 double dist_proj = 0. ;
6485 int projete_modifie = 0 ;
6488 DoubleTab coord_projete_old(dim) ;
6489 coord_projete_old = coord_projete ;
6491 for (
int j = 0; j < dim ; j++)
6493 dist_proj += (x(j) - coord_projete(j))*(x(j) - coord_projete(j)) ;
6495 dist_proj = sqrt( dist_proj ) ;
6496 double dist_proj_old = dist_proj ;
6498 const double facsec = 0.5 ;
6499 const int test0 = ( dist_proj < precision ) ;
6500 const int test1 = ( std::fabs( dist_IBC - dist_proj ) > facsec*dist_IBC ) ;
6501 const int test2 = ( dist_proj > (1.+facsec)*Lref ) ;
6505 if( test0 || (test1 && test2) )
6508 for (
int j = 0; j < dim ; j++)
6510 coord_projete(j) = V(j) ;
6512 projete_modifie = 1 ;
6518 const int test = ( dist_IBC >= (1.-facsec)*Lref ) ;
6520 if( (test0 && test) || (test1 && test2) )
6523 for (
int j = 0; j < dim ; j++)
6525 coord_projete(j) = V(j) ;
6527 projete_modifie = 1 ;
6532 if( projete_modifie == 1 )
6535 for (
int j = 0; j < dim ; j++)
6537 dist_proj += (x(j) - coord_projete(j))*(x(j) - coord_projete(j)) ;
6539 dist_proj = sqrt( dist_proj ) ;
6541 Cerr <<
" ------------------------------------------------------- "<<finl;
6542 Cerr <<
" le projete "<<finl;
6543 for(
int i=0; i<dim; i++)
6545 Cerr <<
" X("<<i<<
") = "<<coord_projete_old(i)<<finl ;
6547 Cerr <<
" situe a une distance "<<dist_proj_old<<
" du barycentre de la face "<<finl ;
6548 for(
int i=0; i<dim; i++)
6550 Cerr <<
" Y("<<i<<
") = "<<x(i)<<finl ;
6552 Cerr <<
" a ete modifie en "<<finl ;
6553 for(
int i=0; i<dim; i++)
6555 Cerr <<
" X("<<i<<
") = "<<coord_projete(i)<<finl ;
6557 Cerr <<
"situe a une distance "<<dist_proj<<
" du barycentre de face "<<finl ;
6558 Cerr <<
"distance a l'IBC : "<<dist_IBC<<finl ;
6559 Cerr <<
"Longueur de reference : "<<Lref<<finl ;
6560 Cerr <<
" ------------------------------------------------------- "<<finl;
6567 const DoubleTab& x,
const DoubleTab& f,
6568 DoubleTab& x_ibc)
const
6573 assert(x_ibc.
dimension(0) == nb_colonnes) ;
6579 double norme1 = 0. ;
6580 double norme_infty = 0. ;
6582 for (
int j=0; j<nb_colonnes; j++)
6584 double norme1_int = 0. ;
6585 for (
int i=0; i<nb_lignes; i++)
6587 norme1_int += std::fabs( C(i,j) ) ;
6589 norme1 = std::max( norme1 , norme1_int ) ;
6592 for (
int i=0; i<nb_lignes; i++)
6594 double norme_infty_int = 0.;
6595 for (
int j=0; j<nb_colonnes; j++)
6597 norme_infty_int += std::fabs( C(i,j) ) ;
6599 norme_infty = std::max( norme_infty , norme_infty_int ) ;
6601 rho = 1. / ( norme1 * norme_infty ) ;
6604 DoubleTab lambda(nb_lignes) ;
6607 double distance = 1. ;
6608 double distance_old = 0. ;
6609 double deplacement_relatif = 1.;
6612 double err_uzawa = variables_internes_->seuil_uzawa ;
6613 int nb_iter_max = variables_internes_->nb_iter_uzawa ;
6626 while( deplacement_relatif > err_uzawa && nb_iter <= nb_iter_max )
6630 distance_old = distance ;
6631 for (
int i=0; i<x_ibc.
dimension(0); i++)
6633 for (
int j=0; j<lambda.
dimension(0); j++)
6635 x_ibc(i) -= 0.5 * C(j,i) * lambda(j) ;
6641 for (
int i=0; i<lambda.
dimension(0); i++)
6646 for (
int j=0; j<x_ibc.
dimension(0); j++)
6648 xx += C(i,j)*x_ibc(j) ;
6654 lambda(i) = std::max( lambda(i) + rho * xx , 0. ) ;
6659 lambda(i) = std::min( lambda(i) + rho * xx , 0. ) ;
6663 distance = ( x_ibc(0) - x[0] ) * ( x_ibc(0) - x[0] )
6664 + ( x_ibc(1) - x[1] ) * ( x_ibc(1) - x[1] ) ;
6666 distance += ( x_ibc(2) - x[2] ) * ( x_ibc(2) - x[2] ) ;
6668 deplacement_relatif = std::fabs(distance - distance_old) ;
6676 const DoubleTab& indicatrice_face,
const DoubleTab& indicatrice,
6677 const DoubleTab& dist_face,
const double t,
const double dt,
6678 DoubleTab& Tab100,DoubleTab& Tab101,DoubleTab& Tab102,DoubleTab& Tab103,
6679 IntTab& Tab12,IntTab& CptFacette,DoubleTab& v_imp,DoubleTab& Vertex,
Parser& parser_x,
Parser& parser_y,
Parser& parser_z )
6687 const IntVect& orientation = domaine_vdf.
orientation();
6689 const IntTab& elem_faces = domaine_vf.
elem_faces();
6691 const int nfaces = faces_elem.
dimension(0) ;
6692 const DoubleTab& xv = domaine_vf.
xv();
6693 const double invalid_test = -1.e30;
6695 Cerr <<
"CALCUL DU PROJETE POUR LES FACES FLUIDE ET d > invalid_test" << finl ;
6697 for (
int i_face = 0; i_face < nfaces; i_face++)
6699 double d = dist_face(i_face) ;
6700 if (indicatrice_face(i_face) == 0. && d > invalid_test)
6702 DoubleList mat_0, mat_1, mat_2, mat_3 ;
6703 ArrOfBit fa7(dimfa7);
6705 int nb_facettes_to_uzawa = 0 ;
6707 for(
int ii=0; ii<2; ii++ )
6709 const int elem = faces_elem(i_face,ii) ;
6714 for (
int face_loc=0 ; face_loc<2*dim ; face_loc++)
6716 const int autre_face = elem_faces(elem,face_loc) ;
6718 if (indicatrice_face(autre_face) != 0. && indicatrice_face(autre_face) != 1.)
6722 for(
int ii2=0 ; ii2<2; ii2++)
6724 const int elem_voisin = faces_elem(autre_face,ii2) ;
6728 if( elem_voisin>-1 && indicatrice(elem_voisin) != 0.
6729 && indicatrice(elem_voisin) != 1. )
6732 for(
int i=0 ; i< CptFacette(elem_voisin) ; i++)
6734 if( !fa7.
testsetbit( Tab12(elem_voisin,i) ) )
6736 nb_facettes_to_uzawa++ ;
6737 mat_0.
add( Tab100(elem_voisin,i) ) ;
6738 mat_1.
add( Tab101(elem_voisin,i) ) ;
6739 mat_2.
add( Tab102(elem_voisin,i) ) ;
6740 mat_3.
add( Tab103(elem_voisin,i) ) ;
6749 if( nb_facettes_to_uzawa > 0 )
6751 DoubleTab matrice_plans(nb_facettes_to_uzawa, dim) ;
6753 DoubleTab secmem(nb_facettes_to_uzawa) ;
6755 for(
int i_facette=0 ; i_facette<nb_facettes_to_uzawa ; i_facette++ )
6757 matrice_plans(i_facette,0) = mat_0[i_facette] ;
6758 matrice_plans(i_facette,1) = mat_1[i_facette] ;
6760 matrice_plans(i_facette,2) = mat_2[i_facette] ;
6761 secmem[i_facette] = -mat_3[i_facette] ;
6764 for (
int i_coord=0; i_coord<dim; i_coord++)
6766 x(i_coord) = xv(i_face,i_coord) ;
6770 DoubleTab coord_projete(dim) ;
6771 coord_projete = 0. ;
6772 uzawa(d,matrice_plans,x,secmem,coord_projete) ;
6777 for (
int i_coord=0; i_coord<dim; i_coord++)
6779 V(i_coord) = Vertex(i_face,i_coord) ;
6783 const int ori = orientation[i_face] ;
6784 const int voisin0 = faces_elem(i_face,0) ;
6785 const int voisin1 = faces_elem(i_face,1) ;
6786 const int monophasique = 1 ;
6788 verifprojete(monophasique,Lref,d,x,V,coord_projete,cpt) ;
6796 for (
int j = 0; j < dim; j++)
6798 parser_x.
setVar(j, coord_projete(j)) ;
6799 parser_y.
setVar(j, coord_projete(j)) ;
6801 parser_z.
setVar(j, coord_projete(j)) ;
6808 v_imp(i_face) = parser_x.
eval() ;
6811 v_imp(i_face) = parser_y.
eval() ;
6814 v_imp(i_face) = parser_z.
eval() ;
6821 ArrOfDouble v(dim) ;
6822 v = variables_internes_->loi_horaire_->vitesse(t+dt,coord_projete);
6832 const DoubleTab& indicatrice,
const DoubleTab& dist_face,
6833 const double t,
const double dt, DoubleTab& Tab100,
6834 DoubleTab& Tab101, DoubleTab& Tab102,DoubleTab& Tab103,
6835 IntTab& Tab12,IntTab& CptFacette, DoubleTab& v_imp,
6845 const IntVect& orientation = domaine_vdf.
orientation();
6847 const IntTab& elem_faces = domaine_vf.
elem_faces();
6849 const int nfaces = faces_elem.
dimension(0) ;
6850 const DoubleTab& xv = domaine_vf.
xv() ;
6851 const double invalid_test = -1.e30;
6853 for(
int i_face=0 ; i_face<nfaces ; i_face++ )
6855 double d = dist_face(i_face) ;
6856 const int diphasique_tmp = (indicatrice_face(i_face) != 0. &&
6857 (indicatrice_face(i_face) != 1. || (indicatrice_face(i_face) == 1. && d<=0. && d> invalid_test)));
6858 const int diphasique = ( variables_internes_->is_extra_diphasique_solide ? diphasique_tmp : (diphasique_tmp && d <= 0.) ) ;
6859 const int fluide = ( indicatrice_face(i_face) == 0. ) ;
6860 const int solide = ( variables_internes_->is_extra_solide && (d>0. || (indicatrice_face(i_face) == 1. && d<-1.e20)) ) ;
6861 const int monophasique = ( fluide || solide ) ;
6863 if (diphasique || (monophasique && d> invalid_test ))
6865 DoubleList mat_0, mat_1, mat_2, mat_3 ;
6866 ArrOfBit fa7(dimfa7);
6868 int nb_facettes_to_uzawa = 0 ;
6871 DoubleTab coord_projete(dim) ;
6872 coord_projete = 0. ;
6873 int projete_calcule = 0 ;
6879 if( indicatrice_face(i_face) == 0.5 )
6882 if( faces_elem(i_face,0) > -1 )
6883 voisin = faces_elem(i_face,0) ;
6885 voisin = faces_elem(i_face,1) ;
6887 if( indicatrice(voisin) == 0. || indicatrice(voisin) == 1. )
6889 for(
int j=0 ; j<dim ; j++)
6891 coord_projete(j)=xv(i_face,j) ;
6893 projete_calcule = 1 ;
6897 if( projete_calcule == 0 )
6901 for(
int j=0 ; j<dim ; j++)
6903 coord_projete(j)=xv(i_face,j) ;
6910 for(
int ii=0 ; ii<2; ii++)
6912 const int elem_voisin = faces_elem(i_face,ii) ;
6923 if (indicatrice(elem_voisin) != 0. && indicatrice(elem_voisin) != 1. )
6925 for(
int i=0 ; i< CptFacette(elem_voisin) ; i++)
6927 if( !fa7.
testsetbit( Tab12(elem_voisin,i) ) )
6929 nb_facettes_to_uzawa++ ;
6930 mat_0.
add( Tab100(elem_voisin,i) ) ;
6931 mat_1.
add( Tab101(elem_voisin,i) ) ;
6932 mat_2.
add( Tab102(elem_voisin,i) ) ;
6933 mat_3.
add( Tab103(elem_voisin,i) ) ;
6939 for(
int j=0 ; j<2*dim ; j++)
6941 const int elem_voisin_voisin = faces_elem(elem_faces(elem_voisin,j),0)
6942 + faces_elem(elem_faces(elem_voisin,j),1) - elem_voisin ;
6944 if( elem_voisin_voisin>-1 && indicatrice(elem_voisin_voisin) != 0.
6945 && indicatrice(elem_voisin_voisin) != 1.
6946 && elem_voisin_voisin!=faces_elem(i_face,0)
6947 && elem_voisin_voisin!=faces_elem(i_face,1) )
6949 for(
int i=0 ; i< CptFacette(elem_voisin_voisin) ; i++)
6951 if( !fa7.
testsetbit( Tab12(elem_voisin_voisin,i) ) )
6953 nb_facettes_to_uzawa++ ;
6954 mat_0.
add( Tab100(elem_voisin_voisin,i) ) ;
6955 mat_1.
add( Tab101(elem_voisin_voisin,i) ) ;
6956 mat_2.
add( Tab102(elem_voisin_voisin,i) ) ;
6957 mat_3.
add( Tab103(elem_voisin_voisin,i) ) ;
6965 if( nb_facettes_to_uzawa == 0 && diphasique )
6967 Cerr <<
"erreur concernant la face : "<< i_face <<
" I_f = "<<indicatrice_face(i_face)<<finl ;
6968 Cerr <<
"aucune facette detecte "<< finl ;
6971 else if (nb_facettes_to_uzawa >0)
6974 DoubleTab matrice_plans(nb_facettes_to_uzawa, dim) ;
6976 DoubleTab secmem(nb_facettes_to_uzawa) ;
6978 for(
int i_facette=0 ; i_facette<nb_facettes_to_uzawa ; i_facette++ )
6980 matrice_plans(i_facette,0) = mat_0[i_facette] ;
6981 matrice_plans(i_facette,1) = mat_1[i_facette] ;
6983 matrice_plans(i_facette,2) = mat_2[i_facette] ;
6984 secmem[i_facette] = -mat_3[i_facette] ;
6986 for (
int i_coord=0; i_coord<dim; i_coord++)
6988 x(i_coord) = xv(i_face,i_coord) ;
6990 coord_projete = 0. ;
6991 uzawa(d,matrice_plans,x,secmem,coord_projete) ;
6996 for (
int i_coord=0; i_coord<dim; i_coord++)
6998 V(i_coord) = Vertex(i_face,i_coord) ;
7002 const int ori = orientation[i_face] ;
7003 const int voisin0 = faces_elem(i_face,0) ;
7004 const int voisin1 = faces_elem(i_face,1) ;
7005 const int monophasique2 = 0 ;
7007 verifprojete(monophasique2,Lref,d,x,V,coord_projete,cpt) ;
7015 const int test = (diphasique || fluide || solide ) ;
7023 for (
int j = 0; j < dim; j++)
7025 parser_x.
setVar(j, coord_projete(j)) ;
7026 parser_y.
setVar(j, coord_projete(j)) ;
7028 parser_z.
setVar(j, coord_projete(j)) ;
7035 v_imp(i_face) = parser_x.
eval() ;
7038 v_imp(i_face) = parser_y.
eval() ;
7041 v_imp(i_face) = parser_z.
eval() ;
7049 v = variables_internes_->loi_horaire_->vitesse(t+dt,coord_projete);
7066static void deplacer_maillage_ft_v_impose(
Noms expression_vitesse,
7069 const double dt = temps - maillage.
temps();
7070 Parser parser_x, parser_y, parser_z;
7071 init_parser_v_impose(expression_vitesse,
7072 parser_x, parser_y, parser_z, temps);
7074 const DoubleTab& sommets = maillage.
sommets();
7076 const int dimension3 = (dim == 3);
7078 const int nb_sommets = maillage.
nb_sommets();
7079 DoubleTabFT deplacement(nb_sommets, dim);
7080 for (i = 0; i < nb_sommets; i++)
7082 const double x0 = sommets(i, 0);
7083 const double y0 = sommets(i, 1);
7084 const double z0 = dimension3 ? sommets(i, 2) : 0.;
7085 double x = x0, y = y0, z = z0;
7086 integrer_vitesse_imposee(parser_x, parser_y, parser_z,
7087 temps, dt, x, y, z);
7088 deplacement(i, 0) = x - x0;
7089 deplacement(i, 1) = y - y0;
7091 deplacement(i, 2) = z - z0;
7107 const IntTab& facettes=maillage.
facettes ();
7108 const int dim3 = (deplacement.
line_size() == 3);
7109 ArrOfInt compo_connexe_facettes(nb_facettes);
7110 int n = search_connex_components_local_FT(maillage, compo_connexe_facettes);
7111 int nb_compo_tot=compute_global_connex_components_FT(maillage, compo_connexe_facettes, n);
7113 DoubleTabFT normale;
7114 calculer_normale_sommets_interface(maillage, normale);
7119 deplacement, Vitesses, Positions);
7125 IntVect compo_sommets;
7129 const int dim = deplacement.
dimension(1);
7130 for (
int iface = 0; iface < nb_facettes; iface++)
7132 const int compo = compo_connexe_facettes[iface];
7133 for (
int j = 0; j < dim; j++)
7134 compo_sommets[facettes(iface, j)] = compo;
7142 const int nb_sommets_tot = compo_sommets.
size_totale();
7143 for (
int som = 0; som < nb_sommets_tot; som++)
7148 deplacement(som, 0) = 100.;
7149 deplacement(som, 1) = 100.;
7153 const int compo = compo_sommets[som];
7171 double nx = 0., ny = 0., nz = 0.;
7172 double vi_x = 0., vi_y = 0., vi_z = 0.;
7174 nx = normale(som, 0);
7175 ny = normale(som, 1);
7177 vi_x = deplacement(som, 0);
7178 vi_y = deplacement(som, 1);
7181 vi_z = deplacement(som, 2);
7186 (vi_x - Vitesses(compo, 0)) * nx
7187 + (vi_y - Vitesses(compo, 1)) * ny;
7189 prodscal += (vi_z - Vitesses(compo, 2)) * nz;
7190 double norme_carre = nx * nx + ny * ny + nz * nz;
7193 if (norme_carre != 0.)
7195 prodscal /= norme_carre;
7196 deplacement(som, 0) = nx * prodscal* (1-is_solid_particle) + Vitesses(compo, 0);
7197 deplacement(som, 1) = ny * prodscal* (1-is_solid_particle) + Vitesses(compo, 1);
7199 deplacement(som, 2) = nz * prodscal* (1-is_solid_particle) + Vitesses(compo, 2);
7211 DoubleTab& deplacement = variables_internes_->deplacement_sommets;
7212 const Navier_Stokes_std& eqn_hydraulique = variables_internes_->refequation_vitesse_transport.valeur();
7235#if DEBUG_CONSERV_VOLUME
7240 DoubleTab norm_deplacement(n);
7243 for (
int i=0; i<n; i++)
7246 for (
int j=0; j<deplacement.
dimension(1); j++)
7247 x += deplacement(i,j)*deplacement(i,j);
7248 norm_deplacement[i] = sqrt(x);
7252 dmin = min_array(norm_deplacement);
7253 dmax = max_array(norm_deplacement);
7260#if DEBUG_CONSERV_VOLUME
7263 double dmean_wpch = 0.;
7264 for (
int i=0; i<n; i++)
7267 for (
int j=0; j<deplacement.
dimension(1); j++)
7268 x += deplacement(i,j)*deplacement(i,j);
7269 norm_deplacement[i] = sqrt(x);
7270 dmean_wpch +=sqrt(x);
7273 double dmin_wpch = min_array(norm_deplacement);
7274 double dmax_wpch = max_array(norm_deplacement);
7275 Cerr <<
"Interfacial_velocity before/after pch [min/mean/max]. Time: "
7277 <<
" Without-pch: " << dmin <<
" " << dmean<<
" " << dmax
7278 <<
" With-pch: " << dmin_wpch <<
" " << dmean_wpch <<
" " << dmax_wpch
7279 <<
" [Warning, values invalid in //] " << finl;
7283 if (interpolation_repere_local_)
7285 DoubleTab Positions,Vitesses;
7292 fout.open(
"composantes_connexes.txt",ios::app);
7293 fout <<
"TEMPS: " << temps << std::endl;
7294 const int nb_bulles=Positions.
dimension(0);
7296 for(
int bulle=0; bulle<nb_bulles; bulle++)
7298 fout <<
"composante " << bulle <<
" position " ;
7299 for (
int i=0; i <Positions.
dimension(1); i++) fout <<
" " << Positions(bulle,i);
7300 fout <<
" vitesse " ;
7301 for (
int i=0; i <Vitesses.
dimension(1); i++) fout <<
" " << Vitesses(bulle,i);
7308 const double delta_t = temps - maillage.
temps();
7310 deplacement *= delta_t;
7313#if DEBUG_CONSERV_VOLUME
7316 if (variables_internes_->VOFlike_correction_volume > 0)
7320 DoubleTabFT position_precedente = maillage.
sommets();
7325 ArrOfDoubleFT var_volume;
7335 const DoubleTab& val_champ_vitesse = (bool)(explicit_u_NS_) ? champ_vitesse.
futur() : champ_vitesse.
valeurs();
7338#if DEBUG_CONSERV_VOLUME
7340 double sum_before_rm = 0.;
7341 double sum_before_rm_dvol = 0.;
7342 for (
int i = 0; i < nb_elem; i++)
7343 sum_before_rm +=-dI_dt[i];
7345 sum_before_rm_dvol +=sum_before_rm*delta_t;
7350 ramasse_miettes(maillage, variables_internes_->tmp_flux->valeurs(), dI_dt);
7353#if DEBUG_CONSERV_VOLUME
7355 for (
int i = 0; i < nb_elem; i++)
7358 const double dvoldt_totale =
remaillage_interface().calculer_somme_dvolume(maillage, var_volume);
7359 const double sum_dvol =sum*delta_t;
7360 Cerr <<
" time " << temps <<
" sum_dI_dt " << sum <<
" sum_dvol " << sum_dvol
7361 <<
" sum_before_rm_dI_dt " << sum_before_rm <<
" sum_before_rm_dvol " << sum_before_rm_dvol
7362 <<
" V_lagrangien= " << dvoldt_totale <<finl;
7369 var_volume *= -delta_t;
7371#if DEBUG_CONSERV_VOLUME
7387 ArrOfDoubleFT var_volume_deplacement;
7389 position_precedente,
7390 var_volume_deplacement);
7391#if DEBUG_CONSERV_VOLUME
7393 double dvol_reel_depl =
remaillage_interface().calculer_somme_dvolume(maillage, var_volume_deplacement);
7394 Cerr <<
"Transport_Interfaces_FT_Disc::calculer_vitesse_repere_local " << finl
7395 <<
" volume avt= " << volume_avt << finl
7396 <<
" apres= " << volume_apres << finl
7397 <<
" dvol_theo_depl= " << dvol_theo_depl << finl
7398 <<
" dvol_reel_depl= " << dvol_reel_depl << finl
7403 var_volume -= var_volume_deplacement;
7406 if (variables_internes_->volume_impose_phase_1 > 0.)
7408 DoubleVect values(3);
7413 if (variables_internes_->nom_domaine_volume_impose_ ==
"??")
7418 const DoubleVect& volumes = domaine_vf.
volumes();
7420 const int nb_elem_sous_domaine = sous_domaine.
nb_elem_tot();
7421 const DoubleTab& indic = indicatrice_->valeurs();
7422 const int nb_elem = domaine_vf.
nb_elem();
7423 for (
int i = 0; i < nb_elem_sous_domaine; i++)
7425 const int elem = sous_domaine[i];
7428 values(0) += indic(elem) * volumes(elem);
7429 values(1) += volumes(elem);
7430 values(2) += (1.-indic(elem)) * volumes(elem);
7437 Cerr <<
"Volume_sous_domaine " << variables_internes_->nom_domaine_volume_impose_ <<
7441 const double erreur_volume = variables_internes_->volume_impose_phase_1 - values(0);
7442 Journal() <<
"Transport_Interfaces_FT_Disc::mettre_a_jour "
7443 <<
"correction_volume_impose_phase_1= " << erreur_volume << finl;
7446 double surface_tot = 0.;
7452 for (i = 0; i < nb_facettes; i++)
7453 surface_tot += surfaces[i];
7455 surface_tot =
mp_sum(surface_tot) * 3.;
7456 if (surface_tot > 0.)
7458 const double inv_surface_tot = 1. / surface_tot;
7459 const IntTab& facettes = maillage.
facettes();
7460 const int nb_som_facette = facettes.
dimension(1);
7461 for (i = 0; i < nb_facettes; i++)
7464 const double dvolume = erreur_volume * surfaces[i] * inv_surface_tot;
7465 for (j = 0; j < nb_som_facette; j++)
7467 int isom = facettes(i,j);
7469 var_volume[isom] -= dvolume;
7481 variables_internes_->nb_lissage_correction_volume);
7483#if DEBUG_CONSERV_VOLUME
7485 double dvol_total_before =
remaillage_interface().calculer_somme_dvolume(maillage, var_volume);
7486 Cerr <<
" dvol_error_before= " << dvol_total_before <<
" Volume_before= " << volume_before <<
" time " << temps << finl;
7492 variables_internes_->nb_iterations_correction_volume);
7493#if DEBUG_CONSERV_VOLUME
7496 Cerr <<
" dvol_error_after= " << dvol_total_after <<
" Volume_after= " << volume_after <<
" time " << temps << finl;
7501#if DEBUG_CONSERV_VOLUME
7503 DoubleTabFT position_precedente = maillage.
sommets();
7510#if DEBUG_CONSERV_VOLUME
7512 ArrOfDoubleFT var_volume_deplacement;
7514 position_precedente,
7515 var_volume_deplacement);
7521 double dvol_reel_depl =
remaillage_interface().calculer_somme_dvolume(maillage, var_volume_deplacement);
7522 Cerr <<
"Transport_Interfaces_FT_Disc::calculer_vitesse_repere_local " << finl
7523 <<
" volume avt= " << volume_avt << finl
7524 <<
" apres= " << volume_apres << finl;
7526 Cerr <<
" dvol_reel_depl= " << dvol_reel_depl << finl
7527 <<
" dt= " << delta_t << finl;
7533 if (collision_model_ && nb_particles_tot_>1)
7541 const Equation_base& eqn_hydraulique = variables_internes_->refequation_vitesse_transport.valeur();
7554 DoubleTabFT d2(deplacement);
7561 la_roue_de_vitesse_a_deja_tournee);
7565 for (
int i = 0; i < n; i++)
7567 for (
int j = 0; j < dim; j++)
7569 const double depl2 = d2(i,j);
7570 deplacement(i,j) -= depl2;
7589 flags_compo_a_supprimer, maillage,
7590 variables_internes_->indicatrice_cache->valeurs());
7596 bool indic_updated =
false;
7602 indic_updated =
true;
7605 ArrOfInt liste_elems_sous_domaine;
7609 const int nb_elems_sous_domaine = sous_domaine.
nb_elem_tot();
7610 for (i = 0; i < nb_elems_sous_domaine; i++)
7612 const int i_elem = sous_domaine[i];
7613 const double indic = indicatrice[i_elem];
7614 if (indic != phase_continue)
7617 const int sz_liste_elems_sous_domaine = liste_elems_sous_domaine.
size_array();
7618 if (
mp_sum(sz_liste_elems_sous_domaine) > 0)
7624 ArrOfInt flags_compo_a_supprimer(nb_compo);
7625 for (i = 0; i < sz_liste_elems_sous_domaine; i++)
7627 const int elem = liste_elems_sous_domaine[i];
7628 const int compo = num_compo[elem];
7629 flags_compo_a_supprimer[compo] = 1;
7634 variables_internes_->indicatrice_cache->valeurs());
7636 if (max_array(flags_compo_a_supprimer))
7639 indic_updated =
true;
7644 for (
int ii = 0; ii < n; ii++)
7651 const int i_phase_continue = (phase_continue < 0.5) ? 0 : 1;
7658 return indic_updated;
7670 Cerr<<
"The method Transport_Interfaces_FT_Disc::integrer_ensemble_lagrange"<<finl;
7671 Cerr<<
"must be overloaded."<<finl;
7677 Process::Journal() <<
"Transport_Interfaces_FT_Disc::mettre_a_jour " <<
le_nom() <<
" temps= "<< temps << finl;
7688 bool update_indic = !is_indic_up_to_date;
7696 return is_indicatrice_up_to_date;
7702 switch (variables_internes_->methode_transport)
7706 deplacer_maillage_ft_v_impose(variables_internes_->expression_vitesse_imposee, maillage, temps);
7717 const double dt = temps - maillage.
temps();
7720 const double t = maillage.
temps() + dt * 0.5;
7723 Cerr <<
"Transport_Interfaces_FT_Disc::mettre_a_jour imposed motion by a time law at t:" << t << finl;
7725 const int nb_sommets = maillage.
nb_sommets();
7728 const DoubleTab& sommets = maillage.
sommets();
7731 ArrOfDouble coord(dim);
7732 ArrOfDouble coord_barycentre(dim);
7733 int nb_sommets_reels=0;
7734 for (
int i = 0; i < nb_sommets; i++)
7737 int est_un_sommet_reel=0;
7741 est_un_sommet_reel=1;
7744 for (
int j = 0; j < dim; j++)
7746 double pos = sommets(i, j);
7748 if (est_un_sommet_reel)
7749 coord_barycentre[j] += pos;
7752 coord = variables_internes_->loi_horaire_->position(t+dt,t,coord);
7755 for (
int j = 0; j < dim; j++)
7756 deplacement(i, j) = coord[j] - sommets(i,j);
7766 for (
int j = 0; j < dim; j++)
7768 coord_barycentre[j]=coord_barycentre[j]/nb_sommets_reels;
7772 variables_internes_->loi_horaire_->imprimer(
schema_temps(),coord_barycentre);
7781 Cerr <<
"Error for the method Transport_Interfaces_FT_Disc::mettre_a_jour :\n"
7782 <<
" the requested transport mehod is not developped" << finl;
7802 return indic_updated;
7807 bool indic_updated =
false;
7809 if (variables_internes_->injection_interfaces_temps_.size_array() > 0)
7811 const ArrOfDouble& tps = variables_internes_->injection_interfaces_temps_;
7813 const Noms& expr = variables_internes_->injection_interfaces_expressions_;
7814 assert(expr.size() == n);
7816 const double last_time = variables_internes_->injection_interfaces_last_time_;
7818 for (i = 0; i < n; i++)
7819 if (tps[i] > last_time)
7821 for (; i < n && tps[i] <= temps; i++)
7823 variables_internes_->injection_interfaces_last_time_=tps[i];
7829 DoubleTab sauvegarde(variables_internes_->indicatrice_cache->valeurs());
7832 variables_internes_->indicatrice_cache->valeurs(),
7834 variables_internes_->distance_interface_sommets);
7838 Cerr <<
"Injection_interface time " << temps <<
" " << expr[i];
7843 indic_updated =
true;
7844 double unused_vol_phase_0 = 0.;
7846 unused_vol_phase_0= 0.;
7848 double volume = volume_phase_1-volume_phase_1_old;
7850 volume*=pow(-1.,1-phase);
7851 Cerr <<
" volume " << volume << finl;
7855 Cerr <<
" failure: collision" << finl;
7856 variables_internes_->indicatrice_cache->valeurs() = sauvegarde;
7857 indic_updated =
false;
7862 Cerr <<
"The injection list is fulfilled !"
7863 <<
" (if it is wished, add an interface at a large time in the file)"
7869 return indic_updated;
7874 bool indic_updated =
false;
7886 const Nom expr =
Nom(
"x^2+(y-") +
Nom(Rc)
7887 +
Nom(
"*cos(") +
Nom(thetac)
7888 +
Nom(
"*pi/180.0))^2-") +
Nom(Rc)
7900 DoubleTab sauvegarde (
7901 variables_internes_->indicatrice_cache.valeur ().valeurs ());
7904 expr, 0., maillage_tmp,
7905 variables_internes_->indicatrice_cache.valeur ().valeurs (), phase,
7906 variables_internes_->distance_interface_sommets);
7908 Cerr <<
"Injection_interface time " << temps <<
" " << expr;
7913 indic_updated =
true;
7914 double unused_vol_phase_0 = 0.;
7916 sauvegarde, unused_vol_phase_0);
7917 unused_vol_phase_0 = 0.;
7919 variables_internes_->indicatrice_cache.valeur ().valeurs (),
7920 unused_vol_phase_0);
7921 double volume = volume_phase_1 - volume_phase_1_old;
7923 volume *= pow (-1., 1 - phase);
7924 Cerr <<
" volume " << volume << finl;
7929 t_injection_ = t_present_;
7933 Cerr <<
" failure: collision" << finl;
7934 variables_internes_->indicatrice_cache.valeur ().valeurs () =
7936 indic_updated =
false;
7943 for (
int ii = 0; ii < n; ii++)
7958 const Nom expr =
Nom(
"x^2+(y-") +
Nom(Rc)
7959 +
Nom(
"*cos(") +
Nom(thetac)
7960 +
Nom(
"*pi/180.0))^2-") +
Nom(Rc)
7972 DoubleTab sauvegarde (
7973 variables_internes_->indicatrice_cache.valeur ().valeurs ());
7976 expr, 0., maillage_tmp,
7977 variables_internes_->indicatrice_cache.valeur ().valeurs (), phase,
7978 variables_internes_->distance_interface_sommets);
7980 Cerr <<
"Injection_interface time " << temps <<
" " << expr;
7985 indic_updated =
true;
7986 double unused_vol_phase_0 = 0.;
7988 sauvegarde, unused_vol_phase_0);
7989 unused_vol_phase_0 = 0.;
7991 variables_internes_->indicatrice_cache.valeur ().valeurs (),
7992 unused_vol_phase_0);
7993 double volume = volume_phase_1 - volume_phase_1_old;
7995 volume *= pow (-1., 1 - phase);
7996 Cerr <<
" volume " << volume << finl;
8007 t_injection_ = t_present_;
8013 Cerr <<
" failure: collision" << finl;
8014 variables_internes_->indicatrice_cache.valeur ().valeurs () =
8016 indic_updated =
false;
8022 return indic_updated;
8026 ,
const bool update_indic)
8036 variables_internes_->distance_interface->changer_temps(temps);
8037 variables_internes_->normale_interface->changer_temps(temps);
8040 variables_internes_->indicatrice_cache->changer_temps(temps);
8041 indicatrice_->changer_temps(temps);
8048 indicatrice_faces_->changer_temps(temps);
8060 double volume_phase_0 = 0.;
8064 Cerr <<
"Volume_phase_0 " <<
Nom(volume_phase_0,
"%20.14g") <<
" time " << temps << finl;
8065 Cerr <<
"Volume_phase_1 " <<
Nom(volume_phase_1,
"%20.14g") <<
" time " << temps << finl;
8079 for (i = 0; i < nb_facettes; i++)
8085 const double s_tot =
mp_sum(s);
8088 Cerr <<
"Surface_Totale_Interface " <<
le_nom()
8089 <<
" t= " << temps <<
" surface= " << s_tot << finl;
8094 const DoubleTab& indic = indicatrice_->valeurs();
8096 const DoubleTab& xp = domaine_vf.
xp();
8098 const DoubleVect& volumes = domaine_vf.
volumes();
8100 const int nb_elem = domaine_vf.
nb_elem();
8102 DoubleTrav values(3,dim);
8113 for (i = 0; i < nb_elem; i++)
8115 const double ind = indic(i);
8116 const double v = volumes(i);
8117 const double v1 = ind * v;
8118 const double v0 = (1. - ind) * v;
8122 for (j = 0; j < dim; j++)
8124 const double x = xp(i,j);
8125 values(0,j) += v0 * x;
8126 values(1,j) += v1 * x;
8130 double mp_g0[3] = {0., 0., 0.};
8131 double mp_g1[3] = {0., 0., 0.};
8144 if (values(2,0) > 0.)
8147 for ( j = 0; j<dim; j++)
8148 mp_g0[j] = values(0,j) / values(2,0);
8150 if (values(2,1) > 0.)
8151 for ( j = 0; j<dim; j++)
8152 mp_g1[j] = values(1,j) / values(2,1);
8158 Cerr <<
"Centre_gravite_phases " <<
le_nom() <<
" t= " << temps;
8159 Cerr <<
" phase0: " << mp_g0[0] <<
" " << mp_g0[1] <<
" " << mp_g0[2] ;
8160 Cerr <<
" phase1: " << mp_g1[0] <<
" " << mp_g1[1] <<
" " << mp_g1[2] << finl;
8167 const ArrOfInt& index_elem = intersections.
index_elem();
8168 DoubleTab& surface = variables_internes_->surface_interface->valeurs();
8169 const int nb_elements = surface.
dimension(0);
8170 for (
int element = 0; element < nb_elements; element++)
8172 int index = index_elem[element];
8173 double surface_totale = 0.;
8181 surface[element] = surface_totale;
8184 variables_internes_->surface_interface->mettre_a_jour(temps);
8193 Equation_base& eqn_hydraulique = variables_internes_->refequation_vitesse_transport.valeur();
8228 DoubleTabFT deplacement(vitesse);
8229 deplacement *= coeff;
8247 variables_internes_->indicatrice_cache->changer_temps(temps);
8248 indicatrice_->changer_temps(temps);
8249 variables_internes_->distance_interface->changer_temps(temps);
8250 variables_internes_->normale_interface->changer_temps(temps);
8264 ref_milieu_ = un_milieu;
8276 Cerr <<
"You forgot to associate the diphasic fluid to the problem named " <<
probleme().
le_nom() << finl;
8279 return ref_milieu_.valeur();
8286 Cerr <<
"You forgot to associate the diphasic fluid to the problem named " <<
probleme().
le_nom() << finl;
8289 return ref_milieu_.valeur();
8313 return indicatrice_;
8318 return indicatrice_;
8323 return variables_internes_->marching_cubes_;
8328 return variables_internes_->marching_cubes_;
8333 return variables_internes_->maillage_interface;
8338 return variables_internes_->maillage_interface;
8342 return variables_internes_->remaillage_interface_;
8347 return variables_internes_->remaillage_interface_;
8351 return variables_internes_->topologie_interface_;
8356 return variables_internes_->topologie_interface_;
8361 return variables_internes_->doubletab_pos;
8366 return variables_internes_->intvect_elements;
8371 return variables_internes_->deplacement_sommets;
8376 return variables_internes_->proprietes_particules_;
8381 return variables_internes_->proprietes_particules_;
8386 return variables_internes_->maillage_inject_;
8391 return variables_internes_->maillage_inject_;
8396 return variables_internes_->proprietes_inject_;
8401 return variables_internes_->proprietes_inject_;
8415 std::vector<YAML_data> dvi = variables_internes_->data_a_sauvegarder();
8416 data.insert(data.end(), dvi.begin(), dvi.end());
8427 int special, afaire;
8430 Nom mon_ident(
"variables_internes_transport");
8431 mon_ident +=
Nom(temps,
"%e");
8432 variables_internes_->set_time(temps);
8437 os << mon_ident << finl;
8438 os << variables_internes_->que_suis_je() << finl;
8443 os << mon_ident << finl;
8444 os << variables_internes_->que_suis_je() << finl;
8446 bytes += variables_internes_->sauvegarder(os);
8449 Cerr <<
"Backup of collision model" << finl;
8450 if (collision_model_) bytes += collision_model_.valeur().sauvegarder(os);
8453 Cerr <<
"Backup of particles position and velocity" << finl;
8468 Cerr <<
"Backup of collision model not supported with pdi format" << finl;
8487 Cerr <<
"Transport_Interfaces_FT_Disc::reprendre" << finl;
8493 is >>
id >> type_name;
8494 if ( (!
id.debute_par(
"variables_internes_transport"))
8495 || type_name != variables_internes_->
que_suis_je())
8497 Cerr <<
"Error for the method Transport_Interfaces_FT_Disc::reprendre" << finl;
8498 Cerr << variables_internes_->que_suis_je() <<
" was expected."<< finl;
8502 variables_internes_->maillage_interface.associer_equation_transport(*
this);
8504 variables_internes_->set_restart_fname(restart_fname);
8505 variables_internes_->set_time(
schema_temps().temps_courant());
8510 if (collision_model_) collision_model_.valeur().reprendre(is);
8511 const int nb_particles_tot=collision_model_.valeur().get_nb_particles_tot();
8527 Cerr <<
"Backup of collision model not supported with pdi format" << finl;
8544 return probleme_base_.valeur();
8549 return variables_internes_->parcours_interface_;
8554 return variables_internes_->connectivite_frontieres_;
8559 return variables_internes_->algorithmes_transport_.valeur();
8564 my_map& map_post, DoubleTab *ftab)
const
8566 const Motcle som =
"sommets";
8567 const Motcle elem =
"elements";
8568 const Motcle bi =
"elements et sommets";
8569 const DoubleTab dummytab;
8765 map_post.emplace(
"heat_transfer", map_element_post_FT(elem, &Transport_Interfaces_FT_Disc::fill_ftab_scalar,
8784 const Motcle som =
"sommets";
8785 const Motcle elem =
"elements";
8786 const Motcle bi =
"elements et sommets";
8788 Transport_Interfaces_FT_Disc::my_map map_post;
8789 fill_map_post_FT(map_post,ftab);
8793 map_element_post_FT::func_type the_function=&Transport_Interfaces_FT_Disc::fill_ftab_scalar;
8794 DoubleTab the_tab_values;
8795 bool elem_found=
false;
8798 for (
const auto& [key, value] : map_post)
8803 the_loc=value.location_;
8804 the_function=value.function_;
8805 the_tab_values=value.values_;
8813 Cerr<<
"The real fields to be post-processed are :"<<finl;
8815 for (
const auto& [key, value] : map_post)
8817 Cerr <<
" Fields("<<i<<
") : " << key <<
", Location : ";
8819 Cerr << value.location_ << finl;
8824 else if (!elem_found)
8826 else if (! (the_loc==bi
8840 (this->*the_function)(ftab,the_tab_values);
8857 const Motcle som =
"sommets";
8858 const Motcle elem =
"elements";
8859 const Motcle bi =
"elements et sommets";
8860 const int nb_champs = 5;
8865 fields[2] =
"numero";
8866 fields[3] =
"pe_local";
8867 fields[4] =
"compo_connexe";
8869 Motcles localisations(nb_champs);
8871 localisations[0] = bi;
8872 localisations[1] = bi;
8873 localisations[2] = bi;
8874 localisations[3] = bi;
8875 localisations[4] = elem;
8878 int rank=fields.
search(champ), i;
8882 Cerr<<
"The integer fields to be post-processed are:"<<finl;
8883 for (i=1 ; i<nb_champs ; i++)
8885 Cerr <<
" Fields("<<i<<
") : "<< fields[i] <<
" # Localisations : " << localisations[i] << finl;
8894 else if (! (localisations[rank]==bi
8923 for (i2=0 ; i2<n ; i2++)
8925 (*itab)(i2) = pe_som[i2];
8932 for (i2=0 ; i2<n ; i2++)
8934 (*itab)(i2) = pe_fac[i2];
8941 for (i2=0 ; i2<n ; i2++)
8957 int n2 = search_connex_components_local_FT(maillage, compo);
8958 compute_global_connex_components_FT(maillage, compo, n2);
8960 for (
int ii = 0; ii < nbf; ii++)
8961 (*itab)[ii] = compo[ii];
8966 Cerr <<
"Transport_Interfaces_FT_Disc::get_champ_post_FT : unexpected case" << finl;
8988 return variables_internes_->n_iterations_distance;
9002 if (tag == variables_internes_->distance_sommets_cache_tag)
9004 return variables_internes_->distance_interface_sommets;
9006 variables_internes_->distance_sommets_cache_tag = tag;
9007 assert(tag == variables_internes_->distance_normale_cache_tag);
9010 DoubleTab& dist_som = variables_internes_->distance_interface_sommets;
9045 const DoubleTab& dist_elem,
9046 const DoubleTab& normale_elem,
9047 DoubleTab& dist_som)
const
9049 static const double distance_sommets_invalides = -1.e30;
9053 const DoubleTab& xp = domaine_vf.
xp();
9057 ArrOfInt ncontrib(nb_sommets);
9062 double centre[3] = {0., 0., 0.};
9063 double normale[3] = {0., 0., 0.};
9067 const int nb_som_elem = elem_som.
line_size();
9069 for (elem = 0; elem < nb_elem_tot; elem++)
9071 const double d1 = dist_elem(elem);
9073 if (d1 < distance_sommets_invalides)
9076 for (i = 0; i < dim; i++)
9078 centre[i] = xp(elem, i);
9079 normale[i] = normale_elem(elem, i);
9082 for (i = 0; i < nb_som_elem; i++)
9084 const int som = elem_som(elem, i);
9087 if (som < nb_sommets)
9091 for (j = 0; j < dim; j++)
9093 double position_sommet = coord_som(som, j);
9094 d2 += (position_sommet - centre[j]) * normale[j];
9096 dist_som(som) += d1 + d2;
9102 const double valeur_invalide = distance_sommets_invalides * 1.1;
9103#ifdef __INTEL_COMPILER
9106 for (i = 0; i < nb_sommets; i++)
9108 const int n = ncontrib[i];
9112 dist_som(i) = valeur_invalide;
9115 Debog::verifier(
"Transport_Interfaces_FT_Disc::calculer_distance_interface_sommets",dist_som);
9135 DoubleTab& distance_elements,
9136 DoubleTab& normale_elements,
9137 const int n_iter)
const
9139 statistics().create_custom_counter(
"Calculer_distance_interface",3,
"FrontTracking");
9140 statistics().begin_count(
"Calculer_distance_interface",statistics().get_last_opened_counter_level()+1);
9142 static const double distance_sommets_invalides = -1.e30;
9147 const DoubleTab& centre_element = domaine_vf.
xp();
9151 distance_elements = distance_sommets_invalides * 1.1;
9152 normale_elements = 0.;
9165 const ArrOfInt& index_elem = intersections.
index_elem();
9168 const IntTab& facettes = maillage.
facettes();
9169 const DoubleTab& sommets = maillage.
sommets();
9171 for (
int elem = 0; elem < nb_elem; elem++)
9173 int index = index_elem[elem];
9175 double normale[3] = {0., 0., 0.};
9177 double centre[3] = {0., 0., 0.};
9179 double surface_tot = 0.;
9187#ifdef AVEC_BUG_SURFACES
9188 const double surface = data.surface_intersection_;
9192 surface_tot += surface;
9193 for (
int i = 0; i < dim; i++)
9195 normale[i] += surface * normale_facettes(num_facette, i);
9198 for (
int j = 0; j < dim; j++)
9200 const int som = facettes(num_facette, j);
9201 const double coord = sommets(som, i);
9203 g_i += coord * coeff;
9205 centre[i] += surface * g_i;
9209 if (surface_tot > 0.)
9214 const double inverse_surface_tot = 1. / surface_tot;
9217 for (j = 0; j < dim; j++)
9219 norme += normale[j] * normale[j];
9220 normale_elements(elem, j) = normale[j];
9221 centre[j] *= inverse_surface_tot;
9225 double i_norme = 1./sqrt(norme);
9226 double distance = 0.;
9227 for (j = 0; j < dim; j++)
9229 double n_j = normale[j] * i_norme;
9230 distance += (centre_element(elem, j) - centre[j]) * n_j;
9232 distance_elements(elem) = distance;
9238 Debog::verifier(
"Transport_Interfaces_FT_Disc::calculer_distance_interface",normale_elements);
9239 Debog::verifier(
"Transport_Interfaces_FT_Disc::calculer_distance_interface",distance_elements);
9242 DoubleTab terme_src(normale_elements);
9243 DoubleTab tmp(normale_elements);
9245 const IntTab& face_voisins = domaine_vf.
face_voisins();
9246 const IntTab& elem_faces = domaine_vf.
elem_faces();
9247 const int nb_elem_voisins = elem_faces.
line_size();
9251 for (iteration = 0; iteration < n_iter; iteration++)
9260 const double un_sur_ncontrib = 1. / (1. + nb_elem_voisins);
9262 for (elem = 0; elem < nb_elem; elem++)
9265 double n[3] = {0., 0., 0.};
9266 for (i = 0; i < dim; i++)
9267 n[i] = normale_elements(elem, i);
9268 for (k = 0; k < nb_elem_voisins; k++)
9271 const int face = elem_faces(elem, k);
9272 const int e_voisin = face_voisins(face, 0) + face_voisins(face, 1) - elem;
9274 for (i = 0; i < dim; i++)
9275 n[i] += normale_elements(e_voisin, i);
9277 for (i = 0; i < dim; i++)
9278 tmp(elem, i) = terme_src(elem, i) + n[i] * un_sur_ncontrib;
9280 normale_elements = tmp;
9285 ArrOfIntFT liste_elements;
9288 for (elem = 0; elem < nb_elem; elem++)
9290 double nx = normale_elements(elem, 0);
9291 double ny = normale_elements(elem, 1);
9292 double nz = (dim==3) ? normale_elements(elem, 2) : 0.;
9293 double norme2 = nx*nx + ny*ny + nz*nz;
9296 double i_norme = 1. / sqrt(norme2);
9297 normale_elements(elem, 0) = nx * i_norme;
9298 normale_elements(elem, 1) = ny * i_norme;
9300 normale_elements(elem, 2) = nz * i_norme;
9307 terme_src = distance_elements;
9308 tmp = distance_elements;
9309 for (iteration = 0; iteration < n_iter; iteration++)
9312 const int liste_elem_size = liste_elements.
size_array();
9313 for (i_elem = 0; i_elem < liste_elem_size; i_elem++)
9315 elem = liste_elements[i_elem];
9316 if (terme_src(elem) > distance_sommets_invalides)
9319 tmp(elem) = distance_elements(elem);
9324 double ncontrib = 0.;
9325 double somme_distances = 0.;
9327 for (k = 0; k < nb_elem_voisins; k++)
9330 const int face = elem_faces(elem, k);
9331 const int e_voisin = face_voisins(face, 0) + face_voisins(face, 1) - elem;
9334 const double distance_voisin = distance_elements(e_voisin);
9335 if (distance_voisin > distance_sommets_invalides)
9338 double nx = normale_elements(elem, 0) + normale_elements(e_voisin, 0);
9339 double ny = normale_elements(elem, 1) + normale_elements(e_voisin, 1);
9340 double nz = (dim==3)
9341 ? normale_elements(elem, 2) + normale_elements(e_voisin, 2) : 0.;
9342 double norm2 = nx*nx + ny*ny + nz*nz;
9345 double i_norm = 1./sqrt(norm2);
9351 double dx = centre_element(elem, 0) - centre_element(e_voisin, 0);
9352 double dy = centre_element(elem, 1) - centre_element(e_voisin, 1);
9353 double dz = (dim==3)
9354 ? centre_element(elem, 2) - centre_element(e_voisin, 2) : 0.;
9355 double d = nx * dx + ny * dy + nz * dz + distance_voisin;
9356 somme_distances += d;
9364 double d = somme_distances / ncontrib;
9369 distance_elements = tmp;
9372 statistics().end_count(
"Calculer_distance_interface");
9376 const ArrOfInt& compo_connexes_facettes,
9377 const int nb_compo_tot,
9378 const DoubleTab& vitesse_sommets,
9379 DoubleTab& vitesses,
9380 DoubleTab& positions)
const
9382 assert(nb_compo_tot == vitesses.
dimension(0));
9383 assert(nb_compo_tot == positions.
dimension(0));
9388 const IntTab& facettes = maillage.
facettes();
9389 const DoubleTab& sommets = maillage.
sommets();
9393 ArrOfDouble surfaces_compo(nb_compo_tot);
9400 for (
int i = 0; i < nb_facettes_tot; i++)
9404 const int compo = compo_connexes_facettes[i];
9405 const double surface = surface_facettes[i];
9406 surfaces_compo[compo] += surface;
9408 for (
int j = 0; j < dim; j++)
9411 const int s = facettes(i, j);
9412 for (
int k = 0; k < dim; k++)
9414 positions(compo, k) += surface * sommets(s, k);
9420 positions *= (1. / dim);
9424 tab_divide_any_shape(positions, s);
9430 for (
int fa7 = 0; fa7 < nb_facettes_tot; fa7++)
9435 const double si=surface_facettes[fa7];
9438 for (
int d=0; d<dim; d++)
9439 som[d]=facettes(fa7,d);
9446 const double n0=normale_facettes(fa7,0);
9447 const double n1=normale_facettes(fa7,1);
9448 ds_dt=-n1*vitesse_sommets(som[0],0)+n0*vitesse_sommets(som[0],1)+n1*vitesse_sommets(som[1],0)-n0*vitesse_sommets(som[1],1);
9452 const double n0=normale_facettes(fa7,0)*0.5;
9453 const double n1=normale_facettes(fa7,1)*0.5;
9454 const double n2=normale_facettes(fa7,2)*0.5;
9455 double s2s1[3],d_surface[3];
9457 for (
int i = 0; i < 3; i++)
9464 const int s0 = som[i];
9465 const int s1 = som[ (i+1)%3 ];
9466 const int s2 = som[ (i+2)%3 ];
9468 s2s1[0] = sommets(s1,0) - sommets(s2,0);
9469 s2s1[1] = sommets(s1,1) - sommets(s2,1);
9470 s2s1[2] = sommets(s1,2) - sommets(s2,2);
9472 d_surface[0] = s2s1[1] * n2 - s2s1[2] * n1;
9473 d_surface[1] = s2s1[2] * n0 - s2s1[0] * n2;
9474 d_surface[2] = s2s1[0] * n1 - s2s1[1] * n0;
9476 ds_dt+=d_surface[0]*vitesse_sommets(s0,0) + d_surface[1]*vitesse_sommets(s0,1) +d_surface[2]*vitesse_sommets(s0,2);
9480 const int compo = compo_connexes_facettes[fa7];
9481 for (
int d=0; d<dim; d++)
9485 for (
int d1=0; d1<dim; d1++)
9487 V+=vitesse_sommets(som[d1],d);
9488 p+=sommets(som[d1],d);
9492 vitesses(compo, d) += si * V + (p-positions(compo, d)) * ds_dt;
9500 tab_divide_any_shape(vitesses, s);
9524 DoubleVect& valeurs)
9530 const Domaine& domaine = domaine_vf.
domaine();
9531 const IntTab& face_voisins = domaine_vf.
face_voisins();
9532 const IntTab& elem_faces = domaine_vf.
elem_faces();
9535 const int nb_elem_tot = domaine.nb_elem_tot();
9539 for (
int i_face = 0; i_face < nb_faces_tot; i_face++)
9542 const int elem0 = face_voisins(i_face, 0);
9543 const int elem1 = face_voisins(i_face, 1);
9544 if (elem0 < 0 || elem1 < 0)
9548 double prod_scal = 0.;
9552 for (
int i = 0; i < dim; i++)
9554 double n0 = normale_interface(elem0, i);
9555 double n1 = normale_interface(elem1, i);
9556 prod_scal += (n0 + n1) * 0.5 * domaine_vf.
face_normales(i_face, i);
9558 if (indic(elem0) == 1. || indic(elem1) == 1.)
9559 prod_scal = - prod_scal;
9562 if (indic(elem0) != 0. && indic(elem0) != 1.)
9564 if (indic(elem1) != 0. && indic(elem1) != 1.)
9580 DoubleVect tmp_flux = flux;
9581 ArrOfInt flag(nb_faces_elem);
9582 for (
int i_elem = 0; i_elem < nb_elem_tot; i_elem++)
9586 for (
int j = 0; j < nb_faces_elem; j++)
9588 int i_face = elem_faces(i_elem, j);
9589 const int elem0 = face_voisins(i_face, 0);
9590 const int elem1 = face_voisins(i_face, 1);
9591 if (elem0 < 0 || elem1 < 0)
9594 const double signe = (elem0 == i_elem) ? 1. : -1.;
9595 const double f = flux(i_face);
9605 const double facteur = valeurs(i_elem) / somme;
9606 for (
int j = 0; j < nb_faces_elem; j++)
9608 int i_face = elem_faces(i_elem, j);
9611 tmp_flux[i_face] = tmp_flux[i_face] * facteur;
9618 for (
int i_elem = 0; i_elem < nb_elem_tot; i_elem++)
9620 for (
int j = 0; j < nb_faces_elem; j++)
9622 int i_face = elem_faces(i_elem, j);
9623 flux[i_face]=tmp_flux[i_face];
9628 for (
int i_elem = 0; i_elem < nb_elem_tot; i_elem++)
9630 double x = valeurs(i_elem);
9631 for (
int j = 0; j < nb_faces_elem; j++)
9633 int i_face = elem_faces(i_elem, j);
9634 double f = flux(i_face);
9635 double signe = (i_elem == face_voisins(i_face, 1)) ? 1. : -1.;
9638 valeurs(i_elem) = x;
9648 const DoubleVect& valeurs_euler,
9649 ArrOfDouble& valeurs_lagrange)
9651 const int nb_sommets = maillage.
nb_sommets();
9654 valeurs_lagrange = 0.;
9657 const IntTab& facettes = maillage.
facettes();
9659 const ArrOfInt& index_elem = intersections.
index_elem();
9661 const int nb_elements = valeurs_euler.
size();
9662 assert(nb_elements == index_elem.
size_array());
9666 for (element = 0; element < nb_elements; element++)
9671 const int index_premiere_intersection = index_elem[element];
9672 if (index_premiere_intersection < 0)
9680 double surface_totale = 0.;
9682 int index = index_premiere_intersection;
9687 const double surface_facette = surface_facettes[facette];
9692 const double valeur_euler = valeurs_euler[element];
9696 index = index_premiere_intersection;
9697 if (surface_totale > 0.)
9699 const double inv_surface_tot = 1. / surface_totale;
9704 const double surface_facette = surface_facettes[facette];
9709 for (i = 0; i < dim; i++)
9711 const int sommet = facettes(facette, i);
9717 const double valeur = coeff * fraction_surface * valeur_euler;
9720 valeurs_lagrange[sommet] += valeur;
9732void Transport_Interfaces_FT_Disc::compute_nb_particles_tot()
9736 ArrOfInt id_number_fa7(nb_facettes);
9737 int n = search_connex_components_local_FT(mesh, id_number_fa7);
9738 nb_particles_tot_= compute_global_connex_components_FT(mesh, id_number_fa7, n);
9749 ArrOfInt id_number_fa7(nb_facettes);
9750 int n = search_connex_components_local_FT(mesh, id_number_fa7);
9751 int nb_particles_tot=compute_global_connex_components_FT(mesh, id_number_fa7, n);
9756 Cerr <<
"WARNING, renumbering particles !" << finl;
9763 const IntTab& facettes = mesh.
facettes();
9764 const DoubleTab& sommets = mesh.
sommets();
9766 DoubleVect particles_surfaces(nb_particles_tot);
9771 for (
int i = 0; i < nb_facettes_tot; i++)
9775 const int id_number_lagrangian_vertice = id_number_fa7[i];
9776 const double surface = surface_facettes[i];
9777 particles_surfaces[id_number_lagrangian_vertice] += surface;
9782 const int s = facettes(i, j);
9785 surface * sommets(s, k);
9798 Equation_base& eqn_hydraulique = variables_internes_->refequation_vitesse_transport.valeur();
9822 Equation_base& eqn_hydraulique = variables_internes_->refequation_vitesse_transport.valeur();
9825 DoubleTab correct_particles_position(nb_particles_tot,
dimension);
9826 DoubleTab correct_particles_velocity(nb_particles_tot,
dimension);
9827 IntVect particles_correct_id_number(nb_particles_tot);
9834 for (
int wrong_id_number=0; wrong_id_number<nb_particles_tot; wrong_id_number++)
9837 int correct_id_number;
9838 if (particle_gravity_center_elem==-1) correct_id_number =-1;
9839 else correct_id_number = particles_eulerian_id_number[particle_gravity_center_elem];
9840 particles_correct_id_number(wrong_id_number)=correct_id_number;
9843 int isduplicateValue = collision_model_.valeur().check_for_duplicates(particles_correct_id_number);
9844 if (isduplicateValue == 1)
Process::exit(
"Transport_Interfaces_FT_Disc::swap_particles_position_velocity "
9845 "ERROR: duplicate value of the particles lagrangian ID number");
9848 for (
int wrong_id_number = 0; wrong_id_number < nb_particles_tot; wrong_id_number++)
9850 int good_id_number = particles_correct_id_number[wrong_id_number];
9864 const DoubleVect& mesh_volumes = domain_vf.
volumes();
9865 const IntTab& elem_faces = domain_vf.
elem_faces();
9866 const DoubleTab& phase_indicator_function = indicatrice_->valeurs();
9867 Equation_base& eqn_hydraulique = variables_internes_->refequation_vitesse_transport.valeur();
9879 for (
int elem=0; elem<domain_vf.
nb_elem(); elem++)
9882 if (phase_indicator_function(elem)==id_solid_phase)
9884 const int particle= particles_eulerian_id_number(elem);
9889 tab_velocity(elem_faces(elem,dim))+
9890 tab_velocity(elem_faces(elem,dim+
dimension))
9891 )*mesh_volumes(elem);
9894 0.5*(tab_velocity(elem_faces(elem,dim))+
9895 tab_velocity(elem_faces(elem,dim+
dimension))
9896 )*mesh_volumes(elem),2);
9904 DoubleVect s_vparticles;
9909 for (
int particle=0; particle<nb_particles_tot_; particle++)
9921void Transport_Interfaces_FT_Disc::fill_ftab_scalar(DoubleTab *ftab,
9922 const ArrOfDouble& values)
const
9926 for (
int fa7=0 ; fa7<nb_fa7 ; fa7++)
9927 (*ftab)(fa7,0) = (
float) values(fa7);
9930void Transport_Interfaces_FT_Disc::fill_ftab_scalar(DoubleTab *ftab,
9931 const DoubleVect& values)
const
9935 for (
int fa7=0 ; fa7<nb_fa7 ; fa7++)
9936 (*ftab)(fa7,0) = (
float) values(fa7);
9939void Transport_Interfaces_FT_Disc::fill_ftab_scalar(DoubleTab *ftab,
9940 const DoubleTab& values)
const
9944 for (
int fa7=0 ; fa7<nb_fa7 ; fa7++)
9945 (*ftab)(fa7,0) = (
float) values(fa7);
9948void Transport_Interfaces_FT_Disc::fill_ftab_vector(DoubleTab *ftab,
const DoubleTab& values)
const
9952 const int nb_compo = values.
dimension(1);
9953 ftab->
resize(nb_fa7, nb_compo);
9954 for (
int fa7=0 ; fa7<nb_fa7 ; fa7++)
9955 for (
int k=0 ; k<nb_compo ; k++)
9956 (*ftab)(fa7,k) = (
float) values(fa7,k);
9967 if (variables_internes_->refequation_vitesse_transport)
9972 const Equation_base& eqn_hydraulique = variables_internes_->refequation_vitesse_transport.valeur();
9991 fill_ftab_vector(ftab,vit);
9996 Cerr <<
"Error : velocity nodes post-processing : to be developped." << finl;
10005 if (variables_internes_->refequation_vitesse_transport)
10008 DoubleTab Positions,Vitesses;
10011 const Equation_base& eqn_hydraulique = variables_internes_->refequation_vitesse_transport.valeur();
10024 fill_ftab_vector(ftab,vit);
10029 Cerr <<
"Error : vitesse_repere_local nodes post-processing : to be developped." << finl;
10056 fill_ftab_vector(ftab,values);
: classe Algorithmes_Transport_FT_Disc
int testsetbit(int_t i) const
Renvoie la valeur du bit e, puis met le bit e a 1.
class Champ_Face_VDF Cette classe sert a representer un champ vectoriel dont on ne calcule
classe Champ_Fonc_Face_VDF
DoubleTab & passe(int i=1) override
Renvoie les valeurs du champs a l'instant t-i.
double changer_temps(const double temps) override
Fixe le temps du champ.
DoubleTab & valeurs() override
Renvoie le tableau des valeurs du champ au temps courant.
virtual DoubleTab & futur(int i=1)
virtual DoubleTab & valeurs()=0
classe Champ_base Cette classe est la base de la hierarchie des champs.
virtual const Domaine_dis_base & domaine_dis_base() const
virtual DoubleTab & valeur_aux_elems_passe(const DoubleTab &positions, const IntVect &les_polys, DoubleTab &tab_valeurs) const
double temps() const
Renvoie le temps du champ.
virtual DoubleTab & valeur_aux_elems(const DoubleTab &positions, const IntVect &les_polys, DoubleTab &valeurs) const
provoque une erreur ! doit etre surchargee par les classes derivees
classe Conds_lim Cette classe represente un vecteur de conditions aux limites.
bool is_reinject_activated() const
bool ready_injection() const
const double & get_Rc_inject() const
const double & get_thetaC() const
virtual void suppression_interfaces(const IntVect &num_compo, const ArrOfInt &flags_compo_a_supprimer, int nouvelle_phase)
Methode appelee par Transport_Interfaces_xxx::test_suppression_interfaces_sous_domaine() lorqu'une in...
static void verifier(const char *const msg, double)
void collecter_espace_virtuel(ArrOfDouble &tab, MD_Vector_tools::Operations_echange op) const
void echange_espace_virtuel(ArrOfDouble &tab) const
classe Discretisation_base Cette classe represente un schema de discretisation en espace,...
void discretiser_champ(const Motcle &directive, const Domaine_dis_base &z, const Nom &nom, const Nom &unite, int nb_comp, int nb_pas_dt, double temps, OWN_PTR(Champ_Inc_base)&champ, const Nom &sous_type=NOM_VIDE) const
void construire_elem_virt_pe_num()
const Sous_Domaine_t & ss_domaine(int i) const
SmallArrOfTID_t & chercher_elements(const DoubleTab &pos, SmallArrOfTID_t &elem, int reel=0) const
Recherche des elements contenant les points dont les coordonnees sont specifiees.
virtual void creer_tableau_elements(Array_base &, RESIZE_OPTIONS opt=RESIZE_OPTIONS::COPY_INIT) const
creation d'un tableau parallele de valeurs aux elements.
DoubleTab_t & les_sommets()
int nb_faces_elem(int=0) const
Renvoie le nombre de face de type i des elements geometriques constituants le domaine.
virtual void creer_tableau_sommets(Array_base &, RESIZE_OPTIONS opt=RESIZE_OPTIONS::COPY_INIT) const
Cree un tableau ayant une "ligne" par sommet du maillage.
int orientation(int) const override
inline DoubleVect& Domaine_VDF::porosite_face() {
double distance_face(int, int, int k) const
virtual const DoubleVect & face_surfaces() const
DoubleVect & volumes_entrelaces()
int nb_faces_tot() const
renvoie le nombre total de faces.
void creer_tableau_faces(Array_base &, RESIZE_OPTIONS opt=RESIZE_OPTIONS::COPY_INIT) const
virtual double face_normales(int face, int comp) const
double xv(int num_face, int k) const
const Joint & joint(int i) const
double volumes(int i) const
void construire_face_virt_pe_num()
Remplissage du tableau face_virt_pe_num_ (voir commentaire dans Domaine_VF.
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
double xp(int num_elem, int k) const
int face_voisins(int num_face, int i) const
renvoie l'element voisin de numface dans la direction i.
int moments_a_imprimer() const
const Nom & le_nom() const override
Donne le nom de l'Objet_U Methode a surcharger : renvoie "neant" dans cette implementation.
classe Domaine_dis_base Cette classe est la base de la hierarchie des domaines discretisees.
virtual IntTab & face_voisins()
const Domaine & domaine() const
Lecture dans un fichier d'objets ecrits au format binaire.
Fichier en lecture Cette classe est a la classe C++ ifstream ce que la classe Entree est a la.
static Nom & get_Output()
Renvoie le mode d'ecriture utilise (pour pouvoir le modifier).
static int is_ecriture_special(int &special, int &a_faire)
indique si le format special a ete demande en lecture active par sauvegarde xyz .
static int is_lecture_special()
indique si le format special a ete demande en lecture active par reprise xyz .
Fichier en lecture Cette classe est a la classe C++ ifstream ce que la classe Entree est a la.
Class defining operators and methods for all reading operation in an input flow (file,...
classe Equation_base Le role d'une equation est le calcul d'un ou plusieurs champs....
virtual void set_param(Param &titi) const override
const Nom & le_nom() const override
Renvoie le nom de l'equation.
int reprendre(Entree &) override
On reprend l'inconnue a partir d'un flot d'entree.
virtual std::vector< YAML_data > data_a_sauvegarder() const
for PDI IO: retrieve name, type and dimensions of the data to save/restore. This has to be overrode f...
virtual void associer_pb_base(const Probleme_base &)
S'associe au Probleme passe en parametre.
const Discretisation_base & discretisation() const
Renvoie la discretisation associee a l'equation.
Solveur_Masse_base & solv_masse()
Renvoie le solveur de masse associe a l'equation.
virtual const Champ_Inc_base & inconnue() const =0
int calculate_time_derivative() const
virtual int nombre_d_operateurs() const =0
int sauvegarder(Sortie &) const override
On sauvegarde l'inconnue, puis les sources sur un flot de sortie.
virtual const Champ_Inc_base & derivee_en_temps() const
Probleme_base & probleme()
Renvoie le probleme associe a l'equation.
int lire_motcle_non_standard(const Motcle &, Entree &) override
Lecture des parametres de type non simple d'un objet_U a partir d'un flot d'entree.
int limpr() const
Demande au schema en temps si il faut effectuer une impression.
Schema_Temps_base & schema_temps()
Renvoie le schema en temps associe a l'equation.
Champs_compris champs_compris_
virtual const Operateur & operateur(int) const =0
Domaine_dis_base & domaine_dis()
Renvoie le domaine discretise associe a l'equation.
const int & get_id_fluid_phase() const
: Classe de postraitement des champs euleriens au format lata
int ecrire_entete(const double temps_courant, const int reprise, const int est_le_premier_post) override
Ouvre le fichier maitre en mode ERASE et ecrit l'entete du fichier lata (sur le processeur maitre seu...
int ecrire_domaine(const Domaine &domaine, const int est_le_premier_post) override
voir Format_Post_base::ecrire_domaine On accepte l'ecriture d'un domaine dans un pas de temps,...
int ecrire_temps(const double temps) override
commence l'ecriture d'un nouveau pas de temps En l'occurence pour le format LATA:
int ecrire_champ(const Domaine &domaine, const Noms &unite_, const Noms &noms_compo, int ncomp, double temps_, const Nom &id_du_champ, const Nom &id_du_domaine, const Nom &localisation, const Nom &nature, const DoubleTab &data) override
voir Format_Post_base::ecrire_champ
int initialize(const Nom &file_basename, const int format, const Nom &option_para) override
static Objet_U & objet(const Nom &)
Voir Interprete_bloc::objet_global() BM: la classe Interprete n'est pas le meilleur endroit pour cett...
double fraction_surface_intersection_
int index_facette_suivante_
: class Intersections_Elem_Facettes
const ArrOfInt & index_elem() const
Renvoie un tableau de taille domaine.
const Intersections_Elem_Facettes_Data & data_intersection(int index) const
Renvoie les donnees de l'intersection stockee a l'indice "index" dans le tableau "data" ( 0 <= index ...
const Joint_Items_t & joint_item(JOINT_ITEM type) const
Renvoie les informations de joint pour le type demande.
const IntTab_t & renum_items_communs() const
Voir renum_items_communs_.
Cette classe implemente les operateurs et les methodes virtuelles de la classe EFichier de la facon s...
int get_sequential_items_flags(ArrOfBit &flags, int line_size=1) const
: class Maillage_FT_Disc Cette classe decrit un maillage:
virtual Entree & lire_param_maillage(Entree &is)
Cette fonction permet de lire les parametres pour le maillage des interfaces.
void mesh_tag_increase() const
int get_mesh_tag() const
return mesh_state_tag_
void set_is_solid_particle(const bool is_solid_particle)
void preparer_tableau_avant_transport(ArrOfDouble &tableau, const Desc_Structure_FT &descripteur) const
Prepare un tableau de donnees aux sommets ou aux facettes pour conserver les valeurs apres transport.
int nb_sommets() const
renvoie le nombre de sommets (reels et virtuels) (egal a sommets().
const DoubleTab & sommets() const
renvoie le tableau des sommets (reels et virtuels) dimension(0) = nombre de sommets,
void associer_equation_transport(const Equation_base &equation) override
on remplit refequation_transport_, schema_comm_domaine_ desc_sommets_.comm_group_ et desc_facettes_....
void update_tableau_apres_transport(ArrOfDouble &tableau, int nb_elements, const Desc_Structure_FT &descripteur) const
Voir preparer_tableau_avant_transport.
int nb_facettes() const
renvoie le nombre de facettes (reelles et virtuelles) (egal a facettes().
int facette_virtuelle(int i) const
Renvoie 0 si la facette m'appartient, 1 sinon.
double temps() const
return temps_physique_ (temps_physique_ ne sert a rien en interne.
const Desc_Structure_FT & desc_sommets() const
renvoie le descripteur des sommets (espace_distant/virtuel)
void completer_maillage()
Complete les structures de donnees du maillage.
virtual void ajouter_maillage(const Maillage_FT_Disc &maillage_tmp, int skip_facettes=0)
double changer_temps(double t)
return temps_physique_ = t
virtual const DoubleTab & get_update_normale_facettes() const
Calcule la grandeur demandee, stocke le resultat dans un tableau interne a la classe et renvoie le re...
void compute_gravity_center_fa7()
void facette_PE_owner(ArrOfInt &facette_pe) const
pour postraitement, remplit le tableau en parametre avec le numero du PE proprietaire de chaque facet...
void parcourir_maillage()
Remplit la structure intersections_elem_facettes_.
const Schema_Comm & get_schema_comm_FT() const
int sommet_virtuel(int i) const
const ArrOfInt & sommet_PE_owner() const
pour postraitement, renvoie le numero du PE proprietaire des sommets
void calcul_indicatrice(DoubleVect &indicatrice, const DoubleVect &indicatrice_precedente) const
Calcul de la fonction indicatrice (on suppose que "indicatrice" a la structure d'un tableau de valeur...
virtual void transporter(const DoubleTab &deplacement)
Deplace les sommets de l'interface d'un vecteur "deplacement" fourni, Change eventuellement les somme...
const Intersections_Elem_Facettes & intersections_elem_facettes() const
void nettoyer_elements_virtuels()
Retire toutes les facettes virtuelles et tous les sommets qui ne sont pas utilises.
virtual const ArrOfDouble & get_update_surface_facettes() const
Calcule la grandeur demandee, stocke le resultat dans un tableau interne a la classe et renvoie le re...
virtual const ArrOfDouble & get_update_courbure_sommets() const
Calcule la grandeur demandee, stocke le resultat dans un tableau interne a la classe et renvoie le re...
void creer_tableau_sommets(Array_base &, RESIZE_OPTIONS opt=RESIZE_OPTIONS::COPY_INIT) const
creation d'un tableau aux sommets du maillage Meme principe que Domaine::creer_tableau_sommets()
const IntTab & facettes() const
renvoie le tableau des facettes (reelles et virtuelles) dimension(0) = nombre de facettes,
const ArrOfInt & sommet_elem() const
pour postraitement, renvoie sommet_elem_
int construire_iso(const DoubleVect &valeurs_sommets, double isovaleur, Maillage_FT_Disc &maillage, DoubleVect &indicatrice_approchee, const Maillage_FT_Disc::AjoutPhase phase, int ignorer_collision=0) const
Construction d'un maillage en segments ou en triangles comme l'isovaleur d'une fonction discretisee a...
void associer_domaine_vf(const Domaine_VF &domaine_vf)
Classe Matrice_Morse Represente une matrice M (creuse), non necessairement carree.
classe Milieu_base Cette classe est la base de la hierarchie des milieux (physiques)
Une chaine de caractere (Nom) en majuscules.
Un tableau d'objets de la classe Motcle.
int search(const Motcle &t) const
virtual const Champ_base * get_delta_vitesse_interface() const
Si le champ de vitesse est discontinu (calcul avec changement de phase), renvoie un pointeur vers le ...
void compute_particles_eulerian_id_number(const OWN_PTR(Collision_Model_FT_base)&collision_model_ptr)
const DoubleTab & get_interfacial_area() const
bool get_new_mass_source() const
bool get_is_solid_particle() const
virtual void calculer_dI_dt(DoubleVect &dI_dt, const DoubleTab &tab_vitesse)
virtual const Fluide_Diphasique & fluide_diphasique() const
const IntTab & get_particles_eulerian_id_number() const
void swap_particles_eulerian_id_number(const ArrOfInt &gravity_center_elem)
classe Navier_Stokes_std Cette classe porte les termes de l'equation de la dynamique
const Champ_Inc_base & inconnue() const override
Renvoie la vitesse (champ inconnue de l'equation) (version const).
class Nom Une chaine de caractere pour nommer les objets de TRUST
virtual int finit_par(const char *const n) const
virtual int find(const char *const n) const
const std::string & getString() const
Un tableau de chaine de caracteres (VECT(Nom)).
classe Objet_U Cette classe est la classe de base des Objets de TRUST
const Nom & que_suis_je() const
renvoie la chaine identifiant la classe.
virtual Entree & readOn(Entree &)
Lecture d'un Objet_U sur un flot d'entree Methode a surcharger.
static double precision_geom
static const Nom & nom_du_cas()
Renvoie une reference constante vers le nom du cas.
virtual int reprendre(Entree &)
Reprise d'un Objet_U sur un flot d'entree Methode a surcharger.
virtual Sortie & printOn(Sortie &) const
Ecriture de l'objet sur un flot de sortie Methode a surcharger.
classe Operateur Classe generique de la hierarchie des operateurs.
double calculer_pas_de_temps() const
Calcule le prochain pas de temps.
Helper class to factorize the readOn method of Objet_U classes.
void ajouter_flag(const char *keyword, const bool *value)
Register a boolean flag whose mere presence switches it to true.
void ajouter(const char *keyword, const int *value, Param::Nature nat=Param::OPTIONAL)
Register an integer parameter.
void ajouter_non_std(const char *keyword, const Objet_U *value, Param::Nature nat=Param::OPTIONAL)
Register a keyword handled by Objet_U::lire_motcle_non_standard.
void associer_connectivite_frontieres(const Connectivite_frontieres &connect)
void associer_domaine_dis(const Domaine_dis_base &domaine_dis)
Remplissage des variables persistantes de la classe (refdomaine_vf_, nb_faces_elem_,...
classe Paroi_FT_disc Condition aux limites d'angle de contact pour l'equation
Representation des donnees de la classe Parser.
void addVar(const char *)
virtual void setNbVar(int nvar)
void setVar(const char *sv, double val)
void setString(const std::string &s)
virtual void parseString()
const DoubleTab & get_dUdz_P2_Stokes_th()
const DoubleTab & get_dVdz_P2_Stokes_th()
const DoubleTab & get_dWdy_P2_Stokes_th()
const DoubleTab & get_dWdx_P2_Stokes_th()
const DoubleTab & get_sigma_yz_fa7_Stokes_th() const
const DoubleTab & get_sigma_xy_fa7_Stokes_th() const
const DoubleTab & get_dUdx_P1_Stokes_th()
const DoubleTab & get_U_P2_Stokes_th() const
const DoubleTab & get_dWdx_P1_Stokes_th()
const DoubleTab & get_sigma_xx_fa7_Stokes_th() const
const DoubleTab & get_sigma_yy_fa7_Stokes_th() const
const DoubleTab & get_sigma_xz_fa7_Stokes_th() const
const DoubleTab & get_dWdz_P1_Stokes_th()
const DoubleTab & get_dWdz_P2_Stokes_th()
const DoubleTab & get_sigma_zz_fa7_Stokes_th() const
const DoubleTab & get_pressure_force_Stokes_th_fa7() const
const DoubleTab & get_dUdx_P2_Stokes_th()
const DoubleTab & get_friction_force_Stokes_th_fa7() const
const DoubleTab & get_dWdy_P1_Stokes_th()
const DoubleTab & get_dVdz_P1_Stokes_th()
const DoubleTab & get_dUdz_P1_Stokes_th()
const DoubleTab & get_U_P1_Stokes_th() const
const DoubleTab & get_dVdy_P2()
const DoubleTab & get_sigma_yz_fa7() const
const DoubleTab & get_sigma_xy_fa7() const
const DoubleTab & get_dUdx_P2()
const DoubleTab & get_dWdz_P1()
const DoubleTab & get_sigma_yy_fa7() const
const DoubleTab & get_dWdy_P1()
const DoubleTab & get_dWdz_P2()
const DoubleTab & get_dVdz_P2()
const DoubleTab & get_sigma_yx_fa7() const
const DoubleTab & get_dWdx_P2()
const DoubleTab & get_U_P2() const
const DoubleTab & get_sigma_xx_fa7() const
const DoubleTab & get_dWdx_P1()
const DoubleTab & get_dUdy_P2()
const DoubleTab & get_sigma_zx_fa7() const
const DoubleTab & get_pressure_force_fa7() const
const DoubleTab & get_dVdx_P2()
const DoubleTab & get_friction_force_fa7() const
const DoubleTab & get_U_P1() const
const DoubleTab & get_sigma_xz_fa7() const
const DoubleTab & get_dVdx_P1()
const DoubleTab & get_sigma_zy_fa7() const
const DoubleTab & get_sigma_zz_fa7() const
const DoubleTab & get_dUdy_P1()
const DoubleTab & get_dVdz_P1()
const DoubleTab & get_dWdy_P2()
const DoubleTab & get_dUdz_P2()
const DoubleTab & get_heat_transfer_fa7() const
const DoubleTab & get_dUdz_P1()
const DoubleTab & get_dUdx_P1()
const DoubleTab & get_dVdy_P1()
static const char *const demande_description
const Triple_Line_Model_FT_Disc & tcl() const
const Nom & le_nom() const override
Donne le nom de l'Objet_U Methode a surcharger : renvoie "neant" dans cette implementation.
classe Probleme_base C'est un Probleme_U qui n'est pas un couplage.
const Nom & checkpoint_filename() const
bool is_dilatable() const
bool & reprise_effectuee()
const Schema_Temps_base & schema_temps() const
Renvoie le schema en temps associe au probleme.
virtual int nombre_d_equations() const =0
virtual const Equation_base & equation(int) const =0
const Nom & restart_filename() const
static void mp_max_for_each_item(TRUSTArray< _TYPE_ > &x, int n=-1)
static double mp_min(double)
static int check_int_overflow(trustIdType)
static void mp_sum_for_each(T &arg1, T &arg2)
C++14 compatible mp_sum_for_each: combine multiple mp_sum calls into one collective operation Usage: ...
static double mp_max(double)
static void mp_sum_for_each_item(TRUSTArray< _TYPE_ > &x, int n=-1)
static Sortie & Journal(int message_level=0)
Renvoie un objet statique de type Sortie qui sert de journal d'evenements.
static void abort()
Routine de sortie de Trio-U sur une erreur abort().
static double mp_sum(double)
Calcule la somme de x sur tous les processeurs du groupe courant.
static void barrier()
Synchronise tous les processeurs du groupe courant (attend que tous les processeurs soient arrives a ...
static int me()
renvoie mon rang dans le groupe de communication courant.
static void exit(int exit_code=-1)
Routine de sortie de TRUST dans une region Kokkos.
static int je_suis_maitre()
renvoie 1 si on est sur le processeur maitre du groupe courant (c'est a dire me() == 0),...
Classe qui porte les proprietes de particules.
void nettoyer(const ArrOfInt &som_utilises)
: class Remaillage_FT Cette classe implemente les procedures de remaillage des interfaces pour le Fro...
int traite_decollement(Maillage_FT_Disc &maillage, const DoubleTab &deplacement) const
Cette fonction permet de gerer le decollement de l'interface de la paroi Si un sommet de bord n'a pas...
int traite_adherence(Maillage_FT_Disc &maillage) const
Cette fonction permet de gerer l'adherence de l'interface a la paroi Si une facette possede 3 sommets...
void associer_domaine(const Domaine_dis_base &domaine_dis)
Cette fonction stocke le domaine_dis dans refdomaine_dis_.
double calculer_variation_volume(const Maillage_FT_Disc &maillage, const DoubleTab &position_initiale, ArrOfDouble &varVolume) const
Cette fonction calcule le volume de phase 0 engendre par chaque sommet lors du deplacement de l'inter...
void set_is_solid_particle(const bool is_solid_particle)
void lisser_dvolume(const Maillage_FT_Disc &maillage, ArrOfDouble &var_volume, const int nb_iterations) const
Regularise le champ scalaire "var_volume" defini aux sommets du "maillage".
void corriger_volume_(Maillage_FT_Disc &maillage, ArrOfDouble &var_volume, const int nb_iter_corrections_vol)
Cette classe est a la classe C++ ofstream ce que la classe Sortie est a la classe C++ ostream Elle re...
static void lire_xyz(Maillage_FT_Disc &mesh, const Domaine_VF *domaine_vf, Entree *fichier, const Domaine *domaine_src)
Remplissage du maillage mesh a partir d'un fichier de sauvegarde xyz (fichier) ou d'un domaine source...
static void init_sequential_domain(Domaine_32_64< _SIZE_ > &dom)
Create parallel descriptors for the vertex and element arrays of the domain (necessary because Scatte...
static void uninit_sequential_domain(Domaine_32_64< _SIZE_ > &dom)
methode utilisee par les interpretes qui modifient le domaine (sequentiel), detruit les descripteurs ...
int nb_impr() const
Renvoie le nombre d'impressions effectuees.
double temps_courant() const
Renvoie le temps courant.
int precision_impr() const
void imprimer_temps_courant(SFichier &) const
double pas_de_temps() const
Renvoie le pas de temps (delta_t) courant.
virtual int nb_valeurs_temporelles() const =0
void update_critere_statio(const DoubleTab &tab_critere, Equation_base &equation)
//Actualisation de stationnaire_atteint_ et residu_ (critere residu_<seuil_statio_)
classe Solveur_Masse_base Represente la matrice de masse d'une equation.
virtual DoubleTab & appliquer(DoubleTab &) const
renvoie appliquer_impl(x/coeffient_temporelle) si on a un coefficient temporel sinon renvoie applique...
virtual int ouvrir(const char *name, IOS_OPEN_MODE mode=ios::out)
void precision(int pre) override
void setf(IOS_FORMAT code) override
Classe de base des flux de sortie.
int_t nb_elem_tot() const
void append_array(_TYPE_ valeur)
_SIZE_ size_array() const
void resize_array(_SIZE_ new_size, RESIZE_OPTIONS opt=RESIZE_OPTIONS::COPY_INIT)
TRUSTList & add(_TYPE_)
insertion en queue
TRUSTList & add_if_not(_TYPE_)
Ajout d'un element a la liste ssi il n'existe pas deja.
int contient(_TYPE_) const
Verifie si un element appartient ou non a la liste.
void resize(_SIZE_ n, RESIZE_OPTIONS opt=RESIZE_OPTIONS::COPY_INIT)
_SIZE_ dimension_tot(int) const override
_SIZE_ dimension(int d) const
_SIZE_ size_totale() const
void ref_array(TRUSTArray< _TYPE_, _SIZE_ > &, _SIZE_ start=0, _SIZE_ sz=-1) override
void resize(_SIZE_, RESIZE_OPTIONS opt=RESIZE_OPTIONS::COPY_INIT)
virtual const MD_Vector & get_md_vector() const
virtual void echange_espace_virtuel(IsExchangeBlocking exchange_type=IsExchangeBlocking::DefaultBlocking, const std::string kernel_name="noname")
classe TRUST_2_PDI Encapsulation of PDI methods (library used for IO operations). See the website pdi...
void read(const std::string &name, void *data)
static void set_PDI_restart(int r)
void get_type(const Nom &name, Nom &type)
Generic method to read the type of a TRUST object in the HDF5 file.
static void set_PDI_checkpoint(int c)
static int is_PDI_checkpoint()
static int is_PDI_restart()
void TRUST_start_sharing(const std::string &name, const void *data)
: class Topologie_Maillage_FT Cette classe implemente les procedures de remaillage des interfaces pou...
virtual int calculer_composantes_connexes_pour_suppression(const Domaine_VF &domaine_vf, const DoubleTab &indicatrice, IntVect &num_compo) const
Computes eulerian connex components of the phase indicator function "indicatrice" according to the ge...
virtual void remailler_interface(const double temps, Maillage_FT_Disc &maillage, Champ_base &indicatrice, Remaillage_FT &algo_remaillage_local)
Remaillage de l'interface: - amelioration petites et grandes facettes,.
virtual double suppression_interfaces(const IntVect &num_compo, const ArrOfInt &flags_compo_a_supprimer, Maillage_FT_Disc &maillage, DoubleTab &indicatrice)
Removes all interfaces contained in eulerian elements marked by the "flags_" array,...
int get_phase_continue() const
Nom maillage_interface_xyz_filename(int restart) const
std::vector< YAML_data > data_a_sauvegarder() const
for PDI IO: retrieve name and type and dimensions of the indicatrice tag
Remaillage_FT remaillage_interface_
int indicatrice_cache_tag
OWN_PTR(Champ_Inc_base) indicatrice_cache
int reprendre(Entree &is) override
Reprise d'un Objet_U sur un flot d'entree Methode a surcharger.
Maillage_FT_Disc maillage_interface
int sauvegarder(Sortie &os) const override
Sauvegarde d'un Objet_U sur un flot de sortie Methode a surcharger.
void set_param(Param &titi) const override
DoubleTab & derivee_en_temps_inco(DoubleTab &derivee) override
Calcul de la derivee en temps de l'inconnue : zero.
const Operateur & operateur(int i) const override
void fill_ftab_pressure_force(DoubleTab *ftab, const DoubleTab &dummytab) const
const bool & get_is_solid_particle() const
void fill_ftab_friction_force(DoubleTab *ftab, const DoubleTab &dummytab) const
virtual void get_expression_vitesse_imposee(DoubleTab &vitesse_imp)
Milieu_base & milieu() override
bool compute_particles_rms_
const DoubleTab & get_mean_particles_volumic_squared_velocity() const
virtual void PPP_face_interface_voisin(const DoubleTab &indicatrice, const DoubleTab &indicatrice_face, DoubleTab &Vertex, DoubleTab &PPP)
int nombre_d_operateurs() const override
void update_indicatrice_normale_distance()
Updates normals and distances to interface, then updates indicatrice.
virtual const Champ_base & get_indicatrice_faces()
void fill_ftab_normal_unit(DoubleTab *ftab, const DoubleTab &dummytab) const
virtual double suppression_interfaces(const IntVect &num_compo, const ArrOfInt &flags_compo_a_supprimer)
virtual void preparer_pas_de_temps()
virtual int get_champ_post_FT(const Motcle &champ, Postraitement_base::Localisation loc, DoubleTab *dtab=0) const
Cherche le champ discret aux interfaces dont le nom est "champ", et verifie qu'il peut etre postraite...
void associer_equation_ns(const Navier_Stokes_FT_Disc &ns)
const DoubleTab & get_mean_particles_volumic_velocity() const
int preparer_calcul() override
Tout ce qui ne depend pas des autres problemes eventuels.
void fill_ftab_vertices_curvature(DoubleTab *ftab, const DoubleTab &dummytab) const
int sauvegarder(Sortie &) const override
On sauvegarde l'inconnue, puis les sources sur un flot de sortie.
virtual void remailler_interface()
Remaillage de l'interface : - amelioration petites et grandes facettes,.
virtual const DoubleTab & get_update_distance_interface_sommets() const
Renvoi de la distance signee entre l'interface et les sommets du maillage eulerien.
const Topologie_Maillage_FT & topologie_interface() const
int impr(Sortie &os) const override
Imprime les operateurs de l'equation sur un flot de sortie, de facon inconditionnelle.
void calcul_source(const DoubleTab &inco_val, const DoubleTab &vpoint, const DoubleTab &rho_faces, DoubleTab &source_val, const DoubleTab &vit_imposee, const DoubleTab &indicatrice_faces, const int is_QC, const double dt, const int is_explicite, const double eta)
void associer_milieu_base(const Milieu_base &milieu) override
DoubleTab particles_purely_solid_mesh_volume_
void update_indicatrice() override
Recalcul du champ variables_internes_->indicatrice_cache a partir de la position des interfaces.
bool injecter_interfaces_par_ajout_phase(double temps)
void impr_effort_fluide_interface(DoubleTab &source_val, DoubleTab &pressure_part, DoubleTab &friction_part, DoubleTab &diff_part)
void calculer_scalaire_interpole(const Champ_base &ch_scal, const Maillage_FT_Disc &, DoubleTab &ch_scal_noeuds, int nv_calc) const
DoubleTab & deplacement_som()
const DoubleTab & get_particles_purely_solid_mesh_volume() const
virtual void transporter_sans_changement_topologie(DoubleTab &vitesse, const double coeff, const double temps)
void calcul_indicatrice_faces(const DoubleTab &indicatrice, const IntTab &face_voisins)
void fill_ftab_Stokes(DoubleTab *ftab, const DoubleTab &values) const
const Proprietes_part_vol & proprietes_particules() const
virtual const Connectivite_frontieres & connectivite_frontieres() const
void init_save_file() override
void mettre_a_jour(double temps) override
La valeur de l'inconnue sur le pas de temps a ete calculee.
void nettoyer_proprietes_particules(const ArrOfInt &som_utilises)
virtual void calcul_nb_traverse(const DoubleTab &xe, const double dx, const int dim, const int ori, Maillage_FT_Disc &maillage, int elem, int &traverse)
int verif_Cl() const override
Methode appelee par Equation_base::readOn On verifie que toutes les cl sont de type Paroi_FT_disc.
void init_particles_position_velocity()
const DoubleTab & get_rms_particles_volumic_velocity() const
virtual void calculer_distance_interface(const Maillage_FT_Disc &maillage, DoubleTab &distance_elements, DoubleTab &normale_elements, const int n_iter) const
Calcul d'un champ scalaire aux elements contenant une distance signee entre le centre de l'element et...
virtual void calcul_tolerance_projete_monophasique(const int i_face, const int ori, const int voisin0, const int voisin1, const DoubleTab &indicatrice_face, const DoubleTab &indicatrice, double &tol)
DoubleTab particles_position_collision_
Remaillage_FT & remaillage_interface()
void modifie_source(DoubleTab &so_modif, const DoubleTab &so_val, const DoubleTab &rho_faces, const int n, const int m, const int is_QC, const DoubleVect &vol_entrelaces, const Solveur_Masse_base &solv_masse)
void compute_particles_rms()
virtual const Parcours_interface & parcours_interface() const
Post_Processing_Hydrodynamic_Forces_Stokes post_process_hydro_forces_Stokes_
int preparer_calcul_anticipated()
virtual const Maillage_FT_Disc & maillage_interface_pour_post() const
Renvoie le maillage stocke specialement pour le postraitement (si on veut postraiter un etat intermed...
const Champ_Inc_base & inconnue() const override
const Champ_base & get_indicatrice() override
getter champ variables_internes_->indicatrice_cache a partir de la position des interfaces.
virtual void calculer_distance_interface_faces(const DoubleTab &dist_elem, const DoubleTab &normale_elem, DoubleTab &dist_faces) const
void check_indicatrice_is_up_to_date() override
Checks if the indicator is up to date.
std::vector< YAML_data > data_a_sauvegarder() const override
for PDI IO: retrieve name and type and dimensions of the indicatrice tag
Transport_Interfaces_FT_Disc()
constructeur par defaut
bool injecter_interfaces_pour_TCL(double temps)
void update_critere_statio()
void discretiser() override
Discretisation des champs: - indicatrice_ : champ scalaire discretise aux elements.
void swap_particles_lagrangian_position_velocity()
WARNING, particles_position_collision_ and particles_velocity_collision_ are not used to transport pa...
const Maillage_FT_Disc & maillage_inject() const
virtual void calculer_vitesse_transport_interpolee(const Champ_base &champ_vitesse, const Maillage_FT_Disc &m, DoubleTab &vitesse_noeuds, int nv_calc) const
void fill_ftab_pressure(DoubleTab *ftab, const DoubleTab &dummytab) const
const Probleme_base & get_probleme_base() const
double calculer_pas_de_temps() const override
Calcul du prochain pas de temps.
void integrer_ensemble_lagrange(const double temps) override
virtual const Champ_base & update_indicatrice_faces()
const Proprietes_part_vol & proprietes_inject() const
ArrOfInt gravity_center_elem_
virtual void PPP_face_interface(Maillage_FT_Disc &maillage, const DoubleTab &indicatrice, const DoubleTab &indicatrice_face, DoubleTab &Vertex)
DoubleTab mean_particles_volumic_velocity_
virtual const Champ_base & get_normale_interface() const
virtual void RenumFa7(DoubleTab &Barycentre, DoubleTab &Tab110, DoubleTab &Tab111, DoubleTab &Tab112, IntTab &Tab12, IntTab &CptFacette, const int nb_facettes, const int nb_facettes_dim)
void fill_ftab_Stokes_pressure_th(DoubleTab *ftab, const DoubleTab &values) const
Post_Processing_Hydrodynamic_Forces post_process_hydro_forces_
virtual void PPP_face_voisin(const DoubleTab &indicatrice, const DoubleTab &indicatrice_face, DoubleTab &PPP)
virtual double calculer_integrale_indicatrice(const DoubleVect &indicatrice, double &v_ph0) const
const OWN_PTR(Collision_Model_FT_base) &get_ptr_collision_model() const
void associer_pb_base(const Probleme_base &probleme) override
S'associe au Probleme passe en parametre.
virtual const Algorithmes_Transport_FT_Disc & algorithmes_transport() const
virtual void plan_facette_existant(Maillage_FT_Disc &maillage, DoubleList A, DoubleList B, DoubleList C, DoubleList D, const int i_facette, int &test_liste)
virtual const Champ_base & get_update_distance_interface_faces() const
virtual void lire_maillage_ft_cao(Entree &is)
const Maillage_FT_Disc & maillage_interface() const
virtual void uzawa(const double d, const DoubleTab &matrice, const DoubleTab &x, const DoubleTab &secmem, DoubleTab &solution) const
virtual const Marching_Cubes & marching_cubes() const
virtual void calcul_eq_plan_facette(Maillage_FT_Disc &maillage, const int i_facette, double &a, double &b, double &c, double &d)
void completer_maillage_et_changer_temps(double temps)
DoubleTab particles_velocity_collision_
bool mettre_a_jour_deplacement(double temps)
virtual void interpoler_vitesse_face(const DoubleTab &distance_interface, const int phase, const int stencil_width, DoubleTab &champ, DoubleTab &gradient, const double t, const double dt)
void fill_ftab_velocity(DoubleTab *ftab, const DoubleTab &dummytab) const
void ramasse_miettes(const Maillage_FT_Disc &maillage, DoubleVect &flux, DoubleVect &valeurs)
int lire_motcle_non_standard(const Motcle &, Entree &) override
Lecture des parametres de type non simple d'un objet_U a partir d'un flot d'entree.
virtual void deplacer_maillage_ft_v_fluide(const double temps)
virtual void calcul_tolerance_projete_diphasique(const int i_face, const int ori, const int voisin0, const int voisin1, const DoubleTab &indicatrice, double &tol)
virtual void calcul_OutElemFa7(Maillage_FT_Disc &maillage, const DoubleTab &indicatrice, const int nb_elem, int &nb_fa7_accepted, IntList &OutElem, IntList &OutFa7)
void verifprojete(const int monophasique, const double Lref, double d, const DoubleTab &x, const DoubleTab &V, DoubleTab &coord_projete, int &cpt)
virtual void calcul_maxfa7(Maillage_FT_Disc &maillage, const DoubleTab &indicatrice, const int nb_elem, int &max_fa7, const int exec_planfa7existan)
void fill_ftab_local_reference_frame_velocity(DoubleTab *ftab, const DoubleTab &dummytab) const
static void transfert_conservatif_eulerien_vers_lagrangien_sommets(const Maillage_FT_Disc &maillage, const DoubleVect &valeurs_euler, ArrOfDouble &valeurs_lagrange)
void deplacer_maillage(double temps)
void close_save_file() override
virtual void StockageFa7(Maillage_FT_Disc &maillage, IntTab &CptFacette, DoubleTab &Tab100, DoubleTab &Tab101, DoubleTab &Tab102, DoubleTab &Tab103, DoubleTab &Tab110, DoubleTab &Tab111, DoubleTab &Tab112, IntTab &Tab12, DoubleTab &Barycentre, const DoubleTab &indicatrice, IntList &OutElem, ArrOfBit &fa7, const int exec_planfa7existant)
void mettre_a_jour_hors_deplacement(double temps, const bool update_statio=true, const bool update_indic=true)
IntVect & vecteur_elements()
Entree & lire_cond_init(Entree &is) override
Lecture des conditions initiales.
int reprendre(Entree &) override
On reprend l'inconnue a partir d'un flot d'entree.
void fill_ftab_Stokes_pressure_interp(DoubleTab *ftab, const DoubleTab &values) const
virtual void projete_point_face_interface(int &nb_proj_modif, const int dim_fa7, const DoubleTab &indicatrice_face, const DoubleTab &indicatrice, const DoubleTab &dist_face, const double t, const double dt, DoubleTab &Tab100, DoubleTab &Tab101, DoubleTab &Tab102, DoubleTab &Tab103, IntTab &Tab12, IntTab &CptFacette, DoubleTab &v_imp, DoubleTab &Vertex, Parser &parser_x, Parser &parser_y, Parser &parser_z)
virtual void BaryFa7(Maillage_FT_Disc &maillage, const int i_facette, DoubleTab &Barycentre)
OBS_PTR(Probleme_base) probleme_base_
virtual void calculer_vmoy_composantes_connexes(const Maillage_FT_Disc &maillage, const ArrOfInt &compo_connexes_facettes, const int nb_compo_tot, const DoubleTab &vitesse_sommets, DoubleTab &vitesses, DoubleTab &positions) const
virtual void calcul_vitesse(DoubleTab &vitesse_imp, const DoubleTab &champ_vitesse, const DoubleTab &vpoint, const double temps, const double dt)
virtual void calculer_vitesse_repere_local(const Maillage_FT_Disc &maillage, DoubleTab &deplacement, DoubleTab &Positions, DoubleTab &Vitesses) const
virtual void calculer_distance_interface_sommets(const DoubleTab &dist_elem, const DoubleTab &normale_elem, DoubleTab &dist_som) const
Calcule dist_som, la distance entre l'interface et les sommets du maillage eulerien a partir de dist_...
void update_normale_distance_interface() const
Calcule la normale et la distance a l'interface, evaluees sur une epaisseur egale a n_iterations_dist...
virtual void projete_point_face_fluide(int &nb_proj_modif, const int dim_fa7, const DoubleTab &indicatrice_face, const DoubleTab &indicatrice, const DoubleTab &dist_face, const double t, const double dt, DoubleTab &Tab100, DoubleTab &Tab101, DoubleTab &Tab102, DoubleTab &Tab103, IntTab &Tab12, IntTab &CptFacette, DoubleTab &v_imp, DoubleTab &Vertex, Parser &parser_x, Parser &parser_y, Parser &parser_z)
void calcul_effort_fluide_interface(const DoubleTab &vpoint, const DoubleTab &rho_faces, DoubleTab &source_val, const int is_explicite, const double eta)
Nom suppression_interfaces_sous_domaine_
virtual int calculer_composantes_connexes_pour_suppression(IntVect &num_compo)
void set_is_solid_particle(const bool is_solid_particle)
void ajouter_contribution_saut_vitesse(DoubleTab &deplacement, const bool la_roue_de_vitesse_a_deja_tournee) const
virtual const Champ_base & get_distance_interface() const
void modifier_vpoint_pour_imposer_vit(const DoubleTab &inco_val, DoubleTab &vpoint0, DoubleTab &vpoint, const DoubleTab &rho_faces, DoubleTab &terme_source, const double temps, const double dt, const int is_explicite, const double eta) override
DoubleTab & tableaux_positions()
virtual const int & get_n_iterations_distance() const
DoubleTab mean_particles_volumic_squared_velocity_
void assembler(Matrice_Morse &mat_morse, const DoubleTab &present, DoubleTab &secmem) override
const int & get_vimp_regul() const
bool test_suppression_interfaces_sous_domaine()
bool injecter_supprimer_interfaces(double temps)
DoubleTab rms_particles_volumic_velocity_
classe Transport_Interfaces_base Cette classe constitue la classe de base des equations de transport ...
bool reinjection_tcl() const
const double & t_injection() const
bool ready_inject_tcl() const
const double & thetaC_tcl() const
const double & Rc_inject() const
classe YAML_data : collection of all needed information for data to save/restore in order to write th...