16#include <Connex_components.h>
17#include <IJK_Field_vector.h>
18#include <Connex_components_FT.h>
19#include <EcrFicPartageBin.h>
20#include <Probleme_FTD_IJK_base.h>
21#include <Probleme_FTD_IJK_cut_cell.h>
22#include <Cut_cell_tools.h>
23#include <IJK_Interfaces.h>
24#include <IJK_Lata_writer.h>
25#include <IJK_Navier_Stokes_tools.h>
26#include <Cut_cell_tools.h>
28#include <LecFicDiffuse_JDD.h>
29#include <Linear_algebra_tools.h>
30#include <Octree_Double.h>
31#include <Ouvrir_fichier.h>
34#include <Domaine_VF.h>
36#include <communications.h>
39#include <Perf_counters.h>
41#include <Transport_Interfaces_FT_Disc.h>
42#include <Navier_Stokes_FTD_IJK.h>
51#define CLASSIC_METHOD 1
69IJK_Interfaces::IJK_Interfaces()
78 int size_fixed_arr[3] = {0,0,0};
79 for (
int c=0; c<3; c++)
80 size_fixed_arr[c] = fixed_arr[c].size_array();
81 assert(size_fixed_arr[0] == size_fixed_arr[1]);
82 assert(size_fixed_arr[0] == size_fixed_arr[2]);
84 tab.
resize(size_fixed_arr[0],3);
85 for (
int c=0; c<3; c++)
86 for (
int l=0; l<size_fixed_arr[0]; l++)
87 tab(l,c) = fixed_arr[c](l);
108 Nom fdsom = prefix +
Nom(
"SOMMETS");
109 Nom fdelem = prefix +
Nom(
"ELEMENTS");
110 Nom fdpelocal = prefix +
Nom(
"FACETTE_PE_LOCAL");
111 Nom fdpeowner = prefix +
Nom(
"FACETTE_PE_OWNER");
112 Nom fd_sommet_reel = prefix +
Nom(
"INDEX_SOMMET_REEL");
117 int nsomtot = nb_som;
118 int nelemtot = nb_elem;
125 FloatTab tmp(nb_som, 3);
126 const DoubleTab& sommets = mesh.
sommets();
128 for (
int i = 0; i < nb_som; i++)
129 for (
int j = 0; j < 3; j++)
130 tmp(i, j) = (float)sommets(i, j);
144 IntTab tmp_facettes_renum(nb_elem, 3);
145 const IntTab& facettes = mesh.
facettes();
146 for (
int i = 0; i < nb_elem; i++)
147 for (
int j = 0; j < 3; j++)
148 tmp_facettes_renum(i, j) = facettes(i, j) + offset_sommet;
156 ArrOfInt tmp2(nb_elem);
158 file.
put(tmp2.
addr(), nb_elem, 1);
165 file.
put(tmp2.
addr(), nb_elem, 1);
172 ArrOfInt indice_global(nb_som);
175 ArrOfInt offset_received(np);
176 ArrOfInt offset_to_send(np);
179 offset_to_send = offset_sommet;
181 envoyer_all_to_all(offset_to_send, offset_received);
188 for (
int i = 0; i < nb_som; i++)
190 const int indice_local_sur_proc_proprietaire = sommet_num_owner[i];
191 const int proc_proprietaire = sommet_PE_owner[i];
192 const int offset_proc_proprietaire = offset_received[proc_proprietaire];
193 indice_global[i] = indice_local_sur_proc_proprietaire + offset_proc_proprietaire;
195 file.
ouvrir(fd_sommet_reel);
201 const Nom format =
"INT32";
205 master_file.
ouvrir(filename, ios::app);
206 master_file <<
"Geom " << meshname <<
" type_elem=TRIANGLE_3D" << finl;
208 master_file <<
"Champ SOMMETS " << basename(fdsom) <<
" geometrie=" << meshname;
209 master_file <<
" size=" << nsomtot <<
" composantes=3" << finl;
211 master_file <<
"Champ ELEMENTS " << basename(fdelem) <<
" geometrie=" << meshname;
212 master_file <<
" size=" << nelemtot <<
" composantes=3 format=" << format << finl;
214 master_file <<
"Champ FACETTE_PE_LOCAL " << basename(fdpelocal) <<
" geometrie=" << meshname;
215 master_file <<
" size=" << nelemtot <<
" composantes=1 format=" << format <<
" localisation=ELEM" << finl;
217 master_file <<
"Champ FACETTE_PE_OWNER " << basename(fdpeowner) <<
" geometrie=" << meshname;
218 master_file <<
" size=" << nelemtot <<
" composantes=1 format=" << format <<
" localisation=ELEM" << finl;
220 master_file <<
"Champ INDEX_SOMMET_REEL " << basename(fd_sommet_reel) <<
" geometrie=" << meshname;
221 master_file <<
" size=" << nsomtot <<
" composantes=1 format=" << format <<
" localisation=SOM" << finl;
226void runge_kutta3_update(
const DoubleTab& dvi, DoubleTab& G, DoubleTab& l,
const int step,
const double dt_tot,
229 const double coeff_a[3] = {0., -5. / 9., -153. / 128.};
231 const double coeff_Gk[3] = {1., 4. / 9., 15. / 32.};
233 const double facteurG = coeff_a[step];
234 const double intermediate_dt = compute_fractionnal_timestep_rk3(dt_tot, step);
235 const double delta_t_divided_by_Gk = intermediate_dt / coeff_Gk[step];
247 for (
int isom = 0; isom < nbsom; isom++)
249 for (
int dir = 0; dir < 3; dir++)
253 G(isom, dir) = 111. / intermediate_dt;
254 l(isom, dir) = 111. / intermediate_dt;
256 double x = dvi(isom, dir);
258 l(isom, dir) += x * delta_t_divided_by_Gk;
264 for (
int isom = 0; isom < nbsom; isom++)
266 for (
int dir = 0; dir < 3; dir++)
270 G(isom, dir) = 111. / intermediate_dt;
271 l(isom, dir) = 111. / intermediate_dt;
273 double x = G(isom, dir) * facteurG + dvi(isom, dir);
275 l(isom, dir) += x * delta_t_divided_by_Gk;
281 for (
int isom = 0; isom < nbsom; isom++)
283 for (
int dir = 0; dir < 3; dir++)
287 G(isom, dir) = 111. / intermediate_dt;
288 l(isom, dir) = 111. / intermediate_dt;
290 double x = G(isom, dir) * facteurG + dvi(isom, dir);
291 l(isom, dir) += x * delta_t_divided_by_Gk;
296 Cerr <<
"Error in runge_kutta_update: wrong step" << finl;
310 os <<
" terme_gravite rho_g \n";
312 os <<
" terme_gravite grad_i \n";
315 os <<
" correction_gradient_potentiel \n";
327 os <<
" active_repulsion_paroi \n";
333 Cerr <<
"IJK_Interfaces::printOn -- couleurs des bulles stockees " << finl;
334 os <<
" follow_colors " <<
"\n"
339 Cerr <<
"IJK_Interfaces::printOn -- Groupes d'appartenance stockes " << finl;
345 if (
parcours_.get_correction_parcours_thomas())
346 os <<
" parcours_interface { correction_parcours_thomas } " <<
"\n";
347 if (
parcours_.get_parcours_sans_tolerance())
348 os <<
" parcours_interface { parcours_sans_tolerance } " <<
"\n";
351 os <<
"use_barycentres_velocity" <<
"\n";
353 os <<
"use_barycentres_velocity" <<
"\n";
355 double max_force_compo = 0.;
358 if (max_force_compo > 1.e-16)
361 Cerr <<
"IJK_Interfaces::printOn -- Storing mean into sauv for future "
362 "restart. (max_force_compo= "
363 << max_force_compo <<
" )." << finl;
407 param.ajouter_flag(
"frozen", &
frozen_);
410 param.ajouter(
"parcours_interface",&
parcours_);
427 param.lire_avec_accolades(is);
429 Cout <<
"IJK_Interfaces::readOn : Option gravite : " <<
terme_gravite_
441 Cout <<
"IJK_Interfaces::readOn : il y a : " <<
nb_groups_
442 <<
" classes/groupes de bulles suivis " << finl;
448 <<
" bulles reelles reprises avec une position de reference. " << finl;
452 Cout <<
"IJK_Interfaces::readOn : " <<
mean_force_.dimension(0)
453 <<
" bulles reelles reprises avec une force_moyenne imposee. " << finl;
460 if (un_mot==
"parcours_interface")
467 Cerr <<
"Unknown Keyword " << un_mot <<
" in IJK_Interfaces::lire_motcle_non_standard" << finl;
486 if (!
parcours_.get_correction_parcours_thomas())
488 Cerr <<
"Warning: The cut-cell method overrides default or user parameters to force the use of correction_parcours_thomas in parcours_interface." << finl;
489 parcours_.set_correction_parcours_thomas();
492 if (!
parcours_.get_parcours_sans_tolerance())
494 Cerr <<
"Warning: The cut-cell method overrides default or user parameters to force the use of parcours_sans_tolerance in parcours_interface." << finl;
511 double min_indicatrice = 1.;
514 for (
int n = 0; n < cut_cell_disc.
get_n_loc(); n++)
516 Int3 ijk = cut_cell_disc.
get_ijk(n);
524 if (indicatrice > 0.)
525 min_indicatrice = std::min(min_indicatrice, indicatrice);
527 if (indicatrice <= .5 && indicatrice > .4)
529 else if (indicatrice <= .4 && indicatrice > .3)
531 else if (indicatrice <= .3 && indicatrice > .2)
533 else if (indicatrice <= .2 && indicatrice > .1)
535 else if (indicatrice <= .1 && indicatrice > .05)
537 else if (indicatrice <= .05 && indicatrice > .01)
539 else if (indicatrice <= .01 && indicatrice > .0)
543 assert(indicatrice == 0.);
548 mp_sum_for_each(count_total, count_40pct, count_30pct, count_20pct, count_10pct);
562 Cerr <<
"Bilan de l'indicatrice:" << finl;
563 Cerr <<
" Valeur minimale: " << min_indicatrice << finl;
564 Cerr <<
" Compte sur les cellules de la structure diphasique:" << finl;
565 Cerr <<
" 40-50% " << count_40pct << finl;
566 Cerr <<
" 30-40% " << count_30pct << finl;
567 Cerr <<
" 20-30% " << count_20pct << finl;
568 Cerr <<
" 10-20% " << count_10pct << finl;
569 Cerr <<
" 5-10% " << count_5pct << finl;
570 Cerr <<
" 1-5% " << count_1pct << finl;
571 Cerr <<
" 0-1% " << count_0pct << finl;
572 Cerr <<
" pure " << count_pure << finl;
573 Cerr <<
" ------ " << finl;
574 Cerr <<
" total " << count_total << finl;
582 int iteration_solver_surface_efficace_face = 0;
584 if (type_surface_efficace_face == TYPE_SURFACE_EFFICACE_FACE::CONSERVATION_VOLUME_ITERATIF)
591 iteration_solver_surface_efficace_face,
603 iteration_solver_surface_efficace_face,
618 (type_surface_efficace_interface == TYPE_SURFACE_EFFICACE_INTERFACE::EXPLICITE),
638 if (type_surface_efficace_interface == TYPE_SURFACE_EFFICACE_INTERFACE::CONSERVATION_VOLUME)
665 (type_surface_efficace_face == TYPE_SURFACE_EFFICACE_FACE::EXPLICITE),
675 (type_surface_efficace_interface == TYPE_SURFACE_EFFICACE_INTERFACE::EXPLICITE),
689 const IJK_Field_vector3_double& velocity_ft = ref_ijk_ft_->eq_ns().get_velocity_ft();
691 const DoubleTab& sommets = mesh.
sommets();
693 ArrOfDouble vinterp_component(nbsom);
695 for (
int direction = 0; direction < 3; direction++)
701 ijk_interpolate_skip_unknown_points(velocity_ft[direction], sommets, vinterp_component,
703 for (
int i = 0; i < nbsom; i++)
704 vinterp_(i, direction) = vinterp_component[i];
726 if (nom==
"CONCENTRATION_INTERFACE")
731 if (nom==
"GRADX_CONCENTRATION_INTERFACE")
733 IJK_Field_double& dconc_interf_x =
scalar_post_fields_.at(
"GRADX_CONCENTRATION_INTERFACE");
736 if (nom==
"GRADY_CONCENTRATION_INTERFACE")
738 IJK_Field_double& dconc_interf_y =
scalar_post_fields_.at(
"GRADY_CONCENTRATION_INTERFACE");
741 if (nom==
"GRADZ_CONCENTRATION_INTERFACE")
743 IJK_Field_double& dconc_interf_z =
scalar_post_fields_.at(
"GRADZ_CONCENTRATION_INTERFACE");
746 if (nom==
"LAPLACIAN_CONCENTRATION_INTERFACE")
748 IJK_Field_double& lapla_interf =
scalar_post_fields_.at(
"LAPLACIAN_CONCENTRATION_INTERFACE");
756 if (nom==
"GRADX_SIGMA")
759 dsigma_x.
data() =
maillage_ft_ijk_.Surfactant_facettes().get_interfacial_source_term_sommets(0);
761 if (nom==
"GRADY_SIGMA")
764 dsigma_y.
data() =
maillage_ft_ijk_.Surfactant_facettes().get_interfacial_source_term_sommets(1);
766 if (nom==
"GRADZ_SIGMA")
769 dsigma_z.
data() =
maillage_ft_ijk_.Surfactant_facettes().get_interfacial_source_term_sommets(2);
771 if (nom==
"DISTANCE_AUTRES_INTERFACES")
773 DoubleTab vr_to_closer;
779 return champs_compris_.get_champ(nom);
781 Cerr <<
"ERROR in IJK_Interfaces::get_IJK_field_vector : " << finl;
782 Cerr <<
"Requested field '" << nom <<
"' is not recognized by IJK_Interfaces::get_IJK_field_vector()." << finl;
789 return champs_compris_.get_champ_vectoriel(nom);
791 Cerr <<
"ERROR in IJK_Interfaces::get_IJK_field_vector : " << finl;
792 Cerr <<
"Requested field '" << nom <<
"' is not recognized by IJK_Interfaces::get_IJK_field_vector()." << finl;
798 std::vector<FieldInfo_t> c =
801 {
"INDICATRICE", Entity::ELEMENT, Nature_du_champ::scalaire,
false },
802 {
"INDICATRICE_FT", Entity::ELEMENT, Nature_du_champ::scalaire,
false },
803 {
"OLD_INDICATRICE", Entity::ELEMENT, Nature_du_champ::scalaire,
false },
804 {
"OLD_INDICATRICE_FT", Entity::ELEMENT, Nature_du_champ::scalaire,
false },
805 {
"REPULSION_FT", Entity::ELEMENT, Nature_du_champ::scalaire,
false },
806 {
"CUT_FIELDS_BARY_L", Entity::ELEMENT, Nature_du_champ::vectoriel,
false },
808 {
"COURBURE", Entity::NODE, Nature_du_champ::scalaire,
true },
809 {
"CONCENTRATION_INTERFACE", Entity::ELEMENT, Nature_du_champ::scalaire,
true },
810 {
"GRADX_CONCENTRATION_INTERFACE", Entity::NODE, Nature_du_champ::scalaire,
true },
811 {
"GRADY_CONCENTRATION_INTERFACE", Entity::NODE, Nature_du_champ::scalaire,
true },
812 {
"GRADZ_CONCENTRATION_INTERFACE", Entity::NODE, Nature_du_champ::scalaire,
true },
813 {
"SIGMA", Entity::NODE, Nature_du_champ::scalaire,
true },
814 {
"GRADX_SIGMA", Entity::NODE, Nature_du_champ::scalaire,
true },
815 {
"GRADY_SIGMA", Entity::NODE, Nature_du_champ::scalaire,
true },
816 {
"GRADZ_SIGMA", Entity::NODE, Nature_du_champ::scalaire,
true },
817 {
"LAPLACIAN_CONCENTRATION_INTERFACE", Entity::ELEMENT, Nature_du_champ::scalaire,
true },
818 {
"DISTANCE_AUTRES_INTERFACES", Entity::NODE, Nature_du_champ::scalaire,
true }
821 chps.insert(chps.end(), c.begin(), c.end());
826 for (
const auto& n : champs_compris_.liste_noms_compris())
828 for (
const auto& n : champs_compris_.liste_noms_compris_vectoriel())
834 Cerr <<
"Reset bubble rising velocity calculations" << finl;
836 Cerr <<
"Compute compo_connex from bounding box" << finl;
838 Cerr <<
"Compute compo_connex from interface compo in mixed cells" << finl;
840 Cerr <<
"Compute rising velocity from compo connex (barycentre calc)" << finl;
847 const int thermal_probes_ghost_cells,
848 const bool compute_vint,
849 const bool is_switch)
851 Cerr <<
"Entree dans IJK_Interfaces::initialize" << finl;
854 ref_domaine_ = domaine_FT;
864 const int nb_ghost_cells = std::max(thermal_probes_ghost_cells, (
int) 4);
914 allocate_cell_vector(groups_indicatrice_ft_test_, domaine_FT, 1);
915 allocate_cell_vector(groups_indicatrice_ft_test_, domaine_FT, 1);
919 for (
int i = 0; i < max_authorized_nb_of_components_; i++)
927 for (
int dir = 0; dir < 3; dir++)
935 for (
int dir = 0; dir < 3; dir++)
940 for (
int dir = 0; dir < 3; dir++)
942 const int idx = i * 3 + dir;
967 for (
int bary_compo = 0; bary_compo < 3; bary_compo++)
980 for (
int face_dir = 0; face_dir < 3; face_dir++)
993 for (
int face_dir = 0; face_dir < 3; face_dir++)
1001 for (
int face_dir = 0; face_dir < 3; face_dir++)
1003 for (
int bary_compo = 0; bary_compo < 2; bary_compo++)
1012 for (
int face_dir = 0; face_dir < 3; face_dir++)
1014 for (
int bary_compo = 0; bary_compo < 2; bary_compo++)
1038 for (
int d = 0; d < 3; d++)
1052 Cut_field_vector3_double& cut_field_deformation_velocity =
static_cast<Cut_field_vector3_double&
>(
deformation_velocity_);
1053 allocate_velocity_ephemere(*ref_ijk_ft_->get_cut_cell_disc(), cut_field_deformation_velocity, domaine_NS, 2);
1076 for (
int direction = 0; direction < 3; direction++)
1078 const double ori = domaine_NS.
get_origin(direction);
1088 const double oriFT = domaine_FT.
get_origin(direction);
1096 double duplicate_stencil_width ;
1099 duplicate_stencil_width =
1103 duplicate_stencil_width =
1183 DoubleTab centre_gravite;
1184 ArrOfDouble volume_bulles;
1202 DoubleTab vr_to_closer;
1222 Cerr <<
"Error in option combination: invalid choice of optimized methods "
1223 <<
"for both the color_function in the forces computation "
1224 <<
"(parser_=0) and the indicator function calculation "
1225 <<
"(recompute_indicator_=0)"
1227 Cerr <<
"Maybe you are not using the forces computation?" << finl;
1236 const Domaine& domaine = domaine_vf.
domaine();
1251 courb.nommer(
"COURBURE");
1252 champs_compris_.ajoute_champ(courb);
1256 conc_interf.nommer(
"CONCENTRATION_INTERFACE");
1257 champs_compris_.ajoute_champ(conc_interf);
1261 dconc_interf_x.nommer(
"GRADX_CONCENTRATION_INTERFACE");
1262 champs_compris_.ajoute_champ(dconc_interf_x);
1266 dconc_interf_y.nommer(
"GRADY_CONCENTRATION_INTERFACE");
1267 champs_compris_.ajoute_champ(dconc_interf_y);
1271 dconc_interf_z.nommer(
"GRADZ_CONCENTRATION_INTERFACE");
1272 champs_compris_.ajoute_champ(dconc_interf_z);
1276 sigma.nommer(
"SIGMA");
1277 champs_compris_.ajoute_champ(sigma);
1281 dsigma_x.nommer(
"GRADX_SIGMA");
1282 champs_compris_.ajoute_champ(dsigma_x);
1286 dsigma_y.nommer(
"GRADY_SIGMA");
1287 champs_compris_.ajoute_champ(dsigma_y);
1291 dsigma_z.nommer(
"GRADZ_SIGMA");
1292 champs_compris_.ajoute_champ(dsigma_z);
1296 lapla_interf.nommer(
"LAPLACIAN_CONCENTRATION_INTERFACE");
1297 champs_compris_.ajoute_champ(lapla_interf);
1301 d.nommer(
"DISTANCE_AUTRES_INTERFACES");
1302 champs_compris_.ajoute_champ(d);
1307 Cerr <<
"IJK_Interfaces::milieu not coded ! access it from NS. " << finl;
1313 Cerr <<
"IJK_Interfaces::milieu not coded ! access it from NS. " << finl;
1321 Cerr <<
"Error for the method IJK_Interfaces::associer_pb_base\n";
1322 Cerr <<
" IJK_Interfaces equation must be associated to\n";
1323 Cerr <<
" a Probleme_FTD_IJK_base problem type\n";
1340 ref_ijk_ft_switch_ = ijk_ft_switch;
1349 liste.add(
"INTERFACES");
1350 liste.add(
"COMPO_CONNEXE");
1351 liste.add(
"DISTANCE_AUTRES_INTERFACES");
1353 liste.add(
"COLOR_Y");
1358 liste.add(
"RK3_STORE_VI");
1365 const char *lata_name,
1366 const int lata_step)
const
1371 if (liste_post_instantanes.
contient_(
"INTERFACES"))
1375 if (liste_post_instantanes.
contient_(
"COLOR_Y"))
1379 Cerr <<
"On demande le post-traitement du champ COLOR_Y alors que le flag follow_colors "
1380 <<
"n'a pas ete active... Ce champ est donc inconnu. Post-traitement refuse. " << finl;
1385 n++, dumplata_ft_field(lata_name,
"INTERFACES",
"COLOR_Y",
"ELEM", color, lata_step);
1387 if (liste_post_instantanes.
contient_(
"VINTERP"))
1389 n++, dumplata_ft_field(lata_name,
"INTERFACES",
"VINTERP",
"SOM",
RK3_G_store_vi_, lata_step);
1391 if (liste_post_instantanes.
contient_(
"VI"))
1400 for (
int dir= 0; dir < 3; dir++)
1404 Cerr <<
"Posttraitement du champ VI vide : on stocke 0." << finl;
1408 for (
int i = 0; i < nbsom; i++)
1412 n++, dumplata_ft_field(lata_name,
"INTERFACES",
"VI_X",
"SOM", vi, lata_step);
1414 n++, dumplata_ft_field(lata_name,
"INTERFACES",
"VI_Y",
"SOM", vi, lata_step);
1416 n++, dumplata_ft_field(lata_name,
"INTERFACES",
"VI_Z",
"SOM", vi, lata_step);
1421 if (liste_post_instantanes.
contient_(
"RK3_STORE_VI"))
1423 ArrOfDoubleFT store_vi_compo;
1430 for (
int dir= 0; dir < 3; dir++)
1435 Cerr <<
"Posttraitement du champ RK3_STORE_VI vide : on stocke 0." << finl;
1436 store_vi_compo = 0.;
1440 for (
int i = 0; i < nbsom; i++)
1443 store_vi_compo[i] = val;
1447 n++, dumplata_ft_field(lata_name,
"INTERFACES",
"RK3_STORE_VI_X",
"SOM", store_vi_compo, lata_step);
1449 n++, dumplata_ft_field(lata_name,
"INTERFACES",
"RK3_STORE_VI_Y",
"SOM", store_vi_compo, lata_step);
1451 n++, dumplata_ft_field(lata_name,
"INTERFACES",
"RK3_STORE_VI_Z",
"SOM", store_vi_compo, lata_step);
1472 const Domaine_IJK& geom_FT = ref_domaine_.valeur();
1473 for (
int direction = 0; direction < 3; direction++)
1479 const double oriFT = geom_FT.
get_origin(direction);
1497 DoubleTab bounding_box;
1503 ArrOfInt masque_delete_pour_compo;
1513 int nb_bulles_reelles_deleted = 0;
1514 ArrOfInt old_to_new_compo;
1519 for (
int i = 0; i < masque_delete_pour_compo.
size_array(); i++)
1521 if (masque_delete_pour_compo[i] != 0)
1524 nb_bulles_reelles_futur -= 1;
1525 nb_bulles_reelles_deleted += 1;
1526 old_to_new_compo[i] = -1;
1530 old_to_new_compo[i] = count;
1534 nb_bulles_reelles_futur = count;
1537 envoyer_broadcast(old_to_new_compo, 0);
1538 envoyer_broadcast(nb_bulles_reelles_futur, 0);
1539 envoyer_broadcast(flag, 0);
1540 envoyer_broadcast(nb_bulles_reelles_deleted, 0);
1544 Cerr << nb_bulles_reelles_deleted <<
" marked for deletion out of " <<
nb_bulles_reelles_
1545 <<
" so that only " << nb_bulles_reelles_futur <<
" will remain." << finl;
1549 for (
int fa7 = 0; fa7 < mesh.
nb_facettes(); fa7++)
1551 icompo = compo_connexe_facettes[fa7];
1552 assert(icompo >= 0);
1563 <<
" to " << nb_bulles_reelles_futur << finl;
1571 int inew = old_to_new_compo[icompo];
1572 assert(inew <= icompo);
1589 Cerr <<
"The table of bubbles groups (compo_to_group_), "
1590 <<
" (as well as possibly through_yminus_ and positions_reference_ "
1591 <<
"if needed)" <<
" have been updated in accordance." << finl;
1624 os <<
"# Impression des couleurs de bulles" << finl;
1625 os <<
"# colonne 1 : temps" << finl;
1626 os <<
"# colonne K : Couleur de la bulle K-2" << finl;
1628 os << current_time <<
" ";
1629 for (
int j = 0; j < n; j++)
1640 for (
int i = 0; i < n; i++)
1646 const int compo = compo_facettes[i];
1669 for (
int i = 0; i < n; i++)
1677 Cerr <<
"Exception dans IJK_Interfaces::get_ghost_number_from_compo(). "
1678 <<
" Compo demandee " << compo <<
" introuvable...";
1689 surfaces.
resize_array(nbulles_reelles + nbulles_ghost, RESIZE_OPTIONS::NOCOPY_NOINIT);
1693 for (
int i = 0; i < n; i++)
1697 int compo = compo_facettes[i];
1704 compo = nbulles_reelles - 1 - idx_ghost;
1707 const double s = surfaces_facettes[i];
1708 surfaces[compo] += s;
1719 ArrOfDouble& out)
const
1727 for (
int i = 0; i < n; i++)
1731 int compo = compo_facettes[i];
1736 const double s = in[i];
1741 for (
int c = 0; c < nbulles; c++)
1742 out[c] /= surfaces[c];
1752 Cerr <<
"interfacial surface and normale are sure to be up-to-date" << finl;
1770 bool tmp_readen_flag =
false;
1771 Nom suffix_tmp =
".old";
1774 if (tmp_readen_flag)
1778 suffix_tmp =
".new";
1781 if (tmp_readen_flag)
1788 if (tmp_readen_flag)
1800 ArrOfDouble& bubbles_rising_vel_mag)
1802 bool is_readen =
false;
1805 bubbles_rising_vel_mag.
reset();
1806 for (
int c=0; c<3; c++)
1808 bubbles_rising_dir[c].reset();
1809 bubbles_rising_vel[c].reset();
1811 if (bubbles_bary_computed)
1813 int line_counter = 0;
1817 const Nom interf_dir = dirname(interf_name);
1819 const Nom suffix_bary =
".sauv.barycentres";
1820 const Nom file_folder = interf_dir + case_name + suffix_bary +
".vel";
1823 read_tmp.open(file_folder);
1824 const bool read_file = read_tmp ? true :
false;
1828 EFichier fic_bary_vel(file_folder);
1829 ifstream& ifstream_bary_vel = fic_bary_vel.
get_ifstream();
1830 const char delimiter =
' ';
1831 Cerr <<
"Read coordinates of bubbles barycentres from: " << file_folder << finl;
1832 while (std::getline(ifstream_bary_vel, line))
1834 std::stringstream ssline(line);
1835 Cerr <<
"Line number: " << line_counter << finl;
1836 while (std::getline(ssline, line, delimiter))
1840 Cerr <<
"Param: " << line << finl;
1848 dir = var_index - 1;
1849 bubbles_rising_vel[dir].append_array(std::stod(line));
1854 dir = var_index - 4;
1855 bubbles_rising_dir[dir].append_array(std::stod(line));
1869 fic_bary_vel.
close();
1878 bool is_readen =
false;
1881 for (
int c=0; c<3; c++)
1883 bubbles_bary[c].reset();
1885 if (bubbles_bary_computed)
1887 int line_counter = 0;
1891 const Nom interf_dir = dirname(interf_name);
1893 const Nom suffix_bary =
".sauv.barycentres";
1894 const Nom file_folder = interf_dir + case_name + suffix_bary + suffix;
1896 read_tmp.open(file_folder);
1897 const bool read_file = read_tmp ? true :
false;
1903 const char delimiter =
' ';
1904 Cerr <<
"Read coordinates of bubbles barycentres from: " << file_folder << finl;
1905 while (std::getline(ifstream_bary_old, line))
1907 std::stringstream ssline(line);
1908 Cerr <<
"Line number: " << line_counter << finl;
1909 while (std::getline(ssline, line, delimiter))
1913 Cerr <<
"Param: " << line << finl;
1921 dir = var_index - 1;
1922 bubbles_bary[dir].append_array(std::stod(line));
1942 const Nom end_space =
" ";
1943 const Nom escape =
"\n";
1948 const Nom interf_dir = dirname(interf_name);
1949 const Nom suffix =
".sauv.barycentres";
1950 const int reset = 1;
1951 const Nom bary_header_old =
"ibubble bary_old_x bary_old_y bary_old_z";
1952 const Nom bary_header_new =
"ibubble bary_new_x bary_new_y bary_new_z";
1953 const Nom bary_vel_header =
"ibubble bary_vel_x bary_vel_y bary_vel_z bary_vect_x bary_vect_y bary_vect_z bary_vel_val";
1955 SFichier fic_bary_old = Open_file_folder(interf_dir, suffix +
".old", bary_header_old, reset, 0);
1956 for (
int ibubble=0; ibubble<nbulles_reelles; ibubble++)
1958 fic_bary_old << ibubble << end_space;
1963 fic_bary_old.
close();
1964 SFichier fic_bary_new = Open_file_folder(interf_dir, suffix +
".new", bary_header_new, reset, 0);
1965 for (
int ibubble=0; ibubble<nbulles_reelles; ibubble++)
1967 fic_bary_new << ibubble << end_space;
1972 fic_bary_new.
close();
1973 SFichier fic_bary_vel = Open_file_folder(interf_dir, suffix +
".vel", bary_vel_header, reset, 0);
1974 for (
int ibubble=0; ibubble<nbulles_reelles; ibubble++)
1976 fic_bary_vel << ibubble << end_space;
1985 fic_bary_vel.
close();
1990 DoubleTab& barycentres,
1991 const int& store_values)
1997 if (ref_ijk_ft_->schema_temps_ijk().get_tstep() == 0)
2010 for (
int i = 0; i < nbulles_reelles; i++)
2013 for (
int dir=0; dir<3; dir++)
2018 const int old_new_up_down = (pos_new - pos_old) < (-ldir / 2.);
2019 const int old_new_down_up = (pos_new - pos_old) > (ldir / 2.);
2020 if (old_new_up_down || old_new_down_up)
2022 if (old_new_up_down)
2035 for (
int dir=0; dir<3; dir++)
2048 DoubleTab& centre_gravite)
const
2054 const int nbulles_tot = nbulles_reelles + nbulles_ghost;
2055 volumes.
resize_array(nbulles_tot, RESIZE_OPTIONS::NOCOPY_NOINIT);
2057 centre_gravite.
resize(nbulles_tot, 3, RESIZE_OPTIONS::NOCOPY_NOINIT);
2058 centre_gravite = 0.;
2061 const IntTab& facettes = mesh.
facettes();
2062 const DoubleTab& sommets = mesh.
sommets();
2064 for (
int i = 0; i < n; i++)
2068 int compo = compo_facettes[i];
2075 compo = nbulles_reelles - 1 - idx_ghost;
2077 const double s = surfaces_facettes[i];
2078 const double normale_scalaire_direction = normales_facettes(i, 0);
2080 const int i0 = facettes(i, 0);
2081 const int i1 = facettes(i, 1);
2082 const int i2 = facettes(i, 2);
2083 const double coord_centre_gravite_i = (sommets(i0, 0) + sommets(i1, 0) + sommets(i2, 0)) / 3.;
2084 const double coord_centre_gravite_j = (sommets(i0, 1) + sommets(i1, 1) + sommets(i2, 1)) / 3.;
2085 const double coord_centre_gravite_k = (sommets(i0, 2) + sommets(i1, 2) + sommets(i2, 2)) / 3.;
2086 const double volume_prisme = coord_centre_gravite_i * s * normale_scalaire_direction;
2088 centre_gravite(compo, 0) += volume_prisme * (coord_centre_gravite_i * 0.5);
2089 centre_gravite(compo, 1) += volume_prisme * coord_centre_gravite_j;
2090 centre_gravite(compo, 2) += volume_prisme * coord_centre_gravite_k;
2091 volumes[compo] += volume_prisme;
2096 for (
int i = 0; i < nbulles_tot; i++)
2099 const double x = (volumes[i] == 0.) ? 0. : 1. / volumes[i];
2100 centre_gravite(i, 0) *= x;
2101 centre_gravite(i, 1) *= x;
2102 centre_gravite(i, 2) *= x;
2113 const int nbulles_tot = nbulles_reelles + nbulles_ghost;
2114 const IntTab& facettes = mesh.
facettes();
2115 const DoubleTab& sommets = mesh.
sommets();
2119 ArrOfDouble volumes;
2120 DoubleTab centre_gravite;
2123 DoubleTab d_max(nbulles_tot);
2125 DoubleTab d_min(nbulles_tot);
2131 for (
int i = 0; i < n; i++)
2135 int compo = compo_facettes[i];
2142 compo = nbulles_reelles - 1 - idx_ghost;
2145 const int i0 = facettes(i, 0);
2146 const int i1 = facettes(i, 1);
2147 const int i2 = facettes(i, 2);
2148 const double d_i_0 = sqrt( (sommets(i0, 0)-centre_gravite(compo,0))*(sommets(i0, 0)-centre_gravite(compo,0)) + (sommets(i0, 1)-centre_gravite(compo,1))*(sommets(i0, 1)-centre_gravite(compo,1)) + (sommets(i0, 2)-centre_gravite(compo,2))*(sommets(i0, 2)-centre_gravite(compo,2)) );
2149 const double d_i_1 = sqrt( (sommets(i1, 0)-centre_gravite(compo,0))*(sommets(i1, 0)-centre_gravite(compo,0)) + (sommets(i1, 1)-centre_gravite(compo,1))*(sommets(i1, 1)-centre_gravite(compo,1)) + (sommets(i1, 2)-centre_gravite(compo,2))*(sommets(i1, 2)-centre_gravite(compo,2)) );
2150 const double d_i_2 = sqrt( (sommets(i2, 0)-centre_gravite(compo,0))*(sommets(i2, 0)-centre_gravite(compo,0)) + (sommets(i2, 1)-centre_gravite(compo,1))*(sommets(i2, 1)-centre_gravite(compo,1)) + (sommets(i2, 2)-centre_gravite(compo,2))*(sommets(i2, 2)-centre_gravite(compo,2)) );
2153 d_imax = std::max(d_i_0,std::max(d_i_1,d_i_2));
2154 d_imin = std::min(d_i_0,std::min(d_i_1,d_i_2));
2157 if (d_imax > d_max[compo])
2158 d_max[compo] = d_imax;
2160 if (d_imin < d_min[compo])
2161 d_min[compo] = d_imin;
2166 for (
int i = 0; i < nbulles_tot; i++)
2168 aspect_ratio[i] = d_max[i]/d_min[i];
2187 for (
int bulle = 0; bulle < nbulles_reelles; bulle++)
2189 surfactant_min(bulle) = 1.e10;
2190 surfactant_max(bulle) = -1.e10;
2192 for (
int i = 0; i < n; i++)
2196 surfactant(compo_facettes(i))+=Surf(i)*Sfa7(i);
2197 if(Surf(i)>surfactant_max(compo_facettes(i)))
2199 surfactant_max(compo_facettes(i)) = Surf(i);
2201 if(Surf(i)<surfactant_min(compo_facettes(i)))
2203 surfactant_min(compo_facettes(i)) = Surf(i);
2215 DoubleTab& poussee)
const
2221 const int nbulles_tot = nbulles_reelles + nbulles_ghost;
2222 poussee.
resize(nbulles_tot, 3, RESIZE_OPTIONS::NOCOPY_NOINIT);
2226 const IntTab& facettes = mesh.
facettes();
2227 const DoubleTab& sommets = mesh.
sommets();
2229 for (
int i = 0; i < n; i++)
2233 int compo = compo_facettes[i];
2240 compo = nbulles_reelles - 1 - idx_ghost;
2242 const double s = surfaces_facettes[i];
2244 const int i0 = facettes(i, 0);
2245 const int i1 = facettes(i, 1);
2246 const int i2 = facettes(i, 2);
2247 const double coord_centre_gravite_i = (sommets(i0, 0) + sommets(i1, 0) + sommets(i2, 0)) / 3.;
2248 const double coord_centre_gravite_j = (sommets(i0, 1) + sommets(i1, 1) + sommets(i2, 1)) / 3.;
2249 const double coord_centre_gravite_k = (sommets(i0, 2) + sommets(i1, 2) + sommets(i2, 2)) / 3.;
2250 const double grav_scalaire_position_fois_s =
2251 (grav(0,0) * coord_centre_gravite_i + grav(0,1) * coord_centre_gravite_j + grav(0,2) * coord_centre_gravite_k) * s;
2252 for (
int dir = 0; dir < 3; dir++)
2253 poussee(compo, dir) += grav_scalaire_position_fois_s * normales_facettes(i, dir);
2270 const double vol = dxi * dxj * dxk;
2275 for (
int fa7 = 0; fa7 < n; fa7++)
2282 const double sf = surface_facettes[fa7];
2293 ai(ijk[0], ijk[1], ijk[2]) += surf / vol;
2312 const double vol = dxi * dxj * dxk;
2317 for (
int fa7 = 0; fa7 < n; fa7++)
2323 int compo_fa7 = compo_connex[fa7];
2330 if (compo_fa7 != compo)
2334 const double sf = surface_facettes[fa7];
2345 ai(ijk[0], ijk[1], ijk[2]) += surf / vol;
2364 const double vol = dxi * dxj * dxk;
2366 for (
int fa7 = 0; fa7 < n; fa7++)
2373 const double sf = surface_facettes[fa7];
2374 const int compo_fa7 = compo_connex[fa7];
2376 if (compo_fa7 != compo)
2389 if (ijk[0]==i_ref and ijk[1]==j_ref and ijk[2]==k_ref)
2401static int encoder_compo(
int num_bulle,
int code_deplacement)
2407 return (num_bulle << 6) | (code_deplacement);
2410static int decoder_numero_bulle(
const int code)
2412 const int num_bulle = code >> 6;
2423 const IntTab& facettes = mesh.
facettes();
2424 const ArrOfDouble& courbure =
maillage_ft_ijk_.get_update_courbure_sommets();
2430 IJK_Field_double SI_ft;
2434 kappa_ft.
data() = 0.;
2436 for (
int fa7 = 0; fa7 < n; fa7++)
2442 const double sf=surface_facettes[fa7];
2450 for (
int isom = 0; isom< 3; isom++)
2452 const int num_som = facettes(fa7, isom);
2453 const double kappa = courbure[num_som];
2454 kappa_ft(ijk[0],ijk[1],ijk[2]) += kappa*surf/3.;
2456 SI_ft(ijk[0],ijk[1],ijk[2]) += surf;
2462 const int nx = kappa_ft.
ni();
2463 const int ny = kappa_ft.
nj();
2464 const int nz = kappa_ft.
nk();
2466 for (
int k = 0; k < nz; k++)
2467 for (
int j = 0; j < ny; j++)
2468 for (
int i = 0; i < nx; i++)
2470 if (kappa_ft(i,j,k)!=0.)
2471 kappa_ft(i,j,k)/=SI_ft(i,j,k);
2488 IJK_Field_vector3_double& normale_cell,
2489 const int igroup)
const
2495 const IntTab& facettes = mesh.
facettes();
2496 const ArrOfDouble& courbure =
maillage_ft_ijk_.get_update_courbure_sommets();
2503 const double vol = dxi * dxj * dxk;
2507 kappa_ai.
data() = 0.;
2508 for (
int dir = 0; dir < 3; dir++)
2509 normale_cell[dir].data() = 0.;
2512 for (
int fa7 = 0; fa7 < n; fa7++)
2521 int icompo = compo_facette[fa7];
2525 icompo = decoder_numero_bulle(-icompo);
2530 const double sf=surface_facettes[fa7];
2539 for (
int dir = 0; dir< 3; dir++)
2541 const double nx = normale_facettes(fa7,dir);
2542 normale_cell[dir](ijk[0],ijk[1],ijk[2]) += nx * surf;
2544 const double fac = surf/(3.*vol);
2545 for (
int isom = 0; isom< 3; isom++)
2547 const int num_som = facettes(fa7, isom);
2548 const double kappa = courbure[num_som];
2549 kappa_ai(ijk[0],ijk[1],ijk[2]) += kappa*fac;
2552 ai(ijk[0],ijk[1],ijk[2]) += surf/vol;
2560 const int nx = ai.
ni();
2561 const int ny = ai.
nj();
2562 const int nz = ai.
nk();
2564 for (
int k = 0; k < nz; k++)
2565 for (
int j = 0; j < ny; j++)
2566 for (
int i = 0; i < nx; i++)
2568 double norme_carre = 0.;
2569 for (
int dir = 0; dir < 3; dir++)
2571 const double dnx = normale_cell[dir](i, j, k);
2572 norme_carre += dnx * dnx;
2575 if (norme_carre > 0.)
2576 for (
int dir = 0; dir < 3; dir++)
2577 normale_cell[dir](i, j, k) *= 1. / sqrt(norme_carre);
2584 const ArrOfDouble& surface_facette,
2585 const ArrOfDouble& surface_par_bulle,
2586 const ArrOfInt& compo_connexes_facettes,
2587 const int nbulles_reelles,
2588 const int nbulles_ghost,
2589 const DoubleTab& vitesse_sommets,
2590 DoubleTab& vitesses_translation_bulles)
const
2592 const int nbulles_tot = nbulles_reelles + nbulles_ghost;
2593 assert(vitesses_translation_bulles.
dimension(0) == nbulles_tot);
2594 assert(surface_par_bulle.
size_array() == nbulles_tot);
2595 assert(vitesses_translation_bulles.
dimension(1) == 3);
2597 const IntTab& facettes = maillage.
facettes();
2599 vitesses_translation_bulles = 0.;
2607 for (
int fa7 = 0; fa7 < maillage.
nb_facettes(); fa7++)
2614 int compo = compo_connexes_facettes[fa7];
2622 compo = nbulles_reelles - 1 - idx_ghost;
2625 assert(compo < nbulles_tot);
2627 const double sf = surface_facette[fa7];
2630 for (
int dir = 0; dir < 3; dir++)
2634 for (
int j = 0; j < 3; j++)
2636 const int isom = facettes(fa7, j);
2637 v += vitesse_sommets(isom, dir);
2639 vitesses_translation_bulles(compo, dir) += v * sf / 3.;
2645 for (
int icompo = 0; icompo < nbulles_tot; icompo++)
2646 if (surface_par_bulle[icompo] > 0.)
2647 for (
int dir = 0; dir < 3; dir++)
2648 vitesses_translation_bulles(icompo, dir) /= surface_par_bulle[icompo];
2652 const ArrOfDouble& surface_facette,
2653 const ArrOfDouble& surface_par_bulle,
2654 const ArrOfInt& compo_connexes_facettes,
2655 const int nbulles_reelles,
2656 const int nbulles_ghost,
2657 const DoubleTab& centre_gravite,
2658 const DoubleTab& vitesse_sommets,
2659 const DoubleTab& vitesse_translation_sommets,
2660 DoubleTab& mean_bubble_rotation_vector)
const
2662 const int nbulles_tot = nbulles_reelles + nbulles_ghost;
2663 assert(mean_bubble_rotation_vector.
dimension(0) == nbulles_tot);
2664 assert(surface_par_bulle.
size_array() == nbulles_tot);
2665 assert(mean_bubble_rotation_vector.
dimension(1) == 3);
2667 const IntTab& facettes = maillage.
facettes();
2668 const DoubleTab& sommets = maillage.
sommets();
2670 mean_bubble_rotation_vector = 0.;
2677 for (
int fa7 = 0; fa7 < maillage.
nb_facettes(); fa7++)
2683 int compo = compo_connexes_facettes[fa7];
2691 compo = nbulles_reelles - 1 - idx_ghost;
2694 assert(compo < nbulles_tot);
2696 const double sf = surface_facette[fa7];
2699 for (
int j = 0; j < 3; j++)
2701 const int isom = facettes(fa7, j);
2703 double vit_x = vitesse_sommets(isom, 0) - vitesse_translation_sommets(compo, 0);
2704 double vit_y = vitesse_sommets(isom, 1) - vitesse_translation_sommets(compo, 1);
2705 double vit_z = vitesse_sommets(isom, 2) - vitesse_translation_sommets(compo, 2);
2706 Vecteur3 vit = {vit_x, vit_y, vit_z};
2708 double x = sommets(isom, 0);
2709 double y = sommets(isom, 1);
2710 double z = sommets(isom, 2);
2713 Vecteur3 coord_centre = {centre_gravite(compo,0), centre_gravite(compo,1), centre_gravite(compo,2)};
2715 Vecteur3 radial_position = coord_sommet - coord_centre;
2721 for (
int dir = 0; dir < 3; dir++)
2725 mean_bubble_rotation_vector(compo, dir) += rotation_vector[dir] * sf / 3.;
2732 for (
int icompo = 0; icompo < nbulles_tot; icompo++)
2733 if (surface_par_bulle[icompo] > 0.)
2734 for (
int dir = 0; dir < 3; dir++)
2735 mean_bubble_rotation_vector(icompo, dir) /= surface_par_bulle[icompo];
2747static void calculer_normale_sommets_interface(
const Maillage_FT_IJK& maillage,
2756 const IntTab& facettes = maillage.
facettes();
2757 const int nsommets_faces = facettes.
dimension(1);
2759 normale.
resize(nsom, dim);
2761 double n[3] = {0., 0. , 0.};
2763 ArrOfDouble surface_environnante(nsom);
2764 surface_environnante = 0.;
2766 for (
int i = 0; i < nfaces; i++)
2773 const double surface = surface_facettes[i];
2775 for (
int k = 0; k < dim; k++)
2776 n[k] = normale_facettes(i, k) * surface;
2778 for (
int j = 0; j < nsommets_faces; j++)
2780 const int sommet = facettes(i, j);
2781 surface_environnante[sommet] += surface;
2782 for (
int k = 0; k < dim; k++)
2783 normale(sommet, k) += n[k];
2794 int count[9] = {0, 0, 0, 0, 0, 0, 0, 0, 0};
2795 for (
int isom = 0; isom < nsom; isom++)
2798 for (
int dir = 0; dir < 3; dir++)
2800 if (surface_environnante[isom])
2801 normale(isom, dir) /= surface_environnante[isom];
2803 normale(isom, dir) = 0.;
2806 norme = normale(isom, 0) * normale(isom, 0) + normale(isom, 1) * normale(isom, 1) +
2807 normale(isom, 2) * normale(isom, 2);
2873 Cerr <<
"IJK_Interfaces.cpp:calculer_normale_sommets_interface : Calcul "
2874 "normale non-unitaire : ";
2876 for (
int i = 0; i < 9; i++)
2879 Cerr <<
" " << count[i] / nsom;
2881 Cerr <<
" " << 1 - sum / nsom <<
" (nsom = " << nsom << finl;
2886 const DoubleTab& vitesses_translation_bulles,
2887 const DoubleTab& mean_bubble_rotation_vector,
2888 const DoubleTab& centre_gravite,
2889 ArrOfDouble& var_volume)
2891 Cut_field_vector3_double& cut_field_deformation_velocity =
static_cast<Cut_field_vector3_double&
>(
deformation_velocity_);
2894 const DoubleTab& sommets = mesh.
sommets();
2899 const int nbulles_tot = nbulles_reelles + nbulles_ghost;
2902 var_volume.
resize(nbsom);
2914 for (
int dir = 0 ; dir < 3 ; dir++)
2919 DoubleTab bounding_box;
2925 for (
int icompo = 0; icompo < nbulles_tot; icompo++)
2935 cut_field_deformation_velocity,
2941 DoubleVect delta_volume_theorique_bilan_ns_vdf;
2945 const Domaine& domaine = domaine_vf.
domaine();
2951 assert(ni ==
I_ft().ni());
2952 assert(nj ==
I_ft().nj());
2953 assert(nk ==
I_ft().nk());
2955 for (
int k = 0; k < nk; k++)
2957 for (
int j = 0; j < nj; j++)
2959 for (
int i = 0; i < ni; i++)
2961 const int num_elem = ref_domaine_->convert_ijk_cell_to_packed(i, j, k);
2972 DoubleTab& mean_bubble_rotation_vector,
2973 DoubleTab& centre_gravite,
2974 const int first_step_interface_smoothing)
2980 if (first_step_interface_smoothing)
2993 const int nbulles_tot = nbulles_reelles + nbulles_ghost;
3001 ArrOfDouble surface_par_bulle;
3005 ArrOfDouble volume_par_bulle(nbulles_tot);
3006 vitesses_translation_bulles.
resize(nbulles_tot, 3);
3007 mean_bubble_rotation_vector.
resize(nbulles_tot, 3);
3008 centre_gravite.
resize(nbulles_tot, 3);
3023 vitesses_translation_bulles);
3038 vitesses_translation_bulles,
3039 mean_bubble_rotation_vector);
3041 if (first_step_interface_smoothing)
3043 vitesses_translation_bulles = 0.;
3044 mean_bubble_rotation_vector = 0.;
3053 for (
int bulle = 0; bulle < nbulles_tot; bulle++)
3055 Cerr <<
"composante " << bulle <<
" vitesse_translation ";
3056 for (
int i = 0; i < 3; i++)
3057 Cerr <<
" " << vitesses_translation_bulles(bulle, i);
3058 Cerr <<
" rotation_vector ";
3059 for (
int i = 0; i < 3; i++)
3060 Cerr <<
" " << mean_bubble_rotation_vector(bulle, i);
3078 const DoubleTab& vitesses_translation_bulles,
3079 const DoubleTab& mean_bubble_rotation_vector,
3080 const DoubleTab& centre_gravite,
3081 const double dt_tot,
3083 const int rk_step = -1,
3084 const int first_step_interface_smoothing)
3088 const DoubleTab& sommets = mesh.
sommets();
3090 DoubleTab deplacement(nbsom, 3);
3100 const int nbulles_tot = nbulles_reelles + nbulles_ghost;
3108 ArrOfIntFT compo_connex_som;
3113 DoubleTabFT normale_sommet;
3114 calculer_normale_sommets_interface(mesh, normale_sommet);
3116 for (
int som = 0; som < nbsom; som++)
3121 deplacement(som, 0) = 100.;
3122 deplacement(som, 1) = 100.;
3123 deplacement(som, 2) = 100.;
3130 int icompo = compo_connex_som[som];
3138 icompo = nbulles_reelles - 1 - idx_ghost;
3140 assert(icompo >= 0);
3141 assert(icompo < nbulles_tot);
3159 double rot_x = mean_bubble_rotation_vector(icompo, 0);
3160 double rot_y = mean_bubble_rotation_vector(icompo, 1);
3161 double rot_z = mean_bubble_rotation_vector(icompo, 2);
3162 Vecteur3 rot = {rot_x, rot_y, rot_z};
3164 double x = sommets(som, 0);
3165 double y = sommets(som, 1);
3166 double z = sommets(som, 2);
3169 Vecteur3 coord_centre = {centre_gravite(icompo,0), centre_gravite(icompo,1), centre_gravite(icompo,2)};
3171 Vecteur3 radial_position = coord_sommet - coord_centre;
3176 double prodscal = 0.;
3177 double norme_carre = 0.;
3178 for (
int direction = 0; direction < 3; direction++)
3180 double vi =
vinterp_(som, direction);
3181 double vcompo = vitesses_translation_bulles(icompo, direction) + tangential_velocity[direction];
3182 double n = normale_sommet(som, direction);
3183 prodscal += (vi - vcompo) * n;
3184 norme_carre += n * n;
3187 prodscal /= norme_carre;
3191 for (
int direction = 0; direction < 3; direction++)
3193 double n = normale_sommet(som, direction);
3197 if (correction_semi_locale_volume_bulle)
3200 vinterp_(som, direction) = n * prodscal;
3205 vinterp_(som, direction) = n * prodscal + vitesses_translation_bulles(icompo, direction);
3234 Journal() <<
"RKSTEP " << rk_step << finl;
3235 for (
int i = 0; i < nbsom; i++)
3238 for (
int direction = 0; direction < 3; direction++)
3240 Journal() << vinterp(i, direction) <<
" ";
3242 Journal() <<
" RK3_G_store_vi_ ";
3243 for (
int direction = 0; direction < 3; direction++)
3265 for (
int som = 0; som < nbsom; som++)
3266 for (
int direction = 0; direction < 3; direction++)
3267 deplacement(som, direction) =
vinterp_(som, direction) * dt_tot;
3274 DoubleTab coord_sommets_avant_deplacement(mesh.
sommets());
3290 dt = compute_fractionnal_timestep_rk3(dt_tot, rk_step);
3334 if (!correction_semi_locale_volume_bulle)
3339 ArrOfDouble var_volume_par_bulle(nbulles_tot);
3344 for (
int isom = 0; isom < nbsom; isom++)
3350 int compo = compo_connex_som[isom];
3358 compo = nbulles_reelles - 1 - idx_ghost;
3361 assert(compo < nbulles_tot);
3364 var_volume_par_bulle[compo] += v;
3370 for (
int icompo = 0; icompo < nbulles_tot; icompo++)
3372 dvol[icompo] += var_volume_par_bulle[icompo];
3384 const DoubleTab& vitesses_translation_bulles,
3385 const DoubleTab& mean_bubble_rotation_vector,
3386 const DoubleTab& centre_gravite,
3389 const int rk_step = -1,
3390 const double temps = -1 )
3393 const DoubleTab& sommets = mesh.
sommets();
3396 if (correction_semi_locale_volume_bulle)
3398 const double fractionnal_timestep = (rk_step == -1) ? dt_tot : compute_fractionnal_timestep_rk3(dt_tot, rk_step);
3417 ArrOfDouble surface_par_bulle;
3433 if ((rk_step == -1) || (rk_step == 2))
3436 const IntTab& facettes = mesh.
facettes();
3437 for (
int i = 0; i < nb_facettes; i++)
3442 int compo = compo_connex_apres_transport[i];
3449 compo = nbulles_reelles - 1 - idx_ghost;
3451 const double var_volume_tot = dvol[compo];
3452 const double surface_bulle = surface_par_bulle[compo];
3453 const double s = surface_facette_apres_transport[i];
3454 const double dv = var_volume_tot * s / (surface_bulle * 3.);
3455 for (
int j = 0; j < 3; j++)
3457 const int isom = facettes(i, j);
3473 if ((rk_step == -1) || (rk_step == 2))
3499 const int nbsom_end_transport = mesh.
nb_sommets();
3501 if (!((nbsom_end_transport == size_store) || (0 == size_store)))
3503 Cerr <<
"Une des tailles de tableau n'est pas bonne... "
3504 <<
" size_store = " << size_store
3505 <<
" nbsom_end_transport = " << nbsom_end_transport
3519 const DoubleTab& vitesses_translation_bulles,
3520 const DoubleTab& mean_bubble_rotation_vector,
3521 const DoubleTab& centre_gravite,
3523 const int first_step_interface_smoothing)
3527 const DoubleTab& sommets = mesh.
sommets();
3529 DoubleTab deplacement(nbsom, 3);
3533 ArrOfIntFT compo_connex_som;
3536 for (
int som = 0; som < nbsom; som++)
3541 deplacement(som, 0) = 100.;
3542 deplacement(som, 1) = 100.;
3543 deplacement(som, 2) = 100.;
3547 int icompo = compo_connex_som[som];
3555 icompo = nbulles_reelles - 1 - idx_ghost;
3557 assert(icompo >= 0);
3560 double rot_x = mean_bubble_rotation_vector(icompo, 0);
3561 double rot_y = mean_bubble_rotation_vector(icompo, 1);
3562 double rot_z = mean_bubble_rotation_vector(icompo, 2);
3563 Vecteur3 rot = {rot_x, rot_y, rot_z};
3565 double x = sommets(som, 0);
3566 double y = sommets(som, 1);
3567 double z = sommets(som, 2);
3570 Vecteur3 coord_centre = {centre_gravite(icompo,0), centre_gravite(icompo,1), centre_gravite(icompo,2)};
3572 Vecteur3 radial_position = coord_sommet - coord_centre;
3577 for (
int direction = 0; direction < 3; direction++)
3579 const double fractionnal_timestep = (rk_step == -1) ? dt_tot : compute_fractionnal_timestep_rk3(dt_tot, rk_step);
3580 deplacement(som, direction) = (vitesses_translation_bulles(icompo, direction) + tangential_velocity[direction]) * fractionnal_timestep;
3590 ArrOfDouble& var_volume,
3609 if (algo_remaillage_local.
a_remailler(temps, maillage))
3621 const DoubleTab& bounding_box_bulles,
3622 const Cut_field_vector3_double& cut_field_velocity,
3623 const DoubleTab& vitesses_translation_bulles,
3624 const DoubleTab& mean_bubble_rotation_vector,
3625 const DoubleTab& positions_bulles)
3627 Cut_field_vector3_double& cut_field_deformation_velocity =
static_cast<Cut_field_vector3_double&
>(
deformation_velocity_);
3628 assert(&cut_field_deformation_velocity[0].get_cut_cell_disc() == &cut_field_deformation_velocity[1].get_cut_cell_disc());
3629 assert(&cut_field_deformation_velocity[0].get_cut_cell_disc() == &cut_field_deformation_velocity[2].get_cut_cell_disc());
3630 const Cut_cell_FT_Disc& cut_cell_disc = cut_field_deformation_velocity[0].get_cut_cell_disc();
3632 const int ni = cut_field_deformation_velocity[0].ni();
3633 const int nj = cut_field_deformation_velocity[0].nj();
3634 const int nk = cut_field_deformation_velocity[0].nk();
3635 const int ghost = cut_field_deformation_velocity[0].
ghost();
3646 double origin_x = geom.
get_origin(DIRECTION_I);
3647 double origin_y = geom.
get_origin(DIRECTION_J);
3648 double origin_z = geom.
get_origin(DIRECTION_K);
3654 int imin = std::max(-ghost, (
int)((bounding_box_bulles(compo, 0, 0) - origin_x)/dx - offset_x - 2));
3655 int imax = std::min(ni + ghost, (
int)((bounding_box_bulles(compo, 0, 1) - origin_x)/dx - offset_x + 2));
3656 int jmin = std::max(-ghost, (
int)((bounding_box_bulles(compo, 1, 0) - origin_y)/dy - offset_y - 2));
3657 int jmax = std::min(nj + ghost, (
int)((bounding_box_bulles(compo, 1, 1) - origin_y)/dy - offset_y + 2));
3658 int kmin = std::max(-ghost, (
int)((bounding_box_bulles(compo, 2, 0) - origin_z)/dz - offset_z - 2));
3659 int kmax = std::min(nk + ghost, (
int)((bounding_box_bulles(compo, 2, 1) - origin_z)/dz - offset_z + 2));
3660 for (
int k = kmin; k < kmax; k++)
3662 for (
int j = jmin; j < jmax; j++)
3664 for (
int i = imin; i < imax; i++)
3666 const double x_centre_cell = (i + offset_x + .5)*dx + origin_x;
3667 const double y_centre_cell = (j + offset_y + .5)*dy + origin_y;
3668 const double z_centre_cell = (k + offset_z + .5)*dz + origin_z;
3670 double rot_x = mean_bubble_rotation_vector(compo, 0);
3671 double rot_y = mean_bubble_rotation_vector(compo, 1);
3672 double rot_z = mean_bubble_rotation_vector(compo, 2);
3673 Vecteur3 rot = {rot_x, rot_y, rot_z};
3675 int n = cut_cell_disc.
get_n(i,j,k);
3678 for (
int phase = 0 ; phase < 2 ; phase++)
3684 Vecteur3 coord = {x_cut_cell, y_cut_cell, z_cut_cell};
3686 Vecteur3 coord_centre = {positions_bulles(compo,0), positions_bulles(compo,1), positions_bulles(compo,2)};
3688 Vecteur3 radial_position = coord - coord_centre;
3693 for (
int direction = 0; direction < 3; direction++)
3695 const DoubleTabFT_cut_cell& diph_velocity = (phase == 0) ? cut_field_velocity[direction].diph_v_ : cut_field_velocity[direction].diph_l_;
3696 DoubleTabFT_cut_cell& deformation_diph_velocity = (phase == 0) ? cut_field_deformation_velocity[direction].diph_v_ : cut_field_deformation_velocity[direction].diph_l_;
3697 deformation_diph_velocity(n) = diph_velocity(n) - vitesses_translation_bulles(compo, direction) - tangential_velocity[direction];
3703 Vecteur3 coord = {x_centre_cell, y_centre_cell, z_centre_cell};
3705 Vecteur3 coord_centre = {positions_bulles(compo,0), positions_bulles(compo,1), positions_bulles(compo,2)};
3707 Vecteur3 radial_position = coord - coord_centre;
3712 for (
int direction = 0; direction < 3; direction++)
3714 cut_field_deformation_velocity[direction].pure_(i,j,k) = cut_field_velocity[direction].pure_(i,j,k) - vitesses_translation_bulles(compo, direction) - tangential_velocity[direction];
3729 bounding_box.
resize(nbulles, 3, 2);
3733 const double inval = 1.e20;
3734 bounding_box = inval;
3738 ArrOfIntFT compo_connex_som;
3740 const DoubleTab& sommets = mesh.
sommets();
3742 ArrOfDouble volume_reel;
3746 ArrOfDouble position_xmax_compo;
3747 ArrOfDouble position_xmin_compo;
3750 position_xmax_compo.
resize_array(nbulles, RESIZE_OPTIONS::NOCOPY_NOINIT);
3751 position_xmin_compo.
resize_array(nbulles, RESIZE_OPTIONS::NOCOPY_NOINIT);
3752 position_xmax_compo = -10000.;
3753 position_xmin_compo = 10000.;
3755 for (
int i_sommet2 = 0; i_sommet2 < nbsom; i_sommet2++)
3757 int iconnex = compo_connex_som[i_sommet2];
3758 double coord = sommets(i_sommet2, 0);
3759 if (coord>position_xmax_compo[iconnex])
3760 position_xmax_compo[iconnex] = coord;
3761 if (coord<position_xmin_compo[iconnex])
3762 position_xmin_compo[iconnex] = coord;
3768 for (
int i_sommet = 0; i_sommet < nbsom; i_sommet++)
3770 for (direction=0; direction<3; direction++)
3772 double coord = sommets(i_sommet, direction);
3773 int iconnex = compo_connex_som[i_sommet];
3778 double pos_ref = position(iconnex,0);
3787 double decallage_bulle_reel_ext_domaine_reel = 1.*(position_xmax_compo[iconnex]-position_xmin_compo[iconnex]);
3790 double pos = std::fmod(std::fmod(pos_ref + offset - decallage_bulle_reel_ext_domaine_reel, Lx) + Lx, Lx) + decallage_bulle_reel_ext_domaine_reel;
3792 coord += (pos - pos_ref);
3799 double mmin = bounding_box(iconnex, direction, 0);
3800 double mmax = -bounding_box(iconnex, direction, 1);
3802 bounding_box(iconnex, direction, 0) = coord;
3804 bounding_box(iconnex, direction, 1) = -coord;
3815 for (
int ibulle = 0; ibulle < nbulles; ibulle++)
3816 for (direction = 0; direction < 3; direction++)
3817 bounding_box(ibulle, direction, 1) = -bounding_box(ibulle, direction, 1);
3826 DoubleTab bounding_box;
3828 ArrOfInt masque_duplicata_pour_compo_reel;
3841 DoubleTab bounding_box_offsetp;
3844 DoubleTab bounding_box_offsetm;
3863static int decoder_deplacement(
const int code,
const int dir,
int& compo_bulle_reel)
3866 const int code_deplacement = code & 63;
3867 int tmp = code_deplacement & (1 << (3 + dir));
3873 int signe = (tmp == 0) ? -1 : 1 ;
3874 tmp = code_deplacement & (1 << dir);
3876 int index = tmp >> dir;
3880 compo_bulle_reel = code >> 6;
3884 return signe * index;
3891 DoubleTab& deplacement,
3892 DoubleTab& bounding_box_NS, DoubleTab position,
3897 deplacement.
resize(nbsom, 3);
3899 ArrOfIntFT compo_sommets;
3903 ArrOfIntFT compo_sommets_init;
3907 ArrOfDouble position_xmax_compo;
3908 ArrOfDouble position_xmin_compo;
3911 position_xmax_compo.
resize_array(nbulles, RESIZE_OPTIONS::NOCOPY_NOINIT);
3912 position_xmin_compo.
resize_array(nbulles, RESIZE_OPTIONS::NOCOPY_NOINIT);
3913 position_xmax_compo = -10000.;
3914 position_xmin_compo = 10000.;
3916 const DoubleTab& sommets = mesh.
sommets();
3918 for (
int i_sommet = 0; i_sommet < nbsomreel; i_sommet++)
3922 int iconnex = compo_sommets_init[i_sommet];
3924 double coord = sommets(i_sommet, 0);
3925 if (coord>position_xmax_compo[iconnex])
3926 position_xmax_compo[iconnex] = coord;
3927 if (coord<position_xmin_compo[iconnex])
3928 position_xmin_compo[iconnex] = coord;
3934 for (
int dir = 0; dir < 3; dir++)
3936 for (
int i_sommet = 0; i_sommet < nbsom; i_sommet++)
3939 const int code = compo_sommets[i_sommet];
3944 int compo_bulle_reel = 0;
3946 int decode = decoder_deplacement(code, dir, compo_bulle_reel);
3948 double depl = decode * (bounding_box_NS(dir, 1) - bounding_box_NS(dir, 0));
3949 deplacement(i_sommet, dir) = depl;
3950 double pos_ref = 0 ;
3952 double decallage_bulle_reel_ext_domaine_reel = 0.;
3959 pos_ref = position(compo_bulle_reel,0);
3961 decallage_bulle_reel_ext_domaine_reel = 1.*(position_xmax_compo[compo_bulle_reel]-position_xmin_compo[compo_bulle_reel]);
3962 pos = std::fmod(std::fmod(deplacement(i_sommet, 0) + pos_ref + offset - decallage_bulle_reel_ext_domaine_reel, Lx) + Lx, Lx) + decallage_bulle_reel_ext_domaine_reel;
3964 deplacement(i_sommet, 0) += (pos - pos_ref);
3976static void calculer_deplacement_from_code_compo_connexe_negatif(
const Maillage_FT_IJK& m,
3977 DoubleTab& deplacement,
3978 DoubleTab& bounding_box_NS)
3982 deplacement.
resize(nbsom, 3);
3983 ArrOfIntFT compo_sommets;
3986 for (
int i_sommet = 0; i_sommet < nbsom; i_sommet++)
3989 int code = compo_sommets[i_sommet];
3994 for (
int dir = 0; dir < 3; dir++)
3995 deplacement(i_sommet, dir) = 0.;
4003 for (
int dir = 0; dir < 3; dir++)
4006 int unused_variable = 0 ;
4007 int decode = decoder_deplacement(code, dir, unused_variable);
4008 double depl = decode * (bounding_box_NS(dir, 1) - bounding_box_NS(dir, 0));
4009 deplacement(i_sommet, dir) = depl;
4019static void calculer_deplacement_from_masque_in_array(
const Maillage_FT_IJK& m,
4020 DoubleTab& deplacement,
4021 const ArrOfInt& masque_array,
4022 DoubleTab& bounding_box_NS, DoubleTab position,
const int nbulles)
4026 deplacement.
resize(nbsom,3);
4027 ArrOfIntFT compo_connexe_som;
4030 ArrOfDouble position_xmax_compo;
4031 ArrOfDouble position_xmin_compo;
4034 position_xmax_compo.
resize_array(nbulles, RESIZE_OPTIONS::NOCOPY_NOINIT);
4035 position_xmin_compo.
resize_array(nbulles, RESIZE_OPTIONS::NOCOPY_NOINIT);
4036 position_xmax_compo = -10000.;
4037 position_xmin_compo = 10000.;
4039 const DoubleTab& sommets = m.
sommets();
4041 for (
int i_sommet = 0; i_sommet < nbsom; i_sommet++)
4043 int iconnex = compo_connexe_som[i_sommet];
4044 double coord = sommets(i_sommet, 0);
4045 if (coord>position_xmax_compo[iconnex])
4046 position_xmax_compo[iconnex] = coord;
4047 if (coord<position_xmin_compo[iconnex])
4048 position_xmin_compo[iconnex] = coord;
4054 for (
int i_sommet = 0; i_sommet < nbsom; i_sommet++)
4057 const int icompo = compo_connexe_som[i_sommet];
4058 const int code = masque_array[icompo];
4060 for (
int dir = 0; dir < 3; dir++)
4063 int compo_bulle_reel = 0;
4065 int decode = decoder_deplacement(code, dir, compo_bulle_reel);
4067 double depl = decode * (bounding_box_NS(dir, 1) - bounding_box_NS(dir, 0));
4068 deplacement(i_sommet, dir) = depl;
4069 double pos_ref = 0 ;
4071 double decallage_bulle_reel_ext_domaine_reel = 0.;
4078 pos_ref = position(compo_bulle_reel,0);
4080 decallage_bulle_reel_ext_domaine_reel = 1.*(position_xmax_compo[compo_bulle_reel]-position_xmin_compo[compo_bulle_reel]);
4081 pos = std::fmod(std::fmod(deplacement(i_sommet, 0) + pos_ref + offset - decallage_bulle_reel_ext_domaine_reel, Lx) + Lx, Lx) + decallage_bulle_reel_ext_domaine_reel;
4083 deplacement(i_sommet, 0) += (pos - pos_ref);
4108 maillage_temporaire.
initialize(ref_domaine_.valeur(), refdomaine_dis_.valeur(),
parcours_);
4119 int icompo, index, signe;
4120 for (
int i_facette = 0; i_facette < a_dupliquer.
nb_facettes(); i_facette++)
4122 icompo = compo_connexe_facettes[i_facette];
4139 ArrOfInt liste_bulles_crees;
4141 int nbulles_crees = 0;
4142 for (
int mon_numero_iteration = 1; ; mon_numero_iteration++)
4146 ArrOfInt index_copie(nbulles);
4147 ArrOfInt index_signe(nbulles);
4150 ArrOfInt liste_facettes_pour_suppression;
4151 int reste_a_faire = 0;
4155 for (icompo = 0; icompo < nbulles; icompo++)
4159 index_copie[icompo] = -1;
4163 const int masque = masque_duplicata_pour_compo[icompo] & 7;
4164 for (index = 1; index <= 7; index++)
4166 if ((masque & index) == index)
4169 if (compteur == mon_numero_iteration)
4173 index_copie[icompo] = index;
4184 index_signe[icompo] = -1;
4185 ArrOfInt index_bulle(7);
4186 ArrOfInt signe_bulle(7);
4188 const int perio_x_reel = masque_duplicata_pour_compo[icompo] & 1;
4189 const int perio_y_reel = masque_duplicata_pour_compo[icompo] & (1 << 1);
4190 const int perio_z_reel = masque_duplicata_pour_compo[icompo] & (1 << 2);
4191 const int perio_x_ghost = masque_duplicata_pour_compo[icompo] & (1 << 3);
4194 int signe_3bit = (masque_duplicata_pour_compo[icompo] & (7 << 4)) >> 1;
4199 int bit_position_x_signe_6bit = 3;
4200 int bit_position_x_ghost_8bit = 7;
4201 int signe_x_ghost = (masque_duplicata_pour_compo[icompo] & (1 << bit_position_x_ghost_8bit)) >> bit_position_x_ghost_8bit;
4202 int mask = 1<<bit_position_x_signe_6bit;
4204 for (
int nb_bulle_duplique = 0; nb_bulle_duplique < 7; nb_bulle_duplique++)
4206 index_bulle[nb_bulle_duplique] = -1;
4207 signe_bulle[nb_bulle_duplique] = signe_3bit;
4229 if (perio_z_reel == 4)
4231 if (perio_y_reel == 2 && perio_x_reel == 1 && perio_x_ghost == 8 )
4239 signe_bulle[4] &= (~mask);
4240 signe_bulle[4] |= (signe_x_ghost << bit_position_x_signe_6bit);
4243 signe_bulle[6] &= (~mask);
4244 signe_bulle[6] |= (signe_x_ghost << bit_position_x_signe_6bit);
4246 else if (perio_y_reel == 2 && perio_x_reel == 1 && perio_x_ghost != 8 )
4254 else if (perio_y_reel == 2 && perio_x_reel != 1 && perio_x_ghost == 8 )
4257 signe_bulle[0] &= (~mask);
4258 signe_bulle[0] |= (signe_x_ghost << bit_position_x_signe_6bit);
4261 signe_bulle[2] &= (~mask);
4262 signe_bulle[2] |= (signe_x_ghost << bit_position_x_signe_6bit);
4266 else if (perio_y_reel != 2 && perio_x_reel == 1 && perio_x_ghost == 8 )
4271 signe_bulle[2] &= (~mask);
4272 signe_bulle[2] |= (signe_x_ghost << bit_position_x_signe_6bit);
4275 else if (perio_y_reel != 2 && perio_x_reel != 1 && perio_x_ghost != 8 )
4279 else if (perio_y_reel != 2 && perio_x_reel != 1 && perio_x_ghost == 8 )
4282 signe_bulle[0] &= (~mask);
4283 signe_bulle[0] |= (signe_x_ghost << bit_position_x_signe_6bit);
4286 else if (perio_y_reel != 2 && perio_x_reel == 1 && perio_x_ghost != 8 )
4291 else if (perio_y_reel == 2 && perio_x_reel != 1 && perio_x_ghost != 8 )
4299 Cerr <<
"Erreur dans la gestion des bulles ghost" << finl;
4305 if (perio_y_reel == 2 && perio_x_reel == 1)
4311 else if (perio_y_reel == 2 && perio_x_reel != 1)
4315 else if (perio_y_reel != 2 && perio_x_reel == 1)
4319 else if (perio_y_reel != 2 && perio_x_reel != 1)
4321 index_bulle[0] = -1;
4325 Cerr <<
"Erreur dans la gestion des bulles ghost" << finl;
4330 if (index_bulle[mon_numero_iteration-1] != -1 && mon_numero_iteration-1 <= 6)
4333 index_copie[icompo] = index_bulle[mon_numero_iteration-1];
4334 index_signe[icompo] = signe_bulle[mon_numero_iteration-1];
4341 if (reste_a_faire == 0)
4349 for (
int i_facette = 0; i_facette < nf; i_facette++)
4351 icompo = decoder_numero_bulle(compo_connexe_facettes[i_facette]);
4353 index = index_copie[icompo];
4355 signe = index_signe[icompo];
4378 signe = masque_duplicata_pour_compo[icompo] & (7 << 3);
4379 const int code_deplacement = signe | index;
4390 Cerr <<
"Comment est-ce possible? " << finl;
4395 liste_facettes_pour_suppression.
append_array(i_facette);
4419 DoubleTab deplacement;
4424 ArrOfDouble volume_reel;
4428 calculer_deplacement_from_code_compo_connexe(maillage_temporaire, split,
4440 for (
int iface = 0; iface < nbf; iface++)
4442 const int code_negatif = -compo_connexe_facettes_tempo[iface];
4452 array_trier_retirer_doublons(liste_bulles_crees);
4463 for (
int p = 1; p < nbproc; p++)
4465 recevoir(prov, p, 0, 2001 + p);
4475 array_trier_retirer_doublons(liste_bulles_crees);
4478 nbulles_crees = liste_bulles_crees.
size_array();
4485 Cerr <<
" On vient de creer des bulles ghost alors qu'il y en avait "
4487 <<
" tableau ghost_compo_converter_ n'est pas de taille nulle! " << finl;
4504 Cerr <<
"IJK_Interfaces.cpp::dupliquer_bulle_perio liste_bulles_crees : "
4505 << liste_bulles_crees << finl;
4509 if (new_size != old_size)
4511 Cerr <<
"On dirait que des doublons etaient presents dans la liste ghost_compo_converter_ ... "
4512 <<
"Ca signifie peut-etre que l'on vient de creer des bulles qui existaient deja..."
4513 <<
"C'est etrange... on s'arrete... ";
4547 const DoubleTab& bounding_box_offsetp,
4548 const DoubleTab& bounding_box_offsetm,
4549 const DoubleTab& authorized_bounding_box,
4550 ArrOfInt& masque_duplicata_pour_compo)
4557 for (
int icompo = 0; icompo < nbulles; icompo++)
4559 int masque_sortie_domaine_reel = 0;
4565 for (
int direction = 0; direction < 3; direction++)
4570 if (bounding_box(icompo, direction, 0) < authorized_bounding_box(direction, 0))
4572 masque_sortie_domaine_reel |= (1 << direction);
4573 masque_sortie_domaine_reel |= (16 << direction);
4580 if (bounding_box_offsetp(icompo, 0, 0) < authorized_bounding_box(0, 0))
4582 masque_sortie_domaine_reel |= (1 << 3);
4583 masque_sortie_domaine_reel |= (16 << 3);
4585 if (bounding_box_offsetp(icompo, 0, 1) > authorized_bounding_box(0, 1))
4587 masque_sortie_domaine_reel |= (1 << 3);
4593 if (bounding_box(icompo, direction, 1) > authorized_bounding_box(direction, 1))
4597 masque_sortie_domaine_reel |= (1 << direction);
4605 if (bounding_box_offsetm(icompo, 0, 0) < authorized_bounding_box(0, 0))
4607 masque_sortie_domaine_reel |= (1 << 3);
4608 masque_sortie_domaine_reel |= (16 << 3);
4610 if (bounding_box_offsetm(icompo, 0, 1) > authorized_bounding_box(0, 1))
4612 masque_sortie_domaine_reel |= (1 << 3);
4622 masque_duplicata_pour_compo[icompo] = masque_sortie_domaine_reel;
4627 envoyer_broadcast(masque_duplicata_pour_compo, 0);
4632 const DoubleTab& authorized_bounding_box,
4633 ArrOfInt& masque_duplicata_pour_compo)
4640 for (
int icompo = 0; icompo < nbulles; icompo++)
4642 int masque_sortie_domaine_reel = 0;
4648 for (
int direction = 0; direction < 3; direction++)
4653 if (bounding_box(icompo, direction, 0) < authorized_bounding_box(direction, 0))
4655 masque_sortie_domaine_reel |= (1 << direction);
4656 masque_sortie_domaine_reel |= (8 << direction);
4659 if (bounding_box(icompo, direction, 1) > authorized_bounding_box(direction, 1))
4661 masque_sortie_domaine_reel |= (1 << direction);
4667 masque_duplicata_pour_compo[icompo] = masque_sortie_domaine_reel;
4673 envoyer_broadcast(masque_duplicata_pour_compo, 0);
4686 ArrOfInt liste_facettes_pour_suppression;
4687 for (
int iface = 0; iface < nbfacettes; iface++)
4689 const int icompo = compo_connex[iface];
4709 DoubleTab bounding_box;
4727 for (
int icompo = 0; icompo < nbulles; icompo++)
4729 const int masque_sortie = masque[icompo];
4730 if (masque_sortie & 2)
4734 if (masque_sortie & 16)
4765 DoubleTab deplacement;
4767 ArrOfDouble volume_reel;
4770 calculer_deplacement_from_masque_in_array(mesh, deplacement, masque_deplacement_par_compo,
bounding_box_NS_domain_, position, nbulles);
4776 for (
int ib = 0; ib < nbulles; ib++)
4778 const int code_deplacement = masque_deplacement_par_compo[ib];
4779 if (code_deplacement != 0)
4781 Journal() <<
"IJK_Interfaces::deplacer_bulle_perio : Un deplacement a eu "
4782 <<
"lieu pour la composante " << ib << finl;
4783 Journal() <<
"Deplacement x y z : ";
4784 for (
int dir = 0; dir < 3; dir++)
4786 int unused_variable = 0 ;
4787 const int decode = decoder_deplacement(code_deplacement, dir, unused_variable);
4820 statistics().create_custom_counter(
"Calcul rho mu indicatrice: calcul de l'indicatrice",2,
"IJK");
4821 statistics().begin_count(
"Calcul rho mu indicatrice: calcul de l'indicatrice",statistics().get_last_opened_counter_level()+1);
4824 const IntTab& elem_faces = domaine_vf.
elem_faces();
4825 const IntTab& faces_voisins = domaine_vf.
face_voisins();
4831 const int ni = indic.
ni();
4832 const int nj = indic.
nj();
4833 const int nk = indic.
nk();
4841 const ArrOfInt& index_elem = intersec.
index_elem();
4844 for (
int k = 0; k < nk; k++)
4846 for (
int j = 0; j < nj; j++)
4848 for (
int i = 0; i < ni; i++)
4855 int index = index_elem[num_elem];
4856 double somme_contrib = 0.;
4865 int nb_increment_somme_contrib = 0;
4866 while (somme_contrib > 1.)
4868 somme_contrib -= 1.;
4869 nb_increment_somme_contrib++;
4871 while (somme_contrib < 0.)
4873 somme_contrib += 1.;
4874 nb_increment_somme_contrib++;
4877 if (nb_increment_somme_contrib > 1)
4879 Cerr <<
"nb_increment_somme_contrib= " << nb_increment_somme_contrib << finl;
4891 if (somme_contrib > 0.)
4896 indic(i, j, k) = somme_contrib;
4904 statistics().create_custom_counter(
"Calcul de l'indicatrice : recherche compo connexes",2,
"IJK");
4905 statistics().begin_count(
"Calcul de l'indicatrice : recherche compo connexes",statistics().get_last_opened_counter_level()+1);
4908 int nb_compo_locales = search_connex_components_local(elem_faces, faces_voisins,
num_compo_);
4915 statistics().end_count(
"Calcul de l'indicatrice : recherche compo connexes");
4922 Cerr <<
"GB-check-interf nb_compo_in_num_compo_ - (nb_reelles + nb_ghost)= " <<
nb_compo_in_num_compo_ <<
" - "
4926 Cerr <<
"IJK_Interfaces::calculer_indicatrice Evaluations drapeau vapeur " << drapeau << finl;
4929#define VERIF_DRAPEAU 0
4944 compute_drapeaux_vapeur_v1(mesh, elem_faces, faces_voisins, num_compo,
true, drapeau1);
4945 compute_drapeaux_vapeur_v1(mesh, elem_faces, faces_voisins, num_compo,
false, drapeau1b);
4954 Cerr <<
"IJK_Interfaces::calculer_indicatrice Evaluations drapeaux" << finl;
4956 Cerr <<
"drapeau v1 towards_liquid=true : " << drapeau1 << finl;
4957 Cerr <<
"drapeau v1 towards_liquid=false : " << drapeau1b << finl;
4958 Cerr <<
"drapeau v2 : " << drapeau2 << finl;
4959 Cerr <<
"drapeau v3 : " << drapeau3 << finl;
4960 Cerr <<
"drapeau v4 : " << drapeau << finl;
4964 int sum = drapeau[idx] + drapeau1[idx] + drapeau1b[idx] + drapeau2[idx] + drapeau3[idx];
4965 if ((drapeau[idx] != drapeau1[idx]) || (drapeau[idx] != drapeau1b[idx]) || (drapeau[idx] != drapeau2[idx]) ||
4966 (drapeau[idx] != drapeau3[idx]))
4967 Cerr <<
"Les methodes different... "
4968 <<
"valeurs donnees par les methodes : " << drapeau[idx] << drapeau1[idx] << drapeau1b[idx]
4969 << drapeau2[idx] << drapeau3[idx] <<
" pour idx = " << idx << finl;
4971 if ((sum != 0) && (sum != 5))
4973 Cerr <<
"Les methodes different... "
4974 <<
"sum = " << sum <<
" pour idx = " << idx << finl;
4978 assert(drapeau == drapeau1);
4979 assert(drapeau == drapeau1b);
4980 assert(drapeau == drapeau2);
4981 assert(drapeau == drapeau3);
4984 check_somme_drapeau(drapeau);
4985 check_somme_drapeau(drapeau1);
4986 check_somme_drapeau(drapeau1b);
4987 check_somme_drapeau(drapeau2);
4988 check_somme_drapeau(drapeau3);
4992#define DEBUG_INDIC 0
4997 for (
int k = 0; k < nk; k++)
4999 for (
int j = 0; j < nj; j++)
5001 for (
int i = 0; i < ni; i++)
5004 int compo = num_compo[num_elem];
5006 indic(i, j, k) = compo;
5015 for (
int k = 0; k < nk; k++)
5017 for (
int j = 0; j < nj; j++)
5019 for (
int i = 0; i < ni; i++)
5024 if (compo >= 0 && drapeau[compo] == 0)
5027 indic(i, j, k) = 0.;
5034 statistics().end_count(
"Calcul rho mu indicatrice: calcul de l'indicatrice");
5039 statistics().create_custom_counter(
"Calcul rho mu indicatrice: calcul de l'indicatrice",2,
"IJK");
5040 statistics().begin_count(
"Calcul rho mu indicatrice: calcul de l'indicatrice",statistics().get_last_opened_counter_level()+1);
5045 const int ni = indic.
ni();
5046 const int nj = indic.
nj();
5047 const int nk = indic.
nk();
5052 const ArrOfInt& index_elem = intersec.
index_elem();
5054 for (
int k = 0; k < nk; k++)
5056 for (
int j = 0; j < nj; j++)
5058 for (
int i = 0; i < ni; i++)
5061 int index = index_elem[num_elem];
5063 if (0. < indic(i, j, k) && indic(i, j, k) < 1.)
5064 indic(i, j, k) = 2.;
5066 double somme_contrib = 0.;
5075 int nb_increment_somme_contrib = 0;
5076 while (somme_contrib > 1.)
5078 somme_contrib -= 1.;
5079 nb_increment_somme_contrib++;
5081 while (somme_contrib < 0.)
5083 somme_contrib += 1.;
5084 nb_increment_somme_contrib++;
5087 if (nb_increment_somme_contrib > 1)
5088 Process::exit(
"Error in IJK_Interfaces::calculer_indicatrice_optim !");
5090 if (somme_contrib > 0.)
5094 indic(i, j, k) = somme_contrib;
5102 statistics().end_count(
"Calcul rho mu indicatrice: calcul de l'indicatrice");
5107 statistics().create_custom_counter(
"Calcul rho mu indicatrice: calcul des indicatrices",2,
"IJK");
5108 statistics().begin_count(
"Calcul rho mu indicatrice: calcul des indicatrices",statistics().get_last_opened_counter_level()+1);
5110 const Domaine& domaine = domaine_vf.
domaine();
5111 const IntTab& elem_faces = domaine_vf.
elem_faces();
5112 const IntTab& faces_voisins = domaine_vf.
face_voisins();
5116 for (
int igroup = 0; igroup <
nb_groups_; igroup++)
5118 indic[igroup].data() = 1.;
5122 Cerr <<
"This choice of options has to be validated. "
5123 <<
"It seems you are using a multi-groups of bubbles "
5124 <<
"and the forcing method with the color_function "
5125 <<
"which has not been tested yet!"
5130 IntVect num_compo[max_authorized_nb_of_groups_];
5135 domaine.creer_tableau_elements(num_compo[i]);
5137 const int ni = indic[0].ni();
5138 const int nj = indic[0].nj();
5139 const int nk = indic[0].nk();
5141 assert(
nb_groups_ <= max_authorized_nb_of_groups_);
5150 const ArrOfInt& index_elem = intersec.
index_elem();
5153 for (
int k = 0; k < nk; k++)
5155 for (
int j = 0; j < nj; j++)
5157 for (
int i = 0; i < ni; i++)
5164 int index = index_elem[num_elem];
5165 double somme_contrib[max_authorized_nb_of_groups_] = {0.};
5171 int icompo = compo_connexe[facette];
5177 m = decoder_numero_bulle(-icompo);
5185 for (
int igroup = 0; igroup <
nb_groups_; igroup++)
5187 while (somme_contrib[igroup] > 1.)
5188 somme_contrib[igroup] -= 1.;
5189 while (somme_contrib[igroup] < 0.)
5190 somme_contrib[igroup] += 1.;
5191 if (somme_contrib[igroup] > 0.)
5195 indic[igroup](i, j, k) = somme_contrib[igroup];
5196 num_compo[igroup][num_elem] = -1;
5204 ArrOfInt drapeau[max_authorized_nb_of_groups_];
5205 statistics().create_custom_counter(
"Calcul de l'indicatrice : recherche compo connexes",2,
"IJK");
5208 statistics().begin_count(
"Calcul de l'indicatrice : recherche compo connexes",statistics().get_last_opened_counter_level()+1);
5211 int nb_compo_locales = search_connex_components_local(elem_faces, faces_voisins, num_compo[i]);
5214 statistics().end_count(
"Calcul de l'indicatrice : recherche compo connexes");
5224 for (
int k = 0; k < nk; k++)
5226 for (
int j = 0; j < nj; j++)
5228 for (
int i = 0; i < ni; i++)
5231 for (
int igroup = 0; igroup <
nb_groups_; igroup++)
5233 int compo = num_compo[igroup][num_elem];
5234 if (compo >= 0 && drapeau[igroup][compo] == 0)
5237 indic[igroup](i, j, k) = 0.;
5245 statistics().end_count(
"Calcul rho mu indicatrice: calcul des indicatrices");
5250 statistics().create_custom_counter(
"Calcul rho mu indicatrice: calcul des indicatrices",2,
"IJK");
5251 statistics().begin_count(
"Calcul rho mu indicatrice: calcul des indicatrices",statistics().get_last_opened_counter_level()+1);
5256 const int ni = indic[0].ni();
5257 const int nj = indic[0].nj();
5258 const int nk = indic[0].nk();
5260 assert(
nb_groups_ <= max_authorized_nb_of_groups_);
5266 const ArrOfInt& index_elem = intersec.
index_elem();
5269 for (
int k = 0; k < nk; k++)
5271 for (
int j = 0; j < nj; j++)
5273 for (
int i = 0; i < ni; i++)
5276 int index = index_elem[num_elem];
5277 double somme_contrib[max_authorized_nb_of_groups_] = {0.};
5279 for (
int igroup = 0; igroup <
nb_groups_; igroup++)
5281 if (0. < indic[igroup](i, j, k) && indic[igroup](i, j, k) < 1.)
5282 indic[igroup](i, j, k) = 2.;
5290 int icompo = compo_connexe[facette];
5296 m = decoder_numero_bulle(-icompo);
5303 for (
int igroup = 0; igroup <
nb_groups_; igroup++)
5305 while (somme_contrib[igroup] > 1.)
5306 somme_contrib[igroup] -= 1.;
5307 while (somme_contrib[igroup] < 0.)
5308 somme_contrib[igroup] += 1.;
5309 if (somme_contrib[igroup] > 0.)
5311 indic[igroup](i, j, k) = somme_contrib[igroup];
5319 for (
int igroup = 0; igroup <
nb_groups_; igroup++)
5323 statistics().end_count(
"Calcul rho mu indicatrice: calcul des indicatrices");
5328 const int ni = indic.
ni();
5329 const int nj = indic.
nj();
5330 const int nk = indic.
nk();
5350 for (
int k = 0; k < nk; k++)
5352 for (
int j = 0; j < nj; j++)
5354 for (
int i = 0; i < ni; i++)
5356 if (indic(i, j, k) == 2.)
5358 ArrOfInt voisins_num(6);
5366 ArrOfDouble voisins_valeur(6);
5367 voisins_valeur[0] = imin + i == 0 ? BORD : indic(i - 1, j, k);
5368 voisins_valeur[1] = jmin + j == 0 ? BORD : indic(i, j - 1, k);
5369 voisins_valeur[2] = kmin + k == 0 ? BORD : indic(i, j, k - 1);
5370 voisins_valeur[3] = imin + i == nitot ? BORD : indic(i + 1, j, k);
5371 voisins_valeur[4] = jmin + j == njtot ? BORD : indic(i, j + 1, k);
5372 voisins_valeur[5] = kmin + k == nktot ? BORD : indic(i, j, k + 1);
5379 for (
int v = 0; v < 6; v++)
5381 double valeur_indic_voisin = voisins_valeur[v];
5382 if (valeur_indic_voisin == 0. || valeur_indic_voisin == 1.)
5384 indic(i, j, k) = valeur_indic_voisin;
5396 for (
int v = 0; v < 6; v++)
5398 double valeur_indic_voisin = voisins_valeur[v];
5399 const int direction = v % 3;
5400 const int face_plus = (v > 2) ? 1 : -1;
5402 if (0. < valeur_indic_voisin && valeur_indic_voisin < 1.)
5404 int num_elem = voisins_num[v];
5405 if (num_elem != BORD)
5410 if (indic(i, j, k) != -1)
5442 statistics().create_custom_counter(
"Calcul rho mu indicatrice: calcul de l'indicatrice surface face",2,
"IJK");
5443 statistics().begin_count(
"Calcul rho mu indicatrice: calcul de l'indicatrice surface face",statistics().get_last_opened_counter_level()+1);
5448 const int ni = indic_surfacique_face[0].ni();
5449 const int nj = indic_surfacique_face[0].nj();
5450 const int nk = indic_surfacique_face[0].nk();
5454 for (
int k = 0; k < nk; k++)
5456 for (
int j = 0; j < nj; j++)
5458 for (
int i = 0; i < ni; i++)
5463 indic_surfacique_face[0](i, j, k) = indic(i, j, k);
5464 indic_surfacique_face[1](i, j, k) = indic(i, j, k);
5465 indic_surfacique_face[2](i, j, k) = indic(i, j, k);
5471 indic_surfacique_face[0](i, j, k) = norme[0](i, j, k) == 0. ? 0. : (norme[0](i, j, k) > 0. ? 0. : 1.);
5472 indic_surfacique_face[1](i, j, k) = norme[1](i, j, k) == 0. ? 0. : (norme[1](i, j, k) > 0. ? 0. : 1.);
5473 indic_surfacique_face[2](i, j, k) = norme[2](i, j, k) == 0. ? 0. : (norme[2](i, j, k) > 0. ? 0. : 1.);
5479 indic_surfacique_face[0](i, j, k) = indic(i-1, j, k);
5483 indic_surfacique_face[1](i, j, k) = indic(i, j-1, k);
5487 indic_surfacique_face[2](i, j, k) = indic(i, j, k-1);
5492 baric_face[0][0](i, j, k) = 0.5;
5493 baric_face[0][1](i, j, k) = 0.5;
5494 baric_face[1][0](i, j, k) = 0.5;
5495 baric_face[1][1](i, j, k) = 0.5;
5496 baric_face[2][0](i, j, k) = 0.5;
5497 baric_face[2][1](i, j, k) = 0.5;
5506 const ArrOfInt& index_elem = intersec.
index_elem();
5509 for (
int k = 0; k < nk; k++)
5511 for (
int j = 0; j < nj; j++)
5513 for (
int i = 0; i < ni; i++)
5525 int index = index_elem[num_elem];
5526 double somme_contrib[3] = {0., 0., 0.};
5527 double somme_contrib_baryc[3][2] = {{0., 0.}, {0., 0.}, {0., 0.}};
5546 for (
int dir=0; dir<3; dir++)
5548 const int dir_i = (dir == 0);
5549 const int dir_j = (dir == 1);
5550 const int dir_k = (dir == 2);
5552 if (
est_pure(indic(i-dir_i,j-dir_j,k-dir_k)))
5562 if (somme_contrib[dir]*somme_contrib[dir] > 0)
5564 while (somme_contrib[dir] > 1.)
5566 somme_contrib[dir] -= 1.;
5567 somme_contrib_baryc[dir][0] -= 1./2.;
5568 somme_contrib_baryc[dir][1] -= 1./2.;
5570 while (somme_contrib[dir] < 0.)
5572 somme_contrib[dir] += 1.;
5573 somme_contrib_baryc[dir][0] += 1./2.;
5574 somme_contrib_baryc[dir][1] += 1./2.;
5577 indic_surfacique_face[dir](i, j, k) = somme_contrib[dir];
5581 baric_face[dir][0](i, j, k) = 1./2.;
5582 baric_face[dir][1](i, j, k) = 1./2.;
5586 baric_face[dir][0](i, j, k) = somme_contrib_baryc[dir][0]/abs(somme_contrib[dir]);
5587 baric_face[dir][1](i, j, k) = somme_contrib_baryc[dir][1]/abs(somme_contrib[dir]);
5591 if ((baric_face[dir][0](i, j, k) <= 0) && (baric_face[dir][0](i, j, k) >= -1e-12 || somme_contrib[dir] < 1e-8))
5593 baric_face[dir][0](i, j, k) = 1e-12;
5595 else if ((baric_face[dir][0](i, j, k) >= 1) && (baric_face[dir][0](i, j, k) <= 1+1e-12 || somme_contrib[dir] < 1e-8))
5597 baric_face[dir][0](i, j, k) = 1 - 1e-12;
5599 if ((baric_face[dir][1](i, j, k) <= 0) && (baric_face[dir][1](i, j, k) >= -1e-12 || somme_contrib[dir] < 1e-8))
5601 baric_face[dir][1](i, j, k) = 1e-12;
5603 else if ((baric_face[dir][1](i, j, k) >= 1) && (baric_face[dir][1](i, j, k) <= 1+1e-12 || somme_contrib[dir] < 1e-8))
5605 baric_face[dir][1](i, j, k) = 1 - 1e-12;
5608 assert((baric_face[dir][0](i, j, k) >= 0) && (baric_face[dir][0](i, j, k) <= 1));
5609 assert((baric_face[dir][1](i, j, k) >= 0) && (baric_face[dir][1](i, j, k) <= 1));
5610 assert(((indic_surfacique_face[dir](i, j, k) < 0.) && (baric_face[dir][0](i, j, k) < 0.)) || ((indic_surfacique_face[dir](i, j, k) > 0.) && (baric_face[dir][0](i, j, k) > 0.)));
5611 assert(((indic_surfacique_face[dir](i, j, k) < 0.) && (baric_face[dir][1](i, j, k) < 0.)) || ((indic_surfacique_face[dir](i, j, k) > 0.) && (baric_face[dir][1](i, j, k) > 0.)));
5621 statistics().end_count(
"Calcul rho mu indicatrice: calcul de l'indicatrice surface face");
5627 statistics().create_custom_counter(
"Calcul rho mu indicatrice: calcul de l'indicatrice surface face",2,
"IJK");
5628 statistics().begin_count(
"Calcul rho mu indicatrice: calcul de l'indicatrice surface face",statistics().get_last_opened_counter_level()+1);
5633 const int ni = indic_surfacique_face[0].ni();
5634 const int nj = indic_surfacique_face[0].nj();
5635 const int nk = indic_surfacique_face[0].nk();
5639 for (
int k = 0; k < nk; k++)
5641 for (
int j = 0; j < nj; j++)
5643 for (
int i = 0; i < ni; i++)
5648 indic_surfacique_face[0](i, j, k) = indic(i, j, k);
5649 indic_surfacique_face[1](i, j, k) = indic(i, j, k);
5650 indic_surfacique_face[2](i, j, k) = indic(i, j, k);
5656 indic_surfacique_face[0](i, j, k) = norme[0](i, j, k) == 0. ? 0. : (norme[0](i, j, k) > 0. ? 0. : 1.);
5657 indic_surfacique_face[1](i, j, k) = norme[1](i, j, k) == 0. ? 0. : (norme[1](i, j, k) > 0. ? 0. : 1.);
5658 indic_surfacique_face[2](i, j, k) = norme[2](i, j, k) == 0. ? 0. : (norme[2](i, j, k) > 0. ? 0. : 1.);
5664 indic_surfacique_face[0](i, j, k) = indic(i-1, j, k);
5668 indic_surfacique_face[1](i, j, k) = indic(i, j-1, k);
5672 indic_surfacique_face[2](i, j, k) = indic(i, j, k-1);
5683 const ArrOfInt& index_elem = intersec.
index_elem();
5686 for (
int k = 0; k < nk; k++)
5688 for (
int j = 0; j < nj; j++)
5690 for (
int i = 0; i < ni; i++)
5702 int index = index_elem[num_elem];
5703 double somme_contrib[3] = {0., 0., 0.};
5715 for (
int dir=0; dir<3; dir++)
5717 const int dir_i = (dir == 0);
5718 const int dir_j = (dir == 1);
5719 const int dir_k = (dir == 2);
5721 if (
est_pure(indic(i-dir_i,j-dir_j,k-dir_k)))
5731 if (somme_contrib[dir]*somme_contrib[dir] > 0)
5733 while (somme_contrib[dir] > 1.)
5735 somme_contrib[dir] -= 1.;
5737 while (somme_contrib[dir] < 0.)
5739 somme_contrib[dir] += 1.;
5742 indic_surfacique_face[dir](i, j, k) = somme_contrib[dir];
5751 statistics().end_count(
"Calcul rho mu indicatrice: calcul de l'indicatrice surface face");
5757 statistics().create_custom_counter(
"Calcul rho mu indicatrice: calcul de la surface interfaciale",2,
"IJK");
5758 statistics().begin_count(
"Calcul rho mu indicatrice: calcul de la surface interfaciale",statistics().get_last_opened_counter_level()+1);
5762 const ArrOfDouble& surface_facettes =
maillage_ft_ijk_.get_update_surface_facettes();
5764 const int ni = surf_interface.
ni();
5765 const int nj = surf_interface.
nj();
5766 const int nk = surf_interface.
nk();
5770 for (
int k = 0; k < nk; k++)
5772 for (
int j = 0; j < nj; j++)
5774 for (
int i = 0; i < ni; i++)
5776 surf_interface(i, j, k) = 0.;
5785 const ArrOfInt& index_elem = intersec.
index_elem();
5788 for (
int k = 0; k < nk; k++)
5790 for (
int j = 0; j < nj; j++)
5792 for (
int i = 0; i < ni; i++)
5803 int index = index_elem[num_elem];
5804 double somme_contrib_surf = 0.;
5815 surf_interface(i, j, k) = somme_contrib_surf;
5816 assert(somme_contrib_surf >= 0.);
5822 statistics().end_count(
"Calcul rho mu indicatrice: calcul de la surface interfaciale");
5827 statistics().create_custom_counter(
"Calcul rho mu indicatrice: calcul du barycentre",2,
"IJK");
5828 statistics().begin_count(
"Calcul rho mu indicatrice: calcul du barycentre",statistics().get_last_opened_counter_level()+1);
5832 const int ni = baric[0].ni();
5833 const int nj = baric[0].nj();
5834 const int nk = baric[0].nk();
5838 for (
int k = 0; k < nk; k++)
5840 for (
int j = 0; j < nj; j++)
5842 for (
int i = 0; i < ni; i++)
5845 baric[0](i, j, k) = 0.5;
5846 baric[1](i, j, k) = 0.5;
5847 baric[2](i, j, k) = 0.5;
5856 const ArrOfInt& index_elem = intersec.
index_elem();
5859 for (
int k = 0; k < nk; k++)
5861 for (
int j = 0; j < nj; j++)
5863 for (
int i = 0; i < ni; i++)
5874 int index = index_elem[num_elem];
5875 double somme_contrib = 0.;
5876 double somme_contrib_baryc[3] = {0., 0., 0.};
5891 while (somme_contrib > 1.)
5893 somme_contrib -= 1.;
5894 somme_contrib_baryc[0] -= 1./2.;
5895 somme_contrib_baryc[1] -= 1./2.;
5896 somme_contrib_baryc[2] -= 1./2.;
5899 while (somme_contrib < 0.)
5901 somme_contrib += 1.;
5902 somme_contrib_baryc[0] += 1./2.;
5903 somme_contrib_baryc[1] += 1./2.;
5904 somme_contrib_baryc[2] += 1./2.;
5909 assert(indic(i, j, k) == somme_contrib);
5911 baric[0](i, j, k) = somme_contrib_baryc[0]/abs(somme_contrib);
5912 baric[1](i, j, k) = somme_contrib_baryc[1]/abs(somme_contrib);
5913 baric[2](i, j, k) = somme_contrib_baryc[2]/abs(somme_contrib);
5917 if ((baric[0](i, j, k) <= 0) && (baric[0](i, j, k) >= -1e-12 || somme_contrib < 1e-8))
5919 baric[0](i, j, k) = 1e-12;
5921 else if ((baric[0](i, j, k) >= 1) && (baric[0](i, j, k) <= 1+1e-12 || somme_contrib < 1e-8))
5923 baric[0](i, j, k) = 1 - 1e-12;
5925 if ((baric[1](i, j, k) <= 0) && (baric[1](i, j, k) >= -1e-12 || somme_contrib < 1e-8))
5927 baric[1](i, j, k) = 1e-12;
5929 else if ((baric[1](i, j, k) >= 1) && (baric[1](i, j, k) <= 1+1e-12 || somme_contrib < 1e-8))
5931 baric[1](i, j, k) = 1 - 1e-12;
5933 if ((baric[2](i, j, k) <= 0) && (baric[2](i, j, k) >= -1e-12 || somme_contrib < 1e-8))
5935 baric[2](i, j, k) = 1e-12;
5937 else if ((baric[2](i, j, k) >= 1) && (baric[2](i, j, k) <= 1+1e-12 || somme_contrib < 1e-8))
5939 baric[2](i, j, k) = 1 - 1e-12;
5942 assert((baric[0](i, j, k) >= 0) && (baric[0](i, j, k) <= 1));
5943 assert((baric[1](i, j, k) >= 0) && (baric[1](i, j, k) <= 1));
5944 assert((baric[2](i, j, k) >= 0) && (baric[2](i, j, k) <= 1));
5945 assert(((indic(i, j, k) < 0.) && (baric[0](i, j, k) < 0.)) || ((indic(i, j, k) > 0.) && (baric[0](i, j, k) > 0.)));
5946 assert(((indic(i, j, k) < 0.) && (baric[1](i, j, k) < 0.)) || ((indic(i, j, k) > 0.) && (baric[1](i, j, k) > 0.)));
5947 assert(((indic(i, j, k) < 0.) && (baric[2](i, j, k) < 0.)) || ((indic(i, j, k) > 0.) && (baric[2](i, j, k) > 0.)));
5953 statistics().end_count(
"Calcul rho mu indicatrice: calcul du barycentre");
5964 const Domaine& domaine = domaine_vf.
domaine();
5967 const int nb_elem = domaine.nb_elem();
5968 const int nb_elem_tot = domaine.nb_elem_tot();
5969 Cerr <<
"nb_elem= " << nb_elem <<
" nb_elem_tot= " << nb_elem_tot;
5984#define NUM_COMPO_INVALID (-2000000000)
5992 IJK_Field_vector3_double& vpoint,
5993 IJK_Field_vector3_double& vrepul,
5994 IJK_Field_vector3_double& vabsrepul
5997 statistics().begin_count(STD_COUNTERS::source_terms,statistics().get_last_opened_counter_level()+1);
6009 double interf1, interf2 ;
6011 Int3 ijk_other_elem;
6012 for (
int icol1 = 0; icol1 < max_authorized_nb_of_components_; icol1++)
6014 for (
int direction = 0; direction < 3; direction++)
6016 for (
int k = 0; k < vpoint[direction].nk(); k++)
6018 for (
int j = 0; j < vpoint[direction].nj(); j++)
6020 for (
int i = 0; i < vpoint[direction].ni(); i++)
6027 ijk_other_elem[0] = i;
6028 ijk_other_elem[1] = j;
6029 ijk_other_elem[2] = k;
6030 ijk_other_elem[direction] -= 1;
6032 s2 =
surf_par_compo_[
old()][icol1](ijk_other_elem[0],ijk_other_elem[1],ijk_other_elem[2]);
6035 vpoint[direction](i, j, k) += 1./2. * (s1 * interf1 + s2 * interf2);
6048 const int nkmax = std::max(vpoint[DIRECTION_I].nk(), std::max(vpoint[DIRECTION_J].nk(), vpoint[DIRECTION_K].nk()));
6049 for (
int k = 0; k < nkmax; k++)
6051 for (
int direction = 0; direction < 3; direction++)
6053 if (k >= vpoint[direction].nk())
6061 for (
int j = 0; j < vpoint[direction].nj(); j++)
6063 for (
int i = 0; i < vpoint[direction].ni(); i++)
6068 const int global_face_position = ijk_face[direction] + offset;
6069 if (!perio && (global_face_position == 0 || global_face_position == nb_items_tot))
6072 Int3 ijk_droite = ijk_face;
6073 Int3 ijk_gauche = ijk_face;
6074 ijk_gauche[direction]--;
6078 for (
int gauche_droite = 0; gauche_droite <= 1; gauche_droite++)
6081 for (
int icol1 = 0; icol1 < max_authorized_nb_of_components_; icol1++)
6083 const Int3& elem1 = gauche_droite ? ijk_droite : ijk_gauche;
6084 const Int3& elem2 = gauche_droite ? ijk_gauche : ijk_droite;
6086 int signe = gauche_droite ? -1 : 1;
6088 if (num_compo == NUM_COMPO_INVALID)
6093 for (icol2 = 0; icol2 < max_authorized_nb_of_components_; icol2++)
6097 if (icol2 < max_authorized_nb_of_components_ && gauche_droite == 1)
6110 elem1[0], elem1[1], elem1[2]
6116 double indic_voisin = 0., phi_voisin = 0.;
6117 double repul_voisin = 0.;
6118 if (icol2 == max_authorized_nb_of_components_)
6149 for (
int dir = 0; dir < 3; dir++)
6151 const int idx = icol1 * 3 + dir;
6158 bary_facettes_dans_elem[dir] =
bary_par_compo_[
old()][idx](elem1[0], elem1[1], elem1[2]);
6165 indic_voisin = (ps > 0 ? 1. : 0.);
6168 repul_voisin = repul;
6182 double gradient_phi_indic = (phi_voisin * indic_voisin - phi * indic) / delta_dir * signe;
6184 vpoint[direction](i, j, k) += gradient_phi_indic;
6185 double gradient_repul_indic = (repul_voisin * indic_voisin - repul * indic) / delta_dir * signe;
6186 vrepul[direction](i, j, k) += gradient_repul_indic;
6187 vabsrepul[direction](i, j, k) += fabs(gradient_repul_indic);
6191 double indic_voisin = 0., phi_face = 0.;
6192 double repul_face = 0.;
6193 if (icol2 == max_authorized_nb_of_components_)
6224 for (
int dir = 0; dir < 3; dir++)
6226 const int idx = icol1 * 3 + dir;
6229 bary_facettes_dans_elem[dir] =
bary_par_compo_[
old()][idx](elem1[0], elem1[1], elem1[2]);
6236 indic_voisin = (ps > 0 ? 1. : 0.);
6254 elem2[0], elem2[1], elem2[2]
6259 assert(indic_voisin > -0.5);
6260 if (est_egal(surface+surface_voisin,0.))
6268 phi_face = (phi * surface + phi_voisin * surface_voisin) / (surface + surface_voisin);
6269 repul_face = (repul * surface + repul_voisin * surface_voisin) / (surface + surface_voisin);
6273 double gradient_indic = (indic_voisin - indic) / delta_dir * signe;
6283 vpoint[direction](i, j, k) += phi_face * gradient_indic ;
6284 vrepul[direction](i, j, k) += repul_face * gradient_indic;
6285 vabsrepul[direction](i, j, k) += fabs(repul_face * gradient_indic);
6294 statistics().end_count(STD_COUNTERS::source_terms);
6306static bool intersection_segment_triangle(
const Vecteur3& A,
6320 const double det1 = determinant(CD, CE, CA);
6321 const double det2 = determinant(CD, CE, CB);
6323 if (det1 * det2 >= 0)
6333 const double det3 = determinant(AC, AD, AB);
6334 const double det4 = determinant(AD, AE, AB);
6335 const double det5 = determinant(AE, AC, AB);
6337 if ((det3 * det4 <= 0) || (det3 * det5 <= 0))
6349 const IntTab& facettes = mesh.
facettes();
6350 const DoubleTab& sommets = mesh.
sommets();
6355 ArrOfInt facettes_traversantes;
6357 const int N = facettes_traversantes.
size_array();
6366 int fa7 = facettes_traversantes[index];
6386 for (
int i = 0; i < 3; i++)
6387 sommets_facette[i] =
Vecteur3(sommets, facettes(fa7, i));
6395 for (
int isom = 0; isom < 3; isom++)
6398 for (
int dir = 0; dir < 3; dir++)
6399 coordA[dir] += bary_som * sommets_facette[isom][dir];
6406 coordB[direction] += cell_size[direction] * face_plus;
6411 for (
int index2 = intersections.
index_elem()[num_elem]; index2 >= 0;
6419 const Vecteur3 s0(sommets, facettes(nvl_fa7, 0));
6420 const Vecteur3 s1(sommets, facettes(nvl_fa7, 1));
6421 const Vecteur3 s2(sommets, facettes(nvl_fa7, 2));
6424 coupe = intersection_segment_triangle(coordA, coordB, s0, s1, s2);
6430 for (
int j = 0; j < N; j++)
6432 if (facettes_traversantes[j] == nvl_fa7)
6447 if (normale_facette[direction] * face_plus > 0.)
6460#define COMPUTE_DRAPEAUX_VALIDATION
6463 ArrOfInt& drapeau_vapeur)
const
6465 statistics().create_custom_counter(
"Calcul de l'indicatrice: calculs des drapeaux",2,
"IJK");
6466 statistics().begin_count(
"Calcul de l'indicatrice: calculs des drapeaux",statistics().get_last_opened_counter_level()+1);
6471 const ArrOfInt& index_elem =
maillage_ft_ijk_.intersections_elem_facettes().index_elem();
6472 const ArrOfInt& index_facette =
maillage_ft_ijk_.intersections_elem_facettes().index_facette();
6475 const IntTab& elem_faces = domaine_vf.
elem_faces();
6476 const IntTab& faces_voisins = domaine_vf.
face_voisins();
6478#ifdef COMPUTE_DRAPEAUX_VALIDATION
6479 ArrOfInt composantes_comptage_phase0(drapeau_vapeur.
size_array());
6480 ArrOfInt composantes_comptage_phase1(drapeau_vapeur.
size_array());
6489 const IntTab& facettes = mesh.
facettes();
6490 const DoubleTab& sommets = mesh.
sommets();
6491 for (
int fa7 = 0; fa7 < nb_facettes; fa7++)
6498 for (
int i = 0; i < 3; i++)
6499 sommets_facette[i] =
Vecteur3(sommets, facettes(fa7, i));
6501 int index = index_facette[fa7];
6529 for (
int iface = 0; iface < 6; iface++)
6531 const int direction = iface % 3;
6532 const int face_plus = (iface > 2) ? 1 : -1;
6534 const int num_face = elem_faces(elem, iface);
6535 const int elem_voisin = faces_voisins(num_face, 0) + faces_voisins(num_face, 1) - elem;
6537 if (elem_voisin < 0)
6543 const int global_compo_voisin = vecteur_composantes[elem_voisin];
6547 if (global_compo_voisin < 0)
6550#ifndef COMPUTE_DRAPEAUX_VALIDATION
6553 if (drapeau_vapeur[global_compo_voisin] == 1)
6558 sommets_facette[2] - sommets_facette[0],
6561#ifndef COMPUTE_DRAPEAUX_VALIDATION
6565 if (normale_facette[direction] * face_plus <= 0.)
6574 for (
int isom = 0; isom < 3; isom++)
6577 for (
int dir = 0; dir < 3; dir++)
6578 coordA[dir] += bary_som * sommets_facette[isom][dir];
6587 coordB[direction] += cell_size[direction] * face_plus;
6596 for (
int index2 = index_elem[elem]; index2 >= 0;
6604 const Vecteur3 s0(sommets, facettes(nvl_fa7, 0));
6605 const Vecteur3 s1(sommets, facettes(nvl_fa7, 1));
6606 const Vecteur3 s2(sommets, facettes(nvl_fa7, 2));
6609 coupe = intersection_segment_triangle(coordA, coordB, s0, s1, s2);
6620 if (normale_facette[direction] * face_plus > 0.)
6621 drapeau_vapeur[global_compo_voisin] = 1;
6622#ifdef COMPUTE_DRAPEAUX_VALIDATION
6624 if (normale_facette[direction] * face_plus > 0.)
6626 composantes_comptage_phase1[global_compo_voisin]++;
6634 if (normale_facette[direction] * face_plus < 0.)
6636 composantes_comptage_phase0[global_compo_voisin]++;
6656#ifdef COMPUTE_DRAPEAUX_VALIDATION
6658 const int n = composantes_comptage_phase0.
size_array();
6664 for (
int i = 0; i < n; i++)
6669 if ((composantes_comptage_phase0[i] > 0 && composantes_comptage_phase1[i] > 0) ||
6670 (composantes_comptage_phase0[i] == 0 && composantes_comptage_phase1[i] == 0))
6675 Cerr <<
"[WARNING-Indic-Vote] compo/#phase0/#phase1 ";
6678 <<
"/" << composantes_comptage_phase0[i]
6679 <<
"/" << composantes_comptage_phase1[i]
6683 if (composantes_comptage_phase0[i] > composantes_comptage_phase1[i])
6684 drapeau_vapeur[i] = 0;
6686 drapeau_vapeur[i] = 1;
6690 Cerr <<
"End." << finl;
6692 envoyer_broadcast(drapeau_vapeur, 0);
6696 statistics().end_count(
"Calcul de l'indicatrice: calculs des drapeaux");
6699static inline double norme_carre(
const Vecteur3& x)
6701 return x[0] * x[0] + x[1] * x[1] + x[2] * x[2];
6706 return Vecteur3(x * v[0], x * v[1], x * v[2]);
6710 return Vecteur3(v[0] + w[0], v[1] + w[1], v[2] + w[2]);
6713static inline double calculer_carre_distance_sommet_facette(
const Vecteur3& coord,
6717 const IntTab& facettes = mesh.
facettes();
6718 const DoubleTab& sommets = mesh.
sommets();
6719 Vecteur3 s0(sommets, facettes(num_facette, 0));
6720 Vecteur3 s1(sommets, facettes(num_facette, 1));
6721 Vecteur3 s2(sommets, facettes(num_facette, 2));
6728 for (
double x = 0.; x < 1.01; x += 0.2)
6729 for (
double y = 0.; y < 1. - x + 0.01; y += 0.2)
6731 Vecteur3 s = s0 + x * (s1 - s0) + y * (s2 - s0);
6732 d = std::min(d, norme_carre(s - coord));
6737static void fill_relative_velocity(
const DoubleTab& vinterp_tmp,
const DoubleTab& vinterp,
const IntTab& facettes,
int id_facette,
int som, DoubleTab& vr_to_other)
6739 const double un_tiers = 1. / 3.;
6740 if (id_facette == -1)
6744 for (
int idir = 0; idir < 3; idir++)
6745 vr_to_other(som, idir) = 0.;
6750 const int isom0 = facettes(id_facette, 0);
6751 const int isom1 = facettes(id_facette, 1);
6752 const int isom2 = facettes(id_facette, 2);
6753 for (
int idir = 0; idir < 3; idir++)
6757 const double velocity_me = vinterp_tmp(som, idir);
6758 const double velocity_other =
6759 un_tiers * (vinterp(isom0, idir) + vinterp(isom1, idir) + vinterp(isom2, idir));
6760 vr_to_other(som, idir) = velocity_me - velocity_other;
6768 const ArrOfInt& compo_connexe_sommets,
6769 const DoubleTab& vinterp_tmp,
6771 ArrOfDouble& distance,
6772 DoubleTab& vr_to_other,
6773 const double distmax)
6776 Octree_Double octree;
6777 const IntTab& facettes = mesh.
facettes();
6782 ArrOfInt liste_facettes;
6785 const int nb_som = sommets_a_tester.
dimension(0);
6786 distance.resize_array(nb_som, RESIZE_OPTIONS::NOCOPY_NOINIT);
6788 vr_to_other.
resize(nb_som, 3, RESIZE_OPTIONS::NOCOPY_NOINIT);
6790 assert(vinterp_tmp.
dimension(0) == nb_som);
6791 for (
int i = 0; i < nb_som; i++)
6793 Vecteur3 coord(sommets_a_tester, i);
6794 const int compo_connexe_som = compo_connexe_sommets[i];
6796 double dmin = distmax * distmax;
6799 octree.
search_elements_box(coord[0] - distmax, coord[1] - distmax, coord[2] - distmax, coord[0] + distmax,
6800 coord[1] + distmax, coord[2] + distmax, liste_facettes);
6804 const int nliste = liste_facettes.
size_array();
6805 int idx_facette = -1;
6806 for (
int j = 0; j < nliste; j++)
6808 const int num_facette = liste_facettes[j];
6809 const int compo = compo_connexe_facettes[num_facette];
6810 if (compo == compo_connexe_som)
6814 const double d = calculer_carre_distance_sommet_facette(coord, mesh, num_facette);
6818 idx_facette = num_facette;
6822 distance[i] = std::min(distance[i], sqrt(dmin));
6823 fill_relative_velocity(vinterp_tmp,
vinterp_, facettes, idx_facette, i, vr_to_other);
6830static void check_neighbouring_layer_in_one_direction(
int dir0,
int dir1,
int dir2,
int n,
6831 const Int3& nb_elem_loc,
const Int3& ijk,
6832 const std::map<std::array<int,3>, std::set<int>>& bary_ijk_loc,
const DoubleTab& bary,
6833 const ArrOfInt& compo_connexe_facettes,
const int compo_connexe_som,
6834 const double x,
const double y,
const double z,
6835 double& distance,
int& id_facette )
6838 for(
int sens=0; sens<2; sens++)
6841 int a = sens == 0 ? std::max(ijk[dir0]-n,zero) : std::min(ijk[dir0]+n,nb_elem_loc[dir0]);
6842 for(
int b=std::max(ijk[dir1]-n,zero); b<std::min(ijk[dir1]+n,nb_elem_loc[dir1]); b++)
6843 for(
int c=std::max(ijk[dir2]-n,zero); c<std::min(ijk[dir2]+n,nb_elem_loc[dir2]); c++)
6845 std::array<int,3> current_ijk;
6849 if(bary_ijk_loc.find(current_ijk)!=bary_ijk_loc.end())
6851 std::set<int> bary_in_current_ijk = bary_ijk_loc.at(current_ijk);
6852 for(
const auto b_fa7: bary_in_current_ijk)
6854 if(compo_connexe_facettes[b_fa7] != compo_connexe_som)
6857 double dist = (bary(b_fa7,0)-x)*(bary(b_fa7,0)-x) + (bary(b_fa7,1)-y)*(bary(b_fa7,1)-y) + (bary(b_fa7,2)-z)*(bary(b_fa7,2)-z);
6858 distance = std::min(sqrt(dist),distance);
6886 const ArrOfInt& compo_connexe_sommets,
6887 const DoubleTab& vinterp_tmp,
6889 ArrOfDouble& distance,
6890 DoubleTab& vr_to_other,
6891 const double distmax)
6893 const IntTab& facettes = mesh.
facettes();
6895 const int nb_som = sommets_a_tester.
dimension(0);
6896 const int nb_fa7 = facettes.
dimension(0);
6899 Int3 ijk_glob, ijk_loc, useless;
6901 for(
int dir=0; dir<3; dir++)
6904 distance.resize_array(nb_som, RESIZE_OPTIONS::NOCOPY_NOINIT);
6906 vr_to_other.
resize(nb_som, 3, RESIZE_OPTIONS::NOCOPY_NOINIT);
6908 assert(vinterp_tmp.
dimension(0) == nb_som);
6910 DoubleTab bary(nb_fa7,3);
6912 std::map<std::array<int,3>, std::set<int>> bary_ijk_loc;
6913 for (
int fa7=0; fa7<nb_fa7; fa7++)
6916 int s0 = facettes(fa7,0), s1 = facettes(fa7,1), s2 = facettes(fa7,2);
6917 for(
int dir=0; dir<3; dir++)
6918 bary(fa7,dir) = (sommets_a_tester(s0,dir) + sommets_a_tester(s1,dir) + sommets_a_tester(s2,dir)) / 3.;
6920 splitting.
search_elem(bary(fa7,0), bary(fa7,1), bary(fa7,2), ijk_glob, ijk_loc, useless);
6921 std::array<int,3> ijk;
6922 for(
int dir=0; dir<3; dir++) ijk[dir] = ijk_loc[dir];
6923 bary_ijk_loc[ijk].insert(fa7);
6926 for (
int som=0; som<nb_som; som++)
6928 double x=sommets_a_tester(som,0), y=sommets_a_tester(som,1), z=sommets_a_tester(som,2);
6929 splitting.
search_elem(x, y, z, ijk_glob, ijk_loc, useless);
6930 const int compo_connexe_som = compo_connexe_sommets[som];
6932 bool local_domain_checked =
false;
6934 int id_facette = -1;
6935 double& dist = distance[som];
6938 while(!local_domain_checked && id_facette==-1)
6941 check_neighbouring_layer_in_one_direction(0 , 1, 2, n_layer,
6942 nb_elem_loc, ijk_loc,
6944 compo_connexe_facettes, compo_connexe_som,
6949 check_neighbouring_layer_in_one_direction(1 , 0, 2, n_layer,
6950 nb_elem_loc, ijk_loc,
6952 compo_connexe_facettes, compo_connexe_som,
6959 check_neighbouring_layer_in_one_direction(2 , 0, 1, n_layer,
6960 nb_elem_loc, ijk_loc,
6962 compo_connexe_facettes, compo_connexe_som,
6968 int i=ijk_loc[0], j=ijk_loc[1], k=ijk_loc[2];
6969 local_domain_checked = i-n_layer <= 0 && i+n_layer>= nb_elem_loc[0]
6970 && j-n_layer <= 0 && j+n_layer>= nb_elem_loc[1]
6971 && k-n_layer <= 0 && k+n_layer>= nb_elem_loc[2];
6978 fill_relative_velocity(vinterp_tmp,
vinterp_, facettes, id_facette, som, vr_to_other);
7005 ArrOfInt& compo_sommet, ArrOfDouble& distance,
7006 DoubleTab& vr_to_other,
double distmax)
7008 const Domaine_IJK& splitting = ref_domaine_.valeur();
7023 double min_coord, max_coord;
7025 if (processor_at_left < 0)
7033 const double left_node = coord_nodes[offset];
7034 min_coord = left_node + distmax;
7035 if(!flags[processor_at_left])
7037 flags[processor_at_left] = 1;
7042 if (processor_at_right < 0)
7051 const double right_node = coord_nodes[i_last_node];
7052 max_coord = right_node - distmax;
7053 if(!flags[processor_at_right])
7055 flags[processor_at_right] = 1;
7062 ArrOfInt index_sent_to_left;
7063 ArrOfInt index_sent_to_right;
7064 if (processor_at_left >= 0 || processor_at_right >= 0)
7067 const int nb_som = coord_sommets.
dimension(0);
7068 for (
int i_som = 0; i_som < nb_som; i_som++)
7070 Vecteur3 coord(coord_sommets, i_som);
7071 if (coord[dir] < min_coord)
7077 schema.
send_buffer(processor_at_left) << compo_sommet[i_som]
7078 << coord_sommets(i_som, 0)
7079 << coord_sommets(i_som, 1)
7080 << coord_sommets(i_som, 2)
7081 << vinterp_tmp(i_som, 0)
7082 << vinterp_tmp(i_som, 1)
7083 << vinterp_tmp(i_som, 2);
7085 if (coord[dir] > max_coord)
7089 schema.
send_buffer(processor_at_right) << compo_sommet[i_som]
7090 << coord_sommets(i_som, 0)
7091 << coord_sommets(i_som, 1)
7092 << coord_sommets(i_som, 2)
7093 << vinterp_tmp(i_som, 0)
7094 << vinterp_tmp(i_som, 1)
7095 << vinterp_tmp(i_som, 2);
7102 const int init_count = coord_sommets.
dimension(0);
7104 int count_right = 0;
7105 for (
int left_right = 0; left_right < 2; left_right++)
7107 const int pe_voisin = (left_right==0) ? processor_at_left : processor_at_right;
7108 int& count = (left_right==0) ? count_left : count_right;
7118 double x, y, z, vx, vy, vz;
7119 buf >> x >> y >> z >> vx >> vy >> vz;
7137 int index = init_count;
7139 for (
int left_right = 0; left_right < 2; left_right++)
7141 const int pe_voisin = (left_right==0) ? processor_at_left : processor_at_right;
7142 const int count = (left_right==0) ? count_left : count_right;
7146 for (
int i = 0; i < count; i++)
7147 buf << distance[index+i]
7148 << vr_to_other(index+i,0)
7149 << vr_to_other(index+i,1)
7150 << vr_to_other(index+i,2);
7155 for (
int left_right = 0; left_right < 2; left_right++)
7157 int pe_voisin = (left_right==0) ? processor_at_left : processor_at_right;
7158 const ArrOfInt indices_sommets = (left_right==0) ? index_sent_to_left : index_sent_to_right;
7162 const int count = indices_sommets.
size_array();
7163 for (
int i = 0; i < count; i++)
7168 buf >> vx >> vy >> vz;
7169 int idx = indices_sommets[i];
7170 distance[idx] = std::min(distance[idx], d);
7171 vr_to_other(idx,0) = vx;
7172 vr_to_other(idx,1) = vy;
7173 vr_to_other(idx,2) = vz;
7180 coord_sommets.
resize(init_count, 3);
7181 vinterp_tmp.
resize(init_count, 3);
7184 distance.resize_array(init_count);
7185 vr_to_other.
resize(init_count, 3);
7194 DoubleTab& v_closer)
7198 ArrOfIntFT compo_connexe_sommets;
7200 DoubleTab tmp_sommets = mesh.
sommets();
7210 compo_connexe_sommets,
7223 f.
ouvrir(
"dump_facette.txt", ios::app);
7228 f.
ouvrir(
"dump_facette.txt");
7230 for (
int i = 0; i < 4; i++)
7232 int som = mesh.
facettes()(fa7, i%3);
7233 f << mesh.
sommets()(som, 0) <<
" " << mesh.
sommets()(som, 1) <<
" " << mesh.
sommets()(som, 2) << finl;
7237void dump_elem(
const Domaine_IJK& split,
int ii,
int jj,
int kk,
int append = 0)
7242 f.
ouvrir(
"dump_elem.txt", ios::app);
7247 f.
ouvrir(
"dump_elem.txt");
7254 for (
int i = 0; i < 3; i++)
7260 for (
int dir = 0; dir < 3; dir++)
7262 for (
int seg = 0; seg < 8; seg++)
7265 int pend = seg + (1 << dir);
7268 f << p[pstart & 1][0] <<
" " << p[(pstart >> 1) & 1][1] <<
" " << p[(pstart >> 2) & 1][2] << finl;
7269 f << p[pend & 1][0] <<
" " << p[(pend >> 1) & 1][1] <<
" " << p[(pend >> 2) & 1][2] << finl;
7276void dump_elem(
const Domaine_VF& domaine,
int elem,
int append = 0)
7281 f.
ouvrir(
"dump_elem.txt", ios::app);
7286 f.
ouvrir(
"dump_elem.txt");
7289 for (
int i = 0; i < 3; i++)
7291 p[0][i] = domaine.xv(domaine.elem_faces(elem, i), i);
7292 p[1][i] = domaine.xv(domaine.elem_faces(elem, i+3), i);
7294 for (
int dir = 0; dir < 3; dir++)
7296 for (
int seg = 0; seg < 8; seg++)
7299 int pend = seg + (1 << dir);
7302 f << p[pstart & 1][0] <<
" " << p[(pstart >> 1) & 1][1] <<
" " << p[(pstart >> 2) & 1][2] << finl;
7303 f << p[pend & 1][0] <<
" " << p[(pend >> 1) & 1][1] <<
" " << p[(pend >> 2) & 1][2] << finl;
7324 if (duplicatas_etaient_presents)
7336 ArrOfInt compo_new(nb_facettes);
7337 int n = search_connex_components_local_FT(maillage, compo_new);
7339 int nb_compo_tot=compute_global_connex_components_FT(maillage, compo_new, n);
7344 ArrOfInt rejeton(nb_rejeton);
7351 for (
int fa7 = 0; fa7 < nb_facettes; fa7++)
7352 rejetonArea(compo_std[fa7], compo_new[fa7])+=surface_facettes[fa7];
7362 for (
int ibul_new = 0 ; ibul_new<nb_compo_tot ; ibul_new++)
7363 if (rejetonArea(ibul, ibul_new)!=0.0 && rejetonArea(ibul, ibul_stay)<rejetonArea(ibul, ibul_new))
7365 ibul_stay=ibul_new ;
7370 for (
int ibul_new = 0 ; ibul_new<nb_compo_tot ; ibul_new++)
7371 if (rejetonArea(ibul, ibul_new)!=0.0 && ibul_new!=ibul_stay)
7372 for (
int ii = 0; ii<nb_rejeton ; ii++)
7373 if (rejeton[ii]==-1)
7375 rejeton[ii]=ibul_new;
7379 Cerr <<
"matrice de surface des compo" << finl;
7380 Cerr << rejetonArea << finl;
7381 Cerr <<
"Les rejetons a supprimer sont" << finl;
7382 Cerr << rejeton << finl;
7385 for (
int fa7 = 0; fa7 < nb_facettes; fa7++)
7386 for (
int ii = 0; ii<nb_rejeton ; ii++)
7387 if (compo_new[fa7] == rejeton[ii])
7395 if (duplicatas_etaient_presents)
7409 IJK_Field_vector3_double& rappel,
7410 const IJK_Field_vector3_double& vitesse,
7411 const IJK_Field_double& indic,
7412 const IJK_Field_double& indic_ft,
7413 const double coef_immo,
7415 const double current_time,
7416 const double coef_ammortissement,
7417 const double coef_rayon_force_rappel,
7418 const double integration_time,
7419 const double coef_mean_force,
7420 const double coef_force_time_n)
7439 const int nbulles_tot = nbulles_reelles + nbulles_ghost;
7443 ArrOfDouble surface_par_bulle;
7446 ArrOfIntFT compo_connex_som;
7449 DoubleTab vitesses_translation_bulles(nbulles_tot,3);
7457 vitesses_translation_bulles);
7462 ArrOfDouble volume_reel;
7466 const double deltat = 1.;
7468 for (
int idir=0; idir < 3; idir++)
7477 while (dx> 0.8*ldom)
7480 while (dx< -0.8*ldom)
7483 individual_forces(ib,idir) = coef_mean_force*
mean_force_(ib,idir);
7484 individual_forces(ib,idir) += (1.-coef_mean_force)*coef_immo*(dx);
7485 individual_forces(ib,idir) += (1.-coef_mean_force)*coef_ammortissement*vitesses_translation_bulles(ib,idir);
7486 individual_forces(ib,idir) += (1.-coef_mean_force)*coef_force_time_n*
force_time_n_(ib,idir);
7498 for (
int idir=0; idir < 3; idir++)
7502 mean_force_(ib,idir)+=individual_forces(ib,idir)*deltat;
7507 envoyer_broadcast(individual_forces, 0);
7533 int reset = (!
reprise_) && (tstep==0);
7534 IOS_OPEN_MODE mode = (reset) ? ios::out : ios::app;
7535 for (
int idir=0; idir < 3; idir++)
7537 snprintf(s, 1000,
"%s_bulles_external_force_every_%d.out", nomcas, idir);
7543 snprintf(s, 1000,
"# Individual forces applied inside chiv[bubble_i]=1.\n# value=1./Vol_bubble \\int F_j(bubble_i) dv");
7547 snprintf(s, 1000,
"%.16e ", current_time);
7551 snprintf(s, 1000,
"%.16e ", individual_forces(ib,idir));
7564 const IJK_Field_double& indic_ns,
7565 const IJK_Field_double& indic_ft,
7566 DoubleTab& individual_forces,
7567 const ArrOfDouble& volume_reel,
7568 const DoubleTab& position)
7573 decodeur_num_compo=NUM_COMPO_INVALID;
7579 Int3 ijk_global, ijk_local, ijk_processeur;
7583 double x = position(ib,0);
7584 double y = position(ib,1);
7585 double z = position(ib,2);
7586 ijk_processeur[0]=0;
7587 ijk_processeur[1]=0;
7588 ijk_processeur[2]=0;
7591 s_ft.
search_elem(x,y,z, ijk_global, ijk_local, ijk_processeur);
7593 if (ijk_processeur[0] * ijk_processeur[1] * ijk_processeur[2] == 1)
7598 const int icompo =
num_compo_[num_elem_zvdf];
7600 decodeur_num_compo[icompo] = ib;
7608 ArrOfInt list_continuous_phase;
7610 int phase_continue1=-1000;
7611 int phase_continue2=-1000;
7612 int phase_continue3=-1000;
7613 int phase_continue4=-1000;
7614 int phase_continue5=-1000;
7615 int phase_continue6=-1000;
7619 if (decodeur_num_compo[i] == NUM_COMPO_INVALID)
7622 if (phase_continue1<0)
7624 phase_continue1 = i;
7626 else if (phase_continue2<0)
7628 phase_continue2 = i;
7630 else if (phase_continue3<0)
7632 phase_continue3 = i;
7634 else if (phase_continue4<0)
7636 phase_continue4 = i;
7638 else if (phase_continue5<0)
7640 phase_continue5 = i;
7642 else if (phase_continue6<0)
7644 phase_continue6 = i;
7649 const int nb_parts_continuous = list_continuous_phase.
size_array();
7650 Cerr <<
"nb parts continuous phase : " << list_continuous_phase.
size_array()
7651 <<
" continuous phases are : " << phase_continue1 <<
" " << phase_continue2 <<
" "
7652 << phase_continue3 <<
" " << phase_continue4 <<
" " << phase_continue5 <<
" " << phase_continue6 << finl;
7653 Cerr <<
"GB-check nb_compo - nb_continuous phase - nb_bulles_tot = "
7658 if (nb_parts_continuous == 0)
7660 Cerr <<
"No continuous phase found. " << finl;
7676 integration_cells_per_bubble=0;
7677 integrated_forces =0.;
7678 for (
int idir=0; idir < 3; idir++)
7680 const int nx = rappel_ft[idir].ni();
7681 const int ny = rappel_ft[idir].nj();
7683 int nz = rappel_ft[idir].nk();
7685 (idir==DIRECTION_K))
7687 force_zero_on_walls(rappel_ft[DIRECTION_K]);
7690 if (kmin + nz == nktot)
7703 for (
int k=nzdeb; k < nz; k++)
7704 for (
int j=0; j< ny; j++)
7705 for (
int i=0; i < nx; i++)
7709 const int icompo =
num_compo_[num_elem_zvdf];
7710 bool continuous =
false;
7711 for (
int c=0; c<nb_parts_continuous; c++)
7713 if (icompo == list_continuous_phase[c])
7719 if ((icompo==-1) || (continuous) )
7721 rappel_ft[idir](i,j,k) = 0.;
7725 int id_bulle = decodeur_num_compo[icompo];
7732 const int ibulle_reelle = decoder_numero_bulle(-ighost);
7735 id_bulle = ibulle_reelle;
7740 assert(id_bulle>=0);
7742 const double f = individual_forces(id_bulle,idir);
7743 integrated_forces(id_bulle,idir) += vol*f;
7744 integration_cells_per_bubble(id_bulle,idir) +=1;
7751 const double f = individual_forces(id_bulle,idir);
7752 rappel_ft[idir](i,j,k) = f;
7759 for (
int idir=0; idir < 3; idir++)
7763 individual_forces(ib, idir) = integrated_forces(ib, idir) / (vol*integration_cells_per_bubble(ib,idir));
7771 const IJK_Field_double& indic,
7772 const DoubleTab& individual_forces,
7773 const ArrOfDouble& volume_reel,
7774 const DoubleTab& position,
7775 const double coef_rayon_force_rappel)
7778 noms_forces.dimensionner_force(3);
7779 for (
int idir=0; idir < 3; idir++)
7781 noms_forces[idir] =
"";
7784 double xir = position(ib,0);
7785 double yir = position(ib,1);
7786 double zir = position(ib,2);
7787 double r2 = coef_rayon_force_rappel*coef_rayon_force_rappel*pow((volume_reel[ib]*3)/(4.*M_PI), 2./3.);
7788 noms_forces[idir] +=
"+((";
7789 noms_forces[idir] +=
Nom(
"(X-(")+
Nom(xir,
"%g")+
Nom(
"))*(X-(")+
Nom(xir,
"%g")+
Nom(
"))");
7790 noms_forces[idir] +=
Nom(
"+(Y-(")+
Nom(yir,
"%g")+
Nom(
"))*(Y-(")+
Nom(yir,
"%g")+
Nom(
"))");
7791 noms_forces[idir] +=
Nom(
"+(Z-(")+
Nom(zir,
"%g")+
Nom(
"))*(Z-(")+
Nom(zir,
"%g")+
Nom(
"))");
7792 noms_forces[idir] +=
Nom(
"-")+
Nom(r2,
"%g");
7796 noms_forces[idir] +=
")_LT_0.)*(1.-ff)*((1.-ff)_GT_0.000001)*(" ;
7797 noms_forces[idir] +=
Nom(individual_forces(ib,idir),
"%g");
7798 noms_forces[idir] +=
")" ;
7802 for (
int idir=0; idir < 3; idir++)
7810 const int ibulle_reelle = decoder_numero_bulle(-ighost);
7815 double r2 = coef_rayon_force_rappel*coef_rayon_force_rappel*pow((volume_reel[ibg]*3)/(4.*M_PI), 2./3.);
7816 noms_forces[idir] +=
"+((";
7817 noms_forces[idir] +=
Nom(
"(X-(")+
Nom(xir,
"%g")+
Nom(
"))*(X-(")+
Nom(xir,
"%g")+
Nom(
"))");
7818 noms_forces[idir] +=
Nom(
"+(Y-(")+
Nom(yir,
"%g")+
Nom(
"))*(Y-(")+
Nom(yir,
"%g")+
Nom(
"))");
7819 noms_forces[idir] +=
Nom(
"+(Z-(")+
Nom(zir,
"%g")+
Nom(
"))*(Z-(")+
Nom(zir,
"%g")+
Nom(
"))");
7820 noms_forces[idir] +=
Nom(
"-")+
Nom(r2,
"%g");
7823 noms_forces[idir] +=
")_LT_0.)*(1.-ff)*((1.-ff)_GT_0.000001)*(" ;
7824 noms_forces[idir] +=
Nom(individual_forces(ibulle_reelle,idir),
"%g");
7825 noms_forces[idir] +=
")" ;
7827 set_field_data(rappel[idir], noms_forces[idir], indic, 0. );
7832 const IJK_Field_double& indic,
7833 const ArrOfDouble& volume_reel,
7834 const DoubleTab& position)
const
7836 Nom nom_indicatrices_np;
7837 nom_indicatrices_np =
"";
7844 double r=pow((volume_reel[ib]*3)/(4.*M_PI), 1./3.);
7852 nom_indicatrices_np +=
"+((";
7853 nom_indicatrices_np +=
Nom(
"((X-(")+
Nom(x0,
"%g")+
Nom(
"))*(X-(")+
Nom(x0,
"%g")+
Nom(
"))/(")+
Nom(a,
"%g")+
Nom(
"*")+
Nom(a,
"%g")+
Nom(
"))");
7854 nom_indicatrices_np +=
Nom(
"+((Y-(")+
Nom(y0,
"%g")+
Nom(
"))*(Y-(")+
Nom(y0,
"%g")+
Nom(
"))/(")+
Nom(b,
"%g")+
Nom(
"*")+
Nom(b,
"%g")+
Nom(
"))");
7855 nom_indicatrices_np +=
Nom(
"+((Z-(")+
Nom(z0,
"%g")+
Nom(
"))*(Z-(")+
Nom(z0,
"%g")+
Nom(
"))/(")+
Nom(c,
"%g")+
Nom(
"*")+
Nom(c,
"%g")+
Nom(
"))");
7856 nom_indicatrices_np +=
Nom(
"-1.0");
7858 nom_indicatrices_np +=
")_LT_0.)*(1.0)" ;
7869 double r=pow((volume_reel[ibg]*3)/(4.*M_PI), 1./3.);
7877 nom_indicatrices_np +=
"+((";
7878 nom_indicatrices_np +=
Nom(
"((X-(")+
Nom(x0,
"%g")+
Nom(
"))*(X-(")+
Nom(x0,
"%g")+
Nom(
"))/(")+
Nom(a,
"%g")+
Nom(
"*")+
Nom(a,
"%g")+
Nom(
"))");
7879 nom_indicatrices_np +=
Nom(
"+((Y-(")+
Nom(y0,
"%g")+
Nom(
"))*(Y-(")+
Nom(y0,
"%g")+
Nom(
"))/(")+
Nom(b,
"%g")+
Nom(
"*")+
Nom(b,
"%g")+
Nom(
"))");
7880 nom_indicatrices_np +=
Nom(
"+((Z-(")+
Nom(z0,
"%g")+
Nom(
"))*(Z-(")+
Nom(z0,
"%g")+
Nom(
"))/(")+
Nom(c,
"%g")+
Nom(
"*")+
Nom(c,
"%g")+
Nom(
"))");
7881 nom_indicatrices_np +=
Nom(
"-1.0");
7883 nom_indicatrices_np +=
")_LT_0.)*(1.0)" ;
7886 Cerr <<
"Setting indicatrice_non_perturbe :" << nom_indicatrices_np << finl;
7888 set_field_data(indic_np, nom_indicatrices_np, indic, 0. );
7897 const DoubleTab& gravite,
7898 const double delta_rho,
7902 const bool parcourir
7911 statistics().create_custom_counter(
"Calcul Indicatrice Next",2,
"IJK");
7912 statistics().begin_count(
"Calcul Indicatrice Next",statistics().get_last_opened_counter_level()+1);
7923 statistics().end_count(
"Calcul Indicatrice Next");
7997 for (
int c=0; c < 3; c++)
8020 for (
int c=0; c<3; c++)
8032 for (
int bary_compo = 0; bary_compo < 3; bary_compo++)
8041 for (
int face_dir = 0; face_dir < 3; face_dir++)
8043 for (
int bary_compo = 0; bary_compo < 2; bary_compo++)
8053 for (
int d = 0; d < 3; d++)
8064 for (
int face_dir = 0; face_dir < 3; face_dir++)
8066 for (
int bary_compo = 0; bary_compo < 2; bary_compo++)
8072 statistics().end_count(
"Calcul Indicatrice Next");
8079 IJK_Field_double& indicatrice_intermediaire_ft,
8080 IJK_Field_double& indicatrice_intermediaire_ns,
8081 IJK_Field_vector3_double& indicatrice_surfacique_intermediaire_face_ft,
8082 IJK_Field_vector3_double& indicatrice_surfacique_intermediaire_face_ns,
8083 const bool parcourir
8089 statistics().create_custom_counter(
"Calcul Indicatrice Next",2,
"IJK");
8090 statistics().begin_count(
"Calcul Indicatrice Next",statistics().get_last_opened_counter_level()+1);
8094 indicatrice_intermediaire_ft.
data() = 1.;
8095 indicatrice_intermediaire_ns.
data() = 1.;
8101 statistics().end_count(
"Calcul Indicatrice Next");
8119 ref_ijk_ft_->eq_ns().redistrib_from_ft_elem().redistribute(
8120 indicatrice_intermediaire_ft, indicatrice_intermediaire_ns);
8127 ref_ijk_ft_->eq_ns().get_redistribute_from_splitting_ft_faces(
8128 indicatrice_surfacique_intermediaire_face_ft,
8129 indicatrice_surfacique_intermediaire_face_ns);
8132 statistics().end_count(
"Calcul Indicatrice Next");
8139 indicatrice_ft_test_.echange_espace_virtuel(indicatrice_ft_test_.ghost());
8144 groups_indicatrice_ft_test_.echange_espace_virtuel();
8146 SChaine group_indic;
8153 for (
int k = 0; k < nk; k++)
8154 for (
int j = 0; j < nj; j++)
8155 for (
int i = 0; i < ni; i++)
8160 indic <<
"(" << i <<
"," << j <<
"," << k <<
") : " <<
indicatrice_ft_[
next()](i, j, k) <<
" VS "
8161 << indicatrice_ft_test_(i, j, k) << finl;
8163 for (
int igroup = 0; igroup < nb_grps; igroup++)
8168 group_indic <<
"groupe = " << igroup <<
" (" << i <<
"," << j <<
"," << k
8170 << groups_indicatrice_ft_next_test_[igroup](i, j, k) << finl;
8178 Cerr <<
"Probleme_FTD_IJK_base:: calcul de l'indicatrice faux ! (iteration " << tstep_ <<
" sur " << nb_timesteps_ <<
" )"
8187 int current_time = next_time ?
next() :
old();
8188 int other_time = next_time ?
old() :
next();
8190 double current_indicatrice = next_time ?
In(i,j,k) :
I(i,j,k);
8191 double other_time_indicatrice = next_time ?
I(i,j,k) :
In(i,j,k);
8193 if ((
I(i,j,k) == 0.) && (
In(i,j,k) == 0.))
8197 else if ((
I(i,j,k) == 1.) && (
In(i,j,k) == 1.))
8201 else if (current_indicatrice == 0.)
8211 if (other_time_bary_compo > .5)
8221 else if (current_indicatrice == 1.)
8231 other_time_bary_compo =
opposing_barycentre(other_time_bary_compo, other_time_indicatrice);
8232 if (other_time_bary_compo > .5)
8250 assert(bary >= 0 && bary <= 1.);
8256 int current_time = next_time ?
next() :
old();
8257 int other_time = next_time ?
old() :
next();
8261 double current_indicatrice_surfacique = next_time ? next_indicatrice_surfacique : old_indicatrice_surfacique;
8262 double other_time_indicatrice_surfacique = next_time ? old_indicatrice_surfacique : next_indicatrice_surfacique;
8264 if (face_dir == bary_compo)
8273 int compo2D = (face_dir == 0) ? ((bary_compo == 2) ? 0 : ((bary_compo == 1) ? 1 : -1)) :
8274 ((face_dir == 1) ? ((bary_compo == 0) ? 0 : ((bary_compo == 2) ? 1 : -1)) :
8275 ((face_dir == 2) ? ((bary_compo == 1) ? 0 : ((bary_compo == 0) ? 1 : -1)) :
8277 assert(compo2D >= 0);
8278 if ((old_indicatrice_surfacique == 0.) && (next_indicatrice_surfacique == 0.))
8282 else if ((old_indicatrice_surfacique == 1.) && (next_indicatrice_surfacique == 1.))
8286 else if (current_indicatrice_surfacique == 0.)
8296 if (other_time_bary_compo > .5)
8306 else if (current_indicatrice_surfacique == 1.)
8316 other_time_bary_compo =
opposing_barycentre(other_time_bary_compo, other_time_indicatrice_surfacique);
8317 if (other_time_bary_compo > .5)
8335 assert(bary >= 0 && bary <= 1.);
8354 champs_compris_.switch_ft_fields();
8369 for (
int bary_compo = 0; bary_compo < 3; bary_compo++)
8378 for (
int face_dir = 0; face_dir < 3; face_dir++)
8380 for (
int bary_compo = 0; bary_compo < 2; bary_compo++)
8394 for (
int c = 0; c < 3; c++)
8410 const DoubleTab& gravite,
8411 const double delta_rho,
8419 ArrOfDouble potentiels_sommets;
8420 ArrOfDouble repulsions_sommets;
8423 potentiels_sommets, phi_par_compo);
8431 ArrOfDouble interfacial_source_term_sommet_dir, unite;
8438 for (
int dir = 0; dir < 3; dir++)
8441 interfacial_source_term_sommet_dir(som)=interfacial_source_term_sommet(som, dir);
8447 repulsions_sommets, repuls_par_compo);
8454 for (
int k = 0; k < nk; k++)
8455 for (
int j = 0; j < nj; j++)
8456 for (
int i = 0; i < ni; i++)
8464 ArrOfDouble& potentiels_sommets,
8465 ArrOfDouble& repulsions_sommets,
8466 const DoubleTab& gravite,
8467 const double delta_rho,
8482 potentiels_sommets = courbure;
8486 const double vol_cell = dxi * dxj * dxk;
8487 const DoubleTab& sommets = mesh.
sommets();
8488 const int nb_som = potentiels_sommets.
size_array();
8494 const ArrOfDouble& sigma_sommets =
maillage_ft_ijk_.Surfactant_facettes().get_sigma_sommets();
8495 for (
int i = 0; i < nb_som; i++)
8497 potentiels_sommets[i] *= sigma_sommets[i];
8502 potentiels_sommets *= sigma;
8513 ArrOfIntFT compo_connexe_sommets;
8518 DoubleTab deplacement;
8520 calculer_deplacement_from_code_compo_connexe_negatif(
8524 for (
int i = 0; i < nb_som; i++)
8526 double correction_potentiel_deplacement = 0.;
8527 if (compo_connexe_sommets[i] < 0)
8530 correction_potentiel_deplacement = -(
8531 gravite(0,0) * deplacement(i, 0) +
8532 gravite(0,1) * deplacement(i, 1) +
8533 gravite(0,2) * deplacement(i, 2)
8538 gravite(0,0) * coord[0] +
8539 gravite(0,1) * coord[1] +
8540 gravite(0,2) * coord[2] +
8541 correction_potentiel_deplacement
8543 potentiels_sommets[i] -= p * delta_rho;
8550 double vrx = DMAXFLOAT;
8551 double vry = DMAXFLOAT;
8552 double vrz = DMAXFLOAT;
8553 DoubleTab vr_to_closer;
8559 const int nznodes = z_grid_nodes.
size_array();
8560 const double zmin = z_grid_nodes[0];
8561 const double zmax = z_grid_nodes[nznodes - 1];
8562 double zsmin = (zmax - zmin) / 2.;
8563 for (
int i = 0; i < nb_som; i++)
8569 vrx = vr_to_closer(i, 0);
8570 vry = vr_to_closer(i, 1);
8571 vrz = vr_to_closer(i, 2);
8578 const double z = sommets(i, 2);
8579 double dzs = std::min(z - zmin, zmax - z);
8580 zsmin = std::min(zsmin, dzs);
8595 d = std::min(d, dzs);
8600 potentiels_sommets[i] -= phi;
8601 repulsions_sommets[i] = -phi * vol_cell;
8603 double dmin_local = dmin;
8609 if (fabs(dmin_local - dmin) < 1e-12)
8617 Cerr <<
"Warning. There were equalities ("
8619 <<
") in dmin computation. "
8620 <<
"2 equalities seems possible if it's between two markers "
8621 "separated by a proc frontieer. "
8636 envoyer_broadcast(vrx, iproc);
8637 envoyer_broadcast(vry, iproc);
8638 envoyer_broadcast(vrz, iproc);
8646 double norm_ur = std::sqrt(vrx * vrx + vry * vry + vrz * vrz);
8647 const double dt = ref_ijk_ft_->schema_temps_ijk().get_timestep();
8648 double cfl = dt * std::min(std::min(std::fabs(vrx) / dx, std::fabs(vry) / dy), std::fabs(vrz) / dz);
8649 int reset = (!
reprise_) && (itstep == 0);
8651 "_dmin.out",
"tstep\ttime\t\tdmin\tzsmin\tur\tCFL_local", reset
8653 fic << itstep <<
" " << time <<
" " << dmin <<
" " << zsmin;
8654 fic <<
" " << norm_ur <<
" " << cfl << finl;
8659 potentiels_sommets *= dxi * dxj * dxk;
8662 double lg_euler = std::max(std::max(dxi, dxj), dxk);
8664 const double sqrt_3 = 1.7320508075688772;
8665 bool ok = (lg_lagrange > lg_euler * sqrt_3);
8666 Cerr <<
"Test de la taille de maille eulerienne : lg_lagrange=" << lg_lagrange <<
" > (lg_euler=" << lg_euler
8667 <<
")*1.7320... ? " << (ok ?
"ok" :
"ko") <<
" (ratio : " << lg_lagrange / (lg_euler * sqrt_3)
8668 << (ok ?
" >" :
" <") <<
"1 )" << finl;
const IJK_Interfaces & get_interfaces() const
Int3 get_ijk(int n) const
int get_n(int i, int j, int k) const
static void calcul_surface_interface_efficace(double timestep, const IJK_Field_vector3_double &velocity, const IJK_Field_double &old_indicatrice_ns, const IJK_Field_double &next_indicatrice_ns, const DoubleTabFT_cut_cell_vector3 &vitesse_deplacement_interface, const DoubleTabFT_cut_cell_vector3 &normale_deplacement_interface, DoubleTabFT_cut_cell_scalar &surface_efficace_interface)
static void calcul_surface_interface_efficace_initiale(int methode_explicite, const IJK_Field_double &old_indicatrice_ns, const IJK_Field_double &next_indicatrice_ns, const IJK_Field_double &surface_interface_ns_old, const IJK_Field_double &surface_interface_ns_next, const IJK_Field_vector3_double &normal_of_interf_ns_old, const IJK_Field_vector3_double &normal_of_interf_ns_next, DoubleTabFT_cut_cell_vector3 &normale_deplacement_interface, DoubleTabFT_cut_cell_scalar &surface_efficace_interface, DoubleTabFT_cut_cell_scalar &surface_efficace_interface_initial)
static void imprimer_informations_surface_efficace_interface(int verbosite_surface_efficace_interface, double timestep, const IJK_Field_vector3_double &velocity, const IJK_Field_double &old_indicatrice_ns, const IJK_Field_double &next_indicatrice_ns, const DoubleTabFT_cut_cell_scalar &surface_efficace_interface, const DoubleTabFT_cut_cell_scalar &surface_efficace_interface_initial, const DoubleTabFT_cut_cell_vector3 &normale_deplacement_interface, const DoubleTabFT_cut_cell_vector3 &vitesse_deplacement_interface)
static void calcul_delta_volume_theorique_bilan(int compo, const DoubleTab &bounding_box_bulles, double timestep, const IJK_Field_double &indicatrice_avant_deformation, const IJK_Field_double &indicatrice_apres_deformation, const IJK_Field_vector3_double &indicatrice_surfacique_efficace_deformation_face, const Cut_field_vector3_double &deformation_velocity, IJK_Field_double &delta_volume_theorique_bilan)
static void imprimer_informations_surface_efficace_face(int verbosite_surface_efficace_face, int iteration_solver_surface_efficace_face, double timestep, const Cut_field_vector3_double &velocity, const IJK_Field_double &old_indicatrice_ns, const IJK_Field_double &next_indicatrice_ns, const DoubleTabFT_cut_cell_vector3 &indicatrice_surfacique_efficace_face, const DoubleTabFT_cut_cell_vector3 &indicatrice_surfacique_efficace_face_initial)
static void calcul_vitesse_interface(const IJK_Field_vector3_double &velocity, const IJK_Field_double &old_indicatrice_ns, const IJK_Field_double &next_indicatrice_ns, const IJK_Field_vector3_double &barycentre_phase1_ns_old, const IJK_Field_vector3_double &barycentre_phase1_ns_next, DoubleTabFT_cut_cell_vector3 &coord_deplacement_interface, DoubleTabFT_cut_cell_vector3 &vitesse_deplacement_interface)
static void calcul_surface_face_efficace_initiale(int methode_explicite, const IJK_Field_vector3_double &old_indicatrice_surfacique_face_ns, const IJK_Field_vector3_double &next_indicatrice_surfacique_face_ns, DoubleTabFT_cut_cell_vector3 &indicatrice_surfacique_efficace_face, DoubleTabFT_cut_cell_vector3 &indicatrice_surfacique_efficace_face_initial)
static void calcul_surface_face_efficace_iteratif(int verbosite_surface_efficace_face, double timestep, const Cut_field_vector3_double &velocity, int &iteration_solver_surface_efficace_face, const IJK_Field_double &old_indicatrice_ns, const IJK_Field_double &next_indicatrice_ns, const IJK_Field_vector3_double &old_indicatrice_surfacique_face_ns, const IJK_Field_vector3_double &next_indicatrice_surfacique_face_ns, DoubleTabFT_cut_cell_vector3 &indicatrice_surfacique_efficace_face, const DoubleTabFT_cut_cell_vector3 &indicatrice_surfacique_efficace_face_initial, DoubleTabFT_cut_cell_vector6 &indicatrice_surfacique_efficace_face_correction, DoubleTabFT_cut_cell_scalar &indicatrice_surfacique_efficace_face_absolute_error)
void set_to_uniform_value(int valeur)
void collecter_espace_virtuel(ArrOfDouble &tab, MD_Vector_tools::Operations_echange op) const
void echange_espace_virtuel(ArrOfDouble &tab) const
virtual void creer_tableau_elements(Array_base &, RESIZE_OPTIONS opt=RESIZE_OPTIONS::COPY_INIT) const
creation d'un tableau parallele de valeurs aux elements.
This class encapsulates all the information related to the eulerian mesh for TrioIJK.
int get_offset_local(int direction) const
Returns the local offset in requested direction.
int get_neighbour_processor(int previous_or_next, int direction) const
Returns the index of the requested neighbour processor (-1 if no neighbour).
bool is_uniform(int direction) const
Method returns true if uniform in this direction.
int get_nb_elem_local(int direction) const
Returns the number of elements owned by this processor in the given direction.
bool get_periodic_flag(int direction) const
Method returns true if periodic in this direction.
double get_domain_length(int direction) const
Returns the length of the entire domain in requested direction.
double get_constant_delta(int direction) const
Returns the size of cells in a direction.
double get_origin(int direction) const
Returns the coordinate of the first node (global) of the mesh in the requested direction.
int get_nb_items_global(Localisation loc, int direction) const
Returns the number of local items (on this processor) for the given localisation in the requested dir...
Vecteur3 get_coords_of_dof(int i, int j, int k, Localisation loc) const
Determines the dof of an element along a localisation.
int convert_ijk_cell_to_packed(const FixedVector< int, 3 > ijk) const
Converts the ijk index of an element to a cell index.
const ArrOfDouble & get_node_coordinates(int direction) const
Returns an array with the coordinates of all nodes in the mesh in requested direction.
FixedVector< int, 3 > convert_packed_to_ijk_cell(int index) const
Convert the local index of an element to a vector with IJK indices.
Localisation
Localisation sub class.
void search_elem(const double &x, const double &y, const double &z, FixedVector< int, 3 > &ijk_global, FixedVector< int, 3 > &ijk_local, FixedVector< int, 3 > &ijk_me) const
Find the element which contains the item's coodirnates.
int elem_faces(int i, int j) const
renvoie le numero de le ieme face de la maille num_elem la facon dont ces faces sont numerotees est
int face_voisins(int num_face, int i) const
renvoie l'element voisin de numface dans la direction i.
classe Domaine_dis_base Cette classe est la base de la hierarchie des domaines discretisees.
const Domaine & domaine() const
Fichier en lecture Cette classe est a la classe C++ ifstream ce que la classe Entree est a la.
Ecriture dans un fichier partage Cette classe derive de Ecr_Fic_Par, en utilisant une sortie en binai...
int put(const unsigned *ob, std::streamsize n, std::streamsize pas) override
int ouvrir(const char *name, IOS_OPEN_MODE mode=ios::out) override
Ouvre le fichier avec les parametres mode et prot donnes Ces parametres sont les parametres de la met...
Sortie & syncfile() override
Provoque l'ecriture sur disque des donnees accumulees sur les differents processeurs depuis le dernie...
void set_64b(bool is64) override
ifstream & get_ifstream()
Class defining operators and methods for all reading operation in an input flow (file,...
classe Equation_base Le role d'une equation est le calcul d'un ou plusieurs champs....
Domaine_dis_base & domaine_dis()
Renvoie le domaine discretise associe a l'equation.
void passer_variable_intensive(const Maillage_FT_IJK &mesh)
ArrOfDouble check_conservation(const Maillage_FT_IJK &mesh)
void correction_conservation_globale(const Maillage_FT_IJK &mesh, const ArrOfDouble &surfactant_avant_remaillage, const ArrOfDouble &surfactant_apres_remaillage)
void passer_variable_extensive(const Maillage_FT_IJK &mesh)
void set_disable_surfactant(bool disable_surfactant)
ArrOfDouble get_FT_field_Array() const
void avancer_en_temps(const Maillage_FT_IJK &mesh, const double time_step)
bool get_disable_surfactant() const
void allocate(const Domaine_IJK &d, Domaine_IJK::Localisation l, int ghost_size, int additional_k_layers=0, int nb_compo=1, const Nom &name=Nom(), bool external_storage=false, int monofluide=0, double rov=0., double rol=0., int use_inv_rho_in_pressure_solver=0)
void echange_espace_virtuel(int ghost)
Exchange data over "ghost" number of cells.
Domaine_IJK::Localisation get_localisation() const
const Domaine_IJK & get_domaine() const
const Domaine_IJK & get_domaine() const
void echange_espace_virtuel()
void reset_flags_and_counters()
Nom lata_interfaces_meshname_
ArrOfInt ghost_compo_converter_
void calculer_kappa_ft(IJK_Field_double &kappa_ft)
void calculer_surface_interface(IJK_Field_double &surf_interface, IJK_Field_double &indic)
FixedVector< FixedVector< IJK_Field_double, max_authorized_nb_of_components_ >, 2 > phi_par_compo_
DoubleTab positions_reference_
ArrOfDoubleFT distance_autres_interfaces_
FixedVector< IJK_Field_vector3_double, 2 > indicatrice_surfacique_face_ft_
IJK_Field_vector3_double deformation_velocity_
static double opposing_barycentre(double initial_barycentre, double initial_area)
SurfaceVapeurIJKComputation surface_vapeur_par_face_computation_
int posttraiter_champs_instantanes(const Motcles &liste_post_instantanes, const char *lata_name, const int lata_step) const
double portee_force_repulsion_
Nom fichier_reprise_interface_
const Milieu_base & milieu() const override
void calculer_surfactant(ArrOfDouble &surfactant, ArrOfDouble &surfactant_min, ArrOfDouble &surfactant_max) const
FixedVector< FixedVector< IJK_Field_double, max_authorized_nb_of_components_ >, 2 > indicatrice_par_compo_
void compute_surface_average_per_bubble(const ArrOfDouble &surfaces, const ArrOfDouble &in, ArrOfDouble &out) const
void ajouter_terme_source_interfaces(IJK_Field_vector3_double &vpoint, IJK_Field_vector3_double &vrepul, IJK_Field_vector3_double &vabsrepul) const
int is_terme_gravite_rhog() const
Probleme_FTD_IJK_base & probleme_ijk()
FixedVector< IJK_Field_double, 2 > indicatrice_ns_
void initialize(const Domaine_IJK &splitting_FT, const Domaine_IJK &splitting_NS, const Domaine_dis_base &domaine_dis, const int thermal_probes_ghost_cells=0, const bool compute_vint=true, const bool is_switch=false)
DoubleTab bubbles_velocities_
Nom fichier_sauvegarde_interface_
void set_param_reprise_pb(Param &)
int get_nb_bulles_ghost(const int print=0) const
void deplacer_bulle_perio(const ArrOfInt &masque_deplacement_par_compo)
double factor_length_duplicata_
void calculer_var_volume_remaillage(double timestep, const DoubleTab &vitesses_translation_bulles, const DoubleTab &mean_bubble_rotation_vector, const DoubleTab ¢re_gravite, ArrOfDouble &var_volume)
void calculer_indicatrices(IJK_Field_vector3_double &indic)
void supprimer_certaines_bulles_reelles()
void posttraiter_tous_champs(Motcles &liste) const
FixedVector< FixedVector< IJK_Field_double, max_authorized_nb_of_components_ >, 2 > surface_par_compo_
void transporter_maillage_deformation(const int correction_semi_locale_volume_bulle, const DoubleTab &vitesses_translation_bulles, const DoubleTab &mean_bubble_rotation_vector, const DoubleTab ¢re_gravite, const double dt_tot, ArrOfDouble &dvol, const int rk_step, const int first_step_interface_smoothing=0)
double seuil_indicatrice_negligeable_
void calculer_indicatrice_optim(IJK_Field_double &indic)
DoubleTab bubbles_bary_old_
const IJK_Field_double & I_ft() const
void compute_drapeaux_vapeur_v3(const Maillage_FT_IJK &mesh, const Domaine_IJK &split, const IntVect &vecteur_composantes, ArrOfInt &drapeau_vapeur) const
DoubleTabFT_cut_cell_vector3 vitesse_deplacement_interface_
void calculer_color(ArrOfInt &color) const
void calculer_vitesse_de_deformation(int compo, const DoubleTab &bounding_box_bulles, const Cut_field_vector3_double &cut_field_velocity, const DoubleTab &vitesses_translation_bulles, const DoubleTab &mean_bubble_rotation_vector, const DoubleTab &positions_bulles)
int lire_motcle_non_standard(const Motcle &un_mot, Entree &is) override
Lecture des parametres de type non simple d'un objet_U a partir d'un flot d'entree.
int timestep_sauvegarde_interface_
FixedVector< IJK_Field_double, 2 > surface_interface_ns_
ArrOfDouble var_volume_correction_globale_
void compute_indicatrice_non_perturbe(IJK_Field_double &indic_np, const IJK_Field_double &indic, const ArrOfDouble &volume_reel, const DoubleTab &position) const
void associer_pb_base(const Probleme_base &) override
S'associe au Probleme passe en parametre.
FixedVector< IJK_Field_vector3_double, 2 > bary_of_interf_ns_
int get_ghost_number_from_compo(const int compo) const
void calculer_aire_interfaciale_for_compo(IJK_Field_double &ai, const int compo) const
bool has_readen_barycentres_prev_
void creer_duplicata_bulles()
DoubleTabFT_cut_cell_scalar surface_efficace_interface_
void calculer_vmoy_translation_composantes_connexes(const Maillage_FT_IJK &maillage, const ArrOfDouble &surface_facette, const ArrOfDouble &surface_par_bulle, const ArrOfInt &compo_connexes_facettes, const int nbulles_reelles, const int nbulles_ghost, const DoubleTab &vitesse_sommets, DoubleTab &vitesses_translation_sommets) const
const IJK_Field_double & I() const
void transporter_maillage_remaillage(int correction_semi_locale_volume_bulle, const DoubleTab &vitesses_translation_bulles, const DoubleTab &mean_bubble_rotation_vector, const DoubleTab ¢re_gravite, double dt_tot, ArrOfDouble &dvol, const int rk_step, const double temps)
int flag_positions_reference_
IJK_Field_double indicatrice_apres_remaillage_ft_
void update_indicatrice_variables_monofluides()
FixedVector< IJK_Field_vector3_double, 2 > indicatrice_surfacique_face_ns_
void compute_drapeaux_vapeur_v2(const IntVect &vecteur_composantes, ArrOfInt &drapeau_liquide) const
DoubleTabFT_cut_cell_vector3 indicatrice_surfacique_efficace_face_initial_
FixedVector< FixedVector< IJK_Field_double, 3 *max_authorized_nb_of_components_ >, 2 > bary_par_compo_
DoubleTab RK3_G_store_vi_
FixedVector< FixedVector< IJK_Field_double, 3 *max_authorized_nb_of_components_ >, 2 > normale_par_compo_
IJK_Field_vector3_double indicatrice_surfacique_avant_remaillage_face_ft_
IJK_Composantes_Connex ijk_compo_connex_
void convert_to_IntVect(const ArrOfInt &in, IntVect &out) const
DoubleTab bubbles_bary_new_
FixedVector< FixedVector< FixedVector< IJK_Field_double, 2 >, 3 >, 2 > barycentre_phase1_face_ns_
void dupliquer_bulle_perio(ArrOfInt &masque_duplicata_pour_compo)
ArrOfDouble var_volume_deformation_
void calculer_barycentre(IJK_Field_vector3_double &baric, IJK_Field_double &indic)
std::map< Motcle, IJK_Field_double > scalar_post_fields_
ArrOfDouble var_volume_remaillage_
Remaillage_FT_IJK remaillage_ft_ijk_
int get_nb_bulles_reelles() const
FixedVector< IJK_Field_int, 2 > nb_compo_traversante_
void compute_compo_connex_from_bounding_box()
void preparer_duplicata_bulles(const DoubleTab &bounding_box_of_bubbles, const DoubleTab &bounding_box_offsetp, const DoubleTab &bounding_box_offsetm, const DoubleTab &authorized_bounding_box, ArrOfInt &masque_duplicata_pour_compo_reel)
int ghost_compo_converter(const int i) const
int dt_impression_bilan_indicatrice_
bool read_bubbles_barycentres_vel(const Nom &interf_name, FixedVector< ArrOfDouble, 3 > &bubbles_rising_dir, FixedVector< ArrOfDouble, 3 > &bubbles_rising_vel, ArrOfDouble &bubbles_rising_vel_mag)
const IJK_Field_vector3_double & get_indicatrice_surfacique_face_old() const
bool has_champ(const Motcle &nom) const override
DoubleTab bubbles_velocities_bary_
FixedVector< IJK_Field_vector3_double, 2 > normal_of_interf_
Parcours_interface parcours_
void calculer_phi_repuls_sommet(ArrOfDouble &potentiels_sommets, ArrOfDouble &repulsions_sommets, const DoubleTab &gravite, const double delta_rho, const double sigma, const double time, const int itstep)
double portee_wall_repulsion_
FixedVector< IJK_Field_vector3_double, 2 > barycentre_phase1_ft_
DoubleTabFT_cut_cell_vector3 indicatrice_surfacique_efficace_face_
bool active_repulsion_paroi_
void postraiter_colors(Sortie &os, const double current_time) const
bool read_bubbles_barycentres(const Nom &interf_name, const Nom &suffix, FixedVector< ArrOfDouble, 3 > &bubbles_bary)
void get_noms_champs_postraitables(Noms &noms, Option opt=NONE) const override
FixedVector< FixedVector< FixedVector< IJK_Field_double, 2 >, 3 >, 2 > barycentre_phase1_face_ft_
void calculer_indicatrice_surfacique_barycentre_face(IJK_Field_vector3_double &indic_surfacique_face, FixedVector< FixedVector< IJK_Field_double, 2 >, 3 > &baric_face, IJK_Field_double &indic, IJK_Field_vector3_double &norme)
int get_recompute_indicator() const
void calculer_aire_interfaciale(IJK_Field_double &ai) const
FixedVector< int, 2 > n_faces_mouilles_
void calculer_indicatrice(IJK_Field_double &indic)
void update_old_intersections()
DoubleTab bounding_box_duplicate_criteria_
void detecter_et_supprimer_rejeton(bool duplicatas_etaient_presents)
FixedVector< FixedVector< IJK_Field_double, max_authorized_nb_of_components_ >, 2 > repuls_par_compo_
void transporter_maillage_rigide(const double dt_tot, const DoubleTab &vitesses_translation_bulles, const DoubleTab &mean_bubble_rotation_vector, const DoubleTab ¢re_gravite, const int rk_step, const int first_step_interface_smoothing=0)
int verbosite_surface_efficace_face_
IJK_Field_double indicatrice_avant_remaillage_ft_
ComputeValParCompoInCell val_par_compo_in_cell_computation_
int compute_cell_phase_with_interface_normal(int num_elem, int direction, int face_plus)
bool compute_distance_autres_interfaces_
FixedVector< IJK_Field_vector3_double, 2 > bary_of_interf_
FixedVector< IJK_Field_double, 2 > surface_interface_ft_
void calculer_distance_autres_compo_connexe_ijk(const DoubleTab &sommets_a_tester, const ArrOfInt &compo_connexe_sommets, const DoubleTab &vinterp_tmp, const Maillage_FT_IJK &mesh, ArrOfDouble &distance, DoubleTab &v_closer, const double distmax)
Maillage_FT_IJK maillage_ft_ijk_
static void Fill_postprocessable_fields(std::vector< FieldInfo_t > &chps)
const IJK_Field_vector3_double & get_indicatrice_surfacique_face_next() const
FixedVector< FixedVector< IJK_Field_int, max_authorized_nb_of_components_ >, 2 > compos_traversantes_
const int & nb_groups() const
int nb_compo_in_num_compo_
FixedVector< IJK_Field_vector3_double, 2 > barycentre_phase1_ns_
IJK_Field_double indicatrice_avant_remaillage_ns_
void store_bubbles_barycentres(const Nom &interf_name)
void calculer_phi_repuls_par_compo(FixedVector< IJK_Field_double, max_authorized_nb_of_components_ > &surf_par_compo, FixedVector< FixedVector< IJK_Field_double, max_authorized_nb_of_components_ >, 3 > &source_interf_par_compo, FixedVector< IJK_Field_double, max_authorized_nb_of_components_ > &phi_par_compo, FixedVector< IJK_Field_double, max_authorized_nb_of_components_ > &repuls_par_compo, const DoubleTab &gravite, const double delta_rho, const double sigma, const double time, const int itstep)
DoubleTabFT_cut_cell_scalar surface_efficace_interface_initial_
int verbosite_surface_efficace_interface_
IJK_Field_double field_repulsion_
const IJK_Field_double & In() const
void preparer_duplicata_bulles_masque_6bit(const DoubleTab &bounding_box, const DoubleTab &authorized_bounding_box, ArrOfInt &masque_duplicata_pour_compo)
static double mean_over_compo(const FixedVector< IJK_Field_double, max_authorized_nb_of_components_ > &field_for_compo, const IJK_Field_int &nb_compo_traversante, const int i, const int j, const int k)
void calcul_surface_efficace_face_initial(TYPE_SURFACE_EFFICACE_FACE type_surface_efficace_face)
void dumplata_ft_mesh(const char *filename, const char *meshname, int step) const
void compute_external_forces_(IJK_Field_vector3_double &rappel_ft, IJK_Field_vector3_double &rappel, const IJK_Field_vector3_double &vitesse, const IJK_Field_double &indic_ns, const IJK_Field_double &indic_ft, const double coef_immo, const int tstep, const double current_time, const double coef_ammortissement, const double coef_rayon_force_rappel, double compteur, double coef_mean_force, double coef_force_time_n)
void associer_switch(const Switch_FT_double &ijk_ft_switch)
void calculer_vecteurs_de_deplacement_rigide(DoubleTab &vitesses_translation_bulles, DoubleTab &mean_bubble_rotation_vector, DoubleTab ¢re_gravite, const int first_step_interface_smoothing=0)
void compute_rising_velocities_from_compo()
int disable_rigid_translation_
FixedVector< IJK_Field_vector3_double, 2 > surface_vapeur_par_face_ns_
FixedVector< IJK_Field_vector< double, max_authorized_nb_of_groups_ >, 2 > groups_indicatrice_ft_
void calcul_surface_efficace_interface_initial(TYPE_SURFACE_EFFICACE_INTERFACE type_surface_efficace_interface)
FixedVector< FixedVector< FixedVector< IJK_Field_double, max_authorized_nb_of_components_ >, 3 >, 2 > source_interf_par_compo_
FixedVector< FixedVector< IJK_Field_double, max_authorized_nb_of_components_ >, 2 > surf_par_compo_
bool use_tryggvason_interfacial_source_
DoubleTabFT_cut_cell_scalar indicatrice_surfacique_efficace_face_absolute_error_
void calcul_surface_efficace_face(TYPE_SURFACE_EFFICACE_FACE type_surface_efficace_face, double timestep, const Cut_field_vector3_double &total_velocity)
int disable_rigid_rotation_
bool read_barycentres_velocity_
void switch_indicatrice_next_old()
int timestep_reprise_interface_
double get_barycentre_face(bool next_time, int face_dir, int bary_compo, int phase, int i, int j, int k) const
bool use_barycentres_velocity_
DoubleTab bounding_box_delete_criteria_
FixedVector< FixedVector< IJK_Field_double, max_authorized_nb_of_components_ >, 2 > courbure_par_compo_
void calculer_indicatrices_optim(IJK_Field_vector3_double &indic)
void calculer_indicatrice_surfacique_face(IJK_Field_vector3_double &indic_surfacique_face, IJK_Field_double &indic, IJK_Field_vector3_double &norme)
ArrOfDouble bubbles_velocities_bary_magnitude_
Connectivite_frontieres connectivite_frontieres_
void calculer_vmoy_rotation_composantes_connexes(const Maillage_FT_IJK &maillage, const ArrOfDouble &surface_facette, const ArrOfDouble &surface_par_bulle, const ArrOfInt &compo_connexes_facettes, const int nbulles_reelles, const int nbulles_ghost, const DoubleTab ¢re_gravite, const DoubleTab &vitesse_sommets, const DoubleTab &vitesse_translation_sommets, DoubleTab &mean_bubble_rotation_vector) const
void calculer_distance_autres_compo_connexe2(ArrOfDouble &distance, DoubleTab &v_closer)
bool has_computed_bubble_barycentres_
void compute_compo_connex_from_interface()
void transferer_bulle_perio()
Intersection_Interface_ijk_cell intersection_ijk_cell_
bool has_champ_vectoriel(const Motcle &nom) const override
void parcourir_maillage()
void calculer_indicatrice_intermediaire(IJK_Field_double &indicatrice_intermediaire_ft_, IJK_Field_double &indicatrice_intermediaire_ns_, IJK_Field_vector3_double &indicatrice_surfacique_intermediaire_face_ft_, IJK_Field_vector3_double &indicatrice_surfacique_intermediaire_face_ns_, const bool parcourir=true)
void calculer_aspect_ratio(ArrOfDouble &aspect_ratio) const
bool correction_gradient_potentiel_
static int est_pure(double indicatrice)
double get_barycentre(bool next_time, int bary_compo, int phase, int i, int j, int k) const
void calcul_surface_efficace_interface(TYPE_SURFACE_EFFICACE_INTERFACE type_surface_efficace_interface, double timestep, const Cut_field_vector3_double &velocity)
void compute_external_forces_color_function(IJK_Field_vector3_double &rappel, const IJK_Field_double &indic_ns, const IJK_Field_double &indic_ft, DoubleTab &individual_forces, const ArrOfDouble &volume_reel, const DoubleTab &position)
Intersection_Interface_ijk_face intersection_ijk_face_
FixedVector< IJK_Field_vector3_double, 2 > surface_vapeur_par_face_
void update_surface_normale() const
IJK_Field_double indicatrice_apres_remaillage_ns_
void calculer_bounding_box_bulles(DoubleTab &bounding_box, int option_shear=0) const
void remailler_interface(const double temps, Maillage_FT_IJK &maillage, ArrOfDouble &var_volume, Remaillage_FT_IJK &algo_remaillage_local)
void calculer_poussee_bulles(const DoubleTab &gravite, DoubleTab &poussee) const
const IJK_Field_double & get_IJK_field(const Motcle &nom) override
double delta_p_max_repulsion_
FixedVector< FixedVector< IJK_Field_vector3_double, 3 >, 2 > barycentre_vapeur_par_face_
void recursive_calcul_distance_chez_voisin(DoubleTab &vinterp_tmp, int dir, const Maillage_FT_IJK &mesh, DoubleTab &coord_sommets, ArrOfInt &compo_sommet, ArrOfDouble &distance, DoubleTab &v_closer, double distmax)
void supprimer_duplicata_bulles()
DoubleTab bubbles_rising_vectors_bary_
void calculer_volume_bulles(ArrOfDouble &volumes, DoubleTab ¢re_gravite) const
FixedVector< IJK_Field_vector3_double, 2 > normal_of_interf_ns_
DoubleTab bounding_box_forbidden_criteria_
void compute_bubbles_volume_and_barycentres(ArrOfDouble &volumes, DoubleTab &barycentres, const int &store_values)
void compute_drapeaux_vapeur_v4(const IntVect &vecteur_composantes, ArrOfInt &drapeau_vapeur) const
IJK_Field_double delta_volume_theorique_bilan_ns_
IJK_Field_vector3_double indicatrice_surfacique_efficace_deformation_face_
IJK_Field_vector3_double indicatrice_surfacique_avant_remaillage_face_ns_
void compute_external_forces_parser(IJK_Field_vector3_double &rappel, const IJK_Field_double &indic_ns, const DoubleTab &individual_forces, const ArrOfDouble &volume_reel, const DoubleTab &position, const double coef_rayon_force_rappel)
int nb_bulles_ghost_before_
FixedVector< IJK_Field_double, 2 > indicatrice_ft_
IJK_Field_vector3_double indicatrice_surfacique_apres_remaillage_face_ns_
void imprime_bilan_indicatrice()
void calculer_distance_autres_compo_connexe_octree(const DoubleTab &sommets_a_tester, const ArrOfInt &compo_connexe_sommets, const DoubleTab &vinterp_tmp, const Maillage_FT_IJK &mesh, ArrOfDouble &distance, DoubleTab &v_closer, const double distmax)
void read_bubbles_barycentres_old_new(const Nom &interf_name)
FixedVector< IJK_Field_vector< double, max_authorized_nb_of_groups_ >, 2 > groups_indicatrice_ns_
const IJK_Field_vector3_double & get_IJK_field_vector(const Motcle &nom) override
double delta_p_wall_max_repulsion_
Maillage_FT_IJK old_maillage_ft_ijk_
int update_indicatrice(IJK_Field_double &indic)
DoubleTabFT_cut_cell_vector3 coord_deplacement_interface_
void set_recompute_indicator(int i)
FixedVector< FixedVector< IJK_Field_vector3_double, 3 >, 2 > barycentre_vapeur_par_face_ns_
void sauvegarder_interfaces(const char *lata_name, const Nom &interf_name="??")
void calculer_indicatrice_next(const DoubleTab &gravite, const double delta_rho, const double sigma, const double time, const int itstep, const bool parcourir=true)
void calculer_normales_et_aires_interfaciales(IJK_Field_double &ai, IJK_Field_double &kappa_ai, IJK_Field_vector3_double &normale_cell, const int igroup) const
void calculer_surface_bulles(ArrOfDouble &surfaces) const
IJK_Field_vector3_double indicatrice_surfacique_apres_remaillage_face_ft_
DoubleTabFT_cut_cell_vector3 normale_deplacement_interface_
DoubleTab bounding_box_NS_domain_
DoubleTabFT_cut_cell_vector6 indicatrice_surfacique_efficace_face_correction_
static double Lx_for_shear_perio
static double shear_x_time_
int index_element_suivant_
double contrib_volume_phase1_
double fraction_surface_intersection_
double contrib_barycentre_phase1_[3]
double contrib_aire_faces_phase1_[3]
double contrib_barycentre_faces_phase1_[3][2]
int index_facette_suivante_
: class Intersections_Elem_Facettes
const ArrOfInt & index_elem() const
Renvoie un tableau de taille domaine.
const ArrOfInt & index_facette() const
Renvoie un tableau de taille "nombre de facettes de l'interface" pour un element 0 <= facette < nb_fa...
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 ...
void get_liste_facettes_traversantes(int num_element, ArrOfInt &liste_facettes) const
: class Maillage_FT_Disc Cette classe decrit un maillage:
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.
const ArrOfInt & sommet_num_owner() const
pour postraitement, renvoie le numero des sommets sur le PE proprietaire des sommets
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 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.
const Desc_Structure_FT & desc_sommets() const
renvoie le descripteur des sommets (espace_distant/virtuel)
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 facette_PE_owner(ArrOfInt &facette_pe) const
pour postraitement, remplit le tableau en parametre avec le numero du PE proprietaire de chaque facet...
int sommet_virtuel(int i) const
const ArrOfInt & sommet_PE_owner() const
pour postraitement, renvoie le numero du PE proprietaire des sommets
const Intersections_Elem_Facettes & intersections_elem_facettes() const
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...
const IntTab & facettes() const
renvoie le tableau des facettes (reelles et virtuelles) dimension(0) = nombre de facettes,
const ArrOfInt & compo_connexe_facettes() const
FT_Field & Surfactant_facettes_non_const()
const Domaine_IJK & get_domaine() const
DoubleTab update_sigma_and_interfacial_source_term_sommet(const Domaine_IJK &splitting, bool compute_interfacial_source, bool use_tryggvason_formulation, const double sigma_const=-1.)
const FT_Field & Surfactant_facettes() const
void transporter(const DoubleTab &deplacement) override
Deplace les sommets de l'interface d'un vecteur "deplacement" fourni, Change eventuellement les somme...
void nettoyer_maillage() override
Retire toutes les facettes virtuelles, toutes les facettes invalides (sommet0 == sommet1) et tous les...
void recopie(const Maillage_FT_Disc &source_mesh, Statut_Maillage niveau_copie) override
Recopie une partie du maillage source dans *this.
void update_gradient_laplacien_Surfactant()
void ajouter_maillage_IJK(const Maillage_FT_IJK &added_mesh)
void supprimer_facettes(const ArrOfInt &liste_facettes)
double minimum_longueur_arrete() const
void calculer_compo_connexe_sommets(ArrOfIntFT &compo_connexe_sommets) const
void initialize(const Domaine_IJK &, const Domaine_dis_base &, const Parcours_interface &, const bool use_tryggvason_interfacial_source=false)
void set_composante_connexe(const int i_facette, const int icompo)
classe Milieu_base Cette classe est la base de la hierarchie des milieux (physiques)
Une chaine de caractere (Nom) en majuscules.
Un tableau d'objets de la classe Motcle.
int contient_(const char *const ch) const
void get_redistribute_from_splitting_ft_faces(const IJK_Field_vector3_double &faces_ft, IJK_Field_vector3_double &faces_ns)
Redistribute_Field & redistrib_from_ft_elem()
class Nom Une chaine de caractere pour nommer les objets de TRUST
Un tableau de chaine de caracteres (VECT(Nom)).
const Nom & que_suis_je() const
renvoie la chaine identifiant la classe.
virtual Entree & readOn(Entree &)
Lecture d'un Objet_U sur un flot d'entree Methode a surcharger.
static const Nom & nom_du_cas()
Renvoie une reference constante vers le nom du cas.
virtual Sortie & printOn(Sortie &) const
Ecriture de l'objet sur un flot de sortie Methode a surcharger.
int_t search_elements_box(double xmin, double ymin, double zmin, double xmax, double ymax, double zmax, ArrOfInt_t &elements) const
cherche tous les elements ou points ayant potentiellement une intersection non vide avec la boite don...
void build_elements(const _TAB_TYPE_ &coords, const IntTab_t &elements, const double epsilon, const bool include_virtual)
Construit un octree a partir d'elements volumiques decrits par des ensembles de sommets.
static bool DISABLE_DIPHASIQUE
Helper class to factorize the readOn method of Objet_U classes.
void ajouter(const char *keyword, const int *value, Param::Nature nat=Param::OPTIONAL)
Register an integer parameter.
const Cut_field_vector3_double & get_cut_field_velocity() const
const Nom & le_nom() const override
Donne le nom de l'Objet_U Methode a surcharger : renvoie "neant" dans cette implementation.
classe Probleme_base C'est un Probleme_U qui n'est pas un couplage.
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 trustIdType mppartial_sum(trustIdType i)
Calul de la somme partielle de i sur les processeurs 0 a me()-1 (renvoie 0 sur le processeur 0).
static double mp_max(double)
static void mp_sum_for_each_item(TRUSTArray< _TYPE_ > &x, int n=-1)
static Sortie & Journal(int message_level=0)
Renvoie un objet statique de type Sortie qui sert de journal d'evenements.
static int nproc()
renvoie le nombre de processeurs dans le groupe courant Voir Comm_Group::nproc() et PE_Groups::curren...
static double mp_sum(double)
Calcule la somme de x sur tous les processeurs du groupe courant.
static void mp_min_for_each(T &arg1, T &arg2)
C++14 compatible mp_min_for_each: combine multiple mp_min calls into one collective operation.
static int me()
renvoie mon rang dans le groupe de communication courant.
static void exit(int exit_code=-1)
Routine de sortie de TRUST dans une region Kokkos.
static int je_suis_maitre()
renvoie 1 si on est sur le processeur maitre du groupe courant (c'est a dire me() == 0),...
static void mp_min_for_each_item(TRUSTArray< _TYPE_ > &x, int n=-1)
void redistribute(const IJK_Field_double &input_field, IJK_Field_double &output_field)
void remaillage_local_interface(double temps, Maillage_FT_IJK &maillage)
void barycentrer_lisser_systematique_ijk(Maillage_FT_IJK &maillage, ArrOfDouble &var_volume)
int a_remailler(double temps, const Maillage_FT_Disc &maillage) const
Cette classe derivee de Sortie empile ce qu'on lui envoie dans une chaine de caracteres.
const char * get_str() const
returns a copy of the string stored by the SChaine
Cette classe est a la classe C++ ofstream ce que la classe Sortie est a la classe C++ ostream Elle re...
void echange_taille_et_messages() const
Cette methode lance l'echange de donnees entre tous les processeurs.
Sortie & send_buffer(int num_PE) const
renvoie le buffer correspondant au processeur num_PE pour y entasser des donnees a envoyer.
void end_comm() const
Vide les buffers et libere les ressources: on a fini de lire les donnees recues dans les buffers.
Entree & recv_buffer(int num_PE) const
renvoie le buffer correspondant au processeur num_PE pour y lire les donnees recues.
void begin_comm() const
Reserve les buffers de comm pour une nouvelle communication.
void set_send_recv_pe_list(const ArrOfInt &send_pe_list, const ArrOfInt &recv_pe_list, const int me_to_me=0)
Definit la liste des processeurs a qui on va envoyer et de qui on va recevoir des donnees.
virtual int ouvrir(const char *name, IOS_OPEN_MODE mode=ios::out)
Classe de base des flux de sortie.
void append_array(_TYPE_ valeur)
_SIZE_ size_array() const
TRUSTArray & inject_array(const TRUSTArray &source, _SIZE_ nb_elements=-1, _SIZE_ first_element_dest=0, _SIZE_ first_element_source=0)
void resize_array(_SIZE_ new_size, RESIZE_OPTIONS opt=RESIZE_OPTIONS::COPY_INIT)
void resize(_SIZE_ new_size, RESIZE_OPTIONS opt=RESIZE_OPTIONS::COPY_INIT)
void resize(_SIZE_ n, RESIZE_OPTIONS opt=RESIZE_OPTIONS::COPY_INIT)
_SIZE_ dimension(int d) const
virtual void echange_espace_virtuel(IsExchangeBlocking exchange_type=IsExchangeBlocking::DefaultBlocking, const std::string kernel_name="noname")
static void transfert_conservatif_eulerien_vers_lagrangien_sommets(const Maillage_FT_Disc &maillage, const DoubleVect &valeurs_euler, ArrOfDouble &valeurs_lagrange)
static double produit_scalaire(const Vecteur3 &x, const Vecteur3 &y)
static void produit_vectoriel(const Vecteur3 &x, const Vecteur3 &y, Vecteur3 &resu)