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>
76static void eval_vitesse(
double x,
double y,
double z,
double t,
78 double& vx,
double& vy,
double& vz)
80 int i0=0, i1=1, i2=2, i3=3;
99static void integrer_vitesse_imposee(
101 double temps,
double dt,
double& x,
double& y,
double& z)
109 parser_vx,parser_vy,parser_vz,vx,vy,vz);
111 eval_vitesse(x + vx * dt * 0.5,
115 parser_vx,parser_vy,parser_vz,vx1,vy1,vz1);
117 eval_vitesse(x + vx1 * dt * 0.5,
121 parser_vx,parser_vy,parser_vz,vx2,vy2,vz2);
123 eval_vitesse(x + vx2 * dt,
127 parser_vx,parser_vy,parser_vz,vx3,vy3,vz3);
129 x += (vx + 2.*vx1 + 2.*vx2 + vx3) / 6. * dt;
130 y += (vy + 2.*vy1 + 2.*vy2 + vy3) / 6. * dt;
131 z += (vz + 2.*vz1 + 2.*vz2 + vz3) / 6. * dt;
143static void calculer_normale_sommets_interface(
const Maillage_FT_Disc& maillage,
151 const IntTab& facettes = maillage.
facettes();
152 const int nsommets_faces = facettes.
line_size();
154 normale.
resize(nsom, dim);
156 double n[3]= {0,0,0};
158 for (
int i = 0; i < nfaces; i++)
165 const double surface = surface_facettes[i];
166 for (
int k = 0; k < dim; k++)
167 n[k] = normale_facettes(i,k) * surface;
169 for (
int j = 0; j < nsommets_faces; j++)
171 const int sommet = facettes(i,j);
172 for (
int k = 0; k < dim; k++)
173 normale(sommet, k) += n[k];
197 int index = fullname.
find(
".");
198 std::string fname = index ==-1 ? fullname.
getString() : fullname.
getString().substr(0, index);
199 Nom nom_fich_xyz(fname);
202 nom_fich_xyz +=
"_maillage_interface";
203 nom_fich_xyz +=
".xyz";
211 Cerr <<
"Transport_Interfaces_FT_Disc::init_save_file Backup file for the FT mesh has already been initialized!" << finl;
218 fic_front_sauv_->ouvrir(nom_fich_xyz);
225 if (!fic_front_sauv_)
227 Cerr <<
"Transport_Interfaces_FT_Disc::close_save_file Backup file for the FT mesh has not been initialized!" << finl;
233 fic_front_sauv_.valeur() <<
Nom(
"fin");
234 (fic_front_sauv_.valeur()).flush();
235 (fic_front_sauv_.valeur()).syncfile();
244 std::vector<YAML_data> data = indicatrice_cache->data_a_sauvegarder();
247 std::string name =
pb_name_.getString() +
"_" +
ident_.getString() +
"_indicatrice_cache_tag";
250 data.push_back(ictag);
260 bytes += indicatrice_cache->sauvegarder(os);
273 std::string name =
pb_name_.getString() +
"_" +
ident_.getString() +
"_indicatrice_cache_tag";
282 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;
290 xyz_os << mon_ident << finl;
299 xyz_os << mon_ident << finl;
316 Cerr <<
"Transport_Interfaces_FT_Disc_interne::reprendre" << finl;
322 Nom name =
pb_name_ +
"_" + indicatrice_cache->le_nom();
332 if (!indicatrice_cache)
335 indicatrice_cache.typer(type);
341 Cerr <<
"Transport_Interfaces_FT_Disc_interne::reprendre : WARNING ! PDI format is not supported for the FT mesh so we use xyz format" << finl;
347 xyz_is->ouvrir(nomfic);
350 Cerr <<
"Error during the opening of the restart file : " << nomfic << finl;
354 Cerr <<
"Reading interface in file " << nomfic << finl;
358 avancer_fichier(xyz_is, mon_ident);
363 avancer_fichier(xyz_is, mon_ident);
385 interpolation_repere_local_ = 0;
393Transport_Interfaces_FT_Disc::~Transport_Interfaces_FT_Disc()
395 delete variables_internes_;
462 param.
ajouter(
"VOFlike_correction_volume",&variables_internes_->VOFlike_correction_volume);
465 param.
ajouter(
"nb_lissage_correction_volume",&variables_internes_->nb_lissage_correction_volume);
467 param.
ajouter(
"nb_iterations_correction_volume",&variables_internes_->nb_iterations_correction_volume);
484 param.
ajouter(
"parcours_interface",&variables_internes_->parcours_interface_);
495 param.
ajouter(
"sous_domaine_volume_impose",&variables_internes_->nom_domaine_volume_impose_);
496 param.
ajouter_flag(
"interpolation_repere_local", &interpolation_repere_local_);
505 param.
ajouter(
"n_iterations_interpolation_ibc",&variables_internes_->n_iterations_interpolation_ibc);
513 param.
ajouter(
"nombre_facettes_retenues_par_cellule",&variables_internes_->nomb_fa7_accepted);
517 param.
ajouter(
"seuil_convergence_uzawa",&variables_internes_->seuil_uzawa) ;
521 param.
ajouter(
"nb_iteration_max_uzawa",&variables_internes_->nb_iter_uzawa) ;
528 param.
ajouter(
"temps_debut",&temps_debut_);
529 param.
ajouter(
"vitesse_imposee_regularisee", &variables_internes_->vimp_regul) ;
535 param.
ajouter(
"collision_model_fpi",&collision_model_);
545 if (un_mot==
"maillage")
551 else if (un_mot==
"methode_transport")
554 motcles2[0] =
"vitesse_imposee";
555 motcles2[1] =
"loi_horaire";
556 motcles2[2] =
"vitesse_interpolee";
559 Cerr << un_mot <<
" " << motlu << finl;
560 int rang = motcles2.
search(motlu);
567 <<
" analytical expressions vx(x,y,[z,]t) vy(...) [vz(...)]" << finl;
569 variables_internes_->methode_transport =
577 Nom& expression = variables_internes_->expression_vitesse_imposee[i];
580 Cerr <<
" Component " << i <<
" of the velocity : " << expression << finl;
586 Cerr <<
" Imposed motion by the time scheme : ";
590 Cerr << nom_loi << finl;
597 Cerr <<
" Transport velocity interpolated from the velocity \n"
598 <<
" of the following Navier_Stokes equation :" << finl;
600 variables_internes_->methode_transport =
605 Cerr << nom_eq << finl;
606 variables_internes_->refequation_vitesse_transport =
608 probleme_base_.valeur()).equation_hydraulique(nom_eq);
610 const Equation_base& eqn_hydraulique = variables_internes_->refequation_vitesse_transport.valeur();
623 Cerr <<
"Transport_Interfaces_FT_Disc::lire\n"
624 <<
"The options for methode_transport are :\n"
631 else if (un_mot==
"n_iterations_distance")
635 variables_internes_->n_iterations_distance = n;
637 Cerr <<
"Iterations number of the smoothing tool \"distance_interface\": " << n << finl;
640 else if (un_mot==
"iterations_correction_volume")
644 variables_internes_->iterations_correction_volume = n;
646 variables_internes_->VOFlike_correction_volume = 1;
647 variables_internes_->nb_lissage_correction_volume = n;
648 variables_internes_->nb_iterations_correction_volume = n;
651 Cerr <<
"Obsolete Keyword is read Iterations_correction_volume : " << n << finl;
652 Cerr <<
"For future compatibility, it is recommended to switch to the new syntax : " << finl;
654 Cerr <<
" VOFlike_correction_volume 1 # a flag (0 or 1) for activation #" << finl;
656 Cerr <<
" VOFlike_correction_volume 0 # a flag (0 or 1) for activation #" << finl;
657 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;
658 Cerr <<
" nb_iterations_correction_volume " << n <<
" # to iterate on the volume correction until seuil is reached. #" << finl;
662 Cerr <<
"The value of nb_iterations_correction_volume is read from bloc remaillage at keyword: nb_iter_correction_volume" << finl;
667 else if (un_mot==
"volume_impose_phase_1")
671 variables_internes_->volume_impose_phase_1 = x;
673 Cerr <<
"volume_impose_phase_1 : " << x << finl;
676 else if (un_mot==
"methode_interpolation_v")
679 motcles2[0] =
"valeur_a_elem";
680 motcles2[1] =
"vdf_lineaire";
681 motcles2[2] =
"mean_volumic_velocity";
685 Cerr << un_mot <<
" " << motlu << finl;
686 int rang = motcles2.
search(motlu);
690 variables_internes_->methode_interpolation_v =
694 variables_internes_->methode_interpolation_v =
698 variables_internes_->methode_interpolation_v =
700 if (interpolation_repere_local_)
702 Cerr <<
"ERROR : interpolation_repere_local should not be employed"
703 " with methode_interpolation_v::MEAN_VOLUMIC_VELOCITY !!!!!!!!!!!" << finl;
709 Cerr <<
"Transport_Interfaces_FT_Disc::lire\n"
710 <<
"The options for " << un_mot <<
" are :\n"
716 else if (un_mot==
"parcours_interface")
718 is >> variables_internes_->parcours_interface_;
721 else if (un_mot==
"injecteur_interfaces")
727 Cerr <<
"Reading of the interfaces to inject in the file " << nom_fichier << finl;
729 ArrOfDouble& temps = variables_internes_->injection_interfaces_temps_;
730 ArrOfInt& phase = variables_internes_->injection_interfaces_phase_;
731 Noms& expr = variables_internes_->injection_interfaces_expressions_;
739 temps.append_array(t);
745 Cerr <<
"Time=" << t <<
" phase=" << ph <<
" expression=" << une_expr << finl;
748 envoyer_broadcast(variables_internes_->injection_interfaces_temps_, 0);
749 envoyer_broadcast(variables_internes_->injection_interfaces_phase_, 0);
750 envoyer_broadcast(variables_internes_->injection_interfaces_expressions_, 0);
754 else if (un_mot==
"interpolation_champ_face")
757 motcles2[0] =
"base";
758 motcles2[1] =
"lineaire";
761 int rang = motcles2.
search(motbis);
765 variables_internes_->interpolation_champ_face =
770 variables_internes_->interpolation_champ_face =
774 mots.add(
"vitesse_fluide_explicite");
775 mots.add(
"extrapolation_diphasique_solide");
776 mots.add(
"extrapolation_solide");
777 mots.add(
"distance_projete_face");
782 Motcle accouverte =
"{" , accfermee =
"}" ;
783 if (mot2 == accouverte)
787 while (mot2 != accfermee)
789 int rang2 = mots.
search(mot2);
794 variables_internes_->vf_explicite = 1 ;
795 Cerr <<
"Lecture du type de vitesse fluide pour le calcul de v_imp : vf_explicite " << finl;
800 variables_internes_->is_extra_diphasique_solide = 1 ;
801 variables_internes_->is_extra_solide = 0 ;
802 Cerr <<
"Lecture du type d'extrapolation considere : diphasique_solide " << finl;
807 variables_internes_->is_extra_solide = 1 ;
808 variables_internes_->is_extra_diphasique_solide = 1 ;
809 Cerr <<
"Lecture du type d'extrapolation considere : solide " << finl;
814 variables_internes_->is_distance_projete_face = 1 ;
815 Cerr <<
"Lecture du type de distance consideree dans l'interpolation : distance reelle entre projete et face " << finl;
824 Cerr <<
"Erreur a la lecture des parametres de l'interpolation lineaire : " << finl;
825 Cerr <<
"On attendait : " << accouverte <<
" , on a eu " << mot2 << finl;
831 Cerr <<
"Transport_Interfaces_FT_Disc::lire\n"
832 <<
" les options de interpolation_champ_face sont :\n"
837 else if (un_mot==
"type_vitesse_imposee")
840 motcles2[0] =
"uniforme";
841 motcles2[1] =
"analytique";
844 int rang = motcles2.
search(motbis);
854 Cerr <<
"Transport_Interfaces_FT_Disc::lire\n"
855 <<
" les options de type_vitesse_imposee sont :\n"
860 else if (un_mot==
"distance_IBC_faces")
863 motcles2[0] =
"initiale";
864 motcles2[1] =
"modifiee";
868 int rang = motcles2.
search(motbis);
880 Cerr <<
"Transport_Interfaces_FT_Disc::lire\n"
881 <<
" les options de distance_IBC_faces sont :\n"
886 else if (un_mot==
"distance_projete_faces")
889 motcles2[0] =
"initiale";
890 motcles2[1] =
"modifiee";
891 motcles2[2] =
"simplifiee";
895 int rang = motcles2.
search(motbis);
908 Cerr <<
"Transport_Interfaces_FT_Disc::lire\n"
909 <<
" les options de distance_projete_faces sont :\n"
914 else if (un_mot==
"type_indic_faces")
917 motcles2[0] =
"standard";
918 motcles2[1] =
"modifiee";
919 motcles2[2] =
"ai_based";
922 Cerr << un_mot <<
" " << motlu << finl;
923 int rang = motcles2.
search(motlu);
930 Cerr <<
" Standard interpolation of indicatrice to faces." << finl;
937 <<
" analytical expressions vx(x,y,[z,]t) vy(...) [vz(...)]" << finl;
946 variables_internes_->modified_indic_faces_position = 0. ;
947 variables_internes_->modified_indic_faces_thickness= 1. ;
949 mots.add(
"position");
950 mots.add(
"thickness");
954 Motcle accouverte =
"{" , accfermee =
"}" ;
955 if (mot == accouverte)
959 while (mot != accfermee)
961 int rang2 = mots.
search(mot);
966 is >> variables_internes_->modified_indic_faces_position ;
967 Cerr <<
"Lecture de la position de l interface modifiee " << finl;
972 is >> variables_internes_->modified_indic_faces_thickness ;
973 if ( variables_internes_->modified_indic_faces_thickness < 0.)
975 Cerr <<
"L epaisseur de l interface doit etre positive ou nulle!!" << finl;
978 Cerr <<
"Lecture de l epaisseur de l interface modifiee " << finl;
987 Cerr <<
"Erreur a la lecture des parametres de l'indicatrice modifiee " << finl;
988 Cerr <<
"On attendait : " << accouverte << finl;
993 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;
1001 Cerr <<
" The interpolation of indicatrice to faces is based on the interfacial area"
1002 <<
" and on the normal to the interface." << finl;
1006 Cerr <<
"Transport_Interfaces_FT_Disc::lire\n"
1007 <<
" les options de type_indic_faces sont :\n"
1024 const Conds_lim& les_cl = le_dom_Cl_dis->les_conditions_limites();
1025 const int n = les_cl.size();
1027 for (i = 0; i < n; i++)
1034 Cerr <<
"Error in Transport_Interfaces_FT_Disc::verif_Cl():\n"
1035 <<
" Boundary conditions for Transport_Interfaces_FT_Disc must be\n"
1036 <<
" of type Paroi_FT_disc (example : \"Paroi_FT_disc symetrie\")" << finl;
1043int facettes_trier_retirer_doublons(IntTab& fa7)
1045 using line = std::array<int, _DIM_>;
1046 line* chunked_fa7 =
reinterpret_cast<line*
>(fa7.
addr());
1049 auto same_fa7 = [&](line A, line B)
1052 std::sort(A.begin(), A.end());
1053 std::sort(B.begin(), B.end());
1055 for (
int i = 0; i < dim; i++)
1065 std::sort(chunked_fa7, chunked_fa7 + nb_facettes, [&](
const line& fa1,
const line& fa2)
1067 return (same_fa7(fa1,fa2) < 0);
1072 for (
int i = 1; i < nb_facettes; i++)
1074 line& fa1 = chunked_fa7[i];
1075 line& fa2 = chunked_fa7[i-1];
1076 if(same_fa7(fa1,fa2)!=0)
1079 for(
int j=0; j<dim; j++)
1080 fa7(count,j)=fa7(i,j);
1092 ArrOfInt phase_specifiee;
1094 int default_phase = -1;
1097 int reverse_normal = 0;
1103 Cerr <<
"Error in Transport_Interfaces_FT_Disc::lire_cond_init, expected { after fichier_geom" << finl;
1108 motcles[0] =
"fichier_geom";
1109 motcles[1] =
"point_phase";
1110 motcles[2] =
"default_phase";
1111 motcles[3] =
"lata_dump";
1112 motcles[4] =
"nom_domaine";
1113 motcles[5] =
"reverse_normal";
1119 const int rang = motcles.
search(motlu);
1124 Cerr <<
" Reading .geom file in following file: " << filename << finl;
1131 is >> phase_specifiee[i];
1132 Cerr <<
" Reading point in phase " << phase_specifiee[i];
1133 if (phase_specifiee[i] != 0 && phase_specifiee[i] != 1)
1135 Cerr <<
" Error reading point_phase : expected 0 or 1" << finl;
1141 Cerr <<
" " << points(i,j);
1147 is >> default_phase;
1148 if (default_phase != 0 && default_phase != 1)
1150 Cerr <<
" Error reading default_phase : expected 0 or 1" << finl;
1153 Cerr <<
" Default phase : " << default_phase << finl;
1157 Cerr <<
" Dumping connex components and indicator function in lata file : "
1158 << lata_file << finl;
1162 Cerr <<
" Reading interface data in domain: " << domain_name << finl;
1166 Cerr <<
" Reverse_normal : swap nodes to reverse the normal vector of the surface mesh" << finl;
1169 Cerr <<
" Unknown keyword " << motlu
1170 <<
"\n Known keywords are " << motcles << finl;
1176 if (filename !=
"??" && domain_name !=
"??")
1178 Cerr <<
"Error: you specified both a .geom file name AND a domain name" << finl;
1181 else if (filename ==
"??" && domain_name ==
"??")
1183 Cerr <<
"Error: you must specify a FICHIER_GEOM or a NOM_DOMAINE" << finl;
1186 if (filename !=
"??")
1188 Cerr <<
"Reading geometry file and building interface mesh" << finl;
1189 if (is_a_binary_file(filename))
1203 Cerr <<
"Reading interface in existing domain: " << domain_name << finl;
1206 Cerr <<
"Error : object " << domain_name <<
" is not a domain" << finl;
1213 Cerr <<
"Reversing mesh normal vectors" << finl;
1214 Domaine& domaine = ref_dom.valeur();
1215 IntTab& elems = domaine.les_elems();
1216 const int nb_elem = elems.dimension(0);
1217 const int nb_som_elem = elems.line_size();
1218 if (nb_som_elem != 2 && nb_som_elem != 3)
1220 Cerr <<
"Error: mesh has wrong dimension (must be segments in 2D, triangles in 3D)" << finl;
1223 const int j0 = nb_som_elem - 2;
1224 for (
int i = 0; i < nb_elem; i++)
1226 const int s0 = elems(i, j0);
1227 const int s1 = elems(i, j0+1);
1229 elems(i, j0+1) = s0;
1238 Domaine& domaine = ref_dom.valeur();
1239 IntTab& fa7 = domaine.les_elems();
1242 const int nb_facettes = fa7.
dimension(0);
1247 count = facettes_trier_retirer_doublons<2>(fa7);
1249 count = facettes_trier_retirer_doublons<3>(fa7);
1252 Cerr<<
"End correction SM 12/08 the removed faces number is: "<< nb_facettes <<
" - "<< count+1 <<
" = " << nb_facettes-count-1 << finl;
1257 Cerr <<
"Building interface mesh" << finl;
1261 Cerr <<
"Extracting connex components and assigning indicator function." << finl;
1270 for (
int i = 0; i < nb_elem; i++)
1273 if (index_elem[i] >= 0)
1279 const IntTab& elem_faces = domaine_vf.
elem_faces();
1281 const int nb_local_connex_components = search_connex_components_local(elem_faces, faces_elem, num_compo);
1282 const int nb_connex_components = compute_global_connex_components(num_compo, nb_local_connex_components);
1285 Cerr <<
" found " << nb_connex_components <<
" connex components" << finl;
1287 ArrOfInt phase_of_component(nb_connex_components);
1288 phase_of_component= default_phase;
1292 const int nb_pts = phase_specifiee.
size_array();
1295 ArrOfInt elem_points;
1297 for (
int i = 0; i < nb_pts; i++)
1299 int composante_connexe = -2;
1300 if (elem_points[i] >= 0)
1301 composante_connexe = num_compo[elem_points[i]];
1304 composante_connexe = (int)
mp_max(composante_connexe);
1305 if (composante_connexe == -2)
1307 Cerr <<
"Error : point " << i <<
" is not in the domain" << finl;
1311 else if (composante_connexe == -1)
1313 Cerr <<
"Error : point " << i <<
" is in connex_component -1"
1314 <<
" (this point is in a mesh cell containing the interface)."
1315 <<
" Use the dump_lata keyword to find a point away from the interface" << finl;
1321 const int phase = phase_specifiee[i];
1322 Cerr <<
"Assigning phase " << phase <<
" to component " << composante_connexe << finl;
1323 phase_of_component[composante_connexe] = phase;
1328 DoubleTab& indic = variables_internes_->indicatrice_cache->valeurs();
1329 for (
int i = 0; i < nb_elem; i++)
1331 const int compo = num_compo[i];
1334 const int phase = phase_of_component[compo];
1341 if (lata_file !=
"??")
1343 Cerr <<
"Writing lata file" << finl;
1345 const Domaine& un_dom = domaine_vf.
domaine();
1346 constexpr double TEMPS = 0.;
1347 constexpr int FIRST_POST = 1;
1352 DoubleTab data(nb_elem);
1353 for (
int i = 0; i < nb_elem; i++)
1354 data(i) = num_compo(i);
1361 Nom nom_champ(
"connex_component");
1362 lata.
ecrire_champ(un_dom, unites, noms_compo, 1, TEMPS, nom_champ, nom_dom,
"elem",
"scalar",
1364 nom_champ =
"indicatrice";
1365 lata.
ecrire_champ(un_dom, unites, noms_compo, 1, TEMPS, nom_champ, nom_dom,
"elem",
"scalar",
1367 nom_champ =
"distance";
1368 lata.
ecrire_champ(un_dom, unites, noms_compo, 1, TEMPS, nom_champ, nom_dom,
"elem",
"scalar",
1375 if (phase_of_component.
size_array() > 0 && min_array(phase_of_component) < 0)
1377 Cerr <<
"Error: some connex components of the indicator function were not initialized\n"
1378 <<
" you must either specify more points or specify a default_phase.\n"
1379 <<
" With the lata_dump keyword, you can watch connex components and find where\n"
1380 <<
" points should be added." << finl;
1446 Cerr <<
"Reading initial condition" << finl;
1448 motcles[0] =
"fonction";
1449 motcles[1] =
"fichier_geom";
1450 motcles[2] =
"fonction_ignorer_collision";
1451 motcles[3] =
"reprise";
1456 Cerr <<
"Erreur dans Transport_Interfaces_FT_Disc::lire_cond_init\n";
1457 Cerr <<
" On attendait {\n On a trouve " << motlu << finl;
1460 const Motcle virgule =
",";
1468 int rang = motcles.
search(motlu);
1484 if (nom_phase==
"ajout_phase1")
1488 else if (nom_phase==
"ajout_phase0")
1494 Cerr <<
"Transport_Interfaces_FT::lire_cond_init :\n";
1495 Cerr <<
" The keyword " << nom_phase <<
" is not understood.\n";
1496 Cerr <<
" Allowed keywords are : ajout_phase0 or ajout_phase1" << finl;
1503 if (
probleme().reprise_effectuee())
1505 Cerr <<
" Interface not build since a restarting is expected." << finl;
1509 Cerr <<
" Interface construction : " << expression << finl;
1514 int ignorer_collision = (rang==2);
1516 variables_internes_->indicatrice_cache->valeurs(),
1518 variables_internes_->distance_interface_sommets,
1522 if (ignorer_collision)
1523 Cerr <<
"Warning: the interface is colliding with an existing interface (ignorer_collision=1)." << finl;
1526 Cerr <<
"Error: the interface is colliding with an existing interface." << finl;
1546 Cerr <<
"A .xyz file is waited to restart the calculation for the interface equation." << finl;
1554 is2 >> tmp >> version;
1560 Cerr <<
"Transport_Interfaces_FT::lire_cond_init :\n";
1561 Cerr <<
" The keyword " << motlu <<
" is not understood here.\n";
1562 Cerr <<
" Allowed keywords are :\n" << motcles << finl;
1566 Cerr<<
"lectureF "<<motlu<<finl;
1568 while (motlu==virgule);
1572 Cerr <<
"Error for method Transport_Interfaces_FT_Disc::lire_cond_init\n";
1573 Cerr <<
" A } was expected\n It has been found " << motlu << finl;
1582 assert(!probleme_base_);
1586 Cerr <<
"Error for method Transport_Interfaces_FT_Disc::associer_pb_base\n";
1587 Cerr <<
" The Transport_Interfaces_FT_Disc equation must be associated to a problem\n";
1588 Cerr <<
" of type Probleme_FT_Disc_gen" << finl;
1589 Cerr <<
" or type Pb_Thermohydraulique_Especes_QC" << finl;
1590 Cerr <<
" or type Pb_Thermohydraulique_Especes_Turbulent_QC" << finl;
1593 probleme_base_ = un_probleme;
1597 variables_internes_->set_pb_name(pb_name);
1611 variables_internes_->set_ident(
le_nom());
1619 fieldname =
"INDICATRICE";
1620 fieldname += suffix;
1623 1 , nb_valeurs_temps,
1626 indicatrice_->associer_eqn(*
this);
1631 fieldname =
"INDICATRICE_CACHE";
1632 fieldname += suffix;
1637 variables_internes_->indicatrice_cache);
1638 variables_internes_->indicatrice_cache->associer_eqn(*
this);
1640 variables_internes_->indicatrice_cache->PDI_save_type(
true);
1644 fieldname =
"INDICATRICE_FACES";
1645 fieldname += suffix;
1648 1 , nb_valeurs_temps,
1650 indicatrice_faces_);
1651 indicatrice_faces_->associer_eqn(*
this);
1654 fieldname =
"VITESSE_FILTREE";
1655 fieldname += suffix;
1660 variables_internes_->vitesse_filtree);
1661 variables_internes_->vitesse_filtree->associer_eqn(*
this);
1665 fieldname =
"FLUX_TMP";
1666 fieldname += suffix;
1671 variables_internes_->tmp_flux);
1675 fieldname =
"INDEX_ELEM";
1680 variables_internes_->index_element);
1683 fieldname =
"NELEM_PAR_DIRECTION";
1688 variables_internes_->nelem_par_direction);
1689 champs_compris_.ajoute_champ(variables_internes_->nelem_par_direction);
1691 fieldname =
"DISTANCE_INTERFACE_ELEM";
1692 fieldname += suffix;
1697 variables_internes_->distance_interface);
1698 champs_compris_.ajoute_champ(variables_internes_->distance_interface);
1700 fieldname =
"DISTANCE_INTERFACE_FACE";
1701 fieldname += suffix;
1706 variables_internes_->distance_interface_faces);
1707 champs_compris_.ajoute_champ(variables_internes_->distance_interface_faces);
1709 fieldname =
"DISTANCE_INTERFACE_FACE_COR";
1710 fieldname += suffix;
1715 variables_internes_->distance_interface_faces_corrigee);
1716 champs_compris_.ajoute_champ(variables_internes_->distance_interface_faces_corrigee);
1718 fieldname =
"DISTANCE_INTERFACE_FACE_DIF";
1719 fieldname += suffix;
1724 variables_internes_->distance_interface_faces_difference);
1725 champs_compris_.ajoute_champ(variables_internes_->distance_interface_faces_difference);
1727 fieldname =
"VITESSE_IMP_INTERP";
1728 fieldname += suffix;
1733 vitesse_imp_interp_);
1740 DoubleTab& d = variables_internes_->distance_interface_sommets;
1745 variables_internes_->distance_interface_sommets = 0.;
1747 fieldname =
"NORMALE_INTERFACE";
1748 fieldname += suffix;
1753 variables_internes_->normale_interface);
1758 fieldname =
"SURFACE_INTERFACE";
1759 fieldname += suffix;
1764 variables_internes_->surface_interface);
1771 fieldname =
"TIME_DERIVATIVE_INTERFACE";
1772 fieldname += suffix;
1791 variables_internes_->connectivite_frontieres_.associer_domaine_vf(domaine_vf);
1805 variables_internes_->algorithmes_transport_.typer(
"Algorithmes_Transport_FT_Disc");
1809 Cerr <<
"Discretizing of the boundary conditions... ";
1812 Cerr <<
" OK " << finl;
1814 le_dom_Cl_dis->associer_eqn(*
this);
1815 le_dom_Cl_dis->associer_inconnue(
inconnue());
1827 Journal() <<
"Transport_Interfaces_FT_Disc::remailler_interface " <<
le_nom() << finl;
1831 Champ_base& indicatrice = variables_internes_->indicatrice_cache.valeur();
1840 if (not(preparer_calcul_anticipated_done_))
1846 Process::Journal()<<
"Transport_Interfaces_FT_Disc::preparer_calcul does nothing because it has been anticipated in the problem"<<finl;
1855 preparer_calcul_anticipated_done_ =
true;
1859 compute_nb_particles_tot();
1874 if (collision_model_)
1879 collision_model_.valeur().preparer_calcul(domain_vdf, nb_particles_tot_,
1880 ns, *
this, schema_comm_FT);
1886 le_dom_Cl_dis->initialiser(temps);
1888 if (
probleme().reprise_effectuee())
1890 Cerr <<
"Calculation restarting : no initial meshing" << finl;
1894 Cerr <<
"Initial condition remeshing." << finl;
1901 indicatrice_->changer_temps(temps);
1902 variables_internes_->indicatrice_cache->changer_temps(temps);
1903 indicatrice_faces_->changer_temps(temps);
1904 variables_internes_->vitesse_filtree->changer_temps(temps);
1905 variables_internes_->tmp_flux->changer_temps(temps);
1906 variables_internes_->index_element->changer_temps(temps);
1907 variables_internes_->nelem_par_direction->changer_temps(temps);
1908 variables_internes_->distance_interface->changer_temps(temps);
1909 variables_internes_->distance_interface_faces->changer_temps(temps);
1910 variables_internes_->distance_interface_faces_corrigee->changer_temps(temps);
1911 variables_internes_->distance_interface_faces_difference->changer_temps(temps);
1912 vitesse_imp_interp_->changer_temps(temps);
1913 variables_internes_->normale_interface->changer_temps(temps);
1914 variables_internes_->surface_interface->changer_temps(temps);
1917 variables_internes_->set_checkpoint_fname(checkpoint_fname);
1927 Cerr <<
"Error for the method Transport_Interfaces_FT_Disc::preparer_calcul\n"
1928 <<
"The transport method has not indicated for the equation\n"
1947 variables_internes_->init_save_file();
1952 variables_internes_->close_save_file();
2003 if (tag == variables_internes_->distance_normale_cache_tag)
2005 Cerr <<
"WARNING : Unneeded call to Transport_Interfaces_FT_Disc::update_normale_distance_interface" << finl;
2006 Cerr <<
" mesh tag: " << tag <<
", cache tag: "<< variables_internes_->distance_normale_cache_tag << finl;
2010 DoubleTab& distance = variables_internes_->distance_interface->valeurs();
2011 DoubleTab& normale = variables_internes_->normale_interface->valeurs();
2015 variables_internes_->n_iterations_distance);
2016 variables_internes_->distance_normale_cache_tag = tag;
2031 if (tag == variables_internes_->indicatrice_cache_tag)
2033 Cerr <<
"WARNING : Unneeded call to Transport_Interfaces_FT_Disc::update_indicatrice. tag = " << tag << finl;
2034 Cerr <<
" indicatrice is already up to date" << finl;
2039 if (tag != variables_internes_->distance_normale_cache_tag)
2041 Cerr <<
"ERROR : in Transport_Interfaces_FT_Disc::update_indicatrice : distance_normale not up to date" << finl;
2042 Cerr <<
" Call update_indicatrice_normale_distance before or update_normale_distance_interface to do both at the same time" << finl;
2043 Cerr <<
" mesh tag: " << tag <<
", cache tag: "<< variables_internes_->distance_normale_cache_tag << finl;
2050 Cerr <<
"ERROR : in Transport_Interfaces_FT_Disc::update_indicatrice : maillage pas parcouru" << finl;
2054 DoubleVect& valeurs_indicatrice = variables_internes_->indicatrice_cache->valeurs();
2058 variables_internes_->indicatrice_cache_tag = tag;
2073 if (tag != variables_internes_->indicatrice_cache_tag)
2075 Cerr <<
"Error: indicatrice not up to date" << finl;
2089 if (tag != variables_internes_->indicatrice_cache_tag)
2091 Cerr <<
"Error: indicatrice not up to date in get_indicatrice. " << finl;
2094 return variables_internes_->indicatrice_cache.valeur();
2100 if (tag != variables_internes_->distance_normale_cache_tag)
2102 Cerr <<
"Error: indicatrice not up to date in get_distance_interface. " << finl;
2105 return variables_internes_->distance_interface.valeur();
2111 if (tag != variables_internes_->distance_normale_cache_tag)
2113 Cerr <<
"Error: indicatrice not up to date in get_normale_interface. " << finl;
2116 return variables_internes_->normale_interface.valeur();
2130void interpoler_vitesse_point_vdf(
const Champ_base& champ_vitesse,
2131 const FTd_vecteur3& coord_som,
2133 FTd_vecteur3& vitesse,
2134 const int future_or_past)
2136 const DoubleTab& valeurs_v = (bool)(future_or_past) ? champ_vitesse.
futur() : champ_vitesse.
valeurs();
2138 const IntTab& elem_faces = domaine_vf.
elem_faces();
2139 const DoubleTab& xv = domaine_vf.
xv();
2140 const DoubleTab& xp = domaine_vf.
xp();
2141 const IntTab& face_voisins = domaine_vf.
face_voisins();
2170 for (compo = 0; compo < dim; compo++)
2179 int elem_voisins[3] = { -1, -1, -1 };
2181 int indice_local_face_voisine[3] = { -1, -1, -1 };
2185 int elem_diagonal = -1;
2189 for (direction = 0; direction < dim; direction++)
2192 const double x = coord_som[direction];
2194 if (direction == compo)
2199 const int face_inf_compo = elem_faces(element, compo);
2200 const int face_sup_compo = elem_faces(element, compo + dim);
2201 const double xmin = xv(face_inf_compo, direction);
2202 const double xmax = xv(face_sup_compo, direction);
2203 a = (xmax - x) / (xmax - xmin);
2211 const double centre_elem = xp(element, direction);
2214 const int i_face_voisine = direction + ((x < centre_elem) ? 0 : dim);
2216 const int face_voisine = elem_faces(element, i_face_voisine);
2218 const int elem_voisin =
2219 face_voisins(face_voisine, 0) + face_voisins(face_voisine, 1) - element;
2220 if (elem_voisin >= 0)
2223 const double centre_voisin = xp(elem_voisin, direction);
2224 a = (centre_voisin - x) / (centre_voisin - centre_elem);
2225 elem_voisins[direction] = elem_voisin;
2226 indice_local_face_voisine[direction] = i_face_voisine;
2239 coef[direction] = a;
2244 const int direction1 = (compo + 1) % 3;
2245 const int direction2 = (compo + 2) % 3;
2246 const int elem_voisin1 = elem_voisins[direction1];
2247 const int elem_voisin2 = elem_voisins[direction2];
2249 if (elem_voisin1 >= 0 && elem_voisin2 >= 0)
2252 int i_face_voisine, face_voisine, elem_diagonal1, elem_diagonal2;
2254 i_face_voisine = indice_local_face_voisine[direction2];
2255 face_voisine = elem_faces(elem_voisin1, i_face_voisine);
2257 elem_diagonal1 = face_voisins(face_voisine, 0)
2258 + face_voisins(face_voisine, 1) - elem_voisin1;
2259 if (elem_diagonal1 >= 0)
2263 i_face_voisine = indice_local_face_voisine[direction1];
2264 face_voisine = elem_faces(elem_voisin2, i_face_voisine);
2266 elem_diagonal2 = face_voisins(face_voisine, 0)
2267 + face_voisins(face_voisine, 1) - elem_voisin2;
2268 if (elem_diagonal1 == elem_diagonal2)
2271 elem_diagonal = elem_diagonal1;
2285 if (elem_diagonal < 0)
2288 elem_voisins[direction1] = -1;
2289 elem_voisins[direction2] = -1;
2301 const int direction1 = 1 - compo;
2302 int element_voisin = elem_voisins[direction1];
2303 if (element_voisin < 0)
2304 element_voisin = element;
2305 const int f00 = elem_faces(element, compo);
2306 const int f10 = elem_faces(element, compo + dim);
2307 const int f01 = elem_faces(element_voisin, compo);
2308 const int f11 = elem_faces(element_voisin, compo + dim);
2310 double ci = coef[compo];
2312 double cj = coef[direction1];
2314 ci * cj * valeurs_v(f00)
2315 + (1.-ci) * cj * valeurs_v(f10)
2316 + ci * (1.-cj) * valeurs_v(f01)
2317 + (1.-ci) * (1.-cj) * valeurs_v(f11);
2319 && (xv(f00,0) <DMINFLOAT)
2320 && ((fabs(valeurs_v(f00)-valeurs_v(f10))>DMINFLOAT) || (fabs(valeurs_v(f01)-valeurs_v(f11))>DMINFLOAT)))
2322 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;
2323 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;
2324 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;
2325 const double x = coord_som[0];
2326 Cerr <<
"Here, the difference is "<< (1.-ci) <<
" vs. " << xv(f10, 0)/x << finl;
2327 Cerr <<
"u1= " << valeurs_v(f10) <<
" direction1 : " << direction1 << finl;
2328 Cerr <<
"interpoler_vitesse_point_vdf of Transport_Interface..cpp not exiting but interpolation is adapted" << finl;
2329 Cerr <<
"Former: " << vitesse[compo];
2330 vitesse[compo] = xv(f10, 0)/x * (
2332 + (1.-cj) * valeurs_v(f11));
2333 Cerr <<
" New: " << vitesse[compo] << finl;
2340 const int direction1 = (compo + 1) % 3;
2341 const int direction2 = (compo + 2) % 3;
2342 int elem00 = element;
2343 int elem10 = elem_voisins[direction1];
2344 int elem01 = elem_voisins[direction2];
2345 int elem11 = elem_diagonal;
2346 if (elem10 < 0) elem10 = element;
2347 if (elem01 < 0) elem01 = element;
2348 if (elem11 < 0) elem11 = element;
2350 const int f000 = elem_faces(elem00, compo);
2351 const int f100 = elem_faces(elem00, compo + dim);
2352 const int f010 = elem_faces(elem10, compo);
2353 const int f110 = elem_faces(elem10, compo + dim);
2354 const int f001 = elem_faces(elem01, compo);
2355 const int f101 = elem_faces(elem01, compo + dim);
2356 const int f011 = elem_faces(elem11, compo);
2357 const int f111 = elem_faces(elem11, compo + dim);
2359 double ci = coef[compo];
2360 double cj = coef[direction1];
2361 double ck = coef[direction2];
2364 ( ci * cj * valeurs_v(f000)
2365 + (1.-ci) * cj * valeurs_v(f100)
2366 + ci * (1.-cj) * valeurs_v(f010)
2367 + (1.-ci) * (1.-cj) * valeurs_v(f110)) * ck
2368 +( ci * cj * valeurs_v(f001)
2369 + (1.-ci) * cj * valeurs_v(f101)
2370 + ci * (1.-cj) * valeurs_v(f011)
2371 + (1.-ci) * (1.-cj) * valeurs_v(f111)) * (1. - ck);
2394void interpoler_simple_vitesse_point_vdf(
const Champ_base& champ_vitesse,
2395 const FTd_vecteur3& coord_som,
2397 FTd_vecteur3& vitesse,
2398 const int future_or_past)
2400 const DoubleTab& valeurs_v = (bool)(future_or_past) ? champ_vitesse.
futur() : champ_vitesse.
valeurs();
2402 const IntTab& elem_faces = domaine_vf.
elem_faces();
2403 const DoubleTab& xv = domaine_vf.
xv();
2425 for (compo = 0; compo < dim; compo++)
2433 int direction = compo;
2435 const double x = coord_som[direction];
2439 const int face_inf_compo = elem_faces(element, compo);
2440 const int face_sup_compo = elem_faces(element, compo + dim);
2441 const double xmin = xv(face_inf_compo, direction);
2442 const double xmax = xv(face_sup_compo, direction);
2443 coef = (xmax - x) / (xmax - xmin);
2450 vitesse[compo] = coef * valeurs_v(face_inf_compo) + (1.-coef) * valeurs_v(face_sup_compo);
2457 const DoubleVect& volumes = domaine_vf.
volumes();
2458 int nd, nb_nd = indicatrice.
size();
2459 assert(nb_nd==domaine_vf.
nb_elem());
2461 double integrale = 0.;
2463 for (nd=0 ; nd<nb_nd ; nd++)
2465 integrale += indicatrice(nd) * volumes(nd);
2466 integrale_ph0 += (1 - indicatrice(nd))*volumes(nd);
2488 DoubleTab& vitesse_noeuds,
2491 const bool la_roue_de_vitesse_a_deja_tournee)
const
2493 switch(variables_internes_->methode_interpolation_v)
2502 Champ_base& champ_filtre = variables_internes_->vitesse_filtree.valeur();
2503 champ_vitesse_interp = & champ_filtre;
2512 const DoubleTab& val_champ_vitesse = (bool)(explicit_u_NS_) ? champ_vitesse.
futur() : champ_vitesse.
valeurs();
2513 champ_filtre.
valeurs() = val_champ_vitesse;
2519 if (explicit_u_NS_ && la_roue_de_vitesse_a_deja_tournee) champ_filtre.
futur() = champ_filtre.
valeurs();
2525 champ_vitesse_interp = &champ_vitesse;
2529 Cerr <<
"Error for the method Transport_Interfaces_FT_Disc::calculer_vitesse_transport_interpolee\n";
2530 Cerr <<
"The interpolation from a field " << champ_vitesse.
que_suis_je();
2531 Cerr <<
"is not developped." << finl;
2538 DoubleTab& les_positions = variables_internes_->doubletab_pos;
2539 IntVect& les_elements = variables_internes_->intvect_elements;
2540 DoubleTab& les_vitesses = variables_internes_->doubletab_vitesses;
2542 const DoubleTab& pos = maillage.
sommets();
2544 const int nb_pos_tot = pos.
dimension(0);
2546 les_positions.
resize(nb_pos_tot, dim);
2547 les_elements.
resize(nb_pos_tot);
2549 int nb_positions = 0;
2550 for (i = 0; i < nb_pos_tot; i++)
2552 const int num_elem = elem[i];
2555 for (j = 0; j < dim; j++)
2556 les_positions(nb_positions, j) = pos(i, j);
2557 les_elements(nb_positions) = num_elem;
2561 les_positions.
resize(nb_positions, dim);
2562 les_elements.
resize(nb_positions);
2563 les_vitesses.
resize(nb_positions, dim);
2565 if (explicit_u_NS_ && la_roue_de_vitesse_a_deja_tournee)
2569 champ_vitesse_interp->
valeur_aux_elems(les_positions, les_elements, les_vitesses);
2572 vitesse_noeuds.
resize(nb_pos_tot, dim);
2574 for (i = 0; i < nb_pos_tot; i++)
2578 for (j = 0; j < dim; j++)
2579 vitesse_noeuds(i, j) = les_vitesses(nb_positions, j);
2591 Cerr <<
"Error for the method Transport_Interfaces_FT_Disc::calculer_vitesse_transport_interpolee\n"
2592 <<
"the interpolation VDF_LINEAIRE is valid only for a VDF discretization with a Champ_face field\n"
2593 <<
" (type for the current field: " << champ_vitesse.
que_suis_je() << finl;
2596 const DoubleTab& pos = maillage.
sommets();
2598 const int nb_pos_tot = pos.
dimension(0);
2602 FTd_vecteur3 vitesse;
2603 vitesse_noeuds.
resize(nb_pos_tot, dim);
2604 for (i = 0; i < nb_pos_tot; i++)
2606 const int element = elem[i];
2610 for (j = 0; j < dim; j++)
2611 coord[j] = pos(i,j);
2615 interpoler_vitesse_point_vdf(champ_vitesse, coord, element, vitesse, explicit_u_NS_ && la_roue_de_vitesse_a_deja_tournee);
2620 interpoler_simple_vitesse_point_vdf(champ_vitesse, coord, element, vitesse, explicit_u_NS_ && la_roue_de_vitesse_a_deja_tournee);
2622 for (j = 0; j < dim; j++)
2623 vitesse_noeuds(i,j) = vitesse[j];
2631 const DoubleTab& vertices = maillage.
sommets();
2632 const int nb_vertices_tot = vertices.
dimension(0);
2635 ArrOfInt id_number_fa7(nb_fa7);
2637 int n = search_connex_components_local_FT(maillage, id_number_fa7);
2638 const int nb_particles_tot=compute_global_connex_components_FT(maillage, id_number_fa7, n);
2639 const DoubleTab& phase_indicator_function = indicatrice_->valeurs();
2640 const ArrOfInt& sommets_elem = maillage.
sommet_elem();
2641 const IntTab& facettes = maillage.
facettes();
2642 DoubleTab Particles_velocity(nb_particles_tot,
dimension);
2643 Particles_velocity=0;
2644 DoubleVect V_compo_elem(nb_particles_tot);
2646 ArrOfDouble Particles_surfaces(nb_particles_tot);
2647 Particles_surfaces=0;
2648 Equation_base& eqn_hydraulique = variables_internes_->refequation_vitesse_transport.valeur();
2650 const DoubleTab& tab_velocity=champ_vitesse.
valeurs();
2653 const DoubleVect& volumes_maille = domain_vdf.
volumes();
2654 const IntTab& elem_faces = domain_vdf.
elem_faces();
2673 for (
int elem=0; elem<domain_vdf.
nb_elem(); elem++)
2677 if (phase_indicator_function(elem)==id_solid_phase)
2679 const int compo= particles_eulerian_id_number(elem);
2680 V_compo_elem(compo)+=volumes_maille(elem)*(1-phase_indicator_function(elem));
2681 for (
int dim=0; dim<
dimension; dim++) Particles_velocity(compo,dim)+=0.5*
2682 (tab_velocity(elem_faces(elem,dim))
2683 +tab_velocity(elem_faces(elem,dim+
dimension))
2684 )*volumes_maille(elem);
2690 DoubleVect s_vcompo;
2692 tab_divide_any_shape(Particles_velocity, s_vcompo);
2696 DoubleTab Positions_compo(nb_particles_tot_,
dimension);
2698 for (
int fa7=0; fa7<nb_fa7; fa7++)
2702 const int compo = id_number_fa7(fa7);
2703 const double s_fa7 = surface_fa7(fa7);
2704 Particles_surfaces(compo)+=s_fa7;
2707 for (
int k = 0; k < vertices.
dimension(1); k++)
2709 int som = facettes(fa7, k);
2710 Positions_compo(compo, dim) += s_fa7 * vertices(som, dim)/
dimension;
2717 DoubleVect s_scompo;
2719 tab_divide_any_shape(Positions_compo, s_scompo);
2725 IntVect id_number_lagrangian_vertices;
2727 RESIZE_OPTIONS::NOCOPY_NOINIT);
2728 id_number_lagrangian_vertices = -1;
2730 const int dim = vitesse_noeuds.
dimension(1);
2731 for (
int iface = 0; iface < nb_fa7; iface++)
2733 const int id_number = id_number_fa7[iface];
2734 for (
int j = 0; j < dim; j++)
2735 id_number_lagrangian_vertices[facettes(iface, j)] = id_number;
2742 for (
int som = 0; som < nb_vertices_tot; som++)
2744 if (sommets_elem[som] >= 0)
2746 for (
int dim = 0; dim <
dimension; dim++)
2747 vitesse_noeuds(som, dim) =
2748 Particles_velocity(id_number_lagrangian_vertices(som), dim);
2757 Cerr <<
"Transport_Interfaces_FT_Disc::calculer_vitesse_transport_interpolee\n"
2758 <<
" interpolation method not developped" << finl;
2768 DoubleTab& val_scal_noeuds,
2772 switch(variables_internes_->methode_interpolation_v)
2778 champ_scal_interp=champ_scal;
2788 DoubleTab scal_filtre_val(champ_scal.
valeurs());
2791 ref_cast(
Champ_P1NC, champ_scal).filtrer_L2(scal_filtre_val);
2794 champ_scal_interp->valeurs() = scal_filtre_val;
2803 Cerr<<
"champ_scal type ="<<champ_scal.
que_suis_je()<<finl;
2804 Cerr <<
"Erreur dans Transport_Interfaces_FT_Disc::calculer_scalaire_interpole\n";
2805 Cerr <<
" L'interpolation a partir d'un champ " << champ_scal.
que_suis_je();
2806 Cerr <<
" n'est pas codee." << finl;
2808 Cerr<<
"champ_scal type ="<<champ_scal.
que_suis_je()<<finl;
2809 Cerr <<
"Error for the method Transport_Interfaces_FT_Disc::calculer_scalaire_interpole\n";
2810 Cerr <<
" The interpolation from a field " << champ_scal.
que_suis_je();
2811 Cerr <<
" is not developped." << finl;
2818 DoubleTab& les_positions = variables_internes_->doubletab_pos;
2819 IntVect& les_elements = variables_internes_->intvect_elements;
2822 const DoubleTab& pos = maillage.
sommets();
2824 const int nb_pos_tot = pos.
dimension(0);
2826 les_positions.
resize(nb_pos_tot, dim);
2827 les_elements.
resize(nb_pos_tot);
2829 int nb_positions = 0;
2830 for (i = 0; i < nb_pos_tot; i++)
2832 const int num_elem = elem[i];
2835 for (j = 0; j < dim; j++)
2836 les_positions(nb_positions, j) = pos(i, j);
2837 les_elements(nb_positions) = num_elem;
2841 les_positions.
resize(nb_positions, dim);
2842 les_elements.
resize(nb_positions);
2844 champ_scal_interp->valeur_aux_elems(les_positions, les_elements,val_scal_noeuds);
2857 Cerr <<
"Transport_Interfaces_FT_Disc::calculer_scalaire_interpole\n"
2858 <<
" the interpolation method is not coded" << finl;
2876 Cerr<<
" On ne doit pas resoudre Transport_Interfaces_FT_Disc en implicite "<<finl;
2880static void init_parser_v_impose(
const Noms& expression_vitesse,
Parser& parser_x,
Parser& parser_y,
Parser& parser_z,
double temps)
2884 std::string sx(expression_vitesse[0]);
2892 parser_x.
setVar(
"z", 0.);
2893 parser_x.
setVar(
"t", temps);
2895 std::string sy(expression_vitesse[1]);
2903 parser_y.
setVar(
"z", 0.);
2904 parser_y.
setVar(
"t", temps);
2906 Nom unused_expr(
"0");
2907 std::string sz(dimension3
2908 ? expression_vitesse[2]
2918 parser_z.
setVar(
"z", 0.);
2919 parser_z.
setVar(
"t", temps);
2941 DoubleTab& vpoint,
const DoubleTab& rho_faces,
2942 DoubleTab& source_val,
const double temps,
2943 const double dt,
const int is_explicite,
2946 if (!(temps>temps_debut_))
2952 DoubleTab vit_imposee;
2958 const IntTab& face_voisins = mon_dom_dis.
face_voisins();
2982 if(variables_internes_->vimp_regul && is_explicite)
2985 DoubleTab vitesse(vpoint0);
2987 vitesse += inco_val ;
2989 for (
int i=0 ; i < n; i++)
2991 if (indicatrice_faces(i) > 0. )
2993 double f = indicatrice_faces(i);
2994 for (
int j = 0; j < m; j++)
2997 vit_imposee(i,j) = f*vit_imposee(i,j) + (1.-f)*vitesse(i,j);
2999 vit_imposee(i,j) = f*vit_imposee(i,j) + (1.-f)*vitesse(i,j)/rho_faces(i);
3001 vitesse_imp_interp_->valeurs()(i,j)= vit_imposee( i,j ) ;
3008 calcul_source(inco_val,vpoint0,rho_faces,source_val,vit_imposee,indicatrice_faces,
3009 is_QC,dt,is_explicite,eta);
3012 const DoubleVect& volumes_entrelaces = ref_cast(
Domaine_VF,mon_dom_dis).volumes_entrelaces();
3016 DoubleTab termes_sources_face(vpoint);
3017 termes_sources_face=0.;
3018 modifie_source(termes_sources_face,source_val,rho_faces,n,m,is_QC,volumes_entrelaces,le_solveur_masse);
3021 for (i = 0; i < n; i++)
3022 for (j = 0; j < m; j++)
3023 vpoint(i, j) += termes_sources_face(i,j);
3031 Cerr <<
"Error for the method Transport_Interfaces_FT_Disc::modifier_vpoint_pour_imposer_vit\n"
3032 <<
" The transport equation is not of type \"methode_transport vitesse_imposee\"."
3033 <<
" or \"methode_transport loi_horaire\"."
3038 Debog::verifier(
"Transport_Interfaces_FT_Disc::calculer_vitesse_imposee vpoint",vpoint);
3043 const IntTab& face_voisins)
3045 DoubleTab& indicatrice_faces = indicatrice_faces_->valeurs();
3047 for (
int i = 0; i < nfaces; i++)
3049 const int elem0 = face_voisins(i, 0);
3050 const int elem1 = face_voisins(i, 1);
3051 indicatrice_faces(i)= 0.;
3054 indicatrice_faces(i) = indicatrice(elem0);
3056 indicatrice_faces(i) += indicatrice(elem1);
3057 if (elem0 >= 0 && elem1 >= 0)
3058 indicatrice_faces(i) *= 0.5;
3059 double bmax=1.-1e-9;
3060 if(indicatrice_faces(i) <= 1.0e-9) indicatrice_faces(i) = 0. ;
3061 else if(indicatrice_faces(i) >= bmax) indicatrice_faces(i) = 1. ;
3064 switch(variables_internes_->type_indic_faces_)
3077 const DoubleVect& face_surfaces = domaine_vdf.
face_surfaces();
3079 double& position = variables_internes_->modified_indic_faces_position;
3080 double& thickness = variables_internes_->modified_indic_faces_thickness;
3081 for (
int i = 0; i < nfaces; i++)
3083 double h=volumes_entrelaces(i)/face_surfaces(i);
3084 if (dist_face(i) > (position+thickness/2.)*h || indicatrice_faces(i)==1.)
3085 indicatrice_faces(i)=1.;
3086 else if (dist_face(i) >= (position-thickness/2.)*h && thickness !=0.)
3087 indicatrice_faces(i) = (dist_face(i) - position*h)/(thickness*h) + 0.5 ;
3089 indicatrice_faces(i)=0.;
3101 const Equation_base& eqn_hydraulique = variables_internes_->refequation_vitesse_transport.valeur();
3110 const int vef = (dim == 2);
3113 Cerr <<
"Code never applied or checked in VEF. You should read the algo first and assess it!" << finl;
3118 for (
int face = 0; face < nfaces; face++)
3120 double indic_face = 0.;
3122 for (v = 0; v < 2; v++)
3124 const int elem = face_voisins(face, v);
3128 const double indic = indicatrice[elem];
3130 if (indic <=5e-3 || indic >= 1.-5e-3)
3138 const double ai= interfacial_area(elem);
3139 if (fabs(ai)>DMINFLOAT)
3144 for (
int j = 0; j < dim; j++)
3147 const double nx = normale_elements(elem, j);
3154 Cerr <<
"Never tested. To be verified. It should depend on a scalar product with the vect (xp-xv)" << finl;
3162 const int dir = orientation[face];
3163 const double nx = normale_elements(elem, dir);
3180 tmp = indicatrice(elem);
3182 double bmax=1.-1e-9;
3183 if(tmp <= 1.0e-9) tmp = 0. ;
3184 else if(tmp >= bmax) tmp = 1. ;
3192 Cerr <<
" WTF, c'est impossible" << finl;
3200 const int elem_voisin = face_voisins(face, 1-v);
3201 const double indic = indicatrice[elem_voisin];
3217 indicatrice_faces_->valeurs()(face) = indic_face;
3222 Cerr <<
"Interpolation option AI_BASED in Transport_Interfaces_FT_Disc is not available "
3223 <<
"in Transport_Interfaces_FT_Disc::calcul_indicatrice_faces if we do not "
3224 <<
"have a Navier_Stokes_FT_Disc equation..." << finl;
3230 Cerr <<
"Transport_Interfaces_FT_Disc::calcul_indicatrice_faces\n"
3231 <<
" unknown case?" << finl;
3236 Debog::verifier(
"Transport_Interfaces_FT_Disc::calcul_indicatrice_faces indicatrice_faces",indicatrice_faces);
3241 return variables_internes_->vimp_regul;
3245 return indicatrice_faces_;
3252 const IntTab& face_voisins = mon_dom_dis.
face_voisins();
3254 return indicatrice_faces_;
3266 const DoubleTab& vpoint,
3267 const DoubleTab& rho_faces,
3268 DoubleTab& source_val,
const DoubleTab& vit_imposee,
const DoubleTab& indicatrice_faces,
3271 const int is_explicite,
3275 DoubleTab terme_explicite(vpoint);
3279 if (indicatrice_faces.
dimension(0) == nfaces )
3283 terme_explicite *= dt;
3287 terme_explicite *= 0.;
3290 for (
int i = 0; i < nfaces; i++)
3293 double indic = (indicatrice_faces(i) > 0. ? 1.0 : 0.0) ;
3296 if (variables_internes_->vimp_regul && !is_explicite ) indic = (indicatrice_faces(i) ==1. ? 1.0 : 0.0);
3298 double rho_face = rho_faces(i);
3301 for (
int j = 0; j < dim; j++)
3303 double increment_inco = inco_val(i,j) + terme_explicite(i,j);
3305 tsource = (vit_imposee(i,j) - increment_inco) / dt * indic * rho_face * c;
3307 tsource = (vit_imposee(i,j) * rho_face - increment_inco) / dt * indic * c;
3309 source_val(i,j) = tsource;
3316 Cerr <<
"Erreur dans Transport_Interfaces_FT_Disc::calcul_source\n" << finl;
3317 Cerr <<
"Dimension de indicatrice_face differente de dimension de vpoint.\n" << finl;
3331 file+=
"_Force_totale_sur_";
3332 else if( type==
"force_totale" )
3334 file+=
"_Friction_totale_sur_";
3335 else if( type==
"Friction" )
3337 file+=
"_Friction_conv_diff_sur_";
3338 else if( type==
"Pressure" )
3340 file+=
"_Friction_Pression_sur_";
3341 else if ( type ==
"moment")
3343 file+=
"_Moment_total_sur_";
3344 else if ( type==
"particles_trajectory" )
3346 file+=
"_particles_trajectory_";
3347 else if ( type==
"mean_rms_particles_velocity")
3349 file+=
"_mean_rms_particles_velocity_";
3350 else if (type==
"particles_data")
3352 file+=
"_particles_data_";
3353 else if (type==
"particle_hydrodynamic_forces")
3356 file+=
"_particle_hydrodynamic_forces_";
3358 else if( type==
"Stokes_theoretical_forces" )
3361 file+=
"_Stokes_theoretical_forces_";
3363 else if (type==
"particle_hydrodynamic_forces_bed")
3366 file+=
"_particle_hydrodynamic_forces_bed_";
3368 else if (type==
"fluid_to_particle_heat_transfer")
3371 file+=
"_fluid_to_particle_heat_transfer_";
3373 else if (type==
"fluid_to_particle_heat_transfer_bed")
3376 file+=
"_fluid_to_particle_heat_transfer_bed_";
3378 else if(type==
"Friction_diff")
3380 file+=
"_Friction_diff_sur_" ;
3383 Cerr <<
"The file " << type <<
" is not understood by Transport_Interfaces_FT_Disc::ouvrir_fichier. "
3387 const int rang=files.search(type);
3402 fic << (
Nom)
"# Printing " << (type==
"moment"?
"of the drag moment exerted":
"of the drag exerted");
3403 fic <<
" by the fluid on the interface " << equation.
le_nom();
3404 fic <<
" " << (type==
"moment"?
"[N.m]":
"[N]") << finl;
3421 if (nb_compo>1) ch+=
"Fx";
3422 if (nb_compo>=2) ch+=espace+
"Fy";
3423 if (nb_compo>=3) ch+=espace+
"Fz";
3430 fic <<
"#########################################" << finl;
3431 fic <<
"# Position - Velocity - Collision force #" << finl;
3432 fic <<
"#########################################" << finl;
3433 fic <<
"# Time [s]" << finl;
3434 fic <<
"# Position of the gravity center of the particle [m] (px py pz)" << finl;
3435 fic <<
"# Velocity of the gravity center of the particle [m/s] (vx vy vz)" << finl;
3436 fic <<
"# Collision force discretized on the particle volume [N] (fcx fcy fcz)" << finl;
3438 fic <<
"# Time" << espace <<
"px py pz" << espace <<
"vx vy vz" << espace <<
"fcx fcy fcz" << finl;
3444 fic <<
"#####################################################" << finl;
3445 fic <<
"# Average velocity - Average velocity squared - RMS #" << finl;
3446 fic <<
"#####################################################" << finl;
3447 fic <<
"# Time [s]" << finl;
3448 fic <<
"# Average velocity of purely solid cells. For each purely solid cell,"
3449 " the velocity at gravity center is computed as the average velocity of"
3450 " the opposing faces weighted by its volume. [m/s] (vx_av vy_av vz_av)" << finl;
3451 fic <<
"# Average velocity squared of purely solid cells. [m^2/s^2] (vx2_av vy2_av vz2_av)" << finl;
3452 fic <<
"# Once the average velocity and the average velocity squared is known,"
3453 " the RMS is computed as sqrt(abs(vi_av^2 - vi2_av)) with i in {x,y,z}."
3454 " [-] (rmsx rmsy rmsz)" << finl;
3456 fic <<
"# Time" << espace <<
"vx_av vy_av vz_av" << espace <<
"vx2_av vy2_av vz2_av"
3457 << espace <<
"rmsx rmsy rmsz" << finl;
3463 fic <<
"############################################################" << finl;
3464 fic <<
"# Position - Velocity - Collision force - Collision number #" << finl;
3465 fic <<
"############################################################" << finl;
3466 fic <<
"# Time [s] - Total number of collision " << finl;
3467 fic <<
"# Position of the gravity center of the particle [m] (px py pz)" << finl;
3468 fic <<
"# Velocity of the gravity center of the particle [m/s] (vx vy vz)" << finl;
3469 fic <<
"# Collision force discretized on the particle volume [N] (fcx fcy fcz)" << finl;
3471 fic <<
"# Time" <<espace <<
"total collision number" << finl;
3472 fic <<
"# particle_id px py pz vx vy vz fcx fcy fcz number_of_particles_in_collision"<< finl;
3477 fic <<
"###################################" << finl;
3478 fic <<
"# Hydrodynamic force computation #" << finl;
3479 fic <<
"###################################" << finl;
3481 fic <<
"# Time [s]"<< espace <<
"Particle surface [m^2]" << espace <<
3482 "Pressure force [N] (fpx fpy fpz)" << espace <<
"Friction force [N] (ffx ffy ffz)" << finl;
3487 fic <<
"#####################################################################" << finl;
3488 fic <<
"# Hydrodynamic force computation - Stokes theoretical configuration #" << finl;
3489 fic <<
"#####################################################################" << finl;
3490 fic <<
"# Time [s]"<< finl;
3491 fic <<
"# Stokes theoretical PRESSURE FORCE computed from the integration, on the lagrangian mesh,"
3492 " of the discretized analytical solution [N] (fpx_th fpy_th fpz_th)" << finl;
3493 fic <<
"# Stokes PRESSURE FORCE computed with the developed method on the theoretical"
3494 " pressure field discretized on the eulerian mesh [N] (fpx_th_interp fpy_th_interp"
3495 " fpz_th_interp)" << finl;
3496 fic <<
"# Stokes theoretical FRICTION FORCE computed from the integration,"
3497 " on the lagrangian mesh, of the discretized analytical solution [N]"
3498 " (ffx_th ffy_th ffz_th)" << finl;
3499 fic <<
"# Stokes FRICTION FORCE computed with the developed method on the theoretical"
3500 " velocity field discretized on the eulerian mesh [N] (ffx_th_interp"
3501 " ffy_th_interp ffz_th_interp)" << finl;
3503 fic <<
"# Time" << espace <<
"fpx_th fpy_th fpz_th" << espace <<
3504 "fpx_th_interp fpy_th_interp fpz_th_interp" << espace <<
"ffx_th ffy_th ffz_th" <<
3505 espace <<
"ffx_th_interp ffy_th_interp ffz_th_interp" << finl;
3510 fic <<
"#########################################################" << finl;
3511 fic <<
"# Hydrodynamic force computation in a particle assembly #" << finl;
3512 fic <<
"#########################################################" << finl;
3514 fic <<
"# Time [s]" << espace <<
"particle_id" << espace <<
"Pressure force [N] (fpx fpy fpz)"
3515 << espace <<
"Friction force [N] (ffx ffy ffz)" << espace <<
3516 "Percentage of facets for which forces were computable" << espace <<
3517 "Average fluid velocity in P2" <<
"Percentage of purely fluid cells in P2"<< finl;
3522 fic <<
"#########################" << finl;
3523 fic <<
"# Heat flux computation #" << finl;
3524 fic <<
"#########################" << finl;
3525 fic <<
"# Time [s]" << finl;
3526 fic <<
"# Computation of the heat flux received by the particle from the surrounding fluid. [W] (phi)" << finl;
3528 fic <<
"# Time" << espace <<
"phi" << finl;
3533 fic <<
"################################################" << finl;
3534 fic <<
"# Heat flux computation in a particle assembly #" << finl;
3535 fic <<
"################################################" << finl;
3536 fic <<
"# Time [s]" << finl;
3537 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;
3538 fic <<
"# Average temperature of purely fluid cells in P2. [K] (T_i)" << finl;
3540 fic <<
"Time" << espace <<
"phi_0 T_0 ... phi_N T_N" << finl;
3547 os.
ouvrir(file,ios::app);
3550 os.
setf(ios::scientific);
3553 const int n,
const int m,
const int is_QC,
3557 for (
int face=0; face<n; face++)
3558 for (
int dim=0; dim<m; dim++)
3559 termes_sources_face(face,dim)=vol_entrelaces(face)*source_val(face,dim);
3562 un_solv_masse.
appliquer(termes_sources_face);
3566 for (
int i = 0; i < n; i++)
3568 const double rho_face = rho_faces(i);
3570 for (
int j = 0; j < m; j++)
3571 termes_sources_face(i,j) = termes_sources_face(i,j) / rho_face;
3580 const int nbdim1 = source_val.
line_size() == 1;
3587 DoubleTab termes_sources_face(source_val);
3588 DoubleTab termes_pressure_face(pressure_part);
3589 DoubleTab termes_friction_face(friction_part);
3590 DoubleTab termes_diff_face(diff_part);
3594 const DoubleVect& vol_entrelaces = ref_cast(
Domaine_VF,mon_dom_dis).volumes_entrelaces();
3596 ArrOfInt sequential_items_flags;
3600 for (
int face=0; face<n; face++)
3602 double indic = (indicatrice_faces(face) > 0. ? 1.0 : 0.0);
3603 double coef = vol_entrelaces(face)*indic;
3604 for (
int dim=0; dim<m; dim++)
3606 termes_sources_face(face,dim)=source_val(face,dim)*coef;
3607 termes_pressure_face(face,dim)=pressure_part(face,dim)*coef;
3608 termes_friction_face(face,dim)=friction_part(face,dim)*coef;
3609 termes_diff_face(face,dim)=diff_part(face,dim)*coef;
3613 if (sequential_items_flags[face])
3619 values(0,j) -= termes_sources_face(face,0);
3620 values(1,j) -= termes_pressure_face(face,0);
3621 values(2,j) -= termes_friction_face(face,0);
3622 values(3,j) -= termes_diff_face(face,0);
3628 values(0,j) -= termes_sources_face(face,j);
3629 values(1,j) -= termes_pressure_face(face,j);
3630 values(2,j) -= termes_friction_face(face,j);
3631 values(3,j) -= termes_diff_face(face,j);
3646 ouvrir_fichier(Force,
"force_totale",1,*
this);
3652 Force << espace << values(0,k);
3654 const int impr_mom = 1 ;
3657 ouvrir_fichier(Pressure,
"Pressure",impr_mom,*
this);
3663 Pressure << espace << values(1,k);
3667 ouvrir_fichier(Friction,
"Friction",impr_mom,*
this);
3673 Friction << espace << values(2,k);
3677 ouvrir_fichier(Friction2,
"Friction_diff",impr_mom,*
this);
3682 Friction2 << espace << values(3,k);
3696 ouvrir_fichier(Force,
"force",1,*
this);
3700 Force << espace << force_[k];
3707 ouvrir_fichier(Moment,
"moment",impr_mom,*
this);
3709 for(
int k=0; k<moment_.size_array(); k++)
3710 Moment << espace << moment_[k];
3718 if (nb_particles_tot_<dim_max_impr)
3721 ouvrir_fichier(Particles_trajectory,
"particles_trajectory",1,*
this);
3723 for (
int particle=0; particle<nb_particles_tot_; particle++)
3725 Particles_trajectory << espace;
3728 Particles_trajectory << espace;
3732 Particles_trajectory << espace;
3734 if (collision_model_)
3736 const DoubleTab& lagrangian_contact_forces=
3737 collision_model_.valeur().get_lagrangian_contact_forces();
3739 Particles_trajectory << espace << lagrangian_contact_forces(particle,dim);
3741 Particles_trajectory << finl;
3747 SFichier Moy_Rms_Vitesse_Particule;
3748 ouvrir_fichier(Moy_Rms_Vitesse_Particule,
"mean_rms_particles_velocity",1,*
this);
3756 for (
int particle=0; particle<nb_particles_tot_; particle++)
3758 Moy_Rms_Vitesse_Particule << particle <<
" ";
3759 Moy_Rms_Vitesse_Particule << espace;
3760 Moy_Rms_Vitesse_Particule << particles_purely_solid_mesh_volume(particle)<< espace;
3763 Moy_Rms_Vitesse_Particule << espace << mean_velocity(particle,dim);
3764 Moy_Rms_Vitesse_Particule << espace;
3767 Moy_Rms_Vitesse_Particule << espace << mean_squared_velocity(particle,dim);
3768 Moy_Rms_Vitesse_Particule << espace;
3771 Moy_Rms_Vitesse_Particule << espace << rms_velocity(particle,dim);
3773 Moy_Rms_Vitesse_Particule << finl;
3778 ouvrir_fichier(Particles_data,
"particles_data",1,*
this);
3779 if (collision_model_)
3781 const int collision_number=collision_model_.valeur().get_collision_number();
3783 collision_number<< finl;
3785 for (
int particle = 0; particle < nb_particles_tot_; particle++)
3787 Particles_data << particle <<
" ";
3794 if (collision_model_)
3796 const DoubleTab& lagrangian_contact_forces=
3797 collision_model_.valeur().get_lagrangian_contact_forces();
3798 const ArrOfDouble particles_collision_number=collision_model_.valeur().get_particles_collision_number();
3801 Particles_data << lagrangian_contact_forces(particle,d) << espace ;
3803 Particles_data << particles_collision_number(particle) << espace ;
3805 Particles_data << finl;
3815 if (nb_particles_tot_<dim_max_impr)
3818 ouvrir_fichier(particle_hydro_forces,
"particle_hydrodynamic_forces",1,mon_eq);
3821 for (
int compo=0; compo<nb_particles_tot_; compo++)
3823 particle_hydro_forces << espace << total_surface_interf(compo) << espace;
3824 for (
int dim=0; dim<
dimension; dim++) particle_hydro_forces <<
3825 espace << total_pressure_force(compo,dim);
3826 particle_hydro_forces << espace;
3827 for (
int dim=0; dim<
dimension; dim++) particle_hydro_forces <<
3828 espace << total_friction_force(compo,dim);
3829 particle_hydro_forces << espace;
3832 particle_hydro_forces << finl;
3836 SFichier particle_hydro_forces_bed;
3837 ouvrir_fichier(particle_hydro_forces_bed,
"particle_hydrodynamic_forces_bed",1,mon_eq);
3838 particle_hydro_forces_bed <<
"TIME ";
3840 particle_hydro_forces_bed << finl;
3844 for (
int particle = 0; particle < nb_particles_tot_; particle++)
3847 particle_hydro_forces_bed << particle ;
3848 particle_hydro_forces_bed << espace << total_surface_interf(particle) << espace;
3849 for (dim=0; dim<
dimension; dim++) particle_hydro_forces_bed << espace <<
3850 total_pressure_force(particle,dim);
3851 particle_hydro_forces_bed << espace;
3852 for (dim=0; dim<
dimension; dim++) particle_hydro_forces_bed << espace <<
3853 total_friction_force(particle,dim);
3854 particle_hydro_forces_bed << espace;
3855 for (dim=0; dim<
dimension; dim++) particle_hydro_forces_bed << espace <<
3856 Nb_fa7_ok_prop(particle,dim);
3857 for (dim=0; dim<
dimension; dim++) particle_hydro_forces_bed << espace <<
3858 U_P2_moy(particle,dim);
3859 particle_hydro_forces_bed << espace << Prop_indic_fluide_P2(particle);
3860 particle_hydro_forces_bed << finl;
3871 get_pressure_force_Stokes_th();
3873 get_friction_force_Stokes_th();
3874 int nb_particles_tot = total_pressure_force.
dimension(0);
3876 if (nb_particles_tot<dim_max_impr)
3879 ouvrir_fichier(Force_Particule_th,
"Stokes_theoretical_forces",1,mon_eq);
3881 for (
int particle=0; particle<nb_particles_tot; particle++)
3883 Force_Particule_th << espace;
3884 for (
int dim=0; dim<
dimension; dim++) Force_Particule_th << espace <<
3885 total_pressure_force_Stokes(particle,dim);
3886 Force_Particule_th << espace;
3887 for (
int dim=0; dim<
dimension; dim++) Force_Particule_th << espace <<
3888 total_pressure_force(particle,dim);
3889 Force_Particule_th << espace;
3890 for (
int dim=0; dim<
dimension; dim++) Force_Particule_th << espace <<
3891 total_friction_force_Stokes(particle,dim);
3892 Force_Particule_th << espace;
3893 for (
int dim=0; dim<
dimension; dim++) Force_Particule_th << espace <<
3894 total_friction_force(particle,dim);
3895 Force_Particule_th << finl;
3903 const DoubleVect& total_heat_transfer =
3905 if (nb_particles_tot_ < dim_max_impr)
3908 ouvrir_fichier(Heat_Transfer,
3909 "fluid_to_particle_heat_transfer", 1, mon_eq);
3911 for (
int particle = 0; particle < nb_particles_tot_;
3914 Heat_Transfer << espace << total_heat_transfer(particle);
3916 Heat_Transfer << finl;
3921 ouvrir_fichier(Heat_Transfer_bed,
3922 "fluid_to_particle_heat_transfer_bed", 1, mon_eq);
3924 for (
int particle = 0; particle < nb_particles_tot_;
3927 Heat_Transfer_bed << espace << total_heat_transfer(particle);
3929 Heat_Transfer_bed << finl;
3946 DoubleTab tab_critere(present);
3948 tab_critere = present;
3949 tab_critere -=passe;
3955 const DoubleTab& rho_faces,DoubleTab& source_val,
3956 const int is_explicite,
const double eta)
3972 if ( !is_explicite )
3974 for (
int i = 0; i<n; i++)
3976 double indic = (indicatrice_faces(i) > 0. ? 1.0 : 0.0);
3977 for (
int j = 0; j < m; j++)
3980 source_val(i,j) -= vpoint(i,j) * indic * rho_faces(i) * c;
3982 source_val(i,j) -= vpoint(i,j) * indic * c;
3987 DoubleTab termes_sources_face(vpoint);
3988 const DoubleVect& vol_entrelaces = ref_cast(
Domaine_VF,mon_dom_dis).volumes_entrelaces();
3990 for (
int face=0; face<n; face++)
3991 for (
int dim=0; dim<m; dim++)
3993 double indic = (indicatrice_faces(face) > 0. ? 1.0 : 0.0);
3994 termes_sources_face(face,dim)=vol_entrelaces(face)*vol_entrelaces(face)*source_val(face,dim)*indic;
3998 le_solveur_masse.
appliquer(termes_sources_face);
4007 const ArrOfDouble& centre_gravite = domaine.cg_moments();
4016 ArrOfInt sequential_items_flags;
4020 for (
int i = 0; i < n; i++)
4024 if (sequential_items_flags[i])
4030 dforce[j] = -termes_sources_face(i);
4031 force_[j] += dforce[j];
4036 dforce[j] = -termes_sources_face(i,j);
4037 force_[j] += dforce[j];
4043 xgr[j] = centre_faces(i,j) - centre_gravite[j];
4046 moment_[0] += dforce[1]*xgr[0] - dforce[0]*xgr[1];
4049 moment_[0] += dforce[2]*xgr[1] - dforce[1]*xgr[2];
4050 moment_[1] += dforce[0]*xgr[2] - dforce[2]*xgr[0];
4051 moment_[2] += dforce[1]*xgr[0] - dforce[0]*xgr[1];
4068 const double temps=le_schema_en_temps->temps_courant();
4070 const int dimension3 = (dim==3);
4071 Parser parser_x, parser_y, parser_z;
4073 const IntTab& face_voisins = domaine_vf.
face_voisins();
4074 const int nfaces = face_voisins.
dimension(0);
4075 const DoubleTab& xv = domaine_vf.
xv();
4081 vit_ibc.
resize(nfaces,1);
4083 vit_ibc.
resize(nfaces,dim);
4085 init_parser_v_impose(variables_internes_->expression_vitesse_imposee,
4086 parser_x, parser_y, parser_z,temps);
4087 for (
int i = 0; i < nfaces; i++)
4089 for (
int j = 0; j < dim; j++)
4091 const double coord = xv(i,j);
4092 parser_x.
setVar(j, coord);
4093 parser_y.
setVar(j, coord);
4095 parser_z.
setVar(j, coord);
4103 vit_ibc(i,0) = parser_x.
eval();
4106 vit_ibc(i,0) = parser_y.
eval();
4109 vit_ibc(i,0) = parser_z.
eval();
4115 vit_ibc(i,0) = parser_x.
eval();
4116 vit_ibc(i,1) = parser_y.
eval();
4118 vit_ibc(i,2) = parser_z.
eval();
4128 const DoubleTab& vitesse,
4129 const DoubleTab& vpoint,
4134 const int dimension3 = (dim==3);
4136 const IntTab& face_voisins = mon_dom_dis.
face_voisins();
4137 const int nfaces = face_voisins.
dimension(0);
4138 const DoubleTab& xv = ref_cast(
Domaine_VF, mon_dom_dis).xv();
4144 vitesse_imp.
resize(nfaces,1);
4146 vitesse_imp.
resize(nfaces,dim);
4150 ArrOfDouble coord(dim);
4151 ArrOfDouble v_imp(dim);
4152 for (
int i = 0; i < nfaces; i++)
4155 for (
int j = 0; j < dim; j++)
4157 v_imp = variables_internes_->loi_horaire_->vitesse(temps+dt,coord);
4162 for (
int j = 0; j < dim; j++)
4163 vitesse_imp(i,j) = v_imp[j];
4168 switch(variables_internes_->interpolation_champ_face)
4172 Parser parser_x, parser_y, parser_z;
4173 init_parser_v_impose(variables_internes_->expression_vitesse_imposee,
4174 parser_x, parser_y, parser_z,temps);
4175 for (
int i = 0; i < nfaces; i++)
4178 for (
int j = 0; j < dim; j++)
4181 const double coord = xv(i,j);
4182 parser_x.
setVar(j, coord);
4183 parser_y.
setVar(j, coord);
4185 parser_z.
setVar(j, coord);
4193 vitesse_imp(i,0) = parser_x.
eval();
4196 vitesse_imp(i,0) = parser_y.
eval();
4199 vitesse_imp(i,0) = parser_z.
eval();
4207 vitesse_imp(i,0) = parser_x.
eval();
4208 vitesse_imp(i,1) = parser_y.
eval();
4210 vitesse_imp(i,2) = parser_z.
eval();
4218 if ( !variables_internes_->vf_explicite )
4225 dt_loc = (1./dt_loc)*(le_schema_en_temps->facteur_securite_pas());
4228 dt_loc = std::min(dt_loc,dt);
4229 Cerr <<
"Transport_Interfaces_FT_Disc::calculer_vitesse_imposee avec dt_loc : "<<dt_loc<<finl;
4234 const int phase = 0;
4235 vitesse_imp = vpoint;
4236 vitesse_imp*= dt_loc;
4237 vitesse_imp += vitesse;
4239 DoubleTab gradient(vitesse);
4240 interpoler_vitesse_face(dist_face,phase,variables_internes_->n_iterations_interpolation_ibc,vitesse_imp,gradient,temps,dt);
4241 Debog::verifier(
"Transport_Interfaces_FT_Disc::calcul_vitesse vitesse_imp apres interp. ",vitesse_imp);
4245 Cerr <<
"Transport_Interfaces_FT_Disc::calculer_vitesse_imposee\n"
4246 <<
" interpolation lineaire non codee en VEF" << finl;
4253 Cerr <<
"Transport_Interfaces_FT_Disc::calculer_vitesse_imposee\n"
4254 <<
" methode d'interpolation non codee" << finl;
4261 Cerr <<
"Error for the method Transport_Interfaces_FT_Disc::calcul_vitesse\n"
4262 <<
" The transport equation is not of type \"methode_transport vitesse_imposee\"."
4263 <<
" or \"methode_transport loi_horaire\"."
4282 const DoubleTab& dist_elem,
4283 const DoubleTab& normale_elem,
4284 DoubleTab& dist_face)
const
4286 statistics().create_custom_counter(
"Calculer_distance_interface",3,
"FrontTracking");
4287 statistics().begin_count(
"Calculer_distance_interface",statistics().get_last_opened_counter_level()+1);
4289 static const double distance_faces_invalides = -1.e30;
4292 const IntTab& elem_faces = domaine_vf.
elem_faces();
4293 const DoubleTab& xv = domaine_vf.
xv();
4294 const DoubleTab& xp = domaine_vf.
xp();
4296 ArrOfInt ncontrib(nb_faces);
4301 double centre[3] = {0., 0., 0.};
4302 double normale[3] = {0., 0., 0.};
4306 const int nb_faces_elem = elem_faces.
line_size();
4308 for (elem = 0; elem < nb_elem_tot; elem++)
4310 const double d1 = dist_elem(elem);
4312 if (d1 < distance_faces_invalides)
4315 for (i = 0; i < dim; i++)
4317 centre[i] = xp(elem, i);
4318 normale[i] = normale_elem(elem, i);
4321 for (i = 0; i < nb_faces_elem; i++)
4323 const int face = elem_faces(elem, i);
4326 if (face < nb_faces)
4330 for (j = 0; j < dim; j++)
4332 double position_centre_face = xv(face, j);
4333 d2 += (position_centre_face - centre[j]) * normale[j];
4335 dist_face(face) += d1 + d2;
4341 const double valeur_invalide = distance_faces_invalides * 1.1;
4342#ifdef __INTEL_COMPILER
4345 for (i = 0; i < nb_faces; i++)
4347 const int n = ncontrib[i];
4351 dist_face(i) = valeur_invalide;
4354 statistics().end_count(
"Calculer_distance_interface");
4362 if (tag == variables_internes_->distance_faces_cache_tag)
4363 return variables_internes_->distance_interface_faces.valeur();
4365 variables_internes_->distance_faces_cache_tag = tag;
4366 assert(tag == variables_internes_->distance_normale_cache_tag);
4369 DoubleTab& dist_face = variables_internes_->distance_interface_faces->valeurs();
4372 return variables_internes_->distance_interface_faces.valeur();
4376inline void check(
const DoubleTab& v_imp,
int& i_face,
double v,
int v_est_initialise)
4378 if (v_est_initialise && v_imp(i_face)!= v)
4380 Cerr <<
"=====================================================================================" << finl;
4381 Cerr <<
"You defined a non-uniform function for the keyword: methode_transport vitesse_imposee" << finl;
4382 Cerr <<
"Please add:" << finl;
4383 Cerr <<
"type_vitesse_imposee analytique" << finl;
4389 v_est_initialise = 1;
4410 const DoubleTab& distance_interface_faces,
4412 const int stencil_width,
4414 DoubleTab& gradient,
4415 const double t,
const double dt)
4424 const IntTab& elem_faces = domaine_vf.
elem_faces();
4426 const int nfaces = faces_elem.
dimension(0);
4427 const int nb_elem = domaine_vf.
nb_elem() ;
4428 const int nb_elem_tot = domaine_vf.
nb_elem_tot() ;
4429 const IntVect& orientation = domaine_vdf.
orientation();
4430 const DoubleTab& xv = domaine_vf.
xv();
4431 const DoubleTab& xp = domaine_vf.
xp();
4432 assert(champ.dimension(0) == nfaces);
4433 assert(gradient.dimension(0) == nfaces);
4434 assert(distance_interface_faces.
dimension(0) == nfaces);
4439 Parser parser_x, parser_y, parser_z;
4440 init_parser_v_impose(variables_internes_->expression_vitesse_imposee,parser_x, parser_y, parser_z, t);
4442 const double invalid_test = -1.e30;
4443 const double invalid_value = -2.e30;
4444 gradient = invalid_value;
4445 const double indic_phase = (phase == 0) ? 0. : 1.;
4451 DoubleTab dist_face_cor(distance_interface_faces) ;
4453 DoubleTab typeface(champ) ;
4464 IntTab trav(nb_elem,dim) ;
4465 if( nb_elem < nb_elem_tot )
4471 for(
int i_elem=0 ; i_elem<nb_elem ; i_elem++)
4473 for(
int i=0 ; i<dim ; i++ )
4474 xe(i)=xp(i_elem,i) ;
4476 if( indicatrice(i_elem) != 0. && indicatrice(i_elem) != 1. )
4478 for(
int dir = 0 ; dir < dim ; dir++ )
4480 int a = elem_faces(i_elem, dir) ;
4481 int b = elem_faces(i_elem, dir+dim) ;
4482 const double pas = std::fabs(domaine_vdf.
distance_face(a,b,dir)) ;
4485 trav(i_elem,dir) = traverse ;
4488 else if( indicatrice(i_elem) == 0. )
4491 for(
int k=0; k<dim; k++)
4492 trav(i_elem,k) = -1 ;
4494 else if( indicatrice(i_elem) == 1. )
4497 for(
int k=0; k<dim; k++)
4498 trav(i_elem,k) = -2 ;
4510 int cpt_halfref = 0 ;
4511 DoubleTab typefacejoint ;
4515 for(
int njoint=0; njoint<nbjoints; njoint++)
4517 const Joint& joint_temp = domaine_vf.
joint(njoint);
4519 const int nb_faces = indices_faces_joint.
dimension(0);
4520 for (
int j = 0; j < nb_faces; j++)
4522 int face_de_joint = indices_faces_joint(j, 1);
4528 for( i_face = 0 ; i_face < nfaces ; i_face++ )
4530 if( indicatrice_face(i_face) == 0. || indicatrice_face(i_face) == 1. )
4532 typeface(i_face) = indicatrice_face(i_face) ;
4536 typefacejoint = typeface ;
4538 for(
int i=0 ; i<FACEJOINT.
size() ; i++ )
4540 typeface(FACEJOINT[i]) = std::max(typefacejoint(FACEJOINT[i]),typeface(FACEJOINT[i])) ;
4541 if( typefacejoint(FACEJOINT[i]) < 0 && typeface(FACEJOINT[i]) > -1 )
4546 for( i_face = 0 ; i_face < nfaces ; i_face++ )
4549 if( indicatrice_face(i_face) != 0. && indicatrice_face(i_face) != 1. && indicatrice_face(i_face) != 0.5 )
4552 int reference_trouvee = 0 ;
4554 while( i_voisin < 2 && reference_trouvee == 0)
4556 const int elem_voisin = faces_elem(i_face,i_voisin) ;
4558 if( elem_voisin > -1 &&
4559 ( indicatrice(elem_voisin) == 0. || indicatrice(elem_voisin) == 1. ) )
4561 reference_trouvee = 1 ;
4562 typeface(i_face) = indicatrice(elem_voisin) ;
4569 typefacejoint = typeface ;
4571 for(
int i=0 ; i<FACEJOINT.
size() ; i++ )
4573 typeface(FACEJOINT[i]) = std::max(typefacejoint(FACEJOINT[i]),typeface(FACEJOINT[i])) ;
4574 if( typefacejoint(FACEJOINT[i]) < 0 && typeface(FACEJOINT[i]) > -1
4575 && indicatrice_face(FACEJOINT[i]) != 0. && indicatrice_face(FACEJOINT[i]) != 1.
4576 && indicatrice_face(FACEJOINT[i]) != 0.5 )
4581 for( i_face = 0 ; i_face < nfaces ; i_face++ )
4583 if( indicatrice_face(i_face) == 0.5 )
4586 if( faces_elem(i_face,0) > -1 )
4587 element = faces_elem(i_face,0) ;
4589 element = faces_elem(i_face,1) ;
4591 if( indicatrice(element) == 0. || indicatrice(element) == 1. )
4593 typeface(i_face) = 0.5 ;
4598 typefacejoint = typeface ;
4600 for(
int i=0 ; i<FACEJOINT.
size() ; i++ )
4602 typeface(FACEJOINT[i]) = std::max(typefacejoint(FACEJOINT[i]),typeface(FACEJOINT[i])) ;
4603 if( typefacejoint(FACEJOINT[i]) < 0 && typeface(FACEJOINT[i]) > -1 && indicatrice_face(FACEJOINT[i]) == 0.5 )
4622 int cpt_rest = nfaces - cpt_ref - cpt_vref - cpt_halfref ;
4626 double proc_stop = -1. ;
4633 if( proc_stop < 0. )
4636 for( i_face = 0 ; i_face<nfaces ; i_face++)
4638 if( typeface( i_face ) < 0 )
4640 const int ori = orientation[i_face] ;
4641 const int elem0 = faces_elem(i_face,0) ;
4642 const int elem1 = faces_elem(i_face,1) ;
4645 if( elem0 > -1 && elem0 < nb_elem_tot )
4647 else if( elem1 > -1 && elem1 < nb_elem_tot )
4654 int nb = trav(elem,ori) ;
4657 int j_face = elem_faces(elem, ori) + elem_faces(elem, ori+dim) - i_face ;
4659 if( typeface(j_face) > -1 )
4662 if( fmod( nb, 2. ) > 0 )
4667 typeface(i_face) = 1. - typeface(j_face) ;
4674 typeface(i_face) = typeface(j_face) ;
4680 elem = elem0 + elem1 - elem ;
4681 if( elem > -1 && elem < nb_elem_tot )
4683 nb = trav(elem,ori) ;
4684 j_face = elem_faces(elem, ori) + elem_faces(elem, ori+dim) - i_face ;
4686 if( typeface(j_face) > -1 )
4689 if( fmod( nb, 2. ) > 0 )
4694 typeface(i_face) = 1. - typeface(j_face) ;
4701 typeface(i_face) = typeface(j_face) ;
4710 typefacejoint = typeface ;
4712 for(
int i=0 ; i<FACEJOINT.
size() ; i++ )
4714 typeface(FACEJOINT[i]) = std::max(typefacejoint(FACEJOINT[i]),typeface(FACEJOINT[i])) ;
4715 if( typefacejoint(FACEJOINT[i]) < 0 && typeface(FACEJOINT[i]) > -1
4716 && indicatrice_face(FACEJOINT[i]) != 0. && indicatrice_face(FACEJOINT[i]) != 1.
4717 && indicatrice_face(FACEJOINT[i]) != 0.5 )
4722 if( cpt_diph == cpt_rest )
4726 stop =
mp_sum( proc_stop ) ;
4730 for( i_face = 0 ; i_face < nfaces ; i_face++ )
4732 if( indicatrice_face(i_face) != 0. && indicatrice_face(i_face) != 1. )
4734 int elem_voisin = 0 ;
4739 const int ori = orientation[i_face] ;
4740 if( faces_elem(i_face,0) > -1 )
4742 j_face = elem_faces(faces_elem(i_face,0), ori) + elem_faces(faces_elem(i_face,0), ori+dim) - i_face ;
4743 dxa = std::fabs(domaine_vdf.
distance_face(i_face,j_face,ori)) ;
4745 if( faces_elem(i_face,1) > -1 )
4747 j_face = elem_faces(faces_elem(i_face,1), ori) + elem_faces(faces_elem(i_face,1), ori+dim) - i_face ;
4748 dxb = std::fabs(domaine_vdf.
distance_face(i_face,j_face,ori)) ;
4752 elem_voisin = faces_elem(i_face,1) ;
4757 elem_voisin = faces_elem(i_face,0) ;
4761 for(
int k=0 ; k<dim ; k++ )
4765 int a = elem_faces(elem_voisin, k) ;
4766 int b = elem_faces(elem_voisin, k+dim) ;
4768 Lref += 0.25*pas*pas ;
4771 const double ratio = std::fabs(distance_interface_faces(i_face))/sqrt(Lref) ;
4772 const double tol_distface = 1e-3 ;
4773 const int test1 = ( ratio < tol_distface ) ;
4775 typeface(i_face) = 0.5 ;
4778 typefacejoint = typeface ;
4780 for(
int i=0 ; i<FACEJOINT.
size() ; i++ )
4782 typeface(FACEJOINT[i]) = typefacejoint(FACEJOINT[i]) ;
4795 for( i_face=0 ; i_face<nfaces ; i_face++ )
4797 if( distance_interface_faces(i_face) > invalid_test )
4798 dist_face_cor(i_face) = (2.*typeface(i_face) - 1.)*std::fabs(distance_interface_faces(i_face)) ;
4800 dist_face_cor(i_face) = distance_interface_faces(i_face) ;
4805 for( i_face=0 ; i_face<nfaces ; i_face++ )
4807 dist_face_cor(i_face) = distance_interface_faces(i_face) ;
4814 for( i_face=0 ; i_face<nfaces ; i_face++ )
4816 variables_internes_->distance_interface_faces_corrigee->valeurs()(i_face) = dist_face_cor(i_face) ;
4817 variables_internes_->distance_interface_faces_difference->valeurs()(i_face) = dist_face_cor(i_face) - distance_interface_faces(i_face) ;
4823 int cpt_face_fluide_signefaux = 0 ;
4824 int cpt_face_solide_signefaux = 0 ;
4825 int cpt_face_diphasique_signefaux = 0 ;
4826 int cpt_face_diphasique_signemodifiee = 0 ;
4827 int cpt_face_diphasique_signemodifiee_dnulle = 0 ;
4828 int nfaces_diphasique = 0 ;
4829 int nfaces_fluide = 0 ;
4830 int nfaces_solide = 0 ;
4832 for( i_face=0 ; i_face<nfaces ; i_face++ )
4834 if( indicatrice_face( i_face ) == 0. )
4836 else if( indicatrice_face( i_face ) == 1. )
4839 nfaces_diphasique++;
4841 if( indicatrice_face( i_face ) == 0. && dist_face_cor( i_face ) >= 0. )
4843 cpt_face_fluide_signefaux++ ;
4845 else if( indicatrice_face( i_face ) == 1. && dist_face_cor( i_face ) <= 0. && dist_face_cor( i_face ) > invalid_test )
4847 cpt_face_solide_signefaux++ ;
4849 else if( indicatrice_face( i_face ) != 1. && indicatrice_face( i_face ) != 0. )
4852 if( faces_elem(i_face,0) > -1 )
4853 element = faces_elem(i_face,0) ;
4855 element = faces_elem(i_face,1) ;
4857 if( indicatrice( element ) == 0. && dist_face_cor( i_face ) >= 0. )
4859 cpt_face_diphasique_signefaux++ ;
4861 else if( indicatrice( element ) == 1. && dist_face_cor( i_face ) <= 0. && dist_face_cor( i_face ) > invalid_test )
4863 cpt_face_diphasique_signefaux++ ;
4867 if( std::fabs(dist_face_cor(i_face) - distance_interface_faces(i_face)) != 0. )
4869 cpt_face_diphasique_signemodifiee++ ;
4871 if( dist_face_cor( i_face ) == 0. )
4873 cpt_face_diphasique_signemodifiee_dnulle++ ;
4878 Cerr <<
" -----------------------------------------------------"<< finl ;
4879 Cerr <<
" -----------------------------------------------------"<< finl ;
4880 Cerr <<
" -----------------------------------------------------"<< finl ;
4881 Cerr <<
" Nombre de faces : "<< nfaces <<finl ;
4882 Cerr <<
" Nombre de faces fluides : "<<nfaces_fluide<<finl;
4883 Cerr <<
" Nombre de faces fluides avec mauvais signe : "<<cpt_face_fluide_signefaux<<finl;
4884 Cerr <<
" Nombre de faces solides : "<<nfaces_solide<<finl;
4885 Cerr <<
" Nombre de faces solides avec mauvais signe : "<<cpt_face_solide_signefaux<<finl;
4886 Cerr <<
" Nombre de faces diphasiques : "<<nfaces_diphasique<<finl;
4887 Cerr <<
" Nombre de faces diphasiques avec mauvais signe : "<<cpt_face_diphasique_signefaux<<finl;
4888 Cerr <<
" Nombre de faces diphasiques avec signe modifie : "<<cpt_face_diphasique_signemodifiee<<finl;
4889 Cerr <<
" Nombre de faces diphasiques avec signe modifie et distance nulle : "<<cpt_face_diphasique_signemodifiee_dnulle<<finl;
4890 Cerr <<
" -----------------------------------------------------"<< finl ;
4891 Cerr <<
" -----------------------------------------------------"<< finl ;
4892 Cerr <<
" -----------------------------------------------------"<< finl ;
4897 DoubleTab v_imp(champ) ;
4899 Cerr <<
" INITIALISATION DE Vimp a Vsolide pour les faces solides "<<finl ;
4900 for ( i_face = 0; i_face < nfaces; i_face++)
4902 const double d = dist_face_cor(i_face) ;
4903 if ( d>0. || (indicatrice_face(i_face) == 1. && d < -1.e20) )
4905 for (
int j = 0; j < dim; j++)
4907 const double coord = xv(i_face,j);
4908 parser_x.
setVar(j, coord) ;
4909 parser_y.
setVar(j, coord) ;
4911 parser_z.
setVar(j, coord) ;
4918 v_imp(i_face) = parser_x.
eval() ;
4921 v_imp(i_face) = parser_y.
eval() ;
4924 v_imp(i_face) = parser_z.
eval() ;
4935 Cerr <<
"Transport_Interfaces_FT_Disc_interne::UNIFORME " << finl;
4936 double vx=0,vy=0,vz=0;
4938 for ( i_face = 0; i_face < nfaces; i_face++)
4943 double d = dist_face_cor(i_face) ;
4944 if (indicatrice_face(i_face) == indic_phase && d > invalid_test)
4947 for (
int j = 0; j < dim; j++)
4949 const double coord = xv(i_face,j);
4950 parser_x.
setVar(j, coord);
4951 parser_y.
setVar(j, coord);
4953 parser_z.
setVar(j, coord);
4963 v_imp(i_face) = parser_x.
eval() ;
4964 check(v_imp,i_face,vx,ix);
4967 v_imp(i_face) = parser_y.
eval() ;
4968 check(v_imp,i_face,vy,iy);
4971 v_imp(i_face) = parser_z.
eval() ;
4972 check(v_imp,i_face,vz,iz);
4975 gradient(i_face) = (champ(i_face) - v_imp(i_face)) / d ;
4985 Cerr <<
"Transport_Interfaces_FT_Disc_interne::ANALYTIQUE "<< finl ;
4989 double indic_pos=variables_internes_->modified_indic_faces_position;
4990 double indic_thi=variables_internes_->modified_indic_faces_thickness;
4991 if (indic_pos-0.5*indic_thi < -1.)
4993 Cerr <<
"-------------------------------------------------------------------------" << finl;
4994 Cerr <<
"-------------------------------------------------------------------------" << finl;
4995 Cerr <<
"The datafile ask to use type_vitesse_imposee=ANALYTIQUE and" << finl;
4996 Cerr <<
"type_indic_faces modifiee with a too large transition to reach zero" << finl;
4997 Cerr <<
"There is a big risk that fluid gradients cannot be computed" << finl;
4998 Cerr <<
"because fluid points (Indic_faces=0) are too far from the interface" << finl;
4999 Cerr <<
"To fix that, just do one of the following things:" << finl;
5000 Cerr <<
"If the imposed velocity field is uniform, use type_vitesse_imposee=UNIFORM" << finl;
5001 Cerr <<
"If not, use \"distance_projete_faces simplifiee\" " << finl;
5002 Cerr <<
"-------------------------------------------------------------------------" << finl;
5003 Cerr <<
"-------------------------------------------------------------------------" << finl;
5011 Cerr <<
"Transport_Interfaces_FT_Disc_interne::PROJETE_SIMPLIFIE, calcul a partir des distances et normales" << finl;
5020 for (i_face = 0; i_face < nfaces; i_face++)
5022 double d = dist_face_cor(i_face) ;
5023 const int diphasique_tmp = (indicatrice_face(i_face) != 0. &&
5024 (indicatrice_face(i_face) != 1. || (indicatrice_face(i_face) == 1. && d<=0. && d> invalid_test)));
5025 const int diphasique = ( variables_internes_->is_extra_diphasique_solide ? diphasique_tmp : (diphasique_tmp && d <= 0.) ) ;
5026 const int fluide = ( indicatrice_face(i_face) == 0. ) ;
5027 const int solide = ( variables_internes_->is_extra_solide && (d>0. || (indicatrice_face(i_face) == 1. && d<-1.e20)) ) ;
5028 const int monophasique = ( fluide || solide ) ;
5030 if (diphasique || (monophasique && d> invalid_test ))
5032 DoubleTab coord_projete(dim) ;
5033 DoubleTab normale_face(dim) ;
5034 DoubleTab coord_face(dim) ;
5035 coord_projete = 0. ;
5038 double cpt_normale = 0.;
5040 for (
int i=0 ; i<dim ; i++)
5041 coord_face(i) = xv(i_face, i) ;
5044 for(
int i=0 ; i<2; i++)
5046 const int elem_voisin = faces_elem(i_face,i) ;
5047 if (elem_voisin > -1)
5049 for (
int j = 0; j < dim; j++)
5050 normale_face(j) += normale_elem(elem_voisin, j);
5054 assert(cpt_normale!=0.);
5055 normale_face /=cpt_normale;
5056 double norme_normale = 0.;
5057 for (
int i=0 ; i<dim ; i++)
5058 norme_normale +=normale_face(i)*normale_face(i);
5059 assert(norme_normale!=0.);
5060 normale_face /= sqrt(norme_normale);
5062 for (
int i=0 ; i<dim ; i++)
5063 coord_projete(i) = coord_face(i) - d*normale_face(i);
5065 for (
int j = 0; j < dim; j++)
5067 parser_x.
setVar(j, coord_projete(j)) ;
5068 parser_y.
setVar(j, coord_projete(j)) ;
5070 parser_z.
setVar(j, coord_projete(j)) ;
5077 v_imp(i_face) = parser_x.
eval() ;
5080 v_imp(i_face) = parser_y.
eval() ;
5083 v_imp(i_face) = parser_z.
eval() ;
5105 const int exec_planfa7existant = 0 ;
5121 int nb_fa7_accepted = variables_internes_->nomb_fa7_accepted;
5122 int max_fa7 = nb_fa7_accepted ;
5127 IntList OutElem, OutFa7 ;
5130 calcul_OutElemFa7(maillage,indicatrice,nb_elem,nb_fa7_accepted,OutElem,OutFa7) ;
5136 calcul_maxfa7(maillage,indicatrice,nb_elem,max_fa7,exec_planfa7existant) ;
5148 DoubleTab Tab100( nb_elem, nb_facettes_in_elem ) ;
5149 DoubleTab Tab101( nb_elem, nb_facettes_in_elem ) ;
5150 DoubleTab Tab102( nb_elem, nb_facettes_in_elem ) ;
5151 DoubleTab Tab103( nb_elem, nb_facettes_in_elem ) ;
5152 DoubleTab Tab110( nb_elem, nb_facettes_in_elem ) ;
5153 DoubleTab Tab111( nb_elem, nb_facettes_in_elem ) ;
5154 DoubleTab Tab112( nb_elem, nb_facettes_in_elem ) ;
5155 IntTab Tab12( nb_elem, nb_facettes_in_elem ) ;
5157 IntTab CptFacette( nb_elem ) ;
5161 DoubleTab Vertex(nfaces,dim) ;
5162 DoubleTab PPP(nfaces,dim) ;
5165 if( nb_elem < nb_elem_tot )
5192 DoubleTab Barycentre(nb_facettes, 3) ;
5195 Cerr <<
"STOKAGE DES FACETTES " << finl ;
5198 StockageFa7(maillage,CptFacette,Tab100,Tab101,Tab102,Tab103,
5199 Tab110,Tab111,Tab112,Tab12,Barycentre,indicatrice,
5200 OutElem,fa7,exec_planfa7existant) ;
5206 if( OutFa7.
size() != OutElem.
size()*nb_fa7_accepted )
5208 Cerr <<
"ERREUR DE DIMENSIONNEMENT " << finl ;
5209 Cerr <<
" OutFa7.size() "<<OutFa7.
size()<<finl ;
5210 Cerr <<
" OutElem.size()*nb_fa7_accepted "<<OutElem.
size()*nb_fa7_accepted<<finl ;
5216 IntTab TabOutFa7(OutElem.
size(),nb_fa7_accepted) ;
5217 for(
int i=0 ; i< OutElem.
size() ; i++)
5220 CptFacette( OutElem[i] ) = nb_fa7_accepted ;
5221 for(
int j=0 ; j< nb_fa7_accepted ; j++)
5223 TabOutFa7(i,j) = OutFa7[j+i*nb_fa7_accepted] ;
5226 StockageFa7(maillage,Tab100,Tab101,Tab102,Tab103,Tab110,Tab111,Tab112,
5227 Tab12,Barycentre,OutElem,TabOutFa7,fa7) ;
5240 if( nb_elem < nb_elem_tot )
5242 Cerr <<
" RENUMEROTATION DES FACETTES VIRTUELLES " << finl ;
5243 RenumFa7(Barycentre,Tab110,Tab111,Tab112,Tab12,CptFacette,
5244 nb_facettes,nb_facettes_dim) ;
5263 Cerr <<
" IDENTIFICATICATION DU VERTEX LE PLUS PROCHE POUR LES FACES CONCERNEES PAR LA PROJECTION " << finl ;
5264 Cerr <<
" CAS DES FACES APPARTENANT A AU MOINS UN ELEMENT TRAVERSE : ETAPE 1 " <<finl ;
5267 Cerr <<
" CAS DES FACES APPARTENANT A AU MOINS UN ELEMENT TRAVERSE : ETAPE 2 " <<finl ;
5270 Cerr <<
" CAS DES FACES MONOPHASIQUE : ETAPE 3 " <<finl ;
5273 Cerr <<
" FIN IDENTIFICATION " << finl ;
5276 int nb_proj_modif = 0 ;
5277 int nb_proj_modif_tot = 0 ;
5278 Cerr <<
"CALCUL DU PROJETE POUR LES FACES QUI APPARTIENNENT " << finl ;
5279 Cerr <<
"A UN ELEMENT TRAVERSE PAR L INTERFACE" << finl;
5281 Tab101,Tab102,Tab103,Tab12,CptFacette,v_imp,PPP,parser_x,parser_y,parser_z) ;
5284 Cerr <<
" Nombre de projete modifie pour les faces diphasiques : "<<nb_proj_modif_tot<<
" au temps "<<
schema_temps().
temps_courant()<<finl ;
5288 Cerr <<
"CALCUL DU PROJETE POUR LES FACES FLUIDE " << finl ;
5290 t,dt,Tab100,Tab101,Tab102,Tab103,Tab12,CptFacette,v_imp,PPP,parser_x,parser_y,parser_z) ;
5294 Cerr <<
" Nombre de projete modifie pour les faces fluides : "<<nb_proj_modif_tot<<
" au temps "<<
schema_temps().
temps_courant()<<finl ;
5295 Debog::verifier(
"Transport_Interfaces_FT_Disc::interpoler_vitesse_face vimp",v_imp);
5299 for (i_face = 0; i_face < nfaces; i_face++)
5301 double d = dist_face_cor(i_face) ;
5302 if (indicatrice_face(i_face) == indic_phase && d > invalid_test)
5304 gradient(i_face) = (champ(i_face) - v_imp(i_face))/d ;
5307 gradient.echange_espace_virtuel();
5308 Debog::verifier(
"Transport_Interfaces_FT_Disc::interpoler_vitesse_face gradient avant etalement",gradient);
5314 DoubleTab gradient_old;
5315 Cerr <<
"Nombre d'iterations de la methode interpoler_vitesse_face = " << stencil_width << finl;
5316 for (
int iteration = 0; iteration < stencil_width; iteration++)
5320 gradient_old = gradient;
5324 for (i_face = 0; i_face < nfaces; i_face++)
5326 const double d = dist_face_cor(i_face) ;
5327 if (indicatrice_face(i_face) != indic_phase && d > invalid_test)
5333 const int ori=orientation(i_face);
5334 const int elem0 = faces_elem(i_face, 0);
5335 const int elem1 = faces_elem(i_face, 1);
5342 for (
int direction=0; direction <
dimension; direction++)
5349 const int voisin0 = elem_faces(elem0, ori) + elem_faces(elem0, ori+
dimension) - i_face;
5350 if( (elem_faces(elem0, ori) != i_face) && (elem_faces(elem0, ori+
dimension) != i_face) )
5352 Cerr <<
"Attention, il y a un probleme sur la numerotation des faces. On calcule en effet" << finl ;
5353 Cerr <<
"elem_faces(elem0, ori) = " <<elem_faces(elem0, ori)<< finl ;
5354 Cerr <<
"elem_faces(elem0, ori+dimension) = " <<elem_faces(elem0, ori+
dimension)<< finl ;
5355 Cerr <<
"i_face = " <<i_face<< finl;
5356 Cerr <<
"tandis qu'au moins un des deux calculs 'elem_faces(elem0, ori)' ou 'elem_faces(elem0, ori+dimension)'"<<finl ;
5357 Cerr <<
"doit correspondre au nuemro de la face i_face"<<finl ;
5358 Cerr <<
"Attention, verifier les conditions aux limites."<<finl ;
5359 Cerr <<
"En effet, le codage suivant ne permet pas de prendre en compte des conditions" << finl ;
5360 Cerr <<
"de periodicite avec une vitesse solide imposee analytique." << finl ;
5363 const double grad0 = gradient_old(voisin0);
5364 if (grad0 > invalid_test)
5372 const int voisin1 = elem_faces(elem1, ori) + elem_faces(elem1, ori+
dimension) - i_face;
5373 if( (elem_faces(elem1, ori) != i_face) && (elem_faces(elem1, ori+
dimension) != i_face) )
5375 Cerr <<
"Attention, il y a un probleme sur la numerotation des faces. On calcule en effet" << finl ;
5376 Cerr <<
"elem_faces(elem1, ori) = " <<elem_faces(elem1, ori)<< finl ;
5377 Cerr <<
"elem_faces(elem1, ori+dimension) = " <<elem_faces(elem1, ori+
dimension)<< finl ;
5378 Cerr <<
"i_face = " <<i_face<< finl;
5379 Cerr <<
"tandis qu'au moins un des deux calculs 'elem_faces(elem1, ori)' ou 'elem_faces(elem1, ori+dimension)'"<<finl ;
5380 Cerr <<
"doit correspondre au nuemro de la face i_face"<<finl ;
5381 Cerr <<
"Attention, verifier les conditions aux limites."<<finl ;
5382 Cerr <<
"En effet, le codage suivant ne permet pas de prendre en compte des conditions" << finl ;
5383 Cerr <<
"de periodicite avec une vitesse solide imposee analytique." << finl ;
5386 const double grad1 = gradient_old(voisin1);
5387 if (grad1 > invalid_test)
5401 const int face0=elem_faces(elem,direction);
5402 const int face1=elem_faces(elem,direction+
dimension);
5403 const int elem_voisin0=faces_elem(face0,0)+faces_elem(face0,1)-elem;
5404 const int elem_voisin1=faces_elem(face1,0)+faces_elem(face1,1)-elem;
5408 double grad_voisin0=0.;
5411 double grad_voisin1=0.;
5413 if (elem_voisin0 >= 0)
5415 for (
int compo = 0; compo<
dimension; compo++)
5417 double c_gravite=xv(i_face,compo);
5418 double c_face_voisin1=xv(elem_faces(elem_voisin0,ori),compo);
5419 double c_face_voisin2=xv(elem_faces(elem_voisin0,ori+
dimension),compo);
5420 a_carre+=(c_gravite-c_face_voisin1)*(c_gravite-c_face_voisin1);
5421 b_carre+=(c_gravite-c_face_voisin2)*(c_gravite-c_face_voisin2);
5423 if (a_carre<b_carre)
5424 grad_voisin0=gradient_old(elem_faces(elem_voisin0,ori));
5426 grad_voisin0=gradient_old(elem_faces(elem_voisin0,ori+
dimension));
5427 if (grad_voisin0 > invalid_test)
5429 somme += grad_voisin0;
5433 if (elem_voisin1 >= 0)
5435 for (
int compo = 0; compo<
dimension; compo++)
5437 double c_gravite=xv(i_face,compo);
5438 double c_face_voisin3=xv(elem_faces(elem_voisin1,ori),compo);
5439 double c_face_voisin4=xv(elem_faces(elem_voisin1,ori+
dimension),compo);
5440 c_carre+=(c_gravite-c_face_voisin3)*(c_gravite-c_face_voisin3);
5441 d_carre+=(c_gravite-c_face_voisin4)*(c_gravite-c_face_voisin4);
5443 if (c_carre<d_carre)
5444 grad_voisin1=gradient_old(elem_faces(elem_voisin1,ori));
5446 grad_voisin1=gradient_old(elem_faces(elem_voisin1,ori+
dimension));
5447 if (grad_voisin1 > invalid_test)
5449 somme += grad_voisin1;
5457 gradient(i_face) = somme / coeff;
5460 gradient.echange_espace_virtuel();
5463 Debog::verifier(
"Transport_Interfaces_FT_Disc::interpoler_vitesse_face gradient",gradient);
5466 for (i_face = 0; i_face < nfaces; i_face++)
5468 if (indicatrice_face(i_face) != indic_phase)
5470 double d = dist_face_cor(i_face) ;
5471 const double grad = gradient[i_face];
5472 const double indic = indicatrice_face(i_face) ;
5474 const int interpolation = ( d < 0. ) ;
5475 const int extra_diphasique = ( variables_internes_->is_extra_diphasique_solide && indic != 0. && indic != 1. && d > 0. ) ;
5476 const int extra_solide = ( variables_internes_->is_extra_solide && indic == 1. && (d>0. || d < -1.e20) ) ;
5477 const int extrapolation = ( extra_diphasique || extra_solide ) ;
5481 for (
int j = 0; j < dim; j++)
5483 const double coord = xv(i_face,j);
5484 parser_x.
setVar(j, coord);
5485 parser_y.
setVar(j, coord);
5487 parser_z.
setVar(j, coord);
5492 v_imp(i_face) = 0. ;
5496 v_imp(i_face) = parser_x.
eval();
5499 v_imp(i_face) = parser_y.
eval();
5502 v_imp(i_face) = parser_z.
eval();
5505 champ( i_face ) = v_imp(i_face) ;
5508 if (grad > invalid_test && d > invalid_test && ( interpolation || extrapolation ) )
5510 champ(i_face) += d*grad ;
5519 champ(i_face) = v_imp(i_face) ;
5520 if (grad > invalid_test && d > invalid_test && ( interpolation || extrapolation ) )
5522 champ(i_face) += d * grad ;
5525 vitesse_imp_interp_->valeurs()(i_face)= champ( i_face ) ;
5529 vitesse_imp_interp_->valeurs()(i_face)=-3e30;
5532 champ.echange_espace_virtuel() ;
5533 Debog::verifier(
"Transport_Interfaces_FT_Disc::interpoler_vitesse_face champ",champ );
5539 const int dim,
const int ori,
5545 const ArrOfInt& index_elem = intersection.
index_elem() ;
5547 const IntTab& facettes = maillage.
facettes() ;
5549 const DoubleTab& sommets = maillage.
sommets() ;
5552 int index = index_elem[elem] ;
5554 const double tol_theta = 1e-4 ;
5555 double I_sur_arete_fa7 = 0. ;
5556 double I_sur_arete_vertex = 0. ;
5557 int I_dans_fa7 = 0 ;
5567 if( std::fabs( normale_facettes(i_facette,ori) ) > 0. )
5569 DoubleTab A(dim), B(dim), C(dim), U(dim), I(dim);
5575 for(
int i = 0 ; i<dim ; i++)
5580 for(
int i=0 ; i<dim ; i++)
5582 U(i) = normale_facettes(i_facette,i) ;
5583 A(i) = sommets(facettes(i_facette, 0),i) ;
5584 B(i) = sommets(facettes(i_facette, 1),i) ;
5586 C(i) = sommets(facettes(i_facette, 2),i) ;
5589 for(
int i =0 ; i<dim ; i++)
5592 I(ori) += (A(i) - I(i))*U(i)/U(ori) ;
5597 double FIscalU = (I(ori)-xe(ori))*U(ori) ;
5598 double FI = std::fabs( I(ori)-xe(ori) ) ;
5599 double theta = acos(FIscalU/FI) ;
5600 double theta_rel = std::fabs( theta - (2.0*atan(1.)) ) / (2.0*atan(1.)) ;
5601 double ecart_elem = std::fabs(xe(ori)-I(ori)) ;
5603 if( theta_rel >= tol_theta && ecart_elem <= 0.5*dx )
5607 double AIscalAB = (I(0)-A(0))*(B(0)-A(0)) + (I(1)-A(1))*(B(1)-A(1)) ;
5608 double ABscalAB = (B(0)-A(0))*(B(0)-A(0)) + (B(1)-A(1))*(B(1)-A(1)) ;
5609 double xx = AIscalAB / ABscalAB ;
5611 const int test1 = ( std::fabs( xx ) <= precision ) ;
5612 const int test2 = ( std::fabs( xx - 1. ) <= precision ) ;
5613 const int test3 = ( xx > precision ) ;
5614 const int test4 = ( xx < 1.-precision ) ;
5615 if( test1 || test2 )
5618 I_sur_arete_vertex += 1./2. ;
5620 else if( test3 && test4 )
5628 DoubleTab ABvectAC(dim) ;
5629 ABvectAC(0) = (B(1)-A(1))*(C(2)-A(2))-(B(2)-A(2))*(C(1)-A(1)) ;
5630 ABvectAC(1) = (B(2)-A(2))*(C(0)-A(0))-(B(0)-A(0))*(C(2)-A(2)) ;
5631 ABvectAC(2) = (B(0)-A(0))*(C(1)-A(1))-(B(1)-A(1))*(C(0)-A(0)) ;
5632 double ABvectACscalU = ABvectAC(0)*U(0) + ABvectAC(1)*U(1) + ABvectAC(2)*U(2) ;
5633 double AireABC = ABvectACscalU / 2. ;
5635 DoubleTab IBvectIC(dim) ;
5636 IBvectIC(0) = (B(1)-I(1))*(C(2)-I(2))-(B(2)-I(2))*(C(1)-I(1)) ;
5637 IBvectIC(1) = (B(2)-I(2))*(C(0)-I(0))-(B(0)-I(0))*(C(2)-I(2)) ;
5638 IBvectIC(2) = (B(0)-I(0))*(C(1)-I(1))-(B(1)-I(1))*(C(0)-I(0)) ;
5639 double IBvectICscalU = IBvectIC(0)*U(0) + IBvectIC(1)*U(1) + IBvectIC(2)*U(2) ;
5640 double AireIBC = IBvectICscalU / 2. ;
5642 DoubleTab ICvectIA(dim) ;
5643 ICvectIA(0) = (C(1)-I(1))*(A(2)-I(2))-(C(2)-I(2))*(A(1)-I(1)) ;
5644 ICvectIA(1) = (C(2)-I(2))*(A(0)-I(0))-(C(0)-I(0))*(A(2)-I(2)) ;
5645 ICvectIA(2) = (C(0)-I(0))*(A(1)-I(1))-(C(1)-I(1))*(A(0)-I(0)) ;
5646 double ICvectIAscalU = ICvectIA(0)*U(0) + ICvectIA(1)*U(1) + ICvectIA(2)*U(2) ;
5647 double AireICA = ICvectIAscalU / 2. ;
5649 DoubleTab IAvectIB(dim) ;
5650 IAvectIB(0) = (A(1)-I(1))*(B(2)-I(2))-(A(2)-I(2))*(B(1)-I(1)) ;
5651 IAvectIB(1) = (A(2)-I(2))*(B(0)-I(0))-(A(0)-I(0))*(B(2)-I(2)) ;
5652 IAvectIB(2) = (A(0)-I(0))*(B(1)-I(1))-(A(1)-I(1))*(B(0)-I(0)) ;
5653 double IAvectIBscalU = IAvectIB(0)*U(0) + IAvectIB(1)*U(1) + IAvectIB(2)*U(2) ;
5654 double AireIAB = IAvectIBscalU / 2. ;
5656 double alpha = AireIBC/AireABC ;
5657 double beta = AireICA/AireABC ;
5658 double gamma = AireIAB/AireABC ;
5660 int test1 = ( std::fabs(alpha-1.)<=precision && std::fabs(beta)<=precision && std::fabs(gamma)<=precision ) ;
5661 int test2 = ( std::fabs(beta-1.)<=precision && std::fabs(alpha)<=precision && std::fabs(gamma)<=precision ) ;
5662 int test3 = ( std::fabs(gamma-1.)<=precision && std::fabs(alpha)<=precision && std::fabs(beta)<=precision ) ;
5663 int test4 = ( std::fabs(alpha)<=precision && precision<beta && beta<1.-precision && precision<gamma && gamma<1.-precision ) ;
5664 int test5 = ( std::fabs(beta)<=precision && precision<alpha && alpha<1.-precision && precision<gamma && gamma<1.-precision ) ;
5665 int test6 = ( std::fabs(gamma)<=precision && precision<beta && beta<1.-precision && precision<alpha && alpha<1.-precision ) ;
5666 int test7 = ( precision<alpha && alpha<1.-precision) ;
5667 int test8 = ( precision<beta && beta<1.-precision) ;
5668 int test9 = ( precision<gamma && gamma<1.-precision) ;
5670 if( test1 || test2 || test3 )
5674 I_sur_arete_vertex += 1e-5 ;
5676 else if( test4 || test5 || test6 )
5679 I_sur_arete_fa7 += 1./2. ;
5681 else if( test7 && test8 && test9 )
5689 Cerr <<
" attention seules les dimension 2 et 3 sont traitees "<<finl ;
5697 traverse = int(ceil(I_sur_arete_vertex)) + int(ceil(I_sur_arete_fa7)) + I_dans_fa7 ;
5703 const int nb_elem,
int& nb_fa7_accepted,
5704 IntList& OutElem, IntList& OutFa7 )
5707 const ArrOfInt& index_elem = intersection.
index_elem() ;
5712 for (
int i_elem = 0; i_elem < nb_elem ; i_elem++)
5715 if( indicatrice(i_elem) != 0. && indicatrice(i_elem) != 1. )
5717 int index = index_elem[i_elem] ;
5731 E.
add(surfaces[i_facette]) ;
5740 if( F.
size() > nb_fa7_accepted )
5749 while( i<nb_fa7_accepted )
5752 double max_data = -1e30 ;
5757 if( E[l]>max_data && !Fa7Elem.
contient(F[l]) )
5766 Fa7Elem.
add( F[ll] ) ;
5767 OutFa7.
add( F[ll] ) ;
5777 const DoubleTab& indicatrice_face, DoubleTab& Vertex )
5784 const int nb_elem = domaine_vf.
nb_elem() ;
5785 const int nfaces = faces_elem.
dimension(0) ;
5787 const ArrOfInt& index_elem = intersection.
index_elem() ;
5788 const IntTab& facettes = maillage.
facettes() ;
5789 const DoubleTabFT& sommets = maillage.
sommets() ;
5790 const DoubleTab& xv = domaine_vf.
xv();
5791 DoubleTab coord(dim) ;
5794 for(
int i_face=0 ; i_face<nfaces ; i_face++ )
5797 double dist_fa7 = 1e+30 ;
5798 int ppp_calcule = 0 ;
5801 if( indicatrice_face(i_face) == 0.5 )
5804 if( faces_elem(i_face,0) > -1 )
5805 voisin = faces_elem(i_face,0) ;
5807 voisin = faces_elem(i_face,1) ;
5811 if( indicatrice(voisin) == 0. || indicatrice(voisin) == 1. )
5813 for(
int j=0; j<dim; j++)
5815 Vertex(i_face,j) = xv(i_face,j) ;
5821 if (indicatrice_face(i_face) != 0. && indicatrice_face(i_face) != 1. && ppp_calcule == 0 )
5824 for(
int i=0 ; i<2; i++)
5826 const int elem_voisin = faces_elem(i_face,i) ;
5828 if( elem_voisin<nb_elem && elem_voisin>-1 && indicatrice(elem_voisin) != 0. && indicatrice(elem_voisin) != 1. )
5830 int index = index_elem[elem_voisin] ;
5837 double dist_min = 1e+30 ;
5839 for(
int som=0; som<facettes.
dimension(1); som++)
5842 double dist_som = 0. ;
5843 for(
int j=0; j<dim; j++)
5846 double sj = sommets(facettes(i_facette, som), j) ;
5848 double xj = xv(i_face, j) ;
5850 dist_som += (xj - sj)*(xj - sj) ;
5853 dist_som = sqrt(dist_som) ;
5856 if( dist_som < dist_min )
5858 dist_min = dist_som ;
5859 for(
int j=0; j<dim; j++)
5861 coord(j) = sommets(facettes(i_facette, som), j) ;
5866 if( dist_min < dist_fa7 )
5868 dist_fa7 = dist_min ;
5869 for(
int j=0; j<dim; j++)
5871 Vertex(i_face,j) = coord(j) ;
5885 DoubleTab& Vertex, DoubleTab& PPP )
5892 const IntTab& elem_faces = domaine_vf.
elem_faces();
5893 const int nfaces = faces_elem.
dimension(0) ;
5894 const DoubleTab& xv = domaine_vf.
xv();
5896 for(
int i_face=0 ; i_face<nfaces ; i_face++ )
5900 int ppp_calcule = 0 ;
5903 if( indicatrice_face(i_face) == 0.5 )
5906 if( faces_elem(i_face,0) > -1 )
5907 voisin = faces_elem(i_face,0) ;
5909 voisin = faces_elem(i_face,1) ;
5914 if( indicatrice(voisin) == 0. || indicatrice(voisin) == 1. )
5916 for(
int j=0; j<dim; j++)
5918 PPP(i_face,j) = xv(i_face,j) ;
5924 if (indicatrice_face(i_face) != 0. && indicatrice_face(i_face) != 1. && ppp_calcule == 0 )
5926 for(
int i=0 ; i<dim ; i++)
5928 d += (Vertex(i_face,i)-xv(i_face,i))*(Vertex(i_face,i)-xv(i_face,i)) ;
5932 for(
int k=0 ; k<dim ; k++)
5934 PPP(i_face,k) = Vertex(i_face,k) ;
5939 for(
int i=0; i<2; i++ )
5941 const int elem = faces_elem(i_face,i) ;
5946 for (
int face_loc=0 ; face_loc<2*dim ; face_loc++)
5948 const int autre_face = elem_faces(elem,face_loc) ;
5949 if (indicatrice_face(autre_face) != 0. && (indicatrice_face(autre_face) != 1.
5950 || (indicatrice_face(autre_face) == 1. && d<0.)))
5953 for(
int j=0 ; j<dim ; j++)
5955 dd += (Vertex(autre_face,j)-xv(i_face,j))*(Vertex(autre_face,j)-xv(i_face,j)) ;
5961 for(
int k=0 ; k<dim ; k++)
5963 PPP(i_face,k) = Vertex(autre_face,k) ;
5982 const IntTab& elem_faces = domaine_vf.
elem_faces();
5983 const int nfaces = faces_elem.
dimension(0) ;
5984 const DoubleTab& xv = domaine_vf.
xv();
5986 for(
int i_face=0 ; i_face<nfaces ; i_face++ )
5989 double dmin = 1.e30 ;
5991 if (indicatrice_face(i_face) == 0. || indicatrice_face(i_face) == 1.)
5996 for(
int ii=0; ii<2; ii++ )
5998 const int elem = faces_elem(i_face,ii) ;
6003 for (
int face_loc=0 ; face_loc<2*dim ; face_loc++)
6005 const int face = elem_faces(elem,face_loc) ;
6007 if (indicatrice_face(face) != 0. && indicatrice_face(face) != 1.)
6010 for(
int i=0 ; i<2; i++)
6012 const int elem_voisin = faces_elem(face,i) ;
6016 if( elem_voisin>-1 && indicatrice(elem_voisin) != 0. && indicatrice(elem_voisin) != 1. )
6019 for (
int autre_face_loc=0 ; autre_face_loc<2*dim ; autre_face_loc++)
6021 const int autre_face = elem_faces(elem_voisin,autre_face_loc) ;
6023 for(
int j=0 ; j<dim ; j++)
6025 dd += (PPP(autre_face,j)-xv(i_face,j))*(PPP(autre_face,j)-xv(i_face,j)) ;
6031 for(
int k=0 ; k<dim ; k++)
6033 PPP(i_face,k) = PPP(autre_face,k) ;
6049 const DoubleTab& indicatrice,
6050 const int nb_elem,
int& max_fa7,
6051 const int exec_planfa7existant)
6054 const ArrOfInt& index_elem = intersection.
index_elem() ;
6058 for (
int i_elem = 0; i_elem < nb_elem ; i_elem++)
6061 if( indicatrice(i_elem) != 0. && indicatrice(i_elem) != 1. )
6063 int index = index_elem[i_elem] ;
6066 DoubleList A, B, C, D ;
6074 int test_liste = 0 ;
6077 if( !(A.
est_vide()) && exec_planfa7existant == 1 )
6084 if( exec_planfa7existant == 1 )
6086 double aa=0., bb=0., cc=0., dd=0. ;
6109 DoubleTab& Tab111, DoubleTab& Tab112,
6110 IntTab& Tab12, IntTab& CptFacette,
6111 const int nb_facettes,
6112 const int nb_facettes_dim )
6119 for (
int i_elem = nn; i_elem < nn_tot ; i_elem++)
6122 for(
int cpt=0 ; cpt < CptFacette(i_elem) ; cpt++ )
6124 int bary_trouve = 0 ;
6126 while( ii < nb_facettes && bary_trouve==0 )
6128 bool test_0 = (Barycentre(ii,0)-Tab110(i_elem,cpt))*(Barycentre(ii,0)-Tab110(i_elem,cpt))<eps ;
6129 bool test_1 = (Barycentre(ii,1)-Tab111(i_elem,cpt))*(Barycentre(ii,1)-Tab111(i_elem,cpt))<eps ;
6130 bool test_2 = (Barycentre(ii,2)-Tab112(i_elem,cpt))*(Barycentre(ii,2)-Tab112(i_elem,cpt))<eps ;
6132 if( test_0 && test_1 && test_2 )
6135 Tab12( i_elem, cpt ) = ii ;
6139 if( bary_trouve==0 && ii==nb_facettes )
6144 if( Tab12( i_elem, cpt ) < nb_facettes )
6145 Tab12( i_elem, cpt )+=nb_facettes_dim ;
6152 IntTab& CptFacette, DoubleTab& Tab100,
6153 DoubleTab& Tab101, DoubleTab& Tab102,
6154 DoubleTab& Tab103, DoubleTab& Tab110,
6155 DoubleTab& Tab111, DoubleTab& Tab112,
6156 IntTab& Tab12, DoubleTab& Barycentre,
6157 const DoubleTab& indicatrice,
6158 IntList& OutElem, ArrOfBit& fa7,
6159 const int exec_planfa7existant )
6162 const ArrOfInt& index_elem = intersection.
index_elem() ;
6168 for (
int i_elem = 0; i_elem < nn ; i_elem++)
6170 if( indicatrice(i_elem) != 0. && indicatrice(i_elem) != 1. && !OutElem.
contient(i_elem) )
6173 int index = index_elem[i_elem] ;
6174 CptFacette(i_elem) = 0 ;
6177 DoubleList A, B, C, D ;
6184 int test_liste = 0 ;
6185 if( !(A.
est_vide()) && exec_planfa7existant == 1 )
6194 Tab102(i_elem,cpt),Tab103(i_elem,cpt)) ;
6195 if( exec_planfa7existant == 1 )
6197 A.
add(Tab100(i_elem,cpt)) ;
6198 B.
add(Tab101(i_elem,cpt)) ;
6199 C.
add(Tab102(i_elem,cpt)) ;
6200 D.
add(Tab103(i_elem,cpt)) ;
6202 Tab12(i_elem,cpt) = i_facette ;
6210 BaryFa7(maillage,i_facette,Barycentre) ;
6213 Tab110(i_elem,cpt) = Barycentre(i_facette,0) ;
6214 Tab111(i_elem,cpt) = Barycentre(i_facette,1) ;
6215 Tab112(i_elem,cpt) = Barycentre(i_facette,2) ;
6222 CptFacette(i_elem) = cpt ;
6225 Cerr<<
"Attention, seulement "<< cpt <<
" facettes ont ete stokees "<< finl ;
6226 Cerr<<
"dans l'element "<< i_elem <<
" qui en contient plus."<< finl ;
6227 Cerr<<
"Regarder la finesse du maillage Lagrangien par rapport au maillage Eulerien" <<finl ;
6235 DoubleTab& Tab101, DoubleTab& Tab102, DoubleTab& Tab103,
6236 DoubleTab& Tab110,DoubleTab& Tab111, DoubleTab& Tab112,
6237 IntTab& Tab12, DoubleTab& Barycentre, IntList& OutElem,
6238 IntTab& TabOutFa7, ArrOfBit& fa7 )
6243 while( cpt_elem < OutElem.
size() )
6245 for(
int k=0 ; k<TabOutFa7.
line_size() ; k++ )
6247 const int i_facette = TabOutFa7(cpt_elem,k) ;
6249 Tab101(OutElem[cpt_elem],k),Tab102(OutElem[cpt_elem],k),
6250 Tab103(OutElem[cpt_elem],k)) ;
6251 Tab12(OutElem[cpt_elem],k) = i_facette ;
6259 BaryFa7(maillage,i_facette,Barycentre) ;
6262 Tab110(OutElem[cpt_elem],k) = Barycentre(i_facette,0) ;
6263 Tab111(OutElem[cpt_elem],k) = Barycentre(i_facette,1) ;
6264 Tab112(OutElem[cpt_elem],k) = Barycentre(i_facette,2) ;
6274 const DoubleTabFT& sommets = maillage.
sommets() ;
6275 const IntTab& facettes = maillage.
facettes() ;
6278 DoubleTab CoordSom( facettes.
line_size(), dim ) ;
6281 for (
int i=0 ; i<facettes.
line_size() ; i ++ )
6283 Bary(i_facette,i) = 0. ;
6284 Som(i) = facettes(i_facette, i) ;
6285 for (
int j = 0 ; j < dim ; j++ )
6287 CoordSom(i,j) = sommets(Som(i),j) ;
6291 for (
int k = 0 ; k< dim ; k++ )
6293 for(
int l = 0 ; l< facettes.
line_size() ; l++ )
6295 Bary(i_facette,k) += CoordSom(l,k) / double(facettes.
line_size()) ;
6299 Bary(i_facette,2) = 0. ;
6304 DoubleList B,DoubleList C,DoubleList D,
6305 const int i_facette,
int& test_liste )
6310 const DoubleTabFT& sommets = maillage.
sommets() ;
6311 const IntTab& facettes = maillage.
facettes() ;
6312 double x1, y1, z1 = 0. ;
6313 double x2, y2, z2 = 0. ;
6314 double x3 = 0., y3 = 0., z3 = 0. ;
6316 const int s1 = facettes(i_facette, 0) ;
6317 x1 = sommets(s1,0) ;
6318 y1 = sommets(s1,1) ;
6319 const int s2 = facettes(i_facette, 1) ;
6320 x2 = sommets(s2,0) ;
6321 y2 = sommets(s2,1) ;
6322 double eps = 1e-20 ;
6328 z1 = sommets(s1,2) ;
6329 z2 = sommets(s2,2) ;
6330 const int s3 = facettes(i_facette, 2) ;
6331 x3 = sommets(s3,0) ;
6332 y3 = sommets(s3,1) ;
6333 z3 = sommets(s3,2) ;
6337 while( i< A.
size() && test_liste==0 )
6339 double plan1 = D[i] + A[i] * x1 + B[i] * y1 + C[i] * z1 ;
6340 double plan2 = D[i] + A[i] * x2 + B[i] * y2 + C[i] * z2 ;
6342 plan3 = D[i] + A[i] * x3 + B[i] * y3 + C[i] * z3 ;
6345 if( std::fabs(plan1)<eps && std::fabs(plan2)<eps && std::fabs(plan3)<eps )
6353 double& a,
double& b,
double& c,
double& d )
6356 const DoubleTabFT& sommets = maillage.
sommets() ;
6357 const IntTab& facettes = maillage.
facettes() ;
6360 double x, y, z = 0. ;
6361 double inverse_norme = 0. ;
6363 const int s1 = facettes(i_facette, 0) ;
6365 a = normale_facettes(i_facette, 0) ;
6366 b = normale_facettes(i_facette, 1) ;
6372 c = normale_facettes(i_facette, 2) ;
6375 inverse_norme = 1. / sqrt( a*a + b*b + c*c ) ;
6377 a *= inverse_norme ;
6378 b *= inverse_norme ;
6379 c *= inverse_norme ;
6380 d = - (a * x + b * y + c * z) ;
6385 double& b,
double& c,
double& d )
6389 const DoubleTabFT& sommets = maillage.
sommets() ;
6390 const IntTab& facettes = maillage.
facettes() ;
6393 double x, y, z = 0. ;
6394 double inverse_norme = 0. ;
6396 const int s1 = facettes(i_facette, 0) ;
6398 a = normale_facettes(i_facette, 0) ;
6399 b = normale_facettes(i_facette, 1) ;
6405 c = normale_facettes(i_facette, 2) ;
6408 inverse_norme = 1. / sqrt( a*a + b*b + c*c ) ;
6410 a *= inverse_norme ;
6411 b *= inverse_norme ;
6412 c *= inverse_norme ;
6413 d = - (a * x + b * y + c * z) ;
6417 const int voisin1,
const DoubleTab& indicatrice,
double& tol )
6423 const IntTab& elem_faces = domaine_vf.
elem_faces();
6428 if( voisin0 > -1 && indicatrice(voisin0) != 0. && indicatrice(voisin0) != 1. )
6431 const int face_voisine = elem_faces(voisin0,ori) + elem_faces(voisin0,ori+dim) - i_face ;
6432 L(ori) = std::fabs(domaine_vdf.
distance_face(i_face,face_voisine,ori)) ;
6434 if( voisin1 > -1 && indicatrice(voisin1) != 0. && indicatrice(voisin1) != 1. )
6437 const int face_voisine = elem_faces(voisin1,ori) + elem_faces(voisin1,ori+dim) - i_face ;
6438 double xx = std::fabs(domaine_vdf.
distance_face(i_face,face_voisine,ori)) ;
6439 L(ori) = std::max( L(ori), xx ) ;
6442 for(
int k = 0 ; k<dim ; k++ )
6446 const int face0 = elem_faces(voisin,k) ;
6447 const int face1 = elem_faces(voisin,k+dim) ;
6448 L(k) = std::fabs(domaine_vdf.
distance_face(face0,face1,k))/2. ;
6451 for(
int k = 0 ; k<dim ; k++ )
6459 const int voisin1,
const DoubleTab& indicatrice_face,
6460 const DoubleTab& indicatrice,
double& tol )
6466 const IntVect& orientation = domaine_vdf.
orientation() ;
6467 const IntTab& elem_faces = domaine_vf.
elem_faces() ;
6479 const int face_voisine = elem_faces(voisin0,ori) + elem_faces(voisin0,ori+dim) - i_face ;
6480 if( indicatrice_face(face_voisine) != 0. && indicatrice_face(face_voisine) != 1. )
6482 L(ori) = std::fabs(domaine_vdf.
distance_face(i_face,face_voisine,ori)) ;
6484 const int voisin_voisin = faces_elem(face_voisine,0) + faces_elem(face_voisine,1) - voisin0 ;
6485 if( voisin_voisin > -1 )
6487 const int face_voisine_voisine = elem_faces(voisin_voisin,ori) + elem_faces(voisin_voisin,ori+dim) - face_voisine ;
6488 L(ori) += std::fabs(domaine_vdf.
distance_face(face_voisine,face_voisine_voisine,ori)) ;
6490 for(
int k = 0 ; k<dim ; k++ )
6494 const int face0 = elem_faces(voisin0,k) ;
6495 const int face1 = elem_faces(voisin0,k+dim) ;
6496 L(k) = std::fabs(domaine_vdf.
distance_face(face0,face1,k))/2. ;
6499 for(
int k = 0 ; k<dim ; k++ )
6510 const int face_voisine = elem_faces(voisin1,ori) + elem_faces(voisin1,ori+dim) -i_face ;
6511 if( indicatrice_face(face_voisine) != 0. && indicatrice_face(face_voisine) != 1. )
6513 L(ori) = std::fabs(domaine_vdf.
distance_face(i_face,face_voisine,ori)) ;
6515 const int voisin_voisin = faces_elem(face_voisine,0) + faces_elem(face_voisine,1) - voisin1 ;
6516 if( voisin_voisin > -1 )
6518 const int face_voisine_voisine = elem_faces(voisin_voisin,ori) + elem_faces(voisin_voisin,ori+dim) - face_voisine ;
6519 L(ori) += std::fabs(domaine_vdf.
distance_face(face_voisine,face_voisine_voisine,ori)) ;
6521 for(
int k = 0 ; k<dim ; k++ )
6525 const int face0 = elem_faces(voisin1,k) ;
6526 const int face1 = elem_faces(voisin1,k+dim) ;
6527 L(k) = std::fabs(domaine_vdf.
distance_face(face0,face1,k))/2. ;
6530 for(
int k = 0 ; k<dim ; k++ )
6537 double tol_long = std::max(tol0,tol1) ;
6549 const int face_voisine = elem_faces(voisin0,ori) + elem_faces(voisin0,ori+dim) - i_face ;
6550 L(ori) = std::fabs(domaine_vdf.
distance_face(i_face,face_voisine,ori)) ;
6555 const int face_voisine = elem_faces(voisin1,ori) + elem_faces(voisin1,ori+dim) - i_face ;
6556 double xx = std::fabs(domaine_vdf.
distance_face(i_face,face_voisine,ori)) ;
6557 L(ori) = std::max( L(ori), xx ) ;
6560 DoubleTab xx_max(dim) ;
6562 for (
int i=0 ; i<2*dim ; i++)
6564 const int autre_face = elem_faces(voisin,i) ;
6565 const int face_voisine = elem_faces(voisin,ori) + elem_faces(voisin,ori+dim) - i_face ;
6567 if( autre_face != i_face && autre_face != face_voisine
6568 && indicatrice_face(autre_face) != 0. && indicatrice_face(autre_face) != 1. )
6570 const int voisin_voisin = faces_elem(autre_face,0) + faces_elem(autre_face,1) - voisin ;
6571 const int ori_face = orientation[autre_face] ;
6572 const int autre_face_voisine = elem_faces(voisin_voisin,ori_face) + elem_faces(voisin_voisin,ori_face+dim) - autre_face ;
6573 double xx = std::fabs(domaine_vdf.
distance_face(autre_face,autre_face_voisine,ori_face)) ;
6574 if( xx_max(ori_face) < xx )
6575 xx_max(ori_face) = xx ;
6579 DoubleTab TOL_TMP(dim-1) ;
6582 for(
int k=0 ; k<dim ; k++)
6586 for(
int l=0 ; l<dim ; l++)
6588 if( l != k && l != ori )
6590 const int face0 = elem_faces(voisin,k) ;
6591 const int face1 = elem_faces(voisin,k+dim) ;
6592 double yy = 0.5*std::fabs(domaine_vdf.
distance_face(face0,face1,k)) ;
6593 double a = 0.5*yy + xx_max(k) ;
6594 const int face3 = elem_faces(voisin,l) ;
6595 const int face4 = elem_faces(voisin,l+dim) ;
6596 double b = 0.5*std::fabs(domaine_vdf.
distance_face(face3,face4,l)) ;
6597 TOL_TMP(cpt) = L(ori)*L(ori) + a*a + b*b ;
6603 double tol_trans = TOL_TMP(0) ;
6605 tol_trans = std::max(tol_trans,TOL_TMP(1)) ;
6606 tol_trans = sqrt( tol_trans ) ;
6610 tol = std::max(tol_long,tol_trans) ;
6615 const DoubleTab& V, DoubleTab& coord_projete,
int& cpt )
6617 const int dim = coord_projete.
dimension(0) ;
6618 const double dist_IBC = std::fabs(d) ;
6620 double dist_proj = 0. ;
6621 int projete_modifie = 0 ;
6624 DoubleTab coord_projete_old(dim) ;
6625 coord_projete_old = coord_projete ;
6627 for (
int j = 0; j < dim ; j++)
6629 dist_proj += (x(j) - coord_projete(j))*(x(j) - coord_projete(j)) ;
6631 dist_proj = sqrt( dist_proj ) ;
6632 double dist_proj_old = dist_proj ;
6634 const double facsec = 0.5 ;
6635 const int test0 = ( dist_proj < precision ) ;
6636 const int test1 = ( std::fabs( dist_IBC - dist_proj ) > facsec*dist_IBC ) ;
6637 const int test2 = ( dist_proj > (1.+facsec)*Lref ) ;
6641 if( test0 || (test1 && test2) )
6644 for (
int j = 0; j < dim ; j++)
6646 coord_projete(j) = V(j) ;
6648 projete_modifie = 1 ;
6654 const int test = ( dist_IBC >= (1.-facsec)*Lref ) ;
6656 if( (test0 && test) || (test1 && test2) )
6659 for (
int j = 0; j < dim ; j++)
6661 coord_projete(j) = V(j) ;
6663 projete_modifie = 1 ;
6668 if( projete_modifie == 1 )
6671 for (
int j = 0; j < dim ; j++)
6673 dist_proj += (x(j) - coord_projete(j))*(x(j) - coord_projete(j)) ;
6675 dist_proj = sqrt( dist_proj ) ;
6677 Cerr <<
" ------------------------------------------------------- "<<finl;
6678 Cerr <<
" le projete "<<finl;
6679 for(
int i=0; i<dim; i++)
6681 Cerr <<
" X("<<i<<
") = "<<coord_projete_old(i)<<finl ;
6683 Cerr <<
" situe a une distance "<<dist_proj_old<<
" du barycentre de la face "<<finl ;
6684 for(
int i=0; i<dim; i++)
6686 Cerr <<
" Y("<<i<<
") = "<<x(i)<<finl ;
6688 Cerr <<
" a ete modifie en "<<finl ;
6689 for(
int i=0; i<dim; i++)
6691 Cerr <<
" X("<<i<<
") = "<<coord_projete(i)<<finl ;
6693 Cerr <<
"situe a une distance "<<dist_proj<<
" du barycentre de face "<<finl ;
6694 Cerr <<
"distance a l'IBC : "<<dist_IBC<<finl ;
6695 Cerr <<
"Longueur de reference : "<<Lref<<finl ;
6696 Cerr <<
" ------------------------------------------------------- "<<finl;
6703 const DoubleTab& x,
const DoubleTab& f,
6704 DoubleTab& x_ibc)
const
6709 assert(x_ibc.
dimension(0) == nb_colonnes) ;
6715 double norme1 = 0. ;
6716 double norme_infty = 0. ;
6718 for (
int j=0; j<nb_colonnes; j++)
6720 double norme1_int = 0. ;
6721 for (
int i=0; i<nb_lignes; i++)
6723 norme1_int += std::fabs( C(i,j) ) ;
6725 norme1 = std::max( norme1 , norme1_int ) ;
6728 for (
int i=0; i<nb_lignes; i++)
6730 double norme_infty_int = 0.;
6731 for (
int j=0; j<nb_colonnes; j++)
6733 norme_infty_int += std::fabs( C(i,j) ) ;
6735 norme_infty = std::max( norme_infty , norme_infty_int ) ;
6737 rho = 1. / ( norme1 * norme_infty ) ;
6740 DoubleTab lambda(nb_lignes) ;
6743 double distance = 1. ;
6744 double distance_old = 0. ;
6745 double deplacement_relatif = 1.;
6748 double err_uzawa = variables_internes_->seuil_uzawa ;
6749 int nb_iter_max = variables_internes_->nb_iter_uzawa ;
6762 while( deplacement_relatif > err_uzawa && nb_iter <= nb_iter_max )
6766 distance_old = distance ;
6767 for (
int i=0; i<x_ibc.
dimension(0); i++)
6769 for (
int j=0; j<lambda.
dimension(0); j++)
6771 x_ibc(i) -= 0.5 * C(j,i) * lambda(j) ;
6777 for (
int i=0; i<lambda.
dimension(0); i++)
6782 for (
int j=0; j<x_ibc.
dimension(0); j++)
6784 xx += C(i,j)*x_ibc(j) ;
6790 lambda(i) = std::max( lambda(i) + rho * xx , 0. ) ;
6795 lambda(i) = std::min( lambda(i) + rho * xx , 0. ) ;
6799 distance = ( x_ibc(0) - x[0] ) * ( x_ibc(0) - x[0] )
6800 + ( x_ibc(1) - x[1] ) * ( x_ibc(1) - x[1] ) ;
6802 distance += ( x_ibc(2) - x[2] ) * ( x_ibc(2) - x[2] ) ;
6804 deplacement_relatif = std::fabs(distance - distance_old) ;
6812 const DoubleTab& indicatrice_face,
const DoubleTab& indicatrice,
6813 const DoubleTab& dist_face,
const double t,
const double dt,
6814 DoubleTab& Tab100,DoubleTab& Tab101,DoubleTab& Tab102,DoubleTab& Tab103,
6815 IntTab& Tab12,IntTab& CptFacette,DoubleTab& v_imp,DoubleTab& Vertex,
Parser& parser_x,
Parser& parser_y,
Parser& parser_z )
6823 const IntVect& orientation = domaine_vdf.
orientation();
6825 const IntTab& elem_faces = domaine_vf.
elem_faces();
6827 const int nfaces = faces_elem.
dimension(0) ;
6828 const DoubleTab& xv = domaine_vf.
xv();
6829 const double invalid_test = -1.e30;
6831 Cerr <<
"CALCUL DU PROJETE POUR LES FACES FLUIDE ET d > invalid_test" << finl ;
6833 for (
int i_face = 0; i_face < nfaces; i_face++)
6835 double d = dist_face(i_face) ;
6836 if (indicatrice_face(i_face) == 0. && d > invalid_test)
6838 DoubleList mat_0, mat_1, mat_2, mat_3 ;
6839 ArrOfBit fa7(dimfa7);
6841 int nb_facettes_to_uzawa = 0 ;
6843 for(
int ii=0; ii<2; ii++ )
6845 const int elem = faces_elem(i_face,ii) ;
6850 for (
int face_loc=0 ; face_loc<2*dim ; face_loc++)
6852 const int autre_face = elem_faces(elem,face_loc) ;
6854 if (indicatrice_face(autre_face) != 0. && indicatrice_face(autre_face) != 1.)
6858 for(
int ii2=0 ; ii2<2; ii2++)
6860 const int elem_voisin = faces_elem(autre_face,ii2) ;
6864 if( elem_voisin>-1 && indicatrice(elem_voisin) != 0.
6865 && indicatrice(elem_voisin) != 1. )
6868 for(
int i=0 ; i< CptFacette(elem_voisin) ; i++)
6870 if( !fa7.
testsetbit( Tab12(elem_voisin,i) ) )
6872 nb_facettes_to_uzawa++ ;
6873 mat_0.
add( Tab100(elem_voisin,i) ) ;
6874 mat_1.
add( Tab101(elem_voisin,i) ) ;
6875 mat_2.
add( Tab102(elem_voisin,i) ) ;
6876 mat_3.
add( Tab103(elem_voisin,i) ) ;
6885 if( nb_facettes_to_uzawa > 0 )
6887 DoubleTab matrice_plans(nb_facettes_to_uzawa, dim) ;
6889 DoubleTab secmem(nb_facettes_to_uzawa) ;
6891 for(
int i_facette=0 ; i_facette<nb_facettes_to_uzawa ; i_facette++ )
6893 matrice_plans(i_facette,0) = mat_0[i_facette] ;
6894 matrice_plans(i_facette,1) = mat_1[i_facette] ;
6896 matrice_plans(i_facette,2) = mat_2[i_facette] ;
6897 secmem[i_facette] = -mat_3[i_facette] ;
6900 for (
int i_coord=0; i_coord<dim; i_coord++)
6902 x(i_coord) = xv(i_face,i_coord) ;
6906 DoubleTab coord_projete(dim) ;
6907 coord_projete = 0. ;
6908 uzawa(d,matrice_plans,x,secmem,coord_projete) ;
6913 for (
int i_coord=0; i_coord<dim; i_coord++)
6915 V(i_coord) = Vertex(i_face,i_coord) ;
6919 const int ori = orientation[i_face] ;
6920 const int voisin0 = faces_elem(i_face,0) ;
6921 const int voisin1 = faces_elem(i_face,1) ;
6922 const int monophasique = 1 ;
6924 verifprojete(monophasique,Lref,d,x,V,coord_projete,cpt) ;
6932 for (
int j = 0; j < dim; j++)
6934 parser_x.
setVar(j, coord_projete(j)) ;
6935 parser_y.
setVar(j, coord_projete(j)) ;
6937 parser_z.
setVar(j, coord_projete(j)) ;
6944 v_imp(i_face) = parser_x.
eval() ;
6947 v_imp(i_face) = parser_y.
eval() ;
6950 v_imp(i_face) = parser_z.
eval() ;
6957 ArrOfDouble v(dim) ;
6958 v = variables_internes_->loi_horaire_->vitesse(t+dt,coord_projete);
6968 const DoubleTab& indicatrice,
const DoubleTab& dist_face,
6969 const double t,
const double dt, DoubleTab& Tab100,
6970 DoubleTab& Tab101, DoubleTab& Tab102,DoubleTab& Tab103,
6971 IntTab& Tab12,IntTab& CptFacette, DoubleTab& v_imp,
6981 const IntVect& orientation = domaine_vdf.
orientation();
6983 const IntTab& elem_faces = domaine_vf.
elem_faces();
6985 const int nfaces = faces_elem.
dimension(0) ;
6986 const DoubleTab& xv = domaine_vf.
xv() ;
6987 const double invalid_test = -1.e30;
6989 for(
int i_face=0 ; i_face<nfaces ; i_face++ )
6991 double d = dist_face(i_face) ;
6992 const int diphasique_tmp = (indicatrice_face(i_face) != 0. &&
6993 (indicatrice_face(i_face) != 1. || (indicatrice_face(i_face) == 1. && d<=0. && d> invalid_test)));
6994 const int diphasique = ( variables_internes_->is_extra_diphasique_solide ? diphasique_tmp : (diphasique_tmp && d <= 0.) ) ;
6995 const int fluide = ( indicatrice_face(i_face) == 0. ) ;
6996 const int solide = ( variables_internes_->is_extra_solide && (d>0. || (indicatrice_face(i_face) == 1. && d<-1.e20)) ) ;
6997 const int monophasique = ( fluide || solide ) ;
6999 if (diphasique || (monophasique && d> invalid_test ))
7001 DoubleList mat_0, mat_1, mat_2, mat_3 ;
7002 ArrOfBit fa7(dimfa7);
7004 int nb_facettes_to_uzawa = 0 ;
7007 DoubleTab coord_projete(dim) ;
7008 coord_projete = 0. ;
7009 int projete_calcule = 0 ;
7015 if( indicatrice_face(i_face) == 0.5 )
7018 if( faces_elem(i_face,0) > -1 )
7019 voisin = faces_elem(i_face,0) ;
7021 voisin = faces_elem(i_face,1) ;
7023 if( indicatrice(voisin) == 0. || indicatrice(voisin) == 1. )
7025 for(
int j=0 ; j<dim ; j++)
7027 coord_projete(j)=xv(i_face,j) ;
7029 projete_calcule = 1 ;
7033 if( projete_calcule == 0 )
7037 for(
int j=0 ; j<dim ; j++)
7039 coord_projete(j)=xv(i_face,j) ;
7046 for(
int ii=0 ; ii<2; ii++)
7048 const int elem_voisin = faces_elem(i_face,ii) ;
7059 if (indicatrice(elem_voisin) != 0. && indicatrice(elem_voisin) != 1. )
7061 for(
int i=0 ; i< CptFacette(elem_voisin) ; i++)
7063 if( !fa7.
testsetbit( Tab12(elem_voisin,i) ) )
7065 nb_facettes_to_uzawa++ ;
7066 mat_0.
add( Tab100(elem_voisin,i) ) ;
7067 mat_1.
add( Tab101(elem_voisin,i) ) ;
7068 mat_2.
add( Tab102(elem_voisin,i) ) ;
7069 mat_3.
add( Tab103(elem_voisin,i) ) ;
7075 for(
int j=0 ; j<2*dim ; j++)
7077 const int elem_voisin_voisin = faces_elem(elem_faces(elem_voisin,j),0)
7078 + faces_elem(elem_faces(elem_voisin,j),1) - elem_voisin ;
7080 if( elem_voisin_voisin>-1 && indicatrice(elem_voisin_voisin) != 0.
7081 && indicatrice(elem_voisin_voisin) != 1.
7082 && elem_voisin_voisin!=faces_elem(i_face,0)
7083 && elem_voisin_voisin!=faces_elem(i_face,1) )
7085 for(
int i=0 ; i< CptFacette(elem_voisin_voisin) ; i++)
7087 if( !fa7.
testsetbit( Tab12(elem_voisin_voisin,i) ) )
7089 nb_facettes_to_uzawa++ ;
7090 mat_0.
add( Tab100(elem_voisin_voisin,i) ) ;
7091 mat_1.
add( Tab101(elem_voisin_voisin,i) ) ;
7092 mat_2.
add( Tab102(elem_voisin_voisin,i) ) ;
7093 mat_3.
add( Tab103(elem_voisin_voisin,i) ) ;
7101 if( nb_facettes_to_uzawa == 0 && diphasique )
7103 Cerr <<
"erreur concernant la face : "<< i_face <<
" I_f = "<<indicatrice_face(i_face)<<finl ;
7104 Cerr <<
"aucune facette detecte "<< finl ;
7107 else if (nb_facettes_to_uzawa >0)
7110 DoubleTab matrice_plans(nb_facettes_to_uzawa, dim) ;
7112 DoubleTab secmem(nb_facettes_to_uzawa) ;
7114 for(
int i_facette=0 ; i_facette<nb_facettes_to_uzawa ; i_facette++ )
7116 matrice_plans(i_facette,0) = mat_0[i_facette] ;
7117 matrice_plans(i_facette,1) = mat_1[i_facette] ;
7119 matrice_plans(i_facette,2) = mat_2[i_facette] ;
7120 secmem[i_facette] = -mat_3[i_facette] ;
7122 for (
int i_coord=0; i_coord<dim; i_coord++)
7124 x(i_coord) = xv(i_face,i_coord) ;
7126 coord_projete = 0. ;
7127 uzawa(d,matrice_plans,x,secmem,coord_projete) ;
7132 for (
int i_coord=0; i_coord<dim; i_coord++)
7134 V(i_coord) = Vertex(i_face,i_coord) ;
7138 const int ori = orientation[i_face] ;
7139 const int voisin0 = faces_elem(i_face,0) ;
7140 const int voisin1 = faces_elem(i_face,1) ;
7141 const int monophasique2 = 0 ;
7143 verifprojete(monophasique2,Lref,d,x,V,coord_projete,cpt) ;
7151 const int test = (diphasique || fluide || solide ) ;
7159 for (
int j = 0; j < dim; j++)
7161 parser_x.
setVar(j, coord_projete(j)) ;
7162 parser_y.
setVar(j, coord_projete(j)) ;
7164 parser_z.
setVar(j, coord_projete(j)) ;
7171 v_imp(i_face) = parser_x.
eval() ;
7174 v_imp(i_face) = parser_y.
eval() ;
7177 v_imp(i_face) = parser_z.
eval() ;
7185 v = variables_internes_->loi_horaire_->vitesse(t+dt,coord_projete);
7202static void deplacer_maillage_ft_v_impose(
Noms expression_vitesse,
7205 const double dt = temps - maillage.
temps();
7206 Parser parser_x, parser_y, parser_z;
7207 init_parser_v_impose(expression_vitesse,
7208 parser_x, parser_y, parser_z, temps);
7210 const DoubleTab& sommets = maillage.
sommets();
7212 const int dimension3 = (dim == 3);
7214 const int nb_sommets = maillage.
nb_sommets();
7215 DoubleTabFT deplacement(nb_sommets, dim);
7216 for (i = 0; i < nb_sommets; i++)
7218 const double x0 = sommets(i, 0);
7219 const double y0 = sommets(i, 1);
7220 const double z0 = dimension3 ? sommets(i, 2) : 0.;
7221 double x = x0, y = y0, z = z0;
7222 integrer_vitesse_imposee(parser_x, parser_y, parser_z,
7223 temps, dt, x, y, z);
7224 deplacement(i, 0) = x - x0;
7225 deplacement(i, 1) = y - y0;
7227 deplacement(i, 2) = z - z0;
7243 const IntTab& facettes=maillage.
facettes ();
7244 const int dim3 = (deplacement.
line_size() == 3);
7245 ArrOfInt compo_connexe_facettes(nb_facettes);
7246 int n = search_connex_components_local_FT(maillage, compo_connexe_facettes);
7247 int nb_compo_tot=compute_global_connex_components_FT(maillage, compo_connexe_facettes, n);
7249 DoubleTabFT normale;
7250 calculer_normale_sommets_interface(maillage, normale);
7255 deplacement, Vitesses, Positions);
7261 IntVect compo_sommets;
7265 const int dim = deplacement.
dimension(1);
7266 for (
int iface = 0; iface < nb_facettes; iface++)
7268 const int compo = compo_connexe_facettes[iface];
7269 for (
int j = 0; j < dim; j++)
7270 compo_sommets[facettes(iface, j)] = compo;
7278 const int nb_sommets_tot = compo_sommets.
size_totale();
7279 for (
int som = 0; som < nb_sommets_tot; som++)
7284 deplacement(som, 0) = 100.;
7285 deplacement(som, 1) = 100.;
7289 const int compo = compo_sommets[som];
7307 double nx = 0., ny = 0., nz = 0.;
7308 double vi_x = 0., vi_y = 0., vi_z = 0.;
7310 nx = normale(som, 0);
7311 ny = normale(som, 1);
7313 vi_x = deplacement(som, 0);
7314 vi_y = deplacement(som, 1);
7317 vi_z = deplacement(som, 2);
7322 (vi_x - Vitesses(compo, 0)) * nx
7323 + (vi_y - Vitesses(compo, 1)) * ny;
7325 prodscal += (vi_z - Vitesses(compo, 2)) * nz;
7326 double norme_carre = nx * nx + ny * ny + nz * nz;
7329 if (norme_carre != 0.)
7331 prodscal /= norme_carre;
7332 deplacement(som, 0) = nx * prodscal* (1-is_solid_particle) + Vitesses(compo, 0);
7333 deplacement(som, 1) = ny * prodscal* (1-is_solid_particle) + Vitesses(compo, 1);
7335 deplacement(som, 2) = nz * prodscal* (1-is_solid_particle) + Vitesses(compo, 2);
7347 DoubleTab& deplacement = variables_internes_->deplacement_sommets;
7348 const Navier_Stokes_std& eqn_hydraulique = variables_internes_->refequation_vitesse_transport.valeur();
7371#if DEBUG_CONSERV_VOLUME
7376 DoubleTab norm_deplacement(n);
7379 for (
int i=0; i<n; i++)
7382 for (
int j=0; j<deplacement.
dimension(1); j++)
7383 x += deplacement(i,j)*deplacement(i,j);
7384 norm_deplacement[i] = sqrt(x);
7388 dmin = min_array(norm_deplacement);
7389 dmax = max_array(norm_deplacement);
7396#if DEBUG_CONSERV_VOLUME
7399 double dmean_wpch = 0.;
7400 for (
int i=0; i<n; i++)
7403 for (
int j=0; j<deplacement.
dimension(1); j++)
7404 x += deplacement(i,j)*deplacement(i,j);
7405 norm_deplacement[i] = sqrt(x);
7406 dmean_wpch +=sqrt(x);
7409 double dmin_wpch = min_array(norm_deplacement);
7410 double dmax_wpch = max_array(norm_deplacement);
7411 Cerr <<
"Interfacial_velocity before/after pch [min/mean/max]. Time: "
7413 <<
" Without-pch: " << dmin <<
" " << dmean<<
" " << dmax
7414 <<
" With-pch: " << dmin_wpch <<
" " << dmean_wpch <<
" " << dmax_wpch
7415 <<
" [Warning, values invalid in //] " << finl;
7419 if (interpolation_repere_local_)
7421 DoubleTab Positions,Vitesses;
7428 fout.open(
"composantes_connexes.txt",ios::app);
7429 fout <<
"TEMPS: " << temps << std::endl;
7430 const int nb_bulles=Positions.
dimension(0);
7432 for(
int bulle=0; bulle<nb_bulles; bulle++)
7434 fout <<
"composante " << bulle <<
" position " ;
7435 for (
int i=0; i <Positions.
dimension(1); i++) fout <<
" " << Positions(bulle,i);
7436 fout <<
" vitesse " ;
7437 for (
int i=0; i <Vitesses.
dimension(1); i++) fout <<
" " << Vitesses(bulle,i);
7444 const double delta_t = temps - maillage.
temps();
7446 deplacement *= delta_t;
7449#if DEBUG_CONSERV_VOLUME
7452 if (variables_internes_->VOFlike_correction_volume > 0)
7456 DoubleTabFT position_precedente = maillage.
sommets();
7461 ArrOfDoubleFT var_volume;
7471 const DoubleTab& val_champ_vitesse = (bool)(explicit_u_NS_) ? champ_vitesse.
futur() : champ_vitesse.
valeurs();
7474#if DEBUG_CONSERV_VOLUME
7476 double sum_before_rm = 0.;
7477 double sum_before_rm_dvol = 0.;
7478 for (
int i = 0; i < nb_elem; i++)
7479 sum_before_rm +=-dI_dt[i];
7481 sum_before_rm_dvol +=sum_before_rm*delta_t;
7486 ramasse_miettes(maillage, variables_internes_->tmp_flux->valeurs(), dI_dt);
7489#if DEBUG_CONSERV_VOLUME
7491 for (
int i = 0; i < nb_elem; i++)
7494 const double dvoldt_totale =
remaillage_interface().calculer_somme_dvolume(maillage, var_volume);
7495 const double sum_dvol =sum*delta_t;
7496 Cerr <<
" time " << temps <<
" sum_dI_dt " << sum <<
" sum_dvol " << sum_dvol
7497 <<
" sum_before_rm_dI_dt " << sum_before_rm <<
" sum_before_rm_dvol " << sum_before_rm_dvol
7498 <<
" V_lagrangien= " << dvoldt_totale <<finl;
7505 var_volume *= -delta_t;
7507#if DEBUG_CONSERV_VOLUME
7523 ArrOfDoubleFT var_volume_deplacement;
7525 position_precedente,
7526 var_volume_deplacement);
7527#if DEBUG_CONSERV_VOLUME
7529 double dvol_reel_depl =
remaillage_interface().calculer_somme_dvolume(maillage, var_volume_deplacement);
7530 Cerr <<
"Transport_Interfaces_FT_Disc::calculer_vitesse_repere_local " << finl
7531 <<
" volume avt= " << volume_avt << finl
7532 <<
" apres= " << volume_apres << finl
7533 <<
" dvol_theo_depl= " << dvol_theo_depl << finl
7534 <<
" dvol_reel_depl= " << dvol_reel_depl << finl
7539 var_volume -= var_volume_deplacement;
7542 if (variables_internes_->volume_impose_phase_1 > 0.)
7544 DoubleVect values(3);
7549 if (variables_internes_->nom_domaine_volume_impose_ ==
"??")
7554 const DoubleVect& volumes = domaine_vf.
volumes();
7556 const int nb_elem_sous_domaine = sous_domaine.
nb_elem_tot();
7557 const DoubleTab& indic = indicatrice_->valeurs();
7558 const int nb_elem = domaine_vf.
nb_elem();
7559 for (
int i = 0; i < nb_elem_sous_domaine; i++)
7561 const int elem = sous_domaine[i];
7564 values(0) += indic(elem) * volumes(elem);
7565 values(1) += volumes(elem);
7566 values(2) += (1.-indic(elem)) * volumes(elem);
7573 Cerr <<
"Volume_sous_domaine " << variables_internes_->nom_domaine_volume_impose_ <<
7577 const double erreur_volume = variables_internes_->volume_impose_phase_1 - values(0);
7578 Journal() <<
"Transport_Interfaces_FT_Disc::mettre_a_jour "
7579 <<
"correction_volume_impose_phase_1= " << erreur_volume << finl;
7582 double surface_tot = 0.;
7588 for (i = 0; i < nb_facettes; i++)
7589 surface_tot += surfaces[i];
7591 surface_tot =
mp_sum(surface_tot) * 3.;
7592 if (surface_tot > 0.)
7594 const double inv_surface_tot = 1. / surface_tot;
7595 const IntTab& facettes = maillage.
facettes();
7596 const int nb_som_facette = facettes.
dimension(1);
7597 for (i = 0; i < nb_facettes; i++)
7600 const double dvolume = erreur_volume * surfaces[i] * inv_surface_tot;
7601 for (j = 0; j < nb_som_facette; j++)
7603 int isom = facettes(i,j);
7605 var_volume[isom] -= dvolume;
7617 variables_internes_->nb_lissage_correction_volume);
7619#if DEBUG_CONSERV_VOLUME
7621 double dvol_total_before =
remaillage_interface().calculer_somme_dvolume(maillage, var_volume);
7622 Cerr <<
" dvol_error_before= " << dvol_total_before <<
" Volume_before= " << volume_before <<
" time " << temps << finl;
7628 variables_internes_->nb_iterations_correction_volume);
7629#if DEBUG_CONSERV_VOLUME
7632 Cerr <<
" dvol_error_after= " << dvol_total_after <<
" Volume_after= " << volume_after <<
" time " << temps << finl;
7637#if DEBUG_CONSERV_VOLUME
7639 DoubleTabFT position_precedente = maillage.
sommets();
7646#if DEBUG_CONSERV_VOLUME
7648 ArrOfDoubleFT var_volume_deplacement;
7650 position_precedente,
7651 var_volume_deplacement);
7657 double dvol_reel_depl =
remaillage_interface().calculer_somme_dvolume(maillage, var_volume_deplacement);
7658 Cerr <<
"Transport_Interfaces_FT_Disc::calculer_vitesse_repere_local " << finl
7659 <<
" volume avt= " << volume_avt << finl
7660 <<
" apres= " << volume_apres << finl;
7662 Cerr <<
" dvol_reel_depl= " << dvol_reel_depl << finl
7663 <<
" dt= " << delta_t << finl;
7669 if (collision_model_ && nb_particles_tot_>1)
7677 const Equation_base& eqn_hydraulique = variables_internes_->refequation_vitesse_transport.valeur();
7690 DoubleTabFT d2(deplacement);
7697 la_roue_de_vitesse_a_deja_tournee);
7701 for (
int i = 0; i < n; i++)
7703 for (
int j = 0; j < dim; j++)
7705 const double depl2 = d2(i,j);
7706 deplacement(i,j) -= depl2;
7725 flags_compo_a_supprimer, maillage,
7726 variables_internes_->indicatrice_cache->valeurs());
7732 bool indic_updated =
false;
7738 indic_updated =
true;
7741 ArrOfInt liste_elems_sous_domaine;
7745 const int nb_elems_sous_domaine = sous_domaine.
nb_elem_tot();
7746 for (i = 0; i < nb_elems_sous_domaine; i++)
7748 const int i_elem = sous_domaine[i];
7749 const double indic = indicatrice[i_elem];
7750 if (indic != phase_continue)
7753 const int sz_liste_elems_sous_domaine = liste_elems_sous_domaine.
size_array();
7754 if (
mp_sum(sz_liste_elems_sous_domaine) > 0)
7760 ArrOfInt flags_compo_a_supprimer(nb_compo);
7761 for (i = 0; i < sz_liste_elems_sous_domaine; i++)
7763 const int elem = liste_elems_sous_domaine[i];
7764 const int compo = num_compo[elem];
7765 flags_compo_a_supprimer[compo] = 1;
7770 variables_internes_->indicatrice_cache->valeurs());
7772 if (max_array(flags_compo_a_supprimer))
7775 indic_updated =
true;
7780 for (
int ii = 0; ii < n; ii++)
7787 const int i_phase_continue = (phase_continue < 0.5) ? 0 : 1;
7794 return indic_updated;
7806 Cerr<<
"The method Transport_Interfaces_FT_Disc::integrer_ensemble_lagrange"<<finl;
7807 Cerr<<
"must be overloaded."<<finl;
7813 Process::Journal() <<
"Transport_Interfaces_FT_Disc::mettre_a_jour " <<
le_nom() <<
" temps= "<< temps << finl;
7824 bool update_indic = !is_indic_up_to_date;
7832 return is_indicatrice_up_to_date;
7838 switch (variables_internes_->methode_transport)
7842 deplacer_maillage_ft_v_impose(variables_internes_->expression_vitesse_imposee, maillage, temps);
7853 const double dt = temps - maillage.
temps();
7856 const double t = maillage.
temps() + dt * 0.5;
7859 Cerr <<
"Transport_Interfaces_FT_Disc::mettre_a_jour imposed motion by a time law at t:" << t << finl;
7861 const int nb_sommets = maillage.
nb_sommets();
7864 const DoubleTab& sommets = maillage.
sommets();
7867 ArrOfDouble coord(dim);
7868 ArrOfDouble coord_barycentre(dim);
7869 int nb_sommets_reels=0;
7870 for (
int i = 0; i < nb_sommets; i++)
7873 int est_un_sommet_reel=0;
7877 est_un_sommet_reel=1;
7880 for (
int j = 0; j < dim; j++)
7882 double pos = sommets(i, j);
7884 if (est_un_sommet_reel)
7885 coord_barycentre[j] += pos;
7888 coord = variables_internes_->loi_horaire_->position(t+dt,t,coord);
7891 for (
int j = 0; j < dim; j++)
7892 deplacement(i, j) = coord[j] - sommets(i,j);
7902 for (
int j = 0; j < dim; j++)
7904 coord_barycentre[j]=coord_barycentre[j]/nb_sommets_reels;
7908 variables_internes_->loi_horaire_->imprimer(
schema_temps(),coord_barycentre);
7917 Cerr <<
"Error for the method Transport_Interfaces_FT_Disc::mettre_a_jour :\n"
7918 <<
" the requested transport mehod is not developped" << finl;
7938 return indic_updated;
7943 bool indic_updated =
false;
7945 if (variables_internes_->injection_interfaces_temps_.size_array() > 0)
7947 const ArrOfDouble& tps = variables_internes_->injection_interfaces_temps_;
7949 const Noms& expr = variables_internes_->injection_interfaces_expressions_;
7950 assert(expr.size() == n);
7952 const double last_time = variables_internes_->injection_interfaces_last_time_;
7954 for (i = 0; i < n; i++)
7955 if (tps[i] > last_time)
7957 for (; i < n && tps[i] <= temps; i++)
7959 variables_internes_->injection_interfaces_last_time_=tps[i];
7965 DoubleTab sauvegarde(variables_internes_->indicatrice_cache->valeurs());
7968 variables_internes_->indicatrice_cache->valeurs(),
7970 variables_internes_->distance_interface_sommets);
7974 Cerr <<
"Injection_interface time " << temps <<
" " << expr[i];
7979 indic_updated =
true;
7980 double unused_vol_phase_0 = 0.;
7982 unused_vol_phase_0= 0.;
7984 double volume = volume_phase_1-volume_phase_1_old;
7986 volume*=pow(-1.,1-phase);
7987 Cerr <<
" volume " << volume << finl;
7991 Cerr <<
" failure: collision" << finl;
7992 variables_internes_->indicatrice_cache->valeurs() = sauvegarde;
7993 indic_updated =
false;
7998 Cerr <<
"The injection list is fulfilled !"
7999 <<
" (if it is wished, add an interface at a large time in the file)"
8005 return indic_updated;
8010 bool indic_updated =
false;
8022 const Nom expr =
Nom(
"x^2+(y-") +
Nom(Rc)
8023 +
Nom(
"*cos(") +
Nom(thetac)
8024 +
Nom(
"*pi/180.0))^2-") +
Nom(Rc)
8036 DoubleTab sauvegarde (
8037 variables_internes_->indicatrice_cache.valeur ().valeurs ());
8040 expr, 0., maillage_tmp,
8041 variables_internes_->indicatrice_cache.valeur ().valeurs (), phase,
8042 variables_internes_->distance_interface_sommets);
8044 Cerr <<
"Injection_interface time " << temps <<
" " << expr;
8049 indic_updated =
true;
8050 double unused_vol_phase_0 = 0.;
8052 sauvegarde, unused_vol_phase_0);
8053 unused_vol_phase_0 = 0.;
8055 variables_internes_->indicatrice_cache.valeur ().valeurs (),
8056 unused_vol_phase_0);
8057 double volume = volume_phase_1 - volume_phase_1_old;
8059 volume *= pow (-1., 1 - phase);
8060 Cerr <<
" volume " << volume << finl;
8065 t_injection_ = t_present_;
8069 Cerr <<
" failure: collision" << finl;
8070 variables_internes_->indicatrice_cache.valeur ().valeurs () =
8072 indic_updated =
false;
8079 for (
int ii = 0; ii < n; ii++)
8094 const Nom expr =
Nom(
"x^2+(y-") +
Nom(Rc)
8095 +
Nom(
"*cos(") +
Nom(thetac)
8096 +
Nom(
"*pi/180.0))^2-") +
Nom(Rc)
8108 DoubleTab sauvegarde (
8109 variables_internes_->indicatrice_cache.valeur ().valeurs ());
8112 expr, 0., maillage_tmp,
8113 variables_internes_->indicatrice_cache.valeur ().valeurs (), phase,
8114 variables_internes_->distance_interface_sommets);
8116 Cerr <<
"Injection_interface time " << temps <<
" " << expr;
8121 indic_updated =
true;
8122 double unused_vol_phase_0 = 0.;
8124 sauvegarde, unused_vol_phase_0);
8125 unused_vol_phase_0 = 0.;
8127 variables_internes_->indicatrice_cache.valeur ().valeurs (),
8128 unused_vol_phase_0);
8129 double volume = volume_phase_1 - volume_phase_1_old;
8131 volume *= pow (-1., 1 - phase);
8132 Cerr <<
" volume " << volume << finl;
8143 t_injection_ = t_present_;
8149 Cerr <<
" failure: collision" << finl;
8150 variables_internes_->indicatrice_cache.valeur ().valeurs () =
8152 indic_updated =
false;
8158 return indic_updated;
8162 ,
const bool update_indic)
8172 variables_internes_->distance_interface->changer_temps(temps);
8173 variables_internes_->normale_interface->changer_temps(temps);
8176 variables_internes_->indicatrice_cache->changer_temps(temps);
8177 indicatrice_->changer_temps(temps);
8184 indicatrice_faces_->changer_temps(temps);
8196 double volume_phase_0 = 0.;
8200 Cerr <<
"Volume_phase_0 " <<
Nom(volume_phase_0,
"%20.14g") <<
" time " << temps << finl;
8201 Cerr <<
"Volume_phase_1 " <<
Nom(volume_phase_1,
"%20.14g") <<
" time " << temps << finl;
8215 for (i = 0; i < nb_facettes; i++)
8221 const double s_tot =
mp_sum(s);
8224 Cerr <<
"Surface_Totale_Interface " <<
le_nom()
8225 <<
" t= " << temps <<
" surface= " << s_tot << finl;
8230 const DoubleTab& indic = indicatrice_->valeurs();
8232 const DoubleTab& xp = domaine_vf.
xp();
8234 const DoubleVect& volumes = domaine_vf.
volumes();
8236 const int nb_elem = domaine_vf.
nb_elem();
8238 DoubleTrav values(3,dim);
8249 for (i = 0; i < nb_elem; i++)
8251 const double ind = indic(i);
8252 const double v = volumes(i);
8253 const double v1 = ind * v;
8254 const double v0 = (1. - ind) * v;
8258 for (j = 0; j < dim; j++)
8260 const double x = xp(i,j);
8261 values(0,j) += v0 * x;
8262 values(1,j) += v1 * x;
8266 double mp_g0[3] = {0., 0., 0.};
8267 double mp_g1[3] = {0., 0., 0.};
8280 if (values(2,0) > 0.)
8283 for ( j = 0; j<dim; j++)
8284 mp_g0[j] = values(0,j) / values(2,0);
8286 if (values(2,1) > 0.)
8287 for ( j = 0; j<dim; j++)
8288 mp_g1[j] = values(1,j) / values(2,1);
8294 Cerr <<
"Centre_gravite_phases " <<
le_nom() <<
" t= " << temps;
8295 Cerr <<
" phase0: " << mp_g0[0] <<
" " << mp_g0[1] <<
" " << mp_g0[2] ;
8296 Cerr <<
" phase1: " << mp_g1[0] <<
" " << mp_g1[1] <<
" " << mp_g1[2] << finl;
8303 const ArrOfInt& index_elem = intersections.
index_elem();
8304 DoubleTab& surface = variables_internes_->surface_interface->valeurs();
8305 const int nb_elements = surface.
dimension(0);
8306 for (
int element = 0; element < nb_elements; element++)
8308 int index = index_elem[element];
8309 double surface_totale = 0.;
8317 surface[element] = surface_totale;
8320 variables_internes_->surface_interface->mettre_a_jour(temps);
8329 Equation_base& eqn_hydraulique = variables_internes_->refequation_vitesse_transport.valeur();
8364 DoubleTabFT deplacement(vitesse);
8365 deplacement *= coeff;
8383 variables_internes_->indicatrice_cache->changer_temps(temps);
8384 indicatrice_->changer_temps(temps);
8385 variables_internes_->distance_interface->changer_temps(temps);
8386 variables_internes_->normale_interface->changer_temps(temps);
8400 ref_milieu_ = un_milieu;
8412 Cerr <<
"You forgot to associate the diphasic fluid to the problem named " <<
probleme().
le_nom() << finl;
8415 return ref_milieu_.valeur();
8422 Cerr <<
"You forgot to associate the diphasic fluid to the problem named " <<
probleme().
le_nom() << finl;
8425 return ref_milieu_.valeur();
8449 return indicatrice_;
8454 return indicatrice_;
8459 return variables_internes_->marching_cubes_;
8464 return variables_internes_->marching_cubes_;
8469 return variables_internes_->maillage_interface;
8474 return variables_internes_->maillage_interface;
8478 return variables_internes_->remaillage_interface_;
8483 return variables_internes_->remaillage_interface_;
8487 return variables_internes_->topologie_interface_;
8492 return variables_internes_->topologie_interface_;
8497 return variables_internes_->doubletab_pos;
8502 return variables_internes_->intvect_elements;
8507 return variables_internes_->deplacement_sommets;
8512 return variables_internes_->proprietes_particules_;
8517 return variables_internes_->proprietes_particules_;
8522 return variables_internes_->maillage_inject_;
8527 return variables_internes_->maillage_inject_;
8532 return variables_internes_->proprietes_inject_;
8537 return variables_internes_->proprietes_inject_;
8551 std::vector<YAML_data> dvi = variables_internes_->data_a_sauvegarder();
8552 data.insert(data.end(), dvi.begin(), dvi.end());
8563 int special, afaire;
8566 Nom mon_ident(
"variables_internes_transport");
8567 mon_ident +=
Nom(temps,
"%e");
8568 variables_internes_->set_time(temps);
8573 os << mon_ident << finl;
8574 os << variables_internes_->que_suis_je() << finl;
8579 os << mon_ident << finl;
8580 os << variables_internes_->que_suis_je() << finl;
8582 bytes += variables_internes_->sauvegarder(os);
8585 Cerr <<
"Backup of collision model" << finl;
8586 if (collision_model_) bytes += collision_model_.valeur().sauvegarder(os);
8589 Cerr <<
"Backup of particles position and velocity" << finl;
8604 Cerr <<
"Backup of collision model not supported with pdi format" << finl;
8623 Cerr <<
"Transport_Interfaces_FT_Disc::reprendre" << finl;
8629 is >>
id >> type_name;
8630 if ( (!
id.debute_par(
"variables_internes_transport"))
8631 || type_name != variables_internes_->
que_suis_je())
8633 Cerr <<
"Error for the method Transport_Interfaces_FT_Disc::reprendre" << finl;
8634 Cerr << variables_internes_->que_suis_je() <<
" was expected."<< finl;
8638 variables_internes_->maillage_interface.associer_equation_transport(*
this);
8640 variables_internes_->set_restart_fname(restart_fname);
8641 variables_internes_->set_time(
schema_temps().temps_courant());
8646 if (collision_model_) collision_model_.valeur().reprendre(is);
8647 const int nb_particles_tot=collision_model_.valeur().get_nb_particles_tot();
8663 Cerr <<
"Backup of collision model not supported with pdi format" << finl;
8680 return probleme_base_.valeur();
8685 return variables_internes_->parcours_interface_;
8690 return variables_internes_->connectivite_frontieres_;
8695 return variables_internes_->algorithmes_transport_.valeur();
8700 my_map& map_post, DoubleTab *ftab)
const
8702 const Motcle som =
"sommets";
8703 const Motcle elem =
"elements";
8704 const Motcle bi =
"elements et sommets";
8705 const DoubleTab dummytab;
8901 map_post.emplace(
"heat_transfer", map_element_post_FT(elem, &Transport_Interfaces_FT_Disc::fill_ftab_scalar,
8920 const Motcle som =
"sommets";
8921 const Motcle elem =
"elements";
8922 const Motcle bi =
"elements et sommets";
8924 Transport_Interfaces_FT_Disc::my_map map_post;
8925 fill_map_post_FT(map_post,ftab);
8929 map_element_post_FT::func_type the_function=&Transport_Interfaces_FT_Disc::fill_ftab_scalar;
8930 DoubleTab the_tab_values;
8931 bool elem_found=
false;
8934 for (
const auto& [key, value] : map_post)
8939 the_loc=value.location_;
8940 the_function=value.function_;
8941 the_tab_values=value.values_;
8949 Cerr<<
"The real fields to be post-processed are :"<<finl;
8951 for (
const auto& [key, value] : map_post)
8953 Cerr <<
" Fields("<<i<<
") : " << key <<
", Location : ";
8955 Cerr << value.location_ << finl;
8960 else if (!elem_found)
8962 else if (! (the_loc==bi
8976 (this->*the_function)(ftab,the_tab_values);
8993 const Motcle som =
"sommets";
8994 const Motcle elem =
"elements";
8995 const Motcle bi =
"elements et sommets";
8996 const int nb_champs = 5;
9001 fields[2] =
"numero";
9002 fields[3] =
"pe_local";
9003 fields[4] =
"compo_connexe";
9005 Motcles localisations(nb_champs);
9007 localisations[0] = bi;
9008 localisations[1] = bi;
9009 localisations[2] = bi;
9010 localisations[3] = bi;
9011 localisations[4] = elem;
9014 int rank=fields.
search(champ), i;
9018 Cerr<<
"The integer fields to be post-processed are:"<<finl;
9019 for (i=1 ; i<nb_champs ; i++)
9021 Cerr <<
" Fields("<<i<<
") : "<< fields[i] <<
" # Localisations : " << localisations[i] << finl;
9030 else if (! (localisations[rank]==bi
9059 for (i2=0 ; i2<n ; i2++)
9061 (*itab)(i2) = pe_som[i2];
9068 for (i2=0 ; i2<n ; i2++)
9070 (*itab)(i2) = pe_fac[i2];
9077 for (i2=0 ; i2<n ; i2++)
9093 int n2 = search_connex_components_local_FT(maillage, compo);
9094 compute_global_connex_components_FT(maillage, compo, n2);
9096 for (
int ii = 0; ii < nbf; ii++)
9097 (*itab)[ii] = compo[ii];
9102 Cerr <<
"Transport_Interfaces_FT_Disc::get_champ_post_FT : unexpected case" << finl;
9124 return variables_internes_->n_iterations_distance;
9138 if (tag == variables_internes_->distance_sommets_cache_tag)
9140 return variables_internes_->distance_interface_sommets;
9142 variables_internes_->distance_sommets_cache_tag = tag;
9143 assert(tag == variables_internes_->distance_normale_cache_tag);
9146 DoubleTab& dist_som = variables_internes_->distance_interface_sommets;
9181 const DoubleTab& dist_elem,
9182 const DoubleTab& normale_elem,
9183 DoubleTab& dist_som)
const
9185 static const double distance_sommets_invalides = -1.e30;
9189 const DoubleTab& xp = domaine_vf.
xp();
9193 ArrOfInt ncontrib(nb_sommets);
9198 double centre[3] = {0., 0., 0.};
9199 double normale[3] = {0., 0., 0.};
9203 const int nb_som_elem = elem_som.
line_size();
9205 for (elem = 0; elem < nb_elem_tot; elem++)
9207 const double d1 = dist_elem(elem);
9209 if (d1 < distance_sommets_invalides)
9212 for (i = 0; i < dim; i++)
9214 centre[i] = xp(elem, i);
9215 normale[i] = normale_elem(elem, i);
9218 for (i = 0; i < nb_som_elem; i++)
9220 const int som = elem_som(elem, i);
9223 if (som < nb_sommets)
9227 for (j = 0; j < dim; j++)
9229 double position_sommet = coord_som(som, j);
9230 d2 += (position_sommet - centre[j]) * normale[j];
9232 dist_som(som) += d1 + d2;
9238 const double valeur_invalide = distance_sommets_invalides * 1.1;
9239#ifdef __INTEL_COMPILER
9242 for (i = 0; i < nb_sommets; i++)
9244 const int n = ncontrib[i];
9248 dist_som(i) = valeur_invalide;
9251 Debog::verifier(
"Transport_Interfaces_FT_Disc::calculer_distance_interface_sommets",dist_som);
9271 DoubleTab& distance_elements,
9272 DoubleTab& normale_elements,
9273 const int n_iter)
const
9275 statistics().create_custom_counter(
"Calculer_distance_interface",3,
"FrontTracking");
9276 statistics().begin_count(
"Calculer_distance_interface",statistics().get_last_opened_counter_level()+1);
9278 static const double distance_sommets_invalides = -1.e30;
9283 const DoubleTab& centre_element = domaine_vf.
xp();
9287 distance_elements = distance_sommets_invalides * 1.1;
9288 normale_elements = 0.;
9301 const ArrOfInt& index_elem = intersections.
index_elem();
9304 const IntTab& facettes = maillage.
facettes();
9305 const DoubleTab& sommets = maillage.
sommets();
9307 for (
int elem = 0; elem < nb_elem; elem++)
9309 int index = index_elem[elem];
9311 double normale[3] = {0., 0., 0.};
9313 double centre[3] = {0., 0., 0.};
9315 double surface_tot = 0.;
9323#ifdef AVEC_BUG_SURFACES
9324 const double surface = data.surface_intersection_;
9328 surface_tot += surface;
9329 for (
int i = 0; i < dim; i++)
9331 normale[i] += surface * normale_facettes(num_facette, i);
9334 for (
int j = 0; j < dim; j++)
9336 const int som = facettes(num_facette, j);
9337 const double coord = sommets(som, i);
9339 g_i += coord * coeff;
9341 centre[i] += surface * g_i;
9345 if (surface_tot > 0.)
9350 const double inverse_surface_tot = 1. / surface_tot;
9353 for (j = 0; j < dim; j++)
9355 norme += normale[j] * normale[j];
9356 normale_elements(elem, j) = normale[j];
9357 centre[j] *= inverse_surface_tot;
9361 double i_norme = 1./sqrt(norme);
9362 double distance = 0.;
9363 for (j = 0; j < dim; j++)
9365 double n_j = normale[j] * i_norme;
9366 distance += (centre_element(elem, j) - centre[j]) * n_j;
9368 distance_elements(elem) = distance;
9374 Debog::verifier(
"Transport_Interfaces_FT_Disc::calculer_distance_interface",normale_elements);
9375 Debog::verifier(
"Transport_Interfaces_FT_Disc::calculer_distance_interface",distance_elements);
9378 DoubleTab terme_src(normale_elements);
9379 DoubleTab tmp(normale_elements);
9381 const IntTab& face_voisins = domaine_vf.
face_voisins();
9382 const IntTab& elem_faces = domaine_vf.
elem_faces();
9383 const int nb_elem_voisins = elem_faces.
line_size();
9387 for (iteration = 0; iteration < n_iter; iteration++)
9396 const double un_sur_ncontrib = 1. / (1. + nb_elem_voisins);
9398 for (elem = 0; elem < nb_elem; elem++)
9401 double n[3] = {0., 0., 0.};
9402 for (i = 0; i < dim; i++)
9403 n[i] = normale_elements(elem, i);
9404 for (k = 0; k < nb_elem_voisins; k++)
9407 const int face = elem_faces(elem, k);
9408 const int e_voisin = face_voisins(face, 0) + face_voisins(face, 1) - elem;
9410 for (i = 0; i < dim; i++)
9411 n[i] += normale_elements(e_voisin, i);
9413 for (i = 0; i < dim; i++)
9414 tmp(elem, i) = terme_src(elem, i) + n[i] * un_sur_ncontrib;
9416 normale_elements = tmp;
9421 ArrOfIntFT liste_elements;
9424 for (elem = 0; elem < nb_elem; elem++)
9426 double nx = normale_elements(elem, 0);
9427 double ny = normale_elements(elem, 1);
9428 double nz = (dim==3) ? normale_elements(elem, 2) : 0.;
9429 double norme2 = nx*nx + ny*ny + nz*nz;
9432 double i_norme = 1. / sqrt(norme2);
9433 normale_elements(elem, 0) = nx * i_norme;
9434 normale_elements(elem, 1) = ny * i_norme;
9436 normale_elements(elem, 2) = nz * i_norme;
9443 terme_src = distance_elements;
9444 tmp = distance_elements;
9445 for (iteration = 0; iteration < n_iter; iteration++)
9448 const int liste_elem_size = liste_elements.
size_array();
9449 for (i_elem = 0; i_elem < liste_elem_size; i_elem++)
9451 elem = liste_elements[i_elem];
9452 if (terme_src(elem) > distance_sommets_invalides)
9455 tmp(elem) = distance_elements(elem);
9460 double ncontrib = 0.;
9461 double somme_distances = 0.;
9463 for (k = 0; k < nb_elem_voisins; k++)
9466 const int face = elem_faces(elem, k);
9467 const int e_voisin = face_voisins(face, 0) + face_voisins(face, 1) - elem;
9470 const double distance_voisin = distance_elements(e_voisin);
9471 if (distance_voisin > distance_sommets_invalides)
9474 double nx = normale_elements(elem, 0) + normale_elements(e_voisin, 0);
9475 double ny = normale_elements(elem, 1) + normale_elements(e_voisin, 1);
9476 double nz = (dim==3)
9477 ? normale_elements(elem, 2) + normale_elements(e_voisin, 2) : 0.;
9478 double norm2 = nx*nx + ny*ny + nz*nz;
9481 double i_norm = 1./sqrt(norm2);
9487 double dx = centre_element(elem, 0) - centre_element(e_voisin, 0);
9488 double dy = centre_element(elem, 1) - centre_element(e_voisin, 1);
9489 double dz = (dim==3)
9490 ? centre_element(elem, 2) - centre_element(e_voisin, 2) : 0.;
9491 double d = nx * dx + ny * dy + nz * dz + distance_voisin;
9492 somme_distances += d;
9500 double d = somme_distances / ncontrib;
9505 distance_elements = tmp;
9508 statistics().end_count(
"Calculer_distance_interface");
9512 const ArrOfInt& compo_connexes_facettes,
9513 const int nb_compo_tot,
9514 const DoubleTab& vitesse_sommets,
9515 DoubleTab& vitesses,
9516 DoubleTab& positions)
const
9518 assert(nb_compo_tot == vitesses.
dimension(0));
9519 assert(nb_compo_tot == positions.
dimension(0));
9524 const IntTab& facettes = maillage.
facettes();
9525 const DoubleTab& sommets = maillage.
sommets();
9529 ArrOfDouble surfaces_compo(nb_compo_tot);
9536 for (
int i = 0; i < nb_facettes_tot; i++)
9540 const int compo = compo_connexes_facettes[i];
9541 const double surface = surface_facettes[i];
9542 surfaces_compo[compo] += surface;
9544 for (
int j = 0; j < dim; j++)
9547 const int s = facettes(i, j);
9548 for (
int k = 0; k < dim; k++)
9550 positions(compo, k) += surface * sommets(s, k);
9556 positions *= (1. / dim);
9560 tab_divide_any_shape(positions, s);
9566 for (
int fa7 = 0; fa7 < nb_facettes_tot; fa7++)
9571 const double si=surface_facettes[fa7];
9574 for (
int d=0; d<dim; d++)
9575 som[d]=facettes(fa7,d);
9582 const double n0=normale_facettes(fa7,0);
9583 const double n1=normale_facettes(fa7,1);
9584 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);
9588 const double n0=normale_facettes(fa7,0)*0.5;
9589 const double n1=normale_facettes(fa7,1)*0.5;
9590 const double n2=normale_facettes(fa7,2)*0.5;
9591 double s2s1[3],d_surface[3];
9593 for (
int i = 0; i < 3; i++)
9600 const int s0 = som[i];
9601 const int s1 = som[ (i+1)%3 ];
9602 const int s2 = som[ (i+2)%3 ];
9604 s2s1[0] = sommets(s1,0) - sommets(s2,0);
9605 s2s1[1] = sommets(s1,1) - sommets(s2,1);
9606 s2s1[2] = sommets(s1,2) - sommets(s2,2);
9608 d_surface[0] = s2s1[1] * n2 - s2s1[2] * n1;
9609 d_surface[1] = s2s1[2] * n0 - s2s1[0] * n2;
9610 d_surface[2] = s2s1[0] * n1 - s2s1[1] * n0;
9612 ds_dt+=d_surface[0]*vitesse_sommets(s0,0) + d_surface[1]*vitesse_sommets(s0,1) +d_surface[2]*vitesse_sommets(s0,2);
9616 const int compo = compo_connexes_facettes[fa7];
9617 for (
int d=0; d<dim; d++)
9621 for (
int d1=0; d1<dim; d1++)
9623 V+=vitesse_sommets(som[d1],d);
9624 p+=sommets(som[d1],d);
9628 vitesses(compo, d) += si * V + (p-positions(compo, d)) * ds_dt;
9636 tab_divide_any_shape(vitesses, s);
9660 DoubleVect& valeurs)
9666 const Domaine& domaine = domaine_vf.
domaine();
9667 const IntTab& face_voisins = domaine_vf.
face_voisins();
9668 const IntTab& elem_faces = domaine_vf.
elem_faces();
9671 const int nb_elem_tot = domaine.nb_elem_tot();
9675 for (
int i_face = 0; i_face < nb_faces_tot; i_face++)
9678 const int elem0 = face_voisins(i_face, 0);
9679 const int elem1 = face_voisins(i_face, 1);
9680 if (elem0 < 0 || elem1 < 0)
9684 double prod_scal = 0.;
9688 for (
int i = 0; i < dim; i++)
9690 double n0 = normale_interface(elem0, i);
9691 double n1 = normale_interface(elem1, i);
9692 prod_scal += (n0 + n1) * 0.5 * domaine_vf.
face_normales(i_face, i);
9694 if (indic(elem0) == 1. || indic(elem1) == 1.)
9695 prod_scal = - prod_scal;
9698 if (indic(elem0) != 0. && indic(elem0) != 1.)
9700 if (indic(elem1) != 0. && indic(elem1) != 1.)
9716 DoubleVect tmp_flux = flux;
9717 ArrOfInt flag(nb_faces_elem);
9718 for (
int i_elem = 0; i_elem < nb_elem_tot; i_elem++)
9722 for (
int j = 0; j < nb_faces_elem; j++)
9724 int i_face = elem_faces(i_elem, j);
9725 const int elem0 = face_voisins(i_face, 0);
9726 const int elem1 = face_voisins(i_face, 1);
9727 if (elem0 < 0 || elem1 < 0)
9730 const double signe = (elem0 == i_elem) ? 1. : -1.;
9731 const double f = flux(i_face);
9741 const double facteur = valeurs(i_elem) / somme;
9742 for (
int j = 0; j < nb_faces_elem; j++)
9744 int i_face = elem_faces(i_elem, j);
9747 tmp_flux[i_face] = tmp_flux[i_face] * facteur;
9754 for (
int i_elem = 0; i_elem < nb_elem_tot; i_elem++)
9756 for (
int j = 0; j < nb_faces_elem; j++)
9758 int i_face = elem_faces(i_elem, j);
9759 flux[i_face]=tmp_flux[i_face];
9764 for (
int i_elem = 0; i_elem < nb_elem_tot; i_elem++)
9766 double x = valeurs(i_elem);
9767 for (
int j = 0; j < nb_faces_elem; j++)
9769 int i_face = elem_faces(i_elem, j);
9770 double f = flux(i_face);
9771 double signe = (i_elem == face_voisins(i_face, 1)) ? 1. : -1.;
9774 valeurs(i_elem) = x;
9784 const DoubleVect& valeurs_euler,
9785 ArrOfDouble& valeurs_lagrange)
9787 const int nb_sommets = maillage.
nb_sommets();
9790 valeurs_lagrange = 0.;
9793 const IntTab& facettes = maillage.
facettes();
9795 const ArrOfInt& index_elem = intersections.
index_elem();
9797 const int nb_elements = valeurs_euler.
size();
9798 assert(nb_elements == index_elem.
size_array());
9802 for (element = 0; element < nb_elements; element++)
9807 const int index_premiere_intersection = index_elem[element];
9808 if (index_premiere_intersection < 0)
9816 double surface_totale = 0.;
9818 int index = index_premiere_intersection;
9823 const double surface_facette = surface_facettes[facette];
9828 const double valeur_euler = valeurs_euler[element];
9832 index = index_premiere_intersection;
9833 if (surface_totale > 0.)
9835 const double inv_surface_tot = 1. / surface_totale;
9840 const double surface_facette = surface_facettes[facette];
9845 for (i = 0; i < dim; i++)
9847 const int sommet = facettes(facette, i);
9853 const double valeur = coeff * fraction_surface * valeur_euler;
9856 valeurs_lagrange[sommet] += valeur;
9868void Transport_Interfaces_FT_Disc::compute_nb_particles_tot()
9872 ArrOfInt id_number_fa7(nb_facettes);
9873 int n = search_connex_components_local_FT(mesh, id_number_fa7);
9874 nb_particles_tot_= compute_global_connex_components_FT(mesh, id_number_fa7, n);
9885 ArrOfInt id_number_fa7(nb_facettes);
9886 int n = search_connex_components_local_FT(mesh, id_number_fa7);
9887 int nb_particles_tot=compute_global_connex_components_FT(mesh, id_number_fa7, n);
9892 Cerr <<
"WARNING, renumbering particles !" << finl;
9899 const IntTab& facettes = mesh.
facettes();
9900 const DoubleTab& sommets = mesh.
sommets();
9902 DoubleVect particles_surfaces(nb_particles_tot);
9907 for (
int i = 0; i < nb_facettes_tot; i++)
9911 const int id_number_lagrangian_vertice = id_number_fa7[i];
9912 const double surface = surface_facettes[i];
9913 particles_surfaces[id_number_lagrangian_vertice] += surface;
9918 const int s = facettes(i, j);
9921 surface * sommets(s, k);
9934 Equation_base& eqn_hydraulique = variables_internes_->refequation_vitesse_transport.valeur();
9958 Equation_base& eqn_hydraulique = variables_internes_->refequation_vitesse_transport.valeur();
9961 DoubleTab correct_particles_position(nb_particles_tot,
dimension);
9962 DoubleTab correct_particles_velocity(nb_particles_tot,
dimension);
9963 IntVect particles_correct_id_number(nb_particles_tot);
9970 for (
int wrong_id_number=0; wrong_id_number<nb_particles_tot; wrong_id_number++)
9973 int correct_id_number;
9974 if (particle_gravity_center_elem==-1) correct_id_number =-1;
9975 else correct_id_number = particles_eulerian_id_number[particle_gravity_center_elem];
9976 particles_correct_id_number(wrong_id_number)=correct_id_number;
9979 int isduplicateValue = collision_model_.valeur().check_for_duplicates(particles_correct_id_number);
9980 if (isduplicateValue == 1)
Process::exit(
"Transport_Interfaces_FT_Disc::swap_particles_position_velocity "
9981 "ERROR: duplicate value of the particles lagrangian ID number");
9984 for (
int wrong_id_number = 0; wrong_id_number < nb_particles_tot; wrong_id_number++)
9986 int good_id_number = particles_correct_id_number[wrong_id_number];
10000 const DoubleVect& mesh_volumes = domain_vf.
volumes();
10001 const IntTab& elem_faces = domain_vf.
elem_faces();
10002 const DoubleTab& phase_indicator_function = indicatrice_->valeurs();
10003 Equation_base& eqn_hydraulique = variables_internes_->refequation_vitesse_transport.valeur();
10015 for (
int elem=0; elem<domain_vf.
nb_elem(); elem++)
10018 if (phase_indicator_function(elem)==id_solid_phase)
10020 const int particle= particles_eulerian_id_number(elem);
10025 tab_velocity(elem_faces(elem,dim))+
10026 tab_velocity(elem_faces(elem,dim+
dimension))
10027 )*mesh_volumes(elem);
10030 0.5*(tab_velocity(elem_faces(elem,dim))+
10031 tab_velocity(elem_faces(elem,dim+
dimension))
10032 )*mesh_volumes(elem),2);
10040 DoubleVect s_vparticles;
10045 for (
int particle=0; particle<nb_particles_tot_; particle++)
10057void Transport_Interfaces_FT_Disc::fill_ftab_scalar(DoubleTab *ftab,
10058 const ArrOfDouble& values)
const
10061 ftab->
resize(nb_fa7, 1);
10062 for (
int fa7=0 ; fa7<nb_fa7 ; fa7++)
10063 (*ftab)(fa7,0) = (
float) values(fa7);
10066void Transport_Interfaces_FT_Disc::fill_ftab_scalar(DoubleTab *ftab,
10067 const DoubleVect& values)
const
10070 ftab->
resize(nb_fa7, 1);
10071 for (
int fa7=0 ; fa7<nb_fa7 ; fa7++)
10072 (*ftab)(fa7,0) = (
float) values(fa7);
10075void Transport_Interfaces_FT_Disc::fill_ftab_scalar(DoubleTab *ftab,
10076 const DoubleTab& values)
const
10078 const int nb_fa7 = values.
dimension(0);
10079 ftab->
resize(nb_fa7, 1);
10080 for (
int fa7=0 ; fa7<nb_fa7 ; fa7++)
10081 (*ftab)(fa7,0) = (
float) values(fa7);
10084void Transport_Interfaces_FT_Disc::fill_ftab_vector(DoubleTab *ftab,
const DoubleTab& values)
const
10087 const int nb_fa7 = values.
dimension(0);
10088 const int nb_compo = values.
dimension(1);
10089 ftab->
resize(nb_fa7, nb_compo);
10090 for (
int fa7=0 ; fa7<nb_fa7 ; fa7++)
10091 for (
int k=0 ; k<nb_compo ; k++)
10092 (*ftab)(fa7,k) = (
float) values(fa7,k);
10103 if (variables_internes_->refequation_vitesse_transport)
10108 const Equation_base& eqn_hydraulique = variables_internes_->refequation_vitesse_transport.valeur();
10127 fill_ftab_vector(ftab,vit);
10132 Cerr <<
"Error : velocity nodes post-processing : to be developped." << finl;
10141 if (variables_internes_->refequation_vitesse_transport)
10144 DoubleTab Positions,Vitesses;
10147 const Equation_base& eqn_hydraulique = variables_internes_->refequation_vitesse_transport.valeur();
10160 fill_ftab_vector(ftab,vit);
10165 Cerr <<
"Error : vitesse_repere_local nodes post-processing : to be developped." << finl;
10192 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...