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
Returns the value of bit e, then sets bit e to 1.
class Champ_Fonc_Face_VDF
DoubleTab & passe(int i=1) override
Returns field values at instant t-i.
double changer_temps(const double temps) override
Sets the time of the field.
DoubleTab & valeurs() override
Returns the array of field values at the current time.
virtual DoubleTab & futur(int i=1)
virtual DoubleTab & valeurs()=0
class Champ_base This class is the base of the fields hierarchy.
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
Returns the time of the field.
virtual DoubleTab & valeur_aux_elems(const DoubleTab &positions, const IntVect &les_polys, DoubleTab &valeurs) const
Causes an error! Must be overridden by derived classes.
class Conds_lim This class represents a vector of boundary conditions.
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
class Discretisation_base This class represents a spatial discretization scheme, which
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
Searches for the elements containing the points whose coordinates are specified.
virtual void creer_tableau_elements(Array_base &, RESIZE_OPTIONS opt=RESIZE_OPTIONS::COPY_INIT) const
Creates a parallel array of values at elements.
DoubleTab_t & les_sommets()
int nb_faces_elem(int=0) const
Returns the number of faces of type i of the geometric elements that make up the domain.
virtual void creer_tableau_sommets(Array_base &, RESIZE_OPTIONS opt=RESIZE_OPTIONS::COPY_INIT) const
Creates an array with one "row" per mesh vertex.
int orientation(int) const override
inline DoubleVect& Domaine_VDF::porosite_face() {
double distance_face(int, int, int k) const
Returns the distance between the centers of two same-orientation faces.
virtual const DoubleVect & face_surfaces() const
DoubleVect & volumes_entrelaces()
int nb_faces_tot() const
Returns the total number of 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()
Fills in the face_virt_pe_num_ array (see comment in Domaine_VF.h).
int elem_faces(int i, int j) const
Returns the index of the i-th face of element num_elem; the face numbering convention is.
double xp(int num_elem, int k) const
int face_voisins(int num_face, int i) const
Returns the neighbouring element of num_face in direction i.
int moments_a_imprimer() const
const Nom & le_nom() const override
Returns the name of the Objet_U. Virtual method to override: returns "neant" in this implementation.
class Domaine_dis_base This class is the base of the hierarchy of discretized domains.
virtual IntTab & face_voisins()
const Domaine & domaine() const
Reading from a file of objects written in binary format.
File for reading. This class is to the C++ ifstream class what the Entree class is to the.
static Nom & get_Output()
Returns the write mode in use (so it can be modified).
static int is_ecriture_special(int &special, int &a_faire)
Indicates whether the special format was requested in active writing by xyz save.
static int is_lecture_special()
Indicates whether the special format was requested in active reading by xyz restart.
File for reading. This class is to the C++ ifstream class what the Entree class is to the.
Class defining operators and methods for all reading operation in an input flow (file,...
class Equation_base The role of an equation is the calculation of one or more fields....
virtual void set_param(Param &titi) const override
const Nom & le_nom() const override
Returns the name of the equation.
int reprendre(Entree &) override
We resume the unknown from an input stream.
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 &)
Associates with the Problem passed as parameter.
const Discretisation_base & discretisation() const
Returns the discretization associated with the equation.
Solveur_Masse_base & solv_masse()
Returns the mass solver associated with the 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
We save the unknown, then the source terms to an output stream.
virtual const Champ_Inc_base & derivee_en_temps() const
Probleme_base & probleme()
Returns the problem associated with the equation.
int lire_motcle_non_standard(const Motcle &, Entree &) override
Reads non-simple-type parameters of an Objet_U from an input stream.
int limpr() const
Ask the time scheme if an output is needed.
Schema_Temps_base & schema_temps()
Returns the time scheme associated with the equation.
Champs_compris champs_compris_
virtual const Operateur & operateur(int) const =0
Domaine_dis_base & domaine_dis()
Returns the discretized domain associated with the equation.
const int & get_id_fluid_phase() const
: Post-processing class for Eulerian fields in LATA format.
int ecrire_entete(const double temps_courant, const int reprise, const int est_le_premier_post) override
Opens the master file in ERASE mode and writes the lata file header (on the master processor only).
int ecrire_domaine(const Domaine &domaine, const int est_le_premier_post) override
See Format_Post_base::ecrire_domaine. Writing a domain within a time step is accepted,...
int ecrire_temps(const double temps) override
Starts writing a new time step. For the LATA format specifically:
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 &)
See Interprete_bloc::objet_global(). BM: the Interprete class is not the best place for this.
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
Returns the joint information for the requested type (read-only).
const IntTab_t & renum_items_communs() const
Returns the renum_items_communs_ array (read-only). See renum_items_communs_.
This class implements the operators and virtual methods of the EFichier class as follows: The file to...
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)
Matrice_Morse class - Represents a (sparse) matrix M, not necessarily square,.
Milieu_base This class is the base of the (physical) medium hierarchy.
A character string (Nom) in uppercase.
An array of Motcle objects.
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)
Navier_Stokes_std This class carries the terms of the momentum equation.
const Champ_Inc_base & inconnue() const override
Returns the velocity (unknown field of the equation) (const version).
class Nom: a character string for naming TRUST objects.
virtual int finit_par(const char *const n) const
virtual int find(const char *const n) const
const std::string & getString() const
An array of character strings (VECT(Nom)).
Base class for TRUST objects (Objet_U).
const Nom & que_suis_je() const
Returns the string identifying the class.
virtual Entree & readOn(Entree &)
Reads an Objet_U from an input stream. Virtual method to override.
static double precision_geom
static const Nom & nom_du_cas()
Returns a constant reference to the case name. This method is static.
virtual int reprendre(Entree &)
Restores an Objet_U from an input stream. Virtual method to override.
virtual Sortie & printOn(Sortie &) const
Writes the object to an output stream. Virtual method to override.
class Operateur Generic class of the operator hierarchy.
double calculer_pas_de_temps() const
Computes the next time step.
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 of data for the Parser class.
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
Returns the name of the Objet_U. Virtual method to override: returns "neant" in this implementation.
class Probleme_base It is a Probleme_U that is not a coupling.
const Nom & checkpoint_filename() const
bool is_dilatable() const
bool & reprise_effectuee()
const Schema_Temps_base & schema_temps() const
Returns the time scheme associated with the problem.
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)
Returns a static Sortie object used as an event journal.
static void abort()
Abort routine for TRUST on a fatal error.
static double mp_sum(double)
Computes the sum of x over all processors in the current group.
static void barrier()
Synchronizes all processors in the current group (waits until all processors have reached the barrier...
static int me()
Returns the rank of the local processor in the current communication group. See Comm_Group::rank() an...
static void exit(int exit_code=-1)
Exit routine for TRUST within a Kokkos region.
static int je_suis_maitre()
Returns 1 if on the master processor of the current group (i.e. me() == 0), 0 otherwise.
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)
SFichier is to the C++ ofstream class what Sortie is to the C++ ostream class.
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)
Method used by interpreters that modify the domain (sequential), destroys the descriptors of vertices...
int nb_impr() const
Returns the number of outputs performed.
double temps_courant() const
Returns the current time.
int precision_impr() const
void imprimer_temps_courant(SFichier &) const
double pas_de_temps() const
Returns the current time step (delta_t).
virtual int nb_valeurs_temporelles() const =0
void update_critere_statio(const DoubleTab &tab_critere, Equation_base &equation)
Updates stationnaire_atteint_ and residu_ (criterion: residu_ < seuil_statio_).
Solveur_Masse_base Represents the mass matrix of an equation.
virtual DoubleTab & appliquer(DoubleTab &) const
Returns appliquer_impl(x/temporal_coefficient) if a temporal coefficient is set, otherwise returns ap...
virtual int ouvrir(const char *name, IOS_OPEN_MODE mode=ios::out)
void precision(int pre) override
void setf(IOS_FORMAT code) override
Base class for output streams.
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")
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
Restores an Objet_U from an input stream. Virtual method to override.
Maillage_FT_Disc maillage_interface
int sauvegarder(Sortie &os) const override
Saves an Objet_U to an output stream. Virtual method to override.
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
Everything that does not depend on other possible problems.
void fill_ftab_vertices_curvature(DoubleTab *ftab, const DoubleTab &dummytab) const
int sauvegarder(Sortie &) const override
We save the unknown, then the source terms to an output stream.
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
Prints the equation operators to an output stream, unconditionally.
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
The value of the unknown at the time step has been calculated.
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
Calculation of the next time step.
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
Associates with the Problem passed as parameter.
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
Reads non-simple-type parameters of an Objet_U from an input stream.
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
We resume the unknown from an input stream.
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_
Transport_Interfaces_base This class is the base class for interface transport equations.
bool reinjection_tcl() const
const double & t_injection() const
bool ready_inject_tcl() const
const double & thetaC_tcl() const
const double & Rc_inject() const
YAML_data class: collection of all needed information for data to save/restore in order to write the ...